How to use cfMesh ? A first tutorial based on the Ahmed body

cfMesh is a mesher developed by Franjo Juretić (https://cfmesh.com/) and incorporated in OpenFOAM® (ESI-group version) since v1712. The version available within OpenFOAM® package is v1.1. The source code of cfMesh as well as the user guide can by found inside the installation directory of OpenFOAM® in modules/cfmesh. cfMesh is a very powerful mesher able to generate millions cells in a very short time.

The Ahmed body [1] is a simplified version of a vehicle geometry. Its shape is simple enough to allow for accurate flow simulation but retains some important practical features relevant to automobile bodies. This tutorial aims to described a simple and easy cfMesh workflow for external flows (aerodynamics, hydrodynamics etc..). Of course other workflows are possible but this one is intended to be as simple as possible.

References :

[1] : https://www.cfd-online.com/Wiki/Ahmed_body

[2] : http://cfmesh.com/wp-content/uploads/2015/09/User_Guide-cfMesh_v1.1.pdf

[3] : https://grabcad.com/library/ahmed-body-4

The details of the geometry can be found in the reference [1]. The CAD model can be downloaded from grabcad [3] thanks to Colby Mazzuca.

The geometry has to be exported to stl format to be used a input data of cfMesh. Generally all CAD software allow to export a geometry as stl. Nevertheless it is not always possible to control the density of triangles on the surface, and cfMesh (as well as snappyHexMesh) do prefers a high density of triangles. If you want to improve the surface quality of your stl geometry you can use for example SALOME and the MEFISTO surface mesher. But don’t worry for a first time if you can’t control the triangle density it will also work !

As described in [2] cfMesh prefers an fms format for the geometry. Thus we are going to use the command surfaceFeatureEdges to convert our stl geometry to a fms geometry. But before proceeding to the conversion we need to generate a bounding box to define our inlet/outlet/top/bottom and side boundaries. cfMesh has an useful utility designed for that : surfaceGenerateBoundingBox. In the case directory you can type the following command to generate a new stl surface named ahmedBodyBox.stl containing the bounding box patches:

[pastacode lang=”cpp” manual=”surfaceGenerateBoundingBox%20ahmedBody.stl%20ahmedBodyBox.stl%203.0%204.0%200%202%202%202″ message=”bounding box generation” highlight=”” provider=”manual”/]

The last 6 numbers stand for the size of the box in each direction as indicated below with the surfaceGenerateBoundingBox -help command.

[pastacode lang=”cpp” manual=”surfaceGenerateBoundingBox%20-help%0A%0AUsage%3A%20surfaceGenerateBoundingBox%20%5BOPTIONS%5D%20%3Cinput%20surface%20file%3E%20%3Coutput%20surface%20file%3E%20%3Cx-neg%3E%20%3Cx-pos%3E%20%3Cy-neg%3E%20%3Cy-pos%3E%20%3Cz-neg%3E%20%3Cz-pos%3E” message=”surfaceGenerateBoundingBox -help” highlight=”” provider=”manual”/]

At this point by inspecting the ahmedBodyBox.stl we can see that 6 new patches have been created and added at the end of the file. The 6 new patches are automatically naturally named :

• xMax and xMin
• yMax and yMin
• zMax and zMin

Once the ahmedBodyBox.stl geometry has been generated the conversion to fms format is done by typing the following command  :

[pastacode lang=”cpp” manual=”surfaceFeatureEdges%20-angle%2015%20ahmedBodyBox.stl%20ahmedBody.fms” message=”fms conversion” highlight=”” provider=”manual”/]

For the features capture we use the -angle option with an angle value of 15°. Note that the default angle for features is 45°

The dictionary for mesh (called meshDict) instruction has to be located in the /system directory. cfMesh required minimal user parameters to work.

First the mesh dictionary contains the name of the surface. We have put the same max/min and boundary size to obtain a uniform far field hexaedral mesh.

[pastacode lang=”cpp” manual=”surfaceFile%20%22ahmedBody.fms%22%3B%0A%0AmaxCellSize%200.25%3B%20%2F%2F%5Bm%5D%0AminCellSize%200.25%3B%20%2F%2F%5Bm%5D%0AboundaryCellSize%200.25%3B%20%2F%2F%5Bm%5D” message=”Basic parameters” highlight=”” provider=”manual”/]

Now we want to refine the cells close to the Ahmed body especially in the wake region. Thus, we are going to use a box refinement :

[pastacode lang=”cpp” manual=”objectRefinements%0A%7B%0A%20%20%20%20mainBox%0A%20%20%20%20%7B%0A%20%20%20%20%20%20%20%20type%20box%3B%0A%20%20%20%20%20%20%20%20cellSize%20%20%200.05%3B%20%2F%2F%5Bm%5D%0A%20%20%20%20%20%20%20%20centre%20%20%20%20(0%200%200)%3B%20%2F%2F%5Bm%5D%0A%20%20%20%20%20%20%20%20lengthX%20%20%20%20%20%203%3B%20%2F%2F%5Bm%5D%0A%20%20%20%20%20%20%20%20lengthY%20%20%20%201.0%3B%20%2F%2F%5Bm%5D%0A%20%20%20%20%20%20%20%20lengthZ%20%20%20%20%20%201.5%3B%20%2F%2F%5Bm%5D%0A%20%20%20%20%7D%0A%7D” message=”Wake zone refinement” highlight=”” provider=”manual”/]

Then we would like to refine the cells in the vicinity of the Ahmed body. In the spirit of snappyHexMesh we use the additionalRefinementLevels keyword. A value of 5 means that cells are cut 5 times (in each Cartesian direction) relatively to the far field cells. We also use the keyword refinementThickness in order to define the thickness of the refinement zone.

[pastacode lang=”cpp” manual=”localRefinement%0A%7B%0A%0A%20%20%20%20ahmedBody%0A%20%20%20%20%7B%0A%20%20%20%20%20%20%20%20additionalRefinementLevels%205%3B%0A%20%20%20%20%20%20%20%20refinementThickness%20%200.05%3B%20%2F%2F%5Bm%5D%0A%20%20%20%20%7D%0A%7D” message=”refinement of the ahmed body” highlight=”” provider=”manual”/]

In order to inflate boundary layer it’s necessary to define a special subdict detailed below. Here 5 layers are inflated with a ratio of 1.2. We also use the keyword optimiseLayer to improve the quality of the layers. The allowDiscontinuity 1 keyword ensures that the layer are not spread over neighboring patches. Note that the default value is 0. So you will often need to play with this parameter.

[pastacode lang=”cpp” manual=”boundaryLayers%0A%7B%0A%20%20%20%20patchBoundaryLayers%0A%20%20%20%20%7B%0A%20%20%20%20%20%20%20%20ahmedBody%0A%20%20%20%20%20%20%20%20%7B%0A%20%20%20%20%20%20%20%20%20%20%20%20%20nLayers%20%20%20%205%3B%0A%20%20%20%20%20%20%20%20%20%20%20%20%20thicknessRatio%20%20%20%20%201.2%3B%0A%20%20%20%20%20%20%20%20%20%20%20%20%20allowDiscontinuity%201%3B%0A%20%20%20%20%20%20%20%20%7D%0A%20%20%20%20%7D%0A%0A%20%20%20%20optimiseLayer%201%3B%0A%0A%20%20%20%20optimisationParameters%0A%20%20%20%20%7B%0A%20%20%20%20%20%20%20%20nSmoothNormals%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%205%3B%0A%20%20%20%20%20%20%20%20relThicknessTol%20%09%09%20%20%20%20%20%20%20%200.4%3B%0A%20%20%20%20%20%20%20%20featureSizeFactor%20%09%090.8%3B%0A%20%20%20%20%20%20%20%20reCalculateNormals%20%09%091%3B%20%0A%20%20%20%20%20%20%20%20maxNumIterations%20%09%095%3B%0A%20%20%20%20%7D%0A%0A%7D” message=”Layer inflation” highlight=”” provider=”manual”/]

Finally, because we have used surfaceGenerateBoundingBox it’s useful to rename the patches xMax/xMin/yMax/yMin/zMax/Zmin in a more conventional manner.

[pastacode lang=”cpp” manual=”renameBoundary%0A%7B%0A%20%20%20%20defaultType%20%20%20%20%20wall%3B%0A%0A%20%20%20%20newPatchNames%0A%20%20%20%20%7B%0A%20%20%20%20%20%20%20%20%22xMax%22%20%7B%20newName%20outlet%20%20%3B%20type%20%20patch%3B%20%7D%0A%20%20%20%20%20%20%20%20%22xMin%22%20%7B%20newName%20inlet%20%3B%20type%20%20patch%3B%20%7D%0A%20%20%20%20%20%20%20%20%22yMax%22%20%7B%20newName%20top%20%20%20%20%3B%20type%20%20patch%3B%20%7D%0A%20%20%20%20%20%20%20%20%22yMin%22%20%7B%20newName%20bottom%20%3B%20type%20%20patch%3B%20%7D%0A%20%20%20%20%20%20%20%20%22zMax%22%20%7B%20newName%20side%20%20%20%3B%20type%20%20patch%3B%20%7D%0A%20%20%20%20%20%20%20%20%22zMin%22%20%7B%20newName%20side%20%20%20%3B%20type%20%20patch%3B%20%7D%0A%20%20%20%20%7D%0A%7D” message=”rename of boundary” highlight=”” provider=”manual”/]

For generating the mesh you only need to type in the terminal :

[pastacode lang=”cpp” manual=”cartesianMesh” message=”” highlight=”” provider=”manual”/]

Note that cartesianMesh is automatically parallelized. Don’t forget to run checkMesh to verify that the mesh quality is compatible with the solver requirements. In the present case we obtain the following checkMesh output (number of cells ~ 576 k):

[pastacode lang=”cpp” manual=”Create%20time%0A%0ACreate%20mesh%20for%20time%20%3D%200%0A%0ATime%20%3D%200%0A%0AMesh%20stats%0A%20%20%20%20points%3A%20%20%20%20%20%20%20%20%20%20%20620703%0A%20%20%20%20faces%3A%20%20%20%20%20%20%20%20%20%20%20%201772737%0A%20%20%20%20internal%20faces%3A%20%20%201724583%0A%20%20%20%20cells%3A%20%20%20%20%20%20%20%20%20%20%20%20576180%0A%20%20%20%20faces%20per%20cell%3A%20%20%206.069839%0A%20%20%20%20boundary%20patches%3A%206%0A%20%20%20%20point%20zones%3A%20%20%20%20%20%200%0A%20%20%20%20face%20zones%3A%20%20%20%20%20%20%200%0A%20%20%20%20cell%20zones%3A%20%20%20%20%20%20%200%0A%0AOverall%20number%20of%20cells%20of%20each%20type%3A%0A%20%20%20%20hexahedra%3A%20%20%20%20%20557018%0A%20%20%20%20prisms%3A%20%20%20%20%20%20%20%20484%0A%20%20%20%20wedges%3A%20%20%20%20%20%20%20%200%0A%20%20%20%20pyramids%3A%20%20%20%20%20%20890%0A%20%20%20%20tet%20wedges%3A%20%20%20%200%0A%20%20%20%20tetrahedra%3A%20%20%20%201348%0A%20%20%20%20polyhedra%3A%20%20%20%20%2016440%0A%20%20%20%20Breakdown%20of%20polyhedra%20by%20number%20of%20faces%3A%0A%20%20%20%20%20%20%20%20faces%20%20%20number%20of%20cells%0A%20%20%20%20%20%20%20%20%20%20%20%206%20%20%201874%0A%20%20%20%20%20%20%20%20%20%20%20%207%20%20%201280%0A%20%20%20%20%20%20%20%20%20%20%20%208%20%20%204%0A%20%20%20%20%20%20%20%20%20%20%20%209%20%20%2012212%0A%20%20%20%20%20%20%20%20%20%20%2010%20%20%2028%0A%20%20%20%20%20%20%20%20%20%20%2012%20%20%201034%0A%20%20%20%20%20%20%20%20%20%20%2014%20%20%202%0A%20%20%20%20%20%20%20%20%20%20%2015%20%20%206%0A%0AChecking%20topology…%0A%20%20%20%20Boundary%20definition%20OK.%0A%20%20%20%20Cell%20to%20face%20addressing%20OK.%0A%20%20%20%20Point%20usage%20OK.%0A%20%20%20%20Upper%20triangular%20ordering%20OK.%0A%20%20%20%20Face%20vertices%20OK.%0A%20%20%20%20Number%20of%20regions%3A%201%20(OK).%0A%0AChecking%20patch%20topology%20for%20multiply%20connected%20surfaces…%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20Patch%20%20%20%20Faces%20%20%20Points%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20Surface%20topology%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20ahmedBody%20%20%20%2027612%20%20%20%2027604%20%20ok%20(non-closed%20singly%20connected)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20outlet%20%20%20%20%20%20240%20%20%20%20%20%20273%20%20ok%20(non-closed%20singly%20connected)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20inlet%20%20%20%20%20%20240%20%20%20%20%20%20273%20%20ok%20(non-closed%20singly%20connected)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20top%20%20%20%20%20%20720%20%20%20%20%20%20777%20%20ok%20(non-closed%20singly%20connected)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20bottom%20%20%20%2018478%20%20%20%2018885%20%20ok%20(non-closed%20singly%20connected)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20side%20%20%20%20%20%20864%20%20%20%20%20%20962%20%20ok%20(non-closed%20singly%20connected)%0A%0AChecking%20faceZone%20topology%20for%20multiply%20connected%20surfaces…%0A%20%20%20%20No%20faceZones%20found.%0A%0AChecking%20basic%20cellZone%20addressing…%0A%20%20%20%20No%20cellZones%20found.%0A%0AChecking%20geometry…%0A%20%20%20%20Overall%20domain%20bounding%20box%20(-4.044%20-1.005916e-06%20-2.1945)%20(4%202.338%202.1945)%0A%20%20%20%20Mesh%20has%203%20geometric%20(non-empty%2Fwedge)%20directions%20(1%201%201)%0A%20%20%20%20Mesh%20has%203%20solution%20(non-empty)%20directions%20(1%201%201)%0A%20%20%20%20Boundary%20openness%20(8.025553e-17%20-4.900931e-16%20-1.255092e-16)%20OK.%0A%20%20%20%20Max%20cell%20openness%20%3D%205.841175e-16%20OK.%0A%20%20%20%20Max%20aspect%20ratio%20%3D%2030.28225%20OK.%0A%20%20%20%20Minimum%20face%20area%20%3D%201.133135e-08.%20Maximum%20face%20area%20%3D%200.06809194.%20%20Face%20area%20magnitudes%20OK.%0A%20%20%20%20Min%20volume%20%3D%208.964937e-10.%20Max%20volume%20%3D%200.01687817.%20%20Total%20volume%20%3D%2082.43248.%20%20Cell%20volumes%20OK.%0A%20%20%20%20Mesh%20non-orthogonality%20Max%3A%2080.2906%20average%3A%204.206543%0A%20%20%20*Number%20of%20severely%20non-orthogonal%20(%3E%2070%20degrees)%20faces%3A%2013.%0A%20%20%20%20Non-orthogonality%20check%20OK.%0A%20%20%3C%3CWriting%2013%20non-orthogonal%20faces%20to%20set%20nonOrthoFaces%0A%20%20%20%20Face%20pyramids%20OK.%0A%20%20%20%20Max%20skewness%20%3D%203.817519%20OK.%0A%20%20%20%20Coupled%20point%20location%20match%20(average%200)%20OK.%0A%0AMesh%20OK.%0A%0AEnd%0A” message=”” highlight=”” provider=”manual”/]

The resolution of the numerical model is done with the simpleFoam solver of OpenFOAM. simpleFoam is a special solver for steady-state incompressible laminar/turbulent flows. The air velocity is 40 m/s. The boundary conditions are very classic :

• inlet : fixedValue for U and turbulence. zeroGradient for p.
• outlet : zeroGradient for U and turbulence. fixedValue for p.
• ahmedBody : noSlip for U, zeroGradient for p. Wall laws for turbulence.
• top/bottom/side : zeroGradient for all fields.

For inlet turbulence parameters settings it’s generally advised to use a calculator. Inside the controlDict file it’s often useful to add runtime post-processing function. For example, in this particular case we can obtain forces, y+, wall shear stress and also residuals. During runtime all this results are written in a special folder named postProcessing in the case directory.

[pastacode lang=”cpp” manual=”functions%0A%7B%0A%0A%20%20%20%20forces%0A%20%20%20%20%7B%0A%20%20%20%20%20%20%20%20type%20%20%20%20%20%20%20%20%20%20%20%20forces%3B%0A%20%20%20%20%20%20%20%20libs%20%20%20%20%20%20%20%20%20%20%20%20(%22libforces.so%22)%3B%0A%20%20%20%20%20%20%20%20patches%20%20%20%20%20%20%20%20%20(ahmedBody)%3B%0A%20%20%20%20%20%20%20%20rho%20%20%20%20%20%20%20%20%20%20%20%20%20rhoInf%3B%0A%20%20%20%20%20%20%20%20rhoInf%20%20%20%20%20%20%20%20%20%201%3B%0A%20%20%20%20%20%20%20%20log%20%20%20%20%20%20%20%20%20%20%20%20%20on%3B%0A%20%20%20%20%20%20%20%20writeControl%20%20%20%20timeStep%3B%0A%20%20%20%20%20%20%20%20writeInterval%20%20%201%3B%0A%20%20%20%20%20%20%20%20CofR%20%20%20%20%20%20%20%20%20%20%20%20(0%200%200)%3B%0A%20%20%20%20%7D%0A%20%20%20%20wallShearStress%0A%20%20%20%20%7B%0A%20%20%20%20%20%20%20%20type%20%20%20%20%20%20%20%20%20%20%20%20wallShearStress%3B%0A%20%20%20%20%20%20%20%20libs%20%20%20%20%20%20%20%20%20%20%20%20(%22libfieldFunctionObjects.so%22)%3B%0A%20%20%20%20%20%20%20%20writeControl%20%20%20%20writeTime%3B%0A%20%20%20%20%20%20%20%20patches%20%20%20%20%20%20%20%20%20(ahmedBody)%3B%0A%20%20%20%20%7D%0A%0A%20%20%20%20yPlus%0A%20%20%20%20%7B%0A%20%20%20%20%20%20%20%20type%20%20%20%20%20%20%20%20%20%20%20%20yPlus%3B%0A%20%20%20%20%20%20%20%20libs%20%20%20%20%20%20%20%20%20%20%20%20(%22libfieldFunctionObjects.so%22)%3B%0A%20%20%20%20%20%20%20%20writeControl%20%20%20%20writeTime%3B%0A%20%20%20%20%20%20%20%20patches%20%20%20%20%20%20%20%20%20(ahmedBody)%3B%0A%20%20%20%20%7D%0A%0A%20%20%20%20residuals%0A%20%20%20%20%7B%0A%20%20%20%20%20%20%20%20type%20%20%20%20%20%20%20%20%20%20%20%20residuals%3B%0A%20%20%20%20%20%20%20%20libs%20%20%20%20%20%20%20%20%20%20%20%20(%22libutilityFunctionObjects.so%22)%3B%0A%20%20%20%20%20%20%20%20fields%20%20%20%20%20%20%20%20%20%20(%22.*%22)%3B%0A%20%20%20%20%7D%0A%7D” message=”” highlight=”” provider=”manual”/]