Higher Order Compact Finite Difference Schemes for Unsteady Boundary Layer Flow Problems

Abstract

We investigate the applicability of the compact finite difference relaxation method (CFDRM) in solving unsteady boundary layer flow problems modelled by nonlinear partial differential equations. The CFDRM utilizes the Gauss-Seidel approach of decoupling algebraic equations to linearize the governing equations and solve the resulting system of ordinary differential equations using compact finite difference schemes. The CFDRM has only been used to solve ordinary differential equations modelling boundary layer problems. This work extends its applications to nonlinear partial differential equations modelling unsteady boundary layer flows. The CFDRM is validated on two examples and the results are compared to results of the Keller-box method.

1. Introduction

The many, important applications associated with boundary layer flow and heat transfer induced by a stretching surface have made them one of the most studied problems in the field of fluid dynamics. These applications include aerodynamic extrusion of plastic sheets, the cooling process of metallic plates in a cooling bath, and the cooling and/or drying of paper and textiles Sakiadis [1] was the first to study boundary layer flows induced by a stretching surface. Since then, due to their industrial importance they have been studied by quite a number of researchers. Initially, steady flows were considered, but since it is not always to maintain steady state conditions, in recent years the focus has shifted towards unsteady flows. In this work, we consider unsteady flows modelled by partial differential equations.

In solving the unsteady boundary layer problems many researchers have employed the Williams III and Rhyne [2] transformation. Williams III and Rhyne introduced self-similar coordinates to transform the equations governing the unsteady flows by converting the infinite time scale of unsteady problems to a finite region of integration. Various analytical and numerical methods have been used to solve the transformed equations. The Keller-box method of Cebeci and Bradshaw [3] has been very popular in solving this kind of problems amongst researchers. These include Seshadri et al. [4] who combined the Keller-box and perturbation series approach for the solution of unsteady mixed convection flow along a heated vertical plate. Nazar et al. [5, 6] solved the unsteady boundary layer flow problem due to an impulsively stretching surface in a rotating fluid by means of the Keller-box numerical method, and they obtained a first-order perturbation approximation of the solution. The limitation of the perturbation approach, as noted by Liao [7], is that it gives solutions that are only valid for small time. Liao proposed the use of the homotopy analysis method (HAM) instead and since then the HAM has been dominantly used in solving the unsteady boundary layer flow problems (see [8–16]). The advantage of the HAM over the perturbation techniques is that it is able to produce solutions valid for all time.

In this work we propose the use of a new method that combines the use of higher order compact finite difference (CFD) schemes and the Gauss-Seidel approach of solving algebraic equations. The method called compact finite difference relaxation method (CFDRM) has only been used once by Dlamini et al. [17] to solve ordinary differential equations. The CFDRM is based on simple decoupling and rearrangement of the governing equations and numerically integrating the resulting equations using the CFD schemes. The use of higher order CFD schemes is motivated by the fact that they produce highly accurate solutions on coarser grids with greater computational efficiency. Various CFD schemes used for applications such as interpolation, filtering, and evaluating high-order derivatives were discussed in detail by Lele [18]. CFD schemes have been used to solve, for example, Burger's equation [19, 20], Navier-Stokes equation [21], Korteweg-de Vries equation [22], Black-Scholes equation [23], and many more [24–26]. In this work we explore their usage in solving unsteady boundary layer flow problems.

We consider two examples, the unsteady one-dimensional MHD laminar boundary layer flow due to an impulsively stretching surface that was previously studied by Srinivasa and Eswara [27], and the unsteady three-dimensional MHD flow and heat transfer over an impulsively stretching plate previously studied by Xu et al. [8]. The main objective of the study is to investigate the applicability of the CFDRM to solve unsteady boundary layer problems modelled by partial differential equations. The results are compared with results of the Keller-box method.

The rest of the paper is organized as follows. In Section 2, we discuss the derivation of the compact finite difference schemes, In Section 3 we discuss the development of the CFDRM for computing the solution of an unsteady MHD laminar boundary layer flow due to an impulsively stretching surface. Section 4 presents CFDRM implementation of an unsteady three-dimensional MHD flow and heat transfer over an impulsively stretching plate. Section 5 contains the results and discussion and the conclusions are given in Section 6.

2. Compact Finite Difference Schemes

