Appendix A

A note on the
Quasilinearization Technique

Quasilinearization technique is widely used to solve different types of non-linear ordinary as

well as partial differential equations arising in the fields of engineering, science, fluid dynamics,
solid mechanics, heat and mass transfer. It can be viewed as a generalization of the Newton-
Raphson approximation technique in functional space. This method provides a sequence of
functions which in general converge rapidly to the solution of original equations, in spite of
using uninspired initial guesses. In this method, an iterative sequence of linear equations is
carefully constructed to approximate the nonlinear equations. The unique feature of quasilin-
earization technique is quadratic convergence and monotonicity, which has been found superior
than the built - in iteration of upwind technique or finite amplitude technique. The efficiency
and accuracy of the method have been illustrated through it’s applications to many boundary
value problems, in the book by Bellman and Kalaba [75]. The existence and uniqueness of
the solution by the quasilinearization technique has been very well established [76]. Since a
majority of the boundary-layer problems are of the two-point, nonlinear, boundary value type,
a second order vector system with one independent variable has been considered in the domain
[a, b] to demonstrate the quasilinearization method. However, extension to the systems with
more independent variables is simple and obvious.

Let us consider the vector equation

Y00 = Φ(x, Y, Y0 ), a≤x≤b (A.1)

169
Quasilinearization 170

subject to the boundary conditions

Y(a) = A, Y(b) = B (A.2)

where Y is a vector composed of n dependent variables in the system and x is the independent
variable. The prime (0 ) denotes the differentiation with respect to x, Φ is a vector function
which gives the derivatives of the dependent variable, and A and B are constant vectors whose
values are known.

Applying the method of quasilinearization to the above system, we obtain a linear system of
equations which, in vector notation, can be written as

n
X
(k+1) 00 (k) ∂ (k) Φ £ (k+1)
¤
Y = Φ+ Yi − (k) Yi
i=1
∂Yi
X n
∂ (k) Φ £(k+1) 0 (k) 0 ¤
+ Yi − Yi (A.3)
i=1
∂Yi

where the subscript i denotes the ith component of the vector, the supercript (k) denotes the k th
∂ (k) Φ
approximate and ∂Yi
denote the derivative expressed in terms of k th approximates. Equation
(B.3) can be rearranged and written precisely as a linear system in terms of the (k + 1)th
approximate, to be determined as

(k+1)
Y00 = (k) P (k+1)
Y0 + (k) Q (k+1) Y + (k) R (A.4)

where the coefficients P ,Q and R are matrices of functions of the k th iterative and are therefore
known.
Appendix B

Outline of
Finite-Difference Methods

The finite difference method (FDM) was first developed by A. Thom [181] in the 1920s under
the title “the method of square” to solve nonlinear hydrodynamic equations. The finite differ-
ence technique are based upon the approximations that permit replacing differential equations
by finite differences.
The method of finite differences are widely used in Computational Fluid Dynamics (CFD),
and these methods are extensively applied to the solution of both linear and nonlinear, ordi-
nary as well as partial differential equations. Indeed, in these methods the region of integration
is divided into a network of computational cells by a generally fixed orthogonal grid. Or-
dinary/partial derivatives of functions in various directions are replaced by finite-differences
such as central or backward differences. These finite difference approximations are in alge-
braic form and the solutions are related to grid points.

Application of the finite difference method to any physical problem consists of following
three steps.

1. Division of spatial domain into an orthogonal computational grid. The grid levels are
conveniently indexed in case of two dimensions by integers i and j which increase mono-
tonically along the x- and y- directions, respectively. The intersections of the grid levels
define a set of node points at which dependent variable is defined.

2. Discretization of the governing equations and boundary conditions in space and time to
derive approximately equivalent algebraic equation for each node, using Taylor series

171
Finite-difference 172

expansion. The spatial discretization of the governing partial differential equation is per-
formed over the finite difference grid, while time discretization is carried out by marching
in time direction.

3. Solution of the discretized (algebraic) equations is obtained by a suitable matrix inversion

along with an iterative technique.

For boundary layer equations which are parabolic in nature, the central differences are
used in η-direction (along the boundary layer thickness) and backward differences are used in
streamwise direction ξ as well as in time direction t∗ with constant step-size 4η, 4ξ and 4t∗
in η- ξ- and t∗ -directions, respectively. The mesh-point diagram for finite difference scheme
has been shown in Fig.C.1. Consequently, the relevant finite-differences are given by

X = ªXm,n,p + (1 − ª)Xm+1,n,p

X 0 = [ª(Xm,n,p+1 − Xm,n,p−1 ) + (1 − ª)(Xm,n+1,p−1 − Xm,n+1,p−1 )] /2∆η

Xt∗ = [ª(Xm,n,p − Xm−1,n,p ) + (1 − ª)(Xm,n+1,p − Xm−1,n+1,p )] /∆t∗

Xξ = [ª(Xm,n,p − Xm,n−1,p ) + (1 − ª)(Xm+1,n,p − Xm+1,n−1,p )] /∆ξ

