3.

6 Computational Fluid Dynamics

N i

E
W P

TL
S

x=L

Figure 3.4: Boundary Condition, at x = L, T = TL .

shortened control volume:

qe · Δy − qw · Δy + qn Δx1 − qs Δx1
(TL − TP ) (TP − TW ) (TN − TP )
= −k Δy + k Δy − k Δx1
(Δx/2) Δx Δy
(TP − TS )
+ k Δx1 = Q Δx1 Δy (3.11)
Δy
3
where Δx1 = Δx.
4
Note that TL has been substituted instead of TE in the above equation, so the
boundary condition is being directly applied. In this fashion, by adjusting the
control volume spacing and the placement of nodes, nodal equations can be
obtained at all nodes and these can be solved simultaneously by the matrix
inversion technique, line-by-line technique or point-by-point technique as dis-
cussed earlier. Having done the above exercise, we may like to look at a more
generalized description of the finite volume method.

3.2 A Generalized Approach

As it has been observed, the Finite Volume method uses an integral form of
the equations to be solved. The computational domain is divided into ele-
mentary volumes and the integration is performed within the these elementary
volumes. The method enables one to handle complex geometry without having
the equation written in curvilinear coordinates. The method also preserves the
 Introduction to Finite Volume Method 3.7

conservative property. The elementary control volumes are described by the

coordinates of the vertices of the quadrilaterals (for 2-D) or hexahedrals (for
3-D).

3.2.1 Equations with First Derivatives

Here the finite volume method will be illustrated for the general first-order
equation

∂E ∂F ∂G
+ + =0 (3.12)
∂t ∂x ∂y

where E, F and G represent the dependent variables. For example, for E =

ρ, F = ρu and G = ρv, Eqn. (3.12) is the two- dimensional continuity equation,
and for E = ρu, F = ρu2 , G = ρuv, it is the inviscid momentum equation in the
x-directions, and so on. In a similar manner, for x direction viscous momentum
equation,
 
2 ∂u ∂u
E = ρ u, F = ρu + p − μ ,G = ρ u v − μ (3.13)
∂x ∂y

Assuming the finite volume (quadrilateral) ABCD shown in Fig. 3.5 is the
representative of the control volume Ωp , we consider the area integral of (3.12)
over Ωp :
  
∂E ∂F ∂G
+ + dx dy = 0 (3.14)
ABCD ∂t ∂x ∂y

Recall the Green’s theorem on transformation between double integral and line
integral

   
∂g ∂f
− dx dy = f dx + gdy (3.15)
ABCD ∂x ∂y c

Applying Green’s theorem, (3.14) becomes

  
d
E dV + H · n ds = 0 (3.16)
dt ABCD

where H = (F, G) and n is the outward unit normal of segment ds (see Figure
3.5). For a segment ds = dxî + dy ĵ. On a counter-clockwise contour, the out-
dy î − dxĵ
ward unit normal n = , where ds = dx2 + dy 2 . For the continuity
ds
equation, F = ρu and G = ρv. In Cartesian coordinates,

H · n ds = F dy − G dx (3.17)
3.8 Computational Fluid Dynamics

Equation (3.16) is a statement of conservation. For the particular choice,

E = ρ, F = ρ u, G = ρ v, in Eqn. (3.16) concides with an integral state-
ment of the conservation of mass. As mentioned, the finite volume method
involves discretization of the governing equation in integral form, in contrast to
the finite difference method, which is usually applied to the governing equation
in differential form.

An approximate evaluation of Eqn. (3.16) would be

d
(A Ei, j ) + (F Δy − G Δx)AB + (F Δy − G Δx)BC
dt
+(F Δy − G Δx)CD + (F Δy − G Δx)DA = 0 (3.18)

where A is the area of the quadrilateral ABCD in Fig. 3.5, and the average

Figure 3.5: Finite Volume Mesh System.

value of E over the quadrilateral is represented by Ei, j and the remaining terms
are approximations for the line integral over segments AB, BC, CD, and DA,
respectively . Further, writing

ΔyAB = yB − yA , ΔxAB = xB − xA and

FAB = 0.5 (Fi, j−1 + Fi, j ), GAB = 0.5 (Gi, j−1 + Gi, j )
 Introduction to Finite Volume Method 3.9

and similar expressions for ΔyBC , etc. if A is not a function of time, Eqn. (3.18)
becomes

