4. Tutorials

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. Note that each tutorial contains information about the problem description, case set-up, solver settings and post-processing details. It is suggested that user follow these steps methodically.

4.1. Laminar Flows

4.1.1. Laminar Flat Plate (only applicable to Applied CCM subscription customers)

In this tutorial, simulation of laminar incompressible flow over a two-dimensional sharp leading-edge flat plate using Caelus 4.10 is introduced here. First, pre-requisites 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.

4.1.1.1. Objectives

With the completion of this tutorial, the user will be familiar with setting up Caelus simulation for steady, laminar, incompressible flow over flat-plates in two-dimensions and subsequently post-processing 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

4.1.1.2. Pre-requisites

It is assumed that the user is familiar with the Linux command line environment using a terminal and 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.

4.1.1.3. Background

The flow over a flat-plate 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 sharp-leading edge plate is solved in a time-accurate 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 non-dimensional shear stress distribution (skin-friction coefficient). For more details, the user is suggested to refer the validation of flat-plate in section Flat-plate (only applicable to Applied CCM subscription customers).

The length of the flat-plate 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 Fig. 66, 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\).

Figure 66: Schematic of the flat-plate flow

In the following table, details of the freestream conditions are provided.

Freestream conditions
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}\)

4.1.1.4. 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 flat-plate. The below (Fig. 67) shows the details of the boundaries that will be used in two-dimensions (\(x-y\) plane). First, the flat-plate, 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 no-slip boundary (\(u, v, w = 0\)). Upstream of the leading edge, a slip boundary will be used to simulate a freestream flow approaching the flat-plate. However, downstream of the plate, it would be ideal to further extend up to three plate lengths (highlighted in green) with no-slip 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 slip-wall at \(x = -0.06~m\) and the outlet at the exit plane of the no-slip wall (green region) at \(x = 1.2192~m\). Both inlet and outlet boundary are between \(0\leq y \leq 0.15~m\). A slip-wall condition is used for the entire top boundary.

Figure 67: Flat-plate computational domain

The 2D structured grid is shown in Fig. 68. 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 one-cell thick. The final 3D grid should be then exported to Caelus format (polyMesh). The two \(x-y\) planes obtained as a result of grid extrusion need boundary conditions to be specified. As the flow over a flat-plate 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.

_images/fp-grid-tutorials.png

Figure 68: Flat-plate computational grid

The flat-plate has a total of 400 cells over the region of interest between \(0 \leq x \leq 0.3048~m\) and 286 cells in the no-slip 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.

4.1.1.5. Problem definition

Several important instructions would be shown here to set-up the flat-plate problem along with the detail of configuration files used. A full working case can be found in:

/tutorials-cases/flat-plate/flat-plate

However, it will be shown here to complete the tutorial assuming no prior set-up or data is available.

Directory Structure

First, let us create a directory named my-flat-plate where we intend to carry out the simulation via a terminal window.

cd flat-plate
mkdir my-flat-plate
cd my-flat-plate/

Note

All commands shown here are entered in a terminal window, unless otherwise mentioned

In order to set-up the problem in my-flat-plate, we need to have few more sub-directories where relevant files will be placed. Caelus requires time, constant and system sub-directories. Since we will begin the simulation at time \(t = 0~s\), the time sub-directory should be just 0. This can be obtained by

cd my-flat-plate/
mkdir 0
mkdir constant
mkdir system

ls
0  constant  system

The 0 sub-directory is where additional two files, p and U for pressure (\(p\)) and velocity (\(U\)) respectively need to be placed. The contents of these two files sets the dimensions, initialisation and boundary conditions to the problem, which also form three principle entries required. Let us begin with creating the two files in 0 by

cd my-flat-plate/
touch 0/p
touch 0/U

cd 0/
ls
p  U

It should be noted that Caelus is case sensitive and therefore the user should set-up the directories, files and the contents identical to what is shown here.

Boundary Conditions and Solver Attributes

Boundary Conditions

First let us look at setting up the boundary conditions. Referring back to Fig. 67, 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
  • No-slip 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

Beginning with the pressure (\(p\)), let us open the file p using a text editor and add the following contents to it.

/*------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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;
    }
    symm-left
    {
        type				empty;
    }
    symm-right
    {
        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. Next we explain 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 non-uniform. 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, let us open the file U and add the following to it

/*------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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;
    }
    symm-left
    {
        type				empty;
    }
    symm-right
    {
        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 \(x-y\) plane obtained due to grid extrusion have been named as symm-left and symm-right with specifying empty boundary conditions forcing Caelus to assume the flow to be in two-dimensions. This completes the setting up of boundary conditions.

Grid file and Physical Properties

The flat-plate grid files that is generated for Caelus format has to be placed in the polyMesh sub-directory which resides in the constant directory. This can be done as follows

cd my-flat-plate/
mkdir constant/polyMesh

The grid files can now be saved in polyMesh sub-directory. In addition, the physical properties are specified in various different files present in the directory constant. Let us navigate to constant and create the following files

touch constant/transportProperties
touch constant/turbulenceProperties

cd constant/
ls
polyMesh  transportProperties  turbulenceProperties

In the first file transportProperties, transport model and kinematic viscosity are specified. The contents of this file are as follows

/*------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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.5896306e-5;

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 final in this class is the turbulenceProperties file, where the type of simulation is specified as

/*------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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.

Controls and Solver Attributes

The files required to control the simulation and specifying the type of discretization method along with the linear solver settings need to present in the system directory. Let us begin with creating these files by navigating to system as

touch system/controlDict
touch system/fvSchemes
touch system/fvSolution

cd system/
ls
controlDict  fvSchemes  fvSolution

Let us first begin adding the input data to controlDict file as below

/*------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
	version				2.0;
	format				ascii;
	class				dictionary;
	location			"system";
	object				controlDict;
}

//-------------------------------------------------------------------------------

application				slim;

startFrom				startTime;

startTime				0;

stopAt					endTime;

endTime					0.3;

deltaT					0.00001;

writeControl				runTime;

writeInterval				0.01;

purgeWrite				0;

writeFormat				ascii;

writePrecision				6;

writeCompression 			uncompressed;

timeFormat				general;

timePrecision				6;

runTimeModifiable			true;

Here, the application slim refers to the slim 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. The simulation is run for 0.3 seconds specifying through the keywords stopAt and endTime. Since SLIM solver is time-accurate, we also need to set the time-step via deltaT. The results are written at every 0.01 seconds via writeControl and writeInterval keywords.

The discretization schemes and parameters are specified through the fvSchemes file, shown below and are the contents and should be copied

/*------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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	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;
}

interpolationSchemes
{
	default				linear;
	interpolate(HbyA)		linear;
}

snGradschemes
{
	default				corrected;
}

fluxRequired
{
	default				no;
	p					;
}

Here, the discretization schemes for finite volume discretization of time-derivative, gradient, divergence and Laplacian terms are specified.

The linear solver controls and tolerances are set in fvSolution as given below

/*------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
	version				2.0;
	format				ascii;
	class				dictionary;
	location			"system";
	object				fvSolution;
}

//------------------------------------------------------------------------------

solvers
{
	p
	{
		solver			PCG;
		preconditioner		DIC;
		tolerance		1e-8;
		relTol			0;
	}
	pFinal
	{
		solver			PCG;
		preconditioner		DIC;
		tolerance		1e-8;
		relTol			0;
	}
	U
	{
		solver			PBiCG;
		preconditioner		DILU;
		tolerance		1e-8;
		relTol			0;
	}
}
SLIM
{
	nNonOrthogonalCorrectors	0;
	pRefCell			0;
	pRefValue			0;
}

Different liner solvers are used here to solve pressure and velocity. Also, since the gird is perfectly orthogonal, the orthogonal correction specified via nNonOrthogonalCorrectors is set to 0.

The set-up of the directory structure along with the relevant files are completed. Let us verify again by typing the following command and the directory tree should be identical to the one shown below

cd my-flat-plate/
tree
.
├── 0
│   ├── p
│   └── U
├── constant
│   ├── polyMesh
│   │   ├── boundary
│   │   ├── cellZones
│   │   ├── faces
│   │   ├── faceZones
│   │   ├── neighbour
│   │   ├── owner
│   │   ├── points
│   │   ├── pointZones
│   │   └── sets
│   ├── transportProperties
│   └── turbulenceProperties
└── system
    ├── controlDict
    ├── fvSchemes
    └── fvSolution

4.1.1.6. Execution of the solver

Before execution of the solver, renumbering of the grid or mesh need to be performed as well as checking the quality of the grid. Renumbering the grid-cell 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

cd my-flat-plate/

renumberMesh -overwrite
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 mesh statistics are shown as below

Checking geometry...
  Overall domain bounding box (-0.06 0 0.03) (1.2192 0.15 0.055)
  Mesh (non-empty, non-wedge) directions (1 1 0)
  Mesh (non-empty) directions (1 1 0)
  All edges aligned with or perpendicular to non-empty directions.
  Boundary openness (5.80542e-19 1.1194e-17 1.1403e-14) OK.
  Max cell openness = 2.2093e-16 OK.
  Max aspect ratio = 55.555 OK.
  Minimum face area = 1e-08. Maximum face area = 0.000138887.  Face area magnitudes OK.
  Min volume = 2.5e-10. Max volume = 2.50831e-07.  Total volume = 0.004797.  Cell volumes OK.
  Mesh non-orthogonality Max: 0 average: 0
  Non-orthogonality check OK.
  Face pyramids OK.
  Max skewness = 9.99996e-06 OK.
  Coupled point location match (average 0) OK.

Mesh OK.

End

The mesh non-orthogonality as reported above is 0 and therefore no no-orthogonal corrections are needed in this case. In the case of mesh non-orthogonality 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 my-flat-plate in our case. From the terminal

cd my-flat-plate/

slim | tee my-flat-plate.log

Now the simulation begins and the output of the solution process as seen on the terminal window is saved in the log file, my-flat-plate.log. Once the simulation completes, the log file (my-flat-plate.log) can be processed further to look at the convergence history. Again, in the terminal window

cd my-flat-plate/
foamLog my-flat-plate.log
cd logs/
xmgrace p_0

Now the convergence of pressure can be seen with respect to time. Other convergence properties can also be viewed similarly in the logs/ sub-directory.

4.1.1.7. Results

A brief qualitative data of the simulated flat-plate results are given here. Since the aim here is to obtained the steady solution, the results therefore represent the final steady state condition. In Fig. 69, 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.

_images/fp-velocity-tutorials.png

Figure 69: Contour of velocity magnitude over the flat-plate

4.1.2. Tee Junction (only applicable to Applied CCM subscription customers)

This tutorial introduces the steady, laminar flow through a two-dimensional \(90^\circ\) tee-junction. Here, we will be using Caelus 4.10 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 tee-junction. The details in running a Caelus simulation for the tee-junction will be shown in sufficient detail so that the user is able to repeat the numerical experiment.

4.1.2.1. Objectives

Through the completion of this tutorial, the user will be able to set-up the Caelus simulation for laminar, incompressible flow through a two-dimensional junction and subsequently estimate the flow-split. 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 tee-junction

4.1.2.2. Pre-requisites

It is assumed that the user is familiar with the Linux command line environment using a terminal and Caelus is installed correctly with appropriate environment variables set. The grid used here was generated using Pointwise and the user is free to use their choice of grid generation tool having exporting capabilities to Caelus grid format.

4.1.2.3. Background

The flow in a tee-junction 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 (only applicable to Applied CCM subscription customers).

A schematic of the tee-junction geometry is shown in Fig. 70. 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\).

Figure 70: Tee-junction geometry

In the table, the flow parameters are provided.

Flow conditions
\(Re\) \(V~(m/s)\) \(p~(m^2/s^2)\) \(\nu~(m^2/s)\)
300 1.0 0 Gauge 0.00333

4.1.2.4. Grid Generation

A hexahedral grid is generated for the tee-junction 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 tee-junction geometry and the details are shown in Fig. 71. Due to viscous nature of the flow, the velocity at the walls is zero, which should be represented through a no-slip 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.

Figure 71: Tee junction Computational Domain

The 2D structured grid is shown in Fig. 72 for a \(x-y\) plane. Caelus is a 3D solver and hence requires a 3D grid although the flow here is assumed to be two-dimensional. The 3D grid was obtained by extruding the 2D grid in the third (\(z\) - direction) dimension by one-cell thick. The two \(x-y\) 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.

Figure 72: Tee junction Structured Grid

A coarse grid with a total of 2025 cells is made for the tee-junction 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.

4.1.2.5. Problem definition

We begin with instructions to set-up the tee-junction problem and subsequently configuring the required input files. A full working case can be found in the following location

/tutorials-cases/tee-junction/tee-junction

In the interest of the user, it will be shown here to complete the tutorial assuming no prior set-up data is available.

Directory Structure

First we begin with creating a directory named my-tee-junction. This is the directory where we intend to carry out the simulation. Here all the required set-up is made through working in a terminal window application.

cd tee-junction/
mkdir my-tee-junction/
cd my-tee-junction/

Note

All commands shown here are entered in a terminal window, unless otherwise mentioned

For setting up the problem in my-tee-junction, few more sub-directories are needed, where relevant files are placed. Caelus require a time, constant and system sub-directories. In this case, the time directory would be just 0 as we will begin the simulation at time \(t = 0~s\). These can be obtained by

cd my-tee-junction/
mkdir 0
mkdir constant
mkdir system

ls
0  constant  system

In the 0 sub-directory, two additional files p and U for pressure (\(p\)) and velocity (\(U\)) are required. The input contents of these two files sets the dimensions, initial and boundary conditions to the problem. These three forms the principle entries required. We will now create the two files by

cd my-tee-junction/
touch 0/p
touch 0/U

cd 0/
ls
p  U

The user should note that Caelus is case sensitive and therefore the setting up of the directories, files, and contents should be identical to shown here.

Boundary Conditions and Solver Attributes

Boundary Conditions

Initially, we will be setting up the boundary conditions. Referring to Fig. 71, 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
  • No-slip wall
    • Velocity: Fixed uniform velocity \(u, v, w = 0\)
    • Pressure: Zero gradient
  • Outlet-1
    • Velocity: Zero gradient velocity
    • Pressure: Fixed uniform gauge pressure \(p = 0\)
  • Outlet-2
    • 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 will be done by opening the file p using a text editor and adding the following contents to it.

/*------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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;
	}
	outlet-1
	{
		type			fixedValue;
		value			uniform	0;
	}
	outlet-2
	{
		type			fixedValue;
		value			uniform 0;
	}
	symm-left-right
	{
		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 non-uniform. 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 manner, the input data for the velocity can be added from below

/*------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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=-(maxX-pos().x)*(pos().x-minX)/(0.25*pow(maxX-minX,2))*normal();";
		valueExpression		"Vmax*para";
		value			uniform (0 1 0);
   	}	
	outlet-1
	{
		type			zeroGradient;
	}
	outlet-2
	{
		type			zeroGradient;
	}
	symm-left-right
	{
		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=-(maxX-pos().x)*(pos().x-minX)/(0.25*pow(maxX-minX,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 grid-generation tool. It should be ensured by the user that their names are identically matched. In addition, the two boundaries in \(x-y\) plane obtained due to grid extrusion are named as symm-left-right with applying empty boundary conditions enforcing the flow to be in two-dimensions. 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 tee-junction.

Grid file and Physical Properties

The grid that has been generated for Caelus format has to be placed in the polyMesh sub-directory of constant. First let us create polyMesh as

cd my-tee-junction/
mkdir constant/polyMesh

Now the grid files can be saved in polyMesh sub-directory. Additionally, the physical properties are specified in three different files, placed in the constant sub-directory. These can be generated as follows

touch constant/transportProperties
touch constant/turbulenceProperties

cd constant/
ls
transportProperties  turbulenceProperties

The first file transportProperties, contains the detail about the transport model for the viscosity and kinematic viscosity. The contents are as follows

/*-----------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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 final file in the constant sub-directory to add is the turbulenceProperties. Here, the type of simulation through the keyword simulationType is added with the option laminar; as shown below

/*-----------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*-----------------------------------------------------------------------------*/
FoamFile
{

	version				2.0;
	format				ascii;
	class				dictionary;
	location			"constant";
	object				turbulenceProperties;
}

