Solution of Navier-Stokes.. 6.

11

by aP φP = aE φE + aW φW + aN φN + aS φS + b. The current values of the

dependent variables (φ′ s) are to be evaluated from the known values of the
coefficients (a’s). The evaluation algorithm will be same for the momentum
equations and the pressure correction equation. On a rectangular grid, the
dependent variables at one point (i, j) may be expressed in terms of its neighbors
(Fig. 6.5) as
ai,j ui,j = bi,j ui+1,j + ci,j ui−1,j + di,j ui,j+1 + ei,j ui,j−1 + fi,j
where bi,j is equivalent to aE , ci,j to aW , di,j to aN , ei,j to aS and fi,j to b. The
evaluation of u can be accomplished in the following ways:
(a) Vertical Sweep Upwards may be written in a pseudo FORTRAN code as
DO 10 J = 2, M - 1
DO 10 I = 2, L - 1

10 ai,j ui,j = bi,j ui+1,j + ci,j ui−1,j + (di,j ūi,j+1 + ei,j ūi,j−1 + fi,j )
Here ū is the currently available value in storage and all the coefficients including
fi,j are known. for each j, we shall get a system of equations if we susbtitute
i = 2......L − 1. In other words, for each j, a tridiagonal matrix is available
which can be solved for all i row of points at that j. Once one complete row is
evaluated for any particular j, the next j will be taken up, and so on.
(b) Horizontal Sweep Forward may be written in pseudo FORTRAN code as
DO 20 I = 2, L - 1
DO 20 J = 2, M - 1

20 ai,j ui,j = di,j ui,j+1 + ei,j ui,j−1 + (bi,j ūi+1,j + ci,j ūi−1,j + fi,j )
Again ū is the currently available value in storage from previous calculations.
For each i, we get a system of equation if we substitute j = 2....M − 1. A
tridiagonal matrix is available for each i. Once one complete column of points
are evaluated for any particular i, the next i will be taken up and so on. The
vertical sweep upward and downward are repeated. Similarly the horizontal
sweep forward and rearward are also repeated until convergence is achieved.
For solving tridiagonal system, the tridiagonal matrix algorithm (TDMA) due
to Thomas (1949) is deployed. The above mentioned evaluation procedure is
known as line-by-line TDMA.

6.4 Solution of the Unsteady Navier-Stokes Equa-

tions
6.4.1 Introduction to the MAC method
The MAC method of Harlow and Welch is one of the earliest and most useful
methods for solving the Navier-Stokes equations. This method necessarily deals
6.12 Computational Fluid Dynamics

M
N i,j+1
M−1
W P i,j E
i−1,j i+1,j
S
2
i,j−1
j
1
1 2 L−1 L

Figure 6.5: Dependent variable at (i, j) is expressed in terms of neigh-

boring points.

with a Poisson equation for the pressure and momentum equations for the
computation of velocity. It was basically developed to solve problems with
free surfaces, but can be applied to any incompressible fluid flow problem. A
modified version of the original MAC method due to Hirt and Cook (1972) has
been used by researchers to solve a variety of flow problems.
The text discusses the modified MAC method and highlights the salient
features of the solution algorithm so that the reader is able to write a com-
puter program with some confidence. The important ideas on which the MAC
algorithm is based are:
1. Unsteady Navier-Stokes equations for incompressible flows in weak con-
servative form and the continuity equation are the governing equations.
2. Description of the problem is elliptic in space and parabolic in time. Solu-
tion will be marched in the time direction. At each time step, a converged
solution in space is obtained but this converged solution at a time step
may not be the solution of the physical problem.
3. If the problem is steady, in its physical sense, then after a finite number
of steps in time direction, two consecutive time steps will show identical
solutions. However, in a machine-computation this is not possible hence a
very small upper bound, say, “STAT” is predefined. Typically, STAT may
be chosen between 10−3 and 10−5 . If the maximum discrepency of any
of the velocity components for two consecutive time steps at any location
over the entire space does not exceed STAT, then it can be said that the
steady solution has evolved.
 Solution of Navier-Stokes.. 6.13

4. If the physical problem is basically unsteady in nature, the aforesaid maxi-

mum discrepancy of any dependant variable for two consecutive time steps
will never be less than STAT. However, for such a situation, a specified
velocity component can be stored over a long duration of time and plot of
the velocity component against time (often called as signal) depicts the
character of the flow. Such a flow may be labelled simply as “unsteady”.
5. With the help of the momentum equations, we compute explicitly a pro-
visional value of the velocity components for the next time step.

Consider the weak conservative form of the nondimensional momentum

equation in the x-direction:
∂u ∂(u2 ) ∂(uv) ∂(uw) ∂p 1 2
+ + + =− + ∇ u
∂t ∂x ∂y ∂z ∂x Re

It is assumed that at t = nth level, we have a converged solution. Then

for the next time step

