Free Convection Nanofluid Flow in the Stagnation-Point Region of a Three-Dimensional Body

Abstract

Analytical results are presented for a steady three-dimensional free convection flow in the stagnation point region over a general curved isothermal surface placed in a nanofluid. The momentum equations in x- and y-directions, energy balance equation, and nanoparticle concentration equation are reduced to a set of four fully coupled nonlinear differential equations under appropriate similarity transformations. The well known technique optimal homotopy analysis method (OHAM) is used to obtain the exact solution explicitly, whose convergence is then checked in detail. Besides, the effects of the physical parameters, such as the Lewis number, the Brownian motion parameter, the thermophoresis parameter, and the buoyancy ratio on the profiles of velocities, temperature, and concentration, are studied and discussed. Furthermore the local skin friction coefficients in x- and y-directions, the local Nusselt number, and the local Sherwood number are examined for various values of the physical parameters.

1. Introduction

The stagnation point is defined as the point on the surface of objects in the flow field where the fluid is brought at rest by the object. The problem of flow and heat transfer at a general three-dimensional stagnation-point region has important applications in many manufacturing processes in petrochemical industries, the aerodynamic of plastic sheet, solar central receivers exposed to wind currents, and so forth. The boundary layer flows near a three-dimensional stagnation point of attachment on an isothermal surface have been examined several times in the past. Howarth [1], Davey and Schofield [2] discussed the solutions of boundary layer flow near the stagnation point on a general three-dimensional surface. Poots [3] formulated the boundary-layer equations for the free convection flow at three-dimensional lower stagnation point on a general curved isothermal surface. Banks [4] concluded that the three-dimensional solution can be exhibited at a two-dimensional stagnation point for a Prandtl number Pr=0.72. Further he investigated at infinitely large Prandtl number Pr=∞ the three-dimensional problem which can be reduced to the two-dimensional case. Banks [5] then extended the results of Poots [3] for negative values of 𝑐 (referred to as dual solution) corresponding to saddle points of attachment. Gorla et al. [6] explored dual solutions by a numerical procedure in the investigation of nonnewtonian power-law type fluids near a three-dimensional stagnation point. Ingham et al. [7] considered unsteady free convection flow near a three-dimensional stagnation point of attachment on an isothermal surface. Pop and Merkin [8] validated that the form of the three-dimensional free convection near a stagnation point on an isothermal surface depends on 𝑐=𝑏/𝑎, the ratio of the principle curvature of the body at a stagnation point. Merkin and Mahmood [9] investigated free convection boundary layer on a vertical plate with prescribed surface heat flux. Takhar et al. [10] studied unsteady free convection boundary layer flow in the forward stagnation-point region of a sphere, which is rotating with time dependent angular velocity in an ambient fluid. Slaouti et al. [11] described the momentum and the thermal boundary layer characteristics for unsteady free convection flow in the stagnation-point region of a heated three-dimensional body. Harris et al. [12] examined the free convection boundary layer flow near the lower stagnation point of a cylindrical body. Xu et al. [13] presented series solutions of unsteady free convection flow in the stagnation-point region of a three-dimensional body. Singh et al. [14] extended the work of Xu et al. [13] to the mass transfer case. Admon et al. [15] presented flow and heat transfer analysis for the unsteady free convection flow near the stagnation point of a three-dimensional body.

Nanoparticles are defined to be an object with at least one dimension smaller than 100 nanometer. When particle is such small and the area of the particle which got the surface is much higher than its weight, that makes particle having new and exciting properties. The addition of nanoparticles in a clear fluid like water and a solution where the suspension of nanoparticles can be stabilized; then, this setup is called nanofluids. Choi [16] is the first among the others to use the term nanofluids. The purpose of the nanofluids is to intensify the heat transfer properties and the development of the nanofluid is potentially for several applications such as cooling of electronics and cooling of data center where a computer eject a lot of heat, coolants, and so forth. The in-depth study on nanofluids can be found in the book by Das et al. [17], in the review paper by Wang and Mujumdar [18], Yang et al. [19], and Jahani et al. [20]. Kuzetsov and Nield [21] used Buogiorno [22] model to investigate a natural convection flow of a nanofluid over a vertical plate. They identified governing parameters for the transport process, namely, a buoyancy-ratio number 𝑁_𝑟, a thermophoresis parameter 𝑁_𝑡, a Brownian number 𝑁_𝑏, and a Lewis number Le. Recently the stagnation-point flow of a nanofluid has drawn considerable attention and a significant amount of research is in progress. Bachok et al. [23] studied the effects of nanoparticle volume fraction parameter and the ratio of gradients of velocities on the flow and heat transfer for the steady three-dimensional stagnation-point flow in a nanofluid. Bachok et al. [24] investigated two-dimensional stagnation-point flow of a nanofluid over a stretching/shrinking sheet. Alsaedi et al. [25] examined the influence of heat generation/absorption on the stagnation-point flow of nanofluids towards a linear stretching surface. Mustafa et al. [26] extended the work of Ibrahim et al. [27] with the inclusion of the combined effects of buoyancy force and convective heating.

In this paper the three-dimensional laminar free convection stagnation-point flow of a nanofluid on a general curved isothermal surface is considered. The conservation equations embodying the mass, the momentum, and the energy as well as the nanoparticle volume fraction are reduced to a set of four fully coupled nonlinear equations. Exact analytical approximations are obtained for the velocity, temperature, and concentration profile as well as for the local skin friction coefficients, the local Nusselt number, and the local Sherwood number by means of the OHAM technique [28]. A new kind of error estimation technique is defined and used to evaluate the convergence of the solution series. It should be noted that two different approaches can be used to choose the optimal values of the convergence-control parameters. One approach is the ₋ℎ-curve method which reflects the variation of ₋ℎ with a fixed unknown functional value (for an example, 𝑓^′′(0)) given by a certain order of HAM approximation. Other technique is the minimum error method [29], in which we can select the “best” value of the convergence-control parameters by the minimum errors of the whole equations. On the other hand, the explicit solutions are given based on the HAM technique (for details see Appendix B). Besides The effects of the physical parameters such as Le, 𝑁_𝑏, 𝑁_𝑡, and 𝑁_𝑟 on profiles of velocities, temperature, and concentration as well as the important quantities such as the local skin friction coefficients, the local Nusselt number, and the local Sherwood number are explored in detail.

2. Problem Formulation

Consider steady, viscous, laminar, incompressible, and free convection nanofluid flow in the stagnation-point region of a three-dimensional body. A Cartesian coordinate system (𝑥,𝑦,𝑧) is chosen with the origin 𝑂 at the stagnation point as shown in the Figure 1, where 𝑥- and 𝑦-coordinates are measured along the body surface while the 𝑧-coordinate is measured normal to the body surface. We assume that the flow is steady, the temperature and concentration are 𝑇_𝑤 and 𝐶_𝑤 at the surface, and the temperature and concentration for the ambient fluid is 𝑇_∞ and 𝐶_∞ and the gravitational acceleration is 𝑔, respectively. Under these assumptions the governing equations of continuity, momentum, energy, and nanoparticle volume fraction are as follows.

Figure 1 

Physical configuration.

