Chapter 2 :

Conservation Laws of Fluid Motion and

Boundary Conditions

Ibrahim Sezai
Department of Mechanical Engineering
Eastern Mediterranean University

Fall 2009-2010

Governing Equations of Fluid Flow and Heat Transfer

The governing equations of fluid flow represent

mathematical statements of the conservation laws of
physics.

The mass of fluid is conserved

The rate of change of momentum equals the sum of the
forces on a fluid particle (Newton’s second law)
The rate of change of energy is equal to the sum of the
rate of heat addition to and the rate of work done on a
fluid particle (first law of thermodynamics).
ME555 : Computational Fluid Dynamics 2 I. Sezai – Eastern Mediterranean University

1
The six faces are labelled
N, S, E, W, T, B

The center of the element is

located at position (x, y, z)

ρ = ρ (x, y, z, t) p = p (x, y, z, t)
T = T (x, y, z, t) u = u (x, y, z, t) Fig. 2-1 Fluid element for
conservation laws.
Fluid properties at faces are approximated by means of the two terms
of the Taylor series
The pressure at the W and E faces, can be expressed as
∂p 1 ∂p 1
p− δx and p+ δx
∂x 2 ∂x 2
ME555 : Computational Fluid Dynamics 3 I. Sezai – Eastern Mediterranean University

Mass Conservation in Three Dimensions

⎛ Rate of increase ⎞ ⎛ Net rate of flow ⎞

⎜ of mass in ⎟ = ⎜ of mass into ⎟
⎜ fluid element ⎟ ⎜ fluid element ⎟
⎝ ⎠ ⎝ ⎠

∂ ∂ρ ⎛ ∂( ρu) 1 ⎞ ⎛ ∂( ρu) 1 ⎞
( ρδ xδ yδ z ) = (δ xδ yδ z ) = ⎜ ρ u − δ x ⎟ δ yδ z − ⎜ ρ u + δ x ⎟ δ yδ z
∂t ∂t ⎝ ∂x 2 ⎠ ⎝ ∂x 2 ⎠
⎛ ∂( ρ v) 1 ⎞ ⎛ ∂ ( ρ v) 1 ⎞
+ ⎜ ρv − δ y ⎟ δ xδ z − ⎜ ρ v + δ y ⎟ δ xδ z
⎝ ∂y 2 ⎠ ⎝ ∂y 2 ⎠
⎛ ∂ ( ρ w) 1 ⎞ ⎛ ∂ ( ρ w) 1 ⎞
+ ⎜ ρw − δ z ⎟ δ xδ y − ⎜ ρ w + δ z ⎟ δ xδ y
⎝ ∂z 2 ⎠ ⎝ ∂z 2 ⎠

∂ρ ∂ ( ρ u ) ∂ ( ρ v) ∂ ( ρ w) Net flow of mass out of the control volume

+ + + =0
∂t ∂x ∂y ∂z
∂ρ (2-4)
Or in more compact vector notation Æ + div( ρ u) = 0
∂t ∂ ( ρ u ) ∂ ( ρ v) ∂ ( ρ w)
For an incompressible fluid ρ = const Æ div( ρ u) = 0 or ∂x
+
∂y
+
∂z
=0

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

2
Rates of change following a fluid particle and for a fluid element

The total or substantial derivative of φ with respect to time following a

fluid particle is

Dφ ∂φ ∂φ dx ∂φ dy ∂φ dz
= + + +
Dt ∂t ∂x dt ∂y dt ∂z dt
A fluid particle follows the flow, so
dx / dt = u
dy / dt = v
dz / dt = w

Hence the substantive derivative of φ is given by

Dφ ∂φ ∂φ ∂φ ∂φ ∂φ
= +u +v +w = + u ⋅ gradφ
Dt ∂t ∂x ∂y ∂z ∂t

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

Dφ/Dt defines the rate of change of property φ per unit mass.

