Research Article | Open Access

# On the Multidomain Bivariate Spectral Local Linearisation Method for Solving Systems of Nonsimilar Boundary Layer Partial Differential Equations

**Academic Editor:**Irena Lasiecka

#### Abstract

In this work, a novel approach for solving systems of nonsimilar boundary layer equations over a large time domain is presented. The method is a multidomain bivariate spectral local linearisation method (MD-BSLLM), Legendre-Gauss-Lobatto grid points, a local linearisation technique, and the spectral collocation method to approximate functions defined as bivariate Lagrange interpolation. The method is developed for a general system of nonlinear partial differential equations. The use of the MD-BSLLM is demonstrated by solving a system of nonlinear partial differential equations that describe a class of nonsimilar boundary layer equations. Numerical experiments are conducted to show applicability and accuracy of the method. Grid independence tests establish the accuracy, convergence, and validity of the method. The solution for the limiting case is used to validate the accuracy of the MD-BSLLM. The proposed numerical method performs better than some existing numerical methods for solving a class of nonsimilar boundary layer equations over large time domains since it converges faster and uses few grid points to achieve accurate results. The proposed method uses minimal computation time and its accuracy does not deteriorate with an increase in time.

#### 1. Introduction

Solutions of nonlinear nonsimilar boundary layer partial differential equations have received considerable attention of research scientist in the past few years. Series solution methods have been used by several researchers. These include asymptotic methods [1, 2], extended series solutions [2], and homotopy analysis method (HAM) [3â€“8]. To achieve accurate results when using series based methods, a large number of terms are required. A complex nonsimilar boundary layer differential equation requires many terms and thus solving the problem becomes tedious. Simplification of a complex model may introduce errors and may lead to inaccurate solutions.

Numerical methods have been used to solve nonsimilar boundary layer equations. Numerical methods like implicit finite difference methods have been used to solve nonsimilar boundary layer equations by Hossain [1, 9â€“11]. Hussain [2] and Yih [12] used the Keller-Box method. Finite difference methods require many grid points to achieve accurate results. The accuracy of the results deteriorates fast as the time domain increases for parabolic nonlinear partial differential equations. Recently, Magagula et al. [13] showed that pseudospectral methods can achieve accurate results with few grid points compared to Keller-Box method. Results were obtained in a fraction of a second.

Pseudospectral methods have been used to solve boundary layer ordinary differential equations [14â€“19], and boundary layer partial differential equations [20â€“22] have been solved using spectral local linearisation method (SLLM) and the bivariate spectral local linearisation method (BSLLM), respectively. Other pseudospectral methods have been used to solve differential equations [23â€“25]. The bivariate spectral local linearisation method was developed and generalized by Motsa [16] for solving systems of boundary layer ordinary differential equations. This method uses quasilinearisation to linearise nonlinear coupled system of differential equations. The result is a system of linear decoupled differential equations which can be solved independently. The SLLM was extended by Motsa [21] to a system of nonlinear coupled boundary layer partial differential equations, where spectral collocation was used in both space and time. The BSLLM method is suitable for nonlinear partial differential equations with small time domain. Increasing the time domain reduces the accuracy of the method. In some cases, the computed numerical solutions fail to converge even when the time step is decreased. Thus, accurate solutions of systems of differential equations are usually determined over short time intervals. The realization that accurate solutions can be obtained in small time intervals have led to a number of numerical methods that uses the domain-decomposition technique. The domain-decomposition approach together with the spectral collocation technique has been used to solve equations arising from various scientific fields. Funaro [26] and Harald [27] solved elliptic differential equations, Kopriva [28] solved hyperbolic systems, Peyret [29] solved stiff problems in fluid mechanics, and recently, it has been used to solve chaotic system of ordinary differential equations by [30â€“33] and Lane-Emden type of equations [34]. These methods split large intervals into smaller subintervals so that the equations are solved over a sequence of nonoverlapping subintervals of the domain. Continuity conditions are used to advance the solution across the nonoverlapping subintervals.

In this paper, we present a multidomain numerical approach for finding solutions of systems of nonlinear coupled nonsimilar boundary layer partial differential equations over a large time interval. The multidomain bivariate spectral local linearisation method (MD-BSLLM) is an extension of the BSLLM. The proposed approach is developed for solving system of coupled nonsimilar boundary layer partial differential equations. The system of equations is solved iteratively using a Jacobi like relaxation approach. The domain is divided into smaller nonoverlapping subintervals on which the Chebyshev spectral collocation method is used to solve the equations. A continuity condition is used to advance the solution across the subintervals. The algorithm is easy to develop and it yields very accurate results using only a few discretization nodes. The accuracy of the method is validated against the series solution for limiting cases. The main aim of the study is to explore the applicability of the multidomain bivariate spectral local linearisation collocation method to systems of nonlinear coupled nonsimilar boundary layer partial differential equations over large time intervals. The results confirm that the method is suitable for solving all types of systems of nonlinear coupled nonsimilar boundary layer partial differential equations over a large time interval.

The article is organized as follows: in Section 2, the multidomain bivariate spectral local linearisation method is described for a general system of partial differential equations. In Section 3, numerical experiments and the series solutions are presented in Section 4. In Section 5, the results are discussed and we conclude in Section 6.

#### 2. Multidomain Bivariate Spectral Local Linearisation Method

