Department of Aeronautical Engineering

Computational Fluid Dynamics

Chapter 4: Basics Aspects Of Discretization

1 Presentation title
☻Introduction

☻Discretization Methods

☻Introduction to Finite Difference Method (FDM)

☻Analysis of Numerical Schemes: Consistency, Stability,

Convergence

☻Introduction to Finite Element Method (FEM)

☻Introduction To Finite Volume Method (FVM)

☻Comparison between Finite Element Method vs Finite

Difference
2 Method vs Content
Finite Volume Method (FVM),
 Introduction
Analytical solutions of partial di fferential equations involve
closed-form expressions which give the variation of the
dependent variables continuously throughout the domain.
In contrast, numerical solutions can give answers at only
discrete points in the domain, called grid points.

• No. of points =  (infinite) • No. of nodes = 700 (finite)

• DOF per points = 6 ( 3 rotation &3 translation) • DOF per nodes = 6
• Total equations =   6 =  • Total equations = 700  6 = 4200
• Solving time =  (impossible) • Solving time = finite
♣ Therefore, the continuous domain of study should be replaced with finite set of points,
and the process is called discretization.
3 Introduction
• The spacing between the grids are
and y, they could be uniform or
variable along the study domain.
• However, the vast majority of CFD
applications involve numerical
solutions on a grid which involves
uniform spacing in each direction,
because this greatly simplifies the
programming of the solution, saves • If the gird points are
storage space and usually results in placed in an irregular
greater accuracy. fashion, it s unstructured
• Each grid points are identified by an gird.
index i which runs in the x- • If the gird as set in regular
direction, and an index j which runs manner it is structured
in the y-direction. gird.

4 Introduction
Structured grids are those whose control volumes can be
indexed by (i,j,k) for i = 1,..., ni, j = 1,.., nj, k = 1,..., nk,
Grids can be Cartesian
or curvilinear (usually
body-fitting).

Unstructured grids can accommodate completely arbitrary

geometries. However, there are significant penalties to be
paid for this flexibility, both in terms of the connectivity
data structures and solution algorithms.

5 Introduction
The purpose of the gird generation is to decompose the flow
domain into control volumes/elements. And will results:
☺ Cell vertices;
☺ Connectivity information: Precisely where the nodes
are relative to the vertices depends on whether the
solver uses, for example, cell-centred or cell-vertex
storage. Further complexity is introduced if a staggered
velocity grid is employed

6 Introduction
The shape of the element can be one of the following
depending on the applications and the computational
resource availability. (a) Triangle,
(b) Tetrahedron,
(c) Quadrilateral,
(d) Hexahedron,
(e) Prism, and
(f) Pyramid

7 Introduction
Similarly, the governing equation needs to be
discretized, i.e. governing equations (either PDE or
Integra form) of fluid dynamics will be replaced
with algebraic quotients, yielding a system of
algebraic equations which can be solved for the
flow-field variables at the specific, discrete grid
points in the flow.

8 Introduction
 Discretization Techniques

Finite Difference Finite Volume Finite Element

Method Method Method

Based Partial Based integral form

differential form of of governing
governing equation equation: finite
volume equations
Finite difference
equations:
truncation errors

Types of solutions:
Stability Analysis
Explicitly and
Implicitly
9 Introduction
 Finite Difference Method

10 Finite Difference Method (FDM)

In FDM, the partial differential equations are replaced by
computationally suitable algebraic difference quotients, i.e.
finite differences.
These finite difference representation of the derivatives are
based on Taylor’s series expansion.
Recall, for continuous function of x, i.e. f(x), with all
derivatives are defined at x, then the value of f at () can be
estimated from Tylor’s series expanded about point x as;

More higher order terms, more accuracy, the

result of the numerical analysis is more
closer to the exact solution.

11 Finite Difference Method (FDM)

Based on the Taylor’s series, the velocity of 2d flow expressed as;

The finite difference representation of derivatives are given by;

The approximate finite difference expression for derivatives will be

(first order accurate );

Or formally, This FD expression is

based on the information
to the right of grid point
(i,j); and is called first
order forward difference.
12 Finite Difference Method (FDM)
Similarly, Taylor’s series expression for