//-----------------------------------------------------------------------------

simulationType			laminar;

Controls and Solver Attributes

The necessary files to control the simulation and specifying solver attributes such as discretization method, linear solver settings are to be placed in the system directory. These can be created as follows

touch system/controlDict
touch system/fvSchemes
touch system/fvSolution

cd system/
ls
controlDict  fvSchemes  fvSolution

First, let us add the input data to the controlDict file with the following information.

/*-----------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*-----------------------------------------------------------------------------*/

FoamFile
{
	version				2.0;
	format				ascii;
	class				dictionary;
	location			"system";
	object				controlDict;
}

libs (
      "libOpenFOAM.so"
      "libsimpleSwakFunctionObjects.so"
      "libswakFunctionObjects.so"
      "libgroovyBC.so"
     );
			
//-------------------------------------------------------------------------------

application				slim;

startFrom				startTime;

startTime				0;

stopAt					endTime;

endTime					100;

deltaT					0.001;

writeControl				runTime;

writeInterval				1;

purgeWrite				0;

writeFormat				ascii;

writePrecision				6;

writeCompression 			uncompressed;

timeFormat				general;

timePrecision				6;

runTimeModifiable			true;

Since groovyBC is used, few relevant libraries are need and they are imported by calling the following soon after dictionary FoamFile

libs
(
"libOpenFOAM.so"
"libsimpleSwakFunctionObjects.so"
"libswakFunctionObjects.so"
"libgroovyBC.so"
);

in the controlDict file as is shown above. Next, the application slim; refers to the slim 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 time of 100 seconds and is specified through stopAt and endTime. In addition, we need to set a value for time-step since SLIM is a time-accurate solver and a value of 0.001; is used. This is followed by setting the time interval of 1; second 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 4.10
Web:   www.caelus-cml.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	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;
}

interpolationSchemes
{
	default				linear;
	interpolate(HbyA)		linear;
}

snGradschemes
{
	default				corrected;
}

fluxRequired
{
	default				no;
	p				;
}

As apparent from the above file, discretization schemes are set for time-derivative, gradient, divergence and Laplacian terms.

In the final file, fvSolution, linear solver settings are made as given below

/*-----------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*-----------------------------------------------------------------------------*/

FoamFile
{
	version				2.0;
	format				ascii;
	class				dictionary;
	location			"system";
	object				fvSolution;
}

//------------------------------------------------------------------------------

solvers
{
	p
	{
		solver			PCG;
		preconditioner		DIC;
		tolerance		1e-10;
		relTol			0;
	}
	pFinal
	{
		solver			PCG;
		preconditioner		DIC;
		tolerance		1e-10;
		relTol			0;
	}
	U
	{
		solver			PBiCG;
		preconditioner		DILU;
		tolerance		1e-15;
		relTol			0;
	}
}

SLIM
{
	nNonOrthogonalCorrectors	0;
	pRefCell			0;
	pRefValue			0;
}

In the above file, different linear solvers can be seen to be used to solve pressure and velocity fields. Further, the grid used here is perfectly orthogonal and therefore the orthogonal correction specified via nNonOrthogonalCorrectors is set to 0.

With these, the set-up 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

cd my-tee-junction/
tree
.
├── 0
│   ├── p
│   └── U
├── constant
│   ├── polyMesh
│   │   ├── boundary
│   │   ├── cellZones
│   │   ├── faces
│   │   ├── faceZones
│   │   ├── neighbour
│   │   ├── owner
│   │   ├── points
│   │   └── pointZones
│   ├── transportProperties
│   └── turbulenceProperties
└── system
    ├── controlDict
    ├── fvSchemes
    └── fvSolution

4.1.2.6. 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

cd my-tee-junction
renumberMesh -overwrite
checkMesh

During the process of renumbering, grid-cell 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

Checking geometry...
  Overall domain bounding box (0 0 -0.072111) (4 6 0.072111)
  Mesh (non-empty, non-wedge) directions (1 1 0)
  Mesh (non-empty) directions (1 1 0)
  All edges aligned with or perpendicular to non-empty directions.
  Boundary openness (1.07982e-18 -2.28423e-18 1.37054e-18) OK.
  Max cell openness = 1.98307e-16 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 non-orthogonality Max: 0.00224404 average: 0.000335877
  Non-orthogonality check OK.Face pyramids OK.
  Max skewness = 2.64853e-05 OK.
  Coupled point location match (average 0) OK.

Mesh OK.
End

From the above information, the mesh non-orthogonality is very small and therefore no non-orthogonal 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 which in our case is the my-tee-junction. In the terminal,

cd my-tee-junction/

slim | tee my-tee-junction.log

The progress of the simulation can now be seen in the terminal window and is simultaneously written to the log file my-tee-junction.log, which can further be processed to get the convergence history. After the simulation is completed, in the terminal window

foamLog my-tee-junction.log
cd logs/
xmgrace p_0

The plot indicates the convergence history for pressure with respect to time and the same is shown in Fig. 73. Similarly, the user can check the properties of other quantities in the logs/ sub-directory.

_images/tj-convergence-tutorials.png

Figure 73: Convergence of pressure with respect to time

4.1.2.7. Results

The solution obtained for the tee-junction at steady state is shown here using qualitative contour plots. In Fig. 74, 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.

_images/tj-velocitypressure-tutorials.png

Figure 74: Velocity magnitude and pressure contour plots within the tee-junction

4.1.3. Circular Cylinder (only applicable to Applied CCM subscription customers)

The simulation of laminar, incompressible flow over a two-dimensional circular cylinder is shown in this tutorial. Caelus 4.10 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 time-dependent simulation of a externally separated flow. Further to this, the flow around the cylinder would be visualised using velocity and pressure contours.

4.1.3.1. Objectives

Through this tutorial, the user would be familiar in setting up a time-dependent Caelus simulation for laminar, incompressible flow in two-dimensions 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

4.1.3.2. Pre-requisites

It is assumed that the user is familiar with the Linux command line environment using a terminal and Caelus is installed correctly with appropriate environmental 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 Caelus grid format.

4.1.3.3. 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 non-dimensional relating to the frequency of oscillations (\(f\)) of vortex shedding, cylinder diameter (\(D\)) and the flow velocity (\(U\)) as

\[St = \frac{fD}{U}\]

For a Reynolds number based on the cylinder diameter of \(Re_D = 100\), the Strouhal number is about \(St\approx 0.16-0.17\) as determined through experiments and is nearly independent of the diameter of the cylinder. More details are discussed in section Circular Cylinder (only applicable to Applied CCM subscription customers)

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 Fig. 75 shows the schematic of the cylinder in the \(x-y\) plane.

Figure 75: Schematic of the circular cylinder in two-dimensions

In the below table a summary of the freestream conditions are given

Freestream conditions
\(Re_D\) \(U~(m/s)\) \(p~(m^2/s^2)\) \(\nu~(m^2/s)\)
100 1.0 \((0)\) Gauge \(0.02\)

4.1.3.4. Grid Generation

A hexahedral grid around the circular cylinder was development with a O-grid 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 two-dimensional cylinder

A rectangular computational domain in the \(x-y\) plane has been constructed surrounding the circular cylinder as shown in Fig. 76. 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 set-up. 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 no-slip boundary condition as highlighted in blue.

Figure 76: Circular cylinder computational domain

The hexahedral grid around the cylinder is shown in Fig. 77 for a \(x-y\) plane. Caelus is a 3D solver and requires the grid to be in 3D. This was obtained by extruding the grid in the \(x-y\) 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 two-dimensionality 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.

_images/cc-grid-tutorials.png

Figure 77: O-grid around the cylinder and structured gird representation

The 2D domain consisted of 9260 cells in total. An O-grid 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 O-grid block. Upstream of the O-grid 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 O-grid region.

4.1.3.5. Problem definition

We first begin with instructions to set-up the circular cylinder case in addition to the configuration files that are needed. A full working case is found in the following directory:

/tutorials-cases/circular-cylinder/circular-cylinder/

Although a full working example is provided, a step-by-step procedure is shown here assuming no prior user experience

Directory Structure

A directory named my-circular-cylinder will be initially created, with all the necessary data. Follow the instructions given by in a terminal window.

cd circular-cylinder
mkdir my-circular-cylinder
cd my-circular-cylinder/

Note

All commands shown here are entered in a terminal window, unless otherwise mentioned

For setting up the problem in my-circular-cylinder, we need to have few more sub-directories, where relevant files should be present. Caelus require time, constant and system sub-directories within the main my-circular-cylinder working directory. Since we start the simulation at time, t = 0 s, a time sub-directory named 0 is required. Using the terminal,

cd my-circular-cylinder/
mkdir 0
mkdir constant
mkdir system

ls
0  constant  system

In the 0 sub-directory, additional files, p and U representing pressure (\(p\)) and velocity (\(U\)) respectively are needed. The contents of these files sets the dimensions, initialisation and boundary conditions to the case, which also forms three principle entries required. We now start with creating these in the 0 sub-directory by

cd my-circular-cylinder/
touch 0/p
touch 0/U

cd 0/
ls
p  U

The user should note that Caelus is case sensitive and therefore the set-up of directories and case files have to identical to shown here.

Boundary Conditions

We can now begin with setting up the boundary conditions. Referring back to Fig. 76, 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
  • No-slip 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\)), let us open the file p using a text editor and add the following contents to it.

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
	version				2.0;
	format				ascii;
	class				volScalarField;
	object				p;
}

//--------------------------------------------------------------------------------

dimensions				[0 2 -2 0 0 0 0];

internalField				uniform 0;

boundaryField
{
	inlet-1
	{
		type			zeroGradient;
	}
	inlet-2
	{
		type			zeroGradient;
	}
	outlet
	{
		type			fixedValue;
		value			uniform 0;
	}
	symmetry
	{
		type			empty;
	}
	wall
	{
		type			zeroGradient;
	}
}

In 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 non-uniform. 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, let us open the file U and add the following contents to it

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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
{
	inlet-1
	{
		type			fixedValue;
		value			uniform (1.0 0 0);
	}
	inlet-2
	{
		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 \(x-y\) plane obtained due to grid extrusion have been named as symm-left and symm-right with specifying empty boundary conditions forcing Caelus to assume the flow to be in two-dimensions. This completes the setting up of boundary conditions.

Grid file and Physical Properties

The circular cylinder grid file need to be placed in the polyMesh sub-directory which resides in the constant directory. This can be performed as follows

cd my-circular-cylinder/
mkdir constant/polyMesh

Now the grid files can be saved in polyMesh sub-directory. Additionally, the physical properties are specified in different files present in the constant directory. We can now navigate to the constant and create the following files

touch constant/transportProperties
touch constant/turbulenceProperties

cd constant/
ls
polyMesh  transportProperties  turbulenceProperties

The first file is, where transport model and the kinematic viscosity are specified. The contents of this file are as follows

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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 class is the turbulenceProperties file, where the type of simulation is specified as

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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, which need to present in the system directory. First, let us begin with creating these files by navigating to system as

touch system/controlDict
touch system/fvSchemes
touch system/fvSolution

cd system/
ls
controlDict  fvSchemes  fvSolution

Starting with the controlDict file as shown below

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
	version				2.0;
	format				ascii;
	class				dictionary;
	location			"system";
	object				controlDict;
}

//--------------------------------------------------------------------------------

application				slim;

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;
        outputControl   		timeStep;
        outputInterval  		50;
     }
}
//------------------------------------------------------------------------------

In the above file, the application slim refers to the slim solver that is used here. 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 SLIM solver is time-accurate, we also need to set the time-step 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 4.10
Web:   www.caelus-cml.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	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;
}

interpolationSchemes
{
	default				linear;
	interpolate(HbyA)		linear;
}

snGradschemes
{
	default				corrected;
}

fluxRequired
{
	default				no;
	p				;
}

In the fvSolution file, linear solver controls and tolerances are set as shown in the below file

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
	version				2.0;
	format				ascii;
	class				dictionary;
	location			"system";
	object				fvSolution;
}

//--------------------------------------------------------------------------------

solvers
{
	p
	{
		solver			PCG;
		preconditioner		DIC;
		tolerance		1e-6;
		relTol			0.05;
	}
	pFinal
	{
		solver			PCG;
		preconditioner		DIC;
		tolerance		1e-7;
		relTol			0;
	}
	U
	{
		solver			PBiCG;
		preconditioner		DILU;
		tolerance		1e-6;
		relTol			0;
	}
}

SLIM
{
	nNonOrthogonalCorrectors	1;
	pRefCell			0;
	pRefValue			0;
}

Note that different liner solvers are used here to solve for pressure and velocity. Also, nNonOrthogonalCorrectors is set to 1, since there is some degree of non-orthogonality in the grid.

At this stage, the directory structure along with relevant files are completed. Let us verify again by typing the following command and the directory tree should be identical to the one shown below

cd my-circular-cylinder/
tree
.
├── 0
│   ├── p
│   └── U
├── constant
│   ├── polyMesh
│   │   ├── boundary
│   │   ├── faces
│   │   ├── neighbour
│   │   ├── owner
│   │   ├── points
│   │   └── sets
│   ├── transportProperties
│   └── turbulenceProperties
└── system
    ├── controlDict
    ├── fvSchemes
    └── fvSolution

4.1.3.6. Execution of the solver

Prior to solver execution, renumbering of the grid or mesh need to be performed as well as checking the quality of the grid. Renumbering the grid-cell 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

cd my-circular-cylinder/

renumberMesh -overwrite
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

Checking geometry...
  Overall domain bounding box (-10 -10 0) (40 10 0.537713)
  Mesh (non-empty, non-wedge) directions (1 1 0)
  Mesh (non-empty) directions (1 1 0)
  All edges aligned with or perpendicular to non-empty directions.
  Boundary openness (-2.57817e-19 1.67414e-19 -4.29222e-16) OK.
  Max cell openness = 2.19645e-16 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 non-orthogonality Max: 14.6136 average: 1.75565
  Non-orthogonality check OK.
  Face pyramids OK.
  Max skewness = 0.206341 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 my-circular-cylinder in this case.

cd my-circular-cylinder/
slim | tee my-circular-cylinder.log

The starting up of the solution can be seen now and the output of the solution process as seen on the terminal window is saved in the log file, my-circular-cylinder.log. After the completion of the simulation, the log file (my-circular-cylinder.log) can be processed further to look at the convergence history. Again, in the terminal window

cd my-circular-cylinder/
foamLog my-circular-cylinder.log
cd logs/
xmgrace p_0

With the above, the convergence of pressure can be seen with respect to time. The same is shown in the Fig. 78 and due to the periodic nature of vortex shedding, oscillatory convergence of pressure is seen. In a similar manner, other convergence properties can be seen in the logs/ sub-directory.

_images/cc-convergence-tutorials.png

Figure 78: Convergence of pressure with respect to time

4.1.3.7. Results

In this section, some qualitative results obtained as a result of steady vortex shedding in the wake of the cylinder is shown. Fig. 79 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 break-up that occurs in the near wake of the cylinder, travels several diameters downstream eventually diffusing into the flow.

_images/cc-velocitypressure-tutorials.png

Figure 79: Velocity magnitude and pressure contour plots for the flow over the 2D cylinder