Continuity Equation 𝜕𝑢+𝜕𝑥𝜕𝑣+𝜕𝑦𝜕𝑤𝜕𝑧=0. (1) Momentum equations for 𝑥- and 𝑦-components 𝑢𝜕𝑢𝜕𝑥+𝑣𝜕𝑢𝜕𝑦+𝑤𝜕𝑢1𝜕𝑧=−𝜌_𝑓∞𝜕𝑝𝜕𝜕𝑥+𝜈²𝑢𝜕𝑥²+𝜕²𝑢𝜕𝑦²+𝜕²𝑢𝜕𝑧²+1−𝐶_∞𝛽𝑇−𝑇_∞−𝜌_𝑝−𝜌_𝑓∞𝜌_𝑓∞𝐶−𝐶_∞𝑢𝑔𝑎𝑥,𝜕𝑣𝜕𝑥+𝑣𝜕𝑣𝜕𝑦+𝑤𝜕𝑣1𝜕𝑧=−𝜌_𝑓∞𝜕𝑝𝜕𝜕𝑦+𝜈²𝑣𝜕𝑥²+𝜕²𝑣𝜕𝑦²+𝜕²𝑣𝜕𝑧²+1−𝐶_∞𝛽𝑇−𝑇_∞−𝜌_𝑝−𝜌_𝑓∞𝜌_𝑓∞𝐶−𝐶_∞𝑔𝑏𝑦. (2) Thermal energy equation 𝑢𝜕𝑇𝜕𝑥+𝑣𝜕𝑇𝜕𝑦+𝑤𝜕𝑇𝜕𝜕𝑧=𝛼²𝑇𝜕𝑥²+𝜕²𝑇𝜕𝑦²+𝜕²𝑇𝜕𝑧²+𝜏𝐷_𝐵𝜕𝐶𝜕𝑥𝜕𝑇+𝜕𝑥𝜕𝐶𝜕𝑦𝜕𝑇+𝜕𝑦𝜕𝐶𝜕𝑧𝜕𝑇+𝐷𝜕𝑧_𝑇𝑇_∞𝜕𝑇𝜕𝑥²+𝜕𝑇𝜕𝑦²+𝜕𝑇𝜕𝑧². (3) Nanoparticle volume fraction equation 𝑢𝜕𝐶𝜕𝑥+𝑣𝜕𝐶𝜕𝑦+𝑤𝜕𝐶𝜕𝑧=𝐷_𝐵𝜕²𝐶𝜕𝑥²+𝜕²𝐶𝜕𝑦²+𝜕²𝐶𝜕𝑧²+𝐷_𝑇𝑇_∞𝜕²𝑇𝜕𝑥²+𝜕²𝑇𝜕𝑦²+𝜕²𝑇𝜕𝑧². (4) In the above equations, 𝑢,  𝑣, and  𝑤 are the 𝑥-component, 𝑦-component, and 𝑧-componen, respectively, of the fluid velocity with origin 𝑂; 𝜈 is the kinematic viscosity; 𝛽 is the coefficient of thermal expansion; 𝐶 is the nanoparticle volume fraction; 𝐷_𝐵 is the Brownian diffusion coefficient; 𝐷_𝑇 is the thermophoretic diffusion coefficient; 𝜏 is the heat capacity ratio. Subscripts 𝑝 and 𝑓, respectively, denote the nanoparticles and the base fluid. The two functions 𝑎 and 𝑏 are the principal curvatures of the body at the stagnation point. The corresponding boundary conditions are 𝑢=𝑣=𝑤=0,𝑇=𝑇_𝑤,𝐶=𝐶_𝑤at𝑧=0,𝑢⟶0,𝑣⟶0,𝑇⟶𝑇_∞,𝐶⟶𝐶_∞as𝜂⟶∞. (5) We now introduce the similarity transformations as follows: 𝜂=Gr^1/4𝑎𝑧,𝑢=𝜈𝑎²Gr^1/2𝑥𝑓^′,𝑣=𝜈𝑎²Gr^1/2𝑐𝑦𝑠^′,𝑤=−𝜈𝑎Gr^1/4𝑓+𝑠,𝜃𝜂=𝑇−𝑇_∞𝑇_𝑤−𝑇_∞,𝜙𝜂=𝐶−𝐶_∞𝐶_𝑤−𝐶_∞. (6) Substituting (6) into (1)–(5), the nondimensionalized governing equations take the forms 𝑓^′′′+𝑓+𝑐𝑠𝑓^′′−𝑓^′2𝑁+𝜃−^𝑟PrGr𝜙=0,𝑠^′′′+𝑓+𝑐𝑠𝑠^′′−𝑐𝑠^′2𝑁+𝜃−^𝑟PrGr𝜙=0,1𝜃Pr^′′+𝑓+𝑐𝑠𝜃^′+𝑁_𝑏𝜃Pr^′𝜙^′+𝑁_𝑡𝜃Pr^′2=0,𝜙^′′+𝑁_𝑡𝑁_𝑏𝜃^′′+LePr𝑓+𝑐𝑠𝜙^′=0. (7)(8)(9)(10) In the above equations, the parameters Pr, 𝑁_𝑏, 𝑁_𝑟, and 𝑁_𝑡 are the Prandtl number, the Brownian motion parameter, the buoyancy parameter, and the thermophoresis parameter, respectively, which are defined by 𝜈Pr=𝛼,𝑁_𝑏=𝜏𝐷_𝐵𝐶_𝑤−𝐶_∞𝛼,𝑁_𝑟=𝜌_𝑝−𝜌_𝑓∞𝐶_𝑤−𝐶_∞𝑔𝜇𝛼𝑎³,𝑁_𝑡𝐷=𝜏^𝑇𝑇_∞𝑇_𝑤−𝑇_∞𝛼, (11) where 𝛼 is the thermal conductivity, Gr and Le are the Grashof number and Lewis number, respectively. They are defined as Gr=1−𝐶_∞𝛽𝑔𝑇_𝑤−𝑇_∞𝜈²𝑎³𝛼,Le=𝐷_𝐵𝑏,𝑐=𝑎. (12) The boundary conditions in nondimensional form are 𝑓(0)=0,𝑓^′(0)=0,𝑓^′(∞)𝑠(0)=0,=0,𝑠^′(0)=0,𝑠^′(∞)𝜃(0)(∞)𝜙(0)(∞)=0,=1,𝜃=0,=1,𝜙=0. (13) The quantities of physical interest such as the local skin friction coefficients 𝐶_𝑓𝑥 and 𝐶_𝑓𝑦, along the 𝑥- and 𝑦-directions respectively and the local Nusselt number and the local Sherwood number, which are defined as 𝐶_𝑓𝑥=𝜏_𝑤𝑥𝜌𝑢²_𝑥,𝐶_𝑓𝑦=𝜏_𝑤𝑦𝜌𝑢²_𝑦,Nu_𝑥=𝑥𝑞𝐾𝑇_𝑤−𝑇_∞,Sh_𝑥=𝑥𝑗𝐷_𝐵𝐶_𝑤−𝐶_∞, (14) where 𝑢_𝑥=𝜈𝑎²Gr^1/2𝑥 and 𝑢_𝑦=𝜈𝑎²Gr^1/2𝑐𝑦 are the local reference velocities, and 𝜏_𝑤𝑥, 𝜏_𝑤𝑦, 𝑞_𝑤, and 𝑗_𝑤 are the local wall skin friction, the local heat flux, and the mass flux from the surface defined by 𝜏_𝑤𝑥=𝜇𝜕𝑢𝜕𝑧_𝑧=0,𝜏_𝑤𝑦=𝜇𝜕𝑣𝜕𝑧_𝑧=0,𝑞_𝑤=−𝐾𝜕𝑇𝜕𝑧_𝑧=0,𝑗_𝑤=−𝐷_𝐵𝜕𝐶𝜕𝑍_𝑧=0. (15) Substituting (6) in (15) the nondimensional skin friction coefficients, the Nusselt number, and the Sherwood number are as follows: Re_𝑥^1/2𝐶_𝑓𝑥=𝑓^′′(0),Re_𝑦^1/2𝐶_𝑓𝑦=𝑠^′′(0),Re_𝑥^−1/2Nu_𝑥=−𝜃^′(0),Re_𝑥^−1/2Sh_𝑥=−𝜙^′(0), (16) where Re_𝑥=𝑢_𝑥𝑥/𝜈 and Re_𝑦=𝑢_𝑦𝑥/𝜈 are the local Reynolds number in 𝑥- and 𝑦-directions, respectively. Note that the local Nusselt number and the local Sherwood number have the same values in both 𝑥- and 𝑦-directions, so we only list them in 𝑥-direction.

3. Convergence Criterion

3.1. Exact Residual Error

