SCHOOL OF BIO AND CHEMICAL ENGINEERING

DEPARTMENT OF CHEMICAL ENGINEERING

UNIT – I – Computational Fluid Dynamics – SCHA1403

1
 I. Overview of CFD and Conservation Laws
Computational fluid dynamics deals with equations that represent a balance process for
mass, momentum, energy and chemical species. Many students have seen these
differential equations before in advanced courses in fluid dynamics or in convective heat
and mass transfer. These notes review the derivation process for the basic balance
equations with a view to their eventual conversion to numerical methods in
computational fluid dynamics.

One of the useful facts about the various balance equations for various quantities is their
similar form. This means that a numerical algorithm developed from one particular
balance quantity, say momentum, can be applied to another balance quantity, say energy.
With this goal in mind, the derivation of the basic differential equations provided in these
notes is aimed at demonstrating the similarity of the various quantities for which we
derive a balance equation.

We represent the various quantities such a mass, momentum, energy and mass of an
individual chemical species by the general symbol, . We then derive a “balance
equation” for , which is valid for any quantity. We then have to consider the particular
physics of the different quantities in the balance equations.

Continuity Equation

We use simple laws. Mass is conserved. Newton’s second law tells us that the rate of
change of momentum equals the applied force. The first law of thermodynamics says
that the rate of energy change equals the heat added plus the work done on a system.
Chemical balances tell us that changes in the mass of individual chemical species are
related to sources or sinks provided by chemical reactions. All of these simple laws,
applied to a flow system, can be stated in terms of a simple balance equation

STORAGE + OUTFLOW - INFLOW = SOURCE [1-1]

2
We apply this general balance equation (remembering that the symbol, can represent
mass, momentum or energy) to a variety of physical quantities such as mass, energy,
momentum, and chemical species.

The various quantities that we will be concerned with are shown in Table 1-1.

Table 1-1 – Examples of Quantities that Satisfy a Balance Equation

S.No mass x momentum y momentum z Energy Species
momentum
1 m Mu mv mw E + mV2/2 m(K)
2 1 U v w e + V2/2 W(K)
In this table, u, v, and w are the x, y and z velocity components, E is the total
thermodynamic internal energy, e is the thermodynamic internal energy per unit mass,
and m(K) is the mass of a chemical species, K, W(K) is the mass fraction of species K.
The other energy term, mV2/2, is the kinetic energy.

We will derive the general balance equation for a differential volume in three-
dimensional Cartesian coordinate space. The velocity components in the (x, y, z)
coordinate directions are denoted as (u, v, w).
The volume of the differential element is Δx Δy Δz. If the density of the fluid is denoted
by the symbol ρ, the mass, m, inside this volume is ρ Δx Δy Δz The rate at which the
quantity, Φ, is stored (or accumulates) in this volume over time is then given by the
following equation.

  (m )  ( xyz )  (  )
Storage = = = = xyz [1-2]
t t t t

The mass flow across a control volume face is the product of velocity, density and cross
sectional area. For example, the x direction flow, with a velocity component u, enters
and exits the control volume through an area given by the product the Δy Δz. Thus, the
mass flow in the x direction is given by the product (ρ u Δy Δz).

3
If the flowing mass has a particular per-unit-mass property, , then the flow of the total
quantity, Φ, across the boundary is given by the product of the per-unit-mass property, ,
times the mass flow rate. This means that the flow of across a control volume face is
the product of four factors: density, ρ, the per-unit-mass value of the quantity considered,
, the cross sectional area of the face, and the velocity normal to the face. Note that this
product has the same dimensions as the storage term, the dimensions of divided by
dimensions of time.

For the coordinate system shown, we define inflows as occurring at a particular face
designated as x, y or z. The outflows then occur at the face in the increasing coordinate
direction. These faces are denoted as x+ Δx, y + Δy, and z + Δz. The inflow in the x
direction would be given by the following expression: ρ u |x Δy Δz. Similar terms
apply in the y and z directions. Summing the terms for each coordinate direction gives
the total inflow.

Inflow =  u  x y z +  v  y x z +  w z y x [1-3]

The expression for the outflow is similar. The only difference is the subscript indicating
that the outflow occurs at a different face of the control volume.

Outflow =  u  x + x y z +  v  y + y x z +  w z + z y x [1-4]

The source term, Sφ, must be defined for each individual physical quantity that satisfies a
balance equation. For the mass balance, the source term is zero. In the balance equation
for the ith chemical species, the source term is the chemical production of species i. The
dimensions of Sφ are (phi dimensions) per unit volume, per unit time. This source term in
the balance equation, with dimensions of (phi dimensions) per unit time, is written as
follows.

Source = S x y z [1-5]

All the terms in the equations [1-2] to [1-5] are to be substituted into equation [1-1]
where they will be added to each other or subtracted from each other. In this case, each
term must have the same dimensions. The storage term in equation [1-2] obviously has
4
the dimensions of -dimensions per unit time. In order to have consistent dimensions
when equations [1-2] to [1-5] are substituted into equation [1-1], all the terms in
equations [1-2] to [1-5] must have these dimensions.

For example, consider the storage term in the energy balance equation. In this equation,
represents the energy E, so the -dimensions would be energy dimensions. In the SI
system, the energy would have units of joules (J), and the units for would be J/kg. The
ρφ product would have units of J/m3. The derivative of the ρ product with respect to
time would have units of J/m3-s. Multiplying this term by ΔxΔyΔz gives the final units
for the storage term in the energy equation as J/s.

A similar analysis can be done for both the inflow and outflow terms. The dimensions of
ρu ΔyΔz are (mass/length3)(length/time)( -dimensions/mass)(length2) which equates
to -dimensions divided by time dimensions. For the energy equation, this would be
joules per second, which is consistent with what we found earlier for the storage (or
transient) term. Thus all terms in equation [1-1] to [1-5] will have dimensions of
divided by dimensions of time. If the source term in equation [1-5] must have
dimensions of divided by dimensions of time, the Sφ term must have dimensions of
divided by dimensions of time and divided by dimensions of volume. In the energy
equation, the Sφ term would then have units of J/s-m3.

Combining all the quantities in equations [1-2] to [1-5], according to the verbal equation
[1-1] and dividing by the volume, ΔxΔyΔz, gives the following result.

  u  −  u  v y + y
−  v  w −  w
+ x + x
+ + z + z
= S [1-6]
x y z

t x y z

Taking the limit as the differentials approach zero gives the definitions of the partial
derivative. This limit results in the following differential equation.

  u  v  w
+ + + = S* [1-7]
t x y z

Where we have defined

5
 Lim
S* = S [1-8]
xyz → 0

Notation to Simplify Multidimensional Equations

The repetition in the directional terms is simplified by notation known as Cartesian tensor
notation (sometimes called the Einstein convention). In this notation, the coordinate
directions and velocity components are given numerical subscripts. The repetition of an
index implies summation over all three coordinate directions (two directions in a two-
dimensional problem.) The usual coordinate system (x, y, z) is written as (x1, x2, x3).
The velocity components that we have written as (u, v, w) could also be written as (ux, uy,
uz). In Cartesian tensor notation, the velocity components are written as (u1, u2, u3).
Using the numerical notation for the spatial coordinates and the velocity components we
can rewrite equation [1-7] as follows.

  u  v  w   u1  u 2  u3

+ + + = + + + = S [1-9]
t x y z t x1 x2 x3

With a conventional summation notation, we could write the second part of equation [1-
9] as follows.

 3
 u i
+ = S [1-10]
t i =1 xi

In the Cartesian tensor notation the summation symbol is not used. Summation over
repeated indices is implied. This we would write equation [1-10] as follows in this
notation.

  u i
+ = S [1-11]
t xi

Equation [1-11] is the same as equations [1-7] and [1-10]; it is just using a more compact
notation to show that there are similar terms in each coordinate direction.

Forces on a Fluid Element – The Source Term for Momentum

6
The source term for the momentum balance equation comes from Newton’s second law
expressed as the statement that the net force is the rate of change of momentum. Thus,
the balance equation for the momentum must have a source term has the dimensions of
force. The analysis of forces on a fluid element considers two kinds of forces. Forces
such as gravity or electromagnetic forces that act over the entire volume of the body are
called body forces. In contrast, forces such as pressure and viscous stresses act on the
surfaces of a fluid element.

The main body force that is considered in fluid dynamics problems is gravity. The
general body force is written as ρB. When written this way, B must have dimensions of
acceleration. This can be seen as follows. The source term for all balance terms, , must
have dimensions of per unit time, per unit volume. Since the dimensions of
momentum are mass times distance over time, the dimensions of the source term for
momentum (momentum per unit volume per unit time) must be (mass times distance over
time) per unit volume per unit time. Thus the source term represented by ρB has
dimensions of mass divided by distance-squared and time-squared. Since the dimensions
of ρ are mass over distance-cubed, the resulting dimensions for B must be distance over
time-squared or acceleration. Since the body force has three directional components we
will write these components as Bx, By, and Bz. Alternatively, using the Cartesian tensor
notation, we will write these components, in general, as Bi, i= 1, 2, 3.