In the derivation of the CFD schemes we consider a one-dimensional uniform mesh on the region [𝑎,𝑏] with nodes 𝑥_𝑖(𝑖=1,2,…,𝑁), where 𝑎=𝑥₁<𝑥₂<⋯<𝑥_𝑁=𝑏, (1) 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 [18] 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ℎ². (2)(3)

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,=𝑦_𝑎,𝑦(𝑏)=𝑦_𝑏, (4) where 𝑓(𝑥,𝑦,𝑦^′) is a nonlinear function, 𝑝(𝑥) and 𝑞(𝑥) are known functions of 𝑥, and 𝑦_𝑎 and 𝑦_𝑏 are known constants. In solving (4) we apply the CFD approximation for the first and second derivatives given by (2) and (3), 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ℎ𝑎₁𝑦₁+𝑎₂𝑦₂+𝑎₃𝑦₃+𝑎₄𝑦₄+𝑎₅𝑦₅+𝑎₆𝑦₆+𝑎₇𝑦₇, (5) and when 𝑖=𝑁−1 we use 13𝑦^′_𝑁−2+𝑦^′_𝑁−1=1ℎ𝑏₁𝑦_𝑁+𝑏₂𝑦_𝑁−1+𝑏₃𝑦_𝑁−3+𝑏₄𝑦_𝑁−4+𝑏₅𝑦_𝑁−5+𝑏₆𝑦_𝑁−6+𝑏₇𝑦_𝑁−7, (6) 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 (5) and (6), 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,𝑎₃=8336,𝑎₄=−119,𝑎₅=23,𝑎₆=−37180,𝑎₇=1,𝑏36_𝑘=−𝑎_𝑘,𝑘=1,2,…,7. (7)

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

Using the above equations, the equations for approximating the first- and second-order derivatives can be expressed as 𝐴₁𝑌^′=𝐵₁𝑌+𝐾₁,𝐴₂𝑌^′′=𝐵₂𝑌+𝐾₂, (11) 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^𝑇. (12)

From (11), 𝑌^′ and 𝑌^′′ can be expressed as 𝑌^′=𝐸₁𝑌+𝐻₁,𝑌^′′=𝐸₂𝑌+𝐻₂, (13) where 𝐸₁=𝐴₁⁻¹𝐵₁,𝐸₂=𝐴₂⁻¹𝐵₂,𝐻₁=𝐴₁⁻¹𝐾₁,𝐻₂=𝐴₂⁻¹𝐾₂. (14)

3. Unsteady MHD Laminar Boundary Layer Flow Due to an Impulsively Stretching Surface

We consider the unsteady, laminar incompressible flow of a vicious electrically conducting fluid over a linear surface. This problem was studied by Srinivasa and Eswara [27]. They used the Keller-box method to solve the problem. After applying Williams III and Rhyne [2] similarity transformation to the original equation governing this flow they obtained the following dimensionless governing equation for the flow 𝑓^′′′+12𝜂(1−𝜉)𝑓^′′+𝜉𝑓𝑓^′′−𝑀𝑓^′−𝑓^′2=𝜉(1−𝜉)𝜕𝑓^′,𝜕𝜉 (15) subject to the boundary conditions 𝑓(0,𝜉)=0,𝑓^′(0,𝜉)=1,𝑓^′(∞,𝜉)=0, (16) where 𝜉∈[0,1] is the dimensionless time scale defined as 𝜉=1−𝑒^−𝜏,𝜏=𝑏𝑡, (17) with 𝑏 as positive constant and 𝑡 as the time variable. In the analysis of boundary layer flow problems, a quantity that is of physical interest is the skin friction which in this model is given [4, 5, 7], in dimensionless form, as 𝐶_𝑓Re_𝑥^1/2=𝜉^−1/2𝑓^′′(𝜉,0), (18) where Re_𝑥 is the local Reynolds number.

The initial unsteady solution at 𝜉=0 (𝜏=0) for the governing equation (15) is obtained as a solution of the equation 𝑓^′′′+12𝜂𝑓^′′()𝑓=0,𝑓0,0=0,_′(0,0)=1,𝑓^′(∞,0)=0, (19) where the primes denote differentiation with respect to 𝜂. Solving (19) gives 𝑓𝜂𝜂,0=𝜂erfc2+2√𝜋−𝜂1−exp²4, (20) where erfc(𝜂) is the standard complementary error function defined by 𝜂=2erfc√𝜋^∞_𝜂exp−𝑧²𝑑𝑧. (21)

