On the Analysis of Variable Thermophysical Properties of Thermophoretic Viscoelastic Fluid Flow past a Vertical Surface with n^th Order of Chemical Reaction

Olubode Kolade Koriko¹, Tosin Oreyeni^1*, Onyekachukwu Anselm Oyem²

¹Department of Mathematical Sciences, Federal University of Technology Akure, Akure, Nigeria.
²Department of Mathematical Sciences, Federal University Lokoja, Lokoja, Nigeria.

DOI: 10.4236/oalib.1104271                

Abstract

The objective of this study is to consider the flow of temperature dependent viscosity and thermal conductivity of free convective heat and mass transfer of viscoelastic fluid over a stretching surface with nth order of chemical reaction and thermophoresis. The effect of the temperature dependent dynamic viscosity and thermal conductivity together with modified thermal and solutal Grashof numbers are properly accounted for in order to enhance the transport phenomenon. Similarity transformations are used to convert and parameterize the non-linear partial differential equation to a system of coupled non-linear ordinary differential equation. The approximate analytical solutions of the corresponding BVP are obtained through Optimal Homotopy Analysis Method (OHAM). The effect of some pertinent parameters is tested on velocity, temperature, concentration profiles. It is observed from the computation that, the thickness of the velocity and thermal boundary layer increases with an increase in temperature dependent variable viscosity and thermal conductivity parameters when modified thermal and solutal Grashof numbers and are less than zero. It is also observed that the concentration layer becomes thinner with increasing thermophoresis parameter when the chemical reaction parameter is greater than zero for both cases of first and second order of chemical reaction i.e. when n=1,2.

Keywords

Viscoelastic Fluid, Thermophoresis, Variable Fluid Properties, Homotopy Analysis Method, Chemical Reaction

1. Introduction

Prandtl [1] introduced the concept of Boundary Layer in 1904, in the lecture titled “Uber Flussigkeitsbewegungen bei sehr kleiner Reibung” translated into English as “On fluid flow with very little friction”. During the course of the lecture, he explained that the viscosity of a fluid plays a role in a (very) thin layer adjacent to the surface. Base on the theory, the fluid flow region around a solid body is divided into two regions; a very thin region in the immediate vicinity of a bounding surface, called the boundary layer, where the effect of viscosity is significant, and the other region outside the boundary layer where the fluid viscosity is negligible. The concept of boundary layers is central to the understanding of convection heat and mass transfer between a surface and a fluid flowing past it. In mass transfer problem involving phase change (evaporation, sublimation, condensation, melting etc.) must also involve heat transfer, and the solution of such problems needs to be analyzed by considering simultaneous heat and mass transfer. Some examples of simultaneous heat and mass problems are drying, evaporating cooling, transpiration cooling, combustion of fuel droplets etc. (Cengel and Ghajar [2] ). Chemical reactions can be classified as either homogeneous or heterogeneous, some mass transfer problems involve chemical reactions that occur within the medium and result in the generation of a species throughout, such reactions that occur within the medium are called homogeneous reactions (Cengel and Ghajar [2] ). The study of heat and mass transfer with chemical reaction is of great practical importance in many branches of science and engineering. The effect of chemical reaction on different geometry of the problem has been examined by many authors. Gangadhar [3] studied chemically reacting mhd boundary layer flow of heat and mass transfer over a moving vertical plate in a porous medium with suction. Ibrahim and Makinde [4] discussed chemically reacting mhd boundary layer flow of heat and mass transfer over a moving vertical plate with suction. Recently, Gireesha et al. [5] investigate effect of chemical reaction on mhd boundary layer flow and melting heat transfer of Williamson nanofluid in porous medium; it was reported that, chemical reaction in the system results in consumption of the chemical and hence, leads to decrease of concentration profile. Many fluids are well known to exhibit non-Newtonian behaviour, in view of this, non-newtonian fluid is defined as a fluid in which shear stress is not directly proportional to deformation rate. An example is paint, it is very thick when stored in the can, but becomes thin when sheared by brushing. It is imperative to note that, non-Newtonian fluids are generally categorized as having time-independent or time-dependent behaviour. It is interesting to note that study of non-Newtonian fluids is further complicated by the fact that the apparent viscosity may be time-dependent. It is a usual occurence that after deformation, some fluids partially return to their original position when the applied stress is released, such fluids are called viscoelastic Pritchard and Leylegian [6] . Some recent studies describing the flows of viscoelastic fluid have been undertaken by Mishra et al. [7] , Olanrewaju et al. [8] Narayana et al. [9] , Das [10] , Choudhury and Das [11] .

Thermophoresis of particles is referred to as a mechanism of movement of small particles in the direction of decreasing thermal gradient. In other words, we can say that thermophoresis makes it possible for the deposition of small particles on the cold surfaces. The effect of this phenomenon was first observed in 1870 by Tyndal [12] , when he observed that a particle free zone around a heated surface appeared in dusty air, and later in 1884, Aitken [13] came up with a prove that the microscopic explanation to the effect was due to the heavier bombardment of the particle from the molecules on the hot region compared cold region. In this phenomenon, the gas molecules migrating from the hot side of the particles have a greater velocity than those migrating from the cold side. The faster moving molecules collide with the particles more forcefully. Therefore, the velocity attained by the particle is referred to as thermophoretic velocity while the force experienced by the suspended particles due to the temperature gradient is referred to as thermophoretic force, and the direction of the force is opposite to the temperature gradient. Stanford Shateyi [14] . Thermophoresis is often experienced in our day-to-day living, and some common experience of this phenomenon is the blackening of glass globe of kerosene lanterns, blackening of the white florescent bulb, chimmeys and industrial furnace walls by carbon particles. The principle of thermophoresis is utilized to manufacture graded index silicon dioxide and germanium oxide optical fiber performs used in the field of communications. Moreso, thermophoresis has many practical applications in aerosol technology, deposition of silicon thin films and radioactive particles in nuclear reactor safety simulations as reported by Hayat and Qasim [15] , Alam et al. [16] .

The objective of this study is to consider the heat and mass of an electrically conducting viscoelastic fluid flow over linearly stretching sheet. The novelty of the study is to investigate the effect of thermophoresis and n^th order of chemical reaction on viscoelastic fluid considering variable thermophysical properties.

2. Mathematical Formulation

The incompressible second-order fluids whose constitutive equation based on the postulate of gradually fading memory was given by Coleman and Noll [17] as

$T=-pl+\mu A_1+{\alpha }_1{A}_2+{\alpha }_2{A}_1^2,$ (1)

where T is the stress tensor, p is the pressure, $\mu \mathrm{,}{\alpha }_1\mathrm{,}{\alpha }_2$ are material constants with ${\alpha }_1\le 0$ , and $A_1$ and $A_2$ are defined as;

