1.

GOVERNING EQUATIONS OF FLUID FLOW

Basic equations of fluid flow include continuity, momentum and energy equation. Depending
on the type of flow compressible/incompressible, laminar/turbulent, steady/unsteady the
governing equations are assumed.

Continuity Equation
Principle of conservation of mass: Within the region, the difference between the rate of
increase of the fluid mass contained within the region must be equal to the difference
between the rate at which the fluid mass enters the region and the rate at which the fluid
leaves the region.

Consider an elementary rectangular parallelepiped with sides of length δx, δy and δz.
u, v, w are the velocity components in x, y and z directions.

Δz

Δy

Mass contained in parallelepiped = ρ Δx Δy Δz Δx

Assuming the net flux is leaving,


Rate of change of mass within control volume = - (ρΔxΔyΔz) .....................................(1)
t
Mass of flux entering from left face/unit time =  x u x yz
Mass of flux entering from right face/unit time =  x  x u x  x yz
Net flux in x-direction =  x  x u x  x yz -  x u x yz .............(2)
Similarly, Net flux in y-direction =  y  y v y  y xz -  y v y xz
Net flux in z-direction =  z z wz z xy -  z wz xy
By Taylor series expansion,
( u ) 1  2 ( u ) 2 1  3 ( u ) 3
( u) x  x  u x  x  x  x  ...........
x 2! x 2 3! x 3

1
Neglecting the higher order terms as Δx, Δy, Δz 0

 ( u )
( u) x x  u x  x
x
 ( u )
Net flux in x-direction = [ u x  x  u x ] yz
x

 ( u )
xyz .................................................... (3)
=
x
 ( v)
Similarly, Net flux in y-direction = xyz .......................... (4)
y
 ( w)
Net flux in z-direction = xyz ............................................. (5)
z
Net flux = Rate of change of mass within control volume
(3)+(4)+(5) = (1)
 ( u )  ( v)  ( w) 
[   ] xyz =  xyz
x y z t

  ( u )  ( v)  ( w)
    0 3-D Compressible Continuity Equation
t x y z
Writing the above equation in vector notation


t
 
   .V  0


t
 
 V .    .V  0 
or
D
Dt
 
  .V  0

D  (  ) (  ) (  )
Where  u v w
Dt t x y z

Case I For steady state

Continuity equation is   .V  0 

2
Case II For steady incompressible flow

Continuity equation is V
.  0 .This equation is valid for uniform cross sectional area.
1-D Continuity Equation

Considering a stream tube AV  AV S
s


AV  AV S
s
S

AV

Mass of fluid entering control volume per unit time = AV


Mass of fluid leaving control volume per unit time = AV  AV S
s

Net mass leaving control volume = AV S ............................................(1)
s
Mass of fluid in the control volume = AS

Rate of decrease with time = - AS ...................................................(2)
t

According to law of conservation (1) = (2)


AV S    AS
s t


A   AV   0 1-D variable area Continuity equation
t s

Case I Steady flow


The continuity equation becomes AV   0
s
or AV  constant

Case II Steady incompressible

AV = constant

3
Acceleration of a Fluid Particle
a = Rate of change of velocity

du dv dw
ax  ; ay  ; ax 
dt dt dt

u = f (x, y, z, t)

du u x u y u z
ax    
dt x t y t z t

u u u u
ax  u v w 
x y z t

u u u u
ax  u v w
t x y z

Similarly

a y  v  u v  v v  w v
t x y z

w w w w
az  u v w
t x y z

a
DV V
Dt

t
 V . V 

Local Conventional
acceleration acceleration

Where V  V ( x, y, z , t )  u ( x, y, z , t )iˆ  v( x, y, z , t ) ˆj  u ( x, y, z , t )kˆ

DV V V V V
 u v w
Dt t x y z

Du Dv Dw
ax  ; ay  and az 
Dt Dt Dt

4
Force:
A Fluid in motion is subjected to several forces which results in the variation of the
acceleration and the energies involved in flow phenomenon in the fluid.

Body forces Gravity

Force α Volume of the body
Others-CG forces, Magnetic
electromagnetic etc.

Force Pressure
Surface forces
Force α Surface area Others-Shear, tangential, force of
compressibility, force due to turbulence etc

Line forces Surface tension

Force α Length
F x  Ma x
Ma  Fg  F p  Fv  Ft  Fs  Fc

Fp (pressure) Exists when pressure gradient exists

Fg (gravity) is due to weight of fluid
Fv (viscous) due to Viscosity of fluid
Ft (turbulent) due to turbulence of the flow. Causes additional stresses called
“Reynolds Stresses”.
Fs (surface tension) due to cohesive property of fluid mass. When depth of the fluid is
extremely small, Fs is considerable.
Fc (compressibility) due to elastic property of fluid

Most of the fluids in motion, surface tension and compressibility forces are not significant.
Neglecting these terms
Ma  Fg  F p  Fv  Ft ........................Reynolds equation of motion.

For laminar of viscous flows, the turbulent forces become less significant. Then
Ma  Fg  F p  Fv ..........................Navier stokes equations.

Further, neglecting viscous forces

Ma  Fg  F p ............................Euler’s equation of motion.

5
Euler’s Equation of motion
Pressure force and fluid weight are assumed to be acting on the mass of fluid in motion.

Total pressure force acting on the left face= PΔyΔz

 P 
Total pressure force acting on the right face=  P  x yz
 x 
Net Pressure force in X-direction:

 p 
F px   xyz
 x 

 P 
PΔyΔz ΔΔ Δy P  x yz
 x 

Δz

Δx
mass x acceleration = body force + pressure force
p
 xyz a x  xxyz   xyz 
x
where x is the body force per unit mass
1 P
ax  x 
 x
u u u u
a
But x   u  v  w
t x y z

u u u u 1 P
u v w x 
t x y z =  x

v v v v 1 P
u v  w y 
t x y z =  y Euler’s Equations

w w w w 1 P
u v w = z 
t x y z  z

6
Bernoulli’s Equation
Consider a stream tube of length

Gravity force = mg

Force balancing in ‘S’ direction

dp
-mg - [ pdA  AS ] ma x
ds

dp
-mg AS ma x .........................................................(1)
ds
But m  AS
z

s
v v
ax  v Δs
t s

dp
Equation (1) becomes: pdA  AS
ds

dz dp AS
-g -  ax θ
ds ds AS

dz dp 1 v v
-g -   v PdA mg sinθ
ds  dS t s
v
For steady flow 0
t

dp
gdz   vv  0

Integrating the above equation

v 22  v12
g Z 2  Z1   
dP
 0
 2

7
GENERAL STRESS SYSTEM IN A DEFORMABLE BODY:
In order to determine the surface forces, let us consider a volume element dV= x.y.z ,whose
front lower left corner has the coordinates x, y, z as shown in figure.
 xz
 xz  dx
x

 xy
 xy  dx
x

 x
x Δz x  dx
x
 xy  xz

Δy
Δx

 x  Normal stress acting on the elementary plane which is perpendicular to X-direction

 y  Normal stress acting on the elementary plane which is perpendicular to Y-direction
 xy  Tangential stresses.
