# Combined Analytical-Numerical Solution for MHD Viscous Flow over a Stretching Sheet

Authors: Fazle Mabood; Waqar A. Khan
Journal: Journal of Computational Engineering
 (2014)

---

## Abstract

We studied the two-dimensional flow of viscous and electrically conducting fluid over a stretching sheet under the influence of constant magnetic field. Approximate analytical solution of governing nonlinear boundary layer equation via optimal homotopy asymptotic method (OHAM) is obtained. For numerical comparison we usedRunge-Kutta-Fehlberg fourth-fifth-order method. The effect of different parameters on fluid flow is analyzed. It is found that the OHAM solution is very close to the numerical solution for different assigned values of parameters; this thus indicates the feasibility of the proposed method (OHAM).


---

## Body

## 1. Introduction

Magnetohydrodynamics (MHD) is the study of interaction of conducting fluids with electromagnetic phenomena. The flow of an electrically conducting fluid in the presence of magnetic field has many applications in engineering. Also the flow influence by a moving boundary is of crucial importance in the extrusion processes in chemical industries [1, 2].Sakiadis [3, 4] is a pioneer in this area who has investigated the boundary layer flow with uniform speed over continuously stretching surface. Later on, Tsou et al. [5] experimentally verified the work of Sakiadis. Crane [6] studied the steady-state two-dimensional boundary layer flow caused by a stretching sheet whose velocity varies linearly with the distance from a fixed point on the sheet. Chiam [7] and Dandapat and Gupta [8] considered the motion of micropolar and power-law fluids, respectively. Most attention has so far been devoted to the analysis of flow of viscoelastic fluids [9–13] and the joint effect of viscoelasticity and magnetic field has been worked out by Ariel [14]. Beside this, Liao [15, 16] has provided solution for impermeable and permeable stretching sheets, which shows that multiple solutions for the stretching surfaces are possible under definite conditions.The optimal homotopy asymptotic method is a powerful approximate analytical technique that is straightforward to use and does not require the existence of any small or large parameter. Optimal homotopy asymptotic method (OHAM) is employed to construct the series solution of the problem. This method is a consistent analytical tool and it has already been applied to a number of nonlinear differential equations arising in science and engineering [17–19]. So far, as we are aware, there have been no solutions for MHD viscous flow over stretching sheet via OHAM.This paper is organized as follows. First, in Section2, we formulate the problem. In Section 3 we present basic principles of OHAM. The OHAM solution for MHD viscous flow is given in Section 4. In Section 5, we analyze the comparison of the solutions using OHAM with numerical method (NM). Section 6 is devoted to the concluding remarks.


## 2. Governing Equation

We consider the steady two-dimensional flow of an incompressible viscous and electrically conducting fluid over a stretching sheet. A uniform magnetic field is applied and further assumes that the induced magnetic field is negligible due to small magnetic Reynolds number fluid. Under such assumption, the MHD boundary layer flow over a flat plate is governed by(1)
∂
u
∂
x
+
∂
v
∂
y
=
0
,
u
∂
u
∂
x
+
v
∂
u
∂
y
=

ν

∂
2
u
∂
y
2
-
σ
B
2
(
x
)
ρ
u
,
where u and v are the velocity components in the x and y directions, respectively, ν is the kinematic viscosity, ρ is fluid density, σ is the electrical conductivity, and B is the magnetic field strength, where B
(
x
)
=
B
0
x
(
n
-
1
)
/
2.The boundary conditions are given below:(2)
u
(
x
,
0
)
=
c
x
n
,
v
(
x
,
0
)
=
0
,
u
(
x
,
∞
)
=
0
.
We introduce dimensionless variables for nondimensionalized form of momentum and energy equations [13]:
(3)
η
=
c
(
n
+
1
)
2
ν
x
(
n
-
1
)
/
2
y
,
u
=
c
x
n
f
′
(
η
)
,
v
=
-
c
ν
(
n
+
1
)
2
x
(
n
-
1
)
/
2
{
f
(
η
)
+
n
-
1
n
+
1
η
f
′
(
η
)
}
.
Using (3), the governing equations can be reduced to nonlinear differential equation, where f is a function of the similarity variable η:
(4)
f
′′′
+
f
f
′′
-
β
f
′
2
-
M
f
′
=
0
,
subject to the boundary conditions
(5)
f
⟶
0
,
f
′
⟶
1
as
η
⟶
0
,
f
′
⟶
0
as
η
⟶
∞
,
as
(6)
β
=
2
n
n
+
1
,
M
=
2
σ
β
0
2
ρ
c
(
1
+
n
)
,
where β is the nondimensional parameter and M is the magnetic parameter.