$A_1=\left(gradv\right)+{\left(gradv\right)}^{\text{T}}\mathrm{,}$ (2)

$A_2=\frac{\text{d}}{\text{d}t}{A}_1+A_1\cdot gradv+{\left(gradv\right)}^{\text{T}}\cdot A_1.$ (3)

Coleman and Noll [17] showed that the model (1) exhibits normal-stress differences in shear flow and is an approximation to a simple fluid in the sense of retardation.

We consider a steady two-dimensional free convective boundary layer flow of an electrically conducting viscoelastic fluid of variable viscosity and thermal conductivity. Keeping the origin fixed, the sheet is then stretched with a velocity $u_w\left(x\right)$ , varying linearly with the distance from the slit. The flow is assumed to flow in x-direction which is along vertical surface and y-axis is normal to it. Fluid suction/injection is imposed at the plate surface. The temperature and concentration of the surface $T_w$ and $C_w$ is held uniform at which is higher than the ambient temperature $T_{\infty }$ and concentration $C_{\infty }$ i.e. $\left(T_w>T_{\infty }\right)$ and $\left(C_w>C_{\infty }\right)$ . The uniform magnetic field of magnitude $B_o$ is applied normal to the plate. Also the magnetic Reynolds number is assumed to be small so that the induced magnetic field is negligible in comparison to the applied magnetic field. Under the foregoing assumptions with the Boussinesq approximation, the governing equations of the MHD free convection flow are:

Continuity Equation

$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0,$ (4)

Momentum Equation

$\begin{array}{c}u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=\frac{1}{\rho }\frac{\partial }{\partial y}\left(\mu \left(T\right)\frac{\partial u}{\partial y}\right)-k_o\left[\frac{\partial }{\partial x}\left(u\frac{{\partial }^{2}u}{\partial y^2}\right)+\frac{\partial u}{\partial y}\frac{{\partial }^{2}v}{\partial y^2}+v\frac{{\partial }^3}{\partial y^3}\right]\\ -\sigma \frac{B_o^{2}u}{\rho }+g\beta \left(T-T_{\infty }\right)+g{\beta }^{*}\left(C-C_{\infty }\right),\end{array}$ (5)

Energy Equation

$\begin{array}{c}u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}=\frac{1}{\rho C_p}\frac{\partial }{\partial y}\left(\kappa \left(T\right)\frac{\partial T}{\partial y}\right)\\ +\frac{\kappa a}{\rho C_p}\vartheta \left[A\left(T_w-T_{\infty }\right){{\displaystyle \exp }}^{-y\sqrt{\frac{a}{\vartheta }}}+B\left(T-T_{\infty }\right)\right],\end{array}$ (6)

Concentration Equation

$u\frac{\partial C}{\partial x}+v\frac{\partial C}{\partial y}=D\frac{{\partial }^{2}C}{\partial y^2}-\frac{\partial }{\partial y}\left(V_{T}C\right)-k_n{C}^n$ (7)

Subject to boundary conditions

$u=u_w\left(x\right)=ax,v=v_w\left(x\right),T=T_w,C=C_w,\text{at}y=0,$ (8)

$u\to 0,T\to T_{\infty },C\to C_{\infty },\text{as}y\to \infty ,$ (9)

In order to justify the variation in the thermo-physical property of the viscoelastic fluid as it flows past a vertical heated surface, classical Boussinesq’s approximation is adopted such that the temperature at the surface is greater than temperature of the fluid at the free stream. It is valid to consider the mathematical model of temperature dependent viscosity model which was developed using the experimental data of Batchelor [18] together with the mathematical model of temperature dependent thermal conductivity model of Charraudeau [19] as;

$\mu \left(T\right)={\mu }^{*}\left[1+b\left(T_w-T\right)\right]\text{and}\kappa \left(T\right)={\kappa }^{*}\left[1+\delta \left(T-T_{\infty }\right)\right]$ (10)

${\mu }^{\mathrm{<mo>\; *\; </mo>}}$ and ${\kappa }^{\mathrm{<mo>\; *\; </mo>}}$ are the constant value of the coefficient of viscosity and thermal conductivity at the free stream respectively. The thermophoretic velocity $V_T$ in Equation (7) can be written in the form [20] as;

$V_T=-\kappa \vartheta \frac{\nabla T}{T_{ref}}=-\kappa \vartheta \frac{1}{T_{ref}}\frac{\partial T}{\partial y},$ (11)

where $k\vartheta$ represents the thermophoretic diffusivity, and $\kappa$ is the thermophoretic coefficient which ranges in value from 0.2 to 1.2 as indicated by Batchelor and Shen [21] and is defined from the theory of Talbot et al. [20] which is given by;

$\kappa =\frac{2C_s\left({\lambda }_g/{\lambda }_p+C_{t}Kn\right)\left[1+Kn\left(C_1+C_2{\exp }^{-C_3/Kn}\right)\right]}{\left(1+3C_{m}Kn\right)\left(1+2{\lambda }_g/{\lambda }_p+2C_{t}Kn\right)}$ (12)

Here, $C_1\mathrm{,}{C}_2\mathrm{,}{C}_3\mathrm{,}{C}_m\mathrm{,}{C}_s\mathrm{,}{C}_t$ are constants, $\lambda$ and ${\lambda }_p$ are the thermal conductivities of the fluid and diffused particles, respectively and $Kn$ is the Knudsen number. A thermophoretic parameter $\tau$ can be defined (see Mills et al. [22] and Tsai [23] ) as;

$\tau =\frac{\kappa \left(T_w-T_{\infty }\right)}{T_{ref}}$ (13)

where u and v are components of velocity in x and y directions respectively, $u_w\left(x\right)$ is the wall shrinking or stretching velocity, ( $a>0$ ) for stretching, ( $a<0$ ) for shrinking and ( $a=0$ ) for static wall, $v_w\left(x\right)$ is the wall mass flux velocity, ρ is the fluid density, $\mu \left(=\vartheta \rho \right)$ is the dynamic viscosity, $\vartheta$ is the kinematic viscosity, σ is the electrical conductivity, T is the fluid temperature in the boundary layer, $T_{\infty }$ is the free stream temperature, β is the thermal expansion

coefficient, $\alpha =\frac{\kappa }{\rho C_p}$ is the thermal diffusivity, $k_o$ is the non-Newtonian

visco-elastic parameter, β is the volumetric coefficient of thermal expansion, ${\beta }^{\mathrm{<mo>\; *\; </mo>}}$ is the volumetric concentration coefficient, A and B are the coefficient of exponentially decaying space and temperature dependent heat source/sink, respectively, D is the mass diffusivity and $k_n$ is the chemical reaction parameter.

