Heat and Mass Transfer in Unsteady Boundary Layer Flow of Williamson Nanofluids

Abstract

In this paper, analytic approximation to the heat and mass transfer characteristics of a two-dimensional time-dependent flow of Williamson nanofluids over a permeable stretching sheet embedded in a porous medium has been presented by considering the effects of magnetic field, thermal radiation, and chemical reaction. The governing partial differential equations along with the boundary conditions were reduced to dimensionless forms by using suitable similarity transformation. The resulting system of ordinary differential equations with the corresponding boundary conditions was solved via the homotopy analysis method. The results of the study show that velocity, temperature, and concentration boundary layer thicknesses generally decrease as we move away from the surface of the stretching sheet and the Williamson parameter was found to retard the velocity but it enhances the temperature and concentration profiles near the surface. It was also found that increasing magnetic field strength, thermal radiation, or rate of chemical reaction speeds up the mass transfer but slows down the heat transfer rates in the boundary layer. The results of this study were compared with some previously published works under some restrictions, and they are found in excellent agreement.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

The term boundary layer flow refers to a kind of flow in a relatively narrow region near a solid surface where the effect of viscosity is significant. The theory of boundary layer flow was first introduced by the German scientist Prandtl in 1904. On the other hand, the field of magnetohydrodynamics (MHD) which deals with the interaction between the velocity of electrically conducting fluids and the electromagnetic fields was first introduced in 1942 by the Swedish physicist and electrical engineer Alfvén [1]. In 1961, Sakiadis [2] initiated the study of boundary layer flows over stretching surfaces and formulated the two-dimensional boundary layer equations. Following these pioneer works, many investigators have reported various useful study results. For instance, Kumar [3] used the Runge–Kutta fourth-order method along with shooting technique to analyze the effect of linear thermal stratification in stable stationary ambient fluid on steady MHD convective flow of a viscous incompressible electrically conducting fluid along a stretching sheet in the presence of mass transfer and magnetic effect. Also, Narasu et al. [4] used perturbation technique to obtain analytic solution for free convective unsteady fluid flow in the presence of thermal diffusion and chemical reaction past a vertical porous plate with heat source and slip effects.

In order to enhance the transfer of heat in thermal systems, Choi and Eastmann [5] introduced the concept of nanofluids in 1995 and showed experimentally that the embedding of nanoparticles (particles having diameters less than 100 nm) with the conventional fluids such as water, oil, or ethylene glycol mixture dramatically increases the thermal conductivity and heat transfer of the fluids. The resulting uniform dispersion and stable suspension of nanoparticles in the base fluid is referred to as a nanofluid. Due to their improved heat transfer properties, nanofluids are working fluids that have been used in many energy saving activities of modern industries, cost saving production of technological devices, and advanced life-saving medical treatments. Consequently, several studies have been conducted to study the boundary layer flow behavior of nanofluids. For instance, Ishak et al. [6] anticipated the boundary layer flow of nanofluids past a porous shrinking sheet. Ibrahim and Shanker [7] reported a numerical study of unsteady laminar boundary layer flow and heat transfer of a viscous incompressible fluid over stretching sheet. Sekhar et al. [8] studied the chemical reaction effects on unsteady MHD oscillatory slip flow in an optically thin fluid through a planer channel filled with saturated porous medium in the presence of a temperature-dependent heat source. Haile and Shankar [9] presented the boundary layer flow of a nanofluid through a porous medium subjected to magnetic field, thermal radiation, viscous dissipation, and chemical reaction effects. Ali et al. [10] applied the idea of Caputo–Fabrizio time fractional derivatives to magnetohydrodynamics free convection flow of generalized Walters’-B fluid over a static vertical plate. The unsteady MHD free convection flows of alumina-water and single-walled carbon nanotube-water nanofluids within a symmetrical wavy trapezoidal enclosure was reported by Job et al. [11]. On the other hand, Sathish Kumar et al. [12] investigated the effect of nonlinear thermal radiation on unsteady MHD flow between parallel plates. Saqib et al. [13] used the Laplace transform technique to obtain exact solutions for natural convection flow model of a hybrid nanofluid in two vertical infinite parallel plates. Other recent studies [14–16] on the effects of various thermophysical parameters on the boundary layer flow of nanoparticles have also been reported.

In 1929, Williamson [17] noticed that there are some fluid flows that show both viscous and elastic properties, now called pseudoplastic fluids. He then proposed an equation, called the Williamson model, and verified it experimentally to exhibit the shear thinning or pseudoplastic properties. The boundary layer flow of pesudoplastic fluids is of great interest due to its application in industry such as extrusion of polymer sheets, emulsion coated sheets like photographic films, solutions, melts of high molecular weight polymers, etc [18]. In 2013, Nadeem et al. [19] presented the modeling of a two-dimensional boundary layer equation for the flow of Williamson fluid past a linear and exponentially stretching sheet. Following these pioneer works, a number of researchers have been reporting useful results on boundary layer flow and transport mechanisms of Williamson nanofluids with various thermophysical effects. For instance, Nadeem and Hussain [20] analyzed the flow and heat transfer mechanisms of Williamson nanofluids. The effects of radiation and chemical reaction on the steady boundary layer flow of MHD Williamson fluid through porous medium toward a horizontal linearly stretching sheet in the presence of nanoparticles were investigated numerically by Krishnamurthy et al. [21]. They also studied these effects on Williamson nanofluid slip flow over a stretching sheet embedded in a porous medium [22]. Recently, Reddy et al. [23] examined the magnetohydrodynamic boundary layer flow with heat and mass transfer of Williamson nanofluid over a stretching sheet with variable thickness and variable thermal conductivity under the radiation effect. Kho et al. [24] analyzed the boundary layer flow of Williamson nanofluids past over a stretching sheet in the presence of thermal radiation effect. Shawky et al. [25] also examined the MHD flow with heat and mass transfer of Williamson nanofluids over stretching sheet through porous medium. Ibrahim and Gamachu [26] utilized the spectral Quasilinearization method to inspect a nonlinear convection flow of Williamson nanofluid past a radially stretching surface under the application of electric field.

However, to the best of the authors’ knowledge, no study has been reported on the analytic solution of time-dependent boundary layer flow of Williamson nanofluids over a permeable stretching sheet embedded in a porous medium with the effects of magnetic field, thermal radiation, and chemical reaction. Thus, motivated by the aforementioned works, the present study attempts to fill the existing gaps in this area. To this end, the work of Bibi et al. [27] has been generalized in a way that the flow model has been extended by considering effects of thermal radiation, chemical reaction and porosity of the medium. Moreover, analytic solutions were obtained by employing an efficient method, namely, the homotopy analysis method.

2. Flow Analysis and Mathematical Formulation

In this study, the effects of magnetic field, thermal radiation, and chemical reaction on a two-dimensional unsteady boundary layer flow of an incompressible and viscous Williamson nanofluid over a heated permeable stretching sheet embedded in a porous medium were examined. It was assumed that the magnetic Reynolds number is small in liquid metals and partially ionized fluids; and the effect of polarization of charges is not taken into account. So, the induced magnetic field is negligible in comparison with the applied magnetic field. It was also assumed that the flow is generated from a slit by stretching of a uniformly permeable and semi-infinite flat plate with one end fixed at the slit and embedded in an optically thick porous medium as shown in Figure 1.

[figure omitted; refer to PDF]

The Cartesian coordinate system has been used in a way that the x-axis is along the stretching sheet, y-axis is normal to the sheet, the origin is located at the slit, and the flow in the region $y\ge 0$ is considered. Now, assuming that the sheet starts stretching at $t=0$ and extends horizontally with nonuniform velocity, $U_w=ax/\left(1-\text{ct}\right)$, where a and c are positive constants with dimension (time)⁻¹, then the basic equations for the balance of mass, momentum, energy, and nanoparticle volume fraction of the flow problem can be expressed in vector form as follows.

2.1. Continuity Equation

The continuity equation for conservation of mass becomes$$\begin{array}{}\text{(1)}& \nabla \cdot \mathbf{V}=0,\end{array}$$where $\nabla$ is the differential operator and V is the flow velocity vector.

2.2. Conservation of Momentum Equation