It should be emphasized that 𝑓_𝑘(𝜂), 𝑠_𝑘(𝜂), 𝜃_𝑘(𝜂), and 𝜙_𝑘(𝜂) as defined in Appendix A contain unknown convergence-control parameters 𝑐^𝑓₀, 𝑐^𝑠₀, 𝑐^𝜃₀, and 𝑐^𝜙₀, which determine the convergence region and rate of the homotopy-series solutions. To find out the optimal values of 𝑐^𝑓₀, 𝑐^𝑠₀, 𝑐^𝜃₀, and 𝑐^𝜙₀, first we define the exact squared residual errors at the 𝑘th-order of approximations as follows: Δ^𝑓_𝑘𝑐^𝑓₀,𝑐^𝑠₀,𝑐^𝜃₀,𝑐^𝜙₀=^∞₀𝒩_𝑓^𝑘_𝑖=0𝑓_𝑖,^𝑘_𝑖=0𝑠_𝑖,^𝑘_𝑖=0𝜃_1,𝑖,^𝑘_𝑖=0𝜙_𝑖²,Δ^𝑠_𝑘𝑐^𝑓₀,𝑐^𝑠₀,𝑐^𝜃₀,𝑐^𝜙₀=^∞₀𝒩_𝑠^𝑘_𝑖=0𝑓_𝑖,^𝑘_𝑖=0𝑠_𝑖,^𝑘_𝑖=0𝜃_1,𝑖,^𝑘_𝑖=0𝜙_𝑖²,Δ^𝜃_𝑘𝑐^𝑓₀,𝑐^𝑠₀,𝑐^𝜃₀,𝑐^𝜙₀=^∞₀𝒩_𝜃^𝑘_𝑖=0𝜃_𝑖,^𝑘_𝑖=0𝜙_𝑖²,Δ^𝜙_𝑘𝑐^𝑓₀,𝑐^𝑠₀,𝑐^𝜃₀,𝑐^𝜙₀=^∞₀𝒩_𝜙^𝑘_𝑖=0𝜃_𝑖,^𝑘_𝑖=0𝜙_𝑖²,Δ^𝑡_𝑘𝑐^𝑓₀,𝑐^𝑠₀,𝑐^𝜃₀,𝑐^𝜙₀=Δ^𝑓_𝑘𝑐^𝑓₀,𝑐^𝑠₀,𝑐^𝜃₀,𝑐^𝜙₀+Δ^𝑠_𝑘𝑐^𝑓₀,𝑐^𝑠₀,𝑐^𝜃₀,𝑐^𝜙₀+Δ^𝜃_𝑘𝑐^𝑓₀,𝑐^𝑠₀,𝑐^𝜃₀,𝑐^𝜙₀+Δ^𝜙_𝑘𝑐^𝑓₀,𝑐^𝑠₀,𝑐^𝜃₀,𝑐^𝜙₀. (17)(18) Note that Δ^𝑡_𝑘 contains all the unknown convergence-control parameters 𝑐^𝑓₀, 𝑐^𝑠₀, 𝑐^𝜃₀, and 𝑐^𝜙₀. Obviously, the more quickly Δ^𝑡_𝑘 decreases to zero, the faster the corresponding homotopy-series solutions converge. Thus, at the given order of approximation 𝑘, the corresponding optimal values of the convergence-control parameters are given by the minimum of Δ^𝑡_𝑘, corresponding to a set of four nonlinear algebraic equations 𝜕Δ^𝑡_𝑘𝜕𝑐^𝑓₀=0,𝜕Δ^𝑡_𝑘𝜕𝑐^𝑠₀=0,𝜕Δ^𝑡_𝑘𝜕𝑐^𝜃₀=0,𝜕Δ^𝑡_𝑘𝜕𝑐^𝜙₀=0. (19) However, Δ^𝑡_𝑘 defined by (18) takes too much CPU time to calculate even if the order of approximation is not very high.

3.2. Average Squared Residual Error

Thus, to greatly decrease the CPU time, it is used here the so-called average squared residual error defined by ℰ^𝑓_𝑘𝑐^𝑓₀,𝑐^𝑠₀,𝑐^𝜃₀,𝑐^𝜙₀=1𝑁+1^𝑁_𝑗=0𝒩_𝑓^𝑘_𝑖=0𝑓_𝑖,^𝑘_𝑖=0𝑠_𝑖,^𝑘_𝑖=0𝜃_𝑖,^𝑘_𝑖=0𝜙_{𝑖𝜂=𝑗𝛿𝜂}²,ℰ^𝑠_𝑘𝑐^𝑓₀,𝑐^𝑠₀,𝑐^𝜃₀,𝑐^𝜙₀=1𝑁+1^𝑁_𝑗=0𝒩_𝑠^𝑘_𝑖=0𝑓_𝑖,^𝑘_𝑖=0𝑠_𝑖,^𝑘_𝑖=0𝜃_𝑖,^𝑘_𝑖=0𝜙_{𝑖𝜂=𝑗𝛿𝜂}²,ℰ^𝜃_𝑘𝑐^𝑓₀,𝑐^𝑠₀,𝑐^𝜃₀,𝑐^𝜙₀=1𝑁+1^𝑁_𝑗=0𝒩_𝜃^𝑘_𝑖=0𝜃_𝑖,^𝑘_𝑖=0𝜙_{𝑖𝜂=𝑗𝛿𝜂}²,ℰ^𝜙_𝑘𝑐^𝑓₀,𝑐^𝑠₀,𝑐^𝜃₀,𝑐^𝜙₀=1𝑁+1^𝑁_𝑗=0𝒩_𝜙^𝑘_𝑖=0𝜃_𝑖,^𝑘_𝑖=0𝜙_{𝑖𝜂=𝑗𝛿𝜂}², (20) where 𝑁 is an integer. It is reasonable to set 𝑁=80 and 𝛿𝜂=0.1. we define ℰ^𝑡_𝑘𝑐^𝑓₀,𝑐^𝑠₀,𝑐^𝜃₀,𝑐^𝜙₀=ℰ^𝑓_𝑘𝑐^𝑓₀,𝑐^𝑠₀,𝑐^𝜃₀,𝑐^𝜙₀+ℰ^𝑠_𝑘𝑐^𝑓₀,𝑐^𝑠₀,𝑐^𝜃₀,𝑐^𝜙₀+ℰ^𝜃_𝑘𝑐^𝑓₀,𝑐^𝑠₀,𝑐^𝜃₀,𝑐^𝜙₀+ℰ^𝜙_𝑘𝑐^𝑓₀,𝑐^𝑠₀,𝑐^𝜃₀,𝑐^𝜙₀, (21) where ℰ^𝑡_𝑘 is the total average squared residual error. Obviously, the smaller the value of ℰ^𝑡_𝑘(𝑐^𝑓₀,𝑐^𝑠₀,𝑐^𝜃₀,𝑐^𝜙₀) for given iteration 𝑘, the better the approximation. The OHAM series solutions which greatly depend on optimal values of convergence-control parameters can be obtained accurately by the proper choice of convergence-control parameters.

3.3. Choice of Convergence-Control Parameters

In this subsection we will explain the procedure to obtain the optimal convergence-control parameters. The optimal values for the convergence-control parameters can be obtained by the minimum of total average squared residual error. The total average squared residual error ℰ^𝑡_𝑘(𝑐^𝑓₀,𝑐^𝑠₀,𝑐^𝜃₀,𝑐^𝜙₀) is minimized by using symbolic computation software 𝑀𝑎𝑡ℎ𝑒𝑚𝑎𝑡𝑖𝑐𝑎. It is clear the smaller the value of total averaged squared residual error for a given 𝑘 the better the approximation. We directly employ the command Minimize at any order of iteration, say “𝑘,” to obtain the corresponding global optimal convergence-control parameters. Let us assume Pr=Le=Gr=1.0, 𝑐=0.5, 𝑁_𝑏=𝑁_𝑡=𝑁_𝑟=0.1, and the convergence-control parameters are unknown. The corresponding squared residual errors at different orders of iterations are shown in Table 1. It can be seen as we increase the order that the total average squared residual error is decreasing. Physically the convergence-control parameters are just the artificial parameters introduced in the high order deformation equations to obtain the convergent results. Since the convergent-control parameters are unknown therefore the CPU takes too much computation time. Therefore it is reasonable to choose optimal values of the convergence-control parameters corresponding to 5-order OHAM iteration. Using fifth order OHAM iteration the OHAM approximations at different orders of iterations are shown in Table 2. The results in Table 2 demonstrate the reliability of OHAM series solutions.

Table 1 

Global optimal convergence control parameters at different orders of approximations for Pr=Gr=Le=1.0, 𝑁_𝑏=𝑁_𝑡=𝑁_𝑟=0.1.

Table 2 

Averaged squared residual errors using 𝑐^𝑓₀=−1.56, 𝑐^𝑠₀=−1.86, 𝑐^𝜃₀=−1.17, and 𝑐^𝜙₀=−1.36 where Pr=Gr=Le=1.0, 𝑁_𝑏=𝑁_𝑡=𝑁_𝑟=0.1, and 𝑐=0.5.

4. Results and Discussion

Analytical approximations with the high precision for velocity profile along 𝑥- and 𝑦-directions, temperature distribution, and concentration profile, the local skin friction coefficients along 𝑥- and 𝑦-directions, the local Nusselt number, and the local sherwood number are plotted in Figures 2–13 for various values of the physical parameters Le, 𝑁_𝑟, 𝑁_𝑡, 𝑁_𝑏, and 𝑐.

Figure 2 

Graphs of skin friction coefficients, Nusselt and Sherwood numbers. 𝑁𝑏=𝑁𝑡=𝑁𝑟=0.1 and Le=1.0.

Figure 3 

Graphs of skin friction coefficient along 𝑥- and 𝑦-direction for varying 𝑁𝑟.

