UNIT - II

CONDUCTION AND ADVECTION

2.1 TWO DIMENSIONAL STEADY HEAT CONDUCTION EQUATION

The methodology used in deriving discretised equations in the 1-D case can be easily extended to 2-D problems. Consider
a stationary fluid where heat is transferred by thermal conduction. Moreover assume that the temperature gradients in the
fluid are changing as a function of time. Now heat governing equation in 2-D form is
. ∂  ∂T ∂  ∂T
ρq + k + k = 0 …(2.1)
∂x ∂x  ∂y ∂y
where T is temperature in K.
.
Where q, heat generation per unit volume. In addition to the east (E) and west (W) neighbours in a general grid
(Fig. 2.1) node P now also has north (N) and south (S) neighbours.

Fig. 2.1 : 2-D Grid

When the above equation is formally integrated over the control volume ∆V, we get

⌠ ∂ k ∂T dx ⋅ dy + ⌠ ∂ k ∂T dx ⋅ dy + ⌠ ρq. dV = 0

⌡ ∂x ∂x ⌡ ∂y ∂y ⌡
∆v ∆v ∆v

Taking area of cell face Ae = Aw = ∆y and An = As = ∆x we obtain

keAe
∂T – k A ∂T + k A ∂T – k A ∂T + q. ⋅∆V = 0 …(2.2)
∂xe w w∂xw n n∂xn s s∂xs
This equation represents the balance of the heat conduction equation in a control volume and the fluxes through its cell
faces. Using the approximation introduced we can write expressions for flux through control volume faces:

Flux across the west face =

kA ∂T = k A TP – TW
 ∂xw w w δxWP

Flux across the east face =

kA ∂T = k A TE – TP
 ∂xe e e δxPE

Flux across the south face =

kA ∂T = k A TP – TS
 ∂xs s s δySP

(2.1)
COMPUTATIONAL FLUID DYNAMICS (BE MECH.) (2.2) CONDUCTION AND ADVECTION

Flux across the north face =

kA ∂T = k A TN – TP
 ∂xn n n δyPN
By substituting the above expressions into equation (2.2), we get
TE – TP TP – TW TN – TP TP – TS .
keAe – kwAw + knAn – ksAs + q ⋅∆V = 0 …(2.3)
δxPE δxWP δyPN δySP

Above equation (2.3) is now written in the general discretised equation form
.
apTp = awTw +aETE + aSTS + aNTN + q …(2.4)

where

aw aE as aN ap
kwAw keAe ksAs knAn aw + aE + as
δxwp δxPE δysp δyPN .
+ aN – q

We obtain the distribution of temperature T in a given 2-D situation by writing discretised equations of the form (2.3) at
each grid node of the sub-divided domain. At the boundaries where the temperatures or fluxes are known the discretised
equations are modified to incorporate boundaries conditions.

2.2 IMPLICIT METHOD FOR 2-D UNSTEADY HEAT CONDUCTION EQUATION

The fully implicit method is recommended for general purpose CFD computations on the ground of its superior stability.
Unsteady Heat Conduction Equation in 2-D
∂T ∂  ∂T ∂  ∂T
ρc + k + k = 0 …(2.5)
∂t ∂x ∂x  ∂y  ∂y

where T is temperature in K and t is time.

A 2-D control volume (Fig. 2.1) is considered for Discretisation. The resulting equation
0 0
apTp = awTw + aETE + aSTS + aNTN + apTp

0 0 ∆V
Where aw + aE + aS + aN + ap = ap and ap = ρc
∆t

The neighbouring of P are aW, aE in 1-D problems and aW, aE, aS, an in 2-D.

The following values for the volume and cell face areas.

1-D 2-D

∆V ∆x ∆x ⋅∆y

aw = ae 1 ∆y

an = as … ∆x

In the fully implicit discretisation approach discussed above for 2-D heat transfer problems, the term arising from temporal
0 0 0
discretisation appears as (i) the contribution of ap to the central coefficient ap and (ii) the contribution of apTpas an
additional source term on the right hand side. The other terms are unaltered and are the same as in the discretised
equations for steady state problems. Using this as a basis the discretised equations for transient heat governing equations
are also simple to obtain. The unsteady transport of temperature T is given by
∂
∂t
–
⎡ 0
(ρT) + div (ρVT) = div ( grad T) + ap …(2.6)
COMPUTATIONAL FLUID DYNAMICS (BE MECH.) (2.3) CONDUCTION AND ADVECTION

