Formulation in compressibleInterFoam

Daniel_Khazaei · May 20, 2019, 06:36

Quote:

Originally Posted by dduque

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.

Hi,

Well I'm interested

Would you please share your findings?

I'm trying to do the exact same thing.

Regards

saeedbd · April 17, 2020, 08:49

Hi foamers,

I am working on CompressibleInterFoam. I want to convert the T equation to a general energy conservation equation to remove some assumptions in the solver. Could you finally manage the energy equation?

Assuming the perfectGas and perfectFluid Eos, we have e=Cv*T for each phase where e is internal energy. So, I am replacing T with et/(Cv_water Y_water + Cv_air Y_air) in the solver (et is the total internal energy with et = e1*Y1 + e2*Y2. But it gives me wrong results for bubble collapse case. So in fact, the T equation in original solver is as:

Code:

fvScalarMatrix TEqn                                                                 
    (                                                                                   
        fvm::ddt(rho, T) + fvm::div(rhoPhi, T) - fvm::Sp(contErr, T)                    
      - fvm::laplacian(turbulence.alphaEff(), T)                                        
      + (                                                                               
            fvc::div(fvc::absolute(phi, U), p)()() // - contErr/rho*p                   
          + (fvc::ddt(rho, K) + fvc::div(rhoPhi, K))()() - contErr*K                    
        )                                                                               
       *(                                                                               
           alpha1()/mixture.thermo1().Cv()()                                            
         + alpha2()/mixture.thermo2().Cv()()                                            
        )                                                                               
     ==                                                                                 
        fvOptions(rho, T)                                                               
    );

and now is converted to

Code:

fvm::ddt(rho, et) + fvm::div(rhoPhi, et) - fvm::Sp(contErr, et)
      - fvm::laplacian(turbulence.alphaEff(), et)
      + (
            fvc::div(fvc::absolute(phi, U), p)()()
           // - contErr/rho*p 
          + (fvc::ddt(rho, K) + fvc::div(rhoPhi, K))()() - contErr*K
        )
      *(
           alpha1()/mixture.thermo1().Cv()() + alpha2()/mixture.thermo2().Cv()()
       )
      *(
         Y1*mixture.thermo1().Cv()()+Y2*mixture.thermo2().Cv()()                        
         )
        ==
        fvOptions(rho, et)

where Y is the mass fraction: et=Y1*e1 + Y2*e2 and e1=Cv1*T , e2=Cv2*T.

Also, I tried:

Code:

fvm::ddt(rho, et) + fvm::div(rhoPhi, et) - fvm::Sp(contErr, et)
      - fvm::laplacian(turbulence.alphaEff(), et)
      + (
            fvc::div(fvc::absolute(phi, U), p)()()
           // - contErr/rho*p 
          + (fvc::ddt(rho, K) + fvc::div(rhoPhi, K))()() - contErr*K
        )
        ==
        fvOptions(rho, et)

But it does not work properly as well. The temperature drops at the interface immediately. In more details, I don't know why et in the new energy equation remains constant while I checked that we have velocity and K is NOT zero! (How can it be possible?) Then, after a very small time step 10^-12, since we have Y_water = 0.1 in the gas side of the interface, Temperature drops suddenly according to T = et/(Cv_water*Y_water + Cv_air*Y_air).

The last point, it gives good results for the shock tube case which is not very dependent on temperature and energy.

Do you know how to solve the issue? Thank you for your time.

Saeed,

Tobi · April 17, 2020, 13:19

Hi,

so I am not familiar with the solver you are using and I am also not familiar with its equation. However, VOF solves only one question for energy/momentum and so on. Just one question, you are solving for the field et, which is what? An energy right, the total internal energy. So, how do you transform the energy back into the temperature? I cannot see anything here. Because, as you are not solving for T anymore, the T field should stay constant.

So actually, what you have to do is:

If possible, introduce a energy equation (maybe the one you have is fine)
Furthermore, you need to have the transformation from T -> e and back from e -> T

Hope this might help.
Tobi

saeedbd · April 17, 2020, 13:44

Hi Tobi,

Thank you very much for giving your feedback. Regarding the Temperature, I am obtaining it within the T = et/(Cv_water*Y_water + Cv_air*Y_air) equation which is based on perfectFluid and perfectGas. So to put it into a nutshell, I initialize the system with et=e1Y1+e2*Y2 and then solve the et equation. Then I update the T, based on the new et using the above equation.

Best,

Tobi · April 17, 2020, 14:31

Okay, got it. Is e = Cv T correct?

saeedbd · April 17, 2020, 14:32

Yes with the considered EoSs.
I found an issue with the defined tolerance for et. So now it is being solved and this problem "I don't know why et in the new energy equation remains constant" is solved.

MrZZQi · April 29, 2020, 04:26

Hi openFoamers,

Thanks for your great discussion, it really helps a lot.

However, when I look through the code, I find this term very confusing to me

Code:

p_rghEqnComp1 & p_rgh

This term is explained as 'p_rgh in p_rghEqnComp1' by Richard, and treated as 'Dp/Dt', so did Daniel do in his derivation.

But as I know & is "inner product", the operator is redefined in fvMatrix.C. p_rghEqnComp1 is a fvMatrix, p_rgh is a scalar field, how can the term be interpreted as Dp/Dt?

Anybody can help? Thanks a lot.

Zuo

Shibi · September 9, 2021, 15:36

Hello to all,

I am also trying to understand the formulation of compressibleInterFoam for OpenFOAM ESI. However, the formulation is not clear to me.

I would like to request a further clarification on the formulation of compressibleInterFoam.

From the paper of Miller et al [1], previously referred to in this post, the pressure equation should be:

$\left ( \frac{\alpha_1 \psi_1}{\rho1}+\frac{\alpha_2 \psi_2}{\rho2} \right )\left [ \frac{\partial P}{\partial t} + \vec{v} \cdot \nabla P\right ] + \nabla \cdot \vec{v}=0$

Or, to better use the divergence theorem:

$\left ( \frac{\alpha_1 \psi_1}{\rho1}+\frac{\alpha_2 \psi_2}{\rho2} \right )\left [ \frac{\partial p}{\partial t} + \nabla \cdot \left ( p \vec{v} \right) - p \left ( \nabla \cdot \vec{v} \right) \right] + \nabla \cdot \vec{v}=0$

With the last term on the lefh hand side, $\nabla \cdot \vec{v}$ , we create a incompressiblePEqn and with the remainder terms a compressiblePEqn.

I took a look into the compressibleInterFoam from foam-exted. The pEqn.H has:

An incomplete CompressiblePEqn written as:

Code:

    if (pimple.transonic())
    {
        pEqnComp =
            (fvm::ddt(p) + fvm::div(phi, p) - fvm::Sp(fvc::div(phi), p));
    }
    else
    {
        pEqnComp =
            (fvm::ddt(p) + fvc::div(phi, p) - fvc::Sp(fvc::div(phi), p));
    }

Which should be: $\left [ \frac{\partial p}{\partial t} + \nabla \cdot \left ( p \vec{v} \right) - p \left ( \nabla \cdot \vec{v} \right) \right]$

The incompressiblePEqn which is derived by substituting the cell center velocity from the momentum equation.

Code:

        fvScalarMatrix pEqnIncomp
        (
            fvc::div(phi)
          - fvm::laplacian(rUAf, p)
        );

Now, the complete pressure Equation is assembled and solved as:

Code:

        solve
        (
            (
                max(alpha1, scalar(0))*(psi1/rho1)
              + max(alpha2, scalar(0))*(psi2/rho2)
            )
           *pEqnComp()
          + pEqnIncomp
        );

$\left ( \frac{\alpha_1 \psi_1}{\rho1}+\frac{\alpha_2 \psi_2}{\rho2} \right )\left [ \frac{\partial p}{\partial t} + \nabla \cdot \left ( p \vec{v} \right) - p \left ( \nabla \cdot \vec{v} \right) \right] + \nabla \cdot \vec{v}=0$

After solving the pressure equation, a dgdt term is updated with:

Code:

            dgdt =
                (pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1))
                *(pEqnComp & p);