The Navier–Stoke’s equation for the balance of linear momentum is given by$$\begin{array}{}\text{(2)}& {\rho }_f\left[\frac{\partial \mathbf{V}}{\partial t}+\left(\mathbf{V}\cdot \nabla \right)\mathbf{V}\right]=\nabla \cdot \mathbf{S}+\mu {▽}^2\mathbf{V}+\left(\mathbf{J}\times \mathbf{B}\right)-\frac{\mu }{K_0}\mathbf{V},\end{array}$$where ${\rho }_f$ is the density of the nanofluid; t is time; S is the Cauchy stress tensor; μ is the coefficient of dynamic viscosity; J is the current density; B is the external magnetic field; $\mathbf{J}\times \mathbf{B}$ is the Lorentz force produced by the interaction of the applied magnetic field with velocity of the fluid; and $K_0$ is permeability of the porous medium.

2.3. Conservation of Energy Equation

In the absence of heat source/sink, viscous dissipation, and Joule heating effects, the conservation of energy for heat transfer is given by$$\begin{array}{}\text{(3)}& {\left(\rho C_{\text{P}}\right)}_f\left[\frac{\partial T}{\partial t}+\left(\mathbf{V}\cdot \nabla \right)T\right]=\kappa {▽}^{2}T+{\left(\rho C_{\text{P}}\right)}_p\left[D_{\text{B}}\nabla C\cdot \nabla T+\frac{D_{\text{T}}}{T_{\infty }}\left(\nabla T\cdot \nabla T\right)\right]-\nabla q_{\text{r}},\end{array}$$where $C_{\text{P}}$ is specific heat at constant pressure; ${\left(\rho C_{\text{P}}\right)}_f$ and ${\left(\rho C_{\text{P}}\right)}_p$ are the specific heat capacity of the nanofluid and the nanoparticles, respectively; T is temperature of the nanofluid; κ is thermal conductivity; $D_{\text{B}}$ is the Brownian diffusion coefficient; $D_{\text{T}}$ is the thermophoresis diffusion coefficient; $T_{\infty }$ is the ambient temperature; and $q_{\text{r}}$ is the radiative heat flux.

2.4. Conservation of Nanoparticle Concentration Equation

For a homogeneous chemical reaction, the concentration equation for nanoparticle volume fraction becomes$$\begin{array}{}\text{(4)}& \frac{\partial C}{\partial t}+\left(\mathbf{V}\cdot \nabla \right)C=D_{\text{B}}{▽}^{2}C+\frac{D_T}{T_{\infty }}\left({▽}^{2}T\right)-K_{\text{r}}\left(C-C_{\infty }\right),\end{array}$$where C is nanoparticle volume fraction, and the rates $K_{\text{r}}>0$ and $K_{\text{r}}<0$ denote destructive and constructive reaction rates, respectively.

In order to reduce the momentum (2), the Williamson constitutive equation is employed for the Cauchy stress tensor as defined in [28]:$$\begin{array}{}\text{(5)}& \mathbf{S}=-p\mathbf{I}+{\tau }_1,\end{array}$$where p is pressure; I is unit tensor; and ${\tau }_1$ is a deviatoric extra stress tensor given by$$\begin{array}{}\text{(6)}& {\tau }_1=\left[{\mu }_{\infty }+\frac{{\mu }_0-{\mu }_{\infty }}{1-\Gamma \dot{\gamma }}\right]A_1,\end{array}$$where ${\mu }_0$ and ${\mu }_{\infty }$ are limiting viscosities at zero and infinite shear rates, respectively; $\Gamma >0$ is the material constant of the Williamson fluid; $A_1$ is first Rivlin–Ericksen tensor; and $\dot{\gamma }=\sqrt{\Pi /2}$ is the shear rate with $\Pi =$ trace $\left(A_1^2\right)$. Here, for psedo-plastic fluids, we consider ${\mu }_{\infty }=0$ and $\Gamma \dot{\gamma }<1$ which gives ${\tau }_1=\left[{\mu }_0/\left(1-\Gamma \dot{\gamma }\right)\right]A_1.$ Using binomial expansion, we get$$\begin{array}{}\text{(7)}& {\tau }_1={\mu }_0\left(1+\Gamma \dot{\gamma }\right)A_1.\end{array}$$

Details of the Williamson fluid model can be found in [17], [29], and [19].

For an optically thick fluid, the energy (3) can also be reduced by applying the Rosseland diffusion approximation for the radiative heat flux given by$$\begin{array}{}\text{(8)}& q_{\text{r}}=-\frac{4{\sigma }^{\ast }}{3k^{\ast }}\frac{\partial T^4}{\partial y},\end{array}$$where $k^{\text{*}}$ is the mean absorption constant and ${\sigma }^{\text{*}}$ is Stefan–Boltzmann constant. If the temperature difference within the flow is small, then we can expand $T^4$ in Taylor series about $T_{\infty }$, which after neglecting higher order terms takes the form$$\begin{array}{}\text{(9)}& T^4\approx 4T_{\infty }^{3}T-3T_{\infty }^4.\end{array}$$

So, the heat flux can be approximated as$$\begin{array}{}\text{(10)}& q_{\text{r}}\approx -\frac{16T_{\infty }^3{\sigma }^{\ast }}{3k^{\ast }}\frac{\partial T}{\partial y}.\end{array}$$

Assume a time-dependent magnetic field $\mathbf{B}\left(t\right)=B_0/\left(\sqrt{1-\text{ct}}\right)$ is applied in the y-direction, where $B_0$ is the initial magnetic field strength and $\text{ct}<1$. Then by using the boundary layer approximation, we rewrite the above conservation laws in two-dimensional Cartesian coordinates by using the equations given by Nadeem and Hussain [20] and Bibi et al. [27] as$$\begin{array}{}\text{(11)}& \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0,\text{(12)}& \frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=\upsilon \frac{{\partial }^{2}u}{\partial y^2}+\sqrt{2}\upsilon \Gamma \frac{\partial u}{\partial y}\frac{{\partial }^{2}u}{\partial y^2}-\left[\frac{\sigma B_0^2}{{\rho }_f\left(1-\text{ct}\right)}+\frac{\upsilon }{K_0}\right]u,\text{(13)}& \frac{\partial T}{\partial t}+u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}={\alpha }_f\frac{{\partial }^{2}T}{\partial y^2}+\tau \left[D_{\text{B}}\frac{\partial C}{\partial y}\frac{\partial T}{\partial y}+\frac{D_T}{T_{\infty }}{\left(\frac{\partial T}{\partial y}\right)}^2\right]+\frac{16{\sigma }^{\ast }{T}_{\infty }^3}{3{\left(\rho C_{\text{p}}\right)}_f{k}^{\ast }}\frac{{\partial }^{2}T}{\partial y^2},\text{(14)}& \frac{\partial C}{\partial t}+u\frac{\partial C}{\partial x}+v\frac{\partial C}{\partial y}=D_{\text{B}}\frac{{\partial }^{2}C}{\partial y^2}+\frac{D_{\text{T}}}{T_{\infty }}\frac{{\partial }^{2}T}{\partial y^2}-K_{\text{r}}\left(C-C_{\infty }\right),\end{array}$$where $\left(u,v\right)$ are the velocity components along x and y directions, respectively; $\upsilon =\mu /{\rho }_f$ is the kinematic viscosity; σ is electric conductivity of the fluid; ${\alpha }_f=\kappa /{\left(\rho C_{\text{p}}\right)}_f$ is thermal diffusivity; and $\tau ={\left(\rho C_{\text{p}}\right)}_p/{\left(\rho C_{\text{p}}\right)}_f$ is ratio of effective heat capacities of nanoparticle and the ordinary fluid.

The boundary conditions for the values of velocity, temperature, and nanoparticle concentration were taken as in Bibi et al. [27].

