On the Comparison between Compact Finite Difference and Pseudospectral Approaches for Solving Similarity Boundary Layer Problems

Abstract

We introduce two methods based on higher order compact finite differences for solving boundary layer problems. The methods called compact finite difference relaxation method (CFD-RM) and compact finite difference quasilinearization method (CFD-QLM) are an alternative form of the spectral relaxation method (SRM) and spectral quasilinearization method (SQLM). The SRM and SQLM are Chebyshev pseudospectral-based methods which have been successfully used to solve boundary layer problems. The main objective of this paper is to give a comparison of the compact finite difference approach against the pseudo-spectral approach in solving similarity boundary layer problems. In particular, we seek to identify the most accurate and computationally efficient method for solving systems of boundary layer equations in fluid mechanics. The results of the two approaches are comparable in terms of accuracy for small systems of equations. For larger systems of equations, the proposed compact finite difference approaches are more accurate than the spectral-method-based approaches.

1. Introduction

Low-order (second or lower) finite difference schemes are not accurate enough for solving many problems in fluid dynamics and other areas. Recently the focus has shifted to higher-order compact finite difference (CFD) schemes. Researchers have found significant improvement to the accuracy of numerical solutions by using fourth or sixth-order CFD schemes compared to the convectional second order finite central difference scheme [1]. Various CFD schemes used for applications such as interpolation, filtering, and evaluating high-order derivatives were discussed in detail by Lele [2]. CFD schemes have largely been applied to solve partial differential equations, for example, Burger’s equation [3, 4], Navier-Stokes equation [5], Korteweg-de Vries equation [6], Black-Scholes equation [7], and many more [8–10]. A limited number of researchers have utilised the CFD schemes in ordinary differential equations. For example in [11], Zhao solved two-point boundary value problems. In [12], Zhao and Corless used the CFD schemes to solve integro-differential equations.

The advantage of the higher-order CFD schemes is that they give high accuracy on coarser grids with greater computational efficiency [13]. The difficulty that comes with the higher order CFD schemes especially near the boundaries makes many researchers shy away from using them in practical computations. To retain the accuracy of the schemes at the boundaries, the schemes are adjusted for the boundary points.

When compared to spectral methods, compact schemes are more flexible in terms of application to complex geometries and boundary conditions. Lele [2] pointed out that the use of spectral methods in turbulent fluid flows is limited to flows in simple domains and simple boundary conditions. Spectral methods become less accurate for problems with complex geometries. Rai [14] gave a comprehensive comparison between results obtained using finite differences and spectral methods for direct simulation of turbulent flows. They used high-order accurate upwind schemes. They concluded that the spectral method is extremely accurate but it has restrictions on the type of geometry and grids that can be efficiently handled. They further state that unlike the spectral methods, the finite difference method can be used efficiently with curvilinear grids.

In recent years Motsa and his coworkers [15–18] have developed successful methods based on the spectral method to solve nonlinear similarity boundary layer problems. The methods include, among others, the spectral relaxation method (SRM) [15, 19, 20], the spectral successive linearisation method [21–23], spectral homotopy analysis method [24–27], and the spectral quasilinearisation method (SQLM) [16]. The SRM is based on simple decoupling and rearrangement of the governing equations and numerically integrating the resulting equations using the Chebyshev spectral collocation method. The SQLM combines the quasilinearization technique developed by Bellman and Kalaba [28] to linearize nonlinear differential equations and solve the resulting equations using the spectral method. In this work we present alternative approach to the implementation of the SRM and SQLM. Instead of using the spectral method in these methods, we use the higher order CFD schemes, and hence the resulting methods are the compact finite difference relaxation method (CFD-RM) and compact finite difference quasilinearization method (CFD-QLM), respectively.

The main objective of this work is to compare the spectral-method-based and the CFD-based methods discussed above. We compare the performance of the methods in terms of accuracy and computational speed when solving nonlinear boundary layer problems in one dimension and three dimensions. We first consider the flow of a viscous incompressible electrically conducting fluid over a continuously shrinking sheet which is governed by a third-order nonlinear differential equation. The known exact analytical solution [29, 30] of this problem is used as a benchmark to validate the accuracy of the proposed algorithm discussed in this work. We also consider a 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 [23, 31].

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 Relaxation Method

The spectral relaxation method (SRM) is a new method that has been introduced recently by Motsa et al. [15, 19, 20] to solve initial and boundary value problems. The method is based on simple decoupling and rearrangement of the governing equations and numerically integrating the resulting equations using the Chebyshev spectral collocation method. The algorithm for the method is summarized as follows.(1)Arrange the governing nonlinear equations in a particular order, placing the equations with the least unknowns at the top of the equations list. (2)Assign the labels 𝑍₁,𝑍₂,𝑍₃,… to the ordered equations obtained in the above step, where each 𝑍_𝑖(𝑖=1,2,3,…) is an unknown function which, in the 𝑖th equation, is identified as the unknown function associated with the highest-order derivative. (3)In the equation for 𝑍₁ (1st equation), the iteration scheme is developed by assuming that only linear terms in 𝑍₁ are to be evaluated at the current iteration level (denoted by 𝑟+1) and all other terms (linear and nonlinear) in 𝑍₂,𝑍₃,… are assumed to be known from the previous iteration (denoted by 𝑟). In addition nonlinear terms in 𝑍₁ are also evaluated at the previous iteration. Furthermore, all derivative terms in 𝑍₁ are assumed to be known from the previous iteration. (4)Similarly, in developing the iteration scheme in the equation for 𝑍₂ (2nd equation), only linear terms in 𝑍₂ are evaluated at the current iteration level (𝑟+1) with all other terms evaluated at the previous level, except 𝑍₁ which is now known from the solution of the first equation. (5)This process is repeated in the 𝑖th equation (𝑖=3,4,…) using the updated solutions for 𝑍_𝑖−1 obtained from the previous 𝑖−1 equations. The resulting iteration scheme is integrated using the Chebyshev spectral method. The region of integration is discretized using the Gauss-Lobatto points defined by 𝜏_𝑗=cos𝜋𝑗𝑁,𝑗=1,2,3,…,𝑁, (1) where 𝑁 is the number of collocation points used. The Chebyshev spectral collocation method is based on the idea of introducing a differentiation matrix 𝐷 which is used to approximate the derivatives of the unknown variables 𝑍_𝑖(𝑥)(𝑖=1,2,3,…) at the collocation points as the matrix vector product 𝑑𝑍_𝑖=𝑑𝑥^𝑁_𝑘=0𝐃_𝑗𝑘𝑍_𝑖=𝐃𝐙_𝑖,𝑗=1,2,…,𝑁, (2) where 𝐙_𝑖=[𝑍_𝑖(𝜏₀),𝑍_𝑖(𝜏₁),…,𝑍_𝑖(𝜏_𝑁)] is the vector function at the collocation points 𝜏_𝑗.