Figure 4 

Graphs of skin friction coefficient along 𝑥- and 𝑦-direction for 𝑁𝑏=0.2 and varying 𝑁𝑡.

Figure 5 

Graphs of skin friction coefficient along 𝑥- and 𝑦-direction for 𝑁𝑡=0.2 and varying 𝑁𝑏.

Figure 6 

Graphs of 𝑓𝜂(𝜂) and 𝑠𝜂(𝜂) for different Le.

Figure 7 

Graphs of 𝜃𝜂(𝜂) and 𝜙𝜂(𝜂) for different Le.

Figure 8 

Graphs of 𝑓𝜂(𝜂) and 𝑠𝜂(𝜂) for different 𝑁𝑟.

Figure 9 

Graphs of 𝜃(𝜂) and 𝜙(𝜂) for different 𝑁𝑟.

Figure 10 

Graphs of 𝑓𝜂(𝜂) and 𝑠𝜂(𝜂) for different 𝑁𝑡.

Figure 11 

Graphs of 𝜃(𝜂) and 𝜙𝜂 for different 𝑁𝑡.

Figure 12 

Graphs of 𝑓𝜂(𝜂) and 𝑆𝜂(𝜂) for different 𝑁𝑏.

Figure 13 

Graphs of 𝜃(𝜂) and 𝜙(𝜂) for different 𝑁𝑏.

Figure 2 depicts the graphs of the local skin friction coefficients, the local Nusselt number, and the local Sherwood number for various values of 𝑐 keeping the other parameters such as Pr=Gr=Le=1.0, 𝑁_𝑏=𝑁_𝑡=𝑁_𝑟=0.1. Since we have assumed that 𝑎 and 𝑏 are positive and most shapes of practical interest lies between cylinder (𝑐=0) and sphere (𝑐=1) therefore it is reasonable to consider the values of 𝑐 between 0 and 1 however the current method can also be applied to generate the solutions for (1≤𝑐≤∞). The increase in the values of 𝑐 increases parameter 𝑏 (the curvature of the body in the 𝑥=0 plane) and decreases the parameter 𝑎 (the curvature of the body in the 𝑦=0 plane) and also the impact of the increase in the values of 𝑐 is to decrease the shear stress because of more streamlined flow. It can be seen in Figure 2 that the skin friction coefficients attains their maximum values in the stagnation point under the boundary layer effects and the increase in the value of 𝑐 reduces the skin friction. Since the increase in the values 𝑐 is due to the increase in the parameter 𝑏 hence the decrease in the skin friction is more pronounced in 𝑦-direction. The analysis of the parameter 𝑐 on the heat and mass transfer shows the increase in the values of 𝑐 will increase the local Nusselt number and the local Sherwood number.

Figure 3 reveals the effects of 𝑁_𝑟 and 𝑁_𝑡 on the local skin friction for both buoyancy assisted Gr=1.0 and buoyancy opposing flow Gr=−1.0. For buoyancy assisted flow the temperature at the stagnation point 𝑇_𝑤 is greater than the ambient fluid temperature 𝑇_∞, that is, 𝑇_𝑤>𝑇_∞. The increase in the values of 𝑁_𝑟 will increase the density and nanoparticles mass differences; thus, the fluid becomes more dense which opposes the transport phenomenon. It can be seen in Figure 3 the skin friction coefficient decreases with the increase in the values of 𝑁_𝑟. For buoyancy opposing case the ambient fluid temperature is higher, that is, 𝑇_∞>𝑇_𝑤 the increase in the values of 𝑁_𝑟 will result in the decrease in the density and mass fraction nanoparticles differences hence the fluid is less dense. For buoyancy opposing case this can be seen in Figure 3 the increase in the values of 𝑁_𝑟 increases the skin friction coefficients. Further it can be seen for buoyancy assisted flow the increase in the values of 𝑁_𝑡 increases the temperature difference which increases the density difference, the ambient fluid is more dense hence the skin friction along 𝑥- and 𝑦-directions decreases while quite the opposite is true for buoyancy opposing flow. Figure 4 describes the effects of 𝑁_𝑡 on the local skin friction coefficients. For buoyancy assisted flow the increase in the value of 𝑁_𝑡 decreases the local skin friction in both directions and for buoyancy opposing flow the increase in the value of 𝑁_𝑡 increases the local skin friction. It is evident from Figures 3 and 4 the effects of buoyancy parameter 𝑁_𝑟 are more prominent than thermophoresis parameter. Figure 5 demonstrates that for both buoyancy assisted and buoyancy opposing flows the increase in the values of 𝑁_𝑏 increases the skin friction in both directions.

Figure 6 depicts the general trend of velocity profiles when 𝑇_𝑤>𝑇_∞ and 𝜌_𝑝>𝜌_𝑓∞. It can be seen that velocity is zero at the stagnation point. Since the temperature is higher slightly away from the wall in the ambient fluid region therefore the particles moves from cold to hot region hence the velocity tends to increase. As it goes further away the temperature difference reduces which results in the decay of the buoyancy force and flow starts decreasing. Eventually away from the wall when temperature differences are negligible the flow asymptotically achieves the free stream velocity. The behavior of increase in the values of Le on the velocity field 𝑓^′(𝜂) can be seen from Figure 6. The increase in the values of Le increases the velocity and momentum boundary layer thickness. Physically it validates the fact that increase in the values of Le increases thermal diffusivity in free convection which is due to a higher temperature therefore higher Lewis number results in higher flow field.

Figure 7 depicts the effects of dimensionless temperature and concentration profiles for buoyancy assisted flow 𝑇_𝑤>𝑇_∞ for various values of Le. At the wall the temperature is maximum; as it moves away it is noticed the temperature reduces and asymptotically it approaches to the free stream temperature. It is observed with the increase in the values of Le the temperature and concentration profiles decreases throughout in the boundary layer. This implies that thermal boundary layer thickness decreases with the increase in the values of Le but the effect of Le is more prominent for concentration profiles. It is due to the fact the increase in the Lewis number reduces the mass diffusivity therefore heat transfer and mass transfer slow down through the fluid because the fluid conducts heat slowly relative to its volumetric heat capacity.

Analysis of Figures 8 and 9 reveals the strong effects of buoyancy coefficient 𝑁_𝑟 on velocity, temperature, and concentration profiles. In the absence of buoyancy parameter the density difference is zero therefore the flow is maximum; as we increase the buoyancy parameter the density difference increases which results in the decrease in the fluid flow. The increase in the buoyancy parameter increases the density difference and mass concentration difference hence it enhances the temperature and concentration profiles which shows that the effects of natural convection are significant. In free convection along with other parameters, thermophoresis parameter plays a significant role in the flow, heat transfer, and mass transfer properties.

Figure 10 illustrates the fluid flow is reduced with the increase in the values of 𝑁_𝑡. It is clear the increase in the values of 𝑁_𝑡 increases the temperature difference which results in the decrease in the fluid flow hence nanoparticles in the fluid opposes the transport phenomena through the parameter 𝑁_𝑡. Further, it can be seen the highest value of the velocity occur in the ambient fluid close to the surface but not at the surface.

Figure 11 illustrates the effects of 𝑁_𝑡 on the fluid temperature and concentration profiles. The increase in the values of 𝑁_𝑡 increases the temperature and concentration profiles hence the nanoparticles enhance heat and mass transfer. It is observed that the effects of 𝑁_𝑡 are more pronounced for concentration profiles.

Figure 12 depicts that the increase in the values of 𝑁_𝑏 increases the flow because the increase in Brownian motion parameter will increase the mass diffusivity which enhances the momentum transfer. Figure 13 depicts that the increase in Brownian motion parameter increases the mass diffusivity which enhances the heat transfer process therefore temperature profiles increases with the increase in the values of 𝑁_𝑏. The increase in the values of 𝑁_𝑏 reduces thermal diffusivity therefore the mass transfer will decrease hence the increase in the values of 𝑁_𝑏 decreases the concentration profiles.

5. Conclusions

