Abstract

We propose a simple and efficient method for solving highly nonlinear systems of boundary layer flow problems with exponentially decaying profiles. The algorithm of the proposed method is based on an innovative idea of linearizing and decoupling the governing systems of equations and reducing them into a sequence of subsystems of differential equations which are solved using spectral collocation methods. The applicability of the proposed method, hereinafter referred to as the spectral local linearization method (SLLM), is tested on some well-known boundary layer flow equations. The numerical results presented in this investigation indicate that the proposed method, despite being easy to develop and numerically implement, is very robust in that it converges rapidly to yield accurate results and is more efficient in solving very large systems of nonlinear boundary value problems of the similarity variable boundary layer type. The accuracy and numerical stability of the SLLM can further be improved by using successive overrelaxation techniques.

1. Introduction

The important properties and applications of the boundary layer flows in many engineering areas are well documented and have been widely researched in the last few decades. The mathematical model of the boundary layer flow system is described, in general, by a set of nonlinear partial differential equations which cannot be solved exactly. A common approach that is employed in the solution of boundary layer flow problems is the reduction of the governing partial differential equations into a set of ordinary differential equations (ODEs) after using appropriate assumptions and a set of transformations based on the so-called similarity variables and using the so-called stream function formulations. The transformed ODEs with the corresponding boundary conditions becomes a system of highly nonlinear boundary value problems (BVP) which, in most cases, does not have closed form solutions. The solution of these nonlinear BVPs continues to fascinate and inspire researchers to develop methods of obtaining solutions to these equations which elucidate the intricate properties of the underlying boundary layer flow problem under different conditions.

The standard way of solving the transformed similarity variable boundary layer equations is the numerical approach based on the shooting algorithm with the Runge-Kutta scheme [1, 2]. Other numerical approaches that have been preferred by some researchers include the finite difference method [3], Keller-Box method [4, 5], spectral homotopy analysis method [6, 7], element free Galerkin method [8], and Newton-Raphson based methods such as the quasilinearization method of [9, 10] and the successive linearization method [1113]. Some analytical approaches have also been found to be very useful in solving boundary layer equations. Examples include the homotopy analysis method [1416], the homotopy perturbation method [1719], the differential transformation method [2022], the Lie-group shooting method [23], the parameter iteration method [24], and the variational iteration method [25], among many other approaches. However, analytical methods are limited in their applications and can only be used in simple systems with few equations. Most real life applications in fluid mechanics involve multiple interacting physical processes and are modelled using complex multiple equation systems. For such problems, it is impractical to use analytical methods because their solution process becomes too cumbersome, and convergence to the true solution can be very slow or not possible at all. For, this reason numerical methods are by far the most practical way of seeking solutions to the boundary layer flow type of highly nonlinear systems. Numerical solution methods also have their own disadvantages such as the stability and convergence issues, difficulty in dealing with singularities, and limiting cases. In addition, some numerical solutions are difficult to interpret such as the occurrence of multiple solutions. The quest for the most optimal method of solving nonlinear problems in fluid mechanics is what drives, ever growing interest in the development of new methods and the modification and improvement of existing analytical and numerical methods.

The prime objective of this paper is to present a new numerical method of solving boundary layer equations that seeks to address some of the aforementioned numerical difficulties. We propose a very simple, yet very accurate and convergent iterative algorithm for solving nonlinear systems of equations that model boundary layer flow problems. The proposed method, hereinafter referred to as spectral local linearization method (SLLM), is based on, decoupling and linearizing systems of equations using a combination of a univariate linearization technique and a spectral collocation discretization. The key feature of the SLLM algorithm is that it breaks down a large coupled system of equations into a sequence of smaller subsystems which can be solved sequentially in a very computationally efficient manner. The applicability of the proposed method is tested on the well-known Blasius boundary layer problem and a three-equation coupled system that models the problem of unsteady free convective heat and mass transfer which are used as test cases to compare the performance of the proposed SLLM when compared to other existing methods. The computed SLLM results demonstrate that the method is easy to develop, accurate, convergent, stable, and very efficient when compared with other existing methods of solving some large systems of boundary value problems.