2.2. Compact Finite Difference Relaxation Method

The compact finite difference relaxation method (CFD-RM) uses the same procedure followed in the SRM. The difference is that instead of using spectral methods to solve the resulting iteration schemes, higher order compact finite difference schemes are used. In this work we use sixth-order CFD schemes.

In the derivation of the CFD schemes we consider a one-dimensional uniform mesh on the region [𝑎,𝑏] with nodes 𝑥_𝑖 (𝑖=1,2,…,𝑁) where 𝑎=𝑥₁<𝑥₂<⋯<𝑥_𝑁=𝑏 (3) and a corresponding function 𝑦_𝑖=𝑦(𝑥_𝑖) at the nodes. The distance between any two successive nodes is a constant ℎ=𝑥_𝑖−𝑥_𝑖−1. Sixth-order approximations of the first, second, and third derivatives at interior nodes can be obtained using the following schemes (see [2] for details): 13𝑦^′_𝑖−1+𝑦^′_𝑖+13𝑦^′_𝑖+1=149𝑦_𝑖+1−𝑦_𝑖−1+12ℎ9𝑦_𝑖+2−𝑦_𝑖−2,24ℎ𝑦11^{′′𝑖−1}+𝑦_𝑖^′′+2𝑦11^′′𝑖+1=12𝑦11^𝑖+1−2𝑦_𝑖+𝑦_𝑖−1ℎ²+3𝑦11^𝑖+2−2𝑦_𝑖+𝑦_𝑖−24ℎ²,7𝑦16^{′′′𝑖−1}+𝑦_𝑖^′′′+7𝑦16^{′′′𝑖+1}𝑦=2^𝑖+2−2𝑦_𝑖+1−2𝑦_𝑖−1−𝑦_𝑖−22ℎ³−18𝑦_𝑖+3−3𝑦_𝑖+1−3𝑦_𝑖−1−𝑦_𝑖−38ℎ³. (4)(5)(6)

For illustrative purposes we describe the application of the CFD schemes to second-order differential equations for 𝑦(𝑥) with known boundary conditions at 𝑦(𝑎) and 𝑦(𝑏). Consider the nonlinear differential equations 𝑦^′′+𝑝(𝑥)𝑦^′+𝑞(𝑥)𝑦(𝑥)+𝑓𝑥,𝑦,𝑦^′𝑦(𝑎)=0,=𝑦_𝑎,𝑦(𝑏)=𝑦_𝑏, (7) where 𝑓(𝑥,𝑦,𝑦^′) is a nonlinear function, 𝑝(𝑥) and 𝑞(𝑥) are known functions of 𝑥, and 𝑦_𝑎 and 𝑦_𝑏 are known constants. In solving (7), we apply the CFD approximation for the first and second derivatives given by (4) and (5), respectively, at the interior nodes (𝑖=2,…,𝑁−1). Since we know boundary conditions at 𝑖=1 and 𝑖=𝑁, the CFD schemes must be adjusted for the nodes near the boundary points. In order to maintain the order 𝑂(ℎ⁶) accuracy at the boundary points as in the interior points and to maintain the same tridiagonal format, we use the following one-sided scheme at 𝑖=2: 𝑦^′₂+13𝑦^′₃=1ℎ𝑎₁𝑦₁+𝑎₂𝑦₂+𝑎₃𝑦₃+𝑎₄𝑦₄+𝑎₅𝑦₅+𝑎₆𝑦₆+𝑎₇𝑦₇, (8) and when 𝑖=𝑁−1, we use 13𝑦^′_𝑁−2+𝑦^′_𝑁−1=1ℎ𝑏₁𝑦_𝑁+𝑏₂𝑦_𝑁−1+𝑏₃𝑦_𝑁−3+𝑏₄𝑦_𝑁−4+𝑏₅𝑦_𝑁−5+𝑏₆𝑦_𝑁−6+𝑏₇𝑦_𝑁−7, (9) where 𝑎_𝑖,𝑏_𝑖(𝑖=1,…,7) are constants to be determined. To obtain a sixth-order accurate scheme, we use Taylor series expansion about 𝑥₁ and 𝑥_𝑁 up to 𝑂(ℎ⁷), in (8) and (9), respectively, and equate terms of order ℎ. In each case, we obtain a system of seven linear algebraic equations in seven unknowns which are solved to give 𝑎₁7=−45,𝑎₂=−1712,𝑎₃=83,𝑎36₄=−119,𝑎₅=23,𝑎₆=−37180,𝑎₇=1,𝑏36_𝑘=−𝑎_𝑘,𝑘=1,2,…,7. (10)

Similarly, for the second derivatives, we use 𝑦₂^′′+2𝑦11^3′′=1ℎ²𝑐₁𝑦₁+𝑐₂𝑦₂+𝑐₃𝑦₃+𝑐₄𝑦₄+𝑐₅𝑦₅+𝑐₆𝑦₆+𝑐₇𝑦₇+𝑐₈𝑦₈, (11) at 𝑖=2 and 2𝑦11^{′′𝑁−2}+𝑦^′′_𝑁−1=1ℎ²𝑑₁𝑦_𝑁+𝑑₂𝑦_𝑁−1+𝑑₃𝑦_𝑁−3+𝑑₄𝑦_𝑁−4+𝑑₅𝑦_𝑁−5+𝑑₆𝑦_𝑁−6+𝑑₇𝑦_𝑁−7+𝑑₈𝑦_𝑁−8 (12) at 𝑖=𝑁−1. The parameters 𝑐_𝑖,𝑑_𝑖 can be determined by expanding (11) and (12) using Taylor series and equating powers of ℎ and subsequently solving the resulting equations. This gives 𝑐₁=3145,𝑐₂=−19110,𝑐₃=−339,𝑐110₄=1933396,𝑐₅=−40,𝑐11₆=9655,𝑐₇=−479990,𝑐₈=13,𝑑220_𝑖=𝑐_𝑖,𝑖=1,2,…,8. (13)