## 3. Basic Principles of OHAM

We review the basic principles of OHAM as developed in [17–20] in the following five steps.(i) Let us consider the following differential equation:(7)
A
[
v
(
x
)
]
+
a
(
x
)
=
0
,
x
∈
Ω
,
where Ω is problem domain, A
(
v
)
=
L
(
v
)
+
N
(
v
), where L, N are linear and nonlinear operator, v
(
x
) is an unknown function, and a
(
x
) is a known function.(ii) Construct an optimal homotopy equation as(8)
(
1
-
p
)
[
L
(
ϕ
(
x
;
p
)
)
+
a
(
x
)
]
-
H
(
p
)
[
A
(
ϕ
(
x
;
p
)
)
+
a
(
x
)
]
=
0
,
where 0
≤
p
≤
1 is an embedding parameter, H
(
p
)
=
∑
k
=
1
m
p
k
C
k is auxiliary function on which the convergence of the solution is greatly dependent. The auxiliary function H
(
p
) also adjusts the convergence domain and controls the convergence region.(iii) Expandϕ
(
x
;
p
,
C
j
) in Taylor’s series about p; one has an approximate solution:
(9)
ϕ
(
x
;
p
,
C
j
)
=
v
0
(
x
)
+
∑
k
=
1
∞
v
k
(
x
,
C
j
)
p
k
,
j
=
1,2
,
3
,
…
.
Many researchers [19, 20] have observed that the convergence of the series (9) depends upon C
j, (
j
=
1,2
,
…
,
m
); if it is convergent, then we obtain
(10)
v
~
=
v
0
(
x
)
+
∑
k
=
1
m
v
k
(
x
;
C
j
)
.(iv) Substituting (10) into (8), we have the following residual:
(11)
R
(
x
;
C
j
)
=
L
(
v
~
(
x
;
C
j
)
)
+
a
(
x
)
+
N
(
v
~
(
x
;
C
j
)
)
.
If R
(
x
;
C
j
)
=
0, then v
~ will be the exact solution. For nonlinear problems, generally this will not be the case. For determining C
j, (
j
=
1,2
,
…
,
m
), Galerkin’s Method, Ritz Method, or the method of least squares can be used.(v) Finally, substitute these constants in (11) and one can get the approximate solution.


## 4. Runge-Kutta-Fehlberg Method

The Runge-Kutta-Fehlberg method (RKF45) is a technique to resolve the below problem:(12)
d
y
d
x
=
f
(
x
,
y
)
,
y
(
x
i
)
=
y
i
.
It is called RKF45 because the fourth-order method with five stages is used together with a fifth-order method with six stages that uses all of the points of the first one. It has a procedure to determine if the proper step size h is being used. At each step, two different approximations for the solution are made and compared. If the two answers are in close agreement, the approximation is accepted. If the two answers do not agree to a specified accuracy, the step size is reduced. If the answers agree to more significant digits than required, the step size is increased. Each step requires the use of the following six values.We need the definitions of the following:(13)
K
1
=
h
f
(
x
k
,
y
k
)
,
K
2
=
h
f
(
x
k
+
1
4
h
,
y
k
+
1
4
K
1
)
,
K
3
=
h
f
(
x
k
+
3
8
h
,
y
k
+
3
32
K
1
+
9
32
K
2
)
,
K
4
=
h
f
(
x
k
+
12
13
h
,
y
k
+
1932
2197
K
1
K
4
h
h
h
h
h
h
-
7200
2197
K
2
+
7296
2197
K
3
)
,
K
5
=
h
f
(
x
k
+
h
,
y
k
+
439
216
K
1
-
8
K
2
h
h
K
5
h
h
h
h
+
3680
513
K
3
-
845
4104
K
4
)
,
K
6
=
h
f
(
x
k
+
1
2
h
,
y
k
-
8
27
K
1
+
2
K
2
h
h
h
K
6
h
h
h
-
3544
2565
K
3
+
1859
4104
K
4
-
11
40
K
5
)
.
Then an approximation to the solution of the initial valued problem (I.V.P) is made using a Runge-Kutta method of order 4:
(14)
y
k
+
1
=
y
k
+
25
216
K
1
+
1408
2565
K
3
+
2197
4101
K
4
-
1
5
K
5
,
where the four function values K
1, K
3, K
4
, and K
5 are used.


