Experiment 1: External Flow Analysis

over NACA 0012 Airfoil

1.1 Objective
Analyze external flow over a NACA 0012 airfoil at different angles of attack using the
simpleFoam solver to understand pressure distribution, lift, and drag characteristics.

1.2 Theory
The NACA 0012 airfoil is defined by the equation:
t  √
0.2969 x − 0.1260x − 0.3516x2 + 0.2843x3 − 0.1015x4 (1.1)

y=±
0.2
where t = 0.12 for the 12% thickness ratio.
The lift coefficient is given by:
L
CL = 1 (1.2)
2
ρU 2 c

where L is lift force, ρ is density, U is free-stream velocity, and c is chord length.

The drag coefficient is:
D
CD = 1 2 (1.3)
2
ρU c
For potential flow over symmetric airfoil at small angles:

CL = 2π sin(α) ≈ 2πα (thin airfoil theory) (1.4)

The pressure coefficient is defined as:

p − p∞
Cp = 1 2
(1.5)
2
ρU∞

1
CFD Lab Manual - OpenFOAM 2406 2

Figure 1.1: Numerical domain to be simulated

CFD Lab Manual - OpenFOAM 2406 3

1.3 Governing Equations

For incompressible, steady flow, the Reynolds-Averaged Navier-Stokes (RANS) equations
are:

∇·U=0 (1.6)
∇ · (UU) = −∇p + ν∇2 U + ∇ · τ turbulent (1.7)

For turbulent flow, the k-ω SST model is used:

 
∂k ∂k ∂ ∂k
+ Uj ∗
= Pk − β kω + (ν + σk νt ) (1.8)
∂t ∂xj ∂xj ∂xj
 
∂ω ∂ω ∂ ∂ω 1 ∂k ∂ω
+ Uj 2
= γS − βω + 2
(ν + σω νt ) + 2(1 − F1 )σω2 (1.9)
∂t ∂xj ∂xj ∂xj ω ∂xj ∂xj

1.4 Detailed Procedure

1.4.1 Step 1: Case Setup
Navigate to OpenFOAM tutorials and copy the airfoil case:
cd $FOAM_RUN
cp -r $FOAM_TUTORIALS / incompressible / simpleFoam / airFoil2D .
cd airFoil2D

1.4.2 Step 2: Geometry and Mesh Generation

The tutorial uses a pre-generated mesh. Examine the mesh:
checkMesh

For custom angle of attack, modify system/blockMeshDict:

// In blockMeshDict , modify the transformation matrix
// For angle of attack alpha ( in radians ) :
// x_new = x * cos ( alpha ) + y * sin ( alpha )
// y_new = -x * sin ( alpha ) + y * cos ( alpha )

1.4.3 Step 3: Boundary Conditions

Modify 0/U for inlet velocity and angle of attack:
inlet
{
type fixedValue ;
value uniform (25.75 0 0) ; // For 0 degree AoA
}

outlet
{
type zeroGradient ;
}
CFD Lab Manual - OpenFOAM 2406 4

wall
{
type noSlip ;
}

frontAndBack
{
type empty ;
}

For angle of attack α: - Ux = U cos(α) = 25.75 cos(α) - Uy = U sin(α) = 25.75 sin(α)

Modify 0/p:
inlet
{
type zeroGradient ;
}

outlet
{
type fixedValue ;
value uniform 0;
}

wall
{
type zeroGradient ;
}

1.4.4 Step 4: Turbulence Properties

Modify 0/k (turbulent kinetic energy):
inlet
{
type fixedValue ;
value uniform 0.375; // k = 1.5*( U * I ) ^2 , I =0.01
}

outlet
{
type zeroGradient ;
}

wall
{
type kqRWallFunction ;
value uniform 0.375;
}

Modify 0/omega (specific dissipation rate):

inlet
{
type fixedValue ;
value uniform 1708; // omega = k ^0.5/(0.09^0.25* L ) , L
=0.07* chord
}
CFD Lab Manual - OpenFOAM 2406 5

outlet
{
type zeroGradient ;
}

wall
{
type ome gaWallFu nction ;
value uniform 1708;
}

1.4.5 Step 5: Transport Properties

Modify constant/transportProperties:
transportModel Newtonian ;
nu [0 2 -1 0 0 0 0] 1.516 e -05; // For Re = 1 e6

1.4.6 Step 6: Turbulence Properties

Ensure constant/momentumTransport contains:
simulationType RAS ;
RAS
{
model kOmegaSST ;
turbulence on ;
printCoeffs on ;
}

1.4.7 Step 7: Solution Control

Modify system/fvSolution:
solvers
{
p
{
solver GAMG ;
tolerance 1e -06;
relTol 0.1;
smoother GaussSeidel ;
}

"( U | k | omega ) "

{
solver smoothSolver ;
smoother symGaussSeidel ;
tolerance 1e -05;
relTol 0.1;
}
}

SIMPLE
CFD Lab Manual - OpenFOAM 2406 6

{
n N o n O r t h o g o n a l C o r r e c t o r s 0;
consistent yes ;
residualControl
{
p 1e -4;
U 1e -5;
"( k | omega ) " 1e -5;
}
}

relax ationFac tors

{
equations
{
U 0.9;
k 0.7;
omega 0.7;
}
}

1.4.8 Step 8: Control Dictionary

Modify system/controlDict:
startFrom startTime ;
startTime 0;
stopAt endTime ;
endTime 2000;
deltaT 1;
writeControl timeStep ;
writeInterval 500;
writeFormat ascii ;
writePrecision 6;
writeCompression off ;
timeFormat general ;
timePrecision 6;
runTi meModifi able true ;

functions
{
forces
{
type forces ;
libs (" libforces . so ") ;
writeControl timeStep ;
writeInterval 10;
patches ( wall ) ;
rho rhoInf ;
rhoInf 1.225;
CofR (0.25 0 0) ;
}

coeffs
{
type forceCoeffs ;
libs (" libforces . so ") ;
CFD Lab Manual - OpenFOAM 2406 7

writeControl timeStep ;
writeInterval 10;
patches ( wall ) ;
rho rhoInf ;
rhoInf 1.225;
liftDir (0 1 0) ;
dragDir (1 0 0) ;
CofR (0.25 0 0) ;
lRef 1;
Aref 1;
}
}

1.4.9 Step 9: Running the Simulation

Execute the following commands:
# Generate mesh
blockMesh

# Check mesh quality

checkMesh

# Initialize fields
simpleFoam

# Monitor residuals
tail -f log . simpleFoam

1.4.10 Step 10: Post-Processing

Calculate forces and coefficients:
# Calculate pressure coefficient
foamCalc - time 2000 Cp

# Extract data along airfoil surface

sample - time 2000

# Launch ParaView
paraview

Create sampling dictionary system/sampleDict:

type sets ;
libs (" libsampling . so ") ;
i nt e rp o l at i on S c he m e cellPoint ;
setFormat raw ;
sets
(
airfoilUpper
{
type face ;
axis xyz ;
patches ( wall ) ;
}
CFD Lab Manual - OpenFOAM 2406 8

);
fields ( p U Cp ) ;

1.5 Expected Results

For NACA 0012 at Re = 1 × 106 :

• At α = 0: CL ≈ 0, CD ≈ 0.008

• At α = 5: CL ≈ 0.5, CD ≈ 0.01

• At α = 10: CL ≈ 1.0, CD ≈ 0.015

• Stall angle: approximately 12-15°

CFD Lab Manual - OpenFOAM 2406 9

1.6 Assignments
1.6.1 Assignment 1
Analyze the NACA 0012 airfoil at angles of attack: 0°, 5°, 10°, 15°. Plot CL vs α curve
and compare with theoretical thin airfoil theory: CL = 2π sin(α).

1.6.2 Assignment 2
Investigate the effect of mesh refinement on solution accuracy. Use three different mesh
densities and analyze convergence of lift and drag coefficients. Calculate Grid Conver-
gence Index (GCI).

1.6.3 Assignment 3
Implement NACA 2412 airfoil geometry and compare aerodynamic characteristics with
NACA 0012. Analyze the effect of camber on pressure distribution and stall characteris-
tics.
Experiment 2: Fluid Flow and Heat
Transfer in Double Pipe Heat
Exchanger

2.1 Objective
Simulate conjugate heat transfer in a double pipe heat exchanger using chtMultiRegionFoam
to analyze temperature distribution and heat transfer coefficients.

2.2 Theory
The overall heat transfer coefficient for a double pipe heat exchanger is:

1 1 ro ln(ro /ri ) ro
= + + (2.10)
U hi kwall ri ho
The effectiveness-NTU relationship for counter-flow configuration:

1 − exp(−N T U (1 − Cr ))
ε= (2.11)
1 − Cr exp(−N T U (1 − Cr ))

where Cr = Cmin /Cmax and N T U = U A/Cmin .

Heat transfer rate:
Q = εCmin (Th,in − Tc,in ) (2.12)
Log mean temperature difference:
∆T1 − ∆T2
LM T D = (2.13)
ln(∆T1 /∆T2 )

Local Nusselt number for internal flow:

N u = 0.023Re0.8 P r0.4 (Dittus-Boelter correlation) (2.14)

10
CFD Lab Manual - OpenFOAM 2406 11

Figure 2.2: Domain to be simulated (simplified demonstration)

2.3 Governing Equations

Energy equation for fluid regions:

∂(ρh)
+ ∇ · (ρUh) = ∇ · (αef f ∇h) + Sh (2.15)
∂t
Heat conduction in solid regions:

∂(ρcp T )
= ∇ · (k∇T ) (2.16)
∂t
Turbulent thermal diffusivity:
νt
αef f = α + (2.17)
P rt

2.4 Detailed Procedure

2.4.1 Step 1: Base Case Setup
Since no direct double pipe tutorial exists, we’ll create from the multiRegion tutorial:
CFD Lab Manual - OpenFOAM 2406 12

cd $FOAM_RUN
cp -r $FOAM_TUTORIALS / heatTransfer / ch tMu lt iR egi on Fo am / multiReg ionHeate r
.
mv multi RegionHe ater doublePipeHX
cd doublePipeHX

2.4.2 Step 2: Geometry Creation

Create system/blockMeshDict for concentric pipes:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object blockMeshDict ;
}

convertToMeters 0.001; // mm to m

vertices
(
// Inner pipe vertices ( hot fluid )
(0 0 0) // 0
(1000 0 0) // 1
(1000 10 0) // 2
(0 10 0) // 3
(0 0 10) // 4
(1000 0 10) // 5
(1000 10 10) // 6
(0 10 10) // 7

// Pipe wall vertices

(0 0 -2) // 8
(1000 0 -2) // 9
(1000 12 -2) // 10
(0 12 -2) // 11
(0 0 12) // 12
(1000 0 12) // 13
(1000 12 12) // 14
(0 12 12) // 15

// Outer annulus vertices ( cold fluid )

(0 0 -5) // 16
(1000 0 -5) // 17
(1000 20 -5) // 18
(0 20 -5) // 19
(0 0 20) // 20
(1000 0 20) // 21
(1000 20 20) // 22
(0 20 20) // 23
);

blocks
(
// Inner fluid region
CFD Lab Manual - OpenFOAM 2406 13

hex (0 1 2 3 4 5 6 7) innerFluid (100 10 1) simpleGrading (1 1 1)

// Pipe wall region

hex (8 9 10 11 12 13 14 15) pipeWall (100 2 1) simpleGrading (1 1
1)

// Outer fluid region

hex (16 17 18 19 20 21 22 23) outerFluid (100 10 1) simpleGrading
(1 1 1)
);

edges () ;

boundary
(
innerInlet
{
type patch ;
faces ((0 4 7 3) ) ;
}
innerOutlet
{
type patch ;
faces ((1 2 6 5) ) ;
}
outerInlet
{
type patch ;
faces ((17 18 22 21) ) ;
}
outerOutlet
{
type patch ;
faces ((16 19 23 20) ) ;
}
frontAndBack
{
type empty ;
faces
(
(0 1 2 3)
(4 5 6 7)
(8 9 10 11)
(12 13 14 15)
(16 17 18 19)
(20 21 22 23)
);
}
);

2.4.3 Step 3: Region Definition

Create system/topoSetDict:
FoamFile
{
version 2.0;
CFD Lab Manual - OpenFOAM 2406 14

format ascii ;
class dictionary ;
object topoSetDict ;
}

actions
(
{
name innerFluid ;
type cellSet ;
action new ;
source boxToCell ;
sourceInfo
{
box (0 0 0) (1000 10 10) ;
}
}

{
name pipeWall ;
type cellSet ;
action new ;
source boxToCell ;
sourceInfo
{
box (0 0 -2) (1000 12 12) ;
}
}

{
name outerFluid ;
type cellSet ;
action new ;
source boxToCell ;
sourceInfo
{
box (0 0 -5) (1000 20 20) ;
}
}
);

2.4.4 Step 4: Material Properties

Create constant/innerFluid/transportProperties:
transportModel Newtonian ;
nu [0 2 -1 0 0 0 0] 1e -06; // Water at 60 C
Pr [0 0 0 0 0 0 0] 4.0; // Prandtl number
Prt [0 0 0 0 0 0 0] 0.85; // Turbulent Prandtl number

Create constant/outerFluid/transportProperties:
transportModel Newtonian ;
nu [0 2 -1 0 0 0 0] 1.5 e -05; // Air at 20 C
Pr [0 0 0 0 0 0 0] 0.71; // Prandtl number
Prt [0 0 0 0 0 0 0] 0.85; // Turbulent Prandtl number

Create constant/pipeWall/transportProperties:
CFD Lab Manual - OpenFOAM 2406 15

kappa [1 1 -3 -1 0 0 0] 45; // Steel thermal

conductivity
Cp [0 2 -2 -1 0 0 0] 500; // Specific heat
rho [1 -3 0 0 0 0 0] 7800; // Density

2.4.5 Step 5: Initial and Boundary Conditions

Create 0/innerFluid/T:
FoamFile
{
version 2.0;
format ascii ;
class volScalarField ;
object T;
}

dimensions [0 0 0 1 0 0 0];
internalField uniform 363; // 90 C

boundaryField
{
innerInlet
{
type fixedValue ;
value uniform 363; // 90 C hot water inlet
}

innerOutlet
{
type zeroGradient ;
}

innerFluid_to_pipeWall
{
type compressible ::
turbulentTemperatureCoupledBaffleMixed ;
Tnbr T;
kappaMethod fluidThermo ;
value uniform 363;
}

frontAndBack
{
type empty ;
}
}

Create 0/outerFluid/T:
FoamFile
{
version 2.0;
format ascii ;
class volScalarField ;
object T;
}
CFD Lab Manual - OpenFOAM 2406 16

dimensions [0 0 0 1 0 0 0];
internalField uniform 293; // 20 C

boundaryField
{
outerInlet
{
type fixedValue ;
value uniform 293; // 20 C cold air inlet
}

outerOutlet
{
type zeroGradient ;
}

outerFluid_to_pipeWall
{
type compressible ::
turbulentTemperatureCoupledBaffleMixed ;
Tnbr T;
kappaMethod fluidThermo ;
value uniform 293;
}

frontAndBack
{
type empty ;
}
}

2.4.6 Step 6: Velocity Boundary Conditions

Create 0/innerFluid/U:
FoamFile
{
version 2.0;
format ascii ;
class volVectorField ;
object U;
}

dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0) ;

boundaryField
{
innerInlet
{
type fixedValue ;
value uniform (2 0 0) ; // 2 m / s inlet velocity
}

innerOutlet
{
CFD Lab Manual - OpenFOAM 2406 17

type zeroGradient ;
}

innerFluid_to_pipeWall
{
type noSlip ;
}

frontAndBack
{
type empty ;
}
}

2.4.7 Step 7: Solver Control

Modify system/fvSolution:
solvers
{
rho
{
solver PCG ;
preconditioner DIC ;
tolerance 1e -05;
relTol 0.1;
}

p rgh
{
solver GAMG ;
tolerance 1e -06;
relTol 0.01;
smoother GaussSeidel ;
}

"( U | h | k | epsilon | omega ) "

{
solver smoothSolver ;
smoother symGaussSeidel ;
tolerance 1e -05;
relTol 0.1;
}
}

PIMPLE
{
nOuterCorrectors 2;
nCorrectors 1;
n N o n O r t h o g o n a l C o r r e c t o r s 0;
}

relax ationFac tors

{
equations
{
".*" 0.9;
CFD Lab Manual - OpenFOAM 2406 18

}
}

2.4.8 Step 8: Running the Simulation

Execute the following commands:
# Generate mesh
blockMesh

# Create regions
topoSet
splitMeshRegions - cellZones

# Generate meshes for each region

blockMesh - region innerFluid
blockMesh - region outerFluid
blockMesh - region pipeWall

# Run simulation
ch tM ul tiR eg io nFo am

# Monitor convergence
tail -f log . ch tMu lt iR egi on Fo am

2.4.9 Step 9: Post-Processing

Calculate heat transfer coefficients:
# Calculate wall heat flux
postProcess - func wallHeatFlux - region innerFluid - time 1000

# Calculate effectiveness
postProcess - func " mag ( T ) " - time 1000

# Launch ParaView for all regions

paraview -- state = system / state . pvsm

2.5 Expected Results

For typical operating conditions:

• Hot fluid outlet temperature: 50-70C

• Cold fluid outlet temperature: 30-40C

• Heat exchanger effectiveness: 0.6-0.8

• Overall heat transfer coefficient: 100-500 W/m2 K

• Pressure drop: 1000-5000 Pa

CFD Lab Manual - OpenFOAM 2406 19

2.6 Assignments
2.6.1 Assignment 1
Vary the mass flow rates of hot and cold fluids and analyze the effect on heat exchanger
effectiveness. Compare results with analytical NTU-effectiveness charts.

2.6.2 Assignment 2
Investigate the effect of pipe wall thermal conductivity (steel vs copper) on overall heat
transfer performance. Calculate the thermal resistance distribution.

2.6.3 Assignment 3
Design an optimal double pipe heat exchanger for given inlet conditions by varying pipe
diameters and length. Perform parametric study to maximize heat transfer while mini-
mizing pressure drop.
Experiment 3: Internal Fluid Flow and
Heat Transfer in Centrifugal Pump

3.1 Objective
Analyze internal flow patterns, pressure distribution, and heat generation in a centrifugal
pump using pimpleFoam with moving reference frame (MRF).

3.2 Theory
The pump head is calculated as:

∆p V22 − V12
H= + + z2 − z1 (3.18)
ρg 2g
Pump efficiency:
Hydraulic Power ρgQH
η= = (3.19)
Shaft Power Pshaf t
The velocity triangle relationships at impeller exit:

Vu = U − Vr cot β (3.20)
U2 Vu2 − U1 Vu1
Hideal = (3.21)
g
where U is blade tip speed, Vr is relative velocity, β is blade angle.
Specific speed: √
n Q
ns = 3/4 (3.22)
H
Affinity laws:
Q2 N2
= (3.23)
Q1 N1
 2
H2 N2
= (3.24)
H1 N1
 3
P2 N2
= (3.25)
P1 N1

20
CFD Lab Manual - OpenFOAM 2406 21

Figure 3.3: Domain to be simulated (students will simplify this as per the instructions)

3.3 Governing Equations

Navier-Stokes equations in rotating reference frame:
∂U
+ (U · ∇)U = −∇p + ν∇2 U − 2Ω × U − Ω × (Ω × r) (3.26)
∂t
CFD Lab Manual - OpenFOAM 2406 22

Additional source terms in rotating frame:

FCoriolis = −2Ω × U (3.27)

Fcentrif ugal = −Ω × (Ω × r) (3.28)

3.4 Detailed Procedure

3.4.1 Step 1: Base Case Setup
Use the MRF tutorial as starting point:
cd $FOAM_RUN
cp -r $FOAM_TUTORIALS / incompressible / pimpleFoam / RAS / propeller .
mv propeller centrifugalPump
cd centrifugalPump

3.4.2 Step 2: Geometry Creation

Create 2D centrifugal pump geometry using system/blockMeshDict:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object blockMeshDict ;
}

convertToMeters 0.001; // mm to m

vertices
(
// Impeller inlet eye
(0 0 0) // 0 - center
(20 0 0) // 1
(0 20 0) // 2
( -20 0 0) // 3
(0 -20 0) // 4

// Impeller exit
(100 0 0) // 5
(0 100 0) // 6
( -100 0 0) // 7
(0 -100 0) // 8

// Volute casing
(150 0 0) // 9
(0 150 0) // 10
( -150 0 0) // 11
(0 -150 0) // 12

// Outlet pipe
(200 -10 0) // 13
(200 10 0) // 14
(300 -10 0) // 15
CFD Lab Manual - OpenFOAM 2406 23

(300 10 0) // 16

// Z - direction (2 D case )
(0 0 5) // 17 - center
(20 0 5) // 18
(0 20 5) // 19
( -20 0 5) // 20
(0 -20 5) // 21
(100 0 5) // 22
(0 100 5) // 23
( -100 0 5) // 24
(0 -100 5) // 25
(150 0 5) // 26
(0 150 5) // 27
( -150 0 5) // 28
(0 -150 5) // 29
(200 -10 5) // 30
(200 10 5) // 31
(300 -10 5) // 32
(300 10 5) // 33
);

blocks
(
// Impeller region ( rotating )
hex (0 1 2 0 17 18 19 17) impeller (20 20 1) simpleGrading (1 1 1)
hex (0 2 3 0 17 19 20 17) impeller (20 20 1) simpleGrading (1 1 1)
hex (0 3 4 0 17 20 21 17) impeller (20 20 1) simpleGrading (1 1 1)
hex (0 4 1 0 17 21 18 17) impeller (20 20 1) simpleGrading (1 1 1)

// Expanded impeller region

hex (1 5 6 2 18 22 23 19) impeller (30 30 1) simpleGrading (1 1 1)
hex (2 6 7 3 19 23 24 20) impeller (30 30 1) simpleGrading (1 1 1)
hex (3 7 8 4 20 24 25 21) impeller (30 30 1) simpleGrading (1 1 1)
hex (4 8 5 1 21 25 22 18) impeller (30 30 1) simpleGrading (1 1 1)

// Volute region ( stationary )

hex (5 9 10 6 22 26 27 23) volute (20 20 1) simpleGrading (1 1 1)
hex (6 10 11 7 23 27 28 24) volute (20 20 1) simpleGrading (1 1 1)
hex (7 11 12 8 24 28 29 25) volute (20 20 1) simpleGrading (1 1 1)
hex (8 12 9 5 25 29 26 22) volute (20 20 1) simpleGrading (1 1 1)

// Outlet pipe
hex (13 15 16 14 30 32 33 31) volute (20 5 1) simpleGrading (1 1 1)
);

edges () ;

