HOMOTOPY ANALYSIS FOR MHD HIEMENZ FLOW IN A POROUS MEDIUM WITH THERMAL RADIATION, VELOCITY AND THERMAL SLIPS EFFECTS

The present study deals with the two dimensional steady laminar forced MHD Hiemenz flow past a flat plate in a porous medium. The effects of thermal radiation and partial slips on the flow field have been investigated under the variable wall temperature condition of the plate. The governing equations have been transformed into a set of coupled non-linear ordinary differential equations (ODEs) by using suitable similarity transformations. These equations have been solved analytically by using homotopy analysis method (HAM). The effects of Prandtl number, suction/blowing parameter, permeability parameter, velocity slip parameter, radiation parameter, magnetic parameter, wall temperature exponent and thermal slip on velocity and temperature profiles have also been studied graphically. Our results have been validated with the help of results that have already been published.


INTRODUCTION
A point on the surface of an object in the flow field where the local fluid velocity is zero is called stagnation point. The stagnation point flow analysis plays a significant role in the study of numerous natural and industrial phenomena, because of its applications in the exploring of flows over the tips of submarines, tip of ships, and aircrafts. Also, this study is important in various engineering disciplines like hydrodynamic processes, cooling of nuclear reactors, cooling of electronic devices by fans, etc. On account of the afore-mentioned applications only, Hiemenz (1911) was the very first researcher to do a pioneering work so as to investigate the viscous fluid motion generated by a two-dimensional stagnation point flow over a flat plate. Due to this reason only, the planar laminar flow of an incompressible viscous fluid in a steady state close to a stagnation point is also called Hiemenz flow. For this flow in a plane, Hiemenz gave similarity solutions of the governing Navier-Stokes equations. Thereafter, this kind of study was carried forward by researchers like Eckert (1942) and Beard et al. (1964).
The process of suction/injection is of special significance with reference to the practical problems related to the boundary control applications e.g. film cooling, fibre coating, and coating of wires. Due to this reason, Schlichting and Bußmann (1943) were the first to analyze the effect of suction on the Hiemenz flow. This problem was further extended by Preston (1948). Ariel (1994) made the analysis of the same problem by way of considering uniform suction. As the magnetohydrodynamic (MHD) stagnation point flow problems are having their theoretical as well as practical applications in manufacturing processes like boundary * Currently in Department of Mathematics, Dr. Babasaheb Ambedkar Technological University, Lonere-402 103, Dist. Raigad(M.S.), India † Corresponding author. Email: snasreenbano@yahoo.in layer along material handling conveyers, blood flow problems, extrusion of plastic sheets, cooling of infinite metallic plate in cooling bath, etc., these problems have attracted in the recent past the attention of many a researchers such as Sparrow et al. (1963); Na (1979), and Ariel (1994). Because of the numerous engineering applications of porous media related heat transfer problems in geothermal energy recovery, crude oil extraction, thermal energy storage etc., Ingham and Pop (1998);Vafai (2005); Nield et al. (2006) and Raptis et al. (1982) analyzed hydromagnetic free convection flows through porous media. Thereafter, Takhar and Ram (1994) and Yih (1998) investigated the problem under different conditions. Yih studied the effect of uniform suction on the flow field by invoking the model of the porous medium developed by Raptis and Takhar (1987). The problem of Yih was revisited by Kechil and Hashim (2009); Rashidi and Erfani (2011);Yildirim and Sezer (2012); Ghasemi et al. (2012); Mabood and Khan (2014); Mohammadi et al. (2016) and Bhatti et al. (2016) by using various numerical and analytical methods. Further, the work of Rashidi and Erfani (2011) was extended by Kudenatti et al. (2017) by way of including the non-Darcy velocity in the governing boundary layer equations.
The fluids exhibiting wall slips properties are very important on account of their technological applications e.g. the polishing of artificial heart valves and internal cavities. So, in order to have a better understanding of the slips phenomena, many a researchers like Mooney (1931); Rao and Rajagopal (1999); Khaled and Vafai (2004); Wang (2002); Wang (2006) and Hayat et al. (2007) examined the effects of slips boundary conditions on fluid flows through different geometries.
Motivated by the work of Yih (1998) and others, we have attempted to investigate the two-dimensional steady laminar forced MHD Hiemenz flow against a flat plate with variable wall temperature in a porous medium with velocity and thermal slips conditions in the presence of thermal radiation using HAM. In this analysis, the system of the coupled nonlinear ODEs governing the afore-mentioned problem is solved by means of the Mathematica package BVPh 2.0 Zhao and Liao (2013) , which is valid for nonlinear boundary-value or eigen-value problems with boundary conditions at multiple points, and is governed by coupled nonlinear ODEs. It is based on the HAM Liao Liao (2003), which is an analytic approximation method for solving highly nonlinear problems. As an analytic tool to solve nonlinear differential equations, the HAM Liao (2003) has been used successfully to investigate a variety of nonlinear problems in science, engineering and finance. Unlike perturbation techniques, the HAM has nothing to do with small/large physical parameters, so that it is essentially a non-perturbation method. Besides, based on the homotopy in topology, the HAM enjoys great freedom in choosing auxiliary linear operators and initial guess. Unlike all other analytic techniques, the HAM especially provides us a convenient way to control and adjust the convergence of solution-series, so that the convergence of solution-series can be guaranteed.
The BVPh 2.0 Zhao and Liao (2013) is a HAM-based Mathematica package, providing us an easy-to-use tool for solving coupled nonlinear ODEs governing the boundary value or eigenvalue problems. BVPh 2.0 is a free-available online (http://numericaltank.sjtu.edu.cn/BVPh.htm) package. Different from numerical packages (such as BVP4c), it is based on the idea-computing numerically with functions instead of numbers. Unlike all the other packages, the convergence of results given by the BVPh 2.0 is especially guaranteed by means of the so-called convergence-control parameter in the frame of the HAM. It is well known that the convergencecontrol parameters control the convergence-rate of HAM solutions-series. As for now, two different techniques have been developed for the choice of the convergence-control parameters. One is called −curve Liao (2003) method and the other technique is minimum error method Liao (2012).
In the normal HAM, −curve method is used in which −curves are plotted to guess the convergence-control parameters on the lines of Guled and Singh (2016) and Xinhui et al. (2011). But in this analysis, the determination of the optimal value of the convergence-control parameter is not easily possible. As a matter of this fact, more number of approximations are required to get an accurate solution-series. On the contrary, the optimal homotopy analysis method (OHAM), on which the present mathematica package BVPh 2.0 is based, takes into account the minimum error method. In OHAM, the optimal value of the convergence-control parameter is obtained by minimizing the total average squared residual error at certain order of HAM approximation. So, we are able to find a more accurate and convergent solution-series after less number of approximations using this mathematica package which is based on OHAM. Moreover, in this package, the optimal value of the convergence-control parameter is obtained with the help of an in-built facility of 'GetOptiVar 'command.
As a consequence, the mathematica package BVPh 2.0 has recently been used successfully by numerous researchers like Farooq and Zhi-Liang (2013); Farooq et al. (2014); Farooq et al. (2015); and many others, to deal with various flow problems pertaining to nano-fluids for different geometries.

Fig. 1 Physical Model
Here we intend to study the two dimensional steady laminar forced convection in MHD Hiemenz flow at the stagnation region against a flat plate through a porous medium with thermal radiation. We choose the cartesian co-ordinate system in such a way that the x−axis is parallel and y−axis is normal to the sheet. The incompressible fluid is electrically conducting. A constant magnetic field of strength B0 is applied in ydirection. In the formulation, the wall temperature has been taken as variable and the surface mass flux as uniform. The induced magnetic field, the Hall effects and the viscous dissipation terms are neglected.
Following the Raptis and Takhar (1987) model for the porous medium and by introducing the boundary-layer approximation, the governing continuity, momentum and energy equations can be written as follows: where u and v are velocity components in x− and y−directions, respectively; P is the pressure; ρ is the fluid density; ν = µ ρ is the kinematic viscosity where µ is the coefficient of fluid viscosity; K is porosity parameter, σ is the electrical conductivity; B0 is the applied magnetic field along y−direction; T is the temperature of the fluid and the porous medium which are in local thermal equilibrium; α = k ρcp is the equivalent thermal diffusivity, where k is coefficient of thermal conductivity; cp is the specific heat at constant pressure, and qr is the radiative heat flux. The boundary conditions are defined as follows: Here vw is the uniform surface mass flux positive (i.e. vw > 0) for blowing and negative (i.e. vw < 0) for suction; u slip is velocity slip, which is proportional to the local wall shear stress and is given by N ν ∂u ∂y , where N is Navier's constant slip length; λ is the exponent of the wall temperature; D is temperature slip factor, and ue = ax is the free stream velocity where 'a'is a positive number.
With the help of free stream, Eq.
(2) takes the new form On solving the Eqs. (2) and (5), we obtain From Eq. (6), it is evident that the free stream velocity affects the flow of a fluid, and hence convection of heat will be affected considerably. That is why, this flow is of forced convection type. Using Rosseland approximation for radiation Brewster (1992), we can write in which σ1 depicts the Stefan-Boltzmann constant and k1 is the absorption coefficient. We presume that the temperature variation within the flow is such that T 4 may be expanded in a Taylor?s series.
On expanding T 4 about T∞ and neglecting higher order terms, we get Therefore, we have from Eqs. (7) and (8) ∂qr Using Eq. (9) in Eq. (3), one gets The stream function satisfying Eq.
(1) is defined as u = ∂ψ ∂y and v = − ∂ψ ∂x . The following similarity variables have been defined: where η is similarity variable and primes denote differentiation with respect to η. The following ordinary differential equations can be obtained by substituting (11) into Eqs. (6) and (10): with the boundary conditions which are transformed to where primes denote differentiation with respect to η; P r = ν α is the Prandtl number; Ω = ν Ka is the permeability parameter; M = σB 2 0 ρa is the magnetic parameter; s = − vw √ aα 1 is the mass transfer parameter where s > 0 for suction, s < 0 for blowing and s = 0 reflects the impermeable surface; β = N ν a α is the dimensionless velocity slip parameter; γ = D a α is the dimensionless thermal slip parame- is the radiation parameter and λ is the exponent of wall temperature.

HOMOTOPY TREATMENT
It was Liao (2012) who first of all in the year 2012 developed the HAMbased BVPh 1.0 Mathematica package to find the solutions of highly nonlinear ODEs governing boundary-value and eigenvalue problems having singularity, multi-point boundary conditions and having multiple solutions to their credit. But, Zhao and Liao (2013) developed a new version of this package to deal with the solutions of coupled nonlinear ODEs in the following year. To this package, Liao gave the name BVPh 2.0. As this package is more efficient than the earlier one, the authors in the present work have used it successfully to find the solutions of Eqs. (12) and (13); together with the boundary conditions (14).
For the detailed study of HAM and the package BVPh 2.0, one can refer to the paper of Liao (2013Liao ( ), (2012.
For the application of the package BVPh 2.0, one needs to define the governing equations of the problem in conjunction with the given boundary conditions. Thereafter, an auxiliary linear operator for each governing equation is chosen in such a way that it defines the equation-type of the corresponding high-order equations, and gives a proper guess for each unknown function. With this input, the package automatically provides us with the approximate analytic solutions at any order one would like. Here it is worth mentioning that there exists a convergence-control parameter for each governing equation. This parameter ensures the convergence of solution series. The optimal values of the convergence-control parameters are determined based on the minimum values of the residual squares of governing equations, and also of boundary conditions in some cases.
The HAM is based on the concept of homotopy in topology. This method converts a nonlinear problem into an infinite number of linear sub-problems without taking into consideration the magnitude of the physical parameters. For the given problem, we can have the following solutions Here fm(η) and θm(η) are evaluated by the high-order deformation equations that are governed by the chosen auxiliary linear operators. Based on the governing Eqs. (12) and (13) and the boundary conditions (14) at infinity, it is clear that the above solutions can also be written in the following forms: Here a m,k and b m,k are coefficients to be determined by the package BVPh 2.0. The Eqs. (17) and (18) represent the solution-expressions, which guide us in choosing the auxiliary linear operators, initial approximations and the auxiliary functions. Thus these solution expressions play a vital role in HAM. As per the solution expressions (17) and (18), we choose the following auxiliary linear operators: The above operators have the following properties: where c1, c2, c3, c4 and c5 are arbitrary constants. Within the framework of HAM, we can also choose the following initial approximations: and the following auxiliary functions The initial approximations given by Eqs. (22) and (23) must satisfy the conditions (14). This is the first and foremost requirement for the successful application of the package BVPh 2.0. The analytic solutions of the nonlinear ODEs (12) and (13) under the boundary conditions (14) can be gained automatically with the help of the auxiliary linear operators given by Eqs. (19) and (20), the initial approximations given by Eqs. (22) and (23) and by the auxiliary functions given by Eq. (24). The solutions f (η) and θ(η) given by the BVPh 2.0 contain two unknown convergence-control parameters c f 0 and c θ 0 . These convergence-control parameters are used to ensure the convergence of the series solutions. These convergence-control parameters play a significant role in the HAM. It is the convergence-control parameter that differentiates the HAM from the other analytic approximation methods. Liao and Zhao (2013) used the following average residual error at the m th -order of approximation in order to decrease the CPU time: Here N is an integer and N1, N2 are the nonlinear differential operators defined as: (27) The total error for the m th ?order of approximation is defined as follows: It may here be emphasized that all the boundary conditions are linear and are being satisfied exactly. Further, the optimal values of c f 0 and c θ 0 at the m th ?order of approximation are determined by the minimum of the total error E t m . This is done simply by using 'GetOptiVar'-a command of the package BVPh 2.0.    The system of nonlinear ODEs Eq. 12 and 13 with the boundary conditions 14 has been solved analytically by employing the HAM approach. For the service of the purpose, the Mathematica BVPh 2.0 package has been used successfully. Without any loss of generality, we have considered the following values of the governing parameters: P r = 7, Ω = M = 2, R = λ = s = 1, β = 0.3 and γ = 0.2; in order to obtain the corresponding optimal convergence-control parameters using the command "GetOptiVar"of the package BVPh 2.0. These values for upto the 6th approximations are listed in Table 1. Here, the total error E t m could be decreased to 8.12084 × 10 −4 with the help of the corresponding optimal convergence-control parameters h f = −0.198367 and h θ = −0.73158.

RESULTS AND DISCUSSION
As is clear from Table 1, the residual error of each governing equation decreases after using these optimal convergence-control parameters gained at the 6th-order of approximation. Thus, we succeed in getting the convergent analytic solution for the considered problem in the case as mentioned above.
On similar lines, the convergent analytic approximations corresponding to different physical parameters have also been obtained by means of BVPh 2.0. The Tables 3 and 4 depict the validity of the package for the considered two-dimensional MHD stagnation-point flow. In Table 3, the values of the local skin friction coefficient (f (0)), obtained in respect of different values of the magnetic parameter and suction/blowing parameters, are in excellent agreement with the corresponding values obtained by Bhatti et al. (2016), Kechil andHashim (2009) andYih (1998). Likewise, the values of the local Nusselt number for different values of the wall temperature exponent in case of suction/blowing are in precise agreement with the corresponding values of Kechil and Hashim (2009) and Yih (1998). This is obvious from Table 4. The Figs. 2-8 elaborate the characteristics of the flow and heat transfer for different values of Prandtl number (P r), permeability parameter (Ω), magnetic parameter (M ), suction/blowing parameter (s), velocity slip parameter (β), radiation parameter (R) , exponent of wall temperature (λ) and thermal slip parameter (γ).
The Fig. 2 illustrates the effects of Pr and s on the velocity and temperature profiles in the presence of slip parameters. The figure shows that the velocity profile decreases and temperature profile increases for increasing values of P r in both the cases of suction and blowing. When the thermal diffusivity dominates (i.e. when P r < 1), the velocity profiles come closer to the wall and there exists free stream velocity throughout the boundary-layer.
Furthermore, the Prandtl number (P r) plays a significant role in controlling the thicknesses of thermal and momemtum boundary-layers. The thermal boundary-layer thickness decreases in case of increasing Prandtl number and suction/blowing parameter. This happens simply due to the fact that the higher values of the Prandtl number correspond to the weaker thermal diffusivity. This results in a thinner thermal boundary layer. From Fig. 2, it is again evident that the temperature gradient at the surface at large values of Prandtl number is less as compared to that corresponding to the small values of the Prandtl number. This, in other words, means that there is less heat transfer at the surface if the blowing is strong enough. As a matter of this fact, we conclude that rate of heat transfer at the surface can be controlled by way of applying sufficiently strong blowing. The Fig. 3 describes graphically the effects of porosity of medium, and suction/blowing (s) on velocity and temperature profiles by considering the presence of velocity and thermal slips. It has been found that both the velocity and the thermal boundary-layer thicknesses show a decreasing trend with the increasing values of the permeability parameter (Ω). Further, as is obvious from Fig. 3, the sensitivity of the velocity profile towards permeability parameter is more as compared to temperature profile. In case of suction, the effect of permeability parameter on temperature profile is almost negligible.
The Fig. 4 illustrates the effects of magnetic parameter (M ) and suction/blowing (s) on the velocity and temperature profiles in the presence of velocity and thermal slips. The velocity profile shows an increasing trend in presence of slips corresponding to the increasing values of magnetic parameter. This occurs because of the fact that the magnetic force increases the fluid motion in the boundary-layer due to the presence of the term ( ue − u) in the momentum equation. This term remains positive in the boundary-layer region. Here, a drag-like force called the Lorentz force is produced by the electrically conducting fluid. The Lorentz force is always associated with the magnetic field that makes the boundarylayer thinner. On the contrary, the temperature profiles show a decreasing trend. Moreover, the effect of magnetic parameter on the temperature profiles is more intense in case of blowing than suction. From the Fig. 5, it is obvious that the velocity profiles exhibit an increasing trend while temperature profiles a decreasing trend with respect to the increasing values of suction/blowing parameter (s). This proves the profundity of the effect of suction/blowing parameter on the boundary-layer thickness. The suction reduces the thermal boundary-layer thickness while blowing thickens it. As a result, the process of suction can be used effectively for fast cooling. As the thermal boundary thickness increases with strong blowing, the heated fluid moves farther from the wall and forms an insulating layer of nearly the same temperature as that of the wall. This results into a decrease in the heat transfer rate from the wall, and hence leads to slower cooling.
The Fig. 6 makes the graphical representation of the effects of velocity slip and suction/blowing on momentum and thermal boundary-layer thicknesses. On the expected lines, the velocity is lower in case of no-slip as compared to that in case of slip. Further, the temperature of the fluid decreases from no-slip condition to slip condition. Thus, an increase in slip parameter, results into a decrease both in the momentum and thermal boundary-layer thicknesses. Also, the effect of slip on temperature profile in case of suction is less pronounced in comparison to that in case of blowing. The fluid velocity remains unaffected of the variations in the values of radiation parameter (R), wall temperature exponent (λ) and thermal slip (γ). This is simply because of the flow problem being uncoupled from the thermal problem.
From the Fig. 7, it is concluded that the thermal boundary layer thickness and the temperature distribution increase with the increasing values of the thermal radiation parameter. This is due to the fact that the divergence of the radiative heat flux (i.e. ∂qr ∂y ) increases along with the decreasing values of the Rosseland radiative absorptivity (k1). This, in turn, shows an increase in the rate of radiative heat transfer to the fluid. This causes the fluid temperature to increase. In view of this fact, the effect of radiation becomes more significant as R → ∞, and the radiation effect is negligible as R → 0.
In the Fig. 8, the effect of wall temperature exponent on thermal boundary layer has been investigated, and it has been found that the temperature shows a decreasing trend as the wall temperature exponent increases. In this case, the thermal boundary layer becomes thin. Further, from the Fig. 8, the increase of thermal slip (γ) results into a decrease in thermal boundary layer thickness.

CONCLUSION
The present study is carried out to analyze the two dimensional steady MHD laminar forced Hiemenz flow and heat transfer over a flat plate in a porous medium with varying wall temperature by considering the radiation effect and partial slip conditions. The analysis has been done with the help of BVPh 2.0 package based on the homotopy analysis method (HAM). The effects of the governing parameters P r, Ω, M, s, β, R, λ and γ on the velocity and temperature profiles are examined in details. The following significant conclusions are drawn from the analysis:  1. An increase in Prandtl number (P r) decreases the fluid velocity and hence enhances momemtum boundary layer thickness.
2. The velocity profile is an increasing function of the parameters Ω, M, s and β.
3. The fluid flow is decelerated by increase in P r.
4. The thermal boundary layer thickness exhibit increasing trend along with P r and R.
5. An increase in parameters Ω, M, s, λ and γ reduces the temperature profile.