# The Spectral Homotopy Analysis Method Extended to Systems of Partial Differential Equations

Authors: S. S. Motsa; F. G. Awad; Z. G. Makukula; P. Sibanda
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

The spectral homotopy analysis method is extended to solutions of systems of nonlinear partial differential equations. The SHAM has previously been successfully used to find solutions of nonlinear ordinary differential equations. We solve the nonlinear system of partial differential equations that model the unsteady nonlinear convective flow caused by an impulsively stretching sheet. The numerical results generated using the spectral homotopy analysis method were compared with those found using the spectral quasilinearisation method (SQLM) and the two results were in good agreement.


---

## Body

## 1. Introduction

The current study serves to demonstrate the extension of the spectral homotopy analysis method (SHAM) to systems of nonlinear partial differential equations. The method, originally introduced by Motsa et al. [1, 2], was used to solve nonlinear ordinary differential equations. A recent modification (see Motsa [3]) allows the method to be used to find solutions of nonlinear PDEs. A nonlinear system of partial differential equations that describe unsteady nonlinear convective flow caused by an impulsively stretched plate is used as the test problem in the study. The spectral quasilinearisation method (SQLM), an adaptation of the quasilinearisation method (Bellman and Kalaba [4]), is used as a benchmark to prove the accuracy of the spectral homotopy analysis method.The test equations for this study are obtained by considering unsteady nonlinear convection flow over an impulsively stretched flat surface. Fluid flow over stretching surfaces is important in many practical applications such as extrusion of plastic sheets, paper production, glass blowing, metal spinning, and drawing plastic films. The quality of the final product depends on the rate of heat transfer at the stretching surface. In a previous study in this field, Kumari et al. [5] used the Keller box method and the Nakamura method to investigate the problem of heat transfer in the unsteady free convection flow over a continuous moving vertical sheet in an ambient fluid. Ishak et al. [6] investigated theoretically the unsteady mixed convection boundary layer flow and heat transfer due to a stretching vertical surface in a quiescent viscous and incompressible fluid. Further, Pop and Na [7] and Wang et al. [8] dealt with the unsteady boundary layer flow due to impulsive flow starting from rest of a stretching sheet in a viscous fluid. A numerical solution for unsteady mixed convection boundary layer nanofluid flow and heat transfer due to a stretching vertical sheet was presented by Mahdy [9]. The equations were solved using the implicit finite difference method. Bachok et al. [10] studied the flow and heat transfer problem due to the unsteady, two-dimensional laminar flow of a viscous nanofluid caused by a permeable stretching/shrinking sheet in a quiescent fluid. Numerical solutions of the transformed governing equations were obtained using a shooting method. Sharma et al. [11] used the finite element method to find numerical solutions of the flow and heat transfer problem over a stretching sheet immersed in a nanofluid, with velocity slip at the boundary.In fluid flow problems involving heat transfer, it may be essential to consider a nonlinear relationship between density and temperature. Thermal stratification and heat released by the viscous dissipation triggers some changes in density gradients. To the best of our knowledge there is no study that has been conducted to discuss the effects of nonlinear density on unsteady nonlinear convection flow over an impulsively stretched flat surface. Accordingly, the main aim of the study is to investigate the nonlinear convective flow over an impulsively stretched flat surface under the nonlinear density-temperature relationship. Using the boundary layer approximation, the unsteady momentum, heat, and mass transfer equations are transferred to nonlinear partial differential equations form and solved using the spectral homotopy analysis method (SHAM). This work is the first attempt at applying the SHAM to system of coupled nonlinear PDEs. The SHAM results presented in this study are validated against results from a quasilinearisation based spectral collocation method.


## 2. Governing Equations