boundary
(
inlet
{
type patch ;
faces
(
(0 17 19 2)
(0 2 20 17)
(0 17 20 3)
CFD Lab Manual - OpenFOAM 2406 24

(0 3 21 17)
(0 17 21 4)
(0 4 18 17)
(0 17 18 1)
(0 1 19 17)
);
}

outlet
{
type patch ;
faces ((15 16 33 32) ) ;
}

impellerWalls
{
type wall ;
faces
(
// Impeller blade surfaces
(1 5 22 18)
(5 6 23 22)
// Add more blade surfaces
);
}

casing
{
type wall ;
faces
(
(9 10 27 26)
(10 11 28 27)
(11 12 29 28)
(12 9 26 29)
);
}

frontAndBack
{
type empty ;
faces () ; // Define all front and back faces
}
);

3.4.3 Step 3: MRF Zone Definition

Create constant/MRFProperties:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object MRFProperties ;
}
CFD Lab Manual - OpenFOAM 2406 25

MRF1
{
cellZone impeller ;
active yes ;

// Rotation axis and origin

no nR ot ati ng Pa tch es () ;
origin (0 0 0) ;
axis (0 0 1) ;
omega 157.08; // rad / s (1500 RPM )
}

3.4.4 Step 4: Cell Zone Definition

Create system/topoSetDict:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object topoSetDict ;
}

actions
(
{
name impeller ;
type cellSet ;
action new ;
source cylinderToCell ;
sourceInfo
{
p1 (0 0 -1) ;
p2 (0 0 10) ;
radius 80;
}
}

{
name impeller ;
type cellZoneSet ;
action new ;
source setToCellZone ;
sourceInfo
{
set impeller ;
}
}
);

3.4.5 Step 5: Boundary Conditions

Create 0/U:
FoamFile
CFD Lab Manual - OpenFOAM 2406 26

{
version 2.0;
format ascii ;
class volVectorField ;
object U;
}

dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0) ;

boundaryField
{
inlet
{
type fixedValue ;
value uniform (0 0 -5) ; // Axial inlet velocity
}

outlet
{
type zeroGradient ;
}

impellerWalls
{
type mov in gW all Ve lo cit y ;
value uniform (0 0 0) ;
}

casing
{
type noSlip ;
}

frontAndBack
{
type empty ;
}
}

Create 0/p:
FoamFile
{
version 2.0;
format ascii ;
class volScalarField ;
object p;
}

dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;

boundaryField
{
inlet
{
type zeroGradient ;
}
CFD Lab Manual - OpenFOAM 2406 27

outlet
{
type fixedValue ;
value uniform 0;
}

impellerWalls
{
type zeroGradient ;
}

casing
{
type zeroGradient ;
}

frontAndBack
{
type empty ;
}
}

3.4.6 Step 6: Turbulence Model Setup

Create 0/k:
FoamFile
{
version 2.0;
format ascii ;
class volScalarField ;
object k;
}

dimensions [0 2 -2 0 0 0 0];
internalField uniform 0.375; // k = 1.5*( U * I ) ^2 , I =0.1

boundaryField
{
inlet
{
type fixedValue ;
value uniform 0.375;
}

outlet
{
type zeroGradient ;
}

impellerWalls
{
type kqRWallFunction ;
value uniform 0.375;
}
CFD Lab Manual - OpenFOAM 2406 28

casing
{
type kqRWallFunction ;
value uniform 0.375;
}

frontAndBack
{
type empty ;
}
}

3.4.7 Step 7: Control Dictionary

Modify system/controlDict:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object controlDict ;
}

startFrom startTime ;
startTime 0;
stopAt endTime ;
endTime 10;
deltaT 0.001;
writeControl timeStep ;
writeInterval 100;

functions
{
forces
{
type forces ;
libs (" libforces . so ") ;
writeControl timeStep ;
writeInterval 10;
patches ( impellerWalls ) ;
rho rhoInf ;
rhoInf 1000;
CofR (0 0 0) ;
}

torque
{
type forces ;
libs (" libforces . so ") ;
writeControl timeStep ;
writeInterval 10;
patches ( impellerWalls ) ;
rho rhoInf ;
rhoInf 1000;
CofR (0 0 0) ;
CFD Lab Manual - OpenFOAM 2406 29

binData
{
nBin 1;
direction (0 0 1) ;
cumulative yes ;
}
}

flowRate
{
type sur faceFiel dValue ;
libs (" l i b f i e l d F u n c t i o n O b j e c t s . so ") ;
writeControl timeStep ;
writeInterval 10;
fields ( phi ) ;
operation sum ;
regionType patch ;
name outlet ;
}
}

3.4.8 Step 8: Solver Settings

Modify system/fvSolution:
solvers
{
p
{
solver GAMG ;
tolerance 1e -06;
relTol 0.01;
smoother GaussSeidel ;
}

pFinal
{
$p ;
relTol 0;
}

"( U | k | omega ) "

{
solver smoothSolver ;
smoother symGaussSeidel ;
tolerance 1e -05;
relTol 0.1;
}

"( U | k | omega ) Final "

{
$U ;
relTol 0;
}
}

PIMPLE
CFD Lab Manual - OpenFOAM 2406 30

{
nOuterCorrectors 2;
nCorrectors 1;
n N o n O r t h o g o n a l C o r r e c t o r s 0;
pRefCell 0;
pRefValue 0;
}

relax ationFac tors

{
equations
{
".*" 0.9;
}
}

3.4.9 Step 9: Running the Simulation

Execute the following commands:
# Generate mesh
blockMesh

# Create cell zones for MRF

topoSet

# Check mesh quality

checkMesh

# Run simulation with MRF

pimpleFoam

# Monitor convergence and forces

tail -f log . pimpleFoam

3.4.10 Step 10: Post-Processing

Analyze pump performance:
# Calculate head rise
foamCalc - time ’5: ’ components p

# Plot pump characteristic curves

gnuplot plotResults . gp

# Visualize flow patterns

paraview

Create [Link]:
set xlabel " Flow Rate ( $m ^3 $ / s ) "
set ylabel " Head ( m ) "
plot " postProcessing / forces /0/ forces . dat " using 1:2 with lines title "
Head "
pause -1
CFD Lab Manual - OpenFOAM 2406 31

3.5 Expected Results

For a typical centrifugal pump at 1500 RPM:

• Best efficiency point: Q = 0.05 m3 /s, H = 20 m

• Maximum head (shut-off): H = 25 m

• Maximum efficiency: η = 75-85%

• Power consumption: 10-15 kW

• Typical performance curve shape with decreasing head vs increasing flow

3.6 Assignments
3.6.1 Assignment 1
Generate pump characteristic curves (H-Q, P-Q, η-Q) by varying flow rates. Compare
with typical centrifugal pump performance curves.

3.6.2 Assignment 2
Analyze the effect of impeller rotational speed on pump performance using similarity
laws. Validate the relationships: Q ∝ N , H ∝ N 2 , P ∝ N 3 .

3.6.3 Assignment 3
Investigate cavitation inception by analyzing pressure distribution at the impeller inlet.
Determine NPSH requirements and critical operating conditions.
Experiment 4: Analysis of Karman
Vortex Street over Cylinder

4.1 Objective
Simulate vortex shedding behind a circular cylinder using pimpleFoam to analyze Strouhal
number and unsteady wake characteristics.

4.2 Theory
The Strouhal number is defined as:
fD
St = (4.29)
U
where f is vortex shedding frequency, D is cylinder diameter, and U is free-stream veloc-
ity.
For Reynolds numbers 40 < Re < 1000:
 
19.7
St ≈ 0.198 1 − (4.30)
Re

For higher Reynolds numbers (Re > 1000):

 
21.2
St ≈ 0.212 1 − (4.31)
Re

Drag coefficient correlations:

1.2 13.6
CD = √ + + 0.4 for Re < 2 × 105 (4.32)
Re Re

The vorticity magnitude:

∂v ∂u
ω= − (4.33)
∂x ∂y
Reynolds number based on cylinder diameter:
UD
Re = (4.34)
ν

32
CFD Lab Manual - OpenFOAM 2406 33

4.3 Governing Equations

Unsteady Navier-Stokes equations:
∂U
+ ∇ · (UU) = −∇p + ν∇2 U (4.35)
∂t
∇·U=0 (4.36)

Vorticity transport equation:

∂ω
+ U · ∇ω = ν∇2 ω + ω · ∇U (4.37)
∂t

4.4 Detailed Procedure

4.4.1 Step 1: Base Case Setup
Use cylinder tutorial as starting point:
cd $FOAM_RUN
cp -r $FOAM_TUTORIALS / incompressible / pimpleFoam / laminar / cylinder2D .
cd cylinder2D

If tutorial doesn’t exist, create new case:

cd $FOAM_RUN
mkdir cylinder2D
cd cylinder2D
mkdir 0 constant system
CFD Lab Manual - OpenFOAM 2406 34

Figure 4.4: Logic for creating mesh around the cylinder

4.4.2 Step 2: Geometry and Mesh Creation

Create system/blockMeshDict with O-grid around cylinder:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object blockMeshDict ;
}

convertToMeters 1;

// Cylinder parameters
radius 0.5; // Cylinder radius
domain 20; // Domain extent
CFD Lab Manual - OpenFOAM 2406 35

vertices
(
// Inner cylinder boundary ( radius = 0.5)
( 0.3536 0.3536 0) // 0
( -0.3536 0.3536 0) // 1
( -0.3536 -0.3536 0) // 2
( 0.3536 -0.3536 0) // 3

// Middle layer ( radius = 2)

( 1.414 1.414 0) // 4
( -1.414 1.414 0) // 5
( -1.414 -1.414 0) // 6
( 1.414 -1.414 0) // 7

// Outer boundary
( 20 10 0) // 8
( -5 10 0) // 9
( -5 -10 0) // 10
( 20 -10 0) // 11

// Z - direction (2 D case )
( 0.3536 0.3536 0.1) // 12
( -0.3536 0.3536 0.1) // 13
( -0.3536 -0.3536 0.1) // 14
( 0.3536 -0.3536 0.1) // 15
( 1.414 1.414 0.1) // 16
( -1.414 1.414 0.1) // 17
( -1.414 -1.414 0.1) // 18
( 1.414 -1.414 0.1) // 19
( 20 10 0.1) // 20
( -5 10 0.1) // 21
( -5 -10 0.1) // 22
( 20 -10 0.1) // 23
);

blocks
(
// Inner O - grid blocks around cylinder
hex (0 4 5 1 12 16 17 13) (20 20 1) simpleGrading (1 1 1)
hex (1 5 6 2 13 17 18 14) (20 20 1) simpleGrading (1 1 1)
hex (2 6 7 3 14 18 19 15) (20 20 1) simpleGrading (1 1 1)
hex (3 7 4 0 15 19 16 12) (20 20 1) simpleGrading (1 1 1)

// Outer blocks extending to far field

hex (4 8 9 5 16 20 21 17) (80 20 1) simpleGrading (2 1 1)
hex (5 9 10 6 17 21 22 18) (30 20 1) simpleGrading (1 1 1)
hex (6 10 11 7 18 22 23 19) (80 20 1) simpleGrading (2 1 1)
hex (7 11 8 4 19 23 20 16) (80 20 1) simpleGrading (2 1 1)
);

edges
(
// Inner circle ( cylinder surface )
arc 0 1 (0 0.5 0)
arc 1 2 ( -0.5 0 0)
arc 2 3 (0 -0.5 0)
arc 3 0 (0.5 0 0)
CFD Lab Manual - OpenFOAM 2406 36

arc 12 13 (0 0.5 0.1)

arc 13 14 ( -0.5 0 0.1)
arc 14 15 (0 -0.5 0.1)
arc 15 12 (0.5 0 0.1)

// Middle circle
arc 4 5 (0 2 0)
arc 5 6 ( -2 0 0)
arc 6 7 (0 -2 0)
arc 7 4 (2 0 0)
arc 16 17 (0 2 0.1)
arc 17 18 ( -2 0 0.1)
arc 18 19 (0 -2 0.1)
arc 19 16 (2 0 0.1)
);

boundary
(
inlet
{
type patch ;
faces
(
(9 10 22 21)
);
}

outlet
{
type patch ;
faces
(
(8 20 23 11)
);
}

top
{
type symmetryPlane ;
faces
(
(8 9 21 20)
);
}

bottom
{
type symmetryPlane ;
faces
(
(10 11 23 22)
);
}

cylinder
{
type wall ;
faces
CFD Lab Manual - OpenFOAM 2406 37

(
(0 12 13 1)
(1 13 14 2)
(2 14 15 3)
(3 15 12 0)
);
}

frontAndBack
{
type empty ;
faces
(
// All front and back faces
(0 1 5 4)
(1 2 6 5)
(2 3 7 6)
(3 0 4 7)
(4 5 9 8)
(5 6 10 9)
(6 7 11 10)
(7 4 8 11)
(12 16 17 13)
(13 17 18 14)
(14 18 19 15)
(15 19 16 12)
(16 20 21 17)
(17 21 22 18)
(18 22 23 19)
(19 23 20 16)
);
}
);

mergePatchPairs () ;

4.4.3 Step 3: Initial and Boundary Conditions

Create 0/U:
FoamFile
{
version 2.0;
format ascii ;
class volVectorField ;
object U;
}

dimensions [0 1 -1 0 0 0 0];
internalField uniform (1 0 0) ;

boundaryField
{
inlet
{
type fixedValue ;
value uniform (1 0 0) ; // Re = 100 for D =1 , nu =0.01
CFD Lab Manual - OpenFOAM 2406 38

outlet
{
type zeroGradient ;
}

top
{
type symmetryPlane ;
}

bottom
{
type symmetryPlane ;
}

cylinder
{
type noSlip ;
}

frontAndBack
{
type empty ;
}
}

Create 0/p:
FoamFile
{
version 2.0;
format ascii ;
class volScalarField ;
object p;
}

dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;

boundaryField
{
inlet
{
type zeroGradient ;
}

outlet
{
type fixedValue ;
value uniform 0;
}

top
{
type symmetryPlane ;
}
CFD Lab Manual - OpenFOAM 2406 39

bottom
{
type symmetryPlane ;
}

cylinder
{
type zeroGradient ;
}

frontAndBack
{
type empty ;
}
}

4.4.4 Step 4: Transport Properties

Create constant/transportProperties:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object tr a ns p o rt P ro p e rt i es ;
}

transportModel Newtonian ;
nu [0 2 -1 0 0 0 0] 0.01; // For Re = 100 with U =1 , D =1

4.4.5 Step 5: Control Dictionary with Probes

Create system/controlDict:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object controlDict ;
}

startFrom startTime ;
startTime 0;
stopAt endTime ;
endTime 200; // Long enough to establish periodic
shedding
deltaT 0.005; // Small time step for accuracy
writeControl timeStep ;
writeInterval 100;
writeFormat ascii ;
writePrecision 8;
writeCompression off ;
timeFormat general ;
timePrecision 6;
CFD Lab Manual - OpenFOAM 2406 40

runTi meModifi able true ;

// Courant number limit

maxCo 0.5;

functions
{
forces
{
type forces ;
libs (" libforces . so ") ;
writeControl timeStep ;
writeInterval 1;
patches ( cylinder ) ;
rho rhoInf ;
rhoInf 1;
CofR (0 0 0) ;
}

forceCoeffs
{
type forceCoeffs ;
libs (" libforces . so ") ;
writeControl timeStep ;
writeInterval 1;
patches ( cylinder ) ;
rho rhoInf ;
rhoInf 1;
liftDir (0 1 0) ;
dragDir (1 0 0) ;
CofR (0 0 0) ;
lRef 1; // Cylinder diameter
Aref 1; // Reference area ( D *1 for 2 D )
}

probes
{
type probes ;
libs (" libsampling . so ") ;
writeControl timeStep ;
writeInterval 1;
fields (U p);
probeLocations
(
(1 0.5 0.05) // Wake monitoring point
(2 0 0.05) // Downstream centerline
(3 0 0.05) // Far wake
);
}

vorticity
{
type vorticity ;
libs (" l i b f i e l d F u n c t i o n O b j e c t s . so ") ;
writeControl writeTime ;
writeInterval 1;
}
}
CFD Lab Manual - OpenFOAM 2406 41

4.4.6 Step 6: Numerical Schemes

Create system/fvSchemes:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object fvSchemes ;
}

ddtSchemes
{
default Euler ;
}

gradSchemes
{
default Gauss linear ;
}

divSchemes
{
default none ;
div ( phi , U ) Gauss linearUpwind grad ( U ) ;
div ( phi , k ) Gauss limitedLinear 1;
div ( phi , omega ) Gauss limitedLinear 1;
div (( nuEff * dev2 ( T ( grad ( U ) ) ) ) ) Gauss linear ;
}

laplacianSchemes
{
default Gauss linear orthogonal ;
}

i n t e r p o l a t i on S c h e m e s
{
default linear ;
}

snGradSchemes
{
default orthogonal ;
}

wallDist
{
method meshWave ;
}

4.4.7 Step 7: Solution Control

Create system/fvSolution:
FoamFile
{
version 2.0;
CFD Lab Manual - OpenFOAM 2406 42

format ascii ;
class dictionary ;
object fvSolution ;
}

solvers
{
p
{
solver GAMG ;
tolerance 1e -06;
relTol 0.01;
smoother GaussSeidel ;
}

pFinal
{
$p ;
relTol 0;
}

U
{
solver smoothSolver ;
smoother symGaussSeidel ;
tolerance 1e -05;
relTol 0.1;
}

UFinal
{
$U ;
relTol 0;
}
}

PIMPLE
{
nOuterCorrectors 2;
nCorrectors 1;
n N o n O r t h o g o n a l C o r r e c t o r s 0;
pRefCell 0;
pRefValue 0;
}

relax ationFac tors

{
equations
{
".*" 1;
}
}

4.4.8 Step 8: Running the Simulation

Execute the following commands:
CFD Lab Manual - OpenFOAM 2406 43

# Generate mesh
blockMesh

# Check mesh quality

checkMesh

# Check Courant number

pimpleFoam - checkMesh

# Run simulation
pimpleFoam > log . pimpleFoam 2 >&1 &

# Monitor residuals and forces

tail -f log . pimpleFoam
gnuplot - persist -e " plot ’ postProcessing / forceCoeffs /0/ forceCoeffs . dat
’ using 1:3 with lines "

4.4.9 Step 9: Post-Processing and Analysis

Calculate Strouhal number:
# Extract force coefficients
cd postProcessing / forceCoeffs /0
tail -n +4 forceCoeffs . dat > cl_data . txt

# Use Python / MATLAB for FFT analysis

python3 << EOF
import numpy as np
import matplotlib . pyplot as plt
from scipy import signal

# Load data
data = np . loadtxt ( ’ cl_data . txt ’)
time = data [: , 0]
cl = data [: , 2] # Lift coefficient

# Remove initial transient ( first 50 time units )

mask = time > 50
time = time [ mask ]
cl = cl [ mask ]

# Calculate sampling frequency

dt = time [1] - time [0]
fs = 1.0 / dt

# Perform FFT
frequencies , psd = signal . welch ( cl , fs , nperseg = len ( cl ) //4)

# Find dominant frequency

peak_idx = np . argmax ( psd [1:]) + 1 # Exclude DC component
f_shed = frequencies [ peak_idx ]

# Calculate Strouhal number ( U =1 , D =1)

St = f_shed * 1.0 / 1.0
print ( f " Shedding frequency : { f_shed :.3 f } Hz ")
print ( f " Strouhal number : { St :.3 f }")
CFD Lab Manual - OpenFOAM 2406 44

# Plot results
plt . figure ( figsize =(12 , 4) )
plt . subplot (1 , 2 , 1)
plt . plot ( time , cl )
plt . xlabel ( ’ Time ’)
plt . ylabel ( ’ Lift Coefficient ’)
plt . title ( ’ Lift Coefficient vs Time ’)

plt . subplot (1 , 2 , 2)
plt . loglog ( frequencies [1:] , psd [1:])
plt . axvline ( f_shed , color = ’r ’ , linestyle = ’ - - ’ , label =f ’ f = { f_shed :.3 f
} ’)
plt . xlabel ( ’ Frequency ( Hz ) ’)
plt . ylabel ( ’ PSD ’)
plt . title ( ’ Power Spectral Density ’)
plt . legend ()
plt . tight_layout ()
plt . savefig ( ’ strouh al_analy sis . png ’)
plt . show ()
EOF

Visualize vorticity field:

# Calculate vorticity
postProcess - func vorticity - time 100:200

# Launch ParaView
paraview
# Load case and visualize vorticity contours
# Create animation of vortex shedding

4.5 Expected Results

For Reynolds number Re = 100:
• Strouhal number: St ≈ 0.165
• Drag coefficient: CD ≈ 1.35 ± 0.05
• Lift coefficient amplitude: CL,rms ≈ 0.25
• Shedding frequency: f ≈ 0.165 Hz (for U = 1 m/s, D = 1 m)
• Wake development length: ≈ 4 − 6 diameters

Figure 4.5: Distinct vortex shedding should be seen behind the cylinder
CFD Lab Manual - OpenFOAM 2406 45

4.6 Assignments
4.6.1 Assignment 1
Investigate the effect of Reynolds number (100, 200, 500) on Strouhal number and com-
pare with experimental correlations. Plot St vs Re relationship.

4.6.2 Assignment 2
Analyze the transition from steady to unsteady flow by gradually increasing Reynolds
number. Determine the critical Reynolds number for vortex shedding onset.

4.6.3 Assignment 3
Study the effect of blockage ratio (cylinder diameter to domain height ratio) on vortex
shedding characteristics and develop correction factors.
Experiment 5: Couette Flow Analysis -
Explicit and Implicit Methods

5.1 Objective
Implement and compare explicit and implicit finite difference schemes for solving the
unsteady Couette flow problem using custom OpenFOAM solvers and analyze numerical
stability and accuracy.

5.2 Theory
The governing equation for unsteady Couette flow between parallel plates:
∂u ∂ 2u
=ν 2 (5.38)
∂t ∂y
Initial condition: u(y, 0) = 0 for 0 ≤ y ≤ h
Boundary conditions:
u(0, t) = 0 (5.39)
u(h, t) = Uw (5.40)
Analytical solution using separation of variables:
∞  2 2 
Uw y X 2Uw  nπy  n π νt
u(y, t) = + n
(−1) sin exp − (5.41)
h n=1
nπ h h2
Steady-state solution:
Uw y
usteady (y) = (5.42)
h
Dimensionless time:
νt
τ= (5.43)
h2

Figure 5.6: Flow conditions that need to be simulated (with proper assumptions made
by the student)

46
CFD Lab Manual - OpenFOAM 2406 47

5.3 Finite Difference Schemes

5.3.1 Explicit Scheme (Forward Euler in Time)
un+1 − uni uni+1 − 2uni + uni−1
i
=ν (5.44)
∆t (∆y)2
Rearranging:
ν∆t n
un+1 = uni + (u − 2uni + uni−1 ) (5.45)
i
(∆y)2 i+1
Stability condition (von Neumann analysis):
ν∆t 1
2
≤ (5.46)
(∆y) 2

5.3.2 Implicit Scheme (Backward Euler in Time)

un+1 − uni un+1 − 2un+1 + un+1
i
= ν i+1 i i−1
(5.47)
∆t (∆y)2
Matrix form:
Aun+1 = un (5.48)
where A is a tridiagonal matrix with:

ai,i−1 = −r (5.49)
ai,i = 1 + 2r (5.50)
ai,i+1 = −r (5.51)

and r = ν∆t
(∆y)2
.

5.4 Detailed Procedure

5.4.1 Step 1: Case Setup
Create new case structure:
cd $FOAM_RUN
mkdir couetteFlow
cd couetteFlow
mkdir 0 constant system

5.4.2 Step 2: Create 1D Mesh

Create system/blockMeshDict for 1D domain:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object blockMeshDict ;
}
CFD Lab Manual - OpenFOAM 2406 48

