Totorials: Incompressible Turbulent Flow¶
Turbulent Flat Plate¶
This tutorial considers the simulation of turbulent incompressible flow over a twodimensional sharp leadingedge flat plate using Caelus 9.04. Some basic steps to start a Caelus simulation for a turbulent flow environment will be shown such as specifying input data to define the boundary conditions, fluid properties, turbulence parameters and discretization/solver settings. Subsequently, the velocity contour over the plate will be visualised to identify the developed boundary layer. It will be further shown in sufficient detail to carry out Caelus simulation so that the user is able to reproduce accurately.
Objectives¶
Through this tutorials the user will be familiarised with setting up the Caelus simulation for steady, turbulent, incompressible flow over a sharp leadingedge flatplate in twodimensions. Further, the user will also be able to visualise the boundary layer. The following steps are carried out in this tutorial
 Background
A brief description about the problem
Geometry and freestream details
 Grid generation
Computational domain and boundary details
Computational grid generation
Exporting grid to Caelus
 Problem definition
Directory structure
Setting up boundary conditions, physical properties and control/solver attributes
 Execution of the solver
Monitoring the convergence
Writing the log files
 Results
Visualisation of turbulent boundary layer
Prerequisites¶
It is assumed that the user is familiar with the Linux command line environment using a terminal or Caelusconsole (for Windows OS) and Caelus is installed correctly with appropriate environment variables set. The grid used here is obtained from Turbulence Modeling Resource in a Plot3D format and is exported to Caelus format using Pointwise. However, the user is free to use their choice of grid generation tool to covert the Plot3D file to Caelus format.
Background¶
Turbulent flow over a flatplate configuration presents an ideal case to introduce the user with the turbulent simulation using Caelus. Here, the steadystate solution to the incompressible flow over the plate will be obtained, which results in a turbulent boundary layer. The shear stress distribution along the length of the wall and the velocity profile across the wall would be used to infer the development of the turbulent boundary layer. The user can look at the validation section for more details at Zero Pressure Gradient Flat Plate.
The flatplate length considered for this tutorial is L = 2.0 m and with a unit Reynolds number of \(5 \times 10^6\). Air is used as a fluid and a temperature of T = 300 K is assumed. Based on the Reynolds number and temperature, kinematic viscosity evaluates to \(\nu = 1.38872\times10^{5}~(m^2/s)\). A freestream velocity of \(U = 69.436113~m/s\) is used. In Figure 26, a schematic of the flatplate is shown. Note that the 2D plane of interest is in \(xz\) directions.
The freestream conditions that would be used is given in the below table
Fluid 
\(L~(m)\) 
\(Re/L~(1/m)\) 
\(U~(m/s)\) 
\(p~(m^2/s^2)\) 
\(T~(K)\) 
\(\nu~(m^2/s)\) 

Air 
0.3048 
\(5 \times 10^6\) 
69.436113 
Gauge (0) 
300 
\(1.38872\times10^{5}\) 
Grid Generation¶
The hexahedral grid used in this tutorial is obtained from Turbulence Modeling Resource that has 137 X 97 cells in \(xz\) directions respectively. The original 3D grid is in Plot3D and is then converted to Caelus compatible polyMesh
format.
The computational domain is a rectangular block that encompasses the flatplate. In Figure 27 below, the details of the boundaries in 2D (\(xz\) plane) that will be used is shown. The region of interest, which is highlighted in blue extends between \(0 \leq x \leq 2.0~m\), where the leadingedge is at \(x=0\). Due to the viscous nature of the flow, the velocity at the wall is zero which is represented through a noslip boundary wherein \(u,v,w = 0\). Upstream of the leading edge, a symmetry boundary at the wall will be used. The inlet boundary is placed at the start of the symmetry boundary and the outlet is placed at the exit of the flatplate the noslip wall. The entire top boundary will be again modelled as a symmetry plane.
The polyMesh
grid as noted earlier is in 3D. However, since the flow over a flatplate is twodimensional, the 2D plane that is considered here is in \(xz\) directions. It would therefore be ideal to have onecell thick in the direction (\(y\)), normal to the 2D plane of interest, where the flow is considered symmetry. The two \(xz\) planes obtained as a result of having 3D grid need boundary conditions to be specified. Since the flow is 2D, we do not need to solve for flow in 3D. This can easily be achieved in Caelus by specifying empty boundary condition for each of the two planes. As a consequence, the flow in \(y\) direction would be symmetry.
Note
A velocity value of \(v=0\) needs to be specified at appropriate boundaries although no flow is solved in the \(y\) direction.
In Figure 28, the 2D grid is shown which has 137 X 97 cells in \(xz\) directions respectively. To capture the turbulent boundary layer accurately, the grids are refined close to the wall and \(y^+\) is estimated to be less than 1. Due to this, no wallfunctions would be used to estimate the velocity gradients near the wall and integration is carried up to the wall.
Problem definition¶
In this section, several key instructions would be provided to setup the turbulent flatplate problem along with details of file configuration. A full working case can be found in:
/tutorials/incompressible/simpleSolver/ras/ACCM_flatPlate2D
However,the user is free to start the case setup from scratch consistent with the directory stucture discussed below.
Directory Structure
Note
All commands shown here are entered in a terminal window, unless otherwise mentioned
For setting up the problem, we need to further have few more subdirectories where relevant files can be created. Caelus requires time
, constant
and system
subdirectories. Since we will begin the simulation at time \(t = 0~s\), the time subdirectory should be just 0
.
In the 0
subdirectory, additional files are required for specifying the boundary properties. The following table lists the necessary files required based on the turbulence model chosen.
Parameter 
File name 

Pressure (\(p\)) 

Velocity (\(U\)) 

Turbulent viscosity (\(\nu\)) 

Turbulence field variable (\(\tilde{\nu}\)) 

Turbulent kinetic energy (\(k\)) 

Turbulent dissipation rate (\(\omega\)) 

As can be noted from the above table, we will be considering two turbulence models namely, SpalartAllmaras (SA) and \(k\omega\)  Shear Stress Transport (\(\rm{SST}\)) in the current exercise. These files set the dimensions, initialisation and boundary conditions to the problem, which also forms the three principle entries required.
The user should take into account that Caelus is case sensitive and therefore where applicable, the directory setup should be identical to what is shown here.
Boundary Conditions and Solver Attributes
Boundary Conditions
Initially, let us setup the boundary conditions. Referring back to Fig. %s:num:tfpdomaintutorials, the following are the boundary conditions that will be specified:
 Inlet
Velocity: Fixed uniform velocity \(u = 69.436113~m/s\) in \(x\) direction
Pressure: Zero gradient
Turbulence:
Spalart–Allmaras (Fixed uniform values of \(\nu_{t~\infty}\) and \(\tilde{\nu}_{\infty}\) as given in Turbulence freestream conditions for SA model)
\(k\omega~\textrm{SST}\) (Fixed uniform values of \(k_{\infty}\), \(\omega_{\infty}\) and \(\nu_{t~\infty}\) as given in Turbulence freestream conditions for k\omega~\rm{SST} model)
 Symmetry
Velocity: Symmetry
Pressure: Symmetry
Turbulence: Symmetry
 Noslip wall
Velocity: Fixed uniform velocity \(u, v, w = 0\)
Pressure: Zero gradient
Turbulence:
Spalart–Allmaras (Fixed uniform values of \(\nu_{t}=0\) and \(\tilde{\nu}=0\))
\(k\omega~\textrm{SST}\) (Zero gradient \(k\) and \(\omega\); Calculated \(\nu_t=0\); )
 Outlet
Velocity: Zero gradient velocity
Pressure: Fixed uniform gauge pressure \(p = 0\)
Turbulence:
Spalart–Allmaras (Calculated \(\nu_{t}=0\) and Zero gradient \(\tilde{\nu}\))
\(k\omega~\textrm{SST}\) (Zero gradient \(k\) and \(\omega\); Calculated \(\nu_t=0\); )
 Initialisation
