Chapter 6: Solution Algorithms for

Pressure-Velocity Coupling in Steady

Flows

Ibrahim Sezai
Department of Mechanical Engineering
Eastern Mediterranean University

Fall 2010-2011

Introduction
The convection of a scalar variable  depends on the
magnitude and direction of the local velocity field.
How to find flow field?
Momentum equations can be derived from the
general transport equation (2.39)
 (  )
 div(  u)  div(  grad  )  S (6.1)
t
by replacing the variable  by u, v and w.
Let us consider the equations governing a two-
dimensional, steady flow:
ME555 : Computational Fluid Dynamics 2 I. Sezai – Eastern Mediterranean University

1
 Introduction
X-momentum equation
    u    u 
(  uu )  (  vu )          su (6.2)
x y x  x  y  y 
Y-momentum equation
    v    v 
(  uv)  (  vv)          sv (6.3)
x y x  x  y  y 
Continuity equation
 
(  u )  (  v)  0 (6.4)
x y
su  p / x for x-momentum equation
sv  p / y for y-momentum equation
The convective terms contain non-linear quantities.
All three equations are intricately coupled.
There is no equation for pressure.
ME555 : Computational Fluid Dynamics 3 I. Sezai – Eastern Mediterranean University

The staggered grid

Where to store the velocities?
If the velocities and the pressures are both defined at the
nodes of an ordinary CV a highly non-uniform pressure
field can act like a uniform field in the discretized
momentum equations.

A checker-board
pressure field

Suppose that the pressure field is oscillatory as shown above

ME555 : Computational Fluid Dynamics 4 I. Sezai – Eastern Mediterranean University

2
  pE  pP   pP  pW 
p pe  pw   
2   2 
 
x x x
p  pW
 E
2 x
p pN  pS

y 2 y
The pressure at the central node (P) does not appear in above
equations. This gives zero pressure gradients at all nodal points
indicating uniform pressure field. Not realistic.
Solution: Use a staggered grid system for the velocity components.
That is: evaluate scalar variables (p, ρ, T) at ordinary nodal points.
But calculate velocity components (u, v) at cell faces which are
staggered relative to nodal points.

ME555 : Computational Fluid Dynamics 5 I. Sezai – Eastern Mediterranean University

CV’s for u and v are different from the

scalar CV’s of p and T
u(i, j) is defined at west face of p(i, j).
v(i, j) is defined at south face of p(i, j).
This is backward staggered system.
For u-control volume:
p pP  pW

x  xu
For v-control volume:

p pP  pS

y  yv

This arrangement gives non-zero

pressure gradient terms. → Gives
realistic behavior for pressure field.

ME555 : Computational Fluid Dynamics 6 I. Sezai – Eastern Mediterranean University

3
 Non-staggered (Collocated) Grid System
The non-staggered grid system is complicated for unstructured or
body-fitted mesh systems.
Also, the storage of u,v,w and Pressure to four different locations is
inefficient.
In non-staggered grid system all variables are stored at the same
location (point P).

The problem of checker-board pressure field is avoided by calculating

the cell face velocities from interpolation using the momentum
equations (momentum interpolation method).
ME555 : Computational Fluid Dynamics 7 I. Sezai – Eastern Mediterranean University

The momentum equations

The discretised momentum equations at location P for a point (i, j) is
aP P  aEE  aW W  a N N  aSS  S (6.5)
Ss body
C VP  S trans
S S
dc pres
,
 P0VP
S trans
a  , a 
o o
P P
o
P (transient terms), sbody = sCbody +sPbodyP (body forces per unit volume in the differential equation)
t
e y  wy n x  s x
aE   max   Fe , 0 , aW   max  Fw , 0 , aN   max   Fn , 0 , aS   max  Fs , 0
xe xw yn ys
PVP
aP  aE  aW  aN  aS   sPbody xy  F ,
t
F  Fe  Fw  Fn  Fs ,
 p
  x V  [( pe  pw ) / x]V  ( pe  pw )y for x-momentum equation
S pres

  p V  [( pn  ps ) / y]V  ( pn  ps )x for y-momentum equation
 y
S dc   max   Fe , 0 e  P   max   Fe , 0 e  E  (6.6)
 max   Fw , 0 w  P   max   Fw , 0 w  W 
 max   Fn , 0 n  P   max   Fn , 0 n  N 
 max   Fs , 0 s  P   max   Fs , 0 s  S 