The problem of steady flow and heat transfer of free convection boundary layer flow in the stagnation-point region in the presence of nanoparticles is studied analytically. The governing equations are nondimensionalized using proper nondimensional quantities. The resulting boundary value problem is solved by optimal homotopy analysis method. A new kind of total averaged squared residual error is used to obtain the optimal values of the convergence-control parameters. The procedure to obtain global convergence control parameters is defined in detail. The effects of the various values of the pertinent parameters 𝑐, 𝑁_𝑏, 𝑁_𝑡, 𝑁_𝑟, and Le on the fluid velocity, temperature profile, and concentration profile are illustrated through graphs. We may extract some important findings from our results.(i)The increase in the values of 𝑐 reduces the skin friction.(ii)The increase in the values of 𝑐 increases the Local Nusselt number and local Sherwood number therefore it enhances heat and mass transfer.(iii)The increase in the values of Le and 𝑁_𝑏 will enhance flow in both 𝑥- and 𝑦-directions.(iv)The effect of the increase in the values of 𝑁_𝑟 and 𝑁_𝑡 is to decelerate the flow.(v)The increase in the values of 𝑁_𝑟, 𝑁_𝑏, and 𝑁_𝑡 increases the temperature while quite the opposite is true for Le.(vi)The effect of the increase in the values of Le and 𝑁_𝑏 is to decrease the concentration profiles.(vii)The effect of the increase in the values of 𝑁_𝑟 and 𝑁_𝑡 is to increase the concentration profiles.

We may conclude the flow and heat transfer properties of free convection boundary layer flow in the stagnation-point region in the presence of nanoparticles can be controlled by changing the quantity of the physical parameters. Hence OHAM is very effective method to solve strongly nonlinear problems.

Appendices

A. Homotopy-Based Approach

We apply the optimal homotopy analysis method to solve the system of nonlinear ordinary differential equations (7)–(10) with boundary conditions (13). In the framework of the HAM the solutions for 𝑓(𝜂), 𝑠(𝜂), 𝜃(𝜂), and 𝜙(𝜂) can be expressed explicitly by an infinite number of subfunctions in the following form: 𝑓𝜂=^+∞_𝑘=0𝑓_𝑘𝜂,𝑠𝜂=^+∞_𝑘=0𝑠_𝑘𝜂,𝜃𝜂=^+∞_𝑘=0𝜃_𝑘𝜂,𝜙𝜂=^+∞_𝑘=0𝜙_𝑘𝜂. (A.1) From physical points of view, it is clear that all the solutions should be in the form as given below; therefore 𝑓𝜂=𝐴^𝑘_0,0+^+∞_𝑘=0^+∞_𝑗=1^+∞_𝑖=0𝐴^𝑘_𝑖,𝑗𝜂^𝑖,𝜃𝜂=exp−𝑗𝜂_{+∞𝑘=0+∞𝑗=1+∞𝑖=0}𝐶^𝑘_𝑖,𝑗𝜂^𝑖,𝑠𝜂exp−𝑗𝜂=𝐵^𝑘_0,0+^+∞_𝑘=0^+∞_𝑗=1^+∞_𝑖=0𝐵^𝑘_𝑖,𝑗𝜂^𝑖,𝜙𝜂=exp−𝑗𝜂_{+∞𝑘=0+∞𝑗=1+∞𝑖=0}𝐷^𝑘_𝑖,𝑗𝜂^𝑖exp−𝑗𝜂, (A.2) where 𝐴^𝑘_𝑖,𝑗, 𝐵^𝑘_𝑖,𝑗, 𝐶^𝑘_𝑖,𝑗, and 𝐷^𝑘_𝑖,𝑗 are constant coefficients to be determined. In the HAM, one has great freedom on the choice of the auxiliary linear operator. Thus, we can choose auxiliary linear operators as ℒ₁𝑢=d𝜂∶𝑞³𝑢d𝜂³−d𝑢,ℒd𝜂₂𝑣=d𝜂∶𝑞²𝑣d𝜂²−𝑣, (A.3) which satisfy the following properties: ℒ₁𝐶₁exp−𝜂+𝐶₂exp𝜂+𝐶₃ℒ=0,₂𝐶₉exp−𝜂+𝐶₁₀exp𝜂=0. (A.4) The corresponding auxiliary linear operator for 𝑓(𝜂) and 𝑠(𝜂) is ℒ₁ and ℒ₂ corresponds to 𝜃(𝜂) and 𝜙(𝜂). The auxiliary function ℋ_𝑖 (𝑖=1--3) is chosen as ℋ_𝑖=exp−𝜂. (A.5) In HAM, we also have freedom to choose the initial approximations. Note that, all the initial approximations need to satisfy the boundary conditions (14). Without loss of generality, we set 𝑓₀𝜂=12−𝑒^−𝜂+12𝑒^−2𝜂,𝑠₀𝜂=1512−𝑒^−𝜂+12𝑒^−2𝜂,𝜃₀𝜂,𝜙=exp−𝜂₀𝜂=exp−𝜂. (A.6) Thus, in the frame of HAM the 𝑘th order approximation 𝑓_𝑘(𝜂), 𝑠_𝑘(𝜂), 𝜃_𝑘(𝜂), and 𝜙_𝑘(𝜂) can be obtained using the following recursive formulae: 𝑓_𝑘𝜂=𝜒_𝑘𝑓_𝑘−1𝜂+𝑓^∗_𝑘𝜂+𝐶^𝑘₁+𝐶^𝑘₂𝑒^𝜂+𝐶^𝑘₃𝑒^−𝜂,𝑠_𝑘𝜂=𝜒_𝑘𝑠_𝑘−1𝜂+𝑠^∗_𝑘𝜂+𝐶^𝑘₄+𝐶^𝑘₅𝑒^𝜂+𝐶^𝑘₆𝑒^−𝜂,𝜃_𝑘𝜂=𝜒_𝑘𝜃_𝑘−1𝜂+𝜃^∗_𝑘𝜂+𝐶^𝑘₇+𝐶^𝑘₈𝑒^−𝜂,𝜙_𝑘𝜂=𝜒_𝑘𝜙_𝑘−1𝜂+𝜙^∗_𝑘𝜂+𝐶^𝑘₉+𝐶^𝑘₁₀𝑒^−𝜂, (A.7) where 𝜒_𝑚=0,𝑚≤11,𝑚>1. (A.8)𝐶^𝑘₁, 𝐶^𝑘₂, 𝐶^𝑘₃, 𝐶^𝑘₄, 𝐶^𝑘₅, 𝐶^𝑘₆, 𝐶^𝑘₇, 𝐶^𝑘₈, 𝐶^𝑘₉, and 𝐶^𝑘₁₀ are integral constants determined as 𝐶^𝑘₁=−𝑓_𝑘^∗′(0)−𝑓^∗_𝑘(0),𝐶^𝑘₂=0,𝐶^𝑘₃=𝑓_𝑘^∗′(0),𝐶^𝑘₄=−𝑠_𝑘^∗′(0)−𝑠^∗_𝑘(0),𝐶^𝑘₅=0,𝐶^𝑘₆=𝑠_𝑘^∗′,𝐶^𝑘₇=0,𝐶^𝑘₈=1,𝐶^𝑘₉=0,𝐶^𝑘₁₀=1. (A.9) and 𝑓^∗_𝑘(𝜂), 𝑠^∗_𝑘(𝜂), 𝜃^∗_𝑘(𝜂), and 𝜙^∗_𝑘(𝜂) are particular solutions defined by 𝑓^∗_𝑘𝜂=ℒ₁⁻¹𝑐^𝑓₀𝑅^𝑓_𝑘,𝑠^∗_𝑘𝜂=ℒ₁⁻¹𝑐^𝑓₀𝑅^𝑠_𝑘,𝜃^∗_𝑘𝜂=ℒ₂⁻¹𝑐^𝜃₀𝑅^𝜃_𝑘,𝜙^∗_𝑘𝜂=ℒ₂⁻¹𝑐^𝜙₀𝑅^𝜙_𝑘. (A.10)𝑐^𝑓₀, 𝑐^𝑠₀, 𝑐^𝜃₀, and 𝑐^𝜙₀ are convergence-control parameters, 𝑅^𝑓_𝑘, 𝑅^𝑠_𝑘, 𝑅^𝜃_𝑘, and 𝑅^𝜙_𝑘 are defined based on (8)–(13) for details see Liao [30]. Consider 𝑅^𝑓_𝑘=𝑓^′′′_𝑘−1+^𝑘−1_𝑖=0𝑓_𝑖^′′𝑓_{𝑘−1−𝑖}+𝑐^𝑘−1_𝑖=0𝑓_𝑖^′′𝑠_{𝑘−1−𝑖}−^𝑘−1_𝑖=0𝑓^′_𝑖𝑓^′_{𝑘−1−𝑖}+𝜃_𝑘−1−𝑁_𝑟𝜙PrGr^𝑘−1,𝑅^𝑠_𝑘=𝑠^′′′_𝑘−1+^𝑘−1_𝑖=0𝑠_𝑖^′′𝑓_{𝑘−1−𝑖}+𝑐^𝑘−1_𝑖=0𝑠_𝑖^′′𝑠_{𝑘−1−𝑖}−𝑐^𝑘−1_𝑖=0𝑓^′_𝑖𝑓^′_{𝑘−1−𝑖}+𝜃_𝑘−1−𝑁_𝑟𝜙PrGr^𝑘−1,𝑅^𝜃_𝑘=1𝜃Pr^{′′𝑘−1}+^𝑘−1_𝑖=0𝜃^′_𝑖𝑓_{𝑘−1−𝑖}+𝑐^𝑘−1_𝑖=0𝜃^′_𝑖𝑠_{𝑘−1−𝑖}+𝑁_𝑏Pr^𝑘−1_𝑖=0𝜃^′_𝑖𝜙^′_{𝑘−1−𝑖}+𝑁_𝑡Pr^𝑘−1_𝑖=0𝜃^′_𝑖𝜃^′_{𝑘−1−𝑖},𝑅^𝜙_𝑘=𝜙^′′_𝑘−1+𝑁_𝑡𝑁_𝑏+LePr^𝑘−1_𝑖=0𝜙^′_𝑖𝑓_{𝑚−1−𝑖}+𝑐^𝑘−1_𝑖=0𝜙^′_𝑖𝑠_{𝑘−1−𝑖}. (A.11) Using the above recursive formulae, we can get the explicit analytical approximations of 𝑓_𝑘(𝜂), 𝑠_𝑘(𝜂), 𝜃_𝑘(𝜂), and 𝜙_𝑘(𝜂) for 𝑘=1,2,3,….

