Lecture Notes EG501V/EG503G – Computational Fluid Dynamics (AY 2023/24)

3. Discretization
The aim of these notes is to:
• Introduce the concept of discretization by an example
• Discuss numerical differentiation
• Discuss accuracy

Heat transfer example

A solid − uranium enriched material − rod of length L and cross sectional area A is fully
insulated, except for its end at x=0 (see Figure 3.1). At this end the rod exchanges heat with
the environment according to a heat transfer coefficient h. This implies that at x = 0 the heat
flux (energy stream per unit area) is φq'' = h (Tenv − Trod ( x = 0)) (note that φq'' is negative: heat
leaves the rod, i.e. streams in negative x-direction).

insulation
T=Tenv h heat production q
xx=0
=0 x=L

Figure 3.1

At some specific moment ( t = 0 ) when the rod and the environment have the same
temperature T0 , the fission reaction starts heating up the rod at a rate q Joules per second and
per meter. The environment is kept at T0 . We want to find out how much the rod heats up,
and how quickly it does so.

An unsteady, one-dimensional, microscopic energy balance over a slice of the rod between x
and dx leads to the following parabolic partial differential equation (PDE):

∂T q λ ∂ 2T
= + (3.1)
∂t ρ c p A ρ c p ∂x 2

with λ, ρ , c p the heat conduction coefficient, the density, and the heat capacity of the rod
material respectively.

Boundary conditions are:

• T ( x, 0) = T0 (initial condition)
∂T
• ( L, t ) = 0 (insulation on the right, Neumann condition)
∂x
∂T
• −λ (0, t ) = h (T0 − T (0, t )) (heat transfer on the left; a mixed boundary condition)
∂x

A typical strategy to solve this problem numerically goes through the following steps:
• From continuous functions to discrete numbers: we “discretize” in time and space.

1
• Differential equations then turn into sets of algebraic equations.
• These algebraic equations can be (efficiently) solved on computers.

i=0 i=1 i=2 ... i=N

x=0 L x=L
∆x =
N
Figure 3.2

Discretization in space is illustrated in Figure 3.2. The way the temperature depends on x is
mimicked by keeping track of temperature values in closely spaced points. Discretization in
time implies taking time steps of size ∆t.

Instead of T ( x, t ) we now write Tij with i representing the x-position ( x = i∆x ,i=0…N), and j
the time step number ( t = j∆t ). If we plan on making M time steps, j=0…M.

Initial condition: Ti , j=0 = T0 i = 0...N . Boundary conditions at x=0 and x=L: see below.

How to turn the PDE (Eq. 3.1) into a set of algebraic equations? For sufficiently small time
steps we write:

∂T Ti , j +1 − Tij
≈ (3.2)
∂t ij ∆t

∂ 2T ∂ 2T ∂  ∂T 
The term with of Eq. 3.1 is a bit more difficult. Realize that =   . Then
∂x 2
∂x 2
∂x  ∂x 

 
∂ 2T 1  ∂T ∂T 1  Ti+1, j − Tij Tij − Ti−1 j  Ti+1, j − 2Tij + Ti−1 j
≈  −  ≈ −  = (3.3)
∂x 2 ∆x  ∂x 1 ∂x 1 
 ∆x  ∆x
 ∆ x 
 ∆x 2
ij  i+ , j
2
i− , j 
2

The discrete version of Eq. 3.1 then reads

Ti , j +1 − Tij q λ  Ti+1, j − 2Tij + Ti−1, j 

= +   (3.4)
∆t ρc p A ρc p  ∆x 2 

This implies that from an “old” temperature distribution (at time j∆t ) we can calculate a new
distribution (at ( j + 1) ∆t ):

q∆t λ∆t  Ti+1, j − 2Tij + Ti−1, j 

Ti , j +1 = Tij + +   (3.5)
ρc p A ρc p  ∆x 2 

This is an update rule that can be easily implemented on a computer. Since for j=0 (the initial
situation) Ti , j=0 = T0 (i = 0...N ) we can calculate Ti , j =1 (i = 0...N ) and so on.