Fe   u e y, Fw   u w y, Fn    v n x, Fs    v  s x, (ue , uw , vn , vs are found by MIM method)
e , w , n , s = face values found from a high order (higher than 1st order) convection scheme such as QUICK or CD

ME555 : Computational Fluid Dynamics 8 I. Sezai – Eastern Mediterranean University

4
 The momentum equations
S momentum source term such as body forces
The coefficients aP, aE, aW aN and aS are calculated
by upwind method.
Sdc is the source term resulting from the adoption of
the deferred correction method when any high order
convection scheme, such as QUICK, is used in
estimating the cell face value f.
The coefficients aE, aW etc. contain:
1) Convective flux per unit mass, F
2) Diffusive conductance, D
at control volume faces.
ME555 : Computational Fluid Dynamics 9 I. Sezai – Eastern Mediterranean University

It can be observed that the coefficients of the discretized x- and y-

momentum equations are the same in collocated grid system, provided
that the diffusion coefficient, Γ, is the same in x and y-momentum
equations.
In order to slow down the changes of dependent variables in
consecutive solutions, an under-relaxation factor is introduced into the
discretized equation (6.5) as follows:

 new    (1   ) n 1 (6.7)

where α = under-relaxation factor

n–1 = value of  from the previous iteration
The under-relaxed form of the general equation is
aP
P   aEE  aW W  aN N  aSS  S  
1    a  n1 (6.8)
  P P

ME555 : Computational Fluid Dynamics 10 I. Sezai – Eastern Mediterranean University

5
 Separating the pressure gradient term from the source term,
S  b  S pres
where
b = source term excluding the pressure gradient term
p
S pres   V  [( pe  pw ) / x]V  ( pe  pw ) y for x-momentum equation
x
p
S pres   V  [( pn  ps ) / y ]V  ( pn  ps )x for y-momentum equation
y
Equation (6.8) becomes
  S pres
P    aEE  aW W  aNN  aSS  bP   1    Pn 1 
aP aP
  S pres (6.9)
  E E
 a   aW W  aN N  aSS  BP  
aP
aP
BP  bP 

1    n 1
aPP
 (6.10)

Note that the term 1    P should be added to the source term in
n 1

equations (6.5) and (6.6) if the equations are relaxed.

ME555 : Computational Fluid Dynamics 11 I. Sezai – Eastern Mediterranean University

Momentum Interpolation Method (MIM)

If  stands for u in Eqn. (6.9), the velocity component at nodes P and
E, can be written as
 u  i ai ui  B p   u y  pe  pw  (6.11)
uP  P
 P

a 
u
P P a 
u
P P

 u  i ai ui  B p  E  u y  pe  pw  E (6.12)
uE  
a 
u
P E a 
u
P E

and for the interface velocity at the cell face e

 u   i ai ui  B p e  u y  pE  pP 
ue   (6.13)
a 
u
P e a  u
P e
where the terms on the right-hand side with subscript e should be
interpolated in an appropriate manner. The interface velocity at cell
faces w, n, and s can be obtained similarly.
In Rhie and Chow’s momentum interpolation, the first term and
1/(ap)e in second term of the Eq. (13) are linearly interpolated from
their counterparts in Eqs. (6.11) and (6.12):
ME555 : Computational Fluid Dynamics 12 I. Sezai – Eastern Mediterranean University

6
 Momentum interpolation method (MIM)
 i ai ui  B p  
i ai ui  B p   
i ai ui  B p 

a u   fe 
a u   1  fe  
aPu
  (6.14)
 P e  P E  P
1 1 1
 f e 
 1  f e 
  aPu
e
a u
P E  
aPu
P
(6.15)

where f e is a linear interpolation factor defined as f e  xP /(2 xe )

Substituting  i aiui  Bp aP  terms from Eq’s (6.11), (6.12) and (6.14)
into Eq. (6.13) we obtain
  u y  pE  pP   u y  pe  pw  E 
  f e 

u
aP   aPu   
 


ue   f e uE  1  f e u P   

 
e

y  e w P
p  p
E
 (6.16)

 1  f  u
 
linear interpolation term 

e
aPu
P

  


correction term
The correction term has the function of smoothing the pressure field
(remove the unrealistic pressure field).
ME555 : Computational Fluid Dynamics 13 I. Sezai – Eastern Mediterranean University

Equations (6.13) and (6.16) are essentially equivalent.

Values of F and D for each of the faces e, w, n and s of the control
volume at location (i, j):

Fe  (  uA)e , Fw (i, j )  Fe (i  1, j )
ue  from MIM, e  f e  E  (1  f e ) f P