Surface forces are represented by the notation σij, which denotes the normal or shear
stress force on a face of the fluid facing direction i. This force is resolved into a
component in each coordinate direction. The j subscript denotes a particular direction.
For example, σyx denotes a force, acting in the x direction, on a control volume surface,
facing the y direction.

By convention, the force or stress is considered positive when it is exerted by the fluid
above an element on the fluid below an element. The net force on a fluid element due to
a single component is the difference between forces on two sides of an element. For
example, the force on the y faces acting in the x direction is shown in the diagram on the
next page.

7
Since stress is the force per unit area, the force that results from the stress, σyx,, is σyx
times the area Δx Δz. The net force in the x direction, due to the stress on the y-facing
faces of the element is given by [σyx|y+Δy - -σyx|y] Δx Δz. If we divide this net force (from
one face in one direction) by the control element volume, Δx Δy Δz, we get the result
shown in equation [1-12] for the net force per unit volume.

Net yFace, xDirection force  yx y + y −  yx  yx y + y

−  yx
= xz =
y y
[1-12]
Unit Volume xyz y

In order to get the correct source term in the differential equation, we have to use the
definition in equation [1-8], which requires that we take the limit as the control volume
approaches zero. In this limit, the last term in equation [1-12] becomes the partial
derivative.

Limit Net yFace, xDirection force  yx

= [1-13]
y → 0 Unit Volume y

There will be similar terms for the net x-direction force from faces in the x and z
directions. Writing out these terms and combining them with the y-face term gives the
net force in the x-direction from surface stresses.

 xx  yx  zx
Net x − direction surface force source = + + [1-14]
x y z

If we repeated the analysis for the net surface-force source term in the y and z coordinate
directions, we would obtain similar terms. These are shown below.

 xy  yy  zy
Net y − direction surface force source = + + [1-15]
x y z

 xz  yz  zz
Net z − direction surface force source = + + [1-16]
x y z

Using numerical subscripts we can write each of these terms in the same form. Further,
using the Cartesian tensor notation we can write the net forces in one direction, say

8
direction j, as a single partial derivative with implied summation over repeated indices.
In this manner, equations [1-14] to [1-16] can be written with the following notation.

 1 j  2 j  3 j iij
Net j − direction surface force source = + + = [1-17]
x1 x 2 x3 xi

In the final term, the repeated index, i, implies summation. The non-repeated index, j,
indicated the direction of the force. In subsequent equations for the work done by the
surface forces we will use a summation over the j index as well.

The Mass Balance Equation

As noted in Table 1-1, = 1 in the balance equation for mass. Since mass is neither
produced nor destroyed – we are ignoring the effects of relativity here – the source term
Sφ = 0. Thus, equation [1-10] can be written as follows for conservation of mass.

  u i
+ =0 [1-18]
t xi

This equation is simple enough so that we can write the full equation without Cartesian
tensor notation in a fairly simple manner as follows:

  u  v  w
+ + + =0 [1-19]
t x y z

This is simply obtained from equation [1-18] by using the definition of the Cartesian
tensor notation and switching to the use of (x,y,z) and (u,v,w) for the spatial and velocity
coordinate systems, respectively. It can also be obtained by setting = 1 and Sφ = 0 in
equation [1-7]. Equation [1-18] (or its equivalent form in equation [1-19]) is known as
the continuity equation.

Using the product rule for differentiation, equation [1-18] may be rewritten as follows

  u
+ ui + i =0 [1-20]
t xi xi

9
Note that both the second and third terms in this equation have repeated indices so that
summation over all coordinate directions is implied for each term. Thus, equation [1-20]
may be written as follows

  u v w    
+  + +  +u +v +w =0 [1-21]
t  x y z  x y z

The four terms with density derivatives in this equation are sometimes written as Dρ/Dt,
which is known as the substantive derivative. For any function, , this derivative is
defined as follows:

D     D  
= +u +v +w or = + ui [1-22]
Dt t x y z Dt t xi

In the second definition above, the Cartesian tensor notation is used and summation is
assumed over the repeated index, i. Physically, the substantive derivative represents the
change in properties with time, following a fixed element of mass as it moves through the
fluid. This is in contrast to the usual partial derivative (with respect to time) that
measures the change in time at a fixed point in space through which the fluid flows.

Using the substantive derivative, the continuity equation may be written as follows.

D  u v w  D u D
+  + + =0 + i =0 +  = 0 [1-23]
Dt  x y z  Dt xi Dt

Here we have introduced a new term, Δ, that is called the dilatation. Its definition is
shown in the equation below. For a constant density fluid the continuity equation can be
written as follows.

u v w ui
 + + =0 or  =0 [1-24]
x y z xi

In a vector definition, the dilatation represents the divergence of velocity. In a physical

sense, the dilatation represents the relative change in density with time as we follow a

10
fixed mass of fluid through the flow. If there is a change in density, there must be a net
inflow or outflow of mass to accommodate this change.

Alternative Form for the General Balance Equation and the Conservation Form

The general balance equation may be written in an alternative form as shown below. We
start with equation [1-11] and apply the chain rule of differentiation to obtain the
following result.

    u i  
 + +  +  ui = S [1-25]
t  t xi  xi

The term in brackets is just the continuity equation as shown in equation [1-18]. Since
the right-hand side of this equation is zero, the bracketed term in equation [1-25] is zero,
and we may write the general balance equation as follows.

 
 +  ui = S [1-26]
t xi

Equation [1-26] is the starting point for many analyses in fluid dynamics and convective
heat and mass transfer. It is a valid balance equation. However, researchers in
computational fluid dynamics have found that the finite-difference schemes derived from
this equation have problems in their solution. In particular, the finite-difference
equations derived from the form of equation [1-26] will not conserve mass. (For any
balance equation with a zero source term, the quantity, , in the balance equation should
be conserved. In this special case, balance equations are sometimes called conservation
equations.) However, finite-difference equations based on the original balance equation
form in equation [1-11] will conserve mass. Because of this, all derivations of finite-
difference equations for computational fluid dynamics start with equation [1-11]. This
equation is sometimes called the conservation form because it conserves mass when
converted to a finite difference formulation.* Using this nomenclature we will call
equation [1-26] the “non-conservation” form.

11
The process of deriving equation [1-26] is a two-way street. If we have an equation in
the form given by equation [1-26] we can reverse the derivation to get the form given by
equation [1-11].

The Momentum Balance Equation

The total net force in direction j is the sum of the body force and the net surface force in
that direction. If we add the body force, ρBj, in direction j, to the net surface force in that
direction, given by equation [1-17] we have the following expression in Cartesian tensor
notation.

 1 j  2 j  3 j  ij
Net j − direction source term = + + + B j = + B j [1-27]
x1 x 2 x3 xi

With this force term, the general balance equation [1-11] for the momentum per unit
mass in direction j ( = uj) can be written as shown below. There are three such
equations, one for each coordinate direction.

u j  u i u j  ij
+ = + B j j = 1, 3 [1-28]
t xi xi

For a Newtonian fluid, the stress, σij, is given by the following equation

  ui uj  2
 ij = − P ij +   +  + ( −  ) ij [1-29]
 x j xi  3

Here P is the usual thermodynamic pressure; δij is called the Kronecker delta. It is 1 if i=j
and zero otherwise. The symbols μ and κ are called the dynamic and bulk viscosities; the
latter is only important in high frequency sound wave problems and will not be
considered further in these notes.If you find this notation confusing, you should write
equation [1-29] for σxy, the net surface force on the x face in the y direction. To do this,
you first write equation [1-29] with i = 1 and j = 2 as subscripts. (What is the value of δ12
according to the definition of δ?) Next, substitute x for i and y for j in the subscripts for
σ. Finally, substitute u and v for u1 and u2, respectively, and substitute x and y for x1 and
x2, respectively. The result will have a notation that is more readable than the notation in
12
equation [1-29]. However, this more obvious notation for directions requires you to write
nine such equations, all of which will have the same form as equation [1-29], to represent
all the surface force terms.

We obtain the momentum balance equation for a Newtonian fluid by substituting the
definition of σij in equation [1-29] into equation [1-28].

u j  u i u j    u uj  
+ = − P ij +   i +  + ( − 2  ) ij  + B j j = 1, 3 [1-30]
t xi xi 
  x j xi  3 

The right-hand side of this equation has in implied summation in the repeated index i.
Since δij=0, unless i = j, the terms multiplied by δij can be written only once. When this is
done equation [1-30] can be written as follows:

u j  ui u j P     ui  u j    2 
+ =− +  + + ( − 3  ) + B j j = 1, 3 [1-31]
t xi x j xi   x j xi  x j  

The reader who is still unsure of the i and j subscripts can write equation [1-31], three
times, once for each coordinate direction. The equation for x-momentum is shown
below. The momentum equations for the other two coordinate directions are left as an
exercise for the reader.

u uu vu wu

+ + + = Bx
t x y z
[1-32]
P    u      v  u      w  u    2 
− + 2    +    +  +  +  + ( −  ) 
x x  x  y   x y  z   x z  x  3 

The original momentum balance in equation [1-28], with σij on the right-hand side, is
valid for any relation between σij and velocity gradients. This equation is starting point
for analysis of non-Newtonian fluids. Before such an analysis can proceed, it is
necessary to determine the relationship between σij and other flow properties. For the
remainder of these notes we will use equation [1-28] and consider only Newtonian fluids.

