# Spectral Relaxation Method and Spectral Quasilinearization Method for Solving Unsteady Boundary Layer Flow Problems

Authors: S. S. Motsa; P. G. Dlamini; M. Khumalo
Journal: Advances in Mathematical Physics
 (2014)

---

## Abstract

Nonlinear partial differential equations (PDEs) modelling unsteady boundary-layer flows are solved by the spectral relaxation method (SRM) and the spectral quasilinearization method (SQLM). The SRM and SQLM are Chebyshev pseudospectral based methods that have been successfully used to solve nonlinear boundary layer flow problems described by systems of ordinary differential equations. In this paper application of these methods is extended, for the first time, to systems of nonlinear PDEs that model unsteady boundary layer flow. The new extension is tested on two problems: boundary layer flow caused by an impulsively stretching plate and a coupled four-equation system that models the problem of unsteady MHD flow and mass transfer in a porous space. Numerous simulation experiments
are conducted to determine the accuracy and compare the computational performance of the proposed methods against the popular Keller-box finite difference
scheme which is widely accepted as being one of the ideal tools for solving nonlinear PDEs that model boundary layer flow problems. The results indicate that the
methods are more efficient in terms of computational accuracy and speed compared with the Keller-box.


---

## Body

## 1. Introduction

Partial differential equations (PDEs) arise in a number of physical problems, such as fluid flow, heat transfer, and biological processes. Finding solutions of the PDEs plays a crucial role in understanding the behaviour of these problems. Mostly, the PDEs modelling real-life problems are nonlinear and complex to solve exactly and hence various analytical and numerical methods have been employed to approximate the solutions of these problems. In recent times, many researchers in fluid mechanics have focused their attention on problems involving boundary layer flows of an incompressible fluid over a stretching surface because of their substantial applications in engineering. A large and growing body of literature has investigated problems involving steady flows. However, in some cases the flow field could be unsteady due to a sudden stretching of the flat sheet. Unsteady flows are mostly defined by systems of nonlinear PDEs and are considerably more difficult to solve than steady flows problems which are often simplified into system nonlinear ODEs using the so-called similarity transformations.The problem of unsteady boundary layer flow due to an impulsively stretching surface in a viscous fluid has been considered by a number of researchers. These studies include the work of Seshadri et al. [1] who used the Keller-box method of Cebeci and Bradshaw [2] and a perturbation series approach for the solution of unsteady mixed convection flow along a heated vertical plate. The Keller-box method was also used by Ali et al. [3] to solve a related problem of unsteady boundary layer flow due to an impulsively stretching surface. Nazar et al. [4, 5] solved the unsteady boundary-layer flow problem due to an impulsively stretching surface in a rotating fluid by means of the Keller-box numerical method, and they obtained a first order perturbation approximation of the solution. Liao [6] noted that a limiting factor of the perturbation approach is that it gives solutions that are only valid for small time. As an alternative approach, Liao [6] suggested the use of the homotopy analysis method (HAM) that was meant to address some of the limitations of the perturbation methods by offering solutions that are uniformly valid for all time. In recent years, there has been an increasing amount of literature that has adopted Liao’s analytic approach in solving unsteady boundary layer flows. However, there are limits to how far analytic approaches can be utilised in nonlinear systems of PDEs involving many equations. Nonlinear systems involving many coupled equations are very difficult to solve analytically.In this work, we apply, for the first time, the spectral relaxation method (SRM) and the spectral quasilinearization method (SQLM) to solve nonlinear PDEs describing unsteady boundary layer flow due to an impulsively stretching surface. The SRM was introduced in [7] for the solution of the nonlinear ODE system model of von Karman flow of a Reiner-Rivlin fluid. The method has also been used in the solution of chaotic and hyperchaotic systems [8, 9]. The SRM is based on simple decoupling and rearrangement of the governing nonlinear equations in a Gauss-Seidel manner. The resulting sequence of equations is integrated using the Chebyshev spectral collocation method. On the other hand, in the SQLM, the governing nonlinear equations are linearised using the Newton-Raphson based quasilinearization method (QLM), developed by Bellman and Kalaba [10], and are then integrated using Chebyshev spectral collocation method. A sizeable body of literature now exists on the use of various finite difference based QLM schemes in boundary layer flows described by both nonlinear ODE and PDE-based systems [11–15]. Spectral method based quasilinearisation schemes have also been successfully applied to a range of fluid mechanics based ODE model problems (see, e.g., [16–18]). For problems with smooth solutions, spectral methods are well known [19–21] to be considerably more accurate than other traditional numerical methods such as finite difference and finite elements. In this investigation we revisit the one-dimensional unsteady boundary layer flow due to an impulsively stretching surface that was previously discussed by [6] using the homotopy analysis method and recently in [3] using the Keller-box method. The problem of unsteady three-dimensional MHD flow and mass transfer in a porous space [22] is also investigated. The main purpose of the study is to investigate the applicability and effectiveness of the new SRM approach to systems of nonlinear PDE-based unsteady boundary layer flows of varying levels of complexity. Numerical simulations are conducted on the sample problems using the SRM, SQLM, and Keller-box method. The three methods are compared in terms of accuracy, computational speed, and easy implementation.The rest of the paper is organized as follows. In Section2, we discuss the development of the SRM and SQLM for the solution of an unsteady boundary-layer flow caused by an impulsively stretching plate. Section 3 presents the SRM and SQLM implementation of an unsteady three-dimensional MHD flow and mass transfer in a porous space. Section 4 contains the results and discussion, and the conclusions are given in Section 5.


## 2. Unsteady Boundary-Layer Flows Caused by an Impulsively Stretching Plate

The governing partial differential equations can be obtained by using the standard stream function formulation in conjunction with the transformations suggested by Williams and Rhyne [23]. The dimensionless governing equation is obtained (see [1, 4, 6] for details) as
(1)
∂
3
f
∂
η
3
+
1
2
(
1
-
ξ
)
η
∂
2
f
∂
η
2
+
ξ
[
f
∂
2
f
∂
η
2
-
(
∂
f
∂
η
)
2
]
=
ξ
(
1
-
ξ
)
∂
2
f
∂
ξ
∂
η
,
subject to the boundary conditions
(2)
f
(
0
,
ξ
)
=
0
,
∂
f
∂
η
|
η
=
0
=
1
,
∂
f
∂
η
|
η
→
+
∞
=
0
,
where the primes denote differentiation with respect to the similarity variable η. f is a nondimensional function that gives the velocity and ξ
∈
[
0,1
] is the dimensionless time-scale defined as
(3)
ξ
=
1
-
e
-
τ
,
τ
=
b
t
,
where b is a positive constant and t is the time variable. In the analysis of boundary layer flow problems, a quantity that is of physical interest in the skin friction in this model is given [1, 4, 6], in dimensionless form, as
(4)
C
f
Re
x
1
/
2
=
ξ
-
1
/
2
f
′′
(
ξ
,
0
)
,
where Re
x is the local Reynolds number.The initial unsteady solution atξ
=
0 (τ
=
0) for the governing equation (1) is obtained as a solution of the equation
(5)
f
′′′
+
1
2
η
f
′′
=
0
,
(6)
f
(
0,0
)
=
0
,
f
′
(
0,0
)
=
1
,
f
′
(
∞
,
0
)
=
0
,
where the primes denote differentiation with respect to η. Solving (5) gives
(7)
f
(
η
,
0
)
=
η
erfc
⁡
(
η
2
)
+
2
π
[
1
-
exp
⁡
(
-
η
2
4
)
]
,
where erfc
⁡
(
η
) is the standard complementary error function defined by
(8)
erfc
⁡
(
η
)
=
2
π
∫
η
∞
exp
⁡
⁡
(
-
z
2
)
d
z
.The steady state solution whenξ
=
1, corresponding to τ
→
+
∞, is obtained from
(9)
f
′′′
+
f
f
′′
-
(
f
′
)
2
=
0
,
f
(
0,1
)
=
0
,
f
′
(
0,1
)
=
1
,
f
′
(
∞
,
1
)
=
0
.
The solution to the above equation is
(10)
f
(
η
,
1
)
=
1
-
exp
⁡
⁡
(
η
)
.


### 2.1. Spectral Relaxation Method(SRM)