2.3 DIRICHLET AND NEUMANN BOUNDARY CONDITIONS

2
• Let us consider an elliptic equation (b – 4ac < 0) in two independent variables x and y. As we know that the
characteristic curves for an elliptic equation are imaginary and the methodology associated with the of characteristics is
of no use. For elliptic equations there are no limited regions of influence or domains of dependence. For example
consider point P in the xy plane as shown in Fig. (2.2); assuming that the domain of the problem is defined as the
rectangle abcd.

Fig. 2.2
• Now we introduce a disturbance at point P. The major mathematical characteristic of elliptic equations is that this
disturbance is felt everywhere throughout the domain. Further point P influences all the points in the domain, then in
turn the solution at point P is influenced by the entire closed boundary abcd.
• Therefore the solution at point P must be carried out simultaneously with the solution at all points in the domain. For
this reason, problems involving elliptic equations are frequently called jury (group) problems, because the solution
within the domain depends on the total boundary domain, boundary conditions must be applied over the entire
boundary abcd. These type of boundary conditions can take the following forms:
1. A specification (detailed description) of the dependent variables u and v along the boundary. This type of
boundary condition is called the Dirichlet Conditions.
∂u
2. A specification of derivatives of the dependent variables, such as , along the boundary. This type of boundary
∂x
condition is called Neumann Boundary Condition.
3. A mix of both is called Dirichlet and Neumann Boundary Conditions.
• In a sub-sonic flow, disturbances (which travel at the speed of sound or even faster) can physically find their way
upstream for as far as they want – theoretically; a finite disturbance in an inviscid flow will propagate to infinity in all
directions.
• In reality an incompressible flow is a limiting case of a subsonic flow wherein the Mach number (M = V/a) goes to zero.
In incompressible flow, the compressibility is zero, hence the speed of sound is infinite. An incompressible inviscid flow
is governed by elliptic equations and such flows are the focus point of elliptic behaviour.

2.4 SOLUTION BY EXPLICIT AND ALTERNATING DIRECTION IMPLICIT METHOD

Numerical tools for solving Computational Fluid Dynamics (CFD) problems use two approaches namely: An explicit
approach and implicit approach. It is appropriate to define these two general approaches now. Let us again consider 2-D
unsteady heat conduction equation, with variables x, y, t.
2 2
∂T  ∂ T  ∂ T = 0
ρc + k 2 +k k 2 …(2.7)
∂t  ∂x   ∂y 
2 2
∂T ∂T ∂T
Here we choose to represent with a forward difference and 2 , 2 with a central difference, leading to the particular
∂t ∂x ∂y
form of the difference equation,
k
where α = .
ρC
COMPUTATIONAL FLUID DYNAMICS (BE MECH.) (2.4) CONDUCTION AND ADVECTION

n+1 n n n n
Ti – Ti α (Ti + 1 – 2Ti + Ti – 1)
= 2 …(2.8)
∆t (∆x)
n+1 n n n n
Tj – Tj α (Tj + 1 – 2Tj + Tj – 1)
And = 2 …(2.9)
∆t (∆y)
This equation (2.9) can be written as
n+1 n ∆t n n n
Ti = Ti + α 2 (T – 2Ti + Ti – 1)
(∆x) i + 1
n+1 n ∆t n n n
And Tj = Tj + α 2 (Tj + 1 – 2Tj + Tj – 1)
(∆y)
The marching variable here is time t. To be more specific, consider the finite difference grid sketched in Fig. (2.3). Here i is
the running index in x-direction and j in y-direction.
Time marching means that T at all grid points at time level n + 1 are calculated from the known values at time level n. Then
the same procedure is used to calculate T at all grid points at time level n + 2, using the known values at level n + 1. In this
fashion the solution is progressively obtained by marching in steps of time.

Fig. 2.3 : Illustration of Time-Marching

What we have just represented in the above paragraph is an example of an explicit approach. In an explicit approach each
difference equation contains only one unknown and therefore can be solved explicitly in a straight-forward manner.