The Energy Balance Equation

13
The source terms in the balance equation for total energy (thermodynamic plus kinetic
plus potential) are the net rate of heat addition plus the net rate at which work is done on
the fulid. The heat flow term is written in terms of the heat flux (heat flow per unit area)
in a particular direction, i. This directional heat flux is given the symbol qi. In the
simplest flows, this equation is given the Fourier Law for heat conduction.
T
qi = −k [1-33]
x i

The net contribution to the source term from the heat flow is found from an analysis that
is similar to that done for the surface stresses in equations [1-12] and [1-13]. The net heat
inflow in the x direction is qx|x – qx|x+Δx. The heat flow at the x+Δx face has a negative
sign because the heat flows out of the element at this face if qx is positive. The x-
direction source term is the net heat flow per unit volume.

Net xDirection heat qx − qx x qx − qx

= − x+x yz = − x + x x
[1-34]
Unit Volume xyz x

In order to get the correct source term in the differential equation, we have to use the
definition in equation [1-8], which requires that we take the limit as the control volume
approaches zero. In this limit, the last term in equation [1-34] becomes the partial
derivative.

Limit Net xDirection heat source q

=− x [1-35]
x → 0 Unit Volume x

This analysis can be repeated in the other two coordinate directions. This gives the total
heat source term as the sum of three partial derivatives that can be represented as one
using the Cartesian tensor notation.

q x q y q x q
Heat Rate = − − − =− i [1-36]
x y x xi

The derivation of the work term requires an analysis of the body and surface forces. The
rate at which work is done at any point in the fluid is the product of the force in a given
direction times the velocity component in that direction. For the body-force terms, the
14
work rate is simply the product of ρBj with the velocity component uj summed over all
three coordinate directions.

Body-force work rate = ρ(uBx + vBy + wBz) =ρuiBi [1-37]

The analysis of the work done by surface forces is more complex. On each face there are
forces acting in all three coordinate directions. These forces must be multiplied by the
appropriate velocity component. In addition, the work at the xi + Δxi face is done on the
fluid and is thus positive; at the xi face, the element does work on the adjacent fluid, so
the work term is negative. The difference between there two is the net work. As a
specific example, consider the work done on the faces in the y direction.

The work term on each face is given by the following equation:

y-face surface force work = (uσyx + vσyy + wσyz)Δx Δz =uiσiy Δx Δz [1-38]

If we follow the same analysis used for the net heat source term above or the net stress
force source term in equations [1-12] and [1-13], we will obtain the following result.

(u yx + v yy + w yz ) ui yi
Net yFace Surface Force Work = = [1-39]
y y

A similar analysis in the other coordinate directions provides analogous terms. The total
work term is given by adding the terms in all coordinate directions.

ui xi ui yi ui zi ui ji

Net Surface Force Work = + + = [1-40]
x y z x j

The final term in equation [1-39] has two repeated indices with an implied summation.
Thus, this single term represents the sum of nine different partial derivatives. The reader
should ensure that she or he can write out all nine terms using the coordinates x, y, and z,
and the corresponding velocity components u, v, and w.

15
The energy balance equation is another form of the general balance formula given by
equation [1-11]. Here the per-unit-mass quantity, φ, in the energy equation is the sum of
the (per-unit-mass) thermodynamic internal energy, e, and the kinetic energy, V2/2. The
source term is the sum of the heat source from equation [1-36], the body force work from
equation [1-37], and the surface force work from equation [1-40]. This gives the
following result for the energy balance equation.

 (e + V 2 / 2)  ui (e + V 2 / 2) q ui ji
+ =− i + + ui Bi [1-41]
t xi xi x j

We could stop here, since we now have an energy balance equation. However, many
different forms of the energy equation are used in practice. These equations include a
separate consideration of the thermodynamic internal energy and kinetic energy balances
and the substitution of other thermodynamic properties – enthalpy and temperature – in
place of the internal energy. The derivation of these equations is presented below.

Alternative Energy Balance Equations

The kinetic energy can be eliminated from the total energy equation as follows. First, we
can use the result of equation [1-26] to cast the energy equation in [1-41] into a non-
conservation form.

(e + V 2 / 2) (e + V 2 / 2) q ui ji

 +  ui =− i + + ui Bi [1-42]
t xi xi x j

Similarly, we can use the result of equation [1-26] to write the momentum balance from
equation [1-28] in a non-conservation form.

u j u j  ij
 +  ui = + B j j = 1, 3 [1-43]
t xi xi

Equation [1-43] represents three different equations for the conservation of momentum,
one in each coordinate direction. The next step in the current derivation is simplified by
the use of the Cartesian tensor notation. If we had three equations for x, y, and z
momentum, we would multiply the x-momentum equation by u, the y-momentum
16
equation by v, the z-momentum equation by w and add the three results. We can obtain
the same result by multiplying equation [1-43] by uj, and applying the summation
convention.

 u j u j  ij 
u j  +  ui = + B j  [1-44]
 t xi xi 

Since each term in equation [1-44] now has a repeated j subscript, we have an implied
sum over this subscript. Thus, equation [1-44] represents the task we outlined above of
multiplying each momentum balance equation by the corresponding velocity component
and adding the results. We can make one further simplification to equation [1-44]. In the
two derivatives on the left-hand side we have terms of the form ujduj. From the
summation convention, we know that this term, with its implied summation, is equal to
udu + vdv + wdw = d(u2/2) + d(v2/2)+d(w2/2) = d(u2 + v2 + w2)/2 = d(V2/2). With this
relationship, equation [1-44] may be written as follows.

 (V 2 / 2)  (V 2 / 2)  ij
 + ui = uj + u j B j [1-45]
t xi xi

This equation can be subtracted from the total energy balance in non-conservation form
given by equation [1-42] to get a balance equation for the thermodynamic internal
energy, e.

e e q ui ji  ij q u
 +  ui =− i + −uj = − i +  ji i [1-46]
t xi xi x j xi xi x j

Of course we can use the result that the left-hand side of this equation can be written in
either the non-conservation form, shown above, or the conservation form. We use the
result that the general non-conservation form in equation [1-26] is equivalent to the
general conservation form in equation [1-11]. Applying this general result to equation [1-
46] gives the equivalent conservation form below.

e  u i e q u
+ = − i +  ji i [1-47]
t xi xi x j

17
Further versions of the thermodynamic energy equation can be derived. Equation [1-46]
is the usual starting point for these derivations. The first step is to introduce the enthalpy,
defined as h = e + P/ρ. Differentiating this enthalpy definition gives the result that dh =
de + dP/ρ - Pdρ/ρ2. We can use this result to substitute the enthalpy for the internal
energy in equation [1-46]. To simplify the derivation, we use the substantive derivative
De/Dt.

e e De  Dh 1 DP P D  qi u
 +  ui = =  − + 2 =− +  ji i [1-48]
t xi Dt  Dt  Dt  Dt  xi x j

From the continuity equation forms shown in equations [1-23] and [1-24], we can write
Dρ/Dt as -ρΔ.

D u
= −  i = −  [1-49]
Dt xi

Substituting this result into equation [1-48] gives:

Dh q u DP
 = − i +  ji i + + P [1-50]
Dt xi x j Dt

We can introduce the temperature, T, by using the general relationship between the
thermodynamic internal energy (or the enthalpy) and other thermodynamic properties:

 T  d
de = cv dT −  P − P  2 [1-51]
 T  

dh = c p dT + 1 − T P 
dP
[1-52]


The quantities βP and κT are fluid properties giving the relative change in density at
constant pressure and temperature, respectively. These are formally defined as follows:

1    1   
P = −   and T =   [1-53]
  T  P   P  T

18
For an ideal gas (P = ρRT), βP = 1/T, and κT = 1/P. In this case, equations [1-51] and [1-
52] simplify to de = cvdT and dh = cpdT.

Substituting equation [1-52] into equation [1-50] gives a differential equation for
temperature that uses the constant pressure heat capacity, cP.

Dh  DT 1 −  PT DP  q u DP
 =   cp +  = − i +  ji i + + P [1-54]
Dt  Dt  Dt  xi x j Dt

Combining the DP/Dt terms gives the following result.

DT q u DP
c p = − i +  ji i +  P T + P [1-54]
Dt xi x j Dt

We can also obtain an equation for the temperature that uses the constant volume heat
capacity, cv. This is done by substituting equation [1-51] into equation [1-46].

e e De  DT 1  T  D  qi u
 +  ui = =  c v + 2  P − P  =− +  ji i [1-55]
t xi Dt  Dt    T  Dt  xi x j

We can use the result that Dρ/Dt = –ρΔ, from equation [1-49] to simplify equation [1-55].

DT q u  T 
cv = − i +  ji i +  P − P  [1-56]
Dt xi x j   T 

We can write equation [1-50] for the enthalpy balance and equations [1-54] and [1-56]
for the temperature in conservation form. These results are given below.

h u i h q u DP
+ = − i +  ji i + + P [1-57]
t xi xi x j Dt

 T u i T  qi u DP