4.1.4. Triangular Cavity (only applicable to Applied CCM subscription customers)

In this tutorial, we carry out laminar, incompressible flow inside a triangular cavity in two-dimensions using Caelus 4.10. Details regarding setting up of the directories, fluid properties, boundary conditions, etc will be discussed. Subsequent to this, post-processing the velocity distribution along the center-line will be shown.

4.1.4.1. Objectives

With the completion of this tutorial, the user would be familiar with setting up a steady-state Caelus simulation for laminar, incompressible flow over lip-driven 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

4.1.4.2. Pre-requisites

The user is assumed to be familiar with the Linux command line environment using a terminal and Caelus is installed correctly. The grid used here is generated using Pointwise and the user can use their choice of grid generation tool having exporting capabilities to Caelus grid format.

4.1.4.3. Background

Flow inside lid-driven 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 center-line of cavity. Details regarding the validation of this case is given in Triangular Cavity (only applicable to Applied CCM subscription customers).

The triangular cavity schematic is shown in Fig. 80. 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 two-dimensional plane considered here is in \(x-z\).

Figure 80: Schematic of a triangular cavity

The below table gives the summary of the freestream conditions used here

Freestream conditions
\(Re_D\) \(U~(m/s)\) \(p~(Pa)\) \(\nu~(m^2/s)\)
800 2.0 \((0)\) Gauge \(0.01\)

4.1.4.4. Grid Generation

A hybrid-grid 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 Fig. 81. The velocity at the cavity walls (high lighted in blue) is zero, represented through a no-slip boundary, wherein \(u, v, w = 0\). Whereas the top wall has a uniform velocity in the x-direction.

Figure 81: Computational domain for a triangular cavity

The hybrid grid is shown in Fig. 82. 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 \(x-y\) 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.

_images/tc-grid-tutorials.png

Figure 82: Hybrid grid representation for a triangular cavity

4.1.4.5. Problem definition

This section provides the case set-up procedures along with the configuration files that are needed. A full working case of this problem is given in the following location:

/tutorials-cases/triangular-cavity/triangular-cavity/

The further description gives a step-by-step procedure assuming no prior set-up experience of the user.

Directory Structure

As a first step, a directory named my-triangular-cavity will be created where all the necessary files required for simulation is placed and this is the root or the working directory.

cd triangular-cavity
mkdir my-triangular-cavity
cd my-triangular-cavity/

Note

All commands shown here are entered in a terminal window, unless otherwise mentioned

A few more sub-directories are needed within the main root directory, where relevant files are created. Caelus require time, constant and system sub-directories within the main my-triangular-cavity working directory. Since we start the simulation at time, t = 0 s, a time sub-directory named 0 is required. Using the terminal,

cd my-triangular-cavity/
mkdir 0
mkdir constant
mkdir system

ls
0  constant  system

Within the 0 sub-directory, additional files such as p and U representing pressure and velocity are required. The contents of these files sets the dimensions, initialisation and boundary conditions to the case, which also forms three principle entries required. These can be created using the following commands

cd my-triangular-cavity/
touch 0/p
touch 0/U

cd 0/
ls
p  U

The user should take precautions in setting the directories and files as Caelus is case sensitive. These have to be identical to shown here.

Boundary Conditions

Here, we start with setting-up of the boundary conditions. Referring back to Fig. 81, 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
  • No-slip 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

Starting with the file p for pressure, we can add the following contents using a text editor

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
	version				2.0;
	format				ascii;
	class				volScalarField;
	object				p;
}

//--------------------------------------------------------------------------------

dimensions				[0 2 -2 0 0 0 0];

internalField				uniform 0;