convertToMeters 0.001; // Convert mm to m

vertices
(
(0 0 0) // 0 - bottom wall
(1 0 0) // 1
(1 10 0) // 2 - top wall
(0 10 0) // 3
(0 0 0.1) // 4
(1 0 0.1) // 5
(1 10 0.1) // 6
(0 10 0.1) // 7
);

blocks
(
hex (0 1 2 3 4 5 6 7) (1 100 1) simpleGrading (1 1 1)
);

edges () ;

boundary
(
bottomWall
{
type wall ;
faces
(
(0 1 5 4)
);
}

topWall
{
type wall ;
faces
(
(2 3 7 6)
);
}

sides
{
type empty ;
faces
(
(0 3 2 1)
(4 5 6 7)
(0 4 7 3)
(1 2 6 5)
);
}
);

mergePatchPairs () ;
CFD Lab Manual - OpenFOAM 2406 49

5.4.3 Step 3: Transport Properties

Create constant/transportProperties:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object tr a ns p o rt P ro p e rt i es ;
}

transportModel Newtonian ;
nu [0 2 -1 0 0 0 0] 1e -06; // Water viscosity

5.4.4 Step 4: Initial Conditions

Create 0/U:
FoamFile
{
version 2.0;
format ascii ;
class volVectorField ;
object U;
}

dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0) ;

boundaryField
{
bottomWall
{
type noSlip ;
}

topWall
{
type fixedValue ;
value uniform (1 0 0) ; // Moving wall velocity
}

sides
{
type empty ;
}
}

Create 0/p:
FoamFile
{
version 2.0;
format ascii ;
class volScalarField ;
object p;
}
CFD Lab Manual - OpenFOAM 2406 50

dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;

boundaryField
{
bottomWall
{
type zeroGradient ;
}

topWall
{
type zeroGradient ;
}

sides
{
type empty ;
}
}

5.4.5 Step 5: Custom Solver for Explicit Method

Create explicitCouetteSolver.C:
# include " fvCFD . H "

int main ( int argc , char * argv [])

{
# include " setRootCase . H "
# include " createTime . H "
# include " createMesh . H "

// Read transport properties

IOdictionary t r an s po r t Pr o pe r ti e s
(
IOobject
(
" t r an s po r t Pr o pe r ti e s " ,
runTime . constant () ,
mesh ,
IOobject :: MUST_READ ,
IOobject :: NO_WRITE
)
);

dimen sionedSc alar nu ( tr a ns p or t P ro p er t i es . lookup (" nu ") ) ;

// Create velocity field

volVectorField U
(
IOobject
(
"U",
runTime . timeName () ,
mesh ,
CFD Lab Manual - OpenFOAM 2406 51

IOobject :: MUST_READ ,
IOobject :: AUTO_WRITE
),
mesh
);

// Calculate stability criterion

scalar deltaY = Foam :: sqrt ( mesh . V () [0]) ; // Cell height
scalar maxDt = 0.5 * deltaY * deltaY / nu . value () ;

Info << " Maximum stable time step : " << maxDt << endl ;
Info << " Current time step : " << runTime . deltaT () . value () << endl ;

if ( runTime . deltaT () . value () > maxDt )

{
FatalError << " Time step too large for stability !" << exit (
FatalError ) ;
}

while ( runTime . loop () )

{
Info << " Time = " << runTime . timeName () << endl ;

// Explicit time stepping for diffusion

volVectorField UOld = U ;

// Calculate second derivative ( Laplacian )

volVectorField d2Udy2 = fvc :: laplacian ( nu , U ) ;

// Forward Euler update

U = UOld + runTime . deltaT () * d2Udy2 ;

// Apply boundary conditions

U . c o r r e c t B o u n d a r y C o n d i t i o n s () ;

runTime . write () ;

Info << " ExecutionTime = " << runTime . elapsedCpuTime () << " s "
<< " ClockTime = " << runTime . elapsedClockTime () << " s "
<< nl << endl ;
}

Info << " End \ n " << endl ;

return 0;
}

5.4.6 Step 6: Custom Solver for Implicit Method

Create implicitCouetteSolver.C:
# include " fvCFD . H "

int main ( int argc , char * argv [])

{
# include " setRootCase . H "
# include " createTime . H "
# include " createMesh . H "
CFD Lab Manual - OpenFOAM 2406 52

// Read transport properties

IOdictionary t r an s po r t Pr o pe r ti e s
(
IOobject
(
" t r an s po r t Pr o pe r ti e s " ,
runTime . constant () ,
mesh ,
IOobject :: MUST_READ ,
IOobject :: NO_WRITE
)
);

dimen sionedSc alar nu ( tr a ns p or t P ro p er t i es . lookup (" nu ") ) ;

// Create velocity field

volVectorField U
(
IOobject
(
"U",
runTime . timeName () ,
mesh ,
IOobject :: MUST_READ ,
IOobject :: AUTO_WRITE
),
mesh
);

while ( runTime . loop () )

{
Info << " Time = " << runTime . timeName () << endl ;

// Implicit time stepping using fvm operators

fvVectorMatrix UEqn
(
fvm :: ddt ( U ) - fvm :: laplacian ( nu , U )
);

UEqn . solve () ;

runTime . write () ;

Info << " ExecutionTime = " << runTime . elapsedCpuTime () << " s "
<< " ClockTime = " << runTime . elapsedClockTime () << " s "
<< nl << endl ;
}

Info << " End \ n " << endl ;

return 0;
}

5.4.7 Step 7: Solver Compilation

Create Make/files:
CFD Lab Manual - OpenFOAM 2406 53

explicitCouetteSolver .C
EXE = $ ( FOAM_USER_APPBIN ) / e x p l i c i t C o u e t t e S o l v e r

implicitCouetteSolver .C
EXE = $ ( FOAM_USER_APPBIN ) / i m p l i c i t C o u e t t e S o l v e r

Create Make/options:
EXE_INC = \
- I$ ( LIB_SRC ) / finiteVolume / lnInclude \
- I$ ( LIB_SRC ) / meshTools / lnInclude

EXE_LIBS = \
- lfiniteVolume \
- lmeshTools

Compile the solvers:

wmake e x p l i c i t C o u e t t e S o l v e r
wmake i m p l i c i t C o u e t t e S o l v e r

5.4.8 Step 8: Control Dictionary for Explicit Method

Create system/[Link]:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object controlDict ;
}

application explicitCouetteSolver ;
startFrom startTime ;
startTime 0;
stopAt endTime ;
endTime 0.1;
deltaT 1e -06; // Small time step for stability
writeControl runTime ;
writeInterval 0.01;
writeFormat ascii ;
writePrecision 8;
writeCompression off ;
timeFormat general ;
timePrecision 6;
runTi meModifi able false ;

5.4.9 Step 9: Control Dictionary for Implicit Method

Create system/[Link]:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
CFD Lab Manual - OpenFOAM 2406 54

object controlDict ;
}

application implicitCouetteSolver ;
startFrom startTime ;
startTime 0;
stopAt endTime ;
endTime 0.1;
deltaT 1e -03; // Larger time step allowed
writeControl runTime ;
writeInterval 0.01;
writeFormat ascii ;
writePrecision 8;
writeCompression off ;
timeFormat general ;
timePrecision 6;
runTi meModifi able false ;

5.4.10 Step 10: Numerical Schemes

Create system/fvSchemes:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object fvSchemes ;
}

ddtSchemes
{
default Euler ;
}

gradSchemes
{
default Gauss linear ;
}

divSchemes
{
default none ;
}

laplacianSchemes
{
default Gauss linear orthogonal ;
}

i n t e r p o l a t i on S c h e m e s
{
default linear ;
}

snGradSchemes
{
CFD Lab Manual - OpenFOAM 2406 55

default orthogonal ;
}

5.4.11 Step 11: Solution Control

Create system/fvSolution:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object fvSolution ;
}

solvers
{
U
{
solver smoothSolver ;
smoother symGaussSeidel ;
tolerance 1e -08;
relTol 0;
}
}

PISO
{
nCorrectors 2;
n N o n O r t h o g o n a l C o r r e c t o r s 0;
}

5.4.12 Step 12: Running Simulations

Execute both methods:
# Generate mesh
blockMesh

# Run explicit method

cp system / controlDict . explicit system / controlDict
e x p l i c i t C o u e t t e S o l v e r > log . explicit 2 >&1

# Save results
cp -r [0 -9]* results_explicit /

# Reset initial conditions

rm - rf [0 -9]*
cp -r 0. orig 0

# Run implicit method

cp system / controlDict . implicit system / controlDict
i m p l i c i t C o u e t t e S o l v e r > log . implicit 2 >&1

# Save results
cp -r [0 -9]* results_implicit /
CFD Lab Manual - OpenFOAM 2406 56

5.4.13 Step 13: Post-Processing and Validation

Create analytical solution script analytical [Link]:
import numpy as np
import matplotlib . pyplot as plt

def a na ly tic al _c oue tt e (y , t , h =0.01 , Uw =1.0 , nu =1 e -6 , n_terms =50) :

""" Analytical solution for unsteady Couette flow """
u = Uw * y / h # Steady - state part

# Transient part ( series solution )

for n in range (1 , n_terms + 1) :
u += (2 * Uw / ( n * np . pi ) ) * ( -1) ** n * \
np . sin ( n * np . pi * y / h ) * \
np . exp ( - n **2 * np . pi **2 * nu * t / h **2)

return u

# Create comparison plots

h = 0.01 # Channel height ( m )
Uw = 1.0 # Wall velocity ( m / s )
nu = 1e -6 # Kinematic viscosity ( $m ^2 $ / s )

y = np . linspace (0 , h , 100)
times = [0.001 , 0.01 , 0.05 , 0.1]

plt . figure ( figsize =(12 , 8) )

for i , t in enumerate ( times ) :
plt . subplot (2 , 2 , i +1)
u_analytical = an aly ti ca l_c ou et te (y , t , h , Uw , nu )

# Load OpenFOAM results ( would need actual data parsing )

# u_explicit = lo ad _o pen fo am _da ta (f ’ results_explicit /{ t }/ U ’)
# u_implicit = lo ad _o pen fo am _da ta (f ’ results_implicit /{ t }/ U ’)

plt . plot ( u_analytical , y , ’k - ’ , label = ’ Analytical ’ , linewidth =2)

# plt . plot ( u_explicit , y , ’r - - ’ , label = ’ Explicit ’)
# plt . plot ( u_implicit , y , ’b : ’ , label = ’ Implicit ’)

plt . xlabel ( ’ Velocity ( m / s ) ’)

plt . ylabel ( ’ Height ( m ) ’)
plt . title (f ’ Time = { t } s ’)
plt . legend ()
plt . grid ( True )

plt . tight_layout ()
plt . savefig ( ’ c oue tt e_ com pa ri son . png ’)
plt . show ()

# Calculate error norms

def calculate_errors () :
""" Calculate L2 and L$ \ infty$ error norms """
# Implementation would compare numerical vs analytical solutions
pass
CFD Lab Manual - OpenFOAM 2406 57

5.5 Expected Results

For typical parameters (h = 10 mm, Uw = 1 m/s, ν = 10−6 m2 /s):

• Explicit stability limit: ∆t < 5 × 10−5 s

• Implicit method: stable for larger time steps (∆t = 10−3 s)

• Steady-state time: t95% ≈ 0.1 s

• Both methods converge to linear velocity profile

• Explicit method: faster per time step, more time steps needed

• Implicit method: slower per time step, fewer time steps needed

5.6 Assignments
5.6.1 Assignment 1
Determine the maximum stable time step for explicit scheme and compare computational
efficiency with implicit scheme for the same accuracy level.

5.6.2 Assignment 2
Implement Crank-Nicolson scheme and compare the numerical diffusion effects with ex-
plicit and implicit methods.

5.6.3 Assignment 3
Extend the analysis to oscillatory Couette flow and investigate the phase lag between
wall motion and fluid response using different numerical schemes.
Experiment 6: Heat Conduction in 2D
Flat Plate

6.1 Objective
Solve the 2D heat conduction equation using explicit and implicit formulations with
laplacianFoam and custom solvers to analyze numerical stability and accuracy.

6.2 Theory
The 2D unsteady heat conduction equation:
 2
∂ 2T

∂T ∂ T
=α + (6.52)
∂t ∂x2 ∂y 2
where α = ρckp is thermal diffusivity.
For a rectangular plate with initial temperature T0 and boundary conditions:
T (0, y, t) = T1 (6.53)
T (L, y, t) = T2 (6.54)
T (x, 0, t) = T3 (6.55)
T (x, H, t) = T4 (6.56)
Analytical solution for steady state:
X∞  nπx   nπy 
Tsteady (x, y) = An sinh sin (6.57)
n=1
H H
For explicit finite difference in 2D:
α∆t α∆t
n+1
Ti,j n
= Ti,j + n
(Ti+1,j n
− 2Ti,j n
+ Ti−1,j )+ (T n − 2Ti,j
n n
+ Ti,j−1 ) (6.58)
(∆x) 2 (∆y)2 i,j+1
Stability condition for 2D explicit scheme:
α∆t α∆t 1
2
+ 2
≤ (6.59)
(∆x) (∆y) 2
For uniform grid (∆x = ∆y = h):
h2
∆t ≤ (6.60)
4α
Fourier number:
α∆t
Fo = (6.61)
h2

58
CFD Lab Manual - OpenFOAM 2406 59

2D Heat Conduction in Rectangular Plate

y
∂T/∂y = 0 (Insulated)

Heat Flow →
H = 5 cm

T₁ = 100°C T₂ = 20°C
Initial: T₀ = 20°C

∂T/∂y = 0 (Insulated) x
L = 10 cm

Heat Equation: Numerical Parameters:

∂T/∂t = α(∂²T/∂x² + ∂²T/∂y²) Grid: 50×25 cells
where α = k/(ρcₚ) = 9.75×10⁻⁵ m²/s Stability: Δt ≤ h²/(4α)
Fourier number: Fo = αΔt/h²
Material: Aluminum

Solution Methods:
1. Explicit Finite Difference (Custom Solver) | 2. Implicit Finite Difference (laplacianFoam) | 3. Analytical Solution (Fourier Series)

Figure 6.7: Schematic for the setup of the problem (already simplified)
CFD Lab Manual - OpenFOAM 2406 60

6.3 Detailed Procedure

6.3.1 Step 1: Case Setup
Create new case from laplacianFoam tutorial:
cd $FOAM_RUN
cp -r $FOAM_TUTORIALS / basic / laplacianFoam / flange .
mv flange heatConduction2D
cd heatConduction2D

If tutorial doesn’t exist, create new case:

cd $FOAM_RUN
mkdir heatConduction2D
cd heatConduction2D
mkdir 0 constant system

6.3.2 Step 2: Create 2D Rectangular Mesh

Create system/blockMeshDict:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object blockMeshDict ;
}

convertToMeters 0.01; // cm to m

vertices
(
(0 0 0) // 0
(10 0 0) // 1 - plate length 10 cm
(10 5 0) // 2 - plate height 5 cm
(0 5 0) // 3
(0 0 0.1) // 4 - thickness 1 mm
(10 0 0.1) // 5
(10 5 0.1) // 6
(0 5 0.1) // 7
);

blocks
(
hex (0 1 2 3 4 5 6 7) (50 25 1) simpleGrading (1 1 1)
);

edges () ;

boundary
(
left
{
type patch ;
faces
(
CFD Lab Manual - OpenFOAM 2406 61

(0 4 7 3)
);
}

right
{
type patch ;
faces
(
(1 2 6 5)
);
}

bottom
{
type patch ;
faces
(
(0 1 5 4)
);
}

top
{
type patch ;
faces
(
(2 3 7 6)
);
}

frontAndBack
{
type empty ;
faces
(
(0 3 2 1)
(4 5 6 7)
);
}
);

mergePatchPairs () ;

6.3.3 Step 3: Material Properties

Create constant/transportProperties:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object tr a ns p o rt P ro p e rt i es ;
}
CFD Lab Manual - OpenFOAM 2406 62

DT [0 2 -1 0 0 0 0] 1e -05; // Thermal diffusivity (

aluminum )

Alternative for custom solver with material properties: Create constant/thermophysicalPropertie

FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object thermophysicalProperties ;
}

// Aluminum properties
rho [1 -3 0 0 0 0 0] 2700;
Cp [0 2 -2 -1 0 0 0] 900;
k [1 1 -3 -1 0 0 0] 237;

// Calculated thermal diffusivity

alpha [0 2 -1 0 0 0 0] 9.75 e -05; // k /( rho * Cp )

6.3.4 Step 4: Initial and Boundary Conditions

Create 0/T:
FoamFile
{
version 2.0;
format ascii ;
class volScalarField ;
object T;
}

dimensions [0 0 0 1 0 0 0];
internalField uniform 293; // Initial temperature 20 C

boundaryField
{
left
{
type fixedValue ;
value uniform 373;
}

right
{
type fixedValue ;
value uniform 293;
}

bottom
{
type zeroGradient ; // Insulated
}

top
{
type zeroGradient ; // Insulated
CFD Lab Manual - OpenFOAM 2406 63

frontAndBack
{
type empty ;
}
}

6.3.5 Step 5: Custom Explicit Heat Solver

Create explicitHeatSolver.C:
# include " fvCFD . H "

int main ( int argc , char * argv [])

{
# include " setRootCase . H "
# include " createTime . H "
# include " createMesh . H "

// Read transport properties

IOdictionary t r an s po r t Pr o pe r ti e s
(
IOobject
(
" t r an s po r t Pr o pe r ti e s " ,
runTime . constant () ,
mesh ,
IOobject :: MUST_READ ,
IOobject :: NO_WRITE
)
);

dimen sionedSc alar DT ( tr a ns p or t P ro p er t i es . lookup (" DT ") ) ;

// Create temperature field

volScalarField T
(
IOobject
(
"T",
runTime . timeName () ,
mesh ,
IOobject :: MUST_READ ,
IOobject :: AUTO_WRITE
),
mesh
);

// Calculate stability criteria

scalar deltaX = Foam :: pow ( mesh . V () [0] , 1.0/3.0) ; // Approximate
cell size
scalar maxDt = deltaX * deltaX / (4.0 * DT . value () ) ;

Info << " Cell size ( approx ) : " << deltaX << " m " << endl ;
Info << " Maximum stable time step : " << maxDt << " s " << endl ;
CFD Lab Manual - OpenFOAM 2406 64

Info << " Current time step : " << runTime . deltaT () . value () << " s "
<< endl ;

if ( runTime . deltaT () . value () > maxDt )

{
WarningIn (" main ")
<< " Time step may be too large for stability !" << endl
<< " Consider reducing deltaT below " << maxDt << " s " <<
endl ;
}

// Create residual monitoring

scalar initialResidual = 1.0;
scalar tolerance = 1e -6;

while ( runTime . loop () )

{
Info << " Time = " << runTime . timeName () << endl ;

volScalarField TOld = T ;

// Explicit forward Euler for heat equation

volScalarField laplacianT = fvc :: laplacian ( DT , T ) ;
T = TOld + runTime . deltaT () * laplacianT ;

// Apply boundary conditions

T . c o r r e c t B o u n d a r y C o n d i t i o n s () ;

// Calculate residual
scalar residual = gMax ( mag ( T - TOld ) () . primitiveField () ) /
runTime . deltaT () . value () ;

if ( runTime . timeIndex () == 1) initialResidual = residual ;

scalar n orm al iz edR es id ual = residual / ( initialResidual + SMALL
);

Info << " Temperature residual : " << residual

<< " , Normalized : " << no rm ali ze dR esi du al << endl ;

// Monitor temperature statistics

scalar Tmin = gMin ( T . primitiveField () ) ;
scalar Tmax = gMax ( T . primitiveField () ) ;
scalar Tavg = gAverage ( T . primitiveField () ) ;

Info << " T stats - Min : " << Tmin << " K , Max : " << Tmax
<< " K , Avg : " << Tavg << " K " << endl ;

runTime . write () ;

Info << " ExecutionTime = " << runTime . elapsedCpuTime () << " s "
<< " ClockTime = " << runTime . elapsedClockTime () << " s "
<< nl << endl ;

// Check for convergence to steady state

if ( n orm al iz edR es id ua l < tolerance && runTime . timeIndex () >
100)
{
Info << " Converged to steady state !" << endl ;
CFD Lab Manual - OpenFOAM 2406 65

break ;
}
}

Info << " End \ n " << endl ;

return 0;
}

6.3.6 Step 6: Standard Implicit Heat Solver

For implicit method, use standard laplacianFoam with modifications.
Create system/[Link]:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object controlDict ;
}

application laplacianFoam ;
startFrom startTime ;
startTime 0;
stopAt endTime ;
endTime 1000;
deltaT 1; // Large time step for implicit
writeControl runTime ;
writeInterval 50;
writeFormat ascii ;
writePrecision 8;
writeCompression off ;
timeFormat general ;
timePrecision 6;
runTi meModifi able true ;

functions
{
temperatureStats
{
type fieldMinMax ;
libs (" l i b f i e l d F u n c t i o n O b j e c t s . so ") ;
writeControl timeStep ;
writeInterval 1;
fields (T);
}

probes
{
type probes ;
libs (" libsampling . so ") ;
writeControl timeStep ;
writeInterval 10;
fields (T);
probeLocations
(
(0.025 0.025 0.05) // Quarter point
CFD Lab Manual - OpenFOAM 2406 66

(0.05 0.025 0.05) // Center - bottom

(0.075 0.025 0.05) // Three - quarter point
);
}
}

6.3.7 Step 7: Control Dictionary for Explicit Method

Create system/[Link]:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object controlDict ;
}

application ex pl ic itH ea tS olv er ;

startFrom startTime ;
startTime 0;
stopAt endTime ;
endTime 1000;
deltaT 0.01; // Small time step for stability
writeControl runTime ;
writeInterval 50;
writeFormat ascii ;
writePrecision 8;
writeCompression off ;
timeFormat general ;
timePrecision 6;
runTi meModifi able false ;

functions
{
temperatureStats
{
type fieldMinMax ;
libs (" l i b f i e l d F u n c t i o n O b j e c t s . so ") ;
writeControl timeStep ;
writeInterval 100;
fields (T);
}
}

6.3.8 Step 8: Numerical Schemes

Create system/fvSchemes:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object fvSchemes ;
}
CFD Lab Manual - OpenFOAM 2406 67

ddtSchemes
{
default Euler ; // First - order in time
// Alternative : backward ; // Second - order in time
}

gradSchemes
{
default Gauss linear ;
}

divSchemes
{
default none ;
}

laplacianSchemes
{
default Gauss linear orthogonal ;
// Alternative : Gauss linear corrected ; // For non - orthogonal
meshes
}