The steady state solution when 𝜉=1, corresponding to 𝜏→+∞, is obtained from 𝑓^′′′+𝑓𝑓^′′−𝑓^′2()𝑓=0,𝑓0,1=0,_′(0,1)=1,𝑓^′(∞,1)=0. (22) The solution to the above equation is 𝑓𝜂,1=1−exp𝜂. (23)

Now we discuss the development of the CFDRM to solve the partial differential equation (15). The CFDRM works for systems of differential equations and hence it is convenient to reduce the order of (15) from three to two. We set 𝑓^′=𝑢, so that (15) becomes 𝑢^′′+12𝜂(1−𝜉)𝑢^′+𝜉𝑓𝑢^′−𝑀𝑢−𝑢²=𝜉(1−𝜉)𝜕𝑢,𝑓𝜕𝜉_′=𝑢. (24)

The CFDRM algorithm uses the same ideas as the SRM described in [28–30]. Both employ the Gauss Seidel idea of decoupling systems of algebraic equations to decouple differential equations. In the framework of the SRM the resulting linear system of differential equations is integrated by the spectral collocation method. Instead of using spectral method, the CFDRM uses the compact finite differences to integrate the linearized system of differential equations. Applying the CFDRM on (24) we obtain the following linear PDEs: 𝑢^′′_𝑟+1+𝑎_1,𝑟𝑢^′_𝑟+1+𝑎_2,𝑟𝑢_𝑟+1+𝑎_3,𝑟=𝜉(1−𝜉)𝜕𝑢_𝑟+1,𝑓𝜕𝜉_′𝑟+1=𝑢_𝑟+1,𝑓_𝑟+1(0,𝜉)=0,𝑢_𝑟+1(0,𝜉)=1,𝑢_𝑟+1(∞,𝜉)=0, (25)(26)(27) where 𝑎_1,𝑟=12𝜂(1−𝜉)+𝜉𝑓_𝑟,𝑎_2,𝑟=−𝜉𝑢²_𝑟. (28)

Here, 𝑟 denotes the iteration level. All linear terms are evaluated at current iteration level (𝑟+1) and nonlinear terms at the previous iteration level (𝑟). The initial approximations for solving (25)-(26) are obtained as the solutions at 𝜉=0. Thus 𝑓₀(𝜂,𝜉) and 𝑢₀(𝜂,𝜉) are given by 𝑓₀𝜂𝜂,𝜉=𝜂erfc2+2√𝜋−𝜂1−exp²4,𝑢₀𝜂𝜂,𝜉=erfc2. (29) Starting from given initial approximations (29), the iteration scheme (25) can be solved iteratively for 𝑢_𝑟+1(𝜂,𝜉) when 𝑟=0,1,2,…. The solution for 𝑢_𝑟+1 is used in (26) which is, in turn, solved for 𝑓_𝑟+1. To solve (25) we use sixth order compact finite difference schemes described in the previous section and use an implicit finite difference method in the 𝜉-direction.

The finite difference scheme with centering about a midpoint halfway between 𝜉^𝑛+1 and 𝜉^𝑛 is applied in the 𝜂 direction. This midpoint is defined as 𝜉^𝑛+(1/2)=(𝜉^𝑛+1+𝜉^𝑛)/2. Applying the centering about 𝜉^𝑛+(1/2) to any function, say 𝑢(𝜂,𝜉) and its associated derivative, we obtain 𝑢𝜂_𝑗,𝜉^𝑛+(1/2)=𝑢_𝑗^𝑛+(1/2)=𝑢_𝑗^𝑛+1+𝑢^𝑛_𝑗2,𝜕𝑢𝜕𝜉^𝑛+(1/2)=𝑢_𝑗^𝑛+1−𝑢^𝑛_𝑗.Δ𝜉 (30)

