# OpenFOAM : Run-time post-processing

Team

ALILOU Youssef

​ER-RAIY Aimad

GRUSS Maxime

#### Supervisor

PEDRONO Annaig

Keywords: OpenFOAM    run-time    post-processing    post-treatment    F-14 Tomcat    CFD     pyFoam  snappyHexMesh    blockMesh    wind tunnel    aerodynamism

# Introduction

## Introduction

OpenFOAM (Open Source Field Operation and Manipulation) CFD toolbox is a free CFD software package writen in C++, originally developed at the Imperial College of London in the late 1980s, in order to have a more powerful and flexible simulation tool than the classical one for fluid dynamics simulations, which is the FORTRAN language. Formerly sold under the name FOAM, it has been released as open source in 2004 and is currently (06/2013) owned by OpenCFD Ltd.

Thanks to the fact that it is free of charge and extensible, it is already used in private Design Offices (for example at Audi, Mitsubishi, Shell...) as well as in academic organizations all over the world and is likely to become a model for future CFD softwares. The sustainability of the software is funded by support, contracted developments and the training courses provided by the foundation.

Its features1 allow to solve problems in most areas of engineering and science, from complex multi-fluid flows with chemical reactions, turbulence and heat transfer, to solid dynamics and electromagnetic, but also offers the user a complete freedom of customization, with the support from OpenCFD if required.

Structure:

The working base is the package available on the OpenFOAM website2 for Linux interfaces-currently the 2.2.0 version-, and contains the following extensible list of objects:

