[fisch]

Hello,

i'm currently reading jasaks thesis and have a look inside the pisofoam solver.
Here is one thing i didnt understand in the following part of the solver:

U = rUA*UEqn.H();
phi = (fvc::interpolate(U) & mesh.Sf())+(maybe something)
...
fvScalarMatrix pEqn
(
fvm::laplacian(rUA, p) == fvc::div(phi)
);

Is the equation for phi representing the hole (2nd) r.h.s. of equation 3.141 in the thesis (page 146)? I thought so, because we are evaluating here the UEqn.H/a_p on the Surface multiplied by the Area.
Why are we then taking the divergence of phi in the laplace equation? This would not match the formula 3.141 or 3.143...


Can anybody help me and/or tell me my missinterpretations?
This would be great.

Thanks for any advice,

rupert

[alberto]

Hi, think to discretizing the momentum equation as

A*U = H - grad(p),

then

U = H/A - 1/A * grad(p)

The face flux is then

phi = U_f . S = (H/A)_f . S - (1/A)_f |S| snGrad(p)

OpenFOAM first computes

phi = (H/A)_f . S,

then solves the pressure equation obtained imposing

div(phi) = 0,

which leads to

laplacian (1/A*p) = div(phi)

Once this equation is solved, the flux is corrected using the second part of the flux definition.

[fisch]

Originally Posted by alberto

Hi, think to discretizing the momentum equation as

A*U = H - grad(p),

then

U = H/A - 1/A * grad(p)

The face flux is then

phi = U_f . S = (H/A)_f . S - (1/A)_f |S| snGrad(p)

OpenFOAM first computes

phi = (H/A)_f . S,

then solves the pressure equation obtained imposing

div(phi) = 0,

which leads to

laplacian (1/A*p) = div(phi)

Once this equation is solved, the flux is corrected using the second part of the flux definition.

Hello Alberto, thank you for the reply.

Maybe i have problems understanding the commands...

I understand that taking the divergence of equation phi = U_f . S = (H/A)_f . S - (1/A)_f |S| snGrad(p) removes the left hand side and the right hand side = 0 is building then the poisson equation... right?

But I thought that (1/A)_f |S| snGrad(p) == div[(1/A) * grad(p)] == laplace (1/A * p)
Herein snGrad(p) is the gradient of p at the surfaces.
is this not right?

Having a look on the equation for pressure on http://openfoamwiki.net/index.php/Th...hm_in_OpenFOAM at chapter 2 (the last equation) this would lead (for me) to a pressure equation like:
laplacian (1/A*p) = (H/A)_f . S = phi
(and not div(phi))

Do you understand my problem? where is my misinterpretation???


Thank you a lot,
rupert

[alberto]

Originally Posted by fisch

But I thought that (1/A)_f |S| snGrad(p) == div[(1/A) * grad(p)] == laplace (1/A * p)
Herein snGrad(p) is the gradient of p at the surfaces.
is this not right?

Yes, since you basically dot with the surface vector.

Having a look on the equation for pressure on http://openfoamwiki.net/index.php/Th...hm_in_OpenFOAM at chapter 2 (the last equation) this would lead (for me) to a pressure equation like:
laplacian (1/A*p) = (H/A)_f . S = phi
(and not div(phi))

If you define

phi = (H/A)_f . S - (1/A)_f |S| snGrad(p)

you have

div(phi) = 0

implies

div((H/A)_f . S) = laplacian(1/A*p),

which is exactly

div(phi) = laplacian(1/A*p)

Best,

[fisch]

Good morning,
i hope that i'm not bugging you, but:

you confirmed that
a) (1/A)_f |S| snGrad(p) == div[(1/A) * grad(p)] == laplace (1/A * p)
right?

If a) works it will make the phi equation from
phi = (H/A)_f - (1/A)_f |S| snGrad(p)
to
b) phi = (H/A)_f - laplace (1/A*p)

div(b) leads to
0 = div((H/A)_f) - div(laplace (1/A*p))
or not?

I think that my problem lies in a)
If (1/A)_f |S| snGrad(p) would be equal to (1/A) * grad(p) this would solve my problems... Is this the case? This would mean that
snGrad(p) == grad(p) (defined at the cellcenter) * normalvector (at the surfaces)

Hope you can finally help me

Thanks,
rupert

[alberto]

Originally Posted by fisch

Good morning,
i hope that i'm not bugging you, but:

you confirmed that
a) (1/A)_f |S| snGrad(p) == div[(1/A) * grad(p)] == laplace (1/A * p)
right?

No. How can

div(1/A grad(p))

be equal to

(1/A)_f |S| snGrad(p) ?

You have to think to the velocity predictor

U = H/A - 1/A grad(p),

interpolate it on cell faces, and obtain U_f. At that point you dot it with the surface vector, and, at that point, you have the snGrad(p).

Best,

[fisch]

Originally Posted by alberto

If you define
phi = (H/A)_f - (1/A)_f |S| snGrad(p)
you have
div(phi) = 0
implies
div((H/A)_f) = laplacian(1/A*p),
which is exactly
div(phi) = laplacian(1/A*p)

Hi,

using your definition for phi: "phi = (H/A)_f - (1/A)_f |S| snGrad(p)"

Div of this equation leads to " 0 = div((H/A)_f) - div( (1/A)_f) |S| snGrad(p))"

How can you say here that div( (1/A)_f) |S| snGrad(p)) is equal to laplacian(1/A*p)???

[alberto]

Originally Posted by fisch

Hi,

using your definition for phi: "phi = (H/A)_f - (1/A)_f |S| snGrad(p)"

Div of this equation leads to " 0 = div((H/A)_f) - div( (1/A)_f) |S| snGrad(p))"

How can you say here that div( (1/A)_f) |S| snGrad(p)) is equal to laplacian(1/A*p)???

By definition, the Laplacian of a function defined on a Euclidean space is the divergence of the gradient of the same function. The surface-normal gradient originates from the dot product of the interpolated velocity on cell faces with the face area vector.

The only element that probably can cause confusion is the norm of the surface, which is not in the code, since operators take care of it.

Best,

[fisch]

thank you very much

[star shower]

I learn a lot! Thank you!
For me, the relationship of div(phi)=0 is doubted. Now I know, phi is a scalar, and its div is zero!

[sharonyue]

Originally Posted by star shower

I learn a lot! Thank you!
For me, the relationship of div(phi)=0 is doubted. Now I know, phi is a scalar, and its div is zero!

I dont think its because phi is a scalar then div(phi) =0, But Im also confused it too. http://www.cfd-online.com/Forums/ope...tml#post440172