cp  + =− +  ji i +  P T + P [1-58]
 t xi  xi x j Dt

 T u i T  qi u  T 

cv  + =− +  ji i +  P − P  [1-59]
 t xi  xi x j   T 
19
Substitutions for Stresses and Heat Flux

Up to this point we have kept the equations completely general. In order to solve an
energy equation we have to have some definition for the heat flux and the stress forces in
terms of other fluid properties. For the stress terms, we will use the relations for a
Newtonian fluid given by equation [1-29]. For the heat flux, we will use the Fourier
relationship given in equation [1-33]. In more complex problems, particularly those
involving high-temperature reacting flows, the usual Fourier heat flux may need to be
augmented by one or more of the following:

• radiation heat flux,

• diffusion-thermo, i.e., a heat flux due to a temperature gradient, and

• enthalpy diffusion (a usually ignored term associated with the diffusion of species
with different enthalpy values,

Using only the Fourier Law heat transfer, the source term involving the heat flux in the
energy balance equation can be written as follows.

qi   T   T  T  T  T
− =−  − k  = k = k + k + k [1-60]
xi xi  xi  xi xi x x y y z z

Substituting equation [1-60] and equation [1-29] for the Newtonian stress relation into
equation [1-47] for the thermodynamic internal energy balance, we obtain the following
result.

e  u i e  T   u u j  
 + ( − 2  ) ij  ui
+ = k + − P ij +   i + [1-61]
t xi xi xi   x 
 j xi  3  x j

From the definition of the Kroenecker delta, (δij = 1 if i = j and zero otherwise), and the
definition of the dilatation, Δ, we can simplify the terms in this equation that involve δij
as follows.

u i u j
 ij = = [1-62]
x j x j

20
With this simplification, equation [1-61] becomes.

e  u i e  T  u uj  u i 2
+ = k − P +   i +  + ( −  )2 [1-63]
t xi xi xi  x  x
 j xi  j 3

The terms that are multiplied by the viscosity can be shown to be a sum of perfect
squares, which must always be positive. These terms represent the dissipation of
mechanical energy into heat. They are usually small except for high Mach number flows.
These terms are defined as the dissipation and are usually given the symbol, Φ. In there
notes we will use the symbol ΦD to represent the dissipation to avoid confusion with the
general quantity in a balance equation. The dissipation is simply defined as all the terms
in equation [1-63] that contain the viscosity.

 u u j  u i 2
ΦD =   i +  + ( −  )2 [1-64]
 x  x
 j xi  j 3

With this definition, equation [1-63] may be simply written as follows.

e  u i e  T
+ = k − P + ΦD [1-65]
t xi xi xi

The temperature gradient in the Fourier law conduction term may also be written as a
gradient of enthalpy or internal energy by using equations [1-51] or [1-52].

T 1 e 1  T P  1 
= +  − P 2 [1-66]
xi cv xi cv   T   xi

T 1 h 1 − T P P
= − [1-67]
xi c p xi c p xi

Substituting equation [1-67] into equation [1-65] gives the following result.

e  u i e  k e  1  T P  1 
+ = − P + ΦD +  − P 2 [1-68]
t xi xi cv xi xi cv   T   xi

21
In the deriving equation [1-68], we started with the energy balance equation for the
thermodynamic internal energy and took the following steps: (1) substituted the Fourier
Law expression for the heat flux, (2) substituted the Newtonian fluid stress relations, (3)
did some simplifications and defined a dissipation term, and (4) substituted an energy
gradient for a temperature gradient. If we started with equation [1-57] for the enthalpy
balance and took the same steps we would obtain the following result.

h u i h  k h  1 − T P  P DP
+ = + ΦD +   + [1-69]
t xi xi c p xi xi  c p  xi Dt

If we started with the temperature equations [1-58] or [1-59], involving cv or cp,

respectively, and repeated the process outlined in the paragraph above equation [1-69] we
would obtain the equations shown below. (In these cases, we keep the temperature
gradient in the final equations instead of substituting an energy or enthalpy gradient.)

 T u i T   T DP
cp  + = k + ΦD +  P T [1-70]
 t xi  xi xi Dt

 T uiT   T T
cv  + = k + ΦD + P  [1-71]
 t xi  xi xi T

The Species Balance Equation

For this equation, φ is the mass fraction of species K, W(K). (The mass fraction symbol
W(K) – for weight fraction, which is the same as mass fraction – to avoid confusion with
the w velocity component and the mass, m. A superscript is used for the species index to
avoid confusion with the coordinate subscripts.) The source term is due to the diffusive
flux and the species production by chemical reaction. In a multicomponent mixture, the
different species will move at different velocities as they mix. (If a fully mixed system,
all species move with the same velocity.) The velocity component of species K in the ith
direction is called the particular velocity of species K (in direction i) and is given the
symbol ui(K). The usual velocity component that was used in our balance equations, uj, is
called the mass-average velocity, when dealing with mixtures. This mass-average
velocity is given by the following equation.
22
 ui = W ( k ) ui( K ) [1-72]
K

The diffusive flux of species K, in the ith direction, has the symbol, ji(K), and is given by
Fick’s Law, shown in the following equation.

ji( K ) = W ( K ) ( ui( K ) − ui ) [1-73]

If we look at the inflow and outflow of species K, with its particular velocity component,
ui(K), across the various faces of a control volume, we can write a balance equation by
noting that the net storage and outflow must be balanced by the chemical production of
species K, r(K). This gives the following balance equation.

W ( K )  W ( K )ui( K )
+ = r (K ) [1-74]
t xi

Substituting equation [1-73] into equation [1-74] gives the following result.

W ( K )  W ( K )ui j ( K )
+ = − i + r(K ) [1-75]
t xi xi

This is the basic species balance equation. In order to proceed further with the solution
of this equation we need to relate the diffusive flux to fluid properties. For an isothermal,
binary system the diffusive flux is given by Fick's Law. This defines a property known
as the binary diffusion coefficient; for the diffusion of two species, A and B, this
coefficient is written as DAB. Fick’s law is then written as follows.

W ( A)
ji( A) = − D AB [1-76]
xi

A full consideration of the relationship between diffusive fluxes and flow properties
would have to consider many additional effects that can contribute to a diffusive flux.
These other effects include thermal diffusion, pressure diffusion, body force diffusion,
and the diffusion of one species due to concentration gradients in other species. An
accurate picture of diffusion in multicomponent mixtures may not be helpful when
turbulent flows are involved. In such cases, a simple picture of the laminar diffusion
23
process will suffice since that process will usually be overwhelmed by the turbulent
diffusion. In this simple picture, an average diffusion coefficient for the individual
species in the mixture will be assumed and the diffusive flux for any individual species
will be assumed to follow Fick's law. For simplification we assume that we can define an
average diffusion coefficient for a species K in a mixture, DK,Mix, by one of the following
equations.

1 W (L)
DK , Mix = W DKL
(L)
or = [1-77]
L DK , Mix L D KL

With this definition of an average diffusion coefficient, we can write the diffusive flux,
for a mixture with any number of components by the following approximate equation.

W ( K )
ji( K ) = − D K , Mix [1-78]
xi

If we substitute this equation for the diffusive flux into the species balance given by
equation [1-75], we obtain the following result.

W ( K )  W ( K )ui  W ( K )
+ = D K ,Mix + r (K ) [1-79]
t xi xi xi

This is our general species balance equation.

The General Equation

All the equations derived in the previous section have the following kinds of terms

• A transient term given by a time derivative

• A set of convection terms involving first order derivatives in the three coordinate
directions and the velocity components in those directions

• A “diffusive” term that involves a second derivative in all space dimensions of the
per-unit-mass quantity in the balance equation. All such terms involve
coefficients such as viscosity, thermal conductivity, or the diffusion coefficient.

24
 They are associated with irreversible processes that tend to smooth out gradients
in flows.

• Other terms, including source terms for the quantity in the balance equation

This can be most clearly seen in the species balance equation that we just derived. We
can identify the four kinds of terms in equation [1-79] as shown below.

W ( K )  W ( K )ui  W ( K )
+ = D K ,Mix + r (K )
t xi xi xi [1-80]
Transient Convective Diffusive Source

Here we have identified the second-derivative terms as representing a diffusive process.

In the other equations that we have seen the second derivative terms represented viscous
stresses and heat conduction. All of these processes have similarities at a molecular level
and we can loosely use the name of diffusive for the second derivative terms in all the
other balance equations that we have seen.

In expanding the nomenclature we have written above to other equations, we will rewrite
equation [1-80] in the more general form shown below.

   u i    ( ) 
c + =  + S ( )
 t xi  xi xi [1-81]
Transient Convective Diffusive Source

In this equation c = 1 in all equations except for the temperature-based energy equations
where c may represent cP or cv; Γ(φ) represents some transport property such as the
viscosity or thermal conductivity divided by heat capacity. It will generally have
dimensions of mass/(length times time). In SI units, the typical Γ(φ) will be measured in
kg/m-s. The “Source” term is a combination of true source terms, such as the production
of species by chemical reactions and other terms that do not fit into the classification of
transient, convection and diffusive.

In this formulation, we can represent the equation as follows.

STORAGE + CONVECTION = "DIFFUSION" + "SOURCE" [1-82]

25
The words transient and convection accurately describe the terms in the general equation.
The second-derivative term is truly diffusion only in the species balance equation.
However, in all the other equations, these terms represent some sort of transport
(momentum, heat, or mass) that is driven by a gradient term. As noted above, the
"Source" term will contain terms that arise from phenomena other that true source terms.
In many cases, these other terms may be negligible, especially in cases where properties
are constant or in the case of ideal gases.

When we derive the numerical analysis equations for the differential equations
representing balance phenomena, we will start with equation [1-81]. The numerical
algorithms that we derive for this equation can be applied to any balance equation. We
do need to keep in mind the actual meaning of the different terms in equation [1-81] for
the various balance equations. These terms are shown in Table 1-2. The derivation of
these source terms is outlined below.

The basic momentum balance in equation in the jth direction is given by equation [1-31],
which is copied below.

u j  ui u j P     ui  u j    2 
+ =− +  + + ( − 3  )  + B j j = 1, 3 [1-31]
t xi x j xi   x j xi  x j  

This equation may be rearranged as follows

u j  ui u j P  u j  ui   2 
+ =− +  +  + ( − 3  ) + B j [1-83]
t xi x j xi xi xi x j x j  

We can identify the transient, convective, and diffusive terms. All the remaining terms
are considered the “Source” term. Thus, we define the “Source” term for the momentum
in the jth direction as follows.

P  u   2 
S ( j) = − +  i+ ( − 3  )  + B j [1-84]
x j xi x j x j  

With this definition, the general form for the momentum balance equation in direction j
can be written as follows.
26
 u j  u i u j  u j
+ =  + S ( j) [1-85]
t xi xi xi

For many numerical methods, it will be useful to address the pressure term in the
momentum as a part of the solution algorithm. We consider the subdivision of the source
term into a pressure gradient term and the remaining terms in equation [1-85] in the text
following Table 1-2.

We can get a special case of the momentum equation for no body force terms and
constant properties. For no body forces, Bj = 0. For constant density, Δ = 0. For
constant viscosity, we can bring μ outside the derivative. When we do all these three
things, we get the equation below that we can manipulate in steps. First, we change the
order of differentiation in the mixed second derivatives that are multiplied by the
viscosity. When we do this, we obtain the dilatation, Δ, which is zero for a constant
density fluid. Thus, for the simplest case considered here (constant density and viscosity
and zero body force terms) the source term in the momentum equation is simply the
pressure gradient.

P  ui P  ui p  P
S ( j) = − + =− + =− + =− [1-86]
x j xi x j x j x j xi x i x j x j

The source terms in the various energy equations are simply found by inspection of
equations [1-68] through [1-71] which are copied below. In each of these equations, the
transient and convective terms are on the left-hand side and the diffusive term is the first
term on the right-hand side. All the remaining terms on the right-hand side, for each of
these, equations, is the source term. The source term for each of these equations is shown
in Table 1-2.

e  u i e  k e  1  T P  1 
+ = − P + Φ D +  − P 2 [1-68]
t x i x i c v x i x i c v   T   x i

h u i h  k h  1 − T P  P DP
+ = + ΦD +   + [1-69]
t xi xi c p xi xi  c p  xi Dt

27
  T u i T   T DP
cp  + = k + ΦD +  P T [1-70]
 t xi  xi xi Dt

 T u i T   T T
cv  + = k + ΦD + P  [1-71]
 t xi  xi xi T

The various forms of the energy equation may be simplified in the following cases:

• For low Mach number flows, the dissipation is very small and can be ignored

• For constant density flows, the dilatation, Δ, is zero as are any other density
derivatives.

• For ideal gases, βP = 1/T, κT = 1/P

For low-Mach-number, constant density flows, the source terms in the energy equation
and the temperature equation involving cv are zero. There is a difference between an
incompressible flow and a constant density flow. A constant density flow is one in which
the density is constant. An incompressible flow is one in which the Mach number is
nearly zero. In an incompressible flow the density may change because of temperature,
but the effect of pressure on density may be ignored. Thus, a furnace may be analyzed as
an incompressible flow. The large temperature variation in the furnace leads to large
density variations. However, the pressure changes in the furnace that are used in the
momentum equation are very small compared to the average pressure. Such a flow may
be analyzed as an incompressible flow, but it is not a constant density flow.

28
Table 1-2 – Identification of Terms in the General Balance Equation
φ c Γ(φ) S(φ)
1 1 0 0
u = ux = P  u  v  w   2 
1 μ − +  +  +  +  ( −  )  + Bx
u1 x x x y x z x x  3 
v = uy = P  u  v  w   2 
1 μ − +  +  +  + ( −  )  + B y
u2 y x y y y z y y  3 
w = uz = P  u  v  w   2 
1 μ − +  +  +  + ( −  )   + Bz
u3 z x z y z z z z  3 

 1  T P  1 
e 1 k/cv − P + ΦD +  − P 2
xi cv   T   xi

 1 − T P  P DP
h 1 k/cP ΦD +   +
xi  c p  xi Dt
DP
T cP k ΦD +  P T
Dt
T P
T cv k ΦD + 
T
W(K) 1 ρDK,Mix r(K)

Pressure in the momentum equations

CFD solutions must determine the pressure as well as the velocity components. Because
of this, it is necessary to explicitly treat the pressure term, which is listed as part of the
source term in Table 1-2. To do this, we define a new source term for the momentum
equations, S*(j), which is the same as the source term in Table 1-2, except that the
pressure term is removed. To do this we rewrite the general momentum equation [1-85]
as follows.

u j  ui u j  u j P  u j
+ =  + S ( j) = − +  + S *( j ) [1-87]
t xi xi xi xi xi xi

The new source terms for equation [1-87] are shown in Table 1-3.

29
Table 1-3 – Identification of Terms in the General Momentum Equation
φ c Γ(φ) S*(φ)
1 1 0 0
u = ux =  u  v  w   2 
1 μ  +  +  + ( −  )  + Bx
u1 x x y x z x x  3 
v = uy =  u  v  w   2 
1 μ +  +  +  + ( −  )  + B y
u2 x y y y z y y  3 
w = uz =  u  v  w   2 
1 μ  +  +  + ( −  )  + Bz
u3 x z y z z z z  3 

Note that the general momentum equation [1-87] has the same form as the general
transport equation [1-81], except for the addition of the pressure gradient term. This
allows differencing schemes developed for the general transport equation to be applied to
the momentum equation with additional steps required to handle the pressure terms.

Time derivative approach

An alternative approach to the treatment of pressure, often used in compressible flows

and in some finite-element approaches, uses first-order, time-derivative equations. In this
approach the various flux terms and their constitutive relations are treated as separate
equations in the numerical algorithms. In these approaches the total stress term, σij, is
written as the sum of the pressure, p δij, and the viscous stress, ij.

 ij = −P ij +  ij [1-88]

Comparing this equation with equation [1-29] for σij shows that the viscous stress term is
given by the following equation.

  ui uj  2
 ij =   +  + ( −  ) ij [1-89]
 x j xi  3

Over overall momentum balance, in of the stress, was given by equation [1-28], which is
copied below.

30
 u j  u i u j  ij
+ = + B j j = 1, 3 [1-28]
t xi xi

Using equation [1-88] we can rewrite this general momentum balance equation as follows

u j  uiu j p  ij
+ =− + + B j j = 1, 3 [1-90]
t xi x j xi

 ij
Recall that the summation convention over repeated indices means that the term is
xi
actually the sum of three different terms. (When equation [1-88] was substituted into
p ij p
equation [1-28] there was a term − . However, this term may be written as − ,
xi x j

since δij = 0 if i is not equal to j.)

Equation [1-90] uses only first order derivatives, but it is necessary to use equation [1-89]
to compute the viscous stress terms, ij, in this equation. (Because the viscous stress
terms are symmetric – i.e., ij = ji – it is in only necessary to compute six of the nine
different ij terms.)

The solution of CFD equations, using differential equations that have only first
derivatives, also uses equations [1-18], [1-41], and [1-75], copied below, for mass,
energy, and species balance.

  u i
+ =0 [1-18]
t xi

 (e + V 2 / 2)  ui (e + V 2 / 2) q u 
+ = − i + i ji + ui Bi [1-41]
t xi xi x j

W ( K )  W ( K )ui j ( K )
+ = − i + r(K ) [1-75]
t xi xi

31
We can use equation [1-88] to replace σij, in the energy balance equation by the pressure
and the viscous stress, ij. Doing this and using the fact that δij = 0 if i is not equal to j,
allows us to rewrite the energy balance equation as follows:

 (e + V 2 / 2)  ui (e + V 2 / 2) q Pu j ui ji
+ =− i + + + ui Bi [1-91]
t xi xi x j x j

Since the subscripts are used as summation indices, we can replace i by j on the
convection term and the heat flux gradient term without changing the result. Doing this
allows us to rewrite equation [1-91] as follows

 (e + V 2 / 2)  u j (e + V / 2) q Pu j ui ji
2

+ =− j + + + ui Bi [1-92]
t x j x j x j x j

CFD practitioners who use the first-derivative approach derive a general algorithm from
writing a vector form of the equation. This form summarizes the continuity equation [1-
18], the three momentum equations implied in equation [1-90], the energy equation of
equation [1-92] and the species balance of equation [1-75]. The equation they use is
written as follows.

U E F G
+ + + =H [1-93]
t x y z

Where the terms U, E, F, G, and H are vectors (one-dimensional column matrices)

defined by the following equations:

   U 1   0   h1 
  U   Bx   
 u   2   h2 
 v  U 3   B y   h3 
U=  = U  H=  = h  [1-94]
w Bz
   4    4
  ( e + V 2
/ 2)  U 5    (uBx + vBy + wB z )  h5 
 W ( K )  U    h 
   6  r (K )   6

32
  u   e1 
 uu + p −  xx   
  e2 
 vu −  xy  e3 
E=  = e  [1-95]
wu −  xz
   4
u[  (e + V / 2) + p ]) − u xx − v xy − w xz − qx  e5 
2

  e 
 uW ( K ) + jx( K )   6

 v   f1 
 uv −  yx   
   f2 
 vv + p −  yy   f3 
