
[Sponsors] 
July 10, 2011, 08:33 
interpolation for pressure

#1 
New Member
Harry
Join Date: Mar 2009
Posts: 27
Rep Power: 8 
Does Openfoam implement bodyforceweighted scheme for pressure interpolation?


July 11, 2011, 02:36 

#2 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26 
What solver are you referring to?


July 11, 2011, 02:44 

#3 
New Member
Harry
Join Date: Mar 2009
Posts: 27
Rep Power: 8 
I am new to OpenFoam. this type of interpolations is necessary in the Eulerian model based on the SIMPLE velocity and pressure coupling method in the solution of gassolid/liquid flows where large source terms present due to the strong interaction between phases. it is necessary associated with the intrinsic flaw in Rhie and Chow momentum interpolation practice. A big CFD issue but very limited publications devoted to the problem. It has puzzled me for a long time.


July 11, 2011, 03:09 

#4 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26 
Yes, bubbleFoam and twoPhaseEulerFoam implement an improved interpolation approach to treat gravity and drag. In OpenFOAM such a practise is implemented in some of the solvers (buoyant, multiphase, ...). You can see how this is managed in bubbleFoam here: http://openfoamwiki.net/index.php/BubbleFoam
Basically the force terms are treated as a source term in the momentum equation (well in bubbleFoam they are literally dropped, see interFoam for a complete example), after faceflux reconstruction, and then are moved to the flux used to construct the pressure equation. You might also want to take a look at how OpenFOAM does RhieChow: http://web.student.chalmers.se/~f98faka/rhiechow.pdf P.S. My area of research is multiphase flow too, and I absolutely agree there is not enough good literature on the topic, and what is available is usually hidden in some conference proceedings! 

July 11, 2011, 06:57 

#5  
New Member
Harry
Join Date: Mar 2009
Posts: 27
Rep Power: 8 
Quote:


July 11, 2011, 10:41 

#7 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26 
Sorry, the link on the wiki (where I took it) to the RhieChow document is broken. Laurence gave you the correct one.


July 12, 2011, 00:14 

#8 
New Member
Harry
Join Date: Mar 2009
Posts: 27
Rep Power: 8 
I have read through the relevant documents but cannot get the main idea about how OpenFoam implement the momentum interpolation under the condition with/without body force as source terms. I appreciate that some guy would like to outline this method or redirect me to some further publications.
Best regards, Shibo 

July 12, 2011, 01:34 

#9  
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26 
Quote:
(1) ddt(U) + div(UU) = div(tau)  grad(p)/rho + F This equation can be written in semidiscrete form as: (2) A*U = H  grad(p)/rho + F where the pressure gradient and the force term we want to include in the momentum interpolation is left explicitly out at this point. In other words, only the first three terms (time derivative, unsteady, divergence of the stress tensor) of Eq. 1 are used to define A and H. For example, in an incompressible code, this could read (it might be different, depending on the solver): Code:
fvm::ddt(U) + fvm::div(phi, U) + turbulence>divDevReff(U) From Eq. 2, interpolating on faces and dotting with the surface area vector S the pressure gradient and the force term, we have:  snGrad(p)*S/rho + F_f . S where F_f is F interpolated on faces. This in OF corresponds to: Code:
 fvc::snGrad(p)*mesh.magSf()/rhoa + fvc::interpolate(F) & mesh.Sf() Code:
solve ( UEqn == fvc::reconstruct (  fvc::snGrad(p)*mesh.magSf()/rhoa + fvc::interpolate(F) & mesh.Sf() ) ) U = H/A where A is updated with the predicted value of U. However, keep in mind that H does not directly include the effect of grad(p) and F. From Eq. 2 we can derive the flux to construct the pressure equation: phi = (H/A)_f . S  (1/A)_f * snGrad(p) S/rho + (1/A)_f*F_f . S and imposing div(phi) = 0, you obtain the pressure equation, which will include for the first time all the effects. The flux is then corrected based on the solution of the pressure equation, and the velocity correction is reconstructed from the flux. Only at this point U will "see" the effect of p and F. If you take a look at VOF solvers (i.e. interFoam), you will notice that the code does not solve for p, but for p_rgh = p + rho*gh. This is another way to treat the gravity to address the weakness of the RhieChow interpolation, but the idea is similar. I hope this helps, but please let us know if you have more questions Best,
__________________
Alberto Passalacqua GeekoCFD  A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as live DVD/USB, hard drive image and virtual image. OpenQBMM  An opensource implementation of quadraturebased moment methods Last edited by alberto; July 12, 2011 at 01:35. Reason: Added [code] tags 

