[openfoammaofnepo]

Hello,

In the very latest openfoam e.g. 3.0.0, the multiphase solver folder has lots of update. In the twoPhaseEulerFoam solver, there are two sets of velocity and pressure source files, i.e. pU and pUf. In twoPhaseEulerFoam.C, they are chosed based on the switch "faceMomentum". I checked the source files for "pUf/UEqns.H" and "pU/UEqns.H". There are some differences, for instance, for the U equation, the time derivative term in pUf/UEqns.H is not considered. Could anybody familar this solver give some explanations about this? This is a little confused. Thank you.

[hillman860]

It is switch between face-based and cell-based formulation.
Have you seen the seen the comment in git?

https://github.com/OpenFOAM/OpenFOAM...777bb10a3e1dda

[openfoammaofnepo]

Dear Sankar,

I also saw that information last week. Thank you so much for your reply.

best regards,
OFFO

[openfoammaofnepo]

Hi,

In the solver twoPhaseEulerFoam in OF300, we have the following two lines in createFields.H:

volScalarField& alpha1 = phase1;
volScalarField& alpha2 = phase2;

I think alpha1 and alpha2 are the phase volume fractions for dispersed and continuous phases. However, I did not find, for phase1 and phase2 can return the volume fraction. phase1 and phase2 are the objects of phaseModel, which are defined in twophaseSystem. But I did not find their forms that shows returning the volume fraction. Anybody knows about this issue? Thank you so much.

OFFO

[GerhardHolzinger]

The phase model class is derived, among others, from volScalarField. I.e. the phaseModel is a scalar field. Guess which one?

The constructor of phaseModel tells us more:

Foam::phaseModel::phaseModel
(
    const phaseSystem& fluid,
    const word& phaseName,
    const label index
)
:
    volScalarField
    (
        IOobject
        (
            IOobject::groupName("alpha", phaseName),
            fluid.mesh().time().timeName(),
            fluid.mesh(),
            IOobject::READ_IF_PRESENT,
            IOobject::AUTO_WRITE
        ),
        fluid.mesh(),
        dimensionedScalar("alpha", dimless, 0)
    ),

    fluid_(fluid),
    // and so on ...

Thus, the code you posted makes perfect sense.

[openfoammaofnepo]

Originally Posted by GerhardHolzinger

The phase model class is derived, among others, from volScalarField. I.e. the phaseModel is a scalar field. Guess which one?

The constructor of phaseModel tells us more:

Foam::phaseModel::phaseModel
(
    const phaseSystem& fluid,
    const word& phaseName,
    const label index
)
:
    volScalarField
    (
        IOobject
        (
            IOobject::groupName("alpha", phaseName),
            fluid.mesh().time().timeName(),
            fluid.mesh(),
            IOobject::READ_IF_PRESENT,
            IOobject::AUTO_WRITE
        ),
        fluid.mesh(),
        dimensionedScalar("alpha", dimless, 0)
    ),

    fluid_(fluid),
    // and so on ...

Thus, the code you posted makes perfect sense.

Dear Gerhard,

Thank you so much for your. clear now!

[fs82]

Originally Posted by openfoammaofnepo

Hello,
the time derivative term in pUf/UEqns.H is not considered. Could anybody familar this solver give some explanations about this? This is a little confused. Thank you.

The reason for faceMomentum is not always obvious, but sometimes, you might have seen some striping or oscillations in your phase velocities. The reason for that is not fully understood, but is very similar to the pressure oscillations which are damped by the RhieChow interpolation. faceMomentum extends this idea to damp the oscillations in the velocities for multiphase flows. However, it is experimental and the whole theory is not clear to me. The experimental status is, to my opinion, also the reason why you miss the time derivative. If you check the solver you will find in pEqn.H

surfaceScalarField alphaRhof10
(
    "alphaRhof10",
    fvc::interpolate
    (
        max(alpha1.oldTime(), phase1.residualAlpha())
       *rho1.oldTime()
    )
);

surfaceScalarField alphaRhof20
(
    "alphaRhof20",
    fvc::interpolate
    (
        max(alpha2.oldTime(), phase2.residualAlpha())
       *rho2.oldTime()
    )
);

and a little bit later this lines

surfaceScalarField rAUf1
(
    IOobject::groupName("rAUf", phase1.name()),
    1.0
   /(
        (alphaRhof10 + Vmf)/runTime.deltaT()
      + fvc::interpolate(U1Eqn.A())
      + Kdf
   )
);

surfaceScalarField rAUf2
(
    IOobject::groupName("rAUf", phase2.name()),
    1.0
   /(
        (alphaRhof20 + Vmf)/runTime.deltaT()
      + fvc::interpolate(U2Eqn.A())
      + Kdf
   )
);