They acquire a double index: the first position indicates to which axis the surface element is
perpendicular, and the second position states in which direction the stress τ is pointing.

  x  xy  xz 
 
Stress Tensor Π =  yx  y  yz  .............................................................(1)
 
 zx  zy  z 

The Momentum equation can be written as

Du    xy  xz 
  FBX   x   
Dt  x y z 

Dv   yx  y  yz 
  FBY     
z 
...........................(2)
Dt  x y

Du    zy  z 
  FBX   zx   
Dt  x y z 
8
Symmetry Tensor:
Let τ = Tangential force which causes rotational motion of the fluid.
ς=Normal stress which causes translatory motion of the fluid.
Let us consider a 2-Dimensional element where Δz=1.


Y  yx   yx x
x


 xy  xy   xy x
x

 yx X

Consider a fixed point ‘C’, now momentum about point “C” is given
d z
I zz   xy y.1x  yx x.1y
dt
d z
Ak 2   xy   yx xy ....................................................................(3)
dt
As the element is small, in the limit k 2 tends to zero;

 xy   yx
 xz   zx and .......................................................(4)
 yz   zy
  x  xy  xz 
 
 
Stress Tensor Π = xy  y  yz  which is symmetric....................................(5)
 
 xz  yz  z 
Considering a symmetric stress tensor ,the Momentum equation can be written as

Du    xy  xz 
  FBX   x   
Dt  x y z 

Dv   yx  y  yz 
  FBY      ..........................(6)
Dt  x y z 

Du    zy  z 
  FBX   zx   
Dt  x y z 
For ideal fluids: friction= 0; Shear stress =0

9
So  x   y   z   P
 x 0 0   P 0 0 
   
Stress Tensor Π =  0  y 0  = 0  P 0 
0 
0 z   0 0  P 

Trace of stress tensor =Arithmetic mean of normal stresses.
 x  y  z  3P
= =  P
3 3
For real fluids:  x   y   z
 u v w  u
 x  P        2
 x y z  x
 u v w  v
 y  P        2 .............................(7)
 x y z  y
 u v w  w
 z  P        2
 x y z  z
2
From Stokes Theorem    
3

  x   P  2   u  v  w   2 u
3  x y z  x
2  u v w  v
 y   P        2 ..........................(8)
3  x y z  y
 u v w  w
 z  P        2
 x y z  z
STRAIN RATE:
The rate at which the fluid element is strained in a flow field is called “Strain Rate”.

Let us consider an element AB at time t=0.

At time t +δt, the element is AB .
Relative motion in X-direction = u + du –u =du A
u u u
du= x  y  z A
x y z

Similarly, u u+du
v v v
Relative motion in Y-direction= dv = x  y  z ........(9) B
x y z
w w w
Relative motion in Z-direction= dw= x  y  z B
x y z
Let us consider a 2 dimensional element ABCD,

10
After time  t , the element is A’B’C’D’.

Displacement of cover A:
u
Linear rate of strain=  x=
x
v
Linear rate of rotation=  1= (Counter clockwise)
x

Similarly displacement of D over A:

v
Linear rate of strain=  =
y
y
u
Linear rate of rotation=  2= (Clock Wise)
y

 Xy= Rotation of diagonal AB = 1 (  2-  1) = 1  v  u 

2 2  x y 
Similarly;
v
y=
y
 z= w
z
 XZ= 1  w  u 
2  z z 

 YZ= 1  w  v  the rate of Strain Tensor which is symmetric is

2  y z 
11
  x xy xz 
  xy y yz (10)
xz yz z 

 u 1  v u  1  w u 
      
 x 2  x y  2  x z 
   xz   1  u  v 
 v 1  w v  
 x
   
xy

Strain Rate Tensor= =  xy  y  yz  =  2  y x  y 2  y z  


 xz  yz  z   1  u  w  1  v w 
  
w 
 2  z x  2  z y  z 
   

(11)
u u
du = x + y + u  z (12)
x y z

du= (  x x +  xy y +  xz  z) + (  y  z-  z y ) (13)

Similarly,

dv= (  xy x +  y y +  yz  z) + (  z  x-  x  z) (14)

dw= (  xz x +  yz y +  z  z) + (  x  y-  y  x) (15)

u
from equations (12) and (13), coefficients of x =
x
u
y =
y
=  
xy z

u
z =
z
=  
xz y

u 1  v u 
Coefficient of y = =   =   - z
y xy z
2  x y 

1  v u 
z=   
2  x y 

Similarly, (16)

1  u w  1  w v 
y=    & x =   
2  z x  2  y z 

12
 , Angular Velocity.  x,  y,  z are components.

1
= curl 
2

Thus the overall motion of the particle occupying a position (x, y, z) at time t is composed of,

1. A pure translation, described by the velocity components u, v, w.

2. A rigid body rotation, described by the angular velocity components  x,  y and  z.
3. Dilatation given by linear strains  x  y and  z.
Deformation given by the three components  xy,  yz and  xz.

NAVIER STOKES EQUATIONS

Let us consider a fluid particle having a shape of parallelepiped with sides dx, dy and dz.

The mass of the particle is  .[Link].

D
Acceleration of the particle is, a=
Dt
^ ^ ^
The body force is F BX =  .[Link]. (x i + y j + z k )
Surface force = P

 x  xy  xz
Px=  
x y z
Equation (8) gives the value of  , normal stress

 v u  w u
And  xy = 2   xy =     ;  xy = 2   xz =     ;
 x y   x z 

Momentum Equation is;


Du
 FBX +   
  x  xy  xz 
   
Dt  x y z 
 

  2 u    v u     w u 
= FBX +
x  P  3  (v)  2 x   y   ( x  y )  z    x  z 
   

Du P   u 2    v u     w u 
  FBX   2    (v)    (  )       (1)
Dt x x  x 3  y  x y  z   x z 

Dv P   v u    v 2     w v 
Similarly,  = FBY     (  )  2   (v)      (2)
Dt y x  x y  y  y 3  z   x z 

13
 P  w u    w v   w 2
= FBZ    (  )       2   (v)
Dw
 (3)
Dt z x  x z  y   y z  z  z 3 
Where,

 u v w 
v =    
 x y z 

These equations are the Navier Stokes equations for the 3 Dimensional UNSTEADY,
COMPRESSIBLE, VISCOUS flow.

CASE 1: N-S Equation for UNSTEADY, INCOMPRESSIBLE AND VISCOUS FLOWS:

For Incompressible flows; (  =constant and  = constant)

 u v w 
Continuity equation for incompressible flow is v = 0,  x  y  z  =0;
 
Du P   u 2    v u     w u 
  FBX   2    (v)    (  )      
Dt x x  x 3  y  x y  z   x z 

P   2u  2u  2u    u v w 
= FBX    2  2  2       
x  x y z  x  x y z 

Du P   2u  2u  2u 
  FBX    2  2  2 
Dt x  x y z 

Du P
  FBX  +  2 u
Dt x

Similarly;
Dv P
 = FBY   2 v
Dt z

Dw P
 = FBZ   2 w
Dt z

Case 2 : N-S Equations for STEADY,INCOMPRESSIBLE and VISCOUS flow:

 u u u  P
 u v  w   FBX    2 u
 x y z  x

