CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM

compressibilty effect in alphaEqn.H

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

Like Tree2Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   September 27, 2011, 07:04
Default compressibilty effect in alphaEqn.H
  #1
Senior Member
 
Nima Sam
Join Date: Sep 2009
Location: Tehran, Iran
Posts: 1,111
Blog Entries: 1
Rep Power: 14
nimasam is on a distinguished road
Send a message via Yahoo to nimasam
i have two compressible flow , i would like to track interface between them
i can not use compressibleInterFoam , because the flow is self expanding due to pressure and temperature , and reaction conversion
now my question is related to VOF method for compressible flow
Code:
1-1) ddt (rho1*alpha1)  + div(rho1*alpha1*U1) = 0
1-2)ddt (rho2*alpha2)  + div(rho2*alpha2*U2) = 0

by simplifying  1-1 we have

ddt (alpha1) + div (alpha1*U1) =alpha1(ddt (rho1) + U & grad (rho1) ) /rho1

now we add and subtract div (rho1 U) in the LS of above equation, so

2) ddt (alpha1) + div (alpha1*U1) = alpha1*div(U1);

now we consider 
U = U1*alpha1* + U2*alpha2;
and
Ur = U1 - U2;
so
U1 = U + Ur*alpha2;
we rearange equation 2 to reach equation 3 :

3) ddt (alpha1) + div (alpha1*U)+div(Ur*alpha1*alpha2) = alpha1*div(U1);

the LS of above euqation is what solved with interFoam,
for interFoam which is for incompressible flows div(U1) is zero! 
so alpha1*div(U1) is the effect of compressibility of phase
now i have two ways in my mind

A) div(U1) = - (ddt(rho1) + u & grad (rho1))/rho1

or

B) div(U1) = div (U) + div(Ur*alpha2)

and none of them work wells,

A) the interface is moving but the alpha gets unbounded after several iteration

B) nothing happens, interface is fixed

any hint or idea how to consider compressibilty effect in alphaEqn ?
linch and dawnrain like this.
nimasam is offline   Reply With Quote

Old   May 16, 2012, 12:05
Default
  #2
Senior Member
 
Illya Shevchuk
Join Date: Aug 2009
Location: Darmstadt, Germany
Posts: 174
Rep Power: 7
linch is on a distinguished road
Hi Nima,

I'm working on the same problem, but I have only thermal expansivity and no compressibility. First I wold like to ask, what is your status? Did you find a proper solution?

Second question, which form of the pressure correction equation do you have?

Generally I see two way to solve the VOF and the pressure correction equation (the similarity is that both of them have to ensure mass conservation): not conservative and conservative form, like I posted here

The conservative form converges poorly and as far as I know becomes unstable for density ratio between both phases.

The non-conservative converges much better, but then there is the question "how to obtain velocity divergence". This is exact the question you asked here.

Mathematically, the RHS of your simplified eqn. 1-1 is equal to:
\frac{\alpha_1}{\rho_1} \left( \frac{\partial{\rho_1}}{\partial{t}} +  \vec{U}_1 \cdot \nabla\rho_1 \right) = \frac{\alpha_1}{\rho_1} \frac{D\rho_1}{Dt}

so I use
Code:
alpha1Su = alpha1/rho1*fvc::DDt(phi1,rho1)
for this term. It works somehow, but even if the mathematical formula is correct, such numerical implementation causes mass imbalance in the interfacial zone. May be you have an idea, how to implement this term consistently and reduce the mass imbalance?
linch is offline   Reply With Quote

Old   May 16, 2012, 13:00
Default
  #3
Senior Member
 
Nima Sam
Join Date: Sep 2009
Location: Tehran, Iran
Posts: 1,111
Blog Entries: 1
Rep Power: 14
nimasam is on a distinguished road
Send a message via Yahoo to nimasam
i calculate the source terms like below:
Code:
vDotP1 = fvc::ddt(rho1) + fvc::div(phi,rho1) - rho1*fvc:: div(phi);
vDotP2 = fvc::ddt(rho2) + fvc::div(phi,rho2)- rho2*fvc::div(phi);