Velocity: Fixed uniform velocity \(u = 69.436113\) in \(x\) direction
Pressure: Zero Gauge pressure
Turbulence:
Spalart–Allmaras (Fixed uniform values of \(\nu_{t~\infty}\) and \(\tilde{\nu}_{\infty}\) as given in Turbulence freestream conditions for SA model)
\(k\omega~\textrm{SST}\) (Fixed uniform values of \(k_{\infty}\), \(\omega_{\infty}\) and \(\nu_{t~\infty}\) as given in Turbulence freestream conditions for k\omega~\rm{SST} model)
Starting with the pressure, let us open p
using a text editor, which has the following contents.
/**\
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 2 0 0 0 0];
internalField uniform 0;
boundaryField
{
bottom
{
type symmetryPlane;
}
inflow
{
type zeroGradient;
}
left
{
type empty;
}
outflow
{
type fixedValue;
value uniform 0;
}
right
{
type empty;
}
top
{
type symmetryPlane;
}
wall
{
type zeroGradient;
}
}
// ************************************************************************* //
As can be seen, the above file begins with a dictionary named FoamFile
which contains the standard set of keywords for version, format, location, class and object names.
dimension
is used to specify the physical dimensions of the pressure field. Here, pressure is defined in terms of kinematic pressure with the units (\(m^2/s^2\)) written as
[0 2 2 0 0 0 0]
internalField
is used to specify the initial conditions. It can be either uniform or nonuniform. Since we have a 0 initial uniform gauge pressure, the entry is
uniform 0;
boundaryField
is used to specify the boundary conditions. In this case its the boundary conditions for pressure at all the boundary patches.
In a similar approach, let us open the file U
.
/**\
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
location "0";
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 1 0 0 0 0];
internalField uniform (69.4361 0 0);
boundaryField
{
bottom
{
type symmetryPlane;
}
inflow
{
type fixedValue;
value uniform (69.4361 0 0);
}
left
{
type empty;
}
outflow
{
type zeroGradient;
}
right
{
type empty;
}
top
{
type symmetryPlane;
}
wall
{
type noSlipWall;
}
}
// ************************************************************************* //
As detailed above, the principle entries for velocity field are self explanatory and the dimensions are typically for that of velocity with the units \(m/s\) ([0 1 1 0 0 0 0]
). Since we initialise the flow with a uniform freestream velocity, we set the internalField
to uniform (69.4361 0 0)
which represents three components of velocity. Similarly, inflow
boundary patch has three velocity components.
Similarly, the turbulent properties needed at the boundaries can be set. We begin with opening the file nut
, which is the turbulent kinematic viscosity and is shown below.
/**\
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object nut;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 1 0 0 0 0];
internalField uniform 2.9224e06;
boundaryField
{
bottom
{
type symmetryPlane;
}
inflow
{
type fixedValue;
value uniform 2.9224023e06;
}
left
{
type empty;
}
outflow
{
type calculated;
value uniform 0;
}
right
{
type empty;
}
top
{
type symmetryPlane;
}
wall
{
type fixedValue;
value uniform 0;
}
}
// ************************************************************************* //
Here, the turbulent viscosity is specified as kinematic and therefore the units are \(m^2/s\) ([0 2 1 0 0 0 0]
). The value of turbulence viscosity at freestream, specified at inflow patch is calculated as detailed in Turbulence freestream conditions for SA model and Turbulence freestream conditions for k\omega~\rm{SST} model for SST models respectively and is specified accordingly. The same value also goes for internalField
. Note that a fixedValue
of 0
is used for the wall
which suggests that on the wall, it is only the molecular (laminar) viscosity that prevails.
We shall now look at nuTilda
which is a turbulence field variable, specific to the SA model and has same units ([0 2 1 0 0 0 0]
) as kinematic turbulent viscosity. The details of which are given in Turbulence freestream conditions for SA model. In the file nuTilda
, the entries specified for the boundaryField
are identical to that of turbulent kinematic viscosity explained above.
/**\
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object nuTilda;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 1 0 0 0 0];
internalField uniform 4.166166e05;
boundaryField
{
bottom
{
type symmetryPlane;
}
inflow
{
type fixedValue;
value uniform 4.166166e05;
}
left
{
type empty;
}
outflow
{
type zeroGradient;
}
right
{
type empty;
}
top
{
type symmetryPlane;
}
wall
{
type fixedValue;
value uniform 0;
}
}
// ************************************************************************* //
We now proceed to files k
and omega
, specific to only \(k\omega~\rm{SST}\) model. As we know, \(k\omega~\rm{SST}\) is a turbulence model which solves for the turbulent kinetic energy and the specific rate of dissipation using two partial differential equations. Caelus therefore requires information about these to be specified when this model is used. Firstly, the file k
with the following contents is needed.
/**\
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object k;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 2 0 0 0 0];
internalField uniform 0.0010848;
boundaryField
{
bottom
{
type symmetryPlane;
}
inflow
{
type fixedValue;
value uniform 0.0010848;
}
left
{
type empty;
}
outflow
{
type zeroGradient;
}
right
{
type empty;
}
top
{
type symmetryPlane;
}
wall
{
type fixedValue;
value uniform 1e10;
}
}
// ************************************************************************* //
The unit of kinetic energy is \(m^2/s^2\) and this is set in dimensions
as [0 2 2 0 0 0 0]
. As with other turbulent quantities discussed above, the value of \(k\) (refer Turbulence freestream conditions for k\omega~\rm{SST} model needs to be specified for internalField
, inflow
and wall
. Please note that for wall boundaryField
with no wallfunction, a small, nonzero fixedValue
is required.
Next, the value for \(\omega\) is evaluated in omega
file as shown below and as detailed in Turbulence freestream conditions for k\omega~\rm{SST} model.
/**\
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object omega;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 0 1 0 0 0 0];
internalField uniform 8679.5135;
boundaryField
{
bottom
{
type symmetryPlane;
}
inflow
{
type fixedValue;
value uniform 8679.5135;
}
left
{
type empty;
}
outflow
{
type zeroGradient;
}
right
{
type empty;
}
top
{
type symmetryPlane;
}
wall
{
type omegaWallFunction;
value uniform 1;
}
}
// ************************************************************************* //
The unit of specific rate of dissipation for \(\omega\) is \(1/s\) which is set in dimensions
as [0 0 1 0 0 0 0]
. The internalField
and inflow
gets a fixedValue
. Note that for wall boundaryField
, we specify omegaWallFunction
and this is a model requirement and sets omega to the correct value near wall based on the \(y^+\). In conjunction, the value that goes with omegaWallFunction
can be anything and for simplicity its set to 1
.
Before setting up other parameters, it is important to ensure that the boundary conditions (inflow
, outflow
, top
, etc) specified in the above files should be the grid boundary patches (surfaces) generated by the grid generation tools and their names are identical. Further, the two boundaries in \(xz\) plane named here as left
and right
have empty
boundary conditions which forces Caelus to assume the flow to be in 2D. With this, the setting up of boundary conditions are completed.
Grid file and Physical Properties
The turbulent flatplate grid files is placed in the constant/polyMesh
subdirectory. Additionally, the physical properties are specified in various different files present in the directory constant
.
As you can see in the constant
directory, three files are listed in addition to the polyMesh
subdirectory. In the first file, RASProperties
, the ReynoldsAverageStress (RAS) model is specified as below. Note that depending on the turbulence model you wish to run with, the line that corresponds to that specific model should be enabled
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object RASProperties;
}
//
// For SpalaratAlamaras Model, enable the below line
RASModel SpalartAllmaras;
// For komega SST Model, enable the below line
// RASModel kOmegaSST;
turbulence on;
printCoeffs on;
Next, we look at the transportProperties
file, where transport model and kinematic viscosity is specified.
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
//
transportModel Newtonian;
nu nu [0 2 1 0 0 0 0] 1.388722E5;
As the viscous behaviour is Newtonian, the transportModel
is given using the keyword Newtonian
and the value of molecular (laminar) kinematic viscosity (nu
) is given having the units \(m^2/s\) ([0 2 1 0 0 0 0]
).
The final file in this class is the turbulenceProperties
file, which sets the simulationType
to RASModel
. Both SA and \(k\omega~\rm{SST}\) are classified as Reynolds Average Stress (RAS) models.
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object turbulenceProperties;
}
//
simulationType RASModel;
Controls and Solver Attributes
In this section, the files required to control the simulation and specifying the type of discretization method along with the linear solver settings are provided. These are placed in the system
directory.
First, we begin with the controlDict
file as below
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
//
application simpleSolver;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 20000;
deltaT 1;
writeControl runTime;
writeInterval 1000;
purgeWrite 0;
writeFormat ascii;
writePrecision 12;
writeCompression true;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
//
// Function Objects to obtain forces
functions
{
forces
{
type forces;
functionObjectLibs ("libforces.so");
patches ( wall );
CofR (0 0 0);
rhoName rhoInf;
rhoInf 1.347049;
writeControl timeStep;
writeInterval 50;
}
}
As can be noted in the above file, simpleSolver
solver is used and the simulation begins at t = 0 s
. This logically explains the need for 0
directory where the data files are read at the beginning of the run. Therefore, the keyword startFrom
is set to startTime
, where startTime
would be 0
. Since the simulation is steadystate we specify the total number of iterations as a keyword for endTime
. Via the writeControl
and writeInternal
keywords, the interval at which the solutions are saved can be specified. Also included is the function object to obtain the force over the wall every 50
iterations. Note that for obtaining the force, the freestream density (rhoInf
) is required and is specified with the value.
The discretization schemes for the finite volume discretization that will be used is set through the fvSchemes
file shown below
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object fvSchemes;
}
//
ddtSchemes
{
default steadyState;
}
gradSchemes
{
default Gauss linear;
grad(p) Gauss linear;
grad(U) Gauss linear;
}
divSchemes
{
default none;
div(phi,U) Gauss linearUpwind grad(U);
div(phi,nuTilda) Gauss upwind; // Will be used for SA model only
div(phi,k) Gauss upwind; // will be used for kepsilon & komega only
div(phi,omega) Gauss upwind; // Will be used for komega model only
div((nuEff*dev(T(grad(U))))) Gauss linear;
div(phi,symm(grad(U))) Gauss linear;
}
laplacianSchemes
{
default none;
laplacian(nu,U) Gauss linear corrected;
laplacian(nuEff,U) Gauss linear corrected;
laplacian(DnuTildaEff,nuTilda) Gauss linear corrected; // Will be used for SA model only
laplacian(DkEff,k) Gauss linear corrected; // Will be used for komega & komega only
laplacian(DomegaEff,omega) Gauss linear corrected; // Will be used for komega model only
laplacian(rAUf,p) Gauss linear corrected;
laplacian(1,p) Gauss linear corrected;
}
interpolationSchemes
{
default linear;
interpolate(HbyA) linear;
}
snGradschemes
{
default corrected;
}
The linear solver controls and tolerances are set in fvSolution
as given below
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSolution;
}
//
solvers
{
p
{
solver PCG;
preconditioner SSGS;
tolerance 1e8;
relTol 0.01;
}
U
{
solver PBiCGStab;
preconditioner USGS;
tolerance 1e7;
relTol 0.01;
}
"(komeganuTilda)"
{
solver PBiCGStab;
preconditioner USGS;
tolerance 1e08;
relTol 0;
}
}
SIMPLE
{
nNonOrthogonalCorrectors 1;
pRefCell 0;
pRefValue 0;
}
// relexation factors
relaxationFactors
{
p 0.3;
U 0.5;
nuTilda 0.5;
k 0.5;
omega 0.5;
}
Here, different linear solvers are used to solve velocity, pressure and turbulence quantities. We also set the nNonOrthogonalCorrectors
to 1 in this case. Further, relaxation is set on the primary and turbulent variables so that the solution is more stable. Furthermore, the relTol
is not set to 0
unlike a timeaccurate setup. This is because we are solving for a steadystate solution and a very low (\(\approx 0\)) tolerance at every iteration is not required as the entire system of equations converges to the global tolerance
set as the simulation progresses to steady state.
Now the setup of the directory structure with all the relevant files the directory tree should appear identical to the one shown below
.
├── 0
│ ├── epsilon
│ ├── k
│ ├── nut
│ ├── nuTilda
│ ├── omega
│ ├── p
│ └── U
├── constant
│ ├── polyMesh
│ │ ├── boundary
│ │ ├── faces
│ │ ├── neighbour
│ │ ├── owner
│ │ └── points
│ ├── RASProperties
│ ├── transportProperties
│ └── turbulenceProperties
└── system
├── controlDict
├── fvSchemes
└── fvSolution
Execution of the solver¶
Prior to execution of solver, renumbering of the grid/mesh needs to be performed in addition to checking the quality of the grid/mesh. Renumbering the gridcell list is vital to reduce the matrix bandwidth while quality check gives us the mesh statistics. Renumbering and mesh quality can be determined by executing the following from the top directory.
caelus run  renumberMesh overwrite
caelus run  checkMesh
At this stage, it is suggested that the user should take note of the matrix bandwidth before and after the mesh renumbering. When the checkMesh
is performed, the mesh statistics are shown as below
/**\
Caelus 8.04
Web: www.caeluscml.com
\**/
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Checking geometry...
Overall domain bounding box (0.06 0 0.03) (1.2192 0.15 0.055)
Mesh (nonempty, nonwedge) directions (1 1 0)
Mesh (nonempty) directions (1 1 0)
All edges aligned with or perpendicular to nonempty directions.
Boundary openness (5.80542e19 1.1194e17 1.1403e14) OK.
Max cell openness = 2.2093e16 OK.
Max aspect ratio = 55.555 OK.
Minimum face area = 1e08. Maximum face area = 0.000138887. Face area magnitudes OK.
Min volume = 2.5e10. Max volume = 2.50831e07. Total volume = 0.004797. Cell volumes OK.
Mesh nonorthogonality Max: 0 average: 0
Nonorthogonality check OK.
Face pyramids OK.
Mesh skewness Max: 3.85044e13 average: 9.40402e15 OK.
Coupled point location match (average 0) OK.
Mesh OK.
End
In the above terminal output, we get Failed 1 mesh checks.
and this is because of the high aspect ratio meshes present immediate to the wall due to very low (\(<< 1~y^+\)). However, Caelus can solve on this mesh. The next step is to execute the solver and monitoring the progress of the solution. The solver is always executed from the top directory.
caelus run l myturbulentflatplate.log  simpleSolver
With the execution of the above command, the simulation begins and the progress of the solution is written to the specified log file (myturbulentflatplate.log
). The log file can be further processed to look at the convergence history and this can be done as follows
caelus logs w myturbulentflatplate.log
This allows you to look at the convergence of different variables with respect to the number of iterations carried out. In Fig. %s:num:tfpconvergencetutorials pressure convergence is shown.
Results¶
The turbulent flow over the flat plate is shown here through velocity magnitude contours for SA model. In Fig. %s:num:tfpvelocitytutorials the boundary layer over the entire flatplate and in the region up to \(x=0.10~m\) is emphasised. The growth of the boundary layer can be seen very clearly. Since the Reynolds number of the flow is reasonably high, the turbulent boundary layer seems thin in comparison to the length of the plate.
Bump in a Channel¶
The simulation of turbulent flow over a twodimensional bump in a channel is considered in this tutorial and will be performed using Caelus 9.04. As with the other tutorials, setting up the directory structure, fluid properties, boundary conditions, turbulence properties etc will be shown. Further to this, visualisation of the solution will be shown to look at the velocity and pressure contours over the bump surface. These steps would be shown in sufficient details so that the user is able to reproduce the tutorial accurately.
Objectives¶
Some of the main objectives of this tutorial would be for the user to get familiarise with setting up the Caelus simulation for steady, turbulent, incompressible flow over a twodimensional bump in a channel and be able to postprocess the desired solution. Following would be some of the steps that would be covered.
 Background
A brief description about the problem
Geometry and freestream details
 Grid generation
Computational domain and boundary details
Computational grid generation
Exporting grid to Caelus
 Problem definition
Directory structure
Setting up boundary conditions, physical properties and control/solver attributes
 Execution of the solver
Monitoring the convergence
Writing the log files
 Results
Visualisation of flow near the bump
Prerequisites¶
The user should be familiar with a Linux command line environment via a terminalor caelusconsole (For Windows OS). It is also assumed that Caelus is installed correctly with appropriate environment variables set. The grid used here is obtained from Turbulence Modeling Resource in a Plot3D format and is exported to Caelus format using Pointwise. However, the user is free to use their choice of grid generation tool to covert the Plot3D file to Caelus format.
Background¶
Turbulent flow over a bump in a channel is quite similar to a flatplate flow, except that due to the curvature effect, a pressure gradient is developed. The flow would be assumed to be steadystate and incompressible. To demonstrate the effect of curvature, the skin friction distribution along the bump will be plotted. For further information on this case, the user can refer to Twodimensional Bump in a Channel.
The bump, as shown in the schematic below in Figure 31 has a upstream and a downstream flatplate region that begins at x = 0 m and x = 1.5 m respectively, which gives a total length of L = 1.5 m. The flow has a unit Reynolds number of \(3 \times 10^6\) and Air is used as a fluid with a temperature of 300 K. Based on these values, kinematic viscosity can be evaluated to \(\nu = 2.314537 \times 10^{5} m^2/s\). To match the required Reynolds number, a freestream velocity of U = 69.436 m/s would be used. Note that the twodimensional plane considered here is in \(xz\) directions.
The freestream conditions that will be used is given in the below table (Freestream conditions)
Fluid 
\(L~(m)\) 
\(Re/L~(1/m)\) 
\(U~(m/s)\) 
\(p~(m^2/s^2)\) 
\(T~(K)\) 
\(\nu~(m^2/s)\) 

Air 
1.5 
\(3 \times 10^6\) 
69.436113 
Gauge (0) 
300 
\(2.314537\times10^{5}\) 
Grid Generation¶
The hexahedral grid used in this tutorial is obtained from Turbulence Modeling Resource that contains 704 X 320 cells in \(xz\) directions respectively. The grid originally is in Plot3D format and is converted to Caelus compatible polyMesh
format.
The computational domain is a rectangular channel encompassing the bump. Figure 32 shown below identifies the boundary conditions in a twodimensional plane (\(xz\)). The bump region, highlighted in blue extends between \(0 \leq x \leq 1.5~m\), where the velocity at the wall is zero, wherein \(u,v,w=0\) represented through a noslip boundary. Upstream and downstream of the bump, a symmetry boundary at the wall is used. The inlet and outlet are placed at the end of the symmetry as depicted in the figure and the top boundary has a symmetry condition.
The polyMesh
grid obtained from the conversion of Plot3D is in 3D. However, the flow over a bump is twodimensional and is solved in a 2D plane with \(xz\) directions. Therefore, ideally we can have cells with onecell thick in the direction (\(y\)), normal to the 2D plane, where the flow can be assumed to be symmetry. The two \(xz\) planes obtained as a result of having a 3D grid require boundary conditions to be specified. As the flow is assumed to be 2D, we do not need to solve the flow in 3D and this can easily be achieved in Caelus by specifying empty
boundary condition for each of the two planes. This consequently allows for the flow in \(y\) direction to be symmetry.
Note
A velocity value of \(v=0\) needs to be specified at appropriate boundaries although no flow is solved in the \(y\) direction.
In Figure 33 above, the 2D grid is shown with a distribution of 704 X 320 in \(xz\) directions respectively. The inset focuses the distribution in the region between \(0 \leq x \leq 1.5~m\). As can be seen, the distribution of the grids in the flow normal direction is finer near the wall to capture the turbulent boundary layer more accurately and it is estimated that \(y^+\) is less than 1 for the chosen grid and therefore, no wallfunction has been used and the integration is carried out up to the wall.
Problem definition¶
This section deals with several key instructions need to setup the turbulent flow over a bump. A full working case of this can be found in:
/tutorials/incompressible/simpleSolver/ras/ACCM_bump2D/
The user is free to start the case setup from scratch consistent with the directory stucture discussed below.
Directory Structure
Note
All commands shown here are entered in a terminal window, unless otherwise mentioned
In order to setup the problem Caelus requires time
, constant
and system
subdirectories within the main working directory. Typically, the simulations are started at time \(t = 0~s\), which requires a time subdirectory to be 0
.
Within the 0
subdirectory, additional files specifying the boundary properties are present. The below table lists the necessary files required based
on the turbulence model chosen
Parameter 
File name 

Pressure (\(p\)) 

Velocity (\(U\)) 

Turbulent viscosity (\(\nu\)) 

Turbulence field variable (\(\tilde{\nu}\)) 

Turbulent kinetic energy (\(k\)) 

Turbulent dissipation rate (\(\omega\)) 

In this tutorial, we will be considering two turbulence models namely, SpalartAllmaras (SA) and \(k\omega\)  Shear Stress Transport (\(\rm{SST}\)). The contents of the files listed above sets the dimensions, initialisation and boundary conditions to the defining problem, which also
forms three principle entries required. Firstly, we begin with looking at these files in the 0
subdirectory
The user should take into account that Caelus is case sensitive and therefore the directory setup should be identical to what is shown here.
Boundary Conditions and Solver Attributes
Boundary Conditions
Referring back to Figure 32, the following are the boundary conditions that will be specified:
 Inlet
Velocity: Fixed uniform velocity \(u = 69.436113~m/s\) in \(x\) direction
Pressure: Zero gradient
Turbulence:
Spalart–Allmaras (Fixed uniform values of \(\nu_{t~\infty}\) and \(\tilde{\nu}_{\infty}\) as given in Turbulence freestream conditions for SA model)
\(k\omega~\rm{SST}\) (Fixed uniform values of \(k_{\infty}\), \(\omega_{\infty}\) and \(\nu_{t~\infty}\) as given in Turbulence freestream conditions for k\omega~\rm{SST} model)
 Symmetry
Velocity: Symmetry
Pressure: Symmetry
Turbulence: Symmetry
 Noslip wall
Velocity: Fixed uniform velocity \(u, v, w = 0\)
Pressure: Zero gradient
Turbulence:
Spalart–Allmaras (Fixed uniform values of \(\nu_{t}=0\) and \(\tilde{\nu}=0\))
\(k\omega~\rm{SST}\) (Fixed uniform values of \(k =<<0\) and \(\nu_t=0\); \(\omega\) = omegaWallFunction)
 Outlet
