• OpenFOAM - Airfoil Calculations
  • Prerequisites
  • 1: Setup of basic case
  • 2: Geometry and mesh
    • Exercise 2a
    • Exercise 2b
    • Exercise 2c
  • 3: Calculation of forces and force coefficients
  • 4: Running simpleFoam
  • 5: Running pimpleFoam
  • 6: Download
  • 7: Alternative meshing

OpenFOAM - Airfoil Calculations

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

1. 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
2. Use the new edge handling feature in snappyHexMesh to better resolve sharp edges (like the trailing edge of an airfoil)
3. To use several different solvers in sequence on the same geometry to (hopefully) produce good results with less computational time
4. 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.

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.:

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:

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:

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

and that under snapControls there is a parameter

//- 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:

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:

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:

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:

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:

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:

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):

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:

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.