CFD Online Discussion Forums - About the alphaEqn.H in interPhaseChangeFoam

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Running, Solving & CFD (https://www.cfd-online.com/Forums/openfoam-solving/)

- - About the alphaEqn.H in interPhaseChangeFoam (https://www.cfd-online.com/Forums/openfoam-solving/215998-about-alphaeqn-h-interphasechangefoam.html)

About the alphaEqn.H in interPhaseChangeFoam

Hello dear foamers,

I'm new to OpenFOAM. I read a tutorial (in Chinese) discussing the interPhaseChangeFoam recently, according to which the alpha equation is,

$\frac{\partial \alpha_l}{\partial t} + \nabla \cdot (\alpha_l U) + \nabla \cdot (\alpha_l \alpha_v U_r ) = \alpha_l \nabla \cdot U - \alpha_l (\frac{\dot{m}}{\rho_l} - \frac{\dot{m}}{\rho_v} )+ \frac{\dot{m}}{\rho_l}$ (1)

The third term is the "compression term" adopted by OpenFOAM, and the derivation seems reasonable there.

However, when I looked into the "alphaEqn.H" in the interPhaseChangeFoam directory (OpenFOAM-v1812) , the equation is defined as,

Code:

fvScalarMatrix alpha1Eqn

        (

            fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)

          + fv::gaussConvectionScheme<scalar>

            (

                mesh,

                phi,

                upwind<scalar>(mesh, phi)

            ).fvmDiv(phi, alpha1)

          - fvm::Sp(divU, alpha1)

         ==

            fvm::Sp(vDotvmcAlphal, alpha1)

          + vDotcAlphal

        );

The "compression term" is missed in the equation. Further, the term "fvm::Sp(vDotvmcAlphal, alpha1)" is exactly the last two terms appeared in the Eqn.(1) which are related to mass transfer.
(

Code:

Foam::phaseChangeTwoPhaseMixture::vDotAlphal() const

{volScalarField alphalCoeff(1.0/rho1() - alpha1_*(1.0/rho1() - 1.0/rho2()));

    Pair<tmp<volScalarField>> mDotAlphal = this->mDotAlphal();



    return Pair<tmp<volScalarField>>

    (

        alphalCoeff*mDotAlphal[0],

        alphalCoeff*mDotAlphal[1]

    );

)
So there seems to be an extra term "vDotcAlphal".

So here's my questions:
1. Is the "compression term" $\nabla \cdot (\alpha_l \alpha_v U_r )$ involved in this solver? Actually I find the term "phir" in the following part of the file, but I'm not sure what is happening.

2. What's the meaning of the term "vDotcAlphal"? Or is there something wrong with the equation (1)?

Any instruction is welcome, and thank you very much in advance.

Hello again dear foamers,

I searched the forums and found some similar questions. I figured out the second question by looking into the code, and following is the explanation. As for the first question, I found somewhere saying that the compression term is not involved in the "prediction step", and appears in the following steps, which is related to the algorithm of MULES(I haven't checked that yet). I hope this will help if someone has similar questions.

So, I was confused about the difference between mDot and mDotAlphal and the term vDotcAlphal in the alpha equation at first. Taking constant model (actually is a modified Lee model provided by interCondensatingEvaporatingFoam) for example, the complete form of mass transfer rate is provided by the function mDot. Here we note the terms defined by mDotAlphal of evaporation and condensation as S_e and S_c, respectively.
The definition of these two functions are:

Code:

Foam::Pair<Foam::tmp<Foam::volScalarField>>

Foam::temperaturePhaseChangeTwoPhaseMixtures::constant::mDot() const

{



    volScalarField limitedAlpha1

    (

        min(max(mixture_.alpha1(), scalar(0)), scalar(1))

    );



    volScalarField limitedAlpha2

    (

        min(max(mixture_.alpha2(), scalar(0)), scalar(1))

    );



    const volScalarField& T = mesh_.lookupObject<volScalarField>("T");



    const twoPhaseMixtureEThermo& thermo =

        refCast<const twoPhaseMixtureEThermo>

        (

            mesh_.lookupObject<basicThermo>(basicThermo::dictName)

        );



    const dimensionedScalar& TSat = thermo.TSat();



    const dimensionedScalar T0(dimTemperature, Zero);



    return Pair<tmp<volScalarField>>

    (

        coeffC_*mixture_.rho2()*limitedAlpha2*max(TSat - T.oldTime(), T0),

        coeffE_*mixture_.rho1()*limitedAlpha1*max(T.oldTime() - TSat, T0)

    );

}



Foam::Pair<Foam::tmp<Foam::volScalarField>>

Foam::temperaturePhaseChangeTwoPhaseMixtures::constant::mDotAlphal() const

{

    const volScalarField& T = mesh_.lookupObject<volScalarField>("T");



    const twoPhaseMixtureEThermo& thermo =

        refCast<const twoPhaseMixtureEThermo>

        (

            mesh_.lookupObject<basicThermo>(basicThermo::dictName)

        );



    const dimensionedScalar& TSat = thermo.TSat();



    const dimensionedScalar T0(dimTemperature, Zero);



    return Pair<tmp<volScalarField>>

    (

        coeffC_*mixture_.rho2()*max(TSat - T.oldTime(), T0),

       -coeffE_*mixture_.rho1()*max(T.oldTime() - TSat, T0)

    );

}

We can see the term mDot equals mDotAlphal multiply the volume fraction of the other phase.

And we can write the RHS of alpha1Eqn in alphaEqn.H as,
fvm::Sp(vDotvmcAlphal, alpha1) + vDotcAlphal
= $[\frac{1}{\rho_l} - \alpha_l (\frac{1}{\rho_l}-\frac{1}{\rho_g})]\cdot (S_g -S_l )\cdot \alpha_l +[\frac{1}{\rho_l} - \alpha_l (\frac{1}{\rho_l}-\frac{1}{\rho_g})]\cdot S_l$
= $[\frac{1}{\rho_l} - \alpha_l (\frac{1}{\rho_l}-\frac{1}{\rho_g})] \cdot [(S_g -S_l)\alpha_l +S_l]$
= $[\frac{1}{\rho_l} - \alpha_l (\frac{1}{\rho_l}-\frac{1}{\rho_g})] \cdot [S_g \alpha_l +S_l \alpha_g]$

Based on the difference between mDotAlphal(S_g or S_l) and mDotP(the real mass transfer), we can find that the RHS is exactly the last two terms in Eqn.(1). (However, I still don't understand the necessity of the implementation here.)

And the difference between the mDotAlphal and mDot(which is called mDotP in interPhaseChangeFoam) is just the volume fraction of the other phase, which is correct and necessary(with current form of code). And there is a minus sign in mDotAlphal but not in mDot.（in interCondensatingEvaporatingFoam, but in interPhaseChangeFoam the signs are opposite. Because I'm not familiar with cavitation models, I cannot tell the reason...).

So, despite simple, but I hope this will help some newbees like me. : D

Please see the following link for a derivation of the mass transfer terms in the alphaEqn:
https://www.tfd.chalmers.se/~hani/ku...Yaquan_Sun.pdf

This explains how the equation is derived. However, the equation could be written with a single mass source term. I believe the reason it is formulated in this way is to linearise the equation for improved solution stability.

In this formulation, "fvm::Sp(vDotvmcAlphal, alpha1)" is a linear function of alpha1, allowing this term to be solved implicitly in the A-matrix. The "vDotcAlphal" term must still be solved explicitly, but is zero in evaporation, only taking a value when condensation occurs.

Don't necessarily take my word for it but this is the understanding I've come to.

C Dunne