We first apply the CFD schemes on (25) and (26) to obtain 𝐄₂+𝐚_1,𝑟𝐄₁𝑈_𝑟+1+𝐇₂+𝐚_1,𝑟𝐇₁+𝐚_2,𝑟=𝜉(1−𝜉)𝑑𝑈_𝑟+1,𝐄𝑑𝜉₁𝐹_𝑟+1+𝐇₁=𝑈_𝑟+1,𝑓_𝑟+1𝜂_𝑁𝑥,𝜉=0,𝑢_𝑟+1𝜂₀𝑢,𝜉=0,_𝑟+1𝜂_𝑁𝑥,𝜉=1, (31)(32)(33) with initial approximations 𝑓_𝑟+1𝜂_𝑗,0=𝜂_𝑗𝜂erfc^𝑗2+2√𝜋−𝜂1−exp^2𝑗4,𝑢_𝑟+1𝜂_𝑗𝜂,0=erfc^𝑗2,𝑗=0,1,2,…,𝑁_𝑥,𝑈_𝑟+1=𝑢_𝑟+1𝜂₀𝑢,𝜉_𝑟+1𝜂₁⋮𝑢,𝜉_𝑟+1𝜂_𝑁𝑥−1𝑢,𝜉_𝑟+1𝜂_𝑁𝑥,𝜉,𝐹_𝑟+1=𝑓_𝑟+1𝜂₀𝑓,𝜉_𝑟+1𝜂₁⋮𝑓,𝜉_𝑟+1𝜂_𝑁𝑥−1𝑓,𝜉_𝑟+1𝜂_𝑁𝑥,𝐚,𝜉_2,𝑟=𝑎_2,𝑟𝜂₀𝑎,𝜉_2,𝑟𝜂₁⋮𝑎,𝜉_2,𝑟𝜂_𝑁𝑥−1𝑎,𝜉_2,𝑟𝜂_𝑁𝑥,𝐚,𝜉_1,𝑟=𝑎_1,𝑟𝜂₀𝑎,𝜉_1,𝑟𝜂₂⋱⋱𝑎,𝜉_1,𝑟𝜂_𝑁𝑥.,𝜉 (34)

Next, we apply the finite difference scheme on (31) in the 𝜉-direction with centering about the midpoint 𝜉^𝑛+(1/2) to obtain 𝐗𝑈^𝑛+1_𝑟+1=𝐘𝑈^𝑛_𝑟+1+𝐆, (35) subject to the following boundary and initial conditions: 𝑢_𝑟+1𝜂₀,𝜉^𝑛=0,𝑢_𝑟+1𝜂_𝑁𝑥,𝜉^𝑛=1,𝑛=0,1,2,…,𝑁_𝑡,𝑢_𝑟+1𝜂_𝑗𝜂,0=erfc^𝑗2,𝑗=0,1,2,…,𝑁_𝑥, (36) where 1𝐗=−2𝐄₂+𝐚^𝑛+(1/2)_1,𝑟𝐄₁+𝜉^𝑛+(1/2)1−𝜉^𝑛+(1/2)1Δ𝜉𝐈,𝐘=2𝐄₂+𝐚^𝑛+(1/2)_1,𝑟𝐄₁+𝜉^𝑛+(1/2)1−𝜉^𝑛+(1/2)Δ𝜉𝐈,𝐆=−𝐇₂−𝐚_1,𝑟𝐇₁−𝐚^𝑛+(1/2)_2,𝑟, (37) where 𝐈 is an (𝑁_𝑥+1)×(𝑁_𝑥+1) identity matrix. The matrices 𝐄₁,𝐄₂,𝐇₁, and 𝐇₂ are as described in Section 3 above. Thus, starting from the initial condition 𝑈⁰_𝑟+1, given by (29), (35) can be solved iteratively to give approximate solutions for 𝑢_𝑟+1(𝜂,𝜉), 𝑟=0,1,2,…, until a solution that converges to within a given accuracy level is obtained. The solution 𝑢_𝑟+1 is used in (32) which is, in turn, solved for 𝑓_𝑟+1.

4. Unsteady Three-Dimensional MHD Flow and Heat Transfer over an Impulsively Stretching Plate

We consider the unsteady, three-dimensional flow of an electrically conducting fluid caused by an impulsive stretching surface in two lateral directions. The fluid is electrically conducting in a traverse magnetic field. The magnetic field is neglected under the assumption of a small magnetic Reynolds number. Xu et al. [8] solved this problem using the homotopy analysis method. Here we use the compact finite difference relaxation method to solve the governing equations. The flow is governed by the following three partial differential equations: 𝑓^′′′+()𝜂1−𝜉2𝑓^′′−𝜉𝜕𝑓^′𝑓𝜕𝜉+𝜉𝑓+𝑠^′′−𝑓^′2−𝑀𝑓^′𝑠=0,_′′′+()𝜂1−𝜉2𝑔^′′−𝜉𝜕𝑠^′𝑠𝜕𝜉+𝜉𝑓+𝑠^′′−𝑠^′2−𝑀𝑠^′𝑔=0,_′′()𝜂+Pr1−𝜉2𝑔^′−𝜉𝜕𝑔𝑔𝜕𝜉+Pr𝜉𝑓+𝑠^′=0, (38)(39)(40) with the following boundary conditions: 𝑓(𝜉,0)=𝑠(𝜉,0)=0,𝑓^′()()𝑠𝜉,0=𝑔𝜉,0=1,_′()𝑓𝜉,0=𝑐,_′(𝜉,∞)=𝑠^′(𝜉,∞)=𝑔(𝜉,∞)=0. (41)

