Journal cover Journal topic
Mechanical Sciences An open-access journal for theoretical and applied mechanics
Journal topic
Mech. Sci., 9, 61–70, 2018
https://doi.org/10.5194/ms-9-61-2018
Mech. Sci., 9, 61–70, 2018
https://doi.org/10.5194/ms-9-61-2018

Research article 14 Feb 2018

Research article | 14 Feb 2018

# Heat transfer and MHD flow of non-newtonian Maxwell fluid through a parallel plate channel: analytical and numerical solution

Heat transfer and MHD flow of non-newtonian Maxwell fluid through a parallel plate channel: analytical and numerical solution
Alireza Rahbari1,2, Morteza Abbasi3, Iman Rahimipetroudi4, Bengt Sundén5, Davood Domiri Ganji3, and Mehdi Gholami6 Alireza Rahbari et al.
• 1Department of Mechanical Engineering, Shahid Rajaee Teacher Training University (SRTTU), Tehran, Iran
• 2Research School of Engineering, The Australian National University, Canberra, ACT 2601, Australia
• 3Department of Mechanical Engineering, Sari Branch, Islamic Azad University, Sari, Iran
• 4Young Researchers Club, Sari Branch, Islamic Azad University, Sari, Iran
• 5Department of Energy Sciences, Division of Heat Transfer, Lund University, Lund, Sweden
• 6Department of Mechanical Engineering, Tehran Science and Research Branch, Islamic Azad University, Damavand, Iran

Correspondence: Alireza Rahbari (alireza.rahbari@anu.edu.au, ar.rahbari@gmail.com)

Abstract

Analytical and numerical analyses have been performed to study the problem of magneto-hydrodynamic (MHD) flow and heat transfer of an upper-convected Maxwell fluid in a parallel plate channel. The governing equations of continuity, momentum and energy are reduced to two ordinary differential equation forms by introducing a similarity transformation. The Homotopy Analysis Method (HAM), Homotopy Perturbation Method (HPM) and fourth-order Runge-Kutta numerical method (NUM) are used to solve this problem. Also, velocity and temperature fields have been computed and shown graphically for various values of the physical parameters. The objectives of the present work are to investigate the effect of the Deborah numbers (De), Hartman electric number (Ha), Reynolds number (Rew) and Prandtl number (Pr) on the velocity and temperature fields. As an important outcome, it is observed that increasing the Hartman number leads to a reduction in the velocity values while increasing the Deborah number has negligible impact on the velocity increment.

1 Introduction

Non-Newtonian fluid flow is a fast-growing field of interest due to its various applications in different fields of engineering (Rivlin and Ericksen, 1955). The problem with this type of fluid is that there is not a single constitutive equation due to various rheological parameters appearing in such fluids. Therefore, different models have been developed in the forms of (i) differential type, (ii) rate type and (iii) integral type (Hayat and Awais, 2011; Hayat et al., 2012). The dynamics of materials with the properties of elasticity and viscosity is a fundamental topic in fluid dynamics. This kind of materials is referred to as Maxwell fluid (Mukhopadhyay, 2012; Adegbie et al., 2015). This model categorized as a subclass of rate type fluids can predict the stress relaxation. The Upper Convected Maxwell (UCM) model is the generalization of the Maxwell material for the case of large deformation using the upper convected time derivative.

