Paper 3A3 The Equations of Fluid Flow and Their

Numerical Solution
Introduction
• A great many fluid flow problems are now solved by use of Computational
Fluid Dynamics (CFD) packages.
• One of the major obstacles to the good use of the resulting numerical solutions
is how to evaluate their accuracy, which in turn means understanding the
assumptions and limitations of the methods used.
• The only real way of doing this is to use the understanding of the basic
principles of fluid flow as covered in the rest of papers 3A1-4, together with
an understanding of the general limitations of numerical methods
• Here there is time only to consider a very brief introduction to the equations
and the methods used to solve them. A more extensive course, based mainly
on coursework, is available in the 4A2 module in the fourth year.

1.1
Outline of Course There may only be time to cover 2 of 4, 5 & 6.

1 Some numerical basics.

2 The equations of motion. Physical interpretation and their formulation as

differential equations or control volume (integral) equations. The problem of
turbulence. Hierarchy of levels of simplification.

3 Numerical approximations. Order of accuracy. Point stability analysis.

Classification as hyperbolic, elliptic and parabolic equations.

4 Example of a parabolic equation – the solution of the 2D boundary layer

equations. A simple turbulence model.

5 Example of an elliptic equation – the solution of the 2D compressible potential

flow equation.

6 Example of an hyperbolic equation – the solution of the 1D Euler equations.

1.2
Books
There are many books on CFD; most of them are far too detailed for this course.
The one perhaps most suitable – but still far too detailed – is:
“Computational Fluid Dynamics – The Basics with Applications”
D Anderson Jr. McGraw-Hill 1995

Other Info
5 lectures, 1 example paper. plus 1 (optional) DPO session (Thur Wk 9, 11-1).
This session replaces one of the 6 scheduled lectures. There will be no lecture on
Thu 8th March (Wk 8). There are a number of computer demos/MATLAB
exercises associated with this course. Occasionally, I will give out the odd
MATLAB script for you to play with. This will be done by downloading from
website:-
http://www.eng.cam.ac.uk/~tph

1.3
1 Some Numerical Basics
Before examining techniques for the partial differential equations applicable to
fluid flow, it is worth examining some simple ordinary differential cases –
problems that we can solve by other (analytic) means so that we can check the
results.
Problem
Solve
d2y
+ Ω2 y = 0 (1)
2
dt
subject to the boundary conditions:-
dy
y ( 0) = 0 ( 0) = Ω
dt
for which the solution is, of course, y ( t ) = sin Ωt .

1.4
Finite Difference Scheme
The natural way to represent a function of y
time numerically is to divide time into yn+1
discrete steps each of length ∆t, and to keep
track of the values of y at the end of each time yn
interval. Denote the value of y at t = n∆t by
yn (n=0,1,2,...). A numerical estimate of the
derivative of y at time n∆t can be found by
considering the slope of the line joining t
successive points.
dy ∆y yn +1 − yn ∆t
≈ =
dt ∆t ∆t

One can think of this as approximating the function over the interval by a simple
function (in this case a straight line) and using this approximation as a means of
estimating the derivative.
1.5
To estimate the value of the second derivative y yn+1
yn
(which is related to curvature) at a given time,
straight line approximations to y in each
interval are not good enough. The next level
of complexity would be to take a local yn−1
quadratic fit to y. A quadratic passing through
consecutive points can be shown to be

= α ( t − n∆t ) + β ( t − n∆t ) + γ
2
yapprox
t
yn+1 − 2 yn + yn −1
where α= ∆t ∆t
2 ∆t 2
yn+1 − yn −1
β= γ = yn
2 ∆t

1.6
This has
d2 yn+1 − 2 yn + yn −1
yapprox = 2α =
2
dt ∆t 2
The original differential equation has now become a difference equation. (A finite
difference approx)
d2y yn+1 − 2 yn + yn−1
+ Ω y
2
= 0 → + Ω 2 yn = 0
2 2
dt ∆t

1.7
Algorithm
The difference equation can be recast as
yn +1 = ( 2 − Ω2∆t 2 ) yn − yn −1 (2)

To solve this needs two initial conditions y0 and y1. When these are known, y2
and all succeeding values can be generated in turn. The initial conditions are
dy
y ( 0 ) = 0 and ( 0 ) = Ω which become
dt
y1 − y0
y0 = 0 and = Ω i.e. y1 = Ω∆t
∆t

1.8
The figure shows the case of Ω = 1 and ∆t = .1 . The circles are generated from
the difference equation, while the solid line is the analytic, exact solution. It is
difficult to see the error.

1.9
We have followed this procedure for a trivial equation. It applies without
modification to much more complicated equations which can not be solved
analytically. Turning the differential equation into a difference one is just as easy
for complicated equations.
E.g.

d2y yn+1 − 2 yn + yn −1
−t
+ e sin y 2
= t 2
→ − e− n∆t sin yn2 = n 2 ∆t 2
dt 2 ∆t 2

