# 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.

1. 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.

2. Governing Equations

We consider the unsteady boundary layer flows with heat and mass transfer, caused by a stretching surface at z = 0 in a rotating fluid. Let u, v, w be the velocity components in the direction of Cartesian axes x, y, z, respectively, with the axes rotating at an angular velocity [OMEGA] in the z direction. When t < 0, the surface rotates at an angular velocity [OMEGA] in the z direction so that the fluid is at rest relative to the surface. Initially (for t = 0), both fluid and plate are stationary and have constant temperature [T.sub.[infinity]] and concentration [C.sub.[infinity]]. At time t = 0, the surface at z = 0 is impulsively stretched in the x-direction. For this time, the surface temperature and concentration vary from [T.sub.[infinity]] to [T.sub.w] ([T.sub.w] > [T.sub.[infinity]]) and [C.sub.[infinity]] to [C.sub.w] ([C.sub.w] > [C.sub.[infinity]]), respectively. Due to the Coriolis force, the fluid motion is three-dimensional and the corresponding governing equations are given as [6, 7].

[mathematical expression not reproducible] (1)

[mathematical expression not reproducible], (2)

[mathematical expression not reproducible], (3)

[mathematical expression not reproducible], (4)

[mathematical expression not reproducible], (5)

[mathematical expression not reproducible], (6)

where p denotes the pressure, [rho] is the density, v is the kinematic viscosity, [[nabla].sup.2] is the three-dimensional Laplacian, T and C are the fluid temperature and concentration, respectively, g is the acceleration due to gravity, [[beta].sub.T] is the thermal expansion coefficient, [[beta].sub.C] is the concentration expansion coefficient, [alpha] is the thermal diffusivity, and D is the diffusion coefficient. The initial and boundary conditions are given as follows:

[mathematical expression not reproducible]. (7)

[mathematical expression not reproducible], (8)

[mathematical expression not reproducible]. (9)

Introducing the following similarity transformations [12] in (2)-(6),

[mathematical expression not reproducible], (10)

gives rise to the coupled system of nonlinear PDEs

[mathematical expression not reproducible], (11)

[mathematical expression not reproducible]. (12)

[mathematical expression not reproducible], (13)

[mathematical expression not reproducible], (14)

with the corresponding boundary conditions

[mathematical expression not reproducible]. (15)

In the above, [lambda] is an important dimensionless parameter signifying the relative importance of rotation rate to stretching rate, [gamma] is the buoyancy parameter, [sigma] is the concentration buoyancy parameter, Pr is the Prandtl number, and Sc the Schmidt number given by

[mathematical expression not reproducible]. (16)

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 x- and y-directions using (10) is given by

[mathematical expression not reproducible], (17)

where [Re.sub.x] = a[x.sup.2]/v 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

[Nu.sub.r] = [Nu.sub.x] [([xi][Re.sub.x]).sup.-1/2] = -[theta]' ([xi],0). (18)

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

[Sh.sub.r] = Sh [([xi][Re.sub.x]).sup.-1/2] = -[theta]' ([xi], 0). (19)

3. 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 [xi] direction. The resulting equations are then integrated using the spectral collocation method. Following Motsa [11], we approximate the exact solutions of f([eta], [xi]), g([eta], [xi]), [theta]([eta], [xi]), and ([eta], [xi]) by Lagrange interpolation polynomials that interpolate the functions at the so-called Gauss-Lobatto collocation points defined as

[[zeta].sub.i] = cos i[pi]/[N.sub.t], i = 0,1,2,...,[N.sub.t], (20)

where [zeta] [member of] [-1,1] and [zeta] = 2[xi] - 1, so that the unknown functions f([eta], [xi]), g([eta], [xi]), [theta]([eta], [xi]), and [phi]([eta],[xi]) are approximated as

[mathematical expression not reproducible], (21)

with [y.sub.j]([eta]) = y([eta], [[xi].sub.j]), y = f, g, [theta] or [phi] and [L.sub.j]([xi]) are the characteristic Lagrange cardinal polynomials given by

[mathematical expression not reproducible]. (22)

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 [[xi.sub.]i] = 0,1,2, ...,[N.sub.t]. An important step in the substitution process is the evaluation of the derivatives of [L.sub.j]([xi]) 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

[mathematical expression not reproducible], (23)

