CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Running, Solving & CFD

PISO Details

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree1Likes
  • 1 Post By hartinger

Reply
 
LinkBack Thread Tools Display Modes
Old   January 1, 2005, 17:02
Default Hi All, Please excuse my ign
  #1
Dr B.M. Smith (Smith)
Guest
 
Posts: n/a
Hi All,
Please excuse my ignorance, but I was looking up the standard solver member functions to try to understand those implementations of PISO (such as done in sonicFoam, mhdFoam, turbFoam, etc.,...). I still cannot make complete sense of all the field manipulations. For example, the functions,

.A() in "volScalarField A = UEqn.A()"
.H() in "U = UEqn.H()/A"

what are these? And what does "U = UEqn.H()/A" therefore turn the field "U" into, is it still a velocity? The Doxygen material says these are the "central coefficients" A(), and the H() returns the "H operation source". Those definitions are completely mysterious to me. Can someone explain? These things appear to have something to do with obtaining a field "phi" which looks like a flux defined on the mesh faces, which then gets used in the pressure correction (as in a textbook version of PISO). But I can't quite translate the foam code into equivalent symbolic math manipulation.

Also, what is the variable,

UphiCoeff

used in PISO? Where does it come from? It appears in icoFoam and sonicFoam but not in turbFoam PISO loops. I could not see it defined in any of the header files.

Finally (for now) why is the feild "pRef" defined in some solver PISO routines and not others? What is the function of this pRef. I can't see where is is needed for the pressure correction.

Thanks for any help,
Blair.

PS. I may have some additional follow up questions.
  Reply With Quote

Old   January 1, 2005, 17:10
Default I would suggest looking into
  #2
Hrvoje Jasak (Hjasak)
Guest
 
Posts: n/a
I would suggest looking into my Thesis - all the basic discretisation tricks are there. Unfortunately, this forum is a rather clumsy place ot explain, so I'll just give you a gist:

.A() is the central coefficient "field", which the matrix stores separately from the off-diagonal.

For the momentum equation (solving for U), H() would be something like:

H() = \sum_N a_N U_N

+ various non-orthogonal correction and other explicit terms. Have a look at the derivation of the pressure equation and it should be clear.

phi is indeed the flux field (before and after correction).