## 5. Solution of the Problem via OHAM

According to the OHAM, applying (8) to (4),
(15)
(
1
-
p
)
(
f
′′′
)
-
H
(
p
)
(
f
′′′
+
f
f
′′
-
β
f
′
2
-
M
f
′
)
=
0
,
where primes denote differentiation with respect to η.We considerf and H
(
p
) as the following:
(16)
f
=
f
0
+
p
f
1
+
p
2
f
2
+
p
3
f
3
,
H
(
p
)
=
p
C
1
+
p
2
C
2
+
p
3
C
3
.
Using (16) in (15) and some simplifying and rearranging of the terms based on the powers of p, the zeroth-, first-, and second-order problems are as follows.


### 5.1. Zeroth-Order Problem

We have(17)
f
0
′′′
(
η
)
=
0
,
with boundary conditions
(18)
f
0
(
0
)
=
0
,
f
0
′
(
0
)
=
1
,
f
0
′
(
∞
)
=
0
.
Its solution is
(19)
f
0
(
η
)
=
1
10
(
10
η
-
η
2
)
.


### 5.2. First-Order Problem

We have(20)
f
1
(
3
)
(
η
,
C
1
)
=
-
(
M
C
1
f
0
′
+
β
C
1
(
f
0
′
)
2
h
h
h
h
h
(
f
0
′
)
2
-
C
1
f
0
f
0
′′
-
f
0
′′′
-
C
1
f
0
′′′
)
,
with boundary conditions
(21)
f
1
(
0
)
=
0
,
f
1
′
(
0
)
=
0
,
f
1
′
(
∞
)
=
0
.
Its solution is
(22)
f
1
(
η
,
C
1
)
=
1
6000
×
(
1875
η
2
C
1
+
5000
M
η
2
C
1
+
3750
β
η
2
C
1
hhhl
-
1000
M
η
3
C
1
-
1000
β
η
3
C
1
-
50
η
4
C
1
hhhl
+
50
M
η
4
C
1
+
100
β
η
4
C
1
+
2
η
5
C
1
hhhl
-
4
β
η
5
C
1
)
.


### 5.3. Second-Order Problem

We have(23)
f
2
(
3
)
(
η
,
C
1
,
C
2
)
=
-
(
M
C
2
f
0
′
+
β
C
2
(
f
0
′
)
2
+
M
C
1
f
1
′
hhhh
hl
+
2
β
C
1
f
0
′
f
1
′
-
C
2
f
0
f
0
′′
-
C
1
f
1
f
0
′′
hhh
h
h
l
(
f
0
′
)
2
-
C
1
f
0
f
1
′′
-
C
2
f
0
′′′
-
f
1
′′′
-
C
1
f
1
′′′
)
,
with boundary conditions
(24)
f
2
(
0
)
=
0
,
f
2
′
(
0
)
=
0
,
f
2
′
(
∞
)
=
0
.
Its solution becomes
(25)
f
2
(
η
,
C
1
,
C
2
)
=
1
504000
×
(
1575000
η
2
C
1
+
4200000
M
η
2
C
1
hh
hlh
+
3150000
β
η
2
C
1
-
840000
β
η
3
C
1
hh
hlh
-
42000
η
4
C
1
+
⋯
+
1680
η
5
c
2
hhh
lh
-
3360
β
η
5
c
2
)
.


