The Spectral Homotopy Analysis Method Extended to Systems of Partial Differential Equations

and Applied Analysis 3 subject to the boundary conditions f (ξ, 0) = 0, f 󸀠 (ξ, 0) = 1 θ (ξ, 0) = 1, φ (ξ, 0) = 1, ξ ≥ 0, f 󸀠 (ξ,∞) = 0, θ (ξ,∞) = 0, φ (ξ,∞) = 0 ξ ≥ 0 (13) where the Rayleigh number Ra x , the nonlinear temperature parameter γ, the buoyancy parameterNr, the Prandtl number Pr, the Schmidt number Sc, the Peclet number Pe x , and the Lewis number Le are defined as Ra x = (1 − φ f∞ ) gβ 0 ρf ∞ ΔTx α m a , γ = β 1 (T w − T ∞ ) β 0 , σ = β 3 (φ w − φ ∞ ) β 2 , Nr = (ρ p − ρ ∞ ) (φ w − φ ∞ ) β 0 ρ f∞ (1 − φ ∞ ) (T w − T ∞ ) , Pr = ] α m , Sc = ] D B , Le = α m


Introduction
The current study serves to demonstrate the extension of the spectral homotopy analysis method (SHAM) to systems of nonlinear partial differential equations.The method, originally introduced by Motsa et al. [1,2], was used to solve nonlinear ordinary differential equations.A recent modification (see Motsa [3]) allows the method to be used to find solutions of nonlinear PDEs.A nonlinear system of partial differential equations that describe unsteady nonlinear convective flow caused by an impulsively stretched plate is used as the test problem in the study.The spectral quasilinearisation method (SQLM), an adaptation of the quasilinearisation method (Bellman and Kalaba [4]), is used as a benchmark to prove the accuracy of the spectral homotopy analysis method.
The test equations for this study are obtained by considering unsteady nonlinear convection flow over an impulsively stretched flat surface.Fluid flow over stretching surfaces is important in many practical applications such as extrusion of plastic sheets, paper production, glass blowing, metal spinning, and drawing plastic films.The quality of the final product depends on the rate of heat transfer at the stretching surface.In a previous study in this field, Kumari et al. [5] used the Keller box method and the Nakamura method to investigate the problem of heat transfer in the unsteady free convection flow over a continuous moving vertical sheet in an ambient fluid.Ishak et al. [6] investigated theoretically the unsteady mixed convection boundary layer flow and heat transfer due to a stretching vertical surface in a quiescent viscous and incompressible fluid.Further, Pop and Na [7] and Wang et al. [8] dealt with the unsteady boundary layer flow due to impulsive flow starting from rest of a stretching sheet in a viscous fluid.A numerical solution for unsteady mixed convection boundary layer nanofluid flow and heat transfer due to a stretching vertical sheet was presented by Mahdy [9].The equations were solved using the implicit finite difference method.Bachok et al. [10] studied the flow and heat transfer problem due to the unsteady, twodimensional laminar flow of a viscous nanofluid caused by a permeable stretching/shrinking sheet in a quiescent fluid.Numerical solutions of the transformed governing equations were obtained using a shooting method.Sharma et al. [11] used the finite element method to find numerical solutions of the flow and heat transfer problem over a stretching sheet immersed in a nanofluid, with velocity slip at the boundary.
In fluid flow problems involving heat transfer, it may be essential to consider a nonlinear relationship between density and temperature.Thermal stratification and heat released by the viscous dissipation triggers some changes in density gradients.To the best of our knowledge there

Governing Equations
Consider the problem of unsteady nonlinear convection of a fluid over a stretching flat plate.Initially ( = 0), both the fluid and stretching plate are kept at a constant temperature   and concentration  ∞ where   >  ∞ is for a heated plate and   <  ∞ corresponding to a cooled plate.We assume that at  = 0 the velocity of the stretching plate is   = , where  is a positive constant.From the Boussinesq approximation, density is related to temperature and the concentration by the equation In the case of thermal stratification and heat released by viscous dissipation, wall jet like profiles induce significant changes in density gradients, and the density depends on the temperature or temperature and concentration in a nonlinear form (see [12][13][14][15][16][17]): Karcher and Müller [18] used the formulation below to define the nonlinearity of the relationship between the density, the temperature, and the concentration: where  ∞ is the constant fluid density,  ∞ and  ∞ are the fluid temperature and solutal concentration, respectively,  0 and  2 are the coefficients of thermal and solutal expansion, and  1 denotes the nonlinear coefficient of thermal expansion.Partha [16] investigated the natural nonlinear convection in a non-Darcy porous medium using a temperatureconcentration-dependent density relation in the form With the usual Boussinesq and the boundary layer approximations, the governing equations are in the form subject to the boundary conditions The initial conditions are  < 0 :  = 0, V = 0,  =  w ,  =   , ∀, , where  and V are the velocity components along  and  directions, respectively,  and  are the local fluid temperature and solute concentration across the boundary layer, respectively, ] is the kinematic viscosity,  is the acceleration due to gravity,  0 and  1 are the thermal expansion coefficients,  2 and  3 are solutal expansion coefficients,   is the effective thermal diffusivity, and  is the mass diffusivity.We introduce the following nondimensional variables: The governing equations ( 6)-( 8) along with the boundary conditions ( 9) can be presented in the form Abstract and Applied Analysis 3 subject to the boundary conditions where the Rayleigh number Ra  , the nonlinear temperature parameter , the buoyancy parameter Nr, the Prandtl number Pr, the Schmidt number Sc, the Peclet number Pe  , and the Lewis number Le are defined as The skin friction and heat and mass transfer coefficients are described by the local skin friction coefficient, Nusselt number Nu, and the Sherwood number Sh defined as The nondimensional form of the skin friction coefficient, the reduced Nusselt number, and the reduced Sherwood number are