F= =f  [1-96]
wv −  yz
   4
v [  (e + V / 2) + p ] − u yx − v yy − w yz − q y   f 5 
2

   
 vW ( K ) + j y( K )   f6 

 w   g1 
 uw −  zx   
   g2 
 vw −  zy   g3 
G=  = g  [1-97]
ww + p −  zz
   4
 w[  (e + V / 2) + p ] − u zx − v zy − w zz − qz   g5 
2

  g 
 wW ( K ) + jz( K )   6

The algorithms that use this format typically solve for the U matrix components. (In
equations [1-94] and [1-98], the Ui terms refer to the components of the U vector, not to
velocity components, ui.)

33
TEXT / REFERENCE BOOKS
1. Anderson, J.D., Computational Fluid Dynamics: The Basics with Applications, 2nd
Edition, McGraw Hill International
Editions, 2012.
2. Fletcher, C.A.J., Computational Techniques for Fluid Dynamics, 3rd Edition Springer
Verlag, 2001.
3. Versteeg, H.K. and Malalasekera, W., An Introduction to Computational Fluid
Dynamics: The Finite Volume Method,
2nd Edition, Pearson Education Ltd., 2007.
4. Chung T.J., Computational Fluid Dynamics, 2nd Edition, Cambridge University Press,
2002.