Velocity: Zero gradient velocity
Pressure: Fixed uniform gauge pressure \(p = 0\)
Turbulence:
Spalart–Allmaras (Calculated \(\nu_{t}=0\) and Zero gradient \(\tilde{\nu}\))
\(k\omega~\rm{SST}\) (Zero gradient \(k\) and \(\omega\); Calculated \(\nu_t=0\); )
 Initialisation
Velocity: Fixed uniform velocity \(u = 69.436113\) in \(x\) direction
Pressure: Zero Gauge pressure
Turbulence:
Spalart–Allmaras (Fixed uniform values of \(\nu_{t~\infty}\) and \(\tilde{\nu}_{\infty}\) as given in Turbulence freestream conditions for SA model)
\(k\omega~\rm{SST}\) (Fixed uniform values of \(k_{\infty}\), \(\omega_{\infty}\) and \(\nu_{t~\infty}\) as given in Turbulence freestream conditions for k\omega~\rm{SST} model)
We begin with p
, the pressure file using a text editor with the following content
/**\
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 2 0 0 0 0];
internalField uniform 0;
boundaryField
{
bottom
{
type symmetryPlane;
}
inflow
{
type zeroGradient;
}
left
{
type empty;
}
outflow
{
type fixedValue;
value uniform 0;
}
right
{
type empty;
}
top
{
type symmetryPlane;
}
wall
{
type zeroGradient;
}
}
// ************************************************************************* //
From the above information, it can be seen that the file begins with a dictionary named FoamFile
which contains the standard set of keywords for version, format, location, class and object names. The explanation of the principle entries are as follows
dimension
is used to specify the physical dimensions of the pressure field. Here, pressure is defined in terms of kinematic pressure with the units (\(m^2/s^2\)) written as
[0 2 2 0 0 0 0]
internalField
is used to specify the initial conditions. It can be either uniform or nonuniform. Since we have a 0 initial uniform gauge pressure, the entry is
uniform 0;
boundaryField
is used to specify the boundary conditions. In this case its the boundary conditions for pressure at all the boundary patches.
Similarly, we have the file U
with the following information
/**\
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
location "0";
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 1 0 0 0 0];
internalField uniform (69.4361 0 0);
boundaryField
{
bottom
{
type symmetryPlane;
}
inflow
{
type fixedValue;
value uniform (69.4361 0 0);
}
left
{
type empty;
}
outflow
{
type zeroGradient;
}
right
{
type empty;
}
top
{
type symmetryPlane;
}
wall
{
type noSlipWall;
}
}
// ************************************************************************* //
As with the pressure, the principle entries for velocity field are selfexplanatory and the dimensions are typically for that of velocity with the units \(m/s\) ([0 1 1 0 0 0 0]
). Since the initialisation of the flow is with a uniform freestream velocity, we should set the internalField
to uniform (69.4361 0 0)
representing three components of velocity. In a similar manner, inflow
boundary patch has three velocity components.
In addition to p
and U
, the turbulent properties are also needed at the boundary patches and these can be set in a similar process. We begin with the file nut
, which corresponds to turbulent kinematic viscosity as shown below.
/**\
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object nut;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 1 0 0 0 0];
internalField uniform 4.87067e06;
boundaryField
{
bottom
{
type symmetryPlane;
}
inflow
{
type fixedValue;
value uniform 4.87067e06;
}
left
{
type empty;
}
outflow
{
type calculated;
value uniform 0;
}
right
{
type empty;
}
top
{
type symmetryPlane;
}
wall
{
type fixedValue;
value uniform 0;
}
}
// ************************************************************************* //
Here, the turbulent viscosity is specified as kinematic and therefore the units are \(m^2/s\) ([0 2 1 0 0 0 0]
). The value of turbulence viscosity at freestream, specified at inflow patch is calculated as detailed in Turbulence freestream conditions for SA model and Turbulence freestream conditions for k\omega~\rm{SST} model for SST models respectively and is specified accordingly. The same value also goes for internalField
. Note that a fixedValue
of 0
is used for the wall
which suggests that on the wall, it is only the molecular (laminar) viscosity that prevails.
The next variable is the nuTilda
which is a turbulence field variable, specific to only SA model and has the same units ([0 2 1 0 0 0 0]
) as kinematic turbulent viscosity. The details of which are given in Turbulence freestream conditions for SA model. The following contents given in the file nuTilda
and the entries specified for the boundaryField
are identical to that of turbulent kinematic viscosity explained above.
/**\
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object nuTilda;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 1 0 0 0 0];
internalField uniform 6.943611e05;
boundaryField
{
bottom
{
type symmetryPlane;
}
inflow
{
type fixedValue;
value uniform 6.943611e05;
}
left
{
type empty;
}
outflow
{
type zeroGradient;
}
right
{
type empty;
}
top
{
type symmetryPlane;
}
wall
{
type fixedValue;
value uniform 0;
}
}
// ************************************************************************* //
We now proceed to files k
and omega
, specific to only \(k\omega~\rm{SST}\) model. As we know, \(k\omega~\rm{SST}\) is a turbulence model which solves for the turbulent kinetic energy and the specific rate of dissipation using two partial differential equations. Caelus therefore requires information about these to be specified at the boundary patches when this model is chosen as shown below.
/**\
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object k;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 2 0 0 0 0];
internalField uniform 0.0010848;
boundaryField
{
bottom
{
type symmetryPlane;
}
inflow
{
type fixedValue;
value uniform 0.0010848;
}
left
{
type empty;
}
outflow
{
type zeroGradient;
}
right
{
type empty;
}
top
{
type symmetryPlane;
}
wall
{
type fixedValue;
value uniform 1e10;
}
}
// ************************************************************************* //
The unit of kinetic energy is \(m^2/s^2\) and this is set in dimensions
as [0 2 2 0 0 0 0]
. As with other turbulent quantities discussed above, the value of \(k\) (refer Turbulence freestream conditions for k\omega~\rm{SST} model) needs to be specified for internalField
, inflow
and wall
. Please note that for wall boundaryField
with no wallfunction, a small, nonzero fixedValue
is required.
We now evaluate the value for \(\omega\) in the omega
file as shown below and as detailed in Turbulence freestream conditions for k\omega~\rm{SST} model.
/**\
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object omega;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 0 1 0 0 0 0];
internalField uniform 5207.65;
boundaryField
{
bottom
{
type symmetryPlane;
}
inflow
{
type fixedValue;
value uniform 5207.65;
}
left
{
type empty;
}
outflow
{
type zeroGradient;
}
right
{
type empty;
}
top
{
type symmetryPlane;
}
wall
{
type omegaWallFunction;
value uniform 1;
}
}
// ************************************************************************* //
The unit of specific rate of dissipation for \(\omega\) is \(1/s\) which is set in dimensions
as [0 0 1 0 0 0 0]
. The internalField
and inflow
gets a fixedValue
. Note that for wall boundaryField
, we specify omegaWallFunction
and this is a model requirement and sets omega to the correct value near wall based on the \(y^+\). In conjunction, the value that goes with omegaWallFunction
can be anything and for simplicity its set to 1
.
Before proceeding further, it is important to ensure that the boundary conditions (inflow
, outflow
, top
, etc) added in the above files should be the grid boundary patches (surfaces) generated by the grid generation tool and their names are identical. Additionally, the two boundaries \(xz\) plane named here as left
and right
have empty
boundary conditions which forces Caelus to assume the flow to be in twodimensions. With this, the setting up of the boundary conditions are complete.
Grid file and Physical Properties
The grid file for the turbulentbump need to be placed in constant/polyMesh
subdirectory. In addition to this, the physical properties are specified in various different files present in the constant
directory. The three files that are listed in addition to the polyMesh
subdirectory set the physical properties. The first one, RASProperties
in which the ReynoldsAverageStress (RAS) is specified, is shown below. Please note that depending on the turbulence model you wish to run with, the line that corresponds to that specific model should be enabled.
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object RASProperties;
}
//
// For SpalaratAlamaras Model, enable the below line
RASModel SpalartAllmaras;
// For komega SST Model, enable the below line
// RASModel kOmegaSST;
turbulence on;
printCoeffs on;
Next, kinematic viscosity is specified in the transportProperties
file, as shown below
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
//
transportModel Newtonian;
nu nu [0 2 1 0 0 0 0] 2.314537E5;
As the viscous behaviour is Newtonian, the transportModel
is given using the keyword Newtonian
and the value of molecular (laminar) kinematic viscosity (nu
) is given having the units \(m^2/s\) ([0 2 1 0 0 0 0]
).
The final file in this class is the turbulenceProperties
file, which sets the simulationType
to RASModel
. Both SA and \(k\omega~\rm{SST}\) are classified as Reynolds Average Stress (RAS) models.
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object turbulenceProperties;
}
//
simulationType RASModel;
Controls and Solver Attributes
This section details the files require to control the simulation and the specifying discretization methods in addition to the linear solver settings. These files are placed in the system
directory.
The controlDict
file contains the following details
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
//
application simpleSolver;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 20000;
deltaT 1;
writeControl runTime;
writeInterval 1000;
purgeWrite 0;
writeFormat ascii;
writePrecision 12;
writeCompression true;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
//
// Function Objects to obtain forces
functions
{
forces
{
type forces;
functionObjectLibs ("libforces.so");
patches ( wall );
CofR (0 0 0);
rhoName rhoInf;
rhoInf 0.80822;
writeControl timeStep;
writeInterval 50;
}
}
Referring to the above information, some explanation is needed. Here, simpleSolver
is used and the simulation begins at t = 0 s
. This now explains the logical need for having a 0
directory where the data files are read at the beginning of the run, which is t = 0 s
in this case. Therefore, the keyword startFrom
is set to startTime
, where startTime
would be 0
. The simulation would be carried out as steadystate and therefore we require to specify the total number of iterations as a keyword for endTime
. Through the writeControl
and writeInterval
keywords, the solution intervals at which they are saved can be specified. Also note that a function object to obtain the force over the wall for every 50 iterations is included. In order to obtain this, a freestream density (rhoInf) need to be specified.
The discretization schemes for the finite volume discretization that will be used should be set through the fvSchemes
file show below
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object fvSchemes;
}
//
ddtSchemes
{
default steadyState;
}
gradSchemes
{
default Gauss linear;
grad(p) Gauss linear;
grad(U) Gauss linear;
}
divSchemes
{
default none;
div(phi,U) Gauss linearUpwind grad(U);
div(phi,nuTilda) Gauss upwind; // Will be used for SA model only
div(phi,k) Gauss upwind; // will be used for kepsilon & komega only
div(phi,omega) Gauss upwind; // Will be used for komega model only
div((nuEff*dev(T(grad(U))))) Gauss linear;
div(phi,symm(grad(U))) Gauss linear;
}
laplacianSchemes
{
default none;
laplacian(nu,U) Gauss linear corrected;
laplacian(nuEff,U) Gauss linear corrected;
laplacian(DnuTildaEff,nuTilda) Gauss linear corrected; // Will be used for SA model only
laplacian(DkEff,k) Gauss linear corrected; // Will be used for komega & komega only
laplacian(DomegaEff,omega) Gauss linear corrected; // Will be used for komega model only
laplacian(rAUf,p) Gauss linear corrected;
laplacian(1,p) Gauss linear corrected;
}
interpolationSchemes
{
default linear;
interpolate(HbyA) linear;
}
snGradschemes
{
default corrected;
}
The linear solver controls and tolerances are set in fvSolution
as given below
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSolution;
}
//
solvers
{
p
{
solver PCG;
preconditioner SSGS;
tolerance 1e8;
relTol 0.01;
}
U
{
solver PBiCGStab;
preconditioner USGS;
tolerance 1e7;
relTol 0.01;
}
"(komeganuTilda)"
{
solver PBiCGStab;
preconditioner USGS;
tolerance 1e08;
relTol 0;
}
}
SIMPLE
{
nNonOrthogonalCorrectors 1;
pRefCell 0;
pRefValue 0;
}
// relexation factors
relaxationFactors
{
p 0.3;
U 0.5;
nuTilda 0.5;
k 0.5;
omega 0.5;
}
The user should note that in the fvSolution
file, different linear solvers are used to solve for velocity, pressure and turbulence quantities. We also set the nNonOrthogonalCorrectors
to 1 for this case. To ensure the stability of the solution, the relaxation is set to primary and turbulent variables. The relTol
is set to 0
unlike a timeaccurate setup as we are solving for a steadystate solution and a very low (\(\approx 0\)) tolerance at every iteration is unnecessary. Since the entire system of equations converge to a global set tolerance
the convergence would occur as the solution progresses to a steady state.
With this, the setup of the directory structure with all the relevant files are complete and the directory tree should appear identical to the one shown below
.
├── 0
│ ├── epsilon
│ ├── k
│ ├── nut
│ ├── nuTilda
│ ├── omega
│ ├── p
│ └── U
├── constant
│ ├── polyMesh
│ │ ├── boundary
│ │ ├── faces
│ │ ├── neighbour
│ │ ├── owner
│ │ └── points
│ ├── RASProperties
│ ├── transportProperties
│ └── turbulenceProperties
└── system
├── controlDict
├── DecomposeParDict
├── fvSchemes
└── fvSolution
Execution of the solver¶
It is always important to renumber and to do a quality check on the grid/mesh before executing the solver. Renumbering reduces the matrix bandwidth whereas the quality check shows the mesh statistics. These two can be performed by executing the following commands from the top working directory
caelus run  renumberMesh overwrite
caelus run  checkMesh
At this stage, it is suggested that the user should take note of the bandwidth before and after the mesh renumbering. When the checkMesh
is performed, the mesh statistics are shown as below
/**\
Caelus 8.04
Web: www.caeluscml.com
\**/
Checking geometry...
Overall domain bounding box (10 10 0) (40 10 0.537713)
Mesh (nonempty, nonwedge) directions (1 1 0)
Mesh (nonempty) directions (1 1 0)
All edges aligned with or perpendicular to nonempty directions.
Boundary openness (2.57817e19 1.67414e19 4.29222e16) OK.
Max cell openness = 2.19645e16 OK.
Max aspect ratio = 3.66844 OK.
Minimum face area = 0.00895343. Maximum face area = 0.586971. Face area magnitudes OK.
Min volume = 0.00481437. Max volume = 0.315622. Total volume = 536.025. Cell volumes OK.
Mesh nonorthogonality Max: 14.6136 average: 1.75565
Nonorthogonality check OK.
Face pyramids OK.
Mesh skewness Max: 0.206341 average: 0.00112274 OK.
Coupled point location match (average 0) OK.
Mesh OK.
The output of the checkMesh
indicates that the mesh check has failed through the final message``Failed 1 mesh checks.`` and this is because of the high aspect ratio meshes present immediate to the wall due to very low (\(<< 1~y^+\)). However, Caelus will solve on this mesh.
In this tutorial, it will be shown further to utilise the multicore capability of CPUs thus performing a parallel computation for large grids, such as the one considered here. At first the grid has to be decomposed into smaller pieces that can be solved by each single CPU core. The number of decomposition should be equal to the number of CPU core available for parallel computing. The decomposition should be carried out through a file decomposeParDict
present in the system
subdirectory with the following content,
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object decomposeParDict;
}
//
numberOfSubdomains 4; // It is suggested that the numberOfSubdomains be increased based on available resources for validation cases and to reduce the computation time.
method simple;
simpleCoeffs
{
n (4 1 1);
delta 0.001;
}
In the above file, the keyword numberOfSubdomains
defines the number of decomposed subdomains. In this case, the grid is partitioned into 4 subdomains. We use simple
as the method of decomposition and n
is used to specify the number of decomposition that should be carried out in x, y and z directions respectively. In this case (4 1 1)
performs 4 decompositions in x direction and 1 in both y and z directions. The execution to decompose the grid is carried out again from the top directory as follows
caelus run  decomposePar
Now the decomposition should begin and the details of which are displayed in the terminal window. Subsequently, 4 processor directories will be generated as shown below
0 constant processor0 processor1 processor2 processor3 system
The solver can now be executed for parallel computation in the host machine from the top directory using the following command
caelus run p l myturbulentbump.log  simpleSolver
Note that here it is assumed that the parallel computing is available in the host machine. With the execution of the above commands, the simulation begins and the progress is written to the specified log file (myturbulentbump.log
).
The log file can be further processed to look at the convergence history and this can be done as follows
caelus logs w myturbulentbump.log
This allows you to look at the convergence of the variables with respect to the number of iterations carried out and the Figure 34 indicates the same for pressure.
NACA Airfoil¶
In this tutorial, the turbulent flow over a twodimensional NACA 0012 airfoil at two angles of attack, namely \(0^\circ\) and \(10^\circ\) will be considered. Caelus 9.04 will be used and the basic steps to setup the directory structure, fluid properties, boundary conditions, turbulence properties etc will be shown. Visualisation of pressure and velocity over the airfoil are also shown. With these, the user should be able to reproduce the tutorial accurately.
Objectives¶
The user will get familiar in setting up Caelus simulation for steady, turbulent, incompressible flow over a twodimensional airfoils at different angles of attack. Alongside, the user will be able to decompose the mesh on several CPUs performing a parallel simulation. Some of the steps that would be detailed are as follows
 Background
