Why is OpenFOAM comparatively more challenging?
While OpenFOAM directory files are primarily written in a plain-text configuration, it's design closely resembles C++. Interestingly, OpenFOAM source code is written entirely in C++. Additionally, the software must be operated using Linux or WSL (if in Windows). This are some of the reasons why it is extremely difficult for beginners to get a grasp of this software.
Reality & pre-requisites
To genuinely use OpenFOAM for running fluid simulations, it is mandatory to get an idea of fluid dynamics from a theoretical perspective. Without it, you will run into a heap of errors and eventually quit out of frustration. You must also get yourself more or less comfortable with basic Linux commands (such as mkdir, cd, cp, rm and their respective extensions) & the nano text editor. Personally, I used Ubuntu through WSL (Windows Sub-System for Linux). This post particularly focuses on external fluid steady-state simulations. The directory will be notably different from internal flow or transient flow simulations.
What to simulate?
For the purpose of my project, which inspired this post, I have simulated three types of wings. This post will contain the baseline configuration rectangular wing case files (To learn more about the project itself, visit the public GitHub repository here). The files are accommodated accordingly, but with a few changes, any solid body can be meshed & simulated. Ensure the solid is exported as a .stl file, as that is best supported by snappyHexMesh.
Folders in an OpenFOAM directory
An external flow simulation directory would look something like this:
caseFiles
.
├── 0/
│ ├── U
│ ├── p
│ ├── k
│ ├── nut
│ └── omega
├── constant/
│ ├── triSurface/
│ │ └── *.stl
│ ├── transportProperties
│ └── turbulenceProperties
└── system/
├── blockMeshDict
├── controlDict
├── snappyHexMeshDict
├── fvSchemes
├── fvSolution
└── surfaceFeatureExtract
For a lot of beginners, this may seem insanely confusing. The complexity most definitely dawned on me when I was first getting started. This may be a good time to say that I myself am not particularly very experienced, but I have ran a few simulations & majorly focused them on external aerodynamics, hence this post. OpenFOAM has built in tutorials which have pre-configured cases & can be of good help. Scroll down to learn how to access those.
The "0" folder
This is the folder which includes the case boundary conditions. Every single file in this folder needs mentioning of every single part of the mesh along with their properties to define them.
In our case, U refers to the velocity conditions. We have a box with a wing placed inside it. The "box" acts as a wind tunnel, and the wing is snapped in so fluid can flow over it. Here's the file contents:
/*--------------------------------*- C++ -*----------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (29.886 0 2.616);
boundaryField
{
inlet
{
type fixedValue;
value uniform (29.886 0 2.616);
}
outlet
{
type inletOutlet;
inletValue uniform (0 0 0);
value uniform (29.886 0 2.616);
}
top
{
type slip;
}
bottom
{
type slip;
}
sideMinusY
{
type slip;
}
sidePlusY
{
type symmetryPlane;
}
wing
{
type noSlip;
}
}
// ************************************************************************* //
For further context, here I have simulated a 5 degree angle of attack, and so I rotated the velocity by using vectorization. In most cases here, all vector notations will be in the form (x, y, z). As for the boundary conditions, slip & noSlip indicates whether a surface has friction or not. fixedValue forces the specific velocity vector at the inlet, and inletOutlet is a hybrid condition which acts as a zeroGradient condition if fluid leaves the domain, but switches to fixedValue if fluid attempts to flow back, preventing backflow & unphysical value results from the simulation. There are many more boundary conditions that is not within the scope of this post, but explanations are widely available for each along with their use cases.
The other files have a similar structure, but omega, nut & k requires a little bit of calculation. omega refers to turbulent specific dissipation rate, nut refers to turbulent kinematic viscosity, & k refers to turbulent kinetic energy. Each of these has wall functions that must be decided upon to place on the mesh. There are plenty information available for each wall functions, and their specific use cases. For this case, the following was used:
omegaWallFunctionnutUSpaldingWallFunctionkqRWallFunction
Some conditions must be decided upon based on your simulation case. For this, we will use the following:
Air Density =
Dynamic Viscosity =
Freestream Velocity =
Now with these conditions decided, we can use the following formulas for each of conditions:
3 % turbulence was simulated
0.09 is the (Turbulent viscosity constant). 0.07 x Chord length is calculating (turbulent length scale)
The "constant" folder
In most OpenFOAM case directories, the constant/ folder contains properties that remain constant (unchanged) throughout the entire simulation. The polyMesh/ subfolder contains the *.stl file, which is then used by snappyHexMesh to snap the shape into the domain. There are two more important configuration files in this directory:
1. turbulenceProperties
/*--------------------------------*- C++ -*----------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object turbulenceProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
simulationType RAS;
RAS
{
RASModel kOmegaSST;
turbulence on;
printCoeffs on;
}
// ************************************************************************* //
Under simulationType, you must enter the turbulence model to be used. Some popular models are laminar, RAS & LES. For this case, RAS was the most appropriate, and so RASModel was later specified with a model-specific coefficient (kOmegaSST in this case).
2. transportProperties
/*--------------------------------*- C++ -*----------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
transportModel Newtonian;
nu nu [0 2 -1 0 0 0 0] 1.4607e-05;
// ************************************************************************* //
This file defines the fluid material properties. In this case, it defines our nu( , Kinematic viscosity) as per the previous calculation. It also includes the option to specify the transport model to be used, Newtonian or Non-Newtonian.
The "system" folder
The system/ folder practically defines how a mesh would be made and how the solver will solve the physics defined in 0/ & constant/. It includes a number of important files.
1. controlDict (Simulation Control Directory)
/*--------------------------------*- C++ -*----------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application simpleFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 500;
deltaT 1;
writeControl timeStep;
writeInterval 50;
purgeWrite 2;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
functions
{
forceCoeffs1
{
type forceCoeffs;
libs ("libforces.so");
writeControl timeStep;
writeInterval 1;
patches (wing);
// Force coefficients are reported relative to the FREESTREAM
// direction (5 deg), not the body/chord axes. Our geometry is aligned with the body axes (0 deg), so we need
// to rotate the freestream vector.
// dragDir = freestream unit vector = (cos5, 0, sin5)
// liftDir = perpendicular, rotated +90 deg = (-sin5, 0, cos5)
dragDir (0.9962 0 0.0872);
liftDir (-0.0872 0 0.9962);
CofR (0.05 -0.3 0);
pitchAxis (0 1 0);
rho rhoInf;
rhoInf 1.225;
magUInf 30;
lRef 0.2; // chord length
Aref 0.12; // planform area = chord x span = 0.2 x 0.6 [m^2]
binData //Lift/drag distributed along the span
{
nBin 20;
direction (0 1 0);
cumulative no;
}
}
{
type yPlus;
libs ("libfieldFunctionObjects.so");
writeControl writeTime;
patches (wing);
}
CpField
{
type pressure;
libs ("libfieldFunctionObjects.so");
Cp yes;
Cref 0;
rho rhoInf;
rhoInf 1.225;
URef 30;
writeControl writeTime;
}
residuals1
{
type residuals;
libs ("libutilityFunctionObjects.so");
writeControl timeStep;
writeInterval 1;
fields (p U k omega);
}
}
// ************************************************************************* //
This directory includes control variables that governs the loop of the solver execution. It also controls at which iterations data is to be saved, the residual limit at which the solver can converge, and a forceCoeffs section which calculates the lift/drag forces and a few other forces for each iteration. This gives a numerical database of the solver run, and is very helpful when analyzing KPIs.
2. fvSchemes (Finite Volume Schemes)
/*--------------------------------*- C++ -*----------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
ddtSchemes
{
default steadyState;
}
gradSchemes
{
default Gauss linear;
grad(U) cellLimited Gauss linear 1;
}
divSchemes
{
default none;
div(phi,U) bounded Gauss linearUpwindV grad(U);
div(phi,k) bounded Gauss upwind;
div(phi,omega) bounded Gauss upwind;
div((nuEff*dev2(T(grad(U))))) Gauss linear;
}
laplacianSchemes
{
default Gauss linear corrected;
}
interpolationSchemes
{
default linear;
}
snGradSchemes
{
default corrected;
}
wallDist
{
method meshWave;
}
// ************************************************************************* //
This file defines the discretization schemes used in the governing equations. Since this case is a RANS (Reynold's Averaged Navier Stokes) simulation, and Navier Stokes equations are partial differentiation equations, this directory converts these equations into algebraic forms for the solver to interpret it and solve the equations.
3. fvSolution (Finite Volume Equation Solvers)
/*--------------------------------*- C++ -*----------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
solvers
{
p
{
solver GAMG;
tolerance 1e-06;
relTol 0.1;
smoother GaussSeidel;
}
"(U|k|omega)"
{
solver smoothSolver;
smoother GaussSeidel;
tolerance 1e-08;
relTol 0.1;
}
}
SIMPLE
{
nNonOrthogonalCorrectors 2;
residualControl
{
p 1e-4;
U 1e-4;
"(k|omega|nut)" 1e-4;
}
}
relaxationFactors
{
fields
{
p 0.3;
}
equations
{
U 0.7;
k 0.7;
omega 0.7;
}
}
// ************************************************************************* //
This directory controls which linear equation solver needs to be used for the simulation, the tolerances for the solver & the algorithms that are used to calculate the matrix equations. For this case, the GAMG solver is selected as the matrix solver. Since it is a steady-state simulation, the SIMPLE algorithm sub-directory is used, which essentially controls the stability of the run. For transient simulations, PISO or PIMPLE can be used.
4. blockMeshDict ("Wind Tunnel")
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM
\\ / O peration |
\\ / A nd |
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// All coordinates below are in millimetres. This factor converts them to
// metres, which is the unit OpenFOAM solvers assume (kg-m-s SI units).
convertToMeters 0.001;
// Domain corner coordinates (mm), derived from chord=200mm, span=600mm:
// x: -1000 (5c upstream) to 2200 (10c downstream after 200mm chord)
// y: -1200 (tip + 3c) to 0 (Root chord at symmetryPlane)
// z: -1000 (5c below) to 1000 (5c above)
vertices
(
(-1000 -1200 -1000) // 0
( 2200 -1200 -1000) // 1
( 2200 0 -1000) // 2
(-1000 0 -1000) // 3
(-1000 -1200 1000) // 4
( 2200 -1200 1000) // 5
( 2200 0 1000) // 6
(-1000 0 1000) // 7
);
// Background cell size ~50 mm in every direction.
// 3200/50=64, 1200/50=24, 2000/50=40
// snappyHexMesh refines locally near the wing, this is just the coarse base.
blocks
(
hex (0 1 2 3 4 5 6 7) (64 24 40) simpleGrading (1 1 1)
);
edges
(
);
// Patch names used later in 0/ boundary conditions and snappyHexMeshDict.
boundary
(
inlet
{
type patch;
faces
(
(0 4 7 3)
);
}
outlet
{
type patch;
faces
(
(1 2 6 5)
);
}
bottom
{
type patch;
faces
(
(0 3 2 1)
);
}
top
{
type patch;
faces
(
(4 5 6 7)
);
}
sideMinusY
{
type patch;
faces
(
(0 1 5 4)
);
}
sidePlusY
{
type symmetryPlane;
faces
(
(3 7 6 2)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
This is OpenFOAM's native block-structured mesh generator. It generates a square mesh, and you can configure it to act as the "Wind Tunnel" for an external flow simulation. It is also possible to CAD model it and use a *.stl file, but since OpenFOAM has a native built-in feature, you can choose which path to take based on your preference.
For aerodynamics cases, a general rule for spacing of the domain is to take the length of the object under simulation, then place 5 lengths upstream, 10 lengths downstream, 3 lengths on each side, & 5 lengths on top/bottom. The blockMeshDict text above is marked with notations to help understand the script, feel free to go through the script & try to gain a cognitive understanding of what each line in the script could mean. Do this with all the configuration files, and any doubts will likely be cleared away.
5. snappyHexMeshDict (Adding Solid to Domain)
/*--------------------------------*- C++ -*----------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object snappyHexMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
castellatedMesh true;
snap true;
addLayers true;
// Geometry — surfaces used for castellation and as refinement volumes
geometry
{
rectangular_600-200.stl
{
type triSurfaceMesh;
name wing;
scale 0.001;
}
// A box around the wing where we want extra background refinement
refinementBox
{
type searchableBox;
min (-0.10 -0.70 -0.20);
max ( 0.60 0.00 0.20);
}
teRefinementBox
{
type searchableBox;
min (0.20 -0.7 -0.05);
max (0.60 0.00 0.05);
}
};
// Castellated Mesh Controls
castellatedMeshControls
{
maxLocalCells 2000000;
maxGlobalCells 8000000;
minRefinementCells 10;
nCellsBetweenLevels 4;
features
(
{
file "rectangular_600-200.eMesh";
level 7;
}
);
refinementSurfaces
{
wing
{
level (6 7);
patchInfo
{
type wall;
}
}
}
resolveFeatureAngle 60;
refinementRegions
{
refinementBox
{
mode inside;
levels ((1e15 5));
}
teRefinementBox
{
mode inside;
levels ((1e15 7));
}
}
locationInMesh (1.5 -0.3 0.5);
allowFreeStandingZoneFaces true;
}
// Snap Controls
snapControls
{
nSmoothPatch 3;
tolerance 4.0;
nSolveIter 50;
nRelaxIter 8;
nFeatureSnapIter 10;
implicitFeatureSnap false;
explicitFeatureSnap true;
multiRegionFeatureSnap false;
}
// Boundary Layer Controls
addLayersControls
{
relativeSizes true;
layers
{
wing
{
nSurfaceLayers 5; // 5 layers per spec
}
}
expansionRatio 1.2; // 5 layers @ 1.20 growth rate
minThickness 5e-06;
firstLayerThickness 0.01041; // Used a y1 sizing script for a y+ target of 1
nGrow 1;
featureAngle 60;
nRelaxIter 20;
nSmoothSurfaceNormals 1;
nSmoothNormals 3;
nSmoothThickness 10;
nMedialAxisIter 10;
maxFaceThicknessRatio 0.5;
maxThicknessToMedialRatio 0.5;
minMedialAxisAngle 80;
nBufferCellsNoExtrude 4;
nLayerIter 40;
}
meshQualityControls
{
maxNonOrtho 65;
maxBoundarySkewness 20;
maxInternalSkewness 4;
maxConcave 80;
minFlatness 0.5;
minVol 1e-13;
minTetQuality 1e-30;
minArea -1;
minTwist 0.02;
minDeterminant 0.001;
minFaceWeight 0.05;
minVolRatio 0.01;
minTriangleTwist -1;
nSmoothScale 4;
errorReduction 0.75;
relaxed
{
maxNonOrtho 75;
}
}
mergeTolerance 1e-6;
//************************************************************************* //
This is perhaps the most important mesh configuration file for a CFD case in OpenFOAM. It defines background mesh snapping, refinement regions where the mesh needs to be denser to calculate accurate fluid flow, boundary layer addition settings & mesh quality controls to make sure enough meshing iterations are passed before meshing completes so that the mesh is not broken or give out inaccurate solver results. It is heavily recommended to go through this file on a line-by-line basis, as it is quite dense with information about the case & will give out a solid idea about how the meshing is being done by snappyHexMesh.
6. surfaceFeatureExtractDict (Sharp Edges Control)
/*--------------------------------*- C++ -*----------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object surfaceFeatureExtractDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
rectangular_600-200.stl
{
extractionMethod extractFromSurface;
extractFromSurfaceCoeffs
{
includedAngle 150;
}
{
writeObj yes;
}
// ************************************************************************* //
This is a simple yet essential configuration file that extracts sharp geometric edges & curves from the *.stl and produces an *.eMesh file at the constant/triSurface directory. Here in this case, the angle threshold for a sharp edge is 150 degrees.
That being all the files, here is an image of what a domain looks like after meshing is complete:
The image is a 2D slice of the mesh to show the refinement. Actual mesh being simulated is in 3D.
How to run this case yourself?
First things first, install and set up Linux using the Ubuntu distro (or any other distros) and install OpenFOAM. Make sure to install ParaView as well. Then, visit the repository here and clone this to your local OpenFOAM directory. Enter the case directory and run the following commands in this order:
blockMesh
surfaceFeatureExtract
snappyHexMesh -overwrite
checkMesh
If checkMesh returns "MESH OK", run touch *.foam to add a .foam file to the case directory. Open the *.foam file using ParaView and visually inspect the mesh using "Surface with Edges". Verify the domain mesh using the Slice filter. It should look something similar to the image above. Then proceed and run:
simpleFoam
After the solver converges (you don't have to wait, you can manually stop the solver at ~100 iterations using Ctrl+C and inspect- for a rough idea of the run), open the previously created *.foam file in ParaView and create some visualizations. Here are some examples:
The built-in case
OpenFOAM has a built in case with pre-configured files which follows a structure that is very similar to the one explained in this post. It is a steady-state incompressible flow simulation over a motorbike. Bring the case into your local directory by using the command cp -r $FOAM_TUTORIALS/incompressible/simpleFoam/motorBike, then paste it in your local directory. Go through all the files in the three folders line-by-line and compare it with this case to gain a much greater depth of understanding about steady-state external aerodynamics simulations in OpenFOAM.
Conclusion
Learning CFD is already challenging enough, and adding OpenFOAM to it only makes it more complex. Regardless of that, starting with this more complex form will be frustrating yet rewarding. It will also help genuinely knowing the physics when switching to a GUI-based CFD software.
Top comments (0)