Using the above equations, the equations for approximating the first- and second-order derivatives can be expressed as 𝐴₂𝑌^′′=𝐵₂𝑌+𝐾₂,𝐴₁𝑌^′=𝐵₁𝑌+𝐾₁, (14) where 𝐴₁=11313113131131⋱⋱⋱3113131_{(𝑁−2)×(𝑁−2)},𝐾₁=1ℎ𝑎₁𝑦₁−𝑦₁0⋮⋮0𝑦36_𝑁36−𝑎₁𝑦_{𝑁(𝑁−2)×1},𝐵₁=1ℎ𝑎₂𝑎₃𝑎₄𝑎₅𝑎₆𝑎₇−79079136−1−7369079136⋱⋱⋱⋱⋱−1−7369079136−1−7369079−𝑎₇−𝑎₆−𝑎₅−𝑎₄−𝑎₃−𝑎_{2(𝑁−2)×(𝑁−2)},𝐴₂=1221112112111211211⋱⋱⋱1211211111_{(𝑁−2)×(𝑁−2)},𝐾₂=1ℎ²𝑐₁𝑦₁3𝑦₁0⋮⋮0443𝑦_𝑁𝑐44₁𝑦_{𝑁(𝑁−2)×1},𝐵₂=1ℎ²𝑐₂𝑐₃𝑐₄𝑐₅𝑐₆𝑐₇𝑐₈12−115122123113444412−11512212311344⋱⋱⋱⋱⋱4412−115122123113444412−11512212𝑐11₈𝑐₇𝑐₆𝑐₅𝑐₄𝑐₃𝑐_{2(𝑁−2)×(𝑁−2)}𝑌^′=𝑦^′₂,𝑦^′₃,…,𝑦^′_𝑁−2,𝑦^′_𝑁−1^𝑇,𝑌^′′=𝑦₂^′′,𝑦₃^′′,…,𝑦^′′_𝑁−2,𝑦^′′_𝑁−1^𝑇. (15)(16)(17)

For the purposes of the examples given in this work, we also illustrate the application of the CFD schemes for the following set of boundary conditions: 𝑦(𝑎)=𝑦_𝑎,𝑦^′(𝑎)=𝑦_𝑑𝑎,𝑦^′(𝑏)=𝑦_𝑑𝑏. (18)

We adjust the schemes near the boundary points as discussed above. In this case, the schemes at 𝑖=𝑁−2 are also adjusted. At 𝑖=2, we use the following one-sided scheme: 𝑦^′₂+13𝑦^′₃=𝑎₁𝑦^′₁+1ℎ𝑎₂𝑦₁+𝑎₃𝑦₂+𝑎₄𝑦₃+𝑎₅𝑦₄+𝑎₆𝑦₅+𝑎₇𝑦₆. (19) At 𝑖=𝑁−1, we use 13𝑦^′_𝑁−2+𝑦^′_𝑁−1=𝑏₁𝑦^′_𝑁+1ℎ𝑏₂𝑦_𝑁−1+𝑏₃𝑦_𝑁−2+𝑏₄𝑦_𝑁−3+𝑏₅𝑦_𝑁−4+𝑏₆𝑦_𝑁−5+𝑏₇𝑦_𝑁−6, (20) and at 𝑖=𝑁−2, we use 13𝑦^′_𝑁−3+𝑦^′_𝑁−2+13𝑦^′_𝑁−1=𝑐₁𝑦^′_𝑁+1ℎ𝑐₂𝑦_𝑁−1+𝑐₃𝑦_𝑁−2+𝑐₄𝑦_𝑁−3+𝑐₅𝑦_𝑁−4+𝑐₆𝑦_𝑁−5+𝑐₇𝑦_𝑁−6 (21) with the constants given by 𝑎₁1=−6,𝑎₂=−203360,𝑎₃5=−,𝑎12₄=1918,𝑎₅1=−9,𝑎₆=1,𝑎24₇1=−,𝑏80₁=463,𝑏₂=15184,𝑏₃=−701,𝑏252₄=311189,𝑏₅=−1921,𝑏₆=71,𝑏252₇=−29,𝑐756₁=5441,𝑐₂=373441,𝑐₃=−25,𝑐294₄=−9291323,𝑐₅=−31441,𝑐₆2=−,𝑐147₇5=−.2646 (22)

Similarly for the second derivative, we use 𝑦₂^′′+2𝑦11^3′′=𝑑₁𝑦^′₁ℎ+1ℎ²𝑑₂𝑦₁+𝑑₃𝑦₂+𝑑₄𝑦₃+𝑑₅𝑦₄+𝑑₆𝑦₅+𝑑₇𝑦₆+𝑑₈𝑦₇, (23) at 𝑖=2. At 𝑖=𝑁−1, we use 2𝑦11^{′′𝑁−2}+𝑦^′′_𝑁−1=𝑒₁𝑦^′_𝑁ℎ+1ℎ²𝑒₂𝑦_𝑁−1+𝑒₃𝑦_𝑁−2+𝑒₄𝑦_𝑁−3+𝑒₅𝑦_𝑁−4+𝑒₆𝑦_𝑁−5+𝑒₇𝑦_𝑁−6+𝑒₈𝑦_𝑁−7, (24) and at 𝑁−3, 2𝑦11^{′′𝑁−3}+𝑦^′′_𝑁−2+2𝑦11^{′′𝑁−1}=𝑓₁𝑦^′_𝑁ℎ+1ℎ²𝑓₂𝑦_𝑁−1+𝑓₃𝑦_𝑁−2+𝑓₄𝑦_𝑁−3+𝑓₅𝑦_𝑁−4+𝑓₆𝑦_𝑁−5+𝑓₇𝑦_𝑁−6+𝑓₈𝑦_𝑁−7. (25)