Now for equation (2.8), writing spatial difference on the right hand side in terms of average properties between time level n
and n + 1.
1 n+1 n 1 n+1 1 n+1 n
n+1
Ti
n
– Ti 2 (
Ti + 1 + Ti + 1 +
2 ) n
–2Ti – + 2Ti +(2 )
Ti – 1 + Ti – 1 ( )
= a 2 …(2.10)
∆t (∆x)

This type of differencing employed above is called Crank-Nicholson form. In CFD, Crank-Nicholson form is used frequently
n+1
for finite –differencing solutions of the boundary layer equations. The unknown Ti is not only expressed in terms of the
n n n
known quantities at time level n, namely Ti + 1, Ti and Ti – 1 but also in terms of other unknown quantities at time level n + 1,
n+1 n+1 n+1
namely Ti + 1 , Ti – 1 . At a given grid point (T) does not stand alone; it cannot by itself result in a solution for Ti for all i
can be solved simultaneously. This is an example of an implicit approach. An implicit approach is one where the unknowns
must be obtained by means of a simultaneous solution of the difference equations applied at all the grid points arrayed at
a given time level. Because of this need to solve large systems of simultaneous algebraic equations, implicit methods are
usually involved with the manipulations of large matrices. It is easy to understand that the implicit approach involves a
more complex set of calculations than the explicit approach.
COMPUTATIONAL FLUID DYNAMICS (BE MECH.) (2.5) CONDUCTION AND ADVECTION

Simplifying the nomenclature by denoting the following quantities by A, B and Ni;

a ⋅∆t
A = 2,
2 (∆x)
α⋅∆t
B = 1+ 2,
(∆x)
n α⋅∆t n n n
Ni = – Ti – (
2 T
2(∆x) i+1
– 2Ti + Ti–1 )
We can write equation (2.7) in the form
n+1 n+1 n+1
ATi–1 – BTi + ATi+1 = Ni …(2.11)
At time level n we apply equation (2.11) at grid point 2
AT1 – BT2 + AT3 = N2 …(2.12)
It is easy to remember that T1, T2 and T3 represent three values at time level n + 1 and N2 is a known number. Because of
the stipulated boundary conditions at grid point 1, T1 is a known number. Now equation (2.12) can be written as
–BT2 + AT3 = N2 – AT1
Substituting
'
N2 – AT1 = N2
'
hence –BT2 + AT3 = N2
Now At time level n, we apply equation (2.12)
At grid point 3
AT2 – BT3 + AT4 = N3
At grid point 4
AT3 – BT4 – AT5 = N4
At grid point 5
AT4 – BT5 + AT6 = N5
At grid point 6
AT5 – BT6 + AT5 = N6

Since grid point 7 is on a boundary, T7 is known from the stipulated boundary condition. For grid point 7,
'
AT5 – BT6 + AT5 – BT6 – AT7 = N6

This system of equations can be written in Matrix form as follows:

'

A  T  N 
–B A 0 0 0 T2 N2

0  T =N 
–B A 0 0 3 3

A –B A 0 …(2.13)
0 A  T  N 
4 4

0 A –B
0 –B   T  N 
5 5

0 0 A 6 6

• The solution of the system of equations denoted by (2.13) involves the manipulation of the tridiagonal arrangement;
such solutions are usually obtained using Thoma’s algorithm. An implicit solution would therefore demand a
simultaneous solution of a large system of non-linear equations an exceptionally difficult task. This is a very big
disadvantage of an implicit approach.
COMPUTATIONAL FLUID DYNAMICS (BE MECH.) (2.6) CONDUCTION AND ADVECTION

