buoyantBoussinesqPisoFoam: A detailed explanation
 

All,

I am beginning to use buoyantBoussinesqPisoFoam. Therefore, I thought it wise to really look at the inner-workings of the code both to understand what it is doing and to get used to OOP in C++. I had to do a lot of digging around, so I wrote up detailed notes of how buoyantBoussinesqPisoFoam works and posted it at http://docs.google.com/fileview?id=0...OTdhNDM4&hl=en. Please note, I was not able to fully understand everything, so if you find any mistakes in my notes, please let me know. I hope this may be helpful.

Matt

Dear mchurchf,
Your explanation of buoyantBoussinesqPisoFoam is indeed detailed. So the RANS and LES equation are nearly the same ,except in RANS using turbulent viscosity nu, LES use sub-grid viscosity nuSgs. Am I right ? This is a little like pisoFoam. But in OpenFOAM-1.6, buoyantBoussinesqPisoFoam can only be used for RANS model, can't be used for LES. Because I am beginner, I try to modify buoyantBoussinesqPisoFoam in the refference of pisoFoam, But I am affraid mistake will happen. Could you send your solver here ? That would be best. Thank you very much.

And if I want to add another passive scalar transport euqtion, for example concentration transport equation, it would be like this:

dC/dt +d/dxj*(C*uj)=(nut/Sct+nu0/Sc)*d/dxk(dC/dxk).

Good job :-)

P.S. Could you add the link to the wiki so the link does not go lost in the forum? http://openfoamwiki.net/index.php/Main_OFSolvers

Very good work... Thanks a lot for documenting.

best,

Dear alberto,

In Buoyant problem, I want to use nonuniform value Temperature in Wall. In my problem, my wall temperature is know, but not uniform. my temperature is changing with the hight of the wall. But when I modied the hotRoom/setHotRoom, and generated my own utility, It seems that boundary function can't be used for wall, only can be used for inlet condition. How I can do this ?
Thank you very much.

Does the wall temperature vary with the wall height according to a known law? If so, you can use GroovyBC

http://openfoamwiki.net/index.php/Contrib_groovyBC

to specify it.

Best,

Yes, alberto. This is utility is very good. It is better than my own modified utility.my own one can add inflow T profile, but can't add wall T profile, I don't khow why, Mybe my mesh was transformed from fluent, a littile not compatible with OpenFOAM. But groovyBC worked well.

Thank you very much.

buoyantBoussinesqPisoFoam Wiki and LES-capability
 

Alberto: I've posted this description of buoyantBoussinesqPisoFoam to http://openfoamwiki.net/index.php/Bu...sinesqPisoFoam.

Panda60: yes, you are correct in that regardless if it is RANS or LES, the equations look the same. In LES you supply nu_SGS and in RANS you supply nu_t. You are also correct that, out-of-the-box, buoyantBoussinesqPisoFoam only works with a RANS model. If you go to the above wiki posting, section 4, I explain how to modifiy it to be LES-capable just like pisoFoam.

Matt

Thanks,mchurchf.
Your explanation is very good. And in the guidence of you, I have modified the solve. But I'm a little confused.buoyant solver seems to converge slowly. and for p boundary condition, It seems very important. My case is a wind engineering problem. in outlet condition , it seems that using " fixedValue" for p is not good, and can't converge. when I change to buoyantPressure, then it slowly converge. I don't khow the what's the meaning of buoyantPressure, if it is only can be used in wall boundary, or also can be used for outlet condition ?

Because it seems not so many people researching on buoyant effect, could you give us some suggestion on buoyant solver ? what we must pay attation to.

Thank you very much.

Jiang,

When you say the solver converges slowly, what exactly do you mean? Do you mean that you are using it with an LES model, and it takes a lot of time for it to come to some sort of steady turbulent state? Or are you using it in RANS mode and trying to come to steady state? If you are using it for RANS, it might be better to use buoyantBoussinesqSimpleFoam which is meant for steady state problems.

My understanding of the bouyantPressure boundary condition is that it is used on walls that are not parallel to the hydrostatic pressure gradient. It appears to then apply dp/dz on the wall that corresponds to the hydrostatic pressure gradient.

What sort of problem are you solving? That may let me better help you.

Matt

Dear mchurchf,
I am very thankfulness for your kind hearted. My problem like this:
I have an rectangular simulation domain, an single building was located in the domain. Wind is coming from left to right, the ground temperature is higher than the air. My solver is buoyantBoussinesqSimpleFoam.