The constants are found to be 𝑑₁=91220,𝑑₂=63413600,𝑑₃=−135,𝑑44₄=11188,𝑑₅=118,𝑑₆3=−,𝑑176₇=91100,𝑑₈1=−,𝑒792₁=1759124352,𝑒₂=5511732670,𝑒₃=−63941,𝑒10890₄=31288739204,𝑒₅=−194753267,𝑒₆=3116,𝑒1089₇=−7780199010,𝑒₈=6341,𝑓65340₁=19256088,𝑓₂=16971331,𝑓₃=−3453,𝑓1331₄=55813993,𝑓₅=−4312662,𝑓₆=147,𝑓1331₇=−2457986,𝑓₈=5.1331 (26)

Combining the schemes for approximating the first and second derivatives at interior points with (19)–(23), the equations for approximating the first and second derivatives are given by 𝐴₁𝑌^′=𝐵₁𝑌+𝐾₁+𝐾^′₁,𝐴₂𝑌^′′=𝐵₂𝑌+𝐾₂+𝐾^′₂, (27) where 𝐴₁ and 𝐴₂ are the same as in (14), 𝐵₁=1ℎ𝑎₃𝑎₄𝑎₅𝑎₆𝑎₇−79079136−1−7369079136⋱⋱⋱⋱⋱−1−73690791𝑐36₇𝑐₆𝑐₅𝑐₄𝑐₃𝑐₁𝑏₇𝑏₆𝑏₅𝑏₄𝑏₃𝑏_{2(𝑁−2)×(𝑁−2)},𝐾₁=1ℎ𝑎₂𝑦₁−𝑦₁0⋮036_{(𝑁−2)×1},𝐾^′₁=𝑎₁𝑦^′₁0⋮𝑐₁𝑦^′_𝑁𝑏₁𝑦^′_{𝑁(𝑁−2)×1},𝐵₂=1ℎ²𝑑₃𝑑₄𝑑₅𝑑₆𝑑₇𝑑₈12−115122123113444412−11512212311344⋱⋱⋱⋱⋱4412−11512212311𝑓44₈𝑓₇𝑓₆𝑓₅𝑓₄𝑓₃𝑓₂𝑒₈𝑒₇𝑒₆𝑒₅𝑒₄𝑒₃𝑒_{2(𝑁−2)×(𝑁−2)},𝐾₂=1ℎ²𝑐₂𝑦₁3𝑦₁0⋮044_{(𝑁−2)×1},𝐾^′₂=1ℎ𝑐₁𝑦^′₁0⋮𝑓₁𝑦^′_𝑁𝑒₁𝑦^′_{𝑁(𝑁−2)×1}. (28)

Similarly, third derivatives can be expressed in the form 𝐴₃𝑌^′′′=𝐵₃𝑌+𝐾₃+𝐾^′₃. (29)

2.3. Spectral Quasilinearization Method

In this section we give a brief description of the spectral quasilinearization method (SQLM). In this method we make use of the quasilinearization method (QLM) and spectral method. The QLM was initially proposed by Bellman and Kalaba [28] to solve nonlinear boundary layer problems. To develop the SQLM, we consider a system of 𝑚 nonlinear differential equations in 𝑚 unknowns 𝑧_𝑖(𝜂)(𝑖=1,2,…,𝑚) where 𝜂 is the independent variable. The system can be written as a sum of its linear 𝐿 and nonlinear components 𝑁 as 𝐿𝑧₁𝜂,𝑧₂𝜂,…,𝑧_𝑚𝜂+𝑁𝑧₁𝜂,𝑧₂𝜂,…,𝑧_𝑚𝜂=𝐻𝜂. (30) Define the vector 𝑍_𝑖 to be the vector of the derivatives of the variable 𝑧_𝑖 with respect to the independent variable 𝜂; that is, 𝑍_𝑖=𝑧_𝑖⁽⁰⁾,𝑧_𝑖⁽¹⁾,…,𝑧^(𝑛𝑖)_𝑖, (31) where 𝑧_𝑖⁽⁰⁾=𝑧_𝑖, 𝑧_𝑖^(𝑝) is the 𝑝th derivative of 𝑧_𝑖 with respect to 𝜂, and 𝑛_𝑖(𝑖=1,2,…,𝑚) 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 𝑖=1,2,…,𝑚. With these definitions (30) can be written as 𝐿_𝑖𝑍₁,𝑍₂,…,𝑍_𝑚+𝑁_𝑖𝑍₁,𝑍₂,…,𝑍_𝑚=^𝑚𝑛_𝑗=1^𝑖_𝑝=0𝛼^[𝑝]_𝑖𝑗𝑍_𝑗^(𝑝)+𝑁_𝑖𝑍₁,𝑍₂,…,𝑍_𝑚=𝐻_𝑖, (32) where 𝛼^[𝑝]_𝑖,𝑗 are the constant coefficients of 𝑧_𝑗^(𝑝), the derivative of 𝑧_𝑗(𝑗=1,2,…,𝑚) that appears in the 𝑖th equation for 𝑖=1,2,…,𝑚. Assume that the solution 𝑧_𝑖(𝜂) of (32) at the (𝑟+1)th iteration is 𝑧_𝑖,𝑟+1. If the solution at the previous iteration 𝑧_𝑖,𝑟(𝜂) is sufficiently close to 𝑧_𝑖,𝑟+1, the nonlinear component 𝑁_𝑖 of (32) can be linearised using one-term Taylor series for multiple variables, so that (32) can be approximated as 𝐿_𝑖𝑍_1,𝑟+1,…,𝑍_𝑚,𝑟+1+𝑁_𝑖[⋯]+^𝑚𝑛_𝑗=1^𝑖_𝑝=0𝑍^(𝑝)_𝑗,𝑟+1−𝑍^(𝑝)_𝑗,𝑟𝜕𝑁_𝑖𝜕𝑍_𝑗^(𝑝)[⋯]=𝐻_𝑖, (33) where [⋯]=𝑍_1,𝑟,𝑍_2,𝑟,…,𝑍_𝑚,𝑟. (34)

