[mcchouffe]

Dear OpenFOAM userbase,

I'm looking somwhat deeper into interFoam.
After reeding the gravity vector g from the environmentalProperties dictionary, a surfaceScalarField ghf is constructed using:

surfaceScalarField ghf("gh", g & mesh.Cf());

What exactly is this constructor doing? What exactly is being stored in the cell-face positions?

thanks in advance for all answers

kind regards,
Laika,
still orbiting
laika@inorbit.com

[olwi]

Hi,

According to the Programmer's Guide, mesh.Cf() is the cartesian coordinates of the face centres (Table 2.1, Sect. 2.3.2), and "&" is the inner product ("dot-product") operator (Table 1.2, Sect. 1.4.1).

I don't know interFoam, so I stop there...

/Ola

[mcchouffe]

Hi,

Thanks Olga.
I just trying to find out where that inner product of the two vectors is coming from. Why is the gravity hidden in such a scalar field, and why is it this scalar field that is being used in the solver (Ueqn.H). I would expect a vectorial Body force to be given to the momentum equation, and not a scalar field...

anyone who knows the heart of interFoam who can help me out?

kind regards,
Laika,
still orbiting
laika@inorbit.com

[hjasak]

Easy. Momentum equation should contain the gravity with the variable density due to free surface. However, this will produce a linear pressure variation even in fluid at rest + cause decoupling problems in a segregated cell-centred solver like interFoam.

The cure is to rip out the gravity from the pressure and solve for piezometric pressure instead. This moves the gravity term in the pressure equation as grad(rho g h), with the squeezed computational molecule (hence the face field instead of cell centres). On the negative side, you lose the actual pressure distribution (have a look at the pressure field) and it is not trivial to recover it.

Enjoy,

Hrv

[mcchouffe]

Aha, thanks for shedding a light on this.
Very interesting!

I'd say this counts for other Body Forces as well, and hence if I want to implement in a VOF calculation a steady growing linear acceleration as a body force, I just need to dynamically adapt the g and construct the ghf with (a+g) rather than g. I don't have to add a body force in the momentum equation.

Now the thing is I do need the pressure field. Am I right when I say I just have to subtract (g+a)*h (as calculated by the inner product to construct ghf) from the calculated pressure field to obtain the static pressure field?

kind regards,
Laika
still orbiting

[hjasak]

Well, you have to integrate the static component properly properly - the issue is that at a point in a domain, absolute pressure depends on the height of the water column above it (with respect to the acceleration vector), which in turn pulls in the gamma distribution above the point.

Assembling the integral is not trivial so you have to be clever about it - play with the equations :-)

Keep orbiting,

Hrv

[maka]

Hi Hrv,

+ cause decoupling problems in a segregated cell-centred solver like interFoam.

Can you say more about that? I'm familiar with pEq.H in interfoam and PISO (Isaa85 and Rusche2002). I studied interFoam to learn about how to include body forces properly but I'm not working with multiphase flow. Thanks.

Best regards,
Maka.

[hjasak]

Are you into understanding precisely that is going on or you just want to do the implementation that will not cause you trouble?

In short, the computational molecule supporting the pressure equation, which means that you should convert the volume force into a div(...) of something. The rest follows naturally, ie. you will add this bit as the flux into the pressure equation, meaning that a component of the pressure has been eliminated.

If you need the complete pressure, some post-processing is necessary... but that's another story.

Enjoy,

Hrv

[maka]

yes, I'm actually interested in understanding what is going on, so I have the freedom to modify it. I'm familiar with finite volume method and fractional step, PISO and have been using OpenFOAM for a while. I say that just to lay some background for the discussion.

(a) I needed understand why the addition of a rho*g source term would cause decoupling problems in a segregated cell-centered solver? Usually one faces a pressure oscillation problem, when a Rie-Chow like interpolation is NOT used with collocated grid. Similar problem may happen for other source terms that are gradients but rho*g does not include any gradient.

(b) another issue, comparing the description of flux in Jasak PhD thesis,
[eq 3.144] F=S.Uf=S.((H(U)/aP)_f - (1/aP)_f (grad p)_f);
with Rusche2002 (attached) page 110
[eq 3.33] ...
"phi flux is a result of the pressure-velocity solution procedure outlined in Section 3.2.4 and is not evaluated by taking the dot product of the face area vector and the face interpolate of
the velocity as this would not obey continuity".
(b1) Is this because the U vector was modified by a predictor step that does not guarantee that it is divergent free.
(b2) What if there was no one and we were using the U_old which should be divergent free, after all that was the objective of the velocity-pressure coupling algorithm?

more, in icoFoam (for simplicity I will not consider interFoam version now) V1.3:

U = rUA*UEqn.H();
phi = (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, U, phi);

putting a side ddtPhiCorr which I do not understand yet, the phi does not include a term with pressure gradient like F in Jasak.
(b3) As a result phi is built with a velocity field that does not satisfy momentum, especially if no predictor was used? why? At least, in Jasak thesis, the Uf is some kind of prediction of U based on momentum equation.

(c) Also in OpenFOAM, the divergence form insists on having the flux as separate variable from velocity vector. I looked why in the forum, the statements I got was related to the flux being divergent free. But as I understand this divergent free flux has to corresponds well to at least an old velocity field that was divergent free. phi of icoFoam in the absence of a predictor step is not, since it does not include any pressure gradient effects (old or updated). Did I miss something here?

Sorry for making this a long post but, but after using OpenFOAM for some time, I know that lots of the things implemented are implemented in this way for a good reason and I just wanted to understand the background so, I feel free to modify the code. Most of the above questions are related to one theme, where in the derivation one can find such information? Thanks for your help.

Best regards,
Maka.

[maka]

the attachment is here.

[maka]

too big to be upload, here is the reference:

Rusche, H. Computational fluid dynamics of dispersed two-phase flows at high phase fractions Imperial College of Science, Technology & Medicine, Department of Mechanical Engineering, 2002.

send me and email if you need the PDF file.

Best regards,
Maka.

[maka]

the link:
Rusche2002