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

Formulation in compressibleInterFoam

Register Blogs Community New Posts Updated Threads Search

Like Tree34Likes

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   June 6, 2013, 03:58
Default
  #41
New Member
 
NaiXian Leslie Lu
Join Date: Jun 2009
Location: France
Posts: 26
Rep Power: 17
LESlie is on a distinguished road
Quote:
Originally Posted by richpaj View Post
so that when you discretize you get something like

Ap * Up = F(Up-1, .. )

and where Ap is the diagonal coefficient.

2) as for adjustPhi(phiHbyA, U, p_rgh)

this is applied in the case of incompressible problems when the pressure
reference (or level) is unknown. If you look in the source code
(src/finiteVolume/cfdTools/general/adjustPhi/adjustPhi.C) you'll find:

if (p.needReference())
{

}
Thanks Richard, it's what I was looking for!
Have a good day.

Leslie
__________________
Cheers,
Leslie LU

LESlie is offline   Reply With Quote

Old   July 2, 2014, 07:08
Default
  #42
New Member
 
Cong Gu
Join Date: Jun 2013
Posts: 10
Rep Power: 13
gucong is on a distinguished road
Quote:
Originally Posted by richpaj View Post
dgdt =
(pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1))
*(p_rghEqnComp & p_rgh);
The code shows that
dgdt = (psi2 / rho2 - psi1 / rho1) *DDt(p - rho * g * h).

But in the derivation, it is
dgdt = (psi2 / rho2 - psi1 / rho1) * DDt(p).

They are not equivalent, are they?
gucong is offline   Reply With Quote

Old   July 2, 2014, 12:55
Default
  #43
Member
 
Richard Kenny
Join Date: Mar 2009
Posts: 64
Rep Power: 18
richpaj is on a distinguished road
Hello Gucong,
the expression looks like a variation on the original (circa OF17) but I think the same reasoning/derivation still applies, but now of course it appears the 'dynamic' and 'static' components of the pressure are clearly indicated. Recalling that the static contribution is written as p0 + rho (g . z), and hence the sign.

I've just scanned the OF23x version of compressibleInterfoam and it seems
a lot more is going on in pEqn these days, possibly with a view to stabilizing
the explicit terms. If the latter fails then and
the average density is varying considerably I think an
alternative approach similar to the old dieselFoam might be more appropriate. But that's getting off-topic.

Regards,

Richard K.
richpaj is offline   Reply With Quote

Old   July 2, 2014, 16:38
Default
  #44
New Member
 
Cong Gu
Join Date: Jun 2013
Posts: 10
Rep Power: 13
gucong is on a distinguished road
Quote:
Originally Posted by richpaj View Post
Hello Gucong,
the expression looks like a variation on the original (circa OF17) but I think the same reasoning/derivation still applies, but now of course it appears the 'dynamic' and 'static' components of the pressure are clearly indicated. Recalling that the static contribution is written as p0 + rho (g . z), and hence the sign.

I've just scanned the OF23x version of compressibleInterfoam and it seems
a lot more is going on in pEqn these days, possibly with a view to stabilizing
the explicit terms. If the latter fails then and
the average density is varying considerably I think an
alternative approach similar to the old dieselFoam might be more appropriate. But that's getting off-topic.

Regards,

Richard K.
Thank you so much for the timely reply. Regarding the definition of dgdt, it is p in OF17x. But in OF21x, it is changed to p_rgh (which means p - rho * (vector g) dot (vector x) ). According to my understanding, the change is not justified unless something else is done to compensate for the difference. But I cannot find such things in the code of OF21x. I have yet to see what the solver is like since OF22x. But here's my current question: Is the solver in OF21x correct or not?

Regards,
gucong is offline   Reply With Quote

Old   July 2, 2014, 22:21
Default
  #45
Member
 
Richard Kenny
Join Date: Mar 2009
Posts: 64
Rep Power: 18
richpaj is on a distinguished road
In p(tot) = pd + p0 + rho (g . z),
p0 = 0 and 'vector' g = (0, 0, -g), which accounts for the sign.
Rgds,
Richard K.
richpaj is offline   Reply With Quote

Old   July 2, 2014, 22:42
Default
  #46
New Member
 
Cong Gu
Join Date: Jun 2013
Posts: 10
Rep Power: 13
gucong is on a distinguished road
Quote:
Originally Posted by richpaj View Post
In p(tot) = pd + p0 + rho (g . z),
p0 = 0 and 'vector' g = (0, 0, -g), which accounts for the sign.
Rgds,
Richard K.
My problem is not with the sign. But does DDt(p) = DDt(p_rgh) in general?
gucong is offline   Reply With Quote

Old   July 2, 2014, 23:49
Default
  #47
Member
 
Richard Kenny
Join Date: Mar 2009
Posts: 64
Rep Power: 18
richpaj is on a distinguished road
Oh right, sorry for the misinterpretation.
But otherwise, DDt(p) != DDt(p_rgh) owing to the spatial cpt of DDt(..),
but presumably this latter is assumed
to be much smaller than the value of ddt(p).

It would be interesting to verify numerically.

RGK
richpaj is offline   Reply With Quote

Old   December 11, 2014, 14:30
Default
  #48
Member
 
james wilson
Join Date: Aug 2014
Location: Orlando, Fl
Posts: 39
Rep Power: 12
jameswilson620 is on a distinguished road
Kenny,

Thanks for your guidance on this topic.

Your reference to:

fvScalarMatrix psiConvectionDiffusion
(
fvm::ddt(rho, psi)
+ fv::gaussConvectionScheme<scalar>(mesh, phi, UDs).fvmDiv(phi, psi)
//- fv::gaussLaplacianScheme<scalar, scalar>(mesh, CDs, snGrads)
//.fvmLaplacian(Dpsif, psi)
- fvm::Sp(Sp, psi)
- Su
);

Is found in "IMULESTemplates.C" for a call to Foam::MULES::implicitSolve.

We however are using a call to Foam::MULES::explicitSolve found in MULESTemplates.C (as youve specified) where the structure of the fvScalarMatrix is less clear:

"
template<class RhoType, class SpType, class SuType>
void Foam::MULES::explicitSolve
(
const RhoType& rho,
volScalarField& psi,
const surfaceScalarField& phi,
surfaceScalarField& phiPsi,
const SpType& Sp,
const SuType& Su,
const scalar psiMax,
const scalar psiMin
)
{
const fvMesh& mesh = psi.mesh();
const scalar rDeltaT = 1.0/mesh.time().deltaTValue();
psi.correctBoundaryConditions();
limit(rDeltaT, rho, psi, phi, phiPsi, Sp, Su, psiMax, psiMin, 3, false);
explicitSolve(rDeltaT, rho, psi, phiPsi, Sp, Su);
}
"

Could you comment on this portion of the code and extend your explanation to cover this specific template the compressibleInterFoam solver uses? Your explanation of the source term implementation using "fvScalarMatrix psiConvectionDiffusion
(...)" makes sense but THIS template is not being used. I am assuming MULES::explicitSolve operates the same way as MULES::implicitSolve in terms of how the source terms are passed using the template where the solution of course is different.

Using this style of source term distribution in openFoam, I have been able to validate analytic Stefan type problems, so I do not doubt it's validity one bit. I would however, like to gain a deeper understanding of the form of the FV equations such that I can get creative with new source terms.

I've modified the compressibleInterFoam solver to output Su, Sp and divU (=alpha1*divU). I have a line plot along the center-line at y and a surface plot of Su and Sp at the instant the line plots are shown in the attached images. I hope this helps with visualization.

Looking at the surface plot of Sp, there are portions of the gaseous phase where Sp is non zero. I would assume that this non zero contribution in the alpha equation would lead to un-bounding of alpha in the gasseous phase. If the solver does take the form described by kenny where Sp --> Sp*alpha1 however, this non zero portion would have no contribution since alpha = 0 in this region and therefore would support Kenny's argument.

One final note, The images shown here correspond to expanding gas.. when the gas cools due to convection and diffusion of thermal energy, the sign of the source terms change since the gas is now contracting and the interface must move inward, and as a result, the assignment of divU should change as well divU --> Su (expand) vs. divU --> Sp (contract) to promote diagonal dominance. It is not treated intuitively in this manner in the default solver. Food for thought ; )

James

Su: Screenshot from 2014-12-12 08:59:22.jpg; Sp: Screenshot from 2014-12-12 09:03:51.jpg

Last edited by jameswilson620; December 12, 2014 at 10:17. Reason: adding result data
jameswilson620 is offline   Reply With Quote

Old   December 28, 2014, 05:32
Default
  #49
Senior Member
 
Dongyue Li
Join Date: Jun 2012
Location: Beijing, China
Posts: 848
Rep Power: 18
sharonyue is on a distinguished road
Quote:
Originally Posted by richpaj View Post

ddt( alpha1) + div( alpha1 * U ) = - ( alpha1 / rho1 ) * DDt( rho1 ) ..........(3)

Then, introducing a compressibility so that rho1 = rho1_0 + psi1 * p yields

ddt( alpha1) + div( alpha1 * U ) = - ( alpha1 * psi1 / rho1) DDt( p ) .......(4)
from my point, DDt( rho1 ) is not the same with (psi1 / rho1) DDt( p ). Maybe its :

DDt( rho1 ) = DDt( rho1_0 ) + (psi1 / rho1) DDt( p )



If Im wrong correct it pls.
sharonyue is offline   Reply With Quote

Old   January 5, 2015, 18:51
Default
  #50
Member
 
james wilson
Join Date: Aug 2014
Location: Orlando, Fl
Posts: 39
Rep Power: 12
jameswilson620 is on a distinguished road
So in case anyone is interested..

a call to MULES::explicitSolve(geometricOneField(), alpha1, phi, tphiAlpha(),Sp,Su, 1, 0);

references this template in MULESTemplates.C for OF2.3

Code:
template<class RhoType, class SpType, class SuType>
void Foam::MULES::explicitSolve
(
    const RhoType& rho,
    volScalarField& psi, // Note the use of the pointer, I believe this is how psi is updated without explicitly defining a "return psi"
    const surfaceScalarField& phi,
    surfaceScalarField& phiPsi,
    const SpType& Sp,
    const SuType& Su,
    const scalar psiMax,
    const scalar psiMin
)
{
    const fvMesh& mesh = psi.mesh();
    const scalar rDeltaT = 1.0/mesh.time().deltaTValue();
    psi.correctBoundaryConditions();
    limit(rDeltaT, rho, psi, phi, phiPsi, Sp, Su, psiMax, psiMin, 3, false); //I haven't looked much into this
    explicitSolve(rDeltaT, rho, psi, phiPsi, Sp, Su); //CALL TO ANOTHER TEMPLATE WHERE PSI IS SOLVED FOR THE NEXT TIME STEP
}
which references the next template (above it in MULESTemplates.C):

Code:
template<class RdeltaTType, class RhoType, class SpType, class SuType>
void Foam::MULES::explicitSolve
(
    const RdeltaTType& rDeltaT,
    const RhoType& rho,
    volScalarField& psi,
    const surfaceScalarField& phiPsi,
    const SpType& Sp,
    const SuType& Su
)
{
    Info<< "MULES: Solving for " << psi.name() << endl;

    const fvMesh& mesh = psi.mesh();

    scalarField& psiIf = psi; // initialize psiIf
    const scalarField& psi0 = psi.oldTime(); // Store old values of psi for solution

    psiIf = 0.0; // Set psiIf to zero
    fvc::surfaceIntegrate(psiIf, phiPsi); // Calculate divergence of phiPsi; update psiIf as result
//NOTE: psiIf is utilized for two purposes here. It is used to advance psi in time (what were after) and also as a temporary dummy variable to store the value of the divergence of phiPsi

    if (mesh.moving())
    {
        psiIf =
        (
            mesh.Vsc0()().field()*rho.oldTime().field()
           *psi0*rDeltaT/mesh.Vsc()().field()
          + Su.field()
          - psiIf
        )/(rho.field()*rDeltaT - Sp.field());
    }
    else
    {
// VOF EQUATION
// (psiIf(unknown) - psi0)/dt + psiIf(known dummy place holder for divergence of phiPsi) = psiIf(unknown)*Sp + Su
// OR
// ddt(psi) + div(phiPsi) = psi*Sp + Su 
// Now solve for psiIf, the unknown value of psiIf at t+dt by re arranging the above expression
        psiIf = 
        (
            rho.oldTime().field()*psi0*rDeltaT
          + Su.field()
          - psiIf  // Divergence
        )/(rho.field()*rDeltaT - Sp.field());
    }//NOTE rho.field() is geometricOneField() (simply = 1)

    psi.correctBoundaryConditions();
}
I hope this is helpful since i've struggled with this for a while.

Can anyone comment on limit(rDeltaT, rho, psi, phi, phiPsi, Sp, Su, psiMax, psiMin, 3, false) from the first template?

James
jameswilson620 is offline   Reply With Quote

Old   March 30, 2015, 03:36
Default
  #51
New Member
 
Sasa Goran
Join Date: Feb 2015
Location: Japan
Posts: 23
Rep Power: 11
Supersale is on a distinguished road
OK, i've been looking at this transportationProperties file for a couple of days now, and i'm none the wiser. I guess i figured out how to set rho0 (initial density, i guess, and equals 0 for compressible fluids) and pMin (i guess it's the minimum pressure, and should be above 0) but this psi thing bother me to death. It's compresibity, but i've never seen such a formulation, all the data i can get for psi is for the [1/MPa] unit, not the [s^2/m^2] as specified in the file.
I would grealty appreciate it if someone could help me out in finding the values for psi (ATM i need water and air) or tell me where a good database is.
Also culd you clarify my guesses about rho0 and pMin, or explain what they are indeed. Thank you guys in advance.

Edit: OK found out. The 1/MPa formulation comes from \beta = ( d_rho / d_p ) / rho

So multiplying beta by rho should give psi

Last edited by Supersale; March 31, 2015 at 02:03.
Supersale is offline   Reply With Quote

Old   May 2, 2015, 12:39
Default
  #52
New Member
 
Jianzhi Li
Join Date: Jul 2013
Location: Shanghai, China
Posts: 27
Rep Power: 13
epi_c is on a distinguished road
Send a message via Skype™ to epi_c
Quote:
Originally Posted by sharonyue View Post
from my point, DDt( rho1 ) is not the same with (psi1 / rho1) DDt( p ). Maybe its :

DDt( rho1 ) = DDt( rho1_0 ) + (psi1 / rho1) DDt( p )



If Im wrong correct it pls.
Dear Dongyue,

I saw three expressions for rho:
rho1 = rho1_0 + psi1 * p
rho1 = rho1_0 + psi1 * (p - p_0)
rho1 = psi1 * p
In all cases, rho1_0 and p_0 are constants, their material derivative should be zero.
From its definition, psi1 is isothermal compressibility, aka, be a constant when temperature fixed.
So, DDt(rho1) = psi1 * DDt(p)
epi_c is offline   Reply With Quote

Old   May 2, 2015, 12:42
Default temperature equation in compressibleInterFoam
  #53
New Member
 
Jianzhi Li
Join Date: Jul 2013
Location: Shanghai, China
Posts: 27
Rep Power: 13
epi_c is on a distinguished road
Send a message via Skype™ to epi_c
Hi guys,

Recently I'm learning the solver compressibleInterFoam in OpenFOAM-2.3.x. From the discussion above, I've already figure out the pressure equation and alpha equation, What bothers me is the temperature equation, which seems to be derived from specific total internal energy equation:

ddt(rho * e) + div(rho * U * e) + ddt(rho * K) + div(rho * U * K) = -div(U * p) + laplacian((Cp/Cv) * (kappa / Cp), e)

where e = Cv * T, e has dimension J/kg, and Cv has dimension J/kg/K.

For a multiphase flow, and in solver's description:
Quote:
The momentum and other fluid properties are of the "mixture"
usually, the mixture properties are taken as volume average of both phases, such as dynamic viscosity mu. From intuition, the mixture fluid will behave more like the phase which has more volume, so I can understand this.

But for specific total internal energy, the equation is derived based on the conservation of energy in a control volume (cell). I think the total energy in a computation cell should be taken as mass average of both phases.

Let us use Y1, Y2 denote mass fraction and alpha1, alpha2 denote the volume fraction of two phases. And Cv denotes the specific heat capacity at constant volume of the mixture, Cv1 and Cv2 for each phase. In my opinion, Cv = Y1*Cv1 + Y2*Cv2, so let the energy equation be devided by Cv, we get

ddt(rho * T) + div(rho * U * T) - laplacian((Cp/Cv) * (kappa / Cp), T)
+ (ddt(rho * K) + div(rho * U * K) + div(U * p))/Cv = 0

but in the TEqn:

Code:
    fvScalarMatrix TEqn
    (
        fvm::ddt(rho, T)
      + fvm::div(rhoPhi, T)
      - fvm::laplacian(twoPhaseProperties.alphaEff(turbulence->mut()), T)
      + (
            fvc::div(fvc::absolute(phi, U), p)
          + fvc::ddt(rho, K) + fvc::div(rhoPhi, K)
        )
       *(
           alpha1/twoPhaseProperties.thermo1().Cv()
         + alpha2/twoPhaseProperties.thermo2().Cv()
        )
    );
it can be interpreted as 1/Cv = alpha1/Cv1 + alpha2/Cv2, this expression is neither volume averaging nor mass averaging.

Also I think the term twoPhaseProperties.alphaEff(turbulence->mut()) should be twoPhaseProperties.alphaEff(turbulence->alphat()) or turbulence->alphaEff(), note that the last two expression are equivalent. The only difference is that in OpenFOAM alphat = mut / Prt. Can anybody help me with these?
epi_c is offline   Reply With Quote

Old   October 10, 2015, 18:24
Default
  #54
New Member
 
Cong Gu
Join Date: Jun 2013
Posts: 10
Rep Power: 13
gucong is on a distinguished road
Quote:
Originally Posted by gucong View Post
The code shows that
dgdt = (psi2 / rho2 - psi1 / rho1) *DDt(p - rho * g * h).

But in the derivation, it is
dgdt = (psi2 / rho2 - psi1 / rho1) * DDt(p).

They are not equivalent, are they?
In this paper, they have a term to correct this, namely the "lambda" defined equation (22)

http://www.sciencedirect.com/science...45793013001229
gucong is offline   Reply With Quote

Old   October 10, 2015, 19:36
Default
  #55
New Member
 
Cong Gu
Join Date: Jun 2013
Posts: 10
Rep Power: 13
gucong is on a distinguished road
Quote:
Originally Posted by epi_c View Post
Hi guys,

Recently I'm learning the solver compressibleInterFoam in OpenFOAM-2.3.x. From the discussion above, I've already figure out the pressure equation and alpha equation, What bothers me is the temperature equation, which seems to be derived from specific total internal energy equation:

ddt(rho * e) + div(rho * U * e) + ddt(rho * K) + div(rho * U * K) = -div(U * p) + laplacian((Cp/Cv) * (kappa / Cp), e)

where e = Cv * T, e has dimension J/kg, and Cv has dimension J/kg/K.

For a multiphase flow, and in solver's description:


usually, the mixture properties are taken as volume average of both phases, such as dynamic viscosity mu. From intuition, the mixture fluid will behave more like the phase which has more volume, so I can understand this.

But for specific total internal energy, the equation is derived based on the conservation of energy in a control volume (cell). I think the total energy in a computation cell should be taken as mass average of both phases.

Let us use Y1, Y2 denote mass fraction and alpha1, alpha2 denote the volume fraction of two phases. And Cv denotes the specific heat capacity at constant volume of the mixture, Cv1 and Cv2 for each phase. In my opinion, Cv = Y1*Cv1 + Y2*Cv2, so let the energy equation be devided by Cv, we get

ddt(rho * T) + div(rho * U * T) - laplacian((Cp/Cv) * (kappa / Cp), T)
+ (ddt(rho * K) + div(rho * U * K) + div(U * p))/Cv = 0

but in the TEqn:

Code:
    fvScalarMatrix TEqn
    (
        fvm::ddt(rho, T)
      + fvm::div(rhoPhi, T)
      - fvm::laplacian(twoPhaseProperties.alphaEff(turbulence->mut()), T)
      + (
            fvc::div(fvc::absolute(phi, U), p)
          + fvc::ddt(rho, K) + fvc::div(rhoPhi, K)
        )
       *(
           alpha1/twoPhaseProperties.thermo1().Cv()
         + alpha2/twoPhaseProperties.thermo2().Cv()
        )
    );
it can be interpreted as 1/Cv = alpha1/Cv1 + alpha2/Cv2, this expression is neither volume averaging nor mass averaging.

Also I think the term twoPhaseProperties.alphaEff(turbulence->mut()) should be twoPhaseProperties.alphaEff(turbulence->alphat()) or turbulence->alphaEff(), note that the last two expression are equivalent. The only difference is that in OpenFOAM alphat = mut / Prt. Can anybody help me with these?
I tried to explain the averaging of Cv in the code. Correct me if I'm wrong. First we have energy equation of each phase

ddt( alpha_1 * rho_1 * (Cv_1 * T + K) ) + div( alpha_1 * rho_1 * U * (Cv_1 * T + K) )
- laplacian(kappa_1 , T ) * alpha_1 * rho_1 / rho + div( p * U ) * alpha_1 * rho_1 / rho
+ alpha_1 * rho_1 * g dot U = 0

Note that the thermal diffusion and mechanical work are distributed to the phases in proportion to the mass. Then expand the first two terms by product rule and apply conservation of phase mass

ddt(alpha_1 * rho1) + div(alpha_1 * rho_1 * U) = 0

Should get

alpha_1 * rho_1 * ddt(Cv_1 * T + K) + alpha_1 * rho_1 * U dot grad(Cv_1 * T + K)

for the first two terms.

And let alpha_Eff = (alpha_1 * kappa_1 / Cv_1) + (alpha_2 * kappa_2 / Cv_2)

Divide the whole equation by (rho_1 * Cv_1), maybe assuming constant Cv_1. Then sum it with the similar equation for phase 2. Finaly multiply the sum with rho. You should arrive at similar things in the code, except some small deviations in the work of gravity. They use p_rgh instead of p, which should take care of work by gravity. But it gets complicated if rho is changing.

Last edited by gucong; October 11, 2015 at 05:33.
gucong is offline   Reply With Quote

Old   November 12, 2015, 03:07
Default
  #56
New Member
 
Jianzhi Li
Join Date: Jul 2013
Location: Shanghai, China
Posts: 27
Rep Power: 13
epi_c is on a distinguished road
Send a message via Skype™ to epi_c
Quote:
Originally Posted by gucong View Post
I tried to explain the averaging of Cv in the code. Correct me if I'm wrong. First we have energy equation of each phase

ddt( alpha_1 * rho_1 * (Cv_1 * T + K) ) + div( alpha_1 * rho_1 * U * (Cv_1 * T + K) )
- laplacian(kappa_1 , T ) * alpha_1 * rho_1 / rho + div( p * U ) * alpha_1 * rho_1 / rho
+ alpha_1 * rho_1 * g dot U = 0

Note that the thermal diffusion and mechanical work are distributed to the phases in proportion to the mass. Then expand the first two terms by product rule and apply conservation of phase mass

ddt(alpha_1 * rho1) + div(alpha_1 * rho_1 * U) = 0

Should get

alpha_1 * rho_1 * ddt(Cv_1 * T + K) + alpha_1 * rho_1 * U dot grad(Cv_1 * T + K)

for the first two terms.

And let alpha_Eff = (alpha_1 * kappa_1 / Cv_1) + (alpha_2 * kappa_2 / Cv_2)

Divide the whole equation by (rho_1 * Cv_1), maybe assuming constant Cv_1. Then sum it with the similar equation for phase 2. Finaly multiply the sum with rho. You should arrive at similar things in the code, except some small deviations in the work of gravity. They use p_rgh instead of p, which should take care of work by gravity. But it gets complicated if rho is changing.
Thanks Cong

It's more clear now. But I don't understand the div( p * U ) * alpha_1 * rho_1 / rho term, should it be div( p * U ) * alpha_1?
pmdelgado2 likes this.
epi_c is offline   Reply With Quote

Old   February 15, 2017, 07:02
Default
  #57
New Member
 
Jan Behrendt
Join Date: Nov 2016
Location: Germany
Posts: 3
Rep Power: 10
JanBeh is on a distinguished road
Hi Foamers,

thanks for the benefical discussion above. Recently i'm also trying to figure out the compressibleInterFoam Solver (Version 2.3.x) in addition, i'm pretty new to OF, so please be kind

I already figured out the alphaEqn (including the correct use of MULES) but im struggling with the pEqn. To be more exact I didn't get that pEqnComp part. I guess piprus attached file could help me a lot, but unfortunately it's not available any more .

may someone can help me understanding the correct derivation of the pEqn or can give me a hint how the main equation looks like?

maybe it's an easy, or even a silly, question but i didn't get that part and in my opinion there is a lack of compressible multiphase flow literature especially for compressibleInterFoam. Hence, i need some help

Thanks in advance
JanBeh is offline   Reply With Quote

Old   August 25, 2017, 07:03
Default
  #58
New Member
 
Jan Behrendt
Join Date: Nov 2016
Location: Germany
Posts: 3
Rep Power: 10
JanBeh is on a distinguished road
Hey gucong,
i got a question about that derivation you've posted here:

Quote:
Originally Posted by gucong View Post
I tried to explain the averaging of Cv in the code. Correct me if I'm wrong. First we have energy equation of each phase

ddt( alpha_1 * rho_1 * (Cv_1 * T + K) ) + div( alpha_1 * rho_1 * U * (Cv_1 * T + K) )
- laplacian(kappa_1 , T ) * alpha_1 * rho_1 / rho + div( p * U ) * alpha_1 * rho_1 / rho
+ alpha_1 * rho_1 * g dot U = 0

Note that the thermal diffusion and mechanical work are distributed to the phases in proportion to the mass. Then expand the first two terms by product rule and apply conservation of phase mass

ddt(alpha_1 * rho1) + div(alpha_1 * rho_1 * U) = 0

Should get

alpha_1 * rho_1 * ddt(Cv_1 * T + K) + alpha_1 * rho_1 * U dot grad(Cv_1 * T + K)

for the first two terms.

And let alpha_Eff = (alpha_1 * kappa_1 / Cv_1) + (alpha_2 * kappa_2 / Cv_2)

Divide the whole equation by (rho_1 * Cv_1), maybe assuming constant Cv_1. Then sum it with the similar equation for phase 2. Finaly multiply the sum with rho. You should arrive at similar things in the code, except some small deviations in the work of gravity. They use p_rgh instead of p, which should take care of work by gravity. But it gets complicated if rho is changing.
i agree with you that one can't interpret alpha1/cv_1+alpha2/cv_2 as 1/cv and that the cvs in the code comes from a derivation similar to yours.
i tried to follow your derivation but i guess i've made a mistakesomewhere. After dividing the whole eqn by (rho_1*Cv_1) we should get:

alpha_1*ddt(T)+alpha1*U dot grad(T) + alpha1/Cv_1*ddt(K) + alpha1/Cv_1*U dot grad(K) + div(p*U)*alpha1/Cv_1/rho - laplacian(kappa_1, T)*alpha1/Cv_1/rho

summing it with the similar eqn for phase 2 and multiply the sum with rho we should get:

rho*ddt(T) + rho*U dot(T) - laplacian(alphaEff , T) ... (+ rate of changes due to mechanical, or kinetic, energy)

the code tells us that the time derivation term should be
ddt(rho, T)

may you can tell me where i did a mistake?
dduque likes this.
JanBeh is offline   Reply With Quote

Old   October 8, 2018, 07:03
Default
  #59
New Member
 
Daniel Duque
Join Date: Oct 2018
Posts: 3
Rep Power: 8
ddcampayo is on a distinguished road
Hello, everyone.

First, thanks to all that have contributed to this thread, which has really helped me understand this solver. Specially, references to the article by Miller, Jasak, Boger, Paterson and Nedungadi (C&F 87, 2013). If I understand correctly, we have equations like:

(1) d/dt ( \rho_1 * \alpha_1) + \nabla\cdot ( \rho_1  \alpha_1 U ) =  0,

which may be written as

(2) \rho_1 [  d/dt (\alpha_1) + \nabla\cdot ( \alpha_1 U )  ] + \alpha_1  D/Dt(\rho_1) = 0,


where DDt is the convective time derivative. But if rho1 is a function of the pressure, then
DDt(rho1) = psi1 * DDt(p), where psi1 is the derivative d(rho1)/d p .

This means the equation may be written as

(3) \rho_1 [  d/dt (\alpha_1) + \nabla\cdot ( \alpha_1 U )   ]  = - \alpha_1 * \psi_1 / \rho_1 D/Dt(\rho_1)  .

Similarly for the other fraction, alpha2. Now, Miller at all start doing combinations with these equations, but on the latest OF 6 code we find this, in pEqn.h


Code:
       p_rghEqnComp1 =
            pos(alpha1)
           *(
                (
                    fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1*rho1f)
                  - (fvOptions(alpha1, mixture.thermo1().rho())&rho1)
                )/rho1
              - fvc::ddt(alpha1) - fvc::div(alphaPhi1)
              + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh))
            );
I interpret this as two chunks. First:
Code:
                (
                    fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1*rho1f)
                  - (fvOptions(alpha1, mixture.thermo1().rho())&rho1)
                )/rho1
which is basically Eq (1), divided by rho1. The term inside () should be zero. Second:
Code:
              - fvc::ddt(alpha1) - fvc::div(alphaPhi1)
              + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh))
again, Eq(3), also to be made zero.

Then, after pressure is solved for, dgdt is computed as :

Code:
            dgdt =
            (
                alpha1*(p_rghEqnComp2 & p_rgh)
              - alpha2*(p_rghEqnComp1 & p_rgh)
            );
Now, I am not really not sure about this, but I guess dgdt is then, after the iterations converge,

(4) dgdt = \alpha_1 (\alpha_2 \psi_2/\rho_2 ddt(p) )  - \alpha_2 (\alpha_1 \psi_1/\rho_1 ddt(p) )
       = \alpha_1  \alpha_2   ( \psi_2/\rho_2 - \psi_1/\rho_1)  ddt(p)

Then, we find on alphaSuSp these two terms, defined if dgdt>0 as:

Code:
        Sp[celli] -= dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
        Su[celli] += dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
Finally, on ../VoF/alphaEqn.H, we find stuff like:

Code:
        fvScalarMatrix alpha1Eqn
        (
            (
                LTS
              ? fv::localEulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
              : fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
            )
          + fv::gaussConvectionScheme<scalar>
            (
                mesh,
                phiCN,
                upwind<scalar>(mesh, phiCN)
            ).fvmDiv(phiCN, alpha1)
       // - fvm::Sp(fvc::ddt(dimensionedScalar("1", dimless, 1), mesh)
       //           + fvc::div(phiCN), alpha1)
         ==
            Su + fvm::Sp(Sp + divU, alpha1)
        );
Now, if Su and Sp are given the expressions above, with the dgdt expression, the final equation turns out to be just Eq. (3) !!! (Not completely straightforward, but it can be shown in a few lines).

So, I guess it all makes sense. I would appreciate comments on this, since in my group
we are trying to add a phase-change term to the phase equations, which would read:

ddt( \rho_1 \alpha_1) + \nabla\cdot( \rho_1  \alpha_1  U ) =   mDot
ddt( \rho_2  \alpha_2) + \nabla\cdot( \rho_2  \alpha_2  U ) =  - mDot  ,

where mDot is a source term. This will modify all these expressions, and I’d like to make sure we get it right.

Btw, I do know there is a solver that includes phase change, but that one is incompressible, and the approach is quite different.
ddcampayo is offline   Reply With Quote

Old   October 30, 2018, 07:24
Default
  #60
New Member
 
Daniel Duque
Join Date: Jan 2011
Location: ETSIN, Madrid
Posts: 28
Rep Power: 15
dduque is on a distinguished road
Quote:
Originally Posted by ddcampayo View Post
Hello, everyone.

Just a quick note to this, dduque, is my usual name. I had to use the previous one, ddcampayo, due to some glitch in the CFD-online user management system.


BTW, I am quite convinced my interpretation is mostly correct, with some technical details that I had wrong. If anyone is interested I can correct those.
Daniel_Khazaei likes this.
dduque is offline   Reply With Quote

Reply

Tags
compressible, compressibleinterfoam, theory


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
Low Reynolds Number k-epsilon formulation CFX 10.0 Chris CFX 4 December 8, 2009 00:51
Immersed Boundary Formulation Rave Main CFD Forum 0 August 11, 2008 15:55
DPM Steady formulation with collisions kulwinder FLUENT 0 May 22, 2004 19:44
energy equation formulation Pedro Phoenics 1 July 5, 2001 13:17
Compressible vs. Incompressible formulations Fernando Velasco Hurtado Main CFD Forum 3 January 7, 2000 17:51


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