vDotAlpha = (pos(alpha2)*(vDotP2/rho2) + pos(alpha1)*(vDotP1/rho1));
but now really i dont remember why i implement like this
nimasam is offline   Reply With Quote

Old   May 16, 2012, 13:41
Default
  #4
Senior Member
 
Illya Shevchuk
Join Date: Aug 2009
Location: Darmstadt, Germany
Posts: 174
Rep Power: 7
linch is on a distinguished road
And you don't experience any problems with mass balance? I'll try your implementation, thanks a lot.
linch is offline   Reply With Quote

Old   May 16, 2012, 13:49
Default
  #5
Senior Member
 
Illya Shevchuk
Join Date: Aug 2009
Location: Darmstadt, Germany
Posts: 174
Rep Power: 7
linch is on a distinguished road
in your implementation, vDotP1 is equal to fvc:Dt(phi,rho1), and vDotP2 to DDt(phi,rho2) respectively.

The third line looks like velocity divergence for entire flow, and not for a single phase, and this confuses me a little. Is your vDotAlpha an explicit source term for alpha1? Could you please tell me, how your MULES looks like?
linch is offline   Reply With Quote

Old   May 16, 2012, 13:58
Default
  #6
Senior Member
 
Nima Sam
Join Date: Sep 2009
Location: Tehran, Iran
Posts: 1,111
Blog Entries: 1
Rep Power: 14
nimasam is on a distinguished road
Send a message via Yahoo to nimasam
i think it was some how poor in continuity specially in cumulative continuity was some how high but golbal and local continuity error was around 1e-05
Attached Files
File Type: h alphaEqns.H (2.3 KB, 58 views)
nimasam is offline   Reply With Quote

Old   May 16, 2012, 14:17
Default
  #7
Senior Member
 
Illya Shevchuk
Join Date: Aug 2009
Location: Darmstadt, Germany
Posts: 174
Rep Power: 7
linch is on a distinguished road
Thanks again.

In you code I find this even more strage:
Code:
Sp[celli] -= vDotAlpha[celli]*alpha1[celli];
Su[celli] += vDotAlpha[celli]*alpha1[celli];
Though they have the same dimensions, their meaning is different. Su = \alpha_1 (\nabla \cdot U_1) while Sp = (\nabla \cdot U_1)
Due to my logic
Code:
Sp[celli] -= vDotAlpha[celli];
Su[celli] += vDotAlpha[celli]*alpha1[celli];
would be the "right" pair. But implementation in compressibleInterFoam looks similar to yours. Hm, I still don't get it.

Last edited by linch; May 18, 2012 at 04:35.
linch is offline   Reply With Quote

Old   May 16, 2012, 15:06
Default
  #8
Senior Member
 
Nima Sam
Join Date: Sep 2009
Location: Tehran, Iran
Posts: 1,111
Blog Entries: 1
Rep Power: 14
nimasam is on a distinguished road
Send a message via Yahoo to nimasam
look at here
nimasam is offline   Reply With Quote

Old   May 21, 2012, 08:38
Default
  #9
Senior Member
 
Illya Shevchuk
Join Date: Aug 2009
Location: Darmstadt, Germany
Posts: 174
Rep Power: 7
linch is on a distinguished road
Hi Nima,

me again here Sry, we had a long weekend here in Germany. Back to the topic:

- I read again the thread you linked, but it is still not 100% clear for me, how can Sp and Su have the same form?

- In the compressibleInterFoam our EOS reads rho1 = rho10 + psi1*p, so we can recast DDt(rho1) into psi1*DDt(p). But since you have temperature & pressure dependency of the density, means your EOS reads rho=rho(p,T), you have to modify the pEqn. May I ask how did you do it? Where do you update the density and what is the RHS of your pressure correction eqn.?

- is your vDotAlpha the same as divU ?

Best,
Illya

Last edited by linch; May 21, 2012 at 12:26.
linch is offline   Reply With Quote

