Different Implementation in solving pEqn
 

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

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

Duffusive equation for pressure (convective correction) is treated explicitly

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

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

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

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

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

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::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);

I was wondering what this lines does there??

regards,
Nishant

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?

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

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

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

Yes, it is alright.

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

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

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

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

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

regards,

Nishant

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

Buoyancy in OpenFOAM
 

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)

i think i did some stupid mistake before.

Hi Matvej

I have advanced a bit ahead with my last problem. The previous error is sorted out. it was because of the difference in meshing and boundary condition.

However I am still wondering with the question of different implementation of p-equation. I mean..
when we use phi in pEqn (), we use
phi += pEqn.flux();
And when we use phid etc then we use
phi = pEqn.flux();

I am still wondering why?? Any hint please?

regards,
Nishant

Hi, I'm very sorry for late reply (too much work).
Anyway, if this answer is still needed:

when we update velocity (phi = velocity*density*area) in diffusion pressure equation, we mean that dU ~ const * grad(dp). If we are using convection-diffusion equation for pressure: ddt(psi,p) + div(phi_d, p) - div(rho/A*grad(p)), we are solving for mass flux, which satisfies continuity equation.

Another words:
- diffusion equation - correction to mass fluxes (to satisfy continuity) - phi+=pEqn.flux()
- convection-diffusion - mass fluxes, which continuity - phi=pEqn.flux()

I hope, this helps you. I'm really sorry for late post

please would you guide me in running a new application on OpenFOAM
 

i am new to OpenFOAM.
i was going through the "OpenFOAM Programming Tutorial" written by Tommaso Lucchini. i modified icoFoam to icoScalarTransportFoam and cavity to cavityScalarTransport in accordance to the directions of the notes.
now i cannot set up the case cavityScalarTransport and cannot run it. would you please guide me with the commands for running this new application.thanking you in advance.