[Anne Lincke]

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= (1-beta(T-T_0))*g

In pEqn.H the main code is settled.

First the velocity is solved without taking into account pressure or density.

U = rAU*UEqn().H();

The flux is the interpolation of this velocity.

phi = fvc::interpolate(U) & mesh.Sf();

Then a buoyancy flux is defined

surfaceScalarField buoyancyPhi(rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());

Here my first question: Why is the normal gradient computed? In the equation there is only a term which includes the density, not the gradient of the density.

The flux is corrected with the buoyancy flux

 phi -= buoyancyPhi;

After the pressure correction with the Laplace equation

fvm::laplacian(rAUf, p_rgh) == fvc::div(phi)

the flux is again corrected by taking into account the pressure minus rgh.

phi -= p_rghEqn.flux();

In the last step the velocity is reconstructed from the flux

U -= rAU*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rAUf);

and p_rgh is updated

 p_rgh = p - rhok*gh;

I do not understand the components of the flux. If they are added we obtain

phi = fvc::interpolate(U) & mesh.Sf() - (rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf()) - p_rghEqn.flux();

So the term including the normal gradient of the density is too much.
Regarding the link
http://openfoamwiki.net/index.php/Bu...sinesqPisoFoam

in the solver BuoyantBoussinesqPisoFoam the flux looks like this

phi =fvc::interpolate(U) & mesh.Sf()) + rUAf*fvc::interpolate(rhok)*(g & mesh.Sf())-pEqn.flux();

which makes much more sense to me.
I can not find any explanation in the forum. Who can help me to understand the code?

[Anne Lincke]

No one for an answer???

[romant]

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.

[Anne Lincke]

Thank you for answering!
Now I understand the implementation. If someone is interested I could post that, too.

[Bernhard]

If it is not on this forum yet, then it can be very interesting for future reference if you post your findings here!

[Anne Lincke]

In the momentum equation we have in z-direction (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

 - fvc :: snGrad (p_rgh)

Computing the derivative of this term with product rule and assuming that rho may change with respect to the z-direction we obtain

- 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

- ghf *fvc::snGrad(rhok)

So the derivative of the density is a correction term which is nonzero for changes in density.

[makaveli_lcf]

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

[Anne Lincke]

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 z-direction. 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 right-hand side of the UEqn

 - fvc::snGrad(p_rgh) - ghf*fvc::snGrad(rhok)

2. I have a lecture script where this method is described, but it is not open for everyone. I also found a notice for this method in Peric and Ferziger's "Computational Methods for Fluid Dynamics" at the very beginning of the book (page 11 in the version I have). Search for "hydrostatic pressure" and you will find a short explanation there.

[makaveli_lcf]

Thank you Anne for the detailed explanation! I will search for the algorithm description in Peric's book.

Got it, thanx)

[siamak 60]

Hi

As anyone knows in buossinesq approximation for density is written as :
rho=rhoRef*(1-beta*(T-Tref))
and in oF rhok=1-beta*(T-Tref)


and in Chalmer university file:
Cite as: Venkatesh, B V.: Tutorial of convective heat transfer in a vertical slot.
In Proceedings of CFD with OpenSource Software, 2016, Edited by Nilsson. H.,
http://www.tfd.chalmers.se/~hani/kurser/OS_CFD_2016
CFD with OpenSource software
A course at Chalmers University of Technology Taught by Håkan Nilsson


we have:


The governing equation for temperature is solved, given as
∇.(ρuT ) − ∇.αef f ∇T = Sradiation + ST (2.4)
Where, αef f = νt /prt + ν/pr , Sradiation and ST are source terms due to radiation and user defined source term respectively. The density is updated as ρ = 1 − β(T − Tref )



and phi is defined as phi=rhok*U

now question is how rhoRef is deleted in above equstion while alpha have it:
α=k/(ρ*Cp)


please help in making it clear.



thanks in advance