A brief description about the problem
Geometry and freestream details
 Grid generation
Computational domain and boundary details
Computational grid generation
Exporting grid to Caelus
 Problem definition
Directory structure
Setting up boundary conditions, physical properties and control/solver attributes
 Execution of the solver
Monitoring the convergence
Writing the log files
 Results
Visualisation of flow over the airfoil
Prerequisites¶
It is understood that the user will be familiar with the Linux command line environment via a terminal or Caelusconsole (For Windows OS) and Caelus is installed corrected with appropriate environment variables set. The grid for this case is obtained from Turbulence Modeling Resource as a Plot3D format and is converted to Caelus using Pointwise. The user is however free to use their choice of grid generation tool to convert the original Plot3D grid to Caelus readable format.
Background¶
Turbulent flow over airfoils is an interesting example to highlight some of the capabilities of Caelus. Here, the flow undergoes rapid expansion due to strong surface curvatures thereby inducing pressure and velocity gradients along the surface. Depending on shape of the curvature, adverse or favourable pressure gradients can exist on either side. These can be examined through surface quantities like pressure and skinfriction distributions. The user can refer to the verification and validation of this case at Twodimensional NACA 0012 Airfoil.
The schematic of NACA 0012 airfoil at two angles of attack are shown in Figure 36 for a twodimensional profile. A chord length (C) of 1.0 m is considered for both and has a Reynolds number of \(6 \times 10^6\). The flow is assumed to be Air with a freestream temperature of 300 K. Considering these values, the freestream velocity can be evaluated to U = 52.077 m/s. Note that the geometric plane considered for twodimensionality is in \(xz\) directions.
The freestream conditions are given in the below table
Fluid 
\(C~(m)\) 
\(Re/L~(1/m)\) 
\(U~(m/s)\) 
\(p~(m^2/s^2)\) 
\(T~(K)\) 
\(\nu~(m^2/s)\) 

Air 
1.0 
\(6 \times 10^6\) 
52.0770 
Gauge (0) 
300 
\(8.6795\times10^{6}\) 
As noted earlier, flow at two angles of attack (\(\alpha\)) will be considered in this tutorial. In order to obtain a freestream velocity of 52.0770 m/s at \(\alpha = 0^\circ\) and \(10^\circ\), the velocity components in \(x\) and \(z\) have to be resolved. The following table provides these values
\(\alpha~\rm{Degrees}\) 
\(u~(m/s)\) 
\(w~(m/s)\) 

\(0^\circ\) 
52.0770 
0.0 
\(10^\circ\) 
51.2858 
9.04307 
Grid Generation¶
The structured grid used for this tutorial can cells obtained from Turbulence Modeling Resource in Plot3D format that contains 512 around the airfoil and 256 cells in the flow normal direction. This should be then converted to polyMesh
format.
The computational domain for the NACA 0012 airfoil is shown in Figure 37 along with the boundary conditions. A large domain exists around the airfoil (highlighted in blue) extending 500 chord lengths in the radial direction and the inlet condition is given for the entire boundary highlighted in green, whereas the outlet is placed at the exit plane which is about \(x \approx 500~m\). The velocity on the airfoil surface is zero, wherein \(u, v, w = 0\) represented through a noslip boundary.
The polyMesh
grid is in threedimensions, however the flow over airfoils can be assumed to be 2D at low angles of attack and is solved here in \(xz\) directions. Therefore, a onecell thick grid normal to the 2D flow plane is sufficient, where the flow can be assumed to be symmetry. The two \(xz\) planes that are prevalent require boundary conditions to be specified. Since the flow is assumed to be 2D, we do not need to solve the flow in the thirddimension and this can be easily achieved in Caelus by specifying empty boundary conditions for each of the two planes. Consequently, the flow will be treated as symmetry in \(y\) direction.
Note
A velocity value of \(v=0\) needs to be specified at appropriate boundaries although no flow is solved in the \(y\) direction.
The 2D airfoil grid in \(xz\) plane is shown in Figure 38 which has a distribution of 512 X 256 cells. The grid in the vicinity of airfoil is shown as an inset and a very fine distribution can be noted very close to the wall. It was estimated that \(y^+\) is less than 1 to capture the turbulent boundary layer accurately and no wallfunction is used.
Problem definition¶
In this section, various steps needed to setup the turbulent flow over an airfoil will be shown. A full working case of this can be found in:
/tutorials/incompressible/simpleSolver/ras/ACCM_airFoil2D
However,the user is free to start the case setup from scratch consistent with the directory stucture discussed below.
Directory Structure
Note
All commands shown here are entered in a terminal window, unless otherwise mentioned
The problem requires time
, constant
and system
subdirectories within the main working directory. Here, the simulation will be started at time t = 0 s, which requires a time subdirectory named 0
.
The 0
subdirectory has files in which boundary properties are specified. In the below table, the list of necessary files are provided based on the turbulence model chosen
Parameter 
File name 

Pressure (\(p\)) 

Velocity (\(U\)) 

Turbulent viscosity (\(\nu\)) 

Turbulence field variable (\(\tilde{\nu}\)) 

Turbulent kinetic energy (\(k\)) 

Turbulent dissipation rate (\(\omega\)) 

We will consider two turbulence models in this tutorial, namely SpalartAllmaras (SA) and \(k\omega\)  Shear Stress Transport (\(\rm{SST}\)). The contents of the files listed above sets the dimensions, initialisation and boundary conditions to the defining problem, which also forms three principle entries required.
The user should note that Caelus is case sensitive and therefore the directory and file setup should be identical to what is shown here.
Boundary Conditions and Solver Attributes
Boundary Conditions
Referring back to Figure 37, the following are the boundary conditions that will be specified:
 Inlet
Velocity:
\(\alpha=0^\circ\): Fixed uniform velocity \(u = 52.0770~m/s\); \(v = w = 0.0~m/s\) in \(x, y\) and \(z\) directions respectively
\(\alpha=10^\circ\): Fixed uniform velocity \(u = 51.2858~m/s\); \(v = 0.0~m/s\) and \(w = 9.04307~m/s\) in \(x, y\) and \(z\) directions respectively
Pressure: Zero gradient
Turbulence:
Spalart–Allmaras (Fixed uniform values of \(\nu_{t~\infty}\) and \(\tilde{\nu}_{\infty}\) as given in Turbulent freestream conditions for SA Model)
\(k\omega~\rm{SST}\) (Fixed uniform values of \(k_{\infty}\), \(\omega_{\infty}\) and \(\nu_{t~\infty}\) as given in Turbulent freestream conditions for SST Model)
 Noslip wall
