# 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

Authors: Sandile S. Motsa; Zodwa G. Makukula
Journal: Journal of Applied Mathematics
 (2017)

---

## Abstract

A bivariate spectral homotopy analysis method (BSHAM) is extended to solutions of systems of nonlinear coupled partial differential equations (PDEs). The method 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.


---

## Body

## 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 Keller-box 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 atz=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 Ω in the z direction. When t<0, the surface rotates at an angular velocity Ω 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∞ and concentration C∞. 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∞ to Tw (Tw>T∞) and C∞ to Cw (Cw>C∞), respectively. Due to the Coriolis force, the fluid motion is three-dimensional and the corresponding governing equations are given as [6, 7].(1)∂u∂x+∂v∂y+∂w∂z=0,(2)∂u∂t+u∂u∂x+v∂u∂y+w∂u∂z-2Ωv=-1ρ∂p∂x+ν∇2u+gβTT-T∞+βCC-C∞,(3)∂v∂t+u∂v∂x+v∂v∂y+w∂v∂z+2Ωu=-1ρ∂p∂y+ν∇2v,(4)∂w∂t+u∂w∂x+v∂w∂y+w∂w∂z=-1ρ∂p∂z+ν∇2w,(5)∂T∂t+u∂T∂x+v∂T∂y+w∂T∂z=α∂2T∂z2,(6)∂C∂t+u∂C∂x+v∂C∂y+w∂C∂z=D∂2C∂z2,where p denotes the pressure, ρ is the density, ν is the kinematic viscosity, ∇2 is the three-dimensional Laplacian, T and C are the fluid temperature and concentration, respectively, g is the acceleration due to gravity, βT is the thermal expansion coefficient, βC is the concentration expansion coefficient, α is the thermal diffusivity, and D is the diffusion coefficient. The initial and boundary conditions are given as follows:t<0:(7)u=v=w=0,T=T∞,C=C∞,for  anyx,y,z.t≥0:(8)u=ax,v=w=0,T=Tw,C=Cw,atz=0,(9)u⟶0,v⟶0,w⟶0,T⟶0,C⟶0,asz⟶∞.Introducing the following similarity transformations [12] in (2)–(6),(10)ξ=1-e-τ,τ=at,η=aν1/2ξ-1/2z,u=axf′ξ,η,v=axgξ,η,w=-aν1/2ξ1/2fξ,η,T=T∞+Tw-T∞θξ,η,C=C∞+Cw-C∞ϕξ,η,gives rise to the coupled system of nonlinear PDEs(11)f′′′+η21-ξf′′+ξff′′-f′2+2λg+γθ+σϕ=ξ1-ξ∂f′∂ξ,(12)g′′+η21-ξg′+ξfg′-f′g-2λf′=ξ1-ξ∂g∂ξ,(13)1Prθ′′+η21-ξθ′+ξfθ′=ξ1-ξ∂θ∂ξ,(14)1Scϕ′′+η21-ξϕ′+ξfϕ′=ξ1-ξ∂ϕ∂ξ,with the corresponding boundary conditions(15)fξ,0=0,f′ξ,0=1,f′ξ,∞=0,gξ,0=0,gξ,∞=0,θξ,0=1,θξ,∞=0,ϕξ,0=1,ϕξ,∞=0.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(16)λ=Ωa,γ=gβTTw-T∞a2x,σ=gβCCw-C∞a2x,Pr⁡=να,Sc=νD.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 thex- and y-directions using (10) is given by(17)Cfx=ξRex-1/2f′′ξ,0,Cfy=ξRex-1/2g′ξ,0,where Rex=ax2/ν 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(18)Nur=NuxξRex-1/2=-θ′ξ,0.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(19)Shr=ShξRex-1/2=-ϕ′ξ,0.