34
UNIT – II – Computational Fluid Dynamics – SCHA1403

II. Turbulence and its Modelling

35
 Turbulence

Whenever turbulence is present in a certain flow it appears to be the

dominant over all other flow phenomena. That is why successful modeling
of turbulence greatly increases the quality of numerical simulations.

All analytical and semi-analytical solutions to simple flow cases were

already known by the end of 1940s. On the other hand there are still many
open questions on modeling turbulence and properties of turbulence it-self.
No universal turbulence model exists yet.

Ideal turbulence model

Solving CFD problem usually consists of four main components: geometry

and grid generation, setting-up a physical model, solving it and post-
processing the computed data. The way geometry and grid are generated,
the set problem is computed and the way acquired data is presented is very
well known. Precise theory is available. Unfortunately, that is not true for
setting-up a physical model for turbulence flows.

The problem is that one tries to model very complex phenomena with a
model as simple as possible.

Therefore an ideal model should introduce the minimum amount of

complexity into the modeling equations, while capturing the essence of the
relevant physics.

Complexity of the turbulence model

Complexity of different turbulence models may vary strongly depends on

the details one wants to observe and investigate by carrying out such
numerical simulations. Complexity is due to the nature of Navier-Stokes
equation (N-S equation). N-S equation is inherently nonlinear, time-
dependent, three-dimensional PDE.

Turbulence could be thought of as instability of laminar flow that occurs at

36
 high Reynolds numbers ( Re ). Such instabilities origin form interactions
between non- linear inertial terms and viscous terms in N-S equation. These
interactions are rotational, fully time-dependent and fully three-
dimensional. Rotational and three- dimensional interactions are mutually
connected via vortex stretching. Vortex stretching is not possible in two
dimensional space. That is also why no satisfactory two-dimensional
approximations for turbulent phenomena are available.
Classification of turbulent models

Nowadays turbulent flows may be computed using several different

approaches. Either by solving the Reynolds-averaged Navier-Stokes
equations with suitable models for turbulent quantities or by computing
them directly. The main approaches are summarized below.

Reynolds-Averaged Navier-Stokes (RANS) Models

Eddy-viscosity models (EVM)

One assumes that the turbulent stress is proportional to the

mean rate of strain. Further more eddy viscosity is derived
from turbulent transport equations (usually k + one other
quantity).

Non-linear eddy-viscosity models (NLEVM)

Turbulent stress is modelled as a non-linear function of mean

velocity gradients. Turbulent scales are determined by solving
transport equations (usually k + one other quantity). Model is
set to mimic response of turbulence to certain important types
of strain.

Differential stress models (DSM)

This category consists of Reynolds-stress transport models

(RSTM) or second-order closure models (SOC). One is
required to solve transport equations for all turbulent stresses.

37
 Computation of fluctuating quantities

Large-eddy simulation (LES)

One computes time-varying flow, but models sub-grid-scale

motions.

Direct numerical simulation (DNS)

No modelling what so ever is applied. One is required to

resolve the smallest scales of the flow as well.

Extend of modelling for certain CFD approach is illustrated in the

following figure Figure 2.1. It is clearly seen, that models computing
fluctuation quantities resolve shorter length scales than models solving
RANS equations. Hence they have the ability to provide better results.
However they have a demand of much greater computer power than those
models applying RANS methods.

Figure 2.1 Extend of modelling for certain types of turbulent models

REYNOLDS-AVERAGED NAVIER-STOKES MODELS

The following chapter deals with the concept of Re

Further on methods how to include these ideas into certain numerical

38
 models are presented.

Equations describing instantaneous fluid motion

For easier understanding of certain mathematical ideas it is convenient to

briefly revise N-S equations describing instantaneous fluid motion at the
beginning. All variables describing instantaneous flow are marked with a
tilde.
Reynolds averaging
The concept of Reynolds averaging was introduced by Reynolds in 1895.
One may consider Reynolds averaging in many different ways. There are
three most common perceptions of this term: time averaging, space
averaging or ensemble averaging.Time averaging is appropriate when
considering a stationary turbulence. That is when the flow does not vary on
the average in time. In such cases time average is defined by:
The closure problem

One could pretend that Reynolds stress is indeed a stress and try to write
constitutive relations similar to those for viscous stress. However there is an
important difference among these two stresses. Viscous stress is a property
of a fluid. That is why separate experiments can be carried out in order to
determine corresponding constitutive relations. These relations are valid
then whenever a flow in that particular fluid is observed. On the other hand
Reynolds stress is a property of the flow. Hence it is dependent on the flow
variables them-selves. That is the reason why it changes from flow to flow
and no general constitutive relations are available.

Laminar flow, infinitesimal fluctuations and superposition

One solution to the closure problem is to treat the flow as a laminar flow
with fluctuations superimposed. One subtracts the averaged momentum
equation from equation describing instantaneous motion. The result for
fluctuating motion reads:
u u (v ) u
p ij U u u
uj uj i
xi
xj 39 x
i
j xj
j
x ij
 i +Uj i
x +

TEXT / REFERENCE BOOKS

1. Anderson, J.D., Computational Fluid Dynamics: The Basics with Applications, 2nd
Edition, McGraw Hill International
Editions, 2012.
2. Fletcher, C.A.J., Computational Techniques for Fluid Dynamics, 3rd Edition Springer
Verlag, 2001.
3. Versteeg, H.K. and Malalasekera, W., An Introduction to Computational Fluid
Dynamics: The Finite Volume Method,
2nd Edition, Pearson Education Ltd., 2007.
4. Chung T.J., Computational Fluid Dynamics, 2nd Edition, Cambridge University Press,
2002.

40
UNIT – III – Computational Fluid Dynamics – SCHA1403

41
 III. Finite Volume Method

Heat Conduction, Convection And Diffusion

CFD provides the solution to the governing equations of the flow subject to a particular
initial and boundary conditions. These equations are highly nonlinear and very difficult to
solve even numerically. In applying these equations to a particular problem, some of the
terms may disappear or be negligible which makes the solution much simpler. Various
numerical techniques are developed for each of the particular application of the general
flow equations and their simplified forms. In order to introduce various computational
techniques we will first consider a simple form of the momentum equation, and then
discretize various forms of that equation. The momentum equation (2) for a 1 -
dimensional, incompressible flow with no body force, and constant properties reduces to

The first term in equation is the transient term, the second is the convective term, the
third is the diffusive term, and the fourth is the pressure term. We will consider various
combinations of the terms in this equation and discuss the methods to solve them.

Transient-Diffusive Terms
Consider only the 1st and the 3rd terms in the above equation and, to further simplify,
assume ν=1:

This is the transient diffusion equation which consists of a first derivative in the time
direction t and a second derivative in the space direction x . This is a parabolic partial