Or

Solving for ,

For this case the information left of grid point , are used and
is called first order rearward(backward) difference.

13 Finite Difference Method (FDM)

In most CFD problems, first order accuracy is not
sufficient, to find a finite difference quotient of second
order accuracy, evaluate the difference of the two
equations

Solving for,

This is second order central difference for the derivative

(∂u/∂x) at grid point (i, j).
In this case the information in both side of the grid point
(i,j) is used.
14 Finite Difference Method (FDM)
To obtain a finite-difference expression for the second
partial derivative , recalling from previous slides;

Substitute this into Taylor’s series expression to find;

This is a second-order central second difference for the derivative at

grid point (i, j).
15 Finite Difference Method (FDM)
Difference expressions for the y-derivatives are obtained in
exactly the same fashion. The results are directly analogous
to the previous equations for the x derivatives.
x y
Forward D

Backward
D

Central D

Central
2nd
difference

16 Finite Difference Method (FDM)

• Note that the central second difference given, can be
interpreted as a forward difference of the first
derivatives, with rearward differences used for the first
derivatives. Dropping the O notation for convenience, we
have

• By the same approached the finite difference quotient for

the mixed derivative () at grid point (i, j). For example,

17 Finite Difference Method (FDM)

Then, write the x-derivative as a central difference of the y-derivatives,
and then cast the y-derivatives also in terms of central differences.

Or

Many other difference approximations can be obtained for the above

derivatives, as well as for derivatives of even higher order. The
philosophy is the same.
18 Finite Difference Method (FDM)
 Finite difference quotient at the boundary
• Only one direction to go However, how do we obtain a result
• We can construct a finite which is of second-order accuracy?
difference approximation Why not with Central difference?
for ∂u/∂y at the boundary. How to drive the second order
• It is easy to construct a accurate difference at the boundary?
forward difference as; Reading assignment
The final output is

This is our desired second-order-

accurate difference quotient at the
boundary.
Is called one-sided differences, because they
express a derivative at a point in terms of
dependent variables on only one side of the
point.
19 Finite Difference Method (FDM)
20 Finite Difference Method (FDM)
21 Finite Difference Method (FDM)
22 Finite Difference Method (FDM)
NB: More higher order accurate quotients can be derived,
however we should note that;
☻ Higher order accurate difference quotients require more
grid points, result in more computer time for each time
and spatial steps.
☺ Higher order schemes may require a smaller number of
total grid points in a flow solution to obtain comparable
accuracy.
☺ Higher order difference schemes may result high quality
solutions such as captured shock waves that are sharper
and distinct.
☼ Degree of accuracy employed for one’s task, should
comprehend these points, and consider the application
and computation resource available.
23 Finite Difference Method (FDM)
To examine some aspects of FDM, Consider the following
unsteady, 1d heat conduction equation, with constant
thermal diffusivity equation as;
or

Let us replace the derivatives with

difference equations(i.e. discretization)

Substitute, these difference quotients,

we will find…. Next slide
24 Analysis of Numerical Schemes (FDM)
Ignoring the truncation error the difference quotient
expressed as;

This equation is a bit different than the original partial differential

reduced by increasing the number of grid points, i.e. ∆x→0 &

equation due to the truncation error, however, the difference can be

∆t→0, the truncation error will approaches to zero, when this the
case the FD representation of PDE is called Consistent.
25 Analysis of Numerical Schemes (FDM)
 Explicitly vs Implicitly
CFD techniques fall under two categories i.e. explicit and
implicit techniques. Consider from the previous discussion

26 Analysis of Numerical Schemes (FDM)

• It tend to marching solution in steps of time.
Assume that we know the dependent variable at all x at some instant
in time, say from given initial conditions.
From the above equation, we see that it contains only one unknown,
namely .
In this fashion, the dependent variable at time (t +Δt) can be obtained
explicitly from the known results at time t, i.e. is obtained directly
from the known values , and .

This is an example of an
explicit finite-difference
solution.

27 Analysis of Numerical Schemes (FDM)

Consider the heat conduction equation once again, with the
spatial derivatives rewritten as an average properties
between the time interval n and n+1, as