The continuity Equation (1) is satisfied by introducing a stream function $\psi$ such that

$u=\frac{\partial \psi }{\partial y},v=-\frac{\partial \psi }{\partial x}.$ (14)

Then, Equations (5), (6), (7), (8) and (9) becomes;

$\begin{array}{l}\frac{\partial \psi }{\partial y}\frac{\partial }{\partial x}\frac{\partial \psi }{\partial y}-\frac{\partial \psi }{\partial x}\frac{\partial }{\partial y}\frac{\partial \psi }{\partial y}\\ =\frac{1}{\rho }\frac{\partial }{\partial y}\left(\mu \left(T\right)\frac{\partial }{\partial y}\left(\frac{\partial \psi }{\partial y}\right)\right)\\ -k_o\left[\frac{\partial }{\partial x}\left(\frac{\partial \psi }{\partial y}\frac{{\partial }^2}{\partial y^2}\frac{\partial \psi }{\partial y}\right)-\frac{\partial }{\partial y}\frac{\partial \psi }{\partial y}\frac{{\partial }^2}{\partial y^2}\frac{\partial \psi }{\partial x}-\frac{\partial \psi }{\partial x}\frac{{\partial }^3}{\partial y^3}\frac{\partial \psi }{\partial y}\right]\end{array}$

$-\frac{\sigma B_o^2}{\rho }\frac{\partial \psi }{\partial y}+g\beta \left(T-T_{\infty }\right)+g{\beta }^{*}\left(C-C_{\infty }\right),$ (15)

$\begin{array}{l}\frac{\partial \psi }{\partial y}\frac{\partial T}{\partial y}-\frac{\partial \psi }{\partial x}\frac{\partial T}{\partial y}\\ =\frac{1}{\rho C_p}\frac{\partial }{\partial y}\left(\kappa \left(T\right)\frac{\partial T}{\partial y}\right)+\frac{\kappa a}{\rho C_p\vartheta }\left[A\left(T_w-T_{\infty }\right){{\displaystyle \exp }}^{-y\sqrt{\frac{a}{\vartheta }}}+B\left(T-T_{\infty }\right)\right],\end{array}$ (16)

$\frac{\partial \psi }{\partial y}\frac{\partial C}{\partial x}-\frac{\partial \psi }{\partial x}\frac{\partial C}{\partial y}=D\frac{{\partial }^{2}C}{\partial y^2}-\frac{\partial }{\partial y}\left(V_{T}C\right)-k_n{C}^n$ (17)

subject to

$\frac{\partial \psi }{\partial y}=u_w\left(x\right)=ax,-\frac{\partial \psi }{\partial x}=v_w\left(x\right),T=T_w,C=C_w,\text{at}y=0,$ (18)

$\frac{\partial \psi }{\partial y}\to \mathrm{0,}T\to T_{\infty }\mathrm{,}C\to C_{\infty }\mathrm{,}\text{as}y\to \infty \mathrm{,}$ (19)

The momentum, energy, and concentration equations can be transformed into the corresponding ordinary differential equations by the following transformation

$\eta =y{\left(\frac{a}{\vartheta }\right)}^{1/2},\psi \left(x,y\right)=x{\left(a\vartheta \right)}^{1/2}f\left(\eta \right),\theta \left(\eta \right)=\frac{T-T_{\infty }}{T_w-T_{\infty }},\varphi \left(\eta \right)=\frac{C-C_{\infty }}{C_w-C_{\infty }}$ (20)

where $\eta$ is the independent dimensionless similarity variable. Thus u and v are given by $u=ax{f}^{\prime }\left(\eta \right)$ , $v=-\sqrt{a\vartheta }f\left(\eta \right)$ , substituting variables (20) into Equations (15)-(19), we obtain the following ordinary differential equations:

$\begin{array}{l}\left[1+\xi -\theta \xi \right]\frac{{\text{d}}^{3}f}{\text{d}{\eta }^3}-{\left(\frac{\text{d}f}{\text{d}\eta }\right)}^2+f\frac{{\text{d}}^{2}f}{\text{d}{\eta }^2}-\xi \frac{\text{d}\theta }{\text{d}\eta }\frac{{\text{d}}^{2}f}{\text{d}{\eta }^2}\\ -R_c\left[2\frac{\text{d}f}{\text{d}\eta }\frac{{\text{d}}^{3}f}{\text{d}{\eta }^3}-{\left(\frac{{\text{d}}^{2}f}{\text{d}{\eta }^2}\right)}^2-f\frac{{\text{d}}^{4}f}{\text{d}{\eta }^4}\right]-M\frac{\text{d}f}{\text{d}\eta }+G_r\theta +G_c\varphi =0\end{array}$ (21)

$\left[1+\theta \epsilon \right]\frac{{\text{d}}^2\theta }{\text{d}{\eta }^2}-P_r{S}_t\frac{\text{d}f}{\text{d}\eta }-P_r\theta \frac{\text{d}f}{\text{d}\eta }+P_{r}f\frac{\text{d}\theta }{\text{d}\eta }+\epsilon {\left(\frac{\text{d}\theta }{\text{d}\eta }\right)}^2+\left(A{\text{e}}^{-\eta }+B\theta \right)=0$ (22)

$\frac{{\text{d}}^2\varphi }{\text{d}{\eta }^2}+S_{c}f\frac{\text{d}\varphi }{\text{d}\eta }-S_c\tau \left(\frac{\text{d}\theta }{\text{d}\eta }\frac{\text{d}\varphi }{\text{d}\eta }+\varphi \frac{{\text{d}}^2\theta }{\text{d}{\eta }^2}\right)-S_c\gamma {\varphi }^n=0$ (23)

The corresponding boundary conditions take the form;

$f\left(0\right)=s,f^{\prime }\left(0\right)=1,\theta \left(0\right)=1,\varphi \left(0\right)=1\text{at}\eta =0,$ (24)

$f^{\prime }\left(\eta \right)\to \mathrm{0,}\theta \left(\eta \right)\to \mathrm{0,}\varphi \left(\eta \right)\to 0\text{as}\eta \to \infty \mathrm{.}$ (25)

In the above equations, primes denote differentiation with respect to $\eta$ . The dimensionless velocity, temperature and concentration are represented as $f\left(\eta \right)$ ,