## 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 ξ direction. The resulting equations are then integrated using the spectral collocation method. Following Motsa [11], we approximate the exact solutions of f(η,ξ), g(η,ξ), θ(η,ξ), and ϕ(η,ξ) by Lagrange interpolation polynomials that interpolate the functions at the so-called Gauss-Lobatto collocation points defined as (20)ζi=cos⁡iπNt,i=0,1,2,…,Nt,where ζ∈[-1,1] and ζ=2ξ-1, so that the unknown functions f(η,ξ), g(η,ξ), θ(η,ξ), and ϕ(η,ξ) are approximated as(21)fη,ξ≈∑j=0Ntfη,ξjLjξ=∑j=0NtfjηLjξ,gη,ξ≈∑j=0Ntgη,ξjLjξ=∑j=0NtgjηLjξ,θη,ξ≈∑j=0Ntθη,ξjLjξ=∑j=0NtθjηLjξ,ϕη,ξ≈∑j=0Ntϕη,ξjLjξ=∑j=0NtϕjηLjξ,with yj(η)=y(η,ξj), y=f,g,θ or ϕ and Lj(ξ) are the characteristic Lagrange cardinal polynomials given by(22)Ljξ=∏j=0j≠kMξ-ξkξj-ξk,Ljξk=δjk=0ifj≠k1ifj=k.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 ξi, i=0,1,2,…,Nt. An important step in the substitution process is the evaluation of the derivatives of Lj(ξ) 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(23)∂f′∂ξξ=ξi=2∑j=0Ntdi,jfj′η,∂g∂ξξ=ξi=2∑j=0Ntdi,jgjη,∂θ∂ξξ=ξi=2∑j=0Ntdi,jθjη,∂ϕ∂ξξ=ξi=2∑j=0Ntdi,jϕjη,i=0,1,2,…,Nt,where di,j(i,j=0,1,…,Nt) are entries of the standard Chebyshev differentiation matrix d=[di,j] of size (Nt+1)×(Nt+1) (see, e.g., [15, 16]).Consequently, substituting (23) in (11)–(14) gives(24)fi′′′+121-ξiηfi′′-2ξi1-ξi∑j=0Ntdi,jfj′+ξififi′′-fi′2+2λgi+γθi+σϕi=0,gi′′+121-ξiηgi′-2ξi1-ξi∑j=0Ntdi,jgj+ξifigi′-fi′gi-2λfi′=0,1Prθi′′+121-ξiηθi′-2ξi1-ξi∑j=0Ntdi,jθj+ξifiθi′=0,1Scϕi′′+121-ξiηϕi′-2ξi1-ξi∑j=0Ntdi,jϕj+ξifiϕi′=0,subject to the boundary conditions(25)fi0=0,fi′0=1,fi′∞=0,gi0=0,gi∞=0,θi0=1,θi∞=0,ϕi0=1,ϕi∞=0,i=0,1,…,Nt.With the initial conditions atξ=0 (which corresponds to ζ=ζNt) known, note that the system (24) forms a system of Nt nonlinear coupled ODEs to be solved for fi, gi, θi, and ϕi,i=0,1,…,Nt-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(26)1-qLfiF0η;q,…,FNt-1η;q-fi,0η=qħNfiF0η;q,…,FNt-1η;q-2ξi1-ξidi,MfNt′,1-qLgiG0η;q,…,GNt-1η;q-gi,0η=qħNgiG0η;q,…,GNt-1η;q-2ξi1-ξidi,MgNt,1-qLθiΘ0η;q,…,ΘNt-1η;q-θi,0η=qħNθiΘ0η;q,…,ΘNt-1η;q-2ξi1-ξidi,MθNt,1-qLϕiΦ0η;q,…,ΦNt-1η;q-ϕi,0η=qħNϕiΦ0η;q,…,ΦNt-1η;q-2ξi1-ξidi,MϕNt,with corresponding boundary conditions(27)Fi0;q=0,Fi′0;q=0,Fi∞;q=0,Gi0;q=0,Gi∞;q=0,Θi0;q=0,Θi∞;q=0,Φi0;q=0,Φi∞;q=0.In the above equations,q∈[0,1] is an embedding parameter and ħ represents a nonzero convergence controlling auxiliary parameter. The initial approximate solutions of fi(η), gi(η), θi(η), and ϕi(η) are given as fi,0(η), gi,0(η), θi,0(η), and ϕi,0(η), respectively, for i=0,1,2,…,Nt-1. The components Lp and Np, p=fi,gi,θi or ϕi, are the linear and nonlinear operators of the unknown functions defined by(28)Lfi=fi′′′+121-ξiηfi′′+ξifi′′-fi′-2ξi1-ξi∑j=0Ntdi,jfj′,Nfi=fi′′′+121-ξiηfi′′-2ξi1-ξi∑j=0Ntdi,jfj′+ξififi′′-fi′2+2λgi+γθi+σϕi,Lgi=gi′′+121-ξiηgi′-2ξi1-ξi∑j=0Ntdi,jgj,Ngi=gi′′+121-ξiηgi′-2ξi1-ξi∑j=0Ntdi,jgj+ξifigi′-fi′gi-2λfi′,Lθi=1Prθi′′+121-ξiηθi′-2ξi1-ξi∑j=0Ntdi,jθj,Nθi=1Prθi′′+121-ξiηθi′-2ξi1-ξi∑j=0Ntdi,jθj+ξifiθi′,Lϕi=1Scϕi′′+121-ξiηϕi′-2ξi1-ξi∑j=0Ntdi,jϕj,Nϕi=1Scϕi′′+121-ξiηϕi′-2ξi1-ξi∑j=0Ntdi,jϕj+ξifiϕi′.The initial approximate solutions are chosen to ensure that(29)Lfifi,0=2ξi1-ξidi,MfNt′,fi,00=0,fi,0′=1,fi,0′∞=0,Lgigi,0=2ξi1-ξidi,MgNt,gi,00=0,gi,0∞=0,Lθiθi,0=2ξi1-ξidi,MθNt,θi,00=0,θi,0∞=0,Lϕiϕi,0=2ξi1-ξidi,MϕNt,ϕi,00=0,ϕi,0∞=0.From the zeroth-order equations (26), it can be shown that as q varies from 0 to 1, the unknown variables Fi(η;q), Gi(η;q), Θi(η;q), and Φi(η;q) vary from their initial approximate solutions fi,0(η), gi,0(η), θi,0(η), and ϕi,0(η) to the corresponding solutions fi(η), gi(η), θi(η), and ϕi(η), respectively, of the governing equations (24). Expanding Fi(η;q), Gi(η;q), Θi(η;q), and Φi(η;q) using Taylor series about q gives(30)Fiη;q=Fiη;0+∑n=1+∞fi,nηqn,Giη;q=Giη;0+∑n=1+∞gi,nηqn,Θiη;q=Θiη;0+∑n=1+∞θi,nηqn,Φiη;q=Φiη;0+∑n=1+∞ϕi,nηqn,where(31)fi,nη=1n!∂nFiη;q∂qnq=0,gi,nη=1n!∂nGiη;q∂qnq=0,θi,nη=1n!∂nΘiη;q∂qnq=0,ϕi,nη=1n!∂nΦiη;q∂qnq=0.Therefore, since Fi(η;1)=fi(η), Gi(η;1)=gi(η), Θi(η;1)=θi(η), Φi(η;1)=ϕi(η) and Fi(η;0)=fi,0(η,ξ), Gi(η;0)=gi,0(η,ξ), Θi(η;0)=θi,0(η,ξ), Φi(η;0)=ϕi,0(η,ξ), we obtain(32)fiη=fi,0η+∑m=1+∞fi,mη,giη=gi,0η+∑m=1+∞gi,mη,θiη=θi,0η+∑m=1+∞θi,mη,ϕiη=ϕi,0η+∑m=1+∞ϕi,mη.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, fi,m, gi,m, θi,m, and ϕi,m, m≥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 m!, and finally setting q=0 given as follows:(33)Lfifi,mη-χmfi,m-1η=ħRi,m-1fi,Lgigi,mη-χmgi,m-1η=ħRi,m-1gi,Lθiθi,mη-χmθi,m-1η=ħRi,m-1θi,Lϕiϕi,mη-χmϕi,m-1η=ħRi,m-1ϕi,with boundary conditions(34)fi,m0=0,fi,m′0=0,fi,m′∞=0,gi,m0=0,gi,m∞=0,θi,m0=0,θi,m∞=0,ϕi,m0=0,ϕi,m∞=0.From (33),(35)Ri,m-1fi=Lfifi,m-1-21-χmξi1-ξidi,MfNt′+ξi2λgi,m-1+γθi,m-1+ξiσϕi,m-1+∑n=0m-1fi,m-1-nfi,n′′-fi,m-1-n′fi,n′,Ri,m-1gi=Lgigi,m-1-21-χmξi1-ξidi,MgNt-2ξλfi,m-1′+ξi∑n=0m-1fi,m-1-ngi,n′-fi,m-1-n′gi,n,Ri,m-1θi=Lθiθi,m-1-21-χmξi1-ξidi,MθNt+ξi∑n=0m-1fi,m-1-nθi,n′,Ri,m-1ϕi=Lϕiϕi,m-1-21-χmξi1-ξidi,MϕNt+ξi∑n=0m-1fi,m-1-nϕi,n′.At this stage, the Chebyshev spectral collocation method is applied independently in the η direction to solve for the initial approximate solutions fi,0, gi,0, θi,0, and ϕi,0 and the higher order deformation equations giving rise to fi,m, gi,m, θi,m, and ϕi,m. 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 fi,m will be given as(36)dpfi,mdηpη=ηr=2ηep∑s=0NDr,spfi,mη¯s=DpFi,m,r=0,1,2,…,N,similar to the other variables. In (36), p is the order of the derivative, D=(2/ηe)[Dr,s] (r,s=0,1,2,…,Nx) with [Dr,s] being an (Nx+1)×(Nx+1) Chebyshev derivative matrix, and the vector Fi,m is defined as(37)Fi,m=fi,mη¯0,fi,mη¯1,…,fi,mη¯Nxt.The derivative formula (36) is a result of a linear transformation used to convert the interval η∈[0,ηe] to η¯∈[-1,1] where ηe is a sufficiently large finite value chosen to approximate the conditions at ∞. This gives rise to the matrix equations:(38)Ai,m-1fiFi,m=ħRi,m-1fi,Ai,m-1giGi,m=ħRi,m-1gi,Ai,m-1θiΘi,m=ħRi,m-1θi,Ai,m-1ϕiΦi,m=ħRi,m-1ϕi,where (39)Ai,ifi=D3+a0D2+a1D-2ξi1-ξidi,iD,Ai,jfi=-2ξi1-ξidi,jD,Ri,m-1fi=ħ+χmD3+a0D2+a1D-2ξi1-ξi∑j=0Nt-1di,jDFi,m-1+ξiγΘi,m-1+σΦi,m-1+2λGi,m-1-2ħξi1-ξidi,NtDFi,Nt+ξiħ∑n=0m-1Fi,m-1-nD2Fi,n-DFi,m-1-nDFi,n,Ai,igi=D2+b0D+-2ξi1-ξidi,iI,Ai,jgi=-2ξi1-ξidi,jI,Ri,m-1gi=ħ+χmD2+b0D-2ξi1-ξi∑j=0Nt-1di,jIGi,m-1-2ħλDFi,m-1-2ħξi1-ξidi,NtGi,Nt+ξiħ∑n=0m-1Fi,m-1-nDGi,n-DFi,m-1-nGi,n.Ai,iθi=D2+c0D+-2ξi1-ξidi,iI,Ai,jθi=-2ξi1-ξidi,jI,Ri,m-1θi=ħ+χmD2+c0D-2ξi1-ξi∑j=0Nt-1di,jIΘi,m-1-2ħξi1-ξidi,NtΘi,Nt+ξiħ∑n=0m-1Fi,m-1-nDΘi,n,Ai,iϕi=D2+d0D+-2ξi1-ξidi,iI,Ai,jϕi=-2ξi1-ξidi,jI,Ri,m-1ϕi=ħ+χmD2+d0D-2ξi1-ξi∑j=0Nt-1di,jIΦi,m-1-2ħξi1-ξidi,NtΦi,Nt+ξiħ∑n=0m-1Fi,m-1-nDΦi,n.