Equation (33) can be rewritten as 𝐿_𝑖𝑍_1,𝑟+1,…,𝑍_𝑚,𝑟+1+^𝑚𝑛_𝑗=1^𝑖_𝑝=0𝑍^(𝑝)_𝑗,𝑟+1𝜕𝑁_𝑖𝜕𝑍_𝑗^(𝑝)[⋯]=𝐻_𝑖+^𝑚𝑛_𝑗=1^𝑖_𝑝=0𝑍^(𝑝)_𝑗,𝑟𝜕𝑁_𝑖𝜕𝑍_𝑗^(𝑝)[⋯]−𝑁_𝑖[⋯]. (35)

To solve the iteration scheme (35), we use the Chebyshev spectral collocation method.

2.4. Compact Finite Difference Quasilinearization Method

Using compact finite difference schemes to solve the QLM scheme (35) results in the compact finite difference quasilinearization method (CFD-QLM). We will use the sixth-order finite difference schemes the same way we used them for the CFD-RM.

3. Examples

In this section, two numerical examples are discussed to compare the CFD results against spectral method results. We consider one-dimensional and three-dimensional problems.

3.1. MHD Boundary Layer Flow over a Shrinking Sheet

We first consider a steady one-dimensional laminar flow of a viscous incompressible electrically conducting fluid over a continuously shrinking sheet. The governing equation is given (see, e.g., [29, 30]) in similarity form as 𝑓^′′′+𝑓𝑓^′′−𝑓^′2−𝑀²𝑓^′=0 (36) subject to 𝑓(0)=0,𝑓^′(0)=−1,𝑓^′(∞)=0, (37) where 𝑓=𝑓(𝜂) and 𝑀 is the magnetic interaction parameter. The analytical solution for (36) is given by 𝑓𝜂=1𝛼𝑒^−𝛼𝜂−1,𝛼=𝑀²−1. (38)

We solve (36) using the SRM and CFD-RM and the SQLM and CFD-QLM. In the next subsections we explain the development of each of the methods for the solution of (36).

3.1.1. SRM and CFD-RM

To apply the SRM and CFD-RM on (36), we first set 𝑓^′(𝜂)=𝑔(𝜂) and write the equation as the following system of equations: 𝑓^′𝑔=𝑔,_′′+𝑓𝑔^′−𝑔²−𝑀²𝑔=0 (39) with the boundary conditions 𝑓(0)=0,𝑔(0)=−1,𝑔(∞)=0. (40)

Applying the SRM and CFD-RM on (36), we obtain the following iteration scheme: 𝑓^′_𝑟+1=𝑔_𝑟,𝑓_𝑟+1(0)𝑔=0,_′′𝑟+1+𝑓_𝑟+1𝑔^′_𝑟+1−𝑀²𝑔_𝑟+1=𝑔²_𝑟,𝑔_𝑟+1(0)=−1,𝑔_𝑟+1(∞)=0. (41)(42) The SRM and CFD-RM schemes that result from applying the spectral method and compact finite difference schemes on (41)-(42), respectively, are given by 𝑋₁𝑓_𝑟+1=𝑌₁,𝑓_𝑟+1𝜏_𝑁𝑋=0,₂𝑔_𝑟+1=𝑌₂,𝑔_𝑟+1𝜏_𝑁𝑔=−1,_𝑟+1𝜏₀=0, (43) where for the SRM, 𝑋₁=𝐃,𝑌₁=𝑔_𝑟,𝑋₂=𝐃²+diag𝑓_𝑟+1𝐃−𝑀²𝐈,𝑌₂=𝑔²_𝑟, (44) where 𝐈 is an (𝑁+1)×(𝑁+1) identity matrix and diag is a function that places a vector [] on the main diagonal of an (𝑁+1)×(𝑁+1) matrix of zeros. For the CFD-RM, we obtain 𝑋₁=𝐸₁,𝑌₁=𝑔_𝑟−𝐴₁⁻¹𝐾₁,𝑋₂=𝐸₂+diag𝑓_𝑟+1𝐸₁−𝑀²𝑌𝐼,₂=𝑔²_𝑟−𝐴₂⁻¹𝐾₂−𝐴₁⁻¹𝐾₁, (45) where 𝐸₁=𝐴₁⁻¹𝐵₁ and 𝐸₂=𝐴₂⁻¹𝐵₂. 𝐼 is an (𝑁−1)×(𝑁−1) identity matrix.

3.1.2. SQLM and CFD-QLM

To solve (36) using the SQLM and CFD-QLM, we first linearize the equation using the quasilinearization method. We first reduce the order of the equation by setting 𝑓^′(𝜂)=𝑔(𝜂) and write it as the following system of equations: 𝑓^′𝑔=𝑔,_′′+𝑓𝑔^′−𝑔²−𝑀²𝑔=0 (46) with the boundary conditions 𝑓(0)=0,𝑔(0)=−1,𝑔(∞)=0. (47)

Applying the quasilinearization method on (46), we obtain 𝑓^′_𝑟+1−𝑔_𝑟+1𝑔=0,_′𝑟𝑓_𝑟+1+𝑔^′′_𝑟+1+𝑓_𝑟𝑔^′_𝑟+1−2𝑔_𝑟+𝑀²𝑔_𝑟+1=𝑓_𝑟𝑔^′_𝑟−𝑔²_𝑟 (48)(49) subject to 𝑓_𝑟+1(0)=0,𝑔_𝑟+1(0)=−1,𝑔_𝑟+1(0)=0. (50)