Method of Solution
The system of ( 12) was solved using the spectral homotopy analysis method.The spectral quasilinearisation method (SQLM) was used to validate the results.The spectral homotopy analysis method (SHAM) was first introduced for the solution of nonlinear ordinary differential equation.The method is a hybrid numerical scheme that combines the underlying ideas of the homotopy analysis method (HAM) and the Chebyshev spectral collocation method.The HAM scheme breaks down a nonlinear differential equation into infinitely many linear ordinary differential equations whose solutions are found analytically.The SHAM provides flexibility when solving linear ODEs through the use of the Chebyshev spectral collocation method.The first attempt to extend the application of the SHAM to PDEs was made by Motsa [3] who solved a nonlinear partial differential equation for an unsteady boundary-layer flow caused by an impulsively stretching plate.The blending of homotopy based methods with other methods has also been considered in [19,20] wherein the homotopy perturbation method was used in conjunction with Laplace transform to solve partial differential equations.In this work, we extend the application to a system of nonlinear PDEs.In the framework of the SHAM, the nonlinear equations are decomposed into their linear and nonlinear parts as follows: with L  and N  ( = , , ) representing the linear and nonlinear operators of the equations, respectively.According to the HAM description, the zeroth order deformation equations are formulated as subject to the boundary conditions The functions F, Θ, and Φ are unknown,  ∈ [0, 1] is the embedding parameter, Λ is the nonzero convergence controlling auxiliary parameter, and the functions  0 (, ),  0 (, ), and  0 (, ) are the initial solutions.From the zeroth order deformation equations above, it can be noted that, as  increases from 0 to 1, the functions F, Θ, and Φ vary from the initial approximations  0 (, ),  0 (, ), and  0 (, ), respectively, to the solutions (, ), (, ), and (, ) of the original equations (12).Using the Taylor series to expand F, Θ, and Φ about  yields The series (23) converge provided the auxiliary parameter Λ and the linear operators are properly chosen.The linear operators are often chosen to be the linear parts of the governing equations.The functions   ,   , and   are found from solutions of higher order deformation equations which are obtained by differentiating the zeroth order deformation equations ( 18)  times with respect to , dividing by ! and setting  = 0.This gives rise to the following equations: with boundary conditions Here the prime denotes differentiation with respect to  and where At this stage, in the HAM framework, the higher order differential equations (24) are solved analytically.However, here, it is not possible since the linear operators contain variable ( Initial approximate solutions are chosen to satisfy the rule of solution expression in (28), the corresponding boundary conditions, and the equations Substituting (28) in (29) and balancing terms of equal order in  give whose solutions are, respectively, The function erfc() is the standard complementary error function defined by Substituting (28) in the higher order deformation equations (24) gives rise to the iteration scheme with boundary conditions The system (33) is a linear sequence of ordinary differential equations that can be solved by any numerical method.In this study the equations are solved using the spectral methods whose underlying ideas will be sketched in Section 3.1.

The Spectral Quasilinearisation Method (SQLM).
In this section we present a spectral quasilinearisation method (SQLM) for solving systems of partial differential equations such as those in (12).Consider first an initial unsteady state solution corresponding to  = 0 and obtained by setting  = 0 in (12).This gives the equations Solving (35) analytically gives ) . ( The quasilinearisation method (QLM) is essentially a generalized Newton-Raphson method that was originally used by Bellman and Kalaba [4] for solving functional equations.The QLM scheme is derived by linearising the nonlinear component of the governing equations using Taylor series expansion with the assumption that the difference between the value of the unknown functions at the current iteration level (denoted by  + 1) and the value at the previous iteration level (denoted by ) is small.Applying the QLM on ( 12) leads to the following system of linear equations: Starting from the initial approximations (36), the iteration schemes (37) can be solved iteratively for  +1 (, ),  +1 (, ), and  +1 (, ) when  = 0, 1, 2, . ... In solving (37), we discretize the equation using the Chebyshev spectral collocation method [21] in the -direction and use an implicit finite difference method in the -direction.The spectral method can be implemented on the interval [−1, 1]; therefore it is necessary to transform the domain of the problem to the interval [−1, 1] before applying the spectral method.For the convenience of the numerical computations, the semi-infinite domain in the space direction is approximated by the truncated domain [0,  ∞ ], where  ∞ is a finite number selected to be large enough to represent the behaviour of the flow properties when  is very large.The transformation  =  ∞ ( + 1)/2 is used to map the interval [0,  ∞ ] to [−1, 1].The basic idea behind the spectral collocation method is the introduction of a differentiation matrix  which is used to estimate the derivatives of the unknown variables (), (), and () at the collocation points (grid points) as the matrix vector product.For example, the vector product for () will be as follows: where   +1 is the number of collocation points, D = 2/ ∞ , and is the vector function at the collocation points.Higher order derivatives are obtained as powers of D; that is, where  is the order of the derivative.We choose the Gauss-Lobatto collocation points to define the nodes in [−1, 1] as The matrix  is of size (  + 1) × (  + 1).The grid points on (, ) are defined as where   + 1 and   + 1 are the total numbers of grid points in the and -directions, respectively, and Δ is the spacing in the -direction.The finite difference scheme is applied with centering about a midpoint halfway between  +1 and   .This midpoint is defined as The derivatives with respect to  are defined in terms of the Chebyshev differentiation matrices.Applying the centering about  +(1/2) to any function, say (, ), and its associated derivative we obtain Thus, applying the spectral method to (37) with the finite difference in  gives where A  , B  (,  = 1, 2, 3) and   ( = 1, 2, 3) are (  + 1) × (  +1) matrices and (+1)×1 vectors, respectively, defined as

Results and Discussion
In this section we present numerical solutions to (12) computed using the spectral homotopy analysis method.Approximate solutions of the skin friction coefficient, the surface heat transfer, and the surface mass transfer rates at different values of flow parameters are presented and compared with the SQLM solutions.Velocity, temperature, and concentration profiles at different flow parameters are also presented.All SHAM results were generated using  = 30,  = 100, and Λ = −1 unless stated otherwise.These values were found to give accurate solutions after a numerical  As can be seen from the table, the results match perfectly well for the set accuracy level.The data from the table also reveals a steady increase of   (, 0), −  (, 0), and −  (, 0) with .
Table 2 shows a variation of   (, 0), −  (, 0), and −  (, 0) with selected flow parameters.It can be seen from the data in Table 2 that all , , Nr, and  contribute positively towards the skin friction coefficient, Nusselt number, and Sherwood number as they are seen to increase with the increase in the flow parameters.All SHAM results in the table are accurate for up to six decimal places.
Figure 1 shows the variation of the skin friction coefficient with  and .The skin friction coefficient is shown to be an increasing function of  for values of  = 0.1 and  = 0.3 while showing a decreasing function for  = 0.5.The skin friction coefficient is shown to decrease with the increase in  for any value of  and the decrease is more prominent for increasing values of .
The effect of  and Nr on the skin friction coefficient is presented in Figure 2. The skin friction coefficient is again seen to be a decreasing function of  for any value of  or Nr.Increasing  or Nr further decreases the skin friction coefficient.
Figure 3 shows the effect of the Prandtl number on the temperature profile.The temperature decreases with an  increase in Prandtl number.This is due to the fact that smaller values of Pr are corresponding to larger values of thermal conductivities and therefore heat is able to diffuse away from the stretching sheet.The effect of the Schmidt number Sc on the concentration profile is shown in Figure 4.It can be seen from the figure that, as the Schmidt number increases, the volume fraction boundary layer decreases across the stretching sheet.A higher Schmidt number implies a lower Brownian diffusion coefficient which will give rise to a shorter penetration depth for the concentration boundary layer.
In Figure 5, we present the solution of concentration profile for various values of the Lewis number (Le).Increasing the Lewis number decreases the concentration profile.Larger values of the Lewis number correspond to weaker molecular diffusivity and thinner boundary layer thickness.

Conclusions
The present study was carried out to extend the application of the spectral homotopy analysis method (SHAM) to systems of nonlinear partial differential equations.The method has been used successfully to find solutions of nonlinear ordinary differential equations.We solve the nonlinear system of partial differential equations governing the unsteady nonlinear convection flow over a stretching flat plate.The spectral quasilinearisation method (SQLM) was used to validate the SHAM results.Numerical values of the skin friction coefficient and the reduced Nusselt and Sherwood numbers were generated using both methods.The results matched well for the desired accuracy.Effects of selected flow parameters on the skin friction coefficient, temperature, and concentration profiles were also presented.The study

Figure 1 :Table 1 :
Figure 1: Effect of  and  on the skin friction coefficient.

Figure 2 :
Figure 2: Effect of  and Nr on the skin friction coefficient.
+1 +  1,   +1 +  2,  +1 +  3, = Pr (1 − )  +1  ,   +1 + 1,   +1 +  2,  +1 +  3, = Le Sc (1 − ) The same  = 30 and  = 100 were used for the computation of the SQLM solutions.The  and  in the tables represent the maximum th and th iteration required to produce converging results.Values of the skin friction coefficient and reduced Nusselt and Sherwood numbers at different values of  are presented in Table1.The table also shows a comparison of the SHAM and SQLM results.