### 5.4. Third-Order Problem

We have(26)
f
3
(
3
)
(
η
,
C
1
,
C
2
,
C
3
)
=
-
(
M
C
3
f
0
′
+
β
C
3
(
f
0
′
)
2
+
M
C
2
f
1
′
+
2
β
C
2
f
0
′
f
1
′
hh
hhh
+
β
C
1
(
f
1
′
)
2
+
M
C
1
f
2
′
+
2
β
C
1
f
0
′
f
2
′
-
C
3
f
0
f
0
′′
hhhh
h
-
C
2
f
1
f
0
′′
-
C
1
f
2
f
0
′′
-
C
2
f
0
f
1
′′
-
C
1
f
1
f
1
′′
hhhh
h
(
f
0
′
)
2
-
C
1
f
0
f
2
′′
-
C
3
f
0
′′′
-
C
2
f
1
′′′
-
f
2
′′′
-
C
1
f
2
′′′
)
,
with boundary conditions
(27)
f
3
(
0
)
=
0
,
f
3
′
(
0
)
=
0
,
f
3
′
(
∞
)
=
0
.
Its solution is a long expression; the few terms are
(28)
f
3
(
η
,
C
1
,
C
2
,
C
3
)
=
1
99792000000
×
(
31185000000
η
2
C
1
+
83160000000
M
η
2
C
1
hh
hhhh
+
62370000000
β
η
2
C
1
+
⋯
+
33264000
η
5
C
3
)
.
We obtain the four-term solution using OHAM for p
=
1:
(29)
f
~
(
η
,
C
1
,
C
2
,
C
3
)
=
f
0
(
η
)
+
f
1
(
η
,
C
1
)
+
f
2
(
η
,
C
1
,
C
2
)
+
f
3
(
η
,
C
1
,
C
2
,
C
3
)
.
We use the method of least squares to obtain the unknown convergent constants C
1, C
2, and  C
3 in (29); for particular case, if M
=
5 and β
=
-
1.5, then the values of   C
1, C
2, and C
3 are given as C
1
=
-
0.058765; C
2
=
-
0.051088; C
3
=
-
0.046104.


## 5.1. Zeroth-Order Problem

We have(17)
f
0
′′′
(
η
)
=
0
,
with boundary conditions
(18)
f
0
(
0
)
=
0
,
f
0
′
(
0
)
=
1
,
f
0
′
(
∞
)
=
0
.
Its solution is
(19)
f
0
(
η
)
=
1
10
(
10
η
-
η
2
)
.


## 5.2. First-Order Problem

We have(20)
f
1
(
3
)
(
η
,
C
1
)
=
-
(
M
C
1
f
0
′
+
β
C
1
(
f
0
′
)
2
h
h
h
h
h
(
f
0
′
)
2
-
C
1
f
0
f
0
′′
-
f
0
′′′
-
C
1
f
0
′′′
)
,
with boundary conditions
(21)
f
1
(
0
)
=
0
,
f
1
′
(
0
)
=
0
,
f
1
′
(
∞
)
=
0
.
Its solution is
(22)
f
1
(
η
,
C
1
)
=
1
6000
×
(
1875
η
2
C
1
+
5000
M
η
2
C
1
+
3750
β
η
2
C
1
hhhl
-
1000
M
η
3
C
1
-
1000
β
η
3
C
1
-
50
η
4
C
1
hhhl
+
50
M
η
4
C
1
+
100
β
η
4
C
1
+
2
η
5
C
1
hhhl
-
4
β
η
5
C
1
)
.


## 5.3. Second-Order Problem