$\theta \left(\eta \right)$ and $\varphi \left(\eta \right)$ respectively, $P_r=\frac{\vartheta }{\alpha }$ is the Prandtl number, $M=\sigma \frac{B_o^2}{\rho a}$ is the magnetic parameter, $G_r=\frac{g\beta \left(T_w-T_{\infty }\right)}{a^{2}x}$ is the Modified Thermal Grashof number, $G_c=\frac{g{\beta }^{*}\left(C_w-C_{\infty }\right)}{a^{2}x}$ is the Modified Solutal Grashof number, $R_c=\frac{k_{o}a}{\vartheta }$ is the viscoelastic parameter, $\gamma =\frac{k_n}{a}$ is the chemical reaction parameter, $S_c=\frac{\vartheta }{D}$ is the Schimdt number.

It is worth mentioning here that the chemical reaction parameter $\gamma$ is a real number ( $\gamma <0$ indicates the generative chemical reaction, $\gamma >0$ denotes the destructive chemical reaction, and $\gamma =0$ for the non-reactive species).

The physical quantities of interest are the skin friction coefficient $C_f$ , the local Nusselt number $N{u}_x$ and the local Sherwood number Sh are defined as;

$C_f=\frac{{\tau }_w}{\rho U^2/2},N{u}_x=\frac{x{q}_w}{\kappa \left(T_w-T_{\infty }\right)}.$ (26)

where the shear stress at the all surface is expressed as $\tau$

${\tau }_w=\mu {\left[\frac{\partial u}{\partial y}\right]}_{y=0}\text{or}{\tau }_w=C_f\frac{\rho U^2}{2}$ (27)

$\tau =\mu {\left[ax\sqrt{\frac{a}{x}}{f}^{″}\left(0\right)\right]}_{\eta =0}$ (28)

2.1. Optimal Homotopy Analysis Solutions

In many cases, by means of analyzing the physical background and the initial/boundary conditions of the nonlinear differential problem, we might know what kinds of base functions are proper to represent the solution, even without solving the given nonlinear problem. In view of the boundary conditions (24) and (25), $f\left(\eta \right)$ , $p\left(\eta \right)$ and $\theta \left(\eta \right)$ can be expressed by the set of base functions in the form

$\langle {\eta }^j\exp \left(-nj\right)|j\ge 0,n\ge 0\rangle$ (29)

The solutions $f\left(\eta \right)$ , $p\left(\eta \right)$ and $\theta \left(\eta \right)$ can be represented in a series form as

$f\left(\eta \right)=a_{0,0}^0+{\displaystyle \underset{n=0}{\overset{\infty }{\sum }}}{\displaystyle \underset{k=0}{\overset{\infty }{\sum }}}{a}_{n,k}^k{\eta }^k\exp \left(-nj\right)$ (30)

$\theta \left(\eta \right)={\displaystyle \underset{n=0}{\overset{\infty }{\sum }}}{\displaystyle \underset{k=0}{\overset{\infty }{\sum }}}{b}_{n,k}^k{\eta }^k\exp \left(-nj\right)$ (31)

$\varphi \left(\eta \right)={\displaystyle \underset{n=0}{\overset{\infty }{\sum }}}{\displaystyle \underset{k=0}{\overset{\infty }{\sum }}}{c}_{n,k}^k{\eta }^k\exp \left(-nj\right)$ (32)

In which $a_{n\mathrm{,}k}^k$ , $b_{n\mathrm{,}k}^k$ and $c_{n\mathrm{,}k}^k$ are the coefficients. As long as such a set of base functions are determined, the auxiliary function $H\left(\eta \right)$ , the initial approximation $f_o\left(\eta \right)$ , ${\theta }_o\left(\eta \right)$ and ${\varphi }_o\left(\eta \right)$ , and the auxiliary linear operator $L_f$ , $L_{\theta }$ and $L_{\varphi }$ must be chosen properly in such a way that all solutions of the corresponding high-order deformation of Equations (52)-(54) exist and can be expressed by this set of base functions. Invoking the rule of solution expressions above for $f\left(\eta \right)\mathrm{,}\theta \left(\eta \right)$ and $\varphi \left(\eta \right)$ on (21)-(23) together with boundary conditions (24) and (25), the initial guesses $f_o\left(\eta \right)$ , ${\theta }_o\left(\eta \right)$ and ${\varphi }_o\left(\eta \right)$ which satisfies both the initial and the boundary conditions (24) and (25) are;

$f_o\left(\eta \right)=1+S-\exp \left(-\eta \right),{\theta }_o\left(\eta \right)=\exp \left(-\eta \right),{\varphi }_o\left(\eta \right)=\exp \left(-\eta \right)$ (33)

Linear operators $L_f$ , $L_p$ and $L_{\theta }$ are

$L_f\left[f\left(\eta ;q\right)\right]=\frac{{\partial }^{3}f\left(\eta ;q\right)}{\partial {\eta }^3}-\frac{\partial f\left(\eta ;q\right)}{\partial \eta }$ (34)

$L_{\theta }\left[\theta \left(\eta \mathrm{<mo>\; ;\; </mo>}q\right)\right]=\frac{{\partial }^2\theta \left(\eta \mathrm{<mo>\; ;\; </mo>}q\right)}{\partial {\eta }^2}-\theta \left(\eta \mathrm{<mo>\; ;\; </mo>}q\right)$ (35)

$L_{\varphi }\left[\varphi \left(\eta \mathrm{<mo>\; ;\; </mo>}q\right)\right]=\frac{{\partial }^2\varphi \left(\eta \mathrm{<mo>\; ;\; </mo>}q\right)}{\partial {\eta }^2}-\varphi \left(\eta \mathrm{<mo>\; ;\; </mo>}q\right)$ (36)

The operators $L_f$ , $L_{\theta }$ and $L_{\varphi }$ have the following properties

$\begin{array}{l}{L}_f\left[C_1+C_2\exp \left(-\eta \right)+C_3\exp \left(-\eta \right)\right]=0,\\ L_{\theta }\left[C_4\exp \left(-\eta \right)+C_5\right]=0,\\ L_{\varphi }\left[C_6\exp \left(-\eta \right)+C_7\right]=0\end{array}$ (37)

In which $C_1,C_2,C_3,C_4,C_5,C_6$ and $C_7$ are constants.

$\left(1-q\right)L_f\left[f\left(\eta ;q\right)-f_o\left(\eta \right)\right]=q{\hslash }_f{H}_f\left(\eta \right)N\left[f\left(\eta ;q\right),\theta \left(\eta ;q\right),\varphi \left(\eta ;q\right)\right]$ (38)

$\left(1-q\right)L_{\theta }\left[\theta \left(\eta \mathrm{<mo>\; ;\; </mo>}q\right)-{\theta }_o\left(\eta \right)\right]=q{\hslash }_{\theta }{H}_{\theta }\left(\eta \right)N\left[f\left(\eta ;q\right),\theta \left(\eta ;q\right),\varphi \left(\eta ;q\right)\right]$ (39)