Fn  (  vA) n , Fs (i, j )  Fn (i, j  1)

vn  from MIM,  n  f n  N  (1  f n ) f P

 e  f e  E  (1  f e ) P  n  f n  N  (1  f n ) P

 e Ae  e Ae
De  , Dn  , Dw (i, j )  De (i  1, j ), Ds (i, j )  Dn (i, j  1)
 xe  xe

ME555 : Computational Fluid Dynamics 14 I. Sezai – Eastern Mediterranean University

7
If a property is unknown at a cell face then a suitable two-
point average is used.
ue,vn, … etc. in fluxes Fe, Fn,…at cell faces are calculated
using MIM.
The variables e, n, …etc. at cell faces in the deferred
correction term Sdc (Eqn.(6.6) ) are calculated using a
convection scheme such as UPWIND or QUICK.
During each iteration the u and v velocity component in F
are those obtained from previous iteration.
Hence, coefficients “ae, an,…” are calculated using the u and
v values from previous iteration.

ME555 : Computational Fluid Dynamics 15 I. Sezai – Eastern Mediterranean University

At each iteration level the values of F are computed

using the u- and v-velocity components resulting
from the previous iteration.
Given a pressure field p, the momentum equations
(6.2) and (6.3) can be written in the discretized form
(6.5) for node P at each location (i, j) and then
solved to obtain the velocity fields.
If the pressure field is correct the resulting velocity
field will satisfy continuity.
As the pressure field is unknown, we need a method
for calculating pressure.
ME555 : Computational Fluid Dynamics 16 I. Sezai – Eastern Mediterranean University

8
 The SIMPLE Algorithm
SIMPLE (Semi-Implicit Method for Pressure-Linked
Equations)
For a guessed pressure field p* the corresponding face
velocity can be written using Eq. (6.13) as
 u  i ai ui*  BPu   u y  pE*  pP* 
u 
* e
  (1   u )uen 1 (6.17)
e
a  u
P e a u
P e

A similar equation can be written for the face velocity vn* .

 v  i ai vi*  BPv   v x  p*N  pP* 

vn*  n
  (1   v )vnn 1 (6.18)
a 
v
P n a 
v
P n

ME555 : Computational Fluid Dynamics 17 I. Sezai – Eastern Mediterranean University

Let p', u', v' be the correction needed to correct the guessed pressure
and velocity fields, i.e.
p  p*  p (6.19)
ue  ue*  ue (6.20)
vn  vn  vn
*
(6.21)
Subtraction of eqn. (6.17) from (6.13) gives
 u  i ai ui  BPu   u y  pE  pP 
ue  e
 (6.22)
a 
u
P e a u
P e

As an approximation, in SIMPLE method the first term in the above

equation is neglected giving
ue  deu  pP  pE  (6.23)
where  u Ae
d eu  , Ae  y (area of CV at face e)
a 
u
P e
ME555 : Computational Fluid Dynamics 18 I. Sezai – Eastern Mediterranean University

9
Similarly
 v An
vn  d nv  pP  pN  d nv 
a 
v
P n
(6.24)
Then the corrected velocities become
ue  ue*  deu  pP  pE  (6.25)
vn  vn*  d nv  pP  pN  (6.26)
Discretizing the continuity equation (6.4) gives

(  ue )y  (  uw )y  (  vn )x  (  vs )x  0 (6.27)

Substituting the corrected face velocities such as that given by Eq’s

(6.24) and (6.25) into Eq. (6.27) gives
aP pP  aW pW  aE pE  aS pS  aN pN  b (6.28)
ME555 : Computational Fluid Dynamics 19 I. Sezai – Eastern Mediterranean University

aP pP  aW pW  aE pE  aS pS  aN pN  b (6.28)

where
aE  (  Ad )e aW  (  Ad ) w aN  (  Ad ) n aS  (  Ad ) s
aP  aw  ae  as  an 6.29)

b  u* A  
w
 
  u * A   v* A   v* A
e
 
s
 n

Note that (u*)w, (v*)s,…, etc are calculated using MIM.

After solving the p' field from Eq. (6.28) the face velocities are
corrected using Eq.'s (6.25)and (6.26) and the pressure field is
corrected by using
p  p*   p p (6.30)

αp = pressure under-relaxation factor (chosen between 0 and

1).
ME555 : Computational Fluid Dynamics 20 I. Sezai – Eastern Mediterranean University