At y = 0,$$\begin{array}{c}\text{(15)}& u=U_{\text{w}}\left(x,t\right)=\frac{ax}{1-\text{ct}},v=V_{\text{w}}\left(t\right)=-\frac{V_0}{\sqrt{1-\text{ct}}},\\ T=T_{\text{w}}\left(x,t\right)=T_{\infty }+\frac{T_0{U}_{\text{w}}x}{\upsilon \sqrt{1-\text{ct}}},\\ C=C_{\text{w}}\left(x,t\right)=C_{\infty }+\frac{C_0{U}_{\text{w}}x}{\upsilon \sqrt{1-\text{ct}}},\end{array}$$and as $y⟶\infty$,$$\begin{array}{}\text{(16)}& u⟶0,T⟶T_{\infty },C⟶C_{\infty },\end{array}$$where $U_{\text{w}},T_{\text{w}},C_{\text{w}}$ denote the surface velocity, temperature, and nanoparticle concentration, respectively; $V_{\text{w}}$ is the mass transmission at the surface of the plate; a and $a/\left(1-\text{ct}\right)$ are the initial and the effective stretching rates of the sheet, respectively; $V_0$ is the constant value of velocity; $T_0,C_0$ are the constant values; and $T_{\infty },C_{\infty }$ are the ambient values of temperature and nanoparticle concentration, respectively.

In order to simplify the mathematical analysis, the following similarity transformations were introduced:$$\begin{array}{}\text{(17)}& \eta =y\sqrt{\frac{U_{\text{w}}}{\upsilon x}},\psi \left(x,y,t\right)=\sqrt{U_{\text{w}}\upsilon x}f\left(\eta \right),\theta \left(\eta \right)=\frac{T-T_{\infty }}{T_{\text{w}}-T_{\infty }},\phi \left(\eta \right)=\frac{C-C_{\infty }}{C_{\text{w}}-C_{\infty }},\end{array}$$where η is the dimensionless similarity variable; $f\left(\eta \right)$ is the dimensionless stream function; and $\theta \left(\eta \right)$ and $\phi \left(\eta \right)$ are the dimensionless temperature and nanoparticles volume fraction, respectively.

If we choose the stream function $\psi \left(x,y\right)$ such that the velocity components are related as $u=\partial \psi /\partial y$ and $v=-\partial \psi /\partial x$, then the continuity (11) is identically satisfied. Computing the required partial derivatives with respect to the new similarity variable η and substituting the values into equations (6)–(8), the governing system of partial differential equations are reduced to the following set of ordinary differential equations:$$\begin{array}{}\text{(18)}& \left(1+\text{We}{f}^{″}\right)f^{‴}+f{f}^{″}-{f^{\prime }}^2-A\left(f^{\prime }+\frac{\eta }{2}{f}^{″}\right)-\left(M^2+K_{\text{p}}\right)f^{\prime }=0,\text{(19)}& \left(1+\frac{4}{3}{R}_{\text{d}}\right){\theta }^{″}+\mathrm{Pr}\left[f{\theta }^{\prime }-2f^{\prime }\theta -\frac{A}{2}\left(3\theta +\eta {\theta }^{\prime }\right)+N_{\text{b}}{\theta }^{\prime }{\phi }^{\prime }+N_{\text{t}}{\left({\theta }^{\prime }\right)}^2\right]=0,\text{(20)}& {\phi }^{″}+\frac{N_{\text{t}}}{N_{\text{b}}}{\theta }^{″}+\text{Sc}\left[f{\phi }^{\prime }-2f^{\prime }\phi -\frac{A}{2}\left(3\phi +\eta {\phi }^{\prime }\right)-\gamma \phi \right]=0,\end{array}$$where the prime′ denotes differentiation with respect to η; $\text{We}=\Gamma x\sqrt{2a^3/\upsilon {\left(1-\text{ct}\right)}^3}$ is the Weissenberg number representing the Williamson parameter; $A=c/a$ is the unsteadiness parameter; $M=\sqrt{\sigma B_0^2/a{\rho }_f}$ is the magnetic parameter; $K_{\text{p}}=\upsilon x/U_{\text{w}}{K}_0$ is the porosity parameter; $R_{\text{d}}=4{\sigma }^{\ast }{T}_{\infty }^3/\kappa k^{\ast }$ is the thermal radiation parameter; $P_{\text{r}}=\upsilon /{\alpha }_m$ is the Prandtl number; $N_{\text{b}}=\tau D_{\text{B}}\left(C_{\text{w}}-C_{\infty }\right)/\upsilon$ is the Brownian motion parameter; $N_{\text{t}}=\tau D_{\text{T}}\left(T_{\text{w}}-T_{\infty }\right)/\upsilon T_{\infty }$ is the thermophoresis parameter; $\text{Sc}=\upsilon /D_{\text{B}}$ is the Schmidt number; and $\gamma =K_{\text{r}}x/U_{\text{w}}$ is the chemical reaction parameter with $\gamma >0$ and $\gamma <0$ denoting destructive and generative chemical reaction rates, respectively. Also using the similarity transformation in equation (17), the boundary conditions in equations (15) and (16) can be reduced as$$\begin{array}{c}\text{(21)}& f\left(0\right)=S,f^{\prime }\left(0\right)=1,\\ \theta \left(0\right)=1,\\ \phi \left(0\right)=1,\\ f^{\prime }\left(\eta \right)⟶0,\\ \theta \left(\eta \right)⟶0,\\ \phi \left(\eta \right)⟶0,\hspace{1em}\text{as\hspace{0.17em}}\eta ⟶\infty ,\end{array}$$where $S=V_0/\sqrt{a\upsilon }$ is suction (for $S<0$) or injection (for $s>0$) parameter.

From engineering point of view, it is useful to examine the impacts of skin friction coefficient $C_{\text{f}}$, local Nusselt number ${\text{Nu}}_x$, and Sherwood number ${\text{Sh}}_x$ on the boundary layer profiles of velocity, temperature, and concentration, respectively.

The skin friction coefficient $C_{\text{f}}$ is defined as$$\begin{array}{}\text{(22)}& c_{\text{f}}=\frac{{\tau }_{\text{w}}}{\left(1/2\right){\rho }_{\text{f}}{U}_{\text{w}}^2},\end{array}$$where ${\tau }_{\text{w}}$ is the shear stress at the permeable surface. For a Williamson fluid, it is given by$$\begin{array}{}\text{(23)}& {\tau }_{\text{w}}=-{\mu }_0{\left[\frac{\partial u}{\partial y}+\frac{\Gamma }{\sqrt{2}}{\left(\frac{\partial u}{\partial y}\right)}^2\right]}_{y=0}.\end{array}$$

Up on substitution, this gives$$\begin{array}{}\text{(24)}& -\frac{1}{2}{c}_{\text{f}}{\mathrm{Re}}_x^{1/2}=f^{″}\left(0\right)+\frac{\text{We}}{2}{f^{″}}^2\left(0\right),\end{array}$$where ${\mathrm{Re}}_x=a{x}^2/\upsilon \left(1-\text{ct}\right)$ is the local Reynolds number.

The local Nusselt number is a dimensionless heat transfer coefficient defined as the ratio of convective and conductive heat transfer rates$$\begin{array}{}\text{(25)}& {\text{Nu}}_x=\frac{x{q}_{\text{w}}}{\kappa \left(T_{\text{w}}-T_{\infty }\right)},\end{array}$$where $q_{\text{w}}=-\left(\kappa +16\sigma \ast T_{\infty }^3/3k^{\ast }\right){\left[\partial T/\partial y\right]}_{y=0}$ is the heat flux that measures the heat transfer at the surface of a permeable stretching sheet. So, by substitution this gives$$\begin{array}{}\text{(26)}& {\text{Nu}}_x=-\left(1+\frac{4}{3}{R}_{\text{d}}\right){\mathrm{Re}}_x^{1/2}{\theta }^{\prime }\left(0\right).\end{array}$$

On the other hand, the local Sherwood number ${\text{Sh}}_x$ is defined as the ratio of convective and diffusive mass transfer rates along a surface is which is given by$$\begin{array}{}\text{(27)}& {\text{Sh}}_x=\frac{x{J}_{\text{w}}}{D_{\text{B}}\left(C_{\text{w}}-C_{\infty }\right)},\end{array}$$where $J_{\text{w}}=-D_{\text{B}}{\left[\partial C/\partial y\right]}_{y=0}$ is the mass flux at the surface. Thus, we get$$\begin{array}{}\text{(28)}& {\mathrm{Re}}_x^{-\left(1/2\right)}{\text{Sh}}_x=-{\phi }^{\prime }\left(0\right).\end{array}$$

3. Method of Solution