2. Description of the Methods of Solution

This section presents a brief description of how the proposed iterative methods of solution are developed for a general system of nonlinear ordinary differential equations in unknown functions.

2.1. Spectral Local Linearization Method (SLLM)

Here, we describe the development of the spectral local linearization method (SLLM). Consider a system of nonlinear ordinary differential equations in unknown functions , where is the dependent variable. The system can be written in terms of as a sum of its linear () and nonlinear components () as

To develop the iteration scheme, we apply local linearization of about (the previous iteration) to the th nonlinear equation assuming that all other are known. Thus, at the th equation, is linearized as follows: Thus, at the current iteration with , (1) becomes where [] denotes and and are the approximations of at the current and the previous iteration, respectively. To obtain a decoupled iteration scheme, we appeal to the Gauss-Seidel approach of decoupling linear algebraic systems in linear algebra applications. We therefore arrange the equations in a particular order and solve them in a chronological order. In seeking the solution of in the current iteration level, , we use updated solutions of () obtained as solutions of the previous equations. Thus, for a system of equations, the local linearization iteration scheme becomes where, at the th equation, .

Thus, starting from an initial approximation , the proposed iterative scheme (4) is then solved as a loop until the system converges at a consistent solution for all the variables. To solve the iteration scheme (4), it is convenient to use the Chebyshev pseudospectral method. For this reason the proposed method is referred to as the spectral local linearization iteration method (SLLM) in this work. Spectral methods are now becoming the preferred tools for solving ordinary and partial differential equations because of their elegance and high accuracy in resolving problems with smooth functions.