The volume fraction equation in [1] is written as:

$\frac{\partial \alpha_1}{\partial t} + \nabla \cdot \left ( \alpha_1 \vec{v} \right) - \alpha_1 \left ( \nabla \cdot \vec{v} \right) = \alpha_1 \alpha_2\left ( \frac{\psi_1\rho_2 - \psi_2 \rho_1}{\alpha_1 \psi_1 \rho_2 + \alpha_2 \psi_2 \rho_1} \right ) \nabla \cdot \vec{v}$

Here are my questions:

1º Is the formulation used in [1], the one in compressibleInterFoam?

2º Is it the same formulation on all branches of OpenFOAM (ESI, ORG and extend)?

3º is dtdg = right and side of the above equation?

4º What is the purpose of doing the dot product (&) between a scalarMatrix and a volScalarField? What is the mathematical representation of this operation?

5º In OpenFOAM(ESI) the compressible part (non-transonic) of the pressure equation we have:

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))
             );

Could anyone explain the lines:

Code:

                    fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1*rho1f)
                  - (fvOptions(alpha1, mixture.thermo1().rho())&rho1)
                )/rho1
               - fvc::ddt(alpha1) - fvc::div(alphaPhi1)

Additionally why is it not using the terms: $\nabla \cdot \left ( p \vec{v} \right) - p \left ( \nabla \cdot \vec{v} \right)$ ?