We have(23)
f
2
(
3
)
(
η
,
C
1
,
C
2
)
=
-
(
M
C
2
f
0
′
+
β
C
2
(
f
0
′
)
2
+
M
C
1
f
1
′
hhhh
hl
+
2
β
C
1
f
0
′
f
1
′
-
C
2
f
0
f
0
′′
-
C
1
f
1
f
0
′′
hhh
h
h
l
(
f
0
′
)
2
-
C
1
f
0
f
1
′′
-
C
2
f
0
′′′
-
f
1
′′′
-
C
1
f
1
′′′
)
,
with boundary conditions
(24)
f
2
(
0
)
=
0
,
f
2
′
(
0
)
=
0
,
f
2
′
(
∞
)
=
0
.
Its solution becomes
(25)
f
2
(
η
,
C
1
,
C
2
)
=
1
504000
×
(
1575000
η
2
C
1
+
4200000
M
η
2
C
1
hh
hlh
+
3150000
β
η
2
C
1
-
840000
β
η
3
C
1
hh
hlh
-
42000
η
4
C
1
+
⋯
+
1680
η
5
c
2
hhh
lh
-
3360
β
η
5
c
2
)
.


## 5.4. Third-Order Problem

We have(26)
f
3
(
3
)
(
η
,
C
1
,
C
2
,
C
3
)
=
-
(
M
C
3
f
0
′
+
β
C
3
(
f
0
′
)
2
+
M
C
2
f
1
′
+
2
β
C
2
f
0
′
f
1
′
hh
hhh
+
β
C
1
(
f
1
′
)
2
+
M
C
1
f
2
′
+
2
β
C
1
f
0
′
f
2
′
-
C
3
f
0
f
0
′′
hhhh
h
-
C
2
f
1
f
0
′′
-
C
1
f
2
f
0
′′
-
C
2
f
0
f
1
′′
-
C
1
f
1
f
1
′′
hhhh
h
(
f
0
′
)
2
-
C
1
f
0
f
2
′′
-
C
3
f
0
′′′
-
C
2
f
1
′′′
-
f
2
′′′
-
C
1
f
2
′′′
)
,
with boundary conditions
(27)
f
3
(
0
)
=
0
,
f
3
′
(
0
)
=
0
,
f
3
′
(
∞
)
=
0
.
Its solution is a long expression; the few terms are
(28)
f
3
(
η
,
C
1
,
C
2
,
C
3
)
=
1
99792000000
×
(
31185000000
η
2
C
1
+
83160000000
M
η
2
C
1
hh
hhhh
+
62370000000
β
η
2
C
1
+
⋯
+
33264000
η
5
C
3
)
.
We obtain the four-term solution using OHAM for p
=
1:
(29)
f
~
(
η
,
C
1
,
C
2
,
C
3
)
=
f
0
(
η
)
+
f
1
(
η
,
C
1
)
+
f
2
(
η
,
C
1
,
C
2
)
+
f
3
(
η
,
C
1
,
C
2
,
C
3
)
.
We use the method of least squares to obtain the unknown convergent constants C
1, C
2, and  C
3 in (29); for particular case, if M
=
5 and β
=
-
1.5, then the values of   C
1, C
2, and C
3 are given as C
1
=
-
0.058765; C
2
=
-
0.051088; C
3
=
-
0.046104.


## 6. Results and Discussion