i n t e r p o l a t i on S c h e m e s
{
default linear ;
}

snGradSchemes
{
default orthogonal ;
// Alternative : corrected ; // For non - orthogonal meshes
}

6.3.9 Step 9: Solution Control

Create system/fvSolution:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object fvSolution ;
}

solvers
{
T
{
solver GAMG ;
tolerance 1e -08;
relTol 0.01;
smoother GaussSeidel ;
nPreSweeps 0;
nPostSweeps 2;
CFD Lab Manual - OpenFOAM 2406 68

nFinestSweeps 2;
ca ch eA ggl om er ati on true ;
n C e l l s I n C o a r s e s t L e v e l 10;
agglomerator faceAreaPair ;
mergeLevels 1;
}

TFinal
{
$T ;
relTol 0;
}
}

SIMPLE
{
n N o n O r t h o g o n a l C o r r e c t o r s 0;
residualControl
{
T 1e -6;
}
}

relax ationFac tors

{
equations
{
T 1;
}
}

6.3.10 Step 10: Compilation and Execution

Compile custom solver:
# Create Make directory and files
mkdir Make

cat > Make / files << EOF

ex pl ic itH ea tS olv er . C
EXE = \ $ ( FOAM_USER_APPBIN ) / e xp lic it He atS ol ve r
EOF

cat > Make / options << EOF

EXE_INC = \\
-I \ $ ( LIB_SRC ) / finiteVolume / lnInclude \\
-I \ $ ( LIB_SRC ) / meshTools / lnInclude

EXE_LIBS = \\
- lfiniteVolume \\
- lmeshTools
EOF

# Compile
wmake

Run simulations:
CFD Lab Manual - OpenFOAM 2406 69

# Generate mesh
blockMesh

# Check mesh quality

checkMesh

# Run explicit method

cp system / controlDict . explicit system / controlDict
ex pl ic itH ea tS olv er > log . explicit 2 >&1

# Save results
mkdir results_explicit
cp -r [0 -9]* results_explicit /

# Reset and run implicit method

rm - rf [0 -9]*
cp -r 0. orig 0
cp system / controlDict . implicit system / controlDict
laplacianFoam > log . implicit 2 >&1

# Save results
mkdir results_implicit
cp -r [0 -9]* results_implicit /

6.3.11 Step 11: Analytical Solution for Validation

Create analytical [Link]:
import numpy as np
import matplotlib . pyplot as plt
from scipy . special import erf

def a na ly tic al _2 d_h ea t (x , y , t , L , H , T_left , T_right , T_init , alpha ,

n_terms =50) :
"""
Analytical solution for 2 D heat conduction with fixed boundaries
"""
# Steady - state solution ( linear in x - direction )
T_steady = T_init + ( T_left - T_init ) * ( L - x ) / L + ( T_right -
T_init ) * x / L

# Transient solution ( series )

T_transient = 0
for n in range (1 , n_terms + 1) :
for m in range (1 , n_terms + 1) :
lambda_nm = np . pi * np . sqrt (( n / L ) **2 + ( m / H ) **2)

A_nm = (16 * ( T_left - T_init ) ) / ( np . pi **2 * n * m ) * \

np . sin ( n * np . pi / 2) * np . sin ( m * np . pi / 2)

T_transient += A_nm * np . sin ( n * np . pi * x / L ) * \

np . sin ( m * np . pi * y / H ) * \
np . exp ( - lambda_nm **2 * alpha * t )

return T_steady + T_transient

CFD Lab Manual - OpenFOAM 2406 70

# Parameters
L = 0.1 # Length ( m )
H = 0.05 # Height ( m )
T_left = 373 # Hot boundary ( K )
T_right = 293 # Cold boundary ( K )
T_init = 293 # Initial temperature ( K )
alpha = 1e -5 # Thermal diffusivity ( $m ^2 $ / s )

# Create meshgrid
x = np . linspace (0 , L , 50)
y = np . linspace (0 , H , 25)
X , Y = np . meshgrid (x , y )

# Time points
times = [10 , 50 , 100 , 500]

fig , axes = plt . subplots (2 , 2 , figsize =(12 , 10) )

axes = axes . flatten ()

for i , t in enumerate ( times ) :

T_analytical = an aly ti ca l_2 d_ he at (X , Y , t , L , H , T_left , T_right ,
T_init , alpha )

im = axes [ i ]. contourf ( X *100 , Y *100 , T_analytical , levels =20 , cmap = ’

hot ’)
axes [ i ]. set_title (f ’ Time = { t } s ’)
axes [ i ]. set_xlabel ( ’ x ( cm ) ’)
axes [ i ]. set_ylabel ( ’ y ( cm ) ’)
plt . colorbar ( im , ax = axes [ i ])

plt . tight_layout ()
plt . savefig ( ’ a n a l y t i c a l _ t e m p e r a t u r e . png ’)
plt . show ()

# Compare centerline temperature

y_center = H / 2
T_centerline = an aly ti ca l_2 d_ he at (x , y_center , 100 , L , H , T_left ,
T_right , T_init , alpha )

plt . figure ( figsize =(10 , 6) )

plt . plot ( x *100 , T_centerline , ’k - ’ , linewidth =2 , label = ’ Analytical ’)
plt . xlabel ( ’ Position ( cm ) ’)
plt . ylabel ( ’ Temperature ( K ) ’)
plt . title ( ’ Centerline Temperature at t = 100 s ’)
plt . grid ( True )
plt . legend ()
plt . savefig ( ’ c e n t e r l i n e _ t e m p e r a t u r e . png ’)
plt . show ()

6.3.12 Step 12: Post-Processing and Comparison

Create comparison script compare [Link]:
import numpy as np
import matplotlib . pyplot as plt

def l o a d _ o p e n f o a m _ t e m p e r a t u r e ( case_dir , time_dir ) :

CFD Lab Manual - OpenFOAM 2406 71

"""
Load temperature field from OpenFOAM case
( This is a placeholder - actual implementation would parse OpenFOAM
files )
"""
# Use foamToVTK and then load VTK files , or
# Parse the OpenFOAM T file directly
pass

def plot_convergence () :
""" Plot convergence comparison between explicit and implicit
methods """

# Load residual data from log files

# This would parse the actual log files

time_explicit = np . linspace (0 , 1000 , 100000) # Many small time

steps
time_implicit = np . linspace (0 , 1000 , 1000) # Fewer large time
steps

# Mock residual data for demonstration

resid ual_expl icit = np . exp ( - time_explicit /100) + 0.01* np . random .
random ( len ( time_explicit ) )
resid ual_impl icit = np . exp ( - time_implicit /100) + 0.01* np . random .
random ( len ( time_implicit ) )

plt . figure ( figsize =(12 , 5) )

plt . subplot (1 , 2 , 1)
plt . semilogy ( time_explicit , residual_explicit , ’r - ’ , alpha =0.7 ,
label = ’ Explicit ’)
plt . semilogy ( time_implicit , residual_implicit , ’b - ’ , linewidth =2 ,
label = ’ Implicit ’)
plt . xlabel ( ’ Time ( s ) ’)
plt . ylabel ( ’ Residual ’)
plt . title ( ’ Convergence Comparison ’)
plt . legend ()
plt . grid ( True )

# Computational efficiency
plt . subplot (1 , 2 , 2)
e ff i ci e n cy _ ex p l ic i t = len ( time_explicit ) / 1000 # Relative cost
e ff i ci e n cy _ im p l ic i t = len ( time_implicit ) / 1000

methods = [ ’ Explicit ’ , ’ Implicit ’]

costs = [ efficiency_explicit , ef f ic i en c y _i m pl i c it ]

plt . bar ( methods , costs , color =[ ’ red ’ , ’ blue ’] , alpha =0.7)

plt . ylabel ( ’ Relative Computational Cost ’)
plt . title ( ’ Computational Efficiency ’)

plt . tight_layout ()
plt . savefig ( ’ method _compari son . png ’)
plt . show ()

# Run comparison
plot_convergence ()
CFD Lab Manual - OpenFOAM 2406 72

6.4 Expected Results

For aluminum plate (α = 1×10−5 m2 /s, 10×5 cm):

• Explicit stability limit: ∆t < 0.005 s (for 2 mm mesh)

• Implicit method: stable for ∆t = 1 s

• Steady-state time: ≈ 500 s (95% of final temperature)

• Maximum temperature gradient: near left boundary

• Both methods converge to same steady-state solution

• Explicit: 100,000 time steps for stability

• Implicit: 1,000 time steps with larger ∆t

6.5 Assignments
6.5.1 Assignment 1
Analyze heat conduction in a plate with sinusoidal boundary conditions and compare
numerical results with analytical Fourier series solution.

6.5.2 Assignment 2
Investigate the effect of mesh refinement on solution accuracy using Richardson extrap-
olation method. Calculate the order of accuracy for both schemes.

6.5.3 Assignment 3
Implement alternating direction implicit (ADI) method and compare its performance
with standard explicit and implicit schemes for large-scale problems.
Experiment 7: 1D Wave Propagation
in Still Lake

7.1 Objective
Simulate 1D wave propagation using the shallow water equations with interFoam and
analyze wave characteristics including dispersion and reflection.

7.2 Theory
The 1D shallow water equations:

∂h ∂(hu)
+ =0 (7.62)
 ∂t ∂x
∂(hu) ∂ 1
+ hu2 + gh2 = 0 (7.63)
∂t ∂x 2

where h is water depth, u is velocity, and g is gravitational acceleration.

Wave speed in shallow water:
(7.64)
p
c = gh
For sinusoidal waves in finite depth:

η(x, t) = A cos(kx − ωt + ϕ) (7.65)

Dispersion relation:
ω 2 = gk tanh(kh) (7.66)
√
For shallow water (kh ≪ 1): ω √= gh · k
For deep water (kh ≫ 1): ω = gk
Phase velocity: r
ω g
cp = = tanh(kh) (7.67)
k k
Group velocity:
r  
dω 1 g 2kh
cg = = tanh(kh) 1 + (7.68)
dk 2 k sinh(2kh)

Wave energy:
1
E = ρgA2 (7.69)
2

73
CFD Lab Manual - OpenFOAM 2406 74

Experiment 7: 1D Wave Propagation in Still Lake

Wave Parameters Mesh Parameters
Height: H = 0.1 m Cells: 400 × 40 × 1
Period:1stT Order
Stokes = 2.0 s Δx = 0.05 m Shallow Water
Length:
Wave λ = 6.25 m
Generation Δt = 0.01 s Wave Absorption
INLET OUTLET
ATMOSPHERE (Patch)

Wave Propagation
H=3m

d = 1.5 m
Gauge 1 Gauge
BOTTOM
2 (Wall) Gauge 3 Gauge 4

L = 20 m

Shallow Water Theory y

Wave Speed: c = √(gh)
Dispersion: ω² = gk tanh(kh) x
Phase Velocity: cₚ = ω/k g
Group Velocity: cₘ = dω/dk

Figure 7.8: Schematic of the problem setup (already simplified)

CFD Lab Manual - OpenFOAM 2406 75

7.3 Detailed Procedure

7.3.1 Step 1: Base Case Setup
Use damBreak tutorial as starting point:
cd $FOAM_RUN
cp -r $FOAM_TUTORIALS / multiphase / interFoam / laminar / damBreak .
mv damBreak w aveProap gation1D
cd waveP roapgati on1D

Clean up existing files:

rm - rf [0 -9]*
mkdir 0

7.3.2 Step 2: Create 1D Wave Tank Geometry

Create system/blockMeshDict:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object blockMeshDict ;
}

convertToMeters 1;

vertices
(
(0 0 0) // 0 - wave tank start
(20 0 0) // 1 - wave tank end (20 m long )
(20 1 0) // 2 - water surface
(0 1 0) // 3
(20 2 0) // 4 - air region top
(0 2 0) // 5
(0 0 0.1) // 6 - back face
(20 0 0.1) // 7
(20 1 0.1) // 8
(0 1 0.1) // 9
(20 2 0.1) // 10
(0 2 0.1) // 11
);

blocks
(
hex (0 1 2 3 6 7 8 9) water (400 20 1) simpleGrading (1 1 1)
// Water region
hex (3 2 4 5 9 8 10 11) air (400 20 1) simpleGrading (1 1 1) //
Air region
);

edges () ;

boundary
(
CFD Lab Manual - OpenFOAM 2406 76

inlet
{
type waveInlet ;
faces
(
(0 3 9 6) // Left boundary ( water )
(3 5 11 9) // Left boundary ( air )
);
}

outlet
{
type waveAbsorption ;
faces
(
(1 7 8 2) // Right boundary ( water )
(2 8 10 4) // Right boundary ( air )
);
}

bottom
{
type wall ;
faces
(
(0 6 7 1) // Tank bottom
);
}

atmosphere
{
type patch ;
faces
(
(5 4 10 11) // Top boundary
);
}

frontAndBack
{
type empty ;
faces
(
(0 1 2 3) // Front water
(3 2 4 5) // Front air
(6 9 8 7) // Back water
(9 11 10 8) // Back air
);
}
);

mergePatchPairs () ;

7.3.3 Step 3: Wave Generation Boundary Conditions

Create 0/[Link]:
CFD Lab Manual - OpenFOAM 2406 77

FoamFile
{
version 2.0;
format ascii ;
class volScalarField ;
object alpha . water ;
}

dimensions [0 0 0 0 0 0 0];
internalField uniform 0;

boundaryField
{
inlet
{
type waveAlpha ;
waveModel StokesFirst ;
nPaddle 1;
waveHeight 0.1; // Wave height ( m )
waveAngle 0.0; // Wave direction
wavePeriod 2.0; // Wave period ( s )
wavePhase 0.0; // Initial phase
seaLevel 1.0; // Still water level
s ea L ev e l As R ef e r en c e true ;
value uniform 0;
}

outlet
{
type waveAbsorption ;
absorptionType shallow ;
value uniform 0;
}

bottom
{
type zeroGradient ;
}

atmosphere
{
type inletOutlet ;
inletValue uniform 0;
value uniform 0;
}

frontAndBack
{
type empty ;
}
}

7.3.4 Step 4: Velocity Boundary Conditions

Create 0/U:
FoamFile
CFD Lab Manual - OpenFOAM 2406 78

{
version 2.0;
format ascii ;
class volVectorField ;
object U;
}

dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0) ;

boundaryField
{
inlet
{
type waveVelocity ;
waveModel StokesFirst ;
nPaddle 1;
waveHeight 0.1;
waveAngle 0.0;
wavePeriod 2.0;
wavePhase 0.0;
seaLevel 1.0;
value uniform (0 0 0) ;
}

outlet
{
type waveAbsorption ;
absorptionType shallow ;
value uniform (0 0 0) ;
}

bottom
{
type noSlip ;
}

atmosphere
{
type pressureInletOutletVelocity ;
value uniform (0 0 0) ;
}

frontAndBack
{
type empty ;
}
}

7.3.5 Step 5: Pressure Boundary Conditions

Create 0/p rgh:
FoamFile
{
version 2.0;
format ascii ;
CFD Lab Manual - OpenFOAM 2406 79

class volScalarField ;
object p rgh ;
}

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

boundaryField
{
inlet
{
type fix edFluxPr essure ;
value uniform 0;
}

outlet
{
type fix edFluxPr essure ;
value uniform 0;
}

bottom
{
type fix edFluxPr essure ;
value uniform 0;
}

atmosphere
{
type totalPressure ;
p0 uniform 0;
value uniform 0;
}

frontAndBack
{
type empty ;
}
}

7.3.6 Step 6: Initialize Water Level

Create system/setFieldsDict:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object setFieldsDict ;
}

de fa ul tFi el dV alu es
(
v ol S ca l a rF i el d V al u e alpha . water 0
);
CFD Lab Manual - OpenFOAM 2406 80

regions
(
boxToCell
{
box (0 0 0) (20 1 0.1) ; // Water region
fieldValues
(
v ol S c al a rF i e ld V al u e alpha . water 1
);
}
);

7.3.7 Step 7: Wave Properties

Create constant/waveProperties:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object waveProperties ;
}

inlet
{
waveModel StokesFirst ;

waveHeight 0.1; // Wave height ( m )

wavePeriod 2.0; // Wave period ( s )
waveLength 6.25; // Calculated from dispersion relation
wavePhase 0.0; // Initial phase ( rad )
waveAngle 0.0; // Wave propagation angle ( deg )

seaLevel 1.0; // Still water level ( m )

s ea L ev e l As R ef e r en c e true ;

// Shallow water parameters

depth 1.0; // Water depth ( m )
direction (1 0 0) ; // Wave direction

// Generation and absorption

genAbs 1; // Generate and absorb
nPaddle 1; // Number of paddles
}

outlet
{
absorptionType shallow ;
nPaddle 1;
}

7.3.8 Step 8: Transport Properties

Create constant/transportProperties:
CFD Lab Manual - OpenFOAM 2406 81

FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object tr a ns p o rt P ro p e rt i es ;
}

phases ( water air ) ;

water
{
transportModel Newtonian ;
nu [0 2 -1 0 0 0 0] 1e -06;
rho [1 -3 0 0 0 0 0] 1000;
}

air
{
transportModel Newtonian ;
nu [0 2 -1 0 0 0 0] 1.48 e -05;
rho [1 -3 0 0 0 0 0] 1;
}

sigma [1 0 -2 0 0 0 0] 0.07; // Surface tension

7.3.9 Step 9: Gravitational Properties

Create constant/g:
FoamFile
{
version 2.0;
format ascii ;
class uniformDimensionedVectorField ;
object g;
}

dimensions [0 1 -2 0 0 0 0];
value (0 -9.81 0) ;

7.3.10 Step 10: Control Dictionary with Wave Monitoring

Create system/controlDict:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object controlDict ;
}

application interFoam ;
startFrom startTime ;
CFD Lab Manual - OpenFOAM 2406 82

startTime 0;
stopAt endTime ;
endTime 20; // Simulate for 10 wave periods
deltaT 0.01; // Small time step for wave accuracy
writeControl runTime ;
writeInterval 0.1;
writeFormat ascii ;
writePrecision 8;
writeCompression off ;
timeFormat general ;
timePrecision 6;
runTi meModifi able true ;

// Courant number control

maxCo 0.5;
maxAlphaCo 0.5;
maxDeltaT 0.05;

functions
{
waveElevation
{
type surfaceElevation ;
libs (" libwaves . so ") ;
writeControl timeStep ;
writeInterval 1;
setFormat raw ;
i nt e rp o l at i on S c he m e cellPoint ;
fields ( alpha . water ) ;

sets
(
gauge1
{
type face ;
axis xyz ;
start (2 0 0.05) ;
end (2 2 0.05) ;
}
gauge2
{
type face ;
axis xyz ;
start (5 0 0.05) ;
end (5 2 0.05) ;
}
gauge3
{
type face ;
axis xyz ;
start (10 0 0.05) ;
end (10 2 0.05) ;
}
gauge4
{
type face ;
axis xyz ;
start (15 0 0.05) ;
CFD Lab Manual - OpenFOAM 2406 83

end (15 2 0.05) ;

}
);
}

forces
{
type forces ;
libs (" libforces . so ") ;
writeControl timeStep ;
writeInterval 10;
patches ( bottom ) ;
rho rhoInf ;
rhoInf 1000;
CofR (10 0 0) ;
}

probes
{
type probes ;
libs (" libsampling . so ") ;
writeControl timeStep ;
writeInterval 5;
fields ( U p rgh alpha . water ) ;
probeLocations
(
(2 0.5 0.05) // Near inlet
(5 0.5 0.05) // Quarter point
(10 0.5 0.05) // Center
(15 0.5 0.05) // Three - quarter point
);
}
}

7.3.11 Step 11: Numerical Schemes

Create system/fvSchemes:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object fvSchemes ;
}

ddtSchemes
{
default Euler ;
}

gradSchemes
{
default Gauss linear ;
}

divSchemes
CFD Lab Manual - OpenFOAM 2406 84

{
div ( rhoPhi , U ) Gauss linearUpwind grad ( U ) ;
div ( phi , alpha ) Gauss vanLeer ;
div ( phirb , alpha ) Gauss i n t e r f a c e C o m pr e s s i o n ;
div ((( rho * nuEff ) * dev2 ( T ( grad ( U ) ) ) ) ) Gauss linear ;
}

laplacianSchemes
{
default Gauss linear corrected ;
}

i n t e r p o l a t i on S c h e m e s
{
default linear ;
}

snGradSchemes
{
default corrected ;
}

7.3.12 Step 12: Solution Control

Create system/fvSolution:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object fvSolution ;
}