Velocity: Fixed uniform velocity \(u, v, w = 0\)
Pressure: Zero gradient
Turbulence:
Spalart–Allmaras (Fixed unifSpalart–Allmaras (Fixed uniform values of \(\nu_{t}=0\) and \(\tilde{\nu}=0\))
\(k\omega~\rm{SST}\) (Fixed uniform values of \(k = <<0\) and \(\nu_t=0\); \(\omega\) = omegaWallFunction)
 Outlet
Velocity: Zero gradient velocity
Pressure: Fixed uniform gauge pressure \(p = 0\)
Turbulence:
Spalart–Allmaras (Calculated \(\nu_{t}=0\) and Zero gradient \(\tilde{\nu}\))
\(k\omega~\rm{SST}\) (Zero gradient \(k\) and \(\omega\); Calculated \(\nu_t=0\); )
 Initialisation
Velocity:
\(\alpha=0^\circ\): Fixed uniform velocity \(u = 52.0770~m/s\); \(v = w = 0.0~m/s\) in \(x, y\) and \(z\) directions respectively
\(\alpha=10^\circ\): Fixed uniform velocity \(u = 51.2858~m/s\); \(v = 0.0~m/s\) and \(w = 9.04307~m/s\) in \(x, y\) and \(z\) directions
Pressure: Zero Gauge pressure
Turbulence:
Spalart–Allmaras (Fixed uniform values of \(\nu_{t~\infty}\) and \(\tilde{\nu}_{\infty}\) as given in Turbulent freestream conditions for SA Model)
\(k\omega~\rm{SST}\) (Fixed uniform values of \(k_{\infty}\), \(\omega_{\infty}\) and \(\nu_{t~\infty}\) as given in Turbulent freestream conditions for SST Model)
First, the pressure file named p
has the following contents
/**\
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 2 0 0 0 0];
internalField uniform 0;
boundaryField
{
inlet
{
type zeroGradient;
}
left
{
type empty;
}
outlet
{
type fixedValue;
value uniform 0;
}
right
{
type empty;
}
wall
{
type zeroGradient;
}
}
// ************************************************************************* //
In the information shown above, it can be seen that the file begins with a dictionary named FoamFile
which contains the standard set of keywords for version, format, location, class and object names. The explanation of the principle entries are as follows
dimension
is used to specify the physical dimensions of the pressure field. Here, pressure is defined in terms of kinematic pressure with the units (\(m^2/s^2\)) written as
[0 2 2 0 0 0 0]
internalField
is used to specify the initial conditions. It can be either uniform or nonuniform. Since we have a 0 initial uniform gauge pressure, the entry is
uniform 0;
boundaryField
is used to specify the boundary conditions. In this case its the boundary conditions for pressure at all the boundary patches.
Similarly, the file U
is defined as follows,
/**\
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
location "0";
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 1 0 0 0 0];
internalField uniform (52.077 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (52.077 0 0);
}
left
{
type empty;
}
outlet
{
type zeroGradient;
}
right
{
type empty;
}
wall
{
type noSlipWall;
}
}
// ************************************************************************* //
The principle entries for velocity field are selfexplanatory and the dimensions are typically for that of velocity with the units \(m/s\) ([0 1 1 0 0 0 0]
). The user should note that appropriate entry has to be enabled for both internalField
and inlet
boundaryField
depending on which angle of attack (AOA/\(\alpha\)) is being simulated. Here, both initialisation and inlet have a uniform flow velocity specified with three velocity components. For example at \(\alpha = 10^\circ\), we specify (51.2858 0 9.04307)
for \(u, v, w\) components respectively. Similarly, the wall boundary patch have three velocity components.
The turbulent properties are also required to be specified at the boundary patches and these can be done similar to p
and U
. The file nut
, defines turbulent kinematic viscosity and its boundary conditions as follows.
/**\
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object nut;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 1 0 0 0 0];
internalField uniform 1.8265e06;
boundaryField
{
inlet
{
type fixedValue;
value uniform 1.8265e06;
}
left
{
type empty;
}
outlet
{
type calculated;
value uniform 0;
}
right
{
type empty;
}
wall
{
type fixedValue;
value uniform 0;
}
}
// ************************************************************************* //
As noted above, the turbulent viscosity is specified as kinematic and therefore the units are in \(m^2/s\) ([0 2 1 0 0 0 0]
). The turbulent viscosity value at freestream, specified at inlet patch is calculated as detailed in Turbulent freestream conditions for SA Model and Turbulent freestream conditions for SST Model for SA and SST models respectively and is specified accordingly. The same value also goes for internalField
. Note that a fixedValue
of 0
is used for the wall
which suggests that on the wall, it is only the molecular (laminar) viscosity that prevails.
The next turbulent property to set is the nuTilda
which is a turbulent field variable, specified to only SA model and has the same units ([0 2 1 0 0 0 0]
) as kinematic turbulent viscosity. Details pertaining to this are given in Turbulent freestream conditions for SA Model. The following are the contents of the file nuTilda
and the entries specified for the boundaryField
are identical to that of turbulent kinematic viscosity explained above.
/**\
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object nuTilda;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 1 0 0 0 0];
internalField uniform 2.60385e05;
boundaryField
{
inlet
{
type fixedValue;
value uniform 2.60385e05;
}
left
{
type empty;
}
outlet
{
type zeroGradient;
}
right
{
type empty;
}
wall
{
type fixedValue;
value uniform 0;
}
}
// ************************************************************************* //
We can now proceed to the files k
and omega
, specific to only \(k\omega~\rm{SST}\) model. \(k\omega~\rm{SST}\) is a turbulence model which solves for turbulent kinetic energy and the specific rate of dissipation using two partial differential equations. Caelus therefore requires information about these to be specified at the boundary patches when this model is chosen. Starting with the file k
,
/**\
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object k;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 2 0 0 0 0];
internalField uniform 0.0010999;
boundaryField
{
inlet
{
type fixedValue;
value uniform 0.0010999;
}
left
{
type empty;
}
outlet
{
type zeroGradient;
}
right
{
type empty;
}
wall
{
type fixedValue;
value uniform 1e10;
}
}
// ************************************************************************* //
The unit of kinetic energy is \(m^2/s^2\) and this is set in dimensions
as [0 2 2 0 0 0 0]
. As with other turbulent quantities discussed above, the value of \(k\) (refer Turbulent freestream conditions for SST Model) needs to be specified for internalField
, inlet
and wall
. Please note that for wall boundaryField
with no wallfunction, a small, nonzero fixedValue
is required.
We now proceed to the file omega
and the value for this is evaluated as detailed in Turbulent freestream conditions for SST Model
/**\
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object omega;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 0 1 0 0 0 0];
internalField uniform 13887.21912;
boundaryField
{
inlet
{
type fixedValue;
value uniform 13887.21912;
}
left
{
type empty;
}
outlet
{
type zeroGradient;
}
right
{
type empty;
}
wall
{
type omegaWallFunction;
value uniform 1;
}
}
// ************************************************************************* //
The unit of specific rate of dissipation for \(\omega\) is \(1/s\) which is set in dimensions
as [0 0 1 0 0 0 0]
. The internalField
and inlet
gets a fixedValue
. Note that for wall boundaryField
, we specify omegaWallFunction
and this is a model requirement and sets omega to the correct value near wall based on the \(y^+\). In conjunction, the value that goes with omegaWallFunction
can be anything and for simplicity its set to 1
.
Before proceeding to the next step, it is vital to ensure that the boundary conditions (inlet
, outlet
, wall
, etc) added in the above files should be the grid boundary patches (surfaces) generated by grid generation tool and their names are identical. Additionally, the two boundaries \(xz\) plane named here as left
and right
have empty
boundary conditions which forces Caelus to assume the flow to be in twodimensions. With this, the setting up of the boundary conditions are complete.
Grid file and Physical Properties
The grid files associated with the NACA airfoil need to be placed in polyMesh
subdirectory, which resides in the constant
directory. Note that for both angles of attack, the identical grid is used. This is because the flow incidence angle is relative to the fixed airfoil and the equivalent velocity components can be easily specified thus simulating the airfoil at the required angle of attack. In addition, the physical properties are specified in various different files present in the constant
directory.
polyMesh RASProperties transportProperties turbulenceProperties
As seen above, the three files are listed in addition to the polyMesh
subdirectory. The first one, RASProperties
in which the ReynoldsAverageStress (RAS) model is specified, which is shown below. Please note that depending on the turbulence model you wish to run with, the line that corresponds to that specific model should be enabled.
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object RASProperties;
}
//
// For SpalaratAlamaras Model, enable the below line
RASModel SpalartAllmaras;
// For komega SST Model, enable the below line
// RASModel kOmegaSST;
turbulence on;
printCoeffs on;
SpalartAllmarasCoeffs
{
curvatureCorrection off;
}
kOmegaSSTCoeffs
{
curvatureCorrection off;
Cscale 1;
frMax 1.25;
}
Second from the list is the transportProperties
file, where kinematic viscosity is specified as shown below
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
//
transportModel Newtonian;
nu nu [0 2 1 0 0 0 0] 8.679514E6;
The viscous behaviour is modelled as Newtonian and hence the keyword Newtonian
is used for the transportModel
and the molecular (laminar) kinematic viscosity (\(nu\)) is given having the units \(m^2/s\) ([0 2 1 0 0 0 0]
).
The final file in this class is the turbulenceProperties
file, which sets the simulationType
to RASModel
. Both SA and \(k\omega~\rm{SST}\) are classified as Reynolds Average Stress (RAS) models.
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object turbulenceProperties;
}
//
simulationType RASModel;
Controls and Solver Attributes
In this section, the files required to control the simulation, specifying discretisation methods and linear solver settings are given. These files are to be placed in the system
directory. First, navigate to the system
directory.
The controlDict
file contains the following details
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
//
application simpleSolver;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 10000; // increase endTime to atleast 50000 for validation cases
deltaT 1;
writeControl runTime;
writeInterval 1000;
purgeWrite 2;
writeFormat ascii;
writePrecision 12;
writeCompression true;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
//
// Function Objects to obtain forces
functions
{
forces
{
type forces;
functionObjectLibs ("libforces.so");
patches ( wall );
CofR (0 0 0);
rhoName rhoInf;
rhoInf 2.15527;
writeControl timeStep;
writeInterval 50;
}
}
Referring to the above file, some explanation is required. Here, simpleSolver
is used and the simulation begins at t = 0 s
. This now explains the logical need for having a 0
directory where the data files are read at the beginning of the run, which is t = 0 s
in this case. Therefore, the keyword startFrom
is set to startTime
, where startTime
is set to 0
. The simulation would be carried out as steadystate and therefore we require to specify the total number of iterations as a keyword for endTime
. Via the writeControl
and writeInterval
keywords, the solution intervals at which they are saved can be specified. Also note that a function object to obtain the force over the wall (airfoil surface) for every 50 iterations is included. In order to obtain this, a freestream density (rhoInf
) need to be specified.
The discretization schemes for the finite volume discretization that will be used should be set through the fvSchemes
file shown below.
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object fvSchemes;
}
//
ddtSchemes
{
default steadyState;
}
gradSchemes
{
default Gauss linear;
}
divSchemes
{
default none;
div(phi,U) Gauss linearUpwindBJ grad(U);
div(phi,nuTilda) Gauss upwind;
div(phi,k) Gauss upwind;
div(phi,omega) Gauss upwind;
div((nuEff*dev(T(grad(U))))) Gauss linear;
div(phi,symm(grad(U))) Gauss linear;
}
laplacianSchemes
{
default none;
laplacian(nu,U) Gauss linear corrected;
laplacian(nuEff,U) Gauss linear corrected;
laplacian(DnuTildaEff,nuTilda) Gauss linear corrected;
laplacian(DkEff,k) Gauss linear corrected;
laplacian(DomegaEff,omega) Gauss linear corrected;
laplacian(rAUf,p) Gauss linear corrected;
laplacian(1,p) Gauss linear corrected;
}
interpolationSchemes
{
default linear;
interpolate(HbyA) linear;
}
snGradschemes
{
default corrected;
}
The linear solver controls and tolerances are set in fvSolution
as given below
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSolution;
}
//
solvers
{
p
{
solver PCG;
preconditioner
{
preconditioner GAMG;
smoother SSGS;
cacheAgglomeration true;
nCellsInCoarsestLevel 10;
agglomerator algebraicPair;
mergeLevels 1;
}
tolerance 1e16;
relTol 0.01;
}
U
{
solver PBiCGStab;
preconditioner USGS;
tolerance 1e16;
relTol 0.1;
}
"(komeganuTilda)"
{
solver PBiCGStab;
preconditioner USGS;
tolerance 1e16;
relTol 0.1;
}
}
SIMPLE
{
nNonOrthogonalCorrectors 1;
pRefCell 0;
pRefValue 0;
consistent yes;
correctionFoam yes;
}
// relexation factors
relaxationFactors
{
p 0.3;
U 0.7;
nuTilda 0.7;
k 0.7;
omega 0.7;
}
The user should now take a note that in the fvSolution
file, different linear solvers are used to solve for velocity, pressure and turbulence quantities We also set the nNonOrthogonalCorrectors
to 1 for this case. To ensure the stability of the solution, the relaxation is set to primary and turbulent variables. The relTol
is not set to 0
unlike a timeaccurate setup as we are solving for a steadystate solution and a very low (\(\approx 0\)) tolerance at every iteration is unnecessary. Since the entire system of equations converge to a global set tolerance
the convergence would occur as the solution progresses to a steady state.
With this, the setup of the directory structure with all the relevant files are complete. This can be verified again by issuing the following command and the directory tree should appear identical to the one shown below
cd myturbulentairfoil/
tree
.
├── 0
│ ├── epsilon
│ ├── k
│ ├── nut
│ ├── nuTilda
│ ├── omega
│ ├── p
│ └── U
├── constant
│ ├── polyMesh
│ │ ├── boundary
│ │ ├── faces
│ │ ├── neighbour
│ │ ├── owner
│ │ └── points
│ ├── RASProperties
│ ├── transportProperties
│ └── turbulenceProperties
└── system
├── controlDict
├── decomposePar
├── fvSchemes
└── fvSolution
Execution of the solver¶
Before executing the solver, it is important to renumber and to carry out a quality check on the grid/mesh. Renumbering reduces the bandwidth whereas the quality check shows the mesh statistics. These two can be performed by executing the following commands from the top working directory
caelus run  renumberMesh overwrite
caelus run  checkMesh
When the renumberMesh
is performed, the user should take note of the bandwidth before and after the mesh renumbering. In a similar manner, when the checkMesh
is performed, the mesh statistics are shown as below
/**\
Caelus 8.04
Web: www.caeluscml.com
\**/
Checking geometry...
Overall domain bounding box (484.456616748 0 507.806809185) (501.000007802 10 507.806808887)
Mesh (nonempty, nonwedge) directions (1 0 1)
Mesh (nonempty) directions (1 0 1)
All edges aligned with or perpendicular to nonempty directions.
Boundary openness (8.28364670182e18 1.55457412622e15 1.3301799376e17) OK.
***High aspect ratio cells found, Max aspect ratio: 31784747.6728, number of cells 29164
<<Writing 29164 cells with high aspect ratio to set highAspectRatioCells
Minimum face area = 9.57960157047e11. Maximum face area = 843.777295738. Face area magnitudes OK.
Min volume = 9.57960157047e10. Max volume = 8437.77295738. Total volume = 8755584.35466. Cell volumes OK.
Mesh nonorthogonality Max: 57.3127287926 average: 1.73980453241
Nonorthogonality check OK.
Face pyramids OK.
Mesh skewness Max: 0.201755632336 average: 0.00254964789178 OK.
Coupled point location match (average 0) OK.
Failed 1 mesh checks.
End
As can be noted above, the output of the checkMesh
indicates that the mesh check has failed reporting in the final message as Failed 1 mesh checks
. This is because of the high aspect ratio meshes present immediate to the wall with very low (\(<< y^+\)) values. Nevertheless, this is just a warning and Caelus will solve on this mesh.
As with the previous tutorial, it will be shown here to utilise the multicore capability of CPUs for performing a parallel computation using MPI technique for large grids, such as the one considered here. Before this can be done, the grid has to be decomposed into smaller domains that can be solved by each single CPU core. The number of decomposition should be equal to the number of CPU core available for parallel computing. The decomposition is carried out through a file decomposeParDict
present in the system
subdirectory as shown below.
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object decomposeParDict;
}
//
numberOfSubdomains 4; // It is suggested that the numberOfSubdomains be increased based on available resources for validation cases and to reduce the computation time.
method scotch;
distributed no;
roots
(
);
In the above file, the the keyword numberOfSubdomains
defines the number of decomposed subdomains needed and 4
is used which partitions the grid into 4 subdomains. We use scotch
as the method of decomposition which automatically divides the gird. The execution to decompose the grid is carried out again from the top directory as follows
caelus run  decomposePar
Now the decomposition should begin and the details of which are displayed in the terminal window. Subsequently, 4 processor directories will be generated as we have requested for 4 divisions of grid as shown below
0 constant processor0 processor1 processor2 processor3 system
The solver can now be executed for parallel computation in the host machine from the top directory using the following command
caelus run p l myturbulentairfoil.log  simpleSolver
Note that here it is assumed that the parallel computing is available in the host machine. With the execution of the above commands, the simulation begins and the progress of the solution is written to the specified log file (myturbulentairfoil.log
).
The log file can be further processed to look at the convergence history and this can be done as follows
caelus logs w myturbulentairfoil.log
The above command allows you to look at the convergence of different variables with respect to the number of iterations carried out. The convergence of pressure is shown in Figure 39.
Results¶
The flow over the airfoil at both \(0^\circ\) and \(10^\circ\) degree angle of attack are presented here. In Figure 40, velocity magnitude and pressure contours can be seen for \(\alpha = 0^\circ\) angle of attack. These result are for the SA model. The suction and the pressure surfaces essentially produce the same flow due to \(0^\circ\) angle of incidence and thus contributes to zero lift. In contrast, at \(0^\circ\) angle of incidence in Figure 41, a low pressure region exists on the upper surface and consequently the velocity increases thus generating some lift.
Convex Curvature¶
Turbulent flow in a constantarea duct having a convex curvature is considered in this tutorial. Caelus 9.04 will be used and the process of settingup of directory structure, fluid properties, boundary conditions, turbulent properties, etc will be explained here. In addition to this, the flow within the duct will be visualised. The steps would be sufficient for the user to reproduce the tutorial accurately.
Objectives¶
Through this tutorial, the user will get familiarise with setting up the Caelus simulation for steady, turbulent, incompressible flow in a twodimensional duct having a convex curvature and subsequently postprocess the results. The steps that would be followed in this tutorial is outlined below
 Background
A brief description about the problem
Geometry and inflow details
 Grid generation
Computational domain and boundary details
Computational grid generation
Exporting grid to Caelus
 Problem definition
Directory structure
Setting up boundary conditions, physical properties and control/solver attributes
 Execution of the solver
Monitoring the convergence
Writing the log files
 Results
Visualising the flow within the convex duct
Prerequisites¶
It will be assumed that the user is comfortable and familiar with the Linux command line environment using a terminal or Caelusconsole (on Windows OS). The user should also make sure that Caelus is installed correctly with appropriate environment variables set. The grid used here is obtained from Turbulence Modeling Resource in Plot 3D format and is exported to Caelus using Pointwise. The user can use their choice of grid generation tool to convert Plot3D file to Caelus format.
Background¶
Turbulent flow in a constantarea duct with a curvature is an interesting case. Here as a result of a curvature, pressure gradients occur in the vicinity of the curvature having localised effect. The flow will be assumed as steadystate and incompressible. Nondimensional shearstress (skinfriction coefficient) will be used to show the influence of curvature on the flow. Validation and verification of this exercise is detailed in section Twodimensional Convex Curvature and the user is suggested to refer for more information.
The inlet of the duct as can be seen from the schematic below in Figure 42 has an inclination of \(\alpha = 30^\circ\). This is followed by a rapid bend at the same angle, \(\alpha\) after a distance of about 1.4 m. The downstream extends to 1.6 meters. The inflow has a Reynolds number of \(2.1 \times 10^6\), with Air as the fluid. The temperature of the inflow is at 293 K and U is the inlet velocity. Based on the Reynolds number, temperature and velocity, the kinematic viscosity is evaluated to \(\nu = 1.519470 \times 10^{5}~m^2/s\). The geometricplane for this case in 2D is \(xz\) plane.
The inflow conditions for this case is given in the below table (Freestream conditions)
Fluid 
\(Re/L~(1/m)\) 
\(U~(m/s)\) 
\(p~(m^2/s^2)\) 
\(T~(K)\) 
\(\nu~(m^2/s)\) 

Air 
\(2.1 \times 10^6\) 
31.9088 
Gauge (0) 
293 
\(1.5194\times10^{6}\) 
To achieve a inlet velocity of U = 31.9088 m/s at \(\alpha = 30^\circ\), the velocity components in \(x\) and \(z\) are resolved. These are given in the below table
\(\alpha\) 
\(u~(m/s)\) 
\(w~(m/s)\) 