X 00 = [ª(Xm,n,p+1 − 2Xm,n,p + Xm,n,p−1 ) + (1 − ª)(Xm,n+1,p+1

− 2Xm,n+1,p + Xm,n+1,p−1 )]/(∆η)2 (B.1)

The subscripts m, n and p denote particular locations corresponding to t∗ , ξ and η respectively.

However, the above finite differences are written with a parameter ª which will give the various
finite difference schemes as indicated below:


 0 Explicit method


ª= 1
Crank-Nicolson method

 2

 1 Implicit method

Usually, in an explicit scheme each equation is allowed the determination of one unknown
quantity in terms of known quantities, while an implicit scheme involves the solution of a set
of simultaneous algebraic equations at each station. The explicit scheme is conditionally stable
and requires a very small mesh size. Hence they are found to be unsuitable for engineering
problems. On the other hand, implicit schemes are unconditionally stable. It is appropriate
Finite-difference 173

to remark here that, a numerical method is considered to be convergent if the solution of the
discretized equation tends to the exact solution of the differential equation as the grid spacing
tends to zero. For initial value problems, the Lax equivalence theorem (Fletcher, 1988) states
that “given a properly posed linear initial value problems and a finite difference approxima-
tion to it that satisfies the consistency condition, stability is necessary and sufficient condition
for convergence”. The consistency aspect of finite difference scheme can be attributed to the
discretisation of governing partial differential equations which becomes exact as the grid spac-
ing tends to zero. The difference between the discretized equation and the exact is actually
introduces the truncation error. It is usually estimated by replacing all the nodal values in the
discrete approximation by a Taylor series expansion about a single point. As a result one recov-
ers the original partial differential equations plus a remainder, which represent the truncation
error.

For a method to be consistent, the truncation error must become zero when the mesh spacing
become very very small. Also, a numerical solution is stable if it ensures a bounded solution,
without magnifying the errors that appear in the course of numerical solution process. The
Lax equivalence theorem is of great importance, since it is relatively easy to show stability of
an algorithm and its consistency with the original partial differential equation, whereas it is
usually very difficult to show convergence of its solution to that of the partial differential equa-
tion. In fact, for non- linear problems, which are strongly influenced by boundary conditions,
the stability and convergence of a method are difficult to demonstrate. Therefore, convergence
is usually checked using numerical experiments, i.e., repeating the calculation as a series of
successive refined grids. If the discretization process are consistent, we usually find that the
solution does converge to a grid independent solutions. In this thesis repeated trials of numer-
ical computations have performed and, only grid independent results have been presented, to
conserve the space.

The quasilinearization technique explained in Appendix B is used along with finite-difference

scheme. Indeed, quasilinearization technique makes the system of nonlinear partial differen-
tial equations to become linear, locally which helps in the faster convergence of the numerical
solution. The difference equations written by using above finite-differences along with quasi-
Finite-difference 174

linearization technique are expressed in the general form of matrix equation as:

An Wn−1 + Bn Wn + Cn Wn+1 = Dn (2 ≤ n ≤ N̄ − 1) (B.2)

where the vectors and the coefficient matrices are of the form
     
F  D1  A11 A12 
Wn =  
  , Dn =  
  , An = 



G D2 A21 A22
n n

and Bn and Cn have similar expressions as An . Also, F and G are dependent variables. The
matrix Eqn.(C.2) can be solved by using Varga’s algorithm [113], to get the required solution.
The algorithm is as follows:

Wn = −En Wn+1 + Jn+1 , 1 ≤ n ≤ N̄ − 1

En = (Bn − An En−1 )−1 Cn

where

Jn = (Bn − An En−1 )−1 (Dn − An Jn−1 ), 2 ≤ n ≤ N̄ − 1

   
F (a) F (b)
E1 = EN̄ = 0, J1 =  
,
 J N̄ = 



G(a) G(b)

The column matrices J1 and JN̄ are known, since they represent the boundary conditions
a and b.
Finite Difference method .. 175

Fig. B.1 Mesh point diagram for finite difference scheme

Appendix C

A Brief Note on
Keller box Method

H. B. Keller [85] introduced a method for obtaining similar and non-similar solutions for
boundary layer flow problems, which became popular as Keller-Box method. This method
is having desirable advantages that make it appropriate for the solution of all parabolic partial
differential equations. The main features of this method are:

1. The method is simple, efficient and easy to program;

2. The method enables one to compute extremely close to the point of boundary layer sepa-
ration with no special precautions;

3. Allows very rapid λ - direction (stream wise coordinate) variations;

4. Only slightly more arithmetic calculations are needed to solve as compared to other finite-
difference schemes;

5. The method is unconditionally stable and has second order accuracy with arbitrary (nonuni-
form) λ- and η- spacings.

The general net rectangle in the λ − η plane is shown in Fig.D.1. The net points are defined as:

λ0 = 0 λn = λn−1 + kn , n = 1, 2, ......., J (C.1)

η0 = 0 ηj = ηj−1 + hj , j = 1, 2, ......., J, ηJ = η∞ (C.2)