A dEi,j /dt + 0.5 (Fi, j−1 + Fi, j ) ΔyAB − 0.5 (Gi, j−1 + Gi, j ) ΔxAB
+ 0.5 (Fi, j + Fi+1, j ) ΔyBC − 0.5 (Gi, j + Gi+1, j ) ΔxBC
+ 0.5 (Fi, j + Fi, j+1 ) ΔyCD − 0.5 (Gi, j + Gi, j+1 ) ΔxCD
+ 0.5 (Fi−1, j + Fi, j ) ΔyDA − 0.5 (Gi−1, j + Gi, j ) ΔxDA = 0 (3.19)

For the irregular grid-mesh (i, j), the finite volume Eqn. (3.19) provides a
discretisation in Cartesian coordinates without introducing generalised coordi-
nates. If the grid-mesh is uniform and coincides with lines of constant x and y,
Eqn. (3.19) becomes
d
Δx Δy Ei, j − 0.5 (Gi, j−1 + Gi, j ) Δx + 0.5 (Fi, j + Fi+1, j ) Δy
dt
+0.5 (Gi, j + Gi, j+1 ) Δx − 0.5 (Fi−1, j + Fi, j ) Δy = 0

or,
d Fi+1, j − Fi−1, j Gi, j+1 − Gi, j−1
Ei, j + + =0 (3.20)
dt 2 Δx 2 Δy
which is same as the central difference representation for the spatial terms in
(3.12). The finite volume method, which is extensively used for both incom-
pressible and compressible flows, has the advantages of conservative property.
Most importantly, it allows complex computational domains to be discretised
in a simple way.

3.2.2 Equations with Second Derivatives

In the previous section, the finite volume method was applied to Eqn. (3.12), in
which only first derivatives appeared, and yielded a relatively straightforward
discretization formula given by Eqn. (3.19). The finite volume method requires
to be modified, when the second derivatives are present in the governing partial
differential equation.
Let us consider Laplace’s equation
∂2 φ ∂2 φ
+ =0 (3.21)
∂x2 ∂y 2
The finite volume method follows the application of the subdomain method to
Eqn. (3.21). For the finite volume ABCD as shown in Fig. 3.6, one can write
the area integral of (3.21) as
  2 
∂ φ ∂2 φ
+ dx dy = H · n ds = 0 (3.22)
ABCD ∂x2 ∂y 2 ABCD

where Green’s theorem has again been used and now

 Introduction to Finite Volume Method 3.9

and similar expressions for ΔyBC , etc. if A is not a function of time, Eqn. (3.18)
becomes

A dEi,j /dt + 0.5 (Fi, j−1 + Fi, j ) ΔyAB − 0.5 (Gi, j−1 + Gi, j ) ΔxAB
+ 0.5 (Fi, j + Fi+1, j ) ΔyBC − 0.5 (Gi, j + Gi+1, j ) ΔxBC
+ 0.5 (Fi, j + Fi, j+1 ) ΔyCD − 0.5 (Gi, j + Gi, j+1 ) ΔxCD
+ 0.5 (Fi−1, j + Fi, j ) ΔyDA − 0.5 (Gi−1, j + Gi, j ) ΔxDA = 0 (3.19)

For the irregular grid-mesh (i, j), the finite volume Eqn. (3.19) provides a
discretisation in Cartesian coordinates without introducing generalised coordi-
nates. If the grid-mesh is uniform and coincides with lines of constant x and y,
Eqn. (3.19) becomes
d
Δx Δy Ei, j − 0.5 (Gi, j−1 + Gi, j ) Δx + 0.5 (Fi, j + Fi+1, j ) Δy
dt
+0.5 (Gi, j + Gi, j+1 ) Δx − 0.5 (Fi−1, j + Fi, j ) Δy = 0

or,
d Fi+1, j − Fi−1, j Gi, j+1 − Gi, j−1
Ei, j + + =0 (3.20)
dt 2 Δx 2 Δy
which is same as the central difference representation for the spatial terms in
(3.12). The finite volume method, which is extensively used for both incom-
pressible and compressible flows, has the advantages of conservative property.
Most importantly, it allows complex computational domains to be discretised
in a simple way.

3.2.2 Equations with Second Derivatives

In the previous section, the finite volume method was applied to Eqn. (3.12), in
which only first derivatives appeared, and yielded a relatively straightforward
discretization formula given by Eqn. (3.19). The finite volume method requires
to be modified, when the second derivatives are present in the governing partial
differential equation.
Let us consider Laplace’s equation
∂2 φ ∂2 φ
+ =0 (3.21)
∂x2 ∂y 2
The finite volume method follows the application of the subdomain method to
Eqn. (3.21). For the finite volume ABCD as shown in Fig. 3.6, one can write
the area integral of (3.21) as
  2 
∂ φ ∂2 φ
+ dx dy = H · n ds = 0 (3.22)
ABCD ∂x2 ∂y 2 ABCD