6º Apart from [1] is there any other literature on the formulation of compressibleInterFoam?

Best Regards!

----
[1] -Miller, S. T., Jasak, H., Boger, D. A., Paterson, E. G., & Nedungadi, A. (2013). A pressure-based, compressible, two-phase flowa finite volume method for underwater explosions. Computers and Fluids, 87, 132–143. https://doi.org/10.1016/j.compfluid.2013.04.002

saicharan662000@gmail.com · December 30, 2021, 01:04

Hi Nat,
I know the post is old. But can you describe the evolution of step 4 mathematically instead of describing it in words?
Thanks in advance.

saicharan662000@gmail.com · December 30, 2021, 01:47

Hello nat,
Can u explain the evolution of step 3 to step 4 using math instead of words.
Thanks in advance.

April 17, 2020, 08:49	Energy equation	#62
saeedbd New Member saeedbd Join Date: Feb 2015 Posts: 6 Rep Power: 9	Hi foamers, I am working on CompressibleInterFoam. I want to convert the T equation to a general energy conservation equation to remove some assumptions in the solver. Could you finally manage the energy equation? Assuming the perfectGas and perfectFluid Eos, we have e=CvT for each phase where e is internal energy. So, I am replacing T with et/(Cv_water Y_water + Cv_air Y_air) in the solver (et is the total internal energy with et = e1Y1 + e2Y2. But it gives me wrong results for bubble collapse case. So in fact, the T equation in original solver is as: Code: fvScalarMatrix TEqn ( fvm::ddt(rho, T) + fvm::div(rhoPhi, T) - fvm::Sp(contErr, T) - fvm::laplacian(turbulence.alphaEff(), T) + ( fvc::div(fvc::absolute(phi, U), p)()() // - contErr/rhop + (fvc::ddt(rho, K) + fvc::div(rhoPhi, K))()() - contErrK ) ( alpha1()/mixture.thermo1().Cv()() + alpha2()/mixture.thermo2().Cv()() ) == fvOptions(rho, T) ); and now is converted to Code: fvm::ddt(rho, et) + fvm::div(rhoPhi, et) - fvm::Sp(contErr, et) - fvm::laplacian(turbulence.alphaEff(), et) + ( fvc::div(fvc::absolute(phi, U), p)()() // - contErr/rhop + (fvc::ddt(rho, K) + fvc::div(rhoPhi, K))()() - contErrK ) ( alpha1()/mixture.thermo1().Cv()() + alpha2()/mixture.thermo2().Cv()() ) ( Y1mixture.thermo1().Cv()()+Y2mixture.thermo2().Cv()() ) == fvOptions(rho, et) where Y is the mass fraction: et=Y1e1 + Y2e2 and e1=Cv1T , e2=Cv2T. Also, I tried: Code: fvm::ddt(rho, et) + fvm::div(rhoPhi, et) - fvm::Sp(contErr, et) - fvm::laplacian(turbulence.alphaEff(), et) + ( fvc::div(fvc::absolute(phi, U), p)()() // - contErr/rhop + (fvc::ddt(rho, K) + fvc::div(rhoPhi, K))()() - contErrK ) == fvOptions(rho, et) But it does not work properly as well. The temperature drops at the interface immediately. In more details, I don't know why et in the new energy equation remains constant while I checked that we have velocity and K is NOT zero! (How can it be possible?) Then, after a very small time step 10^-12, since we have Y_water = 0.1 in the gas side of the interface, Temperature drops suddenly according to T = et/(Cv_waterY_water + Cv_airY_air). The last point, it gives good results for the shock tube case which is not very dependent on temperature and energy. Do you know how to solve the issue? Thank you for your time. Saeed,

April 17, 2020, 13:19		#63
Tobi Super Moderator Tobias Holzmann Join Date: Oct 2010 Location: Tussenhausen Posts: 2,654 Blog Entries: 6 Rep Power: 48	Hi, so I am not familiar with the solver you are using and I am also not familiar with its equation. However, VOF solves only one question for energy/momentum and so on. Just one question, you are solving for the field et, which is what? An energy right, the total internal energy. So, how do you transform the energy back into the temperature? I cannot see anything here. Because, as you are not solving for T anymore, the T field should stay constant. So actually, what you have to do is: If possible, introduce a energy equation (maybe the one you have is fine) Furthermore, you need to have the transformation from T -> e and back from e -> T Hope this might help. Tobi saeedbd likes this. __________________ Keep foaming, Tobias Holzmann

April 17, 2020, 14:31		#65
Tobi Super Moderator Tobias Holzmann Join Date: Oct 2010 Location: Tussenhausen Posts: 2,654 Blog Entries: 6 Rep Power: 48	Okay, got it. Is e = Cv T correct? __________________ Keep foaming, Tobias Holzmann

April 17, 2020, 14:32		#66
saeedbd New Member saeedbd Join Date: Feb 2015 Posts: 6 Rep Power: 9	Yes with the considered EoSs. I found an issue with the defined tolerance for et. So now it is being solved and this problem "I don't know why et in the new energy equation remains constant" is solved. Tobi likes this.

April 29, 2020, 04:26		#67
MrZZQi New Member Zhongqi Zuo Join Date: May 2018 Location: China Posts: 5 Rep Power: 5	Hi openFoamers, Thanks for your great discussion, it really helps a lot. However, when I look through the code, I find this term very confusing to me Code: p_rghEqnComp1 & p_rgh This term is explained as 'p_rgh in p_rghEqnComp1' by Richard, and treated as 'Dp/Dt', so did Daniel do in his derivation. But as I know & is "inner product", the operator is redefined in fvMatrix.C. p_rghEqnComp1 is a fvMatrix, p_rgh is a scalar field, how can the term be interpreted as Dp/Dt? Anybody can help? Thanks a lot. Zuo

April 17, 2020, 13:44		#64
saeedbd New Member saeedbd Join Date: Feb 2015 Posts: 6 Rep Power: 9	Hi Tobi, Thank you very much for giving your feedback. Regarding the Temperature, I am obtaining it within the T = et/(Cv_waterY_water + Cv_airY_air) equation which is based on perfectFluid and perfectGas. So to put it into a nutshell, I initialize the system with et=e1Y1+e2*Y2 and then solve the et equation. Then I update the T, based on the new et using the above equation. Best,

September 9, 2021, 15:36		#68
Shibi Member Zé Join Date: Feb 2020 Posts: 90 Rep Power: 4	Hello to all, I am also trying to understand the formulation of compressibleInterFoam for OpenFOAM ESI. However, the formulation is not clear to me. I would like to request a further clarification on the formulation of compressibleInterFoam. From the paper of Miller et al [1], previously referred to in this post, the pressure equation should be: $\left ( \frac{\alpha_1 \psi_1}{\rho1}+\frac{\alpha_2 \psi_2}{\rho2} \right )\left [ \frac{\partial P}{\partial t} + \vec{v} \cdot \nabla P\right ] + \nabla \cdot \vec{v}=0$ Or, to better use the divergence theorem: $\left ( \frac{\alpha_1 \psi_1}{\rho1}+\frac{\alpha_2 \psi_2}{\rho2} \right )\left [ \frac{\partial p}{\partial t} + \nabla \cdot \left ( p \vec{v} \right) - p \left ( \nabla \cdot \vec{v} \right) \right] + \nabla \cdot \vec{v}=0$ With the last term on the lefh hand side, $\nabla \cdot \vec{v}$ , we create a incompressiblePEqn and with the remainder terms a compressiblePEqn. I took a look into the compressibleInterFoam from foam-exted. The pEqn.H has: An incomplete CompressiblePEqn written as: Code: if (pimple.transonic()) { pEqnComp = (fvm::ddt(p) + fvm::div(phi, p) - fvm::Sp(fvc::div(phi), p)); } else { pEqnComp = (fvm::ddt(p) + fvc::div(phi, p) - fvc::Sp(fvc::div(phi), p)); } Which should be: $\left [ \frac{\partial p}{\partial t} + \nabla \cdot \left ( p \vec{v} \right) - p \left ( \nabla \cdot \vec{v} \right) \right]$ The incompressiblePEqn which is derived by substituting the cell center velocity from the momentum equation. Code: fvScalarMatrix pEqnIncomp ( fvc::div(phi) - fvm::laplacian(rUAf, p) ); Now, the complete pressure Equation is assembled and solved as: Code: solve ( ( max(alpha1, scalar(0))(psi1/rho1) + max(alpha2, scalar(0))(psi2/rho2) ) pEqnComp() + pEqnIncomp ); $\left ( \frac{\alpha_1 \psi_1}{\rho1}+\frac{\alpha_2 \psi_2}{\rho2} \right )\left [ \frac{\partial p}{\partial t} + \nabla \cdot \left ( p \vec{v} \right) - p \left ( \nabla \cdot \vec{v} \right) \right] + \nabla \cdot \vec{v}=0$ After solving the pressure equation, a dgdt term is updated with: Code: dgdt = (pos(alpha2)(psi2/rho2) - pos(alpha1)(psi1/rho1)) (pEqnComp & p); The volume fraction equation in [1] is written as: $\frac{\partial \alpha_1}{\partial t} + \nabla \cdot \left ( \alpha_1 \vec{v} \right) - \alpha_1 \left ( \nabla \cdot \vec{v} \right) = \alpha_1 \alpha_2\left ( \frac{\psi_1\rho_2 - \psi_2 \rho_1}{\alpha_1 \psi_1 \rho_2 + \alpha_2 \psi_2 \rho_1} \right ) \nabla \cdot \vec{v}$ Here are my questions: 1º Is the formulation used in [1], the one in compressibleInterFoam? 2º Is it the same formulation on all branches of OpenFOAM (ESI, ORG and extend)? 3º is dtdg = right and side of the above equation? 4º What is the purpose of doing the dot product (&) between a scalarMatrix and a volScalarField? What is the mathematical representation of this operation? 5º In OpenFOAM(ESI) the compressible part (non-transonic) of the pressure equation we have: Code: p_rghEqnComp1 = pos(alpha1) ( ( fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1rho1f) - (fvOptions(alpha1, mixture.thermo1().rho())&rho1) )/rho1 - fvc::ddt(alpha1) - fvc::div(alphaPhi1) + (alpha1psi1/rho1)correction(fvm::ddt(p_rgh)) ); Could anyone explain the lines: Code: fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1rho1f) - (fvOptions(alpha1, mixture.thermo1().rho())&rho1) )/rho1 - fvc::ddt(alpha1) - fvc::div(alphaPhi1) Additionally why is it not using the terms: $\nabla \cdot \left ( p \vec{v} \right) - p \left ( \nabla \cdot \vec{v} \right)$ ? 6º Apart from [1] is there any other literature on the formulation of compressibleInterFoam? Best Regards! ---- [1] -Miller, S. T., Jasak, H., Boger, D. A., Paterson, E. G., & Nedungadi, A. (2013). A pressure-based, compressible, two-phase flowa finite volume method for underwater explosions. Computers and Fluids*, 87, 132–143. https://doi.org/10.1016/j.compfluid.2013.04.002