if I use {fixed value; value 0} for p in outlet condition, the flow didn't converge at all. in the outlet position, will have reverse flow in upper domain, and have strong velocity in lower domain.

if I use buoyantPressure for outlet condition, the flow seems converge ,but in the corner of domain, will have strong velocity ,the flow field is not right.

I searched the whole forum, somebody sai could use totalPressure or fixed meanPressure, I don't unstand .

So I think maybe the pressure boundary condition I ues is not right.
Could you tell me what pressure boundary condition you use ? especially in outlet ?
Thanks.

This is my pressure boundary condition:
BUILDING
{
type zeroGradient;
}
INLET
{
type zeroGradient;
}
OUTLET
{
type fixedValue;
value uniform 0;
}
SYM_LEFT
{
type symmetryPlane;
}
SYM_RIGHT
{
type symmetryPlane;
}
SYM_TOP
{
type symmetryPlane;
}
GROUND
{
type zeroGradient;
}
GAS_INLET
{
type zeroGradient;
}

This picture is my U_velocity contour.
I have been struggling this problem for several days. I need your help. Thanks.

Jiang,

At time zero, how do you initialize your pressure field? Do you have a constant pressure throughout the domain? If so, that could be a problem.

Pressure should decrease with height. Note that in buoyantBoussinesqSimpleFoam, I think the pressure variable is really pressure divided by a reference density, and we'll call it p~, so p~ = p/rho_ref. (Please check the equations to make sure.) Therefore, if you specify that some reference p~ at the ground, then p~(z) = p~(z=0) - g*z. In my simulations, I write my own code to create an initial p~ field that follows these hydrostatic conditions.

So then you have some choices for the outflow conditions. You need something to keep the flow moving, so you can either specify a constant value of velocity, U, at the top of the domain, which should drive the flow. Or you can specify a pressure gradient by making the outlet pressure lower than the inlet pressure. If you fix the inlet pressure at the hydrostatic conditions defined by p~(z) = p~(z=0)-g*z, then the outlet conditions could be p~(z) = p~(z=0) - g*z - (dp/dx)*L. The (dp/dx)*L is the pressure drop across the length, L, of your domain due to a pressure gradient dp/dx. You'll have to figure out which dp/dx is necessary to drive the flow at the velocity you want.

I think you are getting flow reversal and acceleration at the outflow because the pressure might be constant across the outflow plane. Can you post a contour plot of pressure, just like the plot of velocity that you show?

Yes, I just give a const pressure 0 as initial field.
and I use 290K as refference temperature. in Left ,Right and Top of the domain I use symmetry boundary condition.
This is my pressure field.

Because I can't modify the code, so I can't give hydrostatic initial conditions for pressure. I feel that , in buoyancy driven flow, the outlet condition is very important. In the guidance of you, today I try to use fixed value pressure in outlet like this: p=5-9.8*z. This fixed outlet value is small than internal field, an vary with hight, the flow is a little better, but also not correct.
near ground in outlet, the velocity is also high. And the internal pressure seems to change with the outlet pressure value. My temperature field is completely wrong. In my case , ground temperature is 45 degree, the air is 16 degree.

I try to use constant pressure gradient in outlet, but it doesn't work, have continues problem.

Jiang,

I think it would really help you if you started from initial conditions in which the pressure throughout the field was set to hydrostatic. I think all of your outlet problems are coming from the fact that the flow has a hard time adjusting from a constant 0 pressure to hydrostatic.

You'll have to write your own little code to write out the correct internal field. Once you have that initial condition for p, then you can try out different boundary conditions that you have already tried and see if they work.

If I were you, I would also extend my domain upwind and downwind of the building by a factor of probably 3 times what you already have. This will allow the effects of the building to be minimized at the boundaries.

You also said that your temperatures were 16 and 45 degrees. Do you mean degrees C. Note that OpenFOAM defaults to using Kelvin.

Matt

Dear mchurchf ,
I have initialize my pressure field using funkySetFields utilities.
Because my tempeture near ground is about 25 degree, so I estimate the pressure near ground is p=rho*R*T=1.22*287*(273+25)=104342 Pa,
dp=104342-101325=3017 Pa,

so I initilize my pressure field using p=3017-rho*g*z;
in outlet I also use fixed value 0 condition for p.
Now my simulation have going on for 6 hours, 2000 steps, result seems not good. If using SimpleFoam, 2000 steps can get good result . I will let my simulation going on to tomorrow morning ,and then see if can get good result.
If I should initialize my tempeture field using approximate value like pressure?