• It is obvious, that an implicit approach is more involved (tedious) than an explicit approach. With the complexity of the
implicit approach relative to explicit approach, question arises what is the need of dealing with implicit approach.
Actually, we have yet to mention the most important difference between the explicit and implicit approaches. It is
worth to mention that the increments ∆x and ∆t appeared in all the difference equations. For the explicit approach,
once ∆x is chosen, then ∆t is not an independent, rather ∆t is restricted to be equal to or less than a certain value
prescribed by a stability criterion.
• If ∆t is taken larger than the limit imposed by the stability criterion, the time marching procedure will quickly go
unstable and computer program will quickly shut down as numbers going to infinity or taking the square root of a
negative number. In many cases, ∆t must be very small to maintain the stability; this can result in long computer
running times to make calculations over a given interval of time.
• On the other hand, there are no such stability restrictions on an implicit approach. For most implicit methods, stability
can be maintained over much larger values of ∆t than for a corresponding explicit method. Hence, for an implicit
method, considerable fewer steps are required to cover a given interval in time compared to an explicit method.
Therefore, for some applications, even though the implicit methods require more computations time as per time step
due to its relative complexity, the fact that considerable fewer time steps are required to cover a given interval of time
actually can result in a shorter run time on the computer compared to an explicit approach.

2.5 ROBBIN BOUNDARY CONDITIONS

• The stability of least energy solution in the Robin boundary condition case is quite complicated. We state the following
result which deals with one-dimensional case only
p
2 a
at = ∈ axx – a + q , 0 < x < 1, t > 0,
ξ
1

–s ⌠ ar dx
τξt = – ξ + ξ ⌡
0

a > 0, ∈ ax (0, t) + λa (0, t) = ∈ ax (1, t) + λa (1, t) = 0

hx(0, t) = hx (1, t) = 0
Assume that r = 2, 1 p ≤ 3 or r = p + 1, 1 < p < +∞. en for each λ ∈ 0, 1) the least energy solution is stable for τ < τ1 and
unstable for τ > τ1. At τ1, There is a Hop fabifurcation.
The main idea of the proof is similar to that of Theorem, the study an NLEP on a half line with Robin boundary condition:
∞

⌠ w φ
⌡ x0
p–1 0 p
φ” – φ + pwx0 φ – γ(p – 1) ∞ wx0 = αφ, 0 < y < + ∞ φ'(0) – λφ (0) = 0
⌠ w2
⌡ x 0

p–1
where, wx0 = w(y – x0) with w’ (–x0) = λw(–x0). Let Lx0 (φ) = φ” – φ + pwx0 φ.

Then we need to show that ⌠

–1
⌡ wx0[Lx0 (wx0)] > 0.
0

By some lengthy computations, we can show that the function ⌡

⌠ w [L–1 (w )] is an increasing function in x when p < 3
x0 x 0 x0 0
0

and a constant when p = 3 and an decreasing function when p > 3.

COMPUTATIONAL FLUID DYNAMICS (BE MECH.) (2.7) CONDUCTION AND ADVECTION

2.6 STABILITY CRITERIA

Most hyperbolic problems involve the transport of fluid properties. In the equations of motion, the term describing the
transport process is often called convection or advection.
e.g. the 1-D equation of motion is
du ∂u ∂u
= +u
dt ∂t ∂x
1 ∂p 2
= –ρ + v∇ u.
∂x
∂u
Here the advection term u term is non-linear.
∂x
nd
The classical 2 -order hyperbolic wave equation is
2 2
∂u 2 ∂ u
2 = c 2 .
∂t ∂x
The equation descibes wave propagation at a speed of c in two directions.
st
The 1 -order equation that has properties similar
∂u ∂u
+c = 0, c>0.
∂t ∂x
For a pure initial value problem with initial condition
u(x, 0) = F(x), –∞ < x < ∞,
the exact solution is u(x,t) = F(x-ct), which we have obtained earlier using the method of characteristics. We know that the
solution represents a signal propagating at speed c.

2.7 ALTERNATING DIRECTION IMPLICIT TECHNIQUE (ADI)

Let us consider 1-D heat conduction equation in unsteady state.
2
dT ∂T
= α 2 …(2.14)
dt ∂x
• There exists only one independent variable x, in this equation. As long as we are dealing with linear equations, the
implicit solutions using the Crank-Nicholson scheme are directly obtained using Thoma’s algorithm. When solving a
nonlinear problem by means of an implicit scheme, the matter of linearizing the difference equations is of utmost
importance so that Thoma’s algorithm can be used to expedite the calculations.
• The main focus area here, is concerned with the other aspects that destroy the tridiagonal nature of the difference
equation, namely, multi-dimensionality involving more than one variable in addition to the marching variable, time t.
2 2
dT ∂ T ∂ T
= α 2+ 2 …(2.15)
dt ∂x ∂y 
Using Crank-Nicholson Scheme, finite difference form of the above equation can be written as
1 n+1 n 1 n+1 n
n+1
T –T
i,j 2 (nTi + 1,j + Ti + 1, j) + (–2Ti, j – 2Ti , j )
i,j 2 1 n+1 n
= a 2 + (Ti–1,j + Ti–1,j )
∆t (∆x) 2
1 n+1 n 1 n+1 n 1 n+1 n

