[nishant_hull]

Hi Foamer

Do anyone have got any idea, how these two implementation different:-
1st implementation-
fvScalarMatrix pEqn
(
fvm::ddt(psi, p)
+ fvm::div(phid, p)
- fvm::laplacian(rho*rUA, p)
);
where phid is:-

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

2nd implementation:-
fvScalarMatrix pEqn
(
fvm::ddt(psi, p)
+ fvc::div(phi)
- fvm::laplacian(rho*rUA, p)
);
where phi is:-
phi = fvc::interpolate(rho*U) & mesh.Sf()

Any hint or suggestion will be appreciated.

Apart from that, I am also struggling with modification of sonicFoam solver to work for subsonic flow. I have changed some bits of code and presented it here..
http://www.cfd-online.com/Forums/ope...tml#post227455
could you please take a look on it and give some hints?

regards,

Nishant

[mkraposhin]

1st implementation-
fvScalarMatrix pEqn
(
fvm::ddt(psi, p)
+ fvm::div(phid, p)
- fvm::laplacian(rho*rUA, p)
);
where phid is:-

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

it's a
convective-diffusive implicit form of pressure equation.
rho = rhoOld + rho'
because rho = psi*p -> rho = rho'
ddt(rho) + div(rho*U) = ddt(rho') + div(rho' * U) + div(rho/A*grad(p))

div(rho' * U) = div(psi*p*U), combining psi*U = phid we get as in OpenFOAM:

ddt(psi,p) + div(phid,p) + div(rho/A*grad(p))

2nd implementation:-
fvScalarMatrix pEqn
(
fvm::ddt(psi, p)
+ fvc::div(phi)
- fvm::laplacian(rho*rUA, p)
);
where phi is:-
phi = fvc::interpolate(rho*U) & mesh.Sf()

Duffusive equation for pressure (convective correction) is treated explicitly

[nishant_hull]

Thanks for your reply mkraposhin,

I guess the second implementation is what we need for PISO implementation for subsonic type of flow. and the first implementation is mostly used for supersonic or transonic flows. am i right?

regards,

Nishant

[mkraposhin]

Yes, you are right (at most).

First implementation (convective-diffusive) for flows with Ma >= 0.2-0.3

Second implementation (diffusive) for flows with Ma < 0.2-0.3

[nishant_hull]

Thanks mkraposhin,

I have actually modified the sonicFoam solver for low mach number flow. I am trying to simulate a pulse ( v=0.05*sin(2*pi*t) )across the pipe. As you can see that the velocity fluctuations are very small in this case. (however the pulse velocity is actually speed of sound.)
This solver is actually solving the compressible equations but the pulse is not reaching the ther end of the long pipe and the pulse is dying somewhere in the middle.

have a look on my modified solver. see the attached file.

regards,

Nishant

[mkraposhin]

hi, nishant_hull

i found one wrong line in your code:

phi = pEqn.flux().

for diffusive equation (where fvc::div(phi)), i think, that phi should be corrected in a such way:

phi += pEqn.flux()

also, how did you init phi before solving for p?

phi = (fvc::interpolate(rho*U) & mesh.Sf()) + fvc::ddtPhiCorr(...)?

[mkraposhin]

In attached file you will find corrected lines, marked with red color

[nishant_hull]

Thanks mkraposhin,

I think phi is already defined for compressible solvers and added in some header file as:

phi = fvc::interpolate(rho*U) & mesh.Sf()

This variable is also used by the header rhoEqn.H and declaration of the UEqn.

I think you are right there in updating the phi += pEqn.flux()
Thanks for this suggestion. I was still wondering why this is??

However I have seen a subsonic solver (XiFoam) using the line:
DpDt = fvc:Dt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);

I was wondering what this lines does there??

regards,
Nishant

[mkraposhin]

well,

phi is defined, but it needs to be updated.

when we write

phi = (fvc::interpolate(rho * U) & mesh.Sf) + fvc::ddtPhiCorr(...) we update mass flux before solving pressure equation. This value of phi will be used as starting (unsatisfied to continuity, but satisfied to momentum equation).

After solving pressure equation and applying flux, we got phi, which satisfies continuity equation.

In OpenFOAM, variable phi is a key to all problems - momentum equation is linearized using values of phi from last time step.

fvm::div(phi,U) mean div (rho*U,U), where rho*U = phi

about DpDt:

ddt(rho*e) + div(rho*U*e) = .... (1)

e = u + (w^2)/2 (2)

u = h - p/rho (3)

after neglecting w^2/2 and substituting (3) into (1)

we get (h - enthalpy)

ddt(rho*h) - ddt(p) + div(rho*U*h) - div(U*p) = ...
a b c d

ddt(rho*h) + div(rho*U*h) - (ddt(p) + div(U*p))= ...
e f g

(g) = D/Dt (p)

ddt(rho*h) + div(rho*U*h) - DpDt = ...


is it clear?