July 12, 2011, 07:12 

#10 
New Member
Harry
Join Date: Mar 2009
Posts: 27
Rep Power: 8 
thank you so much, alberto! this is a quite interesting treatment which I am not quite sure that I have completely understood it.
Below is my understanding of the entire procedure, which seems to be achieved through two steps. 1) OpenFoam solves the momentum twice. Firstly, the momentum is solved with pressure gradient and body force. Secondly, the momentum is solved without pressure gradient and body force. 2) the velocity from the second solution of momentum is used to establish pressure correction equation, considering the contributions of pressure gradient and body force. If my understanding is correct, it seems to me that it is not necessary to interpolate of gradient and body force to the cell face and add them to the momentum. Did I missed something? I need to get confirmed before I can raise other questions. BTW: I am completely new to OpenFoam, and cannot thus far discuss this using Openfoam style language. 

July 12, 2011, 12:05 

#11 
Senior Member
Andrea Ferrari
Join Date: Dec 2010
Posts: 269
Rep Power: 7 
Hi all,
i'm writing my problem here, since it is similar. i have a question concerning the solution procedure in OF. I'm struggling with problem of spurious velocities at the interface using interFoam and i've read that basically these spurious currents could come from a sort of errors in the splitting of the equations, because there is a force unbalance between the pressure gradient and the surface force term in the 2 steps predictor and corrector. They have suggested to add another correction in the flux calculation (before solving the pressure), which basically has to take into account the difference between these two terms. here the post of my problem Small time step in interFoam (in the post can be also found the title of the paper) i would like to know what do you think about that and if it would be possibile to implement in OF? thanks andrea 

July 12, 2011, 14:50 

#12  
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26 
Quote:
The estimated value of U from this equation is used to update what OF calls H (see my previous post). Then, at this point the flux is defined consistently with the semidiscrete form of the momentum equation, in order to build the pressure equation. This means you need to compute U = H/A, and add the terms formally left out in the predictor. Quote:
Best, 

July 12, 2011, 21:14 

#13  
New Member
Harry
Join Date: Mar 2009
Posts: 27
Rep Power: 8 
Thanks, alberto! this helps a lot! I have almost got the picture of the method in OF, but still with one unclear point. I appreciate your further helps. My questions are given below:
Quote:


July 12, 2011, 21:56 

#14  
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26 
Quote:
solve ( UEqn == fvc::grad(p) ); This does not work so well in the case of strong density gradients or changes in porosity, for example. It is a known problem, and one of the solutions was proposed in Zhang S., Zhao, X., General formulations for RhieChow interpolation, ASME Heat Transfer/Fluids Engineering Summer Conference, HTFED04, Charlotte, USA, 2004. Such a solution suggests to average the body force term (and other strongly varying terms) on cell faces, and then reconstruct the source term from the faceaveraged value. This is very similar to what OpenFOAM does in the lines of codes you are referring to. Best, 

July 13, 2011, 01:10 

#15 
New Member
Harry
Join Date: Mar 2009
Posts: 27
Rep Power: 8 
Dear alberto,
In the Zhang's paper, the treatment of body force seems to be similar to that done by OF. It is achieved through threes steps: 1. First define body forces primitively at cell centers, FP 2. Average these to cell faces, Ff, using weight factors of one half. 3. Redefine F at cell centers in a manner consistent with the second order formula for nodal pressure gradient. I cannot understand the third step. this is also the reason that I cannot understand the treatment in OF Would you please detail the theory for the third step? In OF, this is done by fvc: reconstruct. 

July 14, 2011, 04:30 