Several researchers have investigated the application of MHD flow in different industrial applications (Recebli et al., 2013, 2015; Selimli et al., 2015; Hayat et al., 2017). Reviewing the literature indicates the importance of MHD flow through parallel channels. Hayat et al. (2006) was among the first who studied the MHD flow of an upper convected Maxwell (UCM) fluid over a porous stretching sheet using homotopy analysis method (HAM). Raftari and Yildirim (2010) proposed a novel analytical model (homotopy perturbation method (HPM) to analyze the same problem investigated by (Hayat et al., 2006). Switching from the non-rotating flow to rotating one was made by Sajid et al. (2011). They advanced the knowledge in this field by considering the effect of rotation parameter on the overall performance of the UCM. They concluded that the rotation parameter can help in controlling the thickness of boundary layer. Compared to the analytical techniques used by (Hayat et al., 2006; Raftari and Yildirim, 2010), Abel et al. (2012) proposed a numerical scheme, fourth order Runge–Kutta method, to examine the MHD flow for UCM Maxwell over a stretching sheet. These studies were focused on the no-slip boundary condition imposed along the flat plate. Therefore, Abbasi and Rahimipetroudi (2013) aimed at improving the models in this field by considering the slip boundary condition of an UCM Maxwell fluid via HPM. More complexity was added to the available models by Nadeem et al. (2014). They looked at the performance of the MHD boundary layer flow of a Maxwell fluid with nanoparticles. Afify and Elgazery (2016) extended the problem by Nadeem et al. (2014) by adding the chemical reaction into their model. In a more recent research, effects of thermal radiation (Hayat et al., 2011b) and mass transfer (Hayat et al., 2011a) were studied on the MHD flow of the UCM fluid with porosity in the channel walls.

In recent decades, many attempts have been made to develop analytical methods for solving such nonlinear equations. Analytical techniques such as HPM (Abbasi and Rahimipetroudi, 2013; He, 2000; Abbasi et al., 2014) and HAM (Hayat et al., 2006; Abbasi et al., 2014, 2016; Turkyilmazoglu, 2011) have been successfully applied to solve many types of nonlinear problems. The aim of this study is to investigate the effect of physical parameters on a steady flow of an upper-convected Maxwell fluid in a channel in the presence of an external magnetic field. In addition, the convergence of the series solution is also explicitly discussed. The obtained results from the analytical solution have been compared with the numerical data and the comparison reveals the capability, effectiveness and convenience of the Homotopy Analysis Method (HAM) and the Homotopy Perturbation Method (HPM) methods in solving the stated problem. Finally, the effect of dominating parameters on the MHD flow and heat transfer of UCM in a parallel plate is discussed.

2 Problem statement and mathematical formulation

Figure 1 illustrates the schematic diagram of the magneto-hydrodynamic (MHD) flow of an incompressible UCM fluid in a parallel plate channel. As seen, the x-axis is taken along the centerline of the channel, parallel to the channel surfaces, and the y-axis is transverse to this. The inter-plate region contains an incompressible two-dimensional steady, laminar UCM viscoelastic fluid flow. The flow is symmetric in X-axis direction. The steady state flow occurs in the channel with the fluid suction and injection taking place at both plates with velocityVw. The fluid injection and suction take place through the walls where Vw>0 stands for suction and Vw<0 stands for injection. The walls of the channel are at ${y}^{\ast }=H$ and ${y}^{\ast }=-H$ (where 2H is the channel height). Here u and v are the velocities components in the x and y directions, respectively.

Figure 1Schematic diagram and the coordinate system for the considered flow.

A uniform magnetic field, B0, is imposed along the y-axis. The constitutive equation for a Maxwell fluid is (Bég and Makinde, 2010):

$\begin{array}{}\text{(1)}& \mathit{\tau }+{\mathit{\lambda }}_{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}\stackrel{\mathrm{^}}{\mathit{\tau }}={\mathit{\mu }}_{\mathrm{0}}\phantom{\rule{0.125em}{0ex}}\mathit{\gamma }\end{array}$

where τ, λ1, μ0 and γ are the extra stress tensor, relaxation time, low-shear viscosity and the rate-of-strain tensor, respectively. The upper convected time derivative of the stress tensor $\stackrel{\mathrm{^}}{\mathit{\tau }}$is formulated by:

$\begin{array}{}\text{(2)}& \stackrel{\mathrm{^}}{\mathit{\tau }}=\frac{\partial \mathit{\tau }}{\partial t}+v\phantom{\rule{0.125em}{0ex}}\cdot \mathrm{\nabla }\mathit{\tau }-{\left(\mathrm{\nabla }v\right)}^{T}\cdot \mathit{\tau }-\mathit{\tau }\cdot \mathrm{\nabla }v\end{array}$

where t, v, (⋅)T and v denote time, velocity vector, transpose of tensor and fluid velocity gradient tensor, respectively. Governing equations take the following form for the UCM fluid by implementing the shear stress strain tensor from Eqs. (1) and (2) (Wehgal and Ashraf, 2012; Hayat and Wang, 2003):

$\begin{array}{ll}\text{(3)}& & \frac{\partial {u}^{\ast }}{\partial {x}^{\ast }}+\frac{\partial {\mathit{\nu }}^{\ast }}{\partial {y}^{\ast }}=\mathrm{0}& {u}^{\ast }\frac{\partial {u}^{\ast }}{\partial {x}^{\ast }}+{v}^{\ast }\frac{\partial {u}^{\ast }}{\partial {y}^{\ast }}+\mathit{\lambda }\left[{u}^{\ast \mathrm{2}}\frac{{\partial }^{\mathrm{2}}{u}^{\ast }}{\partial {x}^{\ast \mathrm{2}}}+{v}^{\ast \mathrm{2}}\frac{{\partial }^{\mathrm{2}}{u}^{\ast }}{\partial {y}^{\ast \mathrm{2}}}\right\\ & \phantom{\rule{1em}{0ex}}+\mathrm{2}\phantom{\rule{0.125em}{0ex}}{u}^{\ast }\phantom{\rule{0.125em}{0ex}}{v}^{\ast }\frac{{\partial }^{\mathrm{2}}{u}^{\ast }}{\partial {x}^{\ast }\partial {y}^{\ast }}]\\ \text{(4)}& & \phantom{\rule{1em}{0ex}}=\mathit{\upsilon }\frac{{\partial }^{\mathrm{2}}{u}^{\ast }}{\partial {y}^{\ast \mathrm{2}}}\phantom{\rule{0.125em}{0ex}}+J×B,\text{(5)}& & \mathrm{\nabla }\cdot B=\mathrm{0},\text{(6)}& & \mathrm{\nabla }×B={\mathit{\mu }}_{\mathrm{m}}\phantom{\rule{0.125em}{0ex}}J,\text{(7)}& & \mathrm{\nabla }×E=\mathrm{0},\text{(8)}& & J={\mathit{\sigma }}_{\mathrm{e}}\left(E+V×B\right)\end{array}$

where ρ, , ϑ, p, J, B, b, μm, E, σe are fluid density, nabla operator, kinematic viscosity, pressure, current density, total magnetic field (so that $B={B}_{\mathrm{0}}+b$, b is the induced magnetic field), magnetic permeability, electric field, electrical conductivity of the fluid, respectively. The uniform constant magnetic field B is imposed along the y-direction. Because the magnetic Reynolds number is considered to be small (Shereliff, 1965), the induced magnetic field b is neglected. It is assumed that the electric field is negligible due to the polarization of charges (i.e., Hall effect). Under these assumptions, the MHD body force occurring in Eq. (2) takes the following form (Rossow, 1958; Ganji et al., 2014):

$\begin{array}{}\text{(9)}& J×B={\mathit{\sigma }}_{\mathrm{e}}\left[\left(V×B\right)×{B}_{\mathrm{0}}\right]=-\mathit{\sigma }{B}_{\mathrm{0}}^{\mathrm{2}}{u}^{\ast }\end{array}$

Equations (3) and (4) are rewritten for an electrically conducting incompressible fluid in the presence of a uniform magnetic field as below:

$\begin{array}{ll}\text{(10)}& & \frac{\partial {u}^{\ast }}{\partial {x}^{\ast }}+\frac{\partial {\mathit{\nu }}^{\ast }}{\partial {y}^{\ast }}=\mathrm{0}& {u}^{\ast }\frac{\partial {u}^{\ast }}{\partial {x}^{\ast }}+{v}^{\ast }\frac{\partial {u}^{\ast }}{\partial {y}^{\ast }}+\mathit{\lambda }\left[{u}^{\ast \mathrm{2}}\frac{{\partial }^{\mathrm{2}}{u}^{\ast }}{\partial {x}^{\ast \mathrm{2}}}+{v}^{\ast \mathrm{2}}\frac{{\partial }^{\mathrm{2}}{u}^{\ast }}{\partial {y}^{\ast \mathrm{2}}}\right\\ & +\mathrm{2}\phantom{\rule{0.125em}{0ex}}{u}^{\ast }\phantom{\rule{0.125em}{0ex}}{v}^{\ast }\frac{{\partial }^{\mathrm{2}}{u}^{\ast }}{\partial {x}^{\ast }\partial {y}^{\ast }}]\\ \text{(11)}& & \phantom{\rule{1em}{0ex}}=\mathit{\upsilon }\frac{{\partial }^{\mathrm{2}}{u}^{\ast }}{\partial {y}^{\ast \mathrm{2}}}\phantom{\rule{0.125em}{0ex}}-\frac{\mathit{\sigma }{B}_{\mathrm{0}}^{\mathrm{2}}}{\mathit{\rho }}{u}^{\ast },\end{array}$

In order to complete the formulation of the problem, boundary conditions need to be specified. The appropriate no-slip boundary conditions are identical to (Bég and Makinde, 2010) and owing to symmetry take the form:

$\begin{array}{}\text{(12)}& & {y}^{\ast }=\mathrm{0}:\frac{\partial {u}^{\ast }}{\partial {y}^{\ast }}=\mathrm{0},\text{(13)}& & {y}^{\ast }=H:{u}^{\ast }=\mathrm{0},\phantom{\rule{1em}{0ex}}{v}^{\ast }={V}_{\mathrm{w}}.\end{array}$

The following dimensionless variables are introduced:

$\begin{array}{}\text{(14)}& x=\frac{{x}^{\ast }}{H};\phantom{\rule{1em}{0ex}}y=\frac{{y}^{\ast }}{H};\phantom{\rule{1em}{0ex}}{u}^{\ast }=-\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{w}}\phantom{\rule{0.125em}{0ex}}x\phantom{\rule{0.125em}{0ex}}{f}^{\prime }\left(y\right)\end{array}$

Using Eqs. (10) and (14) yields:

$\begin{array}{}\text{(15)}& {v}^{\ast }=\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{w}}f\left(y\right)\end{array}$

Substituting the dimensionless parameters into the momentum equation leads to:

$\begin{array}{ll}& {f}^{\prime \prime \prime }-{M}^{\mathrm{2}}{f}^{\prime }+{\mathit{\text{Re}}}_{\mathrm{w}}\phantom{\rule{0.125em}{0ex}}\left({{f}^{\prime }}^{\mathrm{2}}-f\phantom{\rule{0.125em}{0ex}}{f}^{\prime \prime }\right)\\ \text{(16)}& & \phantom{\rule{1em}{0ex}}+\mathit{\text{De}}\phantom{\rule{0.125em}{0ex}}\left(\mathrm{2}f\phantom{\rule{0.125em}{0ex}}{f}^{\prime }\phantom{\rule{0.125em}{0ex}}{f}^{\prime \prime }-{f}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{f}^{\prime \prime \prime }\right)=\mathrm{0}\end{array}$

The boundary conditions (12) and (13) in the dimensionless form are given by:

$\begin{array}{ll}& y=\mathrm{0}:{f}^{\prime \prime }=\mathrm{0};\phantom{\rule{1em}{0ex}}f=\mathrm{0}\\ \text{(17)}& & y=\mathrm{1}:{f}^{\prime }=\mathrm{0};\phantom{\rule{1em}{0ex}}f=\mathrm{1}\end{array}$

Although the proposed differential Eq. (16) is of the third order, four boundary conditions are given in Eq. (17). In order to satisfy all the boundary conditions, the third order differential Eq. (16) is differentiated as follows:

$\begin{array}{ll}& f{}^{\prime \prime \prime \prime }-{M}^{\mathrm{2}}{f}^{\prime \prime }+{\mathit{\text{Re}}}_{\mathrm{w}}\phantom{\rule{0.125em}{0ex}}\left({f}^{\prime }\phantom{\rule{0.125em}{0ex}}{f}^{\prime \prime }-f\phantom{\rule{0.125em}{0ex}}{f}^{\prime \prime \prime }\right)\\ \text{(18)}& & \phantom{\rule{1em}{0ex}}+\mathit{\text{De}}\phantom{\rule{0.125em}{0ex}}\left(\mathrm{2}\phantom{\rule{0.125em}{0ex}}{{f}^{\prime }}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{f}^{\prime \prime }+\mathrm{2}\phantom{\rule{0.125em}{0ex}}f\phantom{\rule{0.125em}{0ex}}{{f}^{\prime \prime }}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}-{f}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}f{}^{\prime \prime \prime \prime }\right)=\mathrm{0}\end{array}$

where Rew, De and M are the Reynolds number, Deborah number and Hartman number, respectively, which are defined as:

$\begin{array}{}\text{(19)}& {\mathit{\text{Re}}}_{\mathrm{w}}=\frac{{V}_{\mathrm{w}}\phantom{\rule{0.125em}{0ex}}H}{\mathit{\upsilon }},\phantom{\rule{1em}{0ex}}\mathit{\text{De}}=\frac{\mathit{\lambda }\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{w}}^{\mathrm{2}}}{\mathit{\upsilon }},\phantom{\rule{1em}{0ex}}M=\sqrt{\frac{\mathit{\sigma }\phantom{\rule{0.125em}{0ex}}B{}_{\mathrm{0}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}H}{\mathit{\mu }}}\end{array}$