where [d.sub.i,j](i, j = 0,1,...,[N.sub.t]) are entries of the standard Chebyshev differentiation matrix d = [[d.sub.i,j]] of size ([N.sub.t] + 1) x ([N.sub.t] + 1) (see, e.g., [15,16]).

Consequently, substituting (23) in (11)-(14) gives

[mathematical expression not reproducible], (24)

subject to the boundary conditions

[mathematical expression not reproducible]. (25)

With the initial conditions at [xi] = 0 (which corresponds to [mathematical expression not reproducible]) known, note that the system (24) forms a system of [N.sub.t] nonlinear coupled ODEs to be solved for [f.sub.i], [g.sub.i], [[theta].sub.i], and [[phi].sub.i], i = 0,1,.... [N.sub.t] - 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

[mathematical expression not reproducible], (26)

with corresponding boundary conditions

[mathematical expression not reproducible]. (27)

In the above equations, q [member of] [0,1] is an embedding parameter and h represents a nonzero convergence controlling auxiliary parameter. The initial approximate solutions of [f.sub.i]([eta]), [g.sub.i]([eta]), [[theta].sub.i]([eta], and are given as [f.sub.i,0]([eta]), [g.sub.i,0]([eta]), [[theta].sub.i,0]([eta]), [[phi].sub.i,0]([eta]), respectively, for i = 0,1,2,...,[N.sub.t] - 1. The components [L.sub.p] and [N.sub.p], p = [f.sub.i], [g.sub.i], [[theta].sub.i] or [[phi].sub.i], are the linear and nonlinear operators of the unknown functions defined by

[mathematical expression not reproducible]. (28)

The initial approximate solutions are chosen to ensure that

[mathematical expression not reproducible]. (29)

From the zeroth-order equations (26), it can be shown that as q varies from 0 to 1, the unknown variables [F.sub.i]([eta]; q), [G.sub.i]([eta]; q), [[THETA].sub.i]([eta]; q), and [[PHI].sub.i]([eta]; q) vary from their initial approximate solutions [f.sub.i,0]([eta]), [g.sub.i,0]([eta]), [[theta].sub.i,0]([eta]), and to the corresponding solutions [f.sub.i]([eta]), [g.sub.i]([eta]), [[theta].sub.i]([eta]), and respectively, of the governing equations (24). Expanding [F.sub.i]([eta]; q), [G.sub.i]([eta]; q), [[THETA].sub.i]([eta]; q), and [[PHI].sub.i]([eta]; q) using Taylor series about q gives

[mathematical expression not reproducible], (30)

where

[mathematical expression not reproducible]. (31)

Therefore, since [F.sub.i]([eta];1) = [f.sub.i]([eta]), [G.sub.i]([eta];1) = [g.sub.i]([eta]), [[THETA].sub.i]([eta];1) = [[theta].sub.i]([eta]), [[PHI].sub.i]([eta];1) = [[phi].sub.i]([eta]) and [F.sub.i]([eta];0) = [f.sub.i,0]([eta], [xi]), [G.sub.i]([eta];0) = [g.sub.i,0([eta],] [xi]) [[THETA].sub.i]([eta], 0) = [[theta].sub.i,0]([eta], [xi]), [[PHI].sub.i]([eta], 0) = [[phi].sub.i,0]([eta], [xi]), we obtain

[mathematical expression not reproducible]. (32)

Convergence of the series (32) depends on a proper choice of the auxiliary parameter h. Next, higher order deformation equations are formed, where we are able to get the terms of the series from their solution, [f.sub.i,m], [g.sub.i,m], [d.sub.i,m], and [[phi].sub.i,m], m [greater than or equal to] 1. The higher order deformation equations are constructed by differentiating the zeroth-order deformation equations (26)

m times with respect to q, then dividing by ml, and finally setting q = 0 given as follows:

[mathematical expression not reproducible], (33)

with boundary conditions

[mathematical expression not reproducible]. (34)

From (33),

[mathematical expression not reproducible]. (35)

At this stage, the Chebyshev spectral collocation method is applied independently in the [eta] direction to solve for the initial approximate solutions [f.sub.i,0], [g.sub.i,0], [[theta].sub.i,0], and [[phi].sub.i,0] and the higher order deformation equations giving rise to [f.sub.i,m], [g.sub.i,m], [0.sub.i,m], and [[phi].sub.i,m]. In the spectral domain, the derivatives of the variables with respect to [eta] are defined in terms of the Chebyshev differentiation matrix; for example, the derivative in [f.sub.i,m] will be given as

[mathematical expression not reproducible], (36)

similar to the other variables. In (36), p is the order of the derivative, D = (2/[[eta].sub.e])[[D.sub.r,s]] (r, s = 0,1,2,..., [N.sub.x]) with [[D.sub.r,s]] being an ([N.sub.x] + 1) x ([N.sub.x] + 1) Chebyshev derivative matrix, and the vector [F.sub.i,m] is defined as

[mathematical expression not reproducible]. (37)

The derivative formula (36) is a result of a linear transformation used to convert the interval [eta] [member of] [0, [[eta].sub.e]] to [bar.[eta]] [member of] [-1,1] where [[eta].sub.e] is a sufficiently large finite value chosen to approximate the conditions at [infinity]. This gives rise to the matrix equations:

[mathematical expression not reproducible], (38)

where

[mathematical expression not reproducible]. (39)

4. 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 ([eta]) and time ([xi]), 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

[mathematical expression not reproducible],

where [[DELTA].sub.f], [[DELTA].sub.g], [[DELTA].sub.[theta]], and [[DELTA].sub.[phi]] represent the governing nonlinear PDEs (11), (12), (13), and (14), respectively, and [F.sub.i], [G.sub.i], [[THETA].sub.i], and [[PHI].sub.i] 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 h that will give optimal results.

Variation of the residual error with time [xi] 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 x-direction, f"(0, [xi]), is shown in Figures 9-13. Figure 9 shows the effect of the rotating ratio parameter [lambda] and it can be seen that the local friction decreases with increase in [lambda]. These findings are in agreement with those of [8, 9]. The effect of [gamma] is shown to increase the local skin-friction coefficient in Figure 10. This is caused by the fact that [gamma] is positive, thus assisting the flow enhancing the fluid velocity.

A similar effect is observed in Figure 11 with [sigma] because positive values of [sigma] 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 f"(0, [xi]) 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 y-direction, g'(0, [xi]). As expected, the skin friction along y-direction will also decrease with increase in the rotating ratio parameter A (presented in Figure 14) as observed earlier with f" (0, [xi]) in Figure 9.

The effects of [gamma] and [sigma] on g'(0, [xi]) are depicted in Figures 15 and 16, respectively. Increasing either [gamma] or [sigma] enhances the flow along the x-direction, thus decreasing that along the y-direction. The momentum boundary layer along y-direction decreases, in turn decreasing the local skin-friction coefficient along y-direction, g'(0, [xi]).

Figures 17 and 18 display the effect of the Prandtl and Schmidt number on g'(0,[xi]), respectively. For both flow parameters, the effect is insignificant for the unsteady-state flow. A slight increase in g'(0,[xi]) with increase in both Pr and Sc is reported towards the steady-state flow.

Figure 19 presents the effect of the rotating parameter A on the local Nusselt number. The local Nusselt number decreases with increase in A with the effect being more momentous towards the steady-state flow. Increasing A 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 y and a 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 [gamma] and [sigma].

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 [xi] is significantly pronounced for all values of [xi].

Variation of -0'(0, [xi]) 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 -[theta]'(0,[xi]) is observed towards the steady-state flow.

The effects of the flow parameters on the mass transfer rate are presented in Figures 24-28. The Sherwood number is shown to decrease with increase in [lambda] and Pr in Figures 24 and 27, respectively. The effects of [gamma], [sigma], and Sc are to enhance mass transfer rate, presented in Figures 25, 26, and 28, respectively. The results in the present study are consistent with those of a similar study carried out by [8].

5. 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 h 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 x-and y-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.

https://doi.org/10.1155/2017/5962073

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

References

[1] A. Ishak, R. Nazar, and I. Pop, "Unsteady mixed convection boundary layer flow due to a stretching vertical surface," The Arabian Journal for Science and Engineering, vol. 31, no. 2, pp. 165-182, 2006.

[2] L. J. Crane, "Flow past a stretching plate," Journal of Applied Mathematics and Physics (ZAMP), vol. 21, no. 4, pp. 645-647, 1970.

[3] C. Y. Wang, "Stretching a surface in a roating fluid," Journal of Applied Mathematics and Physics (ZAMP), vol. 39, pp. 177-185, 1988.

[4] Z. Abbas, T. Javed, M. Sajid, and N. Ali, "Unsteady MHD flow and heat transfer on a stretching sheet in a rotating fluid," Journal of the Taiwan Institute of Chemical Engineers, vol. 41, no. 6, pp. 644-650, 2010.

[5] V. Rajeswari and G. Nath, "Unsteady flow over a stretching surface in a rotating fluid," International Journal of Engineering Science, vol. 30, no. 6, pp. 747-756,1992.

[6] R. Nazar, N. Amin, and I. Pop, "Unsteady boundary layer flow due to a stretching surface in a rotating fluid," Mechanics Research Communications, vol. 31, no. 1, pp. 121-128, 2004.

[7] Y. Tan and S.-J. Liao, "Series solution of three-dimensional unsteady laminar viscous flow due to a stretching surface in a rotating fluid," Journal of Applied Mechanics, vol. 74, no. 5, pp. 1011-1018, 2007.

[8] F. G. Awad, S. Motsa, and M. Khumalo, "Heat and mass transfer in unsteady rotating fluid flow with binary chemical reaction and activation energy," PLoS ONE, vol. 9, no. 9, Article ID e107622, 2014.

[9] K. Govardhan, B. Balaswamy, and N. Kishan, "Unsteady boundary layer flow due to a stretching porous surface in a rotating fluid," Engineering Mechanics, vol. 21, no. 4, pp. 269-277, 2014.

[10] M. E. H. Hafidzuddin and R. Nazar, "Numerical solutions of MHD rotating flow and heat transfer over a permeable shrinking sheet," Science Asia, vol. 40, pp. 58-62, 2014.

[11] S. S. Motsa, "On an interpolation based spectral homotopy analysis method for PDE based unsteady boundary layer flows," Abstract and Applied Analysis, Article ID 848170, 7 pages, 2014.

[12] I. Williams III and T. B. Rhyne, "Boundary layer development on a wedge impulsively set into motion," SIAM Journal on Applied Mathematics, vol. 38, no. 2, pp. 215-224,1980.

[13] S. Liao, "Series solutions of unsteady boundary-layer flows over a stretching flat plate," Studies in Applied Mathematics, vol. 117, no. 3, pp. 239-263, 2006.

[14] S. Liao, "An analytic solution of unsteady boundary-layer flows caused by an impulsively stretching plate," Communications in Nonlinear Science and Numerical Simulation, vol. 11, no. 3, pp. 326-339, 2006.

[15] C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang, Spectral Methods in Fluid Dynamics, Springer, New York, NY, USA, 1988.

[16] L. N. Trefethen, Spectral Methods in MATLAB, vol. 10, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, Pa, USA, 2000.

[17] S. Liao, Beyond Perturbation: Introduction to Homotopy Analysis Method, Chapman & Hall/CRC Press, Boca Raton, Fla, USA, 2003.

[18] S.-J. Liao, Homotopy Analysis Method in Nonlinear Differential Equations, Springer, Berlin, Germany, 2012.

Sandile S. Motsa and Zodwa G. Makukula

University of Swaziland, Private Bag 4, Kwaluseni, Matsapha M201, Swaziland

Correspondence should be addressed to Zodwa G. Makukula; zmzmakukula@gmail.com

Received 1 February 2017; Accepted 16 April 2017; Published 9 May 2017

Caption: FIGURE 1: Residual error in f([eta], [xi]) against [??] at different iteration levels ([lambda] = 1, = [gamma] = 0.5, Pr = 3, Sc = 2, [xi] = 0.45).

Caption: FIGURE 2: Residual error in against h at different iteration levels ([lambda] = 1, [sigma] = [gamma] = 0.5, Pr = 3, Sc = 2, [xi] = 0.45).

Caption: FIGURE 3: Residual error in against h at different iteration levels ([lambda] = 1, [sigma] = [gamma] = 0.5, Pr = 3, Sc = 2, [xi] = 0.45).

Caption: FIGURE 4: Residual error in [phi]([eta],[xi]) against [??] at different iteration levels ([lambda] = 1, [sigma] = [gamma] = 0.5, Pr = 3, Sc = 2, [xi] = 0.45).

Caption: FIGURE 5: Residual error in g([eta],[xi]) against [xi] at different iteration levels ([lambda] = 1, [sigma] = [gamma] = 0.5, Pr = 3, Sc = 2, [??] = -0.5).

Caption: FIGURE 6: Residual error in against [xi] at different iteration levels ([lambda] = 1, [sigma] = [gamma] = 0.5, Pr = 3, Sc = 2, [??] = -0.5).

Caption: FIGURE 7: Residual error in 6[zeta], against [xi] at different iteration levels ([lambda] = 1, [sigma] = [gamma] = 0.5, Pr = 3, Sc = 2, [??] = -0.5).

Caption: FIGURE 8: Residual error in [phi]([eta], [xi]) against [xi], at different iteration levels ([lambda] = 1, [sigma] = [gamma] = 0.5, Pr = 3, Sc = 2, [??] = -0.5).

Caption: FIGURE 9: Effect of [gamma] on f"(0, [xi]) ([sigma] = [gamma] = 0.5, Pr = 3, Sc = 2, [??] = -0.5).

Caption: FIGURE 10: Effect of [gamma] on f"(0,[xi]) ([alpha] = 0.5 [lambda] = 0.5, Pr = 3, Sc = 2, [??] = -0.5).

Caption: FIGURE 11: Effect of [sigma] on f"(0,[xi]) ([lambda] = 1, [gamma] = 0.5, Pr = 3, Sc = 2, h = -0.5).

Caption: FIGURE 12: Effect of Pr on f"(0,[xi]) ([sigma] = [gamma] = 0.5, [lambda] = 1, Sc = 2, [??] = -0.5).

Caption: FIGURE 13: Effect of Sc on f"(0,[xi]) ([sigma] = [gamma] = 0.5, Pr = 3, [lambda] = 1, [??] = -0.5).

Caption: FIGURE 14: Effect of [lambda] on g'(0,[xi]) ([sigma] = [gamma] = 0.5, Pr =3, Sc = 2, [??] = -0.5).

Caption: FIGURE 15: Effect of [gamma] on g'(0,[xi]) ([sigma] = 0.5, [lambda] = 1, Pr = 3, Sc = 2, [??] = -0.5).

Caption: FIGURE 16: Effect of a on g'(0, [xi]), ([lambda] = [gamma] = 0.5, Pr =3, Sc = 2, [??] = -0.5).

Caption: FIGURE 17: Effect of Pr on g'(0,[xi]). ([sigma] = [gamma] = 0.5, [lambda] = 1, Sc = 2, [??] = -0.5).

Caption: FIGURE 18: Effect of Sc on g'(0,[xi]) ([sigma] = [gamma] = 0.5, Pr = 3, [lambda] = 1, [??] = -0.5).

Caption: FIGURE 19: Effect of [lambda] on [theta]'(0,[xi]), ([sigma] = [gamma] = 0.5, Pr =3, Sc = 2, [??] = -0.5).

Caption: FIGURE 20: Effect of [gamma] on [theta]'(0,[xi]) ([sigma] = 0.5, [lambda] = 1, Pr = 3, Sc = 2, [??] = -0.5).

Caption: FIGURE 21: Effect of [sigma] on [theta]'(0,[xi]) ([lambda] = 1, [gamma] = 0.5, Pr = 3, Sc = 2, [??] = -0.5).

Caption: FIGURE 22: Effect of Pr on [theta]'(0,[xi]) ([sigma] = [gamma] = 0.5, [lambda] = 1, Sc = 2, [??] = -0.5).

Caption: FIGURE 23: Effect of Sc on [theta]'(0,[xi]) ([sigma] = [gamma] = 0.5, Pr = 3, [lambda] = 1, [??] = -0.5).

Caption: FIGURE 24: Effect of [lambda] on ([sigma] = [gamma] = 0.5, Pr = 3, Sc = 2, [??] = -0.5).

Caption: FIGURE 25: Effect of [gamma] on [phi]'(0,[xi]) ([sigma] = 0.5, [lambda] = 1, Pr = 3, Sc = 2, [??] = -0.5).

Caption: FIGURE 26: Effect of [sigma] on [phi]'(0,[xi]) ([lambda] = 1, [gamma] = 0.5, Pr = 3, Sc = 2, [??] = -0.5).

Caption: FIGURE 27: Effect of Pr on [phi]'(0, ([sigma] = [gamma] = 0.5, [lambda] = 1, Sc = 2, [??] = -0.5).

Caption: FIGURE 28: Effect of Sc on [phi]'(0,[xi]) ([sigma] = [gamma] = 0.5, Pr = 3, [lambda] = 1, [??] = -0.5).