Similarly;

 v v v  P
 u v  w   FBY    2 v
 x y z  y

 w w w  P
 u v w   FBz    2 w
 x y z  z 14
Further If The Body Forces are zero :

* + [ ]

ρ* + [ ]

ρ* + [ ]

ENERGY EQUATION:
For an incompressible flow the unknowns are the velocity components
u,v,w and the pressure p. To obtain these four quantities we have the four equations in the
form of continuity and the three Navier-Strokes [Link] a compressible flow,
temperature(T) is an additonal unknown; the density(ρ) can be calculated using the equation
of state when p and T are known. The additional equation for T is obtained from first law of
Thermodynamics to a fluid particle. This equation is called Energy Equation.

Let us consider a fluid element of parallelepiped shape with sides

dx,dy,dz. The volume of the element ΔV is dx dy dz and its mass is Δm= ρΔV.
Let an amount of heat dQ be added to the fluid element in time dt and increase in internal
energy be ρΔV, increase in kinetic energy be ρΔV d(u2+v2+w2) and dW be the corresponding
work performed by the particle. Then applying the first law of Thermodynamics we get :

…………………..(1)

where is the substantial derivative and

ρΔV ( )

The addition of heat is assumed to be only due to conduction. According to Fourier’s law, the
heat flux, q (J/m2s) per unit area A and time is proportional to the temperature gradient
normal to the surface.

i.e.,

………………….(2)
where k (J/m s K) is the thermal conductivity of the fluid.

15
Hence the heat transferred into the volume  V * through the surface normal to the X
kT 
direction is equal to -   dydz .the amount leaving the x- direction is
 x 
 T    T  
 k x   x  k x dx  dydz.
    

Thus the total amount of heat added b conduction during the time dt to the volume V * can
be written as;

   T    T    T 
dQ  dtV *  k  
k 
  z  k z .......................3
 y  x  y  y   

To obtain the work done, let us consider the contribution from the Component , xy of the
stress.

We get the work done per unit time as;

  u  xx  
dWxx  dydz uxx   u  dx  xx  dx   V * uxx 
  x  x  x

The total work performed b the normal and shear stresses can now be written as

  
dW  V *  uxx  vxy  wxz   uxy  vyy  wyz    utzx  vtyz  wzz ......4
 x  y z

Substituting the stresses xx,yy,zz ,xy,yz,zx can be written in terms of  and velocity
components.

Substituting the above equations we get after simplification:

 De    T    T    T 
   p.V   k   k   k   
 Dt  x  x  y  y  z  z 
Where  is called dissipation function, is

 u  2  v  2  w  2   v u  2  w v  2  u w  2 2  u v w  2
  2                         
 x   x   x    x y   y z   z x  3  x y z 

TURBULENCE:

The Reynolds number of a flow gives a measure of the relative importance of inertia force and
viscous forces. In experiments on fluid systems it is observed that at values below the so called
Reynolds number Recrlt the flow is smooth and adjacent layers of fluid slide past each other
in an orderly fashion. if the applied boundary conditions do not change with time the flow is
steady. this regime is called “LAMINAR FLOW”.

16
At values of the Reynolds number and above Rexpt a complicated series of events takes place
which eventually leads to a radical change of the flow character. In the final state the flow
behaviour is random and chaotic. The motion becomes intrinsically unsteady even with
constant boundary condition. The velocity and all other flow properties vary in a random and
chaotic way. This regime is called “TURBULENT FLOW”.

Take, for example, the common tap by a domestic sink. Slowly turn the tap on and see that
water rips out of the tap. Open the tap further to increase the flow rate until a steady column
of water comes out of the tap. Notice how smooth the water column is, appearing crystal
clear like glass. Increase the flow rate further and the water column surface begins to move
slowly before the whole column becomes opaque. At this final stage the water flows in a
direction which is generally downwards, but if we look at one point in space in the water
column the fluid seems to move in a random fashion, a so called turbulent motion, which is
superimposed on the general flow. This simple experiment with the flow out of a tap
demonstrates that two main types of flow can be seen with viscous fluids; first a smooth
laminar flow, where the water moves layer over layer giving a clear column of liquid , and a
randomly fluctuating turbulent flow.

AVERAGING:
Mathematically, the instantaneous velocity component u can be described as

u = ū + uˊ where ū is the mean velocity and

uˊis the random fluctuating component.

Similarly,
V =v+vˊ w=w+ wˊ p=p+ pˊ ρ = ρ+ ρˊ T= T +Tˊ

The average is formed as the time average at a fixed point in space, thus for example

u= ∫

17
Where T= time period of sampling which is sufficiently large to include adequate number of
fluctuations
T
to 
2
1
T to T
Similarly  dt
2

The integral is to be taken over a sufficiently large time interval T so that the average is
independent of the time. The time averages of the fluctuating quantities are then zero by
definition.

 '  0;  '  0; w'  0; p '  0;  '  0; T '  0;

T T T T
to  to  to  to 

 
2 2 2 2
1 1 1 1
1)  '  T  u ' dt  T  u  u dt  T  udt  T  udt  u  u  0
T T T T
t o t t o t o o
2 2 2 2

2) u 2  u  u'  2
 u u  2uu 'u '
2

Averaging of u  u  u u  2uu '  u '  u u  u '

2 2 2 2

3) Averaging of uv = uv  u  u ' v  v '   

 uv  uv'  u ' v  u ' v'

 uv  u'v'

4) Averaging of uw = u w  u'w'
5) Averaging of vw = v w  v'w'
u u

6)
 

7)  ud   udn

18
MASS WEIGHTED AVERAGING:
~ u ~ v ~ w ~ h ~ T ~ H
u  v  w  h p T 
 H 

u  u~  u *
v  v~  v *
~  w*
ww
u*  v*  w*  0
For incompressible flows(p=constant); then u*  v*  w*  0

~ ou  'u'  'u'

1) u *  u  u  u  u u  
  

REYNOLDS FORM OF CONTINUITY EQUATION:

 
 uj   0
t xj


   '

 
 
   ' uj  u ' j  0 
t xj

   ' 
t
 
t xj
 uj 

xj
 
u ' j 

xj
 'u j   

xj
 
 ' u ' j   0
Taking time average;

 ' 
t
 
t xj
 uj 

xj
 
u ' j 

xj
 
 ' uj 

xj
 
 'u' j  0  



t xj
 uj   

xj
 'u' j  0  



t xj

 uj   ' u ' j  0 

19
FAVRE Average form of Continuity equation:

( )

[ ̃ ]

( ̃) ̃

Taking time average ,

( ̃) ̃

The fourth and last terms can be rewritten as:

( ) ( ) * + ( )

( ̃)

Reynolds form of Momentum Equation:

( ) ( )

Taking time average,

( ) ( )

where

[( ) ]

20
  l 2 lv ' lw' 
 2 
 ' v' w 
Reynolds stress matrix (or) Apparent stress matrix = π =  lv v ' '

2 
 lw v' w ' w ' 
'
 