yn +1 ( )
= ∆t 2 n 2 ∆t 2 − e− n∆t sin yn2 + 2 yn − yn −1

1.10
Check
The very satisfactory agreement obtained in the above figure is far from the whole
story. Let us examine the solutions of the difference equation in more detail.
Put

yn = λ n ( )
λ 2 − 2 − Ω 2 ∆t 2 λ + 1 = 0
1
Ω 2 ∆t 2 Ω 2 ∆t 2 2
i.e. λ = 1− ± iΩ∆t 1 − when Ω∆t < 2
2 4
1
Ω 2 ∆t 2 Ω 2 ∆t 2 2
and λ = 1− ± Ω∆t −1 when Ω∆t > 2
2 4

leading to a general solution yn = Aλ1n + Bλ2n

1.11
Case (a) Ω∆t < 2

Note that for each of the λ's

2 2
Ω ∆t 2
Ω ∆t
2 2
Re ( λ ) + Im ( λ )
2 2
= 1− + Ω 2 ∆t 2 1 − = 1
2 4

Since they have modulus one, the first root can be written λ1 = eiα and the
second is λ2 = e −iα .
Applying the boundary conditions gives

iα −iα Ω∆t
A + B = 0 and Ae + Be = Ω∆t A = −B = giving
iα −iα
e −e

yn =
e iα
Ω∆t
−e −iα (e inα
−e −inα
) yn =
Ω∆t
sin α
sin nα

1.12
Now for small values of Ω∆t , sin α = Im ( λ1 ) ≈ Ω∆t which means α ≈ Ω∆t . The
form for the solution to the difference equation becomes
yn ≈ sin nΩ∆t = sin Ωt
agreeing with the exact solution of the differential equation.

1.13
As Ω∆t gets larger, yn is still periodic, but the amplitude of oscillation is no
longer unity and the period is different to Ω. The figure shows Ω = 1 and ∆t = .5

1.14
Case (b) Ω∆t > 2

In this case (exercise) that λ1 > 1 & λ2 < 1 . General solution yn = Aλ1n + Bλ2n
→ ∞ as n → ∞ . The figure shows the case Ω = 1 and ∆t = 2.000001.

1.15
Analysis
y
2π
(i) If Ω∆t 2, i.e. ∆t , then yn+1
Ω
yn
clearly
dy yn +1 − yn t
≈
dt ∆t ∆t

d2y yn +1 − 2 yn + yn −1
and ≈
2
dt ∆t 2
are good approximations and the solution
of the difference equation is close to the
solution of the differential equation.

1.16
 2π
If ∆t is comparable with , then y
Ω
dy yn +1 − yn
≈
dt ∆t yn

d2y yn +1 − 2 yn + yn −1 ∆t
and ≈ t
dt 2 ∆t 2
is pure fiction. Clearly ∆t has to be
chosen to be much smaller than the
wavelength (or time over which we have
significant variation) of the solution. In yn+1
this case, we know Ω. In most cases we
do not.

1.17
(ii) This raises the general question of just how accurate is the solution.
The way to assess it is by Taylor's theorem.
dy ∆t 2 d 2 y ∆t 3 d 3 y
y ( t + ∆t ) = y ( t ) + ∆t + + + ... (a)
dt t 2! dt 2 3! dt 3
t t

dy ∆t 2 d 2 y ∆t 3 d 3 y
y ( t − ∆t ) = y ( t ) − ∆t + − + ... (b)
dt t 2! dt 2 3! dt 3
t t
Thus
y ( t + ∆t ) − y ( t ) dy ∆t d 2 y ∆t 2 d 3 y
= + + + ...
∆t dt t 2! dt 2 3! dt 3
t t
and
y ( t + ∆t ) − 2 y ( t ) + y ( t − ∆t ) d2y ∆t 2 d 4 y
= + 2 + ... (c)
∆t 2 2
dt t 4! dt 4
t
1.18
The approximation for
dy yn +1 − yn
≈ is then accurate to first order in the timestep ∆t ,
dt ∆t
d2y yn +1 − 2 yn + yn −1
and that for ≈ is second order accurate.
dt 2 ∆t 2

All of the terms neglected on the right hand side of these equations are referred to
as the truncation error.

A glance at equations (a) and (b) above shows that since the terms on the right
hand side of (b) alternate in sign, while those of (a) all have the same sign, then if
we add or subtract (a) and (b) only alternate terms will survive. This means that,
if we use a combination which has symmetry about the central point (as in case (c)
above), the truncation error will be at least second order. Such symmetric
expressions are referred to as central differences.
1.19
We started with polynomial fits to a function in order to create approximations to
derivatives. Taylor series expansions, like (a) and (b) above, can also be used to
do this.

We can make the idea sketched above – that the size of ∆t has to be small in
comparison with the scale over which the solution to the differential equation
varies – more precise by thinking of the properties of the solutions obtained to the
differential and the difference equations.
The solution to the differential equation is of the form eiΩt giving