This differencing is called Crank Nicolson form

In this expression, the unknown is not only expressed in terms of the known
quantities at time index n, namely, and but also in terms of unknown quantities at
time index n +1, namely and
☺ Rather, the FD equationmust be written at all grid points, resulting in a system
of algebraic equations from which the unknown for all i can be solved
simultaneously. This is an example of an implicit finite-difference solution.
☺ Because they deal with the solution of large systems of simultaneous linear
algebraic equations, implicit methods are usually involved with the
manipulation of large matrices.

28 Analysis of Numerical Schemes (FDM)

 Explicitly Approach Implicitly Approach
☺ Relatively simple to set ☺ Stability can be maintained over much larger
up and program. values of Δt, hence using considerably fewer
☻ In terms of our above time steps to make calculations over a given
example, for a given Δx, interval of t. This results in less computer time.
Δt must be less than some ☻ More complicated to set up and program.
limit imposed by stability ☻ Since massive matrix manipulations are usually
constraints. In many required at each time step, the computer time per
cases, Δt must be very time step is much larger than in the explicit
small to maintain approach.
stability; this can result in ☻ Since large Δt can be taken, the truncation
long computer running error is larger, and the use of implicit methods
times to make to follow the exact transients (time variations of
calculations over a given the independent variable) may not be as accurate
interval of t. as an explicit approach. However, for a time-
dependent solution in which the steady state is
the desired result, this relative time-wise
inaccuracy is not important.

29 Analysis of Numerical Schemes (FDM)

 Errors