#16 
Senior Member
Andrea Ferrari
Join Date: Dec 2010
Posts: 269
Rep Power: 7 
Hi,
if the momentum predictor is turned off, where the solution "can see" the effect of the old pressure? (I'm talking about interFoam) in the flux definition only the surface force and the gravity are explicitly added to U_f=(H/A)_f (where both H and A do not depend on the old pressure...correct me if i'm wrong) So if the momentum predictor is turned off, does the new pressure depend on the old pressure? thanks andrea 

July 14, 2011, 11:09 

#17  
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26 
Quote:
@harry: each face flux will provide you a contribution to the cell centered value, which can be reconstructed then based on the face values you defined. In other words you integrate the contribution of the face values to recover the cell value by means of integration. @Andrea_85: the velocity U sees the pressure immediately after the solution of the pressure equation, when you correct it based on the new pressure gradient (you do this by flux reconstruction in interFoam). Best, 

July 14, 2011, 12:16 

#18 
Senior Member
Andrea Ferrari
Join Date: Dec 2010
Posts: 269
Rep Power: 7 
Hi Alberto and thanks for reply.
What i cannot figure out is the case in which the momentum predictor is turned off (set off in the PISO subdictionary in fvSolution). Which is in this case the predicted velocity? In Ueqn.H OF defines A and H, but H still depends on the velocity (that is unknown in my understanding if the momentum preditor is not solved). The system is: A[U]=[SourceTerm], then you can split in diagonal and offdiagonal part A[U] + Aoffd[U] = [S] and then H is defined as H=[S]aoffd[U] If the momentum predictor is turned to off, H is something like: H=  sum_neigbours(a_n*U_n)+U^0/deltaT (U^0 is the velocity from the last time step) the velocities of the "neighboring" are unknown (or not?...are from the previous time step?) so i dont understand how can you calculate the velocity using something that depends on the velocity itself (U=H(U)/A)... Best Regards andrea 

September 6, 2011, 10:09 

#19 
Senior Member

Hi all!
Thanx for a very interesting discussion! Would you please clarify what does fvc::reconstruct operator actually do? Is it: fvc::reconstruct (f_flux) = Surf_integral (f_face_interpolated . Sf) = Sum_over_all_faces(f_face_interpolated . Sf) ????? Thank you in advance!
__________________
Best regards, Dr. Alexander VAKHRUSHEV Christian Doppler Laboratory for "Advanced Process Simulation of Solidification and Melting" Simulation and Modelling of Metallurgical Processes Department of Metallurgy University of Leoben FranzJosefStr. 18 A  8700 Leoben Österreich / Austria Tel.: +43 3842  402  3125 http://smmp.unileoben.ac.at 

September 6, 2011, 21:56 

#20 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26 
From ~/OpenFOAM/OpenFOAM2.0.x/src/finiteVolume/finiteVolume/fvc/fvcReconstruct.C, you find that the reconstruct method is defined as:
Code:
template<class Type> tmp < GeometricField < typename outerProduct<vector,Type>::type, fvPatchField, volMesh > > reconstruct ( const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf ) { typedef typename outerProduct<vector, Type>::type GradType; const fvMesh& mesh = ssf.mesh(); tmp<GeometricField<GradType, fvPatchField, volMesh> > treconField ( new GeometricField<GradType, fvPatchField, volMesh> ( IOobject ( "volIntegrate("+ssf.name()+')', ssf.instance(), mesh, IOobject::NO_READ, IOobject::NO_WRITE ), inv(surfaceSum(sqr(mesh.Sf())/mesh.magSf())) & surfaceSum((mesh.Sf()/mesh.magSf())*ssf), zeroGradientFvPatchField<GradType>::typeName ) ); treconField().correctBoundaryConditions(); return treconField; } 

Tags 
body force, force 
Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
urgent help needed (rhiechow interpolation problem)  Ardalan  Main CFD Forum  2  March 18, 2011 16:22 
Help on 2D interpolation in StarCCM+  madhuri  CDadapco  0  November 3, 2010 15:21 
Surface interpolation schemes and parallelization  jutta  OpenFOAM Running, Solving & CFD  0  February 25, 2010 15:32 
momentum interpolation for collocated grid  Hadian  Main CFD Forum  4  December 25, 2009 08:25 
spline interpolation  bajjal  Main CFD Forum  0  May 29, 2006 08:27 