Figure 2The 1 – validity for M=0, De= 0.4, Rew= 5 and Pr= 0.9 given by the 4, 5, 6 and 7th-order solution.

Figure 3The 1 – validity for M=0.1, De= 0.1, Rew= 5 and Pr= 0.9 given by the 4, 5, 6 and 7th-order solution.

Figure 4The 2 – validity for M=0, De= 0.3, Rew= 5 and Pr= 4 given by the 4, 5, 6 and 7th-order solution.

## Heat transfer problem

By neglecting viscous and ohmic dissipation, the governing equation of energy obtained is as follows:

$\begin{array}{}\text{(20)}& {u}^{\ast }\frac{\partial T}{\partial {x}^{\ast }}+{v}^{\ast }\frac{\partial T}{\partial {y}^{\ast }}=\frac{k}{\mathit{\rho }\phantom{\rule{0.125em}{0ex}}{c}_{\mathrm{p}}}\left(\frac{{\partial }^{\mathrm{2}}T}{\partial {x}^{\ast \mathrm{2}}}+\frac{{\partial }^{\mathrm{2}}T}{\partial {y}^{\ast \mathrm{2}}}\right),\end{array}$

where T(x,y) is the temperature at any point and k is the thermal conductivity. Similar to the case of channel flow, for the flow of most fluids under practically significant conditions, the following relation holds:

$\begin{array}{}\text{(21)}& \frac{{\partial }^{\mathrm{2}}T}{\partial {x}^{\ast \mathrm{2}}}\ll \frac{{\partial }^{\mathrm{2}}T}{\partial {y}^{\ast \mathrm{2}}},\end{array}$

The first term on the right hand side of Eq. (20) can be ignored compared to the second and the energy equation can be approximated as:

$\begin{array}{}\text{(22)}& {u}^{\ast }\frac{\partial T}{\partial {x}^{\ast }}+{v}^{\ast }\frac{\partial T}{\partial {y}^{\ast }}=\frac{k}{\mathit{\rho }\phantom{\rule{0.125em}{0ex}}{c}_{\mathrm{p}}}\left(\frac{{\partial }^{\mathrm{2}}T}{\partial {y}^{\ast \mathrm{2}}}\right),\end{array}$

The flow is assumed to be symmetric, so the appropriate boundary conditions for the above equation are

$\begin{array}{}\text{(23)}& & {y}^{\ast }=\mathrm{0}\to \frac{\partial T}{\partial {y}^{\ast }}=\mathrm{0}\text{(24)}& & {y}^{\ast }=H\to T={T}_{\mathrm{w}}\end{array}$

The non-dimensional variables have been defined as:

$\begin{array}{}\text{(25)}& y=\frac{{y}^{\ast }}{H};\phantom{\rule{1em}{0ex}}\mathit{\theta }=\frac{T}{{T}_{\mathrm{w}}}\end{array}$

Inserting the dimensionless parameters described in Eqs. (14) and (25) into the Eqs. (22)–(24) results in:

$\begin{array}{}\text{(26)}& & {\mathit{\theta }}^{\prime \prime }-\mathit{\text{Pr}}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}f\phantom{\rule{0.125em}{0ex}}\mathit{\theta }=\mathrm{0},\text{(27)}& & y=\mathrm{1}\to \mathit{\theta }\left(y\right)=\mathrm{1}\text{(28)}& & y=\mathrm{0}\to {\mathit{\theta }}^{\prime }\left(y\right)=\mathrm{0}\end{array}$

where $\mathit{\text{Pr}}\phantom{\rule{0.125em}{0ex}}=\frac{\mathit{\mu }{c}_{\mathrm{p}}}{k}$ represents the Prandtl number.

Figure 5Dimensionless velocities predicted by HAM, HPM and numerical method (NUM) for different Rew number and De at Pr=1, $\mathrm{\hslash }=-\mathrm{1.5}$.

Figure 6Effects of De Deborah numbers on the velocity contours for M=1, Pr= 1, Rew=3.

Figure 7Effects of M Hartman numbers on the velocity contours for De=0.2, Pr= 1, Rew=3.

Figure 8Dimensionless temperature predicted by analytical and numerical method (NUM) for different De Deborah number at M=1, Rew=3, Pr= 1, $\mathrm{\hslash }=-\mathrm{1.5}$.

Figure 9Dimensionless temperature predicted by analytical and numerical method (NUM) for different M Hartman number at De= 0.2, Rew=2, Pr= 1, $\mathrm{\hslash }=-\mathrm{1}$.

Figure 10Dimensionless temperature predicted by analytical and numerical method (NUM) for different Pr Prandtl number at M=1, Rew=3, De= 0.3, $\mathrm{\hslash }=-\mathrm{1}$.

3 Analytical methods

## 3.1 Implementation of the Homotopy Perturbation Method

In this section, we employ HPM to solve Eq. (2) subject to the boundary conditions Eq. (3). One can construct the Homotopy function of Eq. (2) as described in (He, 2000):

$\begin{array}{ll}& H\left(f,p\right)=\left(\mathrm{1}-P\right)\left[f{}^{\prime \prime \prime \prime }\left(\mathit{\eta }\right)-{g}_{\mathrm{0}}\left(y\right)\right]\phantom{\rule{0.125em}{0ex}}+p\left\{{f}^{\prime \prime \prime \prime }-{M}^{\mathrm{2}}{f}^{\prime \prime }\right\\\ & \phantom{\rule{1em}{0ex}}+{\mathit{\text{Re}}}_{\mathrm{w}}\phantom{\rule{0.125em}{0ex}}\left({f}^{\prime }\phantom{\rule{0.125em}{0ex}}{f}^{\prime \prime }-f\phantom{\rule{0.125em}{0ex}}{f}^{\prime \prime \prime }\right)\\ \text{(29)}& & \mathit{\text{De}}\left(\mathrm{2}\phantom{\rule{0.125em}{0ex}}{{f}^{\prime }}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{f}^{\prime \prime }+\mathrm{2}\phantom{\rule{0.125em}{0ex}}f\phantom{\rule{0.125em}{0ex}}{{f}^{\prime \prime }}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}-{f}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{f}^{\prime \prime \prime \prime }\right)}\phantom{\rule{0.125em}{0ex}}=\mathrm{0},& H\left(\mathit{\theta },p\right)=\left(\mathrm{1}-P\right)\left[{\mathit{\theta }}^{\prime \prime }\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\left(y\right)-{g}_{\mathrm{0}}\left(y\right)\right]\\ \text{(30)}& & \phantom{\rule{1em}{0ex}}+p\left[{\mathit{\theta }}^{\prime \prime }-\mathit{\text{Pr}}f\phantom{\rule{0.125em}{0ex}}\mathit{\theta }\right]=\mathrm{0},\end{array}$

where $p\in \left[\mathrm{0}\phantom{\rule{0.125em}{0ex}},\phantom{\rule{0.125em}{0ex}}\mathrm{1}\right]$ is an embedded parameter. For p=0 and p=1 one has:

$\begin{array}{ll}& f\left(y,\mathrm{0}\right)\phantom{\rule{0.125em}{0ex}}={f}_{\mathrm{0}}\left(y\right),\phantom{\rule{1em}{0ex}}f\left(y,\mathrm{1}\right)\phantom{\rule{0.125em}{0ex}}=f\left(y\right)\\ \text{(31)}& & \phantom{\rule{1em}{0ex}}\mathit{\theta }\left(y\phantom{\rule{0.125em}{0ex}},\phantom{\rule{0.125em}{0ex}}\mathrm{0}\right)\phantom{\rule{0.125em}{0ex}}={\mathit{\theta }}_{\mathrm{0}}\left(y\right),\phantom{\rule{1em}{0ex}}\mathit{\theta }\left(y,\mathrm{1}\right)\phantom{\rule{0.125em}{0ex}}=\mathit{\theta }\left(y\right)\end{array}$

Note that when p increases from 0 to 1, f(y  , p), θ(y,p) varies from f0(y), θ0(y) to f(y), θ(y).

