Abstract and Applied Analysis

Volume 2014 (2014), Article ID 848170, 7 pages

http://dx.doi.org/10.1155/2014/848170

## On an Interpolation Based Spectral Homotopy Analysis Method for PDE Based Unsteady Boundary Layer Flows

School of Mathematics, Statistics and Computer Science, University of KwaZulu-Natal, Private Bag X01, Scottsville, Pietermaritzburg 3209, South Africa

Received 12 November 2013; Accepted 22 December 2013; Published 2 January 2014

Academic Editor: Robert A. Van Gorder

Copyright © 2014 S. S. Motsa. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

This work presents a new approach to the application of the spectral homotopy analysis method (SHAM) in solving non-linear partial differential equations (PDEs). The proposed approach is based on an innovative idea of seeking solutions that obey a rule of solution expression that is defined in terms of bivariate Lagrange interpolation polynomials. The applicability and effectiveness of the expanded SHAM approach are tested on a non-linear PDE that models the problem of unsteady boundary layer flow caused by an impulsively stretching plate. Numerical simulations are conducted to generate results for the important flow properties such as the local skin friction. The accuracy of the present results is validated against existing results from the literature and against results generated using the Keller-box method. The preliminary results from the proposed study indicate that the present method is more accurate and computationally efficient than more traditional methods used for solving PDEs that describe nonsimilar boundary layer flow.

#### 1. Introduction

In this work, a new approach for applying the spectral homotopy analysis method (SHAM) is used to solve a partial differential equation (PDE) that models the problem of unsteady boundary layer flow caused by an impulsively stretching plate. The SHAM is a discrete numerical version of the homotopy analysis method (HAM) that has been widely applied to solve a wide variety of nonlinear ordinary and partial differential equations with applications in applied mathematics, physics, nonlinear mechanics, finance, and engineering. A detailed systematic description of the HAM and its applications can be found in two books (and a huge list of references cited therein) [1, 2] by Liao who is credited with developing the method.

The SHAM was introduced in [3, 4] and it uses the principles of the traditional HAM and combines them with the Chebyshev spectral collocation method which is used to solve the so-called higher order deformation equations. One of the advantages of the SHAM is that it can accommodate very complex linear operators in its solution algorithm. In a recent application of the SHAM on solving PDE based problems [5], it was observed that the linear operator that gives the best results is one that is selected as the entire linear part of the governing PDEs. This is in sharp contrast to the traditional HAM which can only admit simple ordinary differential equation based linear differential operators with constant coefficient in its solution of PDEs describing the unsteady boundary layer problems of the type discussed in this paper (see, e.g., [6–17]). In the work of [5], the commonly used ODE based linear operator was identified as a limitation if the governing equation was a nonlinear PDE. Partial derivative linear operators were observed to significantly improve the accuracy and convergence speed of the SHAM approach when solving nonlinear partial differential equations. However, it was also observed that the use of PDE based linear operators leads to a complicated linear sequence of higher order deformation equations whose closed form solution was impossible to find. Consequently, the Chebyshev spectral collocation method [18, 19] coupled with a series form of solution constructed from monomial basis functions () was used to solve the higher order deformation equations. However, as it will be pointed out in this work, this approach is only applicable to cases when and becomes less accurate when is very close to 1. Furthermore, the approach of [5] cannot be used directly to solve evolution equations of the form , where is the solution and is a nonlinear operator which contains all the spatial derivatives of . In this investigation, we propose the use of Lagrange interpolation polynomials in place of the monomial basis functions used in [5] for the solution of the PDEs that describe the unsteady boundary layer flow over an impulsively stretching sheet.

The main objective of this work is to introduce an extended application of the SHAM that uses bivariate Lagrange interpolation polynomials as an alternative method for solving nonlinear PDEs. To validate accuracy of the present SHAM results, the governing PDEs are also solved using the Keller-box implicit finite difference approach. The proposed approach is found to be desirable because it is accurate, computationally efficient, and easy to use.

