On a Bivariate Spectral Homotopy Analysis Method for Unsteady Mixed Convection Boundary Layer Flow, Heat, and Mass Transfer due to a Stretching Surface in a Rotating Fluid

A bivariate spectral homotopy analysis method (BSHAM) is extended to solutions of systems of nonlinear coupled partial differential equations (PDEs).Themethod has been used successfully to solve a nonlinear PDE and is now tested with systems.The method is based on a new idea of finding solutions that obey a rule of solution expression that is defined in terms of the bivariate Lagrange interpolation polynomials. The BSHAM is used to solve a system of coupled nonlinear partial differential equations modeling the unsteady mixed convection boundary layer flow, heat, and mass transfer due to a stretching surface in a rotating fluid, taking into consideration the effect of buoyancy forces. Convergence of the numerical solutions was monitored using the residual error of the PDEs. The effects of the flow parameters on the local skin-friction coefficient, the Nusselt number, and the Sherwood number were presented in graphs.


Introduction
Fluid flow dynamics due to a stretching surface find applications in the production of sheeting material including metal and polymer sheets.Examples include the cooling of an infinite metallic plate in a cooling bath, paper production, the boundary layer along material handling conveyers, the aerodynamic extrusion of plastic sheets, glass blowing, metal spinning, and drawing plastic films, to name just a few of these applications [1].The study of a two-dimensional stretching sheet flow was pioneered by Crane [2], who presented an exact analytical solution for the steady two-dimensional stretching sheet problem in a quiescent fluid.Due the important applications of the flow, many researchers have developed Crane's model to study different aspects of boundary layer flow due to a stretching surface.
The two-dimensional flow in the boundary layer induced by a rotating fluid past a stretching surface was first studied by Wang [3].Such flows have applications in geophysics and astrophysics, in solar physics involved in the sunspot development, and in self-cooled liquid metal blankets in fusion reactors when the container is being rotated [4].Rajeswari and Nath [5] extended Wang's model to study unsteady flow.Other reported studies include that of Nazar et al. [6], who carried out a numerical investigation of the induced unsteady flow due to a stretching surface in a rotating fluid, where the unsteadiness is caused by the suddenly stretched surface.The resulting equations were solved using the Keller-box method.
The homotopy analysis method was used by Tan and Liao [7] to find accurate analytic approximations of the partial differential equations modeling the unsteady viscous flows due to a suddenly stretching surface in a rotating fluid.Abbas et al. [4] carried out a numerical study to investigate the unsteady magnetohydrodynamic (MHD) boundary layer flow and heat transfer in an incompressible rotating viscous fluid over a stretching continuous sheet.The resulting governing partial differential equations were solved using the Kellerbox method.The spectral relaxation method was used to solve a coupled highly nonlinear system of partial differential equations due to an unsteady flow over a stretching surface in an incompressible rotating viscous fluid by Awad et al. [8].The effects of binary chemical reaction and Arrhenius activation energy were taken into consideration.
Govardhan et al. [9] considered the unsteady flow due to a stretching surface in a rotating fluid in the presence of porous media, where the unsteadiness is caused by the suddenly stretched surface.The transformed governing partial differential equations were solved numerically by using the Adams predictor-corrector method for some values of the physically governing parameters.The Keller-box method was used by Hafidzuddin and Nazar [10], to carry out a numerical study for the steady magnetohydrodynamics rotating boundary layer flow and heat transfer of a viscous fluid over a permeable shrinking sheet.
Motivated by the first success of the bivariate spectral homotopy analysis method (BSHAM) in finding the solution of a nonlinear PDE reported in [11], the aim of the study is to extend the application of the BSHAM to solutions of systems of coupled nonlinear PDEs.The test problem is that of a three-dimensional unsteady laminar viscous flow with heat and mass transfer due to a stretching sheet in a rotating fluid.The model is an extension of that considered by Nazar et al. [6] and Tan and Liao [7], to include heat, mass transfer, and buoyancy effects on the flow.Using similarity transformations that were first used by Williams and Rhyne [12], the unsteady Navier-Stokes equations were reduced into a system of coupled nonlinear partial differential equations.We would like to acknowledge the fact that the problem can be solved using other methods like the original homotopy analysis method (HAM) as in the works of Liao [13,14], but not with the type of linear operators used with the BSHAM.
The BSHAM is based on finding solutions that obey a rule of solution expression that is expressed in terms of bivariate Lagrange interpolation polynomials.The efficiency and accuracy of the method will be validated using residual errors of the differential equations.Simulations showing important flow characteristics such as the skin-friction coefficient, heat, and mass transfer rates will be carried out.

