Caelus DocumentationÂ¶
 Version
9.04
 Date
Feb 07, 2020
Caelus is a software library suitable for numerical simulations of problems in Continuum Mechanics (CM), with many applications in Computational Fluid Dynamics (CFD) for a wide range of scientific and engineering applications across commercial and academic environments. Caelus was forked from OpenFOAMÂ®, version 2.1.1 released by OpenFOAMÂ® Foundation which is mainly a collection of libraries written in C++. Caelus is developed and maintained by Applied CCM Pty Ltd.
Being a derivative of OpenFOAMÂ®, it is released under the GPL. You can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. See the GNU General Public License for more details.
This documentation is split into three parts: a tutorial, a theory and a validation manual. New users should start with the tutorial manual that provides examples of usage.
Tutorial ManualÂ¶
Tutorials: IntroductionÂ¶
Caelus is a software library suitable for numerical simulations of problems in Continuum Mechanics (CM), with many applications in Computational Fluid Dynamics (CFD) for a wide range of scientific and engineering applications across commercial and academic environments. Caelus is forked from OpenFOAM, version 2.1.1 released by OpenFOAM Foundation which is mainly a collection of libraries written in C++. Caelus is developed and maintained by Applied CCM Pty Ltd.
Welcome to Caelus tutorials page. Here, various capabilities of Caelus are presented through a range of tutorials. These closely follow the case studies performed for the purpose of code validation in a view that users can repeat those numerical experiments successfully while understanding the procedures involved. Each tutorial contains information about the problem description, case setup, solver settings and postprocessing details. It is suggested that user follow these steps methodically.
Note: In order to use Linux based commands on Windows, as given in the tutorial guide, it is suggest that the user install Caelus with the packaged installer. This guide assumes the user has also installed the Caelus Python Library (CPL).
Tutorials: Incompressible Laminar FlowÂ¶
Laminar Flat PlateÂ¶
In this tutorial, simulation of laminar incompressible flow over a twodimensional sharp leadingedge flat plate using Caelus 9.04 is introduced here. First, prerequisites to begin a Caelus simulation is discussed followed by various dictionary entries defining fluid properties, boundary conditions, solver setting, etc that are needed. Finally, the presence of laminar boundary layer is visualised using velocity contours. Here, the basic procedures of running Caelus is shown in sufficient detail such that the user feels comfortable with the usage.
ObjectivesÂ¶
With the completion of this tutorial, the user will be familiar with setting up Caelus simulation for steady, laminar, incompressible flow over flatplates in twodimensions and subsequently postprocessing the results. Following are some of the steps 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 the laminar 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 that Caelus is installed correctly with appropriate environment variables set. The grid used here is generated using Pointwise and the user is free to use their choice of grid generation tool having exporting capabilities to the Caelus grid format.
BackgroundÂ¶
The flow over a flatplate presents an ideal case where initial steps of a Caelus simulation can be introduced to the user in easy steps. Here, laminar, incompressible flow over a sharpleading edge plate is solved in a timeaccurate approach. This results in the formation of laminar boundary layer which is then compared with the famous Blasius [6] analytical solution in the form of a nondimensional shear stress distribution (skinfriction coefficient). For more details, the user is suggested to refer the validation of flatplate in section Flat plate validation.
The length of the flatplate considered here is \(L = 0.3048~m\) with a Reynolds number based on the total length of \(Re = 200,000\). A schematic of the geometry is shown in Figure 1, wherein \(U\) is the flow velocity in \(x\) direction. An inflow temperature of \(T = 300~K\) can assumed for the fluid air which corresponds to a kinematic viscosity (\(\nu\)) of \(\nu = 1.5896306 \times 10^{5}~m^2/s\). Using the values of \(Re\), \(L\) and \(\nu\), we can evaluate the freestream velocity to \(U = 10.4306~m/s\).
In the following table, details of the freestream conditions are provided.
Fluid 
\(L~(m)\) 
\(Re\) 
\(U~(m/s)\) 
\(p~(m^2/s^2)\) 
\(T~(K)\) 
\(\nu~(m^2/s)\) 

Air 
0.3048 
200,000 
69.436113 
Gauge (0) 
300 
\(1.58963\times10^{5}\) 
Grid GenerationÂ¶
As noted earlier, Pointwise has been used here to generate a hexahedral grid. Specific details pertaining to its usage are not discussed here, rather a more generic discussion is given about the computational domain and boundary conditions. This would facilitate the user to obtain a Caelus compatible grid using their choice of grid generating tool.
The computational domain is a rectangular block encompassing the flatplate. The below (Figure 2) shows the details of the boundaries that will be used in twodimensions (\(xy\) plane). First, the flatplate, which is our region of interest (highlighted in blue) is extended between \(0\leq x \leq 0.3048~m\). Because of viscous nature of the flow, the velocity at the wall is zero which can be represented through a noslip boundary (\(u, v, w = 0\)). Upstream of the leading edge, a slip boundary will be used to simulate a freestream flow approaching the flatplate. However, downstream of the plate, it would be ideal to further extend up to three plate lengths (highlighted in green) with noslip wall. This would then ensure that the boundary layer in the vicinity of the trailing edge is not influenced by outlet boundary. Since the flow is incompressible (subsonic), the disturbance caused by the pressure can propagate both upstream as well as downstream. Therefore, the placement of the inlet and outlet boundaries are to be chosen to have minimal or no effect on the solution. The inlet boundary as shown will be placed at start of the slipwall at \(x = 0.06~m\) and the outlet at the exit plane of the noslip wall (green region) at \(x = 1.2192~m\). Both inlet and outlet boundary are between \(0\leq y \leq 0.15~m\). A slipwall condition is used for the entire top boundary.
The 2D structured grid is shown in Figure 3. Since Caelus is a 3D solver, it necessitates the grid to be in 3D. Therefore, the 3D grid should be obtained through extruding the 2D gird in the \(z\) direction by a minimum of onecell thick. The final 3D grid should be then exported to Caelus format (polyMesh). The two \(xy\) planes obtained as a result of grid extrusion need boundary conditions to be specified. As the flow over a flatplate is generally 2D, we do not need to solve the flow in the third dimension. This can easily be achieved in Caelus by specifying empty
boundary condition for each of those two planes.
Note
A velocity value of \(w=0\) needs to be specified at appropriate boundaries although no flow is solved in the \(z\) direction.
The flatplate has a total of 400 cells over the region of interest between \(0 \leq x \leq 0.3048~m\) and 286 cells in the noslip wall that extends for an additional 3 plate lengths (green region in the above figure). In the wall normal direction, 298 cells are used and sufficient refinement close to the wall was made to ensure that accurate boundary layer is captured.
Problem definitionÂ¶
Several important instructions would be shown here to setup the flatplate problem along with the detail of configuration files used. A full working case can be found in:
/tutorials/incompressible/simpleSolver/laminar/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
In order to setup the problem, three main subdirectories containing all the relevant information are used. Caelus requires time
, constant
and system
subdirectories. Since we begin the simulation at time \(t = 0~s\), the time
subdirectory should be just 0
.
The 0
subdirectory is where additional two files, p
and U
for pressure (\(p\)) and velocity (\(U\)) respectively are kept. The contents of these two files sets the dimensions, initialisation and boundary conditions to the problem, which also form three principle entries required.
It should be noted that Caelus is case sensitive and therefore the user should setup the directories (if applicable), files and the contents identical to what is mentioned here.
Boundary Conditions and Solver Attributes
Boundary Conditions
First let us look at setting up the boundary conditions. Referring back to Figure 2, following are the boundary conditions that need to be specified:
 Inlet
Velocity: Fixed uniform velocity \(u = 10.4306~m/s\) in \(x\) direction
Pressure: Zero gradient
 Slip wall
Velocity: slip
Pressure: slip
 Noslip wall
Velocity: Fixed uniform velocity \(u, v, w = 0\)
Pressure: Zero gradient
 Outlet
Velocity: Zero gradient velocity
Pressure: Fixed uniform gauge pressure \(p = 0\)
 Initialisation