$\left(1-q\right)L_{\varphi }\left[\varphi \left(\eta \mathrm{<mo>\; ;\; </mo>}q\right)-{\varphi }_o\left(\eta \right)\right]=q{\hslash }_{\theta }{H}_{\theta }\left(\eta \right)N\left[f\left(\eta ;q\right),\theta \left(\eta ;q\right),\varphi \left(\eta ;q\right)\right]$ (40)

Subject to boundary conditions

$f\left(\eta =0;q\right)=s,\frac{\partial f\left(\eta =0;q\right)}{\partial \eta }=1,\theta \left(\eta =0;q\right)=1,\varphi \left(\eta =0;q\right)=1$ (41)

$\frac{\partial f\left(\eta \to \infty ;q\right)}{\partial \eta }\to 0,\theta \left(\eta \to \infty \right)\to 0,\varphi \left(\eta \to \infty ;q\right)=0$ (42)

where q is embedding parameters and the nonlinear operators are defined as

$\begin{array}{l}\left(1+\xi -\theta \left(\eta \mathrm{<mo>\; ;\; </mo>}q\right)\xi \right)\frac{{\partial }^{3}f\left(\eta \mathrm{<mo>\; ;\; </mo>}q\right)}{\partial {\eta }^3}-\frac{\partial f\left(\eta \mathrm{<mo>\; ;\; </mo>}q\right)}{\partial \eta }\frac{\partial f\left(\eta \mathrm{<mo>\; ;\; </mo>}q\right)}{\partial \eta }\\ +f\left(\eta \mathrm{<mo>\; ;\; </mo>}q\right)\frac{{\partial }^{2}f\left(\eta \mathrm{<mo>\; ;\; </mo>}q\right)}{\partial {\eta }^2}-\xi \frac{\partial \theta \left(\eta \mathrm{<mo>\; ;\; </mo>}q\right)}{\partial \eta }\frac{{\partial }^{2}f\left(\eta \mathrm{<mo>\; ;\; </mo>}q\right)}{\partial {\eta }^2}\\ -Rc(2\frac{\partial f\left(\eta \mathrm{<mo>\; ;\; </mo>}q\right)}{\partial \eta }\frac{{\partial }^{3}f\left(\eta \mathrm{<mo>\; ;\; </mo>}q\right)}{{\partial }^3\eta }-\frac{{\partial }^{2}f\left(\eta \mathrm{<mo>\; ;\; </mo>}q\right)}{\partial {\eta }^2}\frac{{\partial }^{2}f\left(\eta \mathrm{<mo>\; ;\; </mo>}q\right)}{\partial {\eta }^2}\\ -f\left(\eta \mathrm{<mo>\; ;\; </mo>}q\right)\frac{{\partial }^{4}f\left(\eta \mathrm{<mo>\; ;\; </mo>}q\right)}{\partial {\eta }^4})-M\frac{\partial f\mathrm{<mo\; stretchy="false">\; (\; </mo>}\eta \mathrm{<mo>\; ;\; </mo>}q\mathrm{<mo\; stretchy="false">\; )\; </mo>}}{\partial \eta }+G_r\theta \left(\eta \mathrm{<mo>\; ;\; </mo>}q\right)+G_c\varphi \left(\eta \mathrm{<mo>\; ;\; </mo>}q\right)\end{array}$ (43)

$\begin{array}{l}\left(1+\theta \left(\eta ;q\right)\epsilon \right)\frac{{\partial }^2\theta \left(\eta ;q\right)}{\partial {\eta }^2}+P_{r}f\left(\eta ;q\right)\frac{\partial \theta \left(\eta ;q\right)}{\partial \eta }\\ +\epsilon \frac{\partial \theta \left(\eta ;q\right)}{\partial \eta }\frac{\partial \theta \left(\eta ;q\right)}{\partial \eta }+\left(A{\text{e}}^{-\eta }+B\theta \left(\eta ;q\right)\right)=0\end{array}$ (44)

$\begin{array}{l}\frac{{\partial }^2\varphi \left(\eta ;q\right)}{\partial {\eta }^2}+Scf\left(\eta ;q\right)\frac{\partial \varphi \left(\eta ;q\right)}{\partial \eta }\\ -S_c\tau \left(\frac{\partial \theta \left(\eta ;q\right)}{\partial \eta }\frac{\partial \varphi \left(\eta ;q\right)}{\partial \eta }+\varphi \left(\eta ;q\right)\frac{{\partial }^2\theta \left(\eta ;q\right)}{{\partial }^2}\right)-S_c\gamma \varphi {\left(\eta ;q\right)}^n=0\end{array}$ (45)

Obviously, when $q=0$ and $q=1$ , zero order of deformation equations (38) to (40) leads to With the property

$f\left(\eta ;0\right)=f_o\left(\eta \right),f\left(\eta ;1\right)=f_o\left(\eta \right)$ (46)

$\theta \left(\eta ;0\right)={\theta }_o\left(\eta \right),\theta \left(\eta ;1\right)={\theta }_o\left(\eta \right)$ (47)

$\varphi \left(\eta ;0\right)={\varphi }_o\left(\eta \right),\varphi \left(\eta ;1\right)={\varphi }_o\left(\eta \right)$ (48)

Expanding $f\left(\eta ;q\right)$ , $\theta \left(\eta ;q\right)$ and $\varphi \left(\eta ;q\right)$ in Taylor series with respect to the embedding parameter q,

$f\left(\eta ;q\right)=f_o\left(\eta \right)+{\displaystyle \underset{m=1}{\overset{\infty }{\sum }}}{f}_m\left(\eta \right)q^m\text{where}{f}_m\left(\eta \right)=\frac{1}{m!}{\frac{{\partial }^{m}f\left(\eta ;q\right)}{\partial q^m}|}_{q=0}$ (49)

$\theta \left(\eta ;q\right)={\theta }_o\left(\eta \right)+{\displaystyle \underset{m=1}{\overset{\infty }{\sum }}}{\theta }_m\left(\eta \right)q^m\text{where}{\theta }_m\left(\eta \right)=\frac{1}{m!}{\frac{{\partial }^m\theta \left(\eta ;q\right)}{\partial {\eta }^m}|}_{q=0}$ (50)

$\varphi \left(\eta ;q\right)={\varphi }_o\left(\eta \right)+{\displaystyle \underset{m=1}{\overset{\infty }{\sum }}}{\varphi }_m\left(\eta \right)q^m\text{where}{\varphi }_m\left(\eta \right)=\frac{1}{m!}{\frac{{\partial }^m\varphi \left(\eta ;q\right)}{\partial {\eta }^m}|}_{q=0}$ (51)