•  a large base library, containing the physical models such as transport, thermophysical models and turbulence models -RAS, LES...-, reaction kinetics, etc...
•  many solvers applications (over 80), applied to specific cases, relying on the library, from the simplest such as potentialFoam used for potential flows, to simpleFoam for incompressible and stationary, icoFoam for compressible flows, interFoam for multiphase phase flows, etc...
• utility applications (over 170, performing pre- and post-processing tasks: mesh generation (blockMesh that creates a mesh out of a text configuration file), mesh conversion -from CAD softwares for example-, parallel processing utilities, data visualization and conversion units such as foamToVTK that converts from the native OpenFOAM format to VTK format, in order to be able to visualize the results of a computation on the open-source application Paraview3.

In this project we will focus in particular on Run-time post-processing possibilities, which are post-processing methods, ran in parallel with the computation.

Run-time post processing on OpenFOAM

This part aims to show how to apply a run-time post processing part to a script. The interest of doing so is to allow the user following data during the calculation, indeed the time required can be huge, it can take minutes, hours, days, months, centuries, thus it is crucial to know if the results are coherent, according to parameters such as the residuals; how much time it will need, know the progress of the calculation and even visualize the current results which would allow an insider to make a critical analysis. Here are some examples of available classic functionalities:

fieldAverage - temporal averaging of fields.
forces - calculates pressure/viscous forces and moments.
forceCoeffs - calculates lift, drag and moment coefficients.
isoSurface - generation of an isosurface of given fields in one of the standard sample formats, e.g. VTK.
systemCall - execute any system call, e.g. email you to tell you your job is finished.
StreamLines - generates streamlines in one of the sample formats.

The case study chosen here consists of a F-14 Tomcat aircraft, made with a CAD software.
The making of the case will be rapidly described, then we will focus on the run-time post processing part.

The project files are avalaible to download here.

# Setting of the problem

## Setting of the problem

Wind tunnels are a classic kind of tools often used in aerodynamics to test  the physical models on aircrafts, at a smaller scale, and study shapes. In this kind of installations, engineers can change the conditions of the flow which affects the plane. By measuring the forces applied on the model, it allows to estimate the forces which will occur at a realistic scale and in real conditions. With technical diagnosis, engineers can understand and enhance the aircraft performances.

Even if cheaper than real flight tests, those measures are still expensive, time-consuming and harder to set up than numerical simulations, that's why computational modelisations are often chosen to study aerodynamics.

In this project, we chose to modelize a wind tunnel in which is placed an aircraft (F-14 Tomcat), in order to study its aerodynamics, this will provide a large case which will take a considerable time to be done and will give the opportunity to explore some features of run-time post-processing.

# Case setup

## Case setup

Here we are going to create a similar simulation to the wind tunnel experiment, this will consist in the study of an external flow surrounding an aircraft. According to the velocity, the presence of a wall below the plane and of the landing gear we can assimilate the study to a phase of the flight after the take off or before the landing.
For this case we will solve the problem with the simpleFoam solver, which is used for  incompressible turbulent simulations, using the k-omega model for turbulence, and there is no surface friction on the aircraft.

We chose a Mach number equal to 0.29, to ensure the incompressibility hypothesis is valid (which corresponds to 95.2m/s). The fluid used is air, the parameters used are the standard ones at 25°C.

# Geometry

## Geometry

We create geometry of an F-14 Tomcat with a CAD software (Free 30 days Trial version), with a similar method than the one explained in the tutorial on this link.

The blueprints used in the sketch are available on this link :

Blueprints

At the end we have a geometry which looks like this:

Top view

Rear view                                                                                             Side view

Exploded view

Then we export the geometry as an STL file (STereoLithography also known as Standard Tessellation Language) which describes the model associated raw unstructured triangulated surface, without any representation of color, texture or other common CAD model attributes.

An STL file be represented by either Binary or ASCII .In our case we export it as an ASCII-STL that we put in constant/triSurface subdirectory.

We will build our case on the basis of an existing OpenFOAM tutorial, that we will modify to configure our case .

First we copy the motorbike tutorial from the Tutorials directory and name it Aircraft

cp -r \$FOAM_TUTORIALS/mesh/snappyHexMesh/motorBike Aircraft

cd Aircraft 

# Mesh generation

Mesh generation

OpenFOAM has two native mesh generators: blockMesh and snappyHexMesh. The former is used to generate block-structured grids from vertices defined in blockMeshDict, whereas the later generates 3-dimensional meshes containing hexahedra (hex) and split-hexahedra (split-hex) automatically from triangulated surface geometries in (STL) format.

To generate the mesh with the snappyHexMesh tool we proceed as follows:

First, we open our stl geometry in paraview , in order to check the dimensions of the geometry before setting the block of the wind tunnel (the block generated by blockMesh) dimensions, in order to fit the aircraft into that mesh and suppress numerical errors, caused by a too small mesh:

Now that we have got informed on our geometry bounds, we can specify the vertices that define the background mesh, so that it could contain the geometry. We modify then the vertices on  /constant/blockMeshDict  so that they become:

vertices
(

(-200.319 -266.211 -27.5135)

( 311.603 -266.211 -27.5135)

( 311.603 -125.092 -27.5135)

(-200.319 -125.092 -27.5135)

(-200.319 -266.211 143.683)

( 311.603 -266.211 143.683)

( 311.603 -125.092 143.683)

(-200.319 -125.092 143.683)

);

We set the boundary conditions as they are without changing them , and do the same for blocks, It is not necessary to create a fine background mesh, since it will be refined by snappyHexMesh.

Then we run : blockMesh to generate the hexahedral background mesh.

We could verify that the block fits the geometry by visualizing them both on ParaFoam:

Next step is to modify or create the snappyHexMesh input file. One can create or just modify  snappyHexMeshDict input file in  the system sub-directory of our case. This file includes parameters at the top level that control the various stages of the meshing process.

The process of generating a mesh using snappyHexMesh consists of three stages :

1. Cell splitting at feature edges and surfaces : Start the splitting from locationInMesh feature. This edge must be inside the region to be meshed and must no coincide with a cell, and Splitts the cells around the surface according to refinementSurfaces

2. Snapping to surfaces: involves moving cell vertex points onto surface geometry to remove the jagged castellated surface from the mesh.

3.  Mesh layers addition : involves shrinking the existing mesh from the boundary and inserting layers of cells.

The splitting process begins with cells being selected according to specified edge features first within the domain

Thus, we should start by generating feature edges by editing the surfaceFeatureExtract Dictionnary (for 2.2.x) or by executing the line command.

surfaceFeatureExtract -includedAngle 150 surface.stl features (for former versions)

As we are working on the 2.2.0 of OF we edit the surfaceFeatureDict by substituting the line that corresponds to “motorbike.obj” by “aircraft.stl” so that it looks :

aircraft.stl
{
extractionMethod    extractFromSurface;

extractFromSurfaceCoeffs
{

includedAngle   150;
}

subsetFeatures
{
nonManifoldEdges       no;

openEdges       yes;
}

writeObj                yes;
}

Then we run surfaceFeatureExtract to generate features in .eMesh format.

Now we open and edit the snappyHexMeshDict. We make sure that all stages (castellatedMesh , snap and addLayers) are set to True value.

Then we modify it so that the section after “geometry” keyword becomes :

geometry

{

aircraft.stl //our STL geometry file located in constant/triSurface

{

type triSurfaceMesh;

name aircraft; //object name

}

refinementBox

{

type searchableBox; //the zone at which refinement starts

min (-200 -300 30);

max (50 -200 100);

}

};

Features will be captured from and *.eMesh file located in /constant/extendedFeatureEdgeMesh that will have generated by surfaceFeatureExtract execution .

features

(

{

file "aircraft.eMesh";

level 0;

}

);

Then we specify the surfaces around which the splitting is going to be around, which is here our geometry surface.  
refinementSurfaces

{

aircraft

{

level (5 6);//

patchInfo

{

type wall;

inGroups (aircraft);//aircraft is the first line of our ASCII-STL file,
//it is the patch name

}

}

}

5 level is applied accross the surface an 6 is applied to cells that can see  intersections that form an angle in excess of that specified by resolveFeatureAngle>30, the name corresponding to what we should put in inGroups() corresponds to the first line on our STL file.

The refinementRegions contains entries for refinement of the volume regions specified in the geometry sub-dictionary. A refinement mode is applied to each region which can be :

• inside refines inside the volume region.

• outside refines outside the volume region.

• distance refines according to distance to the surface; and can accommodate different levels at multiple distances with the levels keyword.​

 refinementRegions
{
aircraft
{
mode distance;
levels ((0.1 5) (0.4 4) (1 2));
}
}

Finally, we set locationInMesh to (-120 -180 30);

The next stage is snapping the surface which is done in 3 steps which are:

1. Snap the vertices onto the STL surface.
2. Solve for relaxation of the internal mesh.
3. Iterate using the snapControls in SnappyHexMeshDict.

We set snapControls to their default values  :

// Settings for the snapping.
snapControls
{
//- Number of patch smoothing iterations before finding correspondence
//  to surface
nSmoothPatch 3;

//- Relative distance for points to be attracted by surface feature point
//  or edge. True distance is this factor times local
//  maximum edge length.
tolerance 4.0;

//- Number of mesh displacement relaxation iterations.
nSolveIter 0;

//- Maximum number of snapping relaxation iterations. Should stop
//  before upon reaching a correct mesh.
nRelaxIter 5;

// Feature snapping

//- Number of feature edge snapping iterations.
//  Leave out altogether to disable.
nFeatureSnapIter 0;

//- Detect (geometric only) features by sampling the surface
//  (default=false).
implicitFeatureSnap false;

//- Use castellatedMeshControls::features (default = true)
explicitFeatureSnap true;

//- Detect points on multiple surfaces (only for explicitFeatureSnap)
multiRegionFeatureSnap false;
}

The last stage is adding layer , this boundary layer will be applied on the patch "aircraft"

 layers
{
"aircraft"
{
nSurfaceLayers 1;
}
} 

We are now from one step of mesh generation which is running snappyHexMesh.

With the first over detailed geometry we had some problems, like snappyHexMesh crashing when attempting surface snapping.

But when we ran “admesh” (which is ​is a free software, it is distributed under the terms of the GNU (GPL) for processing triangulated solid meshes.) we had a report on problems that are responsible for this crashing , which are related to some very thin surfaces. Those surfaces are the ones including small details.

admesh could be installed quickly and easily by an aptitude (sudo apt-get install admesh).

and run from terminal using the syntax : admesh [stl_file] [options] , or downloaded from the website .

we run admesh with the following options (see documentation):

admesh --exact --remove-unconnected --fill-holes--normal-directions --normal-values aircraft.stl

and we got the output :

================= Results produced by ADMesh version 0.95 ================
Input file         : aircraft.stl
File type          : ASCII STL file
Header             : solid aircraft
============== Size ==============
Min X = -170.016006, Max X =  176.872009
Min Y = -28.550001, Max Y =  74.518005
Min Z = -348.476013, Max Z =  55.994003
========= Facet Status ========== Original ============ Final ====
Number of facets                 : 90002               90002
Facets with 1 disconnected edge  :     0                   0
Facets with 2 disconnected edges :   234                   0
Facets with 3 disconnected edges :    46                   0
Total disconnected facets        :     0                   0
=== Processing Statistics ===     ===== Other Statistics =====
Number of parts       :    49        Volume   :  379064.187500
Degenerate facets     :     0
Edges fixed           :   280
Facets removed        :     0
Facets added          :     0
Facets reversed       :     0
Backwards edges       :     0
Normals fixed         :     0

We have then modified our geometry by removing all the odd parts and small details, and with the same procedure as below we have got a satisfactory result:

OpenFOAM comes with many mesh utilities like checkMesh which is a utility that checks the validity of a mesh. We can run it simply by typing in a terminal window : checkMesh , where the present working directory is the the one of our case.

Once we ran it, the output wad stating that the mesh was Ok! :

Mesh stats
points:           274946
faces:            609482
internal faces:   516433
cells:            169482
faces per cell:   6.64327
boundary patches: 6
point zones:      0
face zones:       0
cell zones:       0

Overall number of cells of each type:
hexahedra:     120280
prisms:        5386
wedges:        0
pyramids:      0
tet wedges:    188
tetrahedra:    9
polyhedra:     43619
Breakdown of polyhedra by number of faces:
faces   number of cells
4   925
5   4177
6   12180
7   37
8   160
9   15971
10   3
11   25
12   6886
15   2733
18   442
21   78
24   2

Checking topology...
Boundary definition OK.
Cell to face addressing OK.
Point usage OK.
Upper triangular ordering OK.
Face vertices OK.
Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces...
Patch    Faces   Points                  Surface topology
frontAndBack      128      162  ok (non-closed singly connected)
inlet      166      198  ok (non-closed singly connected)
outlet      160      189  ok (non-closed singly connected)
lowerWall      160      189  ok (non-closed singly connected)
upperWall      160      189  ok (non-closed singly connected)
aircraft_aircraft    92275    97199  multiply connected (shared edge)
<<Writing 30 conflicting points to set nonManifoldPoints

Checking geometry...
Overall domain bounding box (-300 -200 -3000) (300 200 1000)
Mesh (non-empty, non-wedge) directions (1 1 1)
Mesh (non-empty) directions (1 1 1)
Boundary openness (6.67487e-19 4.77163e-18 5.47761e-20) OK.
Max cell openness = 3.08091e-16 OK.
Max aspect ratio = 44.7833 OK.
Minimum face area = 0.028553. Maximum face area = 25000.  Face area magnitudes OK.
Min volume = 0.206848. Max volume = 750000.  Total volume = 9.59533e+08.  Cell volumes OK.
Mesh non-orthogonality Max: 79.8456 average: 31.4863
*Number of severely non-orthogonal faces: 85258.
Non-orthogonality check OK.
<<Writing 85258 non-orthogonal faces to set nonOrthoFaces
Face pyramids OK.
***Max skewness = 12.7514, 1367 highly skew faces detected which may impair the quality of the results
<<Writing 1367 skew faces to set skewFaces
Coupled point location match (average 0) OK.

Mesh OK.

End

# Initial and boundary conditions

Initial and boundary conditions

Boundary Conditions :

Our Boundary Conditions are as follows :

- Front , Back and upperwalls are of type slip.

- Inlet of type Fixed Inlet Velocity value.

- Outlet of type ZeroGradient.

- Lower wall of type Wall.

- On the patch corresponding to the geometry (aircraft) we set a wall boundary condition.

Initial conditions:

-Internal velocity uniform field (95, 2 0 0)

-Pressure 46040 Pa.

# Run time post-Processing

Setting run-time post processing functionalities :

In system we create 5 new files, which will control the post processing:

- readFields : In which we specify the fields we want to export for postprocessing.

// ************************************************************************* //

readFields

{functionObjectLibs ("libfieldFunctionObjects.so");

type readFields;

fields (p U k);

}

// ************************************************************************* //

- streamLines: It generates stream lines of the fields mentioned on readFields on several formats , here we chose a vtk outpout format supported by paraview.

// ************************************************************************* //

streamLines

{

type streamLine;

outputControl outputTime;

setFormat vtk; //gnuplot; //xmgr; //raw; //jplot; //csv; //ensight;

UName U;

trackForward true;

fields (p U k);

lifeTime 10000;

nSubCycle 5;

cloudName particleTracks;

seedSampleSet uniform; //cloud;//triSurfaceMeshPointSet;

uniformCoeffs

{

type uniform;

axis x; //distance;

start (-150 -200 40);

end (-15 -170 -100);

nPoints 20;

}

}

// ************************************************************************* //

- cuttingPlane : it controls the surface cut with one of the fields data .

cuttingPlane

{

type surfaces;

functionObjectLibs ("libsampling.so");

outputControl outputTime;

surfaceFormat vtk;

fields ( p U );

interpolationScheme cellPoint;

surfaces

(

yNormal

{

type cuttingPlane;

planeType pointAndNormal;

pointAndNormalDict

{

basePoint (-44.64 -198.82 65.76);

normalVector (-0.03 0.96 -0.25);

}

interpolate true;

}

);

}

// ************************************************************************* //

- forceCoeffs : calculates and generates bin gnuplot files .

// ************************************************************************* //

forceCoeffs1

{

type forceCoeffs;

functionObjectLibs ( "libforces.so" );

outputControl outputTime;

log yes;

patches ( "aircraft_aircraft" );

pName p;

UName U;

rhoName rhoInf; // Indicates incompressible

log true;

rhoInf 1; // Redundant for incompressible

liftDir (0 0 1);

dragDir (1 0 0);

CofR (0.72 0 0); // Axle midpoint on ground

pitchAxis (0 1 0);

magUInf 95.9;

lRef 1.42; // Wheelbase length

Aref 0.75; // Estimated

binData

{

nBin 90; // output data into 20 bins

direction (1 0 0); // bin direction

format gnuplot;

cumulative yes;

}

}

// ************************************************************************* //

- wallBoundedStreamLines : Does the same as streamLines , but on the geometry patch .

near

{

functionObjectLibs ("libfieldFunctionObjects.so");

type nearWallFields;

outputControl outputTime;

fields

(

(U UNear)

);

patches (aircraft);

distance 0.001;

}

wallBoundedStreamLines

{

functionObjectLibs ("libfieldFunctionObjects.so");

type wallBoundedStreamLine;

outputControl outputTime;

setFormat vtk;

UName UNear;

trackForward true;

interpolationScheme cellPoint;

fields (p U k UNear);

lifeTime 100;

cloudName wallBoundedParticleTracks;

seedSampleSet patchSeed;

uniformCoeffs

{

type uniform;

axis x; //distance;

start (0.0035 0.0999 0.0001);

end (0.0035 0.0999 0.0099);

nPoints 20;

}

cloudCoeffs

{

type cloud;

axis x; //distance;

points ((0.351516548679288 -0.0116085375585099 1.24));

}

patchSeedCoeffs

{

type patchSeed;

patches (aircraft);

axis x; //distance;

maxPoints 20000;

}

}

We could do some monitoring to watch residuals , continuity in run time , with the pyFoamPlotWatcher.py Utility (of pyFoam which is ​A python library to control OpenFOAM-runs and manipulate OpenFOAM-data. Comes with a number of utilities that should make your life easier  )to which we specify a custom regular expression to plot the forceCoeffs in runtime . To do so , we create a file named customRegexp in the case directory in which we put the following lines :

Force_Coeffcients

{

expr "Cl = (%f%)";

name Force_Coeffcients;

theTitle "Force Coefficients over simulation time";

titles (

"lift coeffcient"

);

type regular;

grid on;

autoscale;

}

Cd

{

expr "Cd = (%f%)";

name Force_Coeffcients_1;

titles (

drag coefficient

);

type slave;

master Force_Coeffcients;

}

Running the simulation

We run the simulation in parallel (see Allrun in the downloadable version of this case) and we generate log files from which pyFoamPlotWatcher.py that we run in the same time parses values and plot them on gnuplot as follows :

The whole case can be downloaded here (download the 3 parts, and then change the extension towards .rar instead of .zip): part 1 , part 2 , part 3

# Results

## Results

The image above is an animation showing the evolution of streamlines, that we can follow while the calculation is running (it is much faster here).

As the interest of this project is not really the study of the aerodynamics of the F-14, the results are not really relevant; the stress was only put on the fact that they were coherent.

Here are some of the results we can display using paraview (with the aid of the utility paraFoam):

• Velocity magnitude, in a plane cutting the domain:

• Streamlines, colored with a scale of pressure:

• Below we can see streamlines along the surface of the plane (WallboundedLines utility), which gives the path taken by the air when encountering the plane, colored by the value of pressure. It shows where are the most resistive parts of the aircraft.

We can note that there are high values of pressure inside the throats leading to the entry of the turbocompressors, which is not occuring in reality, since this is a study of the aircraft when passive in the airflow (engines off).

• Pressure in a  plane perpendicular to a wing:

Here we can see that the difference in pressure between the upper wing and the lower wing (as well as for the upper part of the empennage and its lower part) which creates the lift, it seems coherent according to the lift coefficient ploted during the calculation. Nevertheless, this kind of modelisation may not be reffined enough to study the airfoil; if we wished to do that we could reffine further the mesh near the wing with snappyHexMesh, however the time of calculation would become even more huge. Besides, the mesh is also to gross to be able to do so.

Which is interesting here is that the post-processing part can be begun before the end of the calculation (actually those results where made that way, and do not show the final time of the computation) and allow to stop it if we see incoherent results.

# Conclusion

Conclusion

In this project we have experienced the interest of using run-time post-treatment, indeed the critical observation of the output of the calculation have allowed us to save time by stopping the calculation several times, as we faced many difficulties before having a proper simulation.
As the residuals and the continuity were satisfactory, the root of the error was not likely to be numerical, but rather from the setup of the case, it was supported by the fact that the lift and drag coefficient had unlikely values. All of that plus the vizualisation of  the odd aspect of the flow, lead us to question the mesh, which happened to be too gross to have proper results.

It is also interesting to predict the time needed by the simulation to end, and optimize the time of actual work.

As a conclusion, we can say that the run-time post-processing is almost necessary to anyone who whishes not to loose huge amounts of time by running false calculations, and by not knowing the exact time when it stops.
To go further, we could think of making a graphic interface with pyFoam, which would allow to have a more practical tool, and would be usable and integrable easily in OpenFoam cases, even by novices.

The whole project can be downloaded on this page.