Abstract
The spectral homotopy analysis method is extended to solutions of systems of nonlinear partial differential equations. The SHAM has previously been successfully used to find solutions of nonlinear ordinary differential equations. We solve the nonlinear system of partial differential equations that model the unsteady nonlinear convective flow caused by an impulsively stretching sheet. The numerical results generated using the spectral homotopy analysis method were compared with those found using the spectral quasilinearisation method (SQLM) and the two results were in good agreement.
1. Introduction
The current study serves to demonstrate the extension of the spectral homotopy analysis method (SHAM) to systems of nonlinear partial differential equations. The method, originally introduced by Motsa et al. [1, 2], was used to solve nonlinear ordinary differential equations. A recent modification (see Motsa [3]) allows the method to be used to find solutions of nonlinear PDEs. A nonlinear system of partial differential equations that describe unsteady nonlinear convective flow caused by an impulsively stretched plate is used as the test problem in the study. The spectral quasilinearisation method (SQLM), an adaptation of the quasilinearisation method (Bellman and Kalaba [4]), is used as a benchmark to prove the accuracy of the spectral homotopy analysis method.
The test equations for this study are obtained by considering unsteady nonlinear convection flow over an impulsively stretched flat surface. Fluid flow over stretching surfaces is important in many practical applications such as extrusion of plastic sheets, paper production, glass blowing, metal spinning, and drawing plastic films. The quality of the final product depends on the rate of heat transfer at the stretching surface. In a previous study in this field, Kumari et al. [5] used the Keller box method and the Nakamura method to investigate the problem of heat transfer in the unsteady free convection flow over a continuous moving vertical sheet in an ambient fluid. Ishak et al. [6] investigated theoretically the unsteady mixed convection boundary layer flow and heat transfer due to a stretching vertical surface in a quiescent viscous and incompressible fluid. Further, Pop and Na [7] and Wang et al. [8] dealt with the unsteady boundary layer flow due to impulsive flow starting from rest of a stretching sheet in a viscous fluid. A numerical solution for unsteady mixed convection boundary layer nanofluid flow and heat transfer due to a stretching vertical sheet was presented by Mahdy [9]. The equations were solved using the implicit finite difference method. Bachok et al. [10] studied the flow and heat transfer problem due to the unsteady, two-dimensional laminar flow of a viscous nanofluid caused by a permeable stretching/shrinking sheet in a quiescent fluid. Numerical solutions of the transformed governing equations were obtained using a shooting method. Sharma et al. [11] used the finite element method to find numerical solutions of the flow and heat transfer problem over a stretching sheet immersed in a nanofluid, with velocity slip at the boundary.
In fluid flow problems involving heat transfer, it may be essential to consider a nonlinear relationship between density and temperature. Thermal stratification and heat released by the viscous dissipation triggers some changes in density gradients. To the best of our knowledge there is no study that has been conducted to discuss the effects of nonlinear density on unsteady nonlinear convection flow over an impulsively stretched flat surface. Accordingly, the main aim of the study is to investigate the nonlinear convective flow over an impulsively stretched flat surface under the nonlinear density-temperature relationship. Using the boundary layer approximation, the unsteady momentum, heat, and mass transfer equations are transferred to nonlinear partial differential equations form and solved using the spectral homotopy analysis method (SHAM). This work is the first attempt at applying the SHAM to system of coupled nonlinear PDEs. The SHAM results presented in this study are validated against results from a quasilinearisation based spectral collocation method.
2. Governing Equations
Consider the problem of unsteady nonlinear convection of a fluid over a stretching flat plate. Initially , both the fluid and stretching plate are kept at a constant temperature and concentration where is for a heated plate and corresponding to a cooled plate. We assume that at the velocity of the stretching plate is , where is a positive constant. From the Boussinesq approximation, density is related to temperature and the concentration by the equation In the case of thermal stratification and heat released by viscous dissipation, wall jet like profiles induce significant changes in density gradients, and the density depends on the temperature or temperature and concentration in a nonlinear form (see [12–17]): Karcher and Müller [18] used the formulation below to define the nonlinearity of the relationship between the density, the temperature, and the concentration: where is the constant fluid density, and are the fluid temperature and solutal concentration, respectively, and are the coefficients of thermal and solutal expansion, and denotes the nonlinear coefficient of thermal expansion. Partha [16] investigated the natural nonlinear convection in a non-Darcy porous medium using a temperature-concentration-dependent density relation in the form With the usual Boussinesq and the boundary layer approximations, the governing equations are in the form subject to the boundary conditions The initial conditions are where and are the velocity components along and directions, respectively, and are the local fluid temperature and solute concentration across the boundary layer, respectively, is the kinematic viscosity, is the acceleration due to gravity, and are the thermal expansion coefficients, and are solutal expansion coefficients, is the effective thermal diffusivity, and is the mass diffusivity.
We introduce the following nondimensional variables: The governing equations (6)–(8) along with the boundary conditions (9) can be presented in the form subject to the boundary conditions where the Rayleigh number , the nonlinear temperature parameter , the buoyancy parameter , the Prandtl number , the Schmidt number , the Peclet number , and the Lewis number are defined as The skin friction and heat and mass transfer coefficients are described by the local skin friction coefficient, Nusselt number , and the Sherwood number defined as The nondimensional form of the skin friction coefficient, the reduced Nusselt number, and the reduced Sherwood number are
3. Method of Solution
The system of (12) was solved using the spectral homotopy analysis method. The spectral quasilinearisation method (SQLM) was used to validate the results. The spectral homotopy analysis method (SHAM) was first introduced for the solution of nonlinear ordinary differential equation. The method is a hybrid numerical scheme that combines the underlying ideas of the homotopy analysis method (HAM) and the Chebyshev spectral collocation method. The HAM scheme breaks down a nonlinear differential equation into infinitely many linear ordinary differential equations whose solutions are found analytically. The SHAM provides flexibility when solving linear ODEs through the use of the Chebyshev spectral collocation method. The first attempt to extend the application of the SHAM to PDEs was made by Motsa [3] who solved a nonlinear partial differential equation for an unsteady boundary-layer flow caused by an impulsively stretching plate. The blending of homotopy based methods with other methods has also been considered in [19, 20] wherein the homotopy perturbation method was used in conjunction with Laplace transform to solve partial differential equations. In this work, we extend the application to a system of nonlinear PDEs. In the framework of the SHAM, the nonlinear equations are decomposed into their linear and nonlinear parts as follows: with and () representing the linear and nonlinear operators of the equations, respectively. According to the HAM description, the zeroth order deformation equations are formulated as subject to the boundary conditions The functions , , and are unknown, is the embedding parameter, is the nonzero convergence controlling auxiliary parameter, and the functions , , and are the initial solutions. From the zeroth order deformation equations above, it can be noted that, as increases from to , the functions , , and vary from the initial approximations , , and , respectively, to the solutions , , and of the original equations (12). Using the Taylor series to expand , , and about yields with Therefore, since we obtain The series (23) converge provided the auxiliary parameter and the linear operators are properly chosen. The linear operators are often chosen to be the linear parts of the governing equations. The functions , , and are found from solutions of higher order deformation equations which are obtained by differentiating the zeroth order deformation equations (18) times with respect to , dividing by and setting . This gives rise to the following equations: with boundary conditions Here the prime denotes differentiation with respect to and where At this stage, in the HAM framework, the higher order differential equations (24) are solved analytically. However, here, it is not possible since the linear operators contain variable coefficients. The rule of solution expression is modified in the space direction and we assume solutions of the form Initial approximate solutions are chosen to satisfy the rule of solution expression in (28), the corresponding boundary conditions, and the equations Substituting (28) in (29) and balancing terms of equal order in give whose solutions are, respectively, The function erfc is the standard complementary error function defined by Substituting (28) in the higher order deformation equations (24) gives rise to the iteration scheme with boundary conditions The system (33) is a linear sequence of ordinary differential equations that can be solved by any numerical method. In this study the equations are solved using the spectral methods whose underlying ideas will be sketched in Section 3.1.
3.1. The Spectral Quasilinearisation Method (SQLM)
In this section we present a spectral quasilinearisation method (SQLM) for solving systems of partial differential equations such as those in (12). Consider first an initial unsteady state solution corresponding to and obtained by setting in (12). This gives the equations Solving (35) analytically gives The quasilinearisation method (QLM) is essentially a generalized Newton-Raphson method that was originally used by Bellman and Kalaba [4] for solving functional equations. The QLM scheme is derived by linearising the nonlinear component of the governing equations using Taylor series expansion with the assumption that the difference between the value of the unknown functions at the current iteration level (denoted by ) and the value at the previous iteration level (denoted by ) is small. Applying the QLM on (12) leads to the following system of linear equations: where the coefficients , , and () are known from previous calculations and are given by Starting from the initial approximations (36), the iteration schemes (37) can be solved iteratively for , , and when .
In solving (37), we discretize the equation using the Chebyshev spectral collocation method [21] in the -direction and use an implicit finite difference method in the -direction. The spectral method can be implemented on the interval ; therefore it is necessary to transform the domain of the problem to the interval before applying the spectral method. For the convenience of the numerical computations, the semi-infinite domain in the space direction is approximated by the truncated domain , where is a finite number selected to be large enough to represent the behaviour of the flow properties when is very large. The transformation is used to map the interval to . The basic idea behind the spectral collocation method is the introduction of a differentiation matrix which is used to estimate the derivatives of the unknown variables , , and at the collocation points (grid points) as the matrix vector product. For example, the vector product for will be as follows: where is the number of collocation points, , and is the vector function at the collocation points. Higher order derivatives are obtained as powers of ; that is, where is the order of the derivative. We choose the Gauss-Lobatto collocation points to define the nodes in as The matrix is of size . The grid points on are defined as where and are the total numbers of grid points in the - and -directions, respectively, and is the spacing in the -direction. The finite difference scheme is applied with centering about a midpoint halfway between and . This midpoint is defined as . The derivatives with respect to are defined in terms of the Chebyshev differentiation matrices. Applying the centering about to any function, say , and its associated derivative we obtain Thus, applying the spectral method to (37) with the finite difference in gives where and () are matrices and vectors, respectively, defined as
4. Results and Discussion
In this section we present numerical solutions to (12) computed using the spectral homotopy analysis method. Approximate solutions of the skin friction coefficient, the surface heat transfer, and the surface mass transfer rates at different values of flow parameters are presented and compared with the SQLM solutions. Velocity, temperature, and concentration profiles at different flow parameters are also presented. All SHAM results were generated using , , and unless stated otherwise. These values were found to give accurate solutions after a numerical experimentation. The same and were used for the computation of the SQLM solutions. The and in the tables represent the maximum th and th iteration required to produce converging results. Values of the skin friction coefficient and reduced Nusselt and Sherwood numbers at different values of are presented in Table 1. The table also shows a comparison of the SHAM and SQLM results. As can be seen from the table, the results match perfectly well for the set accuracy level. The data from the table also reveals a steady increase of , , and with .
Table 2 shows a variation of , , and with selected flow parameters. It can be seen from the data in Table 2 that all , and contribute positively towards the skin friction coefficient, Nusselt number, and Sherwood number as they are seen to increase with the increase in the flow parameters. All SHAM results in the table are accurate for up to six decimal places.
Figure 1 shows the variation of the skin friction coefficient with and . The skin friction coefficient is shown to be an increasing function of for values of and while showing a decreasing function for . The skin friction coefficient is shown to decrease with the increase in for any value of and the decrease is more prominent for increasing values of .
The effect of and on the skin friction coefficient is presented in Figure 2. The skin friction coefficient is again seen to be a decreasing function of for any value of or . Increasing or further decreases the skin friction coefficient.
Figure 3 shows the effect of the Prandtl number on the temperature profile. The temperature decreases with an increase in Prandtl number. This is due to the fact that smaller values of are corresponding to larger values of thermal conductivities and therefore heat is able to diffuse away from the stretching sheet.
The effect of the Schmidt number on the concentration profile is shown in Figure 4. It can be seen from the figure that, as the Schmidt number increases, the volume fraction boundary layer decreases across the stretching sheet. A higher Schmidt number implies a lower Brownian diffusion coefficient which will give rise to a shorter penetration depth for the concentration boundary layer.
In Figure 5, we present the solution of concentration profile for various values of the Lewis number (). Increasing the Lewis number decreases the concentration profile. Larger values of the Lewis number correspond to weaker molecular diffusivity and thinner boundary layer thickness.
5. Conclusions
The present study was carried out to extend the application of the spectral homotopy analysis method (SHAM) to systems of nonlinear partial differential equations. The method has been used successfully to find solutions of nonlinear ordinary differential equations. We solve the nonlinear system of partial differential equations governing the unsteady nonlinear convection flow over a stretching flat plate. The spectral quasilinearisation method (SQLM) was used to validate the SHAM results. Numerical values of the skin friction coefficient and the reduced Nusselt and Sherwood numbers were generated using both methods. The results matched well for the desired accuracy. Effects of selected flow parameters on the skin friction coefficient, temperature, and concentration profiles were also presented. The study reveals yet another successful application of the SHAM. The method is a promising tool for solutions of nonlinear partial differential equations and their systems.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
The work is based on the research supported in part by the National Research Foundation of South Africa (Grant no. 85596).