solvers
{
" alpha . water .*"
{
nAlphaCorr 2;
nAlphaSubCycles 1;
cAlpha 1;

MULESCorr yes ;
nLimiterIter 3;

solver smoothSolver ;
smoother symGaussSeidel ;
tolerance 1e -08;
relTol 0;
}

" pcorr .*"

{
solver GAMG ;
tolerance 1e -05;
relTol 0;
smoother GaussSeidel ;
CFD Lab Manual - OpenFOAM 2406 85

p rgh
{
solver GAMG ;
tolerance 1e -07;
relTol 0.05;
smoother GaussSeidel ;
}

p rghFinal
{
$p rgh ;
relTol 0;
}

U
{
solver smoothSolver ;
smoother symGaussSeidel ;
tolerance 1e -06;
relTol 0;
}
}

PIMPLE
{
momen tumPredi ctor no ;
nOuterCorrectors 1;
nCorrectors 3;
n N o n O r t h o g o n a l C o r r e c t o r s 0;
}

relax ationFac tors

{
equations
{
".*" 1;
}
}

7.3.13 Step 13: Running the Simulation

Execute the simulation:
# Generate mesh
blockMesh

# Check mesh quality

checkMesh

# Initialize water level

setFields

# Run simulation
interFoam > log . interFoam 2 >&1 &
CFD Lab Manual - OpenFOAM 2406 86

# Monitor progress
tail -f log . interFoam

# Monitor wave elevation

gnuplot - persist -e " plot ’ postProcessing / waveElevation /0/
surfaceElevation . dat ’ using 1:2 with lines "

7.3.14 Step 14: Wave Analysis and Post-Processing

Create wave analysis script :
import numpy as np
import matplotlib . pyplot as plt
from scipy . signal import find_peaks
from scipy . fft import fft , fftfreq

def load_wave_data ( filename ) :

""" Load wave elevation data from OpenFOAM output """
data = np . loadtxt ( filename , skiprows =1)
time = data [: , 0]
elevation = data [: , 1]
return time , elevation

def analyze_waves ( time , elevation , target_period =2.0) :

""" Analyze wave characteristics """

# Remove initial transient

mask = time > 5.0 # Remove first 5 seconds
time = time [ mask ]
elevation = elevation [ mask ]

# Calculate wave statistics

wave_height = np . max ( elevation ) - np . min ( elevation )
mean_level = np . mean ( elevation )

# Find peaks and troughs

peaks , _ = find_peaks ( elevation , distance = int ( len ( elevation ) /10) )
troughs , _ = find_peaks ( - elevation , distance = int ( len ( elevation ) /10)
)

# Calculate period from peaks

if len ( peaks ) > 1:
periods = np . diff ( time [ peaks ])
measured_period = np . mean ( periods )
else :
measured_period = 0

# Frequency analysis
dt = time [1] - time [0]
frequencies = fftfreq ( len ( elevation ) , dt )
fft_elevation = fft ( elevation )
power_spectrum = np . abs ( fft_elevation ) **2

# Find dominant frequency

po si ti ve_ fr eq _ma sk = frequencies > 0
domin ant_freq _idx = np . argmax ( power_spectrum [ po si ti ve_ fr eq _ma sk ])
CFD Lab Manual - OpenFOAM 2406 87

do mi na nt_ fr eq uen cy = frequencies [ p os iti ve _f req _m as k ][

domin ant_freq _idx ]
dominant_period = 1.0 / dom in an t_f re qu enc y

print ( f " Wave Analysis Results :")

print ( f " Wave Height : { wave_height :.3 f } m ")
print ( f " Mean Water Level : { mean_level :.3 f } m ")
print ( f " Measured Period ( peaks ) : { measured_period :.3 f } s ")
print ( f " Dominant Period ( FFT ) : { dominant_period :.3 f } s ")
print ( f " Target Period : { target_period :.3 f } s ")
print ( f " Period Error : { abs ( dominant_period - target_period ) /
target_period *100:.1 f }%")

return {
’ height ’: wave_height ,
’ period ’: dominant_period ,
’ frequency ’: dominant_frequency ,
’ mean_level ’: mean_level
}

def p lo t_ wav e_ an aly si s ( time , elevation , analysis_results ) :

""" Create comprehensive wave analysis plots """

fig , axes = plt . subplots (2 , 2 , figsize =(15 , 10) )

# Time series
axes [0 ,0]. plot ( time , elevation , ’b - ’ , linewidth =1)
axes [0 ,0]. axhline ( analysis_results [ ’ mean_level ’] , color = ’r ’ ,
linestyle = ’ - - ’ , label = ’ Mean Level ’)
axes [0 ,0]. set_xlabel ( ’ Time ( s ) ’)
axes [0 ,0]. set_ylabel ( ’ Wave Elevation ( m ) ’)
axes [0 ,0]. set_title ( ’ Wave Time Series ’)
axes [0 ,0]. grid ( True )
axes [0 ,0]. legend ()

# Frequency spectrum
dt = time [1] - time [0]
frequencies = fftfreq ( len ( elevation ) , dt )
fft_elevation = fft ( elevation )
power_spectrum = np . abs ( fft_elevation ) **2

positive_mask = frequencies > 0

axes [0 ,1]. loglog ( frequencies [ positive_mask ] , power_spectrum [
positive_mask ])
axes [0 ,1]. axvline ( analysis_results [ ’ frequency ’] , color = ’r ’ ,
linestyle = ’ - - ’ ,
label = f " f = { analysis_results [ ’ frequency ’]:.3 f }
Hz ")
axes [0 ,1]. set_xlabel ( ’ Frequency ( Hz ) ’)
axes [0 ,1]. set_ylabel ( ’ Power Spectral Density ’)
axes [0 ,1]. set_title ( ’ Frequency Spectrum ’)
axes [0 ,1]. grid ( True )
axes [0 ,1]. legend ()

# Phase analysis ( if multiple gauges available )

axes [1 ,0]. plot ( time , elevation , label = ’ Wave Elevation ’)
axes [1 ,0]. set_xlabel ( ’ Time ( s ) ’)
axes [1 ,0]. set_ylabel ( ’ Elevation ( m ) ’)
CFD Lab Manual - OpenFOAM 2406 88

axes [1 ,0]. set_title ( ’ Wave Propagation ’)

axes [1 ,0]. grid ( True )

# Wave statistics
axes [1 ,1]. text (0.1 , 0.8 , f " Wave Height : { analysis_results [ ’ height
’]:.3 f } m " ,
transform = axes [1 ,1]. transAxes , fontsize =12)
axes [1 ,1]. text (0.1 , 0.6 , f " Period : { analysis_results [ ’ period ’]:.3 f }
s",
transform = axes [1 ,1]. transAxes , fontsize =12)
axes [1 ,1]. text (0.1 , 0.4 , f " Frequency : { analysis_results [ ’ frequency
’]:.3 f } Hz " ,
transform = axes [1 ,1]. transAxes , fontsize =12)
axes [1 ,1]. text (0.1 , 0.2 , f " Mean Level : { analysis_results [ ’
mean_level ’]:.3 f } m " ,
transform = axes [1 ,1]. transAxes , fontsize =12)
axes [1 ,1]. set_title ( ’ Wave Statistics ’)
axes [1 ,1]. set_xlim (0 , 1)
axes [1 ,1]. set_ylim (0 , 1)
axes [1 ,1]. axis ( ’ off ’)

plt . tight_layout ()
plt . savefig ( ’ wave_analysis . png ’ , dpi =300 , bbox_inches = ’ tight ’)
plt . show ()

# Example usage ( would load actual OpenFOAM data )

if __name__ == " __main__ ":
# Generate synthetic data for demonstration
time = np . linspace (0 , 20 , 2000)
elevation = 0.05 * np . sin (2* np . pi * time /2.0) + 1.0 + 0.001* np . random
. randn ( len ( time ) )

results = analyze_waves ( time , elevation , target_period =2.0)

pl ot _w ave _a na lys is ( time , elevation , results )

7.4 Expected Results

For sinusoidal waves with H = 0.1 m, T = 2 s, h = 1 m:

• Wave celerity: c ≈ 3.13 m/s (shallow water approximation)

• Wavelength: λ ≈ 6.25 m

• Wave steepness: H/λ ≈ 0.016 (linear wave theory valid)

• Period preservation during propagation

• Minimal dispersion effects for shallow water

• Wave reflection coefficient at absorbing boundary < 5%

CFD Lab Manual - OpenFOAM 2406 89

7.5 Assignments
7.5.1 Assignment 1
Generate and analyze progressive waves of different frequencies. Validate dispersion re-
lation and compare with linear wave theory.

7.5.2 Assignment 2
Simulate wave reflection from a vertical wall and analyze standing wave patterns. Cal-
culate reflection coefficients.

7.5.3 Assignment 3
Investigate nonlinear effects by increasing wave amplitude and analyze higher-order har-
monics generation using spectral analysis.
Experiment 8: 2D Rayleigh-Taylor
Instability Using VOF Method

8.1 Objective
Simulate Rayleigh-Taylor instability using Volume of Fluid (VOF) method with interFoam
to study density-driven convection and interfacial instabilities.

8.2 Theory
Rayleigh-Taylor instability occurs when a denser fluid is positioned above a lighter fluid
in a gravitational field. The interface becomes unstable and develops into characteristic
bubble and spike structures.
Linear growth rate of Rayleigh-Taylor instability:
s
gk(ρh − ρl ) σk 3
γ= − (8.70)
ρh + ρl ρh + ρl

where k is wavenumber, ρh and ρl are heavy and light fluid densities, g is gravitational
acceleration, and σ is surface tension.
Atwood number:
ρh − ρl
A= (8.71)
ρh + ρl
For negligible surface tension, the growth rate simplifies to:

(8.72)
p
γ = Agk

Critical wavelength (surface tension stabilization):

r
σ
λc = 2π (8.73)
(ρh − ρl )g
Mixing layer thickness growth in nonlinear regime:

h(t) = αAgt2 (8.74)

where α ≈ 0.25 is the mixing parameter.

Bubble penetration depth:
hb (t) = αb Agt2 (8.75)
Spike penetration depth:
hs (t) = αs Agt2 (8.76)

90
CFD Lab Manual - OpenFOAM 2406 91

Typical values: αb ≈ 0.05, αs ≈ 0.8.

Figure 8.9: Typical representation of a Rayleigh - Taylor Instability (device a case ac-
cordingly)

8.3 Detailed Procedure

8.3.1 Step 1: Base Case Setup
Create new case from interFoam tutorial:
cd $FOAM_RUN
cp -r $FOAM_TUTORIALS / multiphase / interFoam / laminar / damBreak .
mv damBreak rayleighTaylor
cd rayleighTaylor
rm - rf [0 -9]*
mkdir 0

8.3.2 Step 2: Create 2D Rectangular Domain

Create system/blockMeshDict:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object blockMeshDict ;
}

convertToMeters 0.001; // mm to m

vertices
(
(0 0 0) // 0 - bottom left
(40 0 0) // 1 - bottom right (4 cm width )
(40 80 0) // 2 - top right (8 cm height )
(0 80 0) // 3 - top left
(0 0 1) // 4 - back bottom left
(40 0 1) // 5 - back bottom right
CFD Lab Manual - OpenFOAM 2406 92

(40 80 1) // 6 - back top right

(0 80 1) // 7 - back top left
);

blocks
(
hex (0 1 2 3 4 5 6 7) (80 160 1) simpleGrading (1 1 1)
);

edges () ;

boundary
(
left
{
type symmetryPlane ;
faces
(
(0 4 7 3)
);
}

right
{
type symmetryPlane ;
faces
(
(1 2 6 5)
);
}

bottom
{
type wall ;
faces
(
(0 1 5 4)
);
}

top
{
type wall ;
faces
(
(2 3 7 6)
);
}

frontAndBack
{
type empty ;
faces
(
(0 3 2 1)
(4 5 6 7)
);
}
CFD Lab Manual - OpenFOAM 2406 93

);

mergePatchPairs () ;

8.3.3 Step 3: Transport Properties with Density Ratio

Create constant/transportProperties:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object tr a ns p o rt P ro p e rt i es ;
}

phases ( heavy light ) ;

heavy
{
transportModel Newtonian ;
nu [0 2 -1 0 0 0 0] 1e -06; // Viscosity ( water -
like )
rho [1 -3 0 0 0 0 0] 1200; // Heavy fluid density
}

light
{
transportModel Newtonian ;
nu [0 2 -1 0 0 0 0] 1.5 e -05; // Viscosity ( air - like )
rho [1 -3 0 0 0 0 0] 1; // Light fluid density
}

sigma [1 0 -2 0 0 0 0] 0.01; // Reduced surface

tension

// Calculated Atwood number : A = (1200 -1) /(1200+1) $ \ approx$ 0.998

8.3.4 Step 4: Initial Interface Setup

Create 0/[Link]:
FoamFile
{
version 2.0;
format ascii ;
class volScalarField ;
object alpha . heavy ;
}

dimensions [0 0 0 0 0 0 0];
internalField uniform 0;

boundaryField
{
left
CFD Lab Manual - OpenFOAM 2406 94

{
type symmetryPlane ;
}

right
{
type symmetryPlane ;
}

bottom
{
type zeroGradient ;
}

top
{
type zeroGradient ;
}

frontAndBack
{
type empty ;
}
}

8.3.5 Step 5: Initial Velocity and Pressure

Create 0/U:
FoamFile
{
version 2.0;
format ascii ;
class volVectorField ;
object U;
}

dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0) ;

boundaryField
{
left
{
type symmetryPlane ;
}

right
{
type symmetryPlane ;
}

bottom
{
type noSlip ;
}
CFD Lab Manual - OpenFOAM 2406 95

top
{
type noSlip ;
}

frontAndBack
{
type empty ;
}
}

Create 0/p_rgh:
FoamFile
{
version 2.0;
format ascii ;
class volScalarField ;
object p rgh ;
}

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

boundaryField
{
left
{
type symmetryPlane ;
}

right
{
type symmetryPlane ;
}

bottom
{
type fix edFluxPr essure ;
value uniform 0;
}

top
{
type fix edFluxPr essure ;
value uniform 0;
}

frontAndBack
{
type empty ;
}
}

8.3.6 Step 6: Interface Initialization with Perturbation

Create system/setFieldsDict:
CFD Lab Manual - OpenFOAM 2406 96

FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object setFieldsDict ;
}

de fa ul tFi el dV alu es
(
v ol S ca l a rF i el d V al u e alpha . heavy 0 // Light fluid everywhere
initially
);

regions
(
// Main heavy fluid region ( upper half )
boxToCell
{
box (0 40 0) (40 80 1) ; // Upper half
fieldValues
(
v ol S c al a rF i e ld V al u e alpha . heavy 1
);
}

// Add small perturbation to interface

cylinderToCell
{
p1 (20 38 0) ; // Cylinder center below
interface
p2 (20 38 1) ;
radius 2; // Small perturbation radius
fieldValues
(
v ol S c al a rF i e ld V al u e alpha . heavy 1
);
}

// Multiple perturbations for 2 D instability

cylinderToCell
{
p1 (10 42 0) ;
p2 (10 42 1) ;
radius 1.5;
fieldValues
(
v ol S c al a rF i e ld V al u e alpha . heavy 0
);
}

cylinderToCell
{
p1 (30 38.5 0) ;
p2 (30 38.5 1) ;
radius 1.5;
fieldValues
CFD Lab Manual - OpenFOAM 2406 97

(
v ol S c al a rF i e ld V al u e alpha . heavy 1
);
}
);

8.3.7 Step 7: Gravitational Field

Create constant/g:
FoamFile
{
version 2.0;
format ascii ;
class uniformDimensionedVectorField ;
object g;
}

dimensions [0 1 -2 0 0 0 0];
value (0 -9.81 0) ;

8.3.8 Step 8: Control Dictionary with Monitoring

Create system/controlDict:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object controlDict ;
}

application interFoam ;
startFrom startTime ;
startTime 0;
stopAt endTime ;
endTime 2; // Simulate for 2 seconds
deltaT 1e -05; // Very small initial time step
writeControl adj ustableR unTime ;
writeInterval 0.05; // Write every 0.05 seconds
writeFormat ascii ;
writePrecision 8;
writeCompression off ;
timeFormat general ;
timePrecision 6;
runTi meModifi able true ;

// Adaptive time stepping

adjustTimeStep yes ;
maxCo 0.3;
maxAlphaCo 0.3;
maxDeltaT 0.001;

functions
{
CFD Lab Manual - OpenFOAM 2406 98

interfaceHeight
{
type interfaceHeight ;
libs (" l i b f i e l d F u n c t i o n O b j e c t s . so ") ;
writeControl timeStep ;
writeInterval 50;
locations ((0.01 0.04 0.0005) (0.02 0.04 0.0005) (0.03
0.04 0.0005) ) ;
alpha alpha . heavy ;
}

mixingHeight
{
type fieldMinMax ;
libs (" l i b f i e l d F u n c t i o n O b j e c t s . so ") ;
writeControl timeStep ;
writeInterval 50;
fields ( alpha . heavy ) ;
}

centerlineProbe
{
type probes ;
libs (" libsampling . so ") ;
writeControl timeStep ;
writeInterval 10;
fields ( alpha . heavy U p rgh ) ;
probeLocations
(
(0.02 0.01 0.0005) // Bottom
(0.02 0.02 0.0005) // Quarter
(0.02 0.04 0.0005) // Interface
(0.02 0.06 0.0005) // Three - quarter
(0.02 0.07 0.0005) // Top
);
}

vorticity
{
type vorticity ;
libs (" l i b f i e l d F u n c t i o n O b j e c t s . so ") ;
writeControl writeTime ;
writeInterval 1;
}

bubbleTracking
{
type surfaces ;
libs (" libsampling . so ") ;
writeControl writeTime ;
writeInterval 1;
surfaceFormat vtk ;
surfaces
(
interface
{
type isoSurface ;
isoField alpha . heavy ;
CFD Lab Manual - OpenFOAM 2406 99

isoValue 0.5;
interpolate true ;
}
);
}
}

8.3.9 Step 9: Numerical Schemes for Interface Capturing

Create system/fvSchemes:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object fvSchemes ;
}

ddtSchemes
{
default Euler ;
}

gradSchemes
{
default Gauss linear ;
}

divSchemes
{
div ( rhoPhi , U ) Gauss linearUpwind grad ( U ) ;
div ( phi , alpha ) Gauss vanLeer ;
div ( phirb , alpha ) Gauss i n t e r f a c e C o m pr e s s i o n ; // Interface
sharpening
div ((( rho * nuEff ) * dev2 ( T ( grad ( U ) ) ) ) ) Gauss linear ;
}

laplacianSchemes
{
default Gauss linear corrected ;
}

i n t e r p o l a t i on S c h e m e s
{
default linear ;
}

snGradSchemes
{
default corrected ;
}

fluxRequired
{
default no ;
p rgh ;
CFD Lab Manual - OpenFOAM 2406 100

pcorr ;
alpha . heavy ;
}

8.3.10 Step 10: Solution Control for Multiphase Flow

Create system/fvSolution:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object fvSolution ;
}

solvers
{
" alpha . heavy .*"
{
nAlphaCorr 2;
nAlphaSubCycles 1;
cAlpha 1; // Interface compression

MULESCorr yes ;
nLimiterIter 5;
al ph aA ppl yP re vCo rr yes ;

solver smoothSolver ;
smoother symGaussSeidel ;
tolerance 1e -08;
relTol 0;
minIter 1;
}

" pcorr .*"

{
solver GAMG ;
tolerance 1e -05;
relTol 0;
smoother GaussSeidel ;
nPreSweeps 0;
nPostSweeps 2;
nFinestSweeps 2;
ca ch eA ggl om er ati on true ;
n C e l l s I n C o a r s e s t L e v e l 10;
agglomerator faceAreaPair ;
mergeLevels 1;
}

p rgh
{
solver GAMG ;
tolerance 1e -08;
relTol 0.01;
smoother GaussSeidel ;
nPreSweeps 0;
CFD Lab Manual - OpenFOAM 2406 101

nPostSweeps 2;
nFinestSweeps 2;
ca ch eA ggl om er ati on true ;
n C e l l s I n C o a r s e s t L e v e l 10;
agglomerator faceAreaPair ;
mergeLevels 1;
}

p rghFinal
{
$p rgh ;
relTol 0;
}

U
{
solver smoothSolver ;
smoother symGaussSeidel ;
tolerance 1e -06;
relTol 0;
minIter 1;
}
}

PIMPLE
{
momen tumPredi ctor no ;
nOuterCorrectors 1;
nCorrectors 3;
n N o n O r t h o g o n a l C o r r e c t o r s 0;
pRefCell 0;
pRefValue 0;
}

relax ationFac tors

{
equations
{
".*" 1;
}
}

8.3.11 Step 11: Running the Simulation

Execute the simulation:
# Generate mesh
blockMesh

# Check mesh quality

checkMesh

# Initialize interface with perturbations

setFields

# Run simulation
interFoam > log . interFoam 2 >&1 &
CFD Lab Manual - OpenFOAM 2406 102

# Monitor progress
tail -f log . interFoam

# Real - time monitoring of interface evolution

watch -n 2 " tail -n 5 log . interFoam "

8.3.12 Step 12: Post-Processing and Analysis

Create analysis script rt [Link]:
import numpy as np
import matplotlib . pyplot as plt
from scipy . signal import find_peaks

def a n a l y z e _ b u b b l e _ p e n e t r a t i o n ( time_data , height_data ) :

""" Analyze bubble penetration according to h = $ \ alpha$Agt2 """

# Theoretical parameters
A = 0.998 # Atwood number
g = 9.81 # Gravity
alpha _theoret ical = 0.05 # Bubble growth coefficient

# Fit to h = $ \ alpha$ A g t2
def t he or eti ca l_ hei gh t (t , alpha ) :
return alpha * A * g * t **2

# Find best fit alpha

from scipy . optimize import curve_fit
popt , pcov = curve_fit ( theoretical_height , time_data , height_data )
alpha_fitted = popt [0]

print ( f " Atwood number : { A :.3 f }")

print ( f " Theoretical $ \ alpha$ : { alph a_theoret ical :.3 f }")
print ( f " Fitted $ \ alpha$ : { alpha_fitted :.3 f }")
print ( f " Error : { abs ( alpha_fitted - alpha_t heoretic al ) /
alpha _theoret ical *100:.1 f }%")

return alpha_fitted

def plo t_rt_evo lution ( time_data , bubble_data , spike_data ) :

""" Plot Rayleigh - Taylor instability evolution """

fig , axes = plt . subplots (2 , 2 , figsize =(15 , 10) )

# Bubble penetration
axes [0 ,0]. plot ( time_data , bubble_data , ’bo - ’ , label = ’ Simulation ’)

# Theoretical curve
A = 0.998
g = 9.81
alpha_b = 0.05
theory_bubble = alpha_b * A * g * time_data **2
axes [0 ,0]. plot ( time_data , theory_bubble , ’r - - ’ , label = ’ Theory ( $ \
alpha$ =0.05) ’)

axes [0 ,0]. set_xlabel ( ’ Time ( s ) ’)

CFD Lab Manual - OpenFOAM 2406 103

axes [0 ,0]. set_ylabel ( ’ Bubble Height ( m ) ’)

axes [0 ,0]. set_title ( ’ Bubble Penetration ’)
axes [0 ,0]. legend ()
axes [0 ,0]. grid ( True )

# Spike penetration
axes [0 ,1]. plot ( time_data , spike_data , ’ro - ’ , label = ’ Simulation ’)

alpha_s = 0.8
theory_spike = alpha_s * A * g * time_data **2
axes [0 ,1]. plot ( time_data , theory_spike , ’b - - ’ , label = ’ Theory ( $ \
alpha$ =0.8) ’)

axes [0 ,1]. set_xlabel ( ’ Time ( s ) ’)

axes [0 ,1]. set_ylabel ( ’ Spike Depth ( m ) ’)
axes [0 ,1]. set_title ( ’ Spike Penetration ’)
axes [0 ,1]. legend ()
axes [0 ,1]. grid ( True )

# Growth rates
if len ( time_data ) > 1:
dt = time_data [1] - time_data [0]
bu bb le _gr ow th _ra te = np . gradient ( bubble_data , dt )
spike _growth_ rate = np . gradient ( spike_data , dt )

axes [1 ,0]. plot ( time_data , bubble_growth_rate , ’b - ’ , label = ’

Bubble ’)
axes [1 ,0]. plot ( time_data , spike_growth_rate , ’r - ’ , label = ’ Spike
’)
axes [1 ,0]. set_xlabel ( ’ Time ( s ) ’)
axes [1 ,0]. set_ylabel ( ’ Growth Rate ( m / s ) ’)
axes [1 ,0]. set_title ( ’ Instantaneous Growth Rates ’)
axes [1 ,0]. legend ()
axes [1 ,0]. grid ( True )

# Mixing statistics
mixing_width = bubble_data + spike_data
axes [1 ,1]. plot ( time_data , mixing_width , ’g - ’ , linewidth =2)
axes [1 ,1]. set_xlabel ( ’ Time ( s ) ’)
axes [1 ,1]. set_ylabel ( ’ Mixing Width ( m ) ’)
axes [1 ,1]. set_title ( ’ Total Mixing Layer Width ’)
axes [1 ,1]. grid ( True )

plt . tight_layout ()
plt . savefig ( ’ rt_analysis . png ’ , dpi =300 , bbox_inches = ’ tight ’)
plt . show ()

def c a l c u l a t e _ g r o w t h _ p a r a m e t e r s () :
""" Calculate theoretical growth parameters """

