CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Running, Solving & CFD

About the alphaEqn.H in interPhaseChangeFoam

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree3Likes
  • 3 Post By MrZZQi

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   March 25, 2019, 04:02
Post About the alphaEqn.H in interPhaseChangeFoam
  #1
New Member
 
Zhongqi Zuo
Join Date: May 2018
Location: China
Posts: 5
Rep Power: 5
MrZZQi is on a distinguished road
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.

Last edited by MrZZQi; March 25, 2019 at 04:04. Reason: add some code
MrZZQi is offline   Reply With Quote

Old   March 26, 2019, 09:32
Default
  #2
New Member
 
Zhongqi Zuo
Join Date: May 2018
Location: China
Posts: 5
Rep Power: 5
MrZZQi is on a distinguished road
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
cyw, wanghongjie and sherlock45 like this.

Last edited by MrZZQi; March 27, 2019 at 22:05.
MrZZQi is offline   Reply With Quote

Reply

Tags
alphaeqn.h, interphasechangefoam

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
extra term in alphaEqn.H interPhaseChangeFoam version 2.3.0 alundilong OpenFOAM Running, Solving & CFD 5 November 25, 2016 04:27
alphaEqn.H in twoPhaseEulerFoam cheng1988sjtu OpenFOAM Bugs 15 May 1, 2016 17:12
Different implementation pEqn.H in pimpleFoam vs interPhaseChangeFoam DanielRCalvete OpenFOAM Programming & Development 1 December 4, 2015 11:37
InterPhaseChangeFoam ERROR shipman OpenFOAM Running, Solving & CFD 37 March 23, 2014 13:43
Bug about MULES::implicitSolve for interPhaseChangeFoam in OF-1.6 chiven OpenFOAM Bugs 18 April 18, 2013 23:56


All times are GMT -4. The time now is 00:22.