Velocity: Fixed uniform velocity \(u = 10.4306~m/s\) in \(x\) direction
Pressure: Zero Gauge pressure
Now let us look at the contents and significance of each file in these subdirectories beginning with the pressure (\(p\)) file.
/**\
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
{
downstream
{
type zeroGradient;
}
inflow
{
type zeroGradient;
}
outflow
{
type fixedValue;
value uniform 0;
}
symmleft
{
type empty;
}
symmright
{
type empty;
}
top
{
type slip;
}
upstream
{
type slip;
}
wall
{
type zeroGradient;
}
}
// ************************************************************************* //
The above file begins with a dictionary named FoamFile
which contains standard set of keywords for version, format, location, class and object names. The following are the principle entries required for this case.
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
, shown below sets the boundary conditions 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 (10.4306 0 0);
boundaryField
{
downstream
{
type fixedValue;
value uniform (0 0 0);
}
inflow
{
type fixedValue;
value uniform (10.4306 0 0);
}
outflow
{
type zeroGradient;
}
symmleft
{
type empty;
}
symmright
{
type empty;
}
top
{
type slip;
}
upstream
{
type slip;
}
wall
{
type fixedValue;
value uniform (0 0 0);
}
}
// ************************************************************************* //
The principle entries for velocity field are self explanatory and the dimension are typical for velocity with 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 (10.43064759 0 0)
representing three components of velocity. In a similar manner, inflow
, wall
and downstream
boundary patches have three velocity components.
At this stage 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 tool and their names are identical. Further, the two boundaries in \(xy\) plane obtained due to grid extrusion have been named as symmleft
and symmright
with specifying empty
boundary conditions forcing Caelus to assume the flow to be in twodimensions. This completes the setting up of boundary conditions.
Grid file and Physical Properties
The flatplate grid files that is generated in the Caelus format resides in the constant/polyMesh
subdirectory. It contains information relating to the points, faces, cells, neighbours and owners of the mesh.
In addition, the physical properties are specified in various different files present in the directory constant
. In the transportProperties
file, transport model and kinematic viscosity are specified. The contents of this file are as follows
/**
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.5896306e5;
As the flow is Newtonian, the transportModel
is specified with Newtonian
keyword and the value of kinematic viscosity (nu
) is given which has the units \(m^2/s\) ([0 2 1 0 0 0 0]
).
The next file is the turbulenceProperties
file, where the type of simulation is specified.
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object turbulenceProperties;
}
//
simulationType laminar;
Since the flow here is laminar, the simulationType
would be laminar
. Similarly, in the RASProperties
file, RASModel
is set to laminar
as shown below.
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object RASProperties;
}
//
RASModel laminar;
turbulence off;
printCoeffs on;
Controls and Solver Attributes
The files required to control the simulation and specifying the type of discretization method along with the linear solver settings are present in the system
directory.
The controlDict
file is shown 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 2000;
deltaT 1;
writeControl runTime;
writeInterval 500;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression uncompressed;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
Here, the application simpleSolver
refers to the SIMPLE solver that will be used. We also begin the simulation at \(t = 0~s\), which logically explains the need for 0
directory where the data files are read at the beginning of the run. Therefore, we need to set the keyword startFrom
to startTime
, where startTime
would be 0
. Since simpleSolver
is a steadystate solver, the keyword endTime
corresponds to the total number of iterations.The interval at which output files are written is controlled by writeControl
and writeInterval
keywords.
The discretization schemes and parameters are specified 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((nuEff*dev(T(grad(U))))) Gauss linear;
}
laplacianSchemes
{
default none;
laplacian(nu,U) Gauss linear corrected;
laplacian(nuEff,U) Gauss linear corrected;
laplacian(p) Gauss linear corrected;
laplacian(rAUf,p) Gauss linear corrected;
laplacian((1A(U)),p) Gauss linear corrected;
}
interpolationSchemes
{
default linear;
interpolate(HbyA) linear;
}
snGradschemes
{
default corrected;
}
Here, the discretization schemes for finite volume discretization of timederivative, gradient, divergence and Laplacian terms are specified.
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 1e08;
relTol 0.01;
maxIter 500;
}
U
{
solver PBiCGStab;
preconditioner USGS;
tolerance 1e7;
relTol 0.1;
}
}
SIMPLE
{
nNonOrthogonalCorrectors 0;
pRefCell 0;
pRefValue 0;
consistent true;
}
// relexation factors
relaxationFactors
{
p 1.0;
U 0.7;
}
Different linear solvers are used here to solve pressure and velocity. The GAMG
preconditioner is applied to improve convergence of the p
solver. Also, by setting the keyword consistent
to true
, SIMPLEC
solver is used and therefore a relaxation factor of 1.0
is applied for p
.
The setup of the directory structure along with the relevant files are completed. Let us verify again by typing the following command (on Linux and Caelusconsole prompt) and the directory tree should be identical to the one shown below
tree
.
â”śâ”€â”€ 0
â”‚Â Â â”śâ”€â”€ p
â”‚Â Â â””â”€â”€ U
â”śâ”€â”€ constant
â”‚Â Â â”śâ”€â”€ polyMesh
â”‚Â Â â”‚Â Â â”śâ”€â”€ boundary
â”‚Â Â â”‚Â Â â”śâ”€â”€ faces
â”‚Â Â â”‚Â Â â”śâ”€â”€ neighbour
â”‚Â Â â”‚Â Â â”śâ”€â”€ owner
â”‚Â Â â”‚Â Â â””â”€â”€ points
â”‚Â Â â”śâ”€â”€ RASProperties
â”‚Â Â â”śâ”€â”€ transportProperties
â”‚Â Â â””â”€â”€ turbulenceProperties
â””â”€â”€ system
â”śâ”€â”€ controlDict
â”śâ”€â”€ fvSchemes
â””â”€â”€ fvSolution
Execution of the solverÂ¶
Before execution of the solver, renumbering of the grid or mesh needs to be performed as well as checking the quality of the grid. Renumbering the gridcell list is vital to reduce the matrix bandwidth while the 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
It is suggested for the user to take note of the bandwidth before and after the mesh renumbering. When the checkMesh
is performed, the following information is given as an output
/**\
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
The mesh nonorthogonality as reported above is 0 and therefore no nonorthogonal corrections are needed in this case. In the case of mesh nonorthogonality being high, certain number of corrections are to be accounted for which can be set in the fvSolution
file with the keyword nNonOrthogonalCorrectors
. The next step is to execute the solver and monitoring the progress of the solution. The solver is always executed from the top directory which is ACCM_flatPlate2D
in our case as below
caelus run l myflatplate.log  simpleSolver
Now the simulation begins and the output of the solution process is written to the log file, myflatplate.log
. In another the terminal window the solver process a ca nbe watched through plotting of the residuals
caelus logs w myflatplate.log
Now the convergence of pressure can be seen with respect to number of iterations along with other convergence properties.
ResultsÂ¶
A brief qualitative data of the simulated flatplate results are given here. Since the aim here is to obtained the steady solution, the results therefore represent the final steady state condition. In Figure 5, the contours of velocity magnitude are shown which highlights the development of the boundary layer. Since the Reynolds number of the flow is high, thickness of the boundary layer is quite thin in comparison to the length of the plate.
Tee JunctionÂ¶
This tutorial introduces the steady, laminar flow through a twodimensional \(90^\circ\) teejunction. Here, we will be using Caelus 9.04 and some of the basic requirements to begin a Caelus simulation will be shown. These include, specifying input data defining boundary conditions, fluid properties and discretization/solver settings. At the end, visualisation is carried out to look at the pressure and velocity contours within the teejunction. The details in running a Caelus simulation for the teejunction will be shown in sufficient detail so that the user is able to repeat the numerical experiment.
ObjectivesÂ¶
Through the completion of this tutorial, the user will be able to setup the Caelus simulation for laminar, incompressible flow through a twodimensional junction and subsequently estimate the flowsplit. Following are the steps that will be carried out in this tutorial
 Background
A brief description about the problem
Geometry and flow conditions
 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 inside the teejunction
PrerequisitesÂ¶
It is assumed that the user is familiar with the Linux command line environment using a terminal or Caelusconsole (for Windows OS) and that Caelus is installed correctly with appropriate environment variables set. The grid used here is generated using Pointwise and the user is free to use their choice of grid generation tool having exporting capabilities to the Caelus grid format.
BackgroundÂ¶
The flow in a teejunction presents with a simple introduction in carrying out separated flow simulation in Caelus. Because of the presence of a side branch, the flow separates forming a recirculating region. This in turn affects the mass flow through main and side branches altering the flow splits. For more details, the user can refer to the validation example show in Tee Junction.
A schematic of the teejunction geometry is shown in Figure 6. Here, \(L = 3.0~m\) and \(W = 1.0~m\) with a Reynolds number of \(Re_w = 300\) based on the side branch width. The velocity (\(V\)) is assumed to be \(1~m/s\) in the \(y\) direction for simplicity. With these flow parameters, the kinematic viscosity can be evaluated to \(\nu = 0.00333~m^2/s\).
In the table, the flow parameters are provided.
\(Re\) 
\(V~(m/s)\) 
\(p~(m^2/s^2)\) 
\(\nu~(m^2/s)\) 

300 
1.0 
0 Gauge 
0.00333 
Grid GenerationÂ¶
A hexahedral grid is generated for the teejunction grid using Pointwise. Specific grid generation details are not discussed here, however information regarding the computational domain and boundary conditions are provided. With this, the user will be able to generate an equivalent grid using their choice of tool.
The computational domain should follow the teejunction geometry and the details are shown in Figure 7. Due to viscous nature of the flow, the velocity at the walls is zero, which should be represented through a noslip boundary condition (\(u, v, w = 0\)) highlighted in blue. A fully developed laminar flow with a parabolic velocity profile will also be applied as a profile boundary at the inlet. This would ensure that the velocity is fully developed before it approaches the side branch, otherwise requiring to have sufficient length in main branch for the flow to develop. As shown in the geometry, the domain will have two outlets, one at the end of the main channel and the other at the end of side branch. Also of further importance is that the exit pressures at the two outlets are set equal.
The 2D structured grid is shown in Figure 8 for a \(xy\) plane. Caelus is a 3D solver and hence requires a 3D grid although the flow here is assumed to be twodimensional. The 3D grid was obtained by extruding the 2D grid in the third (\(z\)  direction) dimension by onecell thick. The two \(xy\) planes obtained as a result of mesh extrusion needs boundary conditions to be specified. As the flow is assumed to be 2D, we do not need to solve the flow in \(z\) direction and this was achieved by specifying empty boundary condition for each of those two planes.
Note
A velocity value of \(w=0\) needs to be specified at appropriate boundaries although no flow is solved in the \(z\) direction.
A coarse grid with a total of 2025 cells is made for the teejunction of which, 90 cells are distributed along the height of the main channel, and 45 along the length of the side branch. The distribution is such that a dimensional length, \(L = 1~m\) has a total of 45 cells and this gives a distribution of \((2/3)45 = 30\) cells for the \((2/3) L\) segment of the main channel. The width \(W\) consists of 15 cells.
Problem definitionÂ¶
We begin with instructions to setup the teejunction problem and subsequently configuring the required input files. A full working case can be found in the following location
/tutorials/incompressible/simpleSolver/laminar/ACCM_teeJunction
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 the following directories are needed:time
, constant
and system
, where relevant files are placed. In this case, the time
directory would be named 0
as we begin the simulation at time \(t = 0~s\).
In the 0
subdirectory, two additional files p
and U
for pressure (\(p\)) and velocity (\(U\)) are present. The input contents of these two files set the dimensions, initial and boundary conditions to the problem. These three forms the principle entries required.
It should be noted that Caelus is case sensitive and therefore the user should setup the directories (if applicable), files and the contents identical to what is mentioned here.
Boundary Conditions and Solver Attributes
Boundary Conditions
Referring to Figure 7, the following boundary conditions will be applied:
 Inlet
Velocity: Parabolic velocity profile; average velocity of \(v = 1.0~m/s\) in \(y\) direction
Pressure: Zero gradient
 Noslip wall
Velocity: Fixed uniform velocity \(u, v, w = 0\)
Pressure: Zero gradient
 Outlet1
Velocity: Zero gradient velocity
Pressure: Fixed uniform gauge pressure \(p = 0\)
 Outlet2
Velocity: Zero gradient velocity
Pressure: Fixed uniform gauge pressure \(p = 0\)
 Initialisation
Velocity: Fixed uniform velocity \(u, v, w = 0\)
Pressure: Zero Gauge pressure
The first quantity to define would be the pressure (\(p\)) and this is done by in file p
, which can be opened using a text editor.
/**
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
{
wall
{
type zeroGradient;
}
inlet
{
type zeroGradient;
}
outlet1
{
type fixedValue;
value uniform 0;
}
outlet2
{
type fixedValue;
value uniform 0;
}
symmleftright
{
type empty;
}
}
As we can see, the above file begins with a dictionary named FoamFile
which contains standard set of keywords such as version, format, location, class and object names. This is followed by the principle elements
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 manner, the input data for the velocity file is shown below
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
object U;
}
//
dimensions [0 1 1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
wall
{
type fixedValue;
value uniform (0 0 0);
}
inlet
{
type groovyBC;
variables "Vmax=1.0;xp=pts().x;minX=min(xp);maxX=max(xp);para=(maxXpos().x)*(pos().xminX)/(0.25*pow(maxXminX,2))*normal();";
valueExpression "Vmax*para";
value uniform (0 1 0);
}
outlet1
{
type zeroGradient;
}
outlet2
{
type zeroGradient;
}
symmleftright
{
type empty;
}
}
As noted above, the principle entries for the velocity filed are self explanatory with the typical dimensional units of \(m/s\) ([0 1 1 0 0 0 0]
). The initialisation of the flow is done at \(0~m/s\) which is set using internalField
to uniform (0 0 0);
which represents three components of velocity.
As discussed previously, a parabolic velocity profile is set for the inlet
. This is done through an external library to Caelus called as groovyBC which allows to specify boundary conditions in terms of an expression. In this case an expression for a parabolic velocity profile in the \(y\) direction is obtained by setting the following expression
"Vmax=1.0;xp=pts().x;minX=min(xp);maxX=max(xp);para=(maxXpos().x)*(pos().xminX)/(0.25*pow(maxXminX,2))*normal();"
and the velocity at the centerline is uniform at \(1~m/s\) represented through uniform (0 1 0);
The boundary conditions (inlet
, outlet
, wall
, etc) specified above should be the grid boundary patches (surfaces) generated by the gridgeneration tool. It should be ensured by the user that their names are identically matched. In addition, the two boundaries in \(xy\) plane obtained due to grid extrusion are named as symmleftright
with applying empty
boundary conditions enforcing the flow to be in twodimensions. It should however be noted that the two planes are grouped together and the empty
patch is applied. This is a capability of Caelus, where similar boundaries can be grouped together and is also used for the wall
boundary, where multiple walls are present in teejunction.
Grid file and Physical Properties
The grid that has been generated for Caelus format is placed in the polyMesh
subdirectory of constant
. Additionally, the physical properties are specified in three different files, placed in the constant
subdirectory. The first file transportProperties
, contains the detail about the transport model for the viscosity and kinematic viscosity. The contents are as follows
/**
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] 0.003333;
We use Newtonian;
keyword as the flow is solved under Newtonian assumption, and a kinematic viscosity (\(nu\)) with the units \(m^2/s\) ([0 2 1 0 0 0 0]
) is specified.
The next file in the constant
subdirectory is the turbulenceProperties
. Here, the type of simulation through the keyword simulationType
is set to be laminar;
as shown below
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object turbulenceProperties;
}
//
simulationType laminar;
Similarly, in the RASProperties
file, RASModel
is set to laminar
.
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object RASProperties;
}
//
RASModel laminar;
turbulence off;
printCoeffs on;
Controls and Solver Attributes
The necessary files to control the simulation and specify solver attributes such as discretization method, linear solver settings can be found in the system
directory. The controlDict
file contains information regarding the simulation as shown 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 3000;
deltaT 1;
writeControl runTime;
writeInterval 1000;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression uncompressed;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
libs (
"libsimpleSwakFunctionObjects.so"
"libswakFunctionObjects.so"
"libgroovyBC.so"
);
Since groovyBC is used, few relevant libraries are imported by calling the following soon at the end of the file.
libs
(
"libsimpleSwakFunctionObjects.so"
"libswakFunctionObjects.so"
"libgroovyBC.so"
);
Next, the application simpleSolver;
referring to the SIMPLE solver is used in this simulation. As we begin the simulation at \(t = 0~s\), we need the boundary condition files to be present in the 0
directory, which has been formerly done. The keywords, startTime
to startTime
is used, where startTime
is set to a value 0
. Since simpleSolver
is a steadystate solver, the keyword endTime
corresponds to the total number of iterations.The interval at which output files are written is controlled by writeControl
and writeInterval
keywords.
The schemes for finite volume discretization are specified through fvSchemes
file with the contents as follows
/**
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 linearUpwindBJ grad(U);
div((nuEff*dev(T(grad(U))))) Gauss linear;
}
laplacianSchemes
{
default none;
laplacian(nu,U) Gauss linear corrected;
laplacian(nuEff,U) Gauss linear corrected;
laplacian(p) Gauss linear corrected;
laplacian(rAUf,p) Gauss linear corrected;
laplacian((1A(U)),p) Gauss linear corrected;
}
interpolationSchemes
{
default linear;
interpolate(HbyA) linear;
}
snGradschemes
{
default corrected;
}
As apparent from the above file, discretization schemes are set for timederivative, gradient, divergence and Laplacian terms.
In the final file, fvSolution
, linear solver settings are made 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;
}
}
SIMPLE
{
nNonOrthogonalCorrectors 0;
pRefCell 0;
pRefValue 0;
consistent true;
}
// relexation factors
relaxationFactors
{
p 1.0;
U 0.7;
}
In the above file, different linear solvers can be seen to be used to solve pressure and velocity fields. By setting the keyword consistent
to true
, SIMPLEC
solver is used and therefore a relaxation factor of 1.0
is applied for p
. Further, the grid used here is perfectly orthogonal and therefore the orthogonal correction specified via nNonOrthogonalCorrectors
is set to 0.
With these, the setup of the relevant directories and files are completed. Let us view the directory structure to ensure all are present. The tree should be identical to the one shown below
tree
.
â”śâ”€â”€ 0
â”‚Â Â â”śâ”€â”€ p
â”‚Â Â â””â”€â”€ U
â”śâ”€â”€ constant
â”‚Â Â â”śâ”€â”€ polyMesh
â”‚Â Â â”‚Â Â â”śâ”€â”€ boundary
â”‚Â Â â”‚Â Â â”śâ”€â”€ faces
â”‚Â Â â”‚Â Â â”śâ”€â”€ neighbour
â”‚Â Â â”‚Â Â â”śâ”€â”€ owner
â”‚Â Â â”‚Â Â â””â”€â”€ points
â”‚Â Â â”śâ”€â”€ RASProperties
â”‚Â Â â”śâ”€â”€ transportProperties
â”‚Â Â â””â”€â”€ turbulenceProperties
â””â”€â”€ system
â”śâ”€â”€ controlDict
â”śâ”€â”€ fvSchemes
â””â”€â”€ fvSolution
Execution of the solverÂ¶
The execution of the solver involves few different steps. The first of which is to renumber the grid or mesh followed by checking the mesh quality. Renumbering reduces the matrix bandwidth while quality check shows the mesh statistics. These can be performed as follows
caelus run  renumberMesh overwrite
caelus run  checkMesh
During the process of renumbering, gridcell bandwidth information before and after renumberMesh
is shown and the user can take a note of this. The mesh statistics are as shown below after invoking checkMesh
/**\
Caelus 8.04
Web: www.caeluscml.com
\**/
Checking geometry...
Overall domain bounding box (0 0 0.072111) (4 6 0.072111)
Mesh (nonempty, nonwedge) directions (1 1 0)
Mesh (nonempty) directions (1 1 0)
All edges aligned with or perpendicular to nonempty directions.
Boundary openness (1.07982e18 2.28423e18 1.37054e18) OK.
Max cell openness = 1.98307e16 OK.
Max aspect ratio = 2.65996 OK.
Minimum face area = 0.00127855. Maximum face area = 0.0131236. Face area magnitudes OK.
Min volume = 0.000184395. Max volume = 0.00119416. Total volume = 1.298. Cell volumes OK.
Mesh nonorthogonality Max: 0.00224404 average: 0.000335877
Nonorthogonality check OK.
Face pyramids OK.
Mesh skewness Max: 2.64853e05 average: 1.09754e06 OK.
Coupled point location match (average 0) OK.
Mesh OK.
End
From the above information, the mesh nonorthogonality is very small and therefore no nonorthogonal corrections are required for the solver to be carried out and we set nNonOrthogonalCorrectors
to 0
in the fvSolution
file. In the next step, we will execute the solver and monitor the progress of the simulation. The solver should be executed from the top level directory
caelus run l myteejunction.log  simpleSolver
The progress of the simulation is written to the log file myteejunction.log
, which can further be processed to get the convergence history. In a separate terminal window use
caelus logs w myteejunction.log
The plot indicates the convergence history for pressure and other variables with respect to number of iterations. The same for pressure is shown in Figure 9.
ResultsÂ¶
The solution obtained for the teejunction at steady state is shown here using qualitative contour plots. In Figure 10, velocity magnitude and pressure contour plots are shown. In addition, streamlines superimposed on the velocity magnitude is given. The change in the flow pattern due to the presence of side branch is quite evident from the velocity magnitude contour. The streamlines particularly facilitate to visualise the flow separation phenomenon which is occurring in this case, just before the flow entering the side branch. Also to note is the velocity profile at the inlet, which is fully developed as expected.
Circular CylinderÂ¶
The simulation of laminar, incompressible flow over a twodimensional circular cylinder is shown in this tutorial. Caelus 9.04 will be used and the details of setting up directory structure, fluid properties, boundary conditions, etc will be shown. This tutorial introduces to the user in carrying out a timedependent simulation of a externally separated flow. Further to this, the flow around the cylinder would be visualised using velocity and pressure contours.
ObjectivesÂ¶
Through this tutorial, the user would be familiar in setting up a timedependent Caelus simulation for laminar, incompressible flow in twodimensions for external separated flows. Following will be some of the steps that will be performed.
 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
Showing the flow structure in near and far wake
PrerequisitesÂ¶
It is assumed that the user is familiar with the Linux command line environment using a terminal or Caelusconsole (for Windows OS) and that Caelus is installed correctly with appropriate environment variables set. The grid used here is generated using Pointwise and the user is free to use their choice of grid generation tool having exporting capabilities to the Caelus grid format.
BackgroundÂ¶
The flow over a circular cylinder is a classical configuration to study separation and its related phenomena. This provides an ideal case for the user to use Caelus for flow over bluff bodies that represents externally separated flow. It has been shown that for low Reynolds number flows (\(40 \leq Re \leq 150\)), period vortex shedding occurs in the wake. This phenomena of vortex shedding behind bluff bodies is referred as the Karman Vortex Street [2]. The frequency associated with the oscillations of vortex streets can be expressed via Strouhal number (\(St\)) which is nondimensional relating to the frequency of oscillations (\(f\)) of vortex shedding, cylinder diameter (\(D\)) and the flow velocity (\(U\)) as
For a Reynolds number based on the cylinder diameter of \(Re_D = 100\), the Strouhal number is about \(St\approx 0.160.17\) as determined through experiments and is nearly independent of the diameter of the cylinder. More details are discussed in section Circular Cylinder
The diameter chosen for the cylinder here is \(D = 2~m\) and is the characteristic length for the Reynolds number, which is (\(Re_D = 100\)). The velocity is assumed to be \(U = 1~m/s\) in the \(x\) direction. Based on these, the kinematic velocity can be estimated as \(\nu = 0.02~m^2/s\). The below Figure 11 shows the schematic of the cylinder in the \(xy\) plane.
In the below table a summary of the freestream conditions are given
\(Re_D\) 
\(U~(m/s)\) 
\(p~(m^2/s^2)\) 
\(\nu~(m^2/s)\) 

100 
1.0 
\((0)\) Gauge 
\(0.02\) 
Grid GenerationÂ¶
A hexahedral grid around the circular cylinder was development with a Ogrid topology using Pointwise. Specific grid generation details are omitted while proving sufficient details regarding computational domain and boundary conditions. With these details the user should be able to recreate the required grid for the twodimensional cylinder
A rectangular computational domain in the \(xy\) plane has been constructed surrounding the circular cylinder as shown in Figure 12. A full cylinder was considered as the vortices developed behind the cylinder are of the periodic nature. Upstream of the cylinder, the domain is extended by 5 cylindrical diameters, whereas, downstream it was extended up to 20. Since the flow here is viscous dominated, sufficient downstream length is required to capture the initial vortex separation from the surface of the cylinder and the subsequent shedding process. In the \(y\) direction, the domain is extended up to 5 cylindrical diameters on either side. From the figure, multiple inlet boundaries to this domain can be seen, one at the far end of the upstream and the other two for the top and bottom boundaries. This type of configuration is particularly needed to appropriately model the inflow, very similar to an undisturbed flow in an experimental setup. It should be noted here that for top and bottom boundaries, the flow is in the \(x\) direction. Outlet boundary to the domain is placed at the downstream end which is at 20 cylindrical diameters. Since the fluid behaviour is viscous, the velocity at the wall is zero (\(u, v, w = 0~m/s\)) represented here through a noslip boundary condition as highlighted in blue.
The hexahedral grid around the cylinder is shown in Figure 13 for a \(xy\) plane. Caelus is a 3D solver and requires the grid to be in 3D. This was obtained by extruding the grid in the \(xy\) plane by one cell thick and subsequently specifying empty boundary conditions to the extruded planes. This enforces that Caelus to solve the flow in two dimensions assuming symmetry in the \(z\) direction as is required in this case due to the twodimensionality of the flow.
Note
A velocity value of \(w=0\) needs to be specified at appropriate boundaries although no flow is solved in the \(z\) direction.
The 2D domain consisted of 9260 cells in total. An Ogrid topology is constructed around the cylinder and extended to a maximum of 1 cylindrical diameter in the radial direction with a distribution of 10 cells. Around the cylinder, 84 cells are used giving 21 cells per each Ogrid block. Upstream of the Ogrid in the \(x\) direction, 31 cells were used and 100 in the downstream. The region of interest is about 10 diameters downstream, where the grids are refined. In the \(y\) direction, both positive and negative axes, 21 cells are used beyond the Ogrid region.
Problem definitionÂ¶
We first begin with instructions to setup the circular cylinder case in addition to the configuration files that are needed. A full working case can be found in the following directory:
/tutorials/incompressible/pimpleSolver/laminar/ACCM_circularCylinder/
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 this problem, Caelus requires time
, constant
and system
subdirectories within the main working directory. Since we start the simulation at time, t = 0 s, a time
subdirectory named 0
is required.
The 0
subdirectory contains files, p
and U
, which describe the dimensions, initialisation and boundary conditions of pressure (\(p\)) and velocity (\(U\)) respectively.
It is to be noted that Caelus is case sensitive and therefore the user should setup the directories (if applicable), files and the contents identical to what is mentioned here.
Boundary Conditions
We now begin with setting up the boundary conditions. Referring back to Figure 12, the following are the boundary conditions that need to be specified:
 Inlet
Velocity: Fixed uniform velocity \(u = 1.0~m/s\) in \(x\) direction
Pressure: Zero gradient
 Noslip wall
Velocity: Fixed uniform velocity \(u, v, w = 0\)
Pressure: Zero gradient
 Outlet
Velocity: Zero gradient velocity
Pressure: Fixed uniform gauge pressure \(p = 0\)
 Initialisation
