
[Sponsors] 
January 1, 2005, 17:02 
Hi All,
Please excuse my ign

#1 
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. 

January 1, 2005, 17:10 
I would suggest looking into

#2 
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 offdiagonal. For the momentum equation (solving for U), H() would be something like: H() = \sum_N a_N U_N + various nonorthogonal 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 fiddlefactor (you'll have to reverseengineer that one) :) Enjoy, Hrv 

January 1, 2005, 17:37 
UphiCoeff is not a "fiddlefa

#3 
Guest
Posts: n/a

UphiCoeff is not a "fiddlefactor", it is a coefficient to mix two alternative formulations for the inclusion of the rateorchange 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 website.
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 meshmotion (still haven't found the best appriach for that) but not in the steadystate codes which do not include rateofchange terms anyway. The reference pressure pRef is needed for incompressible flows in closed domains for which the pressure level is unknown. 

January 1, 2005, 17:52 
The derivation clearly shows

#4 
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 timestep 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 secondorder 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 

January 1, 2005, 18:31 
> The derivation clearly show

#5 
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 timestep 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. 

January 1, 2005, 18:41 
Nope, only one, otherwise you

#6 
Guest
Posts: n/a

Nope, only one, otherwise you get a part of div(ddt U_old) with the interpolated oldtime velocities instead of fluxes, which is not zero for incompressible flows as it should be. That gives you the numerical error which introduces timestep dependence.
> I disagree. Any arguments, explanations, clarifications, background or similar? Or just argumentum ad hominem? Hrv 

January 1, 2005, 18:53 
I do not wish to discuss this

#7 
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 opensource.


June 15, 2007, 10:23 
For piso solver we need to men

#8 
New Member
abhishek k n
Join Date: Mar 2009
Location: Gothenburg, Sweden
Posts: 16
Rep Power: 9 
For piso solver we need to mention nNonOrthogonalCorrectors between 0 and 20
What is basis for this range how it is calculated 

June 15, 2007, 13:17 
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: 703
Rep Power: 13 
Use checkMesh as a guide. It reports the mesh nonorthogonality 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. nonconformal grid). 

June 19, 2007, 23:44 
Hello All,
I am trying to c

#10 
Member
Shaun Cooper
Join Date: Mar 2009
Posts: 54
Rep Power: 9 
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 

June 20, 2007, 04:26 
Hi all,
I've read the thesy

#11 
Member
HoangLam
Join Date: Mar 2009
Location: Lausanne, Switzerland
Posts: 60
Rep Power: 9 
Hi all,
I've read the thesys of Jasak. And in the part 3.8, when he speaks about the semidiscretized 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 

June 20, 2007, 10:39 
Hi Shaun,
the derivation fo

#12 
Senior Member
Markus Hartinger
Join Date: Mar 2009
Posts: 102
Rep Power: 9 
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 

June 22, 2007, 02:10 
Hi All,
Thanks Markus, that

#13 
Member
Shaun Cooper
Join Date: Mar 2009
Posts: 54
Rep Power: 9 
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 

June 22, 2007, 08:27 
Hi Shaun,
the derivation fo

#14 
Senior Member
Markus Hartinger
Join Date: Mar 2009
Posts: 102
Rep Power: 9 
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 

June 22, 2007, 09:22 
Hi Markus,
Thanks for your

#15 
Member
SungEun Kim
Join Date: Mar 2009
Posts: 76
Rep Power: 9 
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); SungEun 

June 22, 2007, 11:10 
look here,
http://openfoamwik

#16 
Senior Member
Markus Hartinger
Join Date: Mar 2009
Posts: 102
Rep Power: 9 

July 2, 2009, 14:41 

#17 
Senior Member
isabel
Join Date: Apr 2009
Location: Spain
Posts: 171
Rep Power: 9 
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)) 

November 17, 2009, 12:31 
pEqn in cavitatingFoam

#18 
Member
Hamed Aghajani
Join Date: Mar 2009
Location: London, UK
Posts: 77
Rep Power: 9 
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 + (1gamma) rhol0 ((gamma * psiv+(1gamma)*psil)psi)*psat, where rhol0 = rholsatpsat+psil Then, the pressure equation should be sth like, ddt(psi*p)(rhol0+(psilpsiv)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 

Thread Tools  
Display Modes  


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 