The rate of change of property φ per unit volume for a fluid particle is
ρDφ /Dt, hence
Dφ ⎛ ∂φ ⎞
ρ = ρ⎜ + u ⋅ gradφ ⎟
Dt ⎝ ∂t ⎠
ρ = mass per unit volume.
Lhs of the mass conservation equation (2-4) is
∂ρ
+ div( ρ u)
∂t
The generalization of these terms for an arbitrary conserved property is
∂ ( ρφ ) (2-9)
+ div( ρφ u)
∂t
⎛ Net rate of flow of φ ⎞
⎛ Rate of increase ⎞ ⎜ ⎟
⎜ ⎟ + ⎜ out of fluid element ⎟
⎝ of φ per unit volume ⎠ ⎜ ⎟
⎝ per unit volume ⎠
ME555 : Computational Fluid Dynamics 6 I. Sezai – Eastern Mediterranean University

3
Rewriting eq. (2-9) = 0 (due to conservation
of mass)
 
∂ ( ρφ ) ⎡ ∂φ ⎤ ⎡ ∂ρ ⎤ Dφ
+ div( ρφ u) = ρ ⎢ + u ⋅ gradφ ⎥ + φ ⎢ + div( ρ u) ⎥ = ρ
∂t ⎣ ∂t ⎦ ⎣ ∂t ⎦ Dt

⎛ Rate of increase ⎞ ⎛ Net rate of flow of φ ⎞ ⎛ Rate of increase ⎞

⎜ ⎟ ⎜ ⎟ ⎜ ⎟
⎜ of φ of ⎟ + ⎜ out of ⎟ = ⎜ of φ for a ⎟
⎜ fluid element ⎟ ⎜ fluid element ⎟ ⎜ fluid particle ⎟
⎝ ⎠ ⎝ ⎠ ⎝ ⎠

x-momentum u Du ∂( ρu)
ρ + div ( ρ uu )
Dt ∂t
y-momentum v Dv ∂ ( ρ v)
ρ + div ( ρ vu )
Dt ∂t
z-momentum w Dw ∂ ( ρ w)
ρ + div ( ρ wu )
Dt ∂t
energy E DE ∂( ρ E )
ρ + div ( ρ Eu )
Dt ∂t
ME555 : Computational Fluid Dynamics 7 I. Sezai – Eastern Mediterranean University

The rates of increase of x-, y-, and z-momentum per unit volume are
Du Dv Dw
ρ ρ ρ
Dt Dt Dt
We distinguish two types of forces on fluid particles:
• surface forces - pressure forces
- viscous forces
• body forces - gravity forces
- centrifugal forces source terms
- Coriolis forces
- electromagnetic force
The pressure, a normal stress, is denoted by p.
Viscous stresses are denoted by τ.
ME555 : Computational Fluid Dynamics 8 I. Sezai – Eastern Mediterranean University

4
Fig. 2-3 Stress components on three
faces of fluid element.
The suffices i and j in τij indicate that
the stress component acts in the j-
direction on a surface normal to the i-
direction.

First we consider the x-

components of the forces
due to pressure p and stress
components τxx, τyx and τzx
shown in Fig. 2-4.

Fig. 2-4 Stress components in

the x-direction.
ME555 : Computational Fluid Dynamics 9 I. Sezai – Eastern Mediterranean University

The net force in the x-direction is the sum of

the force components acting in that direction
on the fluid element.
On the pair of faces (E, W) we have

