compressibilty effect in alphaEqn.H
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 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 ? |
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: so I use Code:
alpha1Su = alpha1/rho1*fvc::DDt(phi1,rho1) |
i calculate the source terms like below:
Code:
vDotP1 = fvc::ddt(rho1) + fvc::div(phi,rho1) - rho1*fvc:: div(phi); |
And you don't experience any problems with mass balance? I'll try your implementation, thanks a lot.
|
in your implementation, vDotP1 is equal to fvc::DDt(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? |
1 Attachment(s)
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
|
look at here
|
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 |
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 |
Hi Nima.
Explicit source term in pEqn means div(U), right? Do you have something like Code:
// update density |
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()) Code:
while (pimple.correctNonOrthogonal()) Code:
for (int aCorr=0; aCorr<nAlphaCorr; aCorr++) Code:
div(phi,rho1) Gauss linear; |
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
|
Hi Nima,
is my procedure the same as yours, or do you have different pressure correction eqn.? Best regards, Illya |
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 |
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 |
i solved rhoEqn once,
when I solve alphaEqn, rhoPhi is getting update, then i calculate total rho from rhoPhi from continuity Eqn (look compressibleInterFoam) |
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 |
Quote:
Quote:
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:
Quote:
Quote:
|
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? |
Hi Nima,
Update: I've added 2 rhoEqns, for rho1 and rho2 at the end of my time step. Now I force the continuity of each phase. The drawback is the the EOS is not exactly fulfilled in some interface cells, that caused problems earlier, but I can live with it, since the continuity is more important. Thanks a lot for the discussion, it helped me a lot. If if you're interested, you can still give me your solver some day, and I could test it with my case. I have a high thermal expansivity of both phases and it may be interesting to know, if your solver still performs well. It must not be now, some day you wish. Best wishes and a nice weeked! Illya |
Any references on your work?
Hi,
first need to thank both of you for the valuable discussion and I'm really sorry as I'm trying to refer your work 6 years ago. I'm trying to simulate high pressure gas atomization problem and molten metal is atomizing due to the high pressure gas. The gas inlet is a pressure boundary condition and pressure is in the range of 15-30 atm. However, I tried compressibleInterFoam in OF40 and it crashes with negative temperatures. I know it's due to the sudden expansion of the gas. I tried different solver conditions but still getting the same error. I assume it's a result of taking density as a function of pressure only. I trying to workout the math to obtain the effect of rho(P,T), following the method described in "A pressure based compressible, two phase flow finite volume method for underwater explotions" by Miller et al. I'm still trying to understand the PEqn.H. However, seems like you to have made some progress on this. Can you please give me some references of your work so that I can follow. Thank you. Kalpana |
All times are GMT -4. The time now is 08:08. |