The fluxes of momentum arise as we have carries the averaging process over non- linear
equations. These are additional unknowns in the equations. To solve the equations for these
quantities, they would involve co relations among fluctuating quantities for higher order and
so on. Thus, the number of unknowns always exceeds the number of equations for reynold’s
averaged quantity. This is called “CLOSURE PROBLEM”.
To obtain a closed set of equations, the higher order correlations need to be expressed
in terms of lower order co relations and their derivatives, such expressions are called
“TURBULANCE MODELS”.

A Turbulent model should:

1) Predict accurately a wide range of flows.

2) Be mathematically simple and involve small number of assumptions and model
constants.
3) Be computationally stable and economical in computational time.

TURBULENT MODELS:

A turbulence model is a computational procedure to close the system of mean flow

equations so that a more or less wide variety of flow problems can be calculated. For
turbulence model to be useful in a general purpose CFD code it must have wide applicability,
be accurate, simple and economical to run.

The common turbulence models are classified as:

1) CLASSICAL MODELS: Based on time-averaged) Reynolds equations

1) Zero equation model- Prandtl mixing length model
2) Two equation model – k – ε model
3) Reynolds stress equation model
4) Algebraic stress model
2) LARGE EDDY SIMULATION : Based on space filtered equations

PRANDTL MIXING LENGTH MODEL

On dimensional grounds we assume that the kinematic turbulent viscosity 

21
Which has dimensions m2/s, can be expressed as a product of a turbulent velocity V1(m/s) and
a length scale lm(m).

where C is a dimensionless constant of proportionality.

And | | ;

| |

| || |

ADVANTGES:
 Easy to implement and cheap in terms of computing resources.
 Good predictions for thin shear layers : jets, mixing layers, wakes and boundary layers.
 Well established.

DISADVANTAGES:
 Completely incapable of describing flows with separation and recirculation.
 Only calculates mean flow properties and turbulent shear stress.

K-ε MODEL:
The instantaneous kinetic energy k(t) of a turbulent flow is the sum of mean
kinetic energy K= and the turbulent kinetic energy ( )
; k(t)=K + k .

Momentum equation,

( ) ( ) { [ ] }

Turbulent viscosity or Eddy viscosity =

22
 The value of k and ε come directly from the differential transport equations for turbulent
kinetic energy and dissipation rate as:

 ( k )  (  u j k )   k 
  k   PK  
t x j x j  x j 

 (  )  (  u j k )     
     C 1 Pk  C 2  
t x j x j  x j  k

t
Where diffusion coefficients k =  
k

t
   


The values of constants in the model are:

Cµ = 0.09, Ct1 = 1.44, Ct2 = 1.92,  k =1,   =1.3, Prt = 0.9

Production rate of turbulent kinetic energy Pk is given by

 U U j  U i 2  U i  U k
Pk   t  i     k  t 
 x j xi  x j 3  xi  xk

ADVANTAGES:

 Simplest turbulence model for which only initial and/ or boundary conditions need
to be supplied.

 Excellent performance for many industrially relevant flows.

 Well established; most widely validated turbulence model.

DISADVANTAGES:

 More expensive to implement than mixing length model (two extra PDE’s).

 Poor performance in a variety of important cases such as:

i. Some confined flows

ii. Flows with large extra Strains (e.g. curved boundary layers, swirling flows).

iii. Rotaint flows.

iv. Fully developed flows in non-circular ducts.

23
 2. CLASSIFICATION OF PARTIAL DIFFERENTIAL EQUATIONS
1) According to the order:

The order of a partial differential equation is determined by the highest order partial
derivative present in that equation.
 
G 0 First order
x y

 2 
 0 Second order
x 2
y

2
  3   2  
 3     0 Third order
 x  xy y
Where  is a dependent variable, x and y are independent variables and G is a constant.

2) According to Linearity:

Partial differential equations are broadly classified into linear equations and non-linear
equations.

If the coefficients are constants or functions of the independent variables then they are
called “Linear Equations”.
 
Example: G 0
x y
 2  2 
a 2 +b +c +d=0
x xy y
If the coefficients are constants or functions of the dependent variables and/or any of its
derivatives of either lower or same order, then they are called “Non-Linear Equations”.
 2 
Example: -  =0
x 2
y
ELLIPTIC, PARABOLIC AND HYPERBOLIC EQUATIONS:
Linear second order partial differential equations in two independent variables are further
classified into three canonical forms.
1)Elliptic
2)Parabolic and
3)Hyperbolic Equations.
The general form of this class of equation is
 2  2  2  
a + b + c +d e  f  g = 0
x 2
xy y 2
x y
Where coefficients are either constants or functions of the independent variables.

24
The three canonical forms are determined by the following criteria:

b2 - 4ac < 0 Elliptical

2
b - 4ac = 0 Parabolic
b2 - 4ac > 0 Hyperbolic
a)ELLIPTIC PDE:
Typical examples are:
 2  2 
+ 2 =0 (Laplace’s Equation)
x 2 y
 2  2 
And + 2 = g (x,y) ( Poisson’s Equation)
x 2 y
In both the cases; a=1, b=0 and c=1.
b2- 4ac = -4, which is <0
In Elliptic problems, the function Φ(x,y) must satisfy, both the differential equation over
a closed domain and the boundary conditions on the closed boundary of the domain.

Boundary conditions prescribed on

entire boundary – usually as a value of

Φ(x,y) or the normal derivative
n
b)PARABOLIC PDE:

Typical example is:

  2
 2 (Heat conduction or diffusion equation)
t x
Where α is a positive and real constant
In this case; a = α, b=0 and c = 0.
b2 - 4ac = 0
In parabolic problems, the solution advances outward indefinitely from known initial
values, always satisfying the known boundary conditions as the solution progresses. This
is also called “Marching type of problems”.
Y

Solution marches outward with

Prescribed boundary time from initial condition Prescribed boundary
condition condition
Solution domain for marching
problems.
Initial condition is specified here

25
C) HYPERBOLIC PDE: Typical example is

 2 2  
2
  (Wave equation)
t 2 x 2

Where  2 is a positive and real constant

In this case; a=  2 ,b=0 and c= -1.
b 2 -4ac = 4  ,which is >0
2

The solution domain of hyperbolic PDE has the same open-ended nature as in parabolic PDE.
However, two initial solutions are required to start the solution of hyperbolic equations in
contrast with parabolic equations where only one initial condition is required.
t

Prescribed specified propagating Prescribed

Boundary outward with time for Boundary
Condition specified initial condition Condition

Initial condition is specified

INITIAL AND BOUNDARY CONDITIONS

The initial and boundary conditions must be specified to obtain unique numerical solution to
partial differential equations.

Let us consider one dimensional Transient heat conduction equation;

T  2T
 2
t x

The mathematical representation of the a physical problem is, the temperature within a large
solid slab having finite thickness changes in the x-direction(thickness direction)as a function of
time till steady state is reached.

There are three categories of boundary conditions;

1) Dirichlet Boundary conditions (First Kind)

2) Neumann conditions (Second Kind)
3) Robbins conditions (Third Kind)

26
Dirichlet Boundary conditions (First Kind)

The values of the dependent variables are specified at the boundaries.