2 ( Ti,j + 1 + Ti, j + 1) + (–2 Ti, j – 2Ti,j) + (Ti,j – 1 + Ti,j – 1)

2 2
+a 2 …(2.16)
(∆y)
Above equation (2.16) is equivalent in xy space to the 1-D form given by equation (2.14), which on reducing the
n+1 n+1 n+1 n+1 n+1
tridiagonal form gives five unknowns, namely, Ti+1,j, Ti,j , Ti–1,j, Ti,j + 1 and Ti,j–1 where the last two unknowns prevent a
COMPUTATIONAL FLUID DYNAMICS (BE MECH.) (2.8) CONDUCTION AND ADVECTION

tridiagonal form. Hence Thoma’s algorithm cannot be used. As a result there is need of developing a scheme that will allow
equation (2.16) to be solved by diagonal form. Such a scheme is known as Alternating Direction Implicit Technique (ADI).
The equation (2.16) is solved by means of a marching technique; i.e. T(t + ∆t) is obtained in some way from the known
values of T(t). Let us achieve the solution of T (t + ∆t) in a two step process, where intermediate values of T are found at an
 ∆t . For equation (2.15)
intermediate time, t +
 2
n + 1/2 n n + 1/2 n + 1/2 n+1/2
n n n
Ti,j – Ti,j Ti+1,j – 2Ti,j + Ti–1,j Ti,j+1 – 2Ti,j + Ti,j–1
= α 2 +α 2 …(2.17)
∆t (∆x) (∆y)
2
Equation (2.17) reduces to diagonal form
n + 1/2 n+1/2 n+1/2
ATi–1,j – BTi,j + ATi+1,j = Ni …(2.18)

α⋅∆t α⋅∆t
Where, A = 2 , B = 1 + 2 ,
2(∆x) 2(∆x)
n α⋅∆t n n n
Ni = – Ti,j – 2 (T – 2Ti,j + Ti,j – 1)
2(∆y) 1,j+1
n+1/2
Equation (2.17) yields a solution for Ti,j for all i, keeping j fixed using Thoma’s algorithm. At a fixed value of j, we sweep
in the x–direction. If there are N grid points in the x–direction, then we sweep from i = 1 to N. This sweep utilizes Thoma’s
algorithm once. This calculation is then repeated at the next row of grid points designated by j + 1. It means replace j in
n+1/2
equation (2.17) by j + 1 and solve for Ti,j+1 for all values of i from 1 to N using Thoma’s algorithm. If there are M grid
points in the y-direction, this process is repeated M times.

Fig. 2.4 : Sweeping in the x-Direction to Obtain T at (t + ∆t/2)

n+1/2
At the end of this step, the values of T at the intermediate time t + ∆t/2 are known at all grid points (i, j); i.e. Ti,j is known
at all (i, j).
The second step of the ADI scheme takes the solution to time t + ∆t, using the known values at time t + ∆t/2. For this
second step, the spatial derivatives in heat governing equation are replaced with central differences, where the y derivative
is treated implicitly.
n+1 n+1/2 n+1/2 n+1/2 n+1/2 n+1 n+1 n+1
Ti,j – Ti,j Ti+1,j – 2Ti,j + Ti–1,j Ti,j+1 – 2Ti,j + Ti,j – 1
= α 2 +α 2 …(2.19)
∆t (∆x) (∆y)
2
COMPUTATIONAL FLUID DYNAMICS (BE MECH.) (2.9) CONDUCTION AND ADVECTION

Equation (2.19) reduces to the tridiagonal form

n+1 n+1 n+1
CTi,j+1 – DTi,j + CTi,j–1 = Lj …(2.20)

