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.

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

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.

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

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

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

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.

[knabhishek]

For piso solver we need to mention nNonOrthogonalCorrectors between 0 and 20
What is basis for this range how it is calculated

[msrinath80]

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).

[coops]

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

[lam]

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

[hartinger]

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

[coops]

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

[hartinger]

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

[sek]

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

[hartinger]

look here,
http://openfoamwiki.net/index.php/IcoFoam

[isabel]

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))

[haghajani]

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