#### 2. Governing Equations

We consider an unsteady boundary layer developed by an impulsively stretching plate in a constant pressure viscous flow. The governing partial differential equations can be obtained by using the standard stream function formulation in conjunction with the transformations suggested by Williams III and Rhyne [20]. The dimensionless governing equation is obtained (see [10, 21, 22] for details) as subject to the boundary conditions where the primes denote differentiation with respect to and is the dimensionless time-scale defined as with a positive constant and is the time variable. In the analysis of boundary layer flow problems, a quantity that is of physical interest is the skin friction which is given in this model [10, 21, 22], in dimensionless form, as where is the local Reynolds number.

The governing equation (1) can be considered to have an initial solution at which is obtained as a solution of the equation Solving (5) gives where is the standard complementary error function defined by

#### 3. Method of Solution

The solution method employed for solving the governing nonlinear PDE (1) is summarised in the following steps.(1)Convert the nonlinear PDE (1) into sequence nonlinear ODEs using Lagrange polynomial interpolation in the variable.(2)Decompose the nonlinear sequence of ODEs into a sequence of linear ODEs using the Homotopy analysis method (HAM) approach.(3)Solve the resulting sequence of linear ODEs using the Chebyshev collocation spectral method.

We approximate the exact solution of by a Lagrange form of polynomial which interpolates at the selected points (called collocation points): The approximation has the form where and are the characteristic Lagrange cardinal polynomials defined as that obey the Kronecker delta equation

The equations for the solution of are obtained by substituting (9) in (1) and requiring that the equation be satisfied exactly at the collocation points , . An important step in the substitution process is the evaluation of the derivatives of with respect to . The derivatives can be evaluated analytically if the interpolating points are chosen as the Chebyshev-Gauss-Lobatto points , such that The values of the derivatives at the Chebyshev-Gauss-Lobatto points are computed as where are entries of the standard Chebyshev differentiation matrix of size (see, e.g., [18, 19]).

Thus, by substituting (9) in (1) and making use of the spectral differentiation matrices, we obtain subject to the boundary conditions Equations (14) form a sequence of nonlinear ordinary differential equations that are to be solved for subject to the boundary conditions (15) for . We remark that, since the solution for is known exactly when (corresponding to ), we need only solve (14) up to . The next step in the solution procedure is the linearisation of nonlinear ODE system (14) and (15). This is achieved by employing the Homotopy analysis method [1, 2].

In the framework of the HAM, we begin by choosing a suitable linear operator which, as suggested in [3, 4], should be taken to be the entire linear component of the given differential equation. Thus, we write the governing equation as: where is a known function that is obtained as the first derivative of (6) and the linear and nonlinear operators are defined as follows:

The next step of the HAM decomposition involves the construction of the so-called zeroth-order deformation equation which, for (14), is defined as subject to the boundary conditions where is an embedding parameter, denotes a nonzero convergence controlling auxiliary parameter, and is the initial approximation of the solution of for . The initial approximation is chosen in such a way that We remark that, as it was proposed in the SHAM application in nonlinear ODEs [3, 4], it is convenient to obtain the initial guess for the SHAM algorithm as the solution of the equation formed by considering only the linear components of the given nonlinear equation. It can be noted from the zeroth-order deformation equation (18) that as increases from 0 to 1, varies from the initial approximation to the solution of the original equation (14). Expanding using Taylor series about gives where Thus, since and we obtain The series (23) converges when the auxiliary parameter is properly chosen.

The functions , are obtained as a solution of higher order deformations equations. These higher order deformation equations are obtained by differentiating the zero-order deformation equations (18) times with respect to , then dividing by , and finally setting . This gives subject to the boundary conditions where

The solution for the initial approximation and higher order deformation equation giving are obtained using the Chebyshev spectral collocation method which is applied independently in the direction using Chebyshev-Gauss-Lobatto points defined as where is a finite value that is chosen to be sufficiently large to approximate the conditions at . The derivatives with respect to are defined in terms of the Chebyshev differentiation matrix as where is the order of the derivative, () with being an Chebyshev derivative matrix, and the vector is defined as Thus substituting (28) in the equation that gives the initial approximation (20) gives the following matrix system: where