where Green’s theorem has again been used and now

3.10 Computational Fluid Dynamics

Figure 3.6: Finite Volume for a Curvilinear Grid.

∂φ ∂φ
H · n ds = dy − dx
∂x ∂y
Following the similar steps as in Sec. 3.2.1, the line integral in equation (3.22)
can be evaluated approximately over the segments AB, BC, CD, and DA, by
∂φ ∂φ
ΔyAB − ΔxAB
∂x i, j−1/2 ∂y i, j−1/2
∂φ ∂φ
+ ΔyBC − ΔxBC
∂x i+1/2, j ∂y i+1/2, j
∂φ ∂φ
+ ΔyCD − ΔxCD
∂x i, j+1/2 ∂y i, j+1/2
∂φ ∂φ
+ ΔyDA − ΔxDA
∂x i−1/2, j ∂y i−1/2, j
= 0 (3.23)
Various techniques for evaluating [∂φ /∂x]i, j−1/2 , etc. are discussed by Fletcher
(1988) and Peyret and Taylor (1983). Here [∂φ /∂x]i, j−1/2 is evaluated as the
    
mean value over the area B BC D AA B in Fig. 3.6. Thus we can write
(Fletcher, 1988).
     
∂φ 1 ∂φ 1
= dx dy = φdy (3.24)
∂x i,j−1/2 S A B  C  D  ∂x S A B  C  D 
 Introduction to Finite Volume Method 3.11

     
∂φ 1 ∂φ 1
= dx dy = − φdx (3.25)
∂y i,j−1/2 S A B  C  D  ∂y S A B  C  D 

where Green’s theorem has been used to obtain the final line integrals on the
right, and

φ dy ≈ φi,j−1 ΔyA B  + φB ΔyB  C  + φi,j ΔyC  D + φA ΔyD A
A B  C  D 

A similar expression can be obtained for


φ dx ≈ φi,j−1 ΔxA B  + φB ΔxB  C  + φi,j ΔxC  D + φA ΔxD A
A B  C  D 

If the mesh is not skewed, one can write

ΔyA B  ≈ − ΔyC  D ≈ ΔyAB and ΔyB  C  ≈ −ΔyD A ≈ Δyj−1, j

ΔxA B  ≈ − ΔxC  D ≈ ΔxAB and ΔxB  C  ≈ −ΔxD A ≈ Δxj−1, j
 
 
SAB ≡ SA B  C  D = A B  × B  C   (3.26)

or,
  
 
SAB =  ΔxA B  i + ΔyA B  
j × ΔxB  C  i + ΔyB  C  
j 
= ΔxA B  ΔyB  C  − ΔyA B  ΔxB  C  (3.27)

or,
SAB = ΔxAB Δyj−1, j − ΔyAB Δxj−1, j (3.28)
Therefore (3.24) becomes

∂φ ΔyAB (φi, j−1 − φi, j ) + Δyj−1, j (φB − φA )

= (3.29)
∂x i, j−1/2 SAB

and (3.25) becomes

∂φ − [ΔxAB (φi, j−1 − φi,j ) + Δxj−1, j (φB − φA )]

= (3.30)
∂y i, j−1/2 SAB

Now, the first line of (3.23) can be evaluated as

(Δx2AB + ΔyAB
2
)(φi, j−1 − φi, j ) (ΔxAB Δxj−1, j + ΔyAB Δyj−1, j )(φB − φA )
+
SAB SAB
(3.31)

In a similar manner, one can evaluate the following

  
∂φ 1 ∂φ 1
= dx dy = φ dy (3.32)
∂x i+1/2, j SA B  C  D ∂x SA B  C  D
3.12 Computational Fluid Dynamics

  
∂φ 1 ∂φ 1
= dx dy = − φ dx (3.33)
∂x i+1/2, j SA B  C  D ∂y SA B C  D
 


φ dy ≈ φB ΔyA B  + φi+1, j ΔyB  C  + φC ΔyC  D + φi, j ΔyD A
A B  C  D

and

φ dx ≈ φB ΔxA B  + φi+1, j ΔxB  C  + φC ΔxC  D + φi, j ΔxD A
A B  C  D

If the mesh is not distorted,

ΔyA B  = −ΔyC  D = Δyi, i+1 and ΔyB  C  = −ΔyD A = ΔyBC
ΔxA B  = −ΔxC  D = Δxi, i+1 and ΔxB  C  = −ΔxD A = ΔxBC

The equivalent expressions for [∂φ/∂x]i, j+1/2 , [∂φ/∂y]i, j+1/2 and other similar
terms in Eqn. (3.23), finally yields