In this section, the multidomain bivariate spectral local linearisation method (MD-BSLLM) for approximating solutions of system of nonlinear partial differential equations that model nonsimilar boundary layer equations is introduced and explained in detail. The multidomain bivariate spectral local linearisation method is developed for a general system of nonlinear partial differential equations. Consider a system of nonlinear partial differential equations of the form,whereThe order of differentiation is denoted by ; the required solutions by for and for are the nonlinear operators containing all the spatial derivatives and time derivatives of . The Legendre-Gauss-Lobatto grid points and the corresponding differentiation matrices are defined in the interval . Therefore, the time interval and space region are transformed to using linear transformations and , respectively. Thus . In order to apply the multidomain Legendre-Gauss-Lobatto bivariate spectral local linearisation method on the system of equations under consideration, the interval is first divided into nonoverlapping subintervals for . Each interval is further divided into not necessarily equal divisions. The system of nonlinear partial differential equations under consideration is solved in each subinterval . The solutions in each subinterval are denoted by , for . In the first interval , the solutions , for are obtained subject to the â€œinitial" conditions , for , respectively. For each , at each subinterval , the continuity conditionsare used to implement the MD-BSQLM over the interval . This procedure is repeated to generate a sequence of solutions , for and . It is assumed that, at each subinterval , the solution can be approximated by a bivariate Lagrange interpolation polynomial of the formfor and . The bivariate Lagrange interpolation polynomial interpolates at selected points in both the and directions, for and . The grid points in time are given by the zeros of the polynomialwhere is the first derivative of the Legendre polynomial given byThe zeros of can be computed using any root finding method. In this work, we use Newtonâ€™s method. Hence given and for , , we haveWe note thatand thusSubstituting (15) into (13), we getSimilarly, the grid points in space are given by the zeros of the polynomialand hence the zeros of can be obtained by using the modified Newton methodThe zeros of (11) and (17) are called Legendre-Gauss-Lobatto points [35, 36]. The function is the characteristic Lagrange cardinal polynomial with the Legendre-Gauss-Lobatto points [35, 36]whereAn alternative convenient expression for (19) is given byThe nonlinear operators , for , are first linearised using the quasilinearisation technique as defined by Bellman and Kalaba [37]. The quasilinearisation method is a technique based on the Taylor series expansion of about some previous iteration. It is assumed that the differences between its previous and current solution and all its derivatives are small. Applying the quasilinearisation method yields the following: where and denote previous and current iterations, respectively, and is a vector of the partial derivatives which is defined as We definewhere the prime denotes differentiation with respect to . The linearised equation (22) can be expressed in a compact form as Equations (28)-(31) form a system of decoupled linear partial differential equations. They are solved iteratively for Equations (28)-(31) can further be expressed as follows:where , , and are variable coefficients of , , and , respectively, in the subinterval of the multidomain approach, for and . These variable coefficients correspond to the equation, for . The constant denotes the order of differentiation. Thus, we haveIn general, the th right hand side is given byEquations (32), (33), and (35) are evaluated at the Legendre-Gauss-Lobatto grid points and . The values of the time derivatives are computed at the Legendre-Gauss-Lobatto points , as (for )where is the standard first derivative Legendre-Gauss-Lobatto based differentiation matrix of size , given by [36]The values of the space derivatives at the Legendre-Gauss-Lobatto points (for ) are similarly computed aswhere is the standard first derivative Legendre-Gauss-Lobatto based differentiation matrix of size as defined in [36]. The entries of the space first derivative matrix with respect to the Legendre-Gauss-Lobatto points are exactly the same as those given by (42). Higher th order derivatives are defined aswhere the vector is defined asand the superscript denotes matrix transpose. Substituting (44) and (45) into (32), (33), and (35) yieldswhereThe diagonal matrices of the corresponding variable coefficients are given byImposing boundary conditions for , (47), (48), and (50) can be expressed as the following matrix systemwhereThe vector is defined as for and . The vector corresponds to the initial boundary condition which is always prescribed.

#### 3. Numerical Experiments

In this section, we present systems of equation for various fluid flow models considered in this study. The MD-BSLLM algorithm is applied to each system of equations as a numerical experiment.

##### 3.1. Steady Free Convection Flow Past a Nonisothermal Vertical Porous Cone

In this subsection, the problem of steady two-dimensional laminar free convection flow past a nonisothermal vertical porous cone with variable surface temperature is also considered for numerical experimentation. The governing nonsimilarity system of partial differential equations is expressed in dimensionless form as (see Hossain et. al. [1] for derivation)where is the Prandtl number, is the dimensionless suction parameter, and and are the dimensionless temperature and stream functions, respectively. The appropriate corresponding boundary conditions areThe skin friction coefficient and the Nusselt number describe the shear-stress and heat flux rate at the surface, respectively, and are defined by [1] asIn this case, the order of differentiation and the number of system of equations . In each subinterval , we have the equationwhereWe note that the equations are in terms of the functions and . In order to apply our method, we let and . Therefore, we have The variable coefficients in (64) are given byand the variable coefficients for (65) are given byEquations (62) and (63) can be expressed respectively as the following matrix system:whereThe vectors and are defined as for .

##### 3.2. Steady Free Convective Flow of a Viscous Incompressible Fluid along a Permeable Vertical Flat Plate

A two-dimensional steady free convective flow of a viscous incompressible fluid along a permeable vertical flat plate in the presence of soluble species is being considered in this subsection. The problem considered here is that of a natural convection boundary layer flow, influenced by the combined buoyancy forces from mass and thermal diffusion from a permeable vertical flat surface with nonuniform surface temperature and nonuniform surface species concentration, but with a uniform rate of suction of fluid through the permeable surface. The transformed boundary layer equations are solved numerically near to and far from the leading edge, using the multidomain bivariate spectral collocation method. The equations were solved by Hussain [2] and areThe appropriate boundary conditions are given byIn this case, the order of differentiation and the number of system of equations . Similarly, in each subinterval , we havewhereThe system of equations is in terms of the functions, , , and and in order to apply our method, we set , , and . This leads to The variable coefficients for (82) are given by