In this section we discuss the development of the spectral relaxation method and its application to solve the partial differential equation (1). It is convenient to reduce the order of (1) from three to two. To this end, we set f
′
=
u, so that (1) becomes
(11)
u
′′
+
1
2
η
(
1
-
ξ
)
u
′
+
ξ
[
f
u
′
-
u
2
]
=
ξ
(
1
-
ξ
)
∂
u
∂
ξ
,
f
′
=
u
.The spectral relaxation method [7] algorithm uses the idea of the Gauss-Seidel method to decouple the governing systems of (11). From the decoupled equations an iteration scheme is developed by evaluating linear terms in the current iteration level (denoted by r
+
1) and nonlinear terms in the previous iteration level (denoted by r). Applying the SRM on (11) gives the following linear partial differential equations:
(12)
u
r
+
1
′′
+
a
1
,
r
u
r
+
1
′
+
a
2
,
r
=
ξ
(
1
-
ξ
)
∂
u
r
+
1
∂
ξ
,
(13)
u
r
+
1
(
0
,
ξ
)
=
1
,
u
r
+
1
(
∞
,
ξ
)
=
0
,
(14)
f
r
+
1
′
=
u
r
+
1
,
f
r
+
1
(
0
,
ξ
)
=
0
,
where
(15)
a
1
,
r
=
1
2
η
(
1
-
ξ
)
+
ξ
f
r
,
a
2
,
r
=
-
ξ
u
r
2
.The initial approximation for solving (12)–(14) is obtained as the solutions at ξ
=
0. Thus f
0
(
η
,
ξ
) and u
0
(
η
,
ξ
) are given by
(16)
f
0
(
η
,
ξ
)
=
η
erfc
⁡
(
η
2
)
+
2
π
[
1
-
exp
⁡
(
-
η
2
4
)
]
,
u
0
(
η
,
ξ
)
=
erfc
⁡
(
η
2
)
.
Starting from given initial approximations (16), the iteration schemes (12) can be solved iteratively for u
r
+
1
(
η
,
ξ
) when r
=
0,1
,
2
,
…. The solution for u
r
+
1 is used in (14) which is, in turn, solved for f
r
+
1. To solve (12) we discretize the equation using the Chebyshev spectral method in the η-direction and use an implicit finite difference method in the ξ-direction. For details of the spectral method, we refer interested readers to [19, 21]. Before applying the spectral method, it is convenient to transform the domain on which the governing equation is defined to the interval [
-
1,1
] where the spectral method can be implemented. For convenience, the semi-infinite domain in the space direction is approximated by the truncated domain [
0
,
η
∞
], where η
∞ is a finite number selected to be large enough to represent the behaviour of the flow properties when η is very large. We use the transformation η
=
η
∞
(
Y
+
1
)
/
2 to map the interval [
0
,
η
∞
] to [
-
1,1
]. The basic idea behind the spectral collocation method is the introduction of a differentiation matrix D which is used to approximate the derivatives of the unknown variables f
(
η
) at the collocation points (grid points) as the matrix vector product
(17)
d
F
d
η
=
∑
k
=
0
N
x
D
j
k
f
(
η
k
)
=
D
F
,
j
=
0,1
,
…
,
N
x
,
where N
x
+
1 is the number of collocation points, D
=
2
D
/
η
∞, and
(18)
F
=
[
f
(
Y
0
)
,
f
(
Y
1
)
,
…
,
f
(
Y
N
x
)
]
T
is the vector function at the collocation points. Higher order derivatives are obtained as powers of D; that is,
(19)
F
(
p
)
=
D
p
F
,
where p is the order of the derivative. We choose the Gauss-Lobatto collocation points to define the nodes in [
-
1,1
] as
(20)
Y
j
=
cos
⁡
⁡
(
π
j
N
x
)
,
j
=
0,1
,
…
,
N
x
.
The matrix D is of size (
N
x
+
1
)
×
(
N
x
+
1
). The grid points on (
η
,
ξ
) are defined as
(21)
Y
j
=
cos
⁡
π
j
N
x
,
ξ
n
=
n
Δ
ξ
,
j
=
0,1
,
…
,
N
x
,
0000000000000000000
n
=
0,1
,
…
,
N
t
,
where N
x
+
1, N
t
+
1 are the total number of grid points in the η and ξ-directions, respectively, and Δ
ξ is the spacing in the ξ-direction. The finite difference scheme is applied with centering about a midpoint halfway between ξ
n
+
1 and ξ
n. This midpoint is defined as ξ
n
+
(
1
/
2
)
=
(
ξ
n
+
1
+
ξ
n
)
/
2. The derivatives with respect to η are defined in terms of the Chebyshev differentiation matrices. Applying the centering about ξ
n
+
(
1
/
2
) to any function, say u
(
η
,
ξ
) and its associated derivative, we obtain
(22)
u
(
η
j
,
ξ
n
+
(
1
/
2
)
)
=
u
j
n
+
(
1
/
2
)
=
u
j
n
+
1
+
u
j
n
2
,
(
∂
u
∂
ξ
)
n
+
(
1
/
2
)
=
u
j
n
+
1
-
u
j
n
Δ
ξ
.Before applying the finite differences, we apply the spectral method on (12) and (14) to obtain
(23)
[
D
2
+
a
1
,
r
D
]
U
r
+
1
+
a
2
,
r
=
ξ
(
1
-
ξ
)
d
U
r
+
1
d
ξ
,
(24)
u
r
+
1
(
η
0
,
ξ
)
=
0
,
u
r
+
1
(
η
N
x
,
ξ
)
=
1
,
(25)
u
r
+
1
(
η
j
,
0
)
=
erfc
⁡
(
η
j
2
)
,
j
=
0,1
,
2
,
…
,
N
x
,
(26)
D
F
r
+
1
=
U
r
+
1
,
f
r
+
1
(
η
N
x
,
ξ
)
=
0
,
(27)
U
r
+
1
=
[
u
r
+
1
(
η
0
,
ξ
)
u
r
+
1
(
η
1
,
ξ
)
⋮
u
r
+
1
(
η
N
x
-
1
,
ξ
)
u
r
+
1
(
η
N
x
,
ξ
)
]
,
F
r
+
1
=
[
f
r
+
1
(
η
0
,
ξ
)
f
r
+
1
(
η
1
,
ξ
)
⋮
f
r
+
1
(
η
N
x
-
1
,
ξ
)
f
r
+
1
(
η
N
x
,
ξ
)
]
,
(28)
a
2
,
r
=
[
a
2
,
r
(
η
0
,
ξ
)
a
2
,
r
(
η
1
,
ξ
)
⋮
a
2
,
r
(
η
N
x
-
1
,
ξ
)
a
2
,
r
(
η
N
x
,
ξ
)
]
,
(29)
a
1
,
r
=
[
a
1
,
r
(
η
0
,
ξ
)
a
1
,
r
(
η
2
,
ξ
)
⋱
⋱
a
1
,
r
(
η
N
x
,
ξ
)
]
.Next, we apply the finite difference scheme on (23) in the ξ-direction with centering about the midpoint ξ
n
+
(
1
/
2
) to obtain
(30)
A
U
r
+
1
n
+
1
=
B
U
r
+
1
n
+
K
,
subject to the following boundary and initial conditions:
(31)
u
r
+
1
(
η
0
,
ξ
n
)
=
0
,
u
r
+
1
(
η
N
x
,
ξ
n
)
=
1
,
00000000000000000
n
=
0,1
,
2
,
…
,
N
t
,
u
r
+
1
(
η
j
,
0
)
=
erfc
⁡
(
η
j
2
)
,
j
=
0,1
,
2
,
…
,
N
x
,
where
(32)
A
=
-
1
2
(
D
2
+
a
1
,
r
n
+
(
1
/
2
)
D
)
+
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
I
,
B
=
1
2
(
D
2
+
a
1
,
r
n
+
(
1
/
2
)
D
)
+
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
I
,
K
=
a
2
,
r
n
+
(
1
/
2
)
,
where I is an (
N
x
+
1
)
×
(
N
x
+
1
).Starting from the initial conditionU
r
+
1
0, given by (16), (30) can be solved iteratively to give approximate solutions for u
r
+
1
(
η
,
ξ
),  r
=
0,1
,
2
,
…, until a solution that converges to within a given accuracy level is obtained. The solution u
r
+
1 is used in (26) which is, in turn, solved for f
r
+
1.


### 2.2. Spectral Quasilinearization Method(SQLM)

In this section we present the spectral quasilinearization method (SQLM) for solving the partial differential equation (1). The quasilinearization technique is essentially a generalized Newton-Raphson Method that was originally used by Bellman and Kalaba [10] for solving functional equations. We first set f
′
=
u, so that (1) becomes
(33)
u
′′
+
1
2
η
(
1
-
ξ
)
u
′
+
ξ
[
f
u
′
-
u
2
]
=
ξ
(
1
-
ξ
)
∂
u
∂
ξ
.Applying the QLM on (33) the nonlinear partial differential equation reduces to the following iterative sequence of linear partial differential equations:
(34)
u
r
+
1
′′
+
a
1
,
r
u
r
+
1
′
+
a
2
,
r
u
r
+
1
+
a
3
,
r
f
r
+
1
+
a
4
,
r
=
ξ
(
1
-
ξ
)
∂
u
r
+
1
∂
ξ
,
(35)
u
r
+
1
(
0
,
ξ
)
=
1
,
u
r
+
1
(
∞
,
ξ
)
=
0
,
(36)
f
r
+
1
′
=
u
r
+
1
,
f
r
+
1
(
0
,
ξ
)
=
0
,
where
(37)
a
1
,
r
=
1
2
η
(
1
-
ξ
)
+
ξ
f
r
,
a
2
,
r
=
-
2
ξ
u
r
,
a
3
,
r
=
ξ
u
r
′
a
4
,
r
=
-
ξ
f
r
u
r
′
+
ξ
u
r
2
.The indicesr and r
+
1 denote the previous and current iteration levels, respectively.Starting from given initial approximations, denoted byf
0
(
η
,
ξ
) and u
0
(
η
,
ξ
), (34)–(36) can be solved iteratively for u
r
+
1
(
η
,
ξ
) and f
r
+
1
(
η
,
ξ
) (r
=
0,1
,
2
,
…). We discretize (34) and (36) using the Chebyshev spectral method in the η-direction and we use the implicit finite difference method in the ξ-direction to discretize (34) as described in the previous section. Applying the spectral method and finite differences on (34) and (36) as described previously, we obtain
(38)
[
A
1,1
A
1,2
A
2,1
A
2,2
]
[
U
r
+
1
n
+
1
F
r
+
1
n
+
1
]
=
[
B
1,1
B
1,2
O
O
]
[
U
r
+
1
n
F
r
+
1
n
]
+
[
K
1
K
2
]
subject to the boundary and initial conditions (31), where
(39)
A
1,1
=
1
2
(
D
2
+
a
1
,
r
n
+
(
1
/
2
)
D
+
a
2
,
r
n
+
(
1
/
2
)
)
-
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
I
,
A
1,2
=
1
2
a
3
,
r
n
+
(
1
/
2
)
,
A
2,1
=
-
I
,
A
2,2
=
-
D
,
B
1,1
=
-
1
2
(
D
2
+
a
1
,
r
n
+
(
1
/
2
)
D
+
a
2
,
r
n
+
(
1
/
2
)
)
-
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
I
,
B
1,2
=
-
1
2
a
3
,
r
n
+
(
1
/
2
)
,
K
1
=
a
4
,
r
n
+
(
1
/
2
)
,
K
2
=
O
¯
,
where I is an (
N
x
+
1
)
×
(
N
x
+
1
) identity matrix, O is an (
N
x
+
1
)
×
(
N
x
+
1
) zero matrix, and O
¯ is an (
N
x
+
1
)
×
1 zero vector. Starting from the initial condition U
r
+
1
0, (38) can be solved iteratively to give approximate solutions for u
r
+
1
(
η
,
ξ
) and f
r
+
1
(
η
,
ξ
), r
=
0,1
,
2
,
…, until a solution that converges to within a given accuracy level is obtained.