42
differential equation that can be used to model the temporal changes in the diffusion of
some quantity through a medium. For instance, the transient diffusion of heat
(conduction) in a solid. We will solve this equation using both a finite difference and a
finite element approach.

First we will describe the domain of the problem. Lets assume the diffusion occurs along
a zone with thickness L. The time is usually started from t=0 and it is extended in the
positive direction. Once we have identified the range of this domain, we place points
throughout this domain.. The spacings in the x and t directions can be the same or they
may be different. Each point is labeled using i for special discretization and n for
temporal discretization.

∆x

n+2

n+1

∆t

This procedure is referred to as the grid generation. Once the grid is generated one of the
differencing scheme can be used to discretize the governing equation of the problem,
equation (38). The type of differencing scheme used depends on the particular problem. It
is mainly through testing that one may find the accuracy and efficiency of one scheme
over another. One simple method to discretized the diffusion equation is to use a forward
difference formula for the time derivative and a central difference formula for the spatial
derivative. The discretized equation will then be
43
Note that the velocity at position i and time n+1 depends on the three values at the time
level n. Thus by knowing the values of u at time level n, its value at the next time level
n+1 can be calculated. Therefore, to start the calculation, values of u in all the domain,
e.g. all the x locations, should be known. These known values at t=0 are known as the
initial conditions.

44
We can generate other differencing equations. For instance, the left hand side of equation
can be discretized based on the next time level n+1:

When a direct computation of the dependent variables can be made in terms of known
quantities, the computation is said to be explicit. Some common explicit methods for
parabolic partial differential equations (e.g., equation 38) are:

(1) The Forward Time/Central Space (FTCS) method which is represented by

equation and it is stable for ∆t/∆x ≤1/2.

(2) The Richardson method, where central difference is used for both time and space
and it is unconditionally unstable with no practical value:

(3) The DuFort-Frankel Method, which also uses central difference for both time and
space, but uin in the diffusion term is replace by (uin+1 + uin-1)/2. This modification
makes the difference equations unconditionally stable.

The truncation error for DuFort-Frankel method is order of O[(∆t)2, (∆x)2 , (∆t/∆x)2]. In
equation (43) the only unknown variable is uin+1, therefore, it is explicit.

In equation , several unknown variables are related to the several known variables. This is
referred to as an implicit equation. When the dependent variables are defined by coupled
sets of equations, and either a matrix or iterative technique is needed to obtain the
solution, the numerical method is said to be implicit. Some common implicit methods for
parabolic partial differential equations are:

45
TEXT / REFERENCE BOOKS
1. Anderson, J.D., Computational Fluid Dynamics: The Basics with Applications, 2nd
Edition, McGraw Hill International
Editions, 2012.
2. Fletcher, C.A.J., Computational Techniques for Fluid Dynamics, 3rd Edition Springer
Verlag, 2001.
3. Versteeg, H.K. and Malalasekera, W., An Introduction to Computational Fluid
Dynamics: The Finite Volume Method,
2nd Edition, Pearson Education Ltd., 2007.
4. Chung T.J., Computational Fluid Dynamics, 2nd Edition, Cambridge University Press,
2002.

46
UNIT – IV – Computational Fluid Dynamics – SCHA1403

47
 IV. CFD Methods
(1) Euler Explicit Method: An explicit differencing of equation results in the following
formulation:

u n+1 −u n u n −u n
I i
+c i+1 i
=0 (1)
∆t ∆x

This is an explicit equation since only one unknown, uin+, appears in the equation. This
method is referred to as Euler Explicit Method and, unfortunately, it is unconditionally
unstable and will not work for solving the wave equation. This method is first-order since
the lowest-order term in the truncation error is first order, i.e., ∆t, and ∆x.

(2) First-Order Upwind Method: The Euler method can be made stable by using a
backward difference instead of a forward difference for a positive wave speed:

un+1 −un un −un

I i
+c i i−1
=0 (2)
∆t ∆x

This method is stable for 0 ≤ c∆t/∆x ≤ 1, where c∆t/∆x is referred to as the CFL
(Courant-Friedrichs-Lewy) number. This method is referred to as the First-Order
Upwind Method.

For discretized transport problems, the CFL number determines how many mesh cells, a
fluid element passes during a timestep. For compressible flow, the definition is different.
Here, the CFL number determines how many cells are passed by a propagating

48
Perturbation. Hence, the wave-speed, i.e., fluid speed plus the sound speed, is employed.
For explicit time-stepping schemes, such as Runge-Kutta, the CFL number must be less
than the stability limit for the actual scheme to converge. For implicit and semi-implicit
schemes, the CFL limit does not constitute a stability limit. On the other hand, the range
of parameters in which these schemes converge may often be characterized by the CFL
number.

(3) Lax Method: Another method of making the Euler equation stable is by using an
average value for uin based on the two neighboring points:

uin+1 − (uin+1 + uin−1 ) / 2 + c uin+1 − uin−1 = 0 (3)

∆t 2∆x

This is referred to the Lax Method which is stable for CFL ≤ 1.

(4) Euler Implict Methods are another way of solving Euler equation:

uin+1 − uin + c uin++11 − uin−+11 =0

∆t 2∆x (4)

These methods are unconditionally stable for all time steps, however, a
system of equations must be solved for each time level.

The above methods are all first-order accurate. More accurate second-order methods are
developed to solve the PDEs describing the flow. The commonly used methods are:

(5) Leap Frog Method

u n+1 − u n−1 u n+1 − u n+1
i i i+1 i−1

2∆t +c 2∆x =0 (5)

(6) Lax-Wendroff Method19

uin+1 − uin uin+1 − uin−1 ∆t

+c = c2 2
(u in+1 − 2uin + uin−1 ) (6)
∆t 2∆x 2(∆x)

(7) MacCormack Method20

This is an explicit, predictor-corrector method which is written in the following form.

∆t n
Predictor: (uin+1 )∗ = uin −c (u −u n ) (7)
∆x i +1 i

n+1 1 n n+1 ∗ ∆t n+1 ∗ n+1 ∗