$\begin{array}{ll}& f\left(y\right)={f}_{\mathrm{0}}\left(y\right)+p\phantom{\rule{0.125em}{0ex}}{f}_{\mathrm{1}}\left(y\right)+{p}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{f}_{\mathrm{2}}\left(y\right)+\mathrm{\dots }\\ & \phantom{\rule{1em}{0ex}}=\sum _{i=\mathrm{0}}^{n}{p}^{i}\phantom{\rule{0.125em}{0ex}}{f}_{i}\left(y\right),{g}_{\mathrm{0}}=\mathrm{0}\\ & \phantom{\rule{1em}{0ex}}\mathit{\theta }\left(y\right)={\mathit{\theta }}_{\mathrm{0}}\left(y\right)+p\phantom{\rule{0.125em}{0ex}}{\mathit{\theta }}_{\mathrm{1}}\left(\mathit{\eta }\right)+{p}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathit{\theta }}_{\mathrm{2}}\left(y\right)+\mathrm{\dots }\\ \text{(32)}& & \phantom{\rule{1em}{0ex}}=\sum _{i=\mathrm{0}}^{n}{p}^{i}\phantom{\rule{0.125em}{0ex}}{\mathit{\theta }}_{i}\left(y\right),{g}_{\mathrm{0}}=\mathrm{0}\end{array}$

By substitution of Eq. (32) into the Eqs. (18) and (26), rearranging based on powers of p-terms and solving these equations result in the following terms:

$\begin{array}{ll}\text{(33)}& & {f}_{\mathrm{0}}\left(y\right)=\mathrm{0.5}\phantom{\rule{0.125em}{0ex}}{y}^{\mathrm{3}}-\mathrm{1.5}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}{y}^{\mathrm{2}},& {f}_{\mathrm{1}}\left(y\right)=\mathrm{0.007440476190}\mathit{\text{De}}\phantom{\rule{0.125em}{0ex}}{y}^{\mathrm{9}}\\ & -\mathrm{0.001785714286}\left(\mathrm{36}.\phantom{\rule{0.125em}{0ex}}+\mathrm{2.0}{\mathit{\text{Re}}}_{\mathrm{w}}\right){y}^{\mathrm{7}}\\ & -\mathrm{0.0125000}\left(\mathrm{2.0}{M}^{\mathrm{2}}-\mathrm{9}.\phantom{\rule{0.125em}{0ex}}\right){y}^{\mathrm{5}}\\ & +\mathrm{0.1666666667}\left(\mathrm{0.300}{M}^{\mathrm{2}}-\mathrm{0.3714285714}\right\\ & +\mathrm{0.06428571429}{\mathit{\text{Re}}}_{\mathrm{w}}){y}^{\mathrm{3}}\\ \text{(34)}& & +\left(-\mathrm{0.025}{M}^{\mathrm{2}}+\mathrm{0.00625}\phantom{\rule{0.125em}{0ex}}-\mathrm{0.007142857143}{\mathit{\text{Re}}}_{\mathrm{w}}\right)y& {f}_{\mathrm{2}}\left(y\right)=-\mathrm{0.000040881}{\mathit{\text{De}}}^{\mathrm{2}}{x}^{\mathrm{15}}-{\mathrm{1.16550116610}}^{-\mathrm{25}}\\ & \left(-{\mathrm{3.013392810}}^{\mathrm{21}}{\mathit{\text{De}}}^{\mathrm{2}}-{\mathrm{2.23214210}}^{\mathrm{20}}{\mathit{\text{DeRe}}}_{\mathrm{w}}\right){x}^{\mathrm{13}}\\ & -{\mathrm{2.52525210}}^{-\mathrm{25}}\left(-{\mathrm{3.8678510}}^{\mathrm{21}}{\mathit{\text{DeM}}}^{\mathrm{2}}-{\mathrm{7.907110}}^{\mathrm{21}}{\mathit{\text{De}}}^{\mathrm{2}}\right\\ & -{\mathrm{3.18214210}}^{\mathrm{21}}{\mathit{\text{DeRe}}}_{\mathrm{w}}-{\mathrm{4.2857110}}^{\mathrm{19}}{\mathit{\text{Re}}}_{\mathrm{w}}^{\mathrm{2}}){x}^{\mathrm{11}}\\ & -{\mathrm{6.61375661410}}^{-\mathrm{25}}\left(\begin{array}{l}{\mathrm{1.035010}}^{\mathrm{22}}{\mathit{\text{DeM}}}^{\mathrm{2}}+{\mathrm{2.2510}}^{\mathrm{20}}{M}^{\mathrm{2}}{\mathit{\text{Re}}}_{\mathrm{w}}\\ +{\mathrm{2.518310}}^{\mathrm{22}}{\mathit{\text{De}}}^{\mathrm{2}}+{\mathrm{1.11810}}^{\mathrm{22}}\mathit{\text{De}}\phantom{\rule{0.125em}{0ex}}{\mathit{\text{Re}}}_{\mathrm{w}}\\ +{\mathrm{4.5010}}^{\mathrm{20}}{\mathit{\text{Re}}}_{\mathrm{w}}^{\mathrm{2}}\end{array}\right){x}^{\mathrm{9}}\\ & -{\mathrm{2.38095230110}}^{-\mathrm{24}}\left(\begin{array}{l}-{\mathrm{4.72510}}^{\mathrm{21}}{\mathit{\text{DeM}}}^{\mathrm{2}}+{\mathrm{4.510}}^{\mathrm{20}}{M}^{\mathrm{2}}{\mathit{\text{Re}}}_{\mathrm{w}}\\ -{\mathrm{3.326710}}^{\mathrm{21}}{\mathit{\text{De}}}^{\mathrm{2}}-{\mathrm{4.28928410}}^{\mathrm{21}}{\mathit{\text{DeRe}}}_{\mathrm{w}}\\ -{\mathrm{6.42857143010}}^{\mathrm{19}}{\mathit{\text{Re}}}^{\mathrm{2}}+{\mathrm{2.510}}^{\mathrm{20}}{M}^{\mathrm{4}}\end{array}\right){x}^{\mathrm{7}}\\ & -{\mathrm{1.6666710}}^{-\mathrm{23}}\left(-{\mathrm{1.5010}}^{\mathrm{20}}{M}^{\mathrm{4}}+{\mathrm{1.0854210}}^{\mathrm{21}}{\mathit{\text{DeM}}}^{\mathrm{2}}\right\\ & -{\mathrm{3.214257710}}^{\mathrm{19}}{M}^{\mathrm{2}}{\mathit{\text{Re}}}_{\mathrm{w}}-{\mathrm{8.91987910}}^{\mathrm{20}}{\mathit{\text{De}}}^{\mathrm{2}}+){x}^{\mathrm{5}}\\ & \mathrm{0.166666666}\left(\begin{array}{l}-\mathrm{0.019285}{M}^{\mathrm{4}}+\mathrm{0.14962}{\mathit{\text{DeM}}}^{\mathrm{2}}+\mathrm{0.0164}{M}^{\mathrm{2}}{\mathit{\text{Re}}}_{\mathrm{w}}\\ +\mathrm{0.007946}{\mathit{\text{De}}}^{\mathrm{2}}+\mathrm{0.0104}{\mathit{\text{DeRe}}}_{\mathrm{w}}\\ +\mathrm{0.0040630}{\mathit{\text{Re}}}_{\mathrm{w}}^{\mathrm{2}}\end{array}\right){x}^{\mathrm{3}}\\ & +\left(\mathrm{0.001309}{M}^{\mathrm{4}}-\mathrm{0.0122}\phantom{\rule{0.125em}{0ex}}\mathit{\text{De}}\phantom{\rule{0.125em}{0ex}}{M}^{\mathrm{2}}\right\\ & -\mathrm{0.002053}{M}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathit{\text{Re}}}_{\mathrm{w}}-\mathrm{0.009762}{\mathit{\text{De}}}^{\mathrm{2}}\\ \text{(35)}& & -\mathrm{0.001906}\mathit{\text{DeRe}}-\mathrm{0.00054}{\mathit{\text{Re}}}_{\mathrm{w}}^{\mathrm{2}})x& \mathrm{⋮}\\ \text{(36)}& & {\mathit{\theta }}_{\mathrm{0}}\left(y\right)=\mathrm{1}\text{(37)}& & {\mathit{\theta }}_{\mathrm{1}}\left(y\right)=\phantom{\rule{0.125em}{0ex}}-\phantom{\rule{0.125em}{0ex}}\mathrm{0.025}\mathit{\text{Pr}}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}{y}^{\mathrm{5}}+\mathrm{0.25}\phantom{\rule{0.125em}{0ex}}\mathit{\text{Pr}}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}{y}^{\mathrm{3}}-\mathrm{0.225}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\mathit{\text{Pr}}& {\mathit{\theta }}_{\mathrm{2}}\left(y\right)=\mathrm{0.0000676}\\ & \mathit{\text{Pr}}\phantom{\rule{0.125em}{0ex}}\mathit{\text{De}}\phantom{\rule{0.125em}{0ex}}{y}^{\mathrm{11}}+\mathrm{0.0001388}{\mathit{\text{Pr}}}^{\mathrm{2}}{y}^{\mathrm{10}}-\mathrm{0.00089283}\phantom{\rule{0.125em}{0ex}}\mathit{\text{Pr}}\phantom{\rule{0.125em}{0ex}}\mathit{\text{De}}\phantom{\rule{0.125em}{0ex}}{y}^{\mathrm{9}}\\ & -\mathrm{0.000049603}\mathit{\text{Pr}}\phantom{\rule{0.125em}{0ex}}{x}^{\mathrm{9}}{\mathit{\text{Re}}}_{\mathrm{w}}\\ & -\mathrm{0.00285714}{y}^{\mathrm{8}}{\mathit{\text{Pr}}}^{\mathrm{2}}-\mathrm{0.000595552}\mathit{\text{Pr}}\phantom{\rule{0.125em}{0ex}}{x}^{\mathrm{7}}{M}^{\mathrm{2}}\\ & +\mathrm{0.002645789}\mathit{\text{Pr}}\phantom{\rule{0.125em}{0ex}}\mathit{\text{De}}\phantom{\rule{0.125em}{0ex}}{y}^{\mathrm{7}}-\mathrm{0.0030959651}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\mathit{\text{Pr}}\phantom{\rule{0.125em}{0ex}}\mathit{\text{De}}\phantom{\rule{0.125em}{0ex}}{y}^{\mathrm{5}}\\ & +\mathrm{0.01250000}{x}^{\mathrm{6}}{\mathit{\text{Pr}}}^{\mathrm{2}}+\mathrm{0.00250000}\mathit{\text{Pr}}\phantom{\rule{0.125em}{0ex}}{x}^{\mathrm{5}}{M}^{\mathrm{2}}\\ & +\mathrm{0.00562544000}{y}^{\mathrm{5}}{\mathit{\text{Pr}}}^{\mathrm{2}}+\mathrm{0.0010467104}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\mathit{\text{Pr}}\phantom{\rule{0.125em}{0ex}}\mathit{\text{De}}\phantom{\rule{0.125em}{0ex}}{y}^{\mathrm{3}}\\ & \mathrm{0.000535}\mathit{\text{Pr}}\phantom{\rule{0.125em}{0ex}}{y}^{\mathrm{5}}{\mathit{\text{Re}}}_{\mathrm{w}}-\mathrm{0.004167}\mathit{\text{Pr}}\phantom{\rule{0.125em}{0ex}}{M}^{\mathrm{2}}{y}^{\mathrm{3}}\\ & -\mathrm{0.05625}{\mathit{\text{Pr}}}^{\mathrm{2}}{y}^{\mathrm{3}}+\mathrm{0.000200507}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\mathit{\text{Pr}}\phantom{\rule{0.125em}{0ex}}\mathit{\text{De}}\\ & +\mathrm{0.0408878}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}{\mathit{\text{Pr}}}^{\mathrm{2}}\\ & -\mathrm{0.001190470}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\mathit{\text{Pr}}\phantom{\rule{0.125em}{0ex}}{\mathit{\text{Re}}}_{\mathrm{w}}{y}^{\mathrm{3}}+\mathrm{0.0022619}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\mathit{\text{Pr}}\phantom{\rule{0.125em}{0ex}}{M}^{\mathrm{2}}\\ \text{(38)}& & +\mathrm{0.000703}\phantom{\rule{0.125em}{0ex}}\mathit{\text{Pr}}\phantom{\rule{0.125em}{0ex}}{\mathit{\text{Re}}}_{\mathrm{w}}& \mathrm{⋮}\end{array}$