B. Formulation of Explicit Solutions

The corresponding coefficients in (A.2) are obtained by means of HAM as 𝑓𝜂=𝐴^𝑘_0,0+^+∞_𝑗=1^+∞_𝑖=0𝐴^𝑘_𝑖,𝑗𝜂^𝑖,𝜃𝜂=exp−𝑗𝜂_{+∞𝑗=1+∞𝑖=0}𝐶^𝑘_𝑖,𝑗𝜂^𝑖,𝑠𝜂exp−𝑗𝜂=𝐵^𝑘_0,0+^+∞_𝑗=1^+∞_𝑖=0𝐵^𝑘_𝑖,𝑗𝜂^𝑖,𝜙𝜂=exp−𝑗𝜂_{+∞𝑗=1+∞𝑖=0}𝐷^𝑘_𝑖,𝑗𝜂^𝑖,𝐴exp−𝑗𝜂_𝑘0,0=𝜒_𝑘𝐴^𝑘−1_0,0+^𝑘+1_𝑗=2^𝑘+1_𝑞=1Γ^𝑘_𝑞,𝑗𝑗−1𝜔^𝑞_𝑗,0−1𝜆𝜔^𝑞_𝑗,1−^𝑘_𝑞=01𝜆Δ^𝑘_𝑞,1𝜔^𝑞_1,1+^𝑘+1_𝑗=2𝑗−1Γ^𝑘_0,𝑗𝜔⁰_𝑗,0,𝐴^𝑘_0,1=𝜒_𝑘𝜎^𝑘−1_0,1𝐴^𝑘−1_0,1+^𝑘+1_𝑗=2^𝑘+1_𝑞=11𝜆Γ^𝑘_𝑞,𝑗𝜔^𝑞_𝑗,1+^𝑘_𝑞=01𝜆Δ^𝑘_𝑞,1𝜔^𝑞_1,1−^𝑘+1_𝑗=2^𝑘+1_𝑞=0𝑗Γ^𝑘_𝑞,𝑗𝜔^𝑞_𝑗,0,𝐴^𝑘_0,𝑗,𝐴=0,2≤𝑗≤𝑘+1_𝑘𝑖,1=𝜒_𝑘𝜎^𝑘−1_𝑖,1𝐴^𝑘−1_0,1+^𝑘+1_𝑞=𝑖Δ^𝑘_𝑞,1𝜔^𝑞_1,𝑖,(),𝐴1≤𝑖≤𝑘+1_{𝑘𝑖,𝑗}=𝜒_𝑘𝜎^𝑘−1_𝑖,𝑗𝐴^𝑘−1_0,1+^𝑘+1_𝑞=𝑖Γ^𝑘_𝑞,𝑗𝜔^𝑞_𝑗,1,,𝐵1≤𝑖≤𝑘+1,2≤𝑗≤𝑘+1_𝑘0,0=𝜒_𝑘𝐵^𝑘−1_0,0+^𝑘+1_𝑗=2^𝑘+1_𝑞=1Θ^𝑘_𝑞,𝑗𝑗−1𝜔^𝑞_𝑗,0−1𝜆𝜔^𝑞_𝑗,1−^𝑘_𝑞=01𝜆Λ^𝑘_𝑞,1𝜔^𝑞_1,1+^𝑘+1_𝑗=2𝑗−1Θ^𝑘_0,𝑗𝜔⁰_𝑗,0,𝐵^𝑘_0,1=𝜒_𝑘𝜎^𝑘−1_0,1𝐵^𝑘−1_0,1+^𝑘+1_𝑗=2^𝑘+1_𝑞=11𝜆Θ^𝑘_𝑞,𝑗𝜔^𝑞_𝑗,1+^𝑘_𝑞=01𝜆Λ^𝑘_𝑞,1𝜔^𝑞_1,1−^𝑘+1_𝑗=2^𝑘+1_𝑞=0𝑗Θ^𝑘_𝑞,𝑗𝜔^𝑞_𝑗,0,𝐵^𝑘_0,𝑗,𝐵=0,2≤𝑗≤𝑘+1_𝑘𝑖,1=𝜒_𝑘𝜎^𝑘−1_𝑖,1𝐵^𝑘−1_0,1+^𝑘+1_𝑞=𝑖Θ^𝑘_𝑞,1𝜔^𝑞_1,𝑖,(),𝐶1≤𝑖≤𝑘+1_𝑘0,1=𝜒_𝑘𝜎^𝑘−1_0,1𝐶^𝑘−1_0,0−^𝑘+1_𝑗=2^𝑘+1_𝑞=0Ψ^𝑘_𝑞,𝑗𝛿^𝑞_𝑗,0,𝐶^𝑘_𝑖,1=𝜒_𝑘𝛿^𝑘−1_𝑖,1𝐶^𝑘−1_𝑖,1−^𝑘+1_{[]𝑞=max0,𝑖−1}Ω^𝑘_𝑞,1𝛿^𝑞_1,𝑖(),𝐶1≤𝑖≤𝑘+1_𝑘0,𝑗,𝐶=0,2≤𝑗≤𝑘+1_{𝑘𝑖,𝑗}=𝜒_𝑘𝛿^𝑘−1_𝑖,𝑗𝐶^𝑘−1_𝑖,𝑗+^𝑘+1_𝑞=𝑖Ω^𝑘_𝑞,𝑗𝛿^𝑞_𝑗,𝑖,𝐷1≤𝑖≤𝑘+1,2≤𝑗≤𝑘+1_𝑘0,1=𝜒_𝑘𝜎^𝑘−1_0,1𝐷^𝑘−1_0,0−^𝑘+1_𝑗=2^𝑘+1_𝑞=0Φ^𝑘_𝑞,𝑗𝛿^𝑞_𝑗,0,𝐷^𝑘_𝑖,1=𝜒_𝑘𝛿^𝑘−1_𝑖,1𝐷^𝑘−1_𝑖,1−^𝑘+1_{[]𝑞=max0,𝑖−1}Υ^𝑘_𝑞,1𝛿^𝑞_1,𝑖(),𝐷1≤𝑖≤𝑘+1_𝑘0,𝑗,𝐷=0,2≤𝑗≤𝑘+1_{𝑘𝑖,𝑗}=𝜒_𝑘𝛿^𝑘−1_𝑖,𝑗𝐷^𝑘−1_𝑖,𝑗+^𝑘+1_𝑞=𝑖Φ^𝑘_𝑞,𝑗𝛿^𝑞_𝑗,𝑖1≤𝑖≤𝑘+1,2≤𝑗≤𝑘+1, (B.1) with 𝜎^𝑘_𝑖,𝑗=𝜔1,0≤𝑖≤𝑘+1,1≤𝑗≤𝑘+10,otherwise,_{𝑞𝑛,𝑗}=^{𝑞+1𝑖=𝑗}𝛿^𝑞_𝑛,𝑗𝜆^{𝑖−𝑗+1},0≤𝑗≤𝑞+1,𝑛=1,_{𝑞𝑖=𝑗}𝛿^𝑞_𝑛,𝑗(𝑛𝜆)^{𝑖−𝑗+1}𝛿,0≤𝑗≤𝑞+1,𝑛>1,0,otherwise,_{𝑞𝑛,𝑗}=𝑞!1𝑗!(2𝜆)^𝑞1−𝑗+2,0≤𝑗≤𝑞,𝑛=1,12𝜆𝑞+1,0≤𝑗≤𝑞+1,𝑛>1,2𝜆𝑞!1𝑗![(𝑛+1)𝜆]^𝑞−−𝑗+1𝑞!1𝑗![(𝑛−1)𝜆]^𝑞,−𝑗+10≤𝑗≤𝑞,𝑛>10,𝑗=𝑞+1,𝑛>1. (B.2) The functions in the above recurrence formulae are given by Γ^𝑘_𝑖,𝑗=𝑐^𝑓₀𝐴₃^𝑘−1_𝑖,𝑗+𝐶^𝑘−1_𝑖,𝑗−𝑁_𝑟𝐷PrGr^{𝑘−1𝑖,𝑗}+Δ^𝑘_1𝑖,𝑗+Δ^𝑘_2𝑖,𝑗+Δ^𝑘_3𝑖,𝑗+Γ^𝑘_1𝑖,𝑗+Γ^𝑘_2𝑖,𝑗,Θ^𝑘_𝑖,𝑗=𝑐^𝑠₀𝐵₃^𝑘−1_𝑖,𝑗+𝐶^𝑘−1_𝑖,𝑗−𝑁_𝑟𝐷PrGr^{𝑘−1𝑖,𝑗}+Δ^𝑘_4𝑖,𝑗+Δ^𝑘_5𝑖,𝑗+Δ^𝑘_6𝑖,𝑗+Γ^𝑘_3𝑖,𝑗+Γ^𝑘_4𝑖,𝑗,Ψ^𝑘_𝑖,𝑗=𝑐^𝜃₀1𝐶Pr^{3𝑘−1𝑖,𝑗}+Δ^𝑘_7𝑖,𝑗+Δ^𝑘_8𝑖,𝑗+Δ^𝑘_9𝑖,𝑗+Δ^𝑘_10𝑖,𝑗+Γ^𝑘_5𝑖,𝑗+Γ^𝑘_6𝑖,𝑗,Φ^𝑘_𝑖,𝑗=𝑐^𝜙₀𝐷^𝑘−1_𝑖,𝑗+𝑁_𝑡𝑁_𝑏𝐶^𝑘−1_𝑖,𝑗+Δ^𝑘_11𝑖,𝑗+Δ^𝑘_12𝑖,𝑗+Γ^𝑘_7𝑖,𝑗+Γ^𝑘_8𝑖,𝑗,Δ^𝑘_𝑖,1=𝑐^𝑓₀𝐴₃^𝑘−1_𝑖,1+𝐶^𝑘−1_𝑖,1−𝑁_𝑟𝐷PrGr^{𝑘−1𝑖,1}+Γ^𝑘_1𝑖,1+Γ^𝑘_2𝑖,1,Λ^𝑘_𝑖,1=𝑐^𝑠₀𝐵₃^𝑘−1_𝑖,1+𝐶^𝑘−1_𝑖,1−𝑁_𝑟𝐷PrGr^{𝑘−1𝑖,1}+Γ^𝑘_3𝑖,1+Γ^𝑘_4𝑖,1,Ω^𝑘_𝑖,1=𝑐^𝜃₀𝐶₃^𝑘−1_𝑖,1+Γ^𝑘_5𝑖,1+Γ^𝑘_6𝑖,1,Υ^𝑘_𝑖,1=𝑐^𝜙₀𝐷₃^𝑘−1_𝑖,1+Γ^𝑘_7𝑖,1+Γ^𝑘_8𝑖,1, (B.3) where 𝑐^𝑓_𝑚, 𝑐^𝑠_𝑚, 𝑐^𝜃_𝑚, and 𝑐^𝜙_𝑚 are the HAM auxiliary parameters, and 𝐴^𝑘_1𝑖,𝑗=(𝑖+1)𝜎^𝑘_𝑖+1,𝑗𝐴^𝑘_𝑖+1,𝑗−𝑗𝜆𝜎^𝑘_𝑖+1,𝑗𝐴^𝑘_𝑖,𝑗,𝐵^𝑘_1𝑖,𝑗=(𝑖+1)𝜎^𝑘_𝑖+1,𝑗𝐵^𝑘_𝑖+1,𝑗−𝑗𝜆𝜎^𝑘_𝑖+1,𝑗𝐵^𝑘_𝑖,𝑗,𝐶^𝑘_1𝑖,𝑗=(𝑖+1)𝜎^𝑘_𝑖+1,𝑗𝐶^𝑘_𝑖+1,𝑗−𝑗𝜆𝜎^𝑘_𝑖+1,𝑗𝐶^𝑘_𝑖,𝑗,𝐷^𝑘_1𝑖,𝑗=(𝑖+1)𝜎^𝑘_𝑖+1,𝑗𝐷^𝑘_𝑖+1,𝑗−𝑗𝜆𝜎^𝑘_𝑖+1,𝑗𝐷^𝑘_𝑖,𝑗,𝐴^𝑘_2𝑖,𝑗=(𝑖+2)(𝑖+1)𝜎^𝑘_𝑖+2,𝑗𝐴^𝑘_𝑖+2,𝑗−2𝑗𝜆(𝑖+1)𝜎^𝑘_𝑖+1,𝑗𝐴^𝑘_𝑖+1,𝑗+𝑗𝜆²(𝑖+1)𝜎^𝑘_𝑖,𝑗𝐴^𝑘_𝑖,𝑗,𝐵^𝑘_2𝑖,𝑗=(𝑖+2)(𝑖+1)𝜎^𝑘_𝑖+2,𝑗𝐵^𝑘_𝑖+2,𝑗−2𝑗𝜆(𝑖+1)𝜎^𝑘_𝑖+1,𝑗𝐵^𝑘_𝑖+1,𝑗+𝑗𝜆²(𝑖+1)𝜎^𝑘_𝑖,𝑗𝐵^𝑘_𝑖,𝑗,𝐶^𝑘_2𝑖,𝑗=(𝑖+2)(𝑖+1)𝜎^𝑘_𝑖+2,𝑗𝐶^𝑘_𝑖+2,𝑗−2𝑗𝜆(𝑖+1)𝜎^𝑘_𝑖+1,𝑗𝐶^𝑘_𝑖+1,𝑗+𝑗𝜆²(𝑖+1)𝜎^𝑘_𝑖,𝑗𝐶^𝑘_𝑖,𝑗,𝐷^𝑘_2𝑖,𝑗=(𝑖+2)(𝑖+1)𝜎^𝑘_𝑖+2,𝑗𝐷^𝑘_𝑖+2,𝑗−2𝑗𝜆(𝑖+1)𝜎^𝑘_𝑖+1,𝑗𝐷^𝑘_𝑖+1,𝑗+𝑗𝜆²(𝑖+1)𝜎^𝑘_𝑖,𝑗𝐷^𝑘_𝑖,𝑗,𝐴^𝑘_3𝑖,𝑗=(𝑖+3)(𝑖+2)(𝑖+1)𝜎^𝑘_𝑖+3,𝑗𝐴^𝑘_𝑖+3,𝑗−3𝑗𝜆(𝑖+2)(𝑖+1)𝜎^𝑘_𝑖+2,𝑗𝐴^𝑘_𝑖+2,𝑗+ß(𝑖+1)𝑗𝜆²𝜎^𝑘_𝑖+1,𝑗𝐴^𝑘_𝑖+1,𝑗−𝑗𝜆³𝜎^𝑘_𝑖,𝑗𝐴^𝑘_𝑖,𝑗,𝐵^𝑘_3𝑖,𝑗=(𝑖+3)(𝑖+2)(𝑖+1)𝜎^𝑘_𝑖+3,𝑗𝐵^𝑘_𝑖+3,𝑗−3𝑗𝜆(𝑖+2)(𝑖+1)𝜎^𝑘_𝑖+2,𝑗𝐵^𝑘_𝑖+2,𝑗+ß(𝑖+1)𝑗𝜆²𝜎^𝑘_𝑖+1,𝑗𝐵^𝑘_𝑖+1,𝑗−𝑗𝜆³𝜎^𝑘_𝑖,𝑗𝐵^𝑘_𝑖,𝑗,Δ^𝑘_1𝑖,𝑗=^𝑘−1_𝑛=0^{min{𝑘−𝑛,𝑗−1}}_{𝑠=max{1,𝑗−𝑛−1}}^{min{𝑘−𝑛,𝑖}}_{𝑟=max{0,𝑖−𝑛−1}}𝜎^{𝑘−1−𝑛}_𝑟,𝑠𝐴^{𝑘−1−𝑛}_𝑟,𝑠𝐴^𝑛_{2𝑖−𝑟,𝑗−𝑠},Δ^𝑘_2𝑖,𝑗=𝑐^𝑘−1_𝑛=0^{min{𝑘−𝑛,𝑗−1}}_{𝑠=max{1,𝑗−𝑛−1}}^{min{𝑘−𝑛,𝑖}}_{𝑟=max{0,𝑖−𝑛−1}}𝜎^{𝑘−1−𝑛}_𝑟,𝑠𝐵^{𝑘−1−𝑛}_𝑟,𝑠𝐴^𝑛_{2𝑖−𝑟,𝑗−𝑠},Δ^𝑘_3𝑖,𝑗=−1^𝑘−1_𝑛=0^{min{𝑘−𝑛,𝑗−1}}_{𝑠=max{1,𝑗−𝑛−1}}^{min{𝑘−𝑛,𝑖}}_{𝑟=max{0,𝑖−𝑛−1}}𝜎^{𝑘−1−𝑛}_𝑟,𝑠𝐴₁^{𝑘−1−𝑛}_𝑟,𝑠𝐴^𝑛_{1𝑖−𝑟,𝑗−𝑠},Δ^𝑘_4𝑖,𝑗=^𝑘−1_𝑛=0^{min{𝑘−𝑛,𝑗−1}}_{𝑠=max{1,𝑗−𝑛−1}}^{min{𝑘−𝑛,𝑖}}_{𝑟=max{0,𝑖−𝑛−1}}𝜎^{𝑘−1−𝑛}_𝑟,𝑠𝐴^{𝑘−1−𝑛}_𝑟,𝑠𝐵^𝑛_{2𝑖−𝑟,𝑗−𝑠},Δ^𝑘_5𝑖,𝑗=𝑐^𝑘−1_𝑛=0^{min{𝑘−𝑛,𝑗−1}}_{𝑠=max{1,𝑗−𝑛−1}}^{min{𝑘−𝑛,𝑖}}_{𝑟=max{0,𝑖−𝑛−1}}𝜎^{𝑘−1−𝑛}_𝑟,𝑠𝐵^{𝑘−1−𝑛}_𝑟,𝑠𝐵^𝑛_{2𝑖−𝑟,𝑗−𝑠},Δ^𝑘_6𝑖,𝑗=−𝑐^𝑘−1_𝑛=0^{min{𝑘−𝑛,𝑗−1}}_{𝑠=max1,𝑗−𝑛−1}^{min{𝑘−𝑛,𝑖}}_{{}𝑟=max0,𝑖−𝑛−1}𝜎^{𝑘−1−𝑛}_𝑟,𝑠𝐵₁^{𝑘−1−𝑛}_𝑟,𝑠𝐵^𝑛_{1𝑖−𝑟,𝑗−𝑠},Δ^𝑘_7𝑖,𝑗=^𝑘−1_𝑛=0^{min{𝑘−𝑛,𝑗−1}}_{𝑠=max{1,𝑗−𝑛−1}}^{min{𝑘−𝑛,𝑖}}_{𝑟=max{0,𝑖−𝑛−1}}𝜎^{𝑘−1−𝑛}_𝑟,𝑠𝐴^{𝑘−1−𝑛}_𝑟,𝑠𝐶^𝑛_{1𝑖−𝑟,𝑗−𝑠},Δ^𝑘_8𝑖,𝑗=𝑐^𝑘−1_𝑛=0^{min{𝑘−𝑛,𝑗−1}}_{𝑠=max{1,𝑗−𝑛−1}}^{min{𝑘−𝑛,𝑖}}_{𝑟=max{0,𝑖−𝑛−1}}𝜎^{𝑘−1−𝑛}_𝑟,𝑠𝐵^{𝑘−1−𝑛}_𝑟,𝑠𝐶^𝑛_{1𝑖−𝑟,𝑗−𝑠},Δ^𝑘_9𝑖,𝑗=𝑁_𝑏Pr^𝑘−1_𝑛=0^{min{𝑘−𝑛,𝑗−1}}_{𝑠=max1,𝑗−𝑛−1}^{min{𝑘−𝑛,𝑖}}_{{}𝑟=max0,𝑖−𝑛−1}𝜎^{𝑘−1−𝑛}_𝑟,𝑠𝐶₁^{𝑘−1−𝑛}_𝑟,𝑠𝐷^𝑛_{1𝑖−𝑟,𝑗−𝑠},Δ^𝑘_10𝑖,𝑗=𝑁_𝑡Pr^𝑘−1_𝑛=0^{min{𝑘−𝑛,𝑗−1}}_{𝑠=max1,𝑗−𝑛−1}^{min{𝑘−𝑛,𝑖}}_{{}𝑟=max0,𝑖−𝑛−1}𝜎^{𝑘−1−𝑛}_𝑟,𝑠𝐶₁^{𝑘−1−𝑛}_𝑟,𝑠𝐶^𝑛_{1𝑖−𝑟,𝑗−𝑠},Δ^𝑘_11𝑖,𝑗=Le⋅Pr^𝑘−1_𝑛=0^{min{𝑘−𝑛,𝑗−1}}_{𝑠=max{1,𝑗−𝑛−1}}^{min{𝑘−𝑛,𝑖}}_{𝑟=max{0,𝑖−𝑛−1}}𝜎^{𝑘−1−𝑛}_𝑟,𝑠𝐴^{𝑘−1−𝑛}_𝑟,𝑠𝐷^𝑛_{1𝑖−𝑟,𝑗−𝑠},Δ^𝑘_12𝑖,𝑗×=Le⋅Pr⋅𝑐_{𝑘−1𝑛=0min{𝑘−𝑛,𝑗−1}𝑠=max{1,𝑗−𝑛−1}min{𝑘−𝑛,𝑖}𝑟=max{0,𝑖−𝑛−1}}𝜎^{𝑘−1−𝑛}_𝑟,𝑠𝐵^{𝑘−1−𝑛}_𝑟,𝑠𝐷^𝑛_{1𝑖−𝑟,𝑗−𝑠},Γ^𝑘_1𝑖,𝑗=^𝑘−1_{𝑛=max{𝑗−1,𝑖−1}}𝐴^{𝑘−1−𝑛}_0,0𝐴^𝑛_2𝑖,𝑗,Γ^𝑘_2𝑖,𝑗=^𝑘−1_{𝑛=max{𝑗−1,𝑖−1}}𝐵^{𝑘−1−𝑛}_0,0𝐴^𝑛_2𝑖,𝑗,Γ^𝑘_3𝑖,𝑗=^𝑘−1_{𝑛=max{𝑗−1,𝑖−1}}𝐴^{𝑘−1−𝑛}_0,0𝐵^𝑛_2𝑖,𝑗,Γ^𝑘_4𝑖,𝑗=^𝑘−1_{𝑛=max{𝑗−1,𝑖−1}}𝐵^{𝑘−1−𝑛}_0,0𝐵^𝑛_2𝑖,𝑗,Γ^𝑘_5𝑖,𝑗=^𝑘−1_{𝑛=max{𝑗−1,𝑖−1}}𝐴^{𝑘−1−𝑛}_0,0𝐶^𝑛_1𝑖,𝑗,Γ^𝑘_6𝑖,𝑗=^𝑘−1_{𝑛=max{𝑗−1,𝑖−1}}𝐵^{𝑘−1−𝑛}_0,0𝐶^𝑛_1𝑖,𝑗,Γ^𝑘_7𝑖,𝑗=^𝑘−1_{𝑛=max{𝑗−1,𝑖−1}}𝐴^{𝑘−1−𝑛}_0,0𝐷^𝑛_1𝑖,𝑗,Γ^𝑘_8𝑖,𝑗=^𝑘−1_{𝑛=max{𝑗−1,𝑖−1}}𝐵^{𝑘−1−𝑛}_0,0𝐷^𝑛_2𝑖,𝑗. (B.4) Therefore by setting 𝐴⁰_0,0, 𝐴⁰_0,1, 𝐵⁰_0,0, 𝐵⁰_0,1, and 𝐶⁰_0,1, all coefficients can be determined successively in the order 𝑘=1,2,3,….

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

This work is supported by the Program for New Century Excellent Talents in University (Grant no. NCET-12-0347), to which the authors are most grateful.