where
α⋅∆t α⋅∆t
C = 2 , D = 1 + 2,
2(∆y) 2(∆y)
n+1/2 α⋅∆t n+1/2 n+12 n+1/2
Li = – Ti,j –
2(∆y) (
2 Ti+1,j – 2Ti,j )
+ Ti–1,j

1
n+ n+1/2
2
We know that T is known at all grid points from the first step. Equation (3.19) yields a solution for Ti,j for all j, keeping
n+1
i fixed, using Thoma’s algorithm. At a fixed value of i, we sweep in the y-direction, using equation (3.16) to solve for Ti,j for
all values of j, where j goes from 1 to M.
At the end of this two-step process, the dependent variable T has been marched a value ∆t in the direction of t. Although,
there are two independent spatial variables x and y in addition to the marching variable t, this marching scheme involves
only tridiagonal forms and the solution has been achieved by the repeated application of Thoma’s algorithm. Because the
scheme involves two steps, one in which the difference equation is implicit in x and the other in which the difference
equation is implicit in y, the name of the scheme Alternating Direction Implicit (ADI).
2 2 2
The ADI scheme is second-order accurate in t, x and y and truncation error is O [(∆t) , (∆x) , (∆y) ]. This scheme is very
useful in many fluid flow problems. It is particularly useful for the solution of problems described by parabolic partial
differential equations.

2.8 APPROACH FOR IRREGULAR BOUNDARY FOR 2-D HEAT CONDUCTION PROBLEMS
• The governing equations govern the flow of a fluid in a flow field. However the flow fields are quite different for
different flows, although the governing equations are the same. Then question arises from where these differences
enter. The answer is through boundary conditions, which are different in each case. The boundary conditions and
sometimes the initial conditions, dictate the particular solutions to be obtained from the governing equations.
• Hence, once we decide the governing flow equations, then the real driver for any particular solution is the boundary
conditions.
• This has particular significance in CFD; any numerical solution of the governing flow equations must be made to see a
strong and compelling numerical representation of the boundary conditions. For flow through a duct of fixed shape,
there are boundary conditions which pertain to the inflow and outflow boundaries, such as at the inlet and exit of the
duct.
• If the problem involves an aerodynamic body immersed in a known free stream, then the boundary conditions applied
at a distance infinitely far upstream and downstream of the body are simply that of the given free stream conditions.
These boundary conditions are called physical boundary conditions imposed by nature. In CFD, we have an additional
concern, namely the proper numerical implementation of these physical boundary conditions.
• In the same sense as the real field is driven by the numerical boundary conditions, the computed flow field is driven by
numerical formulation designed to simulate these boundary conditions.
• A general second order Partial Differential Equation (PDE) in two co-ordinates x and y.
2 2
∂T ∂ T ∂ T
= α 2+ 2 …(2.21)
∂t ∂x ∂y 
2
• System of second order or mixtures of first- and second-order PDEs could be classified as (b – 4ac > 0) the equation is
2 2
hyperbolic; (b – 4ac = 0) the equation is parabolic; (b – 4ac < 0) the equation is elliptic.
• The complicated mixture of elliptic, parabolic and hyperbolic behaviours has implications for the way in which
boundary conditions enter into a flow problem, in particular at locations where flows are bounded by fluid boundaries.
COMPUTATIONAL FLUID DYNAMICS (BE MECH.) (2.10) CONDUCTION AND ADVECTION

Initial Conditions for Unsteady Flows:

Everywhere in the solution region ρ, V and T must be given at time t = 0.
Boundary Conditions for Unsteady and Steady Flows:
• On solid walls u = uw(no slip condition); T = Tw (wall fixed temperature).
• On fluid boundaries ρ, V and T must be given at known as a function of position.
• For an incompressible viscous flow there are no conditions on the density, but all other conditions apply without
modification. Open (far field) boundary conditions give the most serious problems for the designer of general-purpose
CFD codes. Subsonic inviscid compressible flow equations have fewer inlet conditions (normally onlyρ and V are
specified) than viscous flow equations and only one outlet condition (generally pressure) exists.
• Outlet boundary conditions may be used in conjunction with the inlet boundary conditions. If the location of the outlet
is selected far away from geometrical disturbances the flow eventually reaches fully developed state where no changes
occur in the flow direction, in such a region we can safely assume that the gradients of all variables (except pressure)
are zero in flow direction. This gives us opportunity to locate the outlet surface perpendicular to the flow direction and
take gradients in the direction normal to the outlet surface equal to zero.
• Wall Boundary Conditions: The wall is the most common boundary encountered in fluid and heat transfer problems.
The no-lip condition (velocity parallel to wall u = 0, velocity perpendicular to wall v = 0) is the appropriate condition
for the velocity components at the wall. Immediately adjacent to the wall, we have an extremely thin viscous sub-layer
followed by the buffer layer and the turbulent core. Normally, we employ the ’wall function’ to represent the effect of
wall boundary.
• At the wall, the physically proper boundary condition for an inviscid flow is that the flow be tangent to the wall. This is
the only boundary condition at the wall; the other flow properties at the wall must be obtained/considered as part of
the solution.
• Moving Wall: Wall movement in the x-direction is felt by the fluid by a change in the wall shear stress. Its value is
adjusted by replacing velocity u by relative velocity u – uwall.
• Considerable effort has been made by the research community to develop robust methods for general purpose CFD
computations involving complex geometry. But no particular method is valid for all the cases.

2.9 SOLUTION OF TWO DIMENSIONAL STEADY AND UNSTEADY HEAT ADVECTION EQUATION
USING FINITE VOLUME METHOD WITH DIRICHLET BC
Heat transfer is an important problem in many disciplines including science, physics and engineering. Various numerical
techniques have been developed in order to solve and simulate real-world heat transfer problems, among which the FVM.
In this project, the FVM is implemented to compute the temperature distribution over a plate, finned at its bottom edge
and subject to prescribed boundary temperatures. A home-made code was developed to implement the numerical
technique. The system of equations arising from the discretization of the governing equation is solved iteratively by means
of the CGM.

The governing transport equation for a two-dimensional steady-state diffusion problem is given by:
∂  ∂Φ ∂  ∂Φ
Γ + Γ + SΦ = 0
∂x  ∂x  ∂y  ∂y 

where x; y are the space dimensions Γ, is the diffusion coefficient, Φ is the diffusive flux and SΦ is a source term. For the
special case of steady-state heat conduction without volumetric heat generation and the diffusive is the temperature T,
has the following form:
∂  ∂T ∂  ∂T
k + k = 0
∂x  ∂x  ∂y  ∂y

where k is the material thermal conductivity.

COMPUTATIONAL FLUID DYNAMICS (BE MECH.) (2.11) CONDUCTION AND ADVECTION

The geometry of the problem is shown in the following Fig. 2.5:

The following assumptions are made to ensure the two-dimensionality of the problem:

• The plate has a uniform thickness. There is no heat transfer through the thickness.

• The fundamental step in the FVM is the integration of the transport equation over a control volume. The domain is
discretized into cells or control volumes as shown in Fig. 2.5 below to approximate the differential equation. The values
of the desired properties are evaluated at the center of the control volume.

Fig. 2.5

The integration of over a control volume yields:

⌠ ∂ Γ ∂Φ dx⋅dy + ⌠ ∂ Γ ∂Φ dx⋅dy + ⌠ S dV = 0

⌡ ∂x  ∂x  ⌡ ∂y  ∂y  ⌡ Φ
∆V ∆V ∆V

In this case, the grid is uniform so all the cell faces have the same area A. The integrated form of the above equation is:
Γ Α ∂Φ – Γ A ∂Φ  + Γ Α ∂Φ – Γ A ∂Φ  + S∆V –
 e e  ∂x  w w
 ∂x    n n  ∂y  s s
 ∂x  
=0
 e w  n s

where the subscripts e, w, n and s denote East, West, North and South respectively.

The fluxes through the cell faces are approximated by a first-order finite difference formulation as follows:

ΓwAw   = ΓwAw
∂Φ (ΦP – ΦW)
 ∂x w δxWP

ΓeAe   = ΓeAe
∂Φ (ΦE – ΦP)
 ∂x e δxPE

ΓsAs   = ΓsAs
∂Φ (ΦP – ΦS)
 ∂x s δySP

ΓnAn   = ΓnAn
∂Φ (ΦP – ΦS)
 n