where kn is the ∆λ-spacing and hj is the ∆η-spacing. Here n and j are just the sequence
of numbers that indicate the coordinate location, not tensor, indices or exponents. The finite-

176
Keller box Method 177

differences in this scheme are expressed as

1h n i
( )nj−1/2 = ( )j + ( )nj−1 (C.3)
2
1h n i
( )n−1/2
j = ( )j + ( )n−1
j (C.4)
2
µ ¶n−1/2 n−1
∂u unj−1/2 − uj−1/2
= (C.5)
∂λ j−1/2 kn
µ ¶n−1/2 n−1/2 n−1/2
∂u uj − uj−1
= (C.6)
∂η j−1/2 hj

The solution of a differential equations by this method can be obtained by the following four
steps:

1. Reduce the equation or system of equations to a first order system;

2. Write the difference equations using central differences;

3. Linearize the resulting algebraic equations (if they are non-linear) by Newton’s method;

4. And, write them in matrix-vector form and use the block-tridiagonal-elimination tech-
nique to solve the linear system.

Keller box method has been widely used and tested extensively on laminar and turbulent bound-
ary layer flows. This method is found to be much faster, more efficient and more flexible to
use, than other methods.
Keller-Box method 178

Fig. C.1 Net rectangle in η – λ plane for Keller-Box method

Appendix D

Introduction to
Magnetohydrodynamics

Magneto-hydrodynamics (MHD) is the branch of continuum mechanics which deals with the
interaction of electrically conducting fluids and electro magnetic forces. The field of MHD
was initiated by Swedish physicist, Hannes Alfven for which he received in 1970 the Noble
prize. The idea of MHD is that magnetic fields can induce currents in a moving conductive
fluid, which create forces on the fluids, and also change the magnetic field itself. It unifies in
a common frame work the electromagnetic and fluid dynamic theories, to yield a description
of the concurrent effects of the magnetic field on the flow and, the flow on the magnetic field.
The most appropriate name for the phenomena would be Magneto Fluid Mechanics, but the
original name MHD is still generally used for the flow analysis of incompressible electrically
conducting fluids including conducting liquids and gases.

For an incompressible hydromagnetic fluid flow, the basic equations are

DU
ρ = ∇.τi + J × B (D.1)
Dt
∇.U = 0 (D.2)

∇.B = 0 (Magnetic field continuity) (D.3)

ρ
∇.E = (Flux continuity) (D.4)
ε0
1 ∂E
∇ × B = µm J + 2 (Ampere’s law) (D.5)
c0 ∂t

∂B
∇×E = − (Faraday’s law) (D.6)
∂t
179
Magnetohydrodynamics 180

where ρ is the density of the fluid, U the velocity field, D/Dt is the time derivative (in which
∂U/∂t is the local acceleration and (U.∇)U is the advection term), (∇.) is the divergence and
τi = −δI + µγ is the Newtonian stress in which −δI is the indeterminate part of the spherical
stress and
γ = ∇U + (∇U)T (D.7)

is the rate of deformation tensor [182]. Further, ε0 represents the permittivity, c0 is the speed of
light, J is the current density, µm = 1/ε0 c20 the magnetic field permeability, E the total electric
field current, B is the magnetic field.

To precede the discussion the following assumptions are made [183,55,184]: The density ρ,
magnetic permeability µm are considered constant throughout the flow field region. The electric
field conductivity σ is also considered as constant and assumed to finite. Total magnetic field B
is perpendicular to the velocity field U and the induced magnetic field b is negligible compared
with the applied magnetic field B0 so that the magnetic Reynolds number is small [185].

We assume a situation where no energy is added or extracted from the fluid by the electric
field, which implies that there is no electric field present in the fluid flow region.

Under these conditions, we readily obtain the simplified form of the MHD Eqns.(A.3)-(A.6),
that is

J × B = σ(U × B) × B = σB20 U (D.8)

where, (J×B) is actually, the Lorentz force on the conducting fluid produced by the interaction
of the current and the magnetic field. The velocity field is given by

U = [u(y, x), U(x), 0] (D.9)

which follows from the continuity equation. The magnetic field usually applied in the y-
direction normal to the boundary layer flow. The flow is in x-direction whose velocity com-
σB02
ponent is u. Thus J × B, representing Lorentz force becomes − ρ
u. The negative sign in
the term is because of retardation. The Eqn.(A.1) with Eqn.(A.8) in the absence of pressure
Magnetohydrodynamics 181

gradient, becomes
∂u ∂u ∂u ∂ 2 u σB02
+u +v =ν 2 − u (D.10)
∂t ∂x ∂y ∂y ρ
σB02
where ν = µ/ρ is the kinematic viscosity of the fluid and the term ρ
u describes the x-
σB02
component of the magnetic field where ρ
is the magnetic parameter which is the ratio of
the electromagnetic force to the inertial force. In fact, the Eqn.(A.10) represents a typical un-
steady, two-dimensional momentum boundary layer equation with an applied transverse mag-
netic field.