\(30^\circ\) 
27.63313 
15.95443 
Grid Generation¶
The structured grid for this case has been obtained from from Turbulence Modeling Resource in a Plot3 format which contains 512 cells in streamwise direction and 192 in the flow normal direction. The Plot3D file has to be converted to polyMesh
format.
In Figure 43, the computational domain can be seen which as expected follows the geometry. The velocity at the internal walls (highlighted in blue) are zero, wherein \(u, v, w =0\) representing a noslip boundary. The inlet and the outlet are applied to the start and end of the domain respectively.
The grid that is used, polyMesh
format is in threedimensions. However it is assumed that the flow through the duct can be modelled as 2D and is the 2D plane considered for this case is \(xz\) directions. As with all the previous cases, onecell thick normal to the 2D plane is sufficient, assuming the flow to be symmetry. The two \(xz\) planes that are present in the grid require boundary conditions to be specified. An empty boundary condition can be used in Caelus for the two planes that forces the solver not to solve in the thirddimension, essentially treating the flow as symmetry in \(y\) direction.
Note
A velocity value of \(v=0\) needs to be specified at appropriate boundaries although no flow is solved in the \(y\) direction.
The 2D grid in \(xz\) plane is shown in Figure 44 having a distribution of 512 X 192 cells. The inset in the figure highlights the region vicinity of the curvature and very fine distribution of cells can be seen close to the wall. It is estimated that \(y^+\) is less than 1 in order to capture turbulent boundary layer accurately and thus no wallfunction is used.
Problem definition¶
This section details the various steps required to setup the turbulent flow inside a convex duct. A full working case of this can be found in:
/tutorials/incompressible/simpleSolver/ras/ACCM_convexCurvature2D
Directory Structure
Note
All commands shown here are entered in a terminal window, unless otherwise mentioned
A few more subdirectories are needed in the toplevel directory to setup the case. Caelus requires time
, constant
and system
subdirectories within the main myturbulentcurvature
working directory. Since we start the simulation at time, t =0 s, a time
subdirectory named 0
is required. A
The 0
subdirectory cotains the files in which boundary properties are specified. In the below table, the list of necessary files are provided based on the turbulence model chosen
Parameter 
File name 

Pressure (\(p\)) 

Velocity (\(U\)) 

Turbulent viscosity (\(\nu\)) 

Turbulence field variable (\(\tilde{\nu}\)) 

Turbulent kinetic energy (\(k\)) 

Turbulent dissipation rate (\(\omega\)) 

We consider simulating this case with three turbulence models, namely SpalartAllmaras (SA), Spalart–Allmaras Rotational/Curvature (SARC) and \(k\omega\)  Shear Stress Transport (\(\rm{SST}\)). The contents that are entered in these files set the dimensions, initialisation and boundary conditions to the defining problem, which also form three principle entries required.
Caelus is case sensitive and therefore the user should carefully setup the case as shown here if applicable.
Boundary Conditions and Solver Attributes
Boundary Conditions
Referring back to figure Computational domain of the convex curvature, following are the boundary conditions that will be specified:
 Inlet
Velocity: Fixed uniform velocity \(u = 27.63313~m/s\); \(v = 0.0~m/s\) and \(w = 15.95443~m/s\) in \(x, y\) and \(z\) directions respectively
Pressure: Zero gradient
Turbulence:
SA & SARC (Fixed uniform values of \(\nu_{t~i}\) and \(\tilde{\nu}_{i}\) as given in Turbulent freestream conditions for SA Model)
\(k\omega~\rm{SST}\) (Fixed uniform values of \(k_{i}\), \(\omega_{i}\) and \(\nu_{t~i}\) as given in Turbulent freestream conditions for SST Model)
 Noslip wall
Velocity: Fixed uniform velocity \(u, v, w = 0\)
Pressure: Zero gradient
Turbulence:
SA & SARC (Fixed uniform values of \(\nu_{t}=0\) and \(\tilde{\nu} =0\))
\(k\omega~\rm{SST}\) (Fixed uniform values of \(k = <<0\) and \(\nu_t=0\); \(\omega\) = omegaWallFunction)
 Outlet
Velocity: Zero gradient velocity
Pressure: Fixed uniform gauge pressure \(p = 0\)
Turbulence:
SA & SARC (Calculated \(\nu_{t}=0\) and Zero gradient \(\tilde{\nu}\))
\(k\omega~\rm{SST}\) (Zero gradient \(k\) and \(\omega\); Calculated \(\nu_t=0\); )
 Initialisation
Velocity: Fixed uniform velocity \(u = 27.63389~m/s\); \(v = 0.0~m/s\) and \(w = 15.95443~m/s\) in \(x, y\) and \(z\) directions respectively
Pressure: Zero Gauge pressure
Turbulence:
SA & SARC (Fixed uniform values of \(\nu_{t~i}\) and \(\tilde{\nu}_{i}\) as given in Turbulent freestream conditions for SA Model)
\(k\omega~\rm{SST}\) (Fixed uniform values of \(k_{i}\), \(\omega_{i}\) and \(\nu_{t~i}\) as given in the above table) as given in Turbulent freestream conditions for SST Model)
First, we begin with the pressure file, p
and using a text editor with the following content.
/**\
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 2 0 0 0 0];
internalField uniform 0;
boundaryField
{
inlet
{
type zeroGradient;
}
left
{
type empty;
}
outlet
{
type fixedValue;
value uniform 0;
}
right
{
type empty;
}
top
{
type zeroGradient;
}
wall
{
type zeroGradient;
}
}
// ************************************************************************* //
As noted above, the file begins with a dictionary named FoamFile
containing the essential set of keywords for version, format, location, class and object names. The explanation of the principle entries are as follows
dimension
is used to specify the physical dimensions of the pressure field. Here, pressure is defined in terms of kinematic pressure with the units (\(m^2/s^2\)) written as
[0 2 2 0 0 0 0]
internalField
is used to specify the initial conditions. It can be either uniform or nonuniform. Since we have a 0 initial uniform gauge pressure, the entry is
uniform 0;
boundaryField
is used to specify the boundary conditions. In this case its the boundary conditions for pressure at all the boundary patches.
Similarly, the contents for the file U
is as follows
/**\
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
location "0";
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 1 0 0 0 0];
internalField uniform (27.63313 0 15.954);
boundaryField
{
inlet
{
type fixedValue;
value uniform (27.63313 0 15.954);
}
left
{
type empty;
}
outlet
{
type zeroGradient;
}
right
{
type empty;
}
top
{
type fixedValue;
value uniform (0 0 0);
}
wall
{
type noSlipWall;
}
}
// ************************************************************************* //
The principle entries for velocity field are selfexplanatory and the dimensions are typically for that of velocity with the units \(m/s\) ([0 1 1 0 0 0 0]
). Here, both initialisation and inlet have a uniform flow velocity specified with three velocity components. For example at \(\alpha = 30^\circ\), we specify (27.63313 0 15.954)
for \(u, v, w\) components respectively. Similarly, the top and wall boundary patch have three velocity components.
The turbulent properties are also required to be specified at the boundary patches and these can be done similar to p
and U
. First, we start with opening the file nut
, which is turbulent kinematic viscosity having the following details
/**\
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object nut;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 1 0 0 0 0];
internalField uniform 1.36756e07;
boundaryField
{
inlet
{
type fixedValue;
value uniform 1.36756e07;
}
left
{
type empty;
}
outlet
{
type calculated;
value uniform 0;
}
right
{
type empty;
}
top
{
type fixedValue;
value uniform 0;
}
wall
{
type fixedValue;
value uniform 0;
}
}
// ************************************************************************* //
In the above file, the turbulent viscosity is specified as kinematic and therefore the units are in \(m^2/s\) ([0 2 1 0 0 0 0]
). The turbulent viscosity value at inlet, specified at inlet patch is calculated as detailed in Turbulent freestream conditions for SA Model and Turbulent freestream conditions for SST Model for SST models respectively and is specified accordingly. The same value also goes for internalField
. Note that a fixedValue
of 0
is used for the wall
which suggests that on the wall, it is only the molecular (laminar) viscosity that prevails.
The next property to set is the nuTilda
, which is the turbulent field variable, specified to only SA and SARC models and has the same units ([0 2 1 0 0 0 0]
) as kinematic turbulent viscosity. For more information about these, the user can look at Turbulent freestream conditions for SA Model. Following are the contents of the file nuTilda
and the entries specified for the boundaryField
are identical to that of turbulent kinematic viscosity explained above.
/**\
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object nuTilda;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 1 0 0 0 0];
internalField uniform 4.55841e05;
boundaryField
{
inlet
{
type fixedValue;
value uniform 4.55841e05;
}
left
{
type empty;
}
outlet
{
type zeroGradient;
}
right
{
type empty;
}
top
{
type fixedValue;
value uniform 0;
}
wall
{
type fixedValue;
value uniform 0;
}
}
// ************************************************************************* //
We can now set the properties of k
and omega
, specific to only \(k\omega~\rm{SST}\) model. \(k\omega~\rm{SST}\) is a turbulence model which solves for turbulent kinetic energy and the specific rate of dissipation using two partial differential equations. Caelus therefore requires information about these to be specified at the boundary patches when this model is chosen. Starting with the file k
,
/**\
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object k;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 2 0 0 0 0];
internalField uniform 0.001052132;
boundaryField
{
inlet
{
type fixedValue;
value uniform 0.001052132;
}
left
{
type empty;
}
outlet
{
type zeroGradient;
}
right
{
type empty;
}
top
{
type fixedValue;
value uniform 1e10;
}
wall
{
type fixedValue;
value uniform 1e10;
}
}
// ************************************************************************* //
The unit of kinetic energy is \(m^2/s^2\) and this is set in dimensions
as [0 2 2 0 0 0 0]
. As with other turbulent quantities discussed above, the value of \(k\) (refer Turbulent freestream conditions for SST Model) needs to be specified for internalField
, inlet
and wall
. Please note that for wall boundaryField
with no wallfunction, a small, nonzero fixedValue
is required.
Similarly in the file omega
, the value is evaluated as detailed in Turbulent freestream conditions for SST Model
/**\
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object omega;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 0 1 0 0 0 0];
internalField uniform 7747.333;
boundaryField
{
inlet
{
type fixedValue;
value uniform 7747.333;
}
left
{
type empty;
}
outlet
{
type zeroGradient;
}
right
{
type empty;
}
top
{
type omegaWallFunction;
value uniform 1;
}
wall
{
type omegaWallFunction;
value uniform 1;
}
}
// ************************************************************************* //
The unit of specific rate of dissipation for \(\omega\) is \(1/s\) which is set in dimensions
as [0 0 1 0 0 0 0]
. The internalField
and inlet
gets a fixedValue
. Note that for wall boundaryField
, we specify omegaWallFunction
and this is a model requirement and sets omega to the correct value near wall based on the \(y^+\). In conjunction, the value that goes with omegaWallFunction
can be anything and for simplicity its set to 1
.
At this stage it is important to ensure that the boundary conditions (inlet
, outlet
, wall
, etc) added in the above files should be the grid boundary patches (surfaces) generated by grid generation tool and their names are identical. In addition, the two boundaries \(xz\) plane named here as left
and right
have empty
boundary conditions which forces Caelus to assume the flow to be in twodimensions. With this, the setting up of the boundary conditions are complete.
Grid file and Physical Properties
The grid file for the convex curvature has to be placed in polyMesh
subdirectory, which is in the constant
directory. Further to is, the physical properties should be specified in various different files present in the constant
directory.
polyMesh RASProperties transportProperties turbulenceProperties
As noted above, the three files are listed in the constant
subdirectory. The first one is the RASProperties
where, ReynoldsAverageStress (RAS) model is specified as shown below. Please note that depending on the turbulence model you wish to run with, the line that corresponds to that specific model should be enabled.
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object RASProperties;
}
//
// For SpalaratAlamaras Model, enable the below line
//RASModel SpalartAllmaras;
// For komega SST Model, enable the below line
RASModel kOmegaSST;
turbulence on;
printCoeffs on;
SpalartAllmarasCoeffs
{
curvatureCorrection on;
}
kOmegaSSTCoeffs
{
curvatureCorrection off;
Cscale 1;
frMax 1.25;
}
Second from the list is the transportProperties
file, where kinematic viscosity is specified as shown below
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
//
transportModel Newtonian;
nu nu [0 2 1 0 0 0 0] 1.519470E5;
The viscous behaviour is modelled as Newtonian and hence the keyword Newtonian
is used for the transportModel
and the molecular (laminar) kinematic viscosity (\(nu\)) is given having the units \(m^2/s\) ([0 2 1 0 0 0 0]
).
The final file in this class is the turbulenceProperties
file, which sets the simulationType
to RASModel
. Both SA and \(k\omega~\rm{SST}\) are classified as Reynolds Average Stress (RAS) models.
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object turbulenceProperties;
}
//
simulationType RASModel;
Controls and Solver Attributes
This section will provide details and settings required to control the simulation, specifying discretisation methods and linear solver settings. These files should be saved in the system
directory.
The controlDict
file contains the following details
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
//
application simpleSolver;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 3000;
deltaT 1;
writeControl runTime;
writeInterval 1000;
purgeWrite 0;
writeFormat ascii;
writePrecision 12;
writeCompression true;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
//
// Function Objects to obtain forces
functions
{
forces
{
type forces;
functionObjectLibs ("libforces.so");
patches ( wall );
CofR (0 0 0);
rhoName rhoInf;
rhoInf 1.2084;
writeControl timeStep;
writeInterval 50;
}
}
With reference to the above files, some explanation is required. In this case, simpleSolver
solver is used and the simulation begins at t = 0 s
. This now explains the logical need for having a 0
directory where the data files are read at the beginning of the run, which is t = 0 s
for this simulation. Therefore, the keyword startFrom
is set to startTime
, where startTime
is set to 0
. The simulation would be carried out as steadystate and therefore we require to specify the total number of iterations as a keyword for endTime
. Via the writeControl
and writeInterval
keywords, the solution intervals at which they are saved can be specified. Also note that a function object to obtain the force over the wall for every 50 iterations is included. In order to obtain this, a inlet/inflow density (rhoInf) need to be specified.
The discretization schemes for the finite volume discretization that will be used should be set through the fvSchemes
file show below and the contents should be copied
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object fvSchemes;
}
//
ddtSchemes
{
default steadyState;
}
gradSchemes
{
default Gauss linear;
grad(p) Gauss linear;
grad(U) Gauss linear;
}
divSchemes
{
default none;
div(phi,U) Gauss linearUpwind grad(U);
div(phi,nuTilda) Gauss upwind; // Will be used for SA & SARC model only
div(phi,k) Gauss upwind; // will be used for kepsilon & komega only
div(phi,omega) Gauss upwind; // Will be used for komega model only
div((nuEff*dev(T(grad(U))))) Gauss linear;
div(phi,symm(grad(U))) Gauss linear;
}
laplacianSchemes
{
default none;
laplacian(nu,U) Gauss linear corrected;
laplacian(nuEff,U) Gauss linear corrected;
laplacian(DnuTildaEff,nuTilda) Gauss linear corrected; // Will be used for SA & SARC model only
laplacian(DkEff,k) Gauss linear corrected; // Will be used for komega & komega only
laplacian(DomegaEff,omega) Gauss linear corrected; // Will be used for komega model only
laplacian(rAUf,p) Gauss linear corrected;
laplacian(1,p) Gauss linear corrected;
}
interpolationSchemes
{
default linear;
interpolate(HbyA) linear;
}
snGradschemes
{
default corrected;
}
The linear solver controls and tolerances are set in fvSolution
as given below
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSolution;
}
//
solvers
{
p
{
solver PCG;
preconditioner SSGS;
tolerance 1e8;
relTol 0.01;
}
U
{
solver PBiCGStab;
preconditioner USGS;
tolerance 1e7;
relTol 0.01;
}
"(komeganuTilda)"
{
solver PBiCGStab;
preconditioner USGS;
tolerance 1e08;
relTol 0;
}
}
SIMPLE
{
nNonOrthogonalCorrectors 1;
pRefCell 0;
pRefValue 0;
}
// relexation factors
relaxationFactors
{
p 0.3;
U 0.5;
nuTilda 0.5;
k 0.5;
omega 0.5;
}
In the fvSolution
file, different linear solvers are used to solve for velocity, pressure and turbulence quantities and this has to be noted by the user. We also set the nNonOrthogonalCorrectors
to 1 for this case. To ensure the stability of the solution, the relaxation is set to primary and turbulent variables. The relTol
is not set to 0
unlike a timeaccurate setup as we are solving for a steadystate solution and a very low (\(\approx 0\)) tolerance at every iteration is unnecessary. Since the entire system of equations converge to a global set tolerance
the convergence would occur as the solution progresses to a steady state.
The setup of the directory structure with all the relevant files are now complete. This can be verified again as follows and the directory tree should appear identical to the one shown below
cd myturbulentcurvature/
tree
.
├── 0
│ ├── k
│ ├── nut
│ ├── nuTilda
│ ├── omega
│ ├── p
│ └── U
├── constant
│ ├── polyMesh
│ │ ├── boundary
│ │ ├── faces
│ │ ├── neighbour
│ │ ├── owner
│ │ └── points
│ ├── RASProperties
│ ├── transportProperties
│ └── turbulenceProperties
└── system
├── controlDict
├── decomposeParDict
├── fvSchemes
└── fvSolution
Execution of the solver¶
Renumbering and checking the grid/mesh quality is important before the solver is executed. Renumbering reduces the matrix bandwidth whereas the quality check shows the mesh statistics. These two can be performed by executing the following commands from the top working directory
cd myturbulentcurvature/
caelus run  renumberMesh overwrite
caelus run  checkMesh
With the execution of renumberMesh overwrite
, the user should note the reduction in bandwidth after renumbering occurs. Similarly, when the checkMesh
is performed, the mesh statistics are shown as below
/**\
Caelus 8.04
Web: www.caeluscml.com
\**/
Checking geometry...
Overall domain bounding box (1.5 1 0.8097167) (1.6 0 0.127)
Mesh (nonempty, nonwedge) directions (1 0 1)
Mesh (nonempty) directions (1 0 1)
All edges aligned with or perpendicular to nonempty directions.
Boundary openness (2.89884272002e17 2.78435627054e16 2.89924087017e17) OK.
***High aspect ratio cells found, Max aspect ratio: 6762.08993149, number of cells 1283
<<Writing 1283 cells with high aspect ratio to set highAspectRatioCells
Minimum face area = 1.52715761849e08. Maximum face area = 0.0318702444976. Face area magnitudes OK.
Min volume = 1.52715761849e08. Max volume = 0.000180011068556. Total volume = 0.417212754326. Cell volumes OK.
Mesh nonorthogonality Max: 27.7205032543 average: 17.6783210246
Nonorthogonality check OK.
Face pyramids OK.
Mesh skewness Max: 0.0357278029285 average: 0.000574009677514 OK.
Coupled point location match (average 0) OK.
Failed 1 mesh checks.
End
Apparent from the above output, the checkMesh
indicates that the mesh check has failed reporting in the final message as Failed 1 mesh checks
. This is because of the high aspect ratio meshes present immediate to the wall with very low (\(<< y^+\)) values. Nevertheless, this is just a warning and Caelus will solve on this mesh.
We can utilise the multicore capability of CPUs for performing a parallel computation for large grids, such as the one considered for this tutorial. Before this can be done, the grid has to be decomposed into smaller domains that can be solved by each single CPU core. The number of decomposition should be equal to the number of CPU core available for parallel computing. The decomposition should be carried out through a file decomposeParDict
present in the system
subdirectory which is shown below.
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object decomposeParDict;
}
//
numberOfSubdomains 4;// It is suggested that the numberOfSubdomains be increased based on available resources for validation cases and to reduce the computation time.
method scotch;
distributed no;
roots
(
);
In the above file, the the keyword numberOfSubdomains
defines the number of decomposed subdomains needed and 4
is used which partitions the grid into 4 subdomains. We use scotch
as the method of decomposition which automatically divides the grid. The execution to decompose the grid is carried out again from the top directory as follows
caelus run  decomposePar
Now the decomposition should start and the details of which are displayed in the terminal window. Subsequently, 4 processor directories will be generated as we have requested for 4 divisions of grid as shown below
0 constant processor0 processor1 processor2 processor3 system
The solver can now be executed for parallel computation from the top directory using the following command
caelus run p l myturbulentcurvature.log  simpleSolver
Note that here it is assumed that the parallel computing is available in the host machine. With the execution of the above commands, the simulation begins and the progress of the solution is written to the specified log file (myturbulentcurvature.log
).
The log file can be further processed to look at the convergence history and this can be done as follows
caelus logs w myturbulentcurvature.log
The above statement allows you to look at the convergence of variables with respect to the number of iterations carried out as shown in Figure 45 for pressure.
Results¶
The flow visualisation of velocity and pressure within the convex duct is presented here. In Figure 46 velocity magnitude and pressure are shown for SA model. Due to the convex bend, the thinning of the turbulent boundary layer occurs on the lower surface and the pressure decreases. Whereas the trends are opposite on the upper surface.
Backward Facing Step¶
This tutorial focuses on the turbulent flow over a twodimensional backward facing step using Caelus 9.04. In particular this tutorial emphasises the use of wall functions for a high \(y^+\) grid. The basic steps to setup the directory structure, fluid properties, boundary conditions, turbulence properties etc will be discussed. Flow visualisation such as pressure and velocity contours within the separated region will be shown. With these steps in place, the user should be able to reproduce the tutorial accurately.
Objectives¶
Through this tutorial, the user will get familiar in setting up Caelus simulation for steady, turbulent, incompressible flow over a twodimensional backward facing step with wall functions. Further, the velocity and pressure contours within the separated region will be highlighted. Following would be the steps that would be carried out.
 Background