⎡⎛ ∂p 1 ⎞ ⎛ ∂τ xx 1 ⎞ ⎤
⎢⎜ p − ∂x 2 δ x ⎟ − ⎜τ xx − ∂x 2 δ x ⎟ ⎥ δ yδ z
⎣⎝ ⎠ ⎝ ⎠⎦
⎡ ⎛ ∂p 1 ⎞ ⎛ ∂τ xx 1 ⎞ ⎤
+ ⎢− ⎜ p + δ x ⎟ + τ xx + δ x δ yδ z
⎣ ⎝ ∂x 2 ⎠ ⎜⎝ ∂x 2 ⎟⎠ ⎥⎦
⎛ ∂p ∂τ ⎞ (2-12a)
= ⎜ − + xx ⎟ δ xδ yδ z
⎝ ∂x ∂x ⎠
The net force in the x-direction on the pair of faces (N, S) is
⎛ ∂τ ⎞ ⎛ ∂τ ⎞ ∂τ
− ⎜ τ yx − yx 12 δ y ⎟ δ xδ z + ⎜ τ yx + yx 12 δ y ⎟ δ xδ z = yx δ xδ yδ z (2-12b)
⎝ ∂y ⎠ ⎝ ∂y ⎠ ∂y
The net force in the x-direction on the pair of faces (T, B) is
⎛ ∂τ ⎞ ⎛ ∂τ ⎞ ∂τ
− ⎜ τ zx − zx 12 δ z ⎟ δ xδ y + ⎜ τ zx + zx 12 δ z ⎟ δ xδ y = zx δ xδ yδ z (2-12c)
⎝ ∂z ⎠ ⎝ ∂z ⎠ ∂z
ME555 : Computational Fluid Dynamics 10 I. Sezai – Eastern Mediterranean University

5
 The total force per unit volume on the fluid due to these surface
stresses is equal to the sum of (2-12a), (2-12b), (2-12c) divided by the
volume δxδyδz:
∂ (− p + τ xx ) ∂τ yx ∂τ zx (2.13)
+ +
∂x ∂y ∂z
To find x-component of the momentum equation:
⎛ Rate of change of ⎞ ⎛ Total force in x-direction ⎞ ⎛ Total force in x -direction ⎞
⎜ x-momentum of ⎟ = ⎜ on the element due to ⎟ + ⎜ on the element due to ⎟
⎜ fluid particle ⎟ ⎜ surface stresses ⎟ ⎜ body forces ⎟
⎝ 
⎠ ⎝  
⎠ ⎝  
⎠
Eqn.(2.11) Eqn.(2.13) S Mx

Du ∂ (− p + τ xx ) ∂τ yx ∂τ zx (2.14a)
ρ = + + + S Mx
Dt ∂x ∂y ∂z
SMx = Body force on the element per unit volume in x-direction
SMz = –ρg (body force due to gravity per unit volume)
ME555 : Computational Fluid Dynamics 11 I. Sezai – Eastern Mediterranean University

Similarly, y-component of the momentum equation is

Dv ∂τ xy ∂ (− p + τ yy ) ∂τ zy
ρ = + + + S My (2.14b)
Dt ∂x ∂y ∂z
and, z-component of the momentum equation is
Dw ∂τ xz ∂τ yz ∂ (− p + τ zz )
ρ = + + + S Mz (2.14c)
Dt ∂x ∂y ∂z

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

6
 Energy Equation in Three Dimensions

The energy equation is derived from the first law of thermodynamics

which states that

( )
Rate of increase ⎛ Net rate of ⎞ ⎛ Net rate of work ⎞
= ⎜ heat added to ⎟ + ⎜ done on ⎟
of energy of fluid particle ⎜ fluid particle ⎟ ⎜ fluid particle ⎟
  
⎝ ⎠ ⎝ ⎠
DE
ρ
Dt

Work Done by Surface Forces = Fsurface forces × V

V = velocity component in the direction of the force.

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

The surface forces given by (2.12a-c) all act in the x-direction.

The net rate of work done by these forces acting in x-direction is

⎡ ∂ [u (− p + τ xx ) ] ∂ (uτ yx ) ∂ (uτ zx ) ⎤
⎢ + + ⎥ δ xδ yδ z (2.16a)
⎣ ∂x ∂y ∂z ⎦

Similarly, work done by surface stresses in y and z-direction are

⎡ ∂ (vτ xy ) ∂ ⎡v(− p + τ yy ) ⎤ ∂ (vτ zy ) ⎤

⎢ + ⎣ ⎦+ ⎥ δ xδ yδ z (2.16b)
⎢⎣ ∂x ∂y ∂z ⎥
⎦
⎡ ∂ ( wτ xz ) ∂ ( wτ yz ) ∂ [ w(− p + τ zz ) ] ⎤
⎢ + + ⎥ δ xδ yδ z (2.16c)
⎣ ∂x ∂y ∂z ⎦

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