UphiCoeff is the algorithmic fiddle-factor (you'll have to reverse-engineer that one) :-)

Enjoy,

Hrv
  Reply With Quote

Old   January 1, 2005, 17:37
Default UphiCoeff is not a "fiddle-fa
  #3
Henry Weller (Henry)
Guest
 
Posts: n/a
UphiCoeff is not a "fiddle-factor", it is a coefficient to mix two alternative formulations for the inclusion of the rate-or-change of momentum in the flux prediction. We have found that the two alternatives have different merits and it is best to mix the approaches. I am still looking for the ultimate formulation and will release it as soon as I have it but in the meantime this will have to do. I will write this up in more detail as soon as I have time and we will put it on the web-site.

Actually this formulation with UphiCoeff is included in turbFoam in the most current release of OpenFOAM. In fact it is now used in all the transient codes except those incuding mesh-motion (still haven't found the best appriach for that) but not in the steady-state codes which do not include rate-of-change terms anyway.

The reference pressure pRef is needed for incompressible flows in closed domains for which the pressure level is unknown.
  Reply With Quote

Old   January 1, 2005, 17:52
Default The derivation clearly shows
  #4
Hrvoje Jasak (Hjasak)
Guest
 
Posts: n/a
The derivation clearly shows the factor should be equal to one, right? Otherwise, you get a dependence of the solution to the time-step size.

I understand the two approaches have different merits; however, as you cannot derive the "correct" value of UphiCoeff, with your permission I will carry on calling it the fiddle factor.

I think a more pressing and a more general problem is a lack of formal description of all algorithms used in FOAM - I propose we write them up and publish them as soon as we can.

Incidentally, you may be interested that my student in Croatia has derived some very interesting results regarding the use of second-order time discretisation for moving mesh simulation - it turns out that you need to fiddle with the use of swept volume fluxes in the convection term in order to get the correct result. I can only automate this if I introduce a DDt derivative into the codes, which requires unification of ddt and convection. Any ideas?

The Thesis is being translated in English and if you'd like a look, please let me know.

Hrv
  Reply With Quote

Old   January 1, 2005, 18:31
Default > The derivation clearly show
  #5
Henry Weller (Henry)
Guest
 
Posts: n/a
> The derivation clearly shows the factor should be equal to one, right? Otherwise, you get a dependence of the solution to the time-step size.

It may be 0 or 1 or any value in between.

> as you cannot derive the "correct" value of UphiCoeff, with your permission I will carry on calling it the fiddle factor.

I disagree.
  Reply With Quote

Old   January 1, 2005, 18:41
Default Nope, only one, otherwise you
  #6
Hrvoje Jasak (Hjasak)
Guest
 
Posts: n/a
Nope, only one, otherwise you get a part of div(ddt U_old) with the interpolated old-time velocities instead of fluxes, which is not zero for incompressible flows as it should be. That gives you the numerical error which introduces time-step dependence.

> I disagree.

Any arguments, explanations, clarifications, background or similar? Or just argumentum ad hominem?

Hrv
  Reply With Quote

Old   January 1, 2005, 18:53
Default I do not wish to discuss this
  #7
Henry Weller (Henry)
Guest
 
Posts: n/a
I do not wish to discuss this any further. If you wish to release your own versions of the algorithms in OpenFOAM please do so, it's open-source.
  Reply With Quote

Old   June 15, 2007, 10:23
Default For piso solver we need to men
  #8
New Member
 
abhishek k n
Join Date: Mar 2009
Location: Gothenburg, Sweden
Posts: 16
Rep Power: 8
knabhishek is on a distinguished road
For piso solver we need to mention nNonOrthogonalCorrectors between 0 and 20
What is basis for this range how it is calculated
knabhishek is offline   Reply With Quote

Old   June 15, 2007, 13:17
Default Use checkMesh as a guide. It r
  #9
Senior Member
 
Srinath Madhavan (a.k.a pUl|)
Join Date: Mar 2009
Location: Edmonton, AB, Canada
Posts: 700
Rep Power: 12
msrinath80 is on a distinguished road
Use checkMesh as a guide. It reports the mesh non-orthogonality and skewness. Search the forum for nNonOrthogonalCorrectors and you will find useful pieces of advice everywhere. It's all out there. You need to look a little harder.

My personal opinion: Try to stick to a pure orthogonal mesh unless the cell count becomes unimaginable. I personally found that a 4 million cell case featuring a pure orthogonal grid that maintains a maximum cell aspect ratio of 1:4 solves faster than a 2.5 million case with hanging nodes (i.e. non-conformal grid).
msrinath80 is offline   Reply With Quote

Old   June 19, 2007, 23:44
Default Hello All, I am trying to c
  #10
Member
 
Shaun Cooper
Join Date: Mar 2009
Posts: 54
Rep Power: 8
coops is on a distinguished road
Hello All,

I am trying to combine parts of buoyantFoam and sonicFoam to create a new solver. However, I am unsure of some of the differences in the PISO loop of the two solvers. In buoyantFoam the variable phi is defined as:

phi = (fvc::interpolate(rho*U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi)
- fvc::interpolate(rUA*rho*gh) * fvc::snGrad(rho)*mesh.magSf();

However in sonicFoam the variable defined is phid:

surfaceScalarField phid =
(
(fvc::interpolate(rho*U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi)
)/fvc::interpolate(p);

One thing I don't understand is why the sonicFoam variable (phid) is divided through by interpolate(p). When this term is used within pEqn
it appears as:

+ fvm::div(phid, p, "div(phid,p)")

whereas the buoyantFoam phi variable appears in the pdEqn as:

fvc::div(phi)

Can someone please explain to me why this is the case.

Many Thanks,

Shaun
coops is offline   Reply With Quote

Old   June 20, 2007, 04:26
Default Hi all, I've read the thesy
  #11
lam
Member
 
Hoang-Lam
Join Date: Mar 2009
Location: Lausanne, Switzerland
Posts: 60
Rep Power: 8
lam is on a distinguished road
Hi all,

I've read the thesys of Jasak.

And in the part 3.8, when he speaks about the semi-discretized form of the momentum equation, I wonder where the "laplacian(nu,Grad(U))" is.

Because, he says:
a_p*U_p=H(U)-Grad(p), where H(U)=-sum[a_n*U_n]+U0/deltaT
and div(UU)=a_p*U_p + sum[a_n*U_n]

First I guess that
laplacian(nu,Grad(U))=a_p * U_p + U0/deltaT

But I'm quite not sure.
Hope that somebody can tell me where the laplacian(nu,Grad(U)) term is, in the discretized form.

Thanks in advance,

Lam
lam is offline   Reply With Quote

Old   June 20, 2007, 10:39
Default Hi Shaun, the derivation fo
  #12
Senior Member
 
Markus Hartinger
Join Date: Mar 2009
Posts: 102
Rep Power: 8
hartinger is on a distinguished road
Hi Shaun,

the derivation for the pressure equation in sonicFoam:

(1)conti:
ddt(rho) + div(rho*U) = 0

(2)momentum:
U = H/A - 1/A * grad(P)

(3)equation of state:
rho=psi*p

(2)(3)(4) in (1):
ddt(psi,p) + div( rho*(H/A) - rho/A * grad(P))
<=>
ddt(psi,p) + div(rho*(H/A))- laplacian(rho/A, p)
to make everything implicit
div(rho*(H/A)) = fvm::div(rho*(H/A)/p, p)
or you could use the equation of state as in sonicFoamAutoMotion
div(rho*(H/A)) = fvm::div(psi*(H/A), p)

cu
markus
hartinger is offline   Reply With Quote

Old   June 22, 2007, 02:10
Default Hi All, Thanks Markus, that
  #13
Member
 
Shaun Cooper
Join Date: Mar 2009
Posts: 54
Rep Power: 8
coops is on a distinguished road
Hi All,

Thanks Markus, that makes that part of the derivation of the equation fro pressure a lot clearer. However, I am still having a little trouble with the formulation of phid (sonicFoam, phi in buoyantFoam) and also:

phi = pEqn.flux();

and how that is then used to correct the velocity field. I am reading Hrv's thesis and Computational Methods for Fluid Dynamics by Ferziger and Peric so hopefully it will become clear to me this time (I have looked at these texts before ;-) ).

Any assistance would be much appreciated,

Shaun
coops is offline   Reply With Quote

Old   June 22, 2007, 08:27
Default Hi Shaun, the derivation fo
  #14
Senior Member
 
Markus Hartinger
Join Date: Mar 2009
Posts: 102
Rep Power: 8
hartinger is on a distinguished road
Hi Shaun,

the derivation for the pressure equation in buoyant Foam:

(1) conti:
ddt(rho) + div(rho*U) = 0

(2) momentum:
U = H/A - 1/A * grad(P)

(3) equation of state
d(rho) = d(psi(p)*p)

(4) pressure
p = pd + rho*gh + pref

-----------------------------
(3),(4) in (1), time derivative:
ddt(rho) = ddt(psi * (pd + rho*gh + pref))
= ddt(psi, pd) + ddt(psi)*pref) + ddt(psi*rho)*gh

-----------------------------
(2)in (1), convective term:
div(rho*U) = div(rho * H/A) - div(rho/A*grad(p))
with (4)
div(rho/A*grad(p)) = div(rho/A*grad(pd)) + div(rho/A*grad(rho*gh+pref))

div(rho/A*grad(pd)) = laplacian(rho/A, p)

gh and pref are constant:
div(rho/A*grad(rho*gh+pref)) = div(rho/A*grad(rho)) * gh

generally, the laplacian is discretised as:
laplacian(a, b) = div(a_f * magSf * snGrad(b))
so, as in the code:
div(rho/A*grad(rho)) * gh = div(interpolate(rUA*rho*gh)*snGrad(rho)*magSf())

after solving for pd the implicit contribution to phi has to be added
phi += pdEqn.flux()

the velocity is corrected using equation (2) keeping in mind that p = pd + rho*gh + pref, it follows:
U = H/A - 1/A * grad(pd + rho*gh + pref)
pref and gh are constant, so:
U = H/A - 1/A * (grad(pd) + grad(rho)*gh

in sonicFoam
the ".flux()" function returns the influence of all implicit terms on the flux. In this case everything is implicit, so you can assign it like that
phi = pEqn.flux();

cu
markus
cheng1988sjtu likes this.
hartinger is offline   Reply With Quote

Old   June 22, 2007, 09:22
Default Hi Markus, Thanks for your
  #15
sek
Member
 
Sung-Eun Kim
Join Date: Mar 2009
Posts: 76
Rep Power: 8
sek is on a distinguished road
Hi Markus,

Thanks for your explanation. One more question on phid. Where does the term

fvc:ddtPhiCorr(rUA, rho, U, phi)
comes from?

=========================================
surfaceScalarField phid =
(
(fvc::interpolate(rho*U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi)
)/fvc::interpolate(p);

Sung-Eun
sek is offline   Reply With Quote

Old   June 22, 2007, 11:10
Default look here, http://openfoamwik
  #16
Senior Member
 
Markus Hartinger
Join Date: Mar 2009
Posts: 102
Rep Power: 8
hartinger is on a distinguished road
look here,
http://openfoamwiki.net/index.php/IcoFoam
hartinger is offline   Reply With Quote

Old   July 2, 2009, 14:41
Default
  #17
Senior Member
 
isabel
Join Date: Apr 2009
Location: Spain
Posts: 171
Rep Power: 8
isabel is on a distinguished road
Does anybody know how to do the derivation for the pressure equation for InterFoam?
Please correct me if I am wrong:


(1) continuity:
div(rho*U) = 0


(2) momentum:
U = H/A -1/A*grad(pd)+sigma*k*n/A


(3) equation of state
rho = psi*p


(4) pressure
p = pd + rho*gh


------------------------
(3) and (4) in (1):
div(H/A – grad(P)/A + sigma*k*n/A) = 0


div(H/A) – div(1/A* div(pd+rho*gh)) + div(sigma*k*n/A) = 0


div(H/A) – div(1/A* div(pd) + 1/A*div(rho*gh)) + div(sigma*k*n/A) = 0


div(H/A + sigma*k*n/A - div(rho*gh)/A) - laplacian(pd/A) = 0


div(phi) - laplacian(rUAf*grad(pd)) = 0


div(phi) = laplacian(rUAf*grad(pd))




being:


phi = phiu + rUAF*(sigma*k*n + gh*grad(rho))
isabel is offline   Reply With Quote

Old   November 17, 2009, 12:31
Default pEqn in cavitatingFoam
  #18
Member
 
Hamed Aghajani
Join Date: Mar 2009
Location: London, UK
Posts: 77
Rep Power: 8
haghajani is on a distinguished road
Hello All,
Regarding the pEqn in cavitatingFoam,
I think because of the continuity
(1) ddt(rho)+div(rho*U)=0;
(2) ddt (rho, U) + div (phi,U)-laplacian(muf,U) = -grad(p);
(3) rho = p*psi + (1-gamma) rhol0 -((gamma * psiv+(1-gamma)*psil)-psi)*psat, where rhol0 = rholsat-psat+psil

Then, the pressure equation should be sth like,
ddt(psi*p)-(rhol0+(psil-psiv)psat)*ddt(gamma)-psat*ddt(psi)+div(rho*U)=0

The implementation in OF shows two extra terms which I can not understand where they came from

Thank you in advance,
Hamed
haghajani is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Problems understanding some piso details tehache OpenFOAM Running, Solving & CFD 3 July 27, 2007 06:02
details prasad Main CFD Forum 0 March 26, 2007 23:25
Some details styoung317 Main CFD Forum 0 June 24, 2005 04:49
details of cfd madhusudan Main CFD Forum 1 April 26, 2005 09:44
Details of Y+ joseph Main CFD Forum 4 July 12, 2001 15:44


All times are GMT -4. The time now is 15:31.