In the above equations prime denotes the derivative with respect to 𝜂 and 𝑐 the stretching parameter is a positive constant. 𝑀 is the local Hartman number and Pr the Prandtl number. The initial unsteady solution can be found exactly by setting 𝜉=0 in the above equations and solving the resulting equations. The closed form analytical solutions are given by 𝑓𝜂0,𝜂=𝜂erfc2+2√𝜋−𝜂1−exp²4,𝑠𝜂0,𝜂=𝑐𝜂erfc2+2√𝜋−𝜂1−exp²4,𝑔√0,𝜂=erfcPr𝜂2. (42)

Now we discuss the development of the spectral relaxation method to solve the system of partial differential equations (38)–(40). First, we set 𝑓^′=𝑢 and 𝑠^′=𝑣, so that (38) and (39) become 𝑢^′′+()𝜂1−𝜉2𝑢^′−𝜉𝜕𝑢𝑢𝜕𝜉+𝜉𝑓+𝑔^′−𝑢²𝑣−𝑀𝑢=0,_′′+()𝜂1−𝜉2𝑣^′−𝜉𝜕𝑣𝑣𝜕𝜉+𝜉𝑓+𝑔^′−𝑣²−𝑀𝑣=0. (43) Applying the CFDRM on the resulting system of nonlinear partial differential equations gives the following linear partial differential equations: 𝑢^′′_𝑟+1+𝑎_1,𝑟𝑢^′_𝑟+1+𝑎_2,𝑟𝑢_𝑟+1+𝑎_3,𝑟=𝜉(1−𝜉)𝜕𝑢_𝑟+1,𝑓𝜕𝜉_′𝑟+1=𝑢_𝑟+1,𝑣^′′_𝑟+1+𝑏_1,𝑟𝑣^′_𝑟+1+𝑏_2,𝑟𝑣_𝑟+1+𝑏_3,𝑟=𝜉(1−𝜉)𝜕𝑣_𝑟+1,𝑠𝜕𝜉_′𝑟+1=𝑣_𝑟+1,𝑔^′′_𝑟+1+𝑐_1,𝑟𝑔^′_𝑟+1=𝜉(1−𝜉)𝜕𝜃_𝑟+1,𝜕𝜉 (44) with the following boundary conditions: 𝑢_𝑟+1(0,𝜉)=𝑔_𝑟+1()𝑓0,𝜉=1,_𝑟+1(0,𝜉)=𝑠_𝑟+1()𝑣0,𝜉=0,_𝑟+1()𝑢0,𝜉=𝑐,_𝑟+1(∞,𝜉)=𝑣_𝑟+1(∞,𝜉)=𝑔_𝑟+1(∞,𝜉)=0, (45) where 𝑎_1,𝑟=12𝜂(1−𝜉)+𝜉𝑓_𝑟+𝑠_𝑟,𝑎_2,𝑟=−𝑀𝜉,𝑎_3,𝑟=−𝜉𝑢²_𝑟,𝑏_1,𝑟=𝑎_1,𝑟,𝑏_2,𝑟=𝑎_2,𝑟,𝑏_3,𝑟=−𝜉𝑣²_𝑟,𝑐_1,𝑟1=Pr2𝜂(1−𝜉)+𝜉𝑓_𝑟+𝑠_𝑟. (46)