## 2.1. Spectral Relaxation Method(SRM)

In this section we discuss the development of the spectral relaxation method and its application to solve the partial differential equation (1). It is convenient to reduce the order of (1) from three to two. To this end, we set f
′
=
u, so that (1) becomes
(11)
u
′′
+
1
2
η
(
1
-
ξ
)
u
′
+
ξ
[
f
u
′
-
u
2
]
=
ξ
(
1
-
ξ
)
∂
u
∂
ξ
,
f
′
=
u
.The spectral relaxation method [7] algorithm uses the idea of the Gauss-Seidel method to decouple the governing systems of (11). From the decoupled equations an iteration scheme is developed by evaluating linear terms in the current iteration level (denoted by r
+
1) and nonlinear terms in the previous iteration level (denoted by r). Applying the SRM on (11) gives the following linear partial differential equations:
(12)
u
r
+
1
′′
+
a
1
,
r
u
r
+
1
′
+
a
2
,
r
=
ξ
(
1
-
ξ
)
∂
u
r
+
1
∂
ξ
,
(13)
u
r
+
1
(
0
,
ξ
)
=
1
,
u
r
+
1
(
∞
,
ξ
)
=
0
,
(14)
f
r
+
1
′
=
u
r
+
1
,
f
r
+
1
(
0
,
ξ
)
=
0
,
where
(15)
a
1
,
r
=
1
2
η
(
1
-
ξ
)
+
ξ
f
r
,
a
2
,
r
=
-
ξ
u
r
2
.The initial approximation for solving (12)–(14) is obtained as the solutions at ξ
=
0. Thus f
0
(
η
,
ξ
) and u
0
(
η
,
ξ
) are given by
(16)
f
0
(
η
,
ξ
)
=
η
erfc
⁡
(
η
2
)
+
2
π
[
1
-
exp
⁡
(
-
η
2
4
)
]
,
u
0
(
η
,
ξ
)
=
erfc
⁡
(
η
2
)
.
Starting from given initial approximations (16), the iteration schemes (12) can be solved iteratively for u
r
+
1
(
η
,
ξ
) when r
=
0,1
,
2
,
…. The solution for u
r
+
1 is used in (14) which is, in turn, solved for f
r
+
1. To solve (12) we discretize the equation using the Chebyshev spectral method in the η-direction and use an implicit finite difference method in the ξ-direction. For details of the spectral method, we refer interested readers to [19, 21]. Before applying the spectral method, it is convenient to transform the domain on which the governing equation is defined to the interval [
-
1,1
] where the spectral method can be implemented. For convenience, the semi-infinite domain in the space direction is approximated by the truncated domain [
0
,
η
∞
], where η
∞ is a finite number selected to be large enough to represent the behaviour of the flow properties when η is very large. We use the transformation η
=
η
∞
(
Y
+
1
)
/
2 to map the interval [
0
,
η
∞
] to [
-
1,1
]. The basic idea behind the spectral collocation method is the introduction of a differentiation matrix D which is used to approximate the derivatives of the unknown variables f
(
η
) at the collocation points (grid points) as the matrix vector product
(17)
d
F
d
η
=
∑
k
=
0
N
x
D
j
k
f
(
η
k
)
=
D
F
,
j
=
0,1
,
…
,
N
x
,
where N
x
+
1 is the number of collocation points, D
=
2
D
/
η
∞, and
(18)
F
=
[
f
(
Y
0
)
,
f
(
Y
1
)
,
…
,
f
(
Y
N
x
)
]
T
is the vector function at the collocation points. Higher order derivatives are obtained as powers of D; that is,
(19)
F
(
p
)
=
D
p
F
,
where p is the order of the derivative. We choose the Gauss-Lobatto collocation points to define the nodes in [
-
1,1
] as
(20)
Y
j
=
cos
⁡
⁡
(
π
j
N
x
)
,
j
=
0,1
,
…
,
N
x
.
The matrix D is of size (
N
x
+
1
)
×
(
N
x
+
1
). The grid points on (
η
,
ξ
) are defined as
(21)
Y
j
=
cos
⁡
π
j
N
x
,
ξ
n
=
n
Δ
ξ
,
j
=
0,1
,
…
,
N
x
,
0000000000000000000
n
=
0,1
,
…
,
N
t
,
where N
x
+
1, N
t
+
1 are the total number of grid points in the η and ξ-directions, respectively, and Δ
ξ is the spacing in the ξ-direction. The finite difference scheme is applied with centering about a midpoint halfway between ξ
n
+
1 and ξ
n. This midpoint is defined as ξ
n
+
(
1
/
2
)
=
(
ξ
n
+
1
+
ξ
n
)
/
2. The derivatives with respect to η are defined in terms of the Chebyshev differentiation matrices. Applying the centering about ξ
n
+
(
1
/
2
) to any function, say u
(
η
,
ξ
) and its associated derivative, we obtain
(22)
u
(
η
j
,
ξ
n
+
(
1
/
2
)
)
=
u
j
n
+
(
1
/
2
)
=
u
j
n
+
1
+
u
j
n
2
,
(
∂
u
∂
ξ
)
n
+
(
1
/
2
)
=
u
j
n
+
1
-
u
j
n
Δ
ξ
.Before applying the finite differences, we apply the spectral method on (12) and (14) to obtain
(23)
[
D
2
+
a
1
,
r
D
]
U
r
+
1
+
a
2
,
r
=
ξ
(
1
-
ξ
)
d
U
r
+
1
d
ξ
,
(24)
u
r
+
1
(
η
0
,
ξ
)
=
0
,
u
r
+
1
(
η
N
x
,
ξ
)
=
1
,
(25)
u
r
+
1
(
η
j
,
0
)
=
erfc
⁡
(
η
j
2
)
,
j
=
0,1
,
2
,
…
,
N
x
,
(26)
D
F
r
+
1
=
U
r
+
1
,
f
r
+
1
(
η
N
x
,
ξ
)
=
0
,
(27)
U
r
+
1
=
[
u
r
+
1
(
η
0
,
ξ
)
u
r
+
1
(
η
1
,
ξ
)
⋮
u
r
+
1
(
η
N
x
-
1
,
ξ
)
u
r
+
1
(
η
N
x
,
ξ
)
]
,
F
r
+
1
=
[
f
r
+
1
(
η
0
,
ξ
)
f
r
+
1
(
η
1
,
ξ
)
⋮
f
r
+
1
(
η
N
x
-
1
,
ξ
)
f
r
+
1
(
η
N
x
,
ξ
)
]
,
(28)
a
2
,
r
=
[
a
2
,
r
(
η
0
,
ξ
)
a
2
,
r
(
η
1
,
ξ
)
⋮
a
2
,
r
(
η
N
x
-
1
,
ξ
)
a
2
,
r
(
η
N
x
,
ξ
)
]
,
(29)
a
1
,
r
=
[
a
1
,
r
(
η
0
,
ξ
)
a
1
,
r
(
η
2
,
ξ
)
⋱
⋱
a
1
,
r
(
η
N
x
,
ξ
)
]
.Next, we apply the finite difference scheme on (23) in the ξ-direction with centering about the midpoint ξ
n
+
(
1
/
2
) to obtain
(30)
A
U
r
+
1
n
+
1
=
B
U
r
+
1
n
+
K
,
subject to the following boundary and initial conditions:
(31)
u
r
+
1
(
η
0
,
ξ
n
)
=
0
,
u
r
+
1
(
η
N
x
,
ξ
n
)
=
1
,
00000000000000000
n
=
0,1
,
2
,
…
,
N
t
,
u
r
+
1
(
η
j
,
0
)
=
erfc
⁡
(
η
j
2
)
,
j
=
0,1
,
2
,
…
,
N
x
,
where
(32)
A
=
-
1
2
(
D
2
+
a
1
,
r
n
+
(
1
/
2
)
D
)
+
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
I
,
B
=
1
2
(
D
2
+
a
1
,
r
n
+
(
1
/
2
)
D
)
+
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
I
,
K
=
a
2
,
r
n
+
(
1
/
2
)
,
where I is an (
N
x
+
1
)
×
(
N
x
+
1
).Starting from the initial conditionU
r
+
1
0, given by (16), (30) can be solved iteratively to give approximate solutions for u
r
+
1
(
η
,
ξ
),  r
=
0,1
,
2
,
…, until a solution that converges to within a given accuracy level is obtained. The solution u
r
+
1 is used in (26) which is, in turn, solved for f
r
+
1.


## 2.2. Spectral Quasilinearization Method(SQLM)