Governing Equations
We consider the unsteady boundary layer flows with heat and mass transfer, caused by a stretching surface at  = 0 in a rotating fluid.Let , V,  be the velocity components in the direction of Cartesian axes , , , respectively, with the axes rotating at an angular velocity Ω in the  direction.When  < 0, the surface rotates at an angular velocity Ω in the  direction so that the fluid is at rest relative to the surface.Initially (for  = 0), both fluid and plate are stationary and have constant temperature  ∞ and concentration  ∞ .At time  = 0, the surface at  = 0 is impulsively stretched in the -direction.For this time, the surface temperature and concentration vary from  ∞ to   (  >  ∞ ) and  ∞ to   (  >  ∞ ), respectively.Due to the Coriolis force, the fluid motion is three-dimensional and the corresponding governing equations are given as [6,7].
where  denotes the pressure,  is the density, ] is the kinematic viscosity, ∇ 2 is the three-dimensional Laplacian,  and  are the fluid temperature and concentration, respectively,  is the acceleration due to gravity,   is the thermal expansion coefficient,   is the concentration expansion coefficient,  is the thermal diffusivity, and  is the diffusion coefficient.The initial and boundary conditions are given as follows: < 0: for any , , .
≥ 0: → 0, Introducing the following similarity transformations [12] in (2)-( 6), gives rise to the coupled system of nonlinear PDEs with the corresponding boundary conditions In the above,  is an important dimensionless parameter signifying the relative importance of rotation rate to stretching rate,  is the buoyancy parameter,  is the concentration buoyancy parameter, Pr is the Prandtl number, and Sc the Schmidt number given by The important physical parameters for the type of boundary layer flow under consideration are the skin-friction coefficient, Nusselt number, and Sherwood number.From the velocity field, the shearing stress at the surface can be obtained, which in the nondimensional form (skin-friction coefficient) in the and -directions using ( 10) is given by where Re  =  2 /] is the local Reynolds number.Using the temperature field, the heat transfer coefficient at the surface can be obtained, which, in the nondimensional form, in terms of the Nusselt number, is given by From the concentration field, we are able to get the mass transfer coefficient at the plate, which, in the nondimensional form, in terms of the Sherwood number, is given by