The SQLM and CFD-QLM schemes that result from applying the spectral method and compact finite difference schemes on (48)-(49), respectively, are given by Δ_1,1Δ_1,2Δ_2,1Δ_2,2𝑓_𝑟+1𝑔_𝑟+1=Φ_1,𝑟Φ_2,𝑟, (51) where in the framework of the SQLM, Δ_1,1Δ=𝐃,_1,2Δ=−𝐈,_2,1=diag𝑔^′_𝑟,Δ_2,2=𝐃²+diag𝑓_𝑟𝐃−diag2𝑔_𝑟+𝑀²,Φ_1,𝑟Φ=𝐎,_2,𝑟=𝑓_𝑟𝑔^′_𝑟−𝑔²_𝑟, (52) where O is an (𝑁+1)×1 zero vector. For the CFD-QLM, we have Δ_1,1=𝐸₁,Δ_1,2Δ=−𝐼,_2,1=diag𝑔^′_𝑟,Δ_2,2=𝐸₂+diag𝑓_𝑟𝐸₁−diag2𝑔_𝑟+𝑀²,Φ_1,𝑟=−𝐴₁⁻¹𝐾₁,Φ_2,𝑟=𝑓_𝑟𝑔^′_𝑟−𝑔²_𝑟−𝐴₂⁻¹𝐾₂−diag𝑓_𝑟𝐴₁⁻¹𝐾₁. (53)

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 [23, 31] for this problem are given as the following dimensionless system of equations and boundary conditions: 𝑓^′′′+𝑓𝑓^′′−𝑓^′2−𝐾𝑓^′−𝐴𝑓^′+𝜂2𝑓^′′1+Gr𝜃+Gc𝜙=0,𝜃Pr^′′−𝑓^′𝜃+𝑓𝜃^′1−𝐴𝜃+2𝜂𝜃^′1=0,𝜙Sc^′′−𝑓^′𝜙+𝑓𝜙^′1−𝐴𝜙+2𝜂𝜙^′−𝛾𝜙=0 (54) subject to 𝑓(0)=𝑓_𝑤,𝑓^′(0)𝜃(0)(0)𝑓=1,=1,𝜙=1,_′(∞)=0,𝜃(∞)=0,𝜙(∞)=0, (55) where 𝑓(𝜂), 𝜃(𝜂), and 𝜙(𝜂) are, respectively, the dimensionless velocity, temperature, and concentration, 𝑓_𝑤 is the suction/injection parameter, 𝛾 is the chemical reaction constant, Pr is the Prandtl number, Sc is the Schmidt number, 𝐾 is the permeability parameter, and Gr and Gc are the temperature- and concentration-dependent Grashof numbers, respectively.

3.2.1. SRM and CFD-RM

To apply the SRM and CFD-RM on (54), we set 𝑓^′(𝜂)=𝑔(𝜂) and obtain the following iteration scheme: 𝑓^′_𝑟+1=𝑔_𝑟,𝑔^′′𝑟+1+𝑓_𝑟+1𝑔^′_𝑟+1−(𝐴+𝐾)𝑔_𝑟+1+𝜂2𝑔^′_𝑟+1=𝑔²_𝑟−𝐺_𝑟𝜃_𝑟−𝐺_𝑐𝜙_𝑟,1𝑃_𝑟𝜃_𝑟+1−𝑔_𝑟+1𝜃_𝑟+1+𝑓_𝑟+1𝜃^′_𝑟+1−𝐴𝜃_𝑟+1+𝜂2𝜃_𝑟+11=0,𝑆_𝑐𝜙_𝑟+1−𝑔_𝑟+1𝜙_𝑟+1+𝑓_𝑟+1𝜙^′_𝑟+1−𝐴+𝛾𝜙_𝑟+1𝜂−𝐴2𝜙^′_𝑟+1=0 (56) subject to 𝑓_𝑟+1(0)=𝑓_𝑤,𝑔_𝑟+1(0)𝑔=1,_𝑟+1(∞)=0,𝜃_𝑟+1(0)𝜃=1,_𝑟+1(∞)=0,𝜙_𝑟+1(0)𝜙=1,_𝑟+1(∞)=0. (57)

Applying the Chebyshev pseudo-spectral method and finite difference schemes, we obtain 𝑋₁𝑓_𝑟+1=𝑌₁,𝑓_𝑟+1(0)=𝑓_𝑤,𝑋₂𝑔_𝑟+1=𝑌₂,𝑔_𝑟+1𝜏₀=1,𝑔_𝑟+1𝜏_𝑁𝑋=0,₃𝜃_𝑟+1=𝑌₃,𝜃_𝑟+1𝜏₀=1,𝜃𝜏_𝑁𝑋=0,₄𝜙_𝑟+1=𝑌₄,𝜙_𝑟+1𝜏₀=1,𝜙𝜏_𝑁=0. (58)

For the SRM, we have 𝑋₁=𝐃,𝑌₁=𝑔_𝑟,𝑋₂=𝐃²+diag𝑓_𝑟+1𝜂−𝐴2()𝑌𝐃−𝐾+𝐴𝐈,₂=𝑔²_𝑟−𝐺_𝑟𝜃_𝑟−𝐺_𝑐𝜙_𝑟,𝑋₃=1𝑃_𝑟𝐃²+diag𝑓_𝑟+1𝜂−𝐴2𝐃−diag𝑔_𝑟+1𝑌−𝐴𝐈,₃𝑋=𝐎,₄=1𝑆_𝑐𝐃²+diag𝑓_𝑟+1𝜂−𝐴2𝐃−diag𝑔_𝑟+1−𝑌𝐴+𝛾𝐈,₄=𝐎, (59) where O is an (𝑁+1)×1 zero vector.

For the CFD-RM, we have 𝑋₁=𝐸₁,𝑋₂=𝐸₂+diag𝑓_𝑟+1𝜂−𝐴2𝐸₁−()𝑋𝐾+𝐴𝐼,₃=1𝑃_𝑟𝐸₂+diag𝑓_𝑟+1𝜂−𝐴2𝐸₁−diag𝑔_𝑟+1𝑋−𝐴𝐼,₄=1𝑆_𝑐𝐸₂+diag𝑓_𝑟+1𝜂−𝐴2𝐸₁−diag𝑔_𝑟+1−𝑌𝐴+𝛾𝐼,₁=𝑔_𝑟−𝐴₁⁻¹𝐾₁,𝑌₂=𝑔²_𝑟−𝐺_𝑟𝜃_𝑟−𝐺_𝑐𝜙_𝑟−𝐴₂⁻¹𝐾₂−diag𝑓_𝑟+1𝜂−𝐴2𝐴₁⁻¹𝐾₁,𝑌₃1=−𝑃_𝑟𝐴₂⁻¹𝐾₂−diag𝑓_𝑟+1𝜂−𝐴2𝐴₁⁻¹𝐾₁,𝑌₄1=−𝑆_𝑐𝐴₂⁻¹𝐾₂−diag𝑓_𝑟+1𝜂−𝐴2𝐴₁⁻¹𝐾₁. (60)

3.2.2. SQLM and CFD-QLM