Velocity: Fixed uniform velocity \(u = 1.0~m/s\) in \(x\) direction
Pressure: Zero Gauge pressure
Beginning with the pressure (\(p\)), the dictionary begins with FoamFile
containing standard set of keywords for version, format, location, class and object names as shown below.
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object p;
}
//
dimensions [0 2 2 0 0 0 0];
internalField uniform 0;
boundaryField
{
inlet1
{
type zeroGradient;
}
inlet2
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0;
}
symmetry
{
type empty;
}
wall
{
type zeroGradient;
}
}
The following provides the explanation to the principle entries
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 manner, the file U
contains the following entries
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
object U;
}
//
dimensions [0 1 1 0 0 0 0];
internalField uniform (1.0 0 0);
boundaryField
{
inlet1
{
type fixedValue;
value uniform (1.0 0 0);
}
inlet2
{
type fixedValue;
value uniform (1.0 0 0);
}
outlet
{
type zeroGradient;
}
symmetry
{
type empty;
}
wall
{
type fixedValue;
value uniform (0 0 0);
}
}
The principle entries for velocity field are self explanatory and the dimension are typical for 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 (1.0 0 0)
representing three components of velocity. Similarly, inlets
and wall
boundary patches have three velocity components.
Before proceeding further, it is important to ensure that the boundary conditions (inlet
, outlet
, wall
, etc) specified in the above files should be the grid boundary patches (surfaces) generated by the grid generation tool and their names are identical. Further, the two boundaries in \(xy\) plane obtained due to grid extrusion have been named as symmleft
and symmright
with specifying empty
boundary conditions forcing Caelus to assume the flow to be in twodimensions. This completes the setting up of boundary conditions.
Grid file and Physical Properties
The circular cylinder grid file is placed in the constant/polyMesh
subdirectory. Additionally, the physical properties are specified in different files present in the constant
directory.
The first file is transportProperties
, where the transport model and the kinematic viscosity are specified. The contents of this file are as follows
/**
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] 0.02;
Since the flow is Newtonian, the transportModel
is specified with Newtonian
keyword and the value of kinematic viscosity (nu
) is given which has the units \(m^2/s\) ([0 2 1 0 0 0 0]
).
The final in this sundirectory is the turbulenceProperties
file, where the type of simulation is specified as
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object turbulenceProperties;
}
//
simulationType laminar;
As the flow here is laminar, the simulationType
would be laminar
.
Controls and Solver Attributes
This section details the files required to control the simulation, specifying the type of discretization method and linear solver settings. These files can be found in the system
directory.
The first file, controlDict
is shown below
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
//
application pimpleSolver;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 360;
deltaT 0.01;
writeControl runTime;
writeInterval 1;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression uncompressed;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
//
// Function Objects to obtain mean values
functions
{
forces
{
type forces;
functionObjectLibs ("libforces.so");
patches ( wall );
CofR (0 0 0);
rhoName rhoInf;
rhoInf 1;
writeControl timeStep;
writeInterval 50;
}
}
//
In this file, the application pimpleSolver
refers to the PIMPLE solver that is used in this tutorial. We also begin the simulation at t = 0 s, which logically explains the need for 0
directory where the data files are read at the beginning of the run. Therefore, we need to set the keyword startFrom
to startTime
, where startTime
would be 0
. The simulation is run for 360 seconds specifying through the keywords stopAt
and endTime
. Since PIMPLE solver is timeaccurate, we also need to set the timestep via deltaT
. The results are written at every 0.01 seconds via writeControl
and writeInterval
keywords.
The discretization schemes and its parameters are specified in the fvSchemes
file which is shown below
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object fvSchemes;
}
//
ddtSchemes
{
default Euler;
}
gradSchemes
{
default Gauss linear;
grad(p) Gauss linear;
grad(U) Gauss linear;
}
divSchemes
{
default none;
div(phi,U) Gauss linearUpwindBJ grad(U);
div((nuEff*dev(T(grad(U))))) Gauss linear;
}
laplacianSchemes
{
default none;
laplacian(nu,U) Gauss linear corrected;
laplacian(nuEff,U) Gauss linear corrected;
laplacian(p) Gauss linear corrected;
laplacian((1A(U)),p) Gauss linear corrected;
laplacian(rAUf,p) Gauss linear corrected;
}
interpolationSchemes
{
default linear;
interpolate(HbyA) linear;
}
snGradschemes
{
default corrected;
}
In the fvSolution
file, linear solver controls and tolerances are set as shown in the below file
/**
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 1e6;
relTol 0.05;
}
pFinal
{
solver PCG;
preconditioner SSGS;
tolerance 1e7;
relTol 0;
}
UFinal
{
solver PBiCGStab;
preconditioner USGS;
tolerance 1e6;
relTol 0;
}
U
{
solver PBiCGStab;
preconditioner USGS;
tolerance 1e6;
relTol 0;
}
}
PIMPLE
{
nCorrectors 2;
nNonOrthogonalCorrectors 1;
pRefCell 0;
pRefValue 0;
}
Note that different linear solvers are used here to solve for pressure and velocity. Also, nNonOrthogonalCorrectors
is set to 1, since there is some degree of nonorthogonality in the grid.
At this stage, the directory structure should be identical to the one shown below. This can be done by using the tree
command on Linux OS.
tree
.
â”śâ”€â”€ 0
â”‚Â Â â”śâ”€â”€ p
â”‚Â Â â””â”€â”€ U
â”śâ”€â”€ constant
â”‚Â Â â”śâ”€â”€ polyMesh
â”‚Â Â â”‚Â Â â”śâ”€â”€ boundary
â”‚Â Â â”‚Â Â â”śâ”€â”€ faces
â”‚Â Â â”‚Â Â â”śâ”€â”€ neighbour
â”‚Â Â â”‚Â Â â”śâ”€â”€ owner
â”‚Â Â â”‚Â Â â”śâ”€â”€ points
â”‚Â Â â”‚Â Â â””â”€â”€ sets
â”‚Â Â â”śâ”€â”€ transportProperties
â”‚Â Â â””â”€â”€ turbulenceProperties
â””â”€â”€ system
â”śâ”€â”€ controlDict
â”śâ”€â”€ fvSchemes
â””â”€â”€ fvSolution
Execution of the solverÂ¶
Prior to solver execution, renumbering of the grid or mesh needs to be performed as well as checking the quality of the grid. Renumbering the gridcell list is vital to reduce the matrix bandwidth while the 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
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.
In the next step, execution of the solver can be performed while monitoring the progress of the solution. The solver is always executed from the top directory which is ACCM_circularCylinder
in this case.
caelus run l mycircularcylinder.log  pimpleSolver
The output of the solution process is saved in the log file, mycircularcylinder.log
. In a separate terminal window the convergence history can be monitored using
caelus logs w mycircularcylinder.log
With the above, the convergence of pressure along with other variables can be seen with respect to time. The same is shown in the Figure 14 and due to the periodic nature of vortex shedding, oscillatory convergence of pressure is seen.
ResultsÂ¶
In this section, some qualitative results obtained as a result of steady vortex shedding in the wake of the cylinder is shown. Figure 15 shows the velocity magnitude and pressure contour for the flow over the cylinder. The formation of vortex shedding is clearly visible from the velocity contour and the pressure variation due to oscillating vortex in the wake. The vortex breakup that occurs in the near wake of the cylinder, travels several diameters downstream eventually diffusing into the flow.
Triangular CavityÂ¶
In this tutorial, we carry out laminar, incompressible flow inside a triangular cavity in twodimensions using Caelus 9.04. Details regarding setting up of the directories, fluid properties, boundary conditions, etc will be discussed. Subsequent to this, postprocessing the velocity distribution along the centerline will be shown.
ObjectivesÂ¶
With the completion of this tutorial, the user would be familiar with setting up a steadystate Caelus simulation for laminar, incompressible flow over lipdriven cavity. Following are the steps that would be performed
 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
Showing the flow structure inside the cavity
PrerequisitesÂ¶
It is assumed that the user is familiar with the Linux command line environment using a terminal or Caelusconsole (for Windows OS) and that Caelus is installed correctly with appropriate environment variables set. The grid used here is generated using Pointwise and the user is free to use their choice of grid generation tool having exporting capabilities to the Caelus grid format.
BackgroundÂ¶
Flow inside liddriven cavities is a classical case to study cases with flow recirculation. In the present case, the top wall of the cavity moves at constant velocity initiating a recirculation motion with the cavity and as a consequence, a boundary layer develops in the direction of the moving lid. The feature that is of interest is the velocity distribution along the centerline of cavity. Details regarding the validation of this case is given in Triangular Cavity.
The triangular cavity schematic is shown in Figure 16. Here D represents the cavity depth which is 4 m and the top width, W = 2 m. For this configuration, the Reynolds number based on the cavity depth is 800 and the wall velocity is assumed and set to 2 m/s. This give us with a kinematic viscosity of 0.01. Note that the twodimensional plane considered here is in \(xz\).
The below table gives the summary of the freestream conditions used here
\(Re_D\) 
\(U~(m/s)\) 
\(p~(Pa)\) 
\(\nu~(m^2/s)\) 

800 
2.0 
\((0)\) Gauge 
\(0.01\) 
Grid GenerationÂ¶
A hybridgrid consisting of quadrilateral and triangular cells has been generated for this cavity geometry using Pointwise. Details regarding the generation of grid is not covered in this tutorial, however details regarding computational domain and boundary conditions are provided.
The computational domain for the triangular cavity follows the cavity geometry due to internal flow configuration. This is in contrast to other flow configurations here where the flow was over the region of interest. A schematic of the domain is shown in Figure 17. The velocity at the cavity walls (high lighted in blue) is zero, represented through a noslip boundary, wherein \(u, v, w = 0\). Whereas the top wall has a uniform velocity in the xdirection.
The hybrid grid is shown in Figure 18. As can be seen, up to a depth of D = 1.35 m, structured grids are used and after which it is filled with triangular unstructured elements. In the structured domain, 40 X 40 cells are used respectively. In the 2D domain, a total of 5538 cells are present, however the polyMesh
format of Caelus should be in 3D. This was achieved by extruding the grid in the \(xy\) plane by one cell thick and subsequently specifying empty boundary conditions to the extruded planes. This should force Caelus to solve the flow the flow in 2D in the extruded direction, which is \(z\).
Note
A velocity value of \(w=0\) needs to be specified at appropriate boundaries although no flow is solved in the \(z\) direction.
Problem definitionÂ¶
This section provides the case setup procedures along with the configuration files that are needed. A full working case of this problem is given in the following location:
/tutorials/incompressible/pimpleSolver/laminar/ACCM_triangularCavity/
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
Caelus requires time
, constant
and system
subdirectories within the main mytriangularcavity
working directory. Since we start the simulation at time, t = 0 s, a time
subdirectory named 0
is required.
The 0
subdirectory contains the pressure, p
and velocity U
files. The contents of these files set the dimensions, initialisation and boundary conditions to the case, which form the three principle entries required.
If applicable, the user should take precautions in setting the directories and files as Caelus is case sensitive. These have to be identical to the names mentioned here.
Boundary Conditions
Next we start with settingup of the boundary conditions. Referring back to Figure 17, the following are the boundary conditions that need to be specified:
 Moving wall
Velocity: Fixed uniform velocity \(u = 2.0~m/s\) in \(x\) direction
Pressure: Zero gradient
 Noslip wall
Velocity: Fixed uniform velocity \(u, v, w = 0\)
Pressure: Zero gradient
 Initialisation
Velocity: Fixed uniform velocity \(u = 0~m/s\) in \(x, y, z\) directions
Pressure: Zero Gauge pressure
The file p
for pressure contains the following information
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object p;
}
//
dimensions [0 2 2 0 0 0 0];
internalField uniform 0;
boundaryField
{
fixedwalls
{
type zeroGradient;
}
movingwall
{
type zeroGradient;
}
symm
{
type empty;
}
}
As noted from the above file, the dictionary begins with FoamFile
containing standard set of keywords for version, format, location, class and object names. The following provides the explanation to the principle entries
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 shown below
/**
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 (0 0 0);
boundaryField
{
fixedwalls
{
type fixedValue;
value uniform (0 0 0);
}
movingwall
{
type fixedValue;
value uniform (2 0 0);
}
symm
{
type empty;
}
}
The principle entries for velocity field are self explanatory and the dimension are typical for velocity with the units \(m/s\) ([0 1 1 0 0 0 0]
). Here, since the top wall moves with a velocity, we set a uniform (2.0 0 0)
for movingwall
boundary patch. Similarly, the cavity walls (fixedwalls
) have uniform (1.0 0 0)
.
At this stage, the user should ensure that the boundary conditions (fixedwalls
, movingwall
and symm
) specified in the above files should be the grid boundary patches (surfaces) generated by the grid generation tool and their names are identical. Further, the two boundaries in \(xy\) plane obtained due to grid extrusion have been combined and named as symm
with specifying empty
boundary conditions forcing Caelus to assume the flow to be in twodimensions. With this, the setting up of boundary conditions are completed.
Grid file and Physical Properties
The triangular cavity grid is placed in constant/polyMesh
subdirectory. In addition, the physical properties are specified in different files, all present in the constant
directory.
The transport model and the kinematic viscosity are specified in the file transportProperties
. The contents of this file are as follows
/**
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] 0.01;
Since the viscous behaviour is modelled as Newtonian, the transportModel
is specified with the keyword Newtonian
and the value of kinematic viscosity is set with has the units \(m^2/s\) ([0 2 1 0 0 0 0]
).
The final file in this class is the turbulenceProperties
file, in which the type of simulation is specified as
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object turbulenceProperties;
}
//
simulationType laminar;
The flow being laminar, the simulationType
is set to laminar
.
Controls and Solver Attributes
In this section, the files required to control the simulation, setting the discretization parameters and linear solver settings are discussed can be found in the system
directory.
The controlDict
file is shown below.
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
//
application pimpleSolver;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 60;
deltaT 0.01;
writeControl runTime;
writeInterval 1;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression uncompressed;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
In the controlDict
file, pimpleSolver
refers to the application, PIMPLE solver that is used here. The simulation is also started at t = 0 s and this logically explains the need for 0
directory where the data files are read at the beginning of the run. Therefore, the keyword startFrom
to startTime
, where startTime
would be 0
is needed. The simulation is run for 200 seconds specifying through the keywords stopAt
and endTime
. Since PIMPLE solver is timeaccurate, we also need to set the timestep via deltaT
. The results are written at every 0.01 seconds via writeControl
and writeInterval
keywords.
The discretization schemes and its parameters are specified in the fvSchemes
file which is shown below
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object fvSchemes;
}
//
ddtSchemes
{
default Euler;
}
gradSchemes
{
default Gauss linear;
grad(p) Gauss linear;
grad(U) Gauss linear;
}
divSchemes
{
default none;
div(phi,U) Gauss linearUpwindBJ grad(U);
div((nuEff*dev(T(grad(U))))) Gauss linear;
}
laplacianSchemes
{
default none;
laplacian(nu,U) Gauss linear corrected;
laplacian(nuEff,U) Gauss linear corrected;
laplacian(p) Gauss linear corrected;
laplacian((1A(U)),p) Gauss linear corrected;
laplacian(rAUf,p) Gauss linear corrected;
}
interpolationSchemes
{
default linear;
interpolate(HbyA) linear;
}
snGradschemes
{
default corrected;
}
In the following file, which is the fvSolution
the linear solver controls and tolerances are set as shown 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 1e10;
relTol 0;
}
pFinal
{
solver PCG;
preconditioner SSGS;
tolerance 1e10;
relTol 0;
}
UFinal
{
solver PBiCGStab;
preconditioner USGS;
tolerance 1e15;
relTol 0;
}
U
{
solver PBiCGStab;
preconditioner USGS;
tolerance 1e15;
relTol 0;
}
}
PIMPLE
{
nCorrectors 2;
nNonOrthogonalCorrectors 1;
pRefCell 0;
pRefValue 0;
}
It should be noted that different linear solvers are used to solve for pressure a velocity. Since we have used hybrid grids, nNonOrthogonalCorrectors
is set to 1
as there would be some degree of nonorthogonality present due to triangular nature of the geometry.
This completes the setup of directory structure along with all the necessary files. This can be verified by using the following commands and the directory tree should be identical to the one shown below
tree
.
â”śâ”€â”€ 0
â”‚Â Â â”śâ”€â”€ p
â”‚Â Â â””â”€â”€ U
â”śâ”€â”€ constant
â”‚Â Â â”śâ”€â”€ polyMesh
â”‚Â Â â”‚Â Â â”śâ”€â”€ boundary
â”‚Â Â â”‚Â Â â”śâ”€â”€ faces
â”‚Â Â â”‚Â Â â”śâ”€â”€ neighbour
â”‚Â Â â”‚Â Â â”śâ”€â”€ owner
â”‚Â Â â”‚Â Â â””â”€â”€ points
â”‚Â Â â”śâ”€â”€ transportProperties
â”‚Â Â â””â”€â”€ turbulenceProperties
â””â”€â”€ system
â”śâ”€â”€ controlDict
â”śâ”€â”€ fvSchemes
â”śâ”€â”€ fvSolution
Execution of the solverÂ¶
Renumbering and checking the mesh quality is needed before the solver is executed. Renumbering the gridcell list is vital to reduce the matrix bandwidth while the 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
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 (1 4 0) (1 0 0.0447214)
Mesh (nonempty, nonwedge) directions (1 1 0)
Mesh (nonempty) directions (1 1 0)
All edges aligned with or perpendicular to nonempty directions.
Boundary openness (6.25534e18 5.30678e18 6.4712e16) OK.
Max cell openness = 2.16191e16 OK.
Max aspect ratio = 2.97779 OK.
Minimum face area = 8.00174e05. Maximum face area = 0.00322809. Face area magnitudes OK.
Min volume = 3.57849e06. Max volume = 0.000109372. Total volume = 0.178886. Cell volumes OK.
Mesh nonorthogonality Max: 41.8682 average: 9.34065
Nonorthogonality check OK.
Face pyramids OK.
Mesh skewness Max: 0.621726 average: 0.0313562 OK.
Coupled point location match (average 0) OK.
Mesh OK.
End
Solver can now be executed and the progress of the solution can be monitored. The solver is always executed from the top directory which is mytriangularcavity
in this case.
caelus run l mytriangularcavity.log pimpleSolver
The output of the solver progress is saved in the log file, mytriangularcavity.log
. The log file (mytriangularcavity.log
) can be monitored to look at the convergence history. In a separate terminal window use
caelus logs w mytriangularcavity.log
The convergence of the pressure can now be seen with respect to time.
Spherical CavityÂ¶
In this tutorial, we look at simulating buoyant flow inside a spherically heated cavity using Caelus 9.04 which will be a threedimension case. Because of the natural convection process, the fluid initiates a motion due to buoyancy effects. This results in steady isotherms which is of interest in this tutorial and will be compared with the analytical data.
ObjectivesÂ¶
With this tutorial, the user would be familiar in setting up a steadystate Caelus simulation for laminar, buoyant flow within a spherically heated cavity. The steps that would be performed in this tutorial 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
Showing the flow structure inside the cavity
PrerequisitesÂ¶
It is assumed that the user is familiar with the Linux command line environment using a terminal or Caelusconsole (for Windows OS) and that Caelus is installed correctly with appropriate environment variables set. The grid used here is generated using Pointwise and the user is free to use their choice of grid generation tool having exporting capabilities to the Caelus grid format.
BackgroundÂ¶
The flow inside a spherically heated cavity is an interesting case of buoyant flow simulation. Here, the flow is enclosed in a spherical cavity and the wall is heated with a specified temperature gradient. Due to the natural convection process, the fluid initiates a motion as a result of buoyancy effects. The characteristic feature that is of interest is the steady isotherms. Details regarding the validation of this case is given in Spherical Cavity.
The schematic representation of the sphere is shown in Figure 21 and only half of the sphere is considered here having the plane of symmetry in #:math:xy plane at \(z=0\). The radius of the sphere is chosen to be \(r=0.5~m\), such that \(x=0\) at \(r=0\) and \(x=0.5\) at \(r=0.5\).
The thermal boundary condition to the spherical wall was generated by specifying temperature (\(T\)) as a function of distance (\(x\)), which are as follows:
with this,
To obtain a smooth distribution of temperature, sufficient cells in the radial direction is required. For buoyancy driven flows, the nondimensional number, Rayleigh number (\(Ra\)) is often used which relates buoyancy and viscous effects of the flow. In this case, we will be simulating with a (\(Ra\)) number of 2000 and a Prandtl number of 0.7.
The below table summarises the flow conditions that will be used
\(Ra\) 
\(Pr\) 
\(T~(K)\) 
\(p~(m^2/s^2)\) 
\(\nu~(m^2/s)\) 
\(\beta~(1/K)\) 

\(2000\) 
\(0.7\) 
\(T = x\) 
\((0)\) Gauge 
\(3.4 \times 10^{4}\) 
\(3.567\times10^{5}\) 
Grid GenerationÂ¶
A hexahedral grid within the sphere was developed via the Ogrid topology using Pointwise. Grid generation details are not discussed here, however details regarding the computational domain and boundary conditions are provided.
As noted earlier, the computational domain considered is a half sphere with the plane of symmetry in the \(xy\) axes at \(z=0~m\). The temperature will be set as given in the above equations and also the initialisation follows the surface temperature (\(T=x\)). In Figure 22 the computational domain, in particular the applied temperature boundary condition when applied is shown.
The structured Ogrid over the spherical wall and on the plane of symmetry is shown in Figure 23. A total of 18564 cells are present within the domain. Over the plane of symmetry, 32 cells are distributed and, over the spherical surface 35 cells are distributed in the radial direction.
Problem definitionÂ¶
We begin with instructions to setup the spherical cavity case and the required configuration files that are needed. A full working case is found in the following directory:
/tutorials/heatTransfer/buoyantBoussinesqSimpleSolver/laminar/ACCM_sphericalCavity/
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 a time
, constant
and system
subdirectories. Since we will be starting the simulation at time \(t = 0~s\), the time
directory would be just 0
.
In the 0
subdirectory, few additional files p
, p_rgh
, alphat
, T
, and U
for pressure, buoyant pressure, turbulent thermal conductivity, temperature and velocity are set respectively. Note that even for a laminar simulation, alphat
is required, although a value of 0
would be used. The contents of these files set the dimensions, initialisation and boundary conditions to the case.
The user should be aware that Caelus is case sensitive and therefore the setting up of the directories, files and contents (when applicable) should be identical.
Boundary Conditions and Solver Attributes
Boundary Conditions
The boundary conditions would be set as follows
 Wall
Velocity: Fixed uniform velocity \(u, v, w = 0\)
Pressure: Uniform zero Buoyant Pressure
Temperature: Temperature as a function of \(x\) (\(T = x\))
Turbulent thermal conductivity: Fixed uniform value of 0
 Symmetry Plane
Velocity: Symmetry
Pressure: Symmetry
Temperature: Symmetry
Turbulent thermal conductivity: Symmetry
 Initialisation
Velocity: Fixed uniform velocity \(u, v, w = 0\)
Pressure: Uniform zero Buoyant Pressure
Temperature: Temperature as a function of \(x\) (\(T = x\))
Turbulent thermal conductivity: A value of 0
Note that to specify a temperature as a function of \(x\), swak4Foam library was used and funkySetFields
utility was employed. The usage of this would be shown later.
The first quantity to define would be the pressure (\(p\)) and this is done in the file p
, 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
{
symm
{
type symmetryPlane;
}
wall1
{
type zeroGradient;
}
wall2
{
type zeroGradient;
}
}
// ************************************************************************* //
The above file begins with a dictionary named FoamFile
containing standard set of keywords such as version, format, location, class and object names. The principle elements follows next
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 input data for the buoyant pressure is shown below
/**\
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object p_rgh;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 2 0 0 0 0];
internalField uniform 0;
boundaryField
{
symm
{
type symmetryPlane;
}
wall1
{
type fixedFluxPressure;
gradient uniform 0;
rho rhok;
value uniform 0;
}
wall2
{
type fixedFluxPressure;
gradient uniform 0;
rho rhok;
value uniform 0;
}
}
// ************************************************************************* //
In the next file, alphat
is defined as follows.
/**\
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object alphat;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 1 0 0 0 0];
internalField uniform 0;
boundaryField
{
symm
{
type symmetryPlane;
}
wall1
{
type fixedValue;
value uniform 0;
}
wall2
{
type fixedValue;
value uniform 0;
}
}
// ************************************************************************* //
For the temperature, before the funkySetFields
expression is used, we need to have the T
file with the generic boundary conditions. This is because when funkySetFields
is applied, the T
file would be overwritten with the profile values.
The final file in this directory is the velocity and the contents are 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 (0 0 0);
boundaryField
{
symm
{
type symmetryPlane;
}
wall1
{
type fixedValue;
value uniform (0 0 0);
}
wall2
{
type fixedValue;
value uniform (0 0 0);
}
}
// ************************************************************************* //
As noted above, the principle entries for the velocity filed are self explanatory with the typical dimensional units of \(m/s\) ([0 1 1 0 0 0 0]
). The initialisation of the flow is done at \(0~m/s\) which is set using internalField
to uniform (0 0 0);
which represents three components of velocity.
The boundary conditions (symm
, wall1
, wall2
, etc) specified in the above files should be the grid boundary patches (surfaces) that are generated by the gridgeneration tool. The user should ensure that their names are identically matched.
Grid file and Physical Properties
The spherical cavity grid is placed in constant/polyMesh
subdirectory. Additionally, physical properties are specified in separate files, placed in the constant
subdirectory.
The first file is g
,wherein the value of acceleration due to gravity is specified. Since we are simulating a buoyancy driven flow, gravity effects should be considered. The contents of g
are as follows
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class uniformDimensionedVectorField;
location "constant";
object g;
}
// 
dimensions [0 1 2 0 0 0 0];
value ( 0 9.81 0 );
In the next file transportProperties
, details about the transport model for viscosity, Prandtl number, coefficient of thermal expansion are specified. The contents are as follows
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
//
transportModel Newtonian;
// Laminar viscosity
nu nu [0 2 1 0 0 0 0] 3.4e4;
// Thermal expansion coefficient
beta beta [0 0 0 1 0 0 0] 3.5677e5;
// Reference temperature
TRef TRef [0 0 0 1 0 0 0] 0;
// Laminar Prandtl number
Pr Pr [0 0 0 0 0 0 0] 0.7;
// Turbulent Prandtl number
Prt Prt [0 0 0 0 0 0 0] 0.7;
Here, Newtonian;
keyword is used since the flow is under Newtonian assumption and a kinematic viscosity (\(nu\)) with the units \(m^2/s\) ([0 2 1 0 0 0 0]
) is specified. Similarly, the value of beta
, TRef
, Pr
and Prt
are specified. Note that a value of 0.7 is used for laminar Prandtl number.
The type of simulation is specified in turbulenceProperties
to laminar
using the keyword simulationType
as shown below
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object turbulenceProperties;
}
//
simulationType laminar;
In addition to this, RASProperties
file is required since we specify turbulent thermal conductivity property.
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object RASProperties;
}
//
RASModel laminar;
turbulence off;
printCoeffs on;
Controls and Solver Attributes
Here we setup the necessary files to control the simulation and specifying the solver attributes such as discretization method, linear solver setting, etc and these reside in the system
directory.
First, there is the controlDict
file with the following information
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
//
application buoyantBoussinesqSimpleSolver;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 2000;
deltaT 1.0;
writeControl runTime;
writeInterval 10;
purgeWrite 5;
writeFormat ascii;
writePrecision 12;
writeCompression uncompressed;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
libs (
"libsimpleSwakFunctionObjects.so"
"libswakFunctionObjects.so"
"libgroovyBC.so"
);
As we will be using funkySetFields
, few relevant libraries are needed and they are imported by using the following lines in the controlDict
.
libs
(
"libsimpleSwakFunctionObjects.so"
"libswakFunctionObjects.so"
"libgroovyBC.so"
);
The next is the application, buoyantBoussinesqSimpleSolver
which is the buoyancy version of the SIMPLE
solver that will be used in this simulation. As we begin the simulation at \(t = 0~s\), we need the boundary condition files to be present in the 0
directory, which has been formerly done. The keywords, startTime
to startTime
is used, where startTime
is set to a value 0
. The flow is simulated for a total of 2000 iterations and is specified through stopAt
and endTime
. This is followed by setting the time interval of 10
iterations to save the results via writeControl
and writeInterval
keywords.
The schemes for finite volume discretization are specified through fvSchemes
file with the contents as follows
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSchemes;
}
//
ddtSchemes
{
default Euler;
}
gradSchemes
{
default Gauss linear;
grad(p) Gauss linear;
grad(U) Gauss linear;
}
divSchemes
{
default none;
div(phi,U) Gauss linearUpwindBJ grad(U);
div(phi,T) Gauss linearUpwindBJ grad(T);
div((nuEff*dev(T(grad(U))))) Gauss linear;
}
laplacianSchemes
{
default none;
laplacian(rAUf,p_rgh) Gauss linear corrected;
laplacian(nu,U) Gauss linear corrected;
laplacian(nuEff,U) Gauss linear corrected;
laplacian(alphaEff,T) Gauss linear corrected;
laplacian(p_rgh) Gauss linear corrected;
laplacian(kappaEff,T) Gauss linear corrected;
laplacian((1A(U)),p) Gauss linear corrected;
laplacian((1A(U)),p_rgh) Gauss linear corrected;
}
interpolationSchemes
{
default linear;
interpolate(HbyA) linear;
}
snGradschemes
{
default corrected;
}
In the above file, the discretization schemes are set for timederivative, gradient, divergence and Laplacian terms.
The final file is the fvSolution
, where linear solver settings are provided and is as given below
/**
Caelus 9.04
Web: www.caeluscml.com
\**/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSolution;
}
//
solvers
{
p_rgh
{
solver GAMG;
tolerance 1e8;
relTol 0.01;
smoother SSGS;
cacheAgglomeration true;
nCellsInCoarsestLevel 10;
agglomerator algebraicPair;
mergeLevels 1;
}
"(UT)"
{
solver PBiCGStab;
preconditioner USGS;
tolerance 1e8;
relTol 0.01;
}
}
SIMPLE
{
nNonOrthogonalCorrectors 0;
pRefCell 0;
pRefValue 0;
residualControl
{
p_rgh 1e4;
T 1e4;
// possibly check turbulence fields
}
}
relaxationFactors
{
fields
{
p_rgh 0.3;
}
equations
{
U 0.7;
T 0.7;
}
}
// ************************************************************************* //
The above file details the different linear solvers that are used to solve for buoyant pressure and velocity fields. Further, nNonOrthogonalCorrectors
is set to 2;
since there is some degree of nonorthogonality present in the grid.
With these, the setup of the relevant directories and files are completed. Let us view the directory structure to ensure all are present. The tree should be identical to the one shown below
tree
.
â”śâ”€â”€ 0
â”‚ â”śâ”€â”€ alphat
â”‚ â”śâ”€â”€ p
â”‚ â”śâ”€â”€ p_rgh
â”‚ â”śâ”€â”€ T
â”‚ â””â”€â”€ U
â”śâ”€â”€ constant
â”‚ â”śâ”€â”€ polyMesh
â”‚ â”‚ â”śâ”€â”€ boundary
â”‚ â”‚ â”śâ”€â”€ faces
â”‚ â”‚ â”śâ”€â”€ neighbour
â”‚ â”‚ â”śâ”€â”€ owner
â”‚ â”‚ â””â”€â”€ points
â”‚ â”śâ”€â”€ RASProperties
â”‚ â”śâ”€â”€ transportProperties
â”‚ â””â”€â”€ turbulenceProperties
â””â”€â”€ system
â”śâ”€â”€ controlDict
â”śâ”€â”€ fvSchemes
â””â”€â”€ fvSolution
Execution of the solverÂ¶
Before the solver can be executed, we have to apply the temperature profile as noted earlier. This step can be performed as follows
funkySetFields field T expression "pos().x" time 0 keepPatches valuePatches "wall1 wall2"
The above expression rewrites the T
file with the appropriate values of temperature as a function of distance \(x\). Following this, the solver can be executed and it involves few different steps. The first of which is to renumber the grid or mesh followed by checking the mesh quality. Renumbering reduces the matrix bandwidth while quality check shows the mesh statistics. These can be performed as follows
caelus run  renumberMesh overwrite
caelus run  checkMesh
During the process of renumbering, gridcell bandwidth information before and after renumberMesh
is shown and the user can take a note of this. The mesh statistics are as shown below after invoking checkMesh
/**\
Caelus 8.04
Web: www.caeluscml.com
\**/
Checking geometry...
Overall domain bounding box (0.5 0.5 0) (0.5 0.5 0.5)
Mesh (nonempty, nonwedge) directions (1 1 1)
Mesh (nonempty) directions (1 1 1)
Boundary openness (2.69541574409e17 1.48776322194e17 2.09950188835e16) OK.
Max cell openness = 2.71077471766e16 OK.
Max aspect ratio = 5.13298999541 OK.
Minimum face area = 7.26421673269e05. Maximum face area = 0.0145134277912. Face area magnitudes OK.
Min volume = 8.32359098625e07. Max volume = 0.000256180182954. Total volume = 0.259143943239. Cell volumes OK.
Mesh nonorthogonality Max: 64.2010673845 average: 15.1363376112
Nonorthogonality check OK.
Face pyramids OK.
Mesh skewness Max: 0.60153646289 average: 0.0467687348165 OK.
Coupled point location match (average 0) OK.
Mesh OK.
End
The above information gives a maximum mesh nonorthogonality angle of 64.2 and therefore nonorthogonal corrections are needed for the solver. In the next step, we will execute the solver and monitor the progress of the simulation. The solver should be executed from the top level directory using,
caelus run l mysphericalcavity.log buoyantBoussinesqSimpleSolver
The progress of the simulation is written to the log file mysphericalcavity.log
, which can be monitored to view the convergence history. In a separate terminal window use
caelus logs w mysphericalcavity.log
The plot indicates the convergence history for pressure with respect to time and a similar plot is shown in Figure 24. The convergence of other properties can also be found by using the above command.
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.
ResultsÂ¶
The flow visualisation over the bump is shown here through the contours of velocity and pressure for SA model. In Figure 35 the variation of velocity and pressure can be seen as the bump is approached. As expected due to the curvature, flow accelerates while pressure reduces over the bump.
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.
Validation ManualÂ¶
Validation: IntroductionÂ¶
This part contains details of validation case for Caelus.
Validation: Incompressible Laminar FlowsÂ¶
This section details the validation work carried out in Caelus 9.04 for several conditions ranging from attached to separated flows. Both transient and steady state conditions have been considered under laminar and turbulent flows for relevant cases.
FlatplateÂ¶
Laminar, incompressible flow over a twodimensional SharpLeading Edge flatplate
Nomenclature
Symbol 
Definition 
Units (SI) 