7
Summing (2.16a-c) yields the total rate of work done on the fluid
particle by surface stresses:

∂ (uτ xx ) ∂ (uτ yx ) ∂ (uτ zx ) ∂ (vτ xy ) ∂ (vτ yy ) ∂ (vτ zy )

[ −div( pu)] + + + + + +
∂x ∂y ∂z ∂x ∂y ∂z
∂ ( wτ xz ) ∂ ( wτ yz ) ∂ ( wτ zz )
+ + +
∂x ∂y ∂z

where
∂ (up ) ∂ (vp ) ∂ ( wp )
− div( pu) = − − −
∂x ∂y ∂z

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

Energy Flux due to Heat Conduction

The heat flux vector has three
components qx, qy, qz

The net rate of heat transfer to the CV due to heat flow in x-direction is

⎡⎛ ∂qx 1 ⎞ ⎛ ∂qx 1 ⎞ ⎤ ∂qx

⎢⎜ qx − ∂x 2 δ x ⎟ − ⎜ qx + ∂x 2 δ x ⎟ ⎥ δ yδ z = − ∂x δ xδ yδ z (2.18b-c)
⎣⎝ ⎠ ⎝ ⎠⎦

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

8
Similarly, the net rates of heat transfer to the fluid due to heat flows in
the y- and z-direction are
∂q ∂q (2.18b-c)
− y δ xδ yδ z and − z δ xδ yδ z
∂y ∂z
The net rate of heat added to CV per unit volume is the sum of (2.18a-
c) divided by δxδyδz
∂q ∂q ∂q
− x − y − z = −div q (2.19)
∂x ∂y ∂z
∂T ∂T ∂T
qx = −k q y = −k qz = − k
∂x ∂y ∂z
This can be written in vector form as
q = −k grad T
Combining (2.19) and (2.20) yields the rate of heat addition to the
CV due to heat conduction
−div q = div(k grad T )
ME555 : Computational Fluid Dynamics 17 I. Sezai – Eastern Mediterranean University

Energy Equation
sum of the net rate of work done on the CV
by surface stresses (2.17)

⎡ ∂ (uτ xx ) ∂ (uτ yx ) ∂ (uτ zx ) ∂ (vτ xy ) ⎤
⎢ −div( pu) + + + + ⎥
DE ∂ x ∂y ∂z ∂x ⎥
ρ =⎢ (2.22)
Dt ⎢ + ∂ (vτ yy ) + ∂ (vτ zy ) + ∂ ( wτ xz ) + ∂ ( wτ yz ) + ∂ ( wτ zz ) ⎥
⎢⎣ ∂y ∂z ∂x ∂y ∂z ⎥⎦

+ div(k grad T ) + SE
  
N
net rate of heat addition rate of increase of energy
to the fluid (2.21) due to sources
E = i + 12 (u 2 + v 2 + w2 )
 

kinetic energy
i = internal (thermal) energy
SE = source of energy per unit volume per unit time (i.e. effects of
potential energy changes)
ME555 : Computational Fluid Dynamics 18 I. Sezai – Eastern Mediterranean University

9
Multiplying
the x-momentum equation (2.14a) by u
the y-momentum equation (2.14a) by v
the z-momentum equation (2.14a) by w

and adding the results together

