# OpenFOAM - Airfoil Calculations

In this tutorial we will look at incompressible flow over NACA 4-digit airfoils. The goal of this document is:

- To show how simple scripting tools such like Octave can be used together with
*snappyHexMesh*to create advanced meshed, in this case a mesh around a NACA 4-digit arifoil - Use the new edge handling feature in
*snappyHexMesh*to better resolve sharp edges (like the trailing edge of an airfoil) - To use several different solvers in sequence on the same geometry to (hopefully) produce good results with less computational time
- To calculate forces, in this case drag and lift on a airfoil

## Prerequisites

Before you start on this tutorial, please make sure that you have

- A decent laptop or workstation
- A working OpenFOAM 2.x.x installation
- GNU Octave (Matlab will probably also work with minor modifications, although not tested)
- Done the basic
*cavity*and*damBreak*tutorials that is distributed together with OpenFOAM and described in the documentation - Done the spillway tutorial provided here (or in any other way familiarized yourself with
*snappyHexMesh*). - Some basic knowledge about airfoils, drag and lift coefficients etc.

## 1: Setup of basic case

This time we will use the *wingMotion* tutorial as our basic case. This is because this tutorial uses two steps in the solution, first it calculates the stationary solution before it starts the transient solver. The advantage of this procedure is that one easily and effectively gets rid of the effects from non-physical initial conditions and quickly can calculate an averaged stationary flow field. This field is then used as the initial condition for a more computationally demanding transient solver. However, this stationary calculation is not necessarily correct, so one should always use good engineering judgements to assess the stationary solution before proceeding.

Copy the tutorial files over in a folder of your choice, rename the all folders so that their names get the correct meaning. The `airfoil_simpleFoam/system/createPatchDict`

and the `airfoil_pimpleFoam/constant/dynamicMeshDict`

can be deleted, because these functions will not be used here.

**Terminal window**

cp $FOAM_TUTORIALS/incompressible/pimpleDyMFoam/wingMotion ./ mv wingMotion airfoil cd airfoil mv wingMotion2D_snappyHexMesh airfoil_snappyHexMesh mv wingMotion2D_simpleFoam airfoil_simpleFoam mv wingMotion2D_pimpleDyMFoam airfoil_pimpleFoam rm airfoil_simpleFoam/system/createPatchDict rm airfoil_pimpleFoam/constant/dynamicMeshDict

## 2: Geometry and mesh

Before starting with the geometry definitions and mesh generation you might want to read about NACA airfoils at Wikipedia. As the geometry of NACA airfoils is described by polynomial expressions, the coordinates of the airfoil surface can easily be calculated. It is not difficult to create a triangulated surface of this geometry:

To save you for a lot of typing, an Octave-script NACA2STL.m has been made. Download the script into a directory of your choice and make it executable.:

**Terminal window**

chmod +x NACA2STL.m ./NACA2STL.m

These commands will create a STL-file of a NACA0015 airfoil with a chord length of 1 meter. You might need to edit the first line in the script (`#!/usr/bin/octave -qf`

) if your Octave-binary is located in another location than `/usr/bin`

. Put the file in the `airfoil_snappyHexMesh/constant/triSurface`

directory.

The next step is to prepare for the meshing with *snappyHexMesh*. As the airfoil has a sharp edge of significant importance, we must make sure that this is resolved perfectly. To do this we will use the new edge handling feature. The first step is to extract the edges from the STL-file. That is done with the `surfaceFeatureExtract`

-tool:

**Terminal window**

cd airfoil_snappyHexMesh surfaceFeatureExtract -includedAngle 150 -writeObj constant/triSurface/airfoil.stl airfoil

This will automatically detect any edges and produce a file `constant/triSurface/airfoil.eMesh`

containing these. Then *snappyHexMesh* will need to be "informed" about the edges and that it should snap to them. Make sure that the `features`

-section of `system/snappyHexMeshDict`

looks like:

**system/snappyHexMeshDict**

features ( { file "airfoil.eMesh"; level 0; } );

and that under `snapControls`

there is a parameter

**system/snappyHexMeshDict**

//- Highly experimental and wip: number of feature edge snapping // iterations. Leave out altogether to disable. nFeatureSnapIter 10;

### Exercise 2a

The NACA2STL.m script will create a solid for any NACA 4-digit airfoil. All parameters are for simplicity located in the top of the script itself. Try to adjust the parameters and follow the changes. The STL file can always be inspected in ParaView.

### Exercise 2b

Modify the rest of *snappyHexMeshDict* to suit your needs. Change refinement levels, regions, surfaces and layers to create a good mesh. Adjust the mesh resolution to the available computational resources. You will typically create a finer mesh if you have a desktop or workstation available than if you just have a laptop.

### Exercise 2c

Define the boundary conditions for all the necessary patches:

- inlet
- outlet
- topAndBottom
- airfoil

Remove all references to patches not present in this analysis, and make good judgements on the turbulent quantities *k* and *omega*. You can use the CFD-online turbulence properties calculator as an aid. Pick a Reynolds number to use in the simulation and calculate the velocity based on the chord length and kinematic viscosity. The kinematic viscosity in the `wingMotion`

tutorial we started from is 1e-5 m^2/s (if you change this you must change it in both the `airfoil_simpleFoam`

and `airfoil_pimpleFoam`

folders).

For the rest of this tutorial the following values are used, however, you are of course free to use your own choices:

Airfoil |
NACA0015 |
---|---|

Chord length |
1 m |

Span |
1 m |

Angle of attack |
5 degrees |

Reynolds number |
2.00e6 |