The homotopy analysis method (HAM), first proposed by Liao in 1992, is one of the relatively recent and powerful analytic methods that has shown great efficiency over the past few years. Based on the concept of homotopy in topology, the method provides us extremely great freedom to choose different solution expressions from which we can choose the one that approximates our solution more efficiently. In 1997, Liao introduced a very useful parameter, called convergence-control parameter, denoted by $\hslash$ or $c_0$, which provides a more convenient approach to ensure the convergence of HAM solutions and to adjust the rate and region of convergence of the solution.

Now, to implement the HAM in this problem, we choose the initial approximations in the form$$\begin{array}{c}\text{(29)}& f_0\left(\eta \right)=S+1-e^{-\eta },{\theta }_0\left(\eta \right)=e^{-\eta },\\ {\phi }_0\left(\eta \right)=e^{-\eta }.\end{array}$$

We then select the auxiliary linear operators as$$\begin{array}{c}\text{(30)}& L_f=\frac{{\text{d}}^{3}f}{\text{d}{\eta }^3}+\frac{{\text{d}}^{2}f}{\text{d}{\eta }^2},L_{\theta }=\frac{{\text{d}}^2\theta }{\text{d}{\eta }^2}+\frac{\text{d}\theta }{\text{d}\eta },\\ L_{\phi }=\frac{{\text{d}}^2\phi }{\text{d}{\eta }^2}+\frac{\text{d}\phi }{\text{d}\eta },\end{array}$$with the property $L_{\text{f}}\left[C_1+C_2\eta +C_3{e}^{-\eta }\right]=0,L_{\theta }\left[C_4+C_5{e}^{-\eta }\right]=0,L_{\phi }\left[C_6+C_7{e}^{-\eta }\right]=0$; where $C_i\left(i=1-7\right)$ are integral constants to be determined by the boundary conditions.

We also choose the nonzero auxiliary functions as$$\begin{array}{}\text{(31)}& H_f\left(\eta \right)=H_{\theta }\left(\eta \right)=H_{\phi }\left(\eta \right)=e^{-\eta }.\end{array}$$

Finally, based on the ordinary differential equations (12)–(14), the nonlinear operators can be defined as follows:$$\begin{array}{c}\text{(32)}& N_f=\left(1+\text{We}\frac{{\partial }^2{\varphi }_f}{\partial {\eta }^2}\right)\frac{{\partial }^3{\varphi }_f}{\partial {\eta }^3}+{\varphi }_f\frac{{\partial }^2{\varphi }_f}{\partial {\eta }^2}-{\left(\frac{\partial {\varphi }_f}{\partial \eta }\right)}^2-A\left(\frac{\partial {\varphi }_f}{\partial \eta }+\frac{\eta }{2}\frac{{\partial }^2{\varphi }_f}{\partial {\eta }^2}\right)-\left(M^2+K_{\text{p}}\right)\frac{\partial {\varphi }_f}{\partial \eta },N_{\theta }=\left(1+\frac{4}{3}{R}_{\text{d}}\right)\frac{{\partial }^2{\varphi }_{\theta }}{\partial {\eta }^2}+P_{\text{r}}\left[{\varphi }_f\frac{\partial {\varphi }_{\theta }}{\partial \eta }-2{\varphi }_{\theta }\frac{\partial {\varphi }_f}{\partial \eta }-\frac{A}{2}\left(3{\varphi }_{\theta }+\eta \frac{\partial {\varphi }_{\theta }}{\partial \eta }\right)+N_{\text{b}}\frac{\partial {\varphi }_{\theta }}{\partial \eta }\frac{\partial {\varphi }_{\phi }}{\partial \eta }+N_{\text{t}}{\left(\frac{\partial {\varphi }_{\theta }}{\partial \eta }\right)}^2\right],\\ N_{\phi }=\frac{{\partial }^2{\varphi }_{\phi }}{\partial {\eta }^2}+\frac{N_{\text{t}}}{N_{\text{b}}}\frac{{\partial }^2{\varphi }_{\theta }}{\partial {\eta }^2}+S_{\text{c}}\left[{\varphi }_f\frac{\partial {\varphi }_{\phi }}{\partial \eta }-2{\varphi }_{\phi }\frac{\partial {\varphi }_f}{\partial \eta }-\frac{A}{2}\left(3{\varphi }_{\phi }+\eta \frac{\partial {\varphi }_{\phi }}{\partial \eta }\right)-\gamma {\varphi }_{\phi }\right],\end{array}$$where ${\varphi }_f\left(\eta ;q\right),{\varphi }_{\theta }\left(\eta ;q\right)$, and ${\varphi }_{\phi }\left(\eta ;q\right)$ are the homotopy approximations of $f,\theta$, and φ, respectively, satisfying$$\begin{array}{c}\text{(33)}& {\varphi }_f\left(0;q\right)=S,{\varphi }_f^{\prime }\left(0;q\right)=1,\\ {\varphi }_{\theta }\left(0;q\right)=1,\\ {\varphi }_{\phi }\left(0;q\right)=1,\\ {\varphi }_f^{\prime }\left(\eta ;q\right)⟶0,\\ {\varphi }_{\theta }\left(\eta ;q\right)⟶0,\\ {\varphi }_{\phi }\left(\eta ;q\right)⟶0,\hspace{1em}\text{as\hspace{0.17em}}\eta ⟶\infty ,\end{array}$$with $q\in \left[\mathrm{0,1}\right]$ is called the embedding parameter.

According to Liao, the corresponding zeroth-order deformation equations can be constructed as$$\begin{array}{}\text{(34)}& \left(1-q\right)L_f\left[{\varphi }_f-f_0\left(\eta \right)\right]={\hslash }_f{H}_f{N}_f\left[{\varphi }_f,{\varphi }_{\theta },{\varphi }_{\phi }\right],\text{(35)}& \left(1-q\right)L_{\theta }\left[{\varphi }_{\theta }-{\theta }_0\left(\eta \right)\right]={\hslash }_{\theta }{H}_{\theta }{N}_{\theta }\left[{\varphi }_f,{\varphi }_{\theta },{\varphi }_{\phi }\right],\text{(36)}& \left(1-q\right)L_{\phi }\left[{\varphi }_{\phi }-{\phi }_0\left(\eta \right)\right]={\hslash }_{\phi }{H}_{\phi }{N}_{\phi }\left[{\varphi }_f,{\varphi }_{\theta },{\varphi }_{\phi }\right].\end{array}$$

Obviously, when $q=0$, the solutions to equations (25)–(27) are ${\varphi }_f\left(\eta ;0\right)=f_0\left(\eta \right),{\varphi }_{\theta }\left(\eta ;0\right)={\theta }_0\left(\eta \right)$, and ${\varphi }_{\phi }\left(\eta ;0\right)={\phi }_0\left(\eta \right)$. Similarly, when $q=1$, the solutions to equations (25)–(27) are given by ${\varphi }_f\left(\eta ;1\right)=f\left(\eta \right),{\varphi }_{\theta }\left(\eta ;1\right)=\theta \left(\eta \right)$, and ${\varphi }_{\phi }\left(\eta ;1\right)=\phi \left(\eta \right)$. That is, as the embedding parameter q increases from 0 to 1, the homotopy solutions ${\varphi }_f\left(\eta ;q\right),{\varphi }_{\theta }\left(\eta ;q\right)$, and ${\varphi }_{\phi }\left(\eta ;q\right)$ vary continuously from the known initial approximations $f_0\left(\eta \right),{\theta }_0\left(\eta \right)$, and ${\phi }_0\left(\eta \right)$ to the unknown exact solutions $f\left(\eta \right),\theta \left(\eta \right)$, and $\phi \left(\eta \right)$, respectively. Now, if this continuous variation is smooth enough, then we expand ${\varphi }_f\left(\eta ;q\right),{\varphi }_{\theta }\left(\eta ;q\right)$, and ${\varphi }_{\phi }\left(\eta ;q\right)$ in the Taylor series about q as follows:$$\begin{array}{}\text{(37)}& {\varphi }_f\left(\eta ;q\right)={\varphi }_f\left(\eta ;0\right)+\sum _{m=1}^{\infty }\frac{1}{m!}{\left.\frac{{\partial }^m{\varphi }_f\left(\eta ;q\right)}{\partial {\eta }^m}\right|}_{q=0}{q}^m,\text{(38)}& {\varphi }_{\theta }\left(\eta ;q\right)={\varphi }_{\theta }\left(\eta ;0\right)+\sum _{m=1}^{\infty }\frac{1}{m!}{\left.\frac{{\partial }^m{\varphi }_{\theta }\left(\eta ;q\right)}{\partial {\eta }^m}\right|}_{q=0}{q}^m,\text{(39)}& {\varphi }_{\phi }\left(\eta ;q\right)={\varphi }_{\phi }\left(\eta ;0\right)+\sum _{m=1}^{\infty }\frac{1}{m!}{\left.\frac{{\partial }^m{\varphi }_{\phi }\left(\eta ;q\right)}{\partial {\eta }^m}\right|}_{q=0}{q}^m.\end{array}$$

