Problems in understanding BuoyantBoussinesqSimpleFoam
1 Attachment(s)
Dear Foamers,
I still habe problems in understanding the source code of BuoyantBoussinesqSimpleFoam. In the attachment you see the momentum equation which has to be solved where GT= rho g= (1beta(TT_0))*g In pEqn.H the main code is settled. First the velocity is solved without taking into account pressure or density. HTML Code:
U = rAU*UEqn().H(); HTML Code:
phi = fvc::interpolate(U) & mesh.Sf(); HTML Code:
surfaceScalarField buoyancyPhi(rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf()); The flux is corrected with the buoyancy flux HTML Code:
phi = buoyancyPhi; HTML Code:
fvm::laplacian(rAUf, p_rgh) == fvc::div(phi) HTML Code:
phi = p_rghEqn.flux(); HTML Code:
U = rAU*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rAUf); HTML Code:
p_rgh = p  rhok*gh; HTML Code:
Regarding the link http://openfoamwiki.net/index.php/Bu...sinesqPisoFoam in the solver BuoyantBoussinesqPisoFoam the flux looks like this HTML Code:
phi =fvc::interpolate(U) & mesh.Sf()) + rUAf*fvc::interpolate(rhok)*(g & mesh.Sf())pEqn.flux(); I can not find any explanation in the forum. Who can help me to understand the code? 
No one for an answer???

I am not sure, but you might have to look up in a numerical methods book, or somewhere else how the SIMPLE algorithm works, maybe that can give you some hints.

Thank you for answering!
Now I understand the implementation. If someone is interested I could post that, too. 
If it is not on this forum yet, then it can be very interesting for future reference if you post your findings here!

In the momentum equation we have in zdirection (in direction of buoyancy) the terms
 dp / dz + rho*g In OpenFOAM, g is a vector (0 0 9,81) which ensures that the buoyancy is only valid for the right coordinate direction. In order to guarantee, that in the pressure correction of the simple algorithm also the buoyancy term rho*g*z is taken into account, the pressure and the buoyancy are melted together in one term p_rgh = p  rho*g*z Instead of the normal gradient of p, the normal gradient of p_rgh is on the RHS.  d/dn (p_rgh) which equals in OF to HTML Code:
 fvc :: snGrad (p_rgh)  d/dz [ p  rho * g *z] =  dp/dz + rho*g + g*z* d rho/dz So the third term g*z* d rho/dz is "too much" and is therefore substracted via HTML Code:
 ghf *fvc::snGrad(rhok) 
Hallo Anne!
Would you please clarify following points: 1. What do you mean by term "two much" (is it phisical meaning or regards difference in some mathematical formulation). 2. Did you use some reference to understand SIMPLE/PISO treatment of the buoyancy? Thank you in advance! Mfg, Alexander 
Hallo Alexander,
I will try to answer your questions: 1. With "too much" I mean the term g*z* d rho/dz which is the result of taking the derivative of  [ p  rho * g *z] in zdirection. It is an additional term due to the product rule. Basically we have the terms  dp / dz + rho*g So by taking the derivative of both, p and rho*g*z we obtain together with the correction term  g*z* d rho/ dz:  d/dz [ p  rho * g *z]  g*z* d rho/ dz=  dp /dz + rho*g + g*z* d rho/ dz  g*z* d rho/ dz =  dp/ dz + rho*g and this is the orignial term in the buoyancy driven momentum equation. In OpenFOAM code this equals to the righthand side of the UEqn HTML Code:
 fvc::snGrad(p_rgh)  ghf*fvc::snGrad(rhok) 
Thank you Anne for the detailed explanation! I will search for the algorithm description in Peric's book.
Got it, thanx) 
All times are GMT 4. The time now is 11:34. 