ũn+1 n
i,j,k − ui,j,k n
= [CON DIF U − DP DX]i,j,k
∆t
or
n
ũn+1 n
i,j,k = ui,j,k + ∆t [CON DIF U − DP DX]i,j,k (6.19)
n
[CON DIF U − DP DX]i,j,k consists of convective and diffusive terms,
n+1
and the pressure gradient. Similarly, the provisional values for ṽi,j,k and
n+1
w̃i,j,k can be explicitly computed. These explicitly advanced velocity com-
ponents may not constitute a realistic flow field. A divergence-free velocity
field has to exist in order to describe a plausible incompressible flow situa-
tion. Now, with these provisional ũn+1 n+1 n+1
i,j,k , ṽi,j,k and w̃i,j,k values, continuity
equation is evaluated in each cell. If (∇ · V ) produces a nonzero value,
there must be some amount of mass accumulation or annihilation in each
cell which is not physically possible. Therefore the pressure at any cell is
directly linked with the value of the (∇·V ) of that cell. Now, on one hand
the pressure has to be calculated with the help of the nonzero divergence
value and on the other, the velocity components have to be adjusted.
The correction procedure continues through an iterative cycle untill the
divergence-free velocity field is ensured. Details of the procedure will be
discussed in the subsequent section.
6. Boundary conditions are to be applied after each explicit evaluation for
the time step is accomplished. Since the governing equations are elliptic in
space, boundary conditions on all confining surfaces are required. More-
over, the boundary conditions are also to be applied after every pressure-
velocity iteration. The five types of boundary conditions to be considered
6.14 Computational Fluid Dynamics

are rigid no-slip walls, free-slip walls, inflow and outflow boundaries, and
periodic (repeating) boundaries.

6.4.2 MAC Formulation

The region in which computations are to be performed is divided into a set
of small cells having edge lengths δx, δy and δz (Fig. 6.6). With respect to
this set of computational cells, velocity components are located at the centre
of the cell faces to which they are normal and pressure and temperature are
defined at the centre of the cells. Cells are labeled with an index (i,j,k) which
denotes the cell number as counted from the origin in the x, y and z directions
respectively. Also pi,j,k is the pressure at the centre of the cell (i,j,k), while
ui,j,k is the x-direction velocity at the centre of the face between cells (i,j,k) and
(i + 1, j, k) and so on (Fig. 6.7). Because of the staggered grid arrangements,
the velocities are not defined at the nodal points, but whenever required, they
are to be found by interpolation. For example, with uniform grids, we can write
1
ui−1/2,j,k = [ui−1,j,k + ui,j,k ]. Where a product or square of such a quantity
2
appears, it is to be averaged first and then the product to be formed.

Figure 6.6: Discretization of a three-dimensional domain.

Convective terms are discretized using a weighted averaged of second upwind

and space centered scheme (Hirt et al., 1975). Diffusive terms are discretized
 Solution of Navier-Stokes.. 6.15

by a central differencing scheme. Let us consider the discretized terms of the

x-momentum equation (Figure 6.7):
vi,j+1,k

ui,j+1,k
vi,j,k
vi+1,j.k
δy
θi,j.k
p ui,j.k ui+1,j,k
ui−1,j,k i,j.k vi+1,j−1,k
wi,j.k−1 vi,j−1,k
δz

δx ui,j−1,k

Figure 6.7: Three-dimensional staggered grid showing the locations

of the discretized variables.

