1- Introduction

1.1- Mathematical problems and engineering solving

Mathematical model can be broadly defined as a formulation or equation that

expresses the essential features of a physical system or process in mathematical
terms. In a very general sense, it can be represented as a functional relationship of
the form:

Fig.1 The engineering problem solving process

1
1.2- CFD:

Computational fluid dynamics or CFD is the analysis of systems involving fluid

flow, heat transfer and associated phenomena such as chemical reactions by means
of computer-based simulation.
The technique is very powerful and spans a wide range of industrial and non-
industrial application areas, some example are:

- Aerodynamics of aircraft and vehicles: lift and drag hydrodynamics of ships

- power plant: combustion in internal combustion engines and gas turbines
- turbo machine: flows inside rotating passages, diffusers etc

Advantages of CFD:

- Substantial reduction of lead times and costs of new designs

- Ability to study systems where controlled experiments are difficult or
impossible to perform (e.g. very large systems)
- Ability to study systems under hazardous conditions and beyond their
normal performance limits (e.g. safety studies and accident scenarios)

1.3-CFD codes:

CFD codes are structured around the numerical algorithms that can tackle fluid
flow problems. All codes contain three main elements: (i) a pre-processor, (ii) a
solver and (iii) a post-processor.

(i)-Pre-processor:
Pre-processing consists of the input of a flow problem to a CFD program, pre-
processing stage involve:

- Definition of the geometry of the region of interest: the

computational domain

- Grid generation – the sub-division of the domain into a number of smaller,

non-overlapping sub-domains: a grid (or mesh) of cells (or control volumes
or elements)
2
 - Selection of the physical and chemical phenomena that need to be modeled
- Definition of fluid properties
- Specification of appropriate boundary conditions at cells which coincide
with or touch the domain boundary
The solution to a flow problem (velocity, pressure, temperature etc.) is defined at
nodes inside each cell. The accuracy of a CFD solution is governed by the number
of cells in the grid.

(ii)-Solver:
There are three distinct streams of numerical solution techniques: finite difference,
finite element and finite volume methods.
The numerical algorithm consists of the following steps:

- Integration of the governing equations of fluid flow over all the (finite)
control volumes of the domain.
- Discretisation – conversion of the resulting integral equations into a system
of algebraic equations.
- Solution of the algebraic equations.

(iii)-Post-processor

As in pre-processing, a huge amount of development work has recently taken place

in the post-processing field. Due to the increased popularity of engineering
workstations, many of which have outstanding graphics capabilities, the leading
CFD packages are now equipped with versatile data visualization tools. These
include:

- Domain geometry and grid display

- Vector plots
- Line and shaded contour plots
- 2D and 3D surface plots
- Particle tracking
- View manipulation (translation, rotation, scaling etc.)
- Colour PostScript output

3
Three mathematical concepts are useful in determining the success or otherwise of
such algorithm: convergence, consistency and stability.

Convergence: is the property of numerical method to produce solution which

approaches the exact solution as grid spacing reduced to zero.

Consistency: the numerical schemes produces system of algebraic equation

equivalent to the original governing equation as the grid spacing tends to zero.

Stability: is associated with damping of errors as the numerical method proceeds.

1.4 Conservation laws of fluid motion and boundary condition

The governing equations of fluid flow represent mathematical statements of the

conservation laws of physics:

- The mass of a 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)

We consider such a small element of fluid with sides δx, δy and δz

Fig.2 Fluid element of conservation laws

4
The six faces are labelled N, S, E, W, T and B, which stands for North, South, East,
West, Top and Bottom
All fluid properties are functions of space and time so we would strictly need to
write ρ(x, y, z, t), p(x, y, z, t), T(x, y, z, t) and u(x, y, z, t) for the density, pressure,
temperature and the velocity vector respectively

1.4.1 Mass conservation:

The first step in the derivation of the mass conservation equation is to write down a
mass balance for the fluid element:

The rate of increase of mass in the fluid element is:

The net rate of flow of mass into the element across its boundaries is given by:

Flows which are directed into the element produce an increase of mass in the
element and get a positive sign and those flows that are leaving the element are
given a negative sign.

5
 Fig.3 mass flow in and out of fluid element

The rate of increase of mass inside the element is now equated to the net rate of
flow of mass into the element across its faces. All terms of the resulting mass
balance are arranged on the left hand side of the equals sign and the expression is
divided by the element volume. This yield:

or in more compact vector notation:

For an incompressible fluid:

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