e (
iΩ t +∆t )
u (t + ∆t )
= = eiΩ∆t
u (t ) eiΩt

1.20
The corresponding solution to the difference equation is of the form λ n
1
Ω 2 ∆t 2 Ω 2 ∆t 2 2
λ = 1− + iΩ∆t 1 −
2 4

This means that

1
u n+1 Ω 2 ∆t 2 Ω 2 ∆t 2 2
= λ = 1− + iΩ∆t 1 −
n 2 4
u
For small ∆t,
u (t + ∆t ) ( Ω∆t ) 2
( Ω∆t ) 3
= eiΩ∆t = 1 + iΩ∆t − − i + ...
u (t ) 2! 3!
u n+1 Ω 2 ∆t 2 Ω3∆t 3
= λ = 1 + iΩ∆t − + i + ...
n 2 8
u
1.21
The sharp-eyed amongst you will ask about the boundary conditions that we
applied at t = 0. The approximation for the derivative at the origin was only first
order accurate, but this did not seem to be reflected in the accuracy of the solution
there. The answer to this is that we were lucky! In general, a first order accurate
boundary condition will lead to a first order accurate solution. In this case,
however, the truncation error in the derivative at t = 0, was of the form
y ( ∆t ) − y ( 0 ) dy ∆t d 2 y ∆t 2 d 3 y
= + + + ...
∆t dt 0 2! dt 2 3! dt 3
0 0

d2y
and it just happened that = 0
2
dt 0

1.22
(iii) Sometimes the properties of the difference equation bear no relationship to
those of the differential equation. For example, the number of independent
solutions of the difference equation is often not the same as the number of
independent solutions of the differential equation. Sometimes the difference
equation is unstable no matter how small ∆t . A good example of this type of
behaviour appears on the examples paper. See question 1.

Everything that has been said about the solution of ode’s her applies also to
PDE’s. As the course progresses, we shall see that ode’s are, in fact, relatively
benign!

1.23
Simple Heat Conduction Equation.

Let us consider a similar exercise for the scalar diffusion equation (e.g. heat
conduction equation in a bar, with the temperatures of the ends held constant).
∂u ∂ 2u
=ν
∂t ∂x 2
with the boundary conditions u ( 0 ) = u ( L ) = 0
At time t = 0, the initial value of u is given by
x L
10 for 0 < x <
L 2
u ( x,0 ) =
x L
10 1 − for < x < L
L 2

1.24
This time we discretise space and time and denote our
approximation to u at x = i∆x , t = n∆t

If we use a centred difference representation for

∂ 2u ∂u
and a one-sided representation for , we have:
∂x 2 ∂t
uin+1 − uin uin+1 − 2uin + uin−1
=ν (3)
∆t ∆x 2

Using the Taylor series expansions above, this means

uin+1 − uin ∂u ∆t ∂ 2u ∆t 2 ∂ 3u ∂u
= + + + ... = + O ( ∆t )
∆t ∂t i,n 2! ∂t 2 3! ∂t 3 ∂t
i ,n i ,n

uin+1 − 2uin + uin−1

∆x 2
=
∂ 2u
∂x i,n
2
+ 2
∆x 2 ∂ 4u
4! ∂x 4
+ ... =
∂ 2u
∂x 2 ( )
+ O ∆x 2
i ,n

1.25
 ( )
So that the truncation error is now O ∆t , ∆x 2 (or our approximation is first
order accurate in time and second order accurate in space).

uin+1 = uin +
ν∆t
∆x 2 (
uin+1 − 2uin + uin−1 ) for 1 < i < N x

N.B. NOT i = 1 or N x
Just as for the ode, we can solve equation (1) for uin +1 and iterate forward in time.
i.e. if we know uin for all i (i.e. 1 ≤ i ≤ N x at some time t = n∆t , this equation is
solved for each value of i in turn to generate uin +1 for all i (1 < i < N x ). We can
also find values for uin +1 for i = 1 and i = N x using the boundary conditions.

1.26
Figure 1 shows the numerical solution for two choices of timestep given by
ν∆t
(i) = .49
∆x 2
5

4.5

3.5

2.5

1.5

0.5

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

1.27
 6

ν∆t
(ii) = .51 5
∆x 2

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

ν∆t
We will show in a later lecture that the limit for stability is = .5
∆x 2

1.28
The following MATLAB (or OCTAVE) m file produced these plots.

nx=101;
for i=1:nx
x(i) = (i-1)*1./(nx-1);
if i < nx/2
u(i) = 10*x(i);
else
u(i) = 10*(1.-x(i));
end
end
hold on;
plot(x,u,'r','LineWidth',1.5);
al = .49;

1.29
nt = 2000;
nplot=100;
for n=1:nt
for i=2:nx-1
un(i) = u(i)+al*(u(i+1)-2*u(i)+u(i-1));
end
un(1) = 0;
un(nx) = 0;
if rem(n,nplot)==0
plot(x,un,'b','Linewidth',1.5);
end
u = un;
end
hold off;

1.30