Contents I

1 Introduction
Computational Fluid Dynamics Computational Methods
Master of Science in Thermal/Automotive Engineering Partial Differential Equations
Governing Equations of Fluid Mechanics and Heat Transfer
2 Principles of Solution of the Governing Equations
Dr.-Ing. Getachew S. Tibba Discretization
Spatial Discretization
Addis Ababa Science and Technology University Temporal Discretization
Stability Analysis
Boundary and Initial Conditions
3 Modeling Diffusion
Steady State Diffusion
Transient Diffusion
4 Modeling Convection and Diffusion
May 25, 2021 Steady Convection and Diffusion
Central Differencing Scheme
G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 1 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 1 / 160

Contents II Contents III

Upwind Differencing Scheme Introduction
Hybrid Scheme Basic Equations of Turbulence
The Power-Law Scheme
First-Order Closures
The QUICK scheme
Second-Order Closures
Multidimensional Convection and Diffusion
Large-Eddy Simulation
False Diffusion
Direct Numerical Simulation
One-way Coordinate System
Transient Convection and Diffusion
5 Modeling Fluid flow
Difficulties in Fluid Flow Modeling 7 Fluid Flow Problems in I.C. Engines
Pressure-Velocity Coupling Flow Through Manifolds
The SIMPLE Algorithm Elements of Air Motion in Engines
The SIMPLER Algorithm Outline of Fluid Dynamics Models
The SIMPLEC Algorithm
The PISO Algorithm
6 Modeling Turbulence
G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 2 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 3 / 160
 Course Objectives Course Objectives

Course Objectives I Course Objectives II

This course brings together the knowledge gained in fluid

mechanics, thermodynamics, heat transfer and numerical methods
in order to develop computational techniques for the engineering
analysis of heat and fluid flow processes. Student will be able to
This course also intends to provide the students with sufficient 1 Under stand the governing of fluid flow, heat transfer and
background to understand the mathematical representation of the
numerical solution.
governing equations of fluid flow, discretization techniques, grid
generation, transformation equations and to numerically solve the 2 Numerically solve the fluid flow field using some popular CFD
flow field problems. techniques. Model fluid flow problems and heat transfer.
The student will be introduced to the modeling and computational
techniques that are incorporated in current computational fluid
dynamics (CFD) software.
The student will also have the opportunity to use a standard CFD
software package to analyse some complex flow situations.

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 4 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 5 / 160

Introduction Computational Methods Introduction Computational Methods

Simply put, computational fluid dynamics is numerically solving

governing equations of the form,
1 Experimental investigations
Based on direct measurement of physical quantities,
Give more reliable information about a physical process if proper
∂ρφ
∂t + ∇ · (~v ρφ) = ∇ · (Γ∇φ) + ~s (1) measurements have been taken.
Drawbacks:
Relatively expensive,
Time consuming,
Prediction of Processes in Transport Phenomena Difficulties in measuring,
Inaccuracy in reading instruments , and
Inaccuracies of measuring equipment.
Experimental Theoretical 2 Theoretical calculations
Rely on mathematical models,
Relatively cheap and less time consuming,
Have ability of providing complete information,
Analytical Computational
Simulate realistic and ideal conditions.
3 Theoretical calculations- limitations:
Figure 1: Where does Computational Methods Lie? Wrong mathematical models result in wrong output.

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 6 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 7 / 160
 Introduction Computational Methods Introduction Partial Differential Equations

In cases of problems for which adequate mathematical description

has not been obtained, it is impossible to rely on the results of 2 2 2
theoretical calculations.
∂ φ
a ∂∂xφ2 + b ∂x∂y + c ∂∂yφ2 + d ∂φ ∂φ
∂x + e ∂y + f φ = g(x, y)

Since analytical calculations are limited to simple problems, The governing equations of fluid flow and heat transfer are mainly
computational methods are becoming more and more popular methods in the form of partial differential equations.
of predicting physical processes, especially with rapid advancement of The main work of computational methods is, therefore, solving
computing technology. single or set of partial differential equations.
Any valid numerical solution of the governing equations should
Few application areas of computational methods:
exhibit the property of obeying the general mathematical
Automotive and engine applications, properties of partial differential equations.
Manufacturing applications, Eventhough partial differential equations have similar form, they show
Civil engineering applications, different properties that affect their solution procedures and
Air conditioning applications, application of boundary and initial conditions.
Fluid machine applications, · · ·

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 8 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 9 / 160

Introduction Partial Differential Equations Introduction Partial Differential Equations

Partial Differential Equation

The governing equations of marching problems are either parabolic
or hyperbolic partial differential equations.
Physical Classification Mathematical Classification Boundary Boundary
Condition Condition
Equilibrium Problems Marching Problems Parabolic Hyperbolic Elliptic

Condition
Boundary
R
Figure 2: Classification of Partial Differential Equations R

1 Physical Classification of Partial Differential Equations

Initial
Equilibrium problems are those which require closed domain on
Condition
which boundary conditions are prescribed. a) b)
A good example is steady state temperature distribution.
The governing equations of equilibrium problems are elliptic partial Figure 3: Physical Classification of Partial Differential Equations: a)
differential equations. Equilibrium Problem Domain, b) Marching Problem Domain
Marching problems are transient or transient-like problems.
They require an open domain for solution. The solution begins from 2 Mathematical Classification of Partial Differential Equations
some initial value and marches forward satisfying boundary Our discussion focuses on second order partial differential equations,
conditions.
G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 10 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 11 / 160
 Introduction Partial Differential Equations Introduction Partial Differential Equations

From chain rule,

Referring to equation (2), if the coefficients a, b, c, · · · are functions dp dx dy
=u +v (5)
of x, y, dφ dφ
dx , dy , φ the differential equation is called quasi-linear
dt dt dt
partial differential equation. and,
dq dx dy
∂2φ ∂2φ ∂2φ ∂φ ∂φ =v +w (6)
a 2 +b +c 2 +d +e + f φ = g(x, y) (2) dt dt dt
∂x ∂x∂y ∂y ∂x ∂y Equations (4), (5) and (6) can be combined as,
The mathematical classification of partial differential equations     
a b c u H
depends only on the second derivative terms of the equation.  dx dy 0   v  =  dp 
dt dt dt (7)
dy dq
∂ 2φ ∂2φ ∂2φ ∂φ ∂φ 0 dxdt dt
w dt
a 2 +b + c 2 = −d −e − f φ + g(x, y) = H (3)
∂x ∂x∂y ∂y ∂x ∂y
The condition for discontinuity will be,
2 2
∂φ ∂φ ∂ φ ∂ φ
Define, = p(t), = q(t), = u(t), = v(t) and  2
∂x ∂y ∂x2 ∂x∂y dy dy
2
∂ φ
= w(t) a −b +c=0 (8)
∂y 2 dx dx
With these definition, equation(3) can be written as
Solving for dy/dx gives,
au + bv + cw = H (4) √
dy b± b2 − 4ac
= (9)
dx 2a
G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 12 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 13 / 160

Introduction Partial Differential Equations Introduction Partial Differential Equations

The curves y(x) that satisfy equation (9) are called the
characteristics of the partial differential equation.
The characteristics are paths in the solution domain along which
information propagates. The terminology hyperbolic, parabolic, and elliptic chosen to
Discontinuities in derivatives of the dependent variable, if exist, also classify partial differential equations reflects the analogy between
propagate through these characteristic paths. the form of the discriminant, b2 − 4ac, for partial differential
Based on the discreminant D, the PDE can be classified as, equations reflects and the form of the discriminant, b2 − 4ac, which
classifies conic sections given by
D Type of PDE Characteristic Curve Families
ax2 + bxy + cy 2 + dx + ey + f = 0 (10)
b2 − 4ac > 0 Hyperbolic Two real
b2 − 4ac = 0 Parabolic Single real
b2 − 4ac < 0 Elliptic No real
If a, b and c are not constants, the classification may change from
point to point in the domain.

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 14 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 15 / 160
 Introduction Partial Differential Equations Introduction Partial Differential Equations

Two initial and two boundary conditions are needed.

Are open type.
Circle 2 Parabolic partial differential equations
Ellipse Example is the unsteady one dimensional heat conduction equation.
Parabola
Hyperbola ∂2T ∂T
a = (12)
∂x2 ∂t
One initial condition and two boundary conditions are needed to
solve equation(12).
Are open type.

Figure 4: Conic Sections Boundary Condition

Boundary Condition

1 Hyperbolic partial differential equations Initial Condition

Example is the wave equation.
2
Figure 5: Domain of the solution of parabolic partial differential equation
2∂ φ ∂2φ
a = (11)
∂x2 ∂t2 3 Elliptic equations

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 16 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 17 / 160

Introduction Partial Differential Equations Introduction Governing Equations of Fluid Mechanics and Heat Tran

We have three basic governing equations of fluid mechanics and heat

Example is steady two-dimensional heat equation. transfer.
∂2T ∂2T Continuity equation
+ =0 (13)
∂x 2 ∂y 2 Linear Momentum equation
Are closed type. Energy equation
Only boundary conditions must be specified on territories.
Boundary Condition
1 Continuity equation – Conservation of mass

∂ρ ∂ρu ∂ρv ∂ρw

+ + + =0 (14)
∂t ∂x ∂y ∂z
R

Continuity equation in compact form