Hi and thanks a lot :D

I have one note about section 3.1. The pseudo-staggered arrangement is not applied because of the central differencing scheme. It effectively represent an improved version of the Rhie-Chow method (you can take a look at Hannu Karema thesis on multiphase flow to have the details, it is available online at VTT, Finland).
This strategy is very effective when you have rapidly changing source terms, for example.

Best,
Alberto

Dear mchurchf,

I already initialize my T and p field with good value, because I use RANS model, and in both side of the domain I use symmetry boundary condition, but after 5000 steps, the flow field is not symmetry, my flow field always have high asymmetry. Why this would happen ?
Thank you very much.

hi all

Thanks Mr mchurchf for the great work

i m new in openfoam i tried to use boussinesqBuoyantFoam with TransportProperties,
the solver is wrote by hjasak
first when i tried to compile this solver i got some dimension problems, so i made some changes ,the prog work but it dont converge (program aborted in less than 1s)

this is what i got

look at attached files

-> reponse.txt
-> boussinesq_1 : the solver
-> cdc_buoyantSimpleFoam ; my case solved with buoyantSimpleFoam
-> cdc : case i trie to solve

the changes i made : i multiply some euations by rho0 to adjust dimensions


best regards and thank you for your help



Attachment 2151

Attachment 2152

Attachment 2153

Attachment 2154

Dear all, I have strugling for using buoyantBoussinesqSimpleFoam for nearly half a month ,but don't have result. I am now doubt that if buoyantBoussinesqSimpleFoam is correct.

Now I describe my model in detail.
I have a rectangular domain, and a building was located in the domain.
wind come from inlet to outlet. the ground temperature is 318K, the air temperature is about 290K. the side and top domian I use symmetry boundary.

It seems that pressure doundary condition is very important when using buoyant solver. Then I try nearly all kinds of boundary condition for p. But all the result is not right. I am nearly crazy.

This is my model :

the following is my result:

1.The first result I get with the following pressure condition:
BUILDING
{
type zeroGradient;
}
INLET
{
type zeroGradient;
}
OUTLET
{
type fixedValue;
value uniform 0;
}
SYM_LEFT
{
type symmetryPlane;
}
SYM_RIGHT
{
type symmetryPlane;
}
SYM_TOP
{
type symmetryPlane;
}
GROUND
{
type zeroGradient;
}
GAS_INLET
{
type zeroGradient;
}

the iteration have going for more than 10000 steps,
I think ,for kEpsilon model it doesn't need so many iterations.

2.The second result I get with the buoyantPressure and totalPressure condition:
boundaryField
{
INLET
{
type buoyantPressure;
rho rhok;
value uniform 0;
}
OUTLET
{
type totalPressure;
p0 uniform 0;
U U;
phi phi;
rho rhok;
psi none;
gamma 1;
value uniform 0;
}
SYM_LEFT
{
type symmetryPlane;
}
SYM_RIGHT
{
type symmetryPlane;
}
SYM_TOP
{
type symmetryPlane;
}
GROUND
{
type buoyantPressure;
rho rhok;
value uniform 0;
}
BUILDING
{
type buoyantPressure;
rho rhok;
value uniform 0;
}
GAS_INLET
{
type buoyantPressure;
rho rhok;
value uniform 0;
}

3.The third result I get with the buoyantPressure for outlet:
boundaryField
{
INLET
{
type buoyantPressure;
rho rhok;
value uniform 0;
}
OUTLET
{
type buoyantPressure;
rho rhok;
value uniform 0;
}
SYM_LEFT
{
type symmetryPlane;
}
SYM_RIGHT
{
type symmetryPlane;
}
SYM_TOP
{
type symmetryPlane;
}
GROUND
{
type buoyantPressure;
rho rhok;
value uniform 0;
}
BUILDING
{
type buoyantPressure;
rho rhok;
value uniform 0;
}
GAS_INLET
{
type buoyantPressure;
rho rhok;
value uniform 0;
}
}

you could see, all the result is not right.
I was blamed by professor, why so easy problem you use so long time and can't get result.
Yes, It indeed an normal and easy problem when using Fluent. I don't know why it is so difficult to solve when using OpenFOAM .

Because I am already struggling this for a long time. Could anyone give me some help ?
what I should pay attation to when using buoyantBoussinesqFoam ?

what the exact meaning of buoyantPressure ? If it can be used for outlet ?
Thank you very much.

My inlet condition for U, k, epsilon and T was get from experiment, and I check ,no problem. I use upwind scheme for all the advection term.

and I think alphat file is ued for compressible flow. So for boussinesq flow it is not needed. Am I right ?

I think if some experts like Thomas Baumann , can come here give some suggestions , that would be excited.

Dear Jiang,

did you also try a unsteady calculation (e.g. using Piso)? From my experience:
If buoyancy effects really become dominant in the flow ,the flow pattern becomes inherently transient. I think it is the same for bouyantBoussinesq and the bouyant solvers using the state equations. Maybe you can try to your cases with these solvers as well (using ideal gas instead of Boussinesq approach).

Maybe this helps.
Stawrogin

Jiang,

I think your third pressure file should work. Note that if you are using buoyantBoussinesqPisoFoam you need to set a reference temperature in transportProperties. It's easy enough to make this quantity into a field quantity by editing the solver and recompiling. The difference between the temperature and the reference temperature can send a fluid blobs racing toward the top of the model and make the whole thing crash. In my experience there is too much variation in temperature over the scale of the PBL to just use a scalar quantity. This is why I typically use a Tref field or buoyantPisoFoam (even though it's for compressible flows).