\(c_f\) 
Skin friction coefficient 
Nondimensional 
\(e\) (subscript) 
Boundary layer edge conditions 

\(L\) 
Length of the plate 
\(m\) 
\(p\) 
Kinematic pressure 
\(Pa/\rho~(m^2/s^2)\) 
\(Re\) 
Reynolds number 
Nondimensional 
\(T\) 
Temperature 
\(K\) 
\(u\) 
Local velocity in xdirection 
\(m/s\) 
\(U\) 
Freestream velocity in xdirection 
\(m/s\) 
\(x\) 
Distance in xdirection 
\(m\) 
\(y\) 
Distance in ydirection 
\(m\) 
\(\nu\) 
Kinematic viscosity 
\(m^2/s\) 
\(\eta\) 
Blasius parameter 
Nondimensional 
Introduction
In this validation case, steady, incompressible, laminar flow over a twodimensional sharpleading edge flatplate at zero angle of incidence is investigated. The flow generates a laminar boundary layer and the computational results are compared with the Blasius solution for the incompressible flow. Validation of the flow over a flatplate forms the basis of these validation efforts as it is perhaps the most well understood configuration for CFD codevalidation.
Blasius, in his work [6] obtained the solution to the Boundary Layer Equations using a transformation technique. Here, equations of continuity and momentum in twodimensional form are converted into a single ordinary differential equation (ODE). The solution to this ODE can be numerically obtained and is regarded as the exact solution to the boundary layer equations. It is only valid for steady, incompressible, laminar flow over a flatplate. One of the highlights of Blasius solution is the analytical expression for the skin friction coefficient (\(c_f\)) distribution along the flatplate given by
where \(Re_{x}\) is the local Reynolds number defined as
\(U\) is the freestream velocity, \(x\) is the distance starting from the leading edge and \(\nu\) is the kinematic viscosity.
Problem definition
This exercise is based on the validation work carried out by the NASA NPARC alliance for a flow over a flatplate using the same conditions in the incompressible limit. A schematic of the geometric configuration is shown in Figure 52.
The length of the plate is L = 0.3048 m wherein, x = 0 is the leading edge, the Reynolds number of the flow based on the length of the plate is 200,000 and U is the velocity in the xdirection. Assuming the inlet flow is at a temperature of 300 K, the kinematic viscosity can be determined from dynamic viscosity and density of the fluid. The value is given in the table below. The value of dynamic viscosity is obtained from the Sutherland viscosity formulation. Using the Reynolds number, plate length and kinematic viscosity, the freestream velocity evaluates to U = 10.4306 m/s. The following table summarises the freestream conditions used:
Fluid 
\(Re\) 
\(U~(m/s)\) 
\(p~(m^2/s^2)\) 
\(T~(K)\) 
\(\nu~(m^2/s)\) 

Air 
200,000 
10.43064 
\((0)\) Gauge 
300 
\(1.58963\times10^{5}\) 
As we have assumed the flow incompressible, the density (\(\rho\)) remains constant. In addition, since the fluid temperature is not considered, the viscosity remains constant. For incompressible flows, the kinematic forms of pressure and viscosity are always used in Caelus.
Computational Domain and Boundary Conditions
The computational domain is a rectangular block encompassing the flatplate. Figure 53 shows the details of the boundaries used in twodimensions (\(xy\) plane). As indicated in blue, the region of interest extends between \(0\leq x \leq 0.3048~m\) and has a noslip boundary condition. Upstream of the leading edge, a slip boundary is used to simulate freestream uniform flow approaching the flatplate. However, downstream of the plate, there is additional noslip wall a further three plate lengths (highlighted in green). This ensures that the boundary layer in the vicinity of the trailing edge is not influenced by the outlet boundary. Since the flow is subsonic, disturbances cause the pressure to propagate both upstream and downstream. Therefore, placement of the inlet and outlet boundaries were chosen to have minimal effect on the solution. The inlet boundary as shown in the figure below is placed at start of the slipwall (\(x = 0.06~m\)) and the outlet at the end of the second noslip wall (\(x = 1.2192~m\)). Both inlet and outlet boundaries are between \(0\leq y \leq 0.15~m\). A slipwall condition is used for the entire top boundary.
Boundary Conditions and Initialisation
The following are the boundary condition details used for the computational domain:
 Inlet
Velocity: Fixed uniform velocity \(u = 10.4306~m/s\) in \(x\) direction
Pressure: Zero gradient
 Slip wall
Velocity: Slip
Pressure: Slip
 Noslip wall
Velocity: Fixed uniform velocity \(u, v, w = 0\)
Pressure: Zero gradient
 Outlet
Velocity: Zero gradient velocity
Pressure: Fixed uniform gauge pressure \(p = 0\)
 Initialisation
