Agenda of the Session

 Discretisation : Finite differences methods and difference equations

 Explicit and Implicit approaches
 Unsteady Problem -Explicit versus Implicit Scheme
 Errors and stability analysis.
 Time marching and space Marching,
 Reflection boundary condition.
 Relaxation techniques.
 Alternating Direction Implicit Method, (ADI)
 Second order Lax-Wendroff Method,
 Mid-point Leap frog method,
 Upwind scheme
 Numerical viscosity and Artificial viscosity.

1
Computational Fluid Dynamics
Obtain approximate solutions to complex problems by Numerical Method

Need to use a discretization method which approximates the

differential equations by a system for algebraic equations, which can
then be solved on a computer.

Accuracy of numerical solutions quality of discretization

 Components of a numerical solution
method

1. Mathematical Model:
Set of PDEs or integro-differantial eqs. and the corresponding
boundary conditions.
2. Discretization Method:
 Finite difference Method
 Finite volume Method
 Finite element Method
 Spectral (element) methods
 Boundary element

PDE’s (continuous) discrete equations (FDE's)

 [Link] &Basic Vector System

4. Numerical Grid: grid generation

• Structured (regular) grid

• Block structured grid
• Unstructured grid

Discrete locations at which the variables are to be calculated

are defined by the numerical grid, or mesh.
5. Finite Approximations: approx used in discretization process is
selected

e.g. Finite difference: approximations for the

derivatives at the grid points need to be selected

The choice influences:

 Accuracy of approximation and Dimensionality of the problem
 Developing the solution method
 Coding, debugging, speed of code

Compromise between simplicity,easy of implementation,accuracy

and computational efficiency has to be made and also Second
order methods in general are used.
6. Solution Method
Discretization yields a large system of linear/non-linear
algebraic equations.
Linear equations Algebraic equation solvers

Non-linear equations Iteration scheme used

i.e. linearize the equations & resulting linear

systems are solved
by iterative techniques.

Unsteady flows: methods based on

marching in time
7. Convergence criteria (for flows:
Steady iterative procedures)
usually by pseudo-time-
marching or equivalent iteration scheme
Need to set convergence for the iterative method. Accuracy & efficiency is important

Absolute convergence: a  a *  
(tolerance)
aa *
Relative convergence: 
a 
 Discretization methods (Finite Difference)

First step in obtaining a numerical solution is to discretize the

geometric domain to define a numerical grid

Each node has one unknown and need one algebraic equation,
which is a relation between the variable value at that node and
those at some of the neighboring nodes.

The approach is to replace each term of the PDE at the

particular node by a finite-difference approximation.

Numbers of equations and unknowns must be equal

Introduction to Finite Difference Method
 A partial derivative replaced with a suitable algebraic difference quotient is
called finite difference. Most finite-difference representations of derivatives
are based on Taylor’s series expansion.
 Taylor’s series expansion:
Consider a continuous function of x, namely, f(x), with all
derivatives defined at x. Then, the value of f atxalocation
x can be
estimated from a Taylor series expanded about point x, that is,

f 1 2 f 1 3 f 1 n f
f ( x  x)  f ( x)  x  2
x  
2
3
x   ... 
3
n
(x) n  ...
x 2! x 3! x n! x

 In general, to obtain more accuracy, additional higher-order terms must be

included.
 Numerical solutions can give answers at only discrete points in the domain, called grid
points.

 If the PDEs are totally replaced by a system of algebraic equations which can be solved
for the values of the flow-field variables at the discrete points only, then in this sense,
the original PDEs have been discretized.

 Moreover, this method of discretization is called the method of finite differences.

 Forward, Backward and Central Differences:

(1) Forward difference:

Neglecting higher-order terms, we can get

f
f ( xi 1 )  f ( xi )  ( ) i ( xi 1  xi ) 
( xi 1  xi ) 2  2 f
( 2 )i 
xi 1  xi  3 f
3
( 3 ) i  ... 
( xi 1  xi ) n  n f
( n ) i  ...
x 2! x 3! x n! x

f f ( xi 1 )  f ( xi ) f ( xi 1 )  f ( xi )
( )i   ………. (a)
x xi 1  xi xi 1

where xi 1  xi 1  xi
 (2) Backward difference:
Neglecting higher-order terms, we can get
f
f ( xi  1 )  f ( xi )  ( ) i ( xi  xi  1 ) 
( xi  xi  1 ) 2  2 f
( 2 )i 
xi  xi  1  3 f
3
( 3 ) i  ...  ( 1) n ( xi  xi  1 )
n
n f
( n ) i  ...
x 2! x 3! x n! x

f f ( xi )  f ( xi  1 ) f ( xi )  f ( xi  1 )
( )i   (b)
x xi  xi  1 xi ……………..

where xi  xi  xi  1

3) Central difference:
(a)-(b) and neglecting higher-order terms, we can get
f f ( xi 1 )  f ( xi  1 )
( )i 
x xi 1  xi  1
 Truncation error:
The higher-order term neglecting in Eqs. (a), (b), (c) constitute the
truncation error & can be written as

Forward: f f i 1  f i
( )i   o(x)
x x

Backward: f fi  fi 1
( )i   o(x)
x x

f f i 1  f i  1
Central: ( )i   o(x) 2
x 2x
 Second derivatives:
 * Central difference:
If xi xi 1 x , then (a)+(b) becomes

2 f f  2 fi  fi 1 2
( 2 ) i  i 1  o ( x )
x (x) 2
* Forward difference:
2 f f i 2  2 f i 1  f i
( 2 )i  2
 o(x)
x (x)

* Backward difference:
2 f f  2 fi 1  fi 2
( 2 )i  i 2
 o(x)
x (x)
 Mixed Derivatives:
 * Taylor series expansion:

f f (x) 2  2 f (y ) 2  2 f (x)(y )  2 f 3 3

f ( x  x, y  y )  f ( x, y )  x  y    2  o[( x ) , ( y ) ]
x y 2! x 2 2! y 2 2! xy

* Central difference:
 2 f  f  f i  1, j 1  f i 1, j  1  f i  1, j  1
   i 1, j 1  o[(x) 2 , (y ) 2 ]
 xy  i , j 4(x)(y )

 2 f  f  f i , j 1  f i 1, j  f i , j
* Forward difference:    i 1, j 1  o[( x) , (y ) ]
 xy  i , j (x)( y )

 2 f  f  f i  1, j  f i , j  1  f i  1, j  1
   i , j  o[(x) , (y ) ]
* Backward difference:  xy  i, j (  x )(  y )
EXPLICIT METHODS

u 2 u

t

y2
t
Δt
ui,j+1
Δy
ui+1,j+1
Present time
ui-1,j P ui+1,j
Previous time
ui,j

ui,j-1 y

y1 y2 yi ym FTCS
tn = nΔt (uniform time step)
L
  2u 
 u  
ui, j 1  ui, j   t  .... Forward-difference in time
 t  t  2 i, j 2
t
i, j
O(t )

 2ui, j
  2 u 
2
= ui1, j u i1, j Second order central-difference in space
 y y
 i, j
 2

 2ui, j
ui, j 1  ui, j
=
ui1, j u i1, j (A) Difference equation
t
y
After rearragement
t 2
ui, j 1 =ui, j
 i1, j
 2ui, j i1, j  (B)
y u

Difference  u
equation (A) is just an approximation for original PDE due to truncation error.
 2

Note: Truncation error for differential equation is O(t, (x) 2 )

Consistency of finite-difference representation of the PDE as Δx→0 &
Δt→0 differential equation reduces to original differential equations.
 i,j+1 Present time

Δt
Time-marching
Previous time direction
i-1,j i,j i+1,j

Properties at level (j+1) (present time) to be calculated from values at level j (previous
time) Remember that parabolic PDEs lend themselves to a marching solution, here
marching variable is time, t

Eq.(B) allows direct calculation of ui,j+1 from the known values on the RHS of eq.(B)

Explicit approach: each difference eq. contains only one unknown and therefore can
be solved explicitly for this unknown in a straight forward manner.

Comments on this method

• Explicit methods can be very unstable and should be used with caution
•
In general, whether the scheme is unstable or not, depends on the ratio, t / y 
2