As noted in Liao [30], the convergence of the above homotopy series depend on the auxiliary parameters ${\hslash }_f,{\hslash }_{\theta }$, and ${\hslash }_{\phi }$. Assuming that these parameters are properly chosen so that the series in equations (28)–(30) converge at $q=1$. Then using the deformations of ${\varphi }_f\left(\eta ;q\right),{\varphi }_{\theta }\left(\eta ;q\right)$, and ${\varphi }_{\phi }\left(\eta ;q\right)$, the homotopy series solutions become$$\begin{array}{}\text{(40)}& f\left(\eta \right)={\varphi }_f\left(\eta ;1\right)\approx f_0\left(\eta \right)+\sum _{m=1}^{\infty }{f}_m\left(\eta \right),\text{(41)}& \theta \left(\eta \right)={\varphi }_{\theta }\left(\eta ;1\right)\approx {\theta }_0\left(\eta \right)+\sum _{m=1}^{\infty }{\theta }_m\left(\eta \right),\text{(42)}& \phi \left(\eta \right)={\varphi }_{\phi }\left(\eta ;1\right)\approx {\phi }_0\left(\eta \right)+\sum _{m=1}^{\infty }{\phi }_m\left(\eta \right),\end{array}$$where $f_m\left(\eta \right)={\left.{\partial }^m{\varphi }_f\left(\eta ;q\right)/\partial {\eta }^m\right|}_{q=0},{\theta }_m\left(\eta \right)={\left.{\partial }^m{\varphi }_{\theta }\left(\eta ;q\right)/\partial {\eta }^m\right|}_{q=0}$ and ${\phi }_m\left(\eta \right)={\left.{\partial }^m{\varphi }_{\phi }\left(\eta ;q\right)/\partial {\eta }^m\right|}_{q=0}$ are called the mth-order homotopy derivatives. Here, we require these derivatives to satisfy the boundary conditions:$$\begin{array}{c}\text{(43)}& f_m\left(0\right)=S,f_m^{\prime }\left(0\right)=1,\\ {\theta }_m\left(0\right)=1,\\ {\phi }_m\left(0\right)=1,\\ f_m^{\prime }\left(\eta \right)⟶0,\\ {\theta }_m\left(\eta \right)⟶0,\\ {\phi }_m\left(\eta \right)⟶0,\hspace{1em}\text{as\hspace{0.17em}}\eta ⟶\infty .\end{array}$$

It is important to note that the equations in (31)–(33) provide us with the relationship between the exact solutions and the initial approximations.

In order to get the terms of the series in equations (28)–(30), we first define the vectors satisfying the boundary conditions (43) as$$\begin{array}{c}\text{(44)}& {\stackrel{⟶}{f}}_{m-1}\left(\eta \right)=\left\{f_0\left(\eta \right),f_1\left(\eta \right),f_2\left(\eta \right),\dots ,f_{m-1}\left(\eta \right)\right\},{\stackrel{⟶}{\theta }}_{m-1}\left(\eta \right)=\left\{{\theta }_0\left(\eta \right),{\theta }_1\left(\eta \right),{\theta }_2\left(\eta \right),\dots {\theta }_{m-1}\left(\eta \right)\right\},\\ {\stackrel{⟶}{\phi }}_{m-1}\left(\eta \right)=\left\{{\phi }_0\left(\eta \right),{\phi }_1\left(\eta \right),{\phi }_2\left(\eta \right),\dots ,{\phi }_{m-1}\left(\eta \right)\right\}.\end{array}$$

