[cheng1988sjtu]

Hi Foamers,


I've got a little confusion here about the pEqn.flux(). So, what does it really mean, say we want to solve the following momentum equation:

ddt(voidfraction*U) + convection_terms == -voidfraction * grad(p) + Other terms.

As the first step, in OpenFOAM, the following UEqn is solved:

ddt(voidfraction*U) + convection_terms == Other terms.

and discretise this equation into linear system, we could obtain the following form, that is:

A() *[U] = H(). and we can define phi = fvc::interpolate (U) & mesh.Sf(); which is the surface mass flux.

and the "real" velocity field should be :

A() *U = H() -voidfraction * grad(p) .

To satisfy the continuity equation, we have:

ddt(voidfraction) + div(voidfraction*U) = 0

Thus the pressure equation becomes:
pEqn ==
{
laplacian(voidfraction*voidfraction/A(), p) = ddt(voidfraction) + div(voidfraction*phi)
}

then, my question is, pEqn.flux() = SnGrad(p) or pEqn.flux() = (voidfraction*voidfraction/A())*SnGrad(p)

how could we update the phi after we solved the pressure equation, which of the following is right?

(1) phi -= pEqn.flux();

(2) phi -= (1/A()) * (voidfraction*pEqn.flux()/(voidfraction*voidfraction/A()));

this could make a big difference is we want to model the flow through a porous media, and I still don't quite get the exact meaning of pEqn.flux(), though it is stated in the following link that The .flux method includes contributions from all implicit terms of the pEqn (the Laplacian)

http://openfoamwiki.net/index.php/IcoFoam

Could anybody help me with this problem? Thanks!

[santiagomarquezd]

Hi, the pEqn you're assembling is:

fvScalarMatrix pEqn
                (
                    fvm::laplacian(rUA, p) == fvc::div(phi)
                );

which means:



so that the method pEqn.flux() gives you the discrete version of the flux of p, which is



this is not calculated by explicit operators like fvc::grad but directly using the matrix's coefficients.

Regards.

[cheng1988sjtu]

Hi, Santiago,

Thanks very much for your kind reply.

That's great, it really helps me a lot to understand the code.

According to your explaination, then I guess, my option (2) should be the right way to correct the phi, right?

(2) phi -= (1/A()) * (voidfraction*pEqn.flux()/(voidfraction*voidfraction/A()));

Regards,

Charlie

[santiagomarquezd]

Looking at the equation you want to solve it seems to be (2), please check lines 78-93 of

https://github.com/OpenFOAM/OpenFOAM...lerFoam/pEqn.H

which is the code for pEqn.H in twoPhaseEulerFoam. This solver has the kind of momentum equations that you dealing with.

Regards.

[cheng1988sjtu]

Thanks!

My first clue of my doubt also comes from the twoPhaseEulerFoam, I think, right now, I have a clear understanding of how the pressure is corrected, Thanks again!

Regards,

Charlie