To apply the SQLM and CFD-QLM on (54), we first apply the quasilinearization method and obtain the following iteration scheme: 𝑓^′′′_𝑟+1+𝑓_𝑟𝜂−𝐴2𝑓^′′_𝑟+1−2𝑓^′_𝑟+𝐾+𝐴𝑓^′_𝑟+1+𝑓_𝑟^′′𝑓_𝑟+1+𝐺_𝑟𝜃_𝑟+1+𝐺_𝑐𝜙_𝑟+1=𝑓_𝑟𝑓_𝑟^′′−𝑓_𝑟^′2,1𝑃_𝑟𝜃^′′_𝑟+1+𝑓_𝑟𝜂−𝐴2𝜃^′_𝑟+1−𝑓^′_𝑟+𝐴𝜃_𝑟+1−𝜃_𝑟𝑓^′_𝑟+1+𝜃^′_𝑟𝑓_𝑟+1=𝜃^′_𝑟𝑓_𝑟−𝑓^′_𝑟𝜃_𝑟,1𝑆_𝑐𝜙^′′_𝑟+1+𝑓_𝑟𝜂−𝐴2𝜙^′_𝑟+1−𝑓^′_𝑟+𝐴+𝛾𝜙_𝑟+1−𝜙_𝑟𝑓^′_𝑟+1+𝜙^′_𝑟𝑓_𝑟+1=𝜙^′_𝑟𝑓_𝑟−𝑓^′_𝑟𝜙_𝑟 (61) subject to 𝑓_𝑟+1(0)=𝑓_𝑤,𝑓^′_𝑟+1(0)𝜃=1,_𝑟+1(0)=1,𝜙_𝑟+1(0)𝑓=1,_′𝑟+1(∞)=0,𝜃_𝑟+1(∞)=0,𝜙_𝑟+1(∞)=0. (62)

Applying the spectral method and compact finite difference schemes, we obtain the SQLM and CFD-QLM iteration schemes, respectively. In matrix form, the SQLM and CFD-QLM schemes are given by Δ_1,1Δ_1,2Δ_1,3Δ_2,1Δ_2,2Δ_2,3Δ_3,1Δ_3,2Δ_3,3𝑓_𝑟+1𝜃_𝑟+1𝜙_𝑟+1=Φ_1,𝑟Φ_2,𝑟Φ_3,𝑟. (63)

For the SQLM, we have Δ_1,1=𝐃³−diag𝑓_𝑟𝜂−𝐴2𝐃²−diag2𝑓^′_𝑟𝐃𝑓+𝐾+𝐴+diag_𝑟^′′,Δ_1,2=𝐺_𝑟Δ𝐈,_1,3=𝐺_𝑐Δ𝐈,_2,1=−diag𝜃_𝑟𝐃+diag𝜃^′_𝑟,Δ_2,2=1𝑃_𝑟𝐃²+diag𝑓_𝑟𝜂−𝐴2𝐃−𝐴𝐈−diag𝑓^′_𝑟,Δ_2,3Δ=𝐎,_3,1=−diag𝜙_𝑟𝐃+diag𝜙^′_𝑟,Δ_3,2Δ=𝐎,_3,3=1𝑆_𝑐𝐃²+diag𝑓_𝑟𝜂−𝐴2𝐃−𝐴+𝛾𝐈−diag𝑓^′_𝑟,Φ_1,𝑟=𝑓_𝑟𝑓_𝑟^′′−𝑓_𝑟^′2,Φ_2,𝑟=𝑓_𝑟𝜃^′_𝑟−𝑓^′_𝑟𝜃_𝑟,Φ_3,𝑟=𝑓_𝑟𝜙_𝑟−𝑓^′_𝑟𝜙_𝑟, (64) where O is an (𝑁+1)×(𝑁+1) zero matrix, and for the CFD-QLM, we have Δ_1,1=𝐸₃−diag𝑓_𝑟𝜂−𝐴2𝐸₂−diag2𝑓^′_𝑟+𝐾+𝐴𝐸₁+diag𝑓_𝑟^′′,Δ_1,2=𝐺_𝑟Δ𝐼,_1,3=𝐺_𝑐Δ𝐼,_2,1=−diag𝜃_𝑟𝐸₁+diag𝜃^′_𝑟,Δ_2,2=1𝑃_𝑟𝐸₂+diag𝑓_𝑟𝜂−𝐴2𝐸₁−𝐴𝐼−diag𝑓^′_𝑟,Δ_2,3Δ=𝑂,_3,1=−diag𝜙_𝑟𝐸₁+diag𝜙^′_𝑟,Δ_3,2Δ=𝑂,_3,3=1𝑆_𝑐𝐸₂+diag𝑓_𝑟𝜂−𝐴2𝐸₁−𝐴+𝛾𝐼−diag𝑓^′_𝑟,Φ_1,𝑟=𝑓_𝑟𝑓_𝑟^′′−𝑓_𝑟^′2−𝐴₃⁻¹𝐾₃−𝐴₃⁻¹𝐾^′_𝑟−diag𝑓_𝑟𝜂−𝐴2𝐴₂⁻¹𝐾₂−diag𝑓_𝑟𝜂−𝐴2𝐴₂⁻¹𝐾^′₂+diag2𝑓^′_𝑟𝐴₁⁻¹𝐾₁−(𝐴+𝐾)𝐴₁⁻¹𝐾₁+diag2𝑓^′_𝑟𝐴₁⁻¹𝐾^′₁−(𝐴+𝐾)𝐴₁⁻¹𝐾^′₁,Φ_2,𝑟=𝑓_𝑟𝜃^′_𝑟−𝑓^′_𝑟𝜃_𝑟−1𝑃_𝑟𝐴₂⁻¹𝐾₂−diag𝑓_𝑟𝜂−𝐴2𝐴₁⁻¹𝐾₁+diag𝜃_𝑟𝐴₁⁻¹𝐾₁+diag𝜃_𝑟𝐴₁⁻¹𝐾^′₁,Φ_3,𝑟=𝑓_𝑟𝜙_𝑟−𝑓^′_𝑟𝜙_𝑟−1𝑆_𝑐𝐴₂⁻¹𝐾₂−diag𝑓_𝑟𝜂−𝐴2𝐴₁⁻¹𝐾₁+diag𝜙_𝑟𝐴₁⁻¹𝐾₁+diag𝜙_𝑟𝐴₁⁻¹𝐾^′₁, (65) where 𝑂 is an (𝑁−1)×(𝑁−1) zero matrix and 𝐸₃=𝐴₃⁻¹𝐵₃.