The Bivariate Spectral Homotopy Analysis Method (BSHAM)
In this section, we give a description of the method of solution used to solve the governing equations ( 11)-( 14), the bivariate spectral homotopy analysis method (BSHAM).The method is a recent innovation presented by Motsa [11].It is an extension of the spectral homotopy analysis method (SHAM), whereby base functions used are defined in terms of the bivariate Lagrange interpolation polynomials.The first step of the solution algorithm is to reduce the nonlinear system of PDEs into a system of nonlinear ordinary differential equation (ODEs).This is achieved by applying the spectral collocation in the  direction.The resulting equations are then integrated using the spectral collocation method.Following Motsa [11], we approximate the exact solutions of (, ), (, ), (, ), and (, ) by Lagrange interpolation polynomials that interpolate the functions at the so-called Gauss-Lobatto collocation points defined as where  ∈ [−1, 1] and  = 2 − 1, so that the unknown functions (, ), (, ), (, ), and (, ) are approximated as with   () = (,   ),  = , ,  or  and   () are the characteristic Lagrange cardinal polynomials given by Equations ( 21) are then substituted in the governing nonlinear system of PDEs ( 11)-( 14) and collocation is applied; that is, the equations are required to be satisfied exactly at the collocation points   ,  = 0, 1, 2, . . .,   .An important step in the substitution process is the evaluation of the derivatives of   () with respect to .Several formulas exist for computing the derivatives if the collocation points are chosen to be Chebyshev-Gauss-Lobatto points.The values of the derivatives at the Chebyshev-Gauss-Lobatto points are computed as where  , (,  = 0, 1, . . .,   ) are entries of the standard Chebyshev differentiation matrix d = [ , ] of size (  + 1) × (  + 1) (see, e.g., [15,16]).Consequently, substituting (23) in ( 11)- (14) gives subject to the boundary conditions With the initial conditions at  = 0 (which corresponds to  =    ) known, note that the system (24) forms a system of   nonlinear coupled ODEs to be solved for   ,   ,   , and   ,  = 0, 1, . . .,   −1 subject to the corresponding boundary conditions.The next step is to apply the HAM basic principle of decoupling the system of equations into a sequence of linear systems.We refer interested readers about the details on Liao's HAM to [17,18].The decoupling process is made possible through the formation of the so-called zeroth-order deformation equation, defined as with corresponding boundary conditions In the above equations,  ∈ [0,1] is an embedding parameter and ℏ represents a nonzero convergence controlling auxiliary parameter.The initial approximate solutions of   (),   (),   (), and   () are given as  ,0 (),  ,0 (),  ,0 (), and  ,0 (), respectively, for  = 0, 1, 2, . . .,   − 1.The components L  and N  ,  =   ,   ,   or   , are the linear and nonlinear operators of the unknown functions defined by The initial approximate solutions are chosen to ensure that From the zeroth-order equations (26), it can be shown that as  varies from 0 to 1, the unknown variables   (; ),   (; ), Θ  (; ), and Φ  (; ) vary from their initial approximate solutions  ,0 (),  ,0 (),  ,0 (), and  ,0 () to the corresponding solutions   (),   (),   (), and   (), respectively, of the governing equations (24).Expanding   (; ),   (; ), Θ  (; ), and Φ  (; ) using Taylor series about  gives (32) Convergence of the series (32) depends on a proper choice of the auxiliary parameter ℏ.Next, higher order deformation equations are formed, where we are able to get the terms of the series from their solution,  , ,  , ,  , , and  , ,  ≥ 1.
The higher order deformation equations are constructed by differentiating the zeroth-order deformation equations ( 26) times with respect to , then dividing by !, and finally setting  = 0 given as follows: with boundary conditions From (33), At this stage, the Chebyshev spectral collocation method is applied independently in the  direction to solve for the initial approximate solutions  ,0 ,  ,0 ,  ,0 , and  ,0 and the higher order deformation equations giving rise to  , ,  , ,  , , and  , .In the spectral domain, the derivatives of the variables with respect to  are defined in terms of the Chebyshev differentiation matrix; for example, the derivative in  , will be given as The derivative formula (36) is a result of a linear transformation used to convert the interval  ∈ [0,   ] to  ∈ [−1, 1] where   is a sufficiently large finite value chosen to approximate the conditions at ∞.This gives rise to the matrix equations: where (39)