# Material properties
rho_heavy = 1200 # kg / $m ^3 $
rho_light = 1 # kg / $m ^3 $
g = 9.81 # m / $s ^2 $
sigma = 0.01 # N/m

# Atwood number
A = ( rho_heavy - rho_light ) / ( rho_heavy + rho_light )
CFD Lab Manual - OpenFOAM 2406 104

# Critical wavelength
lambda_c = 2 * np . pi * np . sqrt ( sigma / (( rho_heavy - rho_light ) * g
))

# Growth rate for different wavelengths

wavelengths = np . logspace ( -3 , -1 , 100) # 1 mm to 10 cm
k = 2 * np . pi / wavelengths

# Without surface tension

gamma_no_st = np . sqrt ( A * g * k )

# With surface tension

gamma_with_st = np . sqrt ( A * g * k - sigma * k **3 / ( rho_heavy +
rho_light ) )
gamma_with_st = np . real ( gamma_with_st ) # Take only real part

plt . figure ( figsize =(10 , 6) )

plt . loglog ( wavelengths *1000 , gamma_no_st , ’b - ’ , label = ’ Without
surface tension ’)
plt . loglog ( wavelengths *1000 , gamma_with_st , ’r - ’ , label = ’ With
surface tension ’)
plt . axvline ( lambda_c *1000 , color = ’k ’ , linestyle = ’ - - ’ ,
label =f ’ Critical $ \ lambda$ = { lambda_c *1000:.1 f } mm ’)
plt . xlabel ( ’ Wavelength ( mm ) ’)
plt . ylabel ( ’ Growth Rate (1/ s ) ’)
plt . title ( ’ Rayleigh - Taylor Growth Rate vs Wavelength ’)
plt . legend ()
plt . grid ( True )
plt . savefig ( ’ g r o w th _ r a t e _ a n a l y s i s . png ’)
plt . show ()

print ( f " Atwood number : { A :.3 f }")

print ( f " Critical wavelength : { lambda_c *1000:.1 f } mm ")
print ( f " Most unstable wavelength : { lambda_c *1000/ np . sqrt (3) :.1 f } mm
")

# Run analysis
if __name__ == " __main__ ":
c a l c u l a t e _ g r o w t h _ p a r a m e t e r s ()

# Generate synthetic data for demonstration

time = np . linspace (0 , 2 , 100)
A = 0.998
g = 9.81
bubble_height = 0.05 * A * g * time **2 + 0.001* np . random . randn ( len (
time ) )
spike_depth = 0.8 * A * g * time **2 + 0.002* np . random . randn ( len (
time ) )

plot_ rt_evolu tion ( time , bubble_height , spike_depth )

8.4 Expected Results

For density ratio ρ_h/ρ_l = 1200 (A ≈ 0.998):
CFD Lab Manual - OpenFOAM 2406 105

• Initial linear growth phase: exponential development of perturbations

• Nonlinear regime: bubble rise and spike descent

• Bubble growth coefficient: α_b ≈ 0.05

• Spike growth coefficient: α_s ≈ 0.8

• Critical wavelength: λ_c ≈ 3 mm (for σ = 0.01 N/m)

• Mixing layer development with characteristic mushroom structures

• Kelvin-Helmholtz secondary instabilities on spike surfaces

8.5 Assignments
8.5.1 Assignment 1
Vary the Atwood number and analyze its effect on instability growth rate. Compare with
linear stability theory predictions.

8.5.2 Assignment 2
Study the effect of surface tension on small-scale instabilities and determine the critical
wavelength below which surface tension stabilizes the interface.

8.5.3 Assignment 3
Investigate the nonlinear regime and characterize the transition to turbulent mixing.
Analyze the mixing efficiency and entrainment rates.
Experiment 9: Lid-Driven Cavity Flow
Using Vorticity-Streamfunction

9.1 Objective
Analyze lid-driven cavity flow using vorticity-streamfunction formulation with custom
OpenFOAM solvers to study recirculation patterns and corner vortices.

9.2 Theory
The vorticity-streamfunction formulation transforms the Navier-Stokes equations into:
Vorticity definition:
∂v ∂u
ω =∇×U= − (9.77)
∂x ∂y
Streamfunction definition:
∂ψ
u= (9.78)
∂y
∂ψ
v=− (9.79)
∂x
Poisson equation for streamfunction:

∇2 ψ = −ω (9.80)

Vorticity transport equation:

∂ω ∂ω ∂ω
+u +v = ν∇2 ω (9.81)
∂t ∂x ∂y
For steady flow:
∂ω ∂ω
u +v = ν∇2 ω (9.82)
∂x ∂y
Boundary conditions for streamfunction: - On walls: ψ = constant, ∂ψ ∂n
= 0 (no-slip)
- Moving wall: ∂n = Uwall
∂ψ
2
Vorticity boundary conditions: - On stationary walls: ω = − ∂∂nψ2 - On moving walls:
ω = 2(Uwallh−uwall ) where h is grid spacing
Stream function circulation:
I ZZ
Γ = U · dl = ω dA (9.83)

106
CFD Lab Manual - OpenFOAM 2406 107

Primary vortex strength:

Ulid α
ωmax = Re (9.84)
L
where α ≈ 0.6 for moderate Reynolds numbers.

9.3 Detailed Procedure

Figure 9.10: A typical case of lid driven cavity

9.3.1 Step 1: Base Case Setup

Create cavity flow case:
cd $FOAM_RUN
cp -r $FOAM_TUTORIALS / incompressible / icoFoam / cavity .
mv cavity lidDrivenCavity
cd lidDrivenCavity

Or create new case:

cd $FOAM_RUN
mkdir lidDrivenCavity
cd lidDrivenCavity
mkdir 0 constant system
CFD Lab Manual - OpenFOAM 2406 108

9.3.2 Step 2: Create Square Cavity Mesh

Create system/blockMeshDict:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object blockMeshDict ;
}

convertToMeters 0.1; // 10 cm cavity

vertices
(
(0 0 0) // 0 - bottom left
(1 0 0) // 1 - bottom right
(1 1 0) // 2 - top right
(0 1 0) // 3 - top left
(0 0 0.1) // 4 - back bottom left
(1 0 0.1) // 5 - back bottom right
(1 1 0.1) // 6 - back top right
(0 1 0.1) // 7 - back top left
);

blocks
(
hex (0 1 2 3 4 5 6 7) (80 80 1) simpleGrading (1 1 1)
);

edges () ;

boundary
(
movingWall
{
type wall ;
faces
(
(3 7 6 2) // Top wall
);
}

fixedWalls
{
type wall ;
faces
(
(0 4 7 3) // Left wall
(1 2 6 5) // Right wall
(0 1 5 4) // Bottom wall
);
}

frontAndBack
{
type empty ;
faces
CFD Lab Manual - OpenFOAM 2406 109

(
(0 3 2 1) // Front
(4 5 6 7) // Back
);
}
);

mergePatchPairs () ;

9.3.3 Step 3: Custom Vorticity-Streamfunction Solver

Create vorticityStreamFunctionSolver.C:
# include " fvCFD . H "

int main ( int argc , char * argv [])

{
# include " setRootCase . H "
# include " createTime . H "
# include " createMesh . H "

// Read transport properties

IOdictionary t r an s po r t Pr o pe r ti e s
(
IOobject
(
" t r an s po r t Pr o pe r ti e s " ,
runTime . constant () ,
mesh ,
IOobject :: MUST_READ ,
IOobject :: NO_WRITE
)
);

dimen sionedSc alar nu ( tr a ns p or t P ro p er t i es . lookup (" nu ") ) ;

// Create fields
volScalarField psi
(
IOobject
(
" psi " ,
runTime . timeName () ,
mesh ,
IOobject :: MUST_READ ,
IOobject :: AUTO_WRITE
),
mesh
);

volScalarField omega
(
IOobject
(
" omega " ,
runTime . timeName () ,
mesh ,
CFD Lab Manual - OpenFOAM 2406 110

IOobject :: MUST_READ ,
IOobject :: AUTO_WRITE
),
mesh
);

volVectorField U
(
IOobject
(
"U",
runTime . timeName () ,
mesh ,
IOobject :: NO_READ ,
IOobject :: AUTO_WRITE
),
mesh ,
dimen sionedVe ctor (" U " , dimVelocity , vector :: zero )
);

// Convergence criteria
scalar tolerance = 1e -6;
label maxIterations = 1000;

Info << " Starting vorticity - streamfunction solution " << endl ;

for ( label iter = 0; iter < maxIterations ; iter ++)

{
Info << " Iteration " << iter + 1 << endl ;

volScalarField omegaOld = omega ;

volScalarField psiOld = psi ;

// Solve Poisson equation for streamfunction

fvScalarMatrix psiEqn
(
fvm :: laplacian ( psi ) + omega
);

psiEqn . solve () ;

// Calculate velocity from streamfunction

U . component (0) = fvc :: grad ( psi ) . component (1) ; // u = dpsi / dy
U . component (1) = - fvc :: grad ( psi ) . component (0) ; // v = - dpsi / dx
U . component (2) = 0;

// Apply velocity boundary conditions

U . c o r r e c t B o u n d a r y C o n d i t i o n s () ;

// Solve vorticity transport equation

fvScalarMatrix omegaEqn
(
fvm :: div ( phi , omega ) - fvm :: laplacian ( nu , omega )
);

// Calculate phi from velocity

su rf ac eSc al ar Fie ld phi = fvc :: interpolate ( U ) & mesh . Sf () ;
CFD Lab Manual - OpenFOAM 2406 111

omegaEqn . solve () ;

// Update vorticity boundary conditions

updat eVortici tyBC ( omega , U , psi , mesh ) ;

// Check convergence
scalar omegaRes = gMax ( mag ( omega - omegaOld ) ) ;
scalar psiRes = gMax ( mag ( psi - psiOld ) ) ;

Info << " Vorticity residual : " << omegaRes << endl ;
Info << " Streamfunction residual : " << psiRes << endl ;

if ( omegaRes < tolerance && psiRes < tolerance )

{
Info << " Converged in " << iter + 1 << " iterations " <<
endl ;
break ;
}
}

// Calculate derived quantities

volScalarField magU = mag ( U ) ;
volScalarField vorticityMag = mag ( omega ) ;

// Write final fields

runTime . write () ;

Info << " Solution completed " << endl ;

Info << " Maximum velocity : " << gMax ( magU ) << " m / s " << endl ;
Info << " Maximum vorticity : " << gMax ( vorticityMag ) << " 1/ s " <<
endl ;

return 0;
}

void update Vorticit yBC ( volScalarField & omega , const volVectorField & U ,
const volScalarField & psi , const fvMesh & mesh )
{
// Update vorticity boundary conditions
const fvPatchList & patches = mesh . boundary () ;

forAll ( patches , patchI )

{
const fvPatch & patch = patches [ patchI ];

if ( patch . name () == " movingWall ")

{
// Moving wall BC : omega = 2*( U_wall - U_fluid ) / dy
scalar Uwall = 1.0; // Moving wall velocity
scalar dy = mesh . deltaCoeffs () [ patchI ];

scalarField & omegaPatch = omega . boundaryField () [ patchI ];

const scalarField & UPatch = U . boundaryField () [ patchI ].
component (0) ;

omegaPatch = 2.0 * ( Uwall - UPatch ) * dy ;

}
else if ( patch . name () == " fixedWalls ")
CFD Lab Manual - OpenFOAM 2406 112

{
// Stationary wall BC : omega = - d2psi / dn2
const scalarField & psiPatch = psi . boundaryField () [ patchI ];
const scalarField & psiInternal = psi . internalField () ;

// Approximate second derivative

scalarField & omegaPatch = omega . boundaryField () [ patchI ];
// Implementation of second derivative calculation
// This is simplified - actual implementation would be more
complex
omegaPatch = 0.0;
}
}
}

9.3.4 Step 4: Initial and Boundary Conditions

Create 0/U:
FoamFile
{
version 2.0;
format ascii ;
class volVectorField ;
object U;
}

dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0) ;

boundaryField
{
movingWall
{
type fixedValue ;
value uniform (1 0 0) ; // Lid velocity
}

fixedWalls
{
type noSlip ;
}

frontAndBack
{
type empty ;
}
}

Create 0/psi:
FoamFile
{
version 2.0;
format ascii ;
class volScalarField ;
object psi ;
}
CFD Lab Manual - OpenFOAM 2406 113

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

boundaryField
{
movingWall
{
type fixedValue ;
value uniform 0;
}

fixedWalls
{
type fixedValue ;
value uniform 0;
}

frontAndBack
{
type empty ;
}
}

Create 0/omega:
FoamFile
{
version 2.0;
format ascii ;
class volScalarField ;
object omega ;
}

dimensions [0 0 -1 0 0 0 0];
internalField uniform 0;

boundaryField
{
movingWall
{
type calculated ;
value uniform 0;
}

fixedWalls
{
type calculated ;
value uniform 0;
}

frontAndBack
{
type empty ;
}
}
CFD Lab Manual - OpenFOAM 2406 114

9.3.5 Step 5: Compare with Primitive Variable Solution

Run standard icoFoam for comparison:
Create 0/p:
FoamFile
{
version 2.0;
format ascii ;
class volScalarField ;
object p;
}

dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;

boundaryField
{
movingWall
{
type zeroGradient ;
}

fixedWalls
{
type zeroGradient ;
}

frontAndBack
{
type empty ;
}
}

9.3.6 Step 6: Control Dictionary

Create system/controlDict:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object controlDict ;
}

application icoFoam ; // or v o r t i c i t y S t r e a m F u n c t i o n S o l v e r
startFrom startTime ;
startTime 0;
stopAt endTime ;
endTime 10;
deltaT 0.01;
writeControl runTime ;
writeInterval 0.5;
writeFormat ascii ;
writePrecision 8;
writeCompression off ;
CFD Lab Manual - OpenFOAM 2406 115

timeFormat general ;
timePrecision 6;
runTi meModifi able true ;

functions
{
vorticity
{
type vorticity ;
libs (" l i b f i e l d F u n c t i o n O b j e c t s . so ") ;
writeControl writeTime ;
writeInterval 1;
}

streamFunction
{
type streamFunction ;
libs (" l i b f i e l d F u n c t i o n O b j e c t s . so ") ;
writeControl writeTime ;
writeInterval 1;
}

forces
{
type forces ;
libs (" libforces . so ") ;
writeControl timeStep ;
writeInterval 50;
patches ( movingWall fixedWalls ) ;
rho rhoInf ;
rhoInf 1;
CofR (0.5 0.5 0) ;
}

probes
{
type probes ;
libs (" libsampling . so ") ;
writeControl timeStep ;
writeInterval 10;
fields ( U p vorticity ) ;
probeLocations
(
(0.5 0.5 0.05) // Center
(0.75 0.75 0.05) // Primary vortex core
(0.9 0.1 0.05) // Bottom right corner
(0.1 0.9 0.05) // Top left corner
);
}

centerlineData
{
type sets ;
libs (" libsampling . so ") ;
writeControl writeTime ;
writeInterval 1;
i nt e rp o l at i on S c he m e cellPoint ;
setFormat raw ;
CFD Lab Manual - OpenFOAM 2406 116

sets
(
h o r i z o n t a l C en t e r l i n e
{
type uniform ;
axis x;
start (0 0.5 0.05) ;
end (1 0.5 0.05) ;
nPoints 81;
}
ve rt ic alC en te rli ne
{
type uniform ;
axis y;
start (0.5 0 0.05) ;
end (0.5 1 0.05) ;
nPoints 81;
}
);
fields ( U vorticity ) ;
}
}

9.3.7 Step 7: Transport Properties

Create constant/transportProperties:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object tr a ns p o rt P ro p e rt i es ;
}

transportModel Newtonian ;
nu [0 2 -1 0 0 0 0] 0.001; // Re = 1000 for L =0.1 m , U =1 m /
s

9.3.8 Step 8: Numerical Schemes

Create system/fvSchemes:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object fvSchemes ;
}

ddtSchemes
{
default Euler ;
}
CFD Lab Manual - OpenFOAM 2406 117

gradSchemes
{
default Gauss linear ;
}

divSchemes
{
default none ;
div ( phi , U ) Gauss linearUpwind grad ( U ) ;
div ( phi , omega ) Gauss linearUpwind grad ( omega ) ;
}

laplacianSchemes
{
default Gauss linear orthogonal ;
}

i n t e r p o l a t i on S c h e m e s
{
default linear ;
}

snGradSchemes
{
default orthogonal ;
}

9.3.9 Step 9: Solution Control

Create system/fvSolution:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object fvSolution ;
}

solvers
{
p
{
solver GAMG ;
tolerance 1e -06;
relTol 0.05;
smoother GaussSeidel ;
}

pFinal
{
$p ;
relTol 0;
}

"( psi | omega ) "

{
CFD Lab Manual - OpenFOAM 2406 118

solver GAMG ;
tolerance 1e -08;
relTol 0.01;
smoother GaussSeidel ;
}

U
{
solver smoothSolver ;
smoother symGaussSeidel ;
tolerance 1e -05;
relTol 0.1;
}

UFinal
{
$U ;
relTol 0;
}
}

PISO
{
nCorrectors 2;
n N o n O r t h o g o n a l C o r r e c t o r s 0;
pRefCell 0;
pRefValue 0;
}

// For vorticity - streamfunction solver

SIMPLE
{
n N o n O r t h o g o n a l C o r r e c t o r s 0;
residualControl
{
psi 1e -6;
omega 1e -6;
}
}

relax ationFac tors

{
equations
{
".*" 0.9;
}
}

9.3.10 Step 10: Running Simulations

Execute both methods:
# Generate mesh
blockMesh

# Check mesh quality

checkMesh
CFD Lab Manual - OpenFOAM 2406 119

# Run primitive variable formulation

icoFoam > log . icoFoam 2 >&1

# Save results
mkdir res ults_pri mitive
cp -r [0 -9]* re sults_pr imitive /

# Reset and compile custom solver

rm - rf [0 -9]*
cp -r 0. orig 0

# Compile vorticity - streamfunction solver

wmake v o r t i c i t y S t r e a m F u n c t i o n S o l v e r

# Run vorticity - streamfunction formulation

v o r t i c i t y S t r e a m F u n c t i o n S o l v e r > log . vorticity 2 >&1

# Save results
mkdir res ults_vor ticity
cp -r [0 -9]* re sults_vo rticity /

9.3.11 Step 11: Benchmark Data Analysis

Create comparison script
import numpy as np
import matplotlib . pyplot as plt

def lo a d _b e nc h m ar k _d a ta () :
""" Load Ghia et al . (1982) benchmark data for Re =1000"""

# Horizontal centerline velocity

y_ghia = np . array ([0.0000 , 0.0547 , 0.0625 , 0.0703 , 0.1016 , 0.1719 ,
0.2813 , 0.4531 , 0.5000 , 0.6172 , 0.7344 , 0.8516 ,
0.9531 , 0.9609 , 0.9688 , 0.9766 , 1.0000])

u_ghia = np . array ([0.00000 , -0.08186 , -0.09266 , -0.10338 , -0.14612 ,

-0.24299 , -0.32726 , -0.17119 , -0.11477 , 0.02135 ,
0.16256 , 0.29093 , 0.55892 , 0.61756 , 0.68439 ,
0.75837 , 1.00000])

# Vertical centerline velocity

x_ghia = np . array ([0.0000 , 0.0625 , 0.0703 , 0.0781 , 0.0938 , 0.1563 ,
0.2266 , 0.2344 , 0.5000 , 0.8047 , 0.8594 , 0.9063 ,
0.9453 , 0.9531 , 0.9609 , 0.9688 , 1.0000])

v_ghia = np . array ([0.00000 , 0.1836 , 0.19713 , 0.20920 , 0.22965 ,

0.28124 , 0.30203 , 0.30174 , 0.05186 , -0.38598 ,
-0.44993 , -0.23827 , -0.22847 , -0.19254 ,
-0.15663 ,
-0.12146 , 0.00000])

return y_ghia , u_ghia , x_ghia , v_ghia

def an a l yz e _c a v it y _f l ow ( results_dir ) :
""" Analyze cavity flow results and compare with benchmark """
CFD Lab Manual - OpenFOAM 2406 120

# Load benchmark data

y_bench , u_bench , x_bench , v_bench = l oa d _b e n ch m ar k _ da t a ()

# Load OpenFOAM results ( placeholder - actual implementation would

parse files )
# This would load the centerline data from OpenFOAM output

plt . figure ( figsize =(15 , 5) )

# Horizontal centerline velocity

plt . subplot (1 , 3 , 1)
plt . plot ( u_bench , y_bench , ’ko - ’ , label = ’ Ghia et al . (1982) ’,
markersize =6)
# plt . plot ( u_openfoam , y_openfoam , ’r - ’ , label = ’ OpenFOAM ’)
plt . xlabel ( ’ u / U_lid ’)
plt . ylabel ( ’ y /L ’)
plt . title ( ’ Horizontal Centerline Velocity ’)
plt . legend ()
plt . grid ( True )

# Vertical centerline velocity

plt . subplot (1 , 3 , 2)
plt . plot ( x_bench , v_bench , ’ko - ’ , label = ’ Ghia et al . (1982) ’,
markersize =6)
# plt . plot ( x_openfoam , v_openfoam , ’r - ’ , label = ’ OpenFOAM ’)
plt . xlabel ( ’ x /L ’)
plt . ylabel ( ’ v / U_lid ’)
plt . title ( ’ Vertical Centerline Velocity ’)
plt . legend ()
plt . grid ( True )

# Vortex center locations for different Re

Re_values = [100 , 400 , 1000 , 3200 , 5000]
# Primary vortex centers from literature
x_primary = [0.6196 , 0.5608 , 0.5313 , 0.5165 , 0.5117]
y_primary = [0.7344 , 0.6054 , 0.5652 , 0.5477 , 0.5333]

plt . subplot (1 , 3 , 3)
plt . semilogx ( Re_values , x_primary , ’bo - ’ , label = ’x - center ’)
plt . semilogx ( Re_values , y_primary , ’ro - ’ , label = ’y - center ’)
plt . xlabel ( ’ Reynolds Number ’)
plt . ylabel ( ’ Vortex Center Location ’)
plt . title ( ’ Primary Vortex Center vs Re ’)
plt . legend ()
plt . grid ( True )

plt . tight_layout ()
plt . savefig ( ’ c a v i t y _ b e n c h m a r k _ c o m p a r i s o n . png ’)
plt . show ()

def p l o t _ s t r e a m l i n e s _ v o r t i c i t y () :
""" Plot streamlines and vorticity contours """

# Create sample data for visualization

x = np . linspace (0 , 1 , 50)
y = np . linspace (0 , 1 , 50)
X , Y = np . meshgrid (x , y )
CFD Lab Manual - OpenFOAM 2406 121

# Approximate analytical solution for visualization

psi = np . sin ( np . pi * X ) * np . sinh ( np . pi * Y ) / np . sinh ( np . pi )
omega = - np . pi **2 * np . sin ( np . pi * X ) * np . sinh ( np . pi * Y ) / np .
sinh ( np . pi )

fig , axes = plt . subplots (1 , 2 , figsize =(12 , 5) )

