CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (https://www.cfd-online.com/Forums/openfoam-solving/)
-   -   pisoFoam laplace equation (https://www.cfd-online.com/Forums/openfoam-solving/87804-pisofoam-laplace-equation.html)

fisch April 29, 2011 08:05

pisoFoam laplace equation
 
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 April 30, 2011 02:40

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 April 30, 2011 08:37

Quote:

Originally Posted by alberto (Post 305670)
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 April 30, 2011 16:57

Quote:

Originally Posted by fisch (Post 305698)
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.

Quote:

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 May 1, 2011 06:47

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 :D

Thanks,
rupert

alberto May 1, 2011 18:19

Quote:

Originally Posted by fisch (Post 305761)
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 May 2, 2011 02:58

Quote:

Originally Posted by alberto (Post 305726)
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 May 2, 2011 03:22

Quote:

Originally Posted by fisch (Post 305857)
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 May 2, 2011 05:04

thank you very much

star shower October 29, 2011 08:14

learn a lot!
 
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 July 16, 2013 22:13

Quote:

Originally Posted by star shower (Post 329976)
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


All times are GMT -4. The time now is 10:23.