In this section we present the spectral quasilinearization method (SQLM) for solving the partial differential equation (1). The quasilinearization technique is essentially a generalized Newton-Raphson Method that was originally used by Bellman and Kalaba [10] for solving functional equations. We first set f
′
=
u, so that (1) becomes
(33)
u
′′
+
1
2
η
(
1
-
ξ
)
u
′
+
ξ
[
f
u
′
-
u
2
]
=
ξ
(
1
-
ξ
)
∂
u
∂
ξ
.Applying the QLM on (33) the nonlinear partial differential equation reduces to the following iterative sequence of linear partial differential equations:
(34)
u
r
+
1
′′
+
a
1
,
r
u
r
+
1
′
+
a
2
,
r
u
r
+
1
+
a
3
,
r
f
r
+
1
+
a
4
,
r
=
ξ
(
1
-
ξ
)
∂
u
r
+
1
∂
ξ
,
(35)
u
r
+
1
(
0
,
ξ
)
=
1
,
u
r
+
1
(
∞
,
ξ
)
=
0
,
(36)
f
r
+
1
′
=
u
r
+
1
,
f
r
+
1
(
0
,
ξ
)
=
0
,
where
(37)
a
1
,
r
=
1
2
η
(
1
-
ξ
)
+
ξ
f
r
,
a
2
,
r
=
-
2
ξ
u
r
,
a
3
,
r
=
ξ
u
r
′
a
4
,
r
=
-
ξ
f
r
u
r
′
+
ξ
u
r
2
.The indicesr and r
+
1 denote the previous and current iteration levels, respectively.Starting from given initial approximations, denoted byf
0
(
η
,
ξ
) and u
0
(
η
,
ξ
), (34)–(36) can be solved iteratively for u
r
+
1
(
η
,
ξ
) and f
r
+
1
(
η
,
ξ
) (r
=
0,1
,
2
,
…). We discretize (34) and (36) using the Chebyshev spectral method in the η-direction and we use the implicit finite difference method in the ξ-direction to discretize (34) as described in the previous section. Applying the spectral method and finite differences on (34) and (36) as described previously, we obtain
(38)
[
A
1,1
A
1,2
A
2,1
A
2,2
]
[
U
r
+
1
n
+
1
F
r
+
1
n
+
1
]
=
[
B
1,1
B
1,2
O
O
]
[
U
r
+
1
n
F
r
+
1
n
]
+
[
K
1
K
2
]
subject to the boundary and initial conditions (31), where
(39)
A
1,1
=
1
2
(
D
2
+
a
1
,
r
n
+
(
1
/
2
)
D
+
a
2
,
r
n
+
(
1
/
2
)
)
-
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
I
,
A
1,2
=
1
2
a
3
,
r
n
+
(
1
/
2
)
,
A
2,1
=
-
I
,
A
2,2
=
-
D
,
B
1,1
=
-
1
2
(
D
2
+
a
1
,
r
n
+
(
1
/
2
)
D
+
a
2
,
r
n
+
(
1
/
2
)
)
-
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
I
,
B
1,2
=
-
1
2
a
3
,
r
n
+
(
1
/
2
)
,
K
1
=
a
4
,
r
n
+
(
1
/
2
)
,
K
2
=
O
¯
,
where I is an (
N
x
+
1
)
×
(
N
x
+
1
) identity matrix, O is an (
N
x
+
1
)
×
(
N
x
+
1
) zero matrix, and O
¯ is an (
N
x
+
1
)
×
1 zero vector. Starting from the initial condition U
r
+
1
0, (38) can be solved iteratively to give approximate solutions for u
r
+
1
(
η
,
ξ
) and f
r
+
1
(
η
,
ξ
), r
=
0,1
,
2
,
…, until a solution that converges to within a given accuracy level is obtained.


## 3. Unsteady Three-Dimensional MHD Flow and Mass Transfer in a Porous Space

We consider the unsteady and three-dimensional flow of a viscous fluid over a stretching surface investigated by Hayat et al. [22]. The fluid is electrically conducting in the presence of a constant applied magnetic field B
0. The induced magnetic field is neglected under the assumption of a small magnetic Reynolds number. The flow is governed by the following four dimensionless partial differential equations:
(40)
f
′′′
+
(
1
-
ξ
)
(
η
2
f
′′
-
ξ
∂
f
′
∂
ξ
)
+
ξ
[
(
f
+
g
)
f
′′
-
(
f
′
)
2
-
M
2
f
′
-
λ
f
′′
]
=
0
,
(41)
g
′′′
+
(
1
-
ξ
)
(
η
2
g
′′
-
ξ
∂
g
′
∂
ξ
)
+
ξ
[
(
f
+
g
)
g
′′
-
(
g
′
)
2
-
M
2
g
′
-
λ
g
′′
]
=
0
,
(42)
θ
′′
+
Pr
(
1
-
ξ
)
(
η
2
θ
′
-
ξ
∂
θ
∂
ξ
)
+
Pr
ξ
(
f
+
g
)
θ
′
=
0
,
(43)
ϕ
′′
+
Sc
(
1
-
ξ
)
(
η
2
ϕ
′
-
ξ
∂
ϕ
∂
ξ
)
+
Sc
ξ
(
f
+
g
)
ϕ
′
-
γ
Sc
ξ
ϕ
=
0
with the following boundary conditions:
(44)
f
(
ξ
,
0
)
=
g
(
ξ
,
0
)
=
0
,
f
′
(
ξ
,
0
)
=
θ
(
ξ
,
0
)
=
ϕ
(
ξ
,
0
)
=
1
,
f
′
(
ξ
,
∞
)
=
g
′
(
ξ
,
∞
)
=
θ
(
ξ
,
∞
)
=
ϕ
(
ξ
,
∞
)
=
0
,
g
′
(
ξ
,
0
)
=
c
.In the above equations prime denotes the derivative with respect toη and c the stretching parameter is a positive constant. M is the local Hartman number, λ the local porosity parameter, Sc the Schmidt number, Pr the Prandtl number, and γ the chemical reaction parameter. The initial unsteady solution can be found exactly by setting ξ
=
0 in the above equations and solving the resulting equations. The closed form analytical solutions are given by
(45)
f
(
0
,
η
)
=
η
erfc
⁡
(
η
2
)
+
2
π
[
1
-
exp
⁡
⁡
(
-
η
2
4
)
]
,
g
(
0
,
η
)
=
c
(
η
erfc
⁡
(
η
2
)
+
2
π
[
1
-
exp
⁡
⁡
(
-
η
2
4
)
]
)
,
θ
(
0
,
η
)
=
erfc
⁡
(
P
r
η
2
)
,
ϕ
(
0
,
η
)
=
erfc
⁡
(
S
c
η
2
)
.Quantities of physical interest in this problems are the skin friction coefficientsC
f
x and C
f
y in x- and y-directions, local Nusselt number Nu, and local Sherwood number Sh which are given in [22] in dimensionless form as
(46)
Re
x
1
/
2
ξ
1
/
2
C
f
x
=
-
f
′′
(
0
,
ξ
)
,
Re
y
1
/
2
ξ
1
/
2
C
f
y
=
-
g
′′
(
0
,
ξ
)
,
Nu
Re
x
-
1
/
2
ξ
1
/
2
=
-
θ
′
(
0
,
ξ
)
,
Sh
Re
x
-
1
/
2
ξ
1
/
2
=
-
ϕ
′
(
0
,
ξ
)
,
where Re
x and Re
y are the local Reynolds numbers, f
′′
(
0
,
ξ
) and g
′′
(
0
,
ξ
) are the surface shear stresses in x- and y-directions, θ
′
(
0
,
ξ
) is the surface heat transfer parameter, and ϕ
′
(
0
,
ξ
) is the surface mass transfer parameter.


### 3.1. Spectral Relaxation Method Solution