For a given (Δy), Δt must be less than some limit imposed by stability constraints
• Relatively simple to set up and program
IMPLICIT METHODS:
In implicit method information at the boundaries at the same level does not feed
into the computation.

First-order backward difference approximation for time-derivative and

second-order central difference approximation for space-derivative

ui, j 1  ui, j ui1, j 1  2ui, j 1  ui1, j 1 (1)

t =
x  2

i-1,j+1 i,j+1
i+1,j+1 Present time-unknown

xi Previous time-known
i,j
BTCS Method
In equation (1): 3 unknowns ui-1,j+1 , ui,j+1 , ui+1,j+1
Thus, it results in a set of coupled finite difference equations at all grid
points
Rearrange equation (1)
t 
uj t  u j1  j j
2 1
i1  1 2 ux 2 1
i1   ui
x xt2  i i
c d
 i ai
 bi

3 unknowns in each FDE

Algebraic equations
Coefficient
matrix→Tridiagonal→Th
omas algorithm (n-1)
unknowns
Alternating Direction Implicit (ADI) Method

u   2

 2  u2 
u 2
Marching technique (1)
  
 x 
t y

u(t+Δt) will be obtained, in some fashion, from the known values of u(t)

Let’s use two-step process: first treat only x derivative implicitly

Step 1:
 ui,n j  2u n  un
u i,n1/
j 2
= u n1/
i1, 2  2un1/
i, j 2 u n1/ 2
 ui,nj i, j
(2)
j i1, j 1 i, j 1
t 
 
2
x
Equation (2) reduces to the tridiagonal form y 
2

(n+1/2) : intermediate time