Boundary conditions are;
B.C.1 T=f(t) or T1 at x=0 t>0 T=T2
B.C.2 T=T2 at x=L Dirichlet
Initial conditions are; T=f(t) condition
T=f(x) at t=0 0≤x≤L or x
T=T0 T=T1
Neumann conditions (Second Kind) 0 L
The derivative of the dependent variable is given as a constant or as a function of the
independent variable on one boundary.
For example;
T
0 at x=L and t ≥0
x
This condition specifies that the temperature gradient at the right boundary is zero
(insulating condition).
T
0
x
Perfect
T=f(t) Insulation
Or
T=T1 x
0 L

Cauchy conditions:
A problem which combines both Dirichlet and Neumann conditions is called Cauchy
conditions.
Robbins conditions (Third Kind):
The derivative of the dependent variable is given as a function of the dependent variable on
the boundary.
For heat conduction problem, this may correspond to the case of cooling of a large steel slab
of finite thickness ‘L’ by water or oil, the heat transfer coefficient, ‘h’ being finite.

Fluid T>Ta
T
Fluid h, Ta -k  h(T  Ta )
x
T
k  h(T  Ta ) x
x
0 L
27
 3. FINITE DIFFERENCE METHODS
According to Taylor Series Expansion

x  u  x 2   2 u  x 3   3 u 
u i 1  ui         
2!  x 2  i 3!  x 3  i
______________ (1)
1!  x  i

x  u  x 2   2u  x3   2u 
ui 1  ui         
1!  x i 2!  x 2 i 3!  x3 i ______________ (2)

From forward difference eq(1)

 u  ui 1  ui
   O(x)
 x  x
From backward eq(2)

 u  ui  ui 1
   O(x)
 x  x

Central difference eq(1)-eq(2)

ui 1  ui 1  u  x 2   3u 
      ..............
2x  x  i 3!  x 3  i
 u  ui 1  ui 1
   O(x 2 )
 x  2x

Truncation error for forward and backward difference is of first order where as for central
difference it is of second order.

Adding first two Taylor Series Expansions (1)+(2)

ui 1  2ui  ui 1   2 u  x 2   4 u 
  2    4   ..........................
x 2  x i 12  x  i
  2 u  ui 1  2ui  ui 1
 2    O(x 2 )
 x  x 2

28
Euler’s Explicit method

u   2u 
   2 
t  x 
ui
n 1
 ui
n
 u n i 1  2ui n  u n i 1 
   

t  x 2

Euler’s Implicit method

u   2u 
   2 
t  xo 
n 1
ui  ui
n
 u n 1i 1  2ui n 1  u n 1i 1 
   

t  x 2

Crank Nicolson Method

u   2u 
   2 
t  xo 
ui
n 1
 ui
n
  u n 1i 1  2ui n 1  u n 1i 1 u n i 1  2ui n  u n i 1 
   
t 2 x 2 x 2 
COMPARISION OF IMPLICIT/EXPLICIT METHODS

 Explicit method involves pointwise updating & requires no matrix inversion. Implicit
scheme needs matrix inversion.
 Computational time for time set is more for more for implicit method than the explicit.

For stability considerations explicit scheme may require very small time steps and hence
several thousand steps to obtain steady solution. Large time steps can be used in implicit
scheme.

29
General Methods to get Higher order Accuracy
Q1: Derive backward difference equation for second derivative with second order accuracy.
Answer:
Second derivative r=2
Second order m=2
m+r=4 .Therefore four terms ui, ui-1, ui-2, ui-3 are considered for the equation.
x  u  x 2   2u  x 3   3u 
u i 1 ui           .................
1!  x i 2!  x 2 i 3!  x 3 i
(2x)  u  (2x) 2   2u  (2x)3   3u 
u i  2  ui           .................
1!  x i 2!  x 2 i 3!  x 3 i
(3x)  u  (3x) 2   2u  (3x)3   3u 
u i  3  ui           .................
1!  x i 2!  x 2 i 3!  x 3 i
ui 1  ui  2  ui  3 
x  u  x 2   2u  x 3   3u 
(     )ui  (  2  3 )    (  4  9 )  2   (  8  27 )    .............
1!  x i 2!  x i 3!  x 3 i

 u    2u    3u 
setting the coefficients of   ,  2  and  3  to 0,1,0 respectively
 x i  x i  x i
(  2  3 )x =0

x 2
(  4  9 ) =1
2
x 3
(  8  27 ) =0
6
Solving the above equations
α=-5/Δx2 β=4/Δx2 γ=-1/Δx2
Substituting the above values in the equation ui 1  ui 2  ui 3
We get

 2u  5ui 1  4ui  2  ui  3  2ui

  O(x 2 )
x 2
x 2

Assignment:
Q2: Derive central difference equation for third derivative with second order accuracy.
Q3: Derive forward difference equation for third derivative with first order accuracy.
30
Solution of 2nd Order Differential Equation On One Dimensional Domain
y   4 y   2 y  x

BCs are at x=0,y=0 and 0 1

at x=1,y=1 0.25 0.5 0.75 4

1 2 3

2 y y
 4 +2y=x
x 2
x

The difference equation is

y i 1  2 y i  y i 1  y  y i 1 
 4 i 1   2 yi  x
x 2  2x  

y 2  2 y1  y 0  y  y0 
 4 2   2 y1  x
0.252  20.25 
At i=1;

16 y2 -32 y1 +8 y2 +2 y1 =0.25
-30 y1 +24 y2 =0.25 ………………(1)

y 3  2 y 2  y1  y  y1 
 4 3   2 y 2  0.5
At i=2;
0.252  20.25 

16 y 3 -32 y2 +16 y1 +8 y 3 -8 y1 +2 y2 =0.5

8 y1 +-30 y2 +24 y 3 =0.5 ………………(2)

y 4  2 y3  y 2  y  y2 
 4 4   2 y 3  0.75
 20.25 
At i=3;
0.252
16 y4 -32 y 3 +16 y2 +8 y4 -8 y2 +2 y 3 =0.75
8 y2 -30 y 3 +24 y4 =0.75
8 y2 -30 y 3 =-23.25 ………………(3)

Solving eq (1),(2) and (3); we get

y1 = 0.83042
y 2 = 1.04844
y 3 = 1.0545 31
 [Link]-NEUMANN STABILITY ANALYSIS
N r 1 x
1k
Let U(0,x)=  Ak e L

k 0

Where; k=Wave Number A=Amplitude of k th harmonic

rt 1  k i k
U
n
i
= e e Where  k 
L

r = Complex constant

t    1  k i 1  k i
 =  Ui
n rt rt n
U i = e e = e  e e
n t    1  k i
U i = e e
=
n n
U i U i

n 1
=  Ui
n
U i

n
n 1 U i
U i =


U
n
i 1 =e
1  k i
U  n
i

U
n
i 1 =e
 1  k i
U  n
i

U
n 1
i 1 = e
1  k i
U  n
i

U
n 1
i 1 = e
 1  k i
U  n
i

STABILITY ANALYSIS OF EXPLICIT METHOD

U  2U

t x 2

U
n 1
 Ui
n
Un  2UnUn 
i
   i 1 i i 1

t  x 2


 t 
Let r=
x t

U
n 1
i
 Ui  r Ui 1  2Ui  Ui 1