In this section we discuss the development of the spectral relaxation method to solve the system of partial differential equations (40)–(43). First, we set f
′
=
u and g
′
=
v, so that (40) and (41) become
(47)
u
′′
+
(
1
-
ξ
)
(
η
2
u
′
-
ξ
∂
u
∂
ξ
)
+
ξ
[
(
f
+
g
)
u
′
-
u
2
-
M
2
u
-
λ
u
]
=
0
,
v
′′
+
(
1
-
ξ
)
(
η
2
v
′
-
ξ
∂
v
∂
ξ
)
+
ξ
[
(
f
+
g
)
v
′
-
v
2
-
M
2
v
-
λ
v
]
=
0
.
Applying the SRM on the resulting system of nonlinear partial differential equations gives the following linear partial differential equations:
(48)
u
r
+
1
′′
+
a
1
,
r
u
r
+
1
′
+
a
2
,
r
u
r
+
1
+
a
3
,
r
=
ξ
(
1
-
ξ
)
∂
u
r
+
1
∂
ξ
,
(49)
f
r
+
1
′
=
u
r
+
1
,
(50)
v
r
+
1
′′
+
b
1
,
r
v
r
+
1
′
+
b
2
,
r
v
r
+
1
+
b
3
,
r
=
ξ
(
1
-
ξ
)
∂
v
r
+
1
∂
ξ
,
(51)
g
r
+
1
′
=
v
r
+
1
,
(52)
θ
r
+
1
′′
+
c
1
,
r
θ
r
+
1
′
=
ξ
(
1
-
ξ
)
∂
θ
r
+
1
∂
ξ
,
(53)
ϕ
r
+
1
′′
+
d
1
,
r
ϕ
r
+
1
′
+
d
2
,
r
u
r
+
1
=
ξ
(
1
-
ξ
)
∂
ϕ
r
+
1
∂
ξ
,
(54)
u
r
+
1
(
0
,
ξ
)
=
θ
r
+
1
(
0
,
ξ
)
=
ϕ
r
+
1
(
0
,
ξ
)
=
1
,
f
r
+
1
(
0
,
ξ
)
=
g
r
+
1
(
0
,
ξ
)
=
0
,
v
r
+
1
(
0
,
ξ
)
=
c
,
u
r
+
1
(
∞
,
ξ
)
=
v
r
+
1
(
∞
,
ξ
)
=
θ
r
+
1
(
∞
,
ξ
)
=
0
,
ϕ
r
+
1
(
∞
,
ξ
)
=
0
,
where
(55)
a
1
,
r
=
1
2
η
(
1
-
ξ
)
+
ξ
(
f
r
+
g
r
)
,
a
2
,
r
=
-
ξ
(
M
2
+
λ
)
,
a
3
,
r
=
-
ξ
u
r
2
,
b
1
,
r
=
a
1
,
r
,
b
2
,
r
=
a
2
,
r
,
b
3
,
r
=
-
ξ
v
r
2
,
c
1
,
r
=
Pr
(
1
2
η
(
1
-
ξ
)
+
ξ
(
f
r
+
g
r
)
)
,
d
1
,
r
=
Sc
(
1
2
η
(
1
-
ξ
)
+
ξ
(
f
r
+
g
r
)
)
,
d
2
,
r
=
-
γ
Sc
ξ
.Starting from given initial approximations, denoted byu
0
(
η
,
ξ
), f
0
(
η
,
ξ
), v
0
(
η
,
ξ
), g
0
(
η
,
ξ
), θ
0
(
η
,
ξ
), and ϕ
0
(
η
,
ξ
), (48)–(53) can be solved iteratively for the unknown functions. To solve the above decoupled system of differential equations we apply Chebyshev spectral collocation method on the space variable and finite differences in the time variable as described previously and obtain the following system of decoupled equations:
(56)
A
1
U
r
+
1
n
+
1
=
B
1
U
r
+
1
n
+
K
1
,
D
F
r
+
1
n
+
1
=
U
r
+
1
n
+
1
,
A
2
V
r
+
1
n
+
1
=
B
2
V
r
+
1
n
+
K
2
,
D
G
r
+
1
n
+
1
=
V
r
+
1
n
+
1
,
A
3
Θ
r
+
1
n
+
1
=
B
3
Θ
r
+
1
n
+
K
2
,
A
4
Φ
r
+
1
n
+
1
=
B
4
Φ
r
+
1
n
+
K
4
,
u
r
+
1
(
η
0
,
ξ
n
)
=
θ
r
+
1
(
η
0
,
ξ
n
)
=
ϕ
r
+
1
(
η
0
,
ξ
n
)
=
1
,
v
r
+
1
(
η
0
,
ξ
n
)
=
c
,
f
r
+
1
(
η
0
,
ξ
n
)
=
g
r
+
1
(
η
0
,
ξ
n
)
=
0
,
u
r
+
1
(
η
∞
,
ξ
n
)
=
v
r
+
1
(
η
∞
,
ξ
n
)
=
θ
r
+
1
(
η
∞
,
ξ
n
)
=
0
,
ϕ
r
+
1
(
η
∞
,
ξ
n
)
=
0
,
f
r
+
1
(
η
j
,
0
)
=
η
erfc
⁡
(
η
j
2
)
+
2
π
[
1
-
exp
⁡
⁡
(
-
η
j
2
4
)
]
,
g
r
+
1
(
η
j
,
0
)
=
c
(
η
erfc
⁡
(
η
j
2
)
+
2
π
[
1
-
exp
⁡
⁡
(
-
η
j
2
4
)
]
)
,
θ
r
+
1
(
η
j
,
0
)
=
erfc
⁡
(
Pr
η
j
2
)
,
ϕ
r
+
1
(
η
j
,
0
)
=
erfc
⁡
(
Sc
η
j
2
)
,
j
=
0,1
,
2
,
…
,
N
x
,
where the matrices above are defined as
(57)
A
1
=
1
2
(
D
2
+
a
1
,
r
n
+
(
1
/
2
)
D
+
a
2
,
r
)
-
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
I
,
A
2
=
1
2
(
D
2
+
b
1
,
r
n
+
(
1
/
2
)
D
+
b
2
,
r
)
-
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
I
,
A
3
=
1
2
(
D
2
+
c
1
,
r
n
+
(
1
/
2
)
D
)
-
Pr
⁡
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
I
,
A
4
=
1
2
(
D
2
+
d
1
,
r
n
+
(
1
/
2
)
D
+
d
2
,
r
)
-
Sc
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
I
,
B
1
=
-
1
2
(
D
2
+
a
1
,
r
n
+
(
1
/
2
)
D
+
a
2
,
r
)
-
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
I
,
B
2
=
-
1
2
(
D
2
+
b
1
,
r
n
+
(
1
/
2
)
D
+
b
2
,
r
)
-
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
I
,
B
3
=
-
1
2
(
D
2
+
c
1
,
r
n
+
(
1
/
2
)
D
)
-
Pr
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
I
,
B
4
=
-
1
2
(
D
2
+
d
1
,
r
n
+
(
1
/
2
)
D
+
d
2
,
r
)
-
Sc
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
I
,
K
1
=
a
3
,
r
n
+
(
1
/
2
)
,
K
2
=
b
3
,
r
n
+
(
1
/
2
)
,
K
3
=
O
¯
,
K
4
=
O
¯
.Above,U
,
F
,
V
,
G
,
Θ, and Φ are the vectors of the functions u
,
f
,
v
,
g
,
θ, and ϕ, respectively, when evaluated at the grid points η
j
(
j
=
0,1
,
…
,
N
x
).


### 3.2. Spectral Quasilinearization Method

In this section we discuss the development of the spectral quasilinearization method to solve the system of partial differential equations (40)–(43). First, we set f
′
=
u and g
′
=
v, so that equations (40) and (41) become
(58)
u
′′
+
(
1
-
ξ
)
(
η
2
u
′
-
ξ
∂
u
∂
ξ
)
+
ξ
[
(
f
+
g
)
u
′
-
u
2
-
M
2
u
-
λ
u
]
=
0
,
v
′′
+
(
1
-
ξ
)
(
η
2
v
′
-
ξ
∂
v
∂
ξ
)
+
ξ
[
(
f
+
g
)
v
′
-
v
2
-
M
2
v
-
λ
v
]
=
0
.
Applying the SQLM on the resulting system of nonlinear partial differential equations gives the following linear partial differential equations:
(59)
u
r
+
1
′′
+
a
1
,
r
u
r
+
1
′
+
a
2
,
r
u
r
+
1
+
a
3
,
r
f
r
+
1
+
a
4
,
r
g
r
+
1
+
a
5
,
r
=
ξ
(
1
-
ξ
)
∂
u
r
+
1
∂
ξ
,
(60)
f
r
+
1
′
=
u
r
+
1
,
(61)
v
r
+
1
′′
+
b
1
,
r
v
r
+
1
′
+
b
2
,
r
v
r
+
1
+
b
3
,
r
f
r
+
1
+
b
4
,
r
g
r
+
1
+
b
5
,
r
=
ξ
(
1
-
ξ
)
∂
v
r
+
1
∂
ξ
,
(62)
g
r
+
1
′
=
v
r
+
1
,
(63)
θ
r
+
1
′′
+
c
1
,
r
θ
r
+
1
′
+
c
2
,
r
f
r
+
1
+
c
3
,
r
g
r
+
1
+
c
4
,
r
=
ξ
(
1
-
ξ
)
∂
θ
r
+
1
∂
ξ
,
(64)
ϕ
r
+
1
′′
+
d
1
,
r
ϕ
r
+
1
′
+
d
2
,
r
u
r
+
1
+
d
3
,
r
f
r
+
1
+
d
4
,
r
g
r
+
1
+
d
5
,
r
=
ξ
(
1
-
ξ
)
∂
ϕ
r
+
1
∂
ξ
,
(65)
u
r
+
1
(
0
,
ξ
)
=
θ
r
+
1
(
0
,
ξ
)
=
ϕ
r
+
1
(
0
,
ξ
)
=
1
,
f
r
+
1
(
0
,
ξ
)
=
g
r
+
1
(
0
,
ξ
)
=
0
,
v
r
+
1
(
0
,
ξ
)
=
c
,
u
r
+
1
(
∞
,
ξ
)
=
v
r
+
1
(
∞
,
ξ
)
=
θ
r
+
1
(
∞
,
ξ
)
=
0
,
ϕ
r
+
1
(
∞
,
ξ
)
=
0
,
where
(66)
a
1
,
r
=
1
2
η
(
1
-
ξ
)
+
ξ
(
f
r
+
g
r
)
,
a
2
,
r
=
-
ξ
(
2
u
r
+
M
2
+
λ
)
,
a
3
,
r
=
a
4
,
r
=
ξ
u
r
′
,
a
5
,
r
=
ξ
(
u
r
2
-
(
f
r
+
g
r
)
u
r
′
)
,
b
1
,
r
=
a
1
,
r
,
b
2
,
r
=
-
ξ
(
2
v
r
+
M
2
+
λ
)
,
b
3
,
r
=
b
4
,
r
=
ξ
v
r
′
,
b
5
,
r
=
ξ
(
v
r
2
-
(
f
r
+
g
r
)
v
r
′
)
,
c
1
,
r
=
Pr
(
1
2
η
(
1
-
ξ
)
+
ξ
(
f
r
+
g
r
)
)
,
c
2
,
r
=
c
3
,
r
=
ξ
Pr
θ
r
′
,
c
4
,
r
=
ξ
Pr
⁡
(
f
r
+
g
r
)
θ
r
′
,
d
1
,
r
=
Sc
(
1
2
η
(
1
-
ξ
)
+
ξ
(
f
r
+
g
r
)
)
,
d
2
,
r
=
-
γ
Sc
ξ
,
d
3
,
r
=
d
4
,
r
=
ξ
Sc
ϕ
r
′
,
d
5
,
r
=
ξ
Sc
(
f
r
+
g
r
)
ϕ
r
′
.
Starting from given initial approximations, denoted by u
0
(
η
,
ξ
), f
0
(
η
,
ξ
), v
0
(
η
,
ξ
), g
0
(
η
,
ξ
), θ
0
(
η
,
ξ
), and ϕ
0
(
η
,
ξ
), equations (59)–(64) can be solved iteratively for the unknown functions. To solve the above decoupled system of differential equations we apply Chebyshev spectral collocation method on the space variable and finite differences in the time variable as described previously.


## 3.1. Spectral Relaxation Method Solution