Figure 6: Domain of the Solution of Elliptic Partial Differential Equation ∂ρ ~ ~)=0
+ ∇ · (ρV (15)
∂t
2 Linear Momentum equation – Newton’s Second Law of Motion

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 18 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 19 / 160
 Introduction Governing Equations of Fluid Mechanics and Heat Tran Introduction Governing Equations of Fluid Mechanics and Heat Tran

X-momentum equation
du ∂p ∂τxx ∂τyx ∂τzx
ρ =− + + + + ρfx (16) 3 Energy equation – First Law of Thermodynamics
dt ∂x ∂x ∂y ∂z
Y-momentum equation Energy equation in compact form
de ~ ~ ~ · (k ∇T
~ )+∇
~ · (V
~ · τ ) + Se
dv ∂p ∂τxy ∂τyy ∂τzy ρ + ∇ · (V p) = ∇ (20)
ρ =− + + + + ρfy (17) dt
dt ∂y ∂x ∂y ∂z
Z-momentum equation Generalized form of the governing equations:
dw ∂p ∂τxz ∂τyz ∂τzz ∂ ~ Φ) = div(ΓgradΦ) + S
ρ =− + + + + ρfz (18) (ρΦ) + div(ρV (21)
dt ∂z ∂x ∂y ∂z ∂t
Equations (16), (17) and (18) are the scalar forms of the linear The general procedures of computationally solving equation (21) are:
momentum equation in x-, y- and z- coordinates, respectively. 1 Dividing the geometry of the flow and/or heat transfer domain in
Linear momentum equation in compact form to smaller regions-meshing
~
dV
ρ ~ · σ + ρ~b
=∇ (19)
dt
G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 20 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 21 / 160

Introduction Governing Equations of Fluid Mechanics and Heat Tran Introduction Governing Equations of Fluid Mechanics and Heat Tran

(δx)w (δx)e

Δx

Figure 7: Discretization Example 5 Solving the system of algebraic equations,

2 Converting the governing partial differential equations in to 6 Interpreting results.
system of linear algebraic equations-descritization 7 Validating results. Computational solutions can be validated
  against analytical or experimental solutions.
d dT ke (TE − TP ) kw (TP − TW )
k +S =0⇒ − + SΔx = 0
dx dx (δx)e (δx)w
(22)
3 In more compact and standard form

aP TP = aE TE + aW TW + b (23)

4 Applying appropriate boundary and initial conditions,

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 22 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 23 / 160
 Principles of Solution Discretization Principles of Solution Discretization

Principles of Solution of the Governing Equations I Spatial Discretization I

Discretization
is the process of transferring continuous models and equations into
discrete counterparts. Spatial Discretization
results in a set of algebraic equations. Spatial discretization is the process of dividing the computational
Consider the general transport equation below: domain in to smaller sub domains. The most common spatial
discretization methods are:
∂ ~ Φ) = div(ΓgradΦ) + S
(ρΦ) + div(ρV (24) The finite difference method,
∂t
The finite volume method, and
The first term accounts for the variation of φ with time while the
rest terms account for variation of φ with space. The finite element method.
Disretization methods:
Spatial Discretization- with respect to space
Temporal Discretization- with respect to time

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 24 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 25 / 160

Principles of Solution Discretization Principles of Solution Discretization

Spatial Discretization II Spatial Discretization III

The Finite Difference Method

Directly applied to the differential form of the governing equations. The finite difference method is further classified as:
Taylor series expansion for the discretization of the derivatives of First Order Forward Difference
the transport variables.
First Order Backward Difference
Suppose we have a scalar function φ. Then according to Taylor
First Order Central Difference
expansion,
Second Order Central Difference
∂φ δx2 ∂2φ δx3 ∂3φ
φ(x0 + δx) = φ(x0 ) + δx |x + |x + |x + ... (25)
∂x o 2 ∂x2 o 3! ∂x3 o

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 26 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 27 / 160
 Principles of Solution Discretization Principles of Solution Discretization

Spatial Discretization IV Spatial Discretization V

First Order Forward Difference

Suppose δx = Δx in equation (25). Then solving for the first derivative, First Order Backward Difference
Suppose δx = −Δx in equation (25). Then solving for the first
∂φ φ(x0 + Δx) − φ(x0 ) Δx ∂ 2 φ Δx2 ∂ 3 φ derivative,
|xo = − |x − |x − ... (26)
∂x Δx 2 ∂x2 o
3! ∂x3 o
∂φ δx2 ∂ 2 φ δx3 ∂ 3 φ
Or, φ(x0 − δx) = φ(x0 ) − δx |xo + |x − |x + ... (29)
∂x 2 ∂x2 o
3! ∂x3 o
∂φ φ(x0 + Δx) − φ(x0 )
|xo = − O(Δx) (27)
∂x Δx Then rearranging and dropping the truncation error,
Removing the truncation error, ∂φ φ(x0 ) − φ(x0 − Δx)
|xo ≈ (30)
∂φ φ(x0 + Δx) − φ(x0 ) ∂x Δx
|x ≈ (28)
∂x o Δx

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 28 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 29 / 160

Principles of Solution Discretization Principles of Solution Discretization

Spatial Discretization VI Spatial Discretization VII

Second Order Central Difference
Writing the Taylor series expansions for φ(x0 + Δx) and φ(x0 − Δx),
and adding yields,
First Order Central Difference
∂2φ φ(x0 + Δx) − 2φ(x0 ) + φ(x0 − Δx)
If we add equations (28) and (30) we get, 2
|xo = + O(Δx2 ) (32)
∂x Δx2
∂φ φ(x0 + Δx) − φ(x0 − Δx) Or,
|xo ≈ (31)
∂x 2Δx ∂2φ φ(x0 + Δx) − 2φ(x0 ) + φ(x0 − Δx)
2
|xo ≈ (33)
There is a complication with equation (31) because it does not include ∂x Δx2
the value for φ(x0 ).
Example
Determine the temperature distribution in the one dimensional heat
conduction problem with uniform heat generation shown using the
direct and iterative finite difference methods and validate the results
with the analytical solution.

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 30 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 31 / 160
 Principles of Solution Discretization Principles of Solution Discretization

Spatial Discretization VIII Spatial Discretization IX

Given
W
q̇ = 50000 m Three interior nodes
q̇ 3
W
k k = 15 mK 1 o
Ti = 1000 C i 2 3
To = 2000 C

Figure 9: Meshed domain

L = 1.0m
Node 1: T2 −2T1 +Ti
Δx2 + kq̇ = 0 ⇒ 2T1 − T2 = kq̇ ∗ Δx2 + Ti = 308.33
x Node 2: T3 −2T2 +T1
Δx2
+ kq̇ = 0 ⇒ −T1 + 2T2 − T3 = kq̇ ∗ Δx2 = 208.33
Figure 8: Example 1 Node 3: To −2T3 +T2
Δx2
+ kq̇ = 0 ⇒ −T2 + 2T3 = kq̇ ∗ Δx2 + To = 408.33
The system of equations obtained can be solved in direct or iterative
methods.
The one
dimensional conduction heat transfer equation is,
d dT
dx k dx + q̇ = 0
2
Assuming uniform thermal conductivity, ddxT2 + kq̇ = 0

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 32 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 33 / 160

Principles of Solution Discretization Principles of Solution Discretization

Spatial Discretization X Spatial Discretization XI

Iterative Method Solution

In this method, initial value of each temperature will be assumed.
Direct Method Solution: Matrix inversion Then each temperature will be written in terms of other temperatures.
 Value of the temperature will be calculated based on values of the
   
2 −1 0  T1   308.33  
1 q̇ 2
−1 2 −1 T = 208.33 other temperatures. T1 = 2 k Δx + Ti + T2
 2  
0 −1 2 T3 408.33
  
T2 = 21 kq̇ Δx2 + T1 + T3
 
Solution T3 = 12 kq̇ Δx2 + T2 + To
 T1   437.50 C 
   

T2 = 566.670 C After about 50 iterations the following results will be obtained.

T3 487.50 C T1 = 437.4990 C
   
T2 = 566.6660 C
T3 = 487.4990 C
The iterative solution is computationally economical compared to the
direct method.
G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 34 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 35 / 160
 Principles of Solution Discretization Principles of Solution Discretization

Spatial Discretization XII Spatial Discretization XIII

Analytical Solution
2 Comment
q̇ 2
The governing equation: ddxT2 + kq̇ = 0 ⇒ T (x) = − 2k x + c1 x + c2
For this simple problem, there is no big difference between the finite
Inserting the known boundary temperatures, the constants c1 and c2
difference solution and the analytical solution.
are obtained to be: c1 = 1766.67 and c2 = 100.0
Then the analytical solution will be
Advantages of the finite difference method:
T1 = 437.50 C
T2 = 566.670 C Simplicity,
T3 = 487.50 C The possibility to easily obtain high-order approximations.
Solution for 3 interior nodes Solution for 30 interior nodes
Disadvantages of the finite difference method:
Requires a structured grid,
Cannot be directly applied in body-fitted (curvilinear)
coordinates- coordinate transformation needed,
Can be applied only to rather simple geometries.

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 36 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 37 / 160

Principles of Solution Discretization Principles of Solution Discretization

Spatial Discretization XIV Spatial Discretization XV

The Finite Element Method

The finite element method discretization transforms the governing
equations from the differential form into an equivalent integral form in Inserting equation (32) in to equation (31), we get a residual R
two different ways.
R = L(φ) 6= 0 (36)
Variational principle and
Method of weighted residuals or the weak formulation. The whole idea of the weighted residual method is to make R=0 over
The method of weighted residuals will be discussed here as it is more the entire domain by multiplying it with some weighting function W.
convenient for fluid mechanics and heat transfer problems. Consider Z
RWi dx = 0 (37)
div(ρ~v φ) − div(Γgradφ) − S = 0 = L(φ) (34)

Assumed solution of the form

φ(x) = a0 + a1 x + a2 x2 + ... + an xn (35)

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 38 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 39 / 160
 Principles of Solution Discretization Principles of Solution Discretization

Spatial Discretization XVI Spatial Discretization XVII

Types of the weighted residual method
Illustration
Depending on the choice of Wi there are different types of method of
weighted residuals: Consider the one dimensional equation,
Collocation method, d2 φ
+ f (x) = 0 (40)
Sub-domain method, dx2
Least squares method,
1 2 3 4
Body with 4 elements
Galerkin method, · · · 1 2 3 4 5
and 5 nodes

In Galerkin’s method, φ is assumed to vary in the element according to, Figure 10: Finite Element Mesh

φ(e) (x) = N1 φj + N2 φj+1 (38)

Insertion of equation (38) in to equation (40) gives us the residual
Where N1 and N2 are shape functions given by,
d2 φ(e)
R(e) = + f (x) (41)
xj+1 − x x − xj dx2
N1 = , N2 = (39)
xj+1 − xj xj+1 − xj
G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 40 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 41 / 160

Principles of Solution Discretization Principles of Solution Discretization

Spatial Discretization XVIII Spatial Discretization XIX

According to Galerkin method, the weighting function wi and the (e)
shape function Ni are the same. Then, equation (37) becomes Now let u = Ni and v = dφdx .
! Then integration by part and rearrangement of equation (42) gives us,
Z xj+1 Z xj+1
d dφ(e)
Ni dx + Ni f (x)dx = 0 (42) xj+1
dNi dφ(e) xj+1
dφ(e) xj+1
Z Z
xj dx dx xj dx = Ni f (x)dx + Ni |x (43)
xj dx dx xj dx j
To simplify this integral, integration by part is used. According to
integration by part rule, the derivative of product of two variables u For element 1 shown in Figure 25, inserting equation (38) in to
and v is given by, equation (43) yields,
x2 Z x2
dφ(1)
Z  
d(uv) = udv + vdu ⇒ udv = d(uv) − vdu dN1 dN1 dN2
φ1 + φ2 dx = N1 f (x)dx − |x (44a)
x1 dx dx dx x1 dx 1
Rb Rb
If performing integral a udv is difficult than performing a vdu, we can
x2 Z x2
dφ(1)
Z  
make replacement as, dN2 dN1 dN2
φ1 + φ2 dx = N2 f (x)dx + |x (44b)
x1 dx dx dx x1 dx 2
Z b Z b
udv = (uv)|ba − vdu
a a
G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 42 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 43 / 160
 Principles of Solution Discretization Principles of Solution Discretization

Spatial Discretization XX Spatial Discretization XXI

We can write equation (44) in compact form as,

" #  " (1) #
(1) (1) 
k11 k12 φ1 F1
(1) (1) = (1) (45)
k21 k22 φ2 F2 For element 3: " #
(3) (3)   " (3) #
k11 k12 φ3 F3
where R (3) (3) = (3) (47)
(1) x (1) R x2 dN1 dN2 k21 k22 φ4 F4
k11 = x12 dN 1 dN1
dx dx dx, k12 = x1 dx dx dx
(1) Rx (1) R x2 dN2 dN2
For element 4:
k21 = x12 dN 2 dN1
dx dx dx, k22 = x1 dx dx dx " #  " (4) #
(4) (4) 
(1) Rx (1) k11 k12 φ4 F4
F1 = x12 N1 f (x)dx − dφdx |x1 (4) (4) = (5) (48)
(1) Rx (1) k21 k22 φ5 F5
F2 = x12 N2 f (x)dx + dφdx |x2
Similar equations can be written for elements 2, 3 and 4 as follows.
For element 2: " #  " (2) #
(2) (2) 
k11 k12 φ2 F2
(2) (2) = (2) (46)
k21 k22 φ3 F3

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 44 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 45 / 160

Principles of Solution Discretization Principles of Solution Discretization

Spatial Discretization XXII Spatial Discretization XXIII

The next step is to assemble the equations written for the four
elements,
 (1) (1)
  (1)
 Example
k11 k12 0 0 0 φ1  F1
 (1) (1) (2) (2) Determine the temperature distribution in the one dimensional heat
+ 0 0  φ2  F2(1) + F1(2) 
  
k21 k22 k k
 (2)
11
(2)
12
(3) (3)
    (2) (3) 
 conduction problem with uniform heat generation shown using the
 φ3  = F2 + F1 
 0 k21 k22 + k11 k12 0  finite element methods and validate the results with the analytical
  

(3) (3) (4) (4)     (3) (4) 
 0 0 k21 k22 + k11 k12  φ4 F2 + F1  solution.

(4) (4) φ5 (4)
0 0 0 k21 k22 F2
(49)
Just as in the case of the fine difference case, the next step is to apply
boundary/initial conditions and solve equation (49).

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 46 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 47 / 160
 Principles of Solution Discretization Principles of Solution Discretization

Spatial Discretization XXIV Spatial Discretization XXV

1 2 3 4
1 2 3 4 5
Given
W
q̇ q̇ = 50000 m3
W
Ti = 1000 C k k = 15 mK
Figure 12: Example 2
To = 2000 C

We will follow the above procedure to solve this problem. We have four
L = 1.0m elements as shown in Figure 5. Since the temperatures at nodes i and o
x are known, the unknown temperatures are only at nodes 1, 2 and 3.
So, our coefficient matrix in equation (49) will be a 3x3 one.
Figure 11: Example 2 The problem is governed by the one dimensional heat equation with
heat source.
d dT


dx k dx + q̇ = 0
Assuming uniform thermal conductivity the above equation becomes,
d2 T q̇
dx2 + k = 0
G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 48 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 49 / 160

Principles of Solution Discretization Principles of Solution Discretization

Spatial Discretization XXVI Spatial Discretization XXVII

So we have, φ = T and f (x) = kq̇
xj+1 −x x−xj
From equation (39) we have, N1 = xj+1 −xj , N2 = xj+1 −xj
Then dN −1 −1 dN2 1 1 Element 2:  
dx = x2 −x1 = Δx and dx = x2 −x1 = Δx for each
element.
1
  " q̇Δx dT (2) #
Now we can determine the coefficient matrix and force vector for each 1 1 −1 T2 2k − dx |x2
Δx −1 = q̇Δx dT (2)
1 T3
element. 2k + dx |x3
Element Element 3:  
R x1:  " q̇Δx dT (3)
 #
(1) (1) R x2 −1 1
1 −1 T3 2k − dx |x3
−1 −1 1 −1
k11 = x12 Δx Δx dx = Δx , k12 = x1 Δx Δx dx = Δx
1
= q̇Δx
Δx −1 1 T4 dT (3)
2k + dx |x4
(1) R x 1 −1 −1 (1) R x2 1 1 1
k21 = x12 Δx Δx dx = Δx , k22 = x1 Δx Δx dx = Δx
(1) x dT (1) −q̇Δx dT (1)
Element 4:  
2 −x q̇
F1 = x12 xΔx
R
k dx − dx |x1 = 2k − dx |x1
  " q̇Δx dT (4) #
1 1 −1 T4 2k − dx |x4
(1) Rx 1 q̇ dT (1) −q̇Δx dT (1) = q̇Δx
F2 = x12 x−x Δx k dx + dx |x2 = 2k + dx |x2
Δx −1 1 T5 dT (4)
2k + dx |x5
Or in more compact way" we can write # Now we assemble the above equations of the four elements
q̇Δx dT (1)
  
1 1 −1 T1 2k − dx |x1
Δx −1 = q̇Δx dT (1)
1 T2
2k + dx |x2
Similar equations can be written for all the four elements.
G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 50 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 51 / 160
 Principles of Solution Discretization Principles of Solution Discretization

Spatial Discretization XXVIII Spatial Discretization XXIX

  
1 −1 0 0 0 T1 = Ti
−1 2 −1 0 2
0  T2  −T3 + 2T4 − T5 = q̇Δx
1 
   k
Δx  0 −1 2 −1 0     T3 = Or
 2
0 0 −1 2 −1  T4 2T2 − T3 = q̇Δx
k + T1

0 0 0 −1 1 T5 = To 2
   −T2 + 2T3 − T4 = q̇Δxk
q̇Δx dT (1) q̇Δx dT (1)

− | q̇Δx2
 q̇Δx dT (1) 2k dx x1
(2)
2k − dx |x1 −T3 + 2T4 = k + T5
 2k + dx |x2 + q̇Δx dT
2k − dx |x2 
q̇Δx

  or

 k 
  q̇Δx2
(2) (3) q̇Δx
 q̇Δx 
dT q̇Δx dT =
+ | + − | k + T1
   

 2k dx x3 2k dx x3  k  2 −1 0 T2
 q̇Δx dT (3) (4) q̇Δx q̇Δx2
|x4 + q̇Δx − dT |x4  
  
−1 2 −1  T3  = 
+

k

 2k dx 2k dx (4)
 k 
q̇Δx (4) q̇Δx dT 0 −1 2 T4 q̇Δx2
+ dTdx |x5
2k 2k + dx |x5 k + T5
Since T1 and T5 are known, the above five equations can be reduced to Now every term in the force vector are known. Therefore, the equation
three equations by removing the first and last rows of the coefficient can be solved iteratively or by the direct method.
matrix. The solution of these equations is again, like the finite difference
2
−T1 + 2T2 − T3 = q̇Δxk solution, T2 = 437.5◦ C, T3 = 566.6666666◦ C and T4 = 487.5◦ C
2
−T2 + 2T3 − T4 = q̇Δxk
G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 52 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 53 / 160

Principles of Solution Discretization Principles of Solution Discretization

Spatial Discretization XXX Spatial Discretization XXXI

The Finite Volume Method Referring to Figure 13, integration of equation (50) between points
directly utilizes the conservation laws - the integral formulation w and e gives us,
over a differential volume. Z e   Z e
d kdT
divides physical space into a number of arbitrary polyhedral dx + Sdx = 0 (51)
w dx dx w
control volumes
surface integral is approximated by the sum of the fluxes crossing
the individual faces (δx)w (δx)e

Consider the one dimensional heat conduction equation with

W w P e E
uniform heat generation δx
 
d kdT
+S =0 (50) Figure 13: Grid used for the Discretization of the Heat Equation
dx dx

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 54 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 55 / 160
 Principles of Solution Discretization Principles of Solution Discretization

Spatial Discretization XXXII Spatial Discretization XXXIII

Equation (51) can be further modified as

Z e   Z e
kdT
d + Sdx = 0 (52)
w dx w
Equation (56) can be written in more standard form of Finite
    Z e
kdT kdT Volume Method representation
⇒ − + Sdx = 0 (53)
dx e dx w w
⇒ aP TP = aE TE + aW TW + b (57)
k(TE − TP ) k(TP − TW
⇒ − + Sδx = 0 (54)
(δx)e (δx)w
 
k k k k
⇒ TE − + TP + TW + Sδx = 0 (55)
(δx)e (δx)e (δx)w (δx)w
⇒ aE TE − aP TP + aW TW + Sδx = 0 (56)

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 56 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 57 / 160

Diffusion Diffusion Steady State Diffusion

Modeling Diffusion I Steady State Diffusion I

Consider the general transport equation
∂ ~ Φ) = div(ΓgradΦ) + S
(ρΦ) + div(ρV (58)
∂t
Chapter Objectives For steady state diffusion, equation (58) reduces to
After completing the chapter, students should be able to: div(ΓgradΦ) + S = 0 (59)
Discretize the equations governing conduction heat transfer,
and for steady state heat transfer, equation (59) becomes
Monitor the stability of the discretized equations, Apply
appropriate boundary and/or initial conditions to the discretized div(k ∗ grad(T )) + q̇ = 0 (60)
equations, Solve the discretized equations.
where k is thermal conductivity, q̇ is heat generated per unit volume.
Equation (60) can be given in expanded form in rectangular coordinate
system as
     
∂ ∂T ∂ ∂T ∂ ∂T
k + k + k + q̇ = 0 (61)
∂x ∂x ∂y ∂y ∂z ∂z
G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 58 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 59 / 160
 Diffusion Steady State Diffusion Diffusion Steady State Diffusion

Steady State Diffusion II Steady State Diffusion III

N For two-dimensional case, we can refer to the grid in Fig. 14 and

the governing equation can be discretized in the following way.

Δyn
n
In addition to the west (W) and east (E) nodes for
one-dimensional case, the node of interest P has neighboring nodes

Δy
W w P e E in south (S) and in north (N) directions for two-dimensional cases.

Δys
s The integration over the control volume will be,
Z eZ n   Z eZ n   Z eZ n
S ∂ ∂T ∂ ∂T
k dxdy+ k dxdy+ q̇dxdy = 0
w s ∂x ∂x w s ∂y ∂y w s
Δxw Δxe (62)
Δx Similar to the case of the one-dimensional heat conduction, this
integral finally results in,
Figure 14: Grid used for the Discretization of Two-Dimensional Heat ∂T ∂T ∂T ∂T
Conduction Ae ke |e −Aw kw |w +An kn |n −As ks |s +q̇ΔxΔy = 0 (63)
∂x ∂x ∂y ∂y
G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 60 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 61 / 160

Diffusion Steady State Diffusion Diffusion Steady State Diffusion

Steady State Diffusion IV Steady State Diffusion V

⇒ aP TP = aE TE + aW TW + aN TN + aS TS + b (64)
where aP = AΔx
e ke
e
+ AΔxw kw
w
+ AΔx
n kn
n
+ AΔx
s ks
s
, aE = AΔx
e ke
e
, aW = Aw k w
Δxw ,
Δxe
An k n As k s
aN = Δxn , aS = Δxs Δxe− Δxe+
For three-dimensional problems, equation (64) becomes,
P e E
aP TP = aE TE + aW TW + aN TN + aS TS + aB TB + aT TT + b (65)
Figure 15: Approximation of Thermal Conductivity
For materials with variable thermal conductivity, k, thermal
conductivity at faces, for example at face e, is given by

ke = fe kP + (1 − fe )kE (66)
Δxe+
where fe = Δxe

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 62 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 63 / 160
 Diffusion Steady State Diffusion Diffusion Steady State Diffusion

Steady State Diffusion VI Steady State Diffusion VII

T Solution of the Linear Algebraic Equations

The system of algebraic equations resulting from discretisation of the
governing equation can be solved by direct or iterative methods.
t E
N
Direct Method
e Matrix inversion and the Tri-diagonal Matrix algorithm (TDMA) are
n
P s the commonly known direct methods.
w
S Matrix inversion
W b
[A]{T } = {B} (67)
⇒ {T } = [A]−1 {B} (68)
B
The matrix inversion method is limited in that the storage of the
Figure 16: Grid used for the Discretization of Three-Dimensional Heat
Conduction coefficient matrix requires large computer memory.

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 64 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 65 / 160

Diffusion Steady State Diffusion Diffusion Steady State Diffusion

Tri-diagonal Matrix algorithm, TDMA I Tri-diagonal Matrix algorithm, TDMA II

The Tri-diagonal Matrix algorithm (TDMA) is based on forward and ⇒ α2 T2 = β2 T3 + γ2 (P1 T2 + Q1 ) + b (74)
backward substitutions.
⇒ (α2 − γ2 P1 )T2 = β2 T3 + γ2 Q1 + b (75)
The one-dimensional conduction equation, aP TP = aE TE + aW TW + b,
can be rewritten as β2 γ2 Q 1 + b
⇒ T2 = T3 + (76)
α2 − γ2 P1 α2 − γ2 P1
αi Ti = βi Ti+1 + γi Ti−1 + b (69) ⇒ T2 = P2 T3 + Q2 (77)
For known temperatures at both boundaries, equation (69) becomes This can be generalized as,

α1 T1 = β1 T2 + γ1 T0 + b (70) ⇒ Ti = Pi Ti+1 + Qi (78)

β1 γ1 T0 + b where
⇒ T1 = T2 + (71) βi
α1 α1 Pi =
αi − γi Pi−1
⇒ T1 = P1 T2 + Q1 (72)
and
For i=2, γi Qi−1 + b
Qi =
α2 T2 = β2 T3 + γ2 T1 + b (73) αi − γi Pi−1
G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 66 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 67 / 160
 Diffusion Steady State Diffusion Diffusion Steady State Diffusion

Tri-diagonal Matrix algorithm, TDMA III TDMA for Multidimension I

In forward substitution, Pi and Qi are determined.
For i=n-1,
Tn−2 = Pn−2 Tn−1 + Qn−2 (79)
Since Tn−1 is known from boundary condition, the other T values are
determined by back substitution.
y
Chosen line Chosen line
Tn−3 = Pn−3 Tn−2 + Qn−3
x

Tn−4 = Pn−4 Tn−3 + Qn−4 Figure 17: Grid for the Line-by-Line Method
...
In Fig. 17, ’*’ show nodes with assumed temperatures while ’o’ show
...
nodes whose temperature to be calculated.
... For multi-dimensional problems, the TDMA method combined
T2 = P2 T3 + Q2 with an iterative method is employed through what is called the
’Line-by-line’ method, shown in Fig. 17.
G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 68 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 69 / 160

Diffusion Steady State Diffusion Diffusion Steady State Diffusion

TDMA for Multidimension II TDMA for Multidimension III

The simplest iterative method is the Gauss-Seidel method,

P
anb Tnb + b
TP = (80)
aP
This implies that iterative method is needed to come to exact
We can write the discretized equation for node P in the usual way solution.
as, Equation (83) is solved for a given line using TDMA and the
aP TP = aN TN + aS TS + aE TE + aW TW + b (81) updated result will be used for the next line.
As TW and TE are assumed to be known, the above equation can This procedure continues repeatedly for each line until a converged
be modified as, solution is obtained.
aP TP = aN TN + aS TS + b∗ (82)
⇒ αi Ti = βi Ti+1 + γi Ti−1 + b ∗ (83)
Since b∗ contains assumed values of TE and TW , the solution of
equation (83) will nor be exact.

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 70 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 71 / 160
 Diffusion Transient Diffusion Diffusion Transient Diffusion

Transient Diffusion I Transient Diffusion II

Spatial integral of the first left hand side term is performed like in
the steady case,
For transient condition, the heat diffusion equation in rectangular Z e  
coordinates is given by, ∂ ∂T
k dx = (aE TE − aP TP + aW TW )
      w ∂x ∂x
∂ ∂T ∂ ∂T ∂ ∂T ∂T
k + k + k + q̇ = ρc (84) and as temperature is assumed constant between sides w and e,
∂x ∂x ∂y ∂y ∂z ∂z ∂t
right hand side integral becomes
For simplicity, let us start with one dimensional transient Z e
conduction. FVM integration of equation (84) gives us, ∂T ∂Tp
ρc dx = ρc Δx
w ∂t ∂t
Z Δt Z e   Z Δt Z e
∂ ∂T ∂T
k dxdt + q̇ΔxΔt = ρc dxdt (85) Then equation (85) becomes,
0 w ∂x ∂x 0 w ∂t
Z Δt Z Δt
∂Tp
(aE TE − aP TP + aW TW )dt + q̇ΔxΔt = ρc Δxdt (86)
0 0 ∂t

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 72 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 73 / 160

Diffusion Transient Diffusion Diffusion Transient Diffusion

Transient Diffusion III Transient Diffusion IV

Solving for TPΔt gives us,

The first left hand side integral is performed if we know the 
Δx

variations of TW , TP and TE with time. One common method is θaP + ρc T Δt = θ(aE TEΔt + aW TW
Δt
) + (1 − θ)(aE TE0 + aW TW
0
)+
Δt P
assuming  
Δx
T = θT Δt + (1 − θ)T 0 ρc − aP (1 − θ) TP0 + q̇Δx
Δt
where θ is a weighting parameter (0 ≤ θ ≤ 1). (88)
This makes the integral in equation (86) to be, For isotropic material with equal grid size, the above equation
becomes
θ(aE TEΔt − aP TPΔt + aW TW
Δt
) + (1 − θ)(aE TE0 − aP TP0 + aW TW
0
)+
   
Δx Δt 2k Δx Δt k Δt k Δt
q̇Δx = ρc (TP − TP0 ) θ + ρc T =θ T + T +
Δt Δx Δt P Δx E Δx W
(87) 
k 0 k 0
 
Δx 2k

(1 − θ) T + T + ρc − (1 − θ) TP0 + q̇Δx
Δx E Δx W Δt Δx
(89)
G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 74 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 75 / 160
 Diffusion Transient Diffusion Diffusion Transient Diffusion

Transient Diffusion V Explicit Scheme I

For explicit (θ = 0) case, equation (89) is modified to,

 
Δx Δt k 0 k 0 Δx 2k
ρc T = T + T + ρc − TP0 + q̇Δx (91)
Equation (89) can be written in compact form as, Δt P Δx E Δx W Δt Δx

In explicit scheme, the current time temperature TPt depends on

X X
aP TPΔt = Δt
anb [θTnb 0
+(1−θ)Tnb ]+[a0P −(1−θ) anb ]TP0 +b (90)
only previous time neighboring temperatures.
Depending on the value of θ we select, the solution scheme can be To obtain stable solution, all the coefficients of the temperatures
explicit or implicit. should be non-negative. This implies that ρc Δx 2k
Δt − Δx > 0. That is
the Explicit scheme is conditionally stable.
For one dimensional heat conduction the stability condition
requires,
(Δx)2
Δt < ρc (92)
2k

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 76 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 77 / 160

Diffusion Transient Diffusion Diffusion Transient Diffusion

Explicit Scheme II CRANK-NICOLSON Scheme I

For the implicit CRANK-NICOLSON scheme (θ = 0.5) case,

equation (89) is modified to,
   
k Δx Δt k Δt k Δt
+ ρc T = 0.5 T + T +
Similarly, for three dimensional conduction with Δx = Δy = Δz Δx Δt P Δx E Δx W
and isotropic material, the stability condition requires that     (94)
k 0 k 0 Δx k
0.5 TE + TW + ρc − TP0 + q̇Δx
Δx Δx Δt Δx
(Δx)3 6k (Δx)4
ρc − > 0 ⇒ Δt < ρc (93)
Δt Δx 6k Unlike the case of explicit scheme, the current time temperature
TPt in the CRANK-NICOLSON case depends on both the previous
and the current time neighboring temperatures.
The CRANK-NICOLSON is also conditionally stable and its
stability condition requires

(Δx)2
Δt < ρc . (95)
k
G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 78 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 79 / 160
 Diffusion Transient Diffusion Diffusion Transient Diffusion

CRANK-NICOLSON Scheme II Fully-Implicit Scheme I

The discretization equation for fully implicit scheme (θ = 1.0) can

be written as,
   
2k Δx Δt k Δt k Δt
+ ρc T = T + T +
Δx Δt P Δx E Δx W
The current time temperatures can be found through iteration by (96)
Δx 0
the TDMA for one-dimensional problems and by the line-by-line ρc T + q̇Δx
method for two-and three-dimensional cases. Δt P
In the fully implicit scheme the current time temperature TPΔt
depends only on the current time neighboring temperatures.
The fully implicit scheme is unconditionally stable as all
coefficients of the temperatures are positive.
Just like in the case of the Crank-Nicolson scheme, the current
time temperatures can be found by iteration.

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 80 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 81 / 160

Diffusion Transient Diffusion Diffusion Transient Diffusion

Example I Example II
Consider a one-dimensional heat transfer through a 0.9m thick wall
shown below. Taking grid spacing of Δx = 0.15, determine the Solution:
transient temperature distribution in the wall at 10 time points. After As top and bottom sides of the geometry are insulated, the problem
what time will the steady state temperature be reached? Use the three can be considered as one-dimensional. The governing equation of this
transient schemes. problem is
Insulated ∂2T q̇ ρcp ∂T
2
+ =
∂x k k ∂t
Modified version of equation (89) can be used to solve this problem.
W
k=100 mK
ρcp Δx2
 
h=20 mW
TpΔt = θ TEΔt + TW
Δt
+ (1 − θ) TE0 + TW0
2K



W
q=1000 m2 W
q̇ = 3000 m
2θ + +
3 k Δt
ρc = 10000 mJ3 K T∞ = 300K
ρcp Δx2
 
−2 + 2θ + Tp0
k Δt

Insulated
G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 82 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 83 / 160
 Convection and Diffusion Steady Convection and Diffusion Convection and Diffusion Steady Convection and Diffusion

Modeling Convection and Diffusion I Modeling Convection and Diffusion II

For one-dimensional problems with no source term, the above
In problems where there is a flow of fluid, the effect of convection equation reduces to,
can be significant and combined to that of diffusion(conduction).
div(ρuφ) = div(Γgrad(φ)) (99)
The general transport equation in differential form:

∂(ρφ) (δx)w (δx)e

+ div(ρ~v φ) = div(Γgrad(φ)) + Sφ (97)
∂t
W w P e E
where δx
φ is a general property, Γ is the diffusion coefficient, ~v is the
velocity, and Sφ is the source term. Figure 18: One Dimensional Grid
This general governing equation for steady state diffusion and
Integration of equation (99) over the control volume shown in
convection becomes,
Figure 18 gives us,
div(ρ~v φ) = div(Γgrad(φ)) + Sφ (98) 
dφ
 
dφ

(ρuφ)e − (ρuφ)w = Γ − Γ (100)
dx e dx w
G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 84 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 85 / 160

Convection and Diffusion Steady Convection and Diffusion Convection and Diffusion Steady Convection and Diffusion

Modeling Convection and Diffusion III Central Differencing Scheme I

The diffusive terms on the right hand side of equation (100) can The central differencing scheme assumes linear-piecewise profile for
be represented in a similar manner as we saw in chapter 3. φ.
φe = 0.5(φP + φE ), φw = 0.5(φP + φW )
There are different schemes used to represent the convective term
(ρuφ) on faces e and w. Then equation (4.4) can be written as,
These schemes include:
(φE − φP ) (φP − φW )
Central Differencing Scheme 0.5(ρu)e (φP +φE )−0.5(ρu)w (φP +φW ) = Γe −Γw
Upwind Differencing Scheme δxe δxw
(101)
Hybrid Differencing Scheme
Power-Law Scheme 0.5Fe (φP + φE ) − 0.5Fw (φP + φW ) = De (φE − φP ) − Dw (φP − φW )
QUICK Scheme (102)
Where, F=flow(convection) strength and D=diffusion conductance

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 86 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 87 / 160
 Convection and Diffusion Steady Convection and Diffusion Convection and Diffusion Steady Convection and Diffusion

Central Differencing Scheme II Central Differencing Scheme III

The discretization equation can be written in standard form as, Boundedness- when a convergent solution of iterative methods is
obtained, the scheme is said to be bounded.
aP φP = aW φW + aE φE (103)
For the solution to be convergent (bounded scheme) the
where Scarborough condition must be satisfied,
P
aE = De − 0.5Fe , aW = Dw + 0.5Fw , aP = De + 0.5Fe + Dw − 0.5Fw |anb | n≤1 at all nodes
<1 at one node, at least (104)
|aP |
To get realistic solution of convection-diffusion problems, the �
fundamental properties of discretization schemes must be fulfilled, where |a|aPnb |
| is sum of neighboring coefficients, and aP is net
Conservativeness, coefficient of central node P.
Boundedness, and
For converged solution, coefficients should be diagonally dominant.
Transportiveness.
Transportiveness- In cases where convection is more dominant over
Conservativeness- the flux of φ leaving a given cell across a face
diffusion, conditions down stream do not have significant effect on
should be equal to the flux of φ entering adjacent cell across the
upstream conditions.
same cell.

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 88 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 89 / 160

Convection and Diffusion Steady Convection and Diffusion Convection and Diffusion Steady Convection and Diffusion

Central Differencing Scheme IV Central Differencing Scheme V

A dimensionless number used to measure the transportiveness of a
flow problem is the Peclet number given by,
F ρu
Pe = = Γ (105)
D δx

where Pe=0 ⇒ purely diffusive flow while Pe=∞ ⇒ purely The central differencing scheme assumes equal effect of upstream
convective flow. and downstream conditions at a node, so it lacks the
transportiveness property at high Pe.
The central differencing scheme is always conservative.
Referring to the coefficients aW = Dw + 0.5Fw and
aE = De − 0.5Fe , there is a possibility for these coefficients to be
negative, which gives unstable solution.
F
aE = De − 0.5Fe > 0 ⇒ = Pe < 2
D
Central differencing scheme is bounded for Pe<2.
G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 90 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 91 / 160
 Convection and Diffusion Steady Convection and Diffusion Convection and Diffusion Steady Convection and Diffusion

Upwind Differencing Scheme I Upwind Differencing Scheme II

φw = φW , φe = φP
u φW φw φE
φP φe Under the standard discretized equation, aP φP = aW φW + aE φE ,
the coefficients are given by,

W w P e E aE = De + max(0, −Fe )

φw = φP , φe = φE aW = Dw + max(0, Fw )
φW φw φP φe φE
u aP = aE + aW + (Fe − Fw )
The upwind scheme is conservative, bounded and transportive.
w e The accuracy of Upwind scheme is first order.
W P E

Figure 19: Upwind Differencing Scheme

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 92 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 93 / 160

Convection and Diffusion Steady Convection and Diffusion Convection and Diffusion Steady Convection and Diffusion

Hybrid Scheme I The Power-Law Scheme I

The hybrid scheme combines the central differencing and the The power-law scheme gives more accurate results than the hybrid
upwind differencing schemes. scheme.
For Pe<2, the central differencing is employed, while
In power-law scheme,
For Pe>2, the upwind scheme is employed for convection and the
diffusion is made zero. The diffusion is set to zero for Pe>10.
For 0<Pe<10 the flux at each face is evaluated using a polynomial
Under the standard discretized equation, aP φP = aW φW + aE φE , expression.
the coefficients are given by,
Under the standard discretized equation, aP φP = aW φW + aE φE ,
aE = max[0, (De − 0.5Fe ), −Fe ] the coefficients are given by,

aW = max[0, (Dw + 0.5Fw ), Fw ] aE = De max[0, (1 − 0.1|P ee |)5 ] + max[0, −Fe ]

aP = aE + aW + (Fe − Fw )
aW = Dw max[0, (1 − 0.1|P ew |)5 ] + max[0, Fw ]
The hybrid scheme gives physically realistic and stable solutions.
aP = aE + aW + (Fe − Fw )
The hybrid scheme is conservative, unconditionally bounded and
transportive. The power-law scheme is conservative, bounded and transportive.
The accuracy of the hybrid scheme is only first order. It is more accurate for one-dimensional problems.
G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 94 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 95 / 160
 Convection and Diffusion Steady Convection and Diffusion Convection and Diffusion Steady Convection and Diffusion

The QUICK scheme I The QUICK scheme II

The QUICK (Quadratic Upstream Interpolation for Convection where
Kinetics) uses a three-point upstream-weighted interpolation for 6 1 3
aW = Dw + αw Fw + αe Fe + (1 − αw )Fw
face cell values. 8 8 8
Referring to Fig. 20, the face values of φ are obtained as, 1
aW W = − αw Fw
8
aP φP = aW φW + aE φE + aW W φW W + aEE φEE (106)
3 6 1
Piece wise interpolation functions aE = De − αe Fe − (1 − αe )Fe + (1 − αw )Fw
8 8 8
φW
φW W φw φE φEE 1
u φP
φe aEE = − (1 − αe )Fe
8

aP = aE + aW + aEE + aW W + (Fe − Fw )
WW W w P e E EE    
|Fw | |Fe |
Figure 20: QUICK Scheme αw = max 0, , αe = max 0,
Fw Fe
G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 96 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 97 / 160

Convection and Diffusion Steady Convection and Diffusion Convection and Diffusion Steady Convection and Diffusion

The QUICK scheme III The QUICK scheme IV

Example
Water at a temperature of 400K enters a tube and leaves at a
temperature of 300K. Determine the temperature distribution in the
The QUICK scheme is:
tube assuming one-dimensional problem using the Central differencing,
Conservative, Upwind differencing, Hybrid differencing, Power-law, and QUICK
Third order accurate, kg
schemes. For water assume density to be ρ=1000 m3 and thermal
Transportive, and
conductivity to be k=0.6W/mK. Take the velocity to be u=0.001m/s.
Conditionally stable.
Compare the results with analytical solution.
To alleviate the stability problems, continuous modifications are
being made to the standard QUICK scheme. Solution
A program has been written using Python programming language and
the graphical result is given below. By varying the velocity and the
number of volumes see the effect of the Peclet number.

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 98 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 99 / 160
 Convection and Diffusion Steady Convection and Diffusion Convection and Diffusion Steady Convection and Diffusion

The QUICK scheme V The QUICK scheme VI

Peclet number=�0.42 Peclet number=�0.07

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 100 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 101 / 160

Convection and Diffusion Steady Convection and Diffusion Convection and Diffusion Steady Convection and Diffusion

The QUICK scheme VII The QUICK scheme VIII

Peclet number=�0.02

Peclet number=�4.12

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 102 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 103 / 160
 Convection and Diffusion Steady Convection and Diffusion Convection and Diffusion Multidimensional Convection and Diffusion

The QUICK scheme IX Multidimensional Convection and Diffusion I

For two-dimensional steady problem with no source term, the

transport equation can be modified as,
peclet number= 0.70
   
∂(ρuφ) ∂(ρvφ) ∂ ∂φ ∂ ∂φ
+ = Γ + Γ (107)
∂x ∂y ∂x ∂x ∂y ∂y

Referring to Figure 21 we can define the total fluxes at faces as,

∂φ ∂φ
Jx = ρuφ − Γ and Jy = ρvφ − Γ (108)
∂x ∂y

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 104 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 105 / 160

Convection and Diffusion Multidimensional Convection and Diffusion Convection and Diffusion Multidimensional Convection and Diffusion

Multidimensional Convection and Diffusion II Multidimensional Convection and Diffusion III

The continuity equation can be discretized as,

N
Δx
(ρu)e Δy − (ρu)w Δy + (ρv)n Δx − (ρv)s Δx = 0 (110)
n
Jn

⇒ Fe − Fw + Fn − Fs = 0 (111)
Δy

Jw Je
W w e E
In standard form, the discretized equation becomes,
s
aP φP = aE φE + aW φW + aN φN + aS φS (112)
Js

S
where

Figure 21: Control Volume for Two-Dimensional Problem aP = aE + aW + aN + aS + Fe − Fw + Fn − Fs

The discretized form of equation (107) will be, Expressions for aE , aW , aN and aS vary for different schemes.
They can be obtained from one-dimensional cases with slight
Je − Jw + Jn − Js = 0 (109) modifications as follows.
G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 106 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 107 / 160
 Convection and Diffusion Multidimensional Convection and Diffusion Convection and Diffusion Multidimensional Convection and Diffusion

Multidimensional Convection and Diffusion IV Multidimensional Convection and Diffusion V

Central Differencing Scheme: aE = De − F2e , aW = DW + F2w ,

aN = Dn − F2n , aS = Ds + F2s QUICK Scheme: aE = De − 83 αe Fe − 68 (1 − αe )Fe − 81 (1 − αw )Fw ,
Upwind Differencing Scheme: aEE = − 18 (1 − αe )Fe ,
aE = De + max(0, −Fe ), aW = DW + max(0, Fw ), aW = Dw + 86 αw Fw + 18 αe Fe + 38 (1 − αw )Fw ,
aN = Dn + max(0, −Fn ), aS = DS + max(0,
Fs ) aW W = − 81 αw Fw
Hybrid Differencing Scheme: aE = max 0, De − F2e , −Fe , aN = Dn − 38 αn Fn − 68 (1 − αn )Fn − 81 (1 − αs )Fs ,
aW = max 0, Dw + F2w , Fw ,

aN N = − 81 (1 − αn )Fn ,
aS = Ds + 86 αs Fs + 18 αn Fn + 38 (1 − αs )Fs ,
aN = max 0, Dn − F2n , −Fn ,


aSS = − 81 αs Fs
aS = max 0, Ds + F2s , Fs


aP = aE + aW + aN + aS + aEE+ aW W + aN N + aSS +
Power-Law Scheme:
Fe − Fw + Fn − Fs
aE = De max 0, (1 − 0.1kP ee k)5
+ max(0, −Fe ),

aW = Dw max 0, (1 − 0.1kP ew k)5
+ max(0, Fw ), We can follow the same procedure to derive the FVM equations
aN = Dn max 0, (1 − 0.1kP en k)5
+ max(0, −Fn ), for three dimensional problems.
aS = Ds max 0, (1 − 0.1kP es k)5 + max(0, Fs )

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 108 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 109 / 160

Convection and Diffusion False Diffusion Convection and Diffusion False Diffusion

False Diffusion I False Diffusion II

Compare the discretized equations for central differencing and

upwind schemes, T

Fw ΓW (ρu)w Hot Fluid Hot Fluid Hot Fluid

Central Differencing: aW = DW + = +
2 δxw 2 T T

ΓW + (ρu)w
2 δxw (ρu)w Cold Fluid Cold Fluid Cold Fluid
Upwind Differencing: aW = DW + Fw = +
δxw 2
a) b) c)
The diffusion coefficient in upwind scheme is equivalent to the
diffusion coefficient in central differencing scheme plus additional Figure 22: Effect of Diffusion on Temperature Distribution, a) Expected
term (ρu)w
2 δxw .
Result with Γ 6= 0. b) Expected with Γ = 0. c) Computational Result
with Γ = 0. (False Diffusion)
This additional term is called false diffusion and has unrealistic
implication.
The effect of false diffusion is significant in case of
multi-dimensional problems.
G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 110 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 111 / 160
 Convection and Diffusion One-way Coordinate System Convection and Diffusion Transient Convection and Diffusion