[nishant_hull]

Oh, I see. This was very much helpful.
I can understand those DpDt equation as well.
Thanks once again.

However I was having a look on another solver sonicLiquidFoam and I observed that the pressure equation there is solved like this:
surfaceScalarField phid
(
"phid",
psi*
(
(fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi)
)
);

phi = (rhoO/psi)*phid;

fvScalarMatrix pEqn
(
fvm::ddt(psi, p)
+ fvc::div(phi)
+ fvm::div(phid, p)
- fvm::laplacian(rho*rUA, p)
);

I was windering how this application is different from the one I am trying to implement? and would it be a better option for low mach number flow?

regards,

Nishant

[mkraposhin]

as a rule in many OpenFOAM compressible solvers accepted that density and pressure relates linearly: rho = psi * p. But for sonicLiquidFoam pressure and density relates as rho = rho0 + psi*p.

General relationship: psi = d(rho)/d(p)=1/(a^2). For ideal gas: psi = rho/p, for liquids rho = rho0 + psi*p

[nishant_hull]

Thanks Matvej,
Things are pretty clear now. Your comments and suggestions were very helpful.
I was trying to run solver with phi suggested by you as:
phi = (fvc::interpolate(rho*U) & mesh.Sf()) + fvc::ddtPhiCorr(rUA, rho, U, phi); but it was showing some dimensional inconsistency. I therefore opted for this option:

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

Hope it should be fine.

regards,
Nishant

[mkraposhin]

Yes, it is alright.

This due to difference between OF1.4 and OF1.5

[nishant_hull]

Thanks Matvej,

I would like to ask you something about the updating rho value in coodles solver through thermo->rho(). In fact I am dealing with very small pressure fluctuations and therefore I want to modify the coodles solver by using a relative pressure value instead of absolute. I want specifying a reference pressure value and therefore wanted to know how the thermo->rho() is modifying the rho.

I guess, thermo->rho() must have been using the pressure value to modify rho. (if i m not wrong then from u = h - p/rho, the rho = p/h, neglecting u) but I am not sure abt the implementation and where it is done.

regards,

Nishant

[mkraposhin]

Well, you are wrong about thinking way the pressure is updated in OpenFOAM

rho = psi * p, for ideal gas psi = d(rho)/d(p) = 1/(RT)

I think, you need to define constant pressure level (like it is done in buoyantFoam) - pRef, then p = p_0 + pRef.

pRef - is a constant value, and p_0 is deviation from pRef.

If you are not sure, please, post solver code here, and i'll try to fix and explain it

[nishant_hull]

Thanks for the reply Matvej.
Sorry for my late reply. and yes, I have seen the implementation of the thermophysical model recently.
The solver is the same Openfoam coodle solver. (or you can also take any simple compressible solver except sonicfoam, to explain it.) I think it will be difficult to attach that here (because of limit on attachment). but if you want me to send the solver code then I can send you at your email. (and i will need your email for that)

regards,

Nishant

[nishant_hull]

see the attached coodle solver of 1.5.
Guess I can upload bigger files now i was having problems before.

regards,

Nishant

[nishant_hull]

Hi Matvej

I have got the buoyantFoam solver from the old 1.4 version. I have seen the implementation of reference pressure there. I will try to implement this with my code now. hopefully it should work.
I will update you in this regard. Thanks once again,

regards,

Nishant

[mkraposhin]

Hi, nishant_hull

I'm also sorry for late reply - too much current work

I don't know, why are you trying to write your own solver - all needed models are implemented in OpenFOAM.

Continuity equation:

d/dt(rho) + div(rho*U) = 0

Momentum conservation equation:

d/dt(rho U) + div(rho U U) - div(R) = -grad(p) + rho*g

R = mu (grad(U) + grad(U)^T))


d/dt(rho U) + div(rho U U) - div(R) = RHS (Right Hand Side)
and
-grad(p) + rho*g = LHS (Left Hand Side)

we get (for example, g=(0, -9.81, 0))

RHS = - grad(p) + rho*g

p = pd + rho*(g & r) + pRef,
where
g = (gx; gy; gz), r = (X, Y, Z),
& - inner product ( a&b = ax*bx + ay*by + az*bz)

we get
-grad(p) + rho*g = -grad(pd + rho*(g & r) + pRef) + rho*g=
-grad(pd) - grad(rho)*(g&r)

g&r in OpenFOAM is 'gh' field

in OpenFOAM <= 1.4.1 we are solving for field pd = p - rho*(g&r) - pRef,
in OpenFOAM >= 1.5 we are solving for field p.

Take a look at BC's (buoyantPressure) at wall:

RHS == 0 == -grad(p) + rho*g --> grad(p) = rho*g, or, in terms of relative pressure pd grad(pd) = -grad(rho)*(g & r)

I hope, this is clear (if not - due to my lack of knowledge of English, please write)

[nishant_hull]

i think i did some stupid mistake before.