Consider the problem of unsteady nonlinear convection of a fluid over a stretching flat plate. Initially(
t
=
0
), both the fluid and stretching plate are kept at a constant temperature T
w and concentration C
∞ where T
w
>
T
∞ is for a heated plate and T
w
<
T
∞ corresponding to a cooled plate. We assume that at t
=
0 the velocity of the stretching plate is u
w
=
a
x, where  a is a positive constant. From the Boussinesq approximation, density is related to temperature and the concentration by the equation
(1)
ρ
=
ρ
0
[
1
-
β
T
(
T
-
T
∞
)
+
β
C
(
C
-
C
∞
)
]
.
In the case of thermal stratification and heat released by viscous dissipation, wall jet like profiles induce significant changes in density gradients, and the density depends on the temperature or temperature and concentration in a nonlinear form (see [12–17]):
(2)
ρ
=
ρ
∞
[
1
-
β
T
(
T
-
T
∞
)
+
β
1
(
T
-
T
∞
)
2
]
.
Karcher and Müller [18] used the formulation below to define the nonlinearity of the relationship between the density, the temperature, and the concentration:
(3)
ρ
=
ρ
∞
[
1
-
β
0
(
T
-
T
∞
)
-
β
1
(
T
-
T
∞
)
2
-
β
2
(
C
-
C
∞
)
]
,
where ρ
∞ is the constant fluid density, T
∞ and C
∞ are the fluid temperature and solutal concentration, respectively, β
0 and β
2 are the coefficients of thermal and solutal expansion, and β
1 denotes the nonlinear coefficient of thermal expansion. Partha [16] investigated the natural nonlinear convection in a non-Darcy porous medium using a temperature-concentration-dependent density relation in the form
(4)
ρ
=
ρ
∞
[
β
0
(
T
-
T
∞
)
+
β
1
(
T
-
T
∞
)
2
+
β
2
(
C
-
C
∞
)
+
β
3
(
C
-
C
∞
)
2
]
.
With the usual Boussinesq and the boundary layer approximations, the governing equations are in the form
(5)
∂
u
∂
x
+
∂
v
∂
y
=
0
,
(6)
∂
u
∂
t
+
u
∂
u
∂
x
+
v
∂
u
∂
y
=
ν
∂
2
u
∂
y
2
+
[
β
0
(
T
-
T
∞
)
+
β
1
(
T
-
T
∞
)
2
]
+
[
β
2
(
C
-
C
∞
)
+
β
3
(
C
-
C
∞
)
2
]
g
,
(7)
u
∂
T
∂
x
+
v
∂
T
∂
y
=
α
m
∂
2
T
∂
y
2
,
(8)
u
∂
C
∂
x
+
v
∂
C
∂
y
=
D
∂
2
C
∂
y
2
,
subject to the boundary conditions
(9)
u
=
a
x
,
v
=
0
,
T
=
T
w
,
C
=
C
w
,
ssssssssssssssssssssss
at
y
=
0
,
t
≥
0
,
u
⟶
0
,
T
⟶
T
∞
,
C
⟶
C
∞
ssssssssssssssssss
as
y
⟶
∞
,
t
≥
0
.
The initial conditions are
(10)
t
<
0
:
u
=
0
,
v
=
0
,
T
=
T
w
,
C
=
C
w
,
sssssssssssssssssssssssssssssssssssss
∀
x
,
y
,
where u and v are the velocity components along x and y directions, respectively, T and C are the local fluid temperature and solute concentration across the boundary layer, respectively, ν is the kinematic viscosity, g is the acceleration due to gravity, β
0 and β
1 are the thermal expansion coefficients, β
2 and β
3 are solutal expansion coefficients, α
m is the effective thermal diffusivity, and D is the mass diffusivity.We introduce the following nondimensional variables:(11)
η
=
a
ν
ξ
,
ξ
=
1
-
exp
⁡
⁡
(
-
τ
)
,
τ
=
a
t
,
ψ
=
a
ν
ξ
x
f
(
ξ
,
η
)
,
θ
(
ξ
,
η
)
=
T
-
T
∞
T
w
-
T
∞
,
ϕ
(
ξ
,
η
)
=
C
-
C
∞
C
w
-
C
∞
.
The governing equations (6)–(8) along with the boundary conditions (9) can be presented in the form
(12)
f
′′′
+
(
1
-
ξ
)
[
η
2
f
′′
-
ξ
∂
f
′
∂
ξ
]
+
ξ
[
f
f
′′
-
f
′
2
]
+
λ
ξ
[
(
1
+
α
θ
)
θ
+
Nr
(
1
+
σ
ϕ
)
ϕ
]
=
0
,
θ
′′
+
Pr
⁡
(
1
-
ξ
)
[
η
2
θ
′
-
ξ
∂
θ
∂
ξ
]
+
Pr
ξ
f
θ
′
=
0
,
1
Le
ϕ
′′
+
Sc
(
1
-
ξ
)
[
η
2
ϕ
′
-
ξ
∂
ϕ
∂
ξ
]
+
Sc
ξ
f
ϕ
′
=
0
,
subject to the boundary conditions
(13)
f
(
ξ
,
0
)
=
0
,
f
′
(
ξ
,
0
)
=
1
θ
(
ξ
,
0
)
=
1
,
ϕ
(
ξ
,
0
)
=
1
,
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
ξ
≥
0
,
f
′
(
ξ
,
∞
)
=
0
,
θ
(
ξ
,
∞
)
=
0
,
ϕ
(
ξ
,
∞
)
=
0
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
ξ
≥
0
,
where the Rayleigh number R
a
x, the nonlinear temperature parameter γ, the buoyancy parameter Nr, the Prandtl number Pr, the Schmidt number Sc, the Peclet number P
e
x, and the Lewis number Le are defined as
(14)
Ra
x
=
(
1
-
ϕ
^
f
∞
)
g
β
0
ρ
f
∞
Δ
T
x
α
m
a
,
γ
=
β
1
(
T
w
-
T
∞
)
β
0
,
σ
=
β
3
(
ϕ
^
w
-
ϕ
^
∞
)
β
2
,
Nr
=
(
ρ
p
-
ρ
∞
)
(
ϕ
^
w
-
ϕ
^
∞
)
β
0
ρ
f
∞
(
1
-
ϕ
^
∞
)
(
T
w
-
T
∞
)
,
Pr
=
ν
α
m
,
Sc
=
ν
D
B
,
Le
=
α
m
D
B
,
Pe
x
=
u
w
x
α
m
,
λ
=
Ra
x
Pe
x
.
The skin friction and heat and mass transfer coefficients are described by the local skin friction coefficient, Nusselt number Nu, and the Sherwood number Sh defined as
(15)
τ
w
=
μ
∂
u
∂
y
|
y
=
0
,
q
w
=
-
α
m
∂
T
∂
y
|
y
=
0
,
q
m
=
-
D
∂
C
∂
y
|
y
=
0
.
The nondimensional form of the skin friction coefficient, the reduced Nusselt number, and the reduced Sherwood number are
(16)
ξ
1
/
2
Re
x
1
/
2
C
f
=
f
′′
(
ξ
,
0
)
,
Nu
R
=
ξ
1
/
2
Re
x
-
1
/
2
Nu
=
-
θ
′
(
ξ
,
0
)
,
ξ
1
/
2
Re
x
-
1
/
2
Sh
=
-
ϕ
′
(
ξ
,
0
)
.


## 3. Method of Solution

