# Different Implementation in solving pEqn

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

 September 3, 2009, 13:28 Different Implementation in solving pEqn #1 Senior Member   Nishant Join Date: Mar 2009 Location: Glasgow, UK Posts: 165 Rep Power: 10 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

September 7, 2009, 14:33
#2
Senior Member

Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 330
Rep Power: 12
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:

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

 September 9, 2009, 08:45 #3 Senior Member   Nishant Join Date: Mar 2009 Location: Glasgow, UK Posts: 165 Rep Power: 10 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

 September 9, 2009, 14:36 #4 Senior Member     Matvey Kraposhin Join Date: Mar 2009 Location: Moscow, Russian Federation Posts: 330 Rep Power: 12 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

September 9, 2009, 18:18
#5
Senior Member

Nishant
Join Date: Mar 2009
Location: Glasgow, UK
Posts: 165
Rep Power: 10
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
 sonicFoam.doc (28.5 KB, 70 views)

 September 10, 2009, 12:00 #6 Senior Member     Matvey Kraposhin Join Date: Mar 2009 Location: Moscow, Russian Federation Posts: 330 Rep Power: 12 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(...)?

September 10, 2009, 12:07
#7
Senior Member

Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 330
Rep Power: 12
In attached file you will find corrected lines, marked with red color
Attached Files
 sonicFoam.doc (12.5 KB, 72 views)

 September 10, 2009, 15:54 #8 Senior Member   Nishant Join Date: Mar 2009 Location: Glasgow, UK Posts: 165 Rep Power: 10 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.

 September 10, 2009, 16:24 #9 Senior Member     Matvey Kraposhin Join Date: Mar 2009 Location: Moscow, Russian Federation Posts: 330 Rep Power: 12 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? Saturne_User likes this.

 September 10, 2009, 16:50 #10 Senior Member   Nishant Join Date: Mar 2009 Location: Glasgow, UK Posts: 165 Rep Power: 10 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

 September 11, 2009, 02:22 #11 Senior Member     Matvey Kraposhin Join Date: Mar 2009 Location: Moscow, Russian Federation Posts: 330 Rep Power: 12 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

 September 11, 2009, 15:22 #12 Senior Member   Nishant Join Date: Mar 2009 Location: Glasgow, UK Posts: 165 Rep Power: 10 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

 September 12, 2009, 05:31 #13 Senior Member     Matvey Kraposhin Join Date: Mar 2009 Location: Moscow, Russian Federation Posts: 330 Rep Power: 12 Yes, it is alright. This due to difference between OF1.4 and OF1.5

 September 15, 2009, 11:03 #14 Senior Member   Nishant Join Date: Mar 2009 Location: Glasgow, UK Posts: 165 Rep Power: 10 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

 September 23, 2009, 16:15 #15 Senior Member     Matvey Kraposhin Join Date: Mar 2009 Location: Moscow, Russian Federation Posts: 330 Rep Power: 12 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

 September 24, 2009, 14:17 #16 Senior Member   Nishant Join Date: Mar 2009 Location: Glasgow, UK Posts: 165 Rep Power: 10 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

September 25, 2009, 15:55
#17
Senior Member

Nishant
Join Date: Mar 2009
Location: Glasgow, UK
Posts: 165
Rep Power: 10
see the attached coodle solver of 1.5.
Guess I can upload bigger files now i was having problems before.

regards,

Nishant
Attached Files
 coodles.C (2.5 KB, 15 views) createFields.H (1.1 KB, 8 views) hEqn.H (172 Bytes, 11 views) pEqn.H (1.3 KB, 31 views) UEqn.H (202 Bytes, 10 views)

 September 29, 2009, 13:43 #18 Senior Member   Nishant Join Date: Mar 2009 Location: Glasgow, UK Posts: 165 Rep Power: 10 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.

October 12, 2009, 14:07
#20
Senior Member

Nishant
Join Date: Mar 2009
Location: Glasgow, UK
Posts: 165
Rep Power: 10
i think i did some stupid mistake before.
Attached Images
 sonicref1p5.jpg (53.1 KB, 49 views) sonicref1p6.jpg (64.1 KB, 34 views)

Last edited by nishant_hull; October 23, 2009 at 12:48.

 Thread Tools Display Modes Linear Mode

 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 OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post Niklas Wikstrom (Wikstrom) OpenFOAM Running, Solving & CFD 122 June 15, 2014 06:20 carsten OpenFOAM Bugs 11 September 12, 2008 11:16 sivakumar OpenFOAM Pre-Processing 9 September 9, 2008 12:53 msrinath80 OpenFOAM Running, Solving & CFD 9 July 22, 2007 02:58 liugx212 OpenFOAM Running, Solving & CFD 3 January 4, 2006 19:07

All times are GMT -4. The time now is 00:14.