A coordinate system where the future event does not affect the
present or past events is called a one-way coordinate. A steady
one dimensional heat conduction in a plate is an example of
two-way coordinate.
Definitely time is a one-way coordinate.
Space coordinates can be one-way depending on the size of the
Peclet number.
Elliptic problems are always two-way coordinate problems while
parabolic problems are one-way coordinate problems.

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 112 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 113 / 160

Fluid flow Difficulties in Fluid Flow Modeling Fluid flow Difficulties in Fluid Flow Modeling

Difficulties in Fluid Flow Modeling I Difficulties in Fluid Flow Modeling II

In modeling convection and diffusion,

The velocity distribution in the flow field is assumed to be known. For two-dimensional problem, the momentum equations can be
The whole idea was to determine a scalar variable φ. obtained from equation (113) as follows,
However, in practical fluid flow problems, the velocity distribution x-momentum equation
is unknown. ∂ ∂ ∂ ∂ ∂u ∂ ∂u ∂P
(ρu) + (ρuu) + (ρvu) = (µ ) + (µ ) − + Su
The velocity is the main parameter to be determined in most flow ∂t ∂x ∂y ∂x ∂x ∂x ∂y ∂x
problems. (114)
y-momentum equation
As the energy equation was used as a transport equation for
temperature, the momentum equation is used as transport ∂ ∂ ∂ ∂ ∂v ∂ ∂v ∂P
(ρv) + (ρuv) + (ρvv) = (µ ) + (µ ) − + Sv
equation for velocity. ∂t ∂x ∂y ∂x ∂x ∂y ∂y ∂y
(115)
Consider the general transport equation for a quantity φ.
In these equations, S is the body force.
∂ ~ φ) = div(Γgradφ) + Sφ
(ρφ) + div(ρV (113)
∂t

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 114 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 115 / 160
 Fluid flow Difficulties in Fluid Flow Modeling Fluid flow Difficulties in Fluid Flow Modeling

Difficulties in Fluid Flow Modeling III Difficulties in Fluid Flow Modeling IV

For steady flow, equation (114) and (115) will be modified to, Consider the pressure distribution term applied to the
∂ ∂ ∂ ∂u ∂ ∂u ∂P two-dimensional grid shown in figure below.
(ρuu) + (ρvu) = (µ ) + (µ ) − + Su (116)
∂x ∂y ∂x ∂x ∂x ∂y ∂x

∂ ∂ ∂ ∂v ∂ ∂v ∂P
(ρuv) + (ρvv) = (µ ) + (µ ) − + Sv (117)
∂x ∂y ∂x ∂x ∂y ∂y ∂y
In addition to the momentum equations, the flow should satisfy
the continuity equation.
∂ ∂
(ρu) + (ρv) = 0 (118)
∂x ∂y
The above equations indicate some difficulties,
Non-linear quantities, (ρu2 , ρv 2 )
Strong coupling of equations (u, and v in all equations) Figure 23: Scalar and Staggered Control Volumes
No transport equation for pressure.
G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 116 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 117 / 160

Fluid flow Difficulties in Fluid Flow Modeling Fluid flow Difficulties in Fluid Flow Modeling

Difficulties in Fluid Flow Modeling V Difficulties in Fluid Flow Modeling VI

If we apply the finite volume method and integrate the pressure

gradient term in the x-momentum equation we get,
Z e
∂P
dx = Pe − Pw (119)
w ∂x

If Pe and Pw are obtained by linear interpolation from PE , PP and

PW ,
   
PP + PW PE + PP PE − PW
pe − Pw = − = (120)
2 2 2

The pressure in our volume of interest, PP is not there in the above

equation and this results in non-realistic solution, specially for
Figure 24: Checker - Board Pressure Distribution
oscillating type pressure distributions like that shown in Figure 24.

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 118 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 119 / 160
 Fluid flow Difficulties in Fluid Flow Modeling Fluid flow Difficulties in Fluid Flow Modeling

Difficulties in Fluid Flow Modeling VII Difficulties in Fluid Flow Modeling VIII

Integration of the x-momentum equation over the u-control

A remedy to avoid this non-realistic behavior is to use a staggered
volume gives,
grid.
Using staggered grid, ((ρuu)I,J − (ρuu)I−1,J )δy + ((ρvu)i,j+1 − (ρvu)i,j )δx =
Scalar quantities are evaluated at ordinary nodal points.     !     !
Velocity is evaluated at staggered grids centered at cell faces. ∂u ∂u ∂u ∂u
µ − µ δy + µ − µ δx−
∂x I,J ∂x I−1,J ∂y i,j+1 ∂y i,j
Referring to the staggered grid, the pressure gradient at the u
control volume can be written as, (P(I,J) − P(I−1,J) )δy + su δxδy
(123)
∂P PI,J − PI−1,J
= (121)
∂x δx ⇒ (FI,J uI,J − FI,J−1 uI,J−1 )δy + (Fi,j+1 ui,j+1 − Fi,j ui,j )δx =
Similarly the pressure gradient at the v control volume can be (DI,J (ui+1,J − ui,J ) − DI−1,J (ui,J − ui−1,J ))δy+
written as, (Di,j+1 (ui,J+1 − ui,J ) − Di,j (ui,J − ui,J−1 ))δx−
∂P PI,J − PI,J−1
= (122) (P(I,J) − P(I−1,J) )δy + su δxδy
∂y δy
(124)
G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 120 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 121 / 160

Fluid flow Difficulties in Fluid Flow Modeling Fluid flow Difficulties in Fluid Flow Modeling

Difficulties in Fluid Flow Modeling IX Difficulties in Fluid Flow Modeling X

I−1,Jρ I,J +ρ ρ +ρ
Fi,J + Fi+1,J 2 ui,J + I+1,J2 I,J ui+1,J
FI,J = (ρu)I,J = =
The variation of u at the u-control volume faces can be determined 2 2
by any of the differencing schemes (central, upwind, hybrid,· · · ). (ρI−1,J + ρI,J )ui,J + (ρI,J + ρI+1,J )ui+1,J
⇒ FI,J =
F and D in equation (123) can be obtained in the following way(to 4
be averaged at non-scalar nodes). Refer to Figure 25. (125a)
u N u (ρI−1,J + ρI,J )ui−1,J + (ρI−1,J + ρI,J )ui,J
J+1 (125b)
n FI−1,J =
v u v v j+1 4
u
e J
W w P E FI−1,j+1 + FI,j+1
v v s u v j Fi,j+1 = (ρv)i,j+1 = =
u 2
J-1 (125c)
S
v v v j-1 (ρI−1,J+1 + ρI−1,J )vI−1,j+1 + (ρI,J+1 + ρI,J )vI,j+1
u u
4
I-1 i I i+1 I+1 FI−1,j + FI,j
Fi,j+1 = (ρv)i,j = =
Figure 25: u Control Volume
2 (125d)
(ρI−1,J−1 + ρI−1,J )vI−1,j + (ρI,J−1 + ρI,J )vI,j
4
G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 122 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 123 / 160
 Fluid flow Difficulties in Fluid Flow Modeling Fluid flow Difficulties in Fluid Flow Modeling

Difficulties in Fluid Flow Modeling XI Difficulties in Fluid Flow Modeling XII

µI,J
DI,J = (125e) ai,J and anb are coefficients that depend on differencing schemes,
xi+1 − xi unb are velocities at neighboring nodes (i − 1, J), (i + 1, J),
µI−1,J (i, J − 1), and (i, J + 1),
DI−1,J = (125f)
xi − xi−1 Ai,J is the face cell area in x - direction, and
bi,J = Sbx Δx.
µI,J+1 + µI,J + µI−1,J + µI−1,J+1
Di,j+1 = (125g) Similarly for the v-control volume, the y-momentum equation
4(yJ+1 − yJ )
becomes,
µI,J + µI−1,J + µI−1,J−1 + µI,J−1
Di,j = (125h)
4(yJ − yJ−1 )
X
aI,j vI,j = anb vnb − (PI,J − PI,J−1 )AI,j + bI,j (127)
The velocities used in the calculation of F are obtained from
initial guesses or previous iterations. where,
aI,j and anb depend on F and D as in the case of u-control volume.
Finally, equation (124) can be rearranged in the following way,
vnb are velocities at neighboring nodes (I, j − 1), (I, j + 1),
X (I − 1, j), and (I + 1, j).
ai,J ui,J = anb unb − (PI,J − PI−1,J )Ai,J + bi,J (126)
F and D are obtained in the same manner as in u-control volume
where case.
G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 124 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 125 / 160

Fluid flow Pressure-Velocity Coupling Fluid flow The SIMPLE Algorithm

Pressure-Velocity Coupling I The SIMPLE Algorithm I

Acronym for Semi-Implicit Method for Pressure Linked Equations.

Is a guess and correct procedure.
The pressure field is unknown in most applications and the main The first step in SIMPLE algorithm is to assume pressure and
task in fluid flow modeling is how to relate the pressure to the velocity distributions and calculate new ones using equations (126
velocity distribution. and 127).
The main velocity-pressure coupling algorithms will be introduced
in this chapter:
X
ai,J u∗i,J = anb u∗nb − (PI,J
∗ ∗
− PI−1,J )Ai,J + bi,J (128a)
1 The SIMPLE Algorithm,
2 The SIMPLER Algorithm, ∗
X
∗ ∗ ∗
aI,j vI,j = anb vnb − (PI,J − PI,J−1 )AI,j + bI,j (128b)
3 The SIMPLEC Algorithm,
4 The PISO Algorithm. The pressure P ∗ in equations (128a) and (128b) is assumed. The
coefficients a are to be determined using assumed velocities u and
v. Velocities u∗ and v ∗ are to be determined from equations (128a)
and (128b).

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 126 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 127 / 160
 Fluid flow The SIMPLE Algorithm Fluid flow The SIMPLE Algorithm

The SIMPLE Algorithm II The SIMPLE Algorithm III

The next step is to correct the pressure and the velocities, The main feature of SIMPLE algorithm is to drop the summation
terms from equations (130a) and (130b) to simplify the solution
P = P∗ + P′ (129a)
procedure. Then these equations reduce to:
u = u∗ + u′ (129b) ai,J u′i,J = (PI−1,J
′ ′
− PI,J )Ai,J
(131a)
v = v∗ + v′ (129c) ⇒ u′i,J = di,J (PI−1,J
′ ′
− PI,J )
where P, u and v are the correct pressure and velocities and P’, u’, ′ ′ ′
aI,j vI,j = (PI,J−1 − PI,J )AI,j
and v’ are pressure and velocity corrections. (131b)
′ ′ ′
Subtraction of equations (128a) and (128b) from equations (126) ⇒ vI,j = dI,j (PI,J−1 − PI,J )
and (127), respectively, give us, The corrected equations now become,
X
ai,J u′i,J = anb u′nb − (PI,J
′ ′
− PI−1,J )Ai,J (130a) ui,J = u∗i,J + u′i,J = u∗i,J + di,J (PI−1,J
′ ′
− PI,J ) (132a)
∗ ′ ∗ ′ ′
vI,j = vI,j + vI,j = vI,j + dI,j (PI,J−1 − PI,J ) (132b)
X
′ ′ ′ ′
aI,j vI,j = anb vnb − (PI,J − PI,J−1 )AI,j (130b)

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 128 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 129 / 160

Fluid flow The SIMPLE Algorithm Fluid flow The SIMPLE Algorithm

The SIMPLE Algorithm IV The SIMPLE Algorithm V

Substitution for the corrected velocities gives transport equation
For a single scalar cell we have four staggered grids, and the for pressure correction, P ′ .
SIMPLE equations for ui+1,J and vI,j+1 will be, ′
aI,J PI,J ′
= aI+1,J PI+1,J ′
+ aI−1,J PI−1,J +
′ ′
(135)
ui+1,J = u∗i+1,J + u′i+1,J = u∗i+1,J + di+1,J (PI,J
′ ′
− PI+1,J ) (133a) aI,J+1 PI,J+1 + aI,J−1 PI,J−1 + b′I,J

∗ ′ ∗ ′ ′
vI,j+1 = vI,j+1 + vI,j+1 = vI,j+1 + dI,j+1 (PI,J − PI,J+1 ) (133b)
Exercise
The flow should also satisfy the continuity equation. Write expressions for coefficients aI,J , aI+1,J , aI−1,J I, J + 1, aI,J−1 and
The continuity equation is discretized using the scalar grid. for the source term bI,J in equation (135).

(ρuA)i+1,J − (ρuA)i,J + (ρvA)I,j+1 − (ρvA)I,j = 0 Once P ′ is determined from its transport equation, it will be used
⇒ (ρA)i+1,J ui+1,J − (ρA)i,J ui,J + (134) to correct the pressure and velocities.
(ρA)I,j+1 vI,j+1 − (ρA)I,j vI,j = 0 After finding the correct pressure and velocities, other scalar
quantities (such as temperature) will be determined from their
own transport equations.
G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 130 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 131 / 160
 Fluid flow The SIMPLE Algorithm Fluid flow The SIMPLE Algorithm

The SIMPLE Algorithm VI The SIMPLE Algorithm VII

The solution continues iteratively until convergence of all the flow

variables is attained.
The general solution procedure of the SIMPLE algorithm is given
in the following figure.

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 132 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 133 / 160

Fluid flow The SIMPLER Algorithm Fluid flow The SIMPLER Algorithm

THE SIMPLER ALGORITHM I THE SIMPLER ALGORITHM II

The SIMPLER algorithm is an improved version of the SIMPLE

algorithm (SIMPLE Revised).
Define Pseudo - Velocities as
The continuity equation is used to derive discretized equation for P
pressure (not pressure correction). anb unb + bi,J
ûi,J = (137a)
The momentum equations re-arranged, ai,J
P
anb vnb + bI,j
P
anb unb + bi,J Ai,J
ui,J = − (PI,J − PI−1,J ) (136a) v̂I,j = (137b)
ai,J ai,J aI,j
P
P
anb vnb + bI,j AI,j anb unb + bi+1,J
vI,j = − (PI,J − PI,J−1 ) (136b) ûi+1,J = (137c)
aI,j aI,j ai+1,J
P
anb vnb + bI,j+1
P
anb unb + bi+1,J Ai+1,J v̂I,j+1 = (137d)
ui+1,J = − (PI+1,J − PI,J ) (136c) aI,j+1
ai+1,J ai+1,J
P
anb vnb + bI,j+1 AI,j+1
vI,j+1 = − (PI,J+1 − PI,J ) (136d)
aI,j+1 aI,j+1
G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 134 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 135 / 160
 Fluid flow The SIMPLER Algorithm Fluid flow The SIMPLER Algorithm

THE SIMPLER ALGORITHM III THE SIMPLER ALGORITHM IV

The velocities now become,

ui,J = ûi,J − di,J (PI,J − PI−1,J ) (138a)

vI,j = v̂I,j − dI,j (PI,J − PI,J−1 ) (138b)

This pressure is used as an initial guess to be used in the SIMPLE
ui+1,J = ûi+1,J − di+1,J (PI+1,J − PI,J ) (138c) algorithm.
vI,j+1 = v̂I,j+1 − dI,j+1 (PI,J+1 − PI,J ) (138d) The rest solution steps in SIMPLER algorithm are similar to those
of SIMPLE algorithm.
The general solution algorithm is given below.
Substituting equations (138) in to the discretized continuity
equation gives transport equation for pressure as,

aI,J PI,J = aI+1,J PI+1,J + aI−1,J PI−1,J +

(139)
aI,J+1 PI,J+1 + aI,J−1 PI,J−1 + bI,J

The pressure distribution is obtained from equation (139).

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 136 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 137 / 160

Fluid flow The SIMPLER Algorithm Fluid flow The SIMPLEC Algorithm

THE SIMPLER ALGORITHM V SIMPLEC Algorithm

The SIMPLEC (SIMPLE Consistent) algorithm is also a modified

version of the SIMPLE algorithm.
SIMPLEC differs from SIMPLE only in velocity correction.

u′i,J = di,J (PI−1,J

′ ′
− PI,J ) (140a)
′ ′ ′
vI,j = dI,j (PI,J−1 − PI,J ) (140b)
Ai,J
� AI,j
�
where di,J = ai,J − anb and dI,j = aI,j − anb
The other steps of SIMPLEC algorithm are similar to that of
SIMPLE.

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 138 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 139 / 160
 Fluid flow The PISO Algorithm Fluid flow The PISO Algorithm

THE PISO ALGORITHM I THE PISO ALGORITHM II

The corrector step 1 of PISO algorithm is similar to the corrector

PISO stands for Pressure Implicit with Splitting of Operators.
step of SIMPLE.
Involves one predictor step and two corrector steps. It may be P ∗∗ = P ∗ + P ′ (142a)
seen as an extension of SIMPLE, with a further corrector step to
enhance it. u∗∗ = u∗ + u′ (142b)
Predictor step:
v ∗∗ = v ∗ + v ′ (142c)
Velocity components are predicted from assumed pressure values,
The first corrected velocities will be
X
ai,J u∗i,J = anb u∗nb − (PI,J
∗ ∗
− PI−1,J )Ai,J + bi,J (141a)
u∗∗ ∗ ′ ′
i,J = ui,J + di,J (PI−1,J − PI,J ) (143a)
X
∗ ∗ ∗ ∗
aI,j vI,j = anb vnb − (PI,J − PI,J−1 )AI,j + bI,j (141b) ∗∗
vI,j ∗
= vI,j ′
+ dI,j (PI,J−1 ′
− PI,J ) (143b)
Corrector step 1: Corrector step 2:

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 140 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 141 / 160

Fluid flow The PISO Algorithm Fluid flow The PISO Algorithm

THE PISO ALGORITHM III THE PISO ALGORITHM IV

To enhance the SIMPLE procedure PISO performs a second Subtraction of equations (144a) from (145a) and (144b) from
corrector step. (145b), respectively gives us,
X
ai,J (u∗∗∗ ∗∗
i,J − ui,J ) = anb (u∗∗ ∗ ” ”
nb − unb ) + Ai,J (PI−1,J − PI,J ) (146a)
X
ai,J u∗∗
i,J =
∗∗
anb u∗nb + Ai,J (PI−1,J ∗∗
− PI,J ) + bi,J (144a)
X
∗∗∗ ∗∗ ∗∗ ∗ ” ”
∗∗
aI,j vI,j =
X
∗
anb vnb + ∗∗
AI,j (PI,J−1 − ∗∗
PI,J ) + bI,j (144b) aI,j (vI,j − vI,j )= anb (vnb − vnb ) + AI,j (PI,J−1 − PI,J ) (146b)

The velocities can be corrected once more as Substitution of equations (146a) and (146b) in to the continuity
equation gives us the transport equation for the second pressure
correction,
X
ai,J u∗∗∗
i,J = anb u∗∗ ∗∗∗ ∗∗∗
nb + Ai,J (PI−1,J − PI,J ) + bi,J (145a)
” ” ”
∗∗∗
X
∗∗ ∗∗∗ ∗∗∗
aI,J PI,J = aI+1,J PI+1,J + aI−1,J PI−1,J +
aI,j vI,j = anb vnb + AI,j (PI,J−1 − PI,J ) + bI,j (145b) ” ”
(147)
aI,J+1 PI,J+1 + aI,J−1 PI,J−1 + b”I,J

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 142 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 143 / 160
 Fluid flow The PISO Algorithm Fluid flow The PISO Algorithm

THE PISO ALGORITHM V THE PISO ALGORITHM VI

Once the second pressure correction is obtained, the twice

corrected pressure will be,
′
P ∗∗∗ = P ∗∗ + p” = P ∗ + P + P ” (148)

The general solution pressure of PISO is shown below.

G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 144 / 160 G. S. Tibba (AASTU) Computational Fluid Dynamics May 25, 2021 145 / 160