December 30, 2021, 01:04		#69
saicharan662000@gmail.com Member hari charan Join Date: Sep 2021 Location: India,hyderabad Posts: 60 Rep Power: 2	Hi Nat, I know the post is old. But can you describe the evolution of step 4 mathematically instead of describing it in words? Thanks in advance.

December 30, 2021, 01:47		#70
saicharan662000@gmail.com Member hari charan Join Date: Sep 2021 Location: India,hyderabad Posts: 60 Rep Power: 2	Hello nat, Can u explain the evolution of step 3 to step 4 using math instead of words. Thanks in advance.

Thread Tools	Search this Thread
Show Printable Version Email this Page	Search this Thread: Advanced Search
Display Modes
Linear Mode Switch to Hybrid Mode Switch to Threaded Mode

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Low Reynolds Number k-epsilon formulation CFX 10.0	Chris	CFX	4	December 7, 2009 23:51
Immersed Boundary Formulation	Rave	Main CFD Forum	0	August 11, 2008 14:55
DPM Steady formulation with collisions	kulwinder	FLUENT	0	May 22, 2004 18:44
energy equation formulation	Pedro	Phoenics	1	July 5, 2001 12:17
Compressible vs. Incompressible formulations	Fernando Velasco Hurtado	Main CFD Forum	3	January 7, 2000 16:51