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.
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 at in a rotating fluid. Let be the velocity components in the direction of Cartesian axes , respectively, with the axes rotating at an angular velocity in the direction. When , the surface rotates at an angular velocity in the direction so that the fluid is at rest relative to the surface. Initially (for ), both fluid and plate are stationary and have constant temperature and concentration . At time , the surface at is impulsively stretched in the -direction. For this time, the surface temperature and concentration vary from to () and to (), respectively. Due to the Coriolis force, the fluid motion is three-dimensional and the corresponding governing equations are given as [6, 7].where denotes the pressure, is the density, is the kinematic viscosity, is the three-dimensional Laplacian, and are the fluid temperature and concentration, respectively, is the acceleration due to gravity, is the thermal expansion coefficient, is the concentration expansion coefficient, is the thermal diffusivity, and is the diffusion coefficient. The initial and boundary conditions are given as follows: : :
Introducing the following similarity transformations [12] in (2)–(6),gives rise to the coupled system of nonlinear PDEswith the corresponding boundary conditions
In the above, is an important dimensionless parameter signifying the relative importance of rotation rate to stretching rate, is the buoyancy parameter, is the concentration buoyancy parameter, Pr is the Prandtl number, and Sc the Schmidt number given by
The important physical parameters for the type of boundary layer flow under consideration are the skin-friction coefficient, Nusselt number, and Sherwood number. From the velocity field, the shearing stress at the surface can be obtained, which in the nondimensional form (skin-friction coefficient) in the - and -directions using (10) is given bywhere is the local Reynolds number.
Using the temperature field, the heat transfer coefficient at the surface can be obtained, which, in the nondimensional form, in terms of the Nusselt number, is given by
From the concentration field, we are able to get the mass transfer coefficient at the plate, which, in the nondimensional form, in terms of the Sherwood number, is given by
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 , , , and by Lagrange interpolation polynomials that interpolate the functions at the so-called Gauss-Lobatto collocation points defined as where and , so that the unknown functions , , , and are approximated aswith , or and are the characteristic Lagrange cardinal polynomials given by
Equations (21) are then substituted in the governing nonlinear system of PDEs (11)–(14) and collocation is applied; that is, the equations are required to be satisfied exactly at the collocation points , . An important step in the substitution process is the evaluation of the derivatives of with respect to . Several formulas exist for computing the derivatives if the collocation points are chosen to be Chebyshev-Gauss-Lobatto points. The values of the derivatives at the Chebyshev-Gauss-Lobatto points are computed aswhere are entries of the standard Chebyshev differentiation matrix of size (see, e.g., [15, 16]).
Consequently, substituting (23) in (11)–(14) gives
subject to the boundary conditions
With the initial conditions at (which corresponds to ) known, note that the system (24) forms a system of nonlinear coupled ODEs to be solved for , , , and 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 aswith corresponding boundary conditions
In the above equations, is an embedding parameter and represents a nonzero convergence controlling auxiliary parameter. The initial approximate solutions of , , , and are given as , , , and , respectively, for . The components and , or , are the linear and nonlinear operators of the unknown functions defined byThe initial approximate solutions are chosen to ensure that
From the zeroth-order equations (26), it can be shown that as varies from to , the unknown variables , , , and vary from their initial approximate solutions , , , and to the corresponding solutions , , , and , respectively, of the governing equations (24). Expanding , , , and using Taylor series about giveswhereTherefore, since , , , and , , , , we obtainConvergence of the series (32) depends on a proper choice of the auxiliary parameter . Next, higher order deformation equations are formed, where we are able to get the terms of the series from their solution, , , , and , . The higher order deformation equations are constructed by differentiating the zeroth-order deformation equations (26) times with respect to , then dividing by , and finally setting given as follows:with boundary conditions
From (33),At this stage, the Chebyshev spectral collocation method is applied independently in the direction to solve for the initial approximate solutions , , , and and the higher order deformation equations giving rise to , , , and . In the spectral domain, the derivatives of the variables with respect to are defined in terms of the Chebyshev differentiation matrix; for example, the derivative in will be given assimilar to the other variables. In (36), is the order of the derivative, () with being an Chebyshev derivative matrix, and the vector is defined as
The derivative formula (36) is a result of a linear transformation used to convert the interval to where is a sufficiently large finite value chosen to approximate the conditions at . This gives rise to the matrix equations:where
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 and 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 bywhere , , , and represent the governing nonlinear PDEs (11), (12), (13), and (14), respectively, and , , , and are the BSHAM approximate solutions at the time collocation points .
Figures 1–4 present the residual error curves against the convergence controlling parameter , at various iteration levels. The residual errors are seen to decrease with increase in iterations which shows that the numerical scheme converges. From Figures 1–4, we are also able to pick the value of that will give optimal results.
Variation of the residual error with time is presented in Figures 5–8. Again the accuracy of the results is seen to improve with improvement in iterations proving convergence of the scheme. The error is seen to be uniform for values of time.
The effect of the different flow parameters on the local skin friction along the -direction, , 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.
A similar effect is observed in Figure 11 with because positive values of also enhance the fluid velocity, in turn increasing the local skin friction. The effects of the Prandtl number and the Schmidt number are to decrease the local skin-friction coefficient 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 ; that is, increasing implies an increase in the kinematic viscosity.
Figures 14–18 depict the effect of the flow parameters on the local skin-friction coefficient along the -direction, . As expected, the skin friction along -direction will also decrease with increase in the rotating ratio parameter (presented in Figure 14) as observed earlier with in Figure 9.
The effects of and on are depicted in Figures 15 and 16, respectively. Increasing either or enhances the flow along the -direction, thus decreasing that along the -direction. The momentum boundary layer along -direction decreases, in turn decreasing the local skin-friction coefficient along -direction, .
Figures 17 and 18 display the effect of the Prandtl and Schmidt number on , respectively. For both flow parameters, the effect is insignificant for the unsteady-state flow. A slight increase in with increase in both Pr and Sc is reported towards the steady-state flow.
Figure 19 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.
The effects of and on the local Nusselt number are presented in Figures 20 and 21, respectively. Both buoyancy forces are seen to enhance the temperature field and the thermal boundary layer. This explains the increase in the heat transfer rate with increase in both and .
The effect of the Prandtl number on the local Nusselt number is presented in Figure 22. The effect of Pr is to increase the heat transfer rate. Physically, this is caused by the fact that increasing Pr reduces the temperature field and the thermal boundary layer thickness, in turn increasing the heat transfer rate. The increase with is significantly pronounced for all values of .
Variation of 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 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 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].
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 - and -directions in the process decreasing the local skin-friction coefficients in both directions. The buoyancy forces were also reported to increase both heat and mass transfer rates of the flow. Analysis of the other flow parameters proved to be consistent with similar findings reported in the literature. The study proved to be a success, and the bivariate spectral homotopy analysis method can be used to find solutions of coupled nonlinear systems of partial differential equations that are valid for the whole space and time domain.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.