# Streamlines
axes [0]. contour (X , Y , psi , levels =20 , colors = ’ blue ’)
axes [0]. set_xlabel ( ’ x /L ’)
axes [0]. set_ylabel ( ’ y /L ’)
axes [0]. set_title ( ’ Streamlines ’)
axes [0]. set_aspect ( ’ equal ’)

# Vorticity contours
im = axes [1]. contourf (X , Y , omega , levels =20 , cmap = ’ RdBu_r ’)
axes [1]. set_xlabel ( ’ x /L ’)
axes [1]. set_ylabel ( ’ y /L ’)
axes [1]. set_title ( ’ Vorticity Contours ’)
axes [1]. set_aspect ( ’ equal ’)
plt . colorbar ( im , ax = axes [1])

plt . tight_layout ()
plt . savefig ( ’ c a v i t y _ f l o w _ v i s u a l i z a t i o n . png ’)
plt . show ()

# Run analysis
if __name__ == " __main__ ":
a na l yz e _ ca v it y _ fl o w ( ’ results_primitive ’)
p l o t _ s t r e a m l i n e s _ v o r t i c i t y ()

9.4 Expected Results

For Re = 1000 cavity flow:

• Primary vortex center: (0.531, 0.565) in normalized coordinates

• Secondary vortices in bottom corners for Re > 1000

• Maximum streamfunction: ψ_max ≈ 0.1 (normalized)

• Maximum vorticity: ω_max ≈ 3-5 U_lid/L at moving wall

• Excellent agreement with Ghia et al. benchmark data

• Vorticity-streamfunction method: faster convergence for steady flow

• Primitive variables: more memory efficient, easier boundary conditions

CFD Lab Manual - OpenFOAM 2406 122

Figure 9.11: Streamline contours for various Re values

9.5 Assignments
9.5.1 Assignment 1
Investigate the effect of Reynolds number (100, 1000, 5000) on flow patterns and identify
the formation of corner vortices.

9.5.2 Assignment 2
Compare vorticity-streamfunction results with primitive variable solutions using icoFoam.
Analyze computational efficiency and accuracy.

9.5.3 Assignment 3
Extend the analysis to 3D lid-driven cavity and study the development of three-dimensional
instabilities at high Reynolds numbers.
Experiment 10: Turbulence Models for
Pipe Flow

10.1 Objective
Compare different turbulence models (k-ϵ, k-ω, SST) for pipe flow analysis using simpleFoam
at various Reynolds numbers and validate against experimental correlations.

10.2 Theory
Reynolds-averaged momentum equation:
∂ Ūi ∂ Ūi 1 ∂ p̄ ∂ 2 Ūi ∂u′i u′j
+ Ūj =− +ν 2 − (10.85)
∂t ∂xj ρ ∂xi ∂xj ∂xj
The Reynolds stress tensor requires closure models:
 
∂ Ūi ∂ Ūj 2
−u′i u′j = νt + − kδij (10.86)
∂xj ∂xi 3

10.2.1 k-ϵ Model

Transport equations:
  
Dk ∂ νt ∂k
= Pk − ε + ν+ (10.87)
Dt ∂xj σk ∂xj
2
  
Dε ε ε ∂ νt ∂ε
= Cε1 Pk − Cε2 + ν+ (10.88)
Dt k k ∂xj σε ∂xj
2
Turbulent viscosity: νt = Cµ kε
Constants: Cµ = 0.09, Cε1 = 1.44, Cε2 = 1.92, σk = 1.0, σε = 1.3

10.2.2 k-ω Model

Transport equations:
  
Dk ∂ νt ∂k
∗
= Pk − β kω + ν+ (10.89)
Dt ∂xj σk ∂xj
  
Dω ω ∂ νt ∂ω
2
= γ Pk − βω + ν+ (10.90)
Dt k ∂xj σω ∂xj
Turbulent viscosity: νt = k
ω

123
CFD Lab Manual - OpenFOAM 2406 124

10.2.3 k-ω SST Model

Combines k-ω near walls and k-ϵ in free stream using blending function F1 :

ϕ = F1 ϕ1 + (1 − F1 )ϕ2 (10.91)

Cross-diffusion term:
1 ∂k ∂ω
CDkω = 2(1 − F1 )σω2 (10.92)
ω ∂xj ∂xj

10.2.4 Friction Factor Correlations

Blasius equation (smooth pipes, Re < 105 ):

f = 0.316Re−0.25 (10.93)

Colebrook-White equation (smooth pipes, all Re):

 
1 2.51
√ = −2 log10 √ (10.94)
f Re f
Velocity profile in wall coordinates:

u+ = y + (viscous sublayer, y + < 5) (10.95)

1
u+ = ln(y + ) + B (log layer, 30 < y + < 500) (10.96)
κ
where κ = 0.41 (Kármán constant) and B = 5.2.

10.3 Detailed Procedure

10.3.1 Step 1: Base Case Setup
Use pipe flow tutorial or create new case:
cd $FOAM_RUN
cp -r $FOAM_TUTORIALS / incompressible / simpleFoam / pitzDaily .
mv pitzDaily pipeFlow
cd pipeFlow

Or create new case:

cd $FOAM_RUN
mkdir pipeFlow
cd pipeFlow
mkdir 0 constant system

10.3.2 Step 2: Create Axisymmetric Pipe Geometry

Create system/blockMeshDict:
CFD Lab Manual - OpenFOAM 2406 125

FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object blockMeshDict ;
}

convertToMeters 0.001; // mm to m

// Pipe dimensions
diameter 20; // 20 mm diameter
radius 10; // 10 mm radius
length 500; // 500 mm length (25 diameters )

// Wedge geometry for axisymmetric simulation

angle 5; // 5 degree wedge
angleRad # calc " degToRad ( $angle ) ";
y1 # calc " $radius * sin ( $angleRad ) ";
z1 # calc " $radius * cos ( $angleRad ) ";

vertices
(
// Inlet face
(0 0 0) // 0 - centerline start
( $radius 0 0) // 1 - wall start
( $radius $y1 $z1 ) // 2 - wall start ( rotated )
(0 $y1 $z1 ) // 3 - centerline start ( rotated )

// Outlet face
( $length 0 0) // 4 - centerline end
(# calc " $length " $radius 0 0) // 5 - wall end
(# calc " $length " $radius $y1 $z1 ) // 6 - wall end ( rotated )
(# calc " $length " 0 $y1 $z1 ) // 7 - centerline end ( rotated )
);

blocks
(
hex (0 1 2 3 4 5 6 7) (250 20 1) simpleGrading (1 10 1) // Grading
toward wall
);

edges () ;

boundary
(
inlet
{
type patch ;
faces
(
(0 3 2 1)
);
}

outlet
{
CFD Lab Manual - OpenFOAM 2406 126

type patch ;
faces
(
(4 5 6 7)
);
}

wall
{
type wall ;
faces
(
(1 2 6 5)
);
}

centerline
{
type symmetryPlane ;
faces
(
(0 4 7 3)
);
}

wedgeFront
{
type wedge ;
faces
(
(0 1 5 4)
);
}

wedgeBack
{
type wedge ;
faces
(
(2 3 7 6)
);
}
);

mergePatchPairs () ;

10.3.3 Step 3: Velocity Boundary Conditions

Create 0/U:
FoamFile
{
version 2.0;
format ascii ;
class volVectorField ;
object U;
}
CFD Lab Manual - OpenFOAM 2406 127

dimensions [0 1 -1 0 0 0 0];
internalField uniform (5 0 0) ; // 5 m / s bulk velocity

boundaryField
{
inlet
{
type fixedValue ;
value uniform (5 0 0) ;
}

outlet
{
type zeroGradient ;
}

wall
{
type noSlip ;
}

centerline
{
type symmetryPlane ;
}

wedgeFront
{
type wedge ;
}

wedgeBack
{
type wedge ;
}
}

10.3.4 Step 4: Pressure Boundary Conditions

Create 0/p:
FoamFile
{
version 2.0;
format ascii ;
class volScalarField ;
object p;
}

dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;

boundaryField
{
inlet
{
CFD Lab Manual - OpenFOAM 2406 128

type zeroGradient ;
}

outlet
{
type fixedValue ;
value uniform 0;
}

wall
{
type zeroGradient ;
}

centerline
{
type symmetryPlane ;
}

wedgeFront
{
type wedge ;
}

wedgeBack
{
type wedge ;
}
}

10.3.5 Step 5: Turbulence Model Configurations

k-ϵ Model Setup
Create 0/k:
FoamFile
{
version 2.0;
format ascii ;
class volScalarField ;
object k;
}

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

// Calculate k = 1.5*( U * I ) ^2 , where I = 0.05 for turbulent pipe flow

// k = 1.5 * (5 * 0.05) ^2 = 0.09375
internalField uniform 0.09375;

boundaryField
{
inlet
{
type fixedValue ;
value uniform 0.09375;
}
CFD Lab Manual - OpenFOAM 2406 129

outlet
{
type zeroGradient ;
}

wall
{
type kqRWallFunction ;
value uniform 0.09375;
}

centerline
{
type symmetryPlane ;
}

wedgeFront
{
type wedge ;
}

wedgeBack
{
type wedge ;
}
}

Create 0/epsilon:
FoamFile
{
version 2.0;
format ascii ;
class volScalarField ;
object epsilon ;
}

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

// Calculate epsilon = C_mu ^0.75 * k ^1.5 / L

// L = 0.07 * h yd rau li c_ di ame te r = 0.07 * 0.02 = 0.0014 m
// epsilon = 0.09^0.75 * 0.09375^1.5 / 0.0014 = 0.206
internalField uniform 0.206;

boundaryField
{
inlet
{
type fixedValue ;
value uniform 0.206;
}

outlet
{
type zeroGradient ;
}

wall
CFD Lab Manual - OpenFOAM 2406 130

{
type ep s i lo n Wa l l Fu n ct i on ;
value uniform 0.206;
}

centerline
{
type symmetryPlane ;
}

wedgeFront
{
type wedge ;
}

wedgeBack
{
type wedge ;
}
}

k-ω Model Setup

Create 0/omega:
FoamFile
{
version 2.0;
format ascii ;
class volScalarField ;
object omega ;
}

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

// Calculate omega = k ^0.5 / ( C_mu ^0.25 * L )

// omega = 0.09375^0.5 / (0.09^0.25 * 0.0014) = 588
internalField uniform 588;

boundaryField
{
inlet
{
type fixedValue ;
value uniform 588;
}

outlet
{
type zeroGradient ;
}

wall
{
type ome gaWallFu nction ;
value uniform 588;
}
CFD Lab Manual - OpenFOAM 2406 131

centerline
{
type symmetryPlane ;
}

wedgeFront
{
type wedge ;
}

wedgeBack
{
type wedge ;
}
}

10.3.6 Step 6: Turbulence Model Selection

Create constant/[Link]:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object mome ntumTran sport ;
}

simulationType RAS ;

RAS
{
model kEpsilon ;
turbulence on ;
printCoeffs on ;

kEpsilonCoeffs
{
Cmu 0.09;
C1 1.44;
C2 1.92;
sigmaEps 1.3;
sigmaK 1.0;
}
}

Create constant/[Link]:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object mome ntumTran sport ;
}

simulationType RAS ;
CFD Lab Manual - OpenFOAM 2406 132

RAS
{
model kOmega ;
turbulence on ;
printCoeffs on ;

kOmegaCoeffs
{
betaStar 0.09;
gamma 5.0/9.0;
beta 0.075;
sigmaK 2.0;
sigmaOmega 2.0;
}
}

Create constant/[Link]:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object mome ntumTran sport ;
}

simulationType RAS ;

RAS
{
model kOmegaSST ;
turbulence on ;
printCoeffs on ;

kOmegaSSTCoeffs
{
alphaK1 0.85;
alphaK2 1.0;
alphaOmega1 0.5;
alphaOmega2 0.856;
beta1 0.075;
beta2 0.0828;
betaStar 0.09;
gamma1 5.0/9.0;
gamma2 0.44;
a1 0.31;
b1 1.0;
c1 10.0;
F3 no ;
}
}

10.3.7 Step 7: Transport Properties

Create constant/transportProperties:
CFD Lab Manual - OpenFOAM 2406 133

FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object tr a ns p o rt P ro p e rt i es ;
}

transportModel Newtonian ;

// For Re = 50 ,000 with D = 0.02 m , U = 5 m / s

// nu = U * D / Re = 5 * 0.02 / 50000 = 2e -06 $m ^2 $ / s
nu [0 2 -1 0 0 0 0] 2e -06;

10.3.8 Step 8: Control Dictionary with Performance Monitoring

Create system/controlDict:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object controlDict ;
}

application simpleFoam ;
startFrom startTime ;
startTime 0;
stopAt endTime ;
endTime 1000;
deltaT 1;
writeControl runTime ;
writeInterval 100;
writeFormat ascii ;
writePrecision 8;
writeCompression off ;
timeFormat general ;
timePrecision 6;
runTi meModifi able true ;

functions
{
wallShearStress
{
type wallShearStress ;
libs (" l i b f i e l d F u n c t i o n O b j e c t s . so ") ;
writeControl writeTime ;
writeInterval 1;
patches ( wall ) ;
}

yPlus
{
type yPlus ;
libs (" l i b f i e l d F u n c t i o n O b j e c t s . so ") ;
writeControl writeTime ;
CFD Lab Manual - OpenFOAM 2406 134

writeInterval 1;
patches ( wall ) ;
}

pressureDrop
{
type sur faceFiel dValue ;
libs (" l i b f i e l d F u n c t i o n O b j e c t s . so ") ;
writeControl timeStep ;
writeInterval 10;
fields (p);
operation areaAverage ;
regionType patch ;
name inlet ;
}

velocityProfile
{
type sets ;
libs (" libsampling . so ") ;
writeControl writeTime ;
writeInterval 1;
i nt e rp o l at i on S c he m e cellPoint ;
setFormat raw ;
sets
(
radialProfile
{
type uniform ;
axis y;
start (0.4 0 0) ; // 40 cm downstream
end (0.4 0.01 0) ; // Pipe radius
nPoints 50;
}
);
fields ( U k epsilon omega yPlus wallShearStress ) ;
}

frictionFactor
{
type coded ;
libs (" l i b u t i l i t y F u n c t i o n O b j e c t s . so ") ;
name frictionFactor ;
writeControl timeStep ;
writeInterval 50;

code
#{
// Calculate friction factor
const fvMesh & mesh = time () . lookupObject < fvMesh >(" region0 ")
;
const volScalarField & p = mesh . lookupObject < volScalarField
>(" p ") ;

// Get inlet and outlet patches

label inletID = mesh . boundaryMesh () . findPatchID (" inlet ") ;
label outletID = mesh . boundaryMesh () . findPatchID (" outlet ") ;
CFD Lab Manual - OpenFOAM 2406 135

// Calculate pressure drop

scalar pInlet = gSum ( p . boundaryField () [ inletID ] * mesh .
boundary () [ inletID ]. magSf () )
/ gSum ( mesh . boundary () [ inletID ]. magSf () ) ;
scalar pOutlet = gSum ( p . boundaryField () [ outletID ] * mesh .
boundary () [ outletID ]. magSf () )
/ gSum ( mesh . boundary () [ outletID ]. magSf () ) ;

scalar dp = pInlet - pOutlet ;

scalar L = 0.5; // Pipe length
scalar D = 0.02; // Pipe diameter
scalar rho = 1000; // Fluid density
scalar U = 5; // Bulk velocity

scalar f = 2 * dp * D / ( L * rho * U * U ) ;

Info << " Friction factor : " << f << endl ;

// Theoretical friction factor ( Blasius )

scalar Re = U * D / 2e -06;
scalar f_theory = 0.316 * pow ( Re , -0.25) ;

Info << " Reynolds number : " << Re << endl ;

Info << " Theoretical friction factor : " << f_theory << endl
;
Info << " Error : " << mag ( f - f_theory ) / f_theory * 100 <<
"%" << endl ;
#};
}
}

10.3.9 Step 9: Solution Control

Create system/fvSolution:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object fvSolution ;
}

solvers
{
p
{
solver GAMG ;
tolerance 1e -06;
relTol 0.1;
smoother GaussSeidel ;
}

"( U | k | epsilon | omega ) "

{
solver smoothSolver ;
smoother symGaussSeidel ;
CFD Lab Manual - OpenFOAM 2406 136

tolerance 1e -05;
relTol 0.1;
}
}

SIMPLE
{
n N o n O r t h o g o n a l C o r r e c t o r s 0;
consistent yes ;

residualControl
{
p 1e -4;
U 1e -5;
"( k | epsilon | omega ) " 1e -5;
}
}

relax ationFac tors

{
equations
{
U 0.9;
".*" 0.9;
}
}

10.3.10 Step 10: Running Multiple Turbulence Models

Create batch script [Link]:
#!/ bin / bash

# Generate mesh
blockMesh

# Array of turbulence models

models =(" kEpsilon " " kOmega " " kOmegaSST ")

for model in " $ { models [ @ ]}"; do

echo " Running $model model ..."

# Copy appropriate turbulence model file

cp constant / momentum Transport . $model constant / momen tumTransp ort

# Clean previous results

rm - rf [1 -9]*

# Run simulation
simpleFoam > log . $model 2 >&1

# Save results
mkdir results_$model
cp -r [0 -9]* results_$model /
cp log . $model results_$model /

echo " $model simulation completed "

CFD Lab Manual - OpenFOAM 2406 137

done

echo " All simulations completed "

Execute:
chmod + x runAllModels . sh
./ runAllModels . sh

10.3.11 Step 11: Reynolds Number Study

Create parametric study script [Link]:
import numpy as np
import matplotlib . pyplot as plt
import subprocess
import os

def r un _p ipe _f lo w_c as e ( Re , model =" kOmegaSST ") :

""" Run pipe flow case for given Reynolds number """

# Pipe parameters
D = 0.02 # Diameter ( m )
U = 5.0 # Velocity ( m / s )

# Calculate kinematic viscosity

nu = U * D / Re

# Update tr a ns p o rt P ro p e rt i es
with open ( ’ constant / transportProperties ’ , ’w ’) as f :
f . write (f ’ ’ ’ FoamFile
{{
version 2.0;
format ascii ;
class dictionary ;
object tr a ns p o rt P ro p e rt i es ;
}}

transportModel Newtonian ;
nu [0 2 -1 0 0 0 0] { nu :.2 e };
’ ’ ’)

# Copy turbulence model

subprocess . run ([ ’ cp ’ , f ’ constant / mome ntumTran sport .{ model } ’ ,
’ constant / momentumTransport ’])

# Clean case
subprocess . run ([ ’ rm ’ , ’-rf ’] + [ str ( i ) for i in range (1 , 10) ])

# Run simulation
result = subprocess . run ([ ’ simpleFoam ’] , capture_output = True , text =
True )

# Extract friction factor from log

log_lines = result . stdout . split ( ’\n ’)
friction_factor = None
for line in log_lines :
if ’ Friction factor : ’ in line :
CFD Lab Manual - OpenFOAM 2406 138

friction_factor = float ( line . split ( ’: ’) [ -1]. strip () )

break

return friction_factor

def t h e o r e t i c a l _ f r i c t i o n _ f a c t o r ( Re ) :
""" Calculate theoretical friction factor """
if Re < 2300:
return 64 / Re # Laminar flow
elif Re < 1 e5 :
return 0.316 * Re **( -0.25) # Blasius equation
else :
# Colebrook - White ( approximation )
return 0.184 * Re **( -0.2)

def c o m p a r e _ t u r b u l e n c e _ m o d e l s () :
""" Compare different turbulence models """

Re_values = np . array ([1 e4 , 2 e4 , 5 e4 , 1 e5 , 2 e5 ])

models = [ ’ kEpsilon ’ , ’ kOmega ’ , ’ kOmegaSST ’]

plt . figure ( figsize =(12 , 8) )

# Theoretical curves
Re_theory = np . logspace (3 , 6 , 100)
f_laminar = 64 / Re_theory
f_blasius = 0.316 * Re_theory **( -0.25)
f_turbulent = 0.184 * Re_theory **( -0.2)

plt . loglog ( Re_theory , f_laminar , ’k - - ’ , label = ’ Laminar (64/ Re ) ’)

plt . loglog ( Re_theory , f_blasius , ’k - ’ , label = ’ Blasius (0.316/ Re
^0.25) ’)
plt . loglog ( Re_theory , f_turbulent , ’k : ’ , label = ’ Turbulent (0.184/ Re
^0.2) ’)

# Simulation results ( placeholder - would run actual simulations )

colors = [ ’ red ’ , ’ blue ’ , ’ green ’]
markers = [ ’o ’ , ’s ’ , ’^ ’]

for i , model in enumerate ( models ) :

# Placeholder data - actual implementation would run
simulations
f_simulated = t h e o r e t i c a l _ f r i c t i o n _ f a c t o r ( Re_values ) * (1 +
0.1* np . random . randn ( len ( Re_values ) ) )
plt . loglog ( Re_values , f_simulated , markers [ i ] , color = colors [ i ] ,
markersize =8 , label =f ’{ model } model ’)

plt . xlabel ( ’ Reynolds Number ’)

plt . ylabel ( ’ Friction Factor ’)
plt . title ( ’ Friction Factor vs Reynolds Number - Turbulence Model
Comparison ’)
plt . legend ()
plt . grid ( True , alpha =0.3)
plt . xlim (1 e3 , 1 e6 )
plt . ylim (1 e -3 , 1e -1)

plt . tight_layout ()
plt . savefig ( ’ f r i c t i o n _ f a c t o r _ c o m p a r i s o n . png ’)
CFD Lab Manual - OpenFOAM 2406 139

plt . show ()

def p l o t _ v e l o c i t y _ p r o f i l e s () :
""" Plot velocity profiles for different turbulence models """

# Load velocity profile data ( placeholder )

r_normalized = np . linspace (0 , 1 , 50)

# Theoretical profiles
# Law of the wall : u + = (1/ kappa ) * ln ( y +) + B
# Power law : u / u_centerline = (1 - r / R ) ^(1/ n )

u_power_law = (1 - r_normalized ) **(1/7) # 1/7 power law

u_log_law = np . ones_like ( r_normalized ) # Simplified log law

plt . figure ( figsize =(12 , 5) )

# Velocity profile
plt . subplot (1 , 2 , 1)
plt . plot ( u_power_law , r_normalized , ’k - ’ , linewidth =2 , label = ’1/7
Power Law ’)
plt . plot ( u_log_law , r_normalized , ’k - - ’ , label = ’ Log Law ’)

# Placeholder simulation data

models = [ ’k - $ \ epsilon$ ’ , ’k - $ \ omega$ ’ , ’k - $ \ omega$ SST ’]
colors = [ ’ red ’ , ’ blue ’ , ’ green ’]
for i , model in enumerate ( models ) :
u_sim = u_power_law + 0.05* np . random . randn ( len ( r_normalized ) )
plt . plot ( u_sim , r_normalized , colors [ i ] , label =f ’{ model } model
’)

plt . xlabel ( ’ u / u_centerline ’)