For brevity, we omit the details of the spectral methods and refer interested readers to [26, 27]. Before applying the spectral method, it is convenient to transform the domain on which the governing equation is defined to the interval on which the spectral method can be implemented. We use the transformation to map the interval on . The basic idea behind the spectral collocation method is the introduction of a differentiation matrix which is used to approximate the derivatives of the unknown variables at the collocation points as the matrix vector product where is the number of the collocation points (grid 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.

2.2. Spectral Quasilinearization Method (SQLM)

In this section, we describe a quasilinearization method that does not use the decoupling approach described in the previous section. This quasilinearization method (QLM) is a generalisation of the Newton-Raphson method and was first proposed by Bellman and Kalaba [28] for solving nonlinear boundary value problems. To obtain the QLM iteration scheme the nonlinear component of a differential equation is linearized using the multivariable Taylor series expansion as opposed to the local linearization approach used in the previous section. The QLM scheme is solved using the Chebyshev pseudospectral method as described in the previous section. Consequently, the method is referred to as the spectral quasilinearization method (SQLM) in this work.

To develop the SQLM scheme we consider a system of nonlinear ordinary differential equations in unknowns functions , where is the dependent variable. The system can be written as a sum of its linear and nonlinear components as subject to the separated boundary conditions where and are linear operators and and are constants for . Define the vector to be the vector of the derivatives of the variable with respect to the dependent variable , that is, where ,   is the th derivative of with respect to and is the highest derivative order of the variable appearing in the system of equations. In addition, we define and to be the linear and nonlinear operators, respectively, that operate on the for . With these definitions, (7) and (8) can be written as where are the constant coefficients of , and the derivative of () that appears in the th equation for .

Noting that, for each variable , the derivatives in the boundary conditions can at most be one less than the highest derivative of in the governing system (7), we define the vector to be the vector of the derivatives of the variable with respect to the dependent variable from up to , that is, The boundary conditions (8) can be written as where () are the constant coefficients of in the boundary conditions and , are the total number of prescribed boundary conditions at and , respectively. We remark that the sum is equal to the sum of the highest orders of the derivatives corresponding to the dependent variables , that is

Assume that the solution of (10) at the th iteration is . If the solution at the previous iteration is sufficiently close to , the nonlinear component of (10) can be linearized using the one-term Taylor series for multiple variables so that (10) can be approximated as subject to where Equation (14) can be rewritten as

The initial approximation, , required to start the iteration scheme (17) is chosen to be a function that satisfies the boundary conditions (8). Applying the Chebyshev spectral collocation on the recursive iteration scheme (17) gives where , , and are given by respectively, where means denote diagonal matrix. Defining , we can write (18) in matrix form as where , are vectors of size and are matrices. Starting from , the recursive sequence (20) is solved iteratively for .

Thus, the size of the coefficient matrix in (20) is and the column vector on the right-hand side has the dimension .

3. Numerical Experiments

In this section, we discuss the implementation of the SLLM and SQLM approaches on two examples of boundary layer flow systems. We consider the Blasius boundary layer problem and a three-equation system that models the problem of unsteady free convective heat and mass transfer on a stretching surface in a porous medium in the presence of a chemical reaction.

3.1. Blasius Boundary Layer Equation

The Blasius boundary layer flow equation, in a dimensionless form, is expressed in terms of the similarity variable and function as where the primes denote differentiation with respect to . To apply the proposed SLLM iterative method on (21), it is convenient to reduce the order of the governing equations by one through setting . This results in the following system:

We observe that (23) is essentially a differential equation whose unknown is , but the equation requires which can be considered to be a known input function. Thus, assuming that is known at a particular iteration level, the local linearization iteration scheme can be expressed by the following iteration formula;

Applying the Chebyshev pseudospectral method on (24), we obtain the following decoupled system of equations: where where is an vector on which the boundary conditions are imposed and and are the values of functions and , respectively, when evaluated at the collocation points. Equation (25) constitutes the SLLM iteration scheme.

The initial approximation required to start the iteration process can be chosen as a function that satisfies the governing boundary conditions and from known physical considerations of the flow properties. For most boundary layer flow problems it is well known that the flow velocity tends to the mainstream flow exponentially. Thus, a suitable initial approximation for the boundary layer problem (21) is

To obtain the SQLM iterative scheme for the Blasius problem (21), we set Thus, using (17), we obtain subject to

Applying the spectral method, with derivative matrices on the iteration scheme (29) and the corresponding boundary conditions, gives the following matrix system: with the boundary conditions where

3.2. Unsteady Free Convective Heat and Mass Transfer on a Stretching Surface in a Porous Medium with Suction/Injection

In this section we consider a three-equation system that models the problem of unsteady free convective heat and mass transfer on a stretching surface in a porous medium in the presence of a chemical reaction. The governing equations [29, 30] for this problem are given as the following dimensionless system of equations and boundary conditions: where , , and are, respectively, the dimensionless velocity, temperature, and concentration, is the suction/injection parameter, is the chemical reaction constant, is the Prandtl number, is the Schmidt number, is the permeability parameter, and and are the temperature and concentration dependent Grashof numbers, respectively.

To apply the proposed SLLM iterative method on (34), we set to obtain

We observe that the energy equation (38) requires only and to resolve the solution of . The same applies to the mass transfer equation (39) for . Equation (37), on the other hand, requires the solutions of ,   and . In developing the SLLM iteration scheme we solve (38) and (39) for , and first, then use their updated solutions in the (37) for . Thus, assuming that and are known at a particular iteration level, the local linearization iteration algorithm can be expressed as subject to

Applying the Chebyshev pseudospectral method on (40)(41), we obtain the following decoupled matrix system of equations: where where is an identity matrix, is an vector, and , ,  ,  and   are the approximate values of , , ,  and   evaluated at the collocation points. After incorporating the boundary conditions in the matrices in (42), the approximate solutions of ,  ,  , and at each iteration level can be obtained by solving (42) in the order in which they are listed. A suitable initial approximation for the SLLM scheme is

To obtain the SQLM iterative scheme for the system (34), we set

To use formula (17) we define the following variables:

Thus, using (45)(46) in (17), we obtain

Applying the Chebyshev pseudospectral method on (47) we obtain the following SQLM scheme in a matrix form: subject to the boundary conditions where and denotes ; that is, the vector elements are placed on the main diagonal of a matrix whose entries everywhere else are zero. Starting from the initial approximations (44), the approximate SQLM solutions for , and are obtained by solving (49).

4. Numerical Convergence, Error, and Stability Analysis of the Iteration Schemes

The convergence and stability of the iteration schemes can be assessed by considering the norm of the difference in the values of functions between two successive iterations. Thus, for each iteration scheme, we define the following maximum error () at the th iteration: where ;   are the governing unknown functions in the nonlinear system. If the iteration scheme converges, the error is expected to decrease with an increase in the number of iterations. In this paper, the unknowns were calculated, for a given number of collocation points , until the following criteria for convergence was fulfilled at iteration : where is the convergence tolerance level. In this study, the convergence tolerance is set to be . The effect of the number of collocation points was examined in order to select the smallest value of which gives a consistent solution to the error level. This is achieved by repeatedly solving the governing equations using the proposed iteration schemes with different values of until the consistent solution is reached.

5. Improving the Convergence of the Spectral Local Linearization Method (SLLM)

In numerical linear algebra, successive overrelaxation (SOR) methods are used to accelerate the convergence rates of the Gauss-Seidel method in the solution of the linear systems of equation. In this section we propose a similar method to improve the convergence of the spectral relaxation method (SLLM) described in the previous section. If the SLLM scheme for solving the function at the th iteration is then the modified version of the SLLM is defined as where , are matrices and is the convergence controlling relaxation parameter. We observe that when (55) reduces to the original SLLM method. The results in the next section will show that, for certain values of , solving the modified SLLM results in improved efficiency and accuracy.

6. Results and Discussion

We solved the governing systems of (21) and (34) we solved using the spectral local linearization method (SLLM) and spectral quasilinearization method (SQLM) as described in the previous section. In this section, we present the results of the numerical computations for the velocity, temperature, and concentration profiles for various input parameters. Values of other flow properties such as skin friction , surface heat transfer rate at the , and mass transfer rate at the wall are also discussed in order to illustrate some special features of the solution and to compare the accuracy, convergence, and efficiency of the SQLM and SLLM algorithms proposed in this study. The accuracy of the present results was verified by comparing them with other results from literature which have been reported to be accurate to within a certain number of decimal digits. Further validation of the solution was established by comparing the present results with the numerical solutions obtained using the MATLAB inbuilt routine bvp4c which is a finite difference code that implements the three-stage Lobatto III formula [33]. In the calculations presented here, the values of the governing physical parameters were chosen deliberately to match some of the results from published literature in order to enable effective comparison.

In Figure 1 we give a comparison of the SLLM and SQLM for the plot of the Logarithm of the error against the number of iterations in the Blasius boundary layer problem. It can be seen from the figure that the spectral quasi-linearization method converges much faster than the spectral local linearization method. The fast convergence rate of the quasilinearization method was established by many researchers including [9, 10, 3436] who extended the application of the quasilinearization method to a wide variety of nonlinear boundary value problems and established that the method converges quadratically. The convergence rate of the proposed SLLM is seen to be linear in this case. The effect of introducing successive overrelaxation with has the beneficial effect of increasing the convergence speed which, in this case, can be characterized by the steeper line. To obtain the optimal value of that yields the best accuracy, the results obtained using several values of and near were compared. It was found that underrelaxation () slows down the convergence, whereas overrelaxation with improves the convergence of the SLLM iteration scheme. The best results for the Blasius boundary layer problem were obtained with . We remark that the optimal was obtained through numerical experimentation.

To further illustrate how the solution converges with an increase in iterations, we give the SLLM computed results for the skin friction in Tables 1 and 2 (with SOR). The results are compared with the recently reported sixteen decimal digit accurate results of [31, 32]. Ganapol [32] used an algorithm based on a Maclaurin series with a Wynn-epsilon convergence acceleration and analytical continuation to obtain highly accurate skin friction coefficients for the Blasius layer flows. It can be seen from Table 1 that full convergence to 9-digit accurate results is achieved after 17 iterations. The 17-digit accurate results of [31] are matched exactly after 31 iterations. We remark that this value was reached after 52 iterations as reported in the work of  [32]. This shows that the basic SLLM is more efficient than the method used in [32].

Table 2 gives the SOR-improved SLLM results (with ) for the skin friction at different iterations. The table indicates that 9-digit accurate results are obtained after only 13 iterations and 17-digit accurate results are reached after 24 iterations. This clearly shows that using the SOR improves the convergence and efficiency of the basic SLLM.

Table 3 gives the SQLM results for the skin friction of the Blasius flow. It can be seen that the SQLM converges very rapidly compared to the SLLM in this example. Full convergence to the 17-digit accurate results reported in [31] is achieved after only 5 iterations. We remark that the results presented in Tables 13 were generated using and . Thus, under the same conditions, the SQLM is more efficient than the SLLM in this case.

The SLLM Log-error plots for the problem of unsteady free convective heat and mass transfer on a stretching surface are given in Figure 2 for varying input parameters. It can be seen that the Logarithm of the error strictly decreases with an increase in the number of iterations when all the parameters are varied. This demonstrates the convergence and stability of the SLLM in solving the three-equation coupled system (34) for varying physical constants. It can also be seen that the convergence speed increases with an increase in ,  , and but decreases with an increase in .

In Table 4 we give a comparison of the values of and computed using the SLLM and SQLM approaches at fixed values of the governing input parameters. The table shows the effects of the number of the collocation points on the convergence with respect to the benchmark bvp4c numerical results. We remark that the stopping criteria used in both methods are given by (52) where the running of the computations is terminated when the maximum norm of the difference between the values of the governing functions at two successive iterations is less than the prescribed tolerance level (set to be in this work). For the selected values of , the SLLM converged after 8 iterations. Increasing improved the accuracy, and the bvp4c results were matched exactly when for both and . The SQLM algorithm was found converge after only 5 iterations when , but the results did not convergent to the bvp4c results. Increasing the number of collocation points slightly improved the accuracy but full convergence to eight decimal digits was not achieved. Further increase in led to the poor accuracy and convergence. This can be seen in the results corresponding to where convergence was not reached even after 100 iterations. These results indicate that the SQLM is not ideal for solving very large systems of equations. The poor accuracy of the SQLM is caused by the numerical difficulties associated with solving large matrix equations. The advantage of using the SLLM in solving large coupled systems of differential equations is that the SLLM algorithm leads to decoupled systems of sub-matrix equations which can easily be solved using standard direct methods in a sequential manner. With the SQLM the coupled system leads to a coupled iteration scheme which after discretization, using the spectral method, leads to large matrices. The dimension of the coefficient matrices and vectors to be solved in the SQLM algorithm is , where is the number of unknowns in the governing systems. For example, in the case when , the dimension of the coefficient matrix used in the solution for , and is in the SQLM, compared to in the case of the SLLM. Furthermore, the displayed accuracy and stability of the SLLM results for both small and large values of suggest that the efficiency of SLLM is due more to the algorithm behind the iteration scheme than to the specific numerical method used to solve the underlying differential equations of the scheme. The SQLM accuracy and stability are strongly linked to the numerical scheme used to solve the differential equations of the scheme. The Chebyshev differentiation matrix is known to be very susceptible to very large round off errors and ill-conditioning for large numbers of collocation points (see e.g., [27, 37]), that is why the SQLM results in Table 4 seem to be less accurate when becomes very large. Table 4 clearly indicates that the SLLM is more practical to use for large systems of equations than SQLM as it yields more accurate results using only few collocation points, and there is no loss of accuracy when the side of the mesh is refined (when becomes large).

In Figure 3 we give a comparison of the Log-error against iterations plot () between the SLLM and SQLM results for a fixed number of collocation points. From Figure 3 it can be seen that the SQLM rapidly converges but the accuracy does not improve after 5 iterations. After 5 iterations the error plot seems to plateau at a certain level and remains at that level even when considering more iterations. The SLLM results indicate strict convergence with the error progressively reduced with an increase in the number of iterations. This explains why the SLLM gives more accurate results than the SQLM for larger systems as was observed in Table 4.

In Table 5 we give a comparison between the basic SLLM results and the results obtained using the SOR accelerated version of the SLLM for the skin friction and wall heat transfer rate . It can be noted from Table 3 that introducing the convergence accelerating relaxation parameter leads to the improvement of the basic SLLM for some cases. It can be seen that convergence improves with an increase in ,  , and and slows down with an increase in . This result is consistent with the observation made in Figure 2. In cases where convergence is slow, underrelaxation with () accelerates the convergence. We also observe that the effect of increasing the values of results in an increase in the skin friction and in the absolute values of the wall heat transfer rate. The wall heat transfer rate decreases in absolute value as increases. An increase in the unsteadiness parameter and the suction parameter results in an increase in the skin friction and the wall heat transfer rate. We remark that the trends in the effect of the mass transfer rate at the wall were found to be similar to the effects observed in the surface heat transfer rate for all varied input parameters. The trends observed in Table 5 are consistent with the observations made in related studies [29, 30].

Figures 4, 5, 6, and 7 illustrate the effects of the input parameters ,  ,   and , respectively, on the temperature and concentration profiles. It is noted that the velocity decreases with , , and but increases with the buoyancy parameter . The temperature decreases with ,  , and but increases with an increase in . Again, the observed trends are consistent with the observations reported in [29, 30]. The effect of the varied parameters on the concentration was found to be similar to the trends observed in the temperature profiles.

7. Conclusion

In this work we have introduced a new method for solving systems of nonlinear boundary value equations. The proposed method, called spectral local linearization method (SLLM), is based on a simple idea of the decoupling systems of equations using the linearization of the unknown functions in a sequential manner according to the order of the listing of the governing equations of the system. The linearized equations, which now form a sequence of linear differential equations with variable coefficients, were solved using the Chebyshev spectral collocation method. The applicability of the method was examined on similarity boundary layer problems with exponentially decaying profiles. The Blasius boundary layer problem and a three-equation coupled system that models the problem of unsteady, free convective heat and mass transfer were considered as test cases to examine the applicability of the proposed SLLM. The accuracy of the SLLM was confirmed against known accurate results reported in the published literature and against numerical results generated using Matlab’s bvp4c routine for solving boundary value problems. The performance of the SLLM was measured against the spectral quasi-linearization method (SQLM) which uses the Chebyshev collocation method to discretize the quasilinearization method (a Newton-Raphson based iteration method for solving boundary value problem).

From this preliminary study, the following observations were made about the proposed SLLM.(1)The algorithm of the SLLM is very easy to develop and implement, as it is based on a simple univariate linearization of nonlinear functions. (2)The method is numerically efficient, since it results in a series of equations which are solved in a sequential manner by reusing the information from the solution of one equation in the next equation. (3)The method was found to converge rapidly to the expected solutions for all the input parameters considered in this study. (4)The convergence speed of the method can readily be improved by using successive overrelaxation (SRM) techniques. In the case of the Blasius boundary layer, overrelaxation improved convergence, and in the unsteady heat and mass transfer problem, the convergence was improved by underrelaxation. (5)The convergence speed of the SLLM improves with an increase in the values of the parameters ,  and   and slows down with an increase in . (6)The method is more computationally efficient than the SQLM in the sense that it gives accurate results even with a few collocation points (small ). The accuracy of the SQLM was seen to deteriorate when the number of collocation points became very large. This is due to the numerical difficulties associated with solving large matrix systems. Unlike the SQLM, the SLLM does not suffer from the loss of accuracy when the number of collocation points is large.

Because of its demonstrated accuracy, efficiency, and simplicity, it is envisaged that the proposed SLLM method could be used as a viable method for solving some classes of similarity variable boundary layer problems.

Acknowledgment

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