2
To be able to apply Eq. 3.5 for i=0 we need to know a Ti−1, j = T−1, j . For i=N we need to know
a Ti+1, j = TN +1, j . The points −1 and N+1 do not belong to the rod, they are virtual points. We
∂T
use them to apply the boundary conditions. For x=L (i.e. i=N) we have = 0 . This we
∂x
approximate as

∂T TN +1, j − TN −1, j
≈ = 0 → TN +1, j = TN −1, j (3.6)
∂x i= N , j 2∆x

∂T
At x=0 (i=0) the boundary condition is more complicated: −λ = h (T0 − T ) . We
∂x
approximate this as

T1, j − T−1, j 2h∆x

−λ = h (T0 − Ti=0, j ) → T−1, j = T1, j + (T0 − Ti=0, j ) (3.7)
2∆x λ

(beware of a potential source of confusion here: T0 is the initial temperature, Ti=0, j denotes
the temperature at i=0, at time instant j).

Numerical differentiation
After the above example, we now more formally discuss discretization. The basics boil down
to transforming derivatives into differences.

Definition of the derivative of a function f of x at x=x0:

df ( x ) f ( x ) − f ( x0 )
≡ f ′ ( x0 ) ≡ lim .
dx x = x x → x0 x − x0
0

In numerical analysis, the difference x − x0 stays finite (does not really go to zero).
f ( x ) − f ( x0 )
Derivatives are therefore approximated by differences: f ′ ( x0 ) ≈ .
x − x0

An important theoretical result (known as mean-value theorem): in the interval [ x0 , x ] there

f ( x ) − f ( x0 )
exists a point ξ for which exactly f ′ (ξ ) = with ξ generally unknown.
x − x0

Suppose we know the function f at positions xi −1 , xi , xi +1 that are spaced by a distance h (see
Figure 3.3). We distinguish three approximations to the derivative at x = xi :

df f −f
Forward difference: ≈ i +1 i
dx i h

df f − f i −1
Backward difference: ≈ i
dx i h

3
 df f −f
Central difference: ≈ i +1 i −1
dx i 2h

Figure 3.3

It looks like (see the figure) that the central difference approximation is the most accurate.
This we can ‘prove’ by means of Taylor expansions. Let’s write

h2 h3
f i +1 = f i + hf i′ +
2
f i′′ +
6
( )
f i′′′ + O h 4 (3.8)

h2 h3
f i −1 = f i − hf i′ +
2
fi′′ −
6
( )
fi′′′ + O h 4 (3.9)

( )
The “ O h 4 ” indicates that the leading order of the terms not included in the expansion is 4.
f i +1 − f i h ′′ h 2 ′′′ f −f
Equation (3.8) can be rewritten as f i′ =
h
− fi −
2 6
( )
f i − O h3 = i +1 i + O ( h )
h

This means that the forward difference approximation of the first derivative has an error
proportional to h. In numerical jargon: the forward difference approximation of the first
derivative is first-order accurate (accurate to h to the power one).

f i − f i −1 h ′′ h 2 ′′′ f − f i −1
′
Equation (3.9) can be rewritten as f i =
h
+ fi −
2 6
fi + O h3 = i
h
( )
+ O (h)

The backward difference approximation of the first derivative also has first-order accuracy.

h3 ′′′
Subtracting Eq. (3.9) from Eq. (3.8) gives: f i +1 − fi −1 = 2hf i′ + 2
6
( )
f i + O h 4 . This can be

f i +1 − f i −1 h 2 ′′′ f −f
′
rewritten as f i =
2h
−
6
( )
fi + O h 4 = i +1 i −1 + O h 2
2h
( )
Apparently, the central difference approximation of the first derivative has second-order (h2)
accuracy. For small h (which virtually always is the case), second-order accuracy is much
more accurate than first-order accuracy. Central difference being more accurate than forward
or backward difference could already be anticipated from Figure 3.3.

The second derivative (and its accuracy) follow from adding up Equations (3.8) and (3.9):

4
 f + f − 2 fi
( )
f i +1 + f i −1 = 2 f i + h 2 f i′′ + O h 4
→ f i′′ = i +1 i −21
h
( )
+ O h 2 . This equation is essentially

∂ 2T
the same as Eq. 3.3 where we needed a discretization for . The expression is second-
∂x 2
order accurate.