Velocity: Fixed uniform velocity \(u = 10.4306~m/s\) in \(x\) direction
Pressure: Zero Gauge pressure
Computational Grid
A 2D structured grid was generated using Pointwise in the \(xy\) plane. Since Caelus is a 3D computational framework, it necessitates the grid to also be 3D. Therefore, a 3D grid was obtained using Pointwise by extruding the 2D grid in the positive \(z\) direction by one cell. The final 3D grid was then exported to the Caelus format (polyMesh). The two \(xy\) planes obtained as a result of grid extrusion need boundary conditions to be specified. As the flow over a flatplate is generally 2D, we do not need to solve the flow in the third dimension. This is achieved in Caelus by specifying empty boundary condition for each plane. Although, no flow is computed in the \(z\) direction, a velocity of \(w = 0\) has to be specified for the velocity boundary condition as indicated above.
Figure 54 shows the 2D grid in the \(xy\) plane. As can be seen, the grid is refined perpendicular to the wall in order to capture resolve the viscous effects. To ensure that the gradients in boundary layer are well resolved, about 50 grid nodes are placed between the wall and the boundary layer edge. Grid refinement is also added at the leading edge so that the growth of the boundary layer is also well resolved. In this particular case, 400 cells were used in the streamwise (\(x\)) direction (\(x \leq 0 \leq 0.3048~m\)) and 600 in the wall normal (\(y\)) direction. For noslip wall beyond \(x > 0.3048\), a similar distribution is used.
Results and Discussion
A timedependent solution to the twodimensional flatplate was obtained using Caelus 9.04. The SLIM transient solver was used here and the flow was simulated sufficiently long (several plate length flow times) such that steady flow was established. For the discretization of timedependent terms, the firstorder Euler scheme was used. A Gauss linear discretization was used for the pressure and velocity gradients. A linear upwind discretization was for the divergence of velocity and mass flux. A linear corrected scheme was used for Laplacian discretization while celltoface centre interpolation used linear interpolation.
In Figure 55, the skinfriction distribution along the flatplate obtained from the CFD simulation is compared with that of the Blasius analytical solution. Here, the distance \(x\) is normalised with the length of the plate (\(L\)). Excellent agreement is observed along the entire length of the flatplate.
At the exit plane of the flatplate at \(x = 0.3048~m\), velocity data was extracted across the boundary layer and compared with the Blasius analytical solution. This is shown in Figure 56 where the velocity profile is plotted using similarity variables from the Blasius solution. Here, \(\eta\) is the nondimensional distance from the wall to the boundary layer edge and \(U_e\) is the velocity at the boundary layer edge. Similar to skinfriction, the velocity profile also exhibits excellent agreement with the Blasius solution.
Conclusions
The steady, incompressible, laminar flow over a twodimensional flatplate was simulated using Caelus 9.04 utilising the SLIM solver. The results were validated against the Blasius analytical solutions resulting in excellent agreement.
Tee JunctionÂ¶
Laminar, incompressible flow in a 90 degree tee junction
Nomenclature
Symbol 
Definition 
Units (SI) 

\(L\) 
Length of the branch 
\(m\) 
\(p\) 
Kinematic pressure 
\(Pa/\rho~(m^2/s^2)\) 
\(Re_w\) 
Reynolds number based on width 
Nondimensional 
\(V\) 
Freestream velocity in ydirection 
\(m/s\) 
\(W\) 
Width of the branch 
\(m\) 
\(x\) 
Distance in xdirection 
\(m\) 
\(y\) 
Distance in ydirection 
\(m\) 
\(\nu\) 
Kinematic viscosity 
\(m^2/s\) 
Introduction
In this validation case, laminar, incompressible flow through a twodimensional \(90^\circ\) tee junction was investigated. Due to the presence of the side branch, the flow separates and forms a recirculating region. The recirculating regions influences the mass flow through the main and side branches. The numerically computed massflow ratio was calculated and compared with experiment.
A comprehensive study of flow through planar branches has been carried out by Hayes et al. [12] due to its prevelance in the biomechanical industry. Of which, the \(90^\circ\) rightangled tee junction is considered here for the purpose of validation.
Problem definition
The following figure shows the schematic of the teejunction. Here, L = 3.0 m, W = 1.0 m respectively, the Reynolds number based on the width is 300, and V is the velocity in the ydirection. For simplicity, we have assumed the velocity, V = 1 m/s. Using these values the resulting kinematic viscosity was 0.00333 \(m^2/s\).
The summary of the flow properties and geometric details are given in the following table.
\(Re_w\) 
\(L\) 
\(W\) 
\(V~(m/s)\) 
\(p~(m^2/s^2)\) 
\(\nu~(m^2/s)\) 

300 
3.0 
1.0 
1.0 
\((0)\) Gauge 
0.00333 
As we have assumed the flow incompressible, the density (\(\rho\)) remains constant. In addition, since the fluid temperature is not considered, the viscosity remains constant. For incompressible flows, the kinematic forms of pressure and viscosity are always used in Caelus.
Computational Domain and Boundary Conditions
Since this is an internal flow problem, the computational domain is contained within teejunction geometry. The details are shown in Figure 58. As indicated, all teejunction walls have a noslip boundary condition which has been highlighted in blue. At the inlet, a fully developed laminar flow parabolic profile is applied, otherwise a much longer main branch would be required for the flow to develop. The domain has two outlets, one at the end of the main channel and the other at the end of side branch. Note the exit pressures at the two outlets are equal.
Boundary Conditions and Initialisation
The following are the boundary condition details used for the computational domain:
 Inlet
Velocity: Parabolic velocity profile; centerline velocity of \(v = 1.0~m/s\) in \(y\) direction
Pressure: Zero gradient
 Noslip wall
Velocity: Fixed uniform velocity \(u, v, w = 0\)
Pressure: Zero gradient
 Outlet1
Velocity: Zero gradient velocity
Pressure: Fixed uniform gauge pressure \(p = 0\)
 Outlet2
Velocity: Zero gradient velocity
Pressure: Fixed uniform gauge pressure \(p = 0\)
 Initialisation
Velocity: Fixed uniform velocity \(u, v, w = 0\)
Pressure: Zero Gauge pressure
Computational Grid
The 2D structured grid is shown in Figure 59. Since Caelus is a 3D computational framework, it necessitates the grid to also be 3D. Therefore, a 3D grid was obtained by extruding the 2D grid in the positive \(z\) direction by one cell. The final 3D grid was then exported to the Caelus format (polyMesh). The two \(xy\) planes obtained as a result of grid extrusion need boundary conditions to be specified. As the flow over a flatplate is generally 2D, we do not need to solve the flow in the third dimension. This is achieved in Caelus by specifying empty boundary condition for each plane. Although, no flow is computed in the \(z\) direction, a velocity of \(w = 0\) has to be specified for the velocity boundary condition as indicated above.
A total of 2025 cells comprise the teejunction of which, 90 cells are distributed along the height of the main channel, and 45 along the length of the side branch. The distribution is such that a dimensional length of \(L = 1~m\) has a total of 45 cells, giving a distribution of \((2/3)*45 = 30\) cells for the \((2/3) L\) segment of the main channel. The width, \(W\), consists of 15 cells.
Results and Discussion
A timedependent solution to the twodimensional flatplate was obtained using Caelus 9.04. The SLIM transient solver was used here and the flow was simulated sufficiently long such that steady separated flow was established. To ensure this, shearstress distribution was monitored on the lower wall of the side branch. The simulation was stopped once the separation and reattachment locations no longer varied with time.
Mass flow was calculated at the inlet and at the main outlet (outlet1) and the massflow ratio was subsequently calculated. The below table compares the SLIM result with the experimental value. As can be noted, the agreement between the two is excellent.
Experimental 
SLIM 
Percentage Difference 


Flow Split 
\(0.887\) 
\(0.886\) 
\(0.112~\%\) 
Conclusions
The steady, incompressible, twodimensional laminar flow in a rightangled \(90^\circ\) teejunction was simulated using Caelus 9.04 with the SLIM solver and validated against experimental data resulting in excellent agreement.
Circular CylinderÂ¶
Laminar, incompressible flow over a circular cylinder
Nomenclature
Symbol 
Definition 
Units (SI) 

\(D\) 
Diameter of the cylinder 
\(m\) 
\(f\) 
Frequency 
\(hz\) 
\(p\) 
Kinematic pressure 
\(Pa/\rho~(m^2/s^2)\) 
\(Re_D\) 
Reynolds number based on diameter 
Nondimensional 
\(St\) 
Strouhal number 
Nondimensional 
\(U\) 
Freestream velocity in xdirection 
\(m/s\) 
\(x\) 
Distance in xdirection 
\(m\) 
\(y\) 
Distance in ydirection 
\(m\) 
\(\nu\) 
Kinematic viscosity 
\(m^2/s\) 
In this validation study, laminar incompressible flow over a 2D circular cylinder is investigated at a Reynolds number of 100. This classical configuration represents flow over a bluff body dominated by a wake region. For flows having a low Reynolds number (\(40 \leq Re_D \leq 150\)), periodic vortex shedding occurs in the wake. The phenomenon of vortex shedding behind bluff bodies is referred to as the Karman Vortex Street [2] and provides an transient case for CFD code validation.
In his work, Roshko [2] experimentally studied wake development behind twodimensional circular cylinders from Reynolds number ranging from 40 to 10000. For Reynolds numbers of 40 to 150, the so called the stable range [Roshko1954], regular vortex streets are formed with no evidence of turbulence motion in the wake. Therefore, at a Reynolds number of 100, the vortex shedding exhibits smooth, coherent structures making it ideally suited for validating laminar CFD calculations. The frequency associated with the oscillations of the vortex streets can be characterized by the Strouhal Number (\(St\)). The Strouhal Number is a nondimensional number defined as
where, \(f\) is the frequency of oscillations of vortex shedding, \(D\) is the diameter of the cylinder and \(U\) is the freestream velocity of the flow. Experimentally [Roshko1954], it has been determined that for a Reynolds number based on the diameter of the cylinder of 100, the Strouhal number \(St \approx 0.16  0.17\). One of the main objectives in this study was to compare the \(St\) for the CFD calculation to the experimental data of Roshko [Roshko1954]. Provided the cylinder has a sufficient span length, the flow characteristics can be assumed to be twodimensional as the experiments suggest.
Problem Definition
Figure 60 shows the schematic of the twodimensional circular cylinder. Here, the diameter D = 2 m and is the characteristic length for the Reynolds number, which is 100. For simplicity, the freestream velocity was taken to be U = 1 m/s in the x direction. Using these values the kinematic viscosity was calculated to be 0.02 \(m^2/s\).
In the table below, a summary of the freestream conditions are provided
\(Re_D\) 
\(D\) 
\(U~(m/s)\) 
\(p~(m^2/s^2)\) 
\(\nu~(m^2/s)\) 

100 
2.0 
1.0 
\((0)\) Gauge 
0.02 
Here, the flow is assumed to be incompressible and therefore the density is constant. Further, no temperature is evaluated in this calculation and hence the viscosity also remains constant. For incompressible flows, the kinematic forms of pressure and viscosity are always used in Caelus.
Computational Domain and Boundary Conditions
A rectangular computational domain in the \(xy\) plane was constructed surrounding the circular cylinder as shown in Figure 61. A full cylinder must be used due to the oscillatory nature of the shed vortices. Rhe domain extends by 5 diameters of cylinder and 20 diameters downstream. Since the flow here is viscous dominated, sufficient downstream length is required to capture the vortex separation from the surface of the cylinder and the subsequently shedding in the wake. In the \(y\) direction, the domain extends 5 diameters on either side. From the figure, multiple inlet boundaries to this domain can be seen, one at the upstream boundary and the other two for the top and bottom boundaries. This type of configuration is needed to appropriately model the inflow, similar to an undisturbed flow in an experimental setup. It is noted that for top and bottom boundaries, the flow is in the \(x\) direction. The outlet is located at the downstream boundary. The cylindrical wall is a noslip boundary condition.
Boundary Conditions and Initialisation
Following are the details of the boundary conditions used:
 Inlet1
Velocity: Fixed uniform velocity \(u = 1.0~m/s\) in \(x\) direction
Pressure: Zero gradient
 Inlet2
Velocity: Fixed uniform velocity \(u = 1.0~m/s\) in \(x\) direction
Pressure; Zero gradient
 Noslip wall
Velocity: Fixed uniform velocity \(u, v, w = 0\)
Pressure: Zero gradient
 Outlet
Velocity: Zero gradient velocity
Pressure: Fixed uniform gauge pressure \(p = 0\)
 Initialisation
Velocity: Fixed uniform velocity \(u = 1.0~m/s\) in \(x\) direction
Pressure: Zero Gauge pressure
Computational Grid
The computational grid in 2D was generated using Pointwise in the \(xy\) plane. Since Caelus is a 3D computational framework, it necessitates the grid to also be 3D. Therefore, a 3D grid was obtained using Pointwise by extruding the 2D grid in the positive \(z\) direction by one cell. The final 3D grid was then exported to the Caelus format (polyMesh). The two \(xy\) planes obtained as a result of grid extrusion need boundary conditions to be specified. As the flow over a flatplate is generally 2D, we do not need to solve the flow in the third dimension. This is achieved in Caelus by specifying empty boundary condition for each plane. Although, no flow is computed in the \(z\) direction, a velocity of \(w = 0\) has to be specified for the velocity boundary condition as indicated above.
The 2D domain consisted of 9260 cells. An Ogrid topology was constructed around the cylinder (see the right figure) with 10 cells in the radial direction and 84 cells in the circumferential direction. 31 cells were used upstream of the Ogrid, in the \(x\) direction while 100 cells were used downstream. The region of interest is about 10 diameters downstream, where the grids are refined. In the \(y\) direction, 21 cells were used above and below the Ogrid region.
Results and Discussion
A timedependent simulation was carried out using the Caelus 9.04 with the SLIM solver. To capture the transient startup process, the calculation was started from time t = 0 s and was simulated up to t = 360 s, while lift and drag forces over the cylindrical surface were monitored at a frequency of 2 Hz. It was found that the onset of vortex shedding occurred after about t = 90 s which was then followed by a steady shedding process. A Fast Fourier transformation (FFT) was carried out on the lift force data and the peak frequency of vortex shedding occurred at \(f = 0.0888\) Hz. Based on this value, it takes about 7.8 cycles for the shedding to start.
Using the peak frequency value of \(f = 0.0888\) Hz, \(St\) was evaluated. The table below compares the computed value from SLIM with that of the experiment. The agreement is good given that experimental uncertainty can be relatively high at low Reynolds numbers.
Frequency (Hz) 
Strouhal Number 


Experimental 
0.0835 
\(0.167\) 
SLIM 
0.0888 
\(0.177\) 
Conclusions
The transiet, incompressible, twodimensional flow over a circular cylinder was simulated using Caelus 9.04 to estimate the peak frequency of vortex shedding. The value was compared to well known experimental data resulting in good agreement.
Triangular CavityÂ¶
Laminar, incompressible flow inside a lid driven Triangular Cavity
Nomenclature
Symbol 
Definition 
Units (SI) 

\(D\) 
Depth of the cavity 
\(m\) 
\(p\) 
Kinematic pressure 
\(Pa/\rho~(m^2/s^2)\) 
\(Re_D\) 
Reynolds number based on depth 
Nondimensional 
\(U\) 
Freestream velocity in xdirection 
\(m/s\) 
\(W\) 
Width 
\(m\) 
\(x\) 
Distance in xdirection 
\(m\) 
\(y\) 
Distance in ydirection 
\(m\) 
\(\nu\) 
Kinematic viscosity 
\(m^2/s\) 
This validation study concerns the laminar, incompressible flow inside a lid driven triangular cavity. Here, the top wall of the cavity moves at a constant velocity initiating a recirculating motion within the cavity. Due to the viscous nature of the flow, a boundary layer develops in the direction of the moving lid, while flow reversal occurs due to the recirculating flow. The flow feature of interest is the velocity distribution along the centreline of the cavity.
Benchmark experiments on this configuration has been reported in Jyotsna and Venka [11] for a Reynolds number of 800. The main objective of this validation case was to compare the \(x\) velocity distribution against experimental data.
Problem Definition
A schematic of the triangular cavity is presented in Figure 63 where depth of the cavity D = 4 m and the width W = 2 m. The Reynolds number based on the cavity depth is 800 and the wall velocity is U = 2 m/s. Using the Reynolds number, U, and D, kinematic viscosity was calculated to be 0.01 \(m^2/s\).
The table below summaries the flow properties
\(Re_D\) 
\(D~(m)\) 
\(W~(m)\) 
\(U~(m/s)\) 
\(p~(m^2/s^2)\) 
\(\nu~(m^2/s)\) 

800 
4.0 
2.0 
2.0 
\((0)\) Gauge 
\(0.01\) 
The flow in this case is assumed to be incompressible and hence the density remained constant throughout. Further, the temperature field is not accounted into the calculation, and therefore the viscosity of the flow can also remain constant. Since viscosity is constant, it becomes more convenient to specify it as kinematic viscosity. It should be noted that in Caelus for incompressible flows, both pressure and viscosity are always specified as kinematic.
Computational Domain and Boundary Conditions
The computational domain is the triangular cavity shown in Figure 64. Highlighted in blue, the side walls of the cavity have a noslip boundary condition while the top wall, highlighted in green, has a uniform velocity in the \(x\) direction.
Boundary Conditions and Initialisation
 Moving wall
Velocity: Fixed uniform velocity \(u = 2.0~m/s\) in \(x\) direction
Pressure: Zero gradient
 Noslip wall
Velocity: Fixed uniform velocity \(u, v, w = 0\)
Pressure: Zero gradient
 Initialisation
Velocity: Fixed uniform velocity \(u, v, w = 0\)
Pressure: Zero Gauge pressure
Computational Grid
The 2D grid in \(xy\) plane is shown in Figure 65. A hybrid grid is employed for this case with a total of 5538 cells. Up to a depth of D = 1.35 m a structured grid is used while below that value an unstructured triangular grid is used. An unstructured grid is used in the bottom portion because it resulted in lower skewness in this vicinity. For the structured region, 40 cells are distributed across the width of the cavity and 40 along the depth. The cavity walls in the unstructured region have 100 cells along each. The interface of the two regions is node matched and has 40 cells across the width. The grid close to the cavity lid was refined to better capture the shear layer.
The flow characteristics in the cavity can be assumed to be two dimensional and here it has been solved with the same assumption. Since Caelus is a 3D computational framework, it necessitates the grid to also be 3D. Therefore, a 3D grid was obtained using Pointwise by extruding the 2D grid in the positive \(z\) direction by one cell. The final 3D grid was then exported to the Caelus format (polyMesh). The two \(xy\) planes obtained as a result of grid extrusion need boundary conditions to be specified. As the flow over a flatplate is generally 2D, we do not need to solve the flow in the third dimension. This is achieved in Caelus by specifying empty boundary condition for each plane. Although, no flow is computed in the \(z\) direction, a velocity of \(w = 0\) has to be specified for the velocity boundary condition as indicated above.
Results and Discussion
A steady solution to the cavity was obtained using Caelus 9.04 with the SLIM solver. While a timedependent approach was used, the solution was simulated sufficiently long so that steady flow was achieved. To determine when this occured the velocity distribution along the cavity centreline was monitored with respect to time. The simulations was stopped when no appreciable changes were observed.
In Figure 66, the \(x\) velocity distribution along the cavity centreline is compared with that of the benchmark experimental data. The \(y\) distance is normalised with the cavity depth (\(D\)) which gives \(y/d = 0\) at the cavity lid and \(y/d = 1\) at the bottom vertex. Similarly, the \(u\) velocity is normalised with the velocity of the cavity lid (\(u_L\)).
As seen in Figure 66 above, the comparison the experiment is excellent.
Conclusions
A steady, incompressible flow past a twodimensional triangular cavity was simulated using Caelus 9.04 on a hybrid grid. The velocity distribution along the centreline of the cavity was compared with the benchmark experimental data, it was found that the SLIM results compared very favorably.
Spherical CavityÂ¶
Natural Convection in a Spherical Cavity
Nomenclature
Symbol 
Definition 
Units (SI) 