b u n1/ 2  a u n1/ 2  c u n1/ 2  d (2')
i i1, j i i, j i i1, j i

where
t 
bi  c i t2 
2 2 , ai   1
  x 
d  un x t u n  2u n  un 
i i, j i, j i, j
i, j 1 1
2 y
Eq.(2’) yields a solution for ui, j
n1/ 2
for all i, keeping j fixed, via Thomas Algorithm
2
In Eq.(2’) first set j=1, and sweep in x (i=1,...,N) to find u i,n1/
j 2
1

Next, set j=2, and sweep in x (i=1,...,N) to find n1/ 2

u i, j
.... 2
M sweeps in x-direction
Need to use Thomas Algorithm M times
 t
At this place eq.(2’)
gives ui,j n+1/2

n+1/2
y
(1,M) (i,j)=(N,M)
Δt/2
Sweep j+1
direction j

x
(i,j)=(1,1) (i,j)=(N,1)
n
At the end of step1 (after M sweeps), the values of u at the intermediate time (t+Δt/2)
n1/ 2
are known at all grid points: i.e. ui, j is known at all (i,j)
Step 2:
Take the solution to the time (t+Δt), using the known values at time (t+Δt/2)
Again replace spatial derivatives with central differences, but this time treat y
derivative implicitly
 u n1  u n1/ 2
n1/ 2
i, j i, j
= u j
i1,  2u n1/
i, j 2  u n1/ 2
 u i,n1
j  2u n1
i, j  u n1
(3)
i1, j 1 i, j 1
t 
x  y 
2 2

Eq.(3) reduces to the tridiagonal form

 dj (3')
b u n1  a u n1  c u n1
j i, j 1 j i, j
j i, j 1

where 
t t 
bj  c j , aj 
2 2
 2
  2 x 
  1
x2  t
d j  i,n j1/ 
2
ui1, 
n1/2 
j 2ui,n1/2
j
 ui1,
n1/2
j 
u
 x  2
yields a solution for ui,n1 for all j, keeping i fixed, via Thomas Algorithm
j

j=1,....,M i=1
j=1,....,M
i=2
......... :
N times Thomas Algorithm
i=
N
 t
Present time

n+1
y
(1,M) i i+1 (i,j)=(N,M)
Δt/2

Sweep direction
x
(i,j)=(1,1) (i,j)=(N,1)
n+1/2
Remarks:
• Involves only tridiagonal forms
• Alternating direction implicit
• Scheme is second-order accurate
•General class of scheme involving splitting of two or more directions in an implicit
solution of the governing flow equation to obtain tridiagonal forms
• Approximate factorization
•For 3-D, see the scheme in Computational Fluid Dynamics for Engineers Vol.1
Klaus A. Hoffmann & S.T. Chiang pg.90
 Errors Involved in Numerical Solutions

 In the solution of differential equations with finite differences, a variety of

schemes are available for the discretization of derivatives and the solution
of the resulting system of algebraic equations.

 In many situations, questions arise regarding the round-off and truncation

errors involved in the numerical computations, as well as the consistency,
stability and the convergence of the finite difference scheme.
 Round-off errors: Computations are rarely made in exact
arithmetic. This means that real numbers are represented in “floating
point” form and as a result, errors are caused due to the rounding-off of the
real numbers. In extreme cases such errors, called “round-off” errors, can
accumulate and become a main source of error.
Truncation error: In finite difference representation of derivative with
Taylor’s series expansion, the higher order terms are neglected by truncating
the series and the error caused as a result of such truncation is called the
“truncation error”.

The truncation error identifies the difference between the exact solution of a
differential equation and its finite difference solution without round-off error.
Error

Round off
error

Optimum
Truncation

error

Number of ste* ps
Explicit v/s implicit Finite Difference Schemes

The transient (time-dependent) heat equation in 1D

In explicit finite difference schemes, the temperature at time n+1

depends explicitly on the temperature at time n. The explicit finite
difference discretization of above equation is

This can be rearranged in the following manner (with all quantities at

time n+1 on the left-hand-side and quantities at time n on the right-
hand-side)
Since we know 𝑇𝑖+1𝑛, 𝑇𝑖𝑛 and 𝑇𝑖−1𝑛, we can compute 𝑇𝑖,n+1.

The major advantage of explicit finite difference methods is that

they are relatively simple and computationally fast.

However, the main drawback is that stable solutions are obtained only
when

𝜕2y/𝜕𝑥2 are evaluated (at least partially) at the new time step.
In implicit finite difference schemes, the spatial derivatives

The simplest implicit discretization of heat equation in 1D is

This can be rearranged so that unknown terms are on the left and known
terms are on the right

Where.

The main advantage of implicit finite difference methods is that there

are no restrictions on the time step, which is good news if we want to
simulate geological processes at high spatial resolution.

Taking large time steps, however, may result in an inaccurate solution.

 Finite Difference Truncation Error

From Taylor’s Series Expansion

f  2 f  x 2  n f  x n
f x  x   f ( x)  x   2     n 
x  x  i, j 2  x  i, j n!

f ( x) sin 2 x

at : x 0.2 f ( x) 0.9511

f (0.22) ???? x 0.02

f
f x  x   f ( x)  x
x

f (0.22)  f (0.2)  2 cos[ 2 (0.2)](0.02) 0.9899

Exact solution for f (0.22) 0.9823 Error 0.775 percent

4. The Lax-Wendroff method:
The L-W method is derived from Taylor series expansion of the dependent
variable as follows

u t +  2 u2 t
u  x, t  t  =u  x, t t +O t 
3
t
2 2! (1)
+
n1  u t +  2
u2 t (2)
uor interms of
t indicest +O t
Now 2!2 the model eq. 3
u n i consider
i

u u (3)
 a
t x
  u  2 2 u (4)
tu2 a tx
2
u a x  t   a x2

Substituting (3) & (4) into (2) produces


  a u  t + t  2
u
u n1
i  i
22
2
a x2
 x 
 
un
 Use central differencing of second order for the spatial derivatives

n  un
n1
ui n
 a ui1 i1 1 2 2
n
ui1  2u ni 2 u ni1
2x t + a t
i
u 2
x 

O t  , x 
Lax-Wendroff method 2 2

Stability analysis shows  explicit method is stable for C ≤ 1

Implicit Formulations
1. Euler’s BTCS method:

u n1
i  u ni a  i1n1 ui1n1
t   2x  u 
1 n1 n1 1 n1
Cui1  u i  Cui1  uin
2 2
 t, x   TDMA
2

 

2. Crank-Nicolson method:
u n
u n1
i
u i n n1
a  ui1  u n1
i1
u n
   i12x i1 
t 2  2x 
O t  , x  TDMA
2 2


SOLUTION OF EULER’S EQUATIONS
Lax-Wendroff Technique
• Explicit
• Particularly suited to marching solutions: hyperbolic & parabolic eqs.
Example: Time-marching solution of an inviscid flow using unsteady Euler eqs.
For unsteady, 2-D inviscid flow eqs. (HYPERBOLIC IN TIME)

Continuity:     u  u    v  v   (1) no boundary forces

t  
x x y y 
u  u  v u 1 P 
x-mom:   (2)
  u y  x
t  x 
v  v
 v v 1 P 
y-mom:   u  
t (3)  x y  y

Energy: e  e e P u P v  (4)
  u  v  
t  x y  x  y 
Cv: specific heat at constant volume
e: internal energy, we have thermodynamic relation
e=e(T,P)
For perfect gas with constant specific heat e(T)=cvT
Eqs. (1) to (4) are hyperbolic with respect to time
Taylor series expansion in time
n
t
i,n1j    
n n
  2   2
 2!
+ ……….......... (5)
   t +  t
i, j 2
  t  i, j  i, j
If flow field at time
 level n is known,
Eq.(5) gives the new flow field at time (n+1)→t+∆t

n n
if    &     i,n1j
2
are found can be calculated explicitly, from eq.(5)
  t  
 t i, j
2
 i, j
Analogous Taylor series for all other dependent variables can be written as follows:

  2u  t +...
n n
 u  (6)
u n1
    t +  t 2!
 t  i, j 2
 i, j
un i, j 2

  2 v2  t +...
n n
i, j  v  (7)
v n1
 v    t +  t
n
2!
i, j
 t  i, j 
 i, j
i, j
2

e n1
e n
 e 
n
  2 e 
n t +... (8)
i, j   t +  t 2
i, j
 t  i, j 
 i, j
 2!
2
Using spatial derivatives (second-order central dif.) from eq.(1)
 n  un   i1,
n  vi,n j 1   i,n j 1 
  
n
 u i 1, j  n j
i1, j
 n v ni, j 1  vi,n j  i,n j 1
(9)
 i,n j i 1, j
 u ni, i, j 
 t    j 2x 2x  2y 2y 
 i, j
In (9), all quantities on RHS are known.

n
Differentiate eq.(1) with respect to time   2 2  
 t ?
 i, j
n  2 u u   2   u  2 v v   2   v
 2       u    v 
 2  xt
 x t xt x t yt y t yt y t  (10)
 t 
 i, j

2 u
=?
x
t
Differentiate eq.(2) wrt x;

1 P
 2u  2 u  u  2  2u u (11)
 u 2       2 x x
xt x x 
v
v are1 expressed
In eq.(11) all terms on RHS  P xy 2
as second-order, central dif. eqs. at time

level n:
y
 n n
 2u n  un i1, j 2
n u i1,  u i1,  u
  u   u
2 i, j  
i, j j i1, j  j 
 2x 
 xt  i, j
 x 
2

ni1, j  ui1,
 u u  ui,n j 1 i1,  vi1,
n n n n
u j 1 i1, j 1 i1, j
 u ni, j v n j

v i,n j 1 1 1
2y
j
2x
4 x y
 2P n  P n
j i1,  i1, j
n n
1 Pi1,
n i, j 1 n j
Pi1,  P i1,
n

 (12)
j
 j 2 2
i, j 
i1, j 2
 n
i, j
x  2

x  x 

In eq.(12) all terms on RHS are known.
2 
Continuing with the evaluation of eq.(10), a number for is found by
xt
differentiating eq.(1) wrt x & replacing all derivatives on RHS with second-
order central differences, similar to eq.(12).
2 v
differentiate eq.(3) wrt y
yt
2 
differentiate eq.(1) wrt y
y
t
u v central difference; eqs.(2) & (3) , respectively.
t , t
Finally;

n
  2 2  is calculated from eq.(10)
 t
i, j
From eq.(5) i,n1j is known

For the remaining flow-field variables, ui,n 1j , vi,n 1

j , i, repeat the above procedure
en 1 j

ui,n1
j  f un
i, j , , u n
i, j 1 i, j 1 , u
i1, j ,
n
i1, j
un un

Remarks on Lax-Wendroff
•second-order accuracy in both space & time
•Algebra is lengthy
 n1
i
 stable solution
ni 1
ebt  cos(k x)  iC sin(k
x)
m
m t
where C  ax
the stability requirement is ebt  1
t
Ca 1 Courant number (CFL condition) (6)
x
important stability requirement for hyperbolic eqs.

3. Midpoint Leapfrog method (CTCS)

second-order central differencing for both time&space derivatives
 un
u n1
i
 u n1
i
ui1n i1
2t  a O (t) 2 ,(x) 2 
2x (7)Method is stable when C ≤ 1
  Two sets of initial values are required to start the solution,
 a starter scheme is needed (affects the order of accuracy of the method),
 Large increase in computer storage.

In numerical mathematics, relaxation methods are iterative methods for

solving systems of equations, including nonlinear systems

Relaxation methods are important especially in the solution of linear

systems used to model elliptic partial differential equations, such as
Laplace's equation and its generalization, Poisson's equation.

These equations describe boundary-value problems, in which the solution-

function's values are specified on boundary of a domain; the problem is to
compute a solution also on its interior.

Relaxation methods are used to solve the linear equations resulting from a
discretization of the differential equation, for example by finite differences
Iterative relaxation of solutions is commonly dubbed smoothing because
with certain equations, such as Laplace's equation, it resembles repeated
application of a local smoothing filter to the solution vector
 Case Study: 1 Dimensional Unsteady State Conduction
Study the temperature distribution across the engine block
wall which is exposed to external heat source and its
outside temperature kept 900k Assumptions:
Solution • No heat generation inside the wall
(study domain).
• Initial condition: the temperature
at the beginning is 1
900k 900k • One dimensional conduction;
• An explicit approach will be
employed

Governing equation Discretization the domain

(mesh generation)
BC/IC

Numerical problems
. Discretize the governing equation
T ( x ,t +dt ) - T ( x ,t ) ( T ( x- dx, t ) - 2 T ( x, t ) + T ( x+dx ,t ) )
=ᾳ
dt dx × d 𝑥

Time & Spatial

discretization

Write the discretised equation for each node as;

ᾳ ×d t
T ( x , t + dt )=T ( x , t ) + ( T ( x -dx , t ) − 2 T ( x , t ) + T ( x + dx , t ))
dx × d 𝑥

Note dx is spatial step in the domain and dt is time step of the simulation.
Recall the spatial index i, which runs upto N (i.e. the total number of division), then
if the total thickness of the wall is W, then dx = W/N. Similarly, for the time
marching let j be the index which runs upto M(the total number of times the
equation will be solved for each node), if the simulation time considered is T, then
dt=T/M

Numerical problems
Therefore, the nodal temperature at time t+dt (j+1), at node
i is given by;
T ( i , j+ 1 )=T (i , j) + C( T ( i- 1 , j) − 2 T ( i , j ) + T ( i + 1 ,j ))
Where: , i=1:N; j=1:M

This equation will be written for each node in the study domain

Node 1:

Node 2:

…
Node N:

Numerical problems
 MATLAB Program for FDM

clear % We are clearing the work space and command window

clc
clf
% domain descritization
alpha = 0.05; % material thermal conductivity
xmin = 0; % position of the first node
xmax = 0.5; % position of the final node, i;e. the thinness of the bar
is 0.5m
N = 100; % Total number of nodes
dx = (xmax-xmin)/(N-1); % spatial step
x = xmin:dx:xmax; % location of nodes
dt = 4.0812E-5; % time step
tmax = 1; % the total simulation time in seconds
t = 0:dt:tmax; % fractions of intermediate simulation time
% problem initialization
phi0 = ones(1,N)*300; % initial conditions of all the nodes
phiL = 900; % The temperature of the left node/boundary
phiR = phiL; % The temperature of the right
node/boundary
% stability check
r = alpha*dt/(dx^2); % for stability, must be 0.5 or less
% solving the temperature equation for each node
for j = 2:length(t) % for time steps
phi = phi0; % initial temperature value at t=0
for i = 1:N % for space steps
if i == 1 || i == N
phi(i) = phiL; % for both boundary nodes that is i=1

and i =N/11
else
phi(i) = phi(i)+r*(phi(i+1)-2*phi(i)+phi(i-1)); % For the internal gird
points
end
end % ending the calculation
phi0 = phi;
plot(x,phi0)
xlabel('node location')
ylabel ('temp')
shg
pause(0.05)
end
 Shock Wave Discontinuities
 Shock waves are mathematically treated as discontinuities, but it was
quickly recognized this would cause problems for any numerical solution.

 Of course, a shock wave is not a true physical discontinuity, but a very

narrow transition zone whose thickness is on the order of a few molecular
mean-free paths.

 Application of the conservation of mass, momentum, and energy conditions

across a shock wave requires that there be a transformation of kinetic
energy into heat energy. Physically, this transformation can be represented
as a viscous dissipation,

 By introducing an unphysically large value of viscosity, investigators were

able to thicken shock transition zones to where they could be resolved
computationally.

 This artificial increase in the value of viscosity became known as an

artificial viscosity.
Viscosity Dissipation
 If the viscosity isn’t large enough, velocity oscillations about the correct

mean velocity are observed to develop behind a shock. These oscillations

can be interpreted as a macroscopic version of heat energy, i.e., fluctuating
kinetic energy in place of fluctuating molecular energy.

Artificial Viscosity
 The proper formulation and magnitude needed for an artificial viscosity has

undergone many refinements over the years and includes tests to apply this
viscosity only in regions undergoing strong compression and with
magnitudes that are various functions of the first and/or second power of the
compression rate.
 Numerical viscosity arises from discrete approximations to the
momentum advection terms in Eulerian equations, or from re-zoning operations
used in Lagrangian formulations.

The origin of the effect is the use of a homogenizing assumption in the

elements or control volumes underlying the approximation scheme.

For instance, when momentum is exchanged between neighboring

elements through a convective flux, the resulting contributions in a given
element are combined with the momentum already there to arrive at a new
value of average momentum for that element.
This combining or homogenizing process introduces a smoothing
effect. When another step to advance in time is taken, this new value is passed
on to the next element in the direction of flow. Repetition of this smoothing
operation over the many steps needed to carry a solution forward in time
contributes to a “diffusion” of momentum in the direction of flow.
Strictly speaking, the numerical diffusion does not behave like a true viscous
diffusion because it is associated with fluid convection and does not possess the
correct stress-versus-strain-rate dependency associated with a real viscosity.

For example, numerical diffusion does not satisfy Newtonian relativity

because it depends on the choice of computational grid, which is an absolute
reference frame for the numerical approximations.
Also, because the amount of numerical diffusion is proportional to the
velocity of flow through a grid, it does not have the rotational symmetry
possessed by a real viscosity.

Numerical Approximations
The difficulty in developing such schemes is that some smoothing must
always be incorporated into a numerical solution to keep it computationally
stable and to smooth out dispersion errors.
 Dispersion errors are those errors that arise because components of a
solution having different grid resolution requirements may propagate through
the grid with slightly different speeds. Whenever this occurs, unphysical
oscillations develop in the solution where these components reinforce or cancel
one another.

The trick is to develop approximation schemes that remain accurate (i.e., have
a minimum of numerical smoothing)

and

at the same time are robust (i.e., have sufficient numerical smoothing to
compensate for dispersion errors and to remain computationally stable for a
wide range of applications).
 Software and Resources
• CFD software was built upon physics, modeling, numerics.
• Two types of available software
– Commercial (e.g., FLUENT, CFX, Star-CCM, COMSOL)
– Research (e.g., CFDSHIP-IOWA, U2RANS)
• More information on CFD can be got on the following website:
– CFD Online: [Link]
– CFD software
• FLUENT: [Link]
• COMSOL [Link]
• CD-adapco: [Link]
– Grid generation software
• Gridgen: [Link]
• GridPro: [Link]
– Visualization software
• Tecplot: [Link]
• Fieldview: [Link]
 ?
Questions are Welcome
Thank You