In this section we discuss the development of the spectral relaxation method to solve the system of partial differential equations (40)–(43). First, we set f
′
=
u and g
′
=
v, so that (40) and (41) become
(47)
u
′′
+
(
1
-
ξ
)
(
η
2
u
′
-
ξ
∂
u
∂
ξ
)
+
ξ
[
(
f
+
g
)
u
′
-
u
2
-
M
2
u
-
λ
u
]
=
0
,
v
′′
+
(
1
-
ξ
)
(
η
2
v
′
-
ξ
∂
v
∂
ξ
)
+
ξ
[
(
f
+
g
)
v
′
-
v
2
-
M
2
v
-
λ
v
]
=
0
.
Applying the SRM on the resulting system of nonlinear partial differential equations gives the following linear partial differential equations:
(48)
u
r
+
1
′′
+
a
1
,
r
u
r
+
1
′
+
a
2
,
r
u
r
+
1
+
a
3
,
r
=
ξ
(
1
-
ξ
)
∂
u
r
+
1
∂
ξ
,
(49)
f
r
+
1
′
=
u
r
+
1
,
(50)
v
r
+
1
′′
+
b
1
,
r
v
r
+
1
′
+
b
2
,
r
v
r
+
1
+
b
3
,
r
=
ξ
(
1
-
ξ
)
∂
v
r
+
1
∂
ξ
,
(51)
g
r
+
1
′
=
v
r
+
1
,
(52)
θ
r
+
1
′′
+
c
1
,
r
θ
r
+
1
′
=
ξ
(
1
-
ξ
)
∂
θ
r
+
1
∂
ξ
,
(53)
ϕ
r
+
1
′′
+
d
1
,
r
ϕ
r
+
1
′
+
d
2
,
r
u
r
+
1
=
ξ
(
1
-
ξ
)
∂
ϕ
r
+
1
∂
ξ
,
(54)
u
r
+
1
(
0
,
ξ
)
=
θ
r
+
1
(
0
,
ξ
)
=
ϕ
r
+
1
(
0
,
ξ
)
=
1
,
f
r
+
1
(
0
,
ξ
)
=
g
r
+
1
(
0
,
ξ
)
=
0
,
v
r
+
1
(
0
,
ξ
)
=
c
,
u
r
+
1
(
∞
,
ξ
)
=
v
r
+
1
(
∞
,
ξ
)
=
θ
r
+
1
(
∞
,
ξ
)
=
0
,
ϕ
r
+
1
(
∞
,
ξ
)
=
0
,
where
(55)
a
1
,
r
=
1
2
η
(
1
-
ξ
)
+
ξ
(
f
r
+
g
r
)
,
a
2
,
r
=
-
ξ
(
M
2
+
λ
)
,
a
3
,
r
=
-
ξ
u
r
2
,
b
1
,
r
=
a
1
,
r
,
b
2
,
r
=
a
2
,
r
,
b
3
,
r
=
-
ξ
v
r
2
,
c
1
,
r
=
Pr
(
1
2
η
(
1
-
ξ
)
+
ξ
(
f
r
+
g
r
)
)
,
d
1
,
r
=
Sc
(
1
2
η
(
1
-
ξ
)
+
ξ
(
f
r
+
g
r
)
)
,
d
2
,
r
=
-
γ
Sc
ξ
.Starting from given initial approximations, denoted byu
0
(
η
,
ξ
), f
0
(
η
,
ξ
), v
0
(
η
,
ξ
), g
0
(
η
,
ξ
), θ
0
(
η
,
ξ
), and ϕ
0
(
η
,
ξ
), (48)–(53) can be solved iteratively for the unknown functions. To solve the above decoupled system of differential equations we apply Chebyshev spectral collocation method on the space variable and finite differences in the time variable as described previously and obtain the following system of decoupled equations:
(56)
A
1
U
r
+
1
n
+
1
=
B
1
U
r
+
1
n
+
K
1
,
D
F
r
+
1
n
+
1
=
U
r
+
1
n
+
1
,
A
2
V
r
+
1
n
+
1
=
B
2
V
r
+
1
n
+
K
2
,
D
G
r
+
1
n
+
1
=
V
r
+
1
n
+
1
,
A
3
Θ
r
+
1
n
+
1
=
B
3
Θ
r
+
1
n
+
K
2
,
A
4
Φ
r
+
1
n
+
1
=
B
4
Φ
r
+
1
n
+
K
4
,
u
r
+
1
(
η
0
,
ξ
n
)
=
θ
r
+
1
(
η
0
,
ξ
n
)
=
ϕ
r
+
1
(
η
0
,
ξ
n
)
=
1
,
v
r
+
1
(
η
0
,
ξ
n
)
=
c
,
f
r
+
1
(
η
0
,
ξ
n
)
=
g
r
+
1
(
η
0
,
ξ
n
)
=
0
,
u
r
+
1
(
η
∞
,
ξ
n
)
=
v
r
+
1
(
η
∞
,
ξ
n
)
=
θ
r
+
1
(
η
∞
,
ξ
n
)
=
0
,
ϕ
r
+
1
(
η
∞
,
ξ
n
)
=
0
,
f
r
+
1
(
η
j
,
0
)
=
η
erfc
⁡
(
η
j
2
)
+
2
π
[
1
-
exp
⁡
⁡
(
-
η
j
2
4
)
]
,
g
r
+
1
(
η
j
,
0
)
=
c
(
η
erfc
⁡
(
η
j
2
)
+
2
π
[
1
-
exp
⁡
⁡
(
-
η
j
2
4
)
]
)
,
θ
r
+
1
(
η
j
,
0
)
=
erfc
⁡
(
Pr
η
j
2
)
,
ϕ
r
+
1
(
η
j
,
0
)
=
erfc
⁡
(
Sc
η
j
2
)
,
j
=
0,1
,
2
,
…
,
N
x
,
where the matrices above are defined as
(57)
A
1
=
1
2
(
D
2
+
a
1
,
r
n
+
(
1
/
2
)
D
+
a
2
,
r
)
-
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
I
,
A
2
=
1
2
(
D
2
+
b
1
,
r
n
+
(
1
/
2
)
D
+
b
2
,
r
)
-
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
I
,
A
3
=
1
2
(
D
2
+
c
1
,
r
n
+
(
1
/
2
)
D
)
-
Pr
⁡
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
I
,
A
4
=
1
2
(
D
2
+
d
1
,
r
n
+
(
1
/
2
)
D
+
d
2
,
r
)
-
Sc
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
I
,
B
1
=
-
1
2
(
D
2
+
a
1
,
r
n
+
(
1
/
2
)
D
+
a
2
,
r
)
-
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
I
,
B
2
=
-
1
2
(
D
2
+
b
1
,
r
n
+
(
1
/
2
)
D
+
b
2
,
r
)
-
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
I
,
B
3
=
-
1
2
(
D
2
+
c
1
,
r
n
+
(
1
/
2
)
D
)
-
Pr
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
I
,
B
4
=
-
1
2
(
D
2
+
d
1
,
r
n
+
(
1
/
2
)
D
+
d
2
,
r
)
-
Sc
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
I
,
K
1
=
a
3
,
r
n
+
(
1
/
2
)
,
K
2
=
b
3
,
r
n
+
(
1
/
2
)
,
K
3
=
O
¯
,
K
4
=
O
¯
.Above,U
,
F
,
V
,
G
,
Θ, and Φ are the vectors of the functions u
,
f
,
v
,
g
,
θ, and ϕ, respectively, when evaluated at the grid points η
j
(
j
=
0,1
,
…
,
N
x
).


## 3.2. Spectral Quasilinearization Method

In this section we discuss the development of the spectral quasilinearization method to solve the system of partial differential equations (40)–(43). First, we set f
′
=
u and g
′
=
v, so that equations (40) and (41) become
(58)
u
′′
+
(
1
-
ξ
)
(
η
2
u
′
-
ξ
∂
u
∂
ξ
)
+
ξ
[
(
f
+
g
)
u
′
-
u
2
-
M
2
u
-
λ
u
]
=
0
,
v
′′
+
(
1
-
ξ
)
(
η
2
v
′
-
ξ
∂
v
∂
ξ
)
+
ξ
[
(
f
+
g
)
v
′
-
v
2
-
M
2
v
-
λ
v
]
=
0
.
Applying the SQLM on the resulting system of nonlinear partial differential equations gives the following linear partial differential equations:
(59)
u
r
+
1
′′
+
a
1
,
r
u
r
+
1
′
+
a
2
,
r
u
r
+
1
+
a
3
,
r
f
r
+
1
+
a
4
,
r
g
r
+
1
+
a
5
,
r
=
ξ
(
1
-
ξ
)
∂
u
r
+
1
∂
ξ
,
(60)
f
r
+
1
′
=
u
r
+
1
,
(61)
v
r
+
1
′′
+
b
1
,
r
v
r
+
1
′
+
b
2
,
r
v
r
+
1
+
b
3
,
r
f
r
+
1
+
b
4
,
r
g
r
+
1
+
b
5
,
r
=
ξ
(
1
-
ξ
)
∂
v
r
+
1
∂
ξ
,
(62)
g
r
+
1
′
=
v
r
+
1
,
(63)
θ
r
+
1
′′
+
c
1
,
r
θ
r
+
1
′
+
c
2
,
r
f
r
+
1
+
c
3
,
r
g
r
+
1
+
c
4
,
r
=
ξ
(
1
-
ξ
)
∂
θ
r
+
1
∂
ξ
,
(64)
ϕ
r
+
1
′′
+
d
1
,
r
ϕ
r
+
1
′
+
d
2
,
r
u
r
+
1
+
d
3
,
r
f
r
+
1
+
d
4
,
r
g
r
+
1
+
d
5
,
r
=
ξ
(
1
-
ξ
)
∂
ϕ
r
+
1
∂
ξ
,
(65)
u
r
+
1
(
0
,
ξ
)
=
θ
r
+
1
(
0
,
ξ
)
=
ϕ
r
+
1
(
0
,
ξ
)
=
1
,
f
r
+
1
(
0
,
ξ
)
=
g
r
+
1
(
0
,
ξ
)
=
0
,
v
r
+
1
(
0
,
ξ
)
=
c
,
u
r
+
1
(
∞
,
ξ
)
=
v
r
+
1
(
∞
,
ξ
)
=
θ
r
+
1
(
∞
,
ξ
)
=
0
,
ϕ
r
+
1
(
∞
,
ξ
)
=
0
,
where
(66)
a
1
,
r
=
1
2
η
(
1
-
ξ
)
+
ξ
(
f
r
+
g
r
)
,
a
2
,
r
=
-
ξ
(
2
u
r
+
M
2
+
λ
)
,
a
3
,
r
=
a
4
,
r
=
ξ
u
r
′
,
a
5
,
r
=
ξ
(
u
r
2
-
(
f
r
+
g
r
)
u
r
′
)
,
b
1
,
r
=
a
1
,
r
,
b
2
,
r
=
-
ξ
(
2
v
r
+
M
2
+
λ
)
,
b
3
,
r
=
b
4
,
r
=
ξ
v
r
′
,
b
5
,
r
=
ξ
(
v
r
2
-
(
f
r
+
g
r
)
v
r
′
)
,
c
1
,
r
=
Pr
(
1
2
η
(
1
-
ξ
)
+
ξ
(
f
r
+
g
r
)
)
,
c
2
,
r
=
c
3
,
r
=
ξ
Pr
θ
r
′
,
c
4
,
r
=
ξ
Pr
⁡
(
f
r
+
g
r
)
θ
r
′
,
d
1
,
r
=
Sc
(
1
2
η
(
1
-
ξ
)
+
ξ
(
f
r
+
g
r
)
)
,
d
2
,
r
=
-
γ
Sc
ξ
,
d
3
,
r
=
d
4
,
r
=
ξ
Sc
ϕ
r
′
,
d
5
,
r
=
ξ
Sc
(
f
r
+
g
r
)
ϕ
r
′
.
Starting from given initial approximations, denoted by u
0
(
η
,
ξ
), f
0
(
η
,
ξ
), v
0
(
η
,
ξ
), g
0
(
η
,
ξ
), θ
0
(
η
,
ξ
), and ϕ
0
(
η
,
ξ
), equations (59)–(64) can be solved iteratively for the unknown functions. To solve the above decoupled system of differential equations we apply Chebyshev spectral collocation method on the space variable and finite differences in the time variable as described previously.