Starting from given initial approximations, denoted by 𝑢₀(𝜂,𝜉), 𝑓₀(𝜂,𝜉), 𝑣₀(𝜂,𝜉), 𝑠₀(𝜂,𝜉), and 𝑔₀(𝜂,𝜉), (44) can be solved iteratively for the unknown functions. To solve the above decoupled system of differential equations we apply CFD schemes on the space variable and finite differences in the time variable as described previously and obtain the following system of decoupled equations: 𝐗₁𝑈^𝑛+1_𝑟+1=𝐘₁𝑈^𝑛_𝑟+1+𝐺₁,𝐄₁𝐹^𝑛+1_𝑟+1+𝐇₁=𝑈^𝑛+1_𝑟+1,𝐗₂𝑉^𝑛+1_𝑟+1=𝐘₂𝑉^𝑛_𝑟+1+𝐺₂,𝐄₁𝑆^𝑛+1_𝑟+1+𝐇₁=𝑉^𝑛+1_𝑟+1,𝐗₃𝐺^𝑛+1_𝑟+1=𝐘₃𝐺^𝑛_𝑟+1+𝐺₃, (47) with the following boundary conditions: 𝑢_𝑟+1𝜂₀,𝜉^𝑛=𝑔_𝑟+1𝜂₀,𝜉^𝑛𝑣=1,_𝑟+1𝜂₀,𝜉^𝑛=𝑐,𝑓_𝑟+1𝜂₀,𝜉^𝑛=𝑠_𝑟+1𝜂₀,𝜉^𝑛𝑢=0,_𝑟+1𝜂_∞,𝜉^𝑛=𝑣_𝑟+1𝜂_∞,𝜉^𝑛=𝑔_𝑟+1𝜂_∞,𝜉^𝑛=0,𝑛=0,1,2,…,𝑁_𝑡, (48) and initial functions 𝑓_𝑟+1𝜂_𝑗𝜂,0=𝜂erfc^𝑗2+2√𝜋−𝜂1−exp^2𝑗4,𝑠_𝑟+1𝜂_𝑗𝜂,0=𝑐𝜂erfc^𝑗2+2√𝜋−𝜂1−exp^2𝑗4,𝑔_𝑟+1𝜂_𝑗√,0=erfcPr𝜂_𝑗2,𝑗=0,1,2,…,𝑁_𝑥. (49) The matrices above are defined as 𝐗₁=12𝐄₂+𝐚^𝑛+(1/2)_1,𝑟𝐄₁+𝐚_2,𝑟−𝜉^𝑛+(1/2)1−𝜉^𝑛+(1/2)𝐗Δ𝜉𝐈,₂=12𝐄₂+𝐛^𝑛+(1/2)_1,𝑟𝐄₁+𝐛_2,𝑟−𝜉^𝑛+(1/2)1−𝜉^𝑛+(1/2)𝐗Δ𝜉𝐈,₃=12𝐄₂+𝐜^𝑛+(1/2)_1,𝑟𝐄₁−Pr𝜉^𝑛+(1/2)1−𝜉^𝑛+(1/2)𝐘Δ𝜉𝐈,₁1=−2𝐄₂+𝐚^𝑛+(1/2)_1,𝑟𝐄₁+𝐚_2,𝑟−𝜉^𝑛+(1/2)1−𝜉^𝑛+(1/2)𝐘Δ𝜉𝐈,₂1=−2𝐄₂+𝐛^𝑛+(1/2)_1,𝑟𝐄₁+𝐛_2,𝑟−𝜉^𝑛+(1/2)1−𝜉^𝑛+(1/2)𝐘Δ𝜉𝐈,₃1=−2𝐄₂+𝐜^𝑛+(1/2)_1,𝑟𝐄₁−Pr𝜉^𝑛+(1/2)1−𝜉^𝑛+(1/2)𝐆Δ𝜉𝐈,₁=𝐚^𝑛+(1/2)_3,𝑟−𝐇₂−𝐚^𝑛+(1/2)_1,𝑟𝐇₁,𝐆₂=𝐛^𝑛+(1/2)_3,𝑟−𝐇₂−𝐛^𝑛+(1/2)_1,𝑟𝐇₁,𝐆₃=−𝐇₂−𝐜^𝑛+(1/2)_1,𝑟𝐇₁. (50)

The vectors 𝑈, 𝐹, 𝑉, 𝑆, and 𝐺 are the vectors of the functions 𝑢, 𝑓, 𝑣, 𝑠, and 𝑔, respectively, when evaluated at the grid points 𝜂_𝑗  (𝑗=0,1,…,𝑁_𝑥).

5. Results and Discussion

In this section we give the CFDRM and Keller-box method results for the two examples described above. The Keller-box method is an implicit finite difference scheme which is made up of four major steps. First, the governing differential equation is reduced to a system of first-order equations. The resulting first-order equations are discretized using central finite differences. The resulting algebraic equations are linearised by Newton's method and written in matrix vector form. The linear system obtained is then solved by block tridiagonal elimination technique.

