CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

Different Implementation in solving pEqn

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

Reply
 
LinkBack Thread Tools Display Modes
Old   September 3, 2009, 13:28
Default Different Implementation in solving pEqn
  #1
Senior Member
 
Nishant
Join Date: Mar 2009
Location: Glasgow, UK
Posts: 165
Rep Power: 8
nishant_hull is on a distinguished road
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
nishant_hull is offline   Reply With Quote

Old   September 7, 2009, 14:33
Default
  #2
Senior Member
 
mkraposhin's Avatar
 
Matvej Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 172
Rep Power: 8
mkraposhin is on a distinguished road
Quote:
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))

Quote:
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
mkraposhin is offline   Reply With Quote

Old   September 9, 2009, 08:45
Default
  #3
Senior Member
 
Nishant
Join Date: Mar 2009
Location: Glasgow, UK
Posts: 165
Rep Power: 8
nishant_hull is on a distinguished road
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
nishant_hull is offline   Reply With Quote

Old   September 9, 2009, 14:36
Default
  #4
Senior Member
 
mkraposhin's Avatar
 
Matvej Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 172
Rep Power: 8
mkraposhin is on a distinguished road
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
mkraposhin is offline   Reply With Quote

Old   September 9, 2009, 18:18
Default
  #5
Senior Member
 
Nishant
Join Date: Mar 2009
Location: Glasgow, UK
Posts: 165
Rep Power: 8
nishant_hull is on a distinguished road
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
Attached Files
File Type: doc sonicFoam.doc (28.5 KB, 56 views)
nishant_hull is offline   Reply With Quote

Old   September 10, 2009, 12:00
Default
  #6
Senior Member
 
mkraposhin's Avatar
 
Matvej Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 172
Rep Power: 8
mkraposhin is on a distinguished road
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 is offline   Reply With Quote

Old   September 10, 2009, 12:07
Default
  #7
Senior Member
 
mkraposhin's Avatar
 
Matvej Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 172
Rep Power: 8
mkraposhin is on a distinguished road
In attached file you will find corrected lines, marked with red color
Attached Files
File Type: doc sonicFoam.doc (12.5 KB, 58 views)
mkraposhin is offline   Reply With Quote

Old   September 10, 2009, 15:54
Default
  #8
Senior Member
 
Nishant
Join Date: Mar 2009
Location: Glasgow, UK
Posts: 165
Rep Power: 8
nishant_hull is on a distinguished road
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

Last edited by nishant_hull; September 10, 2009 at 16:26.
nishant_hull is offline   Reply With Quote

Old   September 10, 2009, 16:24
Default
  #9
Senior Member
 
mkraposhin's Avatar
 
Matvej Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 172
Rep Power: 8
mkraposhin is on a distinguished road
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?
mkraposhin is offline   Reply With Quote

Old   September 10, 2009, 16:50
Default
  #10
Senior Member
 
Nishant
Join Date: Mar 2009
Location: Glasgow, UK
Posts: 165
Rep Power: 8
nishant_hull is on a distinguished road
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
nishant_hull is offline   Reply With Quote

Old   September 11, 2009, 02:22
Default
  #11
Senior Member
 
mkraposhin's Avatar
 
Matvej Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 172
Rep Power: 8
mkraposhin is on a distinguished road
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
mkraposhin is offline   Reply With Quote

Old   September 11, 2009, 15:22
Default
  #12
Senior Member
 
Nishant
Join Date: Mar 2009
Location: Glasgow, UK
Posts: 165
Rep Power: 8
nishant_hull is on a distinguished road
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
nishant_hull is offline   Reply With Quote

Old   September 12, 2009, 05:31
Default
  #13
Senior Member
 
mkraposhin's Avatar
 
Matvej Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 172
Rep Power: 8
mkraposhin is on a distinguished road
Yes, it is alright.

This due to difference between OF1.4 and OF1.5
mkraposhin is offline   Reply With Quote

Old   September 15, 2009, 11:03
Default
  #14
Senior Member
 
Nishant
Join Date: Mar 2009
Location: Glasgow, UK
Posts: 165
Rep Power: 8
nishant_hull is on a distinguished road
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
nishant_hull is offline   Reply With Quote

Old   September 23, 2009, 16:15
Default
  #15
Senior Member
 
mkraposhin's Avatar
 
Matvej Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 172
Rep Power: 8
mkraposhin is on a distinguished road
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
mkraposhin is offline   Reply With Quote

Old   September 24, 2009, 14:17
Default
  #16
Senior Member
 
Nishant
Join Date: Mar 2009
Location: Glasgow, UK
Posts: 165
Rep Power: 8
nishant_hull is on a distinguished road
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 is offline   Reply With Quote

Old   September 25, 2009, 15:55
Default
  #17
Senior Member
 
Nishant
Join Date: Mar 2009
Location: Glasgow, UK
Posts: 165
Rep Power: 8
nishant_hull is on a distinguished road
see the attached coodle solver of 1.5.
Guess I can upload bigger files now i was having problems before.

regards,

Nishant
Attached Files
File Type: c coodles.C (2.5 KB, 14 views)
File Type: h createFields.H (1.1 KB, 7 views)
File Type: h hEqn.H (172 Bytes, 10 views)
File Type: h pEqn.H (1.3 KB, 22 views)
File Type: h UEqn.H (202 Bytes, 9 views)
nishant_hull is offline   Reply With Quote

Old   September 29, 2009, 13:43
Default
  #18
Senior Member
 
Nishant
Join Date: Mar 2009
Location: Glasgow, UK
Posts: 165
Rep Power: 8
nishant_hull is on a distinguished road
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

Last edited by nishant_hull; October 1, 2009 at 11:58.
nishant_hull is offline   Reply With Quote

Old   October 5, 2009, 14:06
Default Buoyancy in OpenFOAM
  #19
Senior Member
 
mkraposhin's Avatar
 
Matvej Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 172
Rep Power: 8
mkraposhin is on a distinguished road
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)
mkraposhin is offline   Reply With Quote

Old   October 12, 2009, 14:07
Default
  #20
Senior Member
 
Nishant
Join Date: Mar 2009
Location: Glasgow, UK
Posts: 165
Rep Power: 8
nishant_hull is on a distinguished road
i think i did some stupid mistake before.
Attached Images
File Type: jpg sonicref1p5.jpg (53.1 KB, 35 views)
File Type: jpg sonicref1p6.jpg (64.1 KB, 26 views)

Last edited by nishant_hull; October 23, 2009 at 12:48.
nishant_hull 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
Moving mesh Niklas Wikstrom (Wikstrom) OpenFOAM Running, Solving & CFD 122 June 15, 2014 06:20
Differences between serial and parallel runs carsten OpenFOAM Bugs 11 September 12, 2008 11:16
Unknown error sivakumar OpenFOAM Pre-Processing 9 September 9, 2008 12:53
IcoFoam parallel woes msrinath80 OpenFOAM Running, Solving & CFD 9 July 22, 2007 02:58
Could anybody help me see this error and give help liugx212 OpenFOAM Running, Solving & CFD 3 January 4, 2006 19:07


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