Freestream velocity |
2.00 m/s |

Kinematic viscosity |
1e-6 m^2/s |

Fluid density |
1000 kg/m^3 |

This parameters is chosen more or less "randomly" (based on the experience that the NACA0015 airfoil is widely used in research), but you should keep in mind that we are (for the purpose of learning) going to use a stationary solver on the problem to create an initial condition for the transient solver. If you choose a too high angle of attack, vortex shedding and other transient phenomena will occur, and the steady-state solver will never converge.

## 3: Calculation of forces and force coefficients

One of the goals of (nearly) all airfoil calculations are to calculate lift and drag. OpenFOAM handle this very smoothly. Add a section in the bottom of the `controlDict`

-files located in both the `airfoil_simpleFoam/system`

and `airfoil_pimpleFoam/system`

directories:

**system/controlDict**

functions { #include "forceCoeffs" }

If you are familiar with the C/C++ preprocessor `#include`

-statement, you will recognize the syntax used here.

Then create a new file `forceCoeffs`

in both {[airfoil_simpleFoam/system}} and `airfoil_pimpleFoam/system`

with the contents:

**system/forceCoeffs**

forces { type forces; functionObjectLibs ("libforces.so"); outputControl timeStep; outputInterval 1; patches ( "airfoil.*" ); pName p; UName U; rhoName rhoInf; log true; CofR (0.25 0 0); rhoInf 1000; } forceCoeffs { type forceCoeffs; functionObjectLibs ( "libforces.so" ); outputControl timeStep; outputInterval 1; patches ( "airfoil.*" ); pName p; UName U; rhoName rhoInf; log true; liftDir (0 0 1); dragDir (-1 0 0); CofR (0.25 0 0); pitchAxis (0 1 0); magUInf 2.00; rhoInf 1000; lRef 1; Aref 1; }

The contents should be self-explanatory. It specifies the lift and drag direction and the values used when calculating the non-dimensional lift and drag coefficient:

It also calculates a pitching-moment coefficient based on the moment Q around the axis specified by `CofR`

and `pitchAxis`

.

The forces and force coefficients will be recorded every timestep, independent on the `writeInterval`

-setting. The recorded values will be put in two files one `foeceCoeffs/0/forceCoeffs.dat`

and `forces/0/forces.dat`

in the two solver directories.

## 4: Running *simpleFoam*

As previously mentioned *simpleFoam* is a stationary solver. It's job is to find a reasonable initial condition for the transient analysis. Since it is a stationary solver the time steps does not influence the analysis. The default value is a time step of 1 and an end time of 3000. This means that the solver will iterate 3000 times and then stop.

In recent versions of OpenFOAM there is also another way of stopping the *simpleFoam* solver. You can specify the residuals at witch the solver shall stop. This threshold is specified in `system/fvSolution`

:

**system/fvSolution**

SIMPLE { nNonOrthogonalCorrectors 0; pRefCell 0; pRefValue 0; residualControl { p 1e-5; U 1e-5; "(k|omega)" 1e-5; } }

The solver will now run until all the residuals at the start of a new iteration have reached 1e-5 or a maximum of 3000 iterations, whichever comes first. When you specify the residuals you might also want to increase the stop time to a lot higher value, for example 10 000. You should however not set this extremely high, as there is a chance that the simulation never will converge, for example if the problem is unsteady (vortex shedding etc).

You are now ready to start the stationary solver. We won't run this in parallel, as it is a relatively simple and quick calculation:

**Terminal window**

simpleFoam

When the analysis are finished, inspect the results in ParaView. You should also check the force coefficients and make sure they are within reasonable ranges.

Remember that this is a stationary solver only, and the purpose is to find a flow field that we can use as an initial condition for the transient solver. Do not expect accurate results! When the solution is complete, you should verify that the results look reasonable. If they don't, just leave them. *simpleFoam* will not always converge toward the desired solution, especially if the problem in reality is transient with for example vortex shedding.

## 5: Running *pimpleFoam*

When *simpleFoam* has finished, you must copy the mesh and boundary conditions from the `airfoil_simpleFoam`

directory to the `airfoil_pimpleFoam`

directory:

**Terminal window**

cd airfoil_pimpleFoam cp -r ../airfoil_simpleFoam/constant/polyMesh/* constant/polyMesh/ cp -r ../airfoil_simpleFoam/0.org ./ cp 0.org/* 0/

Now here comes the crux of this tutorial: mapping the results from *simpleFoam* into an initial condition for *pimpleFoam* to start on. If the results from *simpleFoam* for any reason does not look physical or in any other way seems unusable, skip this step (*pimpleFoam* will just need a bit longer calculation time):

**Terminal window**

mapFields ../airfoil_simpleFoam -sourceTime latestTime -consistent

Now you can open the case in Paraview, and inspect the initial flow field.

The last and final step is to decompose the case and run *pimpleFoam*:

**Terminal window**

decomposePar mpirun -np 4 pimpleFoam -parallel

That's it! You can now plot the lift and drag coefficients with gnuplot or similar tools.

## 6: Download

You can download all necessary case files in airfoil.tar.gz. Scripts to run and clean the case is also attached.

## 7: Alternative meshing

If you think the *snappyHexMesh* and *extrudeMesh* procedure is too cumbersome, it is actually not very difficult to create a C-mesh around the foil just by using *blockmesh*:

We have again prepared a little octave script for you that will create the necessary `blockMeshDict`

-file for you. You should experiment with expansion rates, resolutions and other settings to get the best possible mesh. This meshing is a general procedure, and will work best with foils at low angles of attack.