plt . ylabel ( ’ r /R ’)
plt . title ( ’ Velocity Profile Comparison ’)
plt . legend ()
plt . grid ( True , alpha =0.3)

# Wall shear stress

plt . subplot (1 , 2 , 2)
y_plus = np . logspace (0 , 3 , 100)
u_plus_viscous = y_plus
u_plus_log = (1/0.41) * np . log ( y_plus ) + 5.2

plt . semilogx ( y_plus , u_plus_viscous , ’k - ’ , label = ’ Viscous sublayer

’)
plt . semilogx ( y_plus , u_plus_log , ’k - - ’ , label = ’ Log layer ’)

for i , model in enumerate ( models ) :

# Placeholder data
u_plus_sim = np . where ( y_plus < 11 , y_plus , (1/0.41) * np . log (
y_plus ) + 5.2)
u_plus_sim += 0.5* np . random . randn ( len ( y_plus ) )
plt . semilogx ( y_plus , u_plus_sim , colors [ i ] , alpha =0.7 , label =f
’{ model } model ’)

plt . xlabel ( ’ y + ’)
plt . ylabel ( ’ u + ’)
plt . title ( ’ Law of the Wall ’)
CFD Lab Manual - OpenFOAM 2406 140

plt . legend ()
plt . grid ( True , alpha =0.3)
plt . xlim (1 , 1000)
plt . ylim (0 , 30)

plt . tight_layout ()
plt . savefig ( ’ veloci ty_profi les . png ’)
plt . show ()

# Run analysis
if __name__ == " __main__ ":
c o m p a r e _ t u r b u l e n c e _ m o d e l s ()
p l o t _ v e l o c i t y _ p r o f i l e s ()

10.4 Expected Results

For pipe flow at Re = 50,000:

• k-ϵ model: Good for far-field, overpredicts near-wall turbulence

• k-ω model: Excellent near-wall behavior, sensitive to freestream conditions

• k-ω SST: Best overall performance, combines advantages of both

• Friction factor: f ≈ 0.0055 (Blasius: f = 0.0056)

• y+ values: 30-100 (appropriate for wall functions)

• Velocity profile: Close agreement with 1/7 power law

• Computational cost: k-ϵ < k-ω SST < k-ω

10.5 Assignments
10.5.1 Assignment 1
Compare predicted friction factors with Moody chart correlations for Reynolds numbers
ranging from 104 to 106 .

10.5.2 Assignment 2
Analyze near-wall treatment using different y+ values and evaluate the performance of
wall functions vs low-Re models.

10.5.3 Assignment 3
Investigate turbulence model performance for transitional flow regime (2000 < Re <
4000) and compare with DNS/experimental data.
Experiment 11: Convection-Diffusion
Transport Equation

11.1 Objective
Simulate generic convection-diffusion transport with forced and natural convection over
flat plate and in pipe using scalarTransportFoam and analyze numerical schemes.

11.2 Theory
General transport equation for scalar quantity ϕ:
∂ϕ
+ ∇ · (Uϕ) = ∇ · (Γ∇ϕ) + Sϕ (11.97)
∂t
where: - ϕ = transported scalar (concentration, temperature, etc.) - U = velocity
field - Γ = diffusion coefficient - Sϕ = source term
Peclet number (ratio of convection to diffusion):
UL
Pe = = Re · Sc (11.98)
D
where Sc = ν/D is Schmidt number (momentum to mass diffusivity ratio).
For heat transfer: P r = ν/α (Prandtl number)

11.2.1 Forced Convection over Flat Plate

Boundary layer equations:
∂u ∂u ∂ 2u
u +v =ν 2 (11.99)
∂x ∂y ∂y
∂T ∂T ∂ 2T
u +v =α 2 (11.100)
∂x ∂y ∂y
Local Nusselt number:
hx x
N ux = = 0.332Rex1/2 P r1/3 (laminar) (11.101)
k
Average Nusselt number:
h̄L 1/2
N uL = = 0.664ReL P r1/3 (11.102)
k
Thermal boundary layer thickness:
δT
= P r−1/3 (11.103)
δ
141
CFD Lab Manual - OpenFOAM 2406 142

11.2.2 Natural Convection

Rayleigh number:
gβ∆T L3
Ra = (11.104)
να
Grashof number:
gβ∆T L3
Gr = (11.105)
ν2
Local Nusselt number for vertical plate:
 1/4
Pr
N ux = 0.508Ra1/4
x (laminar) (11.106)
0.952 + P r

Churchill-Chu correlation (full range):

" #2
1/6
0.387RaL
N uL = 0.825 + (11.107)
[1 + (0.492/P r)9/16 ]8/27

11.2.3 Numerical Schemes for Convection-Diffusion

Upwind differencing (stable but diffusive):

ϕf = ϕP when U · n > 0 (11.108)

Central differencing (accurate but can be unstable):

ϕP + ϕN
ϕf = (11.109)
2
High-resolution schemes (TVD): - Linear upwind - QUICK - MUSCL

11.3 Detailed Procedure

11.3.1 Step 1: Forced Convection over Flat Plate
Case Setup

cd $FOAM_RUN
mkdir c on v ec t io n D if f us i o n
cd co n ve c ti o n Di f fu s i on
mkdir flatPlate pipeFlow
cd flatPlate
mkdir 0 constant system

Flat Plate Geometry

Create system/blockMeshDict:
CFD Lab Manual - OpenFOAM 2406 143

FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object blockMeshDict ;
}

convertToMeters 0.01; // cm to m

vertices
(
// Leading edge region
(0 0 0) // 0 - plate start
(0 2 0) // 1 - far field top
( -2 2 0) // 2 - upstream top
( -2 0 0) // 3 - upstream start

// Main plate region

(10 0 0) // 4 - plate end (10 cm )
(10 2 0) // 5 - downstream top

// Z - direction (2 D case )
(0 0 0.1) // 6
(0 2 0.1) // 7
( -2 2 0.1) // 8
( -2 0 0.1) // 9
(10 0 0.1) // 10
(10 2 0.1) // 11
);

blocks
(
// Upstream block
hex (3 0 1 2 9 6 7 8) upstream (20 40 1) simpleGrading (1 10 1)

// Main plate block

hex (0 4 5 1 6 10 11 7) main (100 40 1) simpleGrading (1 10 1)
);

edges () ;

boundary
(
inlet
{
type patch ;
faces
(
(3 2 8 9)
);
}

outlet
{
type patch ;
faces
CFD Lab Manual - OpenFOAM 2406 144

(
(4 5 11 10)
);
}

top
{
type patch ;
faces
(
(2 1 7 8)
(1 5 11 7)
);
}

plate
{
type wall ;
faces
(
(0 4 10 6)
);
}

upstream
{
type symmetryPlane ;
faces
(
(3 0 6 9)
);
}

frontAndBack
{
type empty ;
faces
(
(3 9 8 2)
(2 8 7 1)
(0 1 7 6)
(0 6 10 4)
(4 10 11 5)
(1 5 11 7)
);
}
);

mergePatchPairs () ;

Velocity Field
Create 0/U:
FoamFile
{
version 2.0;
CFD Lab Manual - OpenFOAM 2406 145

format ascii ;
class volVectorField ;
object U;
}

dimensions [0 1 -1 0 0 0 0];
internalField uniform (2 0 0) ; // 2 m / s freestream

boundaryField
{
inlet
{
type fixedValue ;
value uniform (2 0 0) ;
}

outlet
{
type zeroGradient ;
}

top
{
type fixedValue ;
value uniform (2 0 0) ;
}

plate
{
type noSlip ;
}

upstream
{
type symmetryPlane ;
}

frontAndBack
{
type empty ;
}
}

Temperature/Scalar Field
Create 0/T:
FoamFile
{
version 2.0;
format ascii ;
class volScalarField ;
object T;
}

dimensions [0 0 0 1 0 0 0];
internalField uniform 293; // 20 C ambient
CFD Lab Manual - OpenFOAM 2406 146

boundaryField
{
inlet
{
type fixedValue ;
value uniform 293;
}

outlet
{
type zeroGradient ;
}

top
{
type fixedValue ;
value uniform 293;
}

plate
{
type fixedValue ;
value uniform 373; // 100 C heated plate
}

upstream
{
type symmetryPlane ;
}

frontAndBack
{
type empty ;
}
}

Transport Properties
Create constant/transportProperties:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object tr a ns p o rt P ro p e rt i es ;
}

transportModel Newtonian ;
nu [0 2 -1 0 0 0 0] 1.5 e -05; // Air kinematic viscosity
DT [0 2 -1 0 0 0 0] 2.2 e -05; // Thermal diffusivity ( Pr =
0.7)

// For mass transfer

// D [0 2 -1 0 0 0 0] 1.5 e -05; // Mass diffusivity ( Sc =
1.0)
CFD Lab Manual - OpenFOAM 2406 147

Solve Velocity Field First

Run flow solver:
# Generate mesh
blockMesh

# Solve for velocity field

icoFoam

# Copy final velocity field to 0 directory

cp $ ( foamListTimes | tail -1) / U 0/

Solve Scalar Transport

Create system/[Link]:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object controlDict ;
}

application s c al a rT r a ns p or t F oa m ;
startFrom startTime ;
startTime 0;
stopAt endTime ;
endTime 10;
deltaT 0.001;
writeControl runTime ;
writeInterval 0.5;
writeFormat ascii ;
writePrecision 8;
writeCompression off ;
timeFormat general ;
timePrecision 6;
runTi meModifi able true ;

functions
{
heatFlux
{
type sur faceFiel dValue ;
libs (" l i b f i e l d F u n c t i o n O b j e c t s . so ") ;
writeControl timeStep ;
writeInterval 50;
fields (T);
operation sum ;
regionType patch ;
name plate ;

surfaceFormat none ;

// Calculate heat flux : q = -k * dT / dn

writeFields false ;
}
CFD Lab Manual - OpenFOAM 2406 148

nusseltNumber
{
type coded ;
libs (" l i b u t i l i t y F u n c t i o n O b j e c t s . so ") ;
name nusseltNumber ;
writeControl timeStep ;
writeInterval 50;

code
#{
const fvMesh & mesh = time () . lookupObject < fvMesh >(" region0 ")
;
const volScalarField & T = mesh . lookupObject < volScalarField
>(" T ") ;

// Get plate patch

label plateID = mesh . boundaryMesh () . findPatchID (" plate ") ;
const fvPatch & platePatch = mesh . boundary () [ plateID ];

// Calculate heat transfer coefficient

const scalarField & Twall = T . boundaryField () [ plateID ];
const scalarField & gradTn = T . boundaryField () [ plateID ].
snGrad () ;

scalar Tinf = 293.0; // Ambient temperature

scalar k = 0.026; // Thermal conductivity of air
scalar L = 0.1; // Plate length

// Local heat transfer coefficient

scalarField h = -k * gradTn / ( Twall - Tinf ) ;

// Local Nusselt number

scalarField Nu_local = h * mesh . boundary () [ plateID ]. Cf () .
component (0) / k ;

// Average values
scalar h_avg = gSum ( h * platePatch . magSf () ) / gSum (
platePatch . magSf () ) ;
scalar Nu_avg = h_avg * L / k ;

Info << " Average heat transfer coefficient : " << h_avg << "
W / $m ^2 $K " << endl ;
Info << " Average Nusselt number : " << Nu_avg << endl ;

// Theoretical correlation
scalar Re = 2.0 * L / 1.5 e -05; // Reynolds number
scalar Pr = 0.7; // Prandtl number
scalar Nu_theory = 0.664 * pow ( Re , 0.5) * pow ( Pr , 1.0/3.0) ;

Info << " Theoretical Nusselt number : " << Nu_theory << endl
;
Info << " Error : " << mag ( Nu_avg - Nu_theory ) / Nu_theory *
100 << "%" << endl ;
#};
}

te mp er atu re Pr ofi le
CFD Lab Manual - OpenFOAM 2406 149

{
type sets ;
libs (" libsampling . so ") ;
writeControl writeTime ;
writeInterval 1;
i nt e rp o l at i on S c he m e cellPoint ;
setFormat raw ;
sets
(
x_05cm
{
type uniform ;
axis y;
start (0.005 0 0.05) ;
end (0.005 0.02 0.05) ;
nPoints 40;
}
x_2cm
{
type uniform ;
axis y;
start (0.02 0 0.05) ;
end (0.02 0.02 0.05) ;
nPoints 40;
}
x_5cm
{
type uniform ;
axis y;
start (0.05 0 0.05) ;
end (0.05 0.02 0.05) ;
nPoints 40;
}
);
fields ( T U ) ;
}
}

Run scalar transport:

cp system / controlDict . scalarTransport system / controlDict
s ca l ar T r an s po r t Fo a m

11.3.2 Step 2: Pipe Flow Heat Transfer

Setup Pipe Case

cd ../ pipeFlow
mkdir 0 constant system

Copy pipe mesh from Experiment 10 and modify for heat transfer:
Create 0/T:
FoamFile
{
version 2.0;
format ascii ;
CFD Lab Manual - OpenFOAM 2406 150

class volScalarField ;
object T;
}

dimensions [0 0 0 1 0 0 0];
internalField uniform 293; // Cold inlet

boundaryField
{
inlet
{
type fixedValue ;
value uniform 293; // 20 C inlet
}

outlet
{
type zeroGradient ;
}

wall
{
type fixedValue ;
value uniform 373; // 100 C heated wall
}

centerline
{
type symmetryPlane ;
}

wedgeFront
{
type wedge ;
}

wedgeBack
{
type wedge ;
}
}

11.3.3 Step 3: Natural Convection

Vertical Plate Setup
Create case for natural convection:
cd ..
mkdir nat uralConv ection
cd natur alConvec tion
mkdir 0 constant system

Create system/blockMeshDict for vertical plate:

FoamFile
{
version 2.0;
CFD Lab Manual - OpenFOAM 2406 151

format ascii ;
class dictionary ;
object blockMeshDict ;
}

convertToMeters 0.01; // cm to m

vertices
(
(0 0 0) // 0 - plate bottom
(10 0 0) // 1 - domain width (10 cm )
(10 20 0) // 2 - domain top (20 cm height )
(0 20 0) // 3 - plate top
(0 0 0.1) // 4 - back face
(10 0 0.1) // 5
(10 20 0.1) // 6
(0 20 0.1) // 7
);

blocks
(
hex (0 1 2 3 4 5 6 7) (50 100 1) simpleGrading (5 1 1) // Grading
toward plate
);

edges () ;

boundary
(
plate
{
type wall ;
faces
(
(0 3 7 4)
);
}

right
{
type patch ;
faces
(
(1 2 6 5)
);
}

top
{
type patch ;
faces
(
(3 2 6 7)
);
}

bottom
{
CFD Lab Manual - OpenFOAM 2406 152

type patch ;
faces
(
(0 1 5 4)
);
}

frontAndBack
{
type empty ;
faces
(
(0 4 5 1)
(3 2 6 7)
);
}
);

Use buoyantSimpleFoam for natural convection with Boussinesq approximation.

11.3.4 Step 4: Numerical Schemes Comparison

Create system/[Link]:
FoamFile
{
version 2.0;
format ascii ;
class dictionary ;
object fvSchemes ;
}

ddtSchemes
{
default Euler ;
}

gradSchemes
{
default Gauss linear ;
}

divSchemes
{
default none ;
div ( phi , T ) Gauss upwind ; // First - order upwind
}

laplacianSchemes
{
default Gauss linear orthogonal ;
}

i n t e r p o l a t i on S c h e m e s
{
default linear ;
}
CFD Lab Manual - OpenFOAM 2406 153

snGradSchemes
{
default orthogonal ;
}

Create system/[Link]:
divSchemes
{
default none ;
div ( phi , T ) Gauss linearUpwind grad ( T ) ; // Second - order upwind
}

Create system/[Link]:
divSchemes
{
default none ;
div ( phi , T ) Gauss linear ; // Central differencing
}

11.3.5 Step 5: Peclet Number Study

Create parametric study script [Link]:
import numpy as np
import matplotlib . pyplot as plt
import subprocess

def r u n _ c o n v e c t i o n _ d i f f u s i o n ( Pe , scheme =" upwind ") :

""" Run convection - diffusion case for given Peclet number """

# Fixed velocity and length

U = 2.0 # m/s
L = 0.1 # m

# Calculate diffusion coefficient

D = U * L / Pe

# Update transport properties

with open ( ’ constant / transportProperties ’ , ’w ’) as f :
f . write (f ’ ’ ’ FoamFile
{{
version 2.0;
format ascii ;
class dictionary ;
object tr a ns p o rt P ro p e rt i es ;
}}

transportModel Newtonian ;
nu [0 2 -1 0 0 0 0] 1.5 e -05;
DT [0 2 -1 0 0 0 0] { D :.2 e };
’ ’ ’)

# Copy appropriate scheme

subprocess . run ([ ’ cp ’ , f ’ system / fvSchemes .{ scheme } ’ , ’ system /
fvSchemes ’])
CFD Lab Manual - OpenFOAM 2406 154

# Run simulation
subprocess . run ([ ’ scalarTransportFoam ’] , capture_output = True )

# Extract results ( placeholder )

return 0

def p l o t _ s c h e m e _ c o m p a r i s o n () :
""" Compare numerical schemes for different Peclet numbers """

Pe_values = np . logspace (0 , 3 , 20) # Pe = 1 to 1000

schemes = [ ’ upwind ’ , ’ linearUpwind ’ , ’ central ’]

plt . figure ( figsize =(15 , 5) )

# Numerical diffusion
plt . subplot (1 , 3 , 1)
for scheme in schemes :
# Theoretical numerical diffusion coefficients
if scheme == ’ upwind ’:
num_diff = 0.5 * Pe_values / ( Pe_values + 2)
elif scheme == ’ linearUpwind ’:
num_diff = 0.125 * Pe_values / ( Pe_values + 8)
else : # central
num_diff = np . zeros_like ( Pe_values )

plt . loglog ( Pe_values , num_diff , label = scheme )

plt . xlabel ( ’ Peclet Number ’)

plt . ylabel ( ’ Numerical Diffusion ’)
plt . title ( ’ Numerical Diffusion vs Peclet Number ’)
plt . legend ()
plt . grid ( True , alpha =0.3)

# Stability diagram
plt . subplot (1 , 3 , 2)
Pe_crit = np . array ([2 , 8 , 2]) # Critical Pe for each scheme
stability = [ ’ Stable ’ , ’ Conditionally Stable ’ , ’ Unstable ’]
colors = [ ’ green ’ , ’ orange ’ , ’red ’]

for i , scheme in enumerate ( schemes ) :

plt . axhline ( Pe_crit [ i ] , color = colors [ i ] , linestyle = ’ - - ’ ,
label =f ’{ scheme } ( Pe_crit = { Pe_crit [ i ]}) ’)

plt . xlabel ( ’ Peclet Number ’)

plt . ylabel ( ’ Critical Peclet Number ’)
plt . title ( ’ Stability Limits ’)
plt . legend ()
plt . grid ( True , alpha =0.3)
plt . yscale ( ’ log ’)

# Accuracy comparison
plt . subplot (1 , 3 , 3)
# Placeholder data for accuracy
for scheme in schemes :
if scheme == ’ upwind ’:
accuracy = 1 / Pe_values # First - order
elif scheme == ’ linearUpwind ’:
accuracy = 1 / Pe_values **2 # Second - order
CFD Lab Manual - OpenFOAM 2406 155

else : # central
accuracy = 1 / Pe_values **2 * np . where ( Pe_values < 2 , 1 , np
. inf )

plt . loglog ( Pe_values , np . abs ( accuracy ) , label = scheme )

plt . xlabel ( ’ Peclet Number ’)

plt . ylabel ( ’ Truncation Error ’)
plt . title ( ’ Accuracy Comparison ’)
plt . legend ()
plt . grid ( True , alpha =0.3)

plt . tight_layout ()
plt . savefig ( ’ n u m e r i c a l _ s c h e m e s _ c o m p a r i s o n . png ’)
plt . show ()

def v a l i d a t e _ h e a t _ t r a n s f e r () :
""" Validate heat transfer correlations """

# Experimental data for flat plate ( placeholder )

Re_exp = np . array ([1 e4 , 5 e4 , 1 e5 , 2 e5 , 5 e5 ])
Nu_exp = np . array ([32 , 71 , 100 , 141 , 224])

# Theoretical correlations
Re_theory = np . logspace (3 , 6 , 100)
Pr = 0.7

# Laminar
Nu_laminar = 0.664 * Re_theory **0.5 * Pr **(1/3)

# Turbulent
Nu_turbulent = 0.037 * Re_theory **0.8 * Pr **(1/3)

plt . figure ( figsize =(10 , 6) )

plt . loglog ( Re_theory , Nu_laminar , ’b - ’ , label = ’ Laminar (0.664 Re
^0.5 Pr ^1/3) ’)
plt . loglog ( Re_theory , Nu_turbulent , ’r - ’ , label = ’ Turbulent (0.037
Re ^0.8 Pr ^1/3) ’)
plt . loglog ( Re_exp , Nu_exp , ’ko ’ , markersize =8 , label = ’ Experimental
Data ’)

# Simulation results ( placeholder )

Nu_simulation = Nu_laminar * (1 + 0.1* np . random . randn ( len ( Re_theory
)))
plt . loglog ( Re_theory , Nu_simulation , ’g - - ’ , alpha =0.7 , label = ’
OpenFOAM Simulation ’)

plt . xlabel ( ’ Reynolds Number ’)

plt . ylabel ( ’ Nusselt Number ’)
plt . title ( ’ Heat Transfer Validation - Flat Plate ’)
plt . legend ()
plt . grid ( True , alpha =0.3)
plt . xlim (1 e3 , 1 e6 )
plt . ylim (10 , 1000)

plt . tight_layout ()
plt . savefig ( ’ h e a t _ t r a n s f e r _ v a l i d a t i o n . png ’)
plt . show ()
CFD Lab Manual - OpenFOAM 2406 156

# Run analysis
if __name__ == " __main__ ":
p l o t _ s c h e m e _ c o m p a r i s o n ()
v a l i d a t e _ h e a t _ t r a n s f e r ()

11.4 Expected Results

For forced convection over flat plate at Re = 105 , Pr = 0.7:

• Average Nusselt number: Nu ≃ 200 (theory: 196)

• Thermal boundary layer: δT /δ ≃ P r(−1/3) ≃ 1.2

• Heat transfer coefficient: h ≃ 50 W/m2 K

• Upwind scheme: stable but diffusive (Pe > 2)

• Central scheme: accurate but unstable (Pe > 2)

• LinearUpwind: good compromise (stable up to Pe ≃ 8)

• Natural convection: Ra ≃ 109 , Nu ≃ 100

11.5 Assignments
11.5.1 Assignment 1
Analyze the effect of Prandtl number (0.1, 0.7, 7.0) on thermal boundary layer develop-
ment and heat transfer coefficients.

11.5.2 Assignment 2
Investigate numerical schemes (upwind, central, bounded) for different Peclet numbers
and analyze numerical diffusion effects.

11.5.3 Assignment 3
Extend to natural convection over vertical plate and compare with Churchill-Chu corre-
lations for different Rayleigh numbers.