Differentiating the zeroth-order deformation equations (25)–(27) m times with respect to q, then setting $q=0$ and finally dividing the resulting equations by $m!$, we get the so-called mth-order deformation equations given by$$\begin{array}{}\text{(45)}& L_f\left[f_m-{\chi }_m{f}_{m-1}\right]={\hslash }_f{H}_f{R}_{m-1}^f\left({\stackrel{⟶}{f}}_{m-1},{\stackrel{⟶}{\theta }}_{m-1},{\stackrel{⟶}{\phi }}_{m-1}\right),\text{(46)}& L_{\theta }\left[{\theta }_m-{\chi }_m{\theta }_{m-1}\right]={\hslash }_{\theta }{H}_{\theta }{R}_{m-1}^{\theta }\left({\stackrel{⟶}{f}}_{m-1},{\stackrel{⟶}{\theta }}_{m-1},{\stackrel{⟶}{\phi }}_{m-1}\right),\text{(47)}& L_{\phi }\left[{\phi }_m-{\chi }_m{\phi }_{m-1}\right]={\hslash }_{\phi }{H}_{\phi }{R}_{m-1}^{\phi }\left({\stackrel{⟶}{f}}_{m-1},{\stackrel{⟶}{\theta }}_{m-1},{\stackrel{⟶}{\phi }}_{m-1}\right),\end{array}$$where ${\chi }_m=\left\{\begin{array}{cc}0& \text{if\hspace{0.17em}}m\le 1\\ 1& \text{if\hspace{0.17em}}m>1\end{array}\right.$ is the unit step function,$$\begin{array}{c}\text{(48)}& R_{m-1}^f=f_{m-1}^{‴}+\text{We}\sum _{k=0}^{m-1}{f}_{m-1}^{″}{f}_{m-1-k}^{‴}+\sum _{k=0}^{m-1}{f}_k{f}_{m-1-k}^{″}-\sum _{k=0}^{m-1}{f}_k^{\prime }{f}_{m-1-1}^{\prime }-A\left(f_{m-1}^{\prime }+\frac{\eta }{2}{f}_{m-1}^{″}\right)-\left(M^2+K_p\right)f_{m-1}^{\prime },R_{m-1}^{\theta }=\left(1+\frac{4}{3}R\text{d}\right){\theta }_{m-1}^{″}+\mathrm{Pr}\left[\sum _{k=0}^{m-1}{f}_k{\theta }_{m-1-k}^{\prime }-2\sum _{k=0}^{m-1}{\theta }_k{f}_{m-1-k}^{\prime }-\frac{A}{2}\left(3{\theta }_{m-1}+\eta {\theta }_{m-1}^{\prime }\right)+N_{\text{b}}\sum _{k=0}^{m-1}{\theta }_{m-1}^{\prime }{f}_{m-1-k}^{\prime }\right]+\mathrm{Pr}\left[N_{\text{t}}\sum _{k=0}^{m-1}{\theta }_k^{\prime }{\theta }_{m-1-k}^{\prime }\right],\end{array}$$ $$\begin{array}{}\text{(49)}& R_{m-1}^{\phi }={\phi }_{m-1}^{″}+\frac{N_{\text{t}}}{N_{\text{b}}}{\theta }_{m-1}^{″}+\text{Sc}\left[\sum _{r=0}^{m-1}{f}_{\text{r}}{\phi }_{m-1-r}^{\prime }-2\sum _{r=0}^{m-1}{\phi }_{\text{r}}{f}_{m-1-r}^{\prime }-\frac{A}{2}\left(3{\phi }_{m-1}+\eta {\phi }_{m-1}^{\prime }\right)-\gamma {\phi }_{m-1}\right].\end{array}$$

Taking the inverse of the linear operators on both sides of the high-order deformation equations (35)–(37), we obtain the following iterative formula:$$\begin{array}{c}\text{(50)}& f_m\left(\eta \right)={\chi }_m{f}_{m-1}\left(\eta \right)+{\hslash }_f{L}_f^{-1}\left[H_f\left(\eta \right)R_{m-1}^f\left({\stackrel{⟶}{f}}_{m-1}\left(\eta \right),{\stackrel{⟶}{\theta }}_{m-1}\left(\eta \right),{\stackrel{⟶}{\phi }}_{m-1}\left(\eta \right)\right)\right],{\theta }_m\left(\eta \right)={\chi }_m{\theta }_{m-1}\left(\eta \right)+{\hslash }_{\theta }{L}_{\theta }^{-1}\left[H_{\theta }\left(\eta \right)R_{m-1}^{\theta }\left({\stackrel{⟶}{f}}_{m-1}\left(\eta \right),{\stackrel{⟶}{\theta }}_{m-1}\left(\eta \right),{\stackrel{⟶}{\phi }}_{m-1}\left(\eta \right)\right)\right],\\ {\phi }_m\left(\eta \right)={\chi }_m{\phi }_{m-1}\left(\eta \right)+{\hslash }_{\phi }{L}_{\phi }^{-1}\left[H_{\phi }\left(\eta \right)R_{m-1}^{\phi }\left({\stackrel{⟶}{f}}_{m-1}\left(\eta \right),{\stackrel{⟶}{\theta }}_{m-1}\left(\eta \right),{\stackrel{⟶}{\phi }}_{m-1}\left(\eta \right)\right)\right].\end{array}$$

These equations can be solved recursively one after the other in order as a function of all the previous terms of the series. That is, for each $m=\mathrm{1,2},\dots$ , one can generate the terms of the series and the corresponding partial sums as orders of successive approximations.

Even if the proper values of the initial approximations, the auxiliary linear operators, and the auxiliary functions have been selected, the HAM series solutions still contain the auxiliary parameters ${\hslash }_f,{\hslash }_{\theta }$, and ${\hslash }_{\phi }$ that could influence the accuracy, region of convergence, and rate of approximation for the solutions. That is, we still have great freedom to choose the values of the auxiliary parameters. So it remains to choose appropriate values of ${\hslash }_f,{\hslash }_{\theta }$, and ${\hslash }_{\phi }$ that will ensure the convergence of the HAM solutions. In 2003, Liao [30] suggested that the HAM solutions will converge to the exact solutions in regions for which the graphs of the partial sums of the functions or their first few derivatives at some parameters against $\hslash$, called the $\hslash$-curves, are essentially horizontal. That is, the interval of least change in the physical quantity signifies the best interval for the choice of the convergence-control parameters.

To facilitate the computation, the HAM-based Mathematica package, BVPh 2.0 developed by Zhao and Liao in 2013 has been employed. This package has recently been used successfully by many researchers like Farooq and Zhi-Liang [31]; Demir et al. [32]; and Bano et al. [33]. For details regarding use of the package, we recommend the readers to refer Zhao and Liao [34].

It can be seen from Figure 2 that the $\hslash$-curves are nearly horizontal in the ranges$$\begin{array}{c}\text{(51)}& -1.2<{\hslash }_f<-0.3,-1.2<{\hslash }_{\theta }<-0.2,\\ -1.2<{\hslash }_{\phi }<-\mathrm{0.2.}\end{array}$$

[figure omitted; refer to PDF]

That is, these intervals are the valid regions for which the HAM solutions converge to the exact solutions. So, any value from these intervals may be picked up to give the approximate solutions of our flow problem. By comparing the results at some of these possible values of the parameters, one may choose the optimal values of the parameters for which the approximate solutions will converge in the whole region of η. More precisely, the optimal values of the convergence-control parameters can be obtained by minimizing the average squared residual errors of the governing equations at the 20th order of HAM approximation as indicated in Figure 2(b). Now, using the BVPh 2.0 package at 20th order HAM approximation, the optimal values for the convergence-control parameters were found to be$$\begin{array}{c}\text{(52)}& {\hslash }_f\approx -1.03433,{\hslash }_{\theta }\approx -0.879735,\\ {\hslash }_{\phi }\approx -\mathrm{0.767268.}\end{array}$$

These values and the values $M=0.2,\text{We}=0.2,A=0.3,\mathrm{Pr}=0.72,N_{\text{b}}=0.2,N_{\text{t}}=0.1,\text{Sc}=1,S=0.3,K_{\text{p}}=1,R_{\text{d}}=0.3\text{\hspace{0.17em}and\hspace{0.17em}}\gamma =2$ are utilized throughout this study unless stated otherwise.

It is clear from Table 1 that increasing the number of terms decreases the squared residual errors and leads the series to converge.

Table 1

Convergence of some HAM solutions.

Order of HAM approximation	$-f^{″}\left(0\right)$	$-{\theta }^{\prime }\left(0\right)$	$-{\phi }^{\prime }\left(0\right)$	Squared ${\epsilon }_f$	Residual ${\epsilon }_{\theta }$	Errors ${\epsilon }_{\phi }$
5	1.9352	1.0616	1.3121
10	1.9370	1.0531	1.3064	$1.84\times {10}^{-6}$	0.0000388	0.0000201
15	1.9381	1.0511	1.3059
20	1.9381	1.0499	1.3055	$7.35\times {10}^{-8}$	0.0000192	0.0000104
25	1.9381	1.0492	1.3053
30	1.9381	1.0487	1.3051	$2.28\times {10}^{-8}$	0.0000129	0.0000071
35	1.9381	1.0484	1.3050
40	1.9381	1.0481	1.3049	$9.91\times {10}^{-9}$	0.0000098	0.0000054

4. Results and Discussions

In this section, we present and discuss the main results of our mathematical analysis for the unsteady boundary layer flow of Williamson nanofluids induced by a permeable stretching plate embedded in a porous medium in the presence of magnetic field, thermal radiation, and chemical reaction. In particular, the effects of various thermophysical parameters on velocity, temperature, and concentration profiles were examined and presented as follows.

Since there are many theoretically and practically important flow field and transport phenomena that can be unsteady due to sudden movement of surfaces, we begin our result analysis by presenting the effects of unsteadiness parameter on our flow problem.

The result in Figure 3 explains that the velocity gradient, temperature, and concentration profiles are decreasing functions of the unsteadiness parameter in the boundary layer. This is due to the fact that as the unsteadiness parameter $A=c/a$ increases, the velocity $U_{\text{w}}=ax/1-\text{ct}$ of the stretching sheet also decreases which causes the transfer of less amount of heat and mass from the sheet to the nanofluid in the boundary layer region.

[figure omitted; refer to PDF]

The influence of magnetic field on the unknown functions of the given problem is illustrated in Figure 4. It reveals that the increase in magnetic parameter increases the temperature and concentration profiles but reduces the velocity gradient.

[figure omitted; refer to PDF]

This result agrees well with the fact that the increase in magnetic field strength tends to increase a dragging force which resists the motion of the fluid particles and produces heat in the boundary layer.

The effect of Williamson viscoelastic parameter in the flow field can be described in terms of the Weissenberg number $\text{We}$ as illustrated in Figure 5.

[figure omitted; refer to PDF]

The result in Figure 5 reveals that the velocity gradient decreases near the stretching sheet as the Weissenberg number $\text{We}$ increases. The temperature and concentration profiles are increasing with increasing of $\text{We}$ in the boundary layer.

The effect of porosity $K_{\text{p}}$ in the flow problem has been investigated as presented in Figure 6.

[figure omitted; refer to PDF]

Figure 6 indicates that as porosity of the medium increases, the temperature and concentration profiles also increase in the boundary layer. On the other hand, the velocity gradient is found to decrease near the stretching surface and increase away from the surface with increasing porosity of the medium. This is due to the fact that increasing the value of $K_{\text{p}}$ has the tendency to resist the flow.

The result in Figure 7 displays that the rise in thermal radiation increases the temperature but reduces the concentration profiles. This is because the increase in thermal radiation provides more heat to the nanofluid which causes the increase in temperature and thermal boundary layer thickness. Furthermore, this parameter has no significant effect on the velocity gradient of the fluid flow.

[figure omitted; refer to PDF]

The chemical reaction that takes place in nanofluid flow results the conversion of species in the boundary layer. The effects of chemical reaction parameter on the concentration profile of our problem are illustrated in Figure 8.

The result in Figure 8 shows that the concentration profile increases with the generative chemical reaction and decreases with the destructive chemical reaction. This is true because as generative and destructive chemical reactions take place, the amount of nanoparticles within the fluid is getting larger and smaller, respectively.

[figure omitted; refer to PDF]

We see from Figure 9 that the velocity gradient, temperature profile, and concentration of nanoparticles are decreasing functions of the injection parameter.

[figure omitted; refer to PDF]

It can be observed from Table 2 that all the coefficients of skin friction, Nusselt number, and Sherwood number are increasing functions of A. Also, the skin friction coefficient can be maximized by increasing $M,K_{\text{p}}$, or γ; and the coefficient of Nusselt number is increased by decreasing $K_{\text{p}},R_{\text{d}}$, or the generative chemical reaction parameter γ. It is also shown that increasing $R_{\text{d}}$ or γ or decreasing $K_{\text{p}}$ raises the coefficient of Sherwood number.

Table 2

Coefficients of skin friction, Nusselt number, and Sherwood number.

A	M	K p	R d	γ	$-f^{″}\left(0\right)$	$-{\theta }^{\prime }\left(0\right)$	$-{\phi }^{\prime }\left(0\right)$
0.3					1.93811	1.04995	1.30552
0.6					2.04031	1.16875	1.42042
0.9	0.2				2.14202	1.27181	1.52403
	0.8				2.21440	1.02165	1.27691
	1.4				2.91874	1.21627	1.45894
		1			2.44662	1.25545	1.50555
		4			3.05318	1.21715	1.46277
		7			3.33303	1.46380	1.78493
			0.1		3.41649	1.32197	1.60133
			0.4		3.33312	1.20886	1.88957
			0.7	0.2	3.33277	1.03745	1.95725
				0.5	3.33281	1.04622	2.06716
				0.8	3.33283	1.04553	2.17143
				−0.2	3.33275	1.03952	1.78885
				−0.5	3.33271	1.04150	1.64232
				−0.8	3.33263	1.04410	1.46839

We further validate our solutions by comparing them with previously published works in some limited cases. If $\text{We}=M=R_{\text{d}}=K_{\text{p}}=R_{\text{d}}=\text{Sc}=S=\gamma =0$, then our analytic approximation to the skin friction coefficient, $-f^{″}\left(0\right)$, with some values of the unsteadiness parameter A can be compared with the results of some previous studies as shown in Table 3.

The comparison presented in Table 3 reveals that for the selected values of the parameters, the values of $-f^{″}\left(0\right)$ determined in this study are in excellent agreement with some of the previously published works.

Table 3

Comparisons of results on $-f^{″}\left(0\right)$ for $\text{We}=M=K_{\text{p}}=R_{\text{d}}=\text{Sc}=S=\gamma =0,\mathrm{Pr}=0.72,N_{\text{b}}=0.2,\text{\hspace{0.17em}Nt}=0.1$ against some values of A.

A	Sharidan et al. [35]	Chamkha et al. [36]	Mukhopadhyay and Gorla. [37]	Khan and Azam [38]	Bibi et al. [27]	Present study
0.0	—	—	—	1.00000	1.0005	1.00000
0.2	—	—	—	1.06801	1.0685	1.06874
0.4	—	—	—	1.13469	1.1349	1.13521
0.6	—	—	—	1.19912	1.1992	1.19930
0.8	1.26104	1.261512	1.261479	1.26104	—	1.26099
1.2	1.37772	1.378052	1.377850	1.37772	—	1.37755
2.0	—	—	—	1.58737	—	1.58740

5. Concluding Remarks

In this study, the homotopy analysis method (HAM) has been successfully applied to obtain analytical series solutions for the unsteady boundary layer flow of Williamson nanofluids induced by a permeable stretching sheet embedded in a porous medium in the presence of magnetic field, thermal radiation, and chemical reaction. The convergence of the HAM solutions is ensured by using successive iterations, $\hslash$-curves, and the squared residual error analysis. The validity of the results are also verified by comparing them with previous results under some restricted assumptions and found to be in excellent agreement. The effects of various pertinent thermophysical parameters on the dimensionless velocity, temperature, and concentration profiles are examined, and the results are summarized as follows:

(i)

Since the fluid flow is caused only by the stretching sheet, the velocity, temperature, and concentration boundary layer thickness generally decrease as we move away from the surface of the sheet.

(ii)

The velocity gradient of the flow problem is found to be a decreasing function of the unsteadiness parameter, magnetic field strength, and injection parameter. It also decreases near the stretching surface with the increase in porosity of the medium and Weissenberg number.

(iii)

The temperature profile can be increased in the boundary layer by increasing magnetic field strength, porosity of the medium, thermal radiation, or Weissenberg number. This profile can also be increased by decreasing the injection parameter or the unsteadiness parameter.

(iv)

The concentration profile of nanoparticles can be increased in the boundary layer by increasing the magnetic field strength, Weissenberg number, porosity of the medium, or generative chemical reaction. This profile can also be increased by decreasing the unsteadiness parameter, destructive chemical reaction, thermal radiation or injection parameter;

(v)

The rate of momentum transfer can be enhanced mainly by decreasing magnetic field strength, porosity of the medium, chemical reaction, or unsteadiness parameter; it can also be decreased by reducing the effect of thermal radiation.

(vi)

Heat transfer rate can be speeded up by increasing the unsteadiness parameter or by decreasing the porosity of the medium, thermal radiation, or destructive chemical reaction.

(vii)

Mass transfer rate can also be enhanced by increasing the thermal radiation, constructive chemical reaction, or unsteadiness parameter. It can also be facilitated by decreasing the effect of thermal radiation or destructive chemical reaction.

References

[1] H. Alfvén, "Existence of electromagnetic-hydrodynamic waves," Nature, vol. 150 no. 3805, pp. 405-406, DOI: 10.1038/150405d0, 1942.

[2] B. C. Sakiadis, "Boundary-layer behavior on continuous solid surfaces: I. Boundary-layer equations for two-dimensional and axisymmetric flow," AIChE Journal, vol. 7 no. 1, pp. 26-28, DOI: 10.1002/aic.690070108, 1961.

[3] R. Kumar, "MHD boundary layer flow on heat and mass transfer over a stretching sheet with slip effect," Journal of Naval Architecture and Marine Engineering, vol. 10 no. 2, pp. 109-118, DOI: 10.3329/jname.v10i2.16400, 2013.

[4] N. S. Kumar, R. Kumar, A. G. Vijaya Kumar, "Thermal diffusion and chemical reaction effects on unsteady flow past a vertical porous plate with heat source dependent in slip flow regime," Journal of Naval Architecture and Marine Engineering, vol. 13 no. 1, pp. 51-62, DOI: 10.3329/jname.v13i1.20773, 2016.

[5] S. U. S. Choi, J. A. Eastman, "Enhancing thermal conductivity of fluids with nanoparticles," Proceedings of the 1995 International Mechanical Engineering Congress and Exhibition, .

[6] A. Ishak, Y. Y. Lok, I. Pop, "Stagnation-point flow over a shrinking sheet in a micropolar fluid," Chemical Engineering Communications, vol. 197 no. 11, pp. 1417-1427, DOI: 10.1080/00986441003626169, 2010.

[7] W. Ibrahim, B. Shanker, "Unsteady boundary layer flow and heat transfer due to a stretching sheet by quasilinearization technique," World Journal of Mechanics, vol. 1 no. 6, pp. 288-293, DOI: 10.4236/wjm.2011.16036, 2011.

[8] K. Raja Sekhar, G. V. Ramana Reddy, B. D. C. N. Prasad, "Chemically reacting on MHD oscillatory slip flow in a planer channel with varying temperature and concentration," Advances in Applied Science Research, vol. 3 no. 5, pp. 2652-2659, 2012.

[9] E. Haile, B. Shankar, "Heat and mass transfer through a porous media of MHD flow of nanofluids with thermal radiation, viscous dissipation and chemical reaction effects," American Chemical Science Journal, vol. 4 no. 6, pp. 828-846, DOI: 10.9734/acsj/2014/11082, 2014.

[10] F. Ali, M. Saqib, I. Khan, N. Ahmad Sheikh, "Application of caputo-fabrizio derivatives to MHD free convection flow of generalized Walters’-B fluid model," The European Physical Journal Plus, vol. 131 no. 10,DOI: 10.1140/epjp/i2016-16377-x, 2016.

[11] V. M. Job, S. R. Gunakala, B. Rushi Kumar, R. Sivaraj, "Time-dependent hydromagnetic free convection nanofluid flows within a wavy trapezoidal enclosure," Applied Thermal Engineering, vol. 115, pp. 363-377, DOI: 10.1016/j.applthermaleng.2016.12.084, 2016.

[12] M. Sathish Kumar, N. Sandeep, B. Rushi Kumar, "Effect of nonlinear thermal radiation on unsteady MHD flow between parallel plates," Global Journal of Pure and Applied Mathematics, vol. 12, 2016.

[13] M. Saqib, I. Khan, S. Shafie, "Application of fractional differential equations to heat transfer in hybrid nanofluid: modeling and solution via integral transforms," Advances in Differential Equations, vol. 2019 no. 1,DOI: 10.1186/s13662-019-1988-5, 2019.

[14] M. Saqib, F. Ali, I. Khan, N. A. Sheikh, A. Khan, "Entropy generation in different types of fractionalized nanofluids," Arabian Journal for Science and Engineering, vol. 44 no. 1, pp. 531-540, DOI: 10.1007/s13369-018-3342-8, 2019.

[15] M. Saqib, I. Khan, S. Shafie, "Natural convection channel flow of CMC-based CNTs nanofluid," The European Physical Journal Plus, vol. 133,DOI: 10.1140/epjp/i2018-12340-3, 2018.

[16] M. Saqib, F. Ali, I. Khan, N. Ahmad Sheikh, S. Bin Shafie, "Convection in ethylene glycol-based molybdenum disulfide nanofluid: atangana–baleanu fractional derivative approach," Journal of Thermal Analysis and Calorimetry, vol. 135 no. 1, pp. 523-532, DOI: 10.1007/s10973-018-7054-9, 2019.

[17] R. V. Williamson, "The flow of pseudoplastic materials," Industrial & Engineering Chemistry, vol. 21 no. 11, pp. 1108-1111, DOI: 10.1021/ie50239a035, 1929.

[18] N. C. Sacheti, T. E. Bashir, P. Chandran, "Perturbation analysis of 2-dimensional boundary layer flow of an inelastic fluid using Williamson model," International Journal of Applied Engineering Research, vol. 22, pp. 12728-12734, 2017.

[19] S. Nadeem, S. T. Hussain, C. Lee, "Flow of a Williamson fluid over a stretching sheet," Brazilian Journal of Chemical Engineering, vol. 30 no. 3, pp. 619-625, DOI: 10.1590/s0104-66322013000300019, 2013.

[20] S. Nadeem, S. T. Hussain, "Flow and heat transfer analysis of Williamson nanofluid," Applied Nanoscience, vol. 4 no. 8, pp. 1005-1012, DOI: 10.1007/s13204-013-0282-1, 2014.

[21] M. R. Krishnamurthy, B. C. Prasannakumara, B. J. Gireesha, R. S. R. Gorla, "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, vol. 19 no. 1, pp. 53-61, DOI: 10.1016/j.jestch.2015.06.010, 2016.

[22] B. C. Prasannakumara, B. J. Gireesha, R. S. R. Gorla, M. R. Krishnamurthy, "Effects of chemical reaction and nonlinear thermal radiation on Williamson nanofluid slip flow over a stretching sheet embedded in a porous medium," Journal of Aerospace Engineering, vol. 29 no. 5,DOI: 10.1061/(asce)as.1943-5525.0000578, 2016.

[23] S. Reddy, K. Naikoti, M. M. Rashidi, "MHD flow and heat transfer characteristics of Williamson nanofluid over a stretching sheet with variable thickness and variable thermal conductivity," Transactions of A. Razmadze Mathematical Institute, vol. 171 no. 2, pp. 195-211, DOI: 10.1016/j.trmi.2017.02.004, 2017.

[24] Y. B. Kho, A. Hussanan, M. K. A. Mohamed, N. M. Sarif, Z. Ismail, M. Z. Salleh, "Thermal radiation effect on MHD flow and heat transfer analysis of Williamson nanofluid past over a stretching sheet with constant wall temperature," Journal of Physics: Conference Series, vol. 890,DOI: 10.1088/1742-6596/890/1/012034, 2017.

[25] H. M. Shawky, N. T. M. Eldabe, K. A. Kamel, E. A. Abd-Aziz, MHD Flow with Heat and Mass Transfer of Williamson Nanofluid over Stretching Sheet through Porous Medium, 2018.

[26] W. Ibrahima, D. Gamachu, "Nonlinear convection flow of Williamson nanofluid past a radially stretching surface," AIP Advances, vol. 9 no. 8,DOI: 10.1063/1.5113688, 2019.

[27] M. Bibi, K. U. Rehman, M. Y. Malik, M. Tahir, "Numerical study of Unsteady Williamson fluid flow and heat transfer in the presence of MHD through a permeable stretching surface," The European Physical Journal Plus, vol. 133 no. 4,DOI: 10.1140/epjp/i2018-11991-2, 2018.

[28] I. Dapra, G. Scarpi, "Perturbation solution for pulsatile flow of a non-Newtonian Williamson fluid in a rock fracture," International Journal of Rock Mechanics and Mining Sciences, vol. 44 no. 2, pp. 271-278, DOI: 10.1016/j.ijrmms.2006.07.003, 2007.

[29] C. Vasudev, U. R. Rao, M. V. S. Reddy, G. P. Rao, "Peristaltic pumping of Williamson fluid through a porous medium in a horizontal channel with heat transfer," American Journal of Scientific and Industrial Research, vol. 1 no. 3, pp. 656-666, DOI: 10.5251/ajsir.2010.1.3.656.666, 2010.

[30] S. J. Liao, Beyond Perturbation: Introduction to Homotopy Analysis Method, 2003.

[31] U. Farooq, L. Z. Liang, "Nonlinear heat transfer in a two-layer flow with nanofluids by OHAM," Journal of Heat Transfer, vol. 136 no. 2,DOI: 10.1115/1.4025432, 2014.

[32] M. Ş. Demir, M. S. Dokuz, S. Barış, "Semi-analytical solution for the MHD three dimensional stagnation-point flow of a Newtonian fluid towards a uniformly heated and moving vertical plate," Journal of Environmental Science, Computer Science and Engineering & Technology, vol. 5, pp. 39-51, 2016.

[33] N. Bano, B. B. Singh, S. R. Sayyed, "Homotopy analysis for MHD hiemenz flow in a porous medium with thermal radiation, velocity and thermal slips effects," Frontiers in Heat and Mass Transfer, vol. 10,DOI: 10.5098/hmt.10.14, 2018.

[34] Y. Zhao, S. J. Liao, "HAM-based mathematica package BVPh 2.0 for nonlinear boundary value problems," Advances in the Homotopy Analysis Method, vol. 361, 2013.

[35] S. Sharidan, T. Mahmood, I. Pop, "Similarity solutions for the unsteady boundary layer flow and heat transfer due to a stretching sheet," International Journal of Applied Mechanics and Engineering, vol. 11, pp. 647-654, 2006.

[36] A. J. Chamkha, A. M. Aly, M. A. Mansour, "Similarity solution for unsteady heat and mass transfer from a stretching surface embedded in a porous medium with suction/injection and chemical reaction effects," Chemical Engineering Communications, vol. 197 no. 6, pp. 846-858, DOI: 10.1080/00986440903359087, 2010.

[37] S. Mukhopadhyay, R. S. R. Gorla, "Unsteady MHD boundary layer flow of an upper convected Maxwell fluid past a stretching sheet with first order constructive/destructive chemical reaction," Journal of Naval Architecture and Marine Engineering, vol. 9 no. 2, pp. 123-133, DOI: 10.3329/jname.v9i2.12541, 2012.

[38] M. Khan, M. Azam, "Unsteady heat and mass transfer mechanisms in MHD Carreau nanofluid flow," Journal of Molecular Liquids, vol. 225, pp. 554-562, DOI: 10.1016/j.molliq.2016.11.107, 2017.