## 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 (η) 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(40)Resf=ΔfFi,Gi,Θi,Φi∞,Resg=ΔgFi,Gi∞,Resθ=ΔθFi,Θi∞,Resϕ=ΔϕFi,Φi∞,where Δf, Δg, Δθ, and Δϕ represent the governing nonlinear PDEs (11), (12), (13), and (14), respectively, and Fi, Gi, Θi, and Φi are the BSHAM approximate solutions at the time collocation points ξi.Figures1–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.Figure 1
Residual error inf(η,ξ) against ħ at different iteration levels (λ=1, =γ=0.5, Pr=3, Sc=2, ξ=0.45).Figure 2
Residual error ing(η,ξ) against ħ at different iteration levels (λ=1,σ=γ=0.5,Pr=3,Sc=2,ξ=0.45).Figure 3
Residual error inθ(η,ξ) against ħ at different iteration levels (λ=1,σ=γ=0.5,Pr=3,Sc=2,ξ=0.45).Figure 4
Residual error inϕ(η,ξ) against ħ at different iteration levels (λ=1,σ=γ=0.5,Pr=3,Sc=2,ξ=0.45).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.Figure 5
Residual error inf(η,ξ) against ξ at different iteration levels (λ=1,σ=γ=0.5,Pr=3,Sc=2,ħ=-0.5).Figure 6
Residual error ing(η,ξ) against ξ at different iteration levels (λ=1,σ=γ=0.5,Pr=3,Sc=2,ħ=-0.5).Figure 7
Residual error inθ(η,ξ) against ξ at different iteration levels (λ=1,σ=γ=0.5,Pr=3,Sc=2,ħ=-0.5).Figure 8
Residual error inϕ(η,ξ) against ξ at different iteration levels (λ=1,σ=γ=0.5,Pr=3,Sc=2,ħ=-0.5).The effect of the different flow parameters on the local skin friction along thex-direction, f′′(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 skin-friction coefficient in Figure 10. This is caused by the fact that γ is positive, thus assisting the flow enhancing the fluid velocity.Figure 9
Effect ofλ on f′′(0,ξ) (σ=γ=0.5,Pr=3,Sc=2,ħ=-0.5).Figure 10
Effect ofγ on f′′(0,ξ) (σ=0.5λ=0.5,Pr=3,Sc=2,ħ=-0.5).Figure 11
Effect ofσ on f′′(0,ξ) (λ=1,γ=0.5,Pr=3,Sc=2,ħ=-0.5).Figure 12
Effect of Pr onf′′(0,ξ) (σ=γ=0.5,λ=1,Sc=2,ħ=-0.5).Figure 13
Effect of Sc onf′′(0,ξ) (σ=γ=0.5,Pr=3,λ=1,ħ=-0.5).A similar effect is observed in Figure11 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 f′′(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.Figures14–18 depict the effect of the flow parameters on the local skin-friction coefficient along the y-direction, g′(0,ξ). As expected, the skin friction along y-direction will also decrease with increase in the rotating ratio parameter λ (presented in Figure 14) as observed earlier with f′′(0,ξ) in Figure 9.Figure 14
Effect ofλ on g′(0,ξ). (σ=γ=0.5,Pr=3,Sc=2,ħ=-0.5).Figure 15
Effect ofγ on g′(0,ξ) (σ=0.5,λ=1,Pr=3,Sc=2,ħ=-0.5).Figure 16
Effect ofσ on g′(0,ξ) (λ=γ=0.5,Pr=3,Sc=2,ħ=-0.5).Figure 17
Effect of Pr ong′(0,ξ). (σ=γ=0.5,λ=1,Sc=2,ħ=-0.5).Figure 18
Effect of Sc ong′(0,ξ) (σ=γ=0.5,Pr=3,λ=1,ħ=-0.5).The effects ofγ and σ on g′(0,ξ) are depicted in Figures 15 and 16, respectively. Increasing either γ or σ 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,ξ).Figures17 and 18 display the effect of the Prandtl and Schmidt number on g′(0,ξ), respectively. For both flow parameters, the effect is insignificant for the unsteady-state flow. A slight increase in g′(0,ξ) with increase in both Pr and Sc is reported towards the steady-state flow.Figure19 presents the effect of the rotating parameter λ on the local Nusselt number. The local Nusselt number 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.Figure 19
Effect ofλ on θ′(0,ξ) (σ=γ=0.5,Pr=3,Sc=2,ħ=-0.5).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 σ.Figure 20
Effect ofγ on θ′(0,ξ) (σ=0.5,λ=1,Pr=3,Sc=2,ħ=-0.5).Figure 21
Effect ofσ on θ′(0,ξ) (λ=1,γ=0.5,Pr=3,Sc=2,ħ=-0.5).The effect of the Prandtl number on the local Nusselt number is presented in Figure22. 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 ξ.Figure 22
Effect of Pr onθ′(0,ξ) (σ=γ=0.5,λ=1,Sc=2,ħ=-0.5).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].Figure 23
Effect of Sc onθ′(0,ξ) (σ=γ=0.5,Pr=3,λ=1,ħ=-0.5).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.The effects of the flow parameters on the mass transfer rate are presented in Figures24–28. The Sherwood number is shown to decrease with increase in λ and Pr in Figures 24 and 27, respectively. The effects of γ, σ, 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].Figure 24
Effect ofλ on ϕ′(0,ξ) (σ=γ=0.5,Pr=3,Sc=2,ħ=-0.5).Figure 25
Effect ofγ on ϕ′(0,ξ) (σ=0.5,λ=1,Pr=3,Sc=2,ħ=-0.5).Figure 26
Effect ofσ on ϕ′(0,ξ) (λ=1,γ=0.5,Pr=3,Sc=2,ħ=-0.5).Figure 27
Effect of Pr onϕ′(0,ξ) (σ=γ=0.5,λ=1,Sc=2,ħ=-0.5).Figure 28
Effect of Sc onϕ′(0,ξ) (σ=γ=0.5,Pr=3,λ=1,ħ=-0.5).


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



---

*Source: 5962073-2017-05-09.xml*