# 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?

# 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

# Pressure Driven Nozzle Flow with Shock – rhoCentralFoam

In this post I will go over the set up and solution of a pressure driven nozzle flow with a shock located in the diverging section. This refers to the type of flow problem described by region b in my page covering stationary normal shock-waves.