Solving equation (30) gives the initial approximation which is computed as To obtain the approximate solution for , (for ) the spectral collocation method, with discretisation in the direction, is applied to (24). This gives the following matrix system: where and are as defined in (31) and (32) and

The above procedure can easily be extended to any nonlinear partial differential equation system. In the context of the HAM approach the steps involved can be succinctly stated in terms of the HAM jargon as simply *looking for a solution that obeys the following rule of solution expression*:
where is a Lagrange form of polynomial which interpolates independently at selected points in both the and directions defined as
We remark that both and are defined between and suitable linear transformations are employed to transform the original domain of and to .

#### 4. Results and Discussion

In this section we present the numerical solutions of the governing partial differential equations (1) that were computed using the proposed bivariate spectral homotopy analysis method. Starting from the initial analytical solutions at (corresponding to ), the SHAM was used to generate results up to solutions near the steady state values at (corresponding to ). The accuracy of the computed SHAM approximate results was confirmed against numerical results obtained by using the popular Keller-box implicit finite difference method as described by [23]. The Keller-box method is known to be accurate, fast, and easier to program for boundary layer flow problems of the type discussed in this work. The Keller-box method is still the preferred solution method for most researchers (see, e.g., [24, 25]). The algorithm of the Keller-box method begins with the reduction of the governing nonlinear PDEs into a system of first-order equations that are discretized using central differences. The nonlinear algebraic difference equations are linearised using Newton’s method and written in matrix-vector form. The linear matrix systems are solved in an efficient manner using a block-tridiagonal-elimination technique. The grid spacing in both the - and -direction is carefully selected to ensure that the Keller-box computations yield consistent results for the governing velocity and temperature distributions to a convergence level of at least . On the other hand, the SHAM was implemented with collocation points in the variable. Furthermore, the finite value used to approximate the boundary conditions at infinity was set to be .

Using finite terms of the SHAM series we define th order approximation at the collocation points and (for and ) as Assuming that is the SHAM approximate solution obtained using (40) at the collocation (grid) points, the residual error is defined as where is defined as

A crucial factor in obtaining accurate and converging SHAM series solutions is the selection of suitable convergence controlling parameter . The infinity norm of the residual error at particular values of was used to identify the optimal value of that gives the best accuracy. Figure 1 is an illustration of a typical residual error curve that can be used to calculate the optimal value of when . The optimal values of are chosen to be the clearly defined minimum of the residual curve. It can be seen from Figure 1 that the optimal value lies in the range . We also note that the residual error decreases with an increase in the number of collocation (grid) points () used in interpolating in the direction.

Figure 2 gives the residual error in the whole range of using a fixed value of and different collocation points in the direction. We remark that similar curves are obtained for any fixed value of , but, for illustration purposes, was used in Figure 2. It can be seen from the figure that the error decreases with an increase in the number of collocation points used. Another interesting observation to emerge from Figure 2 is that the residual error appears to be nearly uniform for most parts of the domain of . This result is what makes the present SHAM approach with interpolation better than the SHAM version proposed in [5] for the solution of PDEs. We remark that convergence of the SHAM version used in [5] was seen to significantly slow down when was near 1. It can be clearly seen from Figure 2 that the present approach does not suffer from the same limitation.

In Table 1, the results for the skin friction are given for different values of time . Table 1 also gives the number of collocation points () and the computational time required to obtain a solution that is consistent with at least six decimal places. The SHAM results match the Keller-box numerical results exactly for all values of . It can be observed from this table that converged solutions are reached using small number of collocation points for the selected values of . In addition, from the comparison, of the computational times, it is clear that the proposed SHAM approach is more computationally efficient in terms of the amount of time and it takes the method to give the desired results. Consequently, it can be inferred from Table 1 that the SHAM is more efficient than the Keller-box in computing an accurate solution for (1). It is worth mentioning that the apparent computational speed of the proposed SHAM can be explained by the fact that, unlike the Keller-box and other numerical methods, very few grid points in both the and directions are required to give very accurate results.