Results and Discussion
In this section, numerical results of the governing system of partial differential equations ( 11)-( 14) obtained using the bivariate spectral homotopy analysis method are presented.The results show convergence of the series solutions against the convergence controlling parameter, ℏ, residual error curves and variation of the skin-friction coefficient, and Nusselt and Sherwood numbers with the different flow parameters.The number of collocation points used was 60 and 14 in space () and time (), respectively, if not stated.The residual error measures the extent to which the numerical solutions approximate the true solution of the governing differential equations ( 11)-( 14).It is used in this study to monitor the accuracy of the numerical scheme.For each variable function, the residual error is defined by where Δ  , Δ  , Δ  , and Δ  represent the governing nonlinear PDEs ( 11), ( 12), (13), and ( 14), respectively, and F  , G  , Θ  , and Φ  are the BSHAM approximate solutions at the time collocation points   .Figures 1-4 present the residual error curves against the convergence controlling parameter ℏ, at various iteration levels.The residual errors are seen to decrease with increase in iterations which shows that the numerical scheme converges.From Figures 1-4, we are also able to pick the value of ℏ that will give optimal results.
Variation of the residual error with time  is presented in Figures 5-8.Again the accuracy of the results is seen to improve with improvement in iterations proving convergence of the scheme.The error is seen to be uniform for values of time.The effect of the different flow parameters on the local skin friction along the -direction,   (0, ), is shown in Figures 9-13. Figure 9 shows the effect of the rotating ratio parameter  and it can be seen that the local friction decreases with increase in .These findings are in agreement with those of [8,9].The effect of  is shown to increase the local skinfriction coefficient in Figure 10.This is caused by the fact that  is positive, thus assisting the flow enhancing the fluid velocity.A similar effect is observed in Figure 11 with  because positive values of  also enhance the fluid velocity, in turn increasing the local skin friction.The effects of the Prandtl number and the Schmidt number are to decrease the local skin-friction coefficient   (0, ) as shown in Figures 12 and  13, respectively.An increase in the Prandtl number implies an increase in kinematic viscosity, in turn decreasing the velocity of the fluid flow.The skin friction will decrease with decrease in the velocity profile.The same effect is experienced with Sc; that is, increasing Sc implies an increase in the kinematic viscosity.Figures 14-18 depict the effect of the flow parameters on the local skin-friction coefficient along the -direction,   (0, ).As expected, the skin friction along -direction will also decrease with increase in the rotating ratio parameter  (presented in Figure 14) as observed earlier with   (0, ) in Figure 9.
The effects of  and  on   (0, ) are depicted in Figures 15 and 16, respectively.Increasing either  or  enhances the flow along the -direction, thus decreasing that along  the -direction.The momentum boundary layer along direction decreases, in turn decreasing the local skin-friction coefficient along -direction,   (0, ).
Figures 17 and 18 display the effect of the Prandtl and Schmidt number on   (0, ), respectively.For both flow parameters, the effect is insignificant for the unsteady-state flow.A slight increase in   (0, ) with increase in both Pr and Sc is reported towards the steady-state flow.
Figure 19   decreases with increase in  with the effect being more momentous towards the steady-state flow.Increasing  decreases the thermal boundary layer thickness, as reported in [10].A thinner thermal boundary layer consecutively results in a decrease of the heat transfer rate at the surface.
The effects of  and  on the local Nusselt number are presented in Figures 20 and 21, respectively.Both buoyancy forces are seen to enhance the temperature field and the thermal boundary layer.This explains the increase in the heat transfer rate with increase in both  and .
The effect of the Prandtl number on the local Nusselt number is presented in Figure 22.The effect of Pr is to increase the heat transfer rate.Physically, this is caused by the fact that increasing Pr reduces the temperature field and the thermal boundary layer thickness, in turn increasing the heat  transfer rate.The increase with  is significantly pronounced for all values of .
Variation of −  (0, ) with the Schmidt number is depicted in Figure 23.The heat transfer rate decreases with increase in the Schmidt number.The findings are in agreement with those of Awad et al. [8].For a fixed value of , the local Nusselt number increases with Pr as displayed in Figure 22.These findings are consistent with those of Awad et al. [8], whereby it is shown that the thermal boundary layer decreases with increase in Pr, in turn reducing the heat transfer rate at the surface.The effect of the Schmidt number on the local Nusselt number is insignificant for the unsteady-state flow as displayed in Figure 23.A slight increase in −  (0, ) is observed towards the steady-state flow.

Conclusion
The paper extended the application of a bivariate spectral homotopy analysis method to solutions of systems of nonlinear partial differential equations.The test problem used was that of a three-dimensional unsteady laminar viscous flow with heat and mass transfer due to a stretching sheet in a rotating fluid.The transformed equations gave rise to a coupled nonlinear system of partial differential equations.Error analysis was performed in terms of ℏ curves of the residual error of the solutions.Computations were also carried out to analyze the effect of the different flow parameters on the local skin-friction coefficient, the Nusselt number, and the Sherwood number.The buoyancy forces were found to be   supporting the velocity of the flow in both and -directions in the process decreasing the local skin-friction coefficients in both directions.The buoyancy forces were also reported to increase both heat and mass transfer rates of the flow.Analysis of the other flow parameters proved to be consistent with similar findings reported in the literature.The study proved to be a success, and the bivariate spectral homotopy analysis method can be used to find solutions of coupled nonlinear  systems of partial differential equations that are valid for the whole space and time domain.
presents the effect of the rotating parameter  on the local Nusselt number.The local Nusselt number