The solution of this equation, when p→1, is as follows:

$\begin{array}{}\text{(39)}& & f\left(y\right)=\sum _{i=\mathrm{0}}^{\mathrm{9}}{\mathrm{Lim}}_{p\to \mathrm{1}}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}{p}^{i}\phantom{\rule{0.125em}{0ex}}{f}_{i}\left(y\right)\text{(40)}& & \mathit{\theta }\left(y\right)=\sum _{i=\mathrm{0}}^{\mathrm{9}}{\mathrm{Lim}}_{p\to \mathrm{1}}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}{p}^{i}\phantom{\rule{0.125em}{0ex}}{\mathit{\theta }}_{i}\left(y\right)\end{array}$

## 3.2 Implementation of the Homotopy Analysis Method

For HAM solutions, the auxiliary linear operators are chosen in the following form:

$\begin{array}{}\text{(41)}& & L\left(f\right)={f}^{\prime \prime \prime \prime },\text{(42)}& & L\left(\mathit{\theta }\right)={\mathit{\theta }}^{\prime \prime },\end{array}$

The initially guessed function is obtained by solving the following equations:

$\begin{array}{}\text{(43)}& & L\phantom{\rule{0.125em}{0ex}}\left(\phantom{\rule{0.125em}{0ex}}\frac{\mathrm{1}}{\mathrm{6}}{c}_{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{y}^{\mathrm{3}}+\frac{\mathrm{1}}{\mathrm{2}}{c}_{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{y}^{\mathrm{2}}+{c}_{\mathrm{3}}\phantom{\rule{0.125em}{0ex}}y+{c}_{\mathrm{4}}\right)=\mathrm{0},\text{(44)}& & L\phantom{\rule{0.125em}{0ex}}\left(\phantom{\rule{0.125em}{0ex}}{c}_{\mathrm{5}}\phantom{\rule{0.125em}{0ex}}y+{c}_{\mathrm{6}}\right)=\mathrm{0},\end{array}$

with boundary Eqs. (17), (27) and (28) as:

$\begin{array}{}\text{(45)}& & {f}_{\mathrm{0}}\left(y\right)=-\frac{\mathrm{1}}{\mathrm{2}}{y}^{\mathrm{3}}+\frac{\mathrm{3}}{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}y,\text{(46)}& & {\mathit{\theta }}_{\mathrm{0}}\left(y\right)=\mathrm{1},\end{array}$

Let $P\in \left[\mathrm{0},\phantom{\rule{0.125em}{0ex}}\mathrm{1}\right]$ denote the embedded parameter and indicate non–zero auxiliary parameters. Consequently the following equations are constructed.

Zeroth – order deformation equations

$\begin{array}{ll}\text{(47)}& & \left(\mathrm{1}-P\right)L\left[F\left(y;p\right)-{f}_{\mathrm{0}}\left(y\right)\right]=p\mathrm{\hslash }H\left(y\right)N\left[F\left(y;p\right)\right]\text{(48)}& & \left(\mathrm{1}-P\right)L\left[\mathrm{\Theta }\left(y;p\right)-{\mathit{\theta }}_{\mathrm{0}}\left(y\right)\right]=p\mathrm{\hslash }H\left(y\right)N\left[\mathrm{\Theta }\left(y;p\right)\right]& F\left(\mathrm{0};p\right)=\mathrm{1};\phantom{\rule{1em}{0ex}}{F}^{\prime \prime }\left(\mathrm{0};p\right)=\mathrm{0},\\ \text{(49)}& & \phantom{\rule{2em}{0ex}}F\left(\mathrm{1};p\right)=\mathrm{1},\phantom{\rule{1em}{0ex}}{F}^{\prime }\left(\mathrm{1};p\right)=\mathrm{0}\text{(50)}& & \mathit{\theta }\left(\mathrm{1};p\right)=\mathrm{1};\phantom{\rule{1em}{0ex}}{\mathit{\theta }}^{\prime }\left(\mathrm{0};p\right)=\mathrm{0}& N\left[F\left(y;p\right)\right]=\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\frac{{\mathrm{d}}^{\mathrm{4}}F\left(y;p\right)}{\mathrm{d}{y}^{\mathrm{4}}}\\ & \phantom{\rule{1em}{0ex}}+\phantom{\rule{0.125em}{0ex}}\mathit{\text{Re}}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\left[\frac{\mathrm{d}F\left(y;p\right)}{\mathrm{d}y}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\frac{{\mathrm{d}}^{\mathrm{2}}F\left(y;p\right)}{\mathrm{d}{y}^{\mathrm{2}}}\right\\ & \phantom{\rule{1em}{0ex}}-\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}F\left(y;p\right)\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\frac{{\mathrm{d}}^{\mathrm{3}}F\left(y;p\right)}{\mathrm{d}{y}^{\mathrm{3}}}]\\ & \phantom{\rule{1em}{0ex}}+\mathit{\text{De}}\left[\mathrm{2}{\left(\frac{\mathrm{d}F\left(y;p\right)}{\mathrm{d}y}\right)}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}\frac{{\mathrm{d}}^{\mathrm{2}}F\left(y;p\right)}{\mathrm{d}{y}^{\mathrm{2}}}\right\\ & \phantom{\rule{1em}{0ex}}-\mathrm{2}F\left(y;p\right){\left(\frac{{\mathrm{d}}^{\mathrm{2}}F\left(y;p\right)}{\mathrm{d}{y}^{\mathrm{2}}}\right)}^{\mathrm{2}}-F\left(y;p{\right)}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}\frac{{\mathrm{d}}^{\mathrm{4}}F\left(y;p\right)}{\mathrm{d}{y}^{\mathrm{4}}}]\\ \text{(51)}& & +\left({M}^{\mathrm{2}}\right)\frac{{\mathrm{d}}^{\mathrm{2}}F\left(y;p\right)}{\mathrm{d}{y}^{\mathrm{2}}}\text{(52)}& & N\left[\mathit{\theta }\left(y;p\right)\right]=\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\frac{{\mathrm{d}}^{\mathrm{2}}\mathit{\theta }\left(y;p\right)}{\mathrm{d}{y}^{\mathrm{2}}}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}+\mathit{\text{Pr}}\phantom{\rule{0.125em}{0ex}}F\left(y;p\right)\phantom{\rule{0.125em}{0ex}}\frac{\mathrm{d}\mathit{\theta }\left(y;p\right)}{\mathrm{d}y}\end{array}$

