I would like to account for the pressure gradient due to gravity in a static fluid. I know how to define gravity as an environmental field, and the kinematic viscosity "nu" includes the density of the fluid. Therefore, we should be able to calculate the pressure gradient in a static body of fluid, but I do not see this when using icoFoam or potentialFoam. I tried using a very high gravitational acceleration to magnify the effect, but still no gradient. Is there a solver which will do this, or am I doing something else wrong? Do I need to specify an absolute pressure at the "upper" boundary condition?

More generally, there is not much documentation on the solvers in the User or Programmer's manual. Searching for "icoFoam" in Doxygen turns up 0 results. If there is a place to read more about the solvers, please tell me! I hate to bother you people with simple RTFM questions. Thanks!

Most of the incompressible flow codes don't have gravity included although it's simply a matter of adding a g source term to the momentum equation (note the incompressible codes are usually divided by the constant density so p is actually p/rho). However adding gravity to a constant-density case won't do much apart from generating an underlying constant pressure gradient which is why gravity isn't included except in the interface capturing codes where the difference in the phase densities can drive the motion as in the dam-break case.

OpenFOAM also includes buoyant flow solvers buoyantFoam and buoyantSimpleFoam which include the graviational effect in a manner suitable for highly buoyant problems.

You are quite correct there isn't much documentation on the solvers in OpenFOAM and due lack of funds we cannot afford to employ anyone to continue working on it. Writing documetation is an extremely time consuming business and it seems people are not keen to work for nothing so we are looking for sponsorship for this and many other OpenFOAM development projects.

For the moment, the best way to work around the lack of documentation is to become familiar with the source code. Being C++ it is comparatively easy to understand at the top level, which is close to a symbolic syntax. All solver source can be found in OpenFOAM-1.0.x/applications/solvers and is sorted into categories, making it easier to find what you are looking for.

As Henry mentioned, an example of the gravity source can be found in interFoam under multiphase. The code which you (probably) want to modify is turbFoam under the incompressible category.

[coops]

Hi all,

I have a problem trying to implement hydrostatic equilibrium within my model. I am using an altered version of sonicFoam as the solver. The details are as follows:

I am trying to study some atmospheric dynamics. The model that I am working with extends from 100 km to 400 km vertically. Although the model is a transient 3-D model I will limit this to the vertical direction, the one causing my problems!

Initially I give the model a temperature profile that varies as an eighth order polynomial with height, with the lower boundary having a value of 180K and the top a value of 1100K. This function is then easily integratable and allows for an initial temperature profile to be created. Following this one is able to calculate rho from the ideal gas law (which seems to be the way the model already calculates the density).

On paper the fields that are supplied should satisfy hydrostatic equilibrium, however, when I run the model with these field there is always an acceleration in the vertical direction.

Does anyone have any idea what this might be due to? I am using a mesh that is about 300 km vertically with cells every 1500 m at the moment. Where are the values for T and P specified? Surely they are at the same position for rho. Is it possible that the gradient (dp/dz) is evaluated at a different position to -rho*g and hence they are not in equilibrium at each height?

Any suggestions are much appreciated,

Shaun

[coops]

Hi all,

Just to correct the previous message. The temperature profile is integratable which allows the initial pressure profile to be calculated. Using these two field the ideal gas law should be used to calculate rho. On paper, this formulation goes have dP/dz = -rho*g. However, when the field are calculated they are not consistent!

Any help would be great,

Shaun

[coops]

Hello All,

I am still working on this hydrostatic equilibrium problem and haven't been able to resolve it. Can anyone see that there is an obvious problem with this approach???

I am using a 3-D domain that is 1000 km in the x-direction, 300 in the z-direction and 600 in the y-direction.

I give my model an initial pressure profile that varies approximately exponentially (I use a polynomial fit) with height. The boundary conditions of this field are wallbouyantPressure on all boundaries.

I then give an initial temperature profile that is also approximately exponential. The top and bottom boundaries are given fixed values that are obtained by substituting in the required altitudes. The other boundaries are given zeroGradient boundary conditions.

The temperature and pressure fields are designed such that the ideal gas law should allow the density to be calculated. In the altered version of sonicFoam that I am using rho = psi * p where psi = 1 / (R * T). This is the ideal gas law, yes!?

However, when I plot the initial fields using paraFoam and use the line probe the rho field is not as I would expect. At the boundaries the rho function is not smooth, there is a large positive (greater than the mean) spike followed by a negative (+ve value but less that the mean) spike. I then look at the psi field and this is smooth, I look at the temperature field and it is smooth but the product of the two is not! Why is this????

To test out the equations I made a field called myrho and do the following:



const vectorField& centres = mesh.C();

forAll( centres, celli )
{

T[celli] = ....

p[celli] = .....


myrho[celli] = p[celli] / (Rvalue * T[celli]) ;

}
where Rvalue is the value of R!

When I look at myrho the function is smooth and equal to rho everywhere apart from where rho has the spikes!

Any hint as to the cause of this problem would be appreciated!

Regards,

Shaun

[david_h]

Have you tried extracting rho using the "sample" utility instead of paraFoam ?

[coops]

Hello David,

I have now extracted rho using the sample utility. However, I still have the same problem. Here are some plots of the problem:



This is an image of rho from 100 km to 400 km in altitude. As you can see there are some problems near the 100 km mark!




These are the two fields from which rho is meant to be calculated. Why is there a problem at th boundary?

Here is an image of rho and myrho (calculated from myrho = p[celli] / (R * T[celli]) ) from 100 to 115 km.



As you can see, there is something happening in the calculation of rho. I don't know what is the cause of this.

If I haven't explained myself properly let me know and I will try to explain it in more detail. Also, if someone wanted to look at my code and help me determine the errors of my ways I would be happy to post it, just let me know.

Thanks to all that read this. Please offer any assistance you can think of.

Thanks in advance,

Shaun

[david_h]

Earlier you mentioned that you based your solver on sonicFoam. Have you tried using buoyantFoam ?
sonicFoam in contrast to buoyantFoam does not update the continuity and energy equations inside the pressure correction loop.

Also, you might compare your implementation of the buoyancy force with implementations in existing applications. Each of the following solvers uses a different method to include buoyancy:
a) reactingFoam/dieselEngineFoam
b) buoyantFoam
c) interFoam

[coops]

Hi,

I have used the buoyantFoam solver but there is a strange result if I put an initial pressure gradient. For example, if I make the initial pressure 2e-5 * centres[celli][2] (the height) then after the model has run (adding a rho * g term in which should offset the pressure gradient exactly) there are vertical winds produced. They are alternatively stratified (positive then negative then positive etc.). Overall, it looks like the equlibrium is almost obtained but is not satisfied locally.

Can someone correct me if I am wrong but the values of the gradients are computed at the cell faces from the values at the cells centres. So, for example, if I know the pressure at two cells the gradient that is computed is the gradient at the cell face.

I think that this may be a possible cause of some of the issues I am having with the hydrostatic equilibrium problem.

Thanks,

Shaun

[coops]

Hello,

If the above post is correct then it would be useful for me to be able to generate a density value at the cell face. Is there a way of doing this already? I looked at faceInterpolation but I think this might be operating in the wrong direction, ie interpolating from the face to the cell center.

Thanks

Shaun