\(C_p\) 
Specific heat at constant pressure 
\(J/kg \cdot K\) 
\(D\) 
Diameter of the sphere 
\(m\) 
\(g\) 
Gravitational constant 
\(9.80~m/s^2\) 
\(k\) 
Thermal conductivity 
\(W/m \cdot K\) 
\(Nu\) 
Nusselt number 
Nondimensional 
\(p\) 
Kinematic pressure 
\(Pa/\rho~(m^2/s^2)\) 
\(Pr\) 
Prandtl number 
Nondimensional 
\(r\) 
Radius of the sphere 
\(m\) 
\(Ra\) 
Rayleigh number 
Nondimensional 
\(T\) 
Temperature 
\(K\) 
\(u\) 
velocity in xdirection 
\(m/s\) 
\(x\) 
Distance in xdirection 
\(m\) 
\(y\) 
Distance in ydirection 
\(m\) 
\(z\) 
Distance in zdirection 
\(m\) 
\(\beta\) 
Coefficient of thermal expansion 
\(1/K\) 
\(\nu\) 
Kinematic viscosity 
\(m^2/s\) 
\(\rho\) 
Density 
\(kg/m^3\) 
In this study, laminar flow inside a heated spherical cavity is investigated. A temperature gradient was applied to the surface of the cavity inducing fluid motion due to buoyancy effects. Here, a Rayleigh number (\(Ra\)) of 2000 and a Prandtl number (\(Pr\)) of 0.7 is used to set appropriate thermal conditions to the problem.
Analytical solutions to the heated spherical cavity were investigated by McBain and Stephens [4] at various Rayleigh numbers up to \(Ra = 10000\). However, it was noted that asymptotic heat transfer rate predictions as a function of Nusselt number (\(Nu\)) deviated largely with increase in Rayleigh number. Therefore a value of \(Ra = 2000\) was chosen to avoid these deviations. At this value, the analytical heat transfer rate compares very well with the asymptotic solution. Isotherms from the analytical solution were used for comparision with numerical results.
Problem Definition
The schematic of the spherical cavity is depicted Figure 67 and as can be seen, only half of the sphere is considered here with an \(xy\) symmetry plane. The sphere is located at \(x = 0, y = 0, z = 0\) with a radius \(r = 0.5~m\).
As indicated, a nonuniform temperature profile was used as a thermal boundary condition on the spherical wall. The temperature (\(T\)) was specified as a function of \(x\):
with,
The energy equation in a nondimensional form [4] for a steady state can be expressed as:
where, \(Gr\) is the Grashof number. Since, \(Ra = Gr~Pr\) the above equation can be rewritten along with expanding the Laplacian term (\(\nabla^2\)) as
where \(\alpha\) is the thermometric conductivity. It is reasonable to assume that (\(1/Ra\)) is approximately equal to the thermometric conductivity (\(\alpha\)), which is given as
where, \(k\), \(\rho\), \(C_p\) and \(\nu\) are the thermal conductivity, density, specific heat capacity and kinematic viscosity of the fluid respectively. Using the above relation with a value of \(Ra = 2000\) and \(Pr = 0.7\), the kinematic viscosity was calculated to be \(\nu = 3.4 \times 10^{4}~m^2/s\). The coefficient of thermal expansion (\(\beta\)) needed to model the Boussinesq buoyancy term was evaluated from the following relation
where, \(g\) is the acceleration due to gravity and \(D\) is the diameter of the sphere. \(\beta\) was calculated to \(3.567\times10^{5}~1/K\). In the following table, a summary of the properties are given. Note that gravity acts in \(y\) direction.
\(Ra\) 
\(Pr\) 
\(T~(K)\) 
\(p~(m^2/s^2)\) 
\(\nu~(m^2/s)\) 
\(\beta~(1/K)\) 

\(2000\) 
\(0.7\) 
\(T = x\) 
\((0)\) Gauge 
\(3.4 \times 10^{4}\) 
\(3.567\times10^{5}\) 
Although the temperature is calculated in this simulation, a constant viscosity is used. Since the temperature gradient is very small (\(\mathcal{O}(1)\)), effect of temperature on the viscosity would be insignificant. The kinematic definition of pressure is used here.
Computational Domain and Boundary Conditions
The computational domain was a half sphere with an \(xy\) plane of symmetry at \(z = 0~m\). The surface temperature was prescribed as discussed above. The initialisation of the fluid temperature within the sphere follows that of the surface temperature (\(T = x\)) and is depicted in Figure 68 at the symmetry plane. Note that this figure also aids in providing a clarity of understanding for the temperature variation over the spherical surface.
Boundary Conditions and Initialisation
 Wall
Velocity: Fixed uniform velocity \(u, v, w = 0\)
Pressure: Uniform zero Buoyant Pressure
Temperature: Linear unction of \(x\) (\(T = x\))
 Symmetry Plane
Velocity: Symmetry
Pressure: Symmetry
Temperature: Symmetry
 Initialisation
Velocity: Fixed uniform velocity \(u, v, w = 0\)
Pressure: Uniform zero Buoyant Pressure
Temperature: Linear function of \(x\) (\(T = x\))
Computational Grid
The computational grid for the half sphere was generated using Pointwise. A fully structured grid was constructured with a total of 18564 cells. As seen in Figure 69, an OH topology used where an Hblock is centred within 5 Oblocks.
Results and Discussion
The steady solution to the natural convection in a buoyant sphere was obtained using Caelus 9.04 with the SLIM solver that includes the a buoyancy source term based on the Boussinesq assumption. Since SLIM is inherently timeaccurate, the simulation was run sufficiently long such that a steady state was achieved. In Figure 70, compares the isotherms obtained with SLIM and the analytical isotherms obtained with a first order approximation. Close agreement was observed.
Conclusions
A validation study of a buoyant flow inside a spherically heated cavity was conducted using Caelus 9.04. The isotherms obtained from the CFD results were compared with the first order analytical solution and excellent agreement was observed.
Validation: Incompressible Turbulent FlowsÂ¶
This section discusses the turbulence model verification and validation conducted for Caelus 9.04. The implemented Spalartâ€“Allmaras, Spalartâ€“Allmaras Curvature Correction, \(k\omega~\rm{SST}\) turbulence models are verified with NASAâ€™s CFL3D results. To validate these models, the Caelus results are compared with available experimental data. The cases considered for this exercise are obtained from the Turbulence Modeling Resource database.
Zero Pressure Gradient Flat PlateÂ¶
Turbulent, incompressible flow over a twodimensional SharpLeading Edge flatplate
Nomenclature
Symbol 
Definition 
Units (SI) 

\(a\) 
Speed of sound 
\(m/s\) 
\(c_f\) 
Skin friction coefficient 
Nondimensional 
\(k\) 
Turbulent kinetic energy 
\(m^2/s^2\) 
\(L\) 
Length of the plate 
\(m\) 
\(M_\infty\) 
Freestream Mach number 
Nondimensional 
\(p\) 
Kinematic pressure 
\(Pa/\rho~(m^2/s^2)\) 
\(Re_L\) 
Reynolds number 
\(1/m\) 
\(T\) 
Temperature 
\(K\) 
\(u\) 
Local velocity in xdirection 
\(m/s\) 
\(U\) 
Freestream velocity in xdirection 
\(m/s\) 
\(x\) 
Distance in xdirection 
\(m\) 
\(y\) 
Distance in ydirection 
\(m\) 
\(y^+\) 
Wall distance 
Nondimensional 
\(\mu\) 
Dynamic viscosity 
\(kg/m~s\) 
\(\nu\) 
Kinematic viscosity 
\(m^2/s\) 
\(\tilde{\nu}\) 
Turbulence field variable 
\(m^2/s\) 
\(\rho\) 
Density 
\(kg/m^3\) 
\(\omega\) 
Specific dissipation rate 
\(1/s\) 
\(\infty\) 
Freestream conditions 

\(t\) (subscript) 
Turbulent property 
Introduction
In this case, steady turbulent incompressible flow over a twodimensional sharpleading edge flatplate is considered at zero angle of incidence, which generates a turbulent boundary layer with zeropressure gradient over the flat plate. The simulations are carried over a series of four successively refined grids and the solutions are compared with the CFL3D data. This therefore serves as both grid independence study and verification of the turbulence models. The distribution of skinfriction coefficient (\(c_f\)) along the plate is used to verify the accuracy of the models.
Problem definition
This exercise is based on the Turbulence Modeling Resource case for a flatplate and follows the same conditions used in the incompressible limit. However, note that CFL3D uses a freestream Mach number (\(M_\infty\)) of 0.2 as it is a compressible solver. The schematic of the geometric configuration is shown in Figure 71.
The length of the plate is L = 2.0 m, wherein, x = 0 is at the leading edge and the Reynolds number is \(5 \times 10^6\) and \(U_\infty\) is the freestream velocity in the x direction. The inflow temperature in this case is assumed to be 300 K and for Air as a perfect gas, the speed of sound (\(a\)) can be evaluated to 347.180 m/s. Based on the freestream Mach number and speed of sound, velocity can be calculated, which is \(U_\infty\) = 69.436 m/s. Using velocity and Reynolds number, the kinematic viscosity can be evaluated. The following table summarises the freestream conditions used:
Fluid 
\(Re_L~(1/m)\) 
\(U_\infty~(m/s)\) 
\(p_\infty~(m^2/s^2)\) 
\(T_\infty~(K)\) 
\(\nu_\infty~(m^2/s)\) 

Air 
\(5 \times 10^6\) 
69.436113 
\((0)\) Gauge 
300 
\(1.38872\times10^{5}\) 
Note that in an incompressible flow the density (\(\rho\)) does not change and is reasonable to assume a constant density throughout the calculation. In addition, temperature is not considered and therefore the viscosity can also be held constant. In Caelus, for incompressible flows, pressure and viscosity are always specified as kinematic.
Turbulent Properties for Spalartâ€“Allmaras model
The turbulent inflow boundary conditions used for the Spalartâ€“Allmaras model were calculated as \(\tilde{\nu}_{\infty} = 3 \cdot \nu_\infty\) and subsequently turbulent eddy viscosity was evaluated. The following table provides the values of these used in the current simulations:
\(\tilde{\nu}_\infty~(m^2/s)\) 
\(\nu_{t~\infty}~(m^2/s)\) 

\(4.166166 \times 10^{5}\) 
\(2.9224023 \times 10^{6}\) 
Turbulent Properties for komega SST model
The turbulent inflow boundary conditions used for \(k\omega~\rm{SST}\) were calculated as follows and is as given in Turbulence Modeling Resource
Note that the dynamic viscosity for the above equation is obtained from Sutherland formulation and density is calculated as \(\rho = \mu / \nu\). The below table provides the turbulent properties used in the current simulations
\(k_{\infty}~(m^2/s^2)\) 
\(\omega_{\infty}~(1/s)\) 
\(\nu_{t~\infty}~(m^2/s)\) 

\(1.0848 \times 10^{3}\) 
\(8679.5135\) 
\(1.24985 \times 10^{7}\) 
Computational Domain and Boundary Conditions
The computational domain is a rectangular block encompassing the flatplate. Figure 72 below shows the details of the boundaries used in twodimensions (\(xy\) plane). As can be seen, the region of interest (highlighted in blue) extends between \(0\leq x \leq 2.0~m\) and has a noslip boundary condition. Upstream of the leading edge, a symmetry boundary is used to simulate a freestream flow approaching the flatplate. The inlet boundary as shown in Figure 72 is placed at start of the symmetry at \(x = 0.3333~m\) and the outlet at the exit plane of the noslip wall (blue region) at \(x = 2.0~m\). A symmetry plane condition is used for the entire top boundary.
Boundary Conditions and Initialisation
Following are the boundary condition details used for the computational domain:
 Inlet
Velocity: Fixed uniform velocity \(u = 69.436113~m/s\) in \(x\) direction
Pressure: Zero gradient
Turbulence:
SpalartAllmaras (Fixed uniform values of \(\nu_{t~\infty}\) and \(\tilde{\nu}_{\infty}\) as given in the above table)
\(k\omega~\textrm{SST}\) (Fixed uniform values of \(k_{\infty}\), \(\omega_{\infty}\) and \(\nu_{t~\infty}\) as given in the above table)
 Symmetry
Velocity: Symmetry
Pressure: Symmetry
Turbulence: Symmetry
 Noslip wall
Velocity: Fixed uniform velocity \(u, v, w = 0\)
Pressure: Zero gradient
Turbulence:
SpalartAllmaras (Fixed uniform values of \(\nu_{t}=0\) and \(\tilde{\nu}=0\))
\(k\omega~\textrm{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:
SpalartAllmaras (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~m/s\) in \(x\) direction
Pressure: Zero Gauge pressure
Turbulence:
SpalartAllmaras (Fixed uniform values of \(\nu_{t~\infty}\) and \(\tilde{\nu}_{\infty}\) as given in the above table)
\(k\omega~\textrm{SST}\) (Fixed uniform values of \(k_{\infty}\), \(\omega_{\infty}\) and \(\nu_{t~\infty}\) as given in the above table)
Computational Grid
The 3D structured grid was obtained from Turbulence Modeling Resource as a Plot3D and was converted to Caelus format using Pointwise. It should be noted that the flow normal direction in the Plot3D grids is \(z\) and the twodimensional plane of interest is in \(xz\) directions. Since the flowfield of interest is twodimensional, and simpleSolver being a 3D solver, the two \(xz\) planes are specified with empty boundary conditions. As mentioned earlier, a series of four grids were considered from the original set of five, excluding the coarsest. Details of the different grids used are given below.
Grid 
Cells in \(x\)direction 
Cells in \(z\)direction 
Total 
\(y^+\) 

Grid2 
68 
48 
3264 
0.405 
Grid3 
136 
96 
13,056 
0.203 
Grid4 
272 
192 
52,224 
0.101 
Grid5 
544 
384 
208,896 
0.05 
In Figure 73, the 2D grid in the \(xz\) plane is shown for Grid4. As can be seen, the grid is refined close to the wall in order to capture the turbulent boundary layer accurately. All grids have \(y^+ < 1\) and no wall function is used for the wall boundary in the current verification cases.
Results and Discussion
The steadystate solution of the turbulent flow over a flat plate was obtained using Caelus 9.04. The simpleSolver was used here and the solution was simulated for a sufficient length until the residuals for pressure, velocity and turbulent quantities were less than \(1 \times 10^{6}\). The finite volume discretization of the gradient of pressure and velocity was carried out using the linear approach. Where as the divergence of velocity and mass flux was carried out through the linear upwind method. However, for the divergence of the turbulent quantities, upwind approach was utilised and linear approach for the divergence of the Reynolds stress terms. For the discretization of the Laplacian terms, again linear corrected method was used. For some grids having greater than 50 degree nonorthogonal angle, linear limited with a value of 0.5 was used for the Laplacian of the turbulent stress terms.
Spalartâ€“Allmaras
In Figure 74, the skinfriction distribution along the flatplate obtained from Caelus for different grids is shown. As can be seen, all grids produce the same skinfriction values suggesting a gridindependent solution is achieved.
In Figure 75, the skinfriction distribution obtained from Caelus on Grid5 is compared with CFL3D of the same grid. An excellent agreement is obtained all along the plate.
kOmega SST
The skinfriction distribution for various grids obtained from \(k\omega~\rm{SST}\) model is shown in Figure 76.
The skinfriction comparison between Caelus and CFL3D for Grid5 is shown in Figure 77.
Conclusions
The steady turbulent flow over a twodimensional flatplate was simulated using Caelus 9.04 utilising the simpleSolver. The simulations were carried out with two turbulence models and the obtained solutions were verified against CFL3D data. The results were found to be in good agreement with CFL3D and suggesting the turbulence implementation in Caelus is accurate.
Twodimensional Bump in a ChannelÂ¶
Turbulent, incompressible flow over a twodimensional bump in a channel
Symbol 
Definition 
Units (SI) 

\(a\) 
Speed of sound 
\(m/s\) 
\(c_f\) 
Skin friction coefficient 
Nondimensional 
\(k\) 
Turbulent kinetic energy 
\(m^2/s^2\) 
\(L\) 
Length of the bump 
\(m\) 
\(M_\infty\) 
Freestream Mach number 
Nondimensional 
\(p\) 
Kinematic pressure 
\(Pa/\rho~(m^2/s^2)\) 
\(Re_L\) 
Reynolds number 
\(1/m\) 
\(T\) 
Temperature 
\(K\) 
\(u\) 
Local velocity in xdirection 
\(m/s\) 
\(U\) 
Freestream velocity in xdirection 
\(m/s\) 
\(x\) 
Distance in xdirection 
\(m\) 
\(z\) 
Distance in zdirection 
\(m\) 
\(y^+\) 
Wall distance 
Nondimensional 
\(\mu\) 
Dynamic viscosity 
\(kg/m~s\) 
\(\nu\) 
Kinematic viscosity 
\(m^2/s\) 
\(\tilde{\nu}\) 
Turbulence field variable 
\(m^2/s\) 
\(\rho\) 
Density 
\(kg/m^3\) 
\(\omega\) 
Specific dissipation rate 
\(1/s\) 
\(\infty\) 
Freestream conditions 

\(t\) (subscript) 
Turbulent property 
Introduction
This case covers the verification of turbulent incompressible flow over a twodimensional bump in a channel. The bump acts as a perturbation causing local changes to the velocity and pressure over the bump surface. Since the perturbation is quite small, the flow remains attached to the bump surface. The simulations are carried out over a series of four grids which are successively refined and are studied over two turbulence models, serving as a grid independence study. The distribution of skinfriction (\(c_f\)) is then compared with that of the CFL3D data.
Problem definition
This verification exercise is based on the Turbulence Modeling Resource case for a 2D bump in a channel and follows the same conditions used in the limit of incompressibility. In this case, CFL3D uses a freestream Mach number (\(M_\infty\)) of 0.2. The schematic of the geometric configuration is shown in Figure 78.
The location of the plate upstream of the bump begins at x=0 m and ends at x = 1.5 m, giving a total plate length of L = 1.5 m. The flow has a Reynolds number of \(3 \times 10^6\) with a freestream velocity \(U_\infty\) in the xdirection. The temperature of the inflow is assumed to be 300 K and for Air as a perfect gas, the speed of sound (\(a\)) can be evaluated to 347.180 m/s. With the freestream Mach number and speed of sound, the freestream velocity was calculated to \(U_\infty\) = 69.436 m/s. Kinematic viscosity was then obtained from velocity and Reynolds number. The following table summarises the freestream conditions used for this case.
Fluid 
\(Re_L~(1/m)\) 
\(U_\infty~(m/s)\) 
\(p_\infty~(m^2/s^2)\) 
\(T_\infty~(K)\) 
\(\nu_\infty~(m^2/s)\) 

Air 
\(3 \times 10^6\) 
69.436113 
\((0)\) Gauge 
300 
\(2.314537\times10^{5}\) 
Since the flow is incompressible, the density (\(\rho\)) does not change and therefore, a constant density is assumed throughout the calculation. Further, the temperature is not considered and hence it does not have any influence on viscosity, which is therefore kept constant. Note that in Caelus, pressure and viscosity are always specified as kinematic for a incompressible flow simulation.
Turbulent Properties for Spalartâ€“Allmaras model
The turbulent inflow boundary conditions used for the Spalartâ€“Allmaras model were calculated as \(\tilde{\nu}_{\infty} = 3 \cdot \nu_\infty\) and subsequently turbulent eddy viscosity was evaluated. The following table provides the values of these used in the current simulations:
\(\tilde{\nu}_\infty~(m^2/s)\) 
\(\nu_{t~\infty}~(m^2/s)\) 

\(6.943611 \times 10^{5}\) 
\(4.8706713 \times 10^{6}\) 
Turbulent Properties for komega SST model
The turbulent inflow boundary conditions used for \(k\omega~\rm{SST}\) were calculated as follows and is as given in Turbulence Modeling Resource
Note that the dynamic viscosity for the above equation is obtained from Sutherland formulation and density is calculated as \(\rho = \mu / \nu\). The below table provides the turbulent properties used in the current simulations:
\(k_{\infty}~(m^2/s^2)\) 
\(\omega_{\infty}~(1/s)\) 
\(\nu_{t~\infty}~(m^2/s)\) 

\(1.0848 \times 10^{3}\) 
\(5207.6475\) 
\(2.08310 \times 10^{7}\) 
Computational Domain and Boundary Conditions
The computational domain consists of a rectangular channel encompassing the bump. In Figure 79 , the details of the boundaries used in twodimensions (\(xy\) plane) are shown. The region of interest, which is the bump extends between \(0\leq x \leq 1.5~m\) and has a noslip boundary condition. Upstream and downstream of the bump, the symmetry boundary extends about 17 bump lengths. The inlet boundary is placed at the start of the symmetry at \(x = 25.0~m\) and the outlet is placed at \(x = 26.5~m\). For the entire top boundary, symmetry plane condition is used.
Boundary Conditions and Initialisation
Following are the boundary condition details used for the computational domain:
 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 the above table)
\(k\omega~\rm{SST}\) (Fixed uniform values of \(k_{\infty}\), \(\omega_{\infty}\) and \(\nu_{t~\infty}\) as given in the above table)
 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~m/s\) in \(x\) direction