A brief description about the problem
Geometry and freestream details
 Grid generation
Computational domain and boundary details
Computational grid generation
Exporting grid to Caelus
 Problem definition
Directory structure
Setting up boundary conditions, physical properties and control/solver attributes
 Execution of the solver
Monitoring the convergence
Writing the log files
 Results
Visualisation of flow over within the separated region
Prerequisites¶
By now the user should be familiar with the Linux command line environment via a terminalor caelusconsolu (for Windows OS). The grid for this case can be obtained from the full working case directory shown below and was generated using Pointwise and appropriately converting it to Caelus format.
Background¶
Turbulent flow over a backward facing step is a classical configuration to examine steady separated flows. Here due to the presence of the step, the flow undergoes separation and the subsequent shear layer reattaches downstream forming a recirculation region. The thickness of the boundary layer at the lip of the step and the flow Reynolds number are important parameters that determines the length of separated region. A decrease in pressure or favourable pressure gradient in the immediate vicinity of the step is a classical behaviour that contributes to the increase in drag. The user is suggested to refer to the verification and validation of this case at Twodimensional Backward Facing Step for more detailed analysis.
The schematic of the backward facing step is shown in Figure 47 in twodimensions. A step height (H) of 1.0 m is used with a flow Reynolds number of 36000, which is based on the step height. Air is considered as a fluid with a freestream temperature of 298 K and the freestream velocity corresponds to U = 44.315 m/s. The user should note that the twodimensional geometric plane considered is in \(xz\) directions.
Freestream conditions are detailed in the below table
Fluid 
\(Re_H\) 
\(U~(m/s)\) 
\(p~(m^2/s^2)\) 
\(T~(K)\) 
\(\nu_\infty~(m^2/s)\) 

Air 
\(36000\) 
44.31525 
\((0)\) Gauge 
298.330 
\(1.230979\times10^{3}\) 
Grid Generation¶
A fully structured gird is developed for this geometry and is converted to polyMesh
format. The grid in polyMesh
format can be found in the full working directory of this case at
/tutorials/incompressible/simpleSolver/ras/ACCM_backwardFacingStep2D/constant/polyMesh
The computational domain for the step is shown in Figure 48 which also highlights the boundary conditions used. The upstream flat plate extends for up to 110 step heights so that a fully turbulent boundary layer is ensured prior to the step at which the flow separates. Downstream of the step, the plate extends to 50 step heights giving sufficient length for the flow to reattach. The inlet and the outlet is placed as indicated in the Figure 48 and closely flows the experimental setup. The velocity on the step surfaces is zero, wherein \(u,v,w = 0\) represented through a noslip boundary.
The grid in polyMesh
is in threedimensions, although the flow over the step is twodimensional in nature. A onecell thick grid normal to the flow plane (\(xz\)) is therefore sufficient for Caelus to consider symmetry flow in the \(y\) direction. The two resulting planes that are in (\(xz\)) direction need boundary conditions to be specified. Since the flow is assumed to be 2D, we do not need to solve the flow in the thirddimension and this can be easily achieved in Caelus by specifying empty boundary conditions for each of the two planes. This consequently treats the flow in :math`y` direction as symmetry.
Note
A velocity value of \(v=0\) needs to be specified at appropriate boundaries although no flow is solved in the \(y\) direction.
The 2D step grid in \(xz\) plane is shown in figure %s:numref:tstepgridtutorials which has a total of 189 cells in the streamwise direction and 64 cells in the normal direction, with 20 cells representing the step height. As noted earlier, the grid is developed to use with wall functions and hence the \(y^+ \approx 30\) in this case.
Problem definition¶
This section details various steps needed to setup the turbulent flow over a step. A full working case of this can be found in:
/tutorials/incompressible/simpleSolver/ras/ACCM_backwardFacingStep2D
However,the user is free to start the case setup from scratch consistent with the directory stucture discussed below.
Directory Structure
Note
All commands shown here are entered in a terminal window, unless otherwise mentioned
In order to setup the problem Caelus requires time
, constant
and system
subdirectories within the main working directory. Here, the simulation will be started at time t = 0 s, which requires a time subdirectory named 0
.
The 0
subdirectory contains few additional files in which boundary conditions are specified. In the below table, the list of necessary files are provided based on the turbulence model chosen.
Parameter 
File name 

Pressure (\(p\)) 

Velocity (\(U\)) 

Turbulent viscosity (\(\nu\)) 

Turbulence field variable (\(\tilde{\nu}\)) 

Turbulent kinetic energy (\(k\)) 

Turbulent dissipation rate (\(\omega\)) 

Turbulent dissipation (\(\epsilon\)) 

In this tutorial we consider three turbulence models. They are SpalartAllmaras (SA), \(k\omega\)  Shear Stress Transport (\(\rm{SST}\)) and Realizable \(kepsilon\) models. The content of the files listed above sets the dimensions, initialisation and boundary conditions to the defining problem, which also forms three principle entries required.
The user should note that Caelus is case sensitive and therefore the directory and file setup should be identical to what is shown here.
Boundary Conditions and Solver Attributes
Boundary Conditions
Referring back to Figure 48, the following are the boundary conditions that will be specified:
 Inlet
Velocity: Fixed uniform velocity \(u = 44.31525~m/s\) in \(x\) direction
Pressure: Zero gradient
Turbulence:
Spalart–Allmaras (Fixed uniform values of \(\nu_{t~\infty}\) and \(\tilde{\nu}_{\infty}\) as given in Turbulence freestream conditions for SA model)
\(k\omega~\rm{SST}\) (Fixed uniform values of \(k_{\infty}\), \(\omega_{\infty}\) and \(\nu_{t~\infty}\) as given in Turbulence freestream conditions for SA model)
Realizable \(k\epsilon\) (Fixed uniform value of \(k_{\infty}\), \(\epsilon_{\infty}\) and \(\nu_{t_\infty}\) as given in Turbulent freestream conditions for Realizable k\epsilon Model)
 Symmetry
Velocity: Symmetry
Pressure: Symmetry
Turbulence: Symmetry
 Noslip wall
Velocity: Fixed uniform velocity \(u, v, w = 0\)
Pressure: Zero gradient
Turbulence:
Spalart–Allmaras:
\(\nu_t\): type nutUWallFunction with a uniform value of \(\nu_t=0\)
\(\tilde{\nu}\): type fixedValue with a value of \(\tilde{\nu}=0\)
\(k\omega~\rm{SST}\):
\(k\): type kqRWallFunction with a uniform value of \(k_{\infty}\)
\(\omega\): type omegaWallFunction with a uniform value of \(\omega_{\infty}\)
\(\nu_t\): type nutUWallFunction with a uniform value of \(\nu_t=0\)
Realizable \(k\epsilon\):
\(k\): type kqRWallFunction with a uniform value of \(k_{\infty}\)
\(\epsilon\): type epsilonWallFunction with a uniform value of \(\epsilon=0\)
\(\nu_t\): type nutUWallFunction with a uniform value of \(\nu_t=0\)
 Outlet
Velocity: Zero gradient velocity
Pressure: Fixed uniform gauge pressure \(p = 0\)
Turbulence:
Spalart–Allmaras (Calculated \(\nu_{t}=0\) and Zero gradient \(\tilde{\nu}\))
\(k\omega~\rm{SST}\) (Zero gradient \(k\) and \(\omega\); Calculated \(\nu_t=0\); )
Realizable \(k\epsilon\) (Zero gradient \(k\) and \(\epsilon\); Calculated \(\nu_t=0\); )
 Initialisation