Corrector: = u +(ui ) −c ) −(ui−1 ) (8)
ui
i [ (ui ]
49
2 ∆x

50
Here, (uin+1 )∗ is the predicted value for u at point i and time level n+1. The forward and

backward differencing used in the above equations can be changed depending on the
particular problem.

(8) Second-Order Upwind Method

This is a modification of the MacCormack method where upwind (backward)

differences are used in both predictor and corrector.

∆t n
Predictor: (uin+1 )∗ = uin −c (u −uin−1 ) (9)
∆x i

n+1 1 n n+1 ∗ ∆t n+1 ∗ n+1 ∗ ∆t n n n

u +u
Corrector: ui =2 i +(ui [
) −c ∆x (ui ) −(ui−1 ) ]−c ∆x (u i − 2ui−1 i−2 )
(10)

The fluid dynamics of inviscid flows are governed by Euler equations. These equations
may have different character for various flow regimes. For time-dependent flows, the
equations are hyperbolic for all Mach numbers. Therefore, a time-marching method can
be used to obtain the solution. In steady inviscid flows, the Euler equations are elliptic for
subsonic conditions, and hyperbolic for supersonic conditions. Several simplified form of
the Euler equations are used for inviscid flows. For instance, if the flow is
incompressible, by consider the flow is irrotational as well; a solution to the Laplace’s
equation for the velocity potential or stream function can describe the flow field. The
traditional method of solving hyperbolic PDEs are by the method of characteristics.
Alternatively, there are numerous FDM based solution schemes for such flows.

The following options for the discretization of the convection operator.

(1) Upwind Schemes:

In an upwind (UW) scheme the convection term is formed using a first-order accurate
difference equation equating the velocity derivative to the values at the reference point
and its nearest neighbor taken in the upstream direction. This can give very inaccurate
solutions but they are easy to obtain as they converge readily. For compressible flows,
UW is viewed in a different light. Here, instead of the primitive variables, a set of
characteristic variables is often used. The governing equations for the characteristic
variables are locally hyperbolic. Hence, their solutions are wavelike and upwind
differences are the correct treatment. UW here appears under designations such as flux
splitting, flux difference splitting, fluctuation splitting etc.

(2) Hybrid Schemes:

51
A hybrid scheme, where the upwind scheme is used if the Reynolds number is greater
than two, and central differences are used if the Reynolds number is two or less. This is
more accurate than the upwind scheme but does not converge on some grids of points.

(3) QUICK Upwind Schemes:

The quadratic upstream interpolation for convective kinetics (QUICK) scheme31 is a

quadratic upwind scheme used mainly in the finite volume formulation and is more
accurate than the two schemes described above. This scheme uses a three-point upstream-
weighted quadratic interpolation for cell face values. In the QUICK scheme, one adds
one point in each direction and calculates the derivative using the cubic polynomial
drawn through the four involved points. Local truncation error analysis shows third order
accuracy. The QUICK scheme is unconditionally bounded up to cell Reynold numbers of
5. Beyond this limit, it may become unbounded. The QUICK scheme is normally applied
as a correction to the donor cell scheme. In situations with unboundedness, the correction
may locally be limited, thus reverting to the donor cell scheme. The QUICK scheme has
a somewhat different form in finite volume contexts, since here the differences rather
than the derivatives are of interest.

(4) Power-Law Schemes: Power-law schemes are derivatives of QUICK but are more
accurate.

4.5 Incompressible Navier-Stokes Equations

When considering all the terms in equation a special difficulty arises due to the weak
coupling of the velocity and pressure fields. For the incompressible fluids, the continuity
equation is only function of velocity and not a function of pressure. Only the momentum
equations contain pressure terms. Since most of the terms in the momentum equations are
functions of the velocity components it is natural to use these equations to produce the
solutions for the velocity components. Then, the problem is how to obtain the pressure
solution, since continuity does not contain pressure. A direct method is to discretize all
the equations, i.e., continuity and momentum, and solve them simultaneously. This
results in a very large solution vector that contains all variables and consequently very
large computational effort. There are two commonly used methods to resolve this
problem: (1) pressure-based methods, and (2) methods based on the concept of artificial
compressibility (also known as pseudo-compressibility).

4.5.1 Pressure-Based Methods

In the pressure-based method (PBM), (also known as pressure correction, uncoupled, or

segregated methods) a Poisson equation for pressure corrections is formulated, and then it
is updated for the pressure and velocity fields until a divergence-free velocity field is
obtained. There are numerous variety of this method, some of the more popular ones are
52
the marker-and-cell (MAC) method, SIMPLE and SIMPLER methods, the fractional-step
method, and the pressure-implicit with splitting of operators (PISO) method.

TEXT / REFERENCE BOOKS

1. Anderson, J.D., Computational Fluid Dynamics: The Basics with Applications, 2nd
Edition, McGraw Hill International
Editions, 2012.
2. Fletcher, C.A.J., Computational Techniques for Fluid Dynamics, 3rd Edition Springer
Verlag, 2001.
3. Versteeg, H.K. and Malalasekera, W., An Introduction to Computational Fluid
Dynamics: The Finite Volume Method,
2nd Edition, Pearson Education Ltd., 2007.
4. Chung T.J., Computational Fluid Dynamics, 2nd Edition, Cambridge University Press,
2002.

53
UNIT –V – Computational Fluid Dynamics – SCHA1403

54
 V. Grid Generation
Building a Mesh

One of the most cumbersome and time consuming part of the CFD is the mesh
generation. Although for very simple flows, mesh generation is easy, it becomes very
complex when the problem has many cavities and passages. Mesh generation is basically
the discretization of the computational domain. The mesh in finite difference methods
consists of a set of points, which are called nodes. The finite volume method considers
points that form a set of volumes which are called cells. The finite element methods use
sub-volumes called elements which have nodes where the variables are defined. Values
of the dependent variables, such as velocity, pressure, temperature, etc. will be described
for each element.

Element Form
Various forms of elements can be used. However, the most common type in CFD
programs is a hexahedron with eight nodes, one at each corner, and this is known as a
brick element or volume. For two-dimensional applications the equivalent element is a
four-noded quadrilateral. Some finite volume programs have now been released which
have the ability to use tetrahedral in three dimensions or triangles in two dimensions.
Most finite element CFD codes will allow these elements to be used together with a small
range of other element types. Figure 12 shows some of the commonly used sub-domains.

3-noded triangle 4-noded triangle square

4-noded tetrahedral 6-noded prism 8-noded hexahedral

Typical computational elements.

Before generating the mesh, we should know something about the flow behavior. For
55
instance, where in the flow field we have boundary layers, vortices, large gradients in
pressure or velocity, etc. The mesh size and shape should be such that it can capture the
proper physical conditions that occur in the flow. For regions where large gradients exist,

large number of points within the mesh is needed. This is due to using very simple
variation of the parameter, usually, linear, within the each element. Thus the mesh should
be small enough so that a linear approximation between two points is valid.

This is depicted in Fig. 13, were the variation of function u is given along the coordinate
x. We will use a linear variation between the points of a numerical solution. If a coarse
mesh (∆xL) is used for the numerical calculation of the curve, the results would be far
from the actual variation. However, a fine mesh (∆xS) can produce results which are
close to the actual points. The linear approximation results in large errors where the
gradient of u along x is large.

∆xS ∆xL
∆x

Coarse and fine mesh representation of function u.

One of the main difficulties of mesh generation is that, in many cases we do not know
where the large gradients are. Usually, along the solid surfaces, where the boundary layer
is developed, we need to put more points close to the surface in the direction normal to
the surface. Another example is the large pressure changes close to a shock wave in
compressible flows. Grid refinement is needed to resolve important flow details.
Adaptive grid generation is the solution for complex physical and geometrical problems
in which the location of large gradients is not predictable or varies with time, but this is
beyond the scope of this text. Generally, refinement is needed near walls, stagnation
points, in separation regions, and in wake regions. By increasing number of nodes better
accuracy is achieved. Solution should always (if possible) be based on grid independence
tests with same style and mesh arrangement.

56
Grid generation can be assigned to two distinct categories, structured or unstructured
grids. Relating the mesh structure to the numerical method; finite difference programs
require a mesh to have a regular structure and finite element programs can use a mesh

57
 with an irregular structure. In theory finite volume programs could use a mesh with an
irregular structure, but many implementations insist that the mesh has a regular structure.
When a mesh with a regular structure is used there is an advantage in that the solver
program should run faster than if a mesh with an irregular structure is used. This is due to
the implicit relationship that exists between the number of a cell or a point and the
number of its neighbors in a regular mesh, which enables data to be found easily. No such
relationship occurs for meshes that have an irregular structure and so when trying to find
the values of flow variables in neighboring volumes there must be a computational
overhead. This often takes the form of a look-up table which relates the faces to the cells
or the nodes to the elements.

Structured Grid

The main objective of generating a structured grid is to determine the coordinates

transformation that maps the body fitted non-uniform non-orthogonal physical space
(x,y,z) into the transformed orthogonal computational space (ξ,η,ζ).

Γ*2
Γ2
Γ3 Γ*1
Γ*3
Γ1

y
η

Γ4 Γ*4
x ξ

Physical domain Computational domain

The Transformed Computational Domain

There are two steps in generating a structured grid: a) specification of the boundary point
distribution, b) determination of the interior point distribution. The three popular methods
for generating structured grids are:

58
Conformal mapping method

In a conformal mapping the angles between grid lines in computational and physical
domains are the same. This is the most accurate method, but the application of this
method is limited to two-dimensional problems with simple geometries.

Algebraic method

This is one of the most common methods used in commercial codes appropriate for
several engineering applications. Clustering and stretching of grid elements using
algebraic method can be done by different functions such as: polynomial, trigonometric,
logarithmic, and geometric functions. Using the algebraic grid generation results in a
good control over the grid structure and is relatively simple to apply.

Differential equation method

The partial differential equations used to generate a grid can be of elliptic, parabolic, or
hyperbolic type. The most applied one is the elliptic type. In this case we want to have
control over the followings:
a) Grid point distribution on the boundaries,
b) The angle between the boundaries and the gridlines, and
c) The spacing between the gridlines.

Block-structured method

When the geometry is complex, it is very difficult to generate a single zone grid with
adequate control on the distribution of the mesh points using structured grids. There are
three main types of domain decomposition. These are patched zones, overlapped zones,
and overlaid zones. Patched zones have a common boundary line (see Fig. 15). The mesh
lines across the boundaries may be continuous or discontinuous49. In the second
technique, an overlap region exists between the zones. The extent of that region may be
from one up to several mesh points. In the third technique, which is also known as the
Chimera method. Smaller zones are defined on top of a base grid. Inter-zone data transfer
is accomplished by interpolation.

The application of block-structured grid with an algebraic grid generation for each block
is explained by the following example:

Unstructured grid

Unstructured grids have the advantage of generality in that they can be made to conform
to nearly any desired geometry. This generality, however, comes with a price. The grid
59
generation process is not completely automatic and may require considerable user
interaction to produce grids with acceptable degrees of local resolution while at the same
time having a minimum of element distortion. Unstructured grids require more
information to be stored and recovered than structured grids (e.g., the neighbor
connectivity list), and changing element types and sizes can increase numerical
approximation errors.

A popular type of unstructured grid consists of tetrahedral elements (Fig. 15). These grids
tend to be easier to generate than those composed of hexahedral elements, but they
generally have poorer numerical accuracy. For example, it is difficult to construct
approximations that maintain an accurate propagation of one-dimensional flow
disturbances because tetrahedral grid elements have no parallel faces.

In summary, the best choice for a grid system depends on several factors: convenience in
generation, memory requirements, numerical accuracy, flexibility to conform to complex
geometries and flexibility for localized regions of high or low resolution.

TEXT / REFERENCE BOOKS

1. Anderson, J.D., Computational Fluid Dynamics: The Basics with Applications, 2nd
Edition, McGraw Hill International
Editions, 2012.
2. Fletcher, C.A.J., Computational Techniques for Fluid Dynamics, 3rd Edition Springer
Verlag, 2001.
3. Versteeg, H.K. and Malalasekera, W., An Introduction to Computational Fluid
Dynamics: The Finite Volume Method,
2nd Edition, Pearson Education Ltd., 2007.
4. Chung T.J., Computational Fluid Dynamics, 2nd Edition, Cambridge University Press,
2002.

60