Pressure: Zero Gauge pressure
Turbulence:
Spalartâ€“Allmaras (Fixed uniform values of \(\nu_{t~\infty}\) and \(\tilde{\nu}_{\infty}\) as given in the above table)
\(k\omega~\rm{SST}\) (Fixed uniform values of \(k_{\infty}\), \(\omega_{\infty}\) and \(\nu_{t~\infty}\) as given in the above table)
Computational Grid
The 3D computational grid was obtained from Turbulence Modeling Resource as a Plot3D and was converted to Caelus format using Pointwise. In the Plot3D computational grid, the flow normal direction is \(z\) and thus the twodimensional plane of interest is in \(xz\) directions. Further, since the flowfield is of twodimensional, and the simpleSolver being a 3D solver, the two \(xz\) planes are specified with empty boundary conditions. A series of four grids were considered from the original set of five, excluding the coarsest grid and the following table give its details.
Grid 
Cells in \(x\)direction 
Cells in \(z\)direction 
Total 
\(y^+\) 

Grid2 
176 
80 
14,080 
0.236 
Grid3 
352 
160 
56,320 
0.118 
Grid4 
704 
320 
225,280 
0.059 
Grid5 
1408 
640 
901,120 
0.03 
The 2D grid in \(xz\) plane is shown in Figure 80 for Grid3. As can be noted, the grid is sufficiently refined close to the wall in the normal direction. In addition, the grids are refined in the vicinity of the bump, including both upstream and downstream which can be seen in the inset. All grids have a \(y^+ < 1\) and no wall function is used for the wall boundary in the current verification cases.
Results and Discussion
The steadystate solution of the turbulent flow over a twodimensional bump was obtained using Caelus 9.04. The simpleSolver was used for the calculations and was run for a sufficient length until the residuals for pressure, velocity and turbulent quantities were less than \(1 \times 10^{6}\). The finite volume discretization of the gradient of pressure and velocity was carried out using the linear approach. Where as the divergence of velocity and mass flux was carried out through the linear upwind method. However, for the divergence of the turbulent quantities, upwind approach was utilised and linear approach for the divergence of the Reynolds stress terms. For the discretization of the Laplacian terms, again linear corrected method was used. For some grids having greater than 50 degree nonorthogonal angle, linear limited with a value of 0.5 was used for the Laplacian of the turbulent stress terms.
Spalartâ€“Allmaras
The skinfriction distributions over the 2D bump obtained from Caelus for different grids are shown in Figure 81. There is very little difference in the skinfriction beyond Grid2 suggesting that a gridindependence solution is achieved.
In Figure 82 , the comparison between Caelus and CFL3D is made for Grid5 and as can be seen, a very good agreement is obtained over the entire region of the bump.
kOmega SST
The skinfriction distribution variation for different grids obtained from \(k\omega~\rm{SST}\) model is shown in Figure 83.
In Figure 84 , the skinfriction comparison between Caelus and CFL3D is made for Grid5 and is shown.
Conclusions
The steady turbulent flow simulation over a twodimensional bump was carried out using Caelus 9.04 employing simpleSolver. The solutions were obtained with two turbulence models, implemented inhouse and the results were verified against CFL3D data. The comparison was found to be in good agreement with CFL3D suggesting that the turbulence model implementation is accurate in Caelus.
Twodimensional NACA 0012 AirfoilÂ¶
Nomenclature
Turbulent, incompressible flow over a twodimensional NACA airfoil
Symbol 
Definition 
Units (SI) 

\(a\) 
Speed of sound 
\(m/s\) 
\(c_f\) 
Skin friction coefficient 
Nondimensional 
\(c_p\) 
Pressure coefficient 
Nondimensional 
\(C\) 
Chord length 
\(m\) 
\(I\) 
Turbulent intensity 
Percentage 
\(k\) 
Turbulent kinetic energy 
\(m^2/s^2\) 
\(M_\infty\) 
Freestream Mach number 
Nondimensional 
\(p\) 
Kinematic pressure 
\(Pa/\rho~(m^2/s^2)\) 
\(Re_L\) 
Reynolds number 
\(1/m\) 
\(T\) 
Temperature 
\(K\) 
\(u\) 
Local velocity in xdirection 
\(m/s\) 
\(w\) 
Local velocity in zdirection 
\(m/s\) 
\(U\) 
Freestream velocity 
\(m/s\) 
\(x\) 
Distance in xdirection 
\(m\) 
\(z\) 
Distance in ydirection 
\(m\) 
\(y^+\) 
Wall distance 
Nondimensional 
\(\alpha\) 
Angle of attack 
Degrees 
\(\mu\) 
Dynamic viscosity 
\(kg/m~s\) 
\(\nu\) 
Kinematic viscosity 
\(m^2/s\) 
\(\tilde{\nu}\) 
Turbulence field variable 
\(m^2/s\) 
\(\rho\) 
Density 
\(kg/m^3\) 
\(\omega\) 
Specific dissipation rate 
\(1/s\) 
\(\infty\) 
Freestream conditions 

\(t\) (subscript) 
Turbulent property 
Introduction
This case deals with the steady turbulent incompressible flow over a twodimensional NACA 0012 airfoil. The study is conducted at two angles of attack, \(\alpha = 0^\circ\) and \(\alpha = 10^\circ\) respectively. The simulations are carried out over a series of four grids and compared with CFL3D data and the results are also compared with the experimental data. This exercise therefore verifies and validates the turbulence models used through the distribution of skinfriction coefficient (\(c_f\)) and pressure coefficient (\(c_p\)) over the airfoil surface.
Problem definition
This verification and validation exercise is based on the Turbulence Modeling Resource case for a NACA 0012 airfoil and follows the same flow conditions used in the incompressible limit. However, the numerical code CFL3D uses a freestream Mach number (\(M_\infty\)) of 0.15. In Figure 85 the schematic of the airfoil is shown. Note that the 2D plane in Figure 85 is depicted in \(xz\) directions as the computational grid also follows the same 2D plane.
The length of the airfoil chord is C = 1.0 m, wherein, x = 0 is at the leading edge and the Reynolds number is \(6 \times 10^6\) and \(U_\infty\) is the freestream velocity. For \(\alpha = 10^\circ\), the velocity components are evaluated in order to have a same resultant freestream velocity. The freestream flow temperature in this case is assumed to be 300 K and for Air as a perfect gas, the speed of sound (\(a\)) can be evaluated to 347.180 m/s. Based on the freestream Mach number and speed of sound, freestream velocity can be evaluated to \(U_\infty\) = 52.077 m/s. Using the value of velocity and Reynolds number, kinematic viscosity can be calculated. The following table summarises the freestream conditions used:
Fluid 
\(Re_L~(1/m)\) 
\(U_\infty~(m/s)\) 
\(p_\infty~(m^2/s^2)\) 
\(T_\infty~(K)\) 
\(\nu_\infty~(m^2/s)\) 

Air 
\(6 \times 10^6\) 
52.0770 
\((0)\) Gauge 
300 
\(8.679514\times10^{6}\) 
To get a freestream velocity of \(U_\infty\) = 52.077 m/s at \(\alpha = 10^\circ\), the velocity components in \(x\) and \(z\) are resolved. These are provided in the table below
\(\alpha~\rm{(Degrees)}\) 
\(u~(m/s)\) 
\(w~(m/s)\) 

\(0^\circ\) 
52.0770 
0.0 
\(10^\circ\) 
51.2858 
9.04307 
Note that in an incompressible flow, the density (\(\rho\)) does not vary and a constant density can be assumed throughout the calculation. Further, since temperature is not considered here, the viscosity is also held constant. In Caelus for incompressible flow simulations, pressure and viscosity are always specified as kinematic.
Turbulent Properties for Spalartâ€“Allmaras model
The turbulent inflow boundary conditions used for the Spalartâ€“Allmaras model were calculated as \(\tilde{\nu}_{\infty} = 3 \cdot \nu_\infty\) and subsequently turbulent eddy viscosity was evaluated. The following table provides the values of these used in the current simulations:
\(\tilde{\nu}_\infty~(m^2/s)\) 
\(\nu_{t~\infty}~(m^2/s)\) 

\(2.603854 \times 10^{5}\) 
\(1.8265016 \times 10^{6}\) 
Turbulent Properties for komega SST model
The turbulent inflow boundary conditions used for \(k\omega~\rm{SST}\) were calculated as follows and is as given in Turbulence Modeling Resource
Note that the dynamic viscosity in the above equation is obtained from Sutherland formulation and density is evaluated as \(\rho = \mu / \nu\). In the below table, the turbulent properties used in the current simulations are provided.
\(I\) 
\(k_{\infty}~(m^2/s^2)\) 
\(\omega_{\infty}~(1/s)\) 
\(\nu_{t~\infty}~(m^2/s)\) 

\(0.052\%\) 
\(1.0999 \times 10^{3}\) 
\(13887.219\) 
\(7.811564 \times 10^{8}\) 
Computational Domain and Boundary Conditions
The computational domain used for the airfoil simulations and the details of the boundaries are shown in Figure 86 for a \(xz\) plane. The leading edge and the trailing edge extends between \(0 \leq x \leq 1.0~m\) and the entire airfoil has a noslip boundary condition. The farfield domain extends by about 500 chord lengths in the radial direction and the inlet is placed for the entire boundary highlighted in green. The outlet boundary is placed at the exit plane, which is at \(x \approx 500~m\).
Boundary Conditions and Initialisation
 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 the above table)
\(k\omega~\rm{SST}\) (Fixed uniform values of \(k_{\infty}\), \(\omega_{\infty}\) and \(\nu_{t~\infty}\) as given in the above table)
 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:
\(\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 Gauge pressure
Turbulence:
Spalartâ€“Allmaras (Fixed uniform values of \(\nu_{t~\infty}\) and \(\tilde{\nu}_{\infty}\) as given in the above table)
\(k\omega~\rm{SST}\) (Fixed uniform values of \(k_{\infty}\), \(\omega_{\infty}\) and \(\nu_{t~\infty}\) as given in the above table)
Computational Grid
The 3D computational grid for the NACA 0012 airfoil was obtained from Turbulence Modeling Resource as a Plot3D format. Using Pointwise it was then converted to Caelus format. As indicated earlier, the twodimensional plane of interest in the Plot3D grid is in \(xz\) directions. As the flow is considered here to be twodimensional, and simpleSolver being a 3D solver, the two \(xz\) planes are specified with empty boundary conditions consequently treating as symmetry flow in \(y\) direction. To study the sensitivity of the grid, four grids were considered from the original set of five, in which the coarsest grid was excluded from this study. Details of the different grids used are given in the below table. Not that for both angles of attack, same grid is used.
Grid 
Cells over airfoil 
Cells in normal direction 
Total 
\(y^+\) 

Grid2 
128 
64 
14,336 
0.465 
Grid3 
256 
128 
57,344 
0.209 
Grid4 
512 
256 
229,376 
0.098 
Grid5 
1024 
512 
917,504 
0.047 
The below Figure 87 shows the 2D grid in \(xz\) plane for Grid3 and the refinement around the airfoil is shown in the inset. Sufficient refinement can be seen in the wall normal direction and all the grid have a \(y^+ < 1\) and no wall function is used for the airfoil surface throughout the current verification and validation cases.
Results and Discussion
The solution to the turbulent flow over the NACA 0012 airfoil was obtained using Caelus 9.04. SimpleSolver was used and the solutions were run sufficiently long until the residuals for pressure, velocity and turbulence quantities were less than \(1 \times 10^{6}\). The finite volume discretization of the gradient of pressure and velocity was carried out using the linear approach. Where as the divergence of velocity and mass flux was carried out through the linear upwind method. However, for the divergence of the turbulent quantities, upwind approach was utilised and linear approach for the divergence of the Reynolds stress terms. For the discretization of the Laplacian terms, again linear corrected method was used. For some grids having greater than 50 degree nonorthogonal angle, linear limited with a value of 0.5 was used for the Laplacian of the turbulent stress terms.
The verification results are shown firstly for both angles of attack and is followed by the experimental validation data.
Verification results: Spalartâ€“Allmaras
The following Figure 88 and Figure 89 shows the skinfriction distribution over the upper surface for \(\alpha=0^\circ\) and \(\alpha=10^\circ\) from Caelus for different grids. In both cases, Grid4 and Grid5 essentially produces the same solution suggesting a gridindependence solution is obtained.
In Figure 90 and Fig. #turbulentairfoilcaeluscfl3dsacc10 , the skinfriction is compared with CFL3D on Grid4. As can be seen, a very good agreement between the two codes can be seen.
Verification results: komega SST
The skinfriction distribution obtained from using \(k\omega~\rm{SST}\) turbulence model for \(\alpha=0^\circ\) and \(\alpha=10^\circ\) is shown below for different grids. The gridsensitivity behaviour is very similar to the Spalartâ€“Allmaras turbulence case and no change is seen between Grid4 and Grid5.
The comparison of the skinfriction with CFL3D using \(k\omega~\rm{SST}\) is shown in Figure 94 and Figure 95 for both angle of attacks and similar to the previous case, a very good agreement between the two can be seen.
Experimental validation
Here, the Caelus data is compared with the pressurecoefficient (\(c_p\)) obtained experimentally by Gregory, N. and Oâ€™Reilly, C. L [5] for both angles of attack over the upper surface. In addition, the data obtained from CFL3D is also included for verification. There is a very good agreement with the current Caelus and experiments which indicates that the correct turbulence equations are being solved in both Spalartâ€“Allmaras and \(k\omega~\rm{SST}\) models.
Conclusions
Verification and validation over a twodimensional NACA 0012 airfoil for turbulent inflow conditions were carried out using Caelus 9.04 employing simpleSolver. Two turbulence models that are implemented inhouse were used and the solutions were verified with CFL3D data and subsequently validated with the experimental pressure coefficient values. The results were found to be in very good agreement suggesting that the turbulence modelling implementation is appropriate and solves accurately.
Twodimensional Convex CurvatureÂ¶
Turbulent, incompressible flow in a twodimensional convex curvature channel
Nomenclature
Symbol 
Definition 
Units (SI) 

\(a\) 
Speed of sound 
\(m/s\) 
\(c_f\) 
Skin friction coefficient 
Nondimensional 
\(c_p\) 
Pressure coefficient 
Nondimensional 
\(I\) 
Turbulent intensity 
Percentage 
\(k\) 
Turbulent kinetic energy 
\(m^2/s^2\) 
\(M_i\) 
Inlet Mach number 
Nondimensional 
\(p\) 
Kinematic pressure 
\(Pa/\rho~(m^2/s^2)\) 
\(Re_L\) 
Reynolds number 
\(1/m\) 
\(T\) 
Temperature 
\(K\) 
\(u\) 
Local velocity in xdirection 
\(m/s\) 
\(w\) 
Local velocity in zdirection 
\(m/s\) 
\(U\) 
Inlet velocity 
\(m/s\) 
\(x\) 
Distance in xdirection 
\(m\) 
\(y\) 
Distance in ydirection 
\(m\) 
\(y^+\) 
Wall distance 
Nondimensional 
\(\alpha\) 
Angle of attack 
Degrees 
\(\mu\) 
Dynamic viscosity 
\(kg/m~s\) 
\(\nu\) 
Kinematic viscosity 
\(m^2/s\) 
\(\tilde{\nu}\) 
Turbulence field variable 
\(m^2/s\) 
\(\rho\) 
Density 
\(kg/m^3\) 
\(\omega\) 
Specific dissipation rate 
\(1/s\) 
\(i\) 
Inlet conditions 

\(t\) (subscript) 
Turbulent property 

\(ref\) 
Reference pressure 
\(Pa\) 
Introduction
In this case, the turbulent incompressible flow in a constantarea duct having a convex curvature is investigated as a part of verification and validation. The effect of curvature on the capability of turbulence model to predict the boundary layer accurately is of primary concern in this study. Similar to previous cases, the simulations are carried out over a series of four grids that tests the sensitivity of solution as the grid is refined. As with the earlier cases, SpalartAllmaras (SA) and \(k\omega~\rm{SST}\) turbulence models are tested. However due to the presence of strong curvature, a variant of Spalartâ€“Allmaras turbulence model which accounts for rotational and curvature effects is additionally considered. This model is referred to as â€śSpalartâ€“Allmaras Rotational/Curvatureâ€ť (SARC) [9] . The results are then verified against CFL3D data and validated with experimental pressure distributions.
Problem definition
This verification and validation exercise is based on the Turbulence Modeling Resource case and follows the same conditions used in the incompressible limit. As CFL3D is a compressible CFD code, the original simulation used a inlet Mach number (\(M_i\)) of 0.093. A schematic of the geometric configuration is shown below and is depicted in \(xz\) plane as the computational grid also follows the same 2D plane of flow.
As can be seen above, the 2D duct has a rapid bend of \(\alpha = 30^\circ\) after about a distance of 1.4 m and the downstream extends up to 1.6 m in length. The crosssection of the duct is 0.127 m and the inner radius and outer radius of curvature are 0.127 m and 0.254 m respectively. The flow has a Reynolds number of \(2.1 \times 10^6\) and U is the inlet velocity. To achieve a desired inlet velocity at an angle of 30 degrees, the velocity components are evaluated. The inlet temperature for this case is T = 293 K and for Air as a perfect gas, the speed of sound (\(a\)) can be evaluated to 343.106 m/s. Based on the inlet Mach number and speed of sound, the inlet velocity can be calculated to U = 31.908 m/s. Using velocity and Reynolds number kinematic viscosity can be calculated. The table below summarises the inlet conditions used:
Fluid 
\(Re_L~(1/m)\) 
\(U~(m/s)\) 
\(p~(m^2/s^2)\) 
\(T~~(K)\) 
\(\nu_i~(m^2/s)\) 

Air 
\(2.1 \times 10^6\) 
31.9088 
\((0)\) Gauge 
293 
\(1.519470\times10^{5}\) 
In order 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 below
\(\alpha~\rm{(Degrees)}\) 
\(u~(m/s)\) 
\(w~(m/s)\) 

\(30^\circ\) 
27.63389 
15.95443 
It should be noted that in an incompressible flow, the density (\(\rho\)) does not vary and a constant density assumption is valid throughout the calculation. Further, the temperature field is not solved and therefore its influence on viscosity can be neglected and a constant viscosity can be used. In Caelus for incompressible flow simulations, pressure and viscosity are always specified as kinematic.
Turbulent Properties for Spalartâ€“Allmaras and Spalartâ€“Allmaras Curvature Correction model
The inflow conditions used for turbulent properties in SpalartAllamaras model was calculated as \(\tilde{\nu}_{i} = 3 \cdot \nu_i\) and subsequently turbulent eddy viscosity was evaluated. The following table provides the values of these used in the current simulations:
\(\tilde{\nu}_i~(m^2/s)\) 
\(\nu_{t~i}~(m^2/s)\) 

\(4.558411 \times 10^{5}\) 
\(3.197543 \times 10^{6}\) 
Turbulent Properties for komega SST model
The turbulent inflow boundary conditions used for \(k\omega~\rm{SST}\) were calculated as follows and is as given in Turbulence Modeling Resource
Note that the dynamic viscosity in the above equation is obtained from Sutherland formulation and density is evaluated as \(\rho = \mu / \nu\). In the below table, the turbulent properties used in the current simulations are provided
\(I\) 
\(k_{i}~(m^2/s^2)\) 
\(\omega_{i}~(1/s)\) 
\(\nu_{t~i}~(m^2/s)\) 

\(0.083\%\) 
\(1.0521 \times 10^{3}\) 
\(7747.333\) 
\(1.36756 \times 10^{7}\) 
Computational Domain and Boundary Conditions
The computational domain for the duct is quite simple and follows the geometry as is shown in Figure 99. The walls are modelled as noslip boundary and is highlighted in blue and the outlet is placed at the end of the duct.
Boundary Conditions and Initialisation
 Inlet
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 gradient
Turbulence:
SA & SARC (Fixed uniform values of \(\nu_{t~i}\) and \(\tilde{\nu}_{i}\) as given in the above table)
\(k\omega~\rm{SST}\) (Fixed uniform values of \(k_{i}\), \(\omega_{i}\) and \(\nu_{t~i}\) as given in the above table)
 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 the above table)