10
Similarly the nodal velocities are corrected using
uP  uP*  d Pu  pw  pe  (6.31)
vP  v  d
*
P
v
P  ps  pn  (6.32)
where
 u Ae  v An
d Pu  and d Pv 
a 
u
P P a 
v
P P

The pressure corrections at the cell faces appearing in Eqs. (6.31) and
(6.32) are calculated by linear interpolation from the nodal values as
pw  f w pW  (1  f w ) pP (6.33)

pe  f e pE  (1  f e ) pP (6.34)

ps  f s pS  (1  f s ) pP (6.35)
pn  f n pN  (1  f n ) pP (6.36)
ME555 : Computational Fluid Dynamics 21 I. Sezai – Eastern Mediterranean University

Boundary Conditions for Pressure

Since there is no equation for the pressure, no
boundary conditions are needed for the pressure at
the near boundary points.

The pressure values at the boundaries can be

calculated by linear extrapolation using the two
near-boundary node pressures.

ME555 : Computational Fluid Dynamics 22 I. Sezai – Eastern Mediterranean University

11
 Boundary Conditions for Pressure Correction Equation
When the velocities at the boundaries are known, there is no need to
correct the velocities at the boundaries in the derivation of the
pressure correction equation. For example if the velocity at the west
boundary is known then for a control volume near the west boundary:
ue  ue*  deu  pP  pE  uw  uwall
vn  vn  d n  pP  pN 
* v vs  vs*  d sv  pS  pP 
Substituting above equations into the discretized continuity equation
(6.27) we obtain the following pressure correction equation for a
control volume near the west boundary
aP pP  aW pW  aE pE  aS pS  aN pN  b (6.37)
where aE  (  Ad )e aW  0 aN  (  Ad ) n aS  (  Ad ) s
    
b    uA wall   u * A   v* A   u* A
e sn

(6.38)
This formulation corresponds to Neuman b.c. (∂p´/∂n = 0) where n is
normal to boundary.
ME555 : Computational Fluid Dynamics 23 I. Sezai – Eastern Mediterranean University

Comparing Eq.'s (6.37-6.38) with (6.28-6.29) for a near boundary

control volume the same definition of the coefficients as used for the
interior points can be used for a near boundary control volume by
setting the corresponding coefficient (aw in this case) to zero and using
uwall in the b term.
As a result no value of pressure correction at the boundary ( pw ) is
involved in this formulation.
However, the value of the pressure correction is needed for correcting
the nodal velocities near boundaries.
For example, for correcting the u-velocity at a nodal point P near a
west boundary, pw at the west boundary is needed in accordance with
equation (6.31).
This value can be obtained by using ∂ p΄/∂n = 0 at the boundary, that
is using p΄(1, j) = p΄(2, j).

ME555 : Computational Fluid Dynamics 24 I. Sezai – Eastern Mediterranean University

12
The SIMPLE Algorithm
Step 1: Solve the discretized momentum equations
u
aP
v
(1   ) aP (1   )
u P   i ai ui  bp  ( Pw  Pe ) Ax  aPu u n 1 vP  i ai vi  bp  ( Ps  Pn ) Ay  aPv v n 1
   
Step 2: Calculate interface velocity ue (Eqn’s (6.16)) and similarly calculate vn
  y  pE  pP   u y  pe  pw  E  u y  pe  pw  P 
 
ue   f e u E  1  f e u P    u  f e 
 1  f e 
  aPu
e
 
aPu
E
 
aPu
P 
Use this velocity to find flux terms, Fe, Fw, etc…
However, e corresponding to ue in Sdc term (Eqn. 6.6) is found from a convection
scheme such as upwind or QUICK.
Step3: Solve pressure correction equation (6.28)
aPp ' pP  aWp ' pW  aEp ' pE  aSp ' pS  aNp ' pN  b p '
Step 4: Correct pressure and velocities at points P using Eqn’s (6.30), (6.31), (6.32)
p  p*   p p uP  uP*  d Pu  pw  pe  vP  vP*  d Pv  ps  pn 
Step 5: Correct face velocities using equations (6.25) and (6.26):
ue  ue*  d eu  pP  pE  vn  vn*  d nv  pP  pN 
Step 6: Solve all other discretized transport equations (i.e. temperature)
aP (1   )
P  aEE  aW W  aNN  aS S  bP  aP n1
 
Step 7: Repeat step 1 to 7 until convergence.

ME555 : Computational Fluid Dynamics 25 I. Sezai – Eastern Mediterranean University

13