6
The momentum and energy conservation laws make statements regarding changes
of properties of a fluid particle. This is termed the Lagrangian approach. Each

time t. Let the value of a property per unit mass be denoted by 𝜙. The total or
property of such a particle is a function of the position (x, y, z) of the particle and

substantive derivative of φ with respect to time following a fluid particle, written

as Dφ/Dt, is:

As in the case of the mass conservation equation, we are interested in developing

equations for rates of change per unit volume. of property φ per unit volume for a
fluid particle is given by the product of D𝜙/Dt and density ρ, hence:

The mass conservation equation contains the mass per unit volume (i.e. the density
ρ) as the conserved quantity. The sum of the rate of change of density in time and
the convective term in the mass conservation equation for a fluid element is

The generalization of these terms for an arbitrary conserved property is:

7
The previous equation expresses the rate of change in time of 𝜙per unit volume
plus the net flow of 𝜙 out of the fluid element per unit volume. It is now rewritten
to illustrate its relationship with the substantive derivative of 𝜙:

1.4.2 Conservation of linear momentum:

Newton’s second law states that the rate of change of momentum of a fluid particle
equals the sum of the forces on the particle:

We distinguish two types of forces on fluid particles:

Surface forces:
– Pressure forces
– Viscous forces
Body forces:
– Centrifugal force
– Electromagnetic force
– Gravity force

8
Fig.4 Stress components on three faces of fluid element fig.5 Stress components in the x-
direction

The state of stress of a fluid element is defined in terms of the pressure and the
nine viscous stress components shown in the figure 4,5The pressure, a normal
stress, is denoted by p. Viscous stresses are denoted by τ

The net force in the x-direction is the sum of the force components acting in that
direction on the fluid element.

The total force per unit volume on the fluid due to these surface stresses
is equal to the sum of previous equations divided by the volume δxδyδz

9
Total force in the x-direction on the element due to surface stresses plus the rate of
increase of x-momentum due to sources:

It is not too difficult to verify that the y-component of the momentum equation is
given by:

And the z-component of the momentum equation by:

The sign associated with the pressure is opposite to that associated with
the normal viscous stress, because the usual sign convention takes a
tensile stress to be the positive normal stress so that the pressure, which
is by definition a compressive normal stress, has a minus sign.
the source
Terms SMx, SMy and SMz include contributions due to body forces only. For
example, the body force due to gravity would be modeled by SMx = 0, SMy = 0 and
SMz = −ρg.

1.4.3 Navier stocks equation of Newtonian fluid

10
In a Newtonian fluid the viscous stresses are proportional to the rates
of deformation. The three-dimensional 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, and the second viscosity, λ, to
relate stresses to the volumetric deformation. The nine viscous stress components,
of which six are independent, are:

Substitution of the above shear stresses into the Total force in the x-direction on
the element due to surface stresses plus the rate of increase of x-momentum due to
sources yields the so-called Navier–Stokes equations, named after the two
nineteenth-century scientists who derived them independently:

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

the Navier–Stokes equations can be written in the most useful form for the
development of the finite volume method:

12
1.4.4 Conservation of energy:

The energy equation is derived from the first law of thermodynamics, which states
that the rate of change of energy of a fluid particle is equal to the rate of heat
addition to the fluid particle plus the rate of work done on the particle:

As before, we will be deriving an equation for the rate of increase of energy of a

fluid particle per unit volume, which is given by:

Work done by surface forces:

The rate of work done on the fluid particle in the element by a surface force is
equal to the product of the force and velocity component in the direction of the
force.
The net rate of work done by these surface forces acting in the x-direction is given
By:

Surface stress components in the y- and z-direction also do work on the fluid
particle. A repetition of the above process gives the additional rates of work done
on the fluid particle due to the work done by these surface forces:

13
The total rate of work done per unit volume on the fluid particle by all the surface
forces are given by the sum of the previous equations divided by the volume
δxδyδz. The terms containing pressure can be collected together and written more
compactly in vector form:

This yields the following total rate of work done on the fluid particle by surface
Stresses:

Energy flux due to heat conduction:

The heat flux vector q has three components: qx, qy and qz

14
 Fig.6 Components of the heat flux vector

The net rate of heat transfer to the fluid particle due to heat flow in the x-direction
is given by the difference between the rate of heat input across face W and the rate
of heat loss across face E:

Similarly, the net rates of heat transfer to the fluid due to heat flows in the y- and z-
direction are:

The total rate of heat added to the fluid particle per unit volume due to heat flow
across its boundaries is the sum of the previous equations divided by the volume
δxδyδz:

Fourier’s law of heat conduction relates the heat flux to the local temperature
gradient. So

This can be written in vector form as follows:

15
The final form of the rate of heat addition to the fluid particle due to heat
conduction across element boundaries:

The energy equation is:

1.5 Differential and integral forms of the general transport equations:

We introduce a general variable 𝜙 the conservative form of all fluid flow

equations, including equations for scalar quantities such as temperature and
pollutant concentration etc., can usefully be written in the following form:

The previous equation is called transport equation for property 𝜙. It clearly

highlights the various transport processes: the rate of change term and the
16
convective term on the left hand side and the diffusive term (Γ = diffusion
coefficient) and the source term respectively on the right hand side.

For each of the five PDEs for mass, momentum and energy conservation, the key
step of the finite volume method is the integration of the previous equation over a
three-dimensional control volume (CV):

The volume integrals in the second term on the left hand side, the convective term,
and in the first term on the right hand side, the diffusive term, are rewritten as
integrals over the entire bounding surface of the control volume by using Gauss’s
divergence theorem. For a vector a, this theorem states

The physical interpretation of n.a is the component of vector a in the direction of

the vector n normal to surface element dA. Thus the integral of the divergence of a
vector a over a volume is equal to the component of a in the direction normal to
the surface which bounds the volume summed (integrated) over the entire
bounding surface A.

Applying Gauss’s divergence theorem, the previous equation can be written as

follows:

The order of integration and differentiation has been changed in the first term on

the rate of change of the total amount of fluid property 𝜙 in the control volume.
the left hand side of equation to illustrate its physical meaning. This term signifies

The product n.ρ 𝜙 u expresses the flux component of property 𝜙 due to fluid flow
along the outward normal vector n, so the second term on the left hand side of

17
property 𝜙 of the fluid element due to convection.
equation, the convective term, therefore is the net rate of decrease of fluid

the net rate of increase of fluid property 𝜙 of the fluid element due to
The diffusive term, is thus associated with a flux into the element and represents

increase of property 𝜙 as a result of sources inside the fluid element. In words,

diffusion. The final term on the right hand side of this equation gives the rate of

the previous relationship can be expressed as follows:

In steady state problems the rate of change term of transport equation is equal to
zero. This leads to the integrated form of the steady transport equation:

In time-dependent problems it is also necessary to integrate with respect to time t

over a small interval Δt from, say, t until t + Δt. This yields the most general
integrated form of the transport equation:

1.6.1 Classification of physical behaviors

The initial and boundary conditions that are needed in conjunction with the
equations to construct a well-posed Mathematical model of a fluid flow. First we
distinguish two principal categories of physical behavior:
18
 - Equilibrium problems
- Marching problems

Equilibrium problems:

The problems in the first category are steady state situations, e.g. the steady state
distribution of temperature in a rod of solid material or the equilibrium stress
distribution of a solid object under a given applied load, as well as many steady
fluid flows. These and many other steady state problems are governed by elliptic
equations.
Unique solution to this and all elliptic problems can be obtained by specifying
conditions on the dependent variable (here the temperature or its normal derivative
the heat flux) on all the boundaries of the solution domain. Problems requiring data
over the entire boundary are called boundary-value problems.

Fig.7 Steady state temperature distribution of an insulated rod

Marching problems:
Transient heat transfer, all unsteady flows and wave phenomena are examples of
problems in the second category, the marching or propagation problems. These
problems are governed by parabolic or hyperbolic equations. However, not all
marching problems are unsteady.
Parabolic equations describe time-dependent problems, which involve significant
amounts of diffusion. Examples are unsteady viscous flows and unsteady heat
conduction. This type of problem is termed an initial–boundary-value problem.

19
 Fig.8 Transient distribution of temperature in an insulated rod

Hyperbolic equations dominate the analysis of vibration problems. In general

they appear in time-dependent processes with negligible amounts of energy
dissipation. Hyperbolic problems are also initial–boundary-value problems.

Fig.9 Vibrations of a string under tension

1.6.2 Classification method for simple PDEs:

A practical method of classifying PDEs is developed for a general second order

PDE in two co-ordinates x and y. Consider

20
At first we shall assume that the equation is linear and a, b, c, d, e, f and g are
constants.
Simple we have solutions occur if the next characteristic equation below has two
real roots:

The existence or otherwise of roots of the characteristic equation depends on the

value of discriminate (b2 − 4ac). Table .1 outlines the three cases.

Table.1 Classification of linear second-order PDEs

21