D ⎡⎣ 12 (u 2 + v 2 + w2 ⎤⎦ ⎛ ∂τ ∂τ yx ∂τ zx ⎞
ρ = −u ⋅ grad p + u ⎜ xx + + ⎟
Dt ⎝ ∂x ∂y ∂z ⎠
⎛ ∂τ ∂τ ∂τ ⎞
+ v ⎜ xy + yy + zy ⎟
⎝ ∂x ∂y ∂z ⎠
⎛ ∂τ ∂τ ∂τ ⎞
+ w ⎜ xz + yz + zz ⎟ + u ⋅ S M (2.23)
⎝ ∂x ∂y ∂z ⎠

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

Subtracting (2.23) from (2.22)

Di ∂u ∂u
ρ = − p div u + div( k grad T ) + τ xx + τ yx
Dt ∂x ∂y
∂u ∂v ∂v ∂v ∂w
+ τ zx + τ xy + τ yy + τ zy + τ xz
∂z ∂x ∂y ∂z ∂x
∂w ∂w
+ τ yz + τ zz + Si
∂y ∂z
where Si =SE – u.SM
For an incompressible fluid Æ i = cT and div u = 0 (c = specific heat)
DT ∂u ∂u ∂u
ρc = div(k grad T ) + τ xx + τ yx + τ zx
Dt ∂x ∂y ∂z
∂v ∂v ∂v ∂w
+ τ xy + τ yy + τ zy + τ xz
∂x ∂y ∂z ∂x
∂w ∂w
+ τ yz + τ zz + Si
∂y ∂z
ME555 : Computational Fluid Dynamics 20 I. Sezai – Eastern Mediterranean University

10
 h =i+ p/ρ and ho = h + 12 (u 2 + v 2 + w2 )

Specific enthalphy Specific total enthalphy

Combining these two definitions with the one for specific energy E
ho = i + p / ρ + 12 (u 2 + v 2 + w2 ) = E + p / ρ (2.26)
Substituting of (2.26) into (2.22) yields the (total) enthalphy equation
∂ ( ρ ho )
+ div( ρ ho u) = div(k grad T )
∂t
∂p ∂ (uτ xx ) ∂ (uτ yx ) ∂ (uτ zx )
+ + + +
∂t ∂x ∂y ∂z
∂ (vτ xy ) ∂ (vτ yy ) ∂ (vτ zy )
+ + +
∂x ∂y ∂z
∂ ( wτ xz ) ∂ ( wτ yz ) ∂ ( wτ zz )
+ + + + Sh (2.27)
∂x ∂y ∂z
ME555 : Computational Fluid Dynamics 21 I. Sezai – Eastern Mediterranean University

Equations of State
Thermodynamic variables: ρ, p, i and T.
Relationships between the thermodynamic variables can be obtained
through the assumption of thermodynamic equilibrium.
Equations of state for pressure p and specific internal energy i:
p = p(ρ, T) and i = i(ρ, T)
For a perfect gas equations of state are
p = ρRT and i = CvT
In the flow of compressible fluids the equations of state provide the
linkage between the energy equation and mass conservation and
momentum equations.
Liquids and gases flowing at low speeds behave as incompressible
fluids.
Without density variations there is no linkage between the energy
equation and the mass conservation and momentum equations.

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

11
 Navier-Stokes Equations for a Newtonian Fluid
We need a suitable model for the viscous stresses τij.
Viscous stresses can be expressed as functions of the
local deformation rate (or strain rate).
In 3D flows the local rate of deformation is
composed of the linear deformation rate and the
volumetric deformation rate.
All gases and many liquids are isotropic.
The rate of linear deformation of a fluid element has
nine components in 3D, six of which are
independent in isotropic fluids.
They are denoted by the symbol eij.

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

There are three linear elongating deformation components:

∂u ∂v ∂w
exx = eyy = ezz =
∂x ∂y ∂z
There are also shearing linear deformation components:

⎛ ∂u ∂v ⎞ ⎛ ∂u ∂w ⎞
exy = eyx = 12 ⎜ + ⎟ exz = ezx = 12 ⎜ + ⎟
⎝ ∂y ∂x ⎠ ⎝ ∂z ∂x ⎠
⎛ ∂v ∂w ⎞
eyz = ezy = 12 ⎜ + ⎟
⎝ ∂z ∂y ⎠

The volumetric deformation is given by

∂u ∂v ∂w
+ + = div u
∂x ∂y ∂z
ME555 : Computational Fluid Dynamics 24 I. Sezai – Eastern Mediterranean University

12
 In a Newtonian fluid the viscous stresses are
proportional to the rates of deformation.
The 3D form of Newton’s law of viscosity for
compressible flows involves two constants of
proportionality:
- The (first) dynamic viscosity, μ, to relate stresses
to linear deformations,
- The second viscosity, λ, to relate stresses to the
volumetric deformation.

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

The nine viscous stress components, of which six are independent, are

∂u ∂v ∂w
τ xx = 2 μ + λ div u τ yy = 2μ + λ div u τ zz = 2μ + λ div u
∂x ∂y ∂z
⎛ ∂u ∂v ⎞ ⎛ ∂u ∂w ⎞
τ xy = τ yx = μ ⎜ + ⎟ τ xz = τ zx = μ ⎜ + ⎟
⎝ ∂y ∂x ⎠ ⎝ ∂z ∂x ⎠
(2.31)
⎛ ∂v ∂w ⎞
τ yz = τ zy = μ ⎜ + ⎟
⎝ ∂z ∂y ⎠

Not much is known about the second viscosity λ, because its effect is
small.
For gases a good working approximation is λ = –⅔μ
Liquids are incompressible so the mass conservation equation is
div u = 0
ME555 : Computational Fluid Dynamics 26 I. Sezai – Eastern Mediterranean University

13
Substitution of the above shear stresses (2.31) into (2.14a-c) yields the
Navier-Stokes equations
Du ∂p ∂ ⎡ ∂u ⎤ ∂ ⎡ ∂u ∂v ⎤
ρ = − + ⎢ 2μ + λ div u ⎥ + ⎢ μ + ⎥ (2.32a)
Dt ∂x ∂x ⎣ ∂x ⎦ ∂y ⎣ ∂y ∂x ⎦
∂ ⎡ ∂v ∂w ⎤
+ μ + + S Mx
∂z ⎣⎢ ∂z ∂x ⎦⎥

Dv ∂p ∂ ⎡ ∂u ∂v ⎤ ∂ ⎡ ∂v ⎤ (2.32b)
ρ = − + ⎢μ + + 2μ + λ div u ⎥
Dt ∂y ∂x ⎣ ∂y ∂x ⎥⎦ ∂y ⎢⎣ ∂y ⎦
∂ ⎡ ∂v ∂w ⎤
+ μ + + S My
∂z ⎢⎣ ∂z ∂y ⎥⎦
Dw ∂p ∂ ⎡ ∂u ∂w ⎤ ∂ ⎡ ∂v ∂w ⎤
ρ = − + ⎢μ + + μ + (2.32c)
Dt ∂z ∂x ⎣ ∂z ∂x ⎥⎦ ∂y ⎢⎣ ∂z ∂y ⎥⎦
∂ ⎡ ∂w ⎤
+ 2μ + λ div u ⎥ + S Mz
∂z ⎢⎣ ∂z ⎦
ME555 : Computational Fluid Dynamics 27 I. Sezai – Eastern Mediterranean University

Often it is useful to rearrange the viscous stress terms as follows:

∂ ⎡ ∂u ⎤ ∂ ⎡ ∂u ∂v ⎤ ∂ ⎡ ∂v ∂w ⎤
2μ + λ div u ⎥ + ⎢ μ + ⎥+ μ +
∂x ⎢⎣ ∂x ⎦ ∂y ⎣ ∂y ∂x ⎦ ∂z ⎢⎣ ∂z ∂x ⎥⎦
∂ ⎛ ∂u ⎞ ∂ ⎛ ∂u ⎞ ∂ ⎛ ∂u ⎞ div( μ grad u )
= ⎜μ ⎟+ ⎜μ ⎟+ ⎜μ ⎟
∂x ⎝ ∂x ⎠ ∂y ⎝ ∂y ⎠ ∂z ⎝ ∂z ⎠
⎡ ∂ ⎛ ∂u ⎞ ∂ ⎛ ∂v ⎞ ∂ ⎛ ∂w ⎞ ⎤ ∂
+ ⎢ ⎜μ ⎟+ ⎜μ ⎟+ ⎜μ ⎟ ⎥ + (λ div u) S Mx
⎣ ∂x ⎝ ∂x ⎠ ∂y ⎝ ∂x ⎠ ∂z ⎝ ∂x ⎠ ⎦ ∂x
= div( μ grad u ) + S Mx
The viscous stresses in the y- and z-momentum equations can be re-cast in a similar
manner.
To simplify the momentum equations:
‘hide’ the two smaller contributions to the viscous stress terms in the momentum
source.
Defining a new source by
SM = SM + s M
(2.33)
ME555 : Computational Fluid Dynamics 28 I. Sezai – Eastern Mediterranean University

14
the Navier-Stokes equations can be written in the most useful form for
the development of the finite volume method:
Du ∂p
ρ = − + div( μ grad u ) + S Mx
Dt ∂x (2.34a)

Dv ∂p
ρ = − + div( μ grad v) + S My (2.34b)
Dt ∂y

Dw ∂p
ρ = − + div( μ grad w) + S Mz (2.34c)
Dt ∂z

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

If we use the Newtonian model for viscous stresses in the internal

energy equation (2.24) we obtain
Di (2.35)
ρ = − p div u + div( k grad T ) + Φ + Si
Dt

The dissipation function Φ is

⎧ ⎡⎛ ∂u ⎞2 ⎛ ∂v ⎞ 2 ⎛ ∂w ⎞ 2 ⎤ ⎛ ∂u ∂v ⎞ 2 ⎫
⎪ 2 ⎢⎜ ⎟ + ⎜ ⎟ + ⎜ ⎟ ⎥+⎜ + ⎟ ⎪
⎪ ⎢⎣⎝ ∂x ⎠ ⎝ ∂y ⎠ ⎝ ∂z ⎠ ⎥⎦ ⎝ ∂y ∂x ⎠ ⎪
Φ=μ⎨ ⎬ + λ (div u)
2
(2.36)
2
⎪ ⎛ ∂u ∂w ⎞ ⎛ ∂v ∂w ⎞
2
⎪
⎪+ ⎜ + ⎟ +⎜ + ⎟ ⎪
⎩ ⎝ ∂z ∂x ⎠ ⎝ ∂z ∂y ⎠ ⎭

The dissipation function represents a source of internal energy due

to deformation work on the fluid particle.

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

15
Conservative Form of the Governing Equations of Fluid Flow

The conservative or divergence form of the time-dependent 3-D flow

and energy equations of a compressible Newtonian fluid:
∂ρ (2.4)
Mass + div( ρ u) = 0
∂t
∂( ρu) ∂p
x-momentum + div( ρ uu) = − + div( μ grad u ) + S Mx (2.37a)
∂t ∂x
∂ ( ρ v) ∂p
y -momentum + div( ρ vu) = − + div( μ grad v) + S My (2.37b)
∂t ∂y
∂ ( ρ w) ∂p
z -momentum + div( ρ wu) = − + div( μ grad w) + S Mz (2.37c)
∂t ∂z
∂( ρi)
internal energy + div( ρ iu ) = − p div u + div(k grad T ) + Φ + Si (2.38c)
∂t
equations of state p = p ( ρ , T ) and i = i ( ρ , T ) (2.28)
e.g. for a perfect gas:
p = ρ RT and i = CvT (2.29)
Table 2.1
A system of seven equations with seven unknowns Æ this system is mathematically
closed.
ME555 : Computational Fluid Dynamics 31 I. Sezai – Eastern Mediterranean University

Differential and Integral Forms of the General Transport Equations

Equations in Table 2.1 can usefully written in the following form:

∂ ( ρφ )
+ div( ρφ u) = div( μ grad φ ) + Sφ
∂t (2.39)
⎛ Rate of increase ⎞ ⎛ Net rate of flow ⎞ ⎛ Rate of increase ⎞ ⎛ Rate of increase ⎞
⎜ of φ of fluid ⎟ + ⎜ of φ out of ⎟ = ⎜ of φ due to ⎟ + ⎜ of φ due to ⎟
⎜ element ⎟ ⎜ fluid element ⎟ ⎜ diffusion ⎟ ⎜ sources ⎟
⎝ ⎠ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠

Rate of change term convective term diffusive term source term

Equation (2.39) is used as the starting point in finite volume method.

⎧⎪φ = 1, u , v, w, i ⎫⎪
By setting ⎨Γ = 0, μ , k ⎬ → we obtain equations in Table 2.1.
S
⎩⎪ φ = 0, ( S Mx − ∂p / ∂x ),..., ⎭⎪
ME555 : Computational Fluid Dynamics 32 I. Sezai – Eastern Mediterranean University

16
In the finite volume method (2.39) is integrated over 3-D control
volume yielding

∂ ( ρφ )
∫
CV
∂t
dV + ∫ div( ρφ u)dV = ∫ div(Γ grad φ )dV + ∫ Sφ dV
CV CV CV
(2.40)

For a vector a Gauss’ divergence theorem states

(2.41)
CV
∫ div adV = ∫ n ⋅ adA
A

Applying Gauss’ divergence theorem, equation (2.40) can be written

as
∂⎛ ⎞
⎜ ∫ ρφ dV ⎟ + ∫ n ⋅ ( ρφ u)dA = ∫ n ⋅ (Γ grad φ )dA + ∫ Sφ dV (2.42)
∂t ⎝ CV ⎠ A A CV

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

Equation (2.42) can be expressed as follows:

⎛ Net rate of ⎞ ⎛ Rate of increase ⎞
( Rate of + ⎜
)
decrease of φ due to ⎟ = ⎜ of φ due to ⎟ + Net rate of
increase of φ ⎜ convection across ⎟ ⎜ diffusion across ⎟ creation of φ
⎜ the boundaries ⎟ ⎜ the boundaries ⎟
( )
⎝ ⎠ ⎝ ⎠
In steady state problems the rate of change term of (2.42) is equal to
zero.

∫ n ⋅ ( ρφu)dA =∫ n ⋅ (Γ grad φ )dA + ∫ Sφ dV

A A CV
(2.43)

Integrating (2.42) with respect to time

∂⎛ ⎞
∫ ∂t ⎜⎝ ∫ ρφ dV ⎟⎠dt + ∫ ∫ n ⋅ ( ρφu)dAdt
Δt CV Δt A
= ∫ ∫ n ⋅ (Γ grad φ )dAdt + ∫ ∫ Sφ dVdt (2.44)
Δt A Δt CV

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

17
 Auxiliary Conditions for Viscous Fluid Flow Equations
Table 2-5 Boundary conditions for compressible viscous flow.
Initial conditions for unsteady flows:
• Everywhere in the solution region ρ, u and T must be given at time t = 0
Boundary conditions for unsteady and steady flows:
• On solid walls u = uw (no-slip condition)
T = Tw (fixed temperature) or k∂T/∂n = –qw (fixed heat flux)
• On fluid boundaries inlet: ρ, u and T must be known as a function of position
outlet: –p +μ∂un/∂n =Fn and –p +μ∂ut/∂n =Ft (stress continuity)

Suffices: n → normal direction

to boundary
t → tangential direction
F → given surface stress
For incompressible viscous flows:
Table 2.5 is applicable, except that there are no conditions on the density ρ.
ME555 : Computational Fluid Dynamics 35 I. Sezai – Eastern Mediterranean University

Outflow boundaries:
• High Re flows far from solid objects in an external flow
• Fully developed flow out of a duct.

For these boundaries:

Pressure = specified
∂un/∂n = 0
∂T/∂n = 0

Sources and sinks of mass are placed on the inlet and outlet
boundaries to ensure the correct mass flow into and out of domain.

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

18
Boundary conditions for an internal flow problem.

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

Boundary conditions for an external flow problem.

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

19
Example to symmetry boundary conditions:

∂φ
=0
∂r

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

Example to cyclic boundary conditions:

Cyclic b.c.: φ1 = φ2

1
1 2

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

20