∂(u2 ) 1
= [(ui,j,k + ui+1,j,k )(ui,j,k + ui+1,j,k )
∂x 4δx
+ α|(ui,j,k + ui+1,j,k )|(ui,j,k − ui+1,j,k )
− (ui−1,j,k + ui,j,k )(ui−1,j,k + ui,j,k )
− α|(ui−1,j,k + ui,j,k )|(ui−1,j,k − ui,j,k )
≡ DU U DX

∂(uv) 1
= [(vi,j,k + vi+1,j,k )(ui,j,k + ui,j+1,k )
∂y 4δy
+ α|(vi,j,k + vi+1,j,k )|(ui,j,k − ui,j+1,k )
− (vi,j−1,k + vi+1,j−1,k )(ui,j−1,k + ui,j,k )
− α|(vi,j−1,k + vi+1,j−1,k )|(ui,j−1,k − ui,j,k )
≡ DU V DY
6.16 Computational Fluid Dynamics

∂(uw) 1
= [(wi,j,k + wi+1,j,k )(ui,j,k + ui,j,k+1 )
∂z 4δz
+ α|(wi,j,k + wi+1,j,k )|(ui,j,k − ui,j,k+1 )
− (wi,j,k−1 + wi+1,j,k−1 )(ui,j,k−1 + ui,j,k )
− α|(wi,j,k−1 + wi+1,j,k−1 )|(ui,j,k−1 − ui,j,k )
≡ DU W DZ

∂p pi+1,j,k − pi,j,k
= ≡ DP DX
∂x δx

∂ 2u ui+1,j,k − 2ui,j,k + ui−1,j,k

= 2 ≡ D2U DX2
∂x2 (δx)
∂ 2u ui,j+1,k − 2ui,j,k + ui,j−1,k
= ≡ D2U DY 2
∂y 2 (δy)2
∂ 2u ui,j,k+1 − 2ui,j,k + ui,j,k−1
2
= 2 ≡ D2U DZ2
∂z (δz)

with

α −→ 1 Scheme −→ Second Upwind

α −→ 0 Scheme −→ Space centred

Factor α is chosen in such a way that the differencing scheme retains “some-
thing” of second-order accuracy and the required upwinding is done for the
sake of stability. A typical value of α is between 0.2 and 0.3. As mentioned
earlier, the quantity ũn+1
i,j,k is now evaluated explicitly from the discretized form
of Equation (6.2) as
n
ũn+1 n
i,j,k = ui,j,k +δt [CON DIF U − DP DX]i,j,k

where

[CON DIF U − DP DX]ni,j,k = [(−DU U DX − DU V DY − DU W DZ)

− DP DX + (1/Re)(D2U DX2
+ D2U DY 2 + D2U DZ2)]

Similarly, we evaluate
n+1 n n
ṽi,j,k = vi,j,k +δt [CON DIF V − DP DY ]i,j,k (6.20)
n+1 n n
w̃i,j,k = wi,j,k +δt [CON DIF W − DP DZ]i,j,k (6.21)
 Solution of Navier-Stokes.. 6.17

As discussed earlier, the explicitly advanced tilde velocities may not necessarily
lead to a flow field with zero mass divergence in each cell. This implies that, at
this stage the pressure distribution is not correct, the pressure in each cell will
be corrected in such a way that there is no net mass flow in or out of the cell.
In the original MAC method, the corrected pressures were obtained from the
solution of a Poisson equation for pressure. A related technique developed by
Chorin (1967) involved a simultaneous iteration on pressure and velocity com-
ponents. Vieceli (1971) showed that the two methods as applied to MAC are
equivalent. We shall make use of the iterative correction procedure of Chorin
(1967) in order to obtain a divergence-free velocity field. The mathematical
methodology of this iterative pressure-velocity correction procedure will be dis-
cussed herein.

The relationship between the explicitly advanced velocity component and ve-
locity at the previous time step may be written as

[pni,j,k − pni+1,j,k ]
ũn+1 n
i,j,k = ui,j,k + δt + δt [CON DIF U ]ni,j,k (6.22)
∆x
where [CON DIF U ]ni,j,k is only the contribution from convection and diffusion
terms. On the other hand, the corrected velocity component (unknown) will be
related to the corrected pressure (also unknown) in the following way:

[pn+1 n+1
i,j,k − pi+1,j,k ]
un+1 n
i,j,k = ui,j,k + δt + δt [CON DIF U ]ni,j,k (6.23)
∆x
From Equations (6.22) and (6.23)

[pci,j,k − pci+1,j,k ]
un+1 n+1
i,j,k − ũi,j,k = δt
δx
where the pressure correction may be defined as

pci,j,k = pn+1 n
i,j,k − pi,j,k

Neither the pressure correction nor is the quantity un+1

i,j,k explicitly known at
this stage. Hence, one cannot be calculated without the help of the other.
Calculations are done in an iterative cycle and we write

corrected Estimated Correction

[pci,j,k − pci+1,j,k ]
un+1 n+1
i,j,k −→ ũi,j,k + δt
δx
6.18 Computational Fluid Dynamics

In a similar way, we can formulate the following array:

[pci,j,k − pci+1,j,k ]
un+1 n+1
i,j,k −→ ũi,j,k + δt (6.24)
δx
[pci,j,k − pci−1,j,k ]
un+1 n+1
i−1,j,k −→ ũi−1,j,k − δt (6.25)
δx
n+1 n+1
[pci,j,k − pci,j+1,k ]
vi,j,k −→ ṽi,j,k + δt (6.26)
δy
n+1 n+1
[pi,j,k − pci,j−1,k ]
c
vi,j−1,k −→ ṽi,j−1,k − δt (6.27)
δy
n+1 n+1
[pci,j,k − pci,j,k+1 ]
wi,j,k −→ w̃i,j,k + δt (6.28)
δz
n+1 n+1
[pci,j,k − pci,j,k−1 ]
wi,j,k−1 −→ w̃i,j,k−1 − δt (6.29)
δz
The correction is done through the continuity equation. Plugging-in the above
relationship into the continuity equation (6.1) yields
" n+1
ui,j,k − un+1 n+1 n+1 n+1 n+1
#
i−1,j,k vi,j,k − vi,j−1,k wi,j,k − wi,j,k−1
+ +
δx δy δz
" n+1
ũi,j,k − ũn+1 n+1 n+1 n+1 n+1
#
i−1,j,k ṽi,j,k − ṽi,j−1,k w̃i,j,k − w̃i,j,k−1
= + +
δx δy δz
 c c c
pi+1,j,k − 2pi,j,k + pi−1,j,k pi,j+1,k − 2pi,j,k + pci,j−1,k
c c
−δt 2 + 2
(δx) (δy)
pci,j,k+1 − 2pci,j,k + pci,j,k−1

+ (6.30)
(δz)2

or

un+1 n+1 n+1 n+1 n+1 n+1

" #
i,j,k − ui−1,j,k vi,j,k − vi,j−1,k wi,j,k − wi,j,k−1
+ +
δx δy δz
ũn+1 ũn+1 n+1 n+1 n+1 n+1
" #
i,j,k − i−1,j,k ṽi,j,k − ṽi,j−1,k w̃i,j,k − w̃i,j,k−1
= + +
δx δy δz
2 δt (pci,j,k ) 2 δt (pci,j,k ) 2 δt (pci,j,k )
+ 2 + 2 + (6.31)
δx δy δz 2
In deriving the above expression, it is assumed that the pressure corrections in
the neighboring cells are zero. Back to the calculations, we can write
  
1 1 1
0 = (Div)i,j,k + pci,j,k 2δt + +
δx2 δy 2 δz 2
 Solution of Navier-Stokes.. 6.19

or

−(Div)i,j,k
pci,j,k = h  i (6.32)
1 1 1
2δt δx2
+ δy 2
+ δz 2

In order to accelerate the calculation, the pressure correction equation is mod-

ified as
−ω0 (Div)i,j,k
pci,j,k = h  i (6.33)
2δt δx12 + δy12 + δz12

where ω0 is the overrelaxation factor. A value of ω0 = 1.7 is commonly used.

The value of ω0 giving most rapid convergence, should be determined by nu-
merical experimentation. After calculating pci,j,k , the pressure in the cell (i, j, k)
is adjusted as
pn+1 n c
i,j,k −→ pi,j,k + pi,j,k (6.34)
Now the pressure and velocity components for each cell are corrected through
an iterative procedure in such a way that for the final pressure field, the velocity
divergence in each cell vanishes. The process is continued till a divergence-free
velocity is reached with a prescribed upper bound; here a value of 0.0001 is
recommended.
Finally, we discuss another important observation. If the velocity bound-
ary conditions are correct and a divergence-free converged velocity field has
been obtained, eventually correct pressure will be determined in all the cells at
the boundary. Thus, this method avoids the application of pressure boundary
conditions. This typical feature of modified MAC method has been discussed
in more detail by Peyret and Taylor (1983). However, it was also shown by
Brandt, Dendy and Ruppel (1980) that the aforesaid pressure-velocity itera-
tion procedure of correcting pressure is equivalent to the solution of Poisson
equation for pressure. As such from Eqn. (6.30) we can directly write as

(Div)i,j,k
∇2 (pci,j,k ) = (6.35)
δt
The Eqn. (6.35) can be solved implicitly using appropriate boundary condi-
tions for pc at the confining boundaries.

6.4.3 Boundary Conditions

So far we have not discussed the boundary conditions. However, they are
imposed by setting appropriate velocities in the fictitious cells surrounding the
physical domain (Figure 6.8).
Consider, for example, the bottom boundary of the computational (physical)
mesh. If this boundary is to be a rigid no-slip wall, the normal velocity on the
wall must be zero and the tangential velocity components should also be zero.
Here we consider a stationary wall. With reference to the Figure 6.8, we have
6.20 Computational Fluid Dynamics


vi,1,k = 0
for i = 2 to ire

ui,1,k = −ui,2,k
 and k = 2 to kre
wi,1,k = −wi,2,k
If the left side-wall is a free-slip (vanishing shear) boundary, the normal velocity
must be zero and the tangential velocities should have no normal gradient.

wi,j,1 = 0 
for i = 2 to ire
ui,j,1 = ui,j,2
 and j = 2 to jre
vi,j,1 = vi,j,2
If the front plane is provided with inflow boundary conditions, it should be spec-
ified properly. Any desired functional relationship may be recommended. Gen-
erally, normal velocity components are set to zero and a uniform or parabolic
axial velocity may be deployed. Hence with reference to Fig. 6.8, we can write

v1,j,k = −v2,j,k 

w1,j,k = −w2,j,k for j = 2 to jre

u1,j,k = 1.0 or 
 and k = 2 to kre
u1,j,k = 1.5 1 − ((jm − j)/jm )2
 

where jm is the horizontal midplane.

Continuative or outflow boundaries always pose a problem for low-speed

Figure 6.8: Boundary conditions and fictitious boundary cells.

calculations, because whatever prescription is chosen it can affect the entire

flow upstream. What is needed is a prescription that permits fluid to flow out
 Solution of Navier-Stokes.. 6.21

of the mesh with a minimum of upstream influence. Commonly used conditions

for such a boundary is ∇V · n = 0, where n is the unit normal vector.
The boundary condition that has more generality at the outflow is described
by Orlanski (1976). This condition allows changes inside the flow field to be
transmitted outward, but not vice-versa:
∂Ψ ∂Ψ
+ Uav =0
∂t ∂x
where Uav is the average velocity at the outflow plane and Ψ represents u, v, w
or any dependent variable.

6.4.4 Numerical Stability Considerations

For accuracy, the mesh size must be chosen small enough to resolve the expected
spatial variations in all dependent variables. Once a mesh has been chosen,
the choice of the time increment is governed by two restrictions, namely, the
Courant-Fredrichs-Lewy (CFL) condition and the restriction on the basis of
grid-Fourier numbers. According to the CFL condition, material cannot move
through more than one cell in one time step, because the difference equations
assume fluxes only between the adjacent cells. Therefore, the time increment
must satisfy the inequality.
 
δx δy δz
δt < min , , (6.36)
|u| |v| |w|
where the minimum is with respect to every cell in the mesh. Typically, δt
is chosen equal to one-fourth to one-third of the minimum cell transit time.
When the viscous diffusion terms are more important, the condition necessary
to ensure stability is dictated by the restriction on the Grid-Fourier numbers,
which results in
1 (δx2 δy 2 δz 2 )
νδt < · (6.37)
2 (δx2 + δy 2 + δz 2 )
in dimensional form. After non-dimensionilization, this leads to
1 (δx2 δy 2 δz 2 )
δt < · Re (6.38)
2 (δx2 + δy 2 + δz 2 )
The final δt for each time increment is the minimum of the δt’s obtained from
Equations (6.36) and (6.38)
The last quantity needed to ensure numerical stability is the upwind pa-
rameter α. In general, α should be slightly larger than the maximum value of
|uδt/δx] or |vδt/δy] occurring in the mesh, that is,
 
uδt vδt wδt
max , , ≤α<1 (6.39)
δx δy δz
As a ready prescription, a value between 0.2 and 0.4 can be used for α. If α
is too large, an unnecessary amount of numerical diffusion (artificial viscosity)
will be introduced.
6.22 Computational Fluid Dynamics

6.4.5 Higher-Order Upwind Differencing

More accurate solutions are obtained if the convective terms are discretized by
higher-order schemes. Davis and Moore (1982) use the MAC method with a
multidimensional third-order upwinding scheme. Needless to mention that their
marching algorithm for the momentum equation is explicit and the stability re-
striction concerning the CFL condition [uδt/δx ≤ 1 and vδt/δy ≤ 1] is satisfied.

V i−1,j+1 V i,j+1 V i+1,j+1

U i−1,j+1 φi,j+1 U i,j+1 U i+1,j+1

V i−1,j V i,j V i+1,j

1111111
0000000
φ i−1,j 0000000
1111111
φ φ i+1,j
0000000
1111111
U i−1,j i,j U i,j U i+1,j

v, j, y 0000000
1111111
0000000
1111111
V i−1,j−1 V i,j−1 V i+1,j−1

0000000
1111111
U i−1,j−1 φi,j−1 U i,j−1 U i+1,j−1

u, i, x

Figure 6.9: Dependent variables(u, v and φ) on a rectangular grid.

The multidimensional third-order upwinding is in principle similar to one di-

mensional quadratic upstream interpolation scheme introduced by Leonard
(1979). Consider Fig. 6.9. Let φ be any property which can be convected
∂(uφ)
and diffused. The convective term may be represented as
∂x
∂(uφ) (uφ)i+ 1 ,j − (uφ)i− 1 ,j
= 2 2
(6.40)
∂x δx
where the variables φi+ 12 ,j and φi− 21 ,j are defined as

ξ
φi+ 21 ,j = 0.5(φi,j + φi+1,j ) − (φi−1,j − 2φi,j + φi+1,j ) for ui,j ≤ 0 (6.41)
3
and
ξ
φi− 21 ,j = 0.5(φi,j + φi−1,j ) − (φi−2,j − 2φi−1,j + φi,j ) for ui,j > 0 (6.42)
3
 Solution of Navier-Stokes.. 6.23

The parameter ξ can be chosen to increase the accuracy or to alter the diffusion-
like characteristics. It may be pointed out ξ = 3/8 corresponds to the QUICK
scheme of Leonard(1979).
Let us consider two-dimensional momentum equation in weak conservative form
which is given by

∂u ∂(u2 ) ∂(uv) 1 ∂2u ∂2u

 
∂p
+ + =− + + (6.43)
∂t ∂x ∂y ∂x Re ∂x2 ∂y 2

In non-conservative form this may be written as

1 ∂2u ∂2u
 
∂u ∂u ∂u ∂p
+u +v =− + + 2 (6.44)
∂t ∂x ∂y ∂x Re ∂x2 ∂y

Here we introduce a term transport-velocity. The transport velocities for the

second and third terms on the left hand side are u and v respectively. While
dealing with the equations in the conservative form, we shall keep this in mind.
For example, during discretization of the term ∂(uv)/∂y of Eq. (6.43) we should
remember that v is the transport-velocity associated with this term. It is cus-
tomary to define the transport-velocity at the nodal point where the equation
is being defined. In case of the term ∂(uv)/∂y we have to refer to Fig. 6.10 and
write down the product term uv as

(uv)i,j = 0.25ui,j (vi,j + vi+1,j + vi,j−1 + vi+1,j−1 ) (6.45)

Finally the discretization of the term ∂(uv)/∂y for the x-momentum equation
will be accomplished in the following way:
∂(uv) 1
=0.25
∂y 8δy
[3ui,j+1 (vi,j+1 + vi,j + vi+1,j+1 + vi+1,j )
+ 3ui,j (vi,j + vi,j−1 + vi+1,j + vi+1,j−1 )
− 7ui,j−1 (vi,j−1 + vi,j−2 + vi+1,j−1 + vi+1,j−2 )
+ ui,j−2 (vi,j−2 + vi,j−3 + vi+1,j−2 + vi+1,j−3 )] f or V > 0 (6.46)

∂(uv) 1
=0.25
∂y 8δy
[−ui,j+2 (vi,j+2 + vi,j+1 + vi+1,j+2 + vi+1,j+1 )
+ 7ui,j+1 (vi,j+1 + vi,j + vi+1,j+1 + vi+1,j )
− 3ui,j (vi,j + vi,j−1 + vi+1,j + vi+1,j−1 )
− 3ui,j−1 (vi,j−1 + vi,j−2 + vi+1,j−1 + vi+1,j−2 )] f or V < 0 (6.47)

where
V = vi,j + vi+1,j + vi,j−1 + vi+1,j+1
6.24 Computational Fluid Dynamics

vi,j vi+1,j

ui,j

vi,j−1 vi+1,j−1
v,j,y

u,i,x

Figure 6.10: Definition of the transport velocity at a point where the

momentum equation is being discretized.

6.4.6 Sample Results

For unsteady laminar flow past a rectangular obstacle in a channel, Mukhopad-
hyay, Biswas and Sundararajan (1992) use the MAC algorithm to explicitly
march in time. Their results corroborated with the experimental observation
of Okajima (1982). A typical example of numerical flow visualization depicting
the development of von Kármán vortex street is illustrated in Fig. 6.11. The
cross-stream vorticity contours vectors behind a delta-winglet placed inside a
channel are shown in Fig. 6.12. These results were obtained by Biswas, Torii et
al. (1996) who used MAC to solve for a three-dimensional flow field in a chan-
nel containing delta-winglet as a vortex generator. The MAC algorithm has
been extensively used by the researchers to solve flows in complex geometry.
Braza, Chassaing and Ha-Minh (1986) investigated the dynamic characteris-
tics of the pressure and velocity fields of the unsteady wake behind a circular
cylinder using MAC algorithm. Robichaux, Tafti and Vanka (1992) deployed
MAC algorithm for Large Eddy Simulation (LES) of turbulent channel flows.
Of course, they performed the time integration of the discretized equations by
using a fractional step method (Kim and Moin, 1985). Another investigation by
Kim and Benson (1992) suggests that the MAC method is significantly accurate
and at the same time the computational effort is reasonable.
 Solution of Navier-Stokes.. 6.25

Figure 6.11: Streamlines crossing the cylinder in the duct: ReB =

162

6.5 Solution of Energy Equation

The energy equation for incompressible flows, neglecting mechanical work and
gas radiation, may be written as
 
∂T ∗ ∂T ∗ ∂T ∗ ∂T
ρcp + u + v + w = k∇2 T + µφ∗ (6.48)
∂t∗ ∂x∗ ∂y ∗ ∂z ∗

Figure 6.12: (a) & (b) Vorticity contours behind a cross-plane Xc =

2.5 behind the winglet at β = 15◦ (a) experiment, (b) computation
6.26 Computational Fluid Dynamics

Figure 6.12(c) Formation of complex vortex system due to the winglet

where φ∗ is the viscous dissipation given as

" 2  ∗ 2  ∗ 2 #  ∗ 2
∗ ∂u∗ ∂v ∂w ∂v ∂u∗
φ =2 + + + +
∂x∗ ∂y ∗ ∂z ∗ ∂x∗ ∂y ∗
 ∗ 2  ∗ 2
∂w ∂v ∗ ∂w ∂u∗
+ + ∗ + + ∗
∂y ∗ ∂z ∂x∗ ∂z
Equation (6.48) may be non-dimensionalized in the following way:
u∗ v∗ w∗ T − T∞
u= , v= , w= , θ=
U∞ U∞ U∞ Tw − T∞
x∗ y∗ z∗ t∗
x= , v= , z= , t=
L L L L/U∞
Substituting the above variables in equation (6.48) we obtain
 
ρcp U∞ (Tw − T∞ ) ∂θ ∂θ ∂θ ∂θ
+u +v +w
L ∂t ∂x ∂y ∂z
 2 2 2
 2
(Tw − T∞ )k ∂ θ ∂ θ ∂ θ µU∞
= 2 2
+ 2+ 2 + 2
φ (6.49)
L ∂x ∂y ∂z L
where φ is the nondimensional form of φ∗ . Finally, the normalized energy
equation becomes
1 ∂2θ ∂2θ ∂2θ
 
∂θ ∂θ ∂θ ∂θ Ec
+u +v +w = + + + φ (6.50)
∂t ∂x ∂y ∂z P e ∂x2 ∂y 2 ∂z 2 Re
where Pe, the Peclet number is given as
1 (Tw − T∞ )k L
= ·
Pe L2 ρcp U∞ (Tw − T∞ )
1 k k µ 1 1
or = = · = ·
Pe Lρcp U∞ µcp ρLU∞ P r Re
 Solution of Navier-Stokes.. 6.27

Further, Ec, the Eckert number is

2 2
Ec µU∞ L U∞ 1
= 2
· = ·
Re L ρcp U∞ (Tw − T∞ ) cp (Tw − T∞ ) ρU∞ L/µ

6.5.1 Retention of Dissipation

The dissipation term φ is frequently neglected while solving the energy equation
for incompressible flows. As the Mach number M → 0, Ec → 0. However, even
at a low Mach number, φ can be important if (Tw − T∞ ) is very small. Let us
look at these aspects. Since
2
 
U∞ 1 cp T ∞ T w
Ec = , = 2
− 1
cp (Tw − T∞ ) Ec U∞ T∞

and

cp T ∞ cp γRT∞
2
= 2
U∞ γRU∞

where R is the gas constant = cp − p cv , and γ = cp /cv .

Let the local acoustic velocity C = γRT∞ , and Mach number M∞ = U∞ /C
Then,
   
cp T ∞ cp 1 cp 1 1 1
2
= 2
=   2
= 2
U∞ γ(cp − cv ) M∞ γcp 1 − γ 1 M∞ (γ − 1) M∞

Hence,
 
1 1 Tw
= 2
−1
Ec (γ − 1)M∞ T∞

or
2
(γ − 1)M∞
Ec =
((Tw /T∞ ) − 1)

In general for incompressible flows M∞ < 0.3 and γ ≥ 1. Hence Ec is small.

But for very small temperature difference, i.e., if Tw /T∞ is slightly larger than
1, Ec might assume a large value and importance of dissipation arises.
However, for computing incompressible convective flows, the viscous dissipa-
tion is neglected in this chapter and we continue with the steady state energy
equation.
6.28 Computational Fluid Dynamics

6.6 Solution Procedure

The steady state energy equation, neglecting the dissipation term, may be writ-
ten in the following conservative form as
1 ∂2θ ∂2θ ∂2θ
 
∂(uθ) ∂(vθ) ∂(wθ)
+ + = + 2+ 2 (6.51)
∂x ∂y ∂z P e ∂x2 ∂y ∂z
Equation (6.51) may be written as
m
∇2 θ = P e [CON V T ]i,j,k (6.52)
m
where [CON V T ]i,j,k is the discretized convective terms on the left-hand side of
Equation (6.51) and m stands for the iterative counter. To start with, we can
assume any guess value of θ throughout the flow field. Since u, v, w are known
from the solution of momentum equation hence Equation 6.51 is now a linear
equation. However, from the guess value of θ and known correct values of u,
v and w the left-hand side of Equation 6.51 is evaluated. A weighted average
scheme or QUICK scheme may be adapted for discretization of the convective
terms. After discretizing and evaluating right-hand side of Equation (6.52) we
obtain a Poisson equation for the temperature with a source term on the right
hand side. Now, we shall follow SOR technique for solving Equation (6.52).
Consider a discretized equation as
θi+1,j,k − 2θi,j,k + θi−1,j,k θi,j+1,k − 2θi,j,k + θi,j−1,k
2 +
(δx) (δy)2
θi,j,k+1 − 2θi,j,k + θi,j,k−1
+ 2 = S ∗m
(δz)

where
m
S ∗m ≡ P e [CON V T ]i,j,k

One of the terms of [CONVT] may be expressed as (refer to Fig. 6.7)

∂(uθ) 1
= [ui,j,k (θi,j,k + θi+1,j,k ) + α|ui,j,k | (θi,j,k − θi+1,j,k )
∂x 2δx

− ui−1,j,k (θi−1,j,k + θi,j,k ) − α|ui−1,j,k | (θi−1,j,k − θi,j,k )]

Then
 
2 2 2
A∗m − θi,j,k + 2+ 2 = S ∗m
δx2 δy δz
 Solution of Navier-Stokes.. 6.29

or

A∗m − S ∗m
θi,j,k = h i (6.53)
2 2 2
δx2
+ δy 2
+ δz 2

where
m m m m m m
θi+1,j,k + θi−1,j,k θi,j+1,k + θi,j−1,k θi,j,k+1 + θi,j,k−1
A∗m = + +
(δx)2 (δy)2 (δz)2

θi,j,k in Equation (6.53) may be assumed to be the most recent value and it
m′
may be written as θi,j,k . In order to accelerate the speed of computation we
introduce an overrelaxation factor ω. Thus
h ′ i
m+1 m m m
θi,j,k = θi,j,k + ω θi,j,k − θi,j,k (6.54)

′
m m m+1
where θi,j,k is the previous value, θi,j,k the most recent value and θi,j,k the cal-
culated better guess. The procedure will continue till the required convergence
is achieved. This is equivalent to Gauss-Seidel procedure for solving a system
of linear equations.

References
1. Biswas, G., Torii, K., Fujii, D., and Nishino, K., Numerical and Exper-
imental Determination of Flow Structure and Heat Transfer Effects of
Longitudinal Vortices in a Channel Flow, Int. J. Heat Mass Transfer,
Vol. 39, pp. 3441-3451, 1996.
2. Biswas, G., Breuer, M. and Durst, F., Backward-Facing Step Flows for
Various Expansion Ratios at Low and Moderate Reynolds Numbers, Jour-
nal of Fluids Engineering (ASME), Vol. 126, pp. 362-374, 2004.
3. Brandt, A., Dendy, J.E and Rupppel, H., The Multigrid Method for Semi-
Implicity Hydrodynamics Codes, J. Comput. Phys, Vol. 34, pp. 348-370,
1980.
4. Braza, M., Chassaing P. and Ha Minh, H., Numerical Study and Phys-
ical Analysis of the Pressure and Velocity Fields in the Near-Wake of a
Circular Cylinder, J. Fluid Mech., Vol. 165, pp. 79-130, 1986.
5. Chorin, A.J., Numerical Method for Solving Incompressible Viscous Flow
Problems, J. Comput. Phys., Vol. 2, pp. 12-26, 1967.
6. Davis, R.W. and Moore, E. F., A Numerical Study of Vortex Shedding
from Rectangles, J. Fluid Mech., 116: 475-506, 1982.
6.30 Computational Fluid Dynamics

7. Harlow, F.H. and Welch, J.E., Numerical Calculation of Time-dependent

Viscous Incompressible Flow of Fluid with Free Surfaces, Phys. of Fluids,
Vol. 8, pp. 2182-2188, 1965.
8. Hirt, C.W. and Cook, J.L., Calculating Three-Dimensional Flows Around
Structures and Over Rough Terrain, J. Comput. Phys., Vol. 10, pp. 324-
340, 1972.
9. Hirt, C.W., Nicholas, B.D. and Romero,N.C., Numerical Solution Algo-
rithm for Transient Fluid Flows, LA-5852, Los Alamos Scientific Labora-
tory Report, 1975.
10. Kim, J. and Moin, P., Application of Fractional Step Method to Incom-
pressible Navier-Stokes Equation, J. Comput. Phys., Vol. 59, pp. 308-
323, 1985.
11. Kim, S.W. and Benson, T.J., Comparison of the SMAC, PISO and Itera-
tive Time Advancing Schemes for Unsteady Flows, Computers & Fluids,
Vol. 21, pp. 435-454, 1992.
12. Leonard, B.P., A Stable and accurate Convective Modelling Procedure
Based on Quadratic Upstream Interpolation, Comp. Methods Appl. Mech.
Engr., Vol. 19, pp. 59-98, 1979.
13. Mukhopadhyay, A., Biswas, G., and Sundararajan, T., Numerical Inves-
tigation of Confined Wakes behind a Square Cylinder in a Channel, Int.
J. Numer. Methods Fluids, Vol. 14, pp. 1473-1484, 1992.
14. Okajima, A., Strouhal Numbers of Rectangular Cylinders, J. Fluid Mech.,
Vol. 123, pp. 379-398, 1982.
15. Orlanski, I., A Simple Boundary Condition for Unbounded Hyperbolic
Flows, J. Comput Phys., Vol. 21, pp. 251-269, 1976.
16. Patankar, S.V., A Calculation Procedure for Two-Dimensional Elliptic
Situations, Numerical Heat Transfer, Vol. 4, pp. 409-425, 1981.
17. Patankar, S.V., Numerical Heat Transfer and Fluid Flow, Hemisphere
Publishing Co., 1980.
18. Patankar, S.V. and Spalding, D.B., A Calculation Procedure for the Heat,
Mass and Momentum Transfer in Three-Dimensional Parabolic Flows,
Int. J. Heat and Mass Transfer, Vol. 15, pp. 1787-1805, 1972.
19. Peyret, R. and Taylor, T.D., Computational Methods for Fluid Flow,
Springer Verlag, 1983.
20. Robichaux,J., Tafti, D.K. and Vanka, S.P., Large Eddy Simulation of
Turbulence on the CM-2, Numerical Heat Transfer, Part-B, Vol. 21, pp.
367-388, 1992.
 Solution of Navier-Stokes.. 6.31

21. Thomas, L.H., Elliptic Problems in Linear Difference Equation Over A

Network, Waston Sci. Compt. Lab. Report, Columbia University, New
York, 1949.
22. Van Doormal, J.P. and G. D. Raithby, G.D., Enhancement of the SIM-
PLE Methods for Predicting Incompressible Fluid Flows, Numerical Heat
Transfer, Vol. 7, pp. 147-163, 1984.
23. Vanka, S.P., Chen, B. C. J. and Sha, W.T., A Semi-Implicit Calcula-
tion Procedure for Flows Described in Body-Fitted Coordinate System,
Numerical Heat Transfer, Vol. 3, pp. 1-19, 1980.
24. Viecelli, J.A., Computing Method for Incompressible flows Bounded by
Moving Walls, J. Comput. Phys., Vol. 8, pp. 119-143, 1971.

Problems
1. Develop the MAC and SIMPLE algorithms described in this chapter for
a steady one dimensional flow field u(y), driven by a constant (but as yet
undetermined) pressure gradient. For definiteness, flow in a parallel plate
channel can be considered.
2. For a two-dimensional flow, show that the MAC type iterative correc-
tion of pressure and velocity field through the implicit continuity equa-
tion is equivalent to the solution of Poisson equation for pressure. How
does the procedure avoid the need of directly applying pressure-boundary-
conditions?
3. Compare the solutions obtained by the stream function-vorticity method
with those presented in this chapter for the following problems:
(a) Developing flow in a parallel plates channel;
(b) Flow past a square cylinder (Re < 200);
(c) Flow past a periodic array of square cylinders (Re < 50).

4. Apply MAC method to solve the shear driven cavity flow problem shown
in Figure 6.13. Take a grid size of 51 × 51 and solve the flow equations
for Re = 400.
5. The Navier-Stokes equations in three-dimensional spherical-polar coordi-
nate system in θ direction may be written as
!
∂Vθ ∂Vθ Vθ ∂Vθ Vφ ∂Vθ Vr Vθ Vφ2 cotθ
ρ + Vr + + · + −
∂t ∂r r ∂θ rsinθ ∂φ r r
 2
∂p ∂ Vθ 2 ∂Vθ 1 ∂ 2 Vθ cotθ ∂Vθ
= − +µ + + + 2
r∂θ ∂r2 r ∂r r2 ∂θ2 r ∂θ
2

1 ∂ Vθ 2 ∂Vr Vθ 2cotθ cscθ ∂Vφ
+ · + 2 − 2 2 −
r2 sin2 θ ∂φ2 r ∂θ r sin θ r2 ∂φ