MAB (φi, j−1 − φi, j ) + NAB (φB − φA ) + MBC (φi+1, j − φi, j ) + NBC (φC − φB )

+MCD (φi, − φi, j )+NCD (φD − φC )+MDA (φi−1, j − φi, j )+NDA (φA − φD ) = 0
j+1
(3.34)
where the geometrical parameters are
 
MAB = Δx2AB + ΔyAB 2
/ SAB , NAB = (ΔxAB Δxj−1, j + ΔyAB Δyj−1, j ) / SAB
 2 2

MBC = ΔxBC + ΔyBC / SBC , NBC = (ΔxBC Δxi+1, i + ΔyBC Δyi+1, i ) / SBC
 2 2

MCD = ΔxCD + ΔyCD / SCD , NCD = (ΔxCD Δxj+1, j + ΔyCD Δyj+1, j ) / SCD

and
 
MDA = Δx2DA + ΔyDA
2
/ SDA , NDA = (ΔxDA Δxi−1, i + ΔyDA Δyi−1, i ) / SDA

In Eqn. (3.34) φA , φB , φC , φD are evaluated as the average of the four sur-

rounding nodal values. Thus

φA = 0.25 (φi, j + φi−1, j + φi−1, j−1 + φi, j−1 )

Substitution into Eqn. (3.34) generates the following discretized form of Eqn.
(3.21):

0.25 (NCD − NDA ) φi−1, j+1 + [MCD + 0.25 (NBC − NDA ) ] φi, j+1
+ 0.25 (NBC − NCD ) φi+1, j+1 + [MDA + 0.25 (NCD − NAB ) ] φi−1, j
− (MAB + MBC + MCD + MDA ) φi, j + [MBC + 0.25 (NAB − NCD )] φi+1, j
+ 0.25 (NDA − NAB ) φi−1, j−1 + [MAB + 0.25 (NDA − NBC )] φi, j−1
+ 0.25 (NAB − NBC ) φi+1, j−1 =0 (3.35)
 Introduction to Finite Volume Method 3.13

Eqn.(3.35) may be rearranged as

ai, j φi, j = ai+1, j φi+1, j + ai−1, j φi−1, j + ai, j+1 φi, j+1 +
ai, j−1 φi, j−1 + ai+1, j+1 φi+1, j+1 + ai−1, j+1 φi−1, j+1+
ai+1, j−1 φi+1, j−1 + ai−1, j−1 φi−1, j−1 (3.36)
In Eqn. (3.35) the geometrical quantities such as MAB and NAB , etc. need to be
evaluated only once for a given grid and used for all subsequent calculations.The
coefficients, ai, j , ai+1, j , ai−1, j etc. in Eqn (3.36) can be determined by com-
paring Eqn (3.35) with Eqn (3.36). Equation (3.36) is solved conveniently using
a Successive Over-Relaxation (SOR) technique. Equation (3.36) is manipulated
to give an estimate of φn+1 i, j , thus

φ∗i, j = [ai+1, j φi+1, j + ai−1, j φi−1, j + ai, j+1 φi, j+1 +

ai, j−1 φi, j−1 + ai+1, j+1 φi+1, j+1 + ai−1, j+1 φi−1, j+1 +
ai+1, j−1 φi+1, j−1 + ai−1, j−1 φi−1, j−1 ] /ai, j (3.37)
and the improved better value is
φn+1 n ∗ n
i, j = φi, j + λ (φi, j − φi, j ) (3.38)
where λ is the over relaxation or under relaxation parameter.

One can also implement this in the following way

φ∗p = [aE φE + aW φW + aN φN + aS φS + aN E φN E +
aN W φN W + aSE φSE + aSW φSW +] /ap (3.39)
Finally,
φn+1
p = φnp + λ(φ∗p − φnp ) (3.40)
One attractive feature of the finite volume method is that Neumann (derivative)
boundary conditions can be handled as readily as Dirichlet boundary conditions
by direct substitution into Eqn. (3.23).
The discretised Equation (3.35) reduces to the centred finite difference scheme
on a uniform rectangular grid
φi−1, j − 2φi, j + φi+1, j φi, j−1 − 2φi, j + φi, j+1
+ =0 (3.41)
Δx2 Δy 2

References
1. Fletcher, C. A. J., Computational Techniques for Fluid Dynamics, Vol. 1
(Fundamentals and General Techniques) Springer Verlag, 1988.
2. Patankar, S. V., and Spalding, D. B., A Calculation Procedure for Heat
Mass and Momentum Transfer in Three Dimensional Parabolic Flows, Int.
J. Heat Mass Transfer, Vol. 15, pp. 1787-1805, 1972.