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 for (83) are given byand for (84) are given by Equations (79), (80), and (81) can be expressed, respectively, as the following matrix system:whereThe vectors , and are defined as for .
4. Series Solution for the Limiting Case
In this section, the systems of equations presented above are solved using the series solution approach. These series solutions for the limiting case when is very large are used to validate the results obtained using the multidomain bivariate spectral local linearisation method. Equations of the form (58)-(59) and (75)-(87) can be decomposed to a sequence of coupled linear ordinary differential equations when . The solutions of the resulting systems of equations can be determined using elementary methods for solving differential equations.
Consider the series solution of (58)-(59). For large values of , the dominant terms in (58) are and . Balancing the orders of magnitude of these terms gives . Since , we may infer that is of as . Assuming that has the same order as for large , then . Thus for large values of , we define the following transformations: The transformation equation (98) reduces (58)-(59) to where the primes now denote derivatives with respect to . The corresponding boundary conditions are given byUsing the fact that is large, solutions to (99)-(100) are determined in series form asSubstituting (102) into (99)-(101) and then equating the coefficients of like powers of , we obtain the equations for assubject to the following boundary conditions: The analytical solutions of the zeroth order equations (103)-(105) areFor , we sequentially solve the following system of ordinary differential equations:subject to the following boundary conditionsThe series solution of (75)-(87) is now considered. For large values of , the dominant terms in (75) are and . Balancing terms in (75)-(87), for , we getThis leads to the following transformations:Using the transformations (111)-(112) in (75)-(87), we obtainThe corresponding boundary conditions are given byFor large , we seek solutions to (113)-(116) in series form asSubstituting (117) into (113)-(116) and then equating the coefficients of like powers of , we obtain the equations for assubject to the following boundary conditionsThe sequence of equations obtained when is given as subject to the following boundary conditions:
5. Results and Discussion
In this section, we present and analyze the results obtained using both the multidomain bivariate spectral local linearisation method (MD-BSLLM) and the series solution. The MD-BSLLM was implemented using MATLAB 2013b. The series solutions were obtained using Mathematica’s NDSolve. The parameters to generate all the solutions are given in captions for the tables and graphs. In this article, the physical quantities of interest in the field of fluid dynamics, the skin friction, Sherwood number, and Nusselt number are presented. These physical quantities were obtained using both methods. We present the values of these quantities for different large values of and . The values of these quantities are in excellent agreement for both methods. This in turn implies that the new approach (MD-BSLLM) presented here can be used to solve other nonsimilar boundary layer partial differential equations. The effect of number of intervals in the solution is also presented in this section. Graphs showing the velocity profiles, temperature profiles, and concentration profiles are also presented for various large values of and . These profiles are in excellent agreement with the results presented by Hussain [2].
We first present the solutions of (58) and (59). The physical quantities, skin friction, and the Nusselt number results obtained by solving the governing equations (58)-(59) are given in Table 1. These results were obtained directly by solving the sequences of ordinary differential equations, for the large limiting case, using Mathematica’s DSolve command. The series solutions and the MD-BSLLM solutions are accurate up to at least eight decimal digits as indicated in Table 1. Both solutions are in excellent agreement and thus validating the accuracy of the MD-BSLLM method. We observe that an increase in results in a decrease in both the skin friction and the Nusselt number. It is remarkable that this numerical method was able to produce accurate results for large values of with minimum number of grid points. These results were obtained using , , and converged after five iterations.
Figure 1 shows the velocity profiles for different values of . An increase in leads to a decrease in the velocity. The velocity profile is parabolic and hence cannot be obtained by certain numerical methods like the Gauss-Seiedel based pseudospectral relaxation method. This velocity profile is in excellent agreement with that obtained by Hussain [2]; thus the MD-BSLLM method can be used as a numerical method for solving large parameter partial differential equations.
Figure 2 shows the temperature profile for various values of . The temperature profile behaves like an exponential decay curve, that is, an increase in the pseudosimilarity variable results in a decrease in the temperature. Temperature profiles were also obtained for large values. Decreasing leads to an increase in the temperature profile. Figures 1 and 2 were obtained using and .
For the second problem, the skin friction, Sherwood number, and Nusselt number results obtained by solving the governing equations (75)-(87) are given in Table 2. Similarly, these results were obtained directly by solving the sequences of ordinary differential equations, for the large , using Mathematica’s DSolve command. The series solutions and the MD-BSLLM solutions are accurate up to at least eight decimal digits as indicated in Table 2. The MD-BSLLM and series solutions of (75)-(87) are in excellent agreement and thus validating the accuracy of the MD-BSLLM method. Increasing the transpiration parameter results in a decrease of the skin friction, Sherwood number, and Nusselt number. These results were obtained using few grid points, and .
Table 3 gives a comparison of the results obtained using the MD-BSLLM and the bivariate spectral local linearisation method (BSLLM) for large values of . The main purpose of the results is to confirm that the accuracy of the BSLLM method deteriorates as increases. We observe that the computed values agree up to an average of three digits. This implies that decomposing the main interval into subintervals and solving the nonsimilar partial differential equations in each subinterval improves the accuracy of the BSLLM method. Table 3 in addition displays the computational time of both methods. For , , and , the BSLLM method takes about a fraction of a second to achieve inaccurate results meanwhile the MD-BSLLM method takes about 18 seconds to achieve accurate results.
Table 4 shows a comparison of the BSLLM method and the MD-BSLLM method with different grid points. For the MD-BSLLM method, we use , , and , while for the BSLLM method, we use and to achieve comparable accurate results. Since the BSLLM method requires more grid points to achieve accurate results over a large domain, then it takes more computational time compared to the MD-BSLLM method over the same interval. In Table 4, the MD-BSLLM method takes about seconds to achieve accurate results; meanwhile the BSLLM method takes about seconds to achieve accurate results. This implies that the MD-BSLLM method takes less computational time compared to the BSLLM method to achieve accurate results.
In Table 5, we observe that, for small values of , the results from both methods (BSLLM and MD-BSLLM) are in good agreement. The BSLLM method takes about seconds while the MD-BSLLM method takes about seconds to achieve accurate results. Therefore, we can conclude that it is not necessary to use the MD-BSLLM method to solve the partial differential equations in small subdomains of . The multidomain approach is suitable for solving equations with large values of . The results in Table 5 were obtained with and for the BSLLM method; meanwhile , , and were used for the MD-BSLLM method for .
Table 6 shows the effect of the number of subintervals . Table 6 shows that increasing the number of subintervals yields more accurate results. We note that when , that is, we have one domain, then the MD-BSLLM is behaves exactly as the BSLLM method. For , the results do not match because the MD-BSLLM is exactly the BSLLM method. However, for , only ten subintervals ensure accuracy of the MD-BSLLM method as depicted in Table 6. These results confirm that decomposing the main domain into subintervals improves the accuracy of the BSLLM method.
Figures 3, 4, and 5 show the velocity, temperature, and concentration profiles for various values of , respectively. It is observed that an increase in leads to a decrease in the velocity, temperature, and concentration. The velocity profile is parabolic and hence difficult to be obtained by certain numerical methods like the Gauss-Seiedel based pseudospectral relaxation method. The temperature and concentration profiles decay exponentially. The velocity, temperature, and concentration profiles are in excellent agreement with the results by Hossain [1]. This implies that the MD-BSLLM method can be used as a numerical method for solving large parameter nonsimilar, nonlinear partial differential equations.
Figures 6, 7, and 8 show the convergence error graphs for the velocity, temperature, and concentration profiles for various values of , respectively. It is observed that the method converges after approximately three iterations to a very small error. This validates the accuracy and convergence of the method.
6. Conclusion
In this article, a new numerical approach for solving systems of nonsimilar nonlinear partial differential equations over large domains is presented and implemented. The method is generalized for any th order system of partial differential equations. The method was tested on two different numerical experiments. The results were validated by using series solution of the numerical experiments. The results of the numerical method and the series solution were in excellent agreement and hence validating the accuracy of the method over large domains. The idea of decomposing the main domain into subintervals increases the accuracy of the method. The more subintervals we have, the more accurate the method is. However, for small domains, it is recommended to use and implement the BSLLM method other than the MD-BSLLM method since the BSLLM method achieves accurate results within few seconds as opposed to the MD-BSLLM method. This project has added to literature a novel approach for solving nonsimilar nonlinear system of partial differential equations over a large domain. In future, we will develop and solve systems of nonlinear partial differential equations arising from other disciplines.
Data Availability
All the data used to support the findings of this study are freely available from the corresponding author upon request.
Conflicts of Interest
The author declares that there are no conflicts of interest regarding the publication of this article.
Acknowledgments
The author would like to acknowledge funding from the DST-NRF Centre of Excellence in Mathematical and Statistical Sciences (CoE-MaSS). However, the opinions expressed and conclusions arrived at are those of the author and are not necessarily to be attributed to the Centre of Excellence in Mathematical and Statistical Sciences. Part of this work was presented at the 40th South African Symposium of Numerical and Applied Mathematics, held at Stellenbosch University, from 22nd of March till the 24th of March 2016 and at the Conference of the Southern Africa Mathematical Sciences Association 2016, which was held in Pretoria from the 21st to the 24th of November 2016.