and finally

    // Phase-1 predicted flux
    surfaceScalarField phiHbyA1
    (
        IOobject::groupName("phiHbyA", phase1.name()),
        phi1
    );

    phiHbyA1 =
        rAUf1
       *(
            (alphaRhof10 + Vmf)
           *MRF.absolute(phi1.oldTime())/runTime.deltaT()
          + fvc::flux(U1Eqn.H())
          + Vmf*ddtPhi2
          + Kdf*MRF.absolute(phi2)
          - Ff1()
        );

    // Phase-2 predicted flux
    surfaceScalarField phiHbyA2
    (
        IOobject::groupName("phiHbyA", phase2.name()),
        phi2
    );

    phiHbyA2 =
        rAUf2
       *(
            (alphaRhof20 + Vmf)
           *MRF.absolute(phi2.oldTime())/runTime.deltaT()
          + fvc::flux(U2Eqn.H())
          + Vmf*ddtPhi1
          + Kdf*MRF.absolute(phi1)
          - Ff2()
       );

I marked the lines which are needed for the time derivative bold. If you write this lines (Vmf is not needed, it represents the contribution of the virtual mass) on a paper and transfer it to an equation you will end up with the time derivative

(rho^n*alpha^n*U^(n+1) - rho^n*alpha^n*U^n)/dT

The first (implicit) term of the time derivative goes to rAUf and the second (explicit) one to phiHByA. This works for Euler forward, Euler backward and may be for CrankNicolson, but not for Euler backward second order (backward). Here you see the experimental status of the solver and the option in fvSchemes is not recognized by the solver.

If somebody argue this terms exists in a similar fashion also in the cellCentered version, than it might be right, but there are tiny but important differences

volScalarField rAU1
(
    IOobject::groupName("rAU", phase1.name()),
    1.0
   /(
        U1Eqn.A()
      + max(phase1.residualAlpha() - alpha1, scalar(0))
       *rho1/runTime.deltaT()
    )
);
volScalarField rAU2
(
    IOobject::groupName("rAU", phase2.name()),
    1.0
   /(
        U2Eqn.A()
      + max(phase2.residualAlpha() - alpha2, scalar(0))
       *rho2/runTime.deltaT()
    )
);

and

    volVectorField HbyA1
    (
        IOobject::groupName("HbyA", phase1.name()),
        U1
    );
    HbyA1 =
        rAU1
       *(
            U1Eqn.H()
          + max(phase1.residualAlpha() - alpha1, scalar(0))
           *rho1*U1.oldTime()/runTime.deltaT()
        );

    volVectorField HbyA2
    (
        IOobject::groupName("HbyA", phase2.name()),
        U2
    );
    HbyA2 =
        rAU2
       *(
            U2Eqn.H()
         +  max(phase2.residualAlpha() - alpha2, scalar(0))
           *rho2*U2.oldTime()/runTime.deltaT()
        );

In the cellCentered version these terms are needed to avoid division by zero, or in other words zero equals zero in case a phase fraction gets zero. The important difference in the cellCentered Version is, that this term normally euqals zero (check the max condition) and only in the case of alpha < resisdualAlpha it gets larger than zero. In the faceMomentum version this term is always larger than zero. Hope this helps.

Best regards
Fabian

[openfoammaofnepo]

Hello,

Thank you so much for your reply.

In EEqns.H of reactingTwoPhaseEulerFoam solver, we have the following to update the thermodynamic properties:

fluid.correctThermo();

I am interested in looking at the source files of this function. So I checked the class of twoPhaseSystem.H and twoPhaseSystem.C, and actually this function can be found in phaseSystem.C, as follows:

void Foam::phaseSystem::correctThermo()
{
    forAll(phaseModels_, phasei)
    {
        phaseModels_[phasei].correctThermo();
    }
}

Here I am a little confused about phaseModels_, which is the object of the class :

        //- Phase models

        phaseModelList phaseModels_;

Now where can I find the function of correctThermo()? Anybody can give me some hints? Thank you so much in advance.

I am using OpenFOAM 3.0.1.

best regards, OFFO

[openfoammaofnepo]

For record:

template<class BasePhaseModel>
void Foam::MultiComponentPhaseModel<BasePhaseModel>::correctThermo()
{
    volScalarField Yt
    (
        IOobject
        (
            IOobject::groupName("Yt", this->name()),
            this->fluid().mesh().time().timeName(),
            this->fluid().mesh()
        ),
        this->fluid().mesh(),
        dimensionedScalar("zero", dimless, 0)
    );

    PtrList<volScalarField>& Yi = Y();

    forAll(Yi, i)
    {
        if (i != inertIndex_)
        {
            Yt += Yi[i];
        }
    }

    if (inertIndex_ != -1)
    {
        Yi[inertIndex_] = scalar(1) - Yt;
        Yi[inertIndex_].max(0);
    }
    else
    {
        forAll(Yi, i)
        {
            Yi[i] /= Yt;
            Yi[i].max(0);
        }
    }

    BasePhaseModel::correctThermo();
}

[openfoammaofnepo]

Hello,

In twoPhaseEulerFoam, I found that if I change "nAlphaCorr" for alpha.* in fvSolution file from 1 to 2 or any values > 1, then the computations crash. For my own case and for the tutorial /twoPhaseEulerFoam/RAS/fluidizedBed, this happens. Did you also have this similar problem? What causes this? Thanks!