
[Sponsors] 
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.cfdonline.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:
convectivediffusive 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:


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 (convectivediffusive) for flows with Ma >= 0.20.3 Second implementation (diffusive) for flows with Ma < 0.20.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 

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


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? 

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 

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 5, 2009, 14:06 
Buoyancy in OpenFOAM

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

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.
Last edited by nishant_hull; October 23, 2009 at 12:48. 

Thread Tools  
Display Modes  


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