# Rayleigh–Bénard Convection Using buoyantBoussinesqPimpleFoam

Here is an extremely simple simulation to set up that has a surprisingly beautiful output. In this post, we will simulation the classic Rayleigh–Bénard convection (see Wikipedia) in 3D using the buoyant solver, buoyantBoussinesqPimpleFoam.

buoyantBoussinesqPimpleFoam is a solver for buoyant, turbulent, incompressible flows. The incompressible part of the solver comes from the fact that it uses the Boussinesq approximation for buoyancy which is only suitable for flows with small density changes (aka incompressible). A good source for the equations in this solver is this blog post.

## Simulation Set-up

The basic set-up for this case is simple: a hot bottom surface, a cold top surface, either cyclic or zero-gradient sides, and some initial conditions.

For this example, I used the properties of water, I set the bottom plate at 373 K (don’t ask me why… I know it’s close to boiling point of water), and the top plate at 273 K. For this case, we will not use any turbulent modeling and will simply use a laminar model (this simply means there is only molecular viscosity, there are no simplifications applied to the equations).

### Geometry and Mesh

The geometry is shown below. As the geometry is so simple… I will not go over the blockMesh set up. The mesh discretization that I used was simplegrading (i.e. no inflation), with 200x200x50 cells.

### Constant

For this case, we will simulate water. The transportProperties file should look like:

transportModel Newtonian;

// Laminar viscosity
nu [0 2 -1 0 0 0 0] 1e-06;

// Thermal expansion coefficient
beta [0 0 0 -1 0 0 0] 0.000214;

// Reference temperature
TRef [0 0 0 1 0 0 0] 300;

// Laminar Prandtl number
Pr [0 0 0 0 0 0 0] 7.56;

// Turbulent Prandtl number (not used)
Prt [0 0 0 0 0 0 0] 0.7;

I don’t typically delete unused entries from dictionaries. This makes using previous simulations as templates much easier. Therefore note that the turbulent Prandtl number is in the dictionary but it is not used.

Selecting Reference Temperature TRef for buoyantBoussinesqPimpleFoam. To answer this recall that when the Boussinesq buoyancy approximation is used, the solver does not solve for the density. It solves the relative density using the linear function:

$\frac{\rho}{\rho_0}=1-\beta \left(T-T_0\right)$

Therefore, I think it makes sense that we should choose a temperature for $T_{ref}$ that is somewhere in the range of the simulation. Thus I chose Tref=300. Somewhere in the middle!

And the turbulenceProperties file is:

simulationType laminar;

RAS
{
RASModel laminar;

turbulence off;

printCoeffs off;
}

The g file tells the solver the acceleration due to gravity, as well as the direction:

dimensions [0 1 -2 0 0 0 0];
value (0 -9.81 0);

### Boundary Conditions