Given that the Navier-Stokes, continuity, and convection-diffusion equations as introduced in

Lecture Notes 1 only have first and second-order derivatives, expressions given above for the
discretization of first and second-order derivatives (including their level of accuracy) will
largely cover your needs.

Finite difference and finite volume discretization

The above approach to discretization goes by the name Finite Difference (FD) discretization.
One alternative to FD is Finite Volume (FV) discretization. The commonalities and
differences between FD and FV will be illustrated here for the steady-state version of the
convection diffusion equation (Eq. 1.13): ∇⋅ (uφ ) = Γ∇2φ . Compared to Eq. 1.13 a few
things were changed: (1) Steady state so that the ∂ ∂t was dropped; (2) we do not assume
incompressible flow so that we – in general – cannot write ∇⋅ (uφ ) = u ⋅∇φ ; (3) changed the
dependent variable c to φ . In two dimensions the equation reads

∂ ∂  ∂ 2φ ∂ 2φ 
∂x
( x)
φ u +
∂y
( y )  ∂x2 + ∂y 2 
φ u = Γ (3.10)
 

(which is an elliptic PDE).

In two dimensions and on a uniform grid with grid spacing hx and hy in x and y-direction
respectively (see Figure 3.4), an FD discretization based on central differences would read

(φux )i +1, j − (φux )i −1, j (φux )i, j +1 − (φ ux )i, j −1  φ + φ − 2φi , j φi , j +1 + φi , j −1 − 2φi , j 

+ = Γ  i +1, j i −21, j +  (3.11)
2hx 2hy  hx hy2
 

Figure 3.4

The philosophy behind FV discretization is somewhat different. The nodes i, j are now
placed in the centre of a volume (a surface in two dimensions), see Figure 3.5. We now set up
a steady state balance over the volume. The indices e, w, n, s are for the east, west, north, and

5
south face of the volume (as also indicated in the figure). The fluxes through the faces are of a
∂φ
convective and of a diffusive nature. As an example for the west-side: φ u x w hy − Γ hy .
∂x w

Figure 3.5

The steady state balance thus reads

∂φ ∂φ ∂φ ∂φ
0 = φ u x w hy − Γ hy − φ u x e hy + Γ + φ u y hx − Γ hx − φ u y hx + Γ hx (3.12)
∂x w ∂x e s ∂y s n ∂y n

It makes sense to discretize the differentials in Eq. 3.12 with central differences. Again west
∂φ φi , j − φi −1, j
as an example: −Γ hy = −Γhy . If we do this for east, north and south,
∂x w hx
substitute in Eq. 3.12, and divide by hx hy this results in

φux e φux φuy φu y  φi +1, j + φi −1, j − 2φi , j φi , j +1 + φi , j −1 − 2φi , j 

− w
+ n
− s
= Γ +  (3.13)
hx hx hy hy  h 2
hy2
 x 

We now need to make choices for how to express the terms at the faces of the grey volume
(east, west, north, and south) on the left-hand-side of Eq. 3.13 in the values of velocity and φ
on the nodes of the grid. One choice would be to take the average of the two neighbouring
φ u x w 1 (φ u x )i −1, j + (φ u x )i , j
nodes. With west as the example: =2 . Doing this for all sides will
hx hx
eventually lead to exactly the same result as obtained with the finite difference (FD) method
as given in Eq. 3.11.

One can think, however, of alternative ways of expressing face values into nodal values. One
example is “upwind” interpolation. For instance for west: φ u x w = (φ u x )i −1, j if u x > 0 (i.e. if
the “wind” is coming from the west), and φ u x w
= (φ u x )i , j if u x < 0 (i.e. if the “wind” is
coming from the east). The rationale behind upwind interpolation is that the value of φ in the
upstream direction is more important on the face of the finite volume than the value in the
downstream direction.

At this stage, a good question would be: why discuss FV if it – in the end – gives similar
results as FD. To answer this, it should be noted that the grids we used so far were very
simple: uniform (same x,y spacings everywhere), rectangular grids. If grids become more
complicated (see Figure 1.5 for a few examples) and/or get locally refined, FV methods offer

6
more flexibility than FD methods. This is largely due to FV methods standing closer to the
original balance equations that (eventually) underpin fluid dynamics.