1. Discretization error: The difference between the exact analytical
solutionof the partial differential equation (and the exact (round-off
free) solution of the corresponding difference equation. From our
previous discussion, the discretization error is simply the truncation
error for the difference equation plus any errors introduced by the
numerical treatment of the boundary conditions.
2. Round-off error: The numerical error introduced after a repetitive
number of calculations in which the computer is constantly rounding
the numbers to some significant figure.
Let, A = analytical solution of the partial differential equation, D =
exact solution of the difference equation, N = numerical solution
from a real computer with finite accuracy
then,

we can write:
30 Analysis of Numerical Schemes (FDM)
 Stability Analysis
The numerical solution N must satisfy the di fference
equation; then

By definition, D is the exact solution of the di fference

equation, hence it exactly satisfies:

Evaluating the difference of the two;

31 Analysis of Numerical Schemes (FDM)

We now consider aspects of the stability of the di fference
equation, If errors ε are already present at some stage of the
solution of this equation (as they always are in any real
computer solution), then the solution will be stable if the ’s
shrink, or at best stay the same, as the solution progresses
from step n to n +1; on the other hand, if the ’s grow larger
during the progression of the solution from steps n to n +1,
then the solution is unstable. That is, for a solution to be
stable,

32 Analysis of Numerical Schemes (FDM)

• Depending on the properties of the method, solution
errors may be amplified or damped. An iterative
solution method is unstable or divergent if it fails to
yield a solution to the discrete set.
• When solving unsteady problems, we will use numerical
methods which compute the solution at discrete instants
of time, using the solution at one or more previous time
steps as initial conditions.
• Stability analysis allow us to determine whether errors
in the solution remain bounded as time marching
proceeds.
• It is possible to analyze iterative and time marching
methods using stability analysis. However, this is most
convenient for linear problems, and is usually too
difficult
33 forAnalysis
most ofrealistic
Numericalproblems.
Schemes (FDM)
 Convergence

Consistency

Stability

Please review in detail

34 Analysis of Numerical Schemes (FDM)

Examples: (on board discussion session)
1. Steady, one dimensional conduction
2. Transient, 1d conduction (explicitly and implicitly)
3. 1d convention
4. 2d dimensional Laplace's equation
5. 2d wave equation
6. NS-equations

35 Examples (FDM)
Examples: Mathlab code for-----------------------------Grid-----------------------------------------
transient 1d heat conduction
(Reference mathworks) %Setup grid
N=10; % number of nodes [Link]
clc DELTA_x=L/(N); % distance between adjacent nodes
clear all; [Link](m)
%----------------------------- x=0:DELTA_x:L; % position of each node [Link](m)
%--------------------------Stability
Input---------------------------------------- criterion-----------------------------
% Inputs DELTA_t_crit_N =
L=0.05; % wall thickness(m) DELTA_x*rho*cp/(2*(lambda/DELTA_x+alpha));
lambda=5.0; % conductivity(W/m-K) if DELTA_t > DELTA_t_crit_N
rho=2000; % density(kg/m^3) disp(' ')
cp=200; % specific heat capacity(J/kg-K) disp('Time step exeeds the limit');
return
T_ini=293.2; % initial temperature(K) end
T_inf=273.2; % external temperature(K) %----------------------Initial wall
alpha=500; % heat transfer coefficient(K) temperatures--------------------------- %Initial wall
%-----------------------------Time temperatures T(i,1)
stepping-------------------------------- for i=1:N+1
% Setup time steps T(i,1)=T_ini;
end
M=100; % number of time steps % Step trough time
t=40; for j=1:(M-1)
DELTA_t=t/(M); % time step duration(s) % Heat flux condition(q=n*(-k*dT/dx))[W/m^2]
for j=1:M T(1,j+1)=T(1,j)+2*lambda*(T(2,j)-T(1,j))*DELTA_t/
time(j)=(j-1)*DELTA_t; (rho*cp*DELTA_x^2);
end for i=2:(N)
T(i,j+1)=T(i,j)+lambda*(T(i-1,j)+T(i+1,j)-
2*T(i,j))*DELTA_t/(rho*cp*DELTA_x^2);
end
% Heat flux condition(q=n*(-k*dT/dx))[W/m^2] + heat
transfer coefficient(hout*(Tinf-T))[W/(m^2*K]
T(N+1,j+1)=T(N,j)+(2*lambda*(T(N-1,j)-T(N,j))/
36 Examples (FDM)
(rho*cp*DELTA_x^2)+2*alpha*(T_inf-T(N,j))/
(rho*cp*DELTA_x))*DELTA_t;
Examples: Mathlab code for transient 1d heat conduction
(Reference mathworks)

%-----------------------------
Plot-----------------------------------------
plot(time,T)
figure
plot(x,T)
figure
g=T(:,M);
plot(x,g)

37 Examples (FDM)
 Finite Volume Method
It is based on the approximate solution of an integral form
of governing equations.
The flow field or domain is sub-divided, as in the finite
difference method, into a set of non-overlapping cells that
cover the whole domain.
In the finite volume method (FVM) the term cell is used
instead of the term grid point used in Finite difference
method.
The conservation laws are applied to determine the flow
variables in some discrete points of the cells, called nodes.
These nodes are at typical locations of the cells, such as
cell-centres, cell-vertices or mid-sides.

38 Introduction to Finite Volume Method (FVM)

The fluid field/ study domain will be divided in finite
number of volumes as

2d dimensional domain

39 Introduction to Finite Volume Method (FVM)

After the flow the governing equation will applied for each
small volume.

• Surface and volume integrals approximated

• Interpolation used to express variable values at CV faces
in terms of nodal values
• It results in an algebraic equation per CV

40 Introduction to Finite Volume Method (FVM)

To demonstrate FVM; consider the 2D steady heat
conduction in rectangular geometry with source term. Study
domain is
Step 1: Grid generation

The governing equation for heat conduction

There are two alternative ways of setting up the nodal
equations are the weighted residual approach and the
physical approach.
41 Introduction to Finite Volume Method (FVM)
Using the weighted residual approach, the 2-D heat
conduction equation can be approximately satisfied by:
Integrate this equation
by part, we found

Or

43 Introduction to Finite Volume Method (FVM)

For
. a typical node P with neighbors
E,N,W,S (standing for east, north,
west and south etc.) and
corresponding control volume
boundaries in those directions
denoted by e,n,w,s etc., the heat
balance for the control volume can
be written as follows (for unit depth
in z-direction):

With all q’s are is the heat flux (per

unit area)
44 Introduction to Finite Volume Method (FVM)
In the implementation of the FVM procedure, the heat
fluxes are expressed in terms of the nodal temperatures (T,
etc. at the CV centers) using piecewise interpolation
around the control volume for the field variable
(temperature in this case). Thus, assuming temperature to
have linear variation between points E and P, the heat flux
can be evaluated as follows:
…

Then, the nodal equation for point P becomes:

or

where,
45 Introduction to Finite Volume Method (FVM)
This nodal equation will be written for each node in the
domain, as
Node 1
Node 2 …. to node n
We will have finite number of equations, in matrix form it
will have the following form and ready for computation

46 Introduction to Finite Volume Method (FVM)

The general procedure for FVM

Step 1: Grid generation

Step 2: Discretization

Step 3: Solution of equations

Example:

47 Introduction to Finite Volume Method (FVM)

Example:
Consider the problem of source-free heat conduction in an insulated rod whose ends are
maintained at constant temperatures of 100 °C and 500 °C respectively. The one-
dimensional problem sketched in figure below is governed by equation given below.
Calculate the steady state temperature distribution in the rod using Finite Volume Method
and compare the results with exact analytical solution. Thermal conductivity k equals 1000
W/m/K, cross-sectional area A is 10 x 10-3 m2, use dx = 0.1 m

(1)

48 Introduction to Finite Volume Method (FVM)

Solution
.
• Divide the length of the rod into five equal control volumes as shown in Figure below.
This gives dx = 0.1 m.

• The grid consists of five nodes. For each one of nodes 2, 3 and 4 temperature values to
the east and west are available as nodal values.
• Consequently, discretised equations can be readily written for control volumes
surrounding these nodes:

• The thermal conductivity (ke = kw = k), node spacing (x) and cross-sectional area (Ae =
AW = A) are constants. Therefore, the discretised equation for nodal points 2, 3 and 4 is

49 Introduction to Finite Volume Method (FVM)

.
Solution…
With,

• SU and SP are zero in this case since there is no source term in the governing equation

• Nodes 1 and 5 are boundary nodes, and therefore require special attention. Integration of
equation (1) over the control volume surrounding point 1 gives

(2)

Re-arrange,
(3)

50 Introduction to Finite Volume Method (FVM) By Mesfin Belayneh

. Solution…
• Comparing equation 3 with equation 4 it can be easily identified that the fixed
temperature boundary condition enters the calculation as a source term (S U + SPTP) with
SU = (2kA/x)TA and Sp = -2kA/x and that the link to the (west) boundary side has been
suppressed by setting coefficient aW to zero.

(4)

(3)

• Equation 3 may be cast in the same linear form to yield the discretised equation for
boundary node 1:

With,

51 Introduction to Finite Volume Method (FVM)

.Solution…
• The control volume surrounding node 5 can be treated in a similar manner. Its discretised
equation is given by

Re-arrange,

• The discretised equation for boundary node 5 is

Where,

52 Introduction to Finite Volume Method (FVM)

.Solution…
• The discretization process has yielded one equation for each of the nodal points 1 to 5.
• Substitution of numerical values gives kA/x = 100 and the coefficients of each
discretised equation can easily be worked out. Their values are given in Table

• The resulting set of algebraic equations for this example is

53 Introduction to Finite Volume Method (FVM)

.Solution…
This set of equations can be re-arranged as

• For TA = 100 and TB = 500 the solution can obtained by using, for example, Gaussian
elimination: • The exact solution is a linear distribution
between the specified boundary temperatures:
T = 800x + 100. Figure shows that the exact
solution and the numerical results coincide.

54 Introduction to Finite Volume Method (FVM)

. FDM,FEM,FVM
(i) Finite-Difference Method
Discretize the governing differential equations directly;
e.g.

Although many heat transfer problems may be solved using the

finite difference methods, as soon as irregular geometries or an
unusual specification of boundary conditions are encountered, the
finite difference technique becomes difficult to use.
55 Conclusion
. (ii) Finite-Element Method
Express the solution as a weighted sum of shape
functions , substitute into some form of the governing
equations and solve for the coefficients. e.g., for
velocity,

(ii) Finite-Volume Method

Discretize the governing control-volume equations
directly; e.g.

56 Conclusion
End of the chapter

57 Conclusion