## 4. Results and Discussion

In this section we present the SRM and SQLM results for the two examples described above. Numerical simulations were carried out to obtain approximate numerical values of the quantities of physical interest, namely, the surface shear stresses, surface heat transfer, and the surface mass transfer parameter. In all the spectral method based numerical simulations a finite computational domain of extentη
∞
=
20 was chosen in the η-direction. Through numerical experimentation, this value was found to give accurate results for all the selected governing physical parameters used in the generation of results. Increasing the value of η did not change the results to a significant extent. The number of collocation points used in the spectral method discretization was N
x
=
60 in all cases. We remark that both the SRM and SQLM algorithms are based on the computation of the value of some quantity, say F
r
+
1
n
+
1, at each time step. This is achieved by iterating using the relaxation method or the quasilinearization method using a known value at the previous time step n as initial approximation. The iteration calculations are carried until some desired tolerance level, ε, is attained. In this study, the tolerance level was set to be ε
=
1
0
-
8. The tolerance level is defined as the maximum values of the infinity norm of the difference between the values of the calculated quantities and its first two derivatives at successive iterations. For example, in calculating F
r
+
1
n
+
1, the tolerance level and convergence criteria are defined as
(67)
max
⁡
⁡
{
∥
F
r
+
1
n
+
1
-
F
r
n
+
1
∥
∞
,
∥
U
r
+
1
n
+
1
-
U
r
n
+
1
∥
∞
,
∥
W
r
+
1
n
+
1
-
W
r
n
+
1
∥
∞
}
<
ε
,
where U
=
F
′ and W
=
F
′′. The accuracy of the computed SRM and SQLM approximate results was confirmed against numerical results obtained by using the popular Keller-box implicit finite difference method as described by [2]. The Keller-box method has been reported in literature to be accurate, fast, and easier to program for boundary layer flow problems. The algorithm begins with the reduction of the governing nonlinear PDEs into a system of first order equations that are discretized using central differences. The nonlinear algebraic difference equations are linearised using Newton’s method and written in matrix-vector form. The linear matrix systems are solved in an efficient manner using block-tridiagonal-elimination technique. The grid spacing in both the η-direction and ξ-direction is carefully selected to ensure that the Keller-box computations yield consistent results for the governing velocity and temperature distributions to a convergence level of at least 1
0
-
8.Tables1 and 2 give the approximate numerical values of the skin friction f
′′
(
0
,
ξ
) for various step sizes Δ
ξ, computed using the SRM and SQLM, respectively, for Example 1. The tables also give the total computational time for the integration in the whole time domain to be completed. We remark that the results were computed using the same number of collocation points N
x and the same η
∞. Reducing the step size Δ
ξ improves the accuracy of the results until the results are consistent to within eight decimal digits. The results displayed in the tables are quite revealing in several ways. First, it is clear from the comparison of the computational run times that the SRM takes less computer time than the SQLM. The results also indicate that the SRM converges more rapidly than the SQLM results when the step size Δ
ξ is reduced. Full convergence to within eight decimal digits is reached when Δ
ξ is at least 0.0005 in the SRM compared to Δ
ξ
=
0.0001 in the case of the SQLM. This means that the SRM converges faster than the SQLM with a decrease in the step size Δ
ξ. Furthermore, there is good agreement between the two sets of results when Δ
ξ is very small for all values of ξ. The apparent superiority of the SRM in terms of computational efficiency and accuracy when compared to the SQLM may be explained by the fact that the SRM algorithm reduced a coupled system of equations into smaller sequences of decoupled equations which are solved one after the other. Smaller sized matrices are less susceptible to round-off errors and ill-conditioning and take less computer time to invert.Table 1
Skin frictionf
′
′
(
0
,
ξ
) for various values of ξ computed using SRM.
ξ
Δ
ξ
0.01
0.001
0.0005
0.0002
0.0001
0.1
−0.61046835
−0.61046762
−0.61046761
−0.61046761
−0.61046761
0.3
−0.70126751
−0.70126681
−0.70126680
−0.70126680
−0.70126680
0.5
−0.78982903
−0.78982837
−0.78982837
−0.78982837
−0.78982837
0.7
−0.87626715
−0.87626654
−0.87626653
−0.87626653
−0.87626653
0.9
−0.96053875
−0.96053800
−0.96053800
−0.96053800
−0.96053800
Time
0.66
3.01
5.68
12.54
23.83Table 2
Skin frictionf
′
′
(
0
,
ξ
) for various values of ξ computed using SQLM.
ξ
Δ
ξ
0.01
0.001
0.0005
0.0002
0.0001
0.1
−0.61045544
−0.61046674
−0.61046742
−0.61046758
−0.61046761
0.3
−0.70126943
−0.70126664
−0.70126676
−0.70126679
−0.70126680
0.5
−0.78981759
−0.78982831
−0.78982835
−0.78982836
−0.78982837
0.7
−0.87625663
−0.87626652
−0.87626653
−0.87626653
−0.87626653
0.9
−0.96053069
−0.96053800
−0.96053800
−0.96053800
−0.96053800
Time
0.60
3.29
6.62
20.37
53.92Table3 gives a comparison of the amount of time it takes for each method to generate numerical solutions that have converged to within 1
0
-
8 at selected time levels. As can be seen from Table 3, the SRM is much faster compared to the SQLM in computing the numerical solutions. For very small time steps the SRM appears to be at least twice as fast as the SQLM. The results from Tables 1, 2, and 3 indicate that the SRM is much more computationally efficient and gives more accurate results than the SQLM under the same conditions.Table 3
SRM and SQLM computational times for Example 1.
ξ
∖
Δ
ξ
0.01
0.001
0.0005
0.0002
0.0001
Spectral relaxation method
0.1
0.07
0.22
0.43
0.95
1.93
0.3
0.13
0.65
1.38
3.18
6.01
0.5
0.20
1.14
2.26
5.48
10.69
0.7
0.27
1.72
3.14
7.70
15.08
0.9
0.35
2.30
4.25
10.05
19.57
Spectral qausilinearisation method
0.1
0.03
0.41
0.76
2.27
4.72
0.3
0.11
1.26
2.28
6.78
14.67
0.5
0.20
2.20
3.92
11.01
24.19
0.7
0.29
3.14
5.47
15.25
33.29
0.9
0.39
4.09
7.09
20.11
42.65Figures1 and 2 display the number of iterations required to yield converged results to within the tolerance level of 1
0
-
8 against the time for the SRM and SQLM, respectively. The results are given for different values of the number of grid points N
t. It can be seen from Figure 1 that more iterations are required to give the converged results when the number of grid points is small. For larger values of N
t, between four and six iterations are required in the range 0
≤
0.9. For ξ near 1, the number of iterations required increases. The trends in the results for SQLM, as seen in Figure 2, are similar to those for the SRM. However, in the case of the SQLM, full convergence is achieved with only four iterations in a wider range of ξ and for smaller values of N
t than the SRM. This observation indicates that the actual convergence rates (with an increase in iterations) of the SQLM are greater than those of the SRM.Figure 1
Variation of the SRM iterations with timeξ in Example 1.Figure 2
Variation of the SQLM iterations with timeξ in Example 1.Tables4 and 5 give a comparison of the SRM and SQLM approximate numerical solutions, respectively, against the Keller-box results for the skin frictions, surface heat transfer parameter, and surface mass transfer parameters. We remark that the results reported in Tables 4 and 5 were generated using a tolerance level of 1
0
-
7 in the SRM, SQLM, and Keller-box implementations. It can be noted from Tables 4 and 5 that both the SRM and SQLM converge to the Keller-box results when Δ
ξ is sufficiently small. The convergence to six-decimal-digit accurate results is more or less the same for both SRM and SQLM schemes. We remark that the Keller-box results given in Tables 4 and 5 were calculated using nonuniform step size in the η-direction and a uniform step size Δ
ξ
=
0.0005 in the ξ-direction. Using a nonuniform grid size significantly improves the computation time of the Keller-box method. Thus, to speed up the computation times for the Keller-box method, computations were carried out with an initial step size of Δ
η
0
=
0.001. This was gradually increased by the variable grid parameter (VGP) factor of 1.005 between successive grid points in accordance with the formula η
j
=
η
j
-
1
+
VGP
×
Δ
η
j
-
1 for j
=
1,2
,
…
,
J (where J is the number of grid points in the η direction). The value of η
∞ was fixed at η
∞
=
10 for the Keller-box implementation.Table 4
Comparison between the SRM and Keller-Box approximate numerical values forf
′
′
(
0
,
ξ
), g
′
′
(
0
,
ξ
), θ
′
(
0
,
ξ
), and ϕ
′
(
0
,
ξ
) when λ
=
0.5, M
=
2, c
=
0.5, Sc
=
γ
=
1, and Pr
=
1.5 in Example 2.
Spectral relaxation method
Keller-box
ξ
∖
Δ
ξ
0.01
0.002
0.001
0.0005
0.0005
f
′
′
(
0
,
ξ
)
0.1
−0.851309
−0.851259
−0.851258
−0.851257
−0.851257
0.3
−1.316738
−1.316706
−1.316705
−1.316705
−1.316705
0.5
−1.685327
−1.685307
−1.685306
−1.685306
−1.685306
0.7
−1.992622
−1.992609
−1.992608
−1.992608
−1.992608
0.9
−2.259344
−2.259335
−2.259335
−2.259335
−2.259335
g
′
′
(
0
,
ξ
)
0.1
−0.417173
−0.417151
−0.417151
−0.417150
−0.417150
0.3
−0.639617
−0.639603
−0.639602
−0.639602
−0.639602
0.5
−0.817659
−0.817649
−0.817649
−0.817649
−0.817649
0.7
−0.966610
−0.966604
−0.966603
−0.966603
−0.966604
0.9
−1.095987
−1.095983
−1.095983
−1.095983
−1.095983
θ
′
(
0
,
ξ
)
0.1
−0.710885
−0.710882
−0.710882
−0.710882
−0.710882
0.3
−0.742845
−0.742842
−0.742842
−0.742842
−0.742843
0.5
−0.765247
−0.765244
−0.765244
−0.765244
−0.765244
0.7
−0.777274
−0.777270
−0.777270
−0.777270
−0.777270
0.9
−0.770819
−0.770807
−0.770807
−0.770807
−0.770807
ϕ
′
(
0
,
ξ
)
0.1
−0.634447
−0.634443
−0.634443
−0.634443
−0.634444
0.3
−0.766870
−0.766867
−0.766867
−0.766867
−0.766867
0.5
−0.891209
−0.891207
−0.891207
−0.891207
−0.891207
0.7
−1.010046
−1.010045
−1.010045
−1.010045
−1.010045
0.9
−1.125550
−1.125549
−1.125549
−1.125549
−1.125549Table 5
Comparison between the SRM and Keller-Box approximate numerical values forf
′
′
(
0
,
ξ
), g
′
′
(
0
,
ξ
), θ
′
(
0
,
ξ
), and ϕ
′
(
0
,
ξ
) when λ
=
0.5, M
=
2, c
=
0.5, Sc
=
γ
=
1, and Pr
=
1.5 in Example 2.
Spectral quasilinearisation method
Keller-Box
ξ
∖
Δ
ξ
0.01
0.002
0.001
0.0005
0.0005
f
′
′
(
0
,
ξ
)
0.1
−0.851309
−0.851259
−0.851258
−0.851257
−0.851257
0.3
−1.316738
−1.316706
−1.316705
−1.316705
−1.316705
0.5
−1.685327
−1.685307
−1.685306
−1.685306
−1.685306
0.7
−1.992622
−1.992609
−1.992608
−1.992608
−1.992608
0.9
−2.259344
−2.259335
−2.259335
−2.259335
−2.259335
g
′
′
(
0
,
ξ
)
0.1
−0.417173
−0.417151
−0.417151
−0.417150
−0.417150
0.3
−0.639617
−0.639603
−0.639602
−0.639602
−0.639602
0.5
−0.817659
−0.817649
−0.817649
−0.817649
−0.817649
0.7
−0.966610
−0.966604
−0.966603
−0.966603
−0.966604
0.9
−1.095987
−1.095983
−1.095983
−1.095983
−1.095983
θ
′
(
0
,
ξ
)
0.1
−0.710885
−0.710882
−0.710882
−0.710882
−0.710882
0.3
−0.742845
−0.742842
−0.742842
−0.742842
−0.742843
0.5
−0.765247
−0.765244
−0.765244
−0.765244
−0.765244
0.7
−0.777274
−0.777270
−0.777270
−0.777270
−0.777270
0.9
−0.770819
−0.770807
−0.770807
−0.770807
−0.770807
ϕ
′
(
0
,
ξ
)
0.1
−0.634447
−0.634443
−0.634443
−0.634443
−0.634444
0.3
−0.766870
−0.766867
−0.766867
−0.766867
−0.766867
0.5
−0.891209
−0.891207
−0.891207
−0.891207
−0.891207
0.7
−1.010046
−1.010045
−1.010045
−1.010045
−1.010045
0.9
−1.125550
−1.125549
−1.125549
−1.125549
−1.125549A comparison of the computational times between the SRM, SQLM, and Keller-box method is given in Table6 for the computation of results that converge to within six decimal digits. It can be seen from the table that there is a substantial difference in the computation times of the three methods with the SRM being at least four times faster than the SQLM and the Keller-box being the slowest method. The demonstrated speed of the spectral method based methods is primarily due to the intrinsic property of the spectral collocation method to be able to give accurate approximate results using only a few grid points. Only 60 collocation points were used to generate results that converge to at least six decimal digits of accuracy. On the other hand, the Keller-box required a lot more grid points in the η direction to give the same amount of accuracy. The apparent computational speed of the SRM over the SQLM is in accord with the observation made earlier in the case of the one-equation system. As can be seen in Table 6, the superiority in computational efficiency of the SRM over the SQLM is much more pronounced in the current example that involves a system of four coupled equations. Thus, the SRM is a better alternative method that can be used to obtain numerical solutions of systems of PDEs arising in boundary layer flow problems.Table 6
SRM and SQLM computational times for Example 2.
ξ
∖
Δ
ξ
0.01
0.001
0.0005
0.0002
0.0001
Spectral relaxation method
0.1
0.13
0.91
1.52
3.92
8.52
0.3
0.37
3.02
5.68
12.07
26.13
0.5
0.62
5.24
10.09
21.32
43.58
0.7
0.91
7.40
14.54
32.06
62.84
0.9
1.20
9.97
18.90
42.86
87.48
Spectral qausilinearisation method
0.1
0.57
4.53
8.76
21.19
43.05
0.3
1.50
14.18
27.28
69.15
140.01
0.5
2.43
23.33
45.72
119.98
235.45
0.7
3.39
32.57
64.30
168.09
326.96
0.9
4.66
41.87
83.24
216.00
418.40
Keller-box method
0.1
5.30
47.67
100.59
0.3
17.62
143.90
294.22
0.5
30.27
241.03
486.32
0.7
42.94
337.97
677.87
0.9
55.68
464.56
900.30Figures3 and 4 show the variation of the SRM and SQLM iterations, respectively, over time for different values of N
t. The indicated number of iterations is the total number of iterations required to obtain results that are consistent to within a tolerance level of 1
0
-
6. It can be noted from Figure 3 that the total number of iterations required for the SRM is between 4 and 9 with the total iterations increasing as ξ tends to 1. In contrast, the range of the required number of iterations is 3 to 5 in the case of the SQLM. Thus the convergence rate of the SQLM in terms of iterations is higher than that of the SRM.Figure 3
Variation of the SRM iterations with timeξ in Example 2.Figure 4
Variation of the SQLM iterations with timeξ in Example 2.