Tables1 and 2 show the comparison of OHAM results with numerical RKF45 and RK-4 for different values of parameters. It is noteworthy to mention here that the low error of OHAM is remarkable, while the effectiveness of the proposed method (OHAM) can be concluded from Figures 1(a), 1(b), 1(c), and 1(d), Figures 2(a), 2(b), 2(c), and 2(d), and Figures 3(a) and 3(b), where we have presented the effect of parameters M and β for several assigned values on f
(
η
). The effect of parameters on dimensionless velocity profile is observed in Figure 4, where the boundary layer thickness decreases with the increase in magnetic parameter M keeping fixed the other parameter β.Table 1
Comparison of OHAM results with RKF45 forβ
=
5, M
=
5.
η
f
(
η
)
f
′
(
η
)
OHAM
RKF45
Ab. error
OHAM
RKF45
Ab. error
0.0
0
0
0
1
1
0
0.1
0.086665
0.086668
3.0
×
10
-
6
0.746974
0.747032
5.8
×
10
-
5
0.2
0.151657
0.151663
6.0
×
10
-
6
0.562191
0.562169
2.2
×
10
-
5
0.3
0.200655
0.200658
3.0
×
10
-
6
0.424331
0.424023
3.1
×
10
-
4
0.4
0.237614
0.237598
1.6
×
10
-
5
0.319327
0.319201
1.2
×
10
-
4
0.5
0.265301
0.265277
2.4
×
10
-
5
0.237681
0.237665
1.6
×
10
-
4
0.6
0.285709
0.285688
2.1
×
10
-
5
0.172819
0.172901
8.2
×
10
-
4
0.7
0.300262
0.300252
1.0
×
10
-
5
0.119889
0.120035
1.4
×
10
-
4
0.8
0.309964
0.309969
5.0
×
10
-
5
0.075254
0.075405
1.5
×
10
-
4
0.9
0.315493
0.315512
1.9
×
10
-
5
0.036055
0.036162
1.1
×
10
-
4
1.0
0.317278
0.317303
2.5
×
10
-
5
0
0
0Table 2
Comparison of OHAM results with RK-4 forβ
=
2, M
=
1.
η
f
(
η
)
f
′
(
η
)
OHAM
RK-4
Ab. error
OHAM
RK-4
Ab. error
0.0
0
0
0
1
1
0
0.1
0.091523
0.091539
1.6
×
10
-
5
0.835210
0.835228
1.8
×
10
-
5
0.2
0.167828
0.167837
9.0
×
10
-
6
0.694252
0.694268
1.6
×
10
-
5
0.3
0.231016
0.231020
4.0
×
10
-
6
0.572149
0.572172
2.3
×
10
-
5
0.4
0.282761
0.282770
9.0
×
10
-
6
0.465013
0.465045
3.2
×
10
-
5
0.5
0.324417
0.324424
7.0
×
10
-
6
0.369751
0.369794
4.3
×
10
-
5
0.6
0.357036
0.357042
6.0
×
10
-
6
0.283923
0.283941
1.8
×
10
-
5
0.7
0.381444
0.381459
1.5
×
10
-
5
0.205462
0.205495
3.3
×
10
-
5
0.8
0.398326
0.398334
8.0
×
10
-
6
0.132827
0.132846
1.9
×
10
-
5
0.9
0.408161
0.408178
1.7
×
10
-
5
0.064653
0.064696
4.3
×
10
-
5
1.0
0.411353
0.411388
3.5
×
10
-
5
0
0
0Comparison of the results off
(
η
) using OHAM and numerical method for several values of M at β
=
1.
(a)
(b)
(c)
(d)Comparison of the results off
(
η
) using OHAM and numerical method for several values of M at β
=
1.5.
(a)
(b)
(c)
(d)Comparison of the results off
(
η
) using OHAM and numerical method for different values of β at M
=
50.
(a)
(b)Figure 4
Results off
′
(
η
) using OHAM for different values of M at β
=
-
1.5.


## 7. Conclusion

In this study, we have successfully applied the optimal homotopy asymptotic method for MHD viscous flow over a stretching sheet. Both numerical and approximate analytical results are obtained for the problem. The results obtained via OHAM are presented in tabular and graphical forms. An excellent agreement is observed in OHAM and numerical solution. It was found that OHAM results are matched well with numerical solution. The approximate solution obtained using OHAM is also consistent with the numerical solution for various values ofβ and M. The objective of this work was to show the usefulness of OHAM.The proposed method is simple in applicability, as it does not require discretization like numerical methods. Furthermore, this method provides a convenient way to control the convergence by optimally determining the auxiliary constants. Additionally, this method converges rapidly at lower order of approximations. Therefore, OHAM shows its latent strength and potential for the solution of nonlinear physical problems in real-life applications.



---

*Source: 634328-2014-03-13.xml*