\(k\omega~\rm{SST}\) (Fixed uniform values of \(k_{i}\), \(\omega_{i}\) and \(\nu_{t~i}\) as given in the above table)
Computational Grid
The computational grid in 3D for the convex curvature duct was obtained from Turbulence Modeling Resource as a Plot3D format. The same was used in Caelus by converting it in the required format with the use of Pointwise. Since the flow field is assumed to the twodimensional, the 2D computational plane of interest is in \(xz\) directions. SimpleSolver is a 3D solver, therefore the two additional \(xz\) planes are specified with empty boundary conditions. To look at the effect of grid sensitivity, four grids were considered from the original set of five, while the coarsest grid was not included in this study. The table below gives the details of different grids used.
Grid 
Cells in streamwise direction 
Cells in normal direction 
Total 

Grid2 
128 
48 
6144 
Grid3 
256 
96 
98304 
Grid4 
512 
192 
98,304 
Grid5 
1024 
384 
393,216 
The 2D convex curvature grid is shown in Figure 100 below, for Grid3 in \(xz\) plane and the inset shows the grids in the vicinity of the strong curvature. Grids are refined close to the wall in order to capture the turbulent boundary layer and all grids have a \(y^+ < 1\) and no wall function is used through the validation and verification of this configuration.
Results and Discussion
The turbulent flow inside the convex curvature duct was simulated using Caelus 9.04 through the use of simpleSolver. The solutions were run until the residuals for pressure, velocity and turbulent quantities were less than \(1 \times 10^{6}\). The finite volume discretization of the gradient of pressure and velocity was carried out using the linear approach. Where as the divergence of velocity and mass flux was carried out through the linear upwind method. However, for the divergence of the turbulent quantities, upwind approach was utilised and linear approach for the divergence of the Reynolds stress terms. For the discretization of the Laplacian terms, again linear corrected method was used. For some grids having greater than 50 degree nonorthogonal angle, linear limited with a value of 0.5 was used for the Laplacian of the turbulent stress terms.
The verification results of the turbulence model are discussed first, which is then followed by the experimental validation.
Verification results Spalartâ€“Allmaras (SA)
In Figure 101, the skinfriction distribution obtained over the lower wall of the duct is shown from Caelus for different grids. As can be seen, there is very little difference in skinfriction variation among the different grids. The oscillatory behaviour noted at Grid5 very close to the corner is also apparent in CFL3D data as will be see in Figure 101.
The skinfriction coefficient comparison with CFL3D is shown in Figure 102. It should be noted that the available solution from CFL3D was for Grid4 and hence to be consistent, Grid4 solution from Caelus is used for comparison. In the vicinity of strong curvature region, Caelus compares very well with CFL3D, however both upstream and downstream, there seems to be some difference in the solution.
Verification results Spalartâ€“Allmaras Rotational/Curvature (SARC)
The following Figure 103 shows the grid sensitivity and verification with CFL3D data respectively. The solution that have been used for verification is obtained from Grid4 and the trends are similar to what is noted for the SA model.
Verification results kOmega SST
In Figure 105, the skinfriction sensitivity is shown over the lower wall obtained using Caelus with the \(k\omega~\rm{SST}\) model. After Grid3, not much difference in values can be noted. With Grid5, however some oscillations can be see upstream and in the vicinity of the curvature.
The skinfriction comparison between Caelus and CFL3D is shown in Figure 106. A very good agreement between the two is obtained.
Experimental validation
This section details the experimental validation carried out for Caelus and both skinfriction and pressure coefficients obtained experimentally by Smits, A. J et al. [8] are compared. Further, CFL3D is also included. In Figure 107, skinfriction distribution obtained from Caelus using different turbulence models is compared with the experiments. Both SARC and \(k\omega~\rm{SST}\) has a fair agreement with experiments down stream of the curvature, more so with the SARC model. However upstream they all seem to predict nearly the same values.
Figure 108 shows the comparison of pressurecoefficient (\(c_p\)) distribution with experiments and CFL3D data on the lower surface. Firstly, the solutions obtained from Caelus with the three turbulence models essentially produces the same values and matches exactly with the CFL3D data. In comparison with experiments, the agreement is very good in the upstream, vicinity of the curvature and downstream and, identical to CFL3Dâ€™s behaviour. Note that for obtaining the pressurecoefficient (\(c_p\)) values, a reference pressure (\(p_{ref}\)) is needed. However, this is not specified in Turbulence Modeling Resource for this case and hence a value of 145 Pa has been used.
Experimental skinfriction data [8] is also available over the upper surface of the duct and is used to compare the Caelus results. Figure 109 shows the comparison. Similar to the behaviour noted for the lower surface, the SARC model tends to be closer to the experimental data. In general, all the three turbulence model have similar trends and agrees very closely with the CFL3D data.
Conclusions
A detailed verification and validation of a turbulent flow in a convex curvature duct were carried out using Caelus 9.04 and simpleSolver. Here, three turbulence models used and the solutions were verified against CFL3D data. As a part of validation, Caelus results were compared with the experimental data obtained on both lower and upper surfaces. The comparison was good with both CFL3D as well as with experiments. This suggests that the implementation of the turbulence models is correct and is being solved accurately.
Twodimensional Backward Facing StepÂ¶
Turbulent, incompressible flow over a twodimensional backward facing step
Nomenclature
Symbol 
Definition 
Units (SI) 

\(a\) 
Speed of sound 
\(m/s\) 
\(c_f\) 
Skin friction coefficient 
Nondimensional 
\(c_p\) 
Pressure coefficient 
Nondimensional 
\(H\) 
Step height 
\(m\) 
\(I\) 
Turbulent intensity 
Percentage 
\(k\) 
Turbulent kinetic energy 
\(m^2/s^2\) 
\(M_\infty\) 
Freestream Mach number 
Nondimensional 
\(p\) 
Kinematic pressure 
\(Pa/\rho~(m^2/s^2)\) 
\(Re_H\) 
Reynolds number based on step height 
Nondimensional 
\(T\) 
Temperature 
\(K\) 
\(u\) 
Local velocity in xdirection 
\(m/s\) 
\(u_*\) 
Frictional velocity 
\(m^2/s\) 
\(U\) 
Freestream velocity in xdirection 
\(m/s\) 
\(U_{ref}\) 
Reference velocity 
\(m/s\) 
\(x\) 
Distance in xdirection 
\(m\) 
\(z\) 
Distance in zdirection 
\(m\) 
\(y^+\) 
Wall distance 
Nondimensional 
\(\epsilon\) 
Turbulent dissipation 
\(m^2/s^3\) 
\(\mu\) 
Dynamic viscosity 
\(kg/m~s\) 
\(\nu\) 
Kinematic viscosity 
\(m^2/s\) 
\(\tilde{\nu}\) 
Turbulence field variable 
\(m^2/s\) 
\(\rho\) 
Density 
\(kg/m^3\) 
\(\tau_w\) 
Wall shear stress 
\(kg/m~s^2\) 
\(\omega\) 
Specific dissipation rate 
\(1/s\) 
\(\infty\) 
Freestream conditions 

\(t\) (subscript) 
Turbulent property 
Introduction
This study investigates steady turbulent, incompressible flow over a twodimensional backward facing step at zero angle of incidence. Unlike the previous cases, the efficacy of wall functions in separated flow is evaluated for different turbulent models. The validation of these wall functions are carried out through the comparison of skinfriction coefficient (\(c_f\)) and pressure coefficient (\(c_p\)) downstream of the step with those of the experimental data. In addition to Spalartâ€“Allmaras and \(k  \omega~\rm{SST}\), Realizable \(k\epsilon\) turbulence model was also considered in this exercise.
Problem definition
The backward facing step configuration is obtained from the Turbulence Modeling Resource and is a widely considered case for the purpose of verification and validation. This particular study is based on the experiments carried out by Driver and Seegmiller [15]. The schematic of the step configuration in Figure 110 below as considered in the Turbulence Modeling Resource.
The height of the step (H) is chosen to be 1.0 m and is located at x = 0 m. Upstream of the step, the plate extends to 110 step heights such that a fully developed turbulent boundary layer with proper thickness exists at separation (x = 0 m). This is followed by a downstream plate which extends up to 50 step heights. The flow has a Reynolds number of 36000 based on the step height with a freestream velocity (\(U_\infty\)) in the xdirection. The inflow is assumed to be Air as a perfect gas with a temperature of 298 K, giving the value of speed of sound (\(a\)) at 346.212 m/s. A reference Mach number (\(M_\infty\)) of 0.128 is used to obtain the freestream velocity, which was \(U_\infty\) = 44.315 m/s. The kinematic viscosity was then evaluated based on the velocity and the Reynolds number. The following table summarises the freestream conditions used for the backward facing step.
Fluid 
\(Re_H\) 
\(U_\infty~(m/s)\) 
\(p_\infty~(m^2/s^2)\) 
\(T_\infty~(K)\) 
\(\nu_\infty~(m^2/s)\) 

Air 
\(36000\) 
44.31525 
\((0)\) Gauge 
298.330 
\(1.230979\times10^{3}\) 
As with all the previous cases, the flow is incompressible and therefore the density (\(\rho\)) does not change throughout the simulation. Further, the temperature is not accounted and has no influence on viscosity and is also held constant. However note that in Caelus, pressure and viscosity for incompressible flow are always specified as kinematic.
Turbulent Properties for Spalartâ€“Allmaras model
The turbulent boundary conditions at the freestream used for SpalartAllmaras model were calculated according to \(\tilde{\nu}_{\infty} = 3 \cdot \nu_\infty\) and turbulent eddy viscosity was subsequently evaluated. In the following table, the values used in the current simulation are provided:
\(\tilde{\nu}_\infty~(m^2/s)\) 
\(\nu_{t~\infty}~(m^2/s)\) 

\(3.692937 \times 10^{3}\) 
\(2.590450 \times 10^{4}\) 
Turbulent Properties for komega SST model
The turbulent inflow boundary conditions used for \(k\omega~\rm{SST}\) were calculated as follows and is as given in Turbulence Modeling Resource
The dynamic viscosity in the above equation is obtained from Sutherland formulation and the density is calculated through \(\mu/\nu\). In the table below, the turbulent properties used for the current \(k\omega~\rm{SST}\) simulation are provided
\(I\) 
\(k_{\infty}~(m^2/s^2)\) 
\(\omega_{\infty}~(1/s)\) 
\(\nu_{t~~\infty}~(m^2/s)\) 

\(0.061\%\) 
\(1.0961 \times 10^{3}\) 
\(97.37245\) 
\(1.10787 \times 10^{5}\) 
Turbulent Properties for Realizable kepsilon model
The turbulent inflow properties used for Realizable \(k\epsilon\) model were evaluated as follows
where, \(\lambda\) is the turbulent length scale and is evaluated at 0.22 of the boundary layer thickness at separation (1.5H) and \(C_\mu\) is a model constant with a value 0.09. The following table provides the evaluated turbulent properties. Note that the turbulent intensity was assumed to be 1 % for this particular model.
\(I\) 
\(k_{\infty}~(m^2/s^2)\) 
\(\epsilon_{\infty}~(m^2/s^3)\) 
\(\nu_{t\infty}~(m^2/s)\) 

\(1\%\) 
\(294.57 \times 10^{3}\) 
\(0.079598\) 
\(98.11 \times 10^{3}\) 
Computational Domain and Boundary Conditions
The computational domain simply follows the step geometry for the entire bottom region. In Figure 111 below, the boundary details in twodimensions (\(xz\) plane) are shown. The walls of the upstream plate, step and the downstream plate that extend between \(110 \leq x \leq 50~m\) are modelled as noslip wall boundary condition. Similarly, the top plate is also modelled as a noslip wall. Upstream of the leading edge, that is, \(x \leq 110\) symmetry boundary extends for a length of 20 step heights and the inlet boundary is placed at the start of the symmetry. The outlet is placed at the end of the downstream plate, which is at \(x = 50~m\).
Boundary Conditions and Initialisation
 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 the above table)
\(k\omega~\rm{SST}\) (Fixed uniform values of \(k_{\infty}\), \(\omega_{\infty}\) and \(\nu_{t~\infty}\) as given in the above table)
Realizable \(k\epsilon\) (Fixed uniform value of \(k_{\infty}\), \(\epsilon_{\infty}\) and \(\nu_{t_\infty}\) as given in the above table)
 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 an initial value of \(\nu_t=0\)
\(\tilde{\nu}\): type fixedValue with a value of \(\tilde{\nu}=0\)
\(k\omega~\rm{SST}\):
\(k\): type kqRWallFunction with an initial value of \(k_{\infty}\)
\(\omega\): type omegaWallFunction with an initial value of \(\omega_{\infty}\)
\(\nu_t\): type nutUWallFunction with an initial value of \(\nu_t=0\)
Realizable \(k\epsilon\):
\(k\): type kqRWallFunction with an initial value of \(k_{\infty}\)
\(\epsilon\): type epsilonWallFunction with an initial value of \(\epsilon=0\)
\(\nu_t\): type nutUWallFunction with an initial 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~m/s\) in \(x\) direction
Pressure: Zero Gauge pressure
Turbulence:
Spalartâ€“Allmaras (Fixed uniform values of \(\nu_{t~\infty}\) and \(\tilde{\nu}_{\infty}\) as given in the above table)
\(k\omega~\rm{SST}\) (Fixed uniform values of \(k_{\infty}\), \(\omega_{\infty}\) and \(\nu_{t~\infty}\) as given in the above table)
Realizable \(k\epsilon\) (Fixed uniform values of \(k_{\infty}\), \(\epsilon_{\infty}\) and \(\nu_{t~\infty}\) as given in the above table)
Computational Grid
The 3D computational grid was generated using Pointwise. and was converted to Caelus format. Note that the plane of interest is in \(xz\) directions. Further, since the flow field is assumed to be twodimensional, the two \(xz\) obtained as a result of a 3D grid are specified with empty boundary condition. This forces the 3D solver, simpleSolver to treat the flow in \(y\) direction as symmetry. The following figure shows the 2D grid. As noted earlier, the purpose of this exercise is to validate the wall functions and hence the grid is designed for a \(y^+ = 30\). In order to design such a gird, the first cell height (\(\Delta z\)) from wall in the wall normal (\(z\)) direction was obtained from the following set of equations
where, \(u_*\) is the frictional velocity given by
In the above equation, \(\tau_w\) is the shearstress at the wall and was be estimated using the skinfriction (\(c_f\)) relation, given as
where, \(c_f\) was obtained for the flatplate as given in Schlichting [7] and is shown below
where, \(Re_x\) is the Reynolds number based on the length of