I am interested in your temperature and U BC files. is your model heating up or cooling down? At the moment I'd just switch the model to laminar and focus on the basics.

Good Luck

Scott

Scott,

Can you clarify what you mean by using a Tref field? Do you mean that you set Tref not at just one point, but everywhere? In that way you can set the Tref field to be stratified?

If I understand correctly, if you set Tref to be stratified and initialize T to have the same stratification, then rho_k = 1 everywhere initially, correct?

Thanks,

Matt

That's exactly right, rhok would be zero initially but over time it won't be as blobs are displaced vertically (in this case by a building but also by terrain or turbulence) and advected to areas with a different Tref.

In the real atmosphere you get lots of buoyancy waves on a mesoscale due to stratification. I'm not convinced these things are dominate at the microscale (~100m) but they are present. My feeling is that when you turn on the turbulence a lot of these small scale wave features go away becasue TKE is larger in magnitude but I could be wrong.

Note, while I'm tossing the idea out here I haven't had a lot of luck with this modified solver. More luck than assuming a scalar Tref but I typically use buoyantPisoFoam for the atmospheric PBL. When I first started doing this I had concerns about using a compressible flow solver but from everything I've read I see no reason it won't work. I'd be interested in the groups thoughts regarding this subject (but lets not hijack this thread).

Scott

Dear stawrogin and scotth2o,
Thank you very much.
my inlet condition was got from experiment data, include U, T, k, epsilon.
The ground temperature is about 318K. The air temperature near ground is 295K. the air temperature in top is 282K. All these are from experiment. So in my model the air is a little heat up.
So the temperature differences is not large, dt=295-282=13K. I think this flow is dominanted by inlet velocity, not by buoyancy. The flow field only has buoyancy effect. so BuoyantBoussinesq is better than Buoyant buoyant solver.I have used field vararble TRef, but result not change so much.

This is my velocity boundary field. I use potentialFoam to initialize my velocity, and use experiment profile as inlet condition.

U
OUTLET
{
type zeroGradient;
}
SYM_LEFT
{
type symmetryPlane;
}
SYM_RIGHT
{
type symmetryPlane;
}
SYM_TOP
{
type symmetryPlane;
}
GROUND
{
type fixedValue;
value uniform (0 0 0);
}
BUILDING
{
type fixedValue;
value uniform (0 0 0);
}
GAS_INLET
{
type fixedValue;
value uniform (0 0 0.467);
}

This is my temperature file:

T
OUTLET
{
type zeroGradient;
}
SYM_LEFT
{
type symmetryPlane;
}
SYM_RIGHT
{
type symmetryPlane;
}
SYM_TOP
{
type symmetryPlane;
}
GROUND
{
type fixedValue;
value uniform 318.3;
}
BUILDING
{
type groovyBC;
valueExpression "509.51*pow(pos().z,4)-9.3428*pow(pos().z,3)+128.44*pow(pos().z,2)-51.296*pos().z +44.264+273.15";

}
GAS_INLET
{
type fixedValue;
value uniform 303.4;
}

Now I have been convinced that buoyantPressure is better for p outlet condition, because this result all the things are better, except that in the corner of the domain will have minus velocity after several steps, and then this minus velocity begin an instability to the whole field.Why this would happen ?

This is my picture for 200 steps. I think if there is not minus velocity in the domain corner, this result would be very good.