The auxiliary parameters are so properly chosen that the series (38)-(40) converge at $q=1$ . Hence,

$f\left(\eta ;q\right)=f_o\left(\eta \right)+{\displaystyle \underset{m=1}{\overset{\infty }{\sum }}}{f}_m\left(\eta \right)q^m$ (52)

$\theta \left(\eta ;q\right)={\theta }_o\left(\eta \right)+{\displaystyle \underset{m=1}{\overset{\infty }{\sum }}}{\theta }_m\left(\eta \right)q^m$ (53)

$\varphi \left(\eta ;q\right)={\varphi }_o\left(\eta \right)+{\displaystyle \underset{m=1}{\overset{\infty }{\sum }}}{\theta }_m\left(\eta \right)q^m$ (54)

For the mth order deformation, differentiate (38) to (40) m times with respect to q, divide by $m!$ and set $q=0$ , then we have;

$L_f\left[f_m\left(\eta \right)-{\chi }_m{f}_{m-1}\left(\eta \right)\right]={\hslash }_f{H}_f\left(\eta \right)R_m^f\left(\eta \right)$ (55)

$L_{\theta }\left[{\theta }_m\left(\eta \right)-{\chi }_m{\theta }_{m-1}\left(\eta \right)\right]={\hslash }_{\theta }{H}_{\theta }\left(\eta \right)R_m^{\theta }\left(\eta \right)$ (56)

$L_{\varphi }\left[{\varphi }_m\left(\eta \right)-{\chi }_m{\varphi }_{m-1}\left(\eta \right)\right]={\hslash }_{\theta }{H}_{\varphi }\left(\eta \right)R_m^{\varphi }\left(\eta \right)$ (57)

Subject to

$f_m\left(\eta =0;0\right)=0,\frac{\partial f\left(\eta =0;0\right)}{\partial \eta }=0,{\theta }_m\left(\eta =0\right)=0,{\varphi }_m\left(\eta =0\right)=0$ (58)

where

$\begin{array}{c}{R}_m^f\left(\eta \right)=\left(1+Xi-{\theta }_{m-1}\right)\frac{{\text{d}}^3{f}_{m-1}}{\text{d}{\eta }^3}-{\displaystyle \underset{k=0}{\overset{m-1}{\sum }}}\frac{\text{d}{f}_{m-1-k}}{\text{d}\eta }\frac{\text{d}{f}_k}{\text{d}\eta }\\ +{\displaystyle \underset{k=0}{\overset{m-1}{\sum }}}{f}_k\frac{{\text{d}}^2{f}_{m-1-k}}{\text{d}{\eta }^2}+\xi {\displaystyle \underset{k=0}{\overset{m-1}{\sum }}}\frac{\text{d}{\theta }_{m-1}}{\text{d}\eta }\frac{{\text{d}}^2{f}_{m-1}}{\text{d}{\eta }^2}\\ -R_c{\displaystyle \underset{k=0}{\sum }}m-1\left(2\frac{\text{d}{f}_{m-1-k}}{\text{d}\eta }\frac{{\text{d}}^3{f}_k}{\text{d}{\eta }^3}-\frac{\text{d}{f}_{m-1-k}}{\text{d}\eta }\frac{\text{d}{f}_k}{\text{d}\eta }-f_{m-1-k}\frac{{\text{d}}^4{f}_k}{\text{d}{\eta }^4}\right)\\ -M\frac{\text{d}{f}_{m-1}}{\text{d}\eta }+G_r{\theta }_{m-1}+G_c{\varphi }_{m-1}\end{array}$ (59)

$\begin{array}{c}{R}_m^{\theta }\left(\eta \right)=\left(1+{\theta }_{m-1}\epsilon \right)\frac{{\text{d}}^2{\theta }_{m-1}}{\text{d}{\eta }^2}+P_r{\displaystyle \underset{k=0}{\overset{m-1}{\sum }}}{f}_{m-1-k}\frac{\text{d}{\theta }_k}{\text{d}\eta }\\ +\epsilon {\displaystyle \underset{k=0}{\overset{m-1}{\sum }}}\left(\frac{\text{d}{\theta }_{m-1-k}}{\text{d}\eta }\frac{\text{d}{\theta }_k}{\text{d}\eta }\right)+\left(A{\text{e}}^{-y\sqrt[a]{\vartheta }}+B{\theta }_{m-1}\right)\\ =0\end{array}$ (60)

$\begin{array}{c}{R}_m^{\varphi }\left(\eta \right)=\frac{{\text{d}}^2{\varphi }_{m-1}}{\text{d}{\eta }^2}+S_c{\displaystyle \underset{k=0}{\overset{m-1}{\sum }}}\left(f_{m-1-k}\frac{{\varphi }_k}{d\eta }\right)\\ -S_c\tau {\displaystyle \underset{k=0}{\overset{m-1}{\sum }}}\left(\frac{\text{d}{\theta }_{m-1-k}}{\text{d}\eta }\frac{\text{d}{\varphi }_k}{\text{d}\eta }+{\varphi }_{m-1-k}\frac{{\text{d}}^2{\theta }_k}{\text{d}\eta }\right)-S_c\gamma {\varphi }_{m-1}^n\end{array}$ (61)

And

${\chi }_m=0\text{when}m\le 1$

${\chi }_m=1\text{when}m>1$

The general solutions of equations are given by

$f_m\left(\eta \right)=f_m^{*}+C_1+C_2\exp \left(\eta \right)+C_3\exp (\; -\; \eta \; )$

${\theta }_m\left(\eta \right)={\theta }_m^{*}+C_4+C_5\exp \left(\eta \right)+C_6\exp (\; -\; \eta \; )$

${\varphi }_m\left(\eta \right)={\varphi }_m^{*}\left(\eta \right)+C_7\exp \left(\eta \right)+C_8\exp (\; -\; \eta \; )$

Here, $f_m^{*}\left(\eta \right)$ , ${\theta }_m^{*}\left(\eta \right)$ and ${\varphi }_m^{*}\left(\eta \right)$ are the particular solutions of Equations (52)-(54). Following the rule of solution expression, the rule of coefficient ergodicity and the rule of solution existence as discussed in [24] - [27] we choose auxiliary functions as

$H_f=H_p=H_{\theta }=1$ (62)

2.2. Convergence of the Optimal Homotopy Solutions