4. Results and Discussion

In this section we present results for Examples 1 and 2 and give a comparison between the compact-finite-difference-based methods and the spectral-method-based methods, that is, the SRM and CFD-RM and the SQLM and CFD-QLM. The solution 𝑓(𝜂) and the velocity profiles 𝑓^′(𝜂) for varying values of 𝑀 for Example 1 are shown by Figure 1. The numerical results are compared against the exact analytical solution and good agreement is observed in all cases of varying parameter 𝑀. Tables 1 and 2 show the computed skin friction 𝑓^′′(0) to an accuracy of 10⁻⁹ obtained using the four methods.

Table 1 

Comparison of the SRM and CFD-RM results for the solution of 𝑓^′′(0) in Example 1.

Table 2 

Comparison of the SQLM and CFD-QLM results for the solution of 𝑓^′′(0) in Example 1.

(a)

(b)

Figure 1 

𝑓(𝜂) and the velocity profile 𝑓′(𝜂) for 𝑀=2,5,10.

It can be seen from Table 1 that in terms of computational speed, the SRM is efficient compared to the CFD-RM. This is because quite a large number of grid points are needed for the CFD-RM to give an accuracy of 10⁻⁹ as compared to the SRM. From Table 2 we observe similar results with the SQLM being computationally faster than the CFD-QLM. Again the number of grid points is the main reason for the difference in speed. The spectral method based approaches only require few grid points to give highly accurate results.

Figure 2 shows the errors at different iterations for the CFD-RM, and CFD-QLM respectively. In all cases the error decreases with each iteration which shows convergence of the methods. Increasing the value of the magnetic interaction parameter, 𝑀 increases the convergence of both the CFD-RM and CFD-QLM. This may be explained by the observation that when 𝑀 is very large, the dominant terms in the equation from which the SRM scheme is derived are 𝑔^′′−𝑀²𝑔. Solving this equation gives a solution of the form 𝑔=𝑐𝑒^−𝑀𝜂 (where 𝑐 is a constant). For large values of 𝑀, the 𝑀²−1 appearing in the given exact solution (38) is approximately equal to 𝑀². Consequently, the exact solution reduces to an exponential equation that has the same form as the approximate equation obtained from the SRM scheme at large 𝑀. The CFD-QLM shows a faster convergence than the CFD-RM as shown in Figure 2.

(a)

(b)

Figure 2 

Effect of varying 𝑀 on the error for the CFD-RM and CFD-QLM.

Figures 3 and 4 show the velocity 𝑓^′(𝜂), temperature 𝜃(𝜂), and concentration 𝜙(𝜂) for Example 2 for varying values of 𝐺_𝑟 and 𝑓_𝑤, respectively. We observe that the graphs are qualitatively similar to those reported in [23, 31]. The comparison of the performance of the methods is given by Tables 3 and 4. In Example 2 we compute the solution to an accuracy of 10⁻⁸. The SRM is computationally faster than the CFD-RM since it requires fewer grid points to give the required accuracy. We observe different results for the SQLM in Example 2 compared to Example 1. When comparing the results with the SRM, CFD-RM, and CFD-QLM, we observe that the SQLM fails to give results to the accuracy of 10⁻⁸ as it can be seen in Table 4. In this case the CFD-QLM is much better than the SQLM in terms of accuracy. This is one instance where the advantage of using CFD over spectral approach to integrate linearised equations is highlighted. The spectral approach is less accurate than the CFD approach when large systems of equations are considered.

Table 3 

Comparison of the SRM and CFD-RM results for the solution of 𝑓^′′(0) in Example 2.

Table 4 

Comparison of the SQLM and CFD-QLM results for the solution of 𝑓^′′(0) in Example 2.

(a)

(b)

(c)

Figure 3 

Velocity 𝑓(𝜂), temperature 𝜃(𝜂), and concentration 𝜙(𝜂) profiles in Example 2 for various values of Gr(Gr=0,5,20) when Pr=0.72,Sc=0.6,𝐴=0.5,𝛾=1,Gc=1,𝑓𝑤=1, and 𝐾=1.

(a)

(b)

Figure 4 

Velocity 𝑓(𝜂), temperature 𝜃(𝜂), and concentration 𝜙(𝜂) profiles in Example 2 for various values of 𝑓𝑤(𝑓𝑤=0,2,4) when Pr=0.72,Sc=0.6,𝐴=0.5,𝛾=1,Gc=1,Gr=1, and 𝐾=1.

In terms of convergence between the CFD-RM and CFD-QLM, we observe a faster convergence with the CFD-QLM than the CFD-RM. This is depicted in Figures 5 and 6.

(a)

(b)

Figure 5 

Effect of varying Gr on the errors for the CFD-RM and the CFD-QLM.

(a)

(b)

Figure 6 

Effect of varying 𝑓𝑤 on the errors for the CFD-RM and the CFD-QLM.

5. Conclusion

In this work we have introduced two new methods for solving systems of nonlinear boundary value equations. These methods called the compact finite difference relaxation method (CFD-RM) and compact finite difference quasilinearization method (CFD-QLM) are a modification of the spectral relaxation method (SRM) and the spectral quasilinearization method (SQLM), respectively. We compared the CFD schemes and spectral methods in solving nonlinear similarity boundary layer problems by comparing the CFD-based CFD-RM and CFD-QLM against the spectral-method-based SRM and SQLM. We did the comparison on one dimensional and three-dimensional problems. After comparing these methods we can conclude that they are all highly accurate with the spectral method outperforming the CFD in terms of computational speed. The SQLM though was less accurate for the three-dimensional problems. The CFD-QLM was able to handle the three-dimensional problem. We also observed that the CFD-QLM converges faster than the CFD-RM. When solving nonlinear boundary value problems, we recommend that the spectral-method-based algorithms must be used in problems involving small systems of equations and when computational speed is of importance because they are easy to code and computationally faster. CFD-based algorithms are ideal for solving large systems of equations when high accuracy is required.