For p=0 and p=1, one has:

$\begin{array}{}\text{(53)}& & F\left(y;\mathrm{0}\right)={f}_{\mathrm{0}}\left(y\right)\phantom{\rule{2em}{0ex}}F\left(y;\mathrm{1}\right)=f\left(y\right)\text{(54)}& & \mathit{\theta }\left(y;\mathrm{0}\right)={\mathit{\theta }}_{\mathrm{0}}\left(y\right)\phantom{\rule{2em}{0ex}}\mathit{\theta }\left(y;\mathrm{1}\right)=\mathit{\theta }\left(y\right)\end{array}$

As p increases from 0 to 1 then F(y;p), θ(y;p) varies from f0(y), θ0(y) to f(y), θ(y). By Taylor's theorem, F(y;p), θ(y;p) can be expanded in a power series of p as follows:

$\begin{array}{ll}& F\left(y;p\right)={f}_{\mathrm{0}}\left(y\right)+\sum _{m-\mathrm{1}}^{\mathrm{\infty }}{f}_{\mathrm{m}}\left(y\right){p}^{m}\\ \text{(55)}& & \phantom{\rule{2em}{0ex}}{f}_{\mathrm{m}}\left(y\right)={\frac{\mathrm{1}}{m}\frac{{\partial }^{m}\left(F\left(y;p\right)\right)}{\partial {p}^{m}}|}_{p=\mathrm{0}}& \mathit{\theta }\left(y;p\right)={\mathit{\theta }}_{\mathrm{0}}\left(y\right)+\sum _{m-\mathrm{1}}^{\mathrm{\infty }}{\mathit{\theta }}_{\mathrm{m}}\left(y\right){p}^{m}\\ \text{(56)}& & \phantom{\rule{2em}{0ex}}{\mathit{\theta }}_{\mathrm{m}}\left(y\right)={\frac{\mathrm{1}}{m}\frac{{\partial }^{m}\left(\mathit{\theta }\left(y;p\right)\right)}{\partial {p}^{m}}|}_{p=\mathrm{0}}\end{array}$

there is chosen in such a way that this series is convergent at p=1, therefore through Eqs. (55) and (56), it can be concluded that:

$\begin{array}{}\text{(57)}& & f\left(y\right)={f}_{\mathrm{0}}\left(y\right)+\sum _{m-\mathrm{1}}^{\mathrm{\infty }}{f}_{\mathrm{m}}\left(y\right),\text{(58)}& & \mathit{\theta }\left(y\right)={\mathit{\theta }}_{\mathrm{0}}\left(y\right)+\sum _{m-\mathrm{1}}^{\mathrm{\infty }}{\mathit{\theta }}_{\mathrm{m}}\left(y\right),\end{array}$

mth – order deformation equations

The mth-order deformation equations are written as:

$\begin{array}{ll}\text{(59)}& & L\left[{f}_{\mathrm{m}}\left(y\right)-\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}{\mathit{\chi }}_{\mathrm{m}}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}{f}_{m-\mathrm{1}}\left(y\right)\right]=\mathrm{\hslash }\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}H\left(y\right)\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}{R}_{\mathrm{m}}^{f}\left(y\right)\text{(60)}& & L\left[{\mathit{\theta }}_{\mathrm{m}}\left(y\right)-\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}{\mathit{\chi }}_{\mathrm{m}}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}{\mathit{\theta }}_{m-\mathrm{1}}\left(y\right)\right]=\mathrm{\hslash }\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}H\left(y\right)\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}{R}_{\mathrm{m}}^{\mathit{\theta }}\left(y\right)& {F}_{\mathrm{m}}\left(\mathrm{0};p\right)=\mathrm{0};\phantom{\rule{1em}{0ex}}{F}_{\mathrm{m}}^{\prime \prime }\left(\mathrm{0};p\right)=\mathrm{0},\\ \text{(61)}& & \phantom{\rule{1em}{0ex}}{F}_{\mathrm{m}}\left(\mathrm{1};p\right)=\mathrm{0},\phantom{\rule{1em}{0ex}}{F}_{\mathrm{m}}^{\prime }\left(\mathrm{1};p\right)=\mathrm{0}\text{(62)}& & {\mathit{\theta }}_{\mathrm{m}}\left(\mathrm{1};p\right)=\mathrm{0};\phantom{\rule{1em}{0ex}}{\mathit{\theta }}_{\mathrm{m}}^{\prime }\left(\mathrm{0};p\right)=\mathrm{0}& {R}_{\mathrm{m}}^{f}\left(y\right)={{f}^{\prime \prime }}_{m-\mathrm{1}}^{\prime \prime }+\sum _{k=\mathrm{0}}^{m-\mathrm{1}}\left[\mathit{\text{Re}}\left(\phantom{\rule{0.125em}{0ex}}{f}_{m-\mathrm{1}-k}^{\prime }\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}{f}_{k}^{\prime \prime }-{f}_{m-\mathrm{1}-k}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}{f}_{k}^{\prime \prime \prime }\right)\right\\ & \phantom{\rule{1em}{0ex}}+\mathit{\text{De}}\phantom{\rule{0.125em}{0ex}}{f}_{m-\mathrm{1}-k}^{\prime }\left(\sum _{l=\mathrm{0}}^{k}\phantom{\rule{0.25em}{0ex}}\left(\phantom{\rule{0.125em}{0ex}}\mathrm{2}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}{f}_{k-l}^{\prime }\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}{f}_{l}^{\prime \prime }\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\right)\right)\\ & -\mathit{\text{De}}\phantom{\rule{0.125em}{0ex}}{f}_{m-\mathrm{1}-k}\left(\sum _{l=\mathrm{0}}^{k}\phantom{\rule{0.25em}{0ex}}\left(\phantom{\rule{0.125em}{0ex}}\mathrm{2}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}{f}_{k-l}^{\prime \prime }\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}{f}_{l}^{\prime \prime }\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}+{f}_{k-l}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}{f}_{l}^{\prime \prime \prime \prime }\right)\right)]\\ \text{(63)}& & \phantom{\rule{1em}{0ex}}+\left({M}^{\mathrm{2}}\right)\phantom{\rule{0.125em}{0ex}}{f}_{m-\mathrm{1}}^{\prime \prime }\text{(64)}& & {R}_{\mathrm{m}}^{\mathit{\theta }}\left(y\right)={\mathit{\theta }}_{m-\mathrm{1}}^{\prime \prime }-\sum _{k=\mathrm{0}}^{m-\mathrm{1}}\left[\mathit{\text{Pr}}\left(\phantom{\rule{0.125em}{0ex}}{f}_{m-\mathrm{1}-k}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}{\mathit{\theta }}_{k}^{\prime }\right)\right]\text{(65)}& & {\mathit{\chi }}_{\mathrm{m}}=\left\{\begin{array}{l}\mathrm{0},\phantom{\rule{2em}{0ex}}m\le \mathrm{1}\\ \mathrm{1},\phantom{\rule{2em}{0ex}}m>\mathrm{1}\end{array}\right\\end{array}$

For simplicity, it is assumed that:

$\begin{array}{}\text{(66)}& H\left(y\right)=\mathrm{1}\end{array}$

For different values of m, the solution is obtained by maple analytic solution device. The first deformation is presented as:

$\begin{array}{ll}& {f}_{\mathrm{1}}\left(y\right)=-\frac{\mathrm{5}}{\mathrm{672}}{\mathrm{\hslash }}_{\mathrm{1}}\mathit{\text{De}}{y}^{\mathrm{9}}-\frac{\mathrm{3}}{\mathrm{2}}{\mathrm{\hslash }}_{\mathrm{1}}\left(-\frac{\mathrm{1}}{\mathrm{420}}{\mathit{\text{Re}}}_{\mathrm{w}}-\frac{\mathrm{3}}{\mathrm{70}}\mathit{\text{De}}\right){y}^{\mathrm{7}}\\ & \phantom{\rule{1em}{0ex}}-\frac{\mathrm{3}}{\mathrm{2}}{\mathrm{\hslash }}_{\mathrm{1}}\left(\frac{\mathrm{3}}{\mathrm{40}}B-\frac{\mathrm{1}}{\mathrm{60}}{M}^{\mathrm{2}}\right){y}^{\mathrm{5}}\\ & \phantom{\rule{1em}{0ex}}+\left(\frac{\mathrm{13}}{\mathrm{210}}{\mathrm{\hslash }}_{\mathrm{1}}\mathit{\text{De}}-\frac{\mathrm{3}}{\mathrm{280}}{\mathrm{\hslash }}_{\mathrm{1}}{\mathit{\text{Re}}}_{\mathrm{w}}-\frac{\mathrm{1}}{\mathrm{20}}{\mathrm{\hslash }}_{\mathrm{1}}{M}^{\mathrm{2}}\right){y}^{\mathrm{3}}\\ \text{(67)}& & \phantom{\rule{1em}{0ex}}+\left(-\frac{\mathrm{1}}{\mathrm{160}}{\mathrm{\hslash }}_{\mathrm{1}}\mathit{\text{De}}+\frac{\mathrm{1}}{\mathrm{140}}{\mathrm{\hslash }}_{\mathrm{1}}{\mathit{\text{Re}}}_{\mathrm{w}}+\frac{\mathrm{1}}{\mathrm{40}}{h}_{\mathrm{1}}{M}^{\mathrm{2}}\right)y\text{(68)}& & {\mathit{\theta }}_{\mathrm{1}}\left(y\right)=\frac{\mathrm{1}}{\mathrm{40}}{\mathrm{\hslash }}_{\mathrm{2}}\mathit{\text{Pr}}{y}^{\mathrm{5}}-\frac{\mathrm{1}}{\mathrm{4}}{\mathrm{\hslash }}_{\mathrm{2}}\mathit{\text{Pr}}{y}^{\mathrm{3}}+\frac{\mathrm{9}}{\mathrm{40}}{\mathrm{\hslash }}_{\mathrm{2}}\mathit{\text{Pr}}\end{array}$

The solutions of f(y),θ(y) are too long and will be shown in obtained results.

### 3.2.1 Convergence of the HAM solution

As pointed out by Liao (2012), the convergence region and rate of solution series can be adjusted and controlled by means of the auxiliary parameter . In order to check the convergence of the present solution, the so-called 1-curve of f′′(1) is depicted in Figs. 2 and 3 and the 2-curve of θ(1) is shown in Fig. 4.

The solutions converge for values which correspond to the horizontal line segment in the curve. In order to investigate the range of admissible values of the auxiliary parameter , the -curve of f′′(1) is illustrated for various quantities of De, M and Rew in Fig. 3. Similarly, the 2-curve of θ(1) is shown in Fig. 4. As seen clearly in Figs. 2–4, it can be concluded that ${\mathrm{\hslash }}_{\mathrm{1}}={\mathrm{\hslash }}_{\mathrm{2}}=-\mathrm{1}$ are suitable values for different quantities within $\mathrm{0.1}<\mathit{\text{De}}<\mathrm{0.6}$, $\mathrm{0} and $-\mathrm{5}<{\mathit{\text{Re}}}_{\mathrm{w}}<\mathrm{5}\left(\le T\le \mathrm{6}\right)$.

4 Numerical method

The above system of non-linear ordinary differential Eqs. (17) and (23) along with the boundary conditions (16), (24) and (25) are solved numerically using the algebra package Maple 16.0. The package uses a boundary value (B-V) problem procedure. The algorithm can be used to find moderate accurate solutions for ODE boundary and initial value problems, both with a global error bound. The method uses either Richardson extrapolation or deferred corrections with a base method of either the trapezoid or midpoint method. The trapezoid method is generally efficient for typical problems, while the midpoint method is a powerful method for solving harmless end-point singularities that the trapezoid method cannot. The midpoint method, also known as the fourth-order Runge–Kutta–Fehlberg method, improves the Euler method by adding a midpoint in the step which increases the accuracy by one order. Thus, the midpoint method is used as a suitable numerical technique (Hatami et al., 2013; Aziz, 2006).

5 Results and discussion

In the present study HAM and HPM methods are applied to obtain an explicit analytic solution of the MHD flow of an UCM fluid in a channel (Fig. 1). Comparison between the results obtained by the numerical method (Runge-Kutta-Fehlberg technique), Homotopy Perturbation and Homotopy Analysis Methods for different values of active parameters is shown in Fig. 5. The results are proved to be precise and accurate in solving a wide range of mathematical and engineering problems especially Fluid mechanics cases. This investigation is completed by depicting the effects of some important parameters to determine how these parameters affect the fluid flow and heat transfer. From a physical point of view, Figs. 6–10 are presented in order to highlight the effects of the Deborah number De, Hartman number M and Prandtl number Pr on the velocity and temperature profiles.

Figure 6 shows the effect of Deborah number on the velocity distribution for M=1, Pr=1, Rew=3. This number may be interpreted as the ratio of the relaxation time, and the characteristic time of an experiment or a computer simulation probing the response of the material (Reiner, 1964). Higher De values imply a strongly elastic behavior. The Newtonian fluids possess no relaxation time (i.e., De= 0). An increase in the elastic parameter (i.e. as De rises from 0 to 0.6) results in a decrease of the velocity components at any given point (Hayat et al., 2009; Mostafa and Mahmoud, 2011). Figure 7 shows the effect of a Hartman number M on the velocity components for De=0.2, Pr=1, Rew=3. An increase of the magnetic parameter leads to a decrease in the velocity components and contours at a given point. This is due to the fact that the applied transverse magnetic field produces damping in form of a Lorentz force thereby decreasing the magnitude of the velocity. The drop in velocity as a consequence of an increase in the strength of the magnetic field is observed.

Figures 8 show the variation of temperature for different values of the viscoelastic parameter Deborah number De. From Fig. 8, it is found that the temperature field increases with an increase of the Deborah number De. It is remarkable to note that, since the Deborah number has depended on the relaxation time, the larger values of Deborah number imply to higher relaxation time. It is well known fact that the larger relaxation time fluids leads to a higher temperature (Khan et al., 2016).

Also, the effects of Hartman parameter M on the temperature profile have been illustrated in Fig. 9. It can be inferred from this figure that the temperature increases with an increase in the magnetic field. Application of a magnetic field normal to an electrically conducting fluid has the tendency to produce a resistive type of body force called the Lorentz force which opposes the fluid motion, causing a flow retardation effect. This causes the fluid velocity to decrease. However, this decrease in flow speed is accompanied by corresponding increases in the fluid thermal state level (Mostafa and Mahmoud, 2011). This can be attributed to the fact that, u is small and an increase in Ha increases the Joule dissipation which is proportional to Ha and therefore, the temperature increases.

Figure 10 shows the effect of Prandtl number on the dimensionless temperature profile. As observed, an increase in the Prandtl number leads to a decrease in the temperature. This is in agreement with the physical fact that the thermal boundary layer thickness decreases with increasing Pr.

6 Conclusions

In this investigation, the analytical approaches called Homotopy Analysis Method (HAM), Homotopy Perturbation Method (HPM) have been successfully applied to find the most accurate analytical solution for the velocity and temperature distributions of MHD flow of an upper-convected Maxwell fluid in a channel. Furthermore, the obtained solutions by the proposed method have been compared with the direct numerical solutions using the Runge-Kutta-Fehlberg technique. Effects of different physical parameters, such as, De Deborah number, M Hartman number and Pr Prandtl number on the velocity and temperature profiles have been investigated. The main conclusions are as follows:

• The comparison shows that the HAM and HPM solutions are highly accurate and provide rapid achievement in computing the flow characteristics. According to the previous publications this method is a powerful technique for finding analytical solutions in science and engineering problems.

• The results show that an increase in the M Hartman number is associated with a reduction in the velocity profile. Moreover, it is worthwhile to mention that, the velocity patterns are marginally influenced by the changes in De Deborah number parameter.

• Also it is observed that an increase of Pr Prandtl number results in decreasing the thermal boundary layer thickness and a more uniform temperature distribution across the boundary layer. The reason is that smaller values of Pr are equivalent to increasing the thermal conductivities, and therefore heat is able to diffuse away more rapidly than for higher values of Pr. This means that the heat loss increases for larger Pr as the boundary layer gets thinner.

Data availability
Data availability.

All the data used in this manuscript can be obtained by requesting from the corresponding author.

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Edited by: Ali Konuralp
Reviewed by: two anonymous referees

References

Abbasi, M. and Rahimipetroudi, I.: Analytical solution of an upper- convective Maxwell fluid in porous channel with slip at the boundaries by using the Homotopy Perturbation Method, IJNDES, 5, 7–17, 2013.

Abbasi, M., Ganji, D. D., Rahimipetroudi, I., and Khaki, M.: Comparative Analysis of MHD Boundary-Layer Flow of Viscoelastic Fluid in Permeable Channel with Slip Boundaries by using HAM, VIM, HPM, Walailak, J. Sci. Technol., 11, 551–567, 2014.

Abbasi, M., Khaki, M., Rahbari, A., Ganji, D. D., and Rahimipetroudi, I.: Analysis of MHD flow characteristics of an UCM viscoelastic flow in a permeable channel under slip conditions, J. Braz. Soc. Mech. Sci. Eng., 38, 977–988, 2016.

Abel, M. S., Tawade, J. V., and Nandeppanavar, M. M.: MHD flow and heat transfer for the upper-convected Maxwell fluid over a stretching sheet, Meccanica, 47, 385–393, 2012.

Adegbie, K. S., Omowaye, A.J., Disu, A. B., and Animasaun, I. L.: Heat and Mass Transfer of Upper Convected Maxwell Fluid Flow with Variable Thermo-Physical Properties over a Horizontal Melting Surface, Appl. Math., 6, 1362–1379, 2015.

Afify, A. A. and Elgazery, N. S.: Effect of a chemical reaction on magnetohydrodynamic boundary layer flow of a Maxwell fluid over a stretching sheet with nanoparticles, Particuology, 29, 154–161, 2016.

Anwar Bég, O. and Makinde, O. D.: Viscoelastic flow and species transfer in a Darcian high-permeability channel, J. Pet. Sci. Eng., 76, 93–99, 2010.

Aziz, A.: Heat Conduction with Maple, edited by: Edwards, R. T., Philadelphia (PA), 2006.

Ganji, D. D., Abbasi, M., Rahimi, J., Gholami, M., and Rahimipetroudi, I.: On the MHD Squeeze Flow between Two Parallel Disks with Suction or Injection via HAM and HPM, Appl. Mech. Mater., 3, 270–280, 2014.

Hatami, M., Nouri, R., and Ganji, D. D.: Forced convection analysis for MHD Al2O3 – water nano fluid flow over a horizontal plate, J. Mol. Liq., 187, 294–301, 2013.

Hayat, T. and Awais, M.: Three-dimensional flow of upper-convected Maxwell (UCM) fluid, Int. J. Numer. Meth. Fl., 66, 875–884, 2011.

Hayat, T. and Wang, Y.: Magnetohydrodynamic flow due to noncoaxial rotations of a porous disk and a fourth-grade fluid infinity, Math. Probl. Eng., 2, 47–64, 2003.

Hayat, T., Abbas, Z., and Sajid, M.: Series solution for the upper-convected Maxwell fluid over a porous stretching plate, Phys. Lett. A, 358, 396–403, 2006.

Hayat, T., Abbas, Z., and Sajid, M.: MHD stagnation-point flow of an upper-convected Maxwell fluid over a stretching surface, Chaos Soliton. Fract., 39, 840–848, 2009.

Hayat, T., Awais, M., Qasim, M., Awatif, A., and Hendi, A.: Effects of mass transfer on the stagnation point flow of an upper-convected Maxwell (UCM) fluid, Int. J. Heat Mass Tran., 54, 3777–3782, 2011a.

Hayat, T., Sajjad, R., Abbas, Z., Sajid, M., Awatif, A., and Hendi, A.: Radiation effects on MHD flow of Maxwell fluid in a channel with porous medium, Int. J. Heat Mass Tran., 54, 854–862, 2011b.

Hayat, T., Mustafa, M., Shehzadand, S. A., and Obaidat, S.: Melting heat transfer in the stagnation-point flow of an upper-convected Maxwell (UCM) fluid past a stretching sheet, Int. J. Numer. Meth. Fl., 68, 233–243, 2012.

Hayat, T., Waqas, M., Ijaz Khan, M., and Alsaedi, A.: Impacts of constructive and destructive chemical reactions in magnetohydrodynamic (MHD) flow of Jeffrey liquid due to nonlinear radially stretched surface, J. Mol. Liq., 225, 302–310, 2017.

He, J. H.: A coupling method of homotopy technique and perturbation technique for nonlinear problems, Int. J. Nonlinear Mech., 35, 37–43, 2000.

Khan, Z., Shah, R. A., Islam, S., Jan, B., Imran, M., and Tahir, F.: Steady flow and heat transfer analysis of Phan-Thein-Tanner fluid in double-layer optical fiber coating analysis with Slip Conditions, Sci. Rep., 6, 34593, https://doi.org/10.1038/srep34593, 2016.

Liao, S. J.: Homotopy Analysis Method in Nonlinear Differential Equation, Berlin and Beijing: Springer & Higher Education Press, 2012.

Mostafa, A. and Mahmoud, A.: The effects of variable fluid properties on mhd maxwell fluids over a stretching surface in the presence of heat generation/absorption, Chem. Eng. Commun., 198, 131–146, 2011.

Mukhopadhyay, S.: Upper-Convected Maxwell Fluid Flow over an Unsteady Stretching Surface Embedded in Porous Medium Subjected to Suction/Blowing, Z. Naturforsch., 67, 641–646, 2012.

Nadeem, S., Ul Haq, R., and Khan, Z. H.: Numerical study of MHD boundary layer flow of a Maxwell fluid past a stretching sheet in the presence of nanoparticles, J. Taiwan Inst. Chem. E., 45, 121–126, 2014.

Raftari, B. and Yildirim, A.: The application of homotopy perturbation method for MHD flows of UCM fluids above porous stretching sheets, Comput. Math. Appl., 59, 3328–3337, 2010.

Recebli, Z., Selimli, S., and Gedik, E.: Three dimensional numerical analysis of magnetic field effect on Convective heat transfer during the MHD steady state laminar flow of liquid lithium in a cylindrical pipe, Comput. Fluids, 15, 410–417, 2013.

Recebli, Z., Selimli, S., and Ozkaymak, M.: Theoretical analyses of immiscible MHD pipe flow, Int. J. Hydrogen Energ., 40, 15365–15373, 2015.

Reiner, M.: The Deborah Number, Phys. Today, 17, 62–65, 1964.

Rivlin, R. S. and Ericksen, J. L.: Stress deformation relation for isotropic materials, J. Rat. Mech. Anal., 4, 323–425, 1955.

Rossow, V. J.: On flow of electrically conducting fluids over a flat plate in the presence of a transverse magnetic field, Tech. Report, 1358 NASA 1958.

Sajid, M., Iqbal, Z., Hayat, T., and Obaidat, S.: Series Solution for Rotating Flow of an Upper Convected Maxwell Fluid over a Stretching Sheet, Commun. Theor. Phys., 56, 740–744, 2011.

Selimli, S., Recebli, Z., and Arcaklioglu, E.: Combined effects of magnetic and electrical field on the hydrodynamic and thermophysical parameters of magnetoviscous fluid flow, Int. J. Heat Mass Tran., 86, 426–432, 2015.

Shereliff, J. A.: A text book of magneto-hydrodynamics, Pergoman Press, place City Oxford, 1965.

Turkyilmazoglu, M.: Numerical and analytical solutions for the flow and heat transfer near the equator of an MHD boundary layer over a porous rotating sphere, Int. J. Therm. Sci., 50, 831–842, 2011.

Wehgal, A. R. and Ashraf, M.: MHD asymmetric flow between two porous disks, Punjab University Journal of Mathematics, 44, 9–21, 2012.