It is obvious that the series (64)-(66) consist of the non-zero auxiliary parameters ${\hslash }_f$ , ${\hslash }_{\theta }$ and ${\hslash }_{\varphi }$ which can adjust and control the convergence. The interval on h-axis for which the h-curve becomes parallel to the h-axis is recognized as the set of admissible values of ${\hslash }_f$ , ${\hslash }_{\theta }$ and ${\hslash }_{\varphi }$ for which the solution series converges. These figures show that the ranges for the acceptable values of ${\hslash }_f$ , ${\hslash }_{\theta }$ and ${\hslash }_{\varphi }$ are $-1.7\le {\hslash }_f\le -0.3$ , $-1.2\le {\hslash }_{\theta }\le -0.38$ and $-1.6\le {\hslash }_{\varphi }\le -0.3$ . Obviously, from the h-curves for this problem, we obtained the approximate optimal values of ${\hslash }_f$ , ${\hslash }_{\theta }$ and ${\hslash }_{\varphi }$ at 10th-order of approximation as −1.20967, −0.999414 and −1.46231.

3. Results and Discussion

In order to gain an insight into the behavior of the fluid as it flows, analytic approximate solution of the dimensionless governing equation described in the previous section has been carried out using various values of elastic parameter $R_c$ , heat source parameter A and B, Prandtl number $P_r$ , Magnetic parameter M, thermophoresis parameter $\tau$ , temperature-dependent variable viscosity and thermal conductivity parameters ξ and ε when ( $G_r=G_c<0$ ) which physically means cooling of the fluid or heating of the surface (opposing flow) and when ( $G_r=G_c>0$ ) means heating of the fluid or cooling of the surface (assisting flow). Table 1, Figure 1(a) and Figure 1(b) reveal the influence of temperature dependent viscous and thermal conductivity parameters on velocity and temperature profiles respectively when ( $G_r=G_c<0$ ). It is observed from Figure 1(a) that as ξ and ε increase from 1.0 to 4.0 the velocity distribution increases within the range of $0.5\le \eta \le 5.6$ . At exact value of $\eta =5.8$ all profile converges quickly towards the freestream, and likewise it is noticed from Figure Figure 1(b) that the temperature profile is increased as ξ and ε increases.

Figure 2(a) and Figure 2(b) present the effect of ξ and ε on velocity and temperature profiles when ( $G_r=G_c>0$ ). It is noticed that as values of ξ and ε increase, there is a slight decrease in the magnitude of the velocity profile within the region of $0\le \eta \le 1.2$ , making all the curves to merge together and thereafter at exact value of $\eta =1.3$ there is a conspicuous increase in the strength of the velocity and tends to satisfy the boundary condition. Likewise the temperature profile increases with the increase in the values of ξ and ε. Hence, the boundary layer thickness increases.

The effect of thermophoresis parameter $\tau$ on concentration profiles when $\gamma >0$ for different orders of chemical reaction are presented in Figures 3(a) and Figures 3(b). It can been seen from the plots that for lesser values of Schmidth number $S_c=0.22$ , the magnitude of the concentration profiles decreases and thus, the concentration layer becomes thinner for both cases of $n=1$ and $n=2$ . Physically, thermophoresis is best explained as the migration

Table 1. Values of $-{\theta }^{\prime }\left(0\right)$ for various values $R_c$ when $M=0.5$ , $G_r=G_c=0.5$ , $A=0.4$ , $B=0.3$ , $P_r=0.71$ , $S_c=0.22$ , $\tau =0.4$ , $\gamma =0.5$ , $s=0.3$ , $\xi =0.3$ , $\xi =0.3$ , $s=0.3$ , $n=1$ .

(a) (b)

Figure 1. (a) Effect of viscous and thermal conductivity parameters ξ and ε on velocity profile when $G_r=G_c<0;$ (b) Effect of viscous and thermal conductivity parameters ξ and ε on temperature profile when $G_r=G_c<0$ .

(a) (b)

Figure 2. (a) Effect of viscous and thermal conductivity parameters ξ and ε on velocity profile when $G_r=G_c>0;$ (b) Effect of viscous and thermal conductivity parameters ξ and ε on temperature profile when $G_r=G_c<0$ .

of small sized particles in the direction of decreasing thermal gradient and this can be traced to the fact that, when there is much heat energy supplied in the boundary layer, the small particles suspended in the fluid tends to move to a

(a) (b)

Figure 3. (a) Effect of thermophoresis parameter τ on concentration profile when $n=1;$ (b) Effect of thermophoresis parameter τ on concentration profile when $n=2$ .

region where there is little or low heat energy and thus leading to reduction in the concentration of the species. It is further observed that the concentration of the species quickly reduces faster towards the freestream and asmpotically satisfies the boundary condition at $n=1$ , compared with when the order of the chemical reaction increases to $n=2$ .

Figure 4(a) and Figure 4(b) clearly reveal that increasing Prandtl number $P_r$ causes decrease in velocity, temperature profiles and the associated boundary layer thickness. This can be attributed to the fact that, high $P_r$ for fluid implies low thermal conductivity and therefore the fluid with high $P_r$ attains lower temperature and so the temperature distribution decreases. The effect of space dependent and temperature dependent parameters A and B are shown in Figure 5(a) and Figure 5(b). It is observed that as A and B increase the velocity and temperature profiles increase respectively. This effect is significant because when A and B have increasing positive values, more heat will be produced within the fluid boundary and thus influence the thermal boundary layer (Figure 6, Figure 7).

4.Conclusions

The study considers heat and mass transfer of viscoelastic fluid along vertical surface with thermophoresis and variable fluid properties. Series solutions for velocity, temperature and concentration fields are developed and discussed. We have investigated the effects of various governing parameters. Some key observations are mentioned below

1) The thickness of the velocity and thermal boundary layer increases with an increase in ξ and ε when ( $G_r=G_c<0$ ).

2) The velocity and temperature distributions are increasing functions of ξ and ε when ( $G_r=G_c>0$ ).

3) Concentration boundary layer thickness decreases with an increase in τ when $\gamma >0$ for both cases of $n=1,2$ .

(a) (b)

Figure 4. (a) Effect of Prandtl number $P_r$ on velocity profile; (b) Effect of Prandtl number $P_r$ on temperature profile.

(a) (b)

Figure 5. (a) Effect of heat source parameter A & B on velocity profile; (b) Effect of heat source parameter A & B on temperature profile.

(a) (b)

Figure 6. (a) The h-curve of $f^{″}\left(0\right)$ obtained at 10th-order of approximation; (b) The h-curve of ${\theta }^{\prime }\left(0\right)$ obtained at 10th-order of approximation.

Figure 7. The h-curve of ${\varphi }^{\prime }\left(0\right)$ obtained at 10th-order of approximation.