boundaryField
{
	fixed-walls
	{
		type			zeroGradient;
	}
	moving-wall
	{
		type			zeroGradient;
	}
	{
		symm			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 non-uniform. 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 U can be added as shown below

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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
{
	fixed-walls
	{
		type            	fixedValue;
        value				uniform (0 0 0);
	}
	moving-wall
	{
		type            	fixedValue;
        value				uniform (2 0 0);
	}
	{
		symm			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 moving-wall boundary patch. Similarly, the cavity walls (fixed-walls) have uniform (1.0 0 0).

At this stage, the user should ensure that the boundary conditions (fixed-walls, moving-wall 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 \(x-y\) 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 two-dimensions. With this, the setting up of boundary conditions are completed.

Grid file and Physical Properties

The triangular cavity grid should be placed in polyMesh sub-directory which resides in the constant directory. Using the terminal

cd my-triangular-cavity/
mkdir constant/polyMesh

The grid files can now be saved in the polyMesh sub-directory. In addition, the physical properties are specified in different files, all present in the constant directory. By navigating to the constant directory, create the following files

touch constant/transportProperties
touch constant/turbulenceProperties

cd constant/
ls
polyMesh  transportProperties  turbulenceProperties

We start with the file transportProperties, where transport model and the kinematic viscosity are specified. The contents of this file are as follows

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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 4.10
Web:   www.caelus-cml.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. These files would be created in the system directory as shown below

touch system/controlDict
touch system/fvSchemes
touch system/fvSolution

cd system/
ls
controlDict  fvSchemes  fvSolution

We start with the controlDict file and add the following contents

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
	version				2.0;
	format				ascii;
	class				dictionary;
	location			"system";
	object				controlDict;
}

//-------------------------------------------------------------------------------

application				slim;

startFrom				startTime;

startTime				0;

stopAt					endTime;

endTime					200;

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, slim refers to the application, slim 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 SLIM solver is time-accurate, we also need to set the time-step 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 4.10
Web:   www.caelus-cml.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	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;
}

interpolationSchemes
{
	default				linear;
	interpolate(HbyA)		linear;
}

snGradschemes
{
	default				corrected;
}

fluxRequired
{
	default				no;
	p				;
}

In the following file, which is the fvSolution the linear solver controls and tolerances are set as shown below

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
	version				2.0;
	format				ascii;
	class				dictionary;
	location			"system";
	object				fvSolution;
}

//------------------------------------------------------------------------------

solvers
{
	p
	{
		solver			PCG;
		preconditioner		DIC;
		tolerance		1e-10;
		relTol			0;
	}
	pFinal
	{
		solver			PCG;
		preconditioner		DIC;
		tolerance		1e-10;
		relTol			0;
	}
	U
	{
		solver			PBiCG;
		preconditioner		DILU;
		tolerance		1e-15;
		relTol			0;
	}
}

SLIM
{
	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 non-orthogonality present due to triangular nature of the geometry.

This completes the set-up 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

cd my-triangular-cavity/
tree
.
├── 0
│   ├── p
│   └── U
├── constant
│   ├── polyMesh
│   │   ├── boundary
│   │   ├── faces
│   │   ├── neighbour
│   │   ├── owner
│   │   ├── points
│   │   └── sets
│   ├── transportProperties
│   └── turbulenceProperties
└── system
    ├── controlDict
    ├── fvSchemes
    ├── fvSolution

4.1.4.6. Execution of the solver

Renumbering and checking the mesh quality is needed before the solver is executed. Renumbering the grid-cell 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

cd my-triangular-cavity/

renumberMesh -overwrite
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

Checking geometry...
  Overall domain bounding box (-1 -4 0) (1 0 0.0447214)
  Mesh (non-empty, non-wedge) directions (1 1 0)
  Mesh (non-empty) directions (1 1 0)
  All edges aligned with or perpendicular to non-empty directions.
  Boundary openness (-6.25534e-18 -5.30678e-18 -6.4712e-16) OK.
  Max cell openness = 2.16191e-16 OK.
  Max aspect ratio = 2.97779 OK.
  Minimum face area = 8.00174e-05. Maximum face area = 0.00322809.  Face area magnitudes OK.
  Min volume = 3.57849e-06. Max volume = 0.000109372.  Total volume = 0.178886.  Cell volumes OK.
  Mesh non-orthogonality Max: 41.8682 average: 9.34065
  Non-orthogonality check OK.
  Face pyramids OK.
  Max skewness = 0.621726 OK.
  Coupled point location match (average 0) OK.

Mesh OK.

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 my-triangular-cavity in this case.

cd my-triangular-cavity/

slim | tee my-triangular-cavity.log

The solution start-up process can now be seen and the output of the progress is saved in the log file, my-triangular-cavity.log. With the completion of the simulation, the log file (my-triangular-cavity.log) can be processed further to look at the convergence history. Again, in the terminal window

cd my-triangular-cavity/
foamLog my-triangular-cavity.log
cd logs/
xmgrace p_0

The convergence of the pressure can now be seen with respect to time. Similarly other convergence properties can be seen in logs sub-directory

_images/tc-convergence-tutorials.png

Figure 83: Convergence of pressure with respect to time

4.1.4.7. Results

The flow within the cavity is shown here at steady state condition. Fig. 84 presents the velocity and pressure contour plots. In addition, the streamlines indicate the multiple vortices formed within the cavity.

_images/tc-velocitypressure-tutorials.png

Figure 84: Velocity magnitude and pressure contour plots within the triangular cavity

4.1.5. Spherical Cavity (only applicable to Applied CCM subscription customers)

In this tutorial, we look at simulating buoyant flow inside a spherically heated cavity using Caelus 4.10 which will be a three-dimension 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.

4.1.5.1. Objectives

With this tutorial, the user would be familiar in setting up a steady-state 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

4.1.5.2. Pre-requisites

It is assumed that the user is familiar with the Linux command line environment using a terminal and Caelus is installed correctly with appropriate environmental 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 Caelus grid format.

4.1.5.3. 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 (only applicable to Applied CCM subscription customers).

The schematic representation of the sphere is shown in Fig. 85 and only half of the sphere is considered here having the plane of symmetry in #:math:x-y 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\).

_images/sphere-schematic-tutorials.png

Figure 85: Schematic of a sphere

The thermal boundary condition to the spherical wall was generated by specifying temperature (\(T\)) as a function of distance (\(x\)), which are as follows:

\[T = x\]

with this,

\[T = -0.5~K \quad \text{at} \quad x = -0.5~m\]
\[T = 0.5~K \quad \text{at} \quad x = 0.5~m\]

To obtain a smooth distribution of temperature, sufficient cells in the radial direction is required. For buoyancy driven flows, the non-dimensional 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

Freestream conditions
\(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}\)

4.1.5.4. Grid Generation

A hexahedral grid within the sphere was developed via the O-grid 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 \(x-y\) 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 Fig. 86 the computational domain, in particular the applied temperature boundary condition when applied is shown.

_images/sphere-domain-tutorials.png

Figure 86: Computational domain and temperature boundary condition for a spherical cavity

The structured O-grid over the spherical wall and on the plane of symmetry is shown in Fig. 87. 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.

_images/sphere-grid-tutorials.png

Figure 87: O-grid representation distribution on the wall and on the plane of symmetry

4.1.5.5. Problem definition

We begin with instructions to set-up the spherical cavity case and the required configuration files that are needed. A full working case is found in the following directory:

/tutorials-cases/triangular-cavity/spherical-cavity/

The description that follows here give a step-by-step procedure assuming no prior set-up experience of the user

Directory Structure

To begin with, a directory named my-spherical-cavity will be created and all the necessary files required for the simulation would be placed. This is the root or the working directory. The user should note that the set-up is made through working in a terminal window application.

cd spherical-cavity
mkdir my-spherical-cavity
cd my-spherical-cavity/

Note

All commands shown here are entered in a terminal window, unless otherwise mentioned

In order to set-up the problem in my-spherical-cavity, few more sub-directories are required, where relevant files are placed. Caelus require a time, constant and system sub-directories. Since we will be starting the simulation at time \(t = 0~s\), the time directory would be just 0. This can be obtained by

cd my-spherical-cavity/
mkdir 0
mkdir constant
mkdir system

ls
0  constant  system

In the 0 sub-directory, few additional files p, p_rgh, kappat, T, and U for pressure, buoyant pressure, turbulent thermal conductivity, temperature and velocity are needed respectively. Note that even for a laminar simulation, kappat is required, although a value of 0 would be used. The contents of these files sets the dimensions, initialisation and boundary conditions to the case, which also forms three principle entries required. We begin with creating these files in the 0 subdirectory by

cd my-spherical-cavity/
touch 0/p
touch 0/p_rgh
touch 0/kappat
touch 0/T
touch 0/U

cd 0/
ls kappat  p  p_rgh  T  U

The user should be aware that Caelus is case sensitive and therefore the setting up of the directories, files and contents 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 before the solver is executed.

The first quantity to define would be the pressure (\(p\)) and this will be done by opening the file p using a text editor and adding the following contents to it.

/*------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
	version				2.0;
	format				ascii;
	class				volScalarField;
	object				p;
}

//-------------------------------------------------------------------------------

dimensions				[0 2 -2 0 0 0 0];

internalField				uniform 0;

boundaryField
{
	symm
	{
		type			symmetryPlane;
	}
	wall-1
	{
		type			zeroGradient;
	}
	wall-2
	{
		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 non-uniform. 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 can be added from below

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*-------------------------------------------------------------------------------*/
FoamFile
{
	version     			2.0;
	format      			ascii;
	class       			volScalarField;
	object      			p_rgh;
}
//-------------------------------------------------------------------------------/

dimensions      			[0 2 -2 0 0 0 0];

internalField   			uniform 0;

boundaryField
{
	symm
	{
		type            	symmetryPlane;
	}

	wall-1
	{
		type            	buoyantPressure;
		rho             	rhok;
		value           	uniform 0;
	}
	wall-2
	{
		type            	buoyantPressure;
		rho             	rhok;
		value          		uniform 0;
	}
}

In the next, the contents for kappat will be added.

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/
FoamFile
{
	version     			2.0;
	format      			ascii;
	class       			volScalarField;
	location    			"0";
	object     				kappat;
}
// ------------------------------------------------------------------------------/

dimensions      			[0 2 -1 0 0 0 0];

internalField   			uniform 0;

boundaryField
{
	symm
	{
		type            	symmetryPlane;
	}
	wall-1
	{
		type            	fixedValue;
		value           	uniform 0;
	}
	wall-2
	{
		type            	fixedValue;
		value          		uniform 0;
	}
}

The next would be 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. Add the following to the T file

/*------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
	version				2.0;
	format				ascii;
	class				volScalarField;
	object				T;
}

//-------------------------------------------------------------------------------

dimensions				[0 0 0 1 0 0 0];

internalField				uniform 0;

boundaryField
{
	symm
	{
		type			symmetryPlane;
	}
	wall-1
	{
		type			zeroGradient;
	}
	wall-2
	{
		type			zeroGradient;
	}
}
			

The final file in this directory is the velocity and the contents can be added from below

/*------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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
{
	symm
	{
		type			symmetryPlane;
	}
	wall-1
	{
		type			fixedValue;
		value			uniform (0 0 0);
	}
	wall-2
	{
		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, wall-1, wall-2, etc) specified in the above files should be the grid boundary patches (surfaces) that are generated by the grid-generation tool. The user should ensure that their names are identically matched.

Grid file and Physical Properties

Now the spherical cavity grid should be placed in polyMesh sub-directory which resides in the constant directory. Using the terminal

cd my-spherical-cavity/
mkdir constant/polyMesh

Now the grids can be saved in the polyMesh sub-directory. Additionally, physical properties are specified in separate files, placed in the constant sub-directory. These can be generated as follows

touch constant/transportProperties
touch constant/turbulenceProperties

cd constant/
ls
g  polyMesh  transportProperties  turbulenceProperties

The first file g is where 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 4.10
Web:   www.caelus-cml.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, which is the transportProperties, details about the transport model for viscosity, Prandtl number, coefficient of thermal expansion are specified. The contents are as follows

/*------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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.4e-4;

// Thermal expansion coefficient
beta					beta [0 0 0 -1 0 0 0] 3.5677e-5;

// 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 final file in this directory is the turbulenceProperties. Here, the type of simulation is specified through the keyword simulationType with the option laminar; for a laminar simulation as shown below

/*------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{

	version				2.0;
	format				ascii;
	class				dictionary;
	location			"constant";
	object				turbulenceProperties;
}

//--------------------------------------------------------------------------------

simulationType				laminar;

Controls and Solver Attributes

Here we set-up 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 and can be created as follows

touch system/controlDict
touch system/fvSchemes
touch system/fvSolution

cd system/
ls
controlDict  fvSchemes  fvSolution

First, we begin with adding the input data to the controlDict file with the following information

/*------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
	version				2.0;
	format				ascii;
	class				dictionary;
	location			"system";
	object				controlDict;
}

libs (
      "libOpenFOAM.so"
      "libsimpleSwakFunctionObjects.so"
      "libswakFunctionObjects.so"
      "libgroovyBC.so"
     );

//-------------------------------------------------------------------------------

application				boussinesqSLIM;

startFrom				startTime;

startTime				0;

stopAt					endTime;

endTime					2000;

deltaT					1.0;

writeControl				runTime;

writeInterval				10;

purgeWrite				0;

writeFormat				ascii;

writePrecision				12;

writeCompression 			uncompressed;

timeFormat				general;

timePrecision				6;

runTimeModifiable			true;

As we will be using funkySetFields, few relevant libraries are needed and they are imported by calling soon after the dictionary FoamFile

libs
(
"libOpenFOAM.so"
"libsimpleSwakFunctionObjects.so"
"libswakFunctionObjects.so"
"libgroovyBC.so"
);

in the controlDict file as is shown above. The next is the application, boussinesqSLIM which is the buoyancy version of the SLIM 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 time of 100 seconds and is specified through stopAt and endTime. In addition, we need to set a value for time-step since boussinesqSLIM is a time-accurate solver and a value of 1; is used. This is followed by setting the time interval of 10; seconds 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 4.10
Web:   www.caelus-cml.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	linearUpwind grad(U);
	div(phi,T)      		Gauss 	linearUpwind grad(T);
	div((nuEff*dev(T(grad(U)))))	Gauss   linear;
}

laplacianSchemes
{
	default				none;
	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;
}

interpolationSchemes
{
	default				linear;
	interpolate(HbyA)		linear;
}

snGradschemes
{
	default				corrected;
}

fluxRequired
{
	default				no;
	p_rgh;				;
}

In the above file, the discretization schemes are set for time-derivative, gradient, divergence and Laplacian terms.

The final file is the fvSolution, where linear solver settings are provided and is as given below

/*------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
	version				2.0;
	format				ascii;
	class				dictionary;
	location			"system";
	object				fvSolution;
}

//--------------------------------------------------------------------------------

solvers
{
	p_rgh
	{
		solver			PCG;
		preconditioner		DIC;
		tolerance		1e-10;
		relTol			0;
	}

	p_rghFinal
	{
		$p_rgh;
		relTol          	0;
	}

	"(U)"
	{
		solver          	PBiCG;
		preconditioner  	DILU;
		tolerance       	1e-15;
		relTol          	0;
	}

	"(U)Final"
	{
		$U;
		relTol          	0;
	}

}

SLIM
{
	nNonOrthogonalCorrectors	2;
	pRefCell			0;
	pRefValue			0;
}

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 non-orthogonality present in the grid.

With these, the set-up 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

cd my-tee-junction/
tree
.
├── 0
│   ├── kappat
│   ├── p
│   ├── p_rgh
│   ├── T
│   └── U
├── constant
│   ├── g
│   ├── polyMesh
│   │   ├── boundary
│   │   ├── faces
│   │   ├── neighbour
│   │   ├── owner
│   │   ├── points
│   │   └── sets
│   ├── transportProperties
│   └── turbulenceProperties
└── system
    ├── controlDict
    ├── fvSchemes
    └── fvSolution

4.1.5.6. 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

cd my-spherical-cavity

funkySetFields -field T -expression "pos().x" -time 0 -keepPatches -valuePatches "wall-1 wall-2"

The above expression re-writes 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

cd my-spherical-cavity

renumberMesh -overwrite
checkMesh

During the process of renumbering, grid-cell 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

Checking geometry...
  Overall domain bounding box (-0.5 -0.5 0) (0.5 0.5 0.5)
  Mesh (non-empty, non-wedge) directions (1 1 1)
  Mesh (non-empty) directions (1 1 1)
  Boundary openness (2.69541574409e-17 -1.48776322194e-17 -2.09950188835e-16) OK.
  Max cell openness = 2.60400799961e-16 OK.
  Max aspect ratio = 5.13298999541 OK.
  Minimum face area = 7.26421673269e-05. Maximum face area = 0.0145134277912.  Face area magnitudes OK.
  Min volume = 8.32359098625e-07. Max volume = 0.000256180182954.  Total volume = 0.259143943239.  Cell volumes OK.
  Mesh non-orthogonality Max: 64.2010673845 average: 15.1363376112
  Non-orthogonality check OK.
  Face pyramids OK.
  Max skewness = 0.601536462888 OK.
  Coupled point location match (average 0) OK.

Mesh OK.

End

The above information gives a maximum mesh non-orthogonality angle of 64.2 and therefore non-orthogonal 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 which in our case is the my-spherical-cavity. In the terminal,

cd my-spherical-cavity/

slim | tee my-spherical-cavity.log

The progress of the simulation can now be seen in the terminal window and the information shown is simultaneously written to the log file my-spherical-cavity.log, which can be further processed to get the convergence history. Once the simulation is complete, in the terminal window

foamLog my-spherical-cavity.log
cd logs/
xmgrace p_rgh_0

The plot indicates the convergence history for pressure with respect to time and the same is shown in Fig. 88. The convergence of other properties can be similarly found in logs/ sub-directory.

_images/sphere-convergence-tutorials.png

Figure 88: Convergence of pressure with respect to time

4.1.5.7. Results

Here, the solution obtained at steady state is shown and is compared with the analytical solution. In Fig. 89, the comparison of temperature isotherms is presented. The analytical solution is a first order approximation. A close agreement between the two can be observed.

_images/sphere-isotherms-tutorials.png

Figure 89: Comparison of temperature isotherms between computational and analytical solutions.

4.2. Turbulent Flows

4.2.1. Turbulent Flat Plate

This tutorial considers the simulation of turbulent incompressible flow over a two-dimensional sharp leading-edge flat plate using Caelus 4.10. 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.

4.2.1.1. Objectives

Through this tutorials the user will be familiarised with setting up the Caelus simulation for steady, turbulent, incompressible flow over a sharp leading-edge flat-plate in two-dimensions. 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

4.2.1.2. Pre-requisites

It is assumed that the user is familiar with the Linux command line environment using a terminal 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.

4.2.1.3. Background

Turbulent flow over a flat-plate configuration presents an ideal case to introduce the user with the turbulent simulation using Caelus. Here, the steady-state 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 flat-plate 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 Fig. 90, a schematic of the flat-plate is shown. Note that the 2D plane of interest is in \(x-z\) directions.

Figure 90: Schematic of the flat-plate flow

The freestream conditions that would 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 0.3048 \(5 \times 10^6\) 69.436113 Gauge (0) 300 \(1.38872\times10^{-5}\)

4.2.1.4. Grid Generation

The hexahedral grid used in this tutorial is obtained from Turbulence Modeling Resource that has 137 X 97 cells in \(x-z\) 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 flat-plate. In Fig. 91 below, the details of the boundaries in 2D (\(x-z\) 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 leading-edge is at \(x=0\). Due to the viscous nature of the flow, the velocity at the wall is zero which is represented through a no-slip 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 flat-plate the no-slip wall. The entire top boundary will be again modelled as a symmetry plane.

Figure 91: Flat-plate computational domain

The polyMesh grid as noted earlier is in 3D. However, since the flow over a flat-plate is two-dimensional, the 2D plane that is considered here is in \(x-z\) directions. It would therefore be ideal to have one-cell thick in the direction (\(y\)), normal to the 2D plane of interest, where the flow is considered symmetry. The two \(x-z\) 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.

_images/t-fp-grid-tutorials.png

Figure 92: Flat-plate computational grid in \(x-z\) plane

In Fig. 92, the 2D grid is shown which has 137 X 97 cells in \(x-z\) 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 wall-functions would be used to estimate the velocity gradients near the wall and integration is carried up to the wall.

4.2.1.5. Problem definition

In this section, several key instructions would be provided to set-up the turbulent flat-plate problem along with details of file configuration. A full working case can be found in:

/tutorials-cases/turbulent-flat-plate/turbulent-flat-plate

It will however be shown here to complete the tutorial assuming no prior set-up or data is available to the user.

Directory Structure

Initially, let us create a directory named my-turbulent-flat-plate where the necessary simulation would be carried out using a terminal window.

cd turbulent-flat-plate
mkdir my-turbulent-flat-plate
cd my-turbulent-flat-plate/

Note

All commands shown here are entered in a terminal window, unless otherwise mentioned

For setting up the problem in my-turbulent-flat-plate, we need to further have few more sub-directories where relevant files can be created. Caelus requires time, constant and system sub-directories. Since we will be begin the simulation at time \(t = 0~s\), the time sub-directory should be just 0. This can be done by

cd my-turbulent-flat-plate/
mkdir 0
mkdir constant
mkdir system

ls
0  constant  system

In the 0 sub-directory, 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\)) p
Velocity (\(U\)) U
Turbulent viscosity (\(\nu\)) nut
Turbulence field variable (\(\tilde{\nu}\)) nuTilda (Only for SA model)
Turbulent kinetic energy (\(k\)) k (Only for \(k-\omega~\rm{SST}\) model)
Turbulent dissipation rate (\(\omega\)) omega (Only for \(k-\omega~\rm{SST}\) model)

As can be noted from the above table, we will be considering two turbulence models namely, Spalart-Allmaras (SA) and \(k-\omega\) - Shear Stress Transport (\(\rm{SST}\)) in the current exercise. The contents of the files listed above sets the dimensions, initialisation and boundary conditions to the problem, which also forms the three principle entries required. Let us start with creating these files in 0 by

cd my-turbulent-flat-plate/
touch 0/p
touch 0/U
touch 0/nut
touch 0/nuTilda
touch 0/k
touch 0/omega

cd 0/
ls
k  nut  nuTilda  omega  p  U

The user should take into account that Caelus is case sensitive and therefore the directory set-up should be identical to what is shown here.

Boundary Conditions and Solver Attributes

Boundary Conditions

Initially, let us set-up the boundary conditions. Referring back to Fig. 91, the following are the boundary conditions that will be specified:

  • Inlet
  • Symmetry
    • Velocity: Symmetry
    • Pressure: Symmetry
    • Turbulence: Symmetry
  • No-slip 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

Starting with the pressure, let us open p using a text editor and add the following contents to it.

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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. Next we shall explain 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 non-uniform. 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 and add the following to it

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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				fixedValue;
        value				uniform (0 0 0);
    }
}

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 and wall boundary patches have three velocity components.

Similar to p and U, the turbulent properties are also needed at the boundaries and these can be set in the similar process. We begin with opening the file nut, which is the turbulent kinematic viscosity and add the following

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
    	version				2.0;
    	format				ascii;
    	class				volScalarField;
    	location			"0";
    	object				nut;
}
//--------------------------------------------------------------------------------

dimensions      			[0 2 -1 0 0 0 0];

internalField   			uniform 2.92240e-06; // Enable this line for SA Model
internalField   			uniform 1.24985e-07; // Enable this line for k-omega SST Model

boundaryField
{
    bottom
    {
        type				symmetryPlane;
    }
    inflow
    {
        type				fixedValue;
        //value				uniform 2.9224023e-06; Enable this line for SA Model
		//value				uniform 1.24985e-07; // Enable this line for k-omega SST Model
    }
    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 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 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 should be added to the file nuTilda and the entries specified for the boundaryField are identical to that of turbulent kinematic viscosity explained above.

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
    	version				2.0;
    	format				ascii;
    	class				volScalarField;
    	location			"0";
    	object				nuTilda;
}
//--------------------------------------------------------------------------------

dimensions				[0 2 -1 0 0 0 0];

internalField				uniform 4.166166e-05;

boundaryField
{
    bottom
    {
        type				symmetryPlane;
    }
    inflow
    {
        type				fixedValue;
        value				uniform 4.166166e-05;
    }
    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 below contents should be added to the kinetic energy file, k

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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 1e-10;
    }
}

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 model needs to be specified for internalField, inflow and wall. Please note that for wall boundaryField with no wall-function, a small, non-zero fixedValue is required.

Next, we add the following contents to the file omega and the value for this is evaluated as detailed in Turbulence freestream conditions for model.

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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 \(x-z\) 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 flat-plate grid files has to be placed in the polyMesh sub-directory which resides in the constant directory. This can be carried out as follows

cd my-turbulent-flat-plate/
mkdir constant/polyMesh

The grid that has been generated can now be saved in polyMesh sub-directory. Additionally, the physical properties are specified in various different files present in the directory constant. Let us navigate to constant and create the following files

touch constant/RASProperties
touch constant/transportProperties
touch constant/turbulenceProperties

cd constant/
ls
polyMesh  RASProperties  transportProperties  turbulenceProperties

As you can see in the above terminal output, the three files are listed in addition to the polyMesh sub-directory. In the first file, RASProperties, the Reynolds-Average-Stress (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 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
	version				2.0;
	format				ascii;
	class				dictionary;
	location			"constant";
	object				RASProperties;
}

//--------------------------------------------------------------------------------

// For Spalarat-Alamaras Model, enable the below line
// RASModel				SpalartAllmarasCurvature;

// For k-omega SST Model, enable the below line
// RASModel				koSST;

turbulence				on;

printCoeffs				on;

Next, we look at the transportProperties file, where transport model and kinematic viscosity is specified.

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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.388722E-5;

As the viscous behaviour is Newtonian, the transportModel is Newtonian keyword 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 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
	version				2.0;
	format				ascii;
	class				dictionary;
	location			"constant";
	object				turbulenceProperties;
}

//--------------------------------------------------------------------------------

simulationType				RASModel;

Controls and Solver Attributes

In this section, the files require 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. Navigating to the system directory, we will create the following files

touch system/controlDict
touch system/fvSchemes
touch system/fvSolution

cd system/
ls
controlDict  fvSchemes  fvSolution

First, we begin with adding the input data to the controlDict file as below

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
	version				2.0;
	format				ascii;
	class				dictionary;
	location			"system";
	object				controlDict;
}

// Enable the below line for Spalart-Almaras turbulence Model
// libs ( "libincompressibleSpalartAlmarasCurvature.so" );

// Enable the below line for SST turbulence Model
// libs ("libincompressibleKOSST.so");
			
//-------------------------------------------------------------------------------

application				simpleFOAM;

startFrom				startTime;

startTime				0;

stopAt					endTime;

endTime					20000;

deltaT					1;

writeControl				runTime;

writeInterval				1000;

purgeWrite				0;

writeFormat				ascii;

writePrecision				12;

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.347049;
        outputControl   		timeStep;
        outputInterval  		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 steady-state 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 should be set through the fvSchemes file show below and the contents should be copied

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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 S-A model only
	div(phi,k)			Gauss 	upwind; // will be used for k-epsilon & k-omega only
	div(phi,omega)			Gauss 	upwind;	// Will be used for k-omega model only
	div((nuEff*dev(T(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 S-A model only
	laplacian(DkEff,k)		Gauss 	linear 	corrected; // Will be used for k-omega & k-omega only
	laplacian(DomegaEff,omega)	Gauss 	linear 	corrected; // Will be used for k-omega model only
	laplacian((1|A(U)),p)		Gauss 	linear	corrected;
	laplacian(1,p)			Gauss 	linear 	corrected;
}

interpolationSchemes
{
	default				linear;
	interpolate(HbyA)		linear;
}

snGradschemes
{
	default				corrected;
}

fluxRequired
{
	default				no;
	p				;
}

The linear solver controls and tolerances are set in fvSolution as given below

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/
FoamFile
{
	version				2.0;
	format				ascii;
	class				dictionary;
	location			"system";
	object				fvSolution;
}

//------------------------------------------------------------------------------

solvers
{
	p
	{
		solver			PCG;
		preconditioner		DIC;
		tolerance		1e-8;
		relTol			0.01;
	}
	U
	{
		solver			PBiCG;
		preconditioner		DILU;
		tolerance		1e-7;
		relTol			0.01;
	}

	"(k|omega|nuTilda)"
	{
		solver          	PBiCG;
		preconditioner  	DILU;
		tolerance       	1e-08;
		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. Of further note, the relTol is not set to 0 unlike a time-accurate set-up. This is because we are solving for a steady-state 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 set-up of the directory structure with all the relevant files are completed. We can verify again by typing the following command and the directory tree should appear identical to the one shown below

cd my-turbulent-flat-plate/
tree
  .
  ├── 0
  │   ├── k
  │   ├── nut
  │   ├── nuTilda
  │   ├── omega
  │   ├── p
  │   └── U
  ├── constant
  │   ├── polyMesh
  │   │   ├── boundary
  │   │   ├── cellZones
  │   │   ├── faces
  │   │   ├── faceZones
  │   │   ├── neighbour
  │   │   ├── owner
  │   │   ├── points
  │   │   ├── pointZones
  │   │   └── sets
  │   ├── RASProperties
  │   ├── transportProperties
  │   └── turbulenceProperties
  └── system
        ├── controlDict
        ├── fvSchemes
        └── fvSolution

4.2.1.6. 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 grid-cell 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.

cd my-turbulent-flat-plate/

renumberMesh -overwrite
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

Checking geometry...
  Overall domain bounding box (-0.33333 -1 0) (2 0 1)
  Mesh (non-empty, non-wedge) directions (1 0 1)
  Mesh (non-empty) directions (1 0 1)
  All edges aligned with or perpendicular to non-empty directions.
  Boundary openness (-3.91843881449e-17 2.55923035071e-16 9.79609703622e-18) OK.
  ***High aspect ratio cells found, Max aspect ratio: 21429.3666252, number of cells 2755
  <<Writing 2755 cells with high aspect ratio to set highAspectRatioCells
  Minimum face area = 8.03897309804e-09. Maximum face area = 0.114850939537.  Face area magnitudes OK.
  Min volume = 8.03897309804e-09. Max volume = 0.00493382049012.  Total volume = 2.33333.  Cell volumes OK.
  Mesh non-orthogonality Max: 6.66818534393e-06 average: 0
  Non-orthogonality check OK.
  Face pyramids OK.
  Max skewness = 3.10964734959e-07 OK.
  Coupled point location match (average 0) OK.

Failed 1 mesh checks.

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 should 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 which is my-turbulent-flat-plate in our case. From the terminal,

cd my-turbulent-flat-plate/

simpleSolver | tee my-turbulent-flat-plate.log

With the execution of the above command, the simulation begins and the progress of the solution can be seen in the terminal window and simultaneously it is written to the specified log file (my-turbulent-flat-plate.log). The log file can be further processed to look at the convergence history and this can be done as follows once the simulation is complete

cd my-turbulent-flat-plate/
foamLog my-turbulent-flat-plate.log
cd logs/
xmgrace p_0

The xmgrace p_0 allows you to look at the convergence of pressure with respect to the number of iterations carried out. In Fig. 93 the same is shown. Other convergence variables can also be accessed similarly in the logs/ sub-directory.

_images/t-fp-convergence-tutorials.png

Figure 93: Convergence of pressure with respect to iterations

4.2.1.7. Results

The turbulent flow over the flat plate is shown here through velocity magnitude contours for SA model. In Fig. 94 the boundary layer over the entire flat-plate 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.

_images/t-fp-velocity-tutorials.png

Figure 94: Contour of velocity magnitude over the flat-plate

4.2.2. Bump in a Channel

The simulation of turbulent flow over a two-dimensional bump in a channel is considered in this tutorial and will be performed using Caelus 4.10. 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.

4.2.2.1. 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 two-dimensional bump in a channel and be able to post-process 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

4.2.2.2. Pre-requisites

The user should be familiar with a Linux command line environment via a terminal. 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.

4.2.2.3. Background

Turbulent flow over a bump in a channel is quite similar to a flat-plate flow, except that due to the curvature effect, a pressure gradient is developed. The flow would be assumed to be steady-state 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 Two-dimensional Bump in a Channel.

The bump, as shown in the schematic below in Fig. 95 has a upstream and a downstream flat-plate 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}\). To match the required Reynolds number, a freestream velocity of U = 69.436 m/s would be used. Note that the two-dimensional plane considered here is in \(x-z\) directions.

Figure 95: Schematic of the 2D bump

The freestream conditions that will be used is given in the below table (Freestream conditions)

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}\)

4.2.2.4. Grid Generation

The hexahedral grid used in this tutorial is obtained from Turbulence Modeling Resource that contains 704 X 320 cells in \(x-z\) 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. Fig. 96 shown below identifies the boundary conditions in a two-dimensional plane (\(x-z\)). 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 no-slip 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.

Figure 96: Computational domain of a 2D bump

The polyMesh grid obtained from the conversion of Plot3D is in 3D. However, the flow over a bump is two-dimensional and is solved in a 2D plane with \(x-z\) directions. Therefore, ideally we can have cells with one-cell thick in the direction (\(y\)), normal to the 2D plane, where the flow can be assumed to be symmetry. The two \(x-z\) 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.

_images/t-bump-grid-tutorials.png

Figure 97: Computational grid of a two-dimensional bump in \(x-z\) plane

In Fig. 97 above, the 2D grid is shown with a distribution of 704 X 320 in \(x-z\) 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 wall-function has been used and the integration is carried out up to the wall.

4.2.2.5. Problem definition

This section deals with several key instructions need to set-up the turbulent flow over a bump. A full working case of this can be found in:

/tutorials-cases/turbulent-bump/turbulent-bump

The guide explained here will show to complete this tutorial assuming no prior set-up or data is available to the user.

Directory Structure

The first step to begin this tutorial is to generate a directory named my-turbulent-bump using a terminal, where the necessary simulation would be carried out.

cd turbulent-bump
mkdir my-turbulent-bump
cd my-turbulent-bump/

Note

All commands shown here are entered in a terminal window, unless otherwise mentioned

In order to set-up the problem in my-turbulent-bump, we need to further create few more sub-directories where relevant files will be created. Caelus requires time, constant and system sub-directories within the main my-turbulent-bump working directory. Typically, the simulations are started at time \(t = 0~s\), which requires a time sub-directory to be 0. Using a terminal these can be done by

cd my-turbulent-bump/
mkdir 0
mkdir constant
mkdir system

ls
0  constant  system

Within the 0 sub-directory, additional files are required for specifying the boundary properties. The below table lists the necessary files required based on the turbulence model chosen

Parameter File name
Pressure (\(p\)) p
Velocity (\(U\)) U
Turbulent viscosity (\(\nu\)) nut
Turbulence field variable (\(\tilde{\nu}\)) nuTilda (Only for SA model)
Turbulent kinetic energy (\(k\)) k (Only for \(k-\omega~\rm{SST}\) model)
Turbulent dissipation rate (\(\omega\)) omega (Only for \(k-\omega~\rm{SST}\) model)

In this tutorial, we will be considering two turbulence models namely, Spalart-Allmaras (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 creating these in the 0 sub-directory by

cd my-turbulent-bump/
touch 0/p
touch 0/U
touch 0/nut
touch 0/nuTilda
touch 0/k
touch 0/omega

cd 0/
ls
k  nut  nuTilda  omega  p  U

The user should take into account that Caelus is case sensitive and therefore the directory set-up should be identical to what is shown here.

Boundary Conditions and Solver Attributes

Boundary Conditions

Referring back to Fig. 96, the following are the boundary conditions that will be specified:

  • Inlet
  • Symmetry
    • Velocity: Symmetry
    • Pressure: Symmetry
    • Turbulence: Symmetry
  • No-slip 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

We begin with opening the p, the pressure file using a text editor and add the following contents to it

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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 non-uniform. 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 can open the file U and add the following to it

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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			fixedValue;
		value			uniform (0 0 0);
	}
}

As with the pressure, 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 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 and wall boundary patches have 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 opening the file nut, which corresponds to turbulent kinematic viscosity and add the following

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
    	version				2.0;
    	format				ascii;
    	class				volScalarField;
    	location			"0";
    	object				nut;
}
//--------------------------------------------------------------------------------

dimensions      			[0 2 -1 0 0 0 0];

//internalField   			uniform 4.87067e-06; // Enable this line for SA Model
//internalField   			uniform 2.08310e-07; // Enable this line for k-omega SST Model

boundaryField
{
	bottom
	{
	type				symmetryPlane;
	}
	inflow
	{
		type			fixedValue;
		//value			uniform 4.87067e-06; // Enable this line for SA Model
		//value			uniform 2.08310e-07; // Enable this line for k-omega SST Model
	}
	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 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 should be added to the file nuTilda and the entries specified for the boundaryField are identical to that of turbulent kinematic viscosity explained above.

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
    	version				2.0;
    	format				ascii;
    	class				volScalarField;
    	location			"0";
    	object				nuTilda;
}
//--------------------------------------------------------------------------------

dimensions				[0 2 -1 0 0 0 0];

internalField				uniform 6.943611e-05;

boundaryField
{
	bottom
	{
		type			symmetryPlane;
	}
	inflow
	{
		type			fixedValue;
		value			uniform 6.943611e-05;
	}
	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. At first, the below contents should be added to the kinetic energy file, k

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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 1e-10;
	}
}

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 model) needs to be specified for internalField, inflow and wall. Please note that for wall boundaryField with no wall-function, a small, non-zero fixedValue is required.

We now add the following contents to the file omega and the value for this is evaluated as detailed in Turbulence freestream conditions for model

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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 \(x-z\) plane named here as left and right have empty boundary conditions which forces Caelus to assume the flow to be in two-dimensions. With this, the setting up of the boundary conditions are complete.

Grid file and Physical Properties

The grid file for the turbulent-bump need to be placed in polyMesh sub-directory which resides in the constant directory. This can be carried out as follows:

cd my-turbulent-bump/
mkdir constant/polyMesh

The generated grid can now be saved in polyMesh sub-directory. In addition to this, the physical properties are specified in various different files present in the constant directory. Let us begin creating these files. First, navigate to the constant directory

touch constant/RASProperties
touch constant/transportProperties
touch constant/turbulenceProperties

cd constant/
ls
polyMesh  RASProperties  transportProperties  turbulenceProperties

As seen above, the three files are listed in addition to the polyMesh sub-directory. The first one, RASProperties in which the Reynolds-Average-Stress (RAS) 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 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
	version				2.0;
	format				ascii;
	class				dictionary;
	location			"constant";
	object				RASProperties;
}

//--------------------------------------------------------------------------------

// For Spalarat-Alamaras Model, enable the below line
// RASModel				SpalartAllmarasCurvature;

// For k-omega SST Model, enable the below line
// RASModel				koSST;

turbulence				on;

printCoeffs				on;

Second from the list is the transportProperties file, where kinematic viscosity is specified as shown below

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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.314537E-5;

Since we will be modelling the viscous behaviour as Newtonian, the transportModel is Newtonian keyword 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 4.10
Web:   www.caelus-cml.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 to be placed in the system directory. Let us begin with navigating to the system directory and create the following files

touch system/controlDict
touch system/fvSchemes
touch system/fvSolution

cd system/
ls
controlDict  fvSchemes  fvSolution

In the controlDict file, add the following details

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
	version				2.0;
	format				ascii;
	class				dictionary;
	location			"system";
	object				controlDict;
}

// Enable the below line for Spalart-Almaras turbulence Model
// libs ( "libincompressibleSpalartAlmarasCurvature.so" );

// Enable the below line for SST turbulence Model
// libs ("libincompressibleKOSST.so");
			
//-------------------------------------------------------------------------------

application				simpleFOAM;

startFrom				startTime;

startTime				0;

stopAt					endTime;

endTime					20000;

deltaT					1;

writeControl				runTime;

writeInterval				1000;

purgeWrite				0;

writeFormat				ascii;

writePrecision				12;

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          		0.80822;
        outputControl   		timeStep;
        outputInterval  		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 steady-state 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 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 and the contents should be copied

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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 S-A model only
	div(phi,k)			Gauss 	upwind; // will be used for k-epsilon & k-omega only
	div(phi,omega)			Gauss 	upwind;	// Will be used for k-omega model only
	div((nuEff*dev(T(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 S-A model only
	laplacian(DkEff,k)		Gauss 	linear 	corrected; // Will be used for k-omega & k-omega only
	laplacian(DomegaEff,omega)	Gauss 	linear 	corrected; // Will be used for k-omega model only
	laplacian((1|A(U)),p)		Gauss 	linear	corrected;
	laplacian(1,p)			Gauss 	linear 	corrected;
}

interpolationSchemes
{
	default				linear;
	interpolate(HbyA)		linear;
}

snGradschemes
{
	default				corrected;
}

fluxRequired
{
	default				no;
	p				;
}

The linear solver controls and tolerances are set in fvSolution as given below

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/
FoamFile
{
	version				2.0;
	format				ascii;
	class				dictionary;
	location			"system";
	object				fvSolution;
}

//------------------------------------------------------------------------------

solvers
{
	p
	{
		solver			PCG;
		preconditioner		DIC;
		tolerance		1e-8;
		relTol			0.01;
	}
	U
	{
		solver			PBiCG;
		preconditioner		DILU;
		tolerance		1e-7;
		relTol			0.01;
	}

	"(k|omega|nuTilda)"
	{
		solver          	PBiCG;
		preconditioner  	DILU;
		tolerance       	1e-08;
		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 not set to 0 unlike a time-accurate set-up as we are solving for a steady-state 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 set-up of the directory structure with all the relevant files are complete. We should verify again by issuing the following command and the directory tree should appear identical to the one shown below

cd my-turbulent-bump/
tree
.
├── 0
│   ├── k
│   ├── nut
│   ├── nuTilda
│   ├── omega
│   ├── p
│   └── U
├── constant
│   ├── polyMesh
│   │   ├── boundary
│   │   ├── cellZones
│   │   ├── faces
│   │   ├── faceZones
│   │   ├── neighbour
│   │   ├── owner
│   │   ├── points
│   │   ├── pointZones
│   │   └── sets
│   ├── RASProperties
│   ├── transportProperties
│   └── turbulenceProperties
└── system
    ├── controlDict
    ├── fvSchemes
    └── fvSolution

4.2.2.6. 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

cd my-turbulent-bump/

renumberMesh -overwrite
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

Checking geometry...
  Overall domain bounding box (-25 -1 0) (26.5 0 5)
  Mesh (non-empty, non-wedge) directions (1 0 1)
  Mesh (non-empty) directions (1 0 1)
  All edges aligned with or perpendicular to non-empty directions.
  Boundary openness (7.07176509381e-19 2.67164792684e-15 1.03247823352e-16) OK.
  ***High aspect ratio cells found, Max aspect ratio: 5347.84086655, number of cells 17714
  <<Writing 17714 cells with high aspect ratio to set highAspectRatioCells
  Minimum face area = 1.60659168306e-09. Maximum face area = 0.901907009882.  Face area magnitudes OK.
  Min volume = 1.60659168306e-09. Max volume = 0.0877839160854.  Total volume = 257.483125413.  Cell volumes OK.
  Mesh non-orthogonality Max: 12.7768609922 average: 3.29301890972
  Non-orthogonality check OK.
  Face pyramids OK.
  Max skewness = 0.132470167186 OK.
  Coupled point location match (average 0) OK.

Failed 1 mesh checks.

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 multi-core 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 sub-directory as shown below.

cd my-turbulent-bump/system

touch decomposeParDict

cd system/
ls
controlDict  decomposeParDict  fvSchemes  fvSolution

Now, the following contents should be added to the decomposeParDict file

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
    version     		2.0;
    format      		ascii;
    class       		dictionary;
    object      		decomposeParDict;
}
//--------------------------------------------------------------------------------

numberOfSubdomains 		8;

method          		simple;

simpleCoeffs
{
	n			(8 1 1);
	delta 			0.001;
}

distributed     		no;

roots
(
);

In the above file, the the keyword numberOfSubdomains defines the number of decomposed sub-domains needed and 8 is used which partitions the grid into 8 sub-domains. 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 (8 1 1) performs 8 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

cd my-turbulent-bump/

decomposePar

Now the decomposition should begin and the details of which are displayed in the terminal window. Subsequently, 8 processor directories will be generated as shown below

cd my-turbulent-bump/
ls
0  constant  processor0  processor1  processor2  processor3  processor4  processor5  processor6  processor7  system

The solver can now be executed for parallel computation from the top directory using the following command

cd my-turbulent-bump/
mpirun -np 8 simpleSolver -parallel | tee my-turbulent-bump.log

Note that here it is assumed that the parallel computing is available in the host machine. If the parallel computing is carried out on other hosts, the command is as follows

cd my-turbulent-bump/
mpirun -np 8 -machinefile machines simpleSolver -parallel | tee my-turbulent-bump.log

where, machines is a plain text file located in the top directory with the list of available hosts names. With the execution of the above commands, the simulation begins and the progress of the solution can be seen in the terminal window and simultaneously it is written to the specified log file (my-turbulent-bump.log).

The log file can be further processed to look at the convergence history and this can be done as follows once the simulation is complete

cd my-turbulent-bump/
foamLog my-turbulent-bump.log
cd logs/
xmgrace p_0

The xmgrace p_0 allows you to look at the convergence of pressure with respect to the number of iterations carried out and the Fig. 98 indicates the same. Other convergence variables can also be accessed similarly in the logs/ sub-directory.

_images/t-bump-convergence-tutorials.png

Figure 98: Convergence of pressure with respect to iterations

4.2.2.7. Results

The flow visualisation over the bump is shown here through the contours of velocity and pressure for SA model. In Fig. 99 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.

_images/t-bump-velocitypressure-tutorials.png

Figure 99: Contours of velocity and pressure over the bump surface

4.2.3. NACA Airfoil

In this tutorial, the turbulent flow over a two-dimensional NACA 0012 airfoil at two angles of attack, namely \(0^\circ\) and \(10^\circ\) will be considered. Caelus 4.10 will be used and the basic steps to set-up 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.

4.2.3.1. Objectives

The user will get familiar in setting up Caelus simulation for steady, turbulent, incompressible flow over a two-dimensional 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

4.2.3.2. Pre-requisites

It is understood that the user will be familiar with the Linux command line environment via a terminal 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.

4.2.3.3. 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 skin-friction distributions. The user can refer to the verification and validation of this case at Two-dimensional NACA 0012 Airfoil.

The schematic of NACA 0012 airfoil at two angles of attack are shown in Fig. 100 for a two-dimensional 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 two-dimensionality is in \(x-z\) directions.

Figure 100: Schematic representation of the airfoil

The freestream conditions are given in the below table

Freestream conditions
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 free-stream 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

Velocity components in x and z directions
\(\alpha~\rm{Degrees}\) \(u~(m/s)\) \(w~(m/s)\)
\(0^\circ\) 52.0770 0.0
\(10^\circ\) 51.2858 9.04307

4.2.3.4. 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 Fig. 101 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 is places 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 no-slip boundary.

Figure 101: Computational domain of a 2D airfoil

The polyMesh grid is in three-dimensions, however the flow over airfoils can be assumed to be 2D at low angles of attack and is solved here in \(x-z\) directions. Therefore, a one-cell thick grid normal to the 2D flow plane is sufficient, where the flow can be assumed to be symmetry. The two \(x-z\) 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 third-dimension 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.

_images/t-airfoil-grid-tutorials.png

Figure 102: Computational grid of a 2D airfoil in \(x-z\) plane

The 2D airfoil grid in \(x-z\) plane is shown in Fig. 102 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 wall-function is used.

4.2.3.5. Problem definition

In this section, various steps needed to set-up the turbulent flow over an airfoil will be shown. A full working case of this can be found in:

/tutorials-cases/turbulent-airfoil/turbulent-airfoil

The following sections will provide details to complete this tutorial assuming no prior set-up or data is available to the user.

Directory Structure

As a first step, a top level working directory has to be created where all the necessary simulation will be carried out. For example, let us name this directory as my-turbulent-airfoil.

cd turbulent-airfoil
mkdir my-turbulent-airfoil
cd my-turbulent-airfoil/

Note

All commands shown here are entered in a terminal window, unless otherwise mentioned

To set-up the problem in my-turbulent-airfoil, a few more sub-directories are needed where relevant files will be created. Caelus requires time, constant and system sub-directories within the main my-turbulent-airfoil working directory. Here, the simulation will be started at time t = 0 s, which requires a time sub-directory named 0. With the use of a terminal window,

cd my-turbulent-airfoil/
mkdir 0
mkdir constant
mkdir system

ls
0  constant  system

The 0 sub-directory requires additional 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\)) p
Velocity (\(U\)) U
Turbulent viscosity (\(\nu\)) nut
Turbulence field variable (\(\tilde{\nu}\)) nuTilda (Only for SA model)
Turbulent kinetic energy (\(k\)) k (Only for \(k-\omega~\rm{SST}\) model)
Turbulent dissipation rate (\(\omega\)) omega (Only for \(k-\omega~\rm{SST}\) model)

We will consider two turbulence models in this tutorial, namely Spalart-Allmaras (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 creating these in the 0 sub-directory by

cd my-turbulent-airfoil/
touch 0/p
touch 0/U
touch 0/nut
touch 0/nuTilda
touch 0/k
touch 0/omega

cd 0/
ls
k  nut  nuTilda  omega  p  U

The user should note that Caelus is case sensitive and therefore the directory and file set-up should be identical to what is shown here.

Boundary Conditions and Solver Attributes

Boundary Conditions

Referring back to Fig. 101, 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:
  • No-slip 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:

First, we start with the pressure and using a text editor add the following contents to the file named p

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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 non-uniform. 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, we can open the file U and add the following to it

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
    	version				2.0;
    	format				ascii;
    	class				volVectorField;
    	location			"0";
    	object				U;
}
//--------------------------------------------------------------------------------

dimensions      			[0 1 -1 0 0 0 0];

// Angle of Attack (AOA)
// Enable the below line for AOA = 0 Degree
// internalField   			uniform (52.077 0 0);

// Enable the below line for AOA = 10 Degrees
// internalField			uniform (51.2858 0 9.04307);


boundaryField
{
	inlet
	{
		type			fixedValue;
		value			uniform (52.077 0 0); // Enable this line for AOA = 0 Degree
		value			uniform (51.2858 0 9.04307); // Enable this line for AOA = 10 Degrees
	}
	left
	{
		type			empty;
	}
	outlet
	{
		type			zeroGradient;
	}
	right
	{
		type			empty;
	}
		wall
	{
		type			fixedValue;
		value			uniform (0 0 0);
	}
}

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]). 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. First, we start with opening the file nut, which is turbulent kinematic viscosity and add the following

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
    	version				2.0;
    	format				ascii;
    	class				volScalarField;
    	location			"0";
    	object				nut;
}
//--------------------------------------------------------------------------------

dimensions      			[0 2 -1 0 0 0 0];

//internalField   			uniform 1.8265e-06; // Enable this line for SA Model
//internalField   			uniform 7.8115e-08; // Enable this line for k-omega SST Model

boundaryField
{
	inlet
	{
		type				fixedValue;
		//value				uniform 1.8265e-06; // Enable this line for SA Model
		//value				uniform 7.8115e-08; // Enable this line for k-omega SST Model
	}
	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 contents should be added to the file nuTilda and the entries specified for the boundaryField are identical to that of turbulent kinematic viscosity explained above.

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/
FoamFile
{
    	version				2.0;
    	format				ascii;
    	class				volScalarField;
    	location			"0";
    	object				nuTilda;
}
//--------------------------------------------------------------------------------

dimensions				[0 2 -1 0 0 0 0];

internalField				uniform 2.60385e-05;

boundaryField
{
	inlet
	{
		type				fixedValue;
		value				uniform 2.60385e-05;
	}
	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, add the below contents to it

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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 1e-10;
	}
}

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 wall-function, a small, non-zero 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 4.10
Web:   www.caelus-cml.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 \(x-z\) plane named here as left and right have empty boundary conditions which forces Caelus to assume the flow to be in two-dimensions. 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 sub-directory, 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. The polyMesh directory can be created as follows

cd my-turbulent-airfoil/
mkdir constant/polyMesh

The grid that has been generated can now be saved in polyMesh sub-directory. In addition, the physical properties are specified in various different files present in the constant directory. Let us begin creating these files. First, navigate to the constant directory

touch constant/RASProperties
touch constant/transportProperties
touch constant/turbulenceProperties

cd constant/
ls
polyMesh  RASProperties  transportProperties  turbulenceProperties

As seen above, the three files are listed in addition to the polyMesh sub-directory. The first one, RASProperties in which the Reynolds-Average-Stress (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 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
	version				2.0;
	format				ascii;
	class				dictionary;
	location			"constant";
	object				RASProperties;
}

//--------------------------------------------------------------------------------

// For Spalarat-Alamaras Model, enable the below line
// RASModel				SpalartAllmarasCurvature;

// For k-omega SST Model, enable the below line
// RASModel				koSST;

turbulence				on;

printCoeffs				on;

Second from the list is the transportProperties file, where kinematic viscosity is specified as shown below

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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.679514E-6;

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 4.10
Web:   www.caelus-cml.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 and create the following files

touch system/controlDict
touch system/fvSchemes
touch system/fvSolution

cd system/
ls
controlDict  fvSchemes  fvSolution

In the controlDict file, add the following details

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
	version				2.0;
	format				ascii;
	class				dictionary;
	location			"system";
	object				controlDict;
}

// Enable the below line for Spalart-Almaras turbulence Model
// libs ( "libincompressibleSpalartAlmarasCurvature.so" );

// Enable the below line for SST turbulence Model
// libs ("libincompressibleKOSST.so");
			
//-------------------------------------------------------------------------------

application				simpleFOAM;

startFrom				startTime;

startTime				0;

stopAt					endTime;

endTime					50000;

deltaT					1;

writeControl				runTime;

writeInterval				1000;

purgeWrite				0;

writeFormat				ascii;

writePrecision				12;

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          		2.15527;
        outputControl   		timeStep;
        outputInterval  		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 steady-state 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 show below and the contents should be copied

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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 S-A model only
	div(phi,k)			Gauss 	upwind; // will be used for k-epsilon & k-omega only
	div(phi,omega)			Gauss 	upwind;	// Will be used for k-omega model only
	div((nuEff*dev(T(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 S-A model only
	laplacian(DkEff,k)		Gauss 	linear 	corrected; // Will be used for k-omega & k-omega only
	laplacian(DomegaEff,omega)	Gauss 	linear 	corrected; // Will be used for k-omega model only
	laplacian((1|A(U)),p)		Gauss 	linear	corrected;
	laplacian(1,p)			Gauss 	linear 	corrected;
}

interpolationSchemes
{
	default				linear;
	interpolate(HbyA)		linear;
}

snGradschemes
{
	default				corrected;
}

fluxRequired
{
	default				no;
	p				;
}

The linear solver controls and tolerances are set in fvSolution as given below

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
	version				2.0;
	format				ascii;
	class				dictionary;
	location			"system";
	object				fvSolution;
}

//------------------------------------------------------------------------------

solvers
{
	p
	{
		solver			PCG;
		preconditioner		DIC;
		tolerance		1e-8;
		relTol			0.01;
	}
	U
	{
		solver			PBiCG;
		preconditioner		DILU;
		tolerance		1e-7;
		relTol			0.01;
	}

	"(k|omega|nuTilda)"
	{
		solver          	PBiCG;
		preconditioner  	DILU;
		tolerance       	1e-08;
		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 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 time-accurate set-up as we are solving for a steady-state 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 set-up 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 my-turbulent-airfoil/
tree
.
├── 0
│   ├── k
│   ├── nut
│   ├── nuTilda
│   ├── omega
│   ├── p
│   └── U
├── constant
│   ├── polyMesh
│   │   ├── boundary
│   │   ├── cellZones
│   │   ├── faces
│   │   ├── faceZones
│   │   ├── neighbour
│   │   ├── owner
│   │   ├── points
│   │   ├── pointZones
│   │   └── sets
│   ├── RASProperties
│   ├── transportProperties
│   └── turbulenceProperties
└── system
    ├── controlDict
    ├── fvSchemes
    └── fvSolution

4.2.3.6. 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

cd my-turbulent-airfoil/

renumberMesh -overwrite
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

Checking geometry...
  Overall domain bounding box (-484.456616748 0 -507.806809185) (501.000007802 10 507.806808887)
  Mesh (non-empty, non-wedge) directions (1 0 1)
  Mesh (non-empty) directions (1 0 1)
  All edges aligned with or perpendicular to non-empty directions.
  Boundary openness (8.28364670182e-18 1.55457412622e-15 -1.3301799376e-17) 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.57960157047e-11. Maximum face area = 843.777295738.  Face area magnitudes OK.
  Min volume = 9.57960157047e-10. Max volume = 8437.77295738.  Total volume = 8755584.35466.  Cell volumes OK.
  Mesh non-orthogonality Max: 19.7377051719 average: 1.63498336689
  Non-orthogonality check OK.
  Face pyramids OK.
  Max skewness = 0.201755632336 OK.
  Coupled point location match (average 0) OK.

Failed 1 mesh checks.

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 multi-core 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 should be carried out through a file decomposeParDict present in the system sub-directory as shown below.

cd my-turbulent-airfoil/system

touch decomposeParDict
cd system/
ls
controlDict  decomposeParDict  fvSchemes  fvSolution

Now, the following contents should be added to the decomposeParDict file

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
    version     		2.0;
    format      		ascii;
    class       		dictionary;
    object      		decomposeParDict;
}
//--------------------------------------------------------------------------------

numberOfSubdomains 		8;

method          		scotch;

distributed     		no;

roots
(
);

In the above file, the the keyword numberOfSubdomains defines the number of decomposed sub-domains needed and 8 is used which partitions the grid into 8 sub-domains. 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

cd my-turbulent-airfoil/

decomposePar

Now the decomposition should begin and the details of which are displayed in the terminal window. Subsequently, 8 processor directories will be generated as we have requested for 8 divisions of grid as shown below

cd my-turbulent-airfoil/
ls
0  constant  processor0  processor1  processor2  processor3  processor4  processor5  processor6  processor7  system

The solver can now be executed for parallel computation from the top directory using the following command

cd my-turbulent-airfoil/
mpirun -np 8 simpleSolver -parallel | tee my-turbulent-airfoil.log

Note that here it is assumed that the parallel computing is available in the host machine. If the parallel computing is carried out on other hosts, the command is as follows

cd my-turbulent-airfoil/
mpirun -np 8 -machinefile machine simpleSolver -parallel | tee my-turbulent-airfoil.log

where, machines is a plain text file located in the top directory with the list of available hosts names. With the execution of the above commands, the simulation begins and the progress of the solution can be seen in the terminal window and simultaneously it is written to the specified log file (my-turbulent-airfoil.log).

The log file can be further processed to look at the convergence history and this can be done as follows once the simulation is complete

cd my-turbulent-airfoil/
foamLog my-turbulent-airfoil.log
cd logs/
xmgrace p_0

The xmgrace p_0 allows you to look at the convergence of pressure with respect to the number of iterations carried out and is shown in Fig. 103. Other convergence variables can also be accessed similarly in the logs/ sub-directory.

_images/t-airfoil-convergence-tutorials.png

Figure 103: Convergence of pressure with respect to iterations

4.2.3.7. Results

The flow over the airfoil at both \(0^\circ\) and \(10^\circ\) degree angle of attack are presented here. In Fig. 104, 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 Fig. 105, a low pressure region exists on the upper surface and consequently the velocity increases thus generating some lift.

_images/t-airfoil-velocitypressure-tutorials-0.png

Figure 104: Velocity magnitude and pressure contours for \(\alpha = 0^\circ\) angle of attack

_images/t-airfoil-velocitypressure-tutorials-10.png

Figure 105: Velocity magnitude and pressure contours for \(\alpha = 10^\circ\) angle of attack

4.2.4. Convex Curvature

Turbulent flow in a constant-area duct having a convex curvature is considered in this tutorial. Caelus 4.10 will be used and the process of setting-up 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.

4.2.4.1. Objectives

Through this tutorial, the user will get familiarise with setting up the Caelus simulation for steady, turbulent, incompressible flow in a two-dimensional duct having a convex curvature and subsequently post-process 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

4.2.4.2. Pre-requisites

It will be assumed that the user is comfortable and familiar with the Linux command line environment using a terminal. 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.

4.2.4.3. Background

Turbulent flow in a constant-area 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 steady-state and incompressible. Non-dimensional shear-stress (skin-friction coefficient) will be used to show the influence of curvature on the flow. Validation and verification of this exercise is detailed in section Two-dimensional 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 Fig. 106 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 geometric-plane for this case in 2-D is \(x-z\) plane.

Figure 106: Schematic representation of the 2D curvature geometry

The inflow conditions for this case is given in the below table (Freestream conditions)

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

Velocity components in x and z directions
\(\alpha\) \(u~(m/s)\) \(w~(m/s)\)
\(30^\circ\) 27.63389 15.95443

4.2.4.4. 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 Fig. 107, 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 no-slip boundary. The inlet and the outlet are applied to the start and end of the domain respectively.

Figure 107: Computational domain of the convex curvature

The grid that is used, polyMesh format is in three-dimensions. 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 \(x-z\) directions. As with all the previous cases, one-cell thick normal to the 2D plane is sufficient, assuming the flow to be symmetry. The two \(x-z\) 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 third-dimension, 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.

_images/t-curvature-grid-tutorials.png

Figure 108: Convex curvature grid in two-dimensions in \(x-z\) plane

The 2D grid in \(x-z\) plane is shown in Fig. 108 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 wall-function is used.

4.2.4.5. Problem definition

This section details the various steps required to set-up the turbulent flow inside a convex duct. A full working case of this can be found in:

/tutorials-cases/turbulent-curvature/turbulent-curvature

Directory Structure

A top level directory is firstly needed where all the necessary simulation will be carried out. Let us name this directory as my-turbulent-curvature. Following the procedure below in a terminal window.

cd turbulent-curvature
mkdir my-turbulent-curvature
cd my-turbulent-curvature/

Note

All commands shown here are entered in a terminal window, unless otherwise mentioned

A few more sub-directories are needed in the top-level directory to set-up the case. Caelus requires time, constant and system sub-directories within the main my-turbulent-curvature working directory. Since we start the simulation at time, t =0 s, a time sub-directory named 0 is required. Again, using a terminal

cd my-turbulent-curvature/
mkdir 0
mkdir constant
mkdir system

ls
0  constant  system

The 0 sub-directory requires additional 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\)) p
Velocity (\(U\)) U
Turbulent viscosity (\(\nu\)) nut
Turbulence field variable (\(\tilde{\nu}\)) nuTilda (Only for SA & SA-CC model)
Turbulent kinetic energy (\(k\)) k (Only for \(k-\omega~\rm{SST}\) model)
Turbulent dissipation rate (\(\omega\)) omega (Only for \(k-\omega~\rm{SST}\) model)

We consider simulating this case with three turbulence models, namely Spalart-Allmaras (SA), Spalart–Allmaras Rotational/Curvature (SA-RC) and \(k-\omega\) - Shear Stress Transport (\(\rm{SST}\)). The contents that are entered in 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 creating these in the 0 sub-directory by

cd my-turbulent-curvature/
touch 0/p
touch 0/U
touch 0/nut
touch 0/nuTilda
touch 0/k
touch 0/omega

cd 0/
ls
k  nut  nuTilda  omega  p  U

Caelus is case sensitive and therefore the user should carefully set-up the case as shown here.

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.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:
  • No-slip wall
    • Velocity: Fixed uniform velocity \(u, v, w = 0\)
    • Pressure: Zero gradient
    • Turbulence:
      • SA & SA-RC (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 & SA-RC (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:

First, we begin with the pressure file, p and using a text editor add the following contents

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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 non-uniform. 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 can be added as follows

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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			fixedValue;
		value			uniform (0 0 0);
	}
}

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]). 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 and add the following

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
    	version				2.0;
    	format				ascii;
    	class				volScalarField;
    	location			"0";
    	object				nut;
}
//--------------------------------------------------------------------------------

dimensions      			[0 2 -1 0 0 0 0];

//internalField   			uniform 3.197546e-06; // Enable this line for SA Model
//internalField   			uniform 1.36756e-07; // Enable this line for k-omega SST Model

boundaryField
{
	inlet
	{
		type			fixedValue;
		//value			uniform 3.197546e-06; // Enable this line for SA Model
		//value			uniform 1.36756e-07; // Enable this line for k-omega SST Model
	}
	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 SA-RC 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. The following contents should be added to the file nuTilda and the entries specified for the boundaryField are identical to that of turbulent kinematic viscosity explained above.

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
    	version				2.0;
    	format				ascii;
    	class				volScalarField;
    	location			"0";
    	object				nuTilda;
}
//--------------------------------------------------------------------------------

dimensions				[0 2 -1 0 0 0 0];

internalField				uniform 4.55841e-05;

boundaryField
{
	inlet
	{
		type			fixedValue;
		value			uniform 4.55841e-05;
	}
	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, add the below contents to it

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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 1e-10;
	}
	wall
	{
		type				fixedValue;
		value				uniform 1e-10;
	}
}

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 wall-function, a small, non-zero fixedValue is required.

We can now proceed to the file omega and the value for this is evaluated as detailed in Turbulent freestream conditions for SST Model

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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 \(x-z\) plane named here as left and right have empty boundary conditions which forces Caelus to assume the flow to be in two-dimensions. 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 sub-directory, which is in the constant directory. The polyMesh directory can not be created as follows

cd my-turbulent-curvature/
mkdir constant/polyMesh

The use should now place the grid in the polyMesh sub-directory. Further to is, the physical properties should be specified in various different files present in the constant directory. Let us begin creating these files. First, navigate to the constant directory

touch constant/RASProperties
touch constant/transportProperties
touch constant/turbulenceProperties

cd constant/
ls
polyMesh  RASProperties  transportProperties  turbulenceProperties

As noted above, the three files are listed in the constant sub-directory. The first one is the RASProperties where, Reynolds-Average-Stress (RAS) model is specified as shown below. b-directory. The first one, RASProperties in which the Reynolds-Average-Stress (RAS) 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 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
	version				2.0;
	format				ascii;
	class				dictionary;
	location			"constant";
	object				RASProperties;
}

//--------------------------------------------------------------------------------

// For Spalarat-Alamaras Model, enable the below line
// RASModel				SpalartAllmarasCurvature;

// For k-omega SST Model, enable the below line
// RASModel				koSST;

turbulence				on;

printCoeffs				on;

Second from the list is the transportProperties file, where kinematic viscosity is specified as shown below

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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.519470E-5;

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 4.10
Web:   www.caelus-cml.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. First, navigate to the system directory and create the following files

touch system/controlDict
touch system/fvSchemes
touch system/fvSolution

cd system/
ls
controlDict  fvSchemes  fvSolution

In the controlDict file, add the following details

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
	version				2.0;
	format				ascii;
	class				dictionary;
	location			"system";
	object				controlDict;
}

// Enable the below line for SA and SA-RC turbulence Models
// libs ( "libincompressibleSpalartAlmarasCurvature.so" );

// Enable the below line for SST turbulence Model
// libs ("libincompressibleKOSST.so");
			
//-------------------------------------------------------------------------------

application				simpleFOAM;

startFrom				startTime;

startTime				0;

stopAt					endTime;

endTime					10000;

deltaT					1;

writeControl				runTime;

writeInterval				1000;

purgeWrite				0;

writeFormat				ascii;

writePrecision				12;

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.2084;
        outputControl   		timeStep;
        outputInterval  		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 steady-state 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 4.10
Web:   www.caelus-cml.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 & SA-RC model only
	div(phi,k)			Gauss 	upwind; // will be used for k-epsilon & k-omega only
	div(phi,omega)			Gauss 	upwind;	// Will be used for k-omega model only
	div((nuEff*dev(T(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 & SA-RC model only
	laplacian(DkEff,k)		Gauss 	linear 	corrected; // Will be used for k-omega & k-omega only
	laplacian(DomegaEff,omega)	Gauss 	linear 	corrected; // Will be used for k-omega model only
	laplacian((1|A(U)),p)		Gauss 	linear	corrected;
	laplacian(1,p)			Gauss 	linear 	corrected;
}

interpolationSchemes
{
	default				linear;
	interpolate(HbyA)		linear;
}

snGradschemes
{
	default				corrected;
}

fluxRequired
{
	default				no;
	p				;
}

The linear solver controls and tolerances are set in fvSolution as given below

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
	version				2.0;
	format				ascii;
	class				dictionary;
	location			"system";
	object				fvSolution;
}

//------------------------------------------------------------------------------

solvers
{
	p
	{
		solver			PCG;
		preconditioner		DIC;
		tolerance		1e-8;
		relTol			0.01;
	}
	U
	{
		solver			PBiCG;
		preconditioner		DILU;
		tolerance		1e-7;
		relTol			0.01;
	}

	"(k|omega|nuTilda)"
	{
		solver          	PBiCG;
		preconditioner  	DILU;
		tolerance       	1e-08;
		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 time-accurate set-up as we are solving for a steady-state 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 set-up 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 my-turbulent-curvature/
tree
.
├── 0
│   ├── k
│   ├── nut
│   ├── nuTilda
│   ├── omega
│   ├── p
│   └── U
├── constant
│   ├── polyMesh
│   │   ├── boundary
│   │   ├── cellZones
│   │   ├── faces
│   │   ├── faceZones
│   │   ├── neighbour
│   │   ├── owner
│   │   ├── points
│   │   ├── pointZones
│   │   └── sets
│   ├── RASProperties
│   ├── transportProperties
│   └── turbulenceProperties
└── system
    ├── controlDict
    ├── decomposeParDict
    ├── fvSchemes
    └── fvSolution

4.2.4.6. 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 my-turbulent-curvature/

renumberMesh -overwrite
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

Checking geometry...
  Overall domain bounding box (-1.5 -1 -0.8097167) (1.6 0 0.127)
  Mesh (non-empty, non-wedge) directions (1 0 1)
  Mesh (non-empty) directions (1 0 1)
  All edges aligned with or perpendicular to non-empty directions.
  Boundary openness (2.89884272002e-17 2.78435627054e-16 2.89924087017e-17) 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.52715761849e-08. Maximum face area = 0.0318702444976.  Face area magnitudes OK.
  Min volume = 1.52715761849e-08. Max volume = 0.000180011068556.  Total volume = 0.417212754326.  Cell volumes OK.
  Mesh non-orthogonality Max: 27.7205032543 average: 17.6783210246
  Non-orthogonality check OK.
  Face pyramids OK.
  Max skewness = 0.0357278065799 OK.
  Coupled point location match (average 0) OK.

Failed 1 mesh checks.

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 multi-core 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 sub-directory which is shown below.

cd my-turbulent-curvature/system

touch decomposeParDict

cd system/
ls
controlDict  decomposeParDict  fvSchemes  fvSolution

At this stage, the following contents should be added to the decomposeParDict file

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
    version     		2.0;
    format      		ascii;
    class       		dictionary;
    object      		decomposeParDict;
}
//--------------------------------------------------------------------------------

numberOfSubdomains 		8;

method          		scotch;

distributed     		no;

roots
(
);

In the above file, the the keyword numberOfSubdomains defines the number of decomposed sub-domains needed and 8 is used which partitions the grid into 8 sub-domains. 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

cd my-turbulent-curvature/

decomposePar

Now the decomposition should start and the details of which are displayed in the terminal window. Subsequently, 8 processor directories will be generated as we have requested for 8 divisions of grid as shown below

cd my-turbulent-curvature/
ls
0  constant  processor0  processor1  processor2  processor3  processor4  processor5  processor6  processor7  system

The solver can now be executed for parallel computation from the top directory using the following command

cd my-turbulent-curvature/
mpirun -np 8 simpleSolver -parallel | tee my-turbulent-curvature.log

Note that here it is assumed that the parallel computing is available in the host machine. If the parallel computing is carried out on other hosts, the command is as follows

cd my-turbulent-curvature/
mpirun -np 8 -machinefile machine simpleSolver -parallel | tee my-turbulent-curvature.log

where, machines is a plain text file located in the top directory with the list of available hosts names. With the execution of the above commands, the simulation begins and the progress of the solution can be seen in the terminal window and simultaneously it is written to the specified log file (my-turbulent-curvature.log).

The log file can be further processed to look at the convergence history and this can be done as follows once the simulation is complete

cd my-turbulent-curvature/
foamLog my-turbulent-curvature.log
cd logs/
xmgrace p_0

The xmgrace p_0 allows you to look at the convergence of pressure with respect to the number of iterations carried out as shown in Fig. 109. Other convergence variables can also be accessed similarly in the logs/ sub-directory.

_images/t-curvature-convergence-tutorials.png

Figure 109: Convergence of pressure with respect to iterations

4.2.4.7. Results

The flow visualisation of velocity and pressure within the convex duct is presented here. In Fig. 110 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.

_images/t-curvature-velocitypressure-tutorials.png

Figure 110: Velocity magnitude and pressure contours within the convex duct

4.2.5. Backward Facing Step

This tutorial focuses on the turbulent flow over a two-dimensional backward facing step using Caelus 4.10. In particular this tutorial emphasises the use of wall functions for a high \(y^+\) grid. The basic steps to set-up 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.

4.2.5.1. Objectives

Through this tutorial, the user will get familiar in setting up Caelus simulation for steady, turbulent, incompressible flow over a two-dimensional 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

4.2.5.2. Pre-requisites

By now the user should be familiar with the Linux command line environment via a terminal. 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.

4.2.5.3. 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 Two-dimensional Backward Facing Step for more detailed analysis.

The schematic of the backward facing step is shown in Fig. 111 in two-dimensions. 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 two-dimensional geometric plane considered is in \(x-z\) directions.

Figure 111: Schematic representation of the backward facing step

Freestream conditions are detailed in the below table

Freestream conditions
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}\)

4.2.5.4. 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-cases/turbulent-step/turbulent-step/grid

The computational domain for the step is shown in Fig. 112 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 Fig. 112 and closely flows the experimental set-up. The velocity on the step surfaces is zero, wherein \(u,v,w = 0\) represented through a no-slip boundary.

Figure 112: Computational domain of a 2D backward facing step

The grid in polyMesh is in three-dimensions, although the flow over the step is two-dimensional in nature. A one-cell thick grid normal to the flow plane (\(x-z\)) is therefore sufficient for Caelus to consider symmetry flow in the \(y\) direction. The two resulting planes that are in (\(x-z\)) 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 third-dimension 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.

_images/t-step-grid-tutorials.png

Figure 113: Computational grid of a 2D backward facing step in \(x-z\) plane

The 2D step grid in \(x-z\) plane is shown in Fig. 113 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.

4.2.5.5. Problem definition

This section details various steps needed to set-up the turbulent flow over a step. A full working case of this can be found in:

/tutorials-cases/turbulent-step/turbulent-step

In the following sections, details to complete this tutorial assuming no prior set-up or data is available to the user.

Directory Structure

The first step to begin this tutorial would be to create a top level working directory where all the necessary simulation will be be carried out. For example, let us name this directory as my-turbulent-step.

cd turbulent-step
mkdir my-turbulent-step
cd my-turbulent-step/

Note

All commands shown here are entered in a terminal window, unless otherwise mentioned

In order to set-up the problem in my-turbulent-step, a few more sub-directories are needed where relevant files will be created. Caelus requires time, constant and system sub-directories within the main my-turbulent-step working directory. Here, the simulation will be started at time t = 0 s, which requires a time sub-directory named 0. With the use of a terminal window,

cd my-turbulent-step/
mkdir 0
mkdir constant
mkdir system

ls
0  constant  system

The 0 sub-directory requires 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\)) p
Velocity (\(U\)) U
Turbulent viscosity (\(\nu\)) nut
Turbulence field variable (\(\tilde{\nu}\)) nuTilda (Only for SA model)
Turbulent kinetic energy (\(k\)) k (Only for \(k-\omega~\rm{SST}\) model)
Turbulent dissipation rate (\(\omega\)) omega (Only for \(k-\omega~\rm{SST}\) model)
Turbulent dissipation (\(\epsilon\)) epsilon (Only for R \(k-\epsilon\) model)

In this tutorial we consider three turbulence models. They are Spalart-Allmaras (SA), \(k-\omega\) - Shear Stress Transport (\(\rm{SST}\)) and Realizable \(k-epsilon\) 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. First we begin with creating these files in the 0 sub-directory through

cd my-turbulent-step/
touch 0/p
touch 0/U
touch 0/nut
touch 0/nuTilda
touch 0/k
touch 0/omega
touch 0/epsilon

cd 0/
ls
epsilon k  nut  nuTilda  omega  p  U

The user should note that Caelus is case sensitive and therefore the directory and file set-up should be identical to what is shown here.

Boundary Conditions and Solver Attributes

Boundary Conditions

Referring back to Fig. 112, the following are the boundary conditions that will be specified:

  • Inlet
  • Symmetry
    • Velocity: Symmetry
    • Pressure: Symmetry
    • Turbulence: Symmetry
  • No-slip 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

First, we begin with the pressure and using a text editor add the following contents to the file named p

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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;
	}
	symm-left
	{
		type            empty;
	}
	symm-right
	{
		type            empty;
	}
	top-wall
	{
		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 non-uniform. 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, we can open the file U and add the following to it

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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;
	}
	symm-left
	{
		type            empty;
	}
	symm-right
	{
		type            empty;
	}
	top-wall
	{
		type            fixedValue;
		value           uniform (0 0 0);
	}
	upstream
	{
		type            symmetryPlane;
	}
	wall
	{
		type            fixedValue;
		value           uniform (0 0 0);
	}
}

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]). 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);. 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. First, we start with opening the file nut, which is turbulent kinematic viscosity and add the following

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
	version     	2.0;
	format      	ascii;
	class       	volScalarField;
	location    	"0";
	object      	nut;
}
// -------------------------------------------------------------------------------

dimensions      	[0 2 -1 0 0 0 0];

//internalField   	uniform 0.00025904507; // Enable this line for SA Model

//internalField   	uniform 0; // Enable this line for k-omega SST Model

//internalField   	uniform 0.09811; // Enable this line for Realizable k-epsilon Model

boundaryField
{
	inlet
	{
		type            calculated;
        //value			unifrom 0.00025904507; // Enable this line for SA Model
        //value			uniform 0; // Enable this line for k-omega SST Model			
		//value			uniform 0.09811; //Enable this line for Realizable k-epsilon Model
	}
	outlet
	{
		type            calculated;
		value           uniform 0;
	}
	symm-left
	{
		type            empty;
	}
	symm-right
	{
		type            empty;
	}
	top-wall
	{
		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 Model and Turbulent freestream conditions for Realizable Model for SA, SST and RKE models respectively and is specified accordingly. The same value also goes for internalField.

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 Turbulence freestream conditions for SA model. The following contents should be added to the file nuTilda and the entries specified for the boundaryField are identical to that of turbulent kinematic viscosity explained above.

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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;
	}
	symm-left
	{
		type            empty;
	}
	symm-right
	{
		type            empty;
	}
	top-wall
	{
		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 \(k-omega~\rm{SST}\) and Realizable \(k-epsilon\) models and represents turbulent kinetic energy.

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/
FoamFile
{
	version			2.0;
	format      		ascii;
	class       		volScalarField;
	location    		"0";
	object      		k;
}
// -------------------------------------------------------------------------------

dimensions      		[0 2 -2 0 0 0 0];

//internalField   		uniform 0.00109611; // Enable this line for k-omega SST Model
//internalField   		uniform 0.2945755; // Enable this line for Realizable k-epsilon Model

boundaryField
{
	inlet
	{
		type            fixedValue;
		//value         uniform 0.00109611; // Enable this line for k-omega SST Model
		//value         uniform 0.2945755; // Enable this line for Realizable k-epsilon Model
	}
	outlet
	{
		type            zeroGradient;
	}
	symm-left
	{
		type            empty;
	}
	symm-right
	{
		type            empty;
	}
	top-wall
	{
		type            kqRWallFunction;
		//value         uniform 0.00109611; // Enable this line for k-omega SST Model
		//value         uniform 0.2945755; // Enable this line for Realizable k-epsilon Model
    }
    upstream
	{
		type            symmetryPlane;
	}
	wall
	{
		type			kqRWallFunction;
		//value         uniform 0.00109611; // Enable this line for k-omega SST Model
		//value         uniform 0.2945755; // Enable this line for Realizable k-epsilon Model
	}
}

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 Model and Turbulent freestream conditions for Realizable 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 Model

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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;
	}
	symm-left
	{
		type            empty;
	}
	symm-right
	{
		type            empty;
	}
	top-wall
	{
		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 Model

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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;
	}
	symm-left
	{
		type            empty;
	}
	symm-right
	{
		type            empty;
	}
	top-wall
	{
		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 \(x-z\) plane named here as symm-left and symm-right have empty boundary conditions which forces Caelus to assume the flow to be in two-dimensions. 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 sub-directory, which resides in the constant directory. We use identical grid for all the three turbulence simulation. The polyMesh directory can be created as follows

cd my-turbulent-step/
mkdir constant/polyMesh

The grid files can now be copied into the polyMesh sub-directory. In addition, the physical properties are specified in various different files present in the constant directory. We can start creating these files. First, navigate to the constant directory

touch constant/RASProperties
touch constant/transportProperties
touch constant/turbulenceProperties

cd constant/
ls
polyMesh  RASProperties  transportProperties  turbulenceProperties

From the above information, there are three files listed in addition to the polyMesh sub-directory. The first one, RASProperties in which the Reynolds-Average-Stress (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 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
	version			2.0;
	format			ascii;
	class			dictionary;
	location		"constant";
	object			RASProperties;
}

//--------------------------------------------------------------------------------

// For Spalarat-Alamaras Model, enable the below line
//RASModel			SpalartAllmarasCurvature;

// For k-omega SST Model, enable the below line
//RASModel			koSST;

// For Realizable k-epsilon 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 4.10
Web:   www.caelus-cml.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.230979E-3;

//Cross Power law and Bird-Carreau non-linear viscosity models
// ----------------------------------------------------------------------------
CrossPowerLawCoeffs
{
	nu0             	nu0 [0 2 -1 0 0 0 0] 1e-06;
	nuInf           	nuInf [0 2 -1 0 0 0 0] 1e-06;
	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] 1e-06;
	nuInf           	nuInf [0 2 -1 0 0 0 0] 1e-06;
	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 4.10
Web:   www.caelus-cml.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 and create the following files

touch system/controlDict
touch system/fvSchemes
touch system/fvSolution

cd system/
ls
controlDict  fvSchemes  fvSolution

In the controlDict file, add the following details

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
	version			2.0;
	format			ascii;
	class			dictionary;
	location		"system";
	object			controlDict;
}
			
// Enable the below line for Spalart-Almaras turbulence Model
// libs ( "libincompressibleSpalartAlmarasCurvature.so" );

// Enable the below line for SST turbulence Model
// libs ("libincompressibleKOSST.so");

//--------------------------------------------------------------------------------

application			simpleSolver;

startFrom			latestTime;

startTime			0;

stopAt				endTime;

endTime				10000;

deltaT				1.0;

writeControl			runTime;

writeInterval			500;

purgeWrite			0;

writeFormat			ascii;

writePrecision			12;

writeCompression 		uncompressed;

timeFormat			general;

timePrecision			12;

runTimeModifiable		true;

//-------------------------------------------------------------------------------

// Function Objects to obtain mean values

/*functions
{
	forces
	{
		type		forces;

        functionObjectLibs 	("libforces.so");
        patches     		( wall );
        CofR      		(0 0 0);
        rhoName         	rhoInf;
        rhoInf          	15.13009E-3;
        outputControl   	timeStep;
        outputInterval  	100;
     }
}

With reference 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 steady-state 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 and the contents should be copied

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.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 S-A model only
	div(phi,k)      		Gauss	upwind; 	                    // will be used for k-epsilon & k-omega only
	div(phi,omega)  		Gauss 	upwind;				// Will be used for k-omega 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 S-A-Curvature-corr 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 S-A model only
	laplacian(DkEff,k) 		Gauss 	linear 	corrected;        // Will be used for k-omega & k-omega only
	laplacian(DomegaEff,omega) 	Gauss 	linear 	corrected;	// Will be used for k-omega model only
	laplacian(DepsilonEff,epsilon)	Gauss 	linear 	corrected;
	laplacian((1|A(U)),p)		Gauss 	linear 		corrected;
	laplacian(1,p)  		Gauss 	linear 		corrected;
}

interpolationSchemes
{
	default				linear;
	interpolate(HbyA)		linear;
}

snGradschemes
{
	default				corrected;
}

fluxRequired
{
	default				no;
	p				;
}
//--------------------------------------------------------------------------------

The linear solver controls and tolerances are set in fvSolution as given below

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/
FoamFile
{
	version				2.0;
	format				ascii;
	class				dictionary;
	location			"system";
	object				fvSolution;
}

//------------------------------------------------------------------------------

solvers
{
	p
	{
		solver			PCG;
		preconditioner		DIC;
		tolerance		1e-8;
		relTol			0.01;
	}
	U
	{
		solver			PBiCG;
		preconditioner		DILU;
		tolerance		1e-7;
		relTol			0.01;
	}

	"(k|omega|nuTilda|epsilon)"
    {
        solver          		PBiCG;
        preconditioner  		DILU;
        tolerance       		1e-08;
        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 time-accurate set-up as we are solving for a steady-state 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 set-up 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 my-turbulent-step/
tree
.
├── 0
│   ├── epsilon
│   ├── k
│   ├── nut
│   ├── nuTilda
│   ├── omega
│   ├── p
│   └── U
├── constant
│   ├── polyMesh
│   │   ├── boundary
│   │   ├── cellZones
│   │   ├── faces
│   │   ├── faceZones
│   │   ├── neighbour
│   │   ├── owner
│   │   ├── points
│   │   ├── pointZones
│   │   └── sets
│   ├── RASProperties
│   ├── transportProperties
│   └── turbulenceProperties
└── system
    ├── controlDict
    ├── fvSchemes
    └── fvSolution

4.2.5.6. 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

cd my-turbulent-step/

renumberMesh -overwrite
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

Checking geometry...
  Overall domain bounding box (-130 0 0) (50 1 9)
  Mesh (non-empty, non-wedge) directions (1 0 1)
  Mesh (non-empty) directions (1 0 1)
  All edges aligned with or perpendicular to non-empty directions.
  Boundary openness (7.60426729195e-19 -3.46752522136e-15 6.3479100872e-18) OK.
  Max cell openness = 2.1851506874e-16 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 non-orthogonality Max: 12.1317921302 average: 0.100202600449
  Non-orthogonality check OK.
  Face pyramids OK.
  Max skewness = 0.00821683161155 OK.
  Coupled point location match (average 0) OK.

Similar to the previous tutorial, it will be shown here to utilise the multi-core 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 should be carried out through a file decomposeParDict present in the system sub-directory as shown below.

cd my-turbulent-step/system

touch decomposeParDict
cd system/
ls
controlDict  decomposeParDict  fvSchemes  fvSolution

Now, the following contents should be added to the decomposeParDict file

/*-------------------------------------------------------------------------------*
Caelus 4.10
Web:   www.caelus-cml.com
\*------------------------------------------------------------------------------*/

FoamFile
{
    version     		2.0;
    format      		ascii;
    class       		dictionary;
    object      		decomposeParDict;
}
//--------------------------------------------------------------------------------

numberOfSubdomains 		8;

method          		scotch;

distributed     		no;

roots
(
);

As noted in the above file, the the keyword numberOfSubdomains defines the number of decomposed sub-domains needed and 8 is used which partitions the grid into 8 sub-domains. 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

cd my-turbulent-step/

decomposePar

With this, the decomposition should begin and the details of which are displayed in the terminal window. Subsequently, 8 processor directories will be generated as we have requested for 8 divisions of grid as shown below

cd my-turbulent-step/
ls
0  constant  processor0  processor1  processor2  processor3  processor4  processor5  processor6  processor7  system

At this stage, the solver can now be executed for parallel computation from the top directory using the following command

cd my-turbulent-step/
mpirun -np 8 simpleSolver -parallel | tee my-turbulent-step.log

Note that here it is assumed that the parallel computing is available in the host machine. If the parallel computing is carried out on other hosts, the command is as follows

cd my-turbulent-step/
mpirun -np 8 -machinefile machine simpleSolver -parallel | tee my-turbulent-step.log

where, machines is a plain text file located in the top directory with the list of available hosts names. With the execution of the above commands, the simulation begins and the progress of the solution can be seen in the terminal window and simultaneously it is written to the specified log file (my-turbulent-step.log).

The log file can be further processed to look at the convergence history and this can be done as follows once the simulation is complete

cd my-turbulent-step/
foamLog my-turbulent-step.log
cd logs/
xmgrace p_0

The xmgrace p_0 allows you to look at the convergence of pressure with respect to the number of iterations carried out and is shown in Fig. 114. Other convergence variables can also be accessed similarly in the logs/ sub-directory.

_images/t-step-convergence-tutorials.png

Figure 114: Convergence of pressure with respect to iterations

4.2.5.7. Results

The flow within the separated region is visualised here through the velocity and pressure contours. In Fig. 115 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.

_images/t-step-velocitypressure-tutorials.png

Figure 115: Velocity magnitude and pressure contours over the backward facing step