Velocity: Fixed uniform velocity \(u = 44.31525\) in \(x\) direction
Pressure: Zero Gauge pressure
Turbulence:
Spalart–Allmaras (Fixed uniform values of \(\nu_{t~\infty}\) and \(\tilde{\nu}_{\infty}\) as given in Turbulence freestream conditions for SA model)
\(k\omega~\rm{SST}\) (Fixed uniform values of \(k_{\infty}\), \(\omega_{\infty}\) and \(\nu_{t~\infty}\) as given in Turbulence freestream conditions for SA model)
Realizable \(k\epsilon\) (Fixed uniform values of \(k_{\infty}\), \(\epsilon_{\infty}\) and \(\nu_{t~\infty}\) as given in Turbulent freestream conditions for Realizable k\epsilon Model)
First, the file p
, shown below, contains initial and boundary conditions for pressure.
/**\
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 2 0 0 0 0];
internalField uniform 0;
boundaryField
{
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0;
}
symmleft
{
type empty;
}
symmright
{
type empty;
}
topwall
{
type zeroGradient;
}
upstream
{
type symmetryPlane;
}
wall
{
type zeroGradient;
}
}
// ************************************************************************* //
As can be noted from the above file, it begins with a dictionary named FoamFile
which contains the standard set of keywords for version, format, location, class and object names. The explanation of the principle entries are as follows
dimension
is used to specify the physical dimensions of the pressure field. Here, pressure is defined in terms of kinematic pressure with the units (\(m^2/s^2\)) written as
[0 2 2 0 0 0 0]
internalField
is used to specify the initial conditions. It can be either uniform or nonuniform. Since we have a 0 initial uniform gauge pressure, the entry is
uniform 0;
boundaryField
is used to specify the boundary conditions. In this case its the boundary conditions for pressure at all the boundary patches.
In the similar way, the file U
is defined for velocity
/**\
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
location "0";
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 1 0 0 0 0];
internalField uniform (44.31525 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (44.31525 0 0);
}
outlet
{
type zeroGradient;
}
symmleft
{
type empty;
}
symmright
{
type empty;
}
topwall
{
type noSlipWall;
}
upstream
{
type symmetryPlane;
}
wall
{
type noSlipWall;
}
}
// ************************************************************************* //
The principle entries for velocity field are selfexplanatory and the dimensions are typically for that of velocity with the units \(m/s\) ([0 1 1 0 0 0 0]
). Here, both initialisation and inlet have a uniform flow velocity specified with three velocity components. Therefore these should be set to uniform (44.31525 0 0);
while the wall is given a noslip condition.
The turbulent properties are also required to be specified at the boundary patches and these can be done similar to p
and U
. We have the file nut
, which defines turbulent kinematic viscosity in the domain as follows
/**\
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object nut;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 1 0 0 0 0];
internalField uniform 0.09811;
boundaryField
{
inlet
{
type calculated;
value uniform 0.09811;
}
outlet
{
type calculated;
value uniform 0;
}
symmleft
{
type empty;
}
symmright
{
type empty;
}
topwall
{
type nutUWallFunction;
Cmu 0.09;
kappa 0.41;
E 9.8;
value uniform 0;
}
upstream
{
type symmetryPlane;
}
wall
{
type nutUWallFunction;
Cmu 0.09;
kappa 0.41;
E 9.8;
value uniform 0;
}
}
// ************************************************************************* //
As noted above, the turbulent viscosity is specified as kinematic and therefore the units are in \(m^2/s\) ([0 2 1 0 0 0 0]
). The turbulent viscosity value at freestream, specified at inlet patch is calculated as detailed in Turbulence freestream conditions for SA model, Turbulent freestream conditions for k\omega~\rm{SST} Model and Turbulent freestream conditions for Realizable k\epsilon Model for SA, SST and RKE models respectively and is specified accordingly. The same value also goes for internalField
.
The next turbulent property set is the nuTilda
which is a turbulent field variable, specified to only SA model and has the same units ([0 2 1 0 0 0 0]
) as kinematic turbulent viscosity. Details pertaining to this are given in Turbulence freestream conditions for SA model. The following is the file nuTilda
with the entries specified for the boundaryField
are identical to that of turbulent kinematic viscosity explained above.
/**\
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object nuTilda;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 1 0 0 0 0];
internalField uniform 0.003692937;
boundaryField
{
inlet
{
type fixedValue;
value uniform 0.003692937;
}
outlet
{
type zeroGradient;
}
symmleft
{
type empty;
}
symmright
{
type empty;
}
topwall
{
type fixedValue;
value uniform 0;
}
upstream
{
type symmetryPlane;
}
wall
{
type fixedValue;
value uniform 0;
}
}
// ************************************************************************* //
We can now proceed to the file k
which is specific to both \(komega~\rm{SST}\) and Realizable \(kepsilon\) models and represents turbulent kinetic energy.
/**\
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object k;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 2 0 0 0 0];
internalField uniform 0.2945755;
boundaryField
{
inlet
{
type fixedValue;
value uniform 0.2945755;
}
outlet
{
type zeroGradient;
}
symmleft
{
type empty;
}
symmright
{
type empty;
}
topwall
{
type kqRWallFunction;
value uniform 0.2945755;
}
upstream
{
type symmetryPlane;
}
wall
{
type kqRWallFunction;
value uniform 0.2945755;
}
}
// ************************************************************************* //
The unit of kinetic energy is \(m^2/s^2\) and this is set in dimensions
as [0 2 2 0 0 0 0]
. The value of \(k\) (refer Turbulent freestream conditions for k\omega~\rm{SST} Model and Turbulent freestream conditions for Realizable k\epsilon Model) needs to be specified for internalField
, inlet
and wall
. Please note that for wall kqRWallFunction
with values of freestream are required.
We now proceed to the file omega
, specific to only \(k\omega~\rm{SST}\) model and the value for this is evaluated as detailed in Turbulent freestream conditions for k\omega~\rm{SST} Model
/**\
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object omega;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 0 1 0 0 0 0];
internalField uniform 97.37245;
boundaryField
{
inlet
{
type fixedValue;
value uniform 97.37245;
}
outlet
{
type zeroGradient;
}
symmleft
{
type empty;
}
symmright
{
type empty;
}
topwall
{
type omegaWallFunction;
Cmu 0.09;
kappa 0.41;
E 9.8;
beta1 0.075;
value uniform 97.37245;
}
upstream
{
type symmetryPlane;
}
wall
{
type omegaWallFunction;
Cmu 0.09;
kappa 0.41;
E 9.8;
beta1 0.075;
value uniform 97.37245;
}
}
// ************************************************************************* //
The unit of specific rate of dissipation for \(\omega\) is \(1/s\) which is set in dimensions
as [0 0 1 0 0 0 0]
. The internalField
and inlet
gets a fixedValue
. Note that for wall boundaryField
, we specify omegaWallFunction
required for high \(y^+\).
The final file in this class is the epsilon
, specific to only Realizable \(k\epsilon\) model and the value for this is evaluated as detailed in Turbulent freestream conditions for Realizable k\epsilon Model
/**\
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object epsilon;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 3 0 0 0 0];
internalField uniform 0.079598;
boundaryField
{
inlet
{
type fixedValue;
value uniform 0.079598;
}
outlet
{
type zeroGradient;
}
symmleft
{
type empty;
}
symmright
{
type empty;
}
topwall
{
type epsilonWallFunction;
value uniform 0;
}
upstream
{
type symmetryPlane;
}
wall
{
type epsilonWallFunction;
value uniform 0;
}
}
// ************************************************************************* //
The unit of turbulent dissipation for \(\epsilon\) is \(m^2/s^3\) which is set in dimensions
as [0 2 3 0 0 0 0]
. The internalField
and inlet
gets a fixedValue
similar to \(\omega\). Note that for wall boundaryField
, we specify epsilonWallFunction
required for high \(y^+\).
Before proceeding to the next step, it is vital to ensure that the boundary conditions (inlet
, outlet
, wall
, etc) added in the above files should be the grid boundary patches (surfaces) generated by grid generation tool and their names are identical. Additionally, the two boundaries \(xz\) plane named here as symmleft
and symmright
have empty
boundary conditions which forces Caelus to assume the flow to be in twodimensions. With this, the setting up of the boundary conditions are complete.
Grid file and Physical Properties
The files associated with the grid need to be placed in polyMesh
subdirectory, which resides in the constant
directory. We use identical grid for all the three turbulence simulation. In addition, the physical properties are specified in various different files present in the constant
directory.
There are three files listed in addition to the polyMesh
subdirectory. The first one, RASProperties
in which the ReynoldsAverageStress (RAS) model is specified, which is shown below. Please note that depending on the turbulence model you wish to run with, the line that corresponds to that specific model should be enabled.
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object RASProperties;
}
//
// For SpalaratAlamaras Model, enable the below line
//RASModel SpalartAllmaras;
// For komega SST Model, enable the below line
//RASModel kOmegaSST;
// For Realizable kepsilon Model, enable the below line
RASModel realizableKE;
turbulence on;
printCoeffs on;
Second from the list is the transportProperties
file, where laminar kinematic viscosity is specified as shown below
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
//
transportModel Newtonian;
nu nu [0 2 1 0 0 0 0] 1.230979E3;
//Cross Power law and BirdCarreau nonlinear viscosity models
// 
CrossPowerLawCoeffs
{
nu0 nu0 [0 2 1 0 0 0 0] 1e06;
nuInf nuInf [0 2 1 0 0 0 0] 1e06;
m m [0 0 1 0 0 0 0] 1;
n n [0 0 0 0 0 0 0] 1;
}
BirdCarreauCoeffs
{
nu0 nu0 [0 2 1 0 0 0 0] 1e06;
nuInf nuInf [0 2 1 0 0 0 0] 1e06;
k k [0 0 1 0 0 0 0] 0;
n n [0 0 0 0 0 0 0] 1;
}
//
The viscous behaviour is modelled as Newtonian and hence the keyword Newtonian
is used for the transportModel
and the molecular (laminar) kinematic viscosity (\(nu\)) is given having the units \(m^2/s\) ([0 2 1 0 0 0 0]
).
The final file in this class is the turbulenceProperties
file, which sets the simulationType
to RASModel
. Here, SA, \(k\omega~\rm{SST}\) and Realizable \(k\epsilon\) are classified as Reynolds Average Stress (RAS) models.
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object turbulenceProperties;
}
//
simulationType RASModel;
//
Controls and Solver Attributes
This section details the files required to control the simulation, discretization methods and linear solvers. These files are to be placed in the system
directory. First, navigate to the system
directory
The following is the controlDict
file,
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
//
application simpleSolver;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 10000;
deltaT 1.0;
writeControl runTime;
writeInterval 500;
purgeWrite 0;
writeFormat ascii;
writePrecision 12;
writeCompression true;
timeFormat general;
timePrecision 12;
runTimeModifiable true;
//
// Function Objects to obtain forces
/*functions
{
forces
{
type forces;
functionObjectLibs ("libforces.so");
patches ( wall );
CofR (0 0 0);
rhoName rhoInf;
rhoInf 15.13009E3;
outputControl timeStep;
outputInterval 100;
}
}
Here, simpleSolver
is used and the simulation begins at t = 0 s
. This now explains the logical need for having a 0
directory where the data files are read at the beginning of the run, which is t = 0 s
in this case. Therefore, the keyword startFrom
is set to startTime
, where startTime
is set to 0
. The simulation would be carried out as steadystate and therefore we require to specify the total number of iterations as a keyword for endTime
. Via the writeControl
and writeInterval
keywords, the solution intervals at which they are saved can be specified. Also note that a function object to obtain the force over the wall (step surface) for every 100 iterations is included. In order to obtain this, a freestream density (rhoInf
) need to be specified.
The discretization schemes for the finite volume discretization that will be used should be set through the fvSchemes
file show below
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object fvSchemes;
}
//
ddtSchemes
{
default steadyState;
}
gradSchemes
{
default Gauss linear;
grad(p) Gauss linear;
grad(U) Gauss linear;
}
divSchemes
{
default none;
div(phi,U) Gauss linearUpwind grad(U);
div(phi,nuTilda) Gauss upwind; // Will be used for SA model only
div(phi,k) Gauss upwind; // will be used for kepsilon & komega only
div(phi,omega) Gauss upwind; // Will be used for komega model only
div(phi,epsilon) Gauss upwind;
div((nuEff*dev(T(grad(U))))) Gauss linear;
div(phi,symm(grad(U))) Gauss linear; // Will be used for SACurvaturecorr only
}
laplacianSchemes
{
default none;
laplacian(nu,U) Gauss linear corrected;
laplacian(nuEff,U) Gauss linear corrected;
laplacian(DnuTildaEff,nuTilda) Gauss linear corrected; // Will be used for SA model only
laplacian(DkEff,k) Gauss linear corrected; // Will be used for komega & komega only
laplacian(DomegaEff,omega) Gauss linear corrected; // Will be used for komega model only
laplacian(DepsilonEff,epsilon) Gauss linear corrected;
laplacian(rAUf,p) Gauss linear corrected;
laplacian(1,p) Gauss linear corrected;
}
interpolationSchemes
{
default linear;
interpolate(HbyA) linear;
}
snGradschemes
{
default corrected;
}
//
The linear solver controls and tolerances are set in fvSolution
as given below
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSolution;
}
//
solvers
{
p
{
solver PCG;
preconditioner SSGS;
tolerance 1e8;
relTol 0.01;
}
U
{
solver PBiCGStab;
preconditioner USGS;
tolerance 1e7;
relTol 0.01;
}
"(komeganuTildaepsilon)"
{
solver PBiCGStab;
preconditioner USGS;
tolerance 1e08;
relTol 0;
}
}
SIMPLE
{
nNonOrthogonalCorrectors 1;
pRefCell 0;
pRefValue 0;
}
relaxationFactors
{
p 0.3;
U 0.3;
nuTilda 0.5;
k 0.5;
omega 0.5;
epsilon 0.5;
}
The user should now take a note that in the fvSolution
file, different linear solvers are used to solve for velocity, pressure and turbulence quantities. We also set the nNonOrthogonalCorrectors
to 1 for this case. To ensure the stability of the solution, the relaxation is set to primary and turbulent variables. The relTol
is not set to 0
unlike a timeaccurate setup as we are solving for a steadystate solution and a very low (\(\approx 0\)) tolerance at every iteration is unnecessary. Since the entire system of equations converge to a global set tolerance
the convergence would occur as the solution progresses to a steady state.
With this, the setup of the directory structure with all the relevant files are complete. This can be verified again by issuing the following command and the directory tree should appear identical to the one shown below
cd myturbulentstep/
tree
.
├── 0
│ ├── epsilon
│ ├── k
│ ├── nut
│ ├── nuTilda
│ ├── omega
│ ├── p
│ └── U
├── constant
│ ├── polyMesh
│ │ ├── boundary
│ │ ├── faces
│ │ ├── neighbour
│ │ ├── owner
│ │ └── points
│ ├── RASProperties
│ ├── transportProperties
│ └── turbulenceProperties
└── system
├── controlDict
├── fvSchemes
└── fvSolution
Execution of the solver¶
Prior to the execution of the solver, it is important to renumber and to carry out a quality check on the grid/mesh. Renumbering reduces the bandwidth whereas the quality check shows the mesh statistics. These two can be performed by executing the following commands from the top working directory
caelus run  renumberMesh overwrite
caelus run  checkMesh
When the renumberMesh
is performed, the user should take note of the bandwidth before and after the mesh renumbering. In a similar manner, when the checkMesh
is performed, the mesh statistics are shown as below
/**\
Caelus 8.04
Web: www.caeluscml.com
\**/
Checking geometry...
Overall domain bounding box (130 0 0) (50 1 9)
Mesh (nonempty, nonwedge) directions (1 0 1)
Mesh (nonempty) directions (1 0 1)
All edges aligned with or perpendicular to nonempty directions.
Boundary openness (7.60426729195e19 3.46752522136e15 6.3479100872e18) OK.
Max cell openness = 2.1851506874e16 OK.
Max aspect ratio = 205.17630165 OK.
Minimum face area = 0.00417125. Maximum face area = 10.4153124857. Face area magnitudes OK.
Min volume = 0.00417125. Max volume = 1.987102972. Total volume = 1490. Cell volumes OK.
Mesh nonorthogonality Max: 12.1317921302 average: 0.100202600449
Nonorthogonality check OK.
Face pyramids OK.
Mesh skewness Max: 0.00821683161155 average: 8.01613031841e06 OK.
Coupled point location match (average 0) OK.
Mesh OK.
End
At this stage, the solver can be executed from the top directory using the following command
caelus run l myturbulentstep.log  simpleSolver
The log file can be further processed to look at the convergence history and this can be done as follows
caelus logs w myturbulentstep.log
The caelusplotResiduals l
allows you to look at the convergence of various variables with respect to the number of iterations carried out and is shown in Figure 50.
Results¶
The flow within the separated region is visualised here through the velocity and pressure contours. In Figure 51 velocity magnitude and pressure contours are shown for SA model. Immediate downstream of the step a decrease in pressure is seen followed by an increases as the shear layer reattaches to the step. This is consistent with the formation of a low velocity recirculating region behind the step.