The system of (12) was solved using the spectral homotopy analysis method. The spectral quasilinearisation method (SQLM) was used to validate the results. The spectral homotopy analysis method (SHAM) was first introduced for the solution of nonlinear ordinary differential equation. The method is a hybrid numerical scheme that combines the underlying ideas of the homotopy analysis method (HAM) and the Chebyshev spectral collocation method. The HAM scheme breaks down a nonlinear differential equation into infinitely many linear ordinary differential equations whose solutions are found analytically. The SHAM provides flexibility when solving linear ODEs through the use of the Chebyshev spectral collocation method. The first attempt to extend the application of the SHAM to PDEs was made by Motsa [3] who solved a nonlinear partial differential equation for an unsteady boundary-layer flow caused by an impulsively stretching plate. The blending of homotopy based methods with other methods has also been considered in [19, 20] wherein the homotopy perturbation method was used in conjunction with Laplace transform to solve partial differential equations. In this work, we extend the application to a system of nonlinear PDEs. In the framework of the SHAM, the nonlinear equations are decomposed into their linear and nonlinear parts as follows:
(17)
ℒ
f
[
f
(
η
,
ξ
)
,
θ
(
η
,
ξ
)
,
ϕ
(
η
,
ξ
)
]
=
f
′′′
+
η
2
(
1
-
ξ
)
f
′′
-
ξ
(
1
-
ξ
)
∂
f
′
∂
ξ
+
λ
ξ
(
θ
+
Nr
ϕ
)
,
𝒩
f
[
f
(
η
,
ξ
)
,
θ
(
η
,
ξ
)
,
ϕ
(
η
,
ξ
)
]
=
f
′′′
+
η
2
(
1
-
ξ
)
f
′′
-
ξ
(
1
-
ξ
)
∂
f
′
∂
ξ
+
λ
ξ
(
θ
+
Nr
ϕ
)
+
ξ
(
f
f
′′
-
f
′
2
)
+
λ
ξ
(
α
θ
2
+
Nr
σ
ϕ
2
)
,
ℒ
θ
[
θ
(
η
,
ξ
)
]
=
θ
′′
+
1
2
Pr
(
1
-
ξ
)
η
θ
′
-
Pr
ξ
(
1
-
ξ
)
∂
θ
∂
ξ
,
𝒩
θ
[
f
(
η
,
ξ
)
,
θ
(
η
,
ξ
)
]
=
θ
′′
+
1
2
Pr
(
1
-
ξ
)
η
θ
′
-
Pr
ξ
(
1
-
ξ
)
∂
θ
∂
ξ
+
Pr
ξ
f
θ
′
,
ℒ
ϕ
[
ϕ
(
η
,
ξ
)
]
=
ϕ
′
′
+
Le
[
1
2
Sc
(
1
-
ξ
)
η
ϕ
′
-
Sc
ξ
(
1
-
ξ
)
∂
ϕ
∂
ξ
]
,
𝒩
ϕ
[
f
(
η
,
ξ
)
,
ϕ
(
η
,
ξ
)
]
=
ϕ
′
′
+
Le
[
1
2
Sc
(
1
-
ξ
)
η
ϕ
′
-
Sc
ξ
(
1
-
ξ
)
∂
ϕ
∂
ξ
]
+
Le
Sc
ξ
f
ϕ
′
,
with ℒ
j and 𝒩
j (j
=
f
,
θ
,
ϕ) representing the linear and nonlinear operators of the equations, respectively. According to the HAM description, thezeroth order deformation equations are formulated as
(18)
(
1
-
q
)
ℒ
f
[
F
(
η
,
ξ
;
q
)
-
f
0
(
η
,
ξ
)
]
=
q
Λ
𝒩
f
[
F
(
η
,
ξ
;
q
)
]
,
(
1
-
q
)
ℒ
θ
[
Θ
(
η
,
ξ
;
q
)
-
θ
0
(
η
,
ξ
)
]
=
q
Λ
𝒩
θ
[
Θ
(
η
,
ξ
;
q
)
]
,
(
1
-
q
)
ℒ
ϕ
[
Φ
(
η
,
ξ
;
q
)
-
ϕ
0
(
η
,
ξ
)
]
=
q
Λ
𝒩
ϕ
[
Φ
(
η
,
ξ
;
q
)
]
,
subject to the boundary conditions
(19)
F
(
0
,
ξ
;
q
)
=
0
,
∂
F
(
η
,
ξ
;
q
)
∂
η
|
η
=
0
=
1
,
∂
F
(
η
,
ξ
;
q
)
∂
η
|
η
=
∞
=
0
,
Θ
(
0
,
ξ
;
q
)
=
1
,
Θ
(
∞
,
ξ
;
q
)
=
0
,
Φ
(
0
,
ξ
;
q
)
=
1
,
Φ
(
∞
,
ξ
;
q
)
=
0
.
The functions F, Θ, and Φ are unknown, q
∈
[
0,1
] is the embedding parameter, Λ is the nonzero convergence controlling auxiliary parameter, and the functions f
0
(
η
,
ξ
),  θ
0
(
η
,
ξ
), and ϕ
0
(
η
,
ξ
) are the initial solutions. From thezeroth order deformation equations above, it can be noted that, as q increases from 0 to 1, the functions F,  Θ, and Φ vary from the initial approximations f
0
(
η
,
ξ
),  θ
0
(
η
,
ξ
), and ϕ
0
(
η
,
ξ
), respectively, to the solutions f
(
η
,
ξ
),  θ
(
η
,
ξ
), and ϕ
(
η
,
ξ
) of the original equations (12). Using the Taylor series to expand F,  Θ, and Φ about q yields
(20)
F
(
η
,
ξ
;
q
)
=
F
(
η
,
ξ
;
0
)
+
∑
n
=
1
+
∞
‍
f
n
(
η
,
ξ
)
q
n
,
Θ
(
η
,
ξ
;
q
)
=
Θ
(
η
,
ξ
;
0
)
+
∑
n
=
1
+
∞
‍
θ
n
(
η
,
ξ
)
q
n
,
Φ
(
η
,
ξ
;
q
)
=
Φ
(
η
,
ξ
;
0
)
+
∑
n
=
1
+
∞
‍
ϕ
n
(
η
,
ξ
)
q
n
,
with
(21)
f
(
η
,
ξ
)
=
1
n
!
∂
n
F
(
η
,
ξ
;
q
)
∂
q
n
|
q
=
0
,
θ
(
η
,
ξ
)
=
1
n
!
∂
n
Θ
(
η
,
ξ
;
q
)
∂
q
n
|
q
=
0
,
ϕ
(
η
,
ξ
)
=
1
n
!
∂
n
Φ
(
η
,
ξ
;
q
)
∂
q
n
|
q
=
0
.
Therefore, since
(22)
F
(
η
,
ξ
;
1
)
=
f
(
η
,
ξ
)
,
Θ
(
η
,
ξ
;
1
)
=
θ
(
η
,
ξ
)
,
Φ
(
η
,
ξ
;
1
)
=
ϕ
(
η
,
ξ
)
,
F
(
η
,
ξ
;
0
)
=
f
0
(
η
,
ξ
)
,
Θ
(
η
,
ξ
;
0
)
=
θ
0
(
η
,
ξ
)
,
Φ
(
η
,
ξ
;
0
)
=
ϕ
0
(
η
,
ξ
)
,
we obtain
(23)
f
(
η
,
ξ
)
=
f
0
(
η
,
ξ
)
+
∑
m
=
1
+
∞
‍
f
m
(
η
,
ξ
)
,
θ
(
η
,
ξ
)
=
θ
0
(
η
,
ξ
)
+
∑
m
=
1
+
∞
‍
θ
m
(
η
,
ξ
)
,
ϕ
(
η
,
ξ
)
=
ϕ
0
(
η
,
ξ
)
+
∑
m
=
1
+
∞
‍
ϕ
m
(
η
,
ξ
)
.
The series (23) converge provided the auxiliary parameter Λ and the linear operators are properly chosen. The linear operators are often chosen to be the linear parts of the governing equations. The functions f
m,  θ
m, and ϕ
m are found from solutions of higher order deformation equations which are obtained by differentiating thezeroth order deformation equations (18) m times with respect to q, dividing by m
! and setting q
=
0. This gives rise to the following equations:
(24)
ℒ
f
[
f
m
(
η
,
ξ
)
-
χ
m
f
m
-
1
(
η
,
ξ
)
]
=
Λ
R
m
f
,
ℒ
θ
[
θ
m
(
η
,
ξ
)
-
χ
m
θ
m
-
1
(
η
,
ξ
)
]
=
Λ
R
m
θ
,
ℒ
ϕ
[
ϕ
m
(
η
,
ξ
)
-
χ
m
ϕ
m
-
1
(
η
,
ξ
)
]
=
Λ
R
m
ϕ
,
with boundary conditions
(25)
f
m
(
0
,
ξ
)
=
0
,
f
m
′
(
0
,
ξ
)
=
0
,
f
m
′
(
+
∞
,
ξ
)
=
0
,
θ
m
(
0
,
ξ
)
=
0
,
θ
m
(
+
∞
,
ξ
)
=
0
,
ϕ
m
(
0
,
ξ
)
=
0
,
ϕ
m
(
+
∞
,
ξ
)
=
0
.
Here the prime denotes differentiation with respect to η and
(26)
R
m
f
=
f
m
-
1
′′′
+
1
2
(
1
-
ξ
)
η
f
m
-
1
′′
-
ξ
(
1
-
ξ
)
∂
f
m
-
1
′
∂
ξ
+
λ
ξ
(
θ
m
-
1
+
Nr
ϕ
m
-
1
)
+
ξ
∑
n
=
0
m
-
1
‍
[
f
m
-
1
-
n
f
n
′′
-
f
′
m
-
1
-
n
f
n
′
+
λ
(
γ
θ
m
-
1
-
n
θ
n
+
Nr
σ
ϕ
m
-
1
-
n
ϕ
n
)
f
n
′′
]
,
R
m
θ
=
θ
m
-
1
′′
+
1
2
Pr
(
1
-
ξ
)
η
ϕ
m
-
1
′
-
Pr
ξ
(
1
-
ξ
)
∂
θ
m
-
1
∂
ξ
+
Pr
ξ
∑
n
=
0
m
-
1
‍
f
m
-
1
-
n
ϕ
n
′
,
R
m
ϕ
=
1
Le
ϕ
m
-
1
′
′
+
1
2
Sc
(
1
-
ξ
)
η
ϕ
m
-
1
′
-
Sc
ξ
(
1
-
ξ
)
∂
ϕ
m
-
1
∂
ξ
+
Sc
ξ
∑
n
=
0
m
-
1
‍
f
m
-
1
-
n
ϕ
n
′
,
where
(27)
χ
m
=
{
0
,
m
≤
1
,
1
,
m
>
1
.
At this stage, in the HAM framework, the higher order differential equations (24) are solved analytically. However, here, it is not possible since the linear operators contain variable coefficients. The rule of solution expression is modified in the ξ space direction and we assume solutions of the form
(28)
f
(
η
,
ξ
)
=
∑
j
=
0
+
∞
∑
k
=
0
+
∞
‍
ξ
k
f
j
,
k
(
η
)
with
f
m
(
η
,
ξ
)

=

∑
k
=
0
+
∞
‍
ξ
k
f
m
,
k
(
η
,
ξ
)
,
θ
(
η
,
ξ
)
=
∑
j
=
0
+
∞
∑
k
=
0
+
∞
‍
ξ
k
θ
j
,
k
(
η
)
with
θ
m
(
η
,
ξ
)

=

∑
k
=
0
+
∞
‍
ξ
k
θ
m
,
k
(
η
,
ξ
)
,
ϕ
(
η
,
ξ
)
=
∑
j
=
0
+
∞
∑
k
=
0
+
∞
‍
ξ
k
ϕ
j
,
k
(
η
)
with
ϕ
m
(
η
,
ξ
)

=

∑
k
=
0
+
∞
‍
ξ
k
ϕ
m
,
k
(
η
,
ξ
)
.
Initial approximate solutions are chosen to satisfy the rule of solution expression in (28), the corresponding boundary conditions, and the equations
(29)
ℒ
f
[
f
0
(
η
,
ξ
)
]
=
0
,
f
0
(
0
,
ξ
)
=
0
,
f
0
′
(
0
,
ξ
)
=
1
,
f
0
′
(
∞
,
ξ
)
=
0
,
ℒ
θ
[
θ
0
(
η
,
ξ
)
]
=
0
,
θ
0
(
0
,
ξ
)
=
1
,
θ
0
(
∞
,
ξ
)
=
0
,
ℒ
ϕ
[
ϕ
0
(
η
,
ξ
)
]
=
0
,
ϕ
0
(
0
,
ξ
)
=
1
,
ϕ
0
(
∞
,
ξ
)
=
0
.
Substituting (28) in (29) and balancing terms of equal order in ξ give
(30)
f
0,0
′′′
+
η
2
f
0,0
′′
=
0
,
f
0,0
(
0
)
=
0
,
f
0,0
′
(
0
)
=
1
,
f
0,0
′
(
∞
)
=
0
,
θ
0,0
′′
+
Pr
η
2
θ
0,0
′
=
0
,
θ
0,0
(
0
)
=
1
,
θ
0,0
(
∞
)
=
0
,
1
Le
ϕ
0,0
′
′
+
Sc
η
2
ϕ
0,0
′
=
0
,
ϕ
0,0
(
0
)
=
1
,
ϕ
0,0
(
∞
)
=
0
,
whose solutions are, respectively,
(31)
f
0,0
(
η
)
=
η
erfc
⁡
(
η
2
)
+
2
π
[
1
-
exp
⁡
(
-
η
2
4
)
]
,
θ
0,0
(
η
)
=
erfc
⁡
(
Pr
y
2
)
,
ϕ
0,0
(
η
)
=
erfc
⁡
(
Le
Sc
y
2
)
.
The function erfc(
η
) is the standard complementary error function defined by
(32)
erfc
⁡
(
η
)
=
2
π
∫
η
∞
‍
exp
⁡
(
-
z
2
)
d
z
.
Substituting (28) in the higher order deformation equations (24) gives rise to the iteration scheme
(33)
f
m
,
k
′′′
+
η
2
f
m
,
k
′′
-
k
f
m
,
k
′
=
(
Λ
+
χ
m
)
{
f
m
-
1
,
k
′′′
+
η
2
f
m
-
1
,
k
′′
-
k
f
m
-
1
,
k
′
}
+
(
Λ
+
χ
m
)
{
(
k
-
1
)
f
m
-
1
,
k
-
1
′
-
η
2
f
m
-
1
,
k
-
1
′′
+
λ
θ
m
-
1
,
k
-
1
+
λ
Nr
ϕ
m
-
1
,
k
-
1
-
η
2
f
m
-
1
,
k
-
1
′′
}
+
η
2
f
m
,
k
-
1
′′
-
(
k
-
1
)
f
m
,
k
-
1
′
-
λ
θ
m
,
k
-
1
-
λ
Nr
ϕ
m
,
k
-
1
×
Λ
∑
n
=
0
m
-
1
‍
∑
p
=
0
k
-
1
‍
[
f
n
,
p
f
m
-
1
-
n
,
k
-
1
-
p
′′
-
f
n
,
p
′
f
m
-
1
-
n
,
k
-
1
-
p
′
+
λ
γ
θ
n
,
p
θ
m
-
1
-
n
,
k
-
1
-
p
+
λ
σ
ϕ
n
,
p
ϕ
m
-
1
-
n
,
k
-
1
-
p
f
m
-
1
-
n
,
k
-
1
-
p
′′
]
,
θ
m
,
k
′′
+
η
2
Pr
⁡
θ
m
,
k
′
-
k
θ
m
,
k
=
(
Λ
+
χ
m
)
{
θ
m
-
1
,
k
′′
+
η
2
Pr
⁡
θ
m
-
1
,
k
′
-
k
Pr
θ
m
-
1
,
k
}
×
(
Λ
+
χ
m
)
{
Pr
⁡
(
k
-
1
)
θ
m
-
1
,
k
-
1
-
Pr
η
2
θ
m
-
1
,
k
-
1
′
}
+
η
2
Pr
⁡
θ
′
m
,
k
-
1
-
Pr
⁡
(
k
-
1
)
θ
m
,
k
-
1
+
Λ
Pr
∑
n
=
0
m
-
1
‍
∑
p
=
0
k
-
1
‍
[
f
n
,
p
θ
m
-
1
-
n
,
k
-
1
-
p
′
]
,
ϕ
m
,
k
′
′
+
Le
[
η
2
Sc
ϕ
m
,
k
′
-
k
Sc
ϕ
m
,
k
]
=
(
Λ
+
χ
m
)
{
ϕ
m
-
1
,
k
′
′
+
Le
Sc
η
2
ϕ
m
-
1
,
k
′
-
k
Le
Sc
ϕ
m
-
1
,
k
}
×
(
Λ
+
χ
m
)
{
Le
Sc
(
k
-
1
)
ϕ
m
-
1
,
k
-
1
-
Le
Sc
η
2
ϕ
m
-
1
,
k
-
1
′
}
+
η
2
Le
Sc
ϕ
m
,
k
-
1
′
-
Le
Sc
(
k
-
1
)
ϕ
m
,
k
-
1
+
Λ
Le
Sc
∑
n
=
0
m
-
1
‍
∑
p
=
0
k
-
1
‍
[
f
n
,
p
ϕ
m
-
1
-
n
,
k
-
1
-
p
′
]
,
with boundary conditions
(34)
f
m
,
k
(
0
)
=
0
,
f
m
,
k
′
(
0
)
=
0
,
f
m
,
k
′
(
∞
)
=
0
,
θ
m
,
k
(
0
)
=
0
,
θ
m
,
k
(
∞
)
=
0
,
ϕ
m
,
k
(
0
)
=
0
,
ϕ
m
,
k
(
∞
)
=
0
.
The system (33) is a linear sequence of ordinary differential equations that can be solved by any numerical method. In this study the equations are solved using the spectral methods whose underlying ideas will be sketched in Section 3.1.


### 3.1. The Spectral Quasilinearisation Method (SQLM)

In this section we present a spectral quasilinearisation method (SQLM) for solving systems of partial differential equations such as those in (12). Consider first an initial unsteady state solution corresponding to τ
=
0 and obtained by setting ξ
=
0 in (12). This gives the equations
(35)
f
′′′
+
1
2
η
f
′′
=
0
,
f
(
0,0
)
=
0
,
f
′
(
0,0
)
=
1
,
f
′
(
∞
,
0
)
=
0
,
θ
′′
+
Pr
η
2
θ
′
=
0
,
θ
(
0,0
)
=
1
,
θ
(
∞
,
0
)
=
0
,
1
Le
ϕ
′
′
+
Sc
η
2
ϕ
′
=
0
,
ϕ
(
0,0
)
=
1
,
ϕ
(
∞
,
0
)
=
0
.
Solving (35) analytically gives
(36)
f
(
η
,
0
)
=
η
erfc
⁡
(
η
2
)
+
2
π
[
1
-
exp
⁡
(
-
η
2
4
)
]
,
θ
(
η
,
0
)
=
erfc
⁡
(
Pr
y
2
)
,
ϕ
(
η
,
0
)
=
erfc
⁡
(
Le
Sc
y
2
)
.
The quasilinearisation method (QLM) is essentially a generalized Newton-Raphson method that was originally used by Bellman and Kalaba [4] for solving functional equations. The QLM scheme is derived by linearising the nonlinear component of the governing equations using Taylor series expansion with the assumption that the difference between the value of the unknown functions at the current iteration level (denoted by r
+
1) and the value at the previous iteration level (denoted by r) is small. Applying the QLM on (12) leads to the following system of linear equations:
(37)
f
r
+
1
′′′
+
a
1
,
r
f
r
+
1
′′
+
a
2
,
r
f
r
+
1
′
+
a
3
,
r
f
r
+
1
+
a
4
,
r
θ
r
+
1
+
a
5
,
r
ϕ
r
+
1
+
a
6
,
r
=
ξ
(
1
-
ξ
)
∂
f
r
+
1
′
∂
ξ
,
θ
r
+
1
′′
+
b
1
,
r
θ
r
+
1
′
+
b
2
,
r
f
r
+
1
+
b
3
,
r
=
Pr
ξ
(
1
-
ξ
)
∂
θ
r
+
1
∂
ξ
,
ϕ
r
+
1
′′
+
c
1
,
r
ϕ
r
+
1
′
+
c
2
,
r
f
r
+
1
+
c
3
,
r
=
Le
Sc
ξ
(
1
-
ξ
)
∂
ϕ
r
+
1
∂
ξ
,
where the coefficients a
p
,
r,  b
p
,
r, and c
p
,
r (p
=
1,2
,
…) are known from previous calculations and are given by
(38)
a
1
,
r
=
(
1
-
ξ
)
η
2
+
ξ
f
r
,
a
2
,
r
=
-
2
ξ
f
r
′
,
a
3
,
r
=
ξ
f
r
′′
,
a
4
,
r
=
λ
ξ
+
2
λ
α
ξ
θ
r
,
a
5
,
r
=
λ
ξ
Nr
+
2
λ
σ
Nr
ξ
ϕ
r
,
a
6
,
r
=
ξ
(
f
r
′
)
2
-
ξ
f
r
f
r
′′
-
λ
ξ
(
α
θ
r
2
+
Nr
σ
ϕ
r
2
)
,
b
1
,
r
=
Pr
(
1
-
ξ
)
η
2
+
Pr
ξ
f
r
,
b
2
,
r
=
Pr
ξ
θ
r
′
,
b
3,2
=
-
Pr
ξ
f
r
θ
r
′
,
c
1
,
r
=
Le
Sc
(
1
-
ξ
)
η
2
+
Le
Sc
ξ
f
r
,
(39)
c
2
,
r
=
Le
Sc
ξ
ϕ
r
′
,
c
3
,
r
=
-
Le
Sc
ξ
f
r
ϕ
r
′
.
Starting from the initial approximations (36), the iteration schemes (37) can be solved iteratively for f
r
+
1
(
η
,
ξ
),  θ
r
+
1
(
η
,
ξ
), and ϕ
r
+
1
(
η
,
ξ
) when r
=
0,1
,
2
,
….In solving (37), we discretize the equation using the Chebyshev spectral collocation method [21] in the η-direction and use an implicit finite difference method in the ξ-direction. The spectral method can be implemented on the interval [
-
1,1
]; therefore it is necessary to transform the domain of the problem to the interval [
-
1,1
] before applying the spectral method. For the convenience of the numerical computations, 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. The transformation η
=
η
∞
(
Y
+
1
)
/
2 is used 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 estimate the derivatives of the unknown variables f
(
η
),  θ
(
η
), and ϕ
(
η
) at the collocation points (grid points) as the matrix vector product. For example, the vector product for f
(
η
) will be as follows:
(40)
d
F
(
1
)
d
η
=
∑
k
=
0
N
x
‍
D
j
k
f
(
Y
k
)
=
D
F
,
j
=
0,1
,
…
,
N
x
,
where N
x
+
1 is the number of collocation points, D
=
2
D
/
η
∞, and
(41)
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,
(42)
F
(
s
)
=
D
s
F
,
where s is the order of the derivative. We choose the Gauss-Lobatto collocation points to define the nodes in [
-
1,1
] as
(43)
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
(44)
Y
j
=
cos
⁡
π
j
N
x
,
ξ
n
=
n
Δ
ξ
,
j
=
0,1
,
…
,
N
x
,
n
=
0,1
,
…
,
N
t
,
where N
x
+
1 and N
t
+
1 are the total numbers 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
(45)
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
Δ
ξ
.
Thus, applying the spectral method to (37) with the finite difference in ξ gives
(46)
[
A
11
A
12
A
13
A
21
A
22
A
23
A
31
A
32
A
33
]
[
F
r
+
1
n
+
1
Θ
r
+
1
n
+
1
Φ
r
+
1
n
+
1
]
=
[
B
11
B
12
B
13
B
21
B
22
B
23
B
31
B
32
B
33
]
[
F
r
+
1
n
Θ
r
+
1
n
Φ
r
+
1
n
]
+
[
K
1
K
2
K
3
]
,
where A
i
j
,
B
i
j (
i
,
j
=
1,2
,
3
) and K
i (i
=
1,2
,
3) are (
N
x
+
1
)
×
(
N
x
+
1
) matrices and (
N
x
+
1
)
×
1 vectors, respectively, defined as
(47)
A
11
=
1
2
(
D
3
+
a
1
,
r
n
+
(
1
/
2
)
D
2
+
a
2
,
r
n
+
(
1
/
2
)
D
+
a
3
,
r
n
+
(
1
/
2
)
)
-
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
D
,
A
12
=
1
2
a
4
,
r
n
+
(
1
/
2
)
,
A
13
=
1
2
a
5
,
r
n
+
(
1
/
2
)
,
A
21
=
1
2
b
2
,
r
n
+
(
1
/
2
)
,
A
22
=
1
2
(
D
2
+
b
1
,
r
n
+
(
1
/
2
)
D
)
-
Pr
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
I
,
A
23
=
0
I
,
A
31
=
1
2
c
2
,
r
n
+
(
1
/
2
)
,
A
32
=
0
I
,
A
33
=
1
2
(
D
2
+
c
1
,
r
n
+
(
1
/
2
)
D
)
-
Le
Sc
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
I
,
K
1
=
-
a
6
,
r
n
+
(
1
/
2
)
,
K
2
=
-
b
3
,
r
n
+
(
1
/
2
)
,
K
3
=
-
c
3
,
r
n
+
(
1
/
2
)
,
B
11
=
-
1
2
(
D
3
+
a
1
,
r
n
+
(
1
/
2
)
D
2
+
a
2
,
r
n
+
(
1
/
2
)
D
+
a
3
,
r
n
+
(
1
/
2
)
)
-
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
D
,
B
12
=
-
1
2
a
4
,
r
n
+
(
1
/
2
)
,
B
13
=
-
1
2
a
5
,
r
n
+
(
1
/
2
)
,
B
21
=
-
1
2
b
2
,
r
n
+
(
1
/
2
)
,
B
22
=
-
1
2
(
D
2
+
b
1
,
r
n
+
(
1
/
2
)
D
)
-
Pr
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
I
,
B
23
=
0
I
,
B
31
=
-
1
2
c
2
,
r
n
+
(
1
/
2
)
,
B
32
=
0
I
,
B
33
=
-
1
2
(
D
2
+
c
1
,
r
n
+
(
1
/
2
)
D
)
-
Le
Sc
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
I
.


## 3.1. The Spectral Quasilinearisation Method (SQLM)

In this section we present a spectral quasilinearisation method (SQLM) for solving systems of partial differential equations such as those in (12). Consider first an initial unsteady state solution corresponding to τ
=
0 and obtained by setting ξ
=
0 in (12). This gives the equations
(35)
f
′′′
+
1
2
η
f
′′
=
0
,
f
(
0,0
)
=
0
,
f
′
(
0,0
)
=
1
,
f
′
(
∞
,
0
)
=
0
,
θ
′′
+
Pr
η
2
θ
′
=
0
,
θ
(
0,0
)
=
1
,
θ
(
∞
,
0
)
=
0
,
1
Le
ϕ
′
′
+
Sc
η
2
ϕ
′
=
0
,
ϕ
(
0,0
)
=
1
,
ϕ
(
∞
,
0
)
=
0
.
Solving (35) analytically gives
(36)
f
(
η
,
0
)
=
η
erfc
⁡
(
η
2
)
+
2
π
[
1
-
exp
⁡
(
-
η
2
4
)
]
,
θ
(
η
,
0
)
=
erfc
⁡
(
Pr
y
2
)
,
ϕ
(
η
,
0
)
=
erfc
⁡
(
Le
Sc
y
2
)
.
The quasilinearisation method (QLM) is essentially a generalized Newton-Raphson method that was originally used by Bellman and Kalaba [4] for solving functional equations. The QLM scheme is derived by linearising the nonlinear component of the governing equations using Taylor series expansion with the assumption that the difference between the value of the unknown functions at the current iteration level (denoted by r
+
1) and the value at the previous iteration level (denoted by r) is small. Applying the QLM on (12) leads to the following system of linear equations:
(37)
f
r
+
1
′′′
+
a
1
,
r
f
r
+
1
′′
+
a
2
,
r
f
r
+
1
′
+
a
3
,
r
f
r
+
1
+
a
4
,
r
θ
r
+
1
+
a
5
,
r
ϕ
r
+
1
+
a
6
,
r
=
ξ
(
1
-
ξ
)
∂
f
r
+
1
′
∂
ξ
,
θ
r
+
1
′′
+
b
1
,
r
θ
r
+
1
′
+
b
2
,
r
f
r
+
1
+
b
3
,
r
=
Pr
ξ
(
1
-
ξ
)
∂
θ
r
+
1
∂
ξ
,
ϕ
r
+
1
′′
+
c
1
,
r
ϕ
r
+
1
′
+
c
2
,
r
f
r
+
1
+
c
3
,
r
=
Le
Sc
ξ
(
1
-
ξ
)
∂
ϕ
r
+
1
∂
ξ
,
where the coefficients a
p
,
r,  b
p
,
r, and c
p
,
r (p
=
1,2
,
…) are known from previous calculations and are given by
(38)
a
1
,
r
=
(
1
-
ξ
)
η
2
+
ξ
f
r
,
a
2
,
r
=
-
2
ξ
f
r
′
,
a
3
,
r
=
ξ
f
r
′′
,
a
4
,
r
=
λ
ξ
+
2
λ
α
ξ
θ
r
,
a
5
,
r
=
λ
ξ
Nr
+
2
λ
σ
Nr
ξ
ϕ
r
,
a
6
,
r
=
ξ
(
f
r
′
)
2
-
ξ
f
r
f
r
′′
-
λ
ξ
(
α
θ
r
2
+
Nr
σ
ϕ
r
2
)
,
b
1
,
r
=
Pr
(
1
-
ξ
)
η
2
+
Pr
ξ
f
r
,
b
2
,
r
=
Pr
ξ
θ
r
′
,
b
3,2
=
-
Pr
ξ
f
r
θ
r
′
,
c
1
,
r
=
Le
Sc
(
1
-
ξ
)
η
2
+
Le
Sc
ξ
f
r
,
(39)
c
2
,
r
=
Le
Sc
ξ
ϕ
r
′
,
c
3
,
r
=
-
Le
Sc
ξ
f
r
ϕ
r
′
.
Starting from the initial approximations (36), the iteration schemes (37) can be solved iteratively for f
r
+
1
(
η
,
ξ
),  θ
r
+
1
(
η
,
ξ
), and ϕ
r
+
1
(
η
,
ξ
) when r
=
0,1
,
2
,
….In solving (37), we discretize the equation using the Chebyshev spectral collocation method [21] in the η-direction and use an implicit finite difference method in the ξ-direction. The spectral method can be implemented on the interval [
-
1,1
]; therefore it is necessary to transform the domain of the problem to the interval [
-
1,1
] before applying the spectral method. For the convenience of the numerical computations, 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. The transformation η
=
η
∞
(
Y
+
1
)
/
2 is used 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 estimate the derivatives of the unknown variables f
(
η
),  θ
(
η
), and ϕ
(
η
) at the collocation points (grid points) as the matrix vector product. For example, the vector product for f
(
η
) will be as follows:
(40)
d
F
(
1
)
d
η
=
∑
k
=
0
N
x
‍
D
j
k
f
(
Y
k
)
=
D
F
,
j
=
0,1
,
…
,
N
x
,
where N
x
+
1 is the number of collocation points, D
=
2
D
/
η
∞, and
(41)
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,
(42)
F
(
s
)
=
D
s
F
,
where s is the order of the derivative. We choose the Gauss-Lobatto collocation points to define the nodes in [
-
1,1
] as
(43)
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
(44)
Y
j
=
cos
⁡
π
j
N
x
,
ξ
n
=
n
Δ
ξ
,
j
=
0,1
,
…
,
N
x
,
n
=
0,1
,
…
,
N
t
,
where N
x
+
1 and N
t
+
1 are the total numbers 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
(45)
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
Δ
ξ
.
Thus, applying the spectral method to (37) with the finite difference in ξ gives
(46)
[
A
11
A
12
A
13
A
21
A
22
A
23
A
31
A
32
A
33
]
[
F
r
+
1
n
+
1
Θ
r
+
1
n
+
1
Φ
r
+
1
n
+
1
]
=
[
B
11
B
12
B
13
B
21
B
22
B
23
B
31
B
32
B
33
]
[
F
r
+
1
n
Θ
r
+
1
n
Φ
r
+
1
n
]
+
[
K
1
K
2
K
3
]
,
where A
i
j
,
B
i
j (
i
,
j
=
1,2
,
3
) and K
i (i
=
1,2
,
3) are (
N
x
+
1
)
×
(
N
x
+
1
) matrices and (
N
x
+
1
)
×
1 vectors, respectively, defined as
(47)
A
11
=
1
2
(
D
3
+
a
1
,
r
n
+
(
1
/
2
)
D
2
+
a
2
,
r
n
+
(
1
/
2
)
D
+
a
3
,
r
n
+
(
1
/
2
)
)
-
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
D
,
A
12
=
1
2
a
4
,
r
n
+
(
1
/
2
)
,
A
13
=
1
2
a
5
,
r
n
+
(
1
/
2
)
,
A
21
=
1
2
b
2
,
r
n
+
(
1
/
2
)
,
A
22
=
1
2
(
D
2
+
b
1
,
r
n
+
(
1
/
2
)
D
)
-
Pr
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
I
,
A
23
=
0
I
,
A
31
=
1
2
c
2
,
r
n
+
(
1
/
2
)
,
A
32
=
0
I
,
A
33
=
1
2
(
D
2
+
c
1
,
r
n
+
(
1
/
2
)
D
)
-
Le
Sc
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
I
,
K
1
=
-
a
6
,
r
n
+
(
1
/
2
)
,
K
2
=
-
b
3
,
r
n
+
(
1
/
2
)
,
K
3
=
-
c
3
,
r
n
+
(
1
/
2
)
,
B
11
=
-
1
2
(
D
3
+
a
1
,
r
n
+
(
1
/
2
)
D
2
+
a
2
,
r
n
+
(
1
/
2
)
D
+
a
3
,
r
n
+
(
1
/
2
)
)
-
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
D
,
B
12
=
-
1
2
a
4
,
r
n
+
(
1
/
2
)
,
B
13
=
-
1
2
a
5
,
r
n
+
(
1
/
2
)
,
B
21
=
-
1
2
b
2
,
r
n
+
(
1
/
2
)
,
B
22
=
-
1
2
(
D
2
+
b
1
,
r
n
+
(
1
/
2
)
D
)
-
Pr
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
I
,
B
23
=
0
I
,
B
31
=
-
1
2
c
2
,
r
n
+
(
1
/
2
)
,
B
32
=
0
I
,
B
33
=
-
1
2
(
D
2
+
c
1
,
r
n
+
(
1
/
2
)
D
)
-
Le
Sc
ξ
n
+
(
1
/
2
)
(
1
-
ξ
n
+
(
1
/
2
)
)
Δ
ξ
I
.


## 4. Results and Discussion

In this section we present numerical solutions to (12) computed using the spectral homotopy analysis method. Approximate solutions of the skin friction coefficient, the surface heat transfer, and the surface mass transfer rates at different values of flow parameters are presented and compared with the SQLM solutions. Velocity, temperature, and concentration profiles at different flow parameters are also presented. All SHAM results were generated using L
=
30, N
=
100, and Λ
=
-
1 unless stated otherwise. These values were found to give accurate solutions after a numerical experimentation. The same L
=
30 and N
=
100 were used for the computation of the SQLM solutions. The M and K in the tables represent the maximum mth and kth iteration required to produce converging results. Values of the skin friction coefficient and reduced Nusselt and Sherwood numbers at different values of ξ are presented in Table 1. The table also shows a comparison of the SHAM and SQLM results. As can be seen from the table, the results match perfectly well for the set accuracy level. The data from the table also reveals a steady increase of f
′′
(
ξ
,
0
), -
θ
′
(
ξ
,
0
), and -
ϕ
′
(
ξ
,
0
) with ξ.Table 1
Comparison of the SHAM solutions forf
′
′
(
ξ
,
0
), -
θ
′
(
ξ
,
0
), and -
ϕ
′
(
ξ
,
0
) against those of the SQLM at different values of ξ.M
=
K
=
30, ℏ = −1, L
=
30, N
=
100, λ
=
0.5, Pr
=
0.7, Nr
=
0.5, Le
=
3, Sc
=
1, α
=
1, and σ
=
1.
ξ
f
′
′
(
ξ
,
0
)
-
θ
′
(
ξ
,
0
)
-
ϕ
′
(
ξ
,
0
)
SHAM
SQLM
SHAM
SQLM
SHAM
SQLM
0.1
−0.540183
−0.540183
0.475790
0.475790
0.993124
0.993124
0.2
−0.516888
−0.516888
0.479891
0.479891
1.010491
1.010491
0.3
−0.494398
−0.494398
0.484385
0.484385
1.029550
1.029550
0.4
−0.472821
−0.472821
0.489330
0.489330
1.050616
1.050616
0.5
−0.452295
−0.452295
0.494804
0.494804
1.074110
1.074110
0.6
−0.432994
−0.432994
0.500905
0.500905
1.100615
1.100615
0.7
−0.415145
−0.415145
0.507775
0.507775
1.130989
1.130989
0.8
−0.399066
−0.399066
0.515629
0.515629
1.166621
1.166621
0.9
−0.385221
−0.385221
0.524843
0.524843
1.210167
1.210167
1.0
−0.374286
−0.374286
0.536307
0.536307
1.268622
1.268622Table2 shows a variation of f
′′
(
ξ
,
0
), -
θ
′
(
ξ
,
0
), and -
ϕ
′
(
ξ
,
0
) with selected flow parameters. It can be seen from the data in Table 2 that all λ
,
α
,
Nr, and σ contribute positively towards the skin friction coefficient, Nusselt number, and Sherwood number as they are seen to increase with the increase in the flow parameters. All SHAM results in the table are accurate for up to six decimal places.Table 2
Variation of the SHAM approximate solutions off
′
′
(
ξ
,
0
), -
θ
′
(
ξ
,
0
), and -
ϕ
′
(
ξ
,
0
) for different flow parameters. M
=
K
=
30, ℏ = −1, L
=
30, N
=
100, Pr
=
0.7,  Le
=
3, and  Sc
=
1.
λ
α
Nr
σ
f
′
′
(
ξ
,
0
)
-
θ
′
(
ξ
,
0
)
-
ϕ
′
(
ξ
,
0
)
0.1
1
0.5
1
−0.721619
0.488512
1.063098
0.5
1
0.5
1
−0.452295
0.494804
1.074110
1.0
1
0.5
1
−0.123169
0.502489
1.087615
0.5
0.1
0.5
1
−0.549445
0.492700
1.070271
0.5
0.5
0.5
1
−0.506182
0.493637
1.071980
0.5
1.0
0.5
1
−0.452295
0.494804
1.074110
0.5
1
0.1
1
−0.516058
0.493757
1.071997
0.5
1
0.5
1
−0.452295
0.494804
1.074110
0.5
1
1.0
1
−0.372924
0.496105
1.076742
0.5
1
0.5
0.1
−0.481198
0.494428
1.073293
0.5
1
0.5
0.5
−0.468347
0.494595
1.073656
0.5
1
0.5
1.0
−0.452295
0.494804
1.074110Figure1 shows the variation of the skin friction coefficient with λ and α. The skin friction coefficient is shown to be an increasing function of ξ for values of λ
=
0.1 and λ
=
0.3 while showing a decreasing function for λ
=
0.5. The skin friction coefficient is shown to decrease with the increase in x
i for any value of α and the decrease is more prominent for increasing values of α.Figure 1
Effect ofλ and α on the skin friction coefficient.The effect ofσ and Nr on the skin friction coefficient is presented in Figure 2. The skin friction coefficient is again seen to be a decreasing function of ξ for any value of σ or Nr. Increasing σ or Nr further decreases the skin friction coefficient.Figure 2
Effect ofσ and Nr on the skin friction coefficient.Figure3 shows the effect of the Prandtl number on the temperature profile. The temperature decreases with an increase in Prandtl number. This is due to the fact that smaller values of Pr are corresponding to larger values of thermal conductivities and therefore heat is able to diffuse away from the stretching sheet.Figure 3
Effect of the Prandtl number on the temperature profile.The effect of the Schmidt numberSc on the concentration profile is shown in Figure 4. It can be seen from the figure that, as the Schmidt number increases, the volume fraction boundary layer decreases across the stretching sheet. A higher Schmidt number implies a lower Brownian diffusion coefficient which will give rise to a shorter penetration depth for the concentration boundary layer.Figure 4
Effect of the Schmidt number on the concentration profile.In Figure5, we present the solution of concentration profile for various values of the Lewis number (Le). Increasing the Lewis number decreases the concentration profile. Larger values of the Lewis number correspond to weaker molecular diffusivity and thinner boundary layer thickness.Figure 5
Effect of the Lewis number on the concentration profile.


## 5. Conclusions

The present study was carried out to extend the application of the spectral homotopy analysis method (SHAM) to systems of nonlinear partial differential equations. The method has been used successfully to find solutions of nonlinear ordinary differential equations. We solve the nonlinear system of partial differential equations governing the unsteady nonlinear convection flow over a stretching flat plate. The spectral quasilinearisation method (SQLM) was used to validate the SHAM results. Numerical values of the skin friction coefficient and the reduced Nusselt and Sherwood numbers were generated using both methods. The results matched well for the desired accuracy. Effects of selected flow parameters on the skin friction coefficient, temperature, and concentration profiles were also presented. The study reveals yet another successful application of the SHAM. The method is a promising tool for solutions of nonlinear partial differential equations and their systems.



---

*Source: 241594-2014-04-29.xml*