∂x δyPN

where P denotes the center of the control volume at which the property Φ is to be evaluated and

δxWP denotes the distance between the west face and the cell center.
(ΦE – ΦP) (ΦP – ΦW) (ΦP – ΦS) (ΦP – ΦS) –
ΓeAe – ΓwAw + ΓnAn – ΓsAs + S∆V = 0
δxPE δxWP δyPN δySP
COMPUTATIONAL FLUID DYNAMICS (BE MECH.) (2.12) CONDUCTION AND ADVECTION

The domain is discretized into square cells of side length dx and dy with dy = dx. Since, for the case at hand, the source
term is absent, i.e. , no heat generation by the plate and that the diffusion coefficient is
(TE – TP) (TP – TW) (TN – TP) (TP – TS)
kA – kA + kA – kA =0
dx dx dy dy
The above equation represents the discretization of the governing transport equation at the cell centers, within the domain
boundaries. The discretization of the temperature gradients in requires a special treatment since for the cell adjacent to an
edge boundary, its center is located at dx = 2 or dy = 2 from the boundary in a uniform grid.
2.10 STABILITY CRITERIA
The advection-diffusion equation in flux form,
ut= – (F)x + s
where,
F = au – dux
After applying the finite volume method this becomes (in semi-discrete form),
‘ F3/2 F1/2 –
w1 = – + + s1
h1 h1
Can Dirchlet conditions be applied to the flux term at the left hand side interface, F1/2’, this would imply that we only know
the advection part of the flux, i.e. the au part where u is known at the boundary.
2.11 INTRODUCTION TO FIRST ORDER UPWIND
• This is the simplest numerical scheme. It is the method that we used earlier in the discretization example.
• We assume that the value of φ at the face is the same as the cell centered value in the cell upstream of the face.
• The main advantages are that it is easy to implement and that it results in very stable calculation, but it also very
diffusive. Gradients in the flow field tend to be smeared out, as we will show later.
• This is often the best scheme to start calculations with.

Fig. 2.6
2.11.1 Compact Difference Schemes: Introduction
th
Consider u whose n derivative is to be evaluated. In CD schemes (Pade’) it is evaluated from,
Nr Mr
(n) 1
∑ ak uj+k =
h
n ∑ bk uj+k
k=–Nl k=–Ml
st
For 1 derivatives this can be written as
[A] {u’} = [B] {u}
COMPUTATIONAL FLUID DYNAMICS (BE MECH.) (2.13) CONDUCTION AND ADVECTION

2.12 SECOND ORDER UPWIND

• We determine the value of φ from the cell values in the two cells upstream of the face.
• This is more accurate than the first order upwind scheme, but in regions with strong gradients it can result in face
values that are outside of the range of cell values. It is then necessary to apply limiters to the predicted face values.
• There are many different ways to implement this, but second-order upwind with limiters is one of the more popular
numerical schemes because of its combination of accuracy and stability.

Fig. 2.7

2.13 QUICK CONVECTION SCHEMES

• Quick stands for Quadratic Upwind Interpolation for Convective Kinetics.
• A quadratic curve is fitted through two upstream nodes and one downstream node.
• This is a very accurate scheme, but in regions with strong gradients, overshoots and undershoots can occur. This can
lead to stability problems in the calculation.

Fig. 2.8

EXERCISE
1. Derive discretised equation of 2- Dimensional steady heat conduction equation.
2. Define Dirichlet and Neumann Boundary Conditions. What are their significance?
3. Compare explicit and implicit approach of numerical techniques. Also discuss Alternating Direction Implicit
Technique (ADI) for 1-D heat conduction equation.
COMPUTATIONAL FLUID DYNAMICS (BE MECH.) (2.14) CONDUCTION AND ADVECTION

4. Using Alternating Direction Implicit Technique, analyse numerically 2-D unsteady heat conduction equation, with
variables x, y, t.
2 2
∂T  ∂ T ∂ T2 = 0
ρc + k 2 +k
∂t  ∂x  ∂y 
5. Discuss Approach for Irregular Boundary for 2-D heat conduction problems.
6. Compute the solution of 2-D unsteady heat conduction equation by implicit Method.
7. Develop an ADI method for 2-D flow.