n
 n n n
 32
 Ui Ui  r  U 
1  k i  1 k x
2Ui  e
n n n n n
e U
i i

 1  r  e 1 k x
2e
 1 k x


= r 2Cos k x  2

= 2r Cos k x   1

   k x 
=2r  2rSin 2  
  2 

2   k x 
=4r Sin  
 2 

2   k x 
=1-4r Sin  
 2 

 t  2   k x 
For a solution to be stable: 1  4 2  Sin   1
 x   2 

For this condition there are low possible situation which must hold simultaneously are;

 t  2   k x 
1) 1-4  2  Sin   ≤1
 x   2 
 t  2   k x 
Thus 4  2  Sin   ≥0
 x   2 
 t 
Since  2  is always poisitive; this condition always holds good.
 x 
 t  2   x 
2) 1-4  2  Sin  k  ≥-1
 x   2 
 t  2   k x 
Thus 4  2 
Sin   -1≤1
 x   2 
 t 1
This condition holds good when 
x 2 2

Explicit Method is stable with very small time increments.

33
STABILITY ANALYSIS OF IMPLICIT METHOD:
U  2U

t x 2
n 1
 Ui
n
Un 1 2Un 1  Un 1 
=U i
   i 1 i i 1

t  x 2


 t 
Let r=
x 2
n 1 n 1 n 1 n 1
U i
 Ui =r[
n
Ui1  2Ui  Ui1 ]
   n
i
n
i
=r[  e
1 k x
i
n
 2  i  e  1  xi ]
n

k
n

  1 =[  e 1 k x
 2 +  e
 1 k x
]