For the ground, gas_inlet and outlet try "type inletOutlet" in the temperature file. The model will likely heat up or cool down for a few time steps before reaching a "more" steady state

While I can't remember the energy equation for buoyantBoussinesqPisoFoam in buoyantPisoFoam there is a div(alphat,grad(h) ) term which leads to fluxes of energy out of the top and the bottom. Setting T to inletOutlet seems to let these fluxes flux.

Typically I would set the top temp to inletOutlet too but I have never used the "symmetryPlane" boundary condition. How does this work at the top of the model? Is the temperature on the right side the mirror image of the T on the left?

Good Luck
Scott

Dear scotth2o,

my relaxationFactors like this, am I right ?

relaxationFactors
{
rho 0.9;
p 0.3;
U 0.7;
T 0.7;
k 0.7;
epsilon 0.7;
R 0.7;
}

And I am little confused, why equation like this ?

fvScalarMatrix TEqn
(
fvm::div(phi, T)
- fvm::Sp(fvc::div(phi), T)
- fvm::laplacian(kappaEff, T)
);


tmp<fvVectorMatrix> UEqn
(
fvm::div(phi, U)
- fvm::Sp(fvc::div(phi), U)
+ turbulence->divDevReff(U)
);

why here have the term fvm::Sp(fvc::div(phi), T) ,and fvm::Sp(fvc::div(phi), U) ?
what is the meaning of this term ?

Because buoyantBoussinesaPisoFoam doesn't have this term. Why buoyantBoussinesaSimpleFoam have this ?

I've never changed the relaxationFactors from the default values. They seem to work fine as far as I know.

Scott

Dear all, who can tell me the meaning of alphat ?
If it is only be used for compressible flow ?
If I use buoyantBoussinesSimpleFoam , if I need it ?

Because I saw in buoyantPisoFoam, it write like this:

floor
{
type compressible::alphatWallFunction;
value uniform 0;
}

but in buoyantBoussinesqPisoFoam ,it write like this:

floor
{
type alphatWallFunction;
value uniform 0;
}.

So I think it will be used in BoussinesqFoam, but I can't find which places have used it in the solver , who can tell me ?

Because If using kEpsilon model , we should need wall function for temperature like velocity, we use nutWallFunction can activate the velocity wall function, So I think we should need a variable to activate temperature wall function ?

If ahphat can do this work ?
Thank you very much.

good job! Matthew.

eqns (15)(16)(18) should be d/dt(keff*dT/dx), not keff*d2T/dx2.
cuz keff(x) is not constant.

Poplar,

You are absolutely correct. I've updated the Wiki page on buoyantBoussinesqPisoFoam to reflect your corrections.

Thanks,

Matt

Boundary conditions for pressure
 

Hi,
I have a general question. It is of common knowledge that there is no need for boundary conditions for pressure when solving incompressible flow in enclosures. Then, what are the supposed boundary conditions for p and p_rgh in this case. I tried to put just constant zero value and obtained some solution that looks good, but is it physical?
Thanks, Yuri

Yuri,

If the flow is totally enclosed, you'll need to specify the pressure gradient normal to the surface. For p, static pressure, I recommend using the "buoyantPressure" boundary condition, which will set the static pressure gradient equal to (rho_k*g) dotted with the boundary surface normal, where rho_k is the Boussinesq density. For p_rgh, which is static pressure minus the hydrostatic part, also use "buoyantPressure". The "bouyantPressure" boundary condition looks at whether it is being used for p or p_rgh and applies the proper gradient. For p_rgh, it sets the gradient to -d(rho_k)/dn * g * height, where d/dn is the surface normal gradient.

See src\finiteVolume\fields\fvPatchFields\derived\buoy antPressure for details.

Matt

I started with open Foam to use it for geological problems. I'm trying to get into the buoyantBoussinesqPimpleFoam solver right now. Everything is working fine but I have one question:

I like to apply a linear temperature distribution as boundary conditions for the side walls. Is it possible to apply such in a simple way in openFOAM? In the tutorial there is setHotRoom used in the ./Allrun to define the boundaries. I didn't use that right know because I wasn't able to modify something in the setHotRoom. The only setHotRoom.C are in the tutorial for buoyantPimpleFoam but if I modify that (and recompile) I get no changes if I run setHotRoom with the Boussinesq solver. That sounds strange for me... Does someone have an idea about that??

Maike.

Hi Matthew

How did u deduce eqn (20) from eqn (19)? And what about the approximation order? Is there any references? Thanks.

// Kai