Double Diffusive Magnetohydrodynamic (MHD) Mixed Convective Slip Flow along a Radiating Moving Vertical Flat Plate with Convective Boundary Condition

In this study combined heat and mass transfer by mixed convective flow along a moving vertical flat plate with hydrodynamic slip and thermal convective boundary condition is investigated. Using similarity variables, the governing nonlinear partial differential equations are converted into a system of coupled nonlinear ordinary differential equations. The transformed equations are then solved using a semi-numerical/analytical method called the differential transform method and results are compared with numerical results. Close agreement is found between the present method and the numerical method. Effects of the controlling parameters, including convective heat transfer, magnetic field, buoyancy ratio, hydrodynamic slip, mixed convective, Prandtl number and Schmidt number are investigated on the dimensionless velocity, temperature and concentration profiles. In addition effects of different parameters on the skin friction factor, , local Nusselt number, , and local Sherwood number are shown and explained through tables.


Introduction
Nonlinear equations play an important role in applied mathematics, physics and issues related to engineering due to their role in describing many real world phenomena. The importance of obtaining exact or approximate solutions of nonlinear partial differential equations is still a big problem that compels scientists and engineers to seek different methods for exact or approximate solutions. A variety of numerical and analytical methods have been developed to obtain accurate approximate and analytic solutions for problems. Numerical methods give discontinuous points of a curve and thus it is very time consuming to obtain a complete curve of results. There are also some analytic techniques for nonlinear equations. Some of these analytic methods are Lyapunov's artificial small parameter method [1], d-expansion method [2], perturbation techniques [3,4], variational iteration method (VIM) [5,6] and homotopy analysis method (HAM) [7,8].
In recent years semi-numerical/analytical methods have become popular in magnetofluid dynamics research as they provide an alternative to purely numerical methods and require significantly less computational resources. One such method, the differential transform method (DTM) was first introduced by Zhou [9] in electrical circuit theory for solving both linear and nonlinear initial value problems. Developing this method for partial differential equations and obtaining closed form series solutions for linear and nonlinear initial value problems was carried out by Chen and Ho [10] in 1999, and Ayaz [11] applied DTM to the system of differential equations.
The significant advantage of the differential transform method over numerical methods is that it does not require linearization or discretization to be applied to nonlinear differential equations and therefore is not affected by the related errors. Also, DTM does not require a perturbation parameter and also the validity is independent of whether or not there exist small parameters in the considered equation. This method has been adapted in recent years and successfully applied to simulate many multi-physical transport phenomena problems including magnetic liquid film flows [12], mixed convection flow [13], micropolar convection [14].
Magnetohydrodynamics (MHD) is concerned with the mutual interaction of fluid flow and magnetic fields. The fluids being investigated must be electrically conducting and non-magnetic, which limits the fluids to liquid metals, hot ionized gases (plasmas) and strong electrolytes. The use of an external magnetic field is a very important issue in many industrial applications, especially as a mechanism to control material construction. Some important examples of magnetohydrodynamic flow of an electrically conducting fluid past a heated surface are MHD power generators, plasma studies, petroleum industries, cooling of nuclear reactors, the boundary-layer control in aerodynamics, and crystal growth [15,16]. The goal of the thermal treatment is to cool the material to a desirable temperature before spooling or removing it. As the high temperature material emerges from a furnace or a die, is exposed to the colder ambient, therefore transient conduction process accompanied by surface heat loss is initiated [17]. When high temperatures are encountered in the application areas, the thermal radiation effect becomes very important. High temperature plasmas, cooling of nuclear reactors, liquid metal fluids, and power generation systems are some important applications of radiative heat transfer from a surface plate to conductive fluids. There have been some studies that consider hydromagnetic radiative heat transfer flows. Spreiter and Rizzi [18] studied solar wind radiative magnetohydrodynamics. Nath et al. [19] obtained a set of similarity solutions for radiative-MHD stellar point explosion dynamics using shooting methods. Noor et al. [20] considered MHD free convection thermophoretic flow over a radiate isothermal inclined plate with heat source/sink effect.
Mixed convection flow and heat transfer over a continuously moving surface is applicable to many industrial fields such as hot rolling, paper production, wire drawing, glass fiber production, aerodynamic extrusion of plastic sheets, the boundary-layer along a liquid film, condensation process of metallic plate in a cooling bath and glass, and also in polymer industries [21]. The flow over a continuous material moving through a quiescent fluid is induced by the movement of the solid material and also by thermal buoyancy which will determine the momentum and thermal transport processes [22]. The first study of the flow field due to a surface moving with a constant velocity in a quiescent fluid was undertaken by Sakiadis [23]. Since then, other researchers investigated various aspects of mixed convection problems such as heat and (or) mass transfer, suction/injection, thermal radiation, MHD flow, porous media, slip flows, etc. [24,25].
In some situations such as the spreading of a liquid on a solid substrate, corner flow and the extrusion of polymer melts from a capillary tube, no slip conditions yield unrealistic behavior and must be replaced by slip conditions especially in applications of microfluidics and nanofluidics [26][27][28]. The difference between the fluid velocity at the wall and the velocity of the wall itself is directly proportional to the shear stress. The proportional factor is called the slip length. The corresponding slip boundary condition is u j j wall~l s Lu Ly , where l s is the slip length [29]. For gaseous flow the slip condition of the velocity and the jump condition of l Pr LT Ly , respectively where s v and s T are the tangential momentum coefficient and the temperature accommodation coefficient [30]. Some relevant papers on slip flows are Kim et al. [31], Martin and Boyd [32], Kuznetsov and Nield [33]. The above researchers restricted to either prescribed temperatures or heat flux at the wall or slip. The idea of using convective boundary conditions which are the generalization of isothermal and thermal slip boundary conditions was introduced by Aziz [34] and was followed by Magyari [35] who found an exact solution of Aziz's [34] problem in a compact integral form and Ishak [36] who extended the same problem for a permeable flat plate.
The goal of the present study is to develop similarity transformations via one parameter linear group of transformations and the corresponding similarity solutions for mixed convection flow of viscous incompressible fluid past a moving vertical flat plate with thermal convective and hydrodynamic slip boundary conditions and to solve the transformed coupled ordinary differential equations using the differential transform method. The effects of the Prandtl numberPr , the Schmidt number Sc, the mixed convective parameter Ri, the buoyancy ratio parameter N, the radiation parameter R, the magnetic field parameter M and the slip parameter a on the flow, heat and mass transfer characteristics are investigated numerically.
In section 2 the geometry of the flow along a moving vertical flat plate under consideration in this problem is modeled with hydrodynamic slip and thermal convective boundary condition assumptions. The system with two independent variables is reduced to one variable equations and a system of nonlinear ordinary differential equations is obtained using a linear group of transformations. Sections 3 and 4 provide methods of solution for the governing equations of the problem, the differential transform method and numerical solution, respectively. In section 5 physical reasons are illustrated for the behavior of the graphs and tables of the problem, and finally inferences are made and conclusions are drawn.  Figure 1 shows the geometry assumed in this study, along with the rectangular coordinates, x x and y y, and the corresponding velocity components, u u and v v (where i represents momentum, ii represents thermal and concentration boundary-layers and in general thermal and concentration boundary-layer thickness are not the same). The temperature of the ambient fluid is T ? , the unknown temperature of the plate is T w and the left surface of the plate is heated from a hot fluid of temperature T f (wT ? ) or is cooled from a cooled fluid (T f vT ? ) by the process of convection which yields a heat transfer variable coefficient h f ( x x). It is also assumed that the ambient fluid is of uniform concentration C ? , the unknown concentration of the plate is C w . A transverse magnetic field with variable strength B( x x)~B 0 x x 1=2 is applied parallel to the y y axis, where B 0 is the constant magnetic field [37]. Variable electric conductivity s~s 0 u u is assumed, where s 0 is the constant electric conductivity [37]. The magnetic Reynolds number is assumed to be small and hence the induced magnetic field can be neglected. Fluid properties are invariant except density, which is assumed to vary only in those changes that drive the flow (i.e., the Boussinesq approximation). Under the assumption of boundary-layer approximations, the governing boundarylayer equations in dimensional form are:

Nomenclature
subject to the boundary conditions: v v~0, u u~U w zU slip where T is the temperature, C is the concentration, u is the kinematic viscosity, k is the thermal conductivity, a is the thermal diffusivity, D is the mass diffusivity of species of the fluid, b T is the volumetric thermal coefficient, b C is the volumetric concentration coefficient, g is the acceleration due to gravity, s 1 is the Stefan-Boltzmann constant, k 1 is the Rosseland mean absorption coefficient, a~k r c p is the thermal diffusivity of the fluid, r is the density of the fluid, m is viscosity, h f ( x x) is the heat transfer coefficient and N 1 is the velocity slip factor.

Normalization
The following boundary-layer variables are introduced to express Eqs. (1)-(5) in dimensionless form: where Re~U 0 L=u is the Reynolds number based on the characteristic length L and characteristic velocity U 0 . The stream function y defined as u~L y Ly , v~{ Ly Lx is substituted into Eqs.
(2)-(5) to reduce the number of equations and number of dependent variables. Therefore the following three dimensionless equations are obtained: Ly Ly Pr Ly Ly Sc Ly Ly Here Pr~u =a is the Prandtl number, Sc~u =D is the Schmidt r is the magnetic field parameter, Re 2 is the mixed convective parameter [38], Riw0 is for aiding buoyancy flow and Riv0 is for opposing buoyancy flow and Ri~0 is for purely forced convective flow in which buoyancy effects are not present.
The boundary conditions take the following form:

Application of Linear Group Analysis and Similarity Equations
The transport Eqs. (7)-(10) form a highly coupled nonlinear boundary value problem. Numerical solutions of these equations are complicated and time consuming. In this section the linear group of transformations is imposed on the problem to combine the two independent variables (x,y) into a single independent variable g (similarity variable) and reduce Eqs. (7)-(10) into ordinary differential equations with corresponding boundary conditions. For this purpose all independent and dependent variables are scaled as: where A, a i (i~1,2,:::,6,7) are constants, the values of a i should be chosen such that the form of the Eqs. (7)-(10) is invariant under the transformations. Eqs. (7)-(10) will be invariant if a i are related by a 1~4 a 2 , a 3~3 a 2 , a 4~a5~0 , a 6~{ a 2 , a 7~a2 : ð12Þ It is clear from Eqs. (11) and (12) that This combination of variables is invariant under this group of transformations and consequently, is an absolute invariant which are functions having the same form before and after the transformation. This functional form is denoted using where g is the similarity independent variable. By the same argument, other absolute invariants are where f (g), h(g) and w(g) are the dimensionless velocity, temperature, concentration function, h f À Á 0 is the constant heat transfer coefficient, N 1 ð Þ 0 is the constant hydrodynamic slip factor. Substituting Eqs. (14) and (15) into Eqs. (7)-(9), the following ordinary differential equations are obtained subject to the boundary conditions where primes denote differentiation with respect to g. Here Bi~h f À Á 0 L=k Ra 1=2 is the Biot number and a~N 1 ð Þ 0 u Re 1=2 = L is the hydrodynamic slip parameter. Nondimensionalizing procedure is explained in detail in Appendix A.

Quantities of Engineering Interest
The physical parameters of interest in the present problem are the skin friction factor C f x x , local Nusselt number Nu x x and local Sherwood number Sh x x which may be determined, respectively by the following expressions: , Sh x x~ Using Eqs. (6), (14), (15), we have from Eq. (20)

Re
1=2 where Re x~U x x=u is the local Reynolds number.

The Differential Transform Method
DTM is employed to obtain semi-analytical/numerical solutions to the well-posed two-point boundary value problem defined by Eqs. (16)- (18) and conditions (19). DTM is an extremely strong technique in finding solutions to magnetohydrodynamic and complex material flow problems. It has also been used very effectively in conjunction with Padé approximants. Rashidi et al. [39] studied transient magnetohydrodynamic flow, heat transfer and entropy generation from a spinning disk using DTM-Padé. To provide a summary of the method, the transformation of the k th derivative of a function in one variable is considered which is defined as: where f g ð Þ is the original function and F k ð Þ is the transformed function. The differential inverse transform of F k ð Þ is: The concept of the differential transform is derived from a Taylor series expansion and in actual applications the function f g ð Þ is expressed by a finite series as follows: The value of m is decided by convergence of the series coefficients. The fundamental mathematical operations performed by DTM are listed in Table 1. Taking differential transforms of Eqs. (16)-(18), the following transformed equations are obtained: where F k ð Þ, H k ð Þ and W k ð Þ are the differential transform of f g ð Þ, h g ð Þ and w g ð Þ, respectively, and the transformed boundary conditions are: where a, b, c, d and " are constants which are computed from the boundary conditions. Here Páde approximants are applied to the problem to increase the convergence of a given series. As shown in Fig. 2, without using the Páde approximant, the different orders of the DTM solution, cannot satisfy the boundary conditions at infinity. Therefore, it is necessary to use DTM-Páde to provide an effective way to handle boundary value problems with boundary conditions at infinity. See Appendix B for a description of the Páde approximant method.

Numerical Solution
The system of nonlinear differential equations (16)-(18) is solved under the boundary conditions (19). The initial boundary conditions for f and h in (19) are unknown in comparison with the case of no-slip and no jump condition of the temperature boundary conditions. Hence, the solution of the system cannot

Transformed function
Original function :::(kzr)G(kzr) doi:10.1371/journal.pone.0109404.t001 proceed numerically using any standard integration routine. Here, following [40][41], a second order numerical technique is adopted. This technique combines the features of the finite difference method and the shooting method and is accurate because it uses central differences.
The semi-infinite integration domain g [ 0,? ½ Þ is replaced by a finite domain g [ 0,g ? ½ Þand g ? should be chosen sufficiently large so that the numerical solution closely approximates the terminal boundary conditions (19). Here a large enough finite value has been substituted for g ? , the numerical infinity, to ensure that the solutions are not affected by imposing the asymptotic conditions at a finite distance. The value of g ? has been kept invariant during the run of the program. Now a mesh is defined by g i~i h (i~0,1,:::,n), with h being the mesh size, and Eqs. (16)-(18) are discretized using central difference approximations for the derivatives, the following equations are obtained at the i th mesh point: Note that Eqs. (16)- (18), which are written at the i th mesh point, the first, second and third derivatives are approximated by central differences centered at i th mesh point. This scheme ensures that the accuracy of O h 2 À Á is preserved in the discretization. Eqs. (31)-(33) are term recurrence relations in F , h and w. So, in order to start the recursion, besides the values of F 0 , h 0 and w 0 , the values of F 1 , F 2 , h 1 and w 1 are also required. These values can be obtained by Taylor series expansion near g~0 for F , h and w, with initial assumptions for the dimensionless functions of F and h.
The values of F 0 ð Þ, F ' 0 ð Þ, h 0 ð Þ and w 0 ð Þ are given as boundary conditions in (19). The values of F ''' 0 ð Þ, h'' 0 ð Þ and w'' 0 ð Þ can be obtained directly from Eqs. (10)-(12) and using the initial assumptions. After obtaining the values of F 1 , h 1 and w 1 , at the next cycle F 2 , h 2 and w 2 are obtained. The order indicated above is followed for the subsequent cycles. The integration is carried out until the values of F , h and w are obtained at all the mesh points.
The three asymptotic boundary conditions (13) and (14) must be satisfied. Initial assumptions are found by applying a shooting method along with the fourth order Runge-Kutta method so as to fulfill the free boundary conditions at g~g ? in (19). The guesses can be improved by a suitable zero-finding algorithm, including Newton's method, Broyden's, etc. [42,43].

Results and Discussion
A linear group of transformations is used to reduce the two independent variables into one and reduce the governing equations into a system of nonlinear ordinary differential equations with associated boundary conditions. Equations (16)-(18) with boundary conditions (19) were solved analytically using the differential transform method and compared with numerical results. Figure 2 shows the results of comparison and great agreement is seen. Typically, the natural convection is negligible when Riv0:1, forced convection is negligible when Riw10, and neither is negligible when 0:1v Riv10. It is useful to note that usually the forced convection is large relative to natural convection except in the case of extremely low forced flow velocities. Here Ri w 0 is chosen to have aiding buoyancy flow.
For Biot number smaller than 0.1 the heat conduction inside the body is quicker than the heat convection away from its surface, and temperature gradients are negligible inside of it. Having a Biot number smaller than 0.1 labels a substance as thermally thin, and temperature can be assumed to be constant throughout the materials volume. The opposite is also true: A Biot number greater than 0.1 (a ''thermally thick'' substance) indicates that one cannot make this assumption, and more complicated heat transfer equations for ''transient heat conduction'' will be required to describe the time-varying and non-spatially-uniform temperature field within the material body. In this research the case of cooling of the plate Gr w 0 is assumed and also the value of 0.5 is chosen for the Biot number and the thermal radiation number is equal to 1 in all diagrams as the control parameters.
In Tables 2-6 effects of some of the parameters on the Skin friction factor, C f x x , local Nusselt number, Nu x x , and local Sherwood number, Sh x x are shown. According to Table 2, increasing the buoyancy ratio parameter, N, causes Re 1=2 x represents the relative magnitude of species buoyancy and thermal buoyancy forces. For N~1 these two forces are of the same magnitude. For Nw 0, the species buoyancy force is dominant and vice versa for N v 1. The weaker contribution of the thermal buoyancy force for N v 1 results in a depletion in the local Nusselt number. The stronger contribution of the species buoyancy force for N w 1 induces enhancement in the local Sherwood number. The momentum field is coupled to the concentration (species) field via the linear species buoyancy force, Ri N w, in the momentum Eq. (16). The species field is coupled to the momentum field via the nonlinear term, 3Scf w = =4, in the species diffusion Eq. (18). Both terms are assistive. The influence of N on both species diffusion and momentum diffusion is therefore very strong as observed in the very large magnitudes of skin friction and local Sherwood number in Table 2.
In Table 3, increasing the mixed convective parameter is observed to strongly increase skin friction and accelerate the boundary-layer flow. Only Ri.0 is considered corresponding to buoyancy-assisted flow. In the range 0:1 v Ri v 10 both free (natural) and forced convection modes contribute. For Ri.10, forced convection effects are negated and for Ri v 0:1 free convection effects vanish. Evidently at low Ri values, buoyancy has a lesser influence on the flow characteristics and skin friction is found to be negative (flow reversal). With Ri exceeding unity any flow reversal is eliminated and the flow is strongly accelerated. Stronger buoyancy effects therefore act to stabilize the flow and aid momentum development. This simultaneously encourages more efficient diffusion of heat and species resulting in a marked increase in heat and mass transfer rates (Re In Table 4 increasing the value of hydrodynamic slip parameter, a, decreases Re 1=2 x x C f x x but weakly increases Re Sh x x . Momentum slip is simulated in the wall velocity boundary condition given in Eq. (19). Increasing momentum slip causes a reduction in the penetration of the stagnant surface through the boundary-layer. This serves to enhance momentum boundary-layer thickness since the flow is decelerated with increasing slip so that skin friction is lowered.
Bi arises in the wall temperature gradient boundary condition in Eq. (19). As Bi increases from v 0:1 (thermally thin case) to w 0:1 (thermally thick case) the rate of thermal conduction heat transfer inside the plate becomes dramatically lower than the heat convection away from its surface, and temperature gradients are increased at the plate. The influence on the flow is to accelerate it. Skin friction is elevated and therefore momentum boundary-layer thickness decreased. As shown in Table 5, an increase in the value of Biot number, Bi, increases Re 1=2 Sh x x . Schmidt number is the ratio of viscous diffusion to molecular (species) diffusion. For Sc v 1, molecular diffusion rate exceeds the momentum diffusion rate and vice versa for Sc w1. Sub-unity values of Schmidt number will therefore result in a deceleration in the flow (reduced skin friction), which will also decrease thermal diffusion rates. Conversely mass transfer will be accentuated in the regime with increasing Sc values.
In Fig. 2 results using different orders of DTM and DTM-Páde are compared with those of the numerical method. Very good agreement is observed for DTM-Páde and the numerical method. In figures 3. a-3. f effects of different parameters are investigated on the flow regime and the following results are observed: N Fig. 3. a shows how the flow responds to change in the magnetic field. Increasing magnetic interaction number M from purely hydrodynamic case M~0 to higher values of M, gives rise to a strong deceleration in the flow. Presence of a magnetic field in an electrically conducting fluid introduces a Lorentz force which acts against the flow in the case that magnetic field is applied in the normal direction as considered in the present problem. The described type of resistive force tends to slow down the flow field.
N A positive rise in N induces an increase in the flow along the plate as seen in Fig. 3

. b.
N There is a clear decrease in the velocity values at the wall accompanying a rise in Prandtl number because the flow is decelerated. Fluids with higher Prandtl numbers will therefore possess higher viscosities (and lower thermal conductivities) which means that such fluids flow slower than lower Pr fluids (Fig. 3. c). As a result the velocity will be decreased substantially with increasing Prandtl number.   situation of a~0 and towards full slip as a??. In the limiting case the frictional resistance between the cooling liquid and the moving plate is eliminated, and the moving plate no longer imposes any motion of the cooling liquid.

Conclusions
In this study, combined heat and mass transfer of the flow along a moving vertical flat plate with hydrodynamic slip and thermal convective boundary condition was considered. In order to reduce the two independent variables into one and hence to reduce the governing equations into a system of nonlinear ordinary differential equations, a linear group of transformations was used. The obtained equations were solved analytically using the differential transform method. The results were verified by results taken from the numerical method and excellent agreement was observed. The effects of different parameters on the skin friction factor, C f x x , local Nusselt number, Nu x x , and local Sherwood number Sh x x were shown and explained through tables and also changes of dimensionless flow and heat and mass transfer rates due to changes in some parameters were analyzed and presented graphically.

APPENDIX A
The governing boundary-layer equations in dimensional form are: subject to the boundary conditions: v v~0, u u~U w zU slip~U0 , C~C w at y y~0, u u?0, T?T ? , C?C ? as y y??: Using the following boundary-layer variables: the following equations are obtained for Eqs. (a2)-(a4): After simplifying the equations and dividing Eq. (a7) by ð Þ =L, and using the definitions for Re~U 0 L=u, Ri~Gr r, R~16 s 1 T 3 ? =3 k 1 k, Pr~u =a, and Sc~u=D following equations are obtained: Pr The stream function y defined as u~Ly=Ly, v~{Ly=Lx is substituted into Eqs. (a10)-(a12) to reduce the number of equations and number of dependent variables, therefore the following three dimensionless equations are obtained: Ly Ly Pr Ly Ly Sc Ly Ly The boundary conditions take the following form: ?0, h?0, w?0 as y??: All independent and dependent variables are scaled as: x Ã~x A a 1 , y Ã~y A a 2 , y Ã~y A a 3 , where A, a i (i~1,2,:::,6,7) are constants. The values of a i should be chosen such that the form of the Eqs. (a13)-(a15) is invariant under the transformations by substituting the above variables: A a 1 za 2 A a 3 za 4 Sc Boundary conditions: A a 6 za 4 k Re 1=2 h Ã ,w Ã~1 at y~0 , L y Ã Ly Ã ?0, h Ã ?0, w Ã ?0 as y??: The construction of L=M ½ approximants involves only algebraic operations [44,45]. Each choice of L, degree of the numerator and M, degree of the denominator, leads to an approximant. How to direct the choice in order to obtain the best approximant is the major difficulty in applying the technique, which necessitates the need for a criterion for the choice depending on the shape of the solution. A criterion which has worked well here is the choice of L=M ½ approximants such that L~M . The approximants are constructed using MATHEMA-TICA software.