4) The velocity and temperature distribution of the transport phenomenon are decreasing properties of the flow.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Prandtl, L. (1905) Uber Flussigkeigsbewegung bei sehr kleiner Reibung. In: Verhandlungen des dritten internationalen Mathematischen Kongresses, Heidelberg, Teubner Verlag, Leipzig, 484-491.
[2]	Cengel, Y.A. and Ghajar, A.J. (2015) Heat and Mass Transfer: Fundamentals and Applications. McGraw-Hill Education, New York.
[3]	Gangadhar, K. and Reddy, N.B. (2013) Chemically Reacting MHD Boundary Layer Flow of Heat and Mass Transfer over a Moving Vertical Plate in a Porous Medium with Suction. Journal of Applied Fluid Mechanics, 6, 107-114.
[4]	Ibrahim, S.Y. and Makinde, O.D. (2010) Chemically Reacting MHD Boundary Layer Flow of Heat and Mass Transfer over a Moving Vertical Plate with Suction. Scientific Research and Essays, 5, 2875-2882.
[5]	Krishnamurthy, M.R., Prasannakumara, B.C. and Gireesha, B.J. (2016) Effect of Chemical Reaction on MHD Boundary Layer Flow and Melting Heat Transfer of Williamson Nanofluid in Porous Medium. Engineering Science and Technology, an International Journal, 19, 53-61. https://doi.org/10.1016/j.jestch.2015.06.010
[6]	Pritchard, P.J. and Leylegian, J.C. (2011) Fox and McDonald's Introduction to Fluid Mechanics. 8th Edition, John Wiley and Sons, Inc., Hoboken.
[7]	Bhukta, D., Dash, G.C. and Mishra, S.R. (2014) Heat and Mass Transfer on MHD Flow of a Viscoelastic Flow through Porous Media over a Shrinking Sheet. International Scholarly Research Notices, 11 p. https://doi.org/10.1155/2014/572162
[8]	Gbadeyan, J.A., Idowu, A.S., Ogunsola, A.W., Agboola, O.O. and Olanrewaju, P.O. (2011) Heat and Mass Transfer for Soret and Dufour's Effect on Mixed Convectioin Boundary Layer Flow over a Stretching Vertical Surface in a Porous Medium Filled with a Visco-elastic Fluid in the Presence of Magnetic Field. Global Journal of Science Frontier Research, 11, 97-114.
[9]	Venkateswarlu, B. and Narayana, P.V.S. (2015) Mhd Visco-Elastic Fluid Flow over a Continuously Moving Vertical Surface with Chemical Reaction. Walailak Journal of Science and Technology, 12, 775-783.
[10]	Das, U.J. (2014) Soret and Dufour Effects on Steady Free Convective Mhd Viscoelastic Fluid Flow Confined between a Long Vertical Wavy Wall and Parallet Flat Wall. Thammasat International Journal of Science and Technology, 19, 9-21.
[11]	Choudhury, R. and Das, S.K. (2014) Visco-Elastic Mhd Free Convective Flow through Porous Media in Presence of Radiation and Chemical Reaction with Heat and Mass Transfer. Journal of Applied Fluid Mechanics, 7, 603-609.
[12]	Tyndal, J. (1870) On Dust and Disease. Proceedings of the Royal Institution, 6, 1-14.
[13]	Aitken, J. (1884) On the Formation of Small Clear Spaces in Dusty Air. Transactions of the Royal Society of Edinburgh, 32, 239-272.
[14]	Standford, S. (2013) A New Numerical Approach to MHD Flow of a Maxwell Fluid Past a Vertical Stretching Sheet in the Presence of Thermophoresis and Chemical Reaction. Boundary Value Problems, 2013, 196.
[15]	Hayat, T. and Qasim, M. (2010) Influence of Thermal Radiation and Joule Heating on MHD Flow of a Maxwell Fluid in the Presence of Thermophoresis. International Journal of Heat and Mass Transfer, 53, 4780-4788. https://doi.org/10.1016/j.ijheatmasstransfer.2010.06.014
[16]	Alam, M.S., Rahman, M.M. and Sattar, M.A. (2007) Similarity Solutions for Hydromagnetic Free Convective Heat and Mass Transfer Flow along a Semi-Infinite Permeable Inclined Flat Plate with Heat Generation and Thermophoresis. Nonlinear Analysis: Modelliing and Control, 12, 433-445.
[17]	Coleman, B.D. and Noll, W. (1990) An Approximation Theorem for Functionals with Application in Continuum Mechanics. Archive for Rational Mechanics and Analysis, 6, 355.
[18]	Batchelor, G.K. (1987) An Introduction to Fluid Dynamics. Cambridge University Press, London.
[19]	Charraudeau, J. (1975) Influence de gradents de properties physiques en convection force application au cas du tube. International Journal of Heat and Mass Transfer, 18, 87-95. https://doi.org/10.1016/0017-9310(75)90011-3
[20]	Talbot, L., Cheng, R.K. and Scheffer, R.W. (1980) Thermophoresis of Particles in a Heated Boundary Layer. Journal of Fluid Mechanics, 101, 737-758. https://doi.org/10.1017/S0022112080001905
[21]	Batchelor, G.K. and Shen, C. (1985) Thermophoretic Deposition of Particles in Gas Flowing over Cold Surface. Journal of Colloid and Interface Science, 107, 21-37.
[22]	Mills, A.F., Hang, X. and Ayazi, F. (1984) The Effect of Wall Suction and Thermophoresis on Aerosol-Particle Deposition from a Laminar Boundary Layer on a Flat Plate. International Journal of Heat and Mass Transfer, 27, 1110-1114. https://doi.org/10.1016/0017-9310(84)90127-3
[23]	Tsai, R. (1999) A Simple Approach for Evaluating the Effect of Wall Suction and Thermophoresis on Aerosol Particle Deposition from a Laminar Flow over a Flat Plate. International Communications in Heat and Mass Transfer, 26, 249-257. https://doi.org/10.1016/S0735-1933(99)00011-1
[24]	Liao, S.J. (2003) Beyond Perturbation: Introduction to Homotopy Analysis Method. Chapman/Hall, CRC Press, Boca Raton.
[25]	Koriko, O.K., Oreyeni, T., Omowaye, A.J. and Animasaun, I.L. (2016) Homotopy Analysis of MHD Free Convective Micropolar Fluid Flow along a Vertical Surface Embedded in Non-Darcian Thermally Stratified Medium. Open Journal of Fluid Dynamics, 6, 198-221.
[26]	Hilton, P.J. (1953) An Introduction to Homotopy Theory. Cambridge University Press, Cambridge.
[27]	Oyem, O.A. and Oreyeni, T. (2017) Homotopy Analysis of Magnetic Field Effect on Free Convection Flow Past a Semi-Infinite Flat Plate. Annals of Faculty Engineering Hunedoara—International Journal of Engineering, 4, 103-110.