#### 5. Conclusion

In this paper, a new approach based on the spectral homotopy analysis method with bivariate Lagrange interpolating polynomials was introduced for the solution of partial differential equations. The applicability of the proposed method was tested on the problem of unsteady boundary layer flows caused by an impulsively stretching sheet. A residual error analysis was conducted in order to assess the accuracy of the present method. Computational efficiency of the method was demonstrated by comparing with results obtained using the Keller-Box implicit finite difference method. It was found that the proposed SHAM approach was significantly much faster than the Keller-box method. The numerical results presented in this study clearly demonstrate the potential of the proposed bivariate SHAM approach for the simulation of PDEs with high efficiency and accuracy. In future studies it would be interesting to explore the use of this method in other classes of nonlinear PDEs.

#### Conflict of Interests

The author declares that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgment

This work is based on the research supported in part by the National Research Foundation of South Africa (Grant no. 85596).

#### References

- S. Liao,
*Beyond Perturbation: Introduction to the Homotopy Analysis Method*, vol. 2 of*CRC Series: Modern Mechanics and Mathematics*, Chapman & Hall/CRC Press, Boca Raton, Fla, USA, 2004. View at MathSciNet - S. J. Liao,
*Homotopy Analysis Method in Nonlinear Differential Equations*, Springer, Berlin, Germany, 2012. - S. S. Motsa, P. Sibanda, and S. Shateyi, “A new spectral-homotopy analysis method for solving a nonlinear second order BVP,”
*Communications in Nonlinear Science and Numerical Simulation*, vol. 15, no. 9, pp. 2293–2302, 2010. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - S. S. Motsa, P. Sibanda, F. G. Awad, and S. Shateyi, “A new spectral-homotopy analysis method for the MHD Jeffery-Hamel problem,”
*Computers & Fluids*, vol. 39, no. 7, pp. 1219–1225, 2010. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - S. S. Motsa, “On the practical use of the spectral homotopy analysis method and local linearisation method for unsteady boundary-layer flows caused by an impulsively stretching plate,”
*Numerical Algorithms*, 2013. View at Publisher · View at Google Scholar - A. Ali and A. Mehmood, “Homotopy analysis of unsteady boundary layer flow adjacent to permeable stretching surface in a porous medium,”
*Communications in Nonlinear Science and Numerical Simulation*, vol. 13, no. 2, pp. 340–349, 2008. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - I. Ahmad, M. Sajid, T. Hayat, and M. Ayub, “Unsteady axisymmetric flow of a second-grade fluid over a radially stretching sheet,”
*Computers & Mathematics with Applications*, vol. 56, no. 5, pp. 1351–1357, 2008. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - M. Kumari and G. Nath, “Analytical solution of unsteady three-dimensional MHD boundary layer flow and heat transfer due to impulsively stretched plane surface,”
*Communications in Nonlinear Science and Numerical Simulation*, vol. 14, no. 8, pp. 3339–3350, 2009. View at Publisher · View at Google Scholar · View at Scopus - T. Hayat, M. Qasim, and Z. Abbas, “Homotopy solution for the unsteady three-dimensional MHD flow and mass transfer in a porous space,”
*Communications in Nonlinear Science and Numerical Simulation*, vol. 15, no. 9, pp. 2375–2387, 2010. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - S. Liao, “An analytic solution of unsteady boundary-layer flows caused by an impulsively stretching plate,”
*Communications in Nonlinear Science and Numerical Simulation*, vol. 11, no. 3, pp. 326–339, 2006. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - A. Mehmood, A. Ali, and T. Shah, “Heat transfer analysis of unsteady boundary layer flow by homotopy analysis method,”
*Communications in Nonlinear Science and Numerical Simulation*, vol. 13, no. 5, pp. 902–912, 2008. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - A. Mehmood, A. Ali, H. S. Takhar, and T. Shah, “Unsteady three-dimensional MHD boundary-layer flow due to the impulsive motion of a stretching surface,”
*Acta Mechanica*, vol. 199, no. 1–4, pp. 241–249, 2008. View at Publisher · View at Google Scholar · View at Scopus - M. Sajid, I. Ahmad, T. Hayat, and M. Ayub, “Series solution for unsteady axisymmetric flow and heat transfer over a radially stretching sheet,”
*Communications in Nonlinear Science and Numerical Simulation*, vol. 13, no. 10, pp. 2193–2202, 2008. View at Publisher · View at Google Scholar · View at Scopus - M. Sajid, I. Ahmad, T. Hayat, and M. Ayub, “Unsteady flow and heat transfer of a second grade fluid over a stretching sheet,”
*Communications in Nonlinear Science and Numerical Simulation*, vol. 14, no. 1, pp. 96–108, 2009. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - Y. Tan and S.-J. Liao, “Series solution of three-dimensional unsteady laminar viscous flow due to a stretching surface in a rotating fluid,”
*Journal of Applied Mechanics, Transactions ASME*, vol. 74, no. 5, pp. 1011–1018, 2007. View at Publisher · View at Google Scholar · View at Scopus - H. Xu, S. J. Liao, and I. Pop, “Series solutions of unsteady boundary layer flow of a micropolar fluid near the forward stagnation point of a plane surface,”
*Acta Mechanica*, vol. 184, no. 1–4, pp. 87–101, 2006. View at Publisher · View at Google Scholar · View at Scopus - H. Xu and I. Pop, “Homotopy analysis of unsteady boundary-layer flow started impulsively from rest along a symmetric wedge,”
*Zeitschrift für Angewandte Mathematik und Mechanik*, vol. 88, no. 6, pp. 507–514, 2008. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang,
*Spectral Methods in Fluid Dynamics*, Springer Series in Computational Physics, Springer, New York, NY, USA, 1988. View at MathSciNet - L. N. Trefethen,
*Spectral Methods in MATLAB*, vol. 10 of*Software, Environments, and Tools*, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, Pa, USA, 2000. View at Publisher · View at Google Scholar · View at MathSciNet - J. C. Williams III and T. B. Rhyne, “Boundary layer development on a wedge impulsively set into motion,”
*SIAM Journal on Applied Mathematics*, vol. 38, no. 2, pp. 215–224, 1980. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - R. Nazar, N. Amin, and I. Pop, “Unsteady boundary layer flow due to a stretching surface in a rotating fluid,”
*Mechanics Research Communications*, vol. 31, no. 1, pp. 121–128, 2004. View at Publisher · View at Google Scholar · View at Scopus - R. Seshadri, N. Sreeshylan, and G. Nath, “Unsteady mixed convection flow in the stagnation region of a heated vertical plate due to impulsive motion,”
*International Journal of Heat and Mass Transfer*, vol. 45, no. 6, pp. 1345–1352, 2002. View at Publisher · View at Google Scholar · View at Scopus - T. Cebeci and P. Bradshaw,
*Physical and Computational Aspects of Convective Heat Transfer*, Springer, New York, NY, USA, 1984. View at MathSciNet - A. J. Chamkha, S. Abbasbandy, A. M. Rashad, and K. Vajravelu, “Radiation effects on mixed convection over a wedge embedded in a porous medium filled with a nanofluid,”
*Transport in Porous Media*, vol. 91, no. 1, pp. 261–279, 2012. View at Publisher · View at Google Scholar · View at MathSciNet - A. J. Chamkha, S. Abbasbandy, A. M. Rashad, and K. Vajravelu, “Radiation effects on mixed convection about a cone embedded in a porous medium filled with a nanofluid,”
*Meccanica*, vol. 48, no. 2, pp. 275–285, 2013. View at Publisher · View at Google Scholar · View at MathSciNet