=r[  [2Cos(  k x )-2  ]

=2r  [Cos(  k x )-1]

=2r  [-2Sin 2 (
 k
x
)]
2

 x
=  4r Sin ( )
2 k

 x
 = 1  4r Sin ( 2 k

2
)

 [1  4r Sin 2 (
 k
x
)]=1
2

1

 2  k x 
1  4 rSin ( )
2 
Unconditionally stable for all values of r

34
 STABILITY ANALYSIS OF CRANK – NICOLSON METHOD
U  2U

t x 2

 Ui 1  2Ut  Ui 1  2Ui  Ui 1 

n 1 n 1 n 1 n 1
 Ui
n n n n
U i
  U i 1

t 2  x 2 x 2 

 t 
Let r=
x 2

U
n 1
i
 Ui 
n r
2
U n 1
i 1
n 1 n 1
 2Ui  Ui 1  Ui 1  2Ui  Ui 1
n n n

U n 1
i
 Ui
n

r


2 e
1 k x
U  2U   e
n
i
n
i
 1 k x
U e
n
i
1 k x
U  2U  e
n
i
n
i
 1 k n
U
i

 1 =
r

2 e
 1 k x
 2   e
 1 k x
e
1 k x
2e
 1 k x

  1 = r  2Cos k x   2  2Cos( k x   2 
2

r 2   k x  2   k x 
  1= 2   4rSin    4 Sin  
  2   2 

  k x  2rSin 2   k x 
  2rSin   =1- 2

 2   2 
  x 
1  2rSin 2  k 
 2 

  x 
1  2rSin 2  k 
 2 

Unconditionally stable for all values of r.

35
 [Link] METHODS
GAUSSIAN ELIMINATION METHOD:
The algebraic equations are to be written in the form of AX=B. Row equations are performed
to obtain lower triangular elements as zero to obtain the solution.
Example :
Let us consider a linear system of n equations with n unknowns x1,x2,x3…..xn
8x2+2x3=-7
3x1+5x2+2x3=8
3x1+x2+4x3=13

[ ] [ ]=[ ]

R2-R3

[ ] [ ]=[ ]

R3-2R2

[ ] [ ]=[ ]

6x3=3 ; x3=0.5
4x2-2x3=-5 ; x2=-1
3x1+x2+4x3=13 ; x1=4
x1=4; x2=-1 ; x3=0.5

POINT JACOBI METHOD:

An initial guess for the unknown elements are assumed for starting the iterations. The solution
obtained for an iteration will be the guess for next iteration. Iterations are to be performed till
the solution converges.
Example:
5x1+x2+2x3=19
x1+4x2-2x3=-2
2x1+3x2+9x3=39
Initial values (or) guess values;
X1=5; x2=5; and x3=5
5x1+ x2+ 2x3=19 x1= - x2 - x3+ …………………..(1)
x1+ 4x2- 2x3 = -2 x2=- x1+ x3- …………………......(2)
2x1+ 3x2 + 8x3 = 39 x3= - x1- x2+ …………………....(3)
36
At (5,5,5);

x1= - (5) - (5)+ = 0.8

x2=- (5)+ (5)- = 0.75
x3= - (5)- + = 1.75

At(0.8,0.75,1.75);

x1= - (0.75) - (1.75)+ = 2.95

x2=- (0.8)+ (1.75)- = 0.175
x3= - (0.8)- + = 4.39375

At(2.95,0.175,4.39375)

x1= - (0.175) - (4.39375)+ = 2.0075

x2=- (2.95)+ (4.39375)- = 0.959375
x3= - (2.95)- + = 4.071875

At(2.0075,0.959375,4.071875)

x1= - (0.959375) - (4.071875)+ = 1.979375

x2=- (2.0075)+ (4.071875)- = 1.0340625
x3= - (2.0075)- + = 4.013359

At(1.979375,1.0340625,4.013359)

x1= - (1.0340625) - (4.013359)+ = 1.9878439

x2=- (1.979375)+ (4.013359)- = 1.0118357
x3= - (1.979375)- + = 3.99238

Iteration No X1 X2 X3
1 5 5 5
2 0.8 0.75 1.75
3 2.95 0.175 4.39375
4 2.0075 0.959375 4.071875
5 1.979375 1.0340625 4.013359
6 1.9878439 1.0118357 3.99238

37
GAUSS-SEIDAL ITERATION METHOD :

An initial guess for the unknown elements are assumed for starting the iteration.
Unlike Jacobi method, the latest value obtained is taken for the iteration calculations.
Iterations are to be performed till the solution converges.
Example:
5x1+x2+2x3=19
x1+4x2-2x3=-2
2x1+3x2+8x3=39
Initial values or guess values; x1=5 ; x2 = 5 ; x3=5

5x1+ x2+ 2x3=19 x1= - x2 - x3+ …………………..(1)

x1+ 4x2- 2x3 = -2 x2=- x1+ x3- …………………......(2)

2x1+ 3x2 + 8x3 = 39 x3= - x1- x2+ …………………....(3)

At (5,5,5);
x1= - (5) - (5)+ = 0.8
x2=- (0.8)+ (5)- = 1.8
x3= - (0.8)- + =4
At (0.8,1.8,4);
x1= - (1.8) - (4)+ = 1.84
x2=- (1.84)+ (4)- = 1.04
x3= - (1.84)- + = 4.025

At (1.84, 1.04, 4.025);

x1= - (1.04) - (4.025)+ = 1.982
x2=- (1.982)+ (4.025)- = 1.017
x3= - (1.982)- + = 3.998125

At (1.982, 1.017, 3.998125);

x1= - (1.017) - (3.998125)+ = 1.99735
x2=- (1.99735)+ (3.998125)- = 0.99975
x3= - (1.99735)- + = 4.000765

38
Iteration No X1 X2 X3
1 5 5 5
2 0.8 1.8 1.8
3 1.84 1.04 1.04
4 1.982 1.017 1.017
5 1.99735 0.999725 0.999725

39
 [Link] OF FINITE DIFFERENCE METHOD :
Finite Difference Method of 1-D Heat Conduction

B.C. are At x=0, T1=30

At x=L, TL=100

1 2 3 4 5 6
T1  30 0 T6  100 0

Finite Difference Equation

At i=1; T1=30
At i=2;

At i=3;

At i=4;

At i=5;
At i=6; T6=100

[ ]x[ ]=

[ ]
The above matrix is a tri-diagonal matrix and can be solved using Thomas Algorithm.

Finite Difference Method 1-D Heat Conduction with Flux

Type BC

B.C. are At x=0, T1=30

At x=L,

1 2 3 4 5 6

T
Finite Difference equation T1  30 0 0
x

40
 Q (x)
2

T i 1
 2T i  T i 1 
k

At i=1; T 1  30
Q (x) Q (x)
2 2

At i=2; T 3  2T 2  T 1  T  2T
3 2
  30
k k
Q (x)
2

At i=3; T 4  2T 3  T 2  5 6 7
k
Q (x)
2

At i=4; T 5  2T 4  T 3 
k
Q (x)
2

At i=5; T 6  2T 5  T 4 
k
Q (x)
2

At i=6; T 7  2T 6  T 5 
k

Image point method:

Imagine 7th node

 T   T  T
  2(x) 7 5
 0 that implies T T5
 x 
7

 Q (x)
2

  30
 2 1 0 0 0  T   
   2
 
k

2 Q (x)
2
 1 1 0 0  T 3  
 0 
1 2 1 0   T    k 
   4
  
Q (x)
2

 0 2 1  T  
5
0 1 
 0  
 0 0 2  2  T 6
  k 
Q (x)
2
 
  
 k 
 Q (x) 
2

  
 k 

The above tridiagonal matrix can be solved using Thoma’s algorithm.

FINITE DIFFERENCE METHOD OF TWO-DIMENSIONAL DOMAIN

SOLUTION TO ELLIPTIC EQUATIONS

LAPLACE EQUATION 4
(1,3) (2,3) (3,3)
 u  u 0
2 2

3
x y
2 2
(1,2) (2,2) (3,2)
U  2U  U
i 1, j i, j i 1, j

U i , j 1
 2U i , j  U 1, j 1
0 2 (1,1) (2,1) (3,1)
x 2
x 2

1
41 U=0
 For uniform grid; x  y U=0

U i 1, j
 4U i , j  U i 1, j  U i , j 1  U i , j 1  0

At (1,1) U  4U  U  U  U
2 ,1 1,1 0 ,1 1, 2 1, 0
0
At (2,1) U  4U  U  U  U
3,1 2 ,1 1,1 2, 2 2, 0
0

At (3,1) U 4 ,1
 4U 3,1  U 2,1  U 3, 2  U 3, 0  0

At (1,2) U 2, 2
 4U 1, 2  U 0, 2  U 1,3  U 1,1  0

At (2,2) U 3, 2
 4U 2, 2  U 1, 2  U 2,3  U 2,1  0

At (3,2) U 4, 2
 4U 3, 2  U 2, 2  U 3,3  U 3,1  0

At (1,3) U 2,3
 4U 1,3  U 0,3  U 1, 4  U 1, 2  0

At (2,3) U 3, 3
 4U 2,3  U 1,3  U 2, 4  U 2, 2  0

At (3,3) U 4,3
 4U 3,3  U 2,3  U 3, 4  U 3, 2  0

 4 1 0 1 0 0 0 0 0  U 11   0 
 1 4 1     
 0 1 0 0 0 0  U 21  0 
0 1 4 0 0 1 0 0 0  U 31   1
     
 1 0 0  4 1 0 1 0 0  U 12   0 
0  
1 0 1 4 1 0 1 0  U 22   0 
      1
 0 0 1 0 1  4 0 0 1  U 32   
0 0 0 1 0 0 4 1 0  U 13    1
     
 0 0 0 0 1 0 1  4 1  U 23   1
0 1  4 9 X 9 U 33  2
 0 0 0 0 1 0

42
Then

 
U i 1, j  2U i , j  U i 1, j   2 U i , j 1  2U i , j  U i , j 1  0
U i 1, j  U i 1, j   2U i , j 1   2U i , j 1  2(1   2 )U i , j  0
U i 1, j  U i 1, j   2U i , j 1   2U i , j 1  U i , j  0

Where   2(1   2 )
At (1,1); U 2,1  U 0,1   2U 1, 2   2U 1,0  U 1,1  0
At(2,1); U 3,1  U 1,1   2U 2, 2   2U 2,0  U 2,1  0 (1,3) (2,3) (3,3)
At(3,1); U 4,1  U 2,1   2U 3, 2   2U 3,0  U 3,1  0
(1,2) (2,2) (3,2)
At(1,2); U 2, 2  U 0, 2   2U 1,3   2U 1,1  U 1, 2  0 U  U1
At(2,2); U 3, 2  U 1, 2   2U 2,3   2U 2,1  U 2, 2  0 (1,1) (2,1) (3,1)
At(3,2); U 4, 2  U 2, 2   2U 3,3   2U 3,1  U 3, 2  0
At(1,3); U 2,3  U 0,3   2U 1, 4   2U 1, 2  U 1,3  0
At(2,3); U 3,3  U 1,3   2U 2, 4   2U 2, 2  U 2,3  0
At(3,3); U 4,3  U 2,3   2U 3, 4   2U 3, 2  U 3,3  0 U  U1

 1 0 2 0 0 0 0 0  U 11    U 01   2U10 
     
1  1 0 2 0 0 0 0  U 21    2
U 20 
0 1  0 0 2 0 0 0  U 31    U10  U 41 
2

 2     
 0 0  1 0 2 0 0  U 12    U 
 
02
0 2 0 1  1 0 2 0  
U 22 
 0 
   
0 0 2 0 1  0 0  2  U 32    U 42 
0 0 0 2 0 0  1 0  U 13   U 03   U14 
 2
     
0 0 0 0 2 0 1  1  U 23    U 24  2

   
0 0 0 0 0 2 0 1   9 X 9 U 33  U 43   2U 34 

43
CONVECTIVE DIFFUSION EQUATION

dT d 2T
The governing equation u  2
dx dx

CENTRAL DIFFERENCING

The difference equation is

u {(Ti+1 + Ti-1)/2Δx} = α { (Ti+1 + Ti-1 – 2Ti ) / Δx2}

2Ti–(1-Pee/2)Ti+1 - (1+Pee/2)Ti-1 = 0

Where Pee = u Δx/α

Cell Pe< 2for spatial stability, when central difference isused

UPWIND DIFFERENCING
u{(Ti - Ti-1)/Δx} = α{(Ti+1+ Ti-1– 2Ti)/Δx2}

FOR U>0

2(2+ Pee)Ti - Ti+1-(1+ Pee) Ti-1 = 0

FOR U<0
2(2+ |Pee|)Ti - (1+ |Pee|) Ti+1 -Ti-1 = 0

The method is unconditionally stable,

Upwinding can be done with higher order accuracy, For node I , we can consider the nodes i-2,
i-1 and i to get second order accurate expression for convective term. Even nodes i-2, i-1, i
and i+I can be taken for third order accuracy.

44
ARTIFICIAL DIFFUSION

Artificial diffusion can be used in flow direction for high speed flows to avoid numerical
oscillations ; need not be used in cross flow direction it can be used to smoothen the solution
at shocks & high gradient regions

u (dT/dx) = α (d2T /dx2) + α0u (d2T /dx2)

The last term on the right is the artificial diffusion term

(uΔx/2α)( Ti+1 -Ti-1) = (1 + α0/α )(Ti+1 + Ti-1 – 2Ti )

By setting (α0/α) = Pee/2, one can get the upwind form from central difference from for
artificial diffusion , 2nd order , 4th order or 6th order expressions etc. Can be used

HIGHER ORDER ARTIFICIAL DIFFUSION

u (dT/dx) = α (d2T /dx2) + α0u (d2T /dx2) + α0iv (d4T/dx4) + α0vi (d6T/dx6)

Where d2T /dx2 = (Ti+1 + Ti-1 – 2Ti ) / Δx2

d4T/dx4= ( Ti+2 – 4Ti+1 + 6Ti – 4Ti-1 + Ti-2)/ Δx4

INCOMPRESSIBLE FLOW FORMULATIONS: Incompressible flow formulations can be solved by

primitive variables – Velocity components and pressure or stream function/ vorticity or stream
function velocity/ vorticity equations.

INCOMPRESSIBLE FLOW EQUATIONS

Two dimensional steady state flow equations are
(∂u/∂x) + (∂v/∂y) = 0 Continuity Equation
u (∂u/∂x) + v (∂v/∂y) = - (1/ρ)( ∂ρ/∂x) + (µ/ρ){(∂2u/∂x2) + (∂2u/∂y2)}x- Momentum Equation

u (∂v/∂x) + v (∂v/∂y) = - (1/ρ)( ∂ρ/∂y) + (µ/ρ){( ∂2v/∂x2) + (∂2v/∂y2)}Y - Momentum Equation

VORTICITY STREAM FUNCTION FORMULATION

Cross differentiating the x - and y – momentum equations, and subtracting one from the
other, we get
Vorticity Transport Equation u (∂σ/ ∂x) + v (∂σ / ∂y) = µ/ρ {(∂2σ/ ∂x2 )+( ∂2σ / ∂y2)}

45
 v u
Where vorticity is defined as     V  (  )
x y
For 2 – d incompressible flow , we can define stream function y as

u = ∂ψ/∂y ; v = - ∂ψ/∂x

In terms of y the equations for vorticity definition and vorticity transport become

  2  2 
 = -  2  2 
 x y 
       2  2 
    
x y y x   x 2 y 2 

LID DRIVEN CAVITY FLOW

A cavity on the flow surface is shown in the following figure.
U  u1 , v  0

Descretisation of the cavity into 16 elements as shown below.

I,j+1
U = u1, v = 0

i-1,j I,j i+1,j

u=0 u=0 i,j-1

v=0 v=0

u=0,v=0

46
FOR INTERIOR POINTS THE DIFFERENCE EQUATION IS
 i 1, j   i 1, j  2 i , j  i , j 1   i , j 1  2 i , j
   i , j
x 2 y 2
 i 1, j   i 1, j   2  i , j 1   i , j 1   21   2  i , j  x 2 i , j
 i , j 1   i , j 1  i , j   i 1, j   i 1, j   i 1, j  i , j   i , j 1     i 1, j   i 1, j  2 i , j  i , j 1   i , j 1  2 i , j 
    
2y x 2x y   x 2 y 2 
Where β = Δx/Δy

APPLICATION OF BOUNDARY CONDITIONS

For stream function:

Any surface having normal velocity = 0must be a y = constant. Surface. Hence, all the
boundaries of the cavity have same y value, which can be conveniently set as zero

47
7.3 STAGGERD MESH PROCEDURE

Pressure nodes are taken as the main nodes.

X-velocity (u) nodes are shifted by x/2 with reference to pressure nodes and y-velocity
(v) nodes

are shifted by y’/2 with reference to pressure nodes.

Such a staggered mesh avoids odd-even decoupling (checquer- board configuration)

between

velocities & pressures

V-VELOCITY

U-VELOCITY

PRESSURE

EQUATIONS OF COMPRESSIBLE FLOW:

Continuity equation

x-momentum , -

y-momentum , -

energy equation ( ) ( )

COMPRESSIBLE FLOW SIMULATION:

For compressible flow the mass, momentum & energy balance equations are integrated
with respect to time to obtain and e values at new time level.

From these variables, the velocity, temperature are derived.

New pressure is obtained using the state equation.

50
[Link]:

Convergence is the property of a numerical method to produce a solution which

approaches the exact solution as the grid spacing, control volume size or element size is
reduced to zero.

[Link]:

Consistent numerical schemes produce system of algebraic equations which can be

demonstrated to be equivalent to the original governing equation as the grid spacing tends to
zero.

[Link]:

Stability is associated with damping of error as the numerical methods proceeds.

If a given numerical error is amplified from one time step to the next, then the calculation will
become unstable.

If the error does not grow; if it decreases from one step to another, then the calculation has a
stable behavior.

[Link]:

Numerical solution of the equation is influenced by three sources of errors:

I. Discretization Error.
II. Truncation Error and
III. Round-Off Error.

I. DISCRETIZATION ERROR:

This is the difference between the exact analytical solution of the partial differential equation
and the exact solution of the corresponding finite difference equation.

The discretization error is simply the truncation error for the finite difference equation plus
any errors introduced by the numerical treatment of the boundary conditions.

II. TRUNCATION ERROR:

Truncation error is due to replacement of an exact mathematical expression by a numerical

approximation.

This is the difference between an exact expression and the corresponding truncation form
used in numerical solution.

51
III. ROUND-OFF ERROR:

This is the numerical error introduced for a repetitive number of calculations in which the
computer is constantly rounding the numbers to some decimal points.

If

A=Analytical solution of the partial difference equation.

D=Exact solution of the finite difference equation.

N=Numerical solution from a real computer with finite accuracy.

Then Discretization Error = A-D . . . . . . . . . (1)

Round-Off Error = =N-D . . . . . . . . . (2)

Equation (2) can be written as; N=D+

The numerical solution N must satisfy the finite difference equation.

Let us consider 1D heat conduction equation . . . . . . . . . (3)

* + . . . . . . . . . (4)
From the definition; D is the exact solution of the finite difference equation; hence it exactly
satisfies the differential equation.

D * + . . . . . . . . (5)

Subtracting eq(5) from eq(4);

* +. . . . . . (6)

From eq(6); we see that the error also satisfies the difference equation.

If the error , grow larger during the progression of the solution from step n to (N+1) then the
solution is unstable.

For a solution to be stable; the mandatory condition is;| |

52