In the “zero” folder, we need the following files: p, p_rgh, T, U, and alphat (this file needs to be present… however it is not used given the laminar simulationType.

T:

dimensions [0 0 0 1 0 0 0];

internalField uniform 273;

boundaryField
{
floor
{
type fixedValue;
value uniform 373;
}
ceiling
{
type fixedValue;
value uniform 273;
}
fixedWalls
{
}
}

p_rgh:

dimensions [0 2 -2 0 0 0 0];

internalField uniform 0;

boundaryField
{
floor
{
type fixedFluxPressure;
rho rhok;
value uniform 0;
}

ceiling
{
type fixedFluxPressure;
rho rhok;
value uniform 0;
}

fixedWalls
{
type fixedFluxPressure;
rho rhok;
value uniform 0;
}
}

p:

dimensions [0 2 -2 0 0 0 0];

internalField uniform 0;

boundaryField
{
floor
{
type calculated;
value $internalField; } ceiling { type calculated; value$internalField;
}

fixedWalls
{
type calculated;
value \$internalField;
}
}

U:

dimensions [0 1 -1 0 0 0 0];

internalField uniform (0 0 0);

boundaryField
{
floor
{
type noSlip;
}

ceiling
{
type noSlip;
}

fixedWalls
{
type noSlip;
}
}

## Simulation Results

The results (as always) are the best part. But especially for this case since they are so nice to look at! I have made a couple animations of temperature fields and contours. Enjoy.

## Conclusion

This case demonstrated the simple set up of a case using buoyantBoussinesqPimpleFoam. The case (Rayleigh-Bénard convection) was simulated in 3D on a fine grid.

Comments and questions are welcome! Keep in mind I set this case up very quickly.

# Time-Varying Cylinder Motion in Cross-flow: timeVaryingFixedUniformValue

This post is a simple demonstration of the timeVaryingFixedUniformValue boundary condition. This boundary condition allows a Dirichlet-type boundary condition to be varied in time. To demonstrate, we will modify the oscillating cylinder case.

## Set-Up

Instead of using the oscillating boundary condition for point displacement. We will have the cylinder do two things:

• Move in a circular motion
• Move in a sinusoidal decay motion

The basics of this boundary condition are extremely simple. Keep in mind that although (here) we are modifying the pointDisplacement boundary condition for the cylinder, the basics of this BC would be the same if you were doing a time varying boundary condition for say pressure or velocity.

In the pointDisplacement file:

 cylinder
{
type timeVaryingUniformFixedValue;
fileName "prescribedMotion";
outOfBounds clamp;
}

fileName points to the file where the time varying boundary condition is defined. Here we used a file called prescribedMotion however you can name it whatever you want. The outOfBounds variable dictates what the simulation should do if the simulation time progresses outside of the time domain defined in the file.

The additional file containing the desired motion prescribedMotion is formatted in the following way:

(
( 0 (0 0 0))
( 0.005 (-0.0000308418795848531 0.00392682932795517 0))
( 0.01 (-0.0001233599085671 0.00785268976953207 0))
( 0.015 (-0.000277531259507496 0.0117766126774107 0))
( 0.02 (-0.00049331789293211 0.0156976298823283 0))
...
( 9.99 (-0.0001233599085671 -0.00785268976953189 0))
( 9.995 (-0.0000308418795848531 -0.00392682932795479 0))
( 10 (0 -3.06161699786838E-016 0))
)

The first column is the time in seconds, and the vector defines the point displacement. In the present tutorial, these points were calculated in libreOffice and then exported into the text file.  I arbitrarily made up the motions purely for the sake of making this blog post.

The circular motion was defined as:

$x=0.25\cos\left(\pi t\right)-0.25$ and $y=0.25\sin\left(\pi t\right)$

Decaying sinusoidal motion was:

$y=\sin(\pi t) \exp\left(-t/2\right)$

The rest of the set-up is identical to the set-up in the oscillating cylinder example. The solver pimpleDyMFoam is then run.

## Conclusions

This post demonstrated how a more complicated motion can be prescribed by using a little math and the timeVaryingUniformFixedValue boundary condition. Always like to hear questions and comments! Has anybody else done something like this?

# Equations for Steady 1D Isentropic Flow

The equations used to describe steady 1D isentropic flow are derived from conservation of mass, momentum, and energy, as well as an equation of state (typically the ideal gas law).

These equations are typically described as ratios between the local static properties (p, T, $\rho$) and their stagnation property as a function of Mach number and the ratio of specific heats, $\gamma$. Recall that Mach number is the ratio between the velocity and the speed of sound, a.

These ratios are given here:

Temperature: $T_o/T = \left(1+\frac{\gamma -1}{2} M^2\right)$

Pressure: $P_o/P = \left(1+\frac{\gamma -1}{2} M^2\right)^{\frac{\gamma}{\gamma-1}}$

Density: $\rho_o/\rho = \left(1+\frac{\gamma -1}{2} M^2\right)^{\frac{1}{\gamma-1}}$

In addition to the relationships between static and stagnation properties, 1D nozzle flow offers an equation regarding the choked cross-sectional flow area (recall that the flow is choked when M=1.)

$A/A^* = \frac{1}{M}\left(\left(\frac{2}{\gamma+1}\right)\left(1+\frac{\gamma -1}{2} M^2\right)\right)^{\frac{\gamma+1}{2\left(\gamma-1\right)}}$

Some excellent references for these equations are:

• Gas Dynamics Vol. I – Zucrow and Hoffman – 1976
• Gas Dynamics – John and Keith – 2nd Ed. – 2006

# Establishing Grid Convergence

Establishing grid convergence is a necessity in any numerical study. It is essential to verify that the equations are being solved correctly and that the solution is insensitive to the grid resolution. Some great background and details can be found from NASA:

https://www.grc.nasa.gov/WWW/wind/valid/tutorial/spatconv.html

First, here is a summary the equations and steps discussed here (in case you don’t want to read the whole example):

1. Complete at least 3 simulations (Coarse, medium, fine) with a constant refinement ratio, r, between them (in our example we use r=2)
2. Choose a parameter indicative of grid convergence. In most cases, this should be the parameter you are studying. ie if you are studying drag, you would use drag.
3. Calculate the order of convergence, p, using:
• $p=\ln(\frac{(f_3-f_2)}{(f_2-f_1)}) / \ln(r)$
4. Perform a Richardson extrapolation to predict the value at h=0
• $f_{h=0}=f_{fine}+\frac{f_1-f_2}{r^p-1}$
5. Calculate grid convergence index (GCI) for the medium and fine refinement levels
• $GCI=\frac{F_s |e|}{r^p-1}$
6. Ensure that grids are in the asymptotic range of convergence by checking:
• $\frac{GCI_{2,3}}{r^p \times GCI_{1,2}} \approxeq 1$

So what is a grid convergence study? Well, the gist of it is that you refine the mesh several times and compare the solutions to estimate the error from discretization. There are several strategies to do this. However, I have always been a fan of the following method: Create a very fine grid and simulate the flow problem. Then reduce the grid density twice, creating a medium grid, and coarse grid.

There are several strategies to do this. However, I have always been a fan of the following method: Create a very fine grid and simulate the flow problem. Then reduce the grid density twice, creating a medium grid, and coarse grid. To keep the process simple, the ratio of refinement should be the same for each step. ie. if you reduce the grid spacing by 2 in each direction between the fine and medium grid, you should reduce it again by 2 between the medium and coarse grid. In the current example, I generated three grids for the cavity problem with a refinement ratio of 2:

• Fine grid: 80 cells in each direction – (6400 cells)
• Medium grid: 40 cells in each direction – (1600 cells)
• Coarse grid: 20 cells in each direction – (400 cells)

Velocity contour plots are shown in the following figures:

We can see from the figures that the quality of the simulation improves as the grid is refined. However, the point of a grid convergence study is to quantify this improvement and to provide insight into the actual quality of the fine grid.

The accuracy of the fine grid is then examined by calculating the effective order of convergence, performing a Richardson extrapolation, and calculating the grid convergence index. As well, (as stated in the article from NASA), it is helpful to ensure that you are in the asymptotic range of convergence.

What are we examining?

It is very important at the start of a CFD study to know what you are going to do with the result. This is because different parameters will converge differently. For example, if you are studying a higher order parameter such as local wall friction, your grid requirements will probably be more strict than if you are studying an integral (and hence lower order) parameter such as coefficient of drag. You only need to ensure that the property or parameter that you are studying is grid independent. For example, if you were studying the pressure increase across a shockwave, you would not check that wall friction somewhere else in the simulation was converged (unless you were studying wall friction as well).  If you want to be able to analyze any property in a simulation, both high and low order, then you should do a very rigorous grid convergence study of primitive, integrated and derived variables.

In our test case lets pretend that in our research or engineering project that we are interested in the centerline pressure and velocity. In particular, let’s say we are interested in the profiles of these variables along the centerline as well as their peak values. The centerline profiles of velocity and pressure are shown in the following figures:

Calculate the effective order of convergence

From our simulations, we have generated the following data for minimum pressure along the centerline, and maximum velocity along the centerline:

The order of convergence, p, is calculated with the following equation:

$p=\ln(\frac{(f_3-f_2)}{(f_2-f_1)}) / \ln(r)$

where r is the ratio of refinement, and f1 to f3 are the results from each grid level.

Using the data we found, p = 1.84 for the minimum pressure and p=1.81 for the maximum velocity.

Perform Richardson extrapolation of the results

Once we have an effective order, p, we can do a Richardson extrapolation. This is an estimate of the true value of the parameter we are examining based on our order of convergence. The extrapolation can be performed with the following equation:

$f_{h=0}=f_{fine}+\frac{f_1-f_2}{r^p-1}$

recall that r is the refinement ratio and is $h_2/h_1$ which in this case is 2.

Using this equation we get the Richardson extrapolated results:

• $P_{min}$ at h=0  -> -0.029941
• $V_{max}$ at h=o  -> 0.2954332

The results are plotted here:

Calculate the Grid Convergence Index (GCI)

Grid convergence index is a standardized way to report grid convergence quality. It is calculated at refinement steps. Thus we will calculate a GCI for steps from grids 3 to 2, and from 2 to 1.

The equation to compute grid convergence index is:

$GCI=\frac{F_s |e|}{r^p-1}$

where e is the error between the two grids and $F_s$ is an optional (but always recommended) safety factor.

Now we can calculate the grid convergence indices for the minimum pressure and maximum velocity.

Minimum pressure

• $GCI_{2,3} = 1.25 \times |\frac{-0.028836-(-0.025987)}{-0.028836}|/(2^{1.84}-1) \times 100 \% = 4.788 \%$
• $GCI_{1,2} = 1.25 \times |\frac{-0.029632-(-0.028836)}{-0.029632}|/(2^{1.84}-1) \times 100 \% = 1.302 \%$

Max velocity

• $GCI_{2,3} =1.25 \times | \frac{0.2892-0.27359}{0.2892}|/(2^{1.81}-1) \times 100 \% =2.69187 \%$
• $GCI_{1,2} = 1.25 \times |\frac{0.29365-0.2892}{0.29365}|/(2^{1.84}-1) \times 100 \% = 0.7559 \%$

Check that we are in the asymptotic range of convergence

It is also necessary to check that we are examing grid converegence within the asymptotic range of convergence. If we are not in the asymtotic range this means that we are not asymptotically approaching a converged answer and thus our solution is definitely not grid indipendent.

With three grids, this can be checked with the following relationship:

$\frac{GCI_{2,3}}{r^p \times GCI_{1,2}} \approxeq 1$

If we are in the asymptotic range then the left-hand side of the above equation should be approximately equal to 1.

In our example we get:

Minimum Pressure

$1.0276 \approxeq 1$

Minimum Velocity

$1.0154 \approxeq 1$

Applying Richardson extrapolation to a range of data

Alternatively to choosing a single value like minimum pressure or maximum velocity. Richardson extrapolation can be applied to a range of data. For example, we can use the equation for Richardson extrapolation to estimate the entire profile of pressure and velocity along the centerline at h=0.

This is shown here:

In this post, we used the cavity tutorial from OpenFOAM to do a simple grid convergence study. We established an order of convergence, performed Richardson extrapolation, calculated grid convergence indices (GCI) and checked for the asymptotic range of convergence.

As I said before, the NASA resource is very helpful and covers a similar example:

As well the papers by Roache are excellent reading for anybody doing numerical analysis in fluids:

# The Ahmed Body

The Ahmed body is a geometric shape first proposed by Ahmed and Ramm in 1984. The shape provides a model to study geometric effects on the wakes of ground vehicles (like cars).

Image highlights:

In this post, I will use simpleFoam to simulate the Ahmed body at a Reynolds number of 10^6 using the k-omega SST turbulence model. The geometry was meshed using cfMesh which I will briefly discuss as well. Here is a breakdown of this post:

1. Geometry Definition
2. Meshing with cfMesh
3. Boundary Conditions
4. Results

Note: I ran this case on my computer with 6 Intel – i7 (3.2 GHz) cores and 32 Gb of RAM. Throughout the simulation about 20-ish gigs of RAM were used.

## Geometry Definition

### STL Creation

The meshing utility cfMesh is similar to snappyHexMesh in that it depends on a geometry file of some type (.stl etc) to create the mesh. But it is different in that the entire domain must be part of the definition.

The ahmed body geometry can be found: http://www.cfd-online.com/Wiki/File:Ahmed.gif

For this simulation, I generated the geometry using SolidWorks. But this wasn’t for any particular reason other than that it was quick since I am familiar with it.

### Preparation for Meshing

Once you have an STL file, you could go straight ahead to meshing it with cfMesh. However, some simple preparations to the STL geometry can improve the quality of the mesh created, and make setting up the case easier.

In particular, when you create the STL file in SolidWorks (or your 3D modeller of choice)  it contains no information about the boundaries and patches. As well, cfMesh works best if the geometry is defined using a .fms or .ftr file format.

Use surfaceFeatureEdge utility to extract edge information and create a .ftr file. Firstly let’s extract edge and face information from our STL file. We also define an angle. This angle tells cfMesh that any angle change large than this (in our case I chose 20 degrees) is a feature edge that must be matched.

surfaceFeatureEdge volume.stl volume.ftr -angle 20

After we run this, the new file volume.ftr contains a bunch of face and patch information. 13 surfaces with feature edges were extracted. The first 6 (volume_0 to volume_5) are the boundaries of the simulation (inlet, outlet, ground, front, back, and top).

mergeSurfacePatches volume.ftr ahmed -patchNames '(volume_6 volume_7 volume_8 volume_9 volume_10 volume_11 volume_12)'

After running this command, volume.ftr now contains 7 patches. We are now ready to move on to setting up cfMesh.

## Meshing with cfMesh

### Set up meshDict file in the system folder

Similar to snappyHexMesh and blockMesh, cfMesh using a dictionary file to control its parameters; meshDict. In this dictionary file we will be modifying a few parameters.

Tell cfMesh what file is to be meshed:

surfaceFile "volume_transformed.ftr";

Set the default grid size:

maxCellSize 0.2;

Set up refinement zones:

We want to set up two refinement zones; a larger one to capture most of the flow further away from the body (including the far wake), and a smaller more refined one to capture the near wake and the flow very close to the Ahmed body.

objectRefinements
{
box1
{
cellSize 25e-3;
type box;
centre (2.4 0 6.5);
lengthX 1;
lengthY 1;
lengthZ 3.5;
}

box2
{
cellSize 5e-3;
type box;
centre (2.4 0 7);
lengthX 0.5;
lengthY 1;
lengthZ 2;
}

}

Set up boundaries to be renamed:

renameBoundary
{
newPatchNames
{
volume_0
{
newName ground;
type wall;
}
volume_1
{
newName back;
type wall;
}
volume_2
{
newName inlet;
type patch;
}
volume_3
{
newName front;
type wall;
}
volume_4
{
newName outlet;
type wall;
}
volume_5
{
newName top;
type patch;
}
ahmed
{
newName ahmed;
type wall;
}
}

Set up boundary layering:

We require boundary layer on both the Ahmed body, as well as the volume_0 patch. Recall that the Ahmed body is surface mounted!

boundaryLayers
{
patchBoundaryLayers
{
ahmed
{
nLayers 10;
thicknessRatio 1.1;
maxFirstLayerThickness 5e-3;
}
volume_0
{
nLayers 10;
thicknessRatio 1.05;
maxFirstLayerThickness 10e-3;
}

}
}

### Run cfMesh

We want to create a hex-dominant grid. This means that the 3D grid will consist primarily of hex cells. To achieve this we will use the cartesianMesh solver from cfMesh.

The results are shown below. The final mesh consisted of approximately 16.9 million cells. The majority of the cells were hexahedra (approximately 99%).

## Boundary Conditions for the Solver

For this case, we are going to run a steady-state RANS simulation using the kwSST model and the solver simpleFOAM. This is simply to demonstrate the running of the solver.

The boundary conditions used are summarized in the following table:

As you can see I have used wall functions for the wall boundary conditions. This is due to the very small cell requirements that would be required to resolve the boundary layer on the ground, as well as on the Ahmed body which is at a Reynolds number of one million.

## Simulation Results

The simulation took about a day and a half on 6 cores. Throughout the simulation, about 20 gb of RAM was used.

## Conclusions

In this post, we meshed and simulated a surface-mounted Ahmed body at a Reynolds number of one million. We meshed it using the open-source meshing add-on cfMesh. We then solved it as a steady-state RANS simulation using the kwSST turbulence model, and the simpleFOAM solver.

The results gave some nice figures and a qualitatively correct result! And it was pretty fun. cfMesh was extremely easy to use and required much less user input than its OpenCFD counterpart snappyHexMesh.

## Some references:

http://www.cfd-online.com/Wiki/Ahmed_body

Some papers studying the Ahmed body:

-See the reference on the above CFD Online page!

As usual please comment and let me know what you think!

Cheers,

curiosityFluids

# Oscillating Cylinder in Laminar Crossflow – pimpleDyMFoam

In this post I am going to simulate an oscillating cylinder in a cross-flow… just for fun… and to provide an additional tutorial case for those wishing to use some of the dynamic meshing features of OpenFOAM.

The case I am going to simulate is a cylinder in a Reynolds number 200 cross-flow (U=2 m/s, D=1 m, nu = 0.01 m^2/s), oscillating at a rate of 0.2 hz.

## Tutorial Files

Please let me know if the download does not work or if there is a problem running the tutorial files. Note: I ran this case in parallel on a relatively fast 6-core computer. It will take a long time if you try to run it as single core.

## Mesh Generation

For this simulation, I built a simple two dimensional O-grid around the cylinder which joined to a series of blocks making a rectangular domain. I built this using the native OpenFOAM meshing utility blockMesh (which I like a lot).

If you are wondering about the blockMesh set up… I intend to a blockMesh tutorial post… not that it’s all that necessary since the OpenFOAM manual covers it pretty well.

## Case Set-up

### dynamicMeshDict

When running a dynamic mesh case a solver runs as part of the solution and solves for the new grid at each timestep. Several are available in OpenFOAM and I am not really trying to do a full post on that right now. So I’ll just tell you that in this case I decided to use the displacementLaplacian solver.

Along with the solver one must define the coefficients that go with the solver. For most of the solvers this means setting the diffusivity for the Laplace solver. Since I am relatively new to these types of solutions I did what we all do… and turned to the great CFD online forum! A discussion in this post (http://www.cfd-online.com/Forums/openfoam-meshing/97299-diffusivity-dynamicmeshdict.htm) made me think that a good starting point would be the inverseDistance model for the mesh diffusivity.

In the end my dynamicMeshDict (which belongs in the constant folder) looked like:

dynamicFvMesh dynamicMotionSolverFvMesh;

motionSolverLibs ( "libfvMotionSolvers.so" );

solver displacementLaplacian;

displacementLaplacianCoeffs
{
diffusivity inverseDistance 1(cylinder);
}

The solver produces the grid movements through time as the simulation is performed. Here is what the mesh looks like while it is moving! I added the stationary mesh and the moving mesh. Click on them to get a clearer picture 🙂

### Boundary Conditions

The set up for this case is simple. An inlet boundary (on the far left) where I specified the incoming velocity and pressure (both are uniform fixedValue). The top and bottom boundaries are slip boundaries (but could also be freestream depending on your preferred set-up style).

The cylinder itself requires some special treatment because of the moving mesh. It must obey the no-slip condition right? So do we set it as uniform (0 0 0) ? No! The cylinder is moving! Luckily there is a handy boundary condition for this in the U file:

 cylinder
{
type movingWallVelocity;
value uniform (0 0 0);
}

For a typical simulation using pimpleFoam a p and U file would be all that are required. However, since we are doing a moving mesh simulation there is another parameter that must be solved for and requires boundary conditions… pointDisplacement.

For the pointDisplacement boundary conditions, we know that all of the outer edges should NOT move. Therefore they are all fixed with a type of fixedValue and  a value of uniform (0 0 0).

The cylinder however is moving and requires a definition. In this simulation we are simulating and oscillating cylinder. Since we are using the displacement solver the type is oscillatingDisplacement. We input and omega (rad/s) and an amplitude (m) in the following way:

 cylinder
{
type oscillatingDisplacement;
omega 1.256;
amplitude (0 0.5 0); // Max piston stroke
value uniform (0 0 0);
}

## Results

Yay! Now its time to look at the results. Well since I am not doing this for any particular scientific study… let’s look at some pretty pictures!

Here is an animation of vorticity:

Looks pretty nice! I personally have a big nerd love for vortex shedding…. I don’t know why.

Obviously if you intend to do any scientific or engineering work with this type of problem you would need to think very carefully about the grid resolution, diffusivity model, Reynolds number, oscillation frequency  etc. All of these were arbitrarily selected here to facilitate the blog post and to provide a nice tutorial example!

## Conclusion

In this post I briefly covered the set-up of this type of dynamic meshing problem. The main difference for running a dynamic mesh case is that you require a dynamic mesh solver (you must specify in the dynamicMeshDict), and you also require boundary conditions for that solver.

Let me know if there are any problems with this blog post or with the tutorial files provided.