## 5. Conclusion

In this paper, we investigated the application of the spectral relaxation method (SRM) and spectral quasilinearisation method (SQLM) in the solution of unsteady boundary layer flows that are described by systems of coupled nonlinear partial differential equations. We considered the model problems of unsteady boundary layer flow caused by an impulsively stretching sheet and the unsteady MHD flow and mass transfer in a porous space. The purpose of this study was to establish the applicability of the SRM, for the first time, to systems of PDEs that model unsteady boundary layer flows. The investigation also sought to assess the accuracy and efficiency of the SRM compared to the SQLM and Keller-box method.The most obvious finding to emerge from this study is that the SRM is significantly more computationally efficient than the SQLM which, in turn, is faster than the Keller-box method. For sufficiently small step-sizes, all the three methods yield results that are consistent to within a given tolerance level. The SRM was observed to convergence faster than the SQLM with a reduction of the step size. It is this feature that makes the SRM computationally efficient as accurate results are obtained using fewer grid points in the time direction. In addition, the SRM algorithm involves the solution of a sequence of smaller sized matrix equations compared to the SQLM. The numerical results presented in this study clearly demonstrate the potential of the SRM scheme for the simulation of numerical solutions of the class of unsteady boundary layer flows equations related to the model equations discussed in this study. The evidence of the accuracy and efficiency of the SRM from this study suggests that the method can be used as a more practical tool for solving unsteady boundary layer flows and for validating the results generated using other numerical methods in the solution of similar boundary layer flow equations. The presented SRM approach adds to a growing body of literature on practical numerical methods for solving complex nonlinear PDEs in some fluid mechanics applications.



---

*Source: 341964-2014-06-18.xml*