Example 1 was solved for 𝑀=1. The value of 𝜂_∞, which correspond to 𝜂=∞, was approximated to be equal to 10. This value was found to be large enough to give consistent results. For the CFDRM, uniform stepsizes Δ𝜂=0.1 in the 𝜂-direction and Δ𝜉=0.0025 in the 𝜉-direction were used to get results accurate to six decimal places. For the Keller-box, stepsizes Δ𝜂=0.005 and Δ𝜉=0.004 were used. Table 1 shows the skin friction 𝑓^′′(0,𝜉) for Example 1 obtained using both CFDRM and Keller-box method. The times taken to compute the solution at different values of 𝜉 are also given in Table 1. It can be seen that the CFDRM takes a very short period of time to compute the solution compared to the Keller-box method. The Keller-box required quite a small Δ𝜂 to get a solution accurate to six decimal places. This means that more points are needed in the Keller-box compared to the CFDRM. Figure 1 shows the effects of varying the magnetic field parameter (𝑀) on the skin friction 𝑓^′′(0,𝜉). The parameters 𝑀=1, 𝑐=0.5, and Pr=0.7 were used to solve the governing equations in Example 2. The value of 𝜂_∞ was approximated to be equal to 10. For the CFDRM, uniform stepsizes Δ𝜂=0.1 and Δ𝜉=0.0025 to get results accurate to six decimal places. For the Keller-box, stepsizes Δ𝜂=0.005 and Δ𝜉=0.002 were used. The quantities of physical interest (𝑓^′′(0,𝜉), 𝑠^′′(0,𝜉), and 𝑔^′(0,𝜉)) are tabulated for both CFDRM and Keller-box and are shown in Table 2. The results are in good agreement with six decimal places as depicted in the table. Table 3 shows the times taken by each method at various values of 𝜉. It is evident from both Tables 2 and 3 that the CFDRM computes the solution in a very short period of time compared to the Keller-box. The main reason is that a very small Δ𝜂 and hence a large number of points are required in the Keller-box method for the required accuracy. Graphical solutions of 𝑓^′′(0,𝜉), 𝑠^′′(0,𝜉), and 𝑔^′(0,𝜉) in the entire range 0≤𝜉≤1 for Example 2 are shown in Figures 2, 3, and 4. The solution is plotted for different values of 𝑀 and 𝑐. Comparing the plots with those of [8], it can be seen that the results of the CFDRM are comparable with their homotopy analysis method results. Xu et al. [8] used the homotopy analysis method to get series solutions which are accurate and uniformly valid for all dimensionless time 0≤𝜉≤1. From these figures it is clear that the CFDRM is able to give results which are accurate and uniformly valid for all dimensionless time 0≤𝜉≤1.

Table 1 

CFDRM and Keller-box results for the skin friction 𝑓^′′(0,𝜉) in Example 1 when 𝑀=1.

Table 2 

CFDRM and Keller-box results for 𝑓^′′(0,𝜉),𝑠^′′(0,𝜉), and 𝑔^′(0,𝜉) in Example 2 when 𝑀=1, 𝑐=0.5, and Pr=0.7.

Table 3 

Computational times for Example 2 at various values of 𝜉.

Figure 1 

CFDRM solution of −𝑓′′(0,𝜉) for different values of 𝑀 for Example 1.

Figure 2 

CFDRM solution of −𝑓′′(0,𝜉) and −𝑠′′(0,𝜉) for different values of 𝑀 when 𝑐=0.5 and Pr=0.7.

Figure 3 

CFDRM solution of −𝑓′′(0,𝜉) and −𝑠′′(0,𝜉) for different values of 𝑐 when 𝑀=1 and Pr=0.7.

Figure 4 

CFDRM solution of −𝑔′(0,𝜉) for different values of 𝑀 and 𝑐, respectively, when Pr=0.7.

6. Conclusion

In this paper, two unsteady boundary layer flow problems have been investigated. The nonlinear partial differential equations governing the unsteady boundary layer flow of the two examples considered are solved using the compact finite difference relaxation method (CFDRM). The CFDRM had previously been applied only to ordinary differential equations modelling boundary layer problems [17]. This work has successfully extended its application to partial differential equations. The results obtained are in good agreement with results obtained using the Keller-box. The CFDRM is computationally faster than the Keller-box method. This is because the CFDRM gives highly accurate solutions on coarser grids.