Old   May 22, 2012, 01:28
Default
  #10
Senior Member
 
Nima Sam
Join Date: Sep 2009
Location: Tehran, Iran
Posts: 1,111
Blog Entries: 1
Rep Power: 14
nimasam is on a distinguished road
Send a message via Yahoo to nimasam
Hi Illya
psi is somehow (1/RT), so after temperature calculation, i update rho,
then i calculate source term in pEqn from rho not p!, so it became an explicit source term for pressure equation
nimasam is offline   Reply With Quote

Old   May 22, 2012, 04:13
Default
  #11
Senior Member
 
Illya Shevchuk
Join Date: Aug 2009
Location: Darmstadt, Germany
Posts: 174
Rep Power: 7
linch is on a distinguished road
Hi Nima.

Explicit source term in pEqn means div(U), right? Do you have something like
Code:
// update density
// update vDotP1 & vDotP2

divU = alpha1*vDotP1/rho1 + alpha2*vDotP2/rho2;

fvScalarMatrix p_rghEqnIncomp
(
    divU
 + fvc::div(phi)
  - fvm::laplacian(rUAf, p_rgh)
);

Last edited by linch; May 22, 2012 at 04:36.
linch is offline   Reply With Quote

Old   May 22, 2012, 08:40
Default
  #12
Senior Member
 
Illya Shevchuk
Join Date: Aug 2009
Location: Darmstadt, Germany
Posts: 174
Rep Power: 7
linch is on a distinguished road
I've implemented the pressure correction like above, but the simple loop doesn't converge. Could you take a look and say, what do you think about it?

In expandableInterFoam.C:
Code:
        while (pimple.loop())
        {
            #include "updateSinglePhaseProperties.H" // rho1, rho2 updated here

            Drho1Dt   = fvc::ddt(rho1) + fvc::div(phi,rho1) - fvc::div(phi)*rho1;
            Drho2Dt   = fvc::ddt(rho2) + fvc::div(phi,rho2) - fvc::div(phi)*rho2;
            rho     == alpha1*rho1 + (scalar(1) - alpha1)*rho2;    // update rho-field (changes due to temperature transport)

            divU = (alpha1*Drho1Dt/rho1 + (scalar(1)-alpha1)*Drho2Dt/rho2);
            
            #include "UEqn.H"
            
            // --- Pressure corrector loop
            while (pimple.correct())
            {
                #include "pEqn.H"
            }
            

            #include "alphaEqnSubCycle.H" // rho, alpha1, phiAlpha updated here
            
            #include "updateFluxes.H" // rhoPhi1, rhoPhi, rhoCpPhi updated here
            
            #include "TEqn.H"
        }
        {
        /*************** continuity check ***************************/
            volScalarField rhoAlpha1 = rho1*alpha1;
            rhoAlpha1.oldTime() = rho1.oldTime()*alpha1.oldTime();
            
            volScalarField rhoAlpha2 = rho2*(scalar(1)-alpha1);
            rhoAlpha2.oldTime() = rho2.oldTime()*(scalar(1)-alpha1.oldTime());

            contErr   = fvc::ddt(rho)+fvc::div(rhoPhi);
            contErr1 = fvc::ddt(rhoAlpha1)+fvc::div(rhoPhi1);
        /***********************************************************/
        }
in pEqn.H:
Code:
    while (pimple.correctNonOrthogonal())
    {
        fvScalarMatrix p_rghEqn
        (
            fvm::laplacian(rAUf, p_rgh) 
            == 
            fvc::div(phi) + divU
        );

        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));

        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));

        if (pimple.finalNonOrthogonalIter())
        {
            phi -= p_rghEqn.flux();
        }
    }
and in alphaEqn.H:
Code:
    for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
    {
        volScalarField dgdt
        (
            IOobject
            (
                "dgdt",
                runTime.timeName(),
                mesh
            ),
            (pos(scalar(1.0) - alpha1)*(Drho2Dt/rho2) - pos(alpha1)*(Drho1Dt/rho1))
        );
        
        divU = fvc::div(phi);

        alpha1Sp = dimensionedScalar("Sp",dimensionSet(0,0,-1,0,0,0,0),0);
        alpha1Su = divU*alpha1;
        

        forAll(dgdt, celli)
        {
            if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0)
            {
                alpha1Sp[celli] -= dgdt[celli]*alpha1[celli];
                alpha1Su[celli] += dgdt[celli]*alpha1[celli];
            }
            else if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0)
            {
                alpha1Sp[celli] += dgdt[celli]*(1.0 - alpha1[celli]);
            }
        }
        
        
        phiAlpha =
            fvc::flux
            (
                phi,
                alpha1,
                alphaScheme
            )
          + fvc::flux
            (
                -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
                alpha1,
                alpharScheme
            );
        MULES::explicitSolve(geometricOneField(), alpha1, phi, phiAlpha, alpha1Sp, alpha1Su, 1, 0);
    }
For Drho1Dt & Drho2Dt I use Gauss linear div-schemes:
Code:
div(phi,rho1)               Gauss linear;
div(phi,rho2)               Gauss linear;
What are significant differences from your implementation? alphaEqn looks similar, Drho1Dt & Drho2Dt also (you called them vDotP1 & vDotP2 respectively). Equation solution order should not matter, if the solution is converged. The only possible difference I see is the pressure correction eqn. and the source term in it.
linch is offline   Reply With Quote

Old   May 22, 2012, 09:03
Default
  #13
Senior Member
 
Illya Shevchuk
Join Date: Aug 2009
Location: Darmstadt, Germany
Posts: 174
Rep Power: 7
linch is on a distinguished road
update: I had a small typo, it converges. But the mass conservation is still poor. I get local error in the order of 1e+02 and not 1e-05
linch is offline   Reply With Quote

Old   May 24, 2012, 07:40
Default
  #14
Senior Member
 
Illya Shevchuk
Join Date: Aug 2009
Location: Darmstadt, Germany
Posts: 174
Rep Power: 7
linch is on a distinguished road
Hi Nima,

is my procedure the same as yours, or do you have different pressure correction eqn.?

Best regards,
Illya
linch is offline   Reply With Quote

Old   May 24, 2012, 10:16
Default
  #15
Senior Member
 
Nima Sam
Join Date: Sep 2009
Location: Tehran, Iran
Posts: 1,111
Blog Entries: 1
Rep Power: 14
nimasam is on a distinguished road
Send a message via Yahoo to nimasam
hi dear Illya
sorry4 late answer, take a quick look on your code, it seems look a like, ofcourse i implemented it for PISO loop in OpenFOAM-1.6, however i suggest you play with algorithm
i prefer for example in each loop first update alpha1, then i have rhoPhi, then update rho from rhoEqn (rhoPhi), and so on ....
keep the structure of compressibleInterFoam
nimasam is offline   Reply With Quote

Old   May 24, 2012, 11:18
Default
  #16
Senior Member
 
Illya Shevchuk
Join Date: Aug 2009
Location: Darmstadt, Germany
Posts: 174
Rep Power: 7
linch is on a distinguished road
Hi Nima,

thanks for the reply. On my opinion the order of equations doesn't really matter for me, since I have an outer (SIMPLE) loop over all equations till the solution converges. The reason to solve the alphaEqn after the pressure correction is following: alphaEqn needs volumetric fluxes "phi" coming from the momentum equation, so the fluxes of the alpha1-phase "phiAlpha" are consistent with the volumetric fluxes "phi" in this way.

What really wracks my brain is the mass imbalance. If I understood you write, you additionally solve the rhoEqn to ensure the continuity is satisfied I think this could be the solution of my problem.

I'm sorry to terrorize you with my endless questions :-), but I have a couple of them again:
1) Do you solve the rhoEqn twice, for each phase?
2) Do you solve the rhoEqn after the PISO loop, means at the end of your time step?

A lot of thanks and best regards,
Illya
linch is offline   Reply With Quote

Old   May 24, 2012, 14:45
Default
  #17
Senior Member
 
Nima Sam
Join Date: Sep 2009
Location: Tehran, Iran
Posts: 1,111
Blog Entries: 1
Rep Power: 14
nimasam is on a distinguished road
Send a message via Yahoo to nimasam
i solved rhoEqn once,
when I solve alphaEqn, rhoPhi is getting update,
then i calculate total rho from rhoPhi from continuity Eqn (look compressibleInterFoam)
nimasam is offline   Reply With Quote

Old   May 25, 2012, 04:52
Default
  #18
Senior Member
 
Illya Shevchuk
Join Date: Aug 2009
Location: Darmstadt, Germany
Posts: 174
Rep Power: 7
linch is on a distinguished road
Thank you Nima,

apparently I'm to dumb to get it

1) solving the continuity eqn. for rho could ensure the overall continuity, but not the continuity of each phase

2) solving rhoEqn in between doesn't ensure the mass balance, because of the density update following in the end of the timestep

3) we also implicitly include conti in alphaEqn and pEqn, but in a non-conservative form, what means only if the source terms in these equations are derived consistently with other equations and their discretization, mass conservation can be achieved

4) I use the same source term as you. Your mass balance is good, my not. Why?

Do you use ideal gas law for the gas phase? Which EOS do you use for the liquid? How high are your temperature caused density variations beneath each phase (in %)? There must be a reasonable reason

Best regards,
Illya
linch is offline   Reply With Quote

Old   May 25, 2012, 05:18
Default
  #19
Senior Member
 
Nima Sam
Join Date: Sep 2009
Location: Tehran, Iran
Posts: 1,111
Blog Entries: 1
Rep Power: 14
nimasam is on a distinguished road
Send a message via Yahoo to nimasam
Quote:
Originally Posted by linch View Post
1) solving the continuity eqn. for rho could ensure the overall continuity, but not the continuity of each phase
you are right!

Quote:
2) solving rhoEqn in between doesn't ensure the mass balance, because of the density update following in the end of the timestep
why?
in each iteration, update rho total, based on continuity
in end of iteration update rho based on mixture equation ( i guess this help each phase mass conservation)

Quote:
3) we also implicitly include conti in alphaEqn and pEqn, but in a non-conservative form, what means only if the source terms in these equations are derived consistently with other equations and their discretization, mass conservation can be achieved
????

Quote:
4) I use the same source term as you. Your mass balance is good, my not. Why?
i really dont know

Quote:
Do you use ideal gas law for the gas phase? Which EOS do you use for the liquid? How high are your temperature caused density variations beneath each phase (in %)? There must be a reasonable reason
actually i had two phase one phase was air (perfect gas) and the other one was polyurethane foam which density changes based on temperature, chemical conversion and ...)
nimasam is offline   Reply With Quote

Old   May 25, 2012, 05:28
Default
  #20
Senior Member
 
Illya Shevchuk
Join Date: Aug 2009
Location: Darmstadt, Germany
Posts: 174
Rep Power: 7
linch is on a distinguished road
3) I mean, the source terms in alphaEqn and pEqn are based on single phase continuity equations.

If you have already published your results, would you mind giving me the whole solver to examine it carefully?
linch is offline   Reply With Quote

Reply

Tags
alpha1, compressibility, continuty, interfoam, mass balance, vof

Thread Tools
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 On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Can Flow-3D solve the Joule-Thomson effect? seasoul FLOW-3D 8 December 21, 2011 01:01
alphaEqn.H in twoPhaseEulerFoam cheng1988sjtu OpenFOAM Bugs 14 September 18, 2011 13:46
Can 'shock waves' occur in viscous fluid flows? diaw Main CFD Forum 105 November 20, 2009 05:19
VOF and strong surface tension effect phsieh2005 Main CFD Forum 1 August 30, 2005 06:56
Iclude particle diameter effect in wall? yueroo FLUENT 0 April 14, 2001 03:56


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