A Novel Numerical Technique for Two-Dimensional Laminar Flow between Two Moving Porous Walls

Abstract

We investigate the steady two-dimensional flow of a viscous incompressible fluid in a rectangular domain that is bounded by two permeable surfaces. The governing fourth-order nonlinear differential equation is solved by applying the spectral-homotopy analysis method and a novel successive linearisation method. Semianalytical results are obtained and the convergence rate of the solution series was compared with numerical approximations and with earlier results where the homotopy analysis and homotopy perturbation methods were used. We show that both the spectral-homotopy analysis method and successive linearisation method are computationally efficient and accurate in finding solutions of nonlinear boundary value problems.

1. Introduction

Laminar viscous flow in tubes that allow seepage across contracting or expanding permeable walls is encountered in the transport of biological fluids such as blood and filtration in kidneys and lungs. Such flows have many other practical applications such as in binary gas diffusion, chromatography, ion exchange, and ground water movement [1–6]. In addition, flow in channels with permeable walls provides a good starting point for the study of flow in multichannel filtration systems such as the wall flow monolith filter used to reduce emissions from diesel engines introduced by Oxarango et al. in [7]. Consequently, in the past four decades a considerable amount of research effort has been expended in the study of laminar flows in rectangular domains that are bounded by permeable walls [8–15].

The equations governing such flows are generally nonlinear and in the past asymptotic techniques, and numerical methods have been used to analyze such flows and to solve the equations; for example, in the pioneering study by Berman [8] asymptotic methods were used to solve the steady flow problem for small suction. In the study by Uchida and Aoki [16], numerical methods were used to solve the governing nonlinear equations and to explain the flow characteristics. Majdalani and Roh [4] and Majdalani [3] studied the oscillatory channel flow with wall injection, and the resulting governing equations were solved using asymptotic formulations (WKB and multiple-scale techniques). The multiple-scale solution was found to be advantageous over the others in that its leading-order term is simpler and more accurate than other formulations, and it displayed clearly the relationship between the physical parameters that control the final motion. It also provided means of quantifying important flow features such as corresponding vortical wave amplitude, rotational depth of penetration, and near wall velocity overshoot to mention a few. Jankowski and Majdalani [12] used the same approach and drew similar conclusions about the multiple-scale solution for oscillatory channel flow with arbitrary suction.An analytical solution by means of the Liouville-Green transformation was developed for laminar flow in a porous channel with large wall suction and a weakly oscillatory pressure by Jankowski and Majdalani [13]. The scope of the problem had many limitations, for example, the study did not consider variations in thermostatic properties and the oscillatory pressure amplitude was taken to be small in comparison with the stagnation pressure. Zhou and Majdalani [17] investigated the mean-flow for slab rocket motors with regressing walls. The transformed governing equation was solved numerically, using finite differences, and asymptotically, using variation of parameters and small parameter perturbations in the blowing Reynolds number. The results from the two methods were compared for different Reynolds numbers Re and the wall regression rate 𝛼, and it was observed that accuracy of the analytical solution deteriorates for small Re and large 𝛼. A good agreement between the solutions was observed for large values of Re. A similar analysis was done by Majdalani and Zhou [6] for moderate-to-large injection and suction driven channel flows with expanding or contracting walls.

In recent years, the use of nonperturbation techniques such as the Adomian decomposition method [18, 19]. He's homotopy perturbation method [20, 21], and the homotopy analysis method [22, 23] has been increasingly preferred to solve nonlinear differential equations that arise in science and engineering. Dinarvand et al. [2] solved Berman’s model of two-dimensional viscous flow in porous channels with wall suction or injection using both the HAM and the homotopy perturbation method (HPM). They concluded that the HPM solution is not valid for large Reynolds numbers, a weakness earlier observed in the case of other perturbation techniques. Using the homotopy analysis method, Xu et al. [24] developed highly accurate series approximations for two-dimensional viscous flow between two moving porous walls and obtained multiple solutions associated with this problem. The multiple solutions associated with this problem were also reported by Zarturska et al. [25]. Although the homotopy analysis method is a reliable and efficient semi-analytical technique, it however suffers from a number of limiting assumptions such as the requirements that the solution ought to conform to the so-called rule of solution expression and the rule of coefficient ergodicity. A modification of the homotopy analysis method, see Motsa et al. [26, 27], seeks to produce a more efficient method while also addressing the limitations of the HAM. In this paper, we use the spectral homotopy analysis method to solve the nonlinear differential equation that governs the flow of a viscous incompressible fluid in a rectangular domain bounded by two permeable walls. The problem is also solved using a new and highly efficient technique, the successive linearisation method (see [28, 29]) so as to independently corroborate and validate the SHAM results. The results are also compared with numerical approximations and the recent results reported in Xu et al. [24] and Dinarvand and Rashidi [30].

2. Governing Equations

Consider two-dimensional laminar, isothermal, and incompressible viscous fluid flow in a rectangular domain bounded by two permeable surfaces that enable the fluid to enter or exit during successive expansions or contractions. The walls are placed at a separation 2𝑎 and contract or expand uniformly at a time-dependent rate 𝑎(̇𝑡). The governing Navier-Stokes equations are given in Majdalani et al. [31] as 𝜕𝑢𝜕𝑥+𝜕𝑣𝜕𝑦𝜕𝑢=0,+𝑢𝜕𝑢𝜕𝑡𝜕𝑥+𝑣𝜕𝑢𝜕𝑦1=−𝜌𝜕𝑝𝜕𝑥+𝜈∇²𝜕𝑣𝑢,+𝑢𝜕𝑣𝜕𝑡𝜕𝑥+𝑣𝜕𝑣𝜕𝑦1=−𝜌𝜕𝑝𝜕𝑦+𝜈∇²𝑣, (2.1)(2.2)(2.3) where 𝑢 and 𝑣 are the velocity components in the 𝑥 and 𝑦 directions, respectively, 𝑝, 𝜌, 𝜈 and 𝑡 are the dimensional pressure, density, kinematic viscosity, and time, respectively. Assuming that inflow or outflow velocity is 𝑣_𝑤, then the boundary conditions are𝑢(𝑥,𝑎)=0,𝑣(𝑎)=−𝑣_𝑤̇𝜕𝑢=−𝑎/𝑐,𝜕𝑦(𝑥,0)=0,𝑣(0)=0,𝑢(0,𝑦)=0, (2.4) where 𝑐(=𝑎/̇𝑣_𝑤) is the injection or suction coefficient. Introducing the stream function 𝜓=𝜈𝑥𝐹(𝑦,𝑡)/𝑎 and the transformations 𝜓𝜓=𝑎𝑎̇𝑢,𝑢=𝑎̇𝑣,𝑣=𝑎̇𝑥,𝑥=𝑎𝑦,𝑦=𝑎𝐹,𝐹=,Re (2.5) Majdalani et al. [31] and Dinarvand and Rashidi [30] showed that (2.1)–(2.3) reduce to the normalized nonlinear differential equation𝐹^𝐼𝑉+𝛼(𝑦𝐹^′′′+3𝐹^′′)+Re(𝐹𝐹^′′′−𝐹^′𝐹^′′)=0, (2.6) subject to the boundary conditions𝐹=0,𝐹^′′=0,at𝑦=0,𝐹=1,𝐹^′=0,at𝑦=1, (2.7)(2.8) where 𝛼(𝑡)=𝑎𝑎̇/𝜈 is the nondimensional wall dilation rate defined to be positive for expansion and negative for contraction, and Re=𝑎𝑣_𝑤/𝜈 is the filtration Reynolds number defined positive for injection and negative for suction through the walls. Equation (2.6) is strongly nonlinear and not easy to solve analytically, and most researchers have studied the classic Berman formula [8]; that is, when 𝛼=0. In this paper, we seek to solve (2.6) subject to the boundary conditions (2.7) and (2.8) using a novel spectral modification of the homotopy analysis method and the successive linearisation method. By comparison with the numerical approximations and previously obtained results, we show that these new techniques are accurate and more efficient than the standard homotopy analysis method.

3. Spectral Homotopy Analysis Method Solution

In applying the spectral-homotopy analysis method, it is convenient to first transform the domain of the problem from [0,1] to [−1,1] and make the governing boundary conditions homogeneous by using the transformations𝑦=𝜉+12,𝑈(𝜉)=𝐹(𝑦)−𝐹₀(𝑦),𝐹₀3(𝑦)=21𝑦−2𝑦³. (3.1) Substituting (3.1) in the governing equation and boundary conditions (2.6)–(2.8) gives16𝑈^1𝑉+8𝑎₁𝑈^′′′+4𝑎₂𝑈^′′+2𝑎₃𝑈^′−3Re𝑈+8Re(𝑈𝑈^′′′−𝑈^′𝑈^′′)=𝜙(𝑦), (3.2) subject to𝑈=0,𝑈^′′=0,𝜉=−1,𝑈=0,𝑈^′=0,𝜉=1, (3.3) where the primes denote differentiation with respect to 𝜉 and𝑎₁3=𝛼𝑦+Re(21𝑦−2𝑦³),𝑎₂3=3𝛼−2Re(1−𝑦²𝑎),₃=3𝑦𝑅e,𝜙(𝑦)=12𝛼𝑦+3Re𝑦³. (3.4)

The initial approximation is taken to be the solution of the nonhomogeneous linear part of the governing equations (3.2) given by16𝑈₀^1𝑉+8𝑎₁𝑈₀^′′′+4𝑎₂𝑈₀^′′+2𝑎₃𝑈^′₀−3Re𝑈₀=𝜙(𝑦), (3.5) subject to𝑈₀=0,𝑈₀^′′𝑈=0,𝜉=−1,₀=0,𝑈^′₀=0,𝜉=1. (3.6) We use the Chebyshev pseudospectral method to solve (3.5)–(3.6). The unknown function 𝑈₀(𝜉) is approximated as a truncated series of Chebyshev polynomials of the form𝑈₀(𝜉)≈U^𝑁₀(𝜉_𝑗)=^𝑁_𝑘=0𝑈_𝑘𝑇_𝑘(𝜉_𝑗),𝑗=0,1,…,𝑁, (3.7) where 𝑇_𝑘 is the 𝑘th Chebyshev polynomial, 𝑈_𝑘, are coefficients and 𝜉₀,𝜉₁,…,𝜉_𝑁 are Gauss-Lobatto collocation points (see [32]) defined by𝜉_𝑗=cos𝜋𝑗𝑁,𝑗=0,1,…,𝑁. (3.8) Derivatives of the functions 𝑈₀(𝜉) at the collocation points are represented as𝑑^𝑟𝑈₀𝑑𝜉^𝑟=^𝑁_𝑘=0𝒟^𝑟_𝑘𝑗𝑈₀(𝜉_𝑗), (3.9) where 𝑟 is the order of differentiation and 𝒟 is the Chebyshev spectral differentiation matrix ([32, 33]). Substituting (3.7)–(3.9) in (3.5)–(3.6) yields𝐀𝐔₀=Φ, (3.10) subject to theboundary conditions 𝑈₀(𝜉₀)=0,𝑈₀(𝜉_𝑁)=0,_𝑁𝑘=0𝒟²_𝑁𝑘𝑈₀(𝜉_𝑘)=0,^𝑁_𝑘=0𝒟_0𝑘𝑈₀(𝜉_𝑘)=0, (3.11)(3.12) where𝐴=16𝒟⁴+8𝐚₁𝒟³+4𝐚₂𝒟²+2𝐚₃𝐔𝒟−3Re𝐈,₀=[𝑈₀(𝜉₀),𝑈₀(𝜉₁),…,𝑈₀(𝜉_𝑁)]^𝑇,Φ=[𝜙(𝑦₀),𝜙(𝑦₁),…,𝜙(𝑦_𝑁)]^𝑇,𝐚_𝑠=diag([𝑎_𝑠(𝑦₀),𝑎_𝑠(𝑦₁),…,𝑎_𝑠(𝑦_𝑁−1),𝑎_𝑠(𝑦_𝑁)]),𝑠=1,2,3. (3.13)(3.14) In the above definitions, the superscript 𝑇 denotes transpose, diag is a diagonal matrix and 𝐈 is an identity matrix of size (𝑁+1)×(𝑁+1).

To implement the boundary conditions (3.11), we delete the first and the last rows and columns of 𝐴 and delete the first and last rows of 𝐔₀ and Φ. The boundary conditions (3.12) are imposed on the resulting first and last rows of the modified matrix 𝐀 and setting the resulting first and last rows of the modified matrix Φ to be zero. The values of [𝑈₀(𝜉₁),𝑈₀(𝜉₂),…,𝑈₀(𝜉_𝑁−1)] are then determined from 𝐔₀=𝐀⁻¹Φ. (3.15) To find the SHAM solutions of (3.2) we begin by defining the following linear operator:𝜕ℒ[𝑈(𝜉;𝑞)]=16⁴𝑈𝜕𝜉⁴+8𝑎₁𝜕³𝑈𝜕𝜉³+4𝑎₂𝜕²𝑈𝜕𝜉²+2𝑎₃𝜕𝑈𝜕𝜉−3Re𝑈, (3.16) where 𝑞∈[0,1] is the embedding parameter, and 𝑈(𝜉;𝑞) is an unknown function.

The zeroth order deformation equation is given by(1−𝑞)ℒ[𝑈(𝜉;𝑞)−𝑈₀(𝜉)]=𝑞ℏ{𝒩[𝑈(𝜉;𝑞)]−Φ}, (3.17) where ℏ is the nonzero convergence controlling auxiliary parameter and 𝒩 is a nonlinear operator given by𝜕𝒩[𝑈(𝜉;𝑞)]=16⁴𝑈𝜕𝜉⁴+8𝑎₁𝜕³𝑈𝜕𝜉³+4𝑎₂𝜕²𝑈𝜕𝜉²+2𝑎₃𝜕𝑈𝜕𝜕𝜉−3Re𝑈+8Re(𝑈³𝑈𝜕𝜉³−𝜕𝑈𝜕𝜕𝜉²𝑈𝜕𝜉²). (3.18)

Differentiating (3.17) 𝑚 times with respect to 𝑞 and then setting 𝑞=0 and finally dividing the resulting equations by 𝑚! yields the 𝑚th order deformation equationsℒ[𝑈_𝑚(𝜉)−𝜒_𝑚𝑈_𝑚−1(𝜉)]=ℏ𝑅_𝑚, (3.19) subject to the boundary conditions𝑈_𝑚(−1)=𝑈_𝑚(1)=𝑈_𝑚^′′(−1)=𝑈^′_𝑚(1)=0, (3.20) where 𝑅_𝑚(𝜉)=16𝑈^1𝑉_𝑚−1+8𝑎₁𝑈^′′′_𝑚−1+4𝑎₂𝑈^′′_𝑚−1+2𝑎₃𝑈^′_𝑚−1−3Re𝑈_𝑚−1+8Re^𝑚−1_𝑛=0(𝑈_𝑛𝑈^′′′_{𝑚−1−𝑛}−𝑈^′_𝑛𝑈^′′_{𝑚−1−𝑛})−𝜙(𝑦)(1−𝜒_𝑚𝜒),_𝑚=0,𝑚≤11,𝑚>1. (3.21)(3.22) Applying the Chebyshev pseudospectral transformation on (3.19)–(3.21) gives𝐀𝐔_𝑚=(𝜒_𝑚+ℏ)𝐀𝐔_𝑚−1−ℏ(1−𝜒_𝑚)Φ+ℏ𝑃_𝑚−1 (3.23) subject to the boundary conditions 𝑈_𝑚(𝜉₀)=0,𝑈_𝑚(𝜉_𝑁)=0,_𝑁𝑘=0𝒟²_𝑁𝑘𝑈_𝑚(𝜉_𝑘)=0,^𝑁_𝑘=0𝒟_0𝑘𝑈_𝑚(𝜉_𝑘)=0, (3.24)(3.25) where 𝐴 and Φ, are as defined in (3.13) and𝐔_𝑚=[𝑈_𝑚(𝜉₀),𝑈_𝑚(𝜉₁),…,𝑈_𝑚(𝜉_𝑁)]^𝑇,𝑃_𝑚−1=8Re^𝑚−1_𝑛=0[𝑈_𝑛(𝒟³𝑈_{𝑚−1−𝑛})−(𝒟𝑈_𝑛)(𝒟²𝑈_{𝑚−1−𝑛})]. (3.26) To implement the boundary conditions (3.24) we delete the first and last rows of 𝑃_𝑚−1 and Φ and delete the first and last rows and first and last columns of 𝐀 in (3.23). The boundary conditions (3.25) are imposed on the first and last row of the modified 𝐀 matrix on the left side of the equal sign in (3.23). The first and last rows of the modified 𝐀 matrix on the right of the equal sign in (3.23) are the set to be zero. This results in the following recursive formula for 𝑚≥1:𝐔_𝑚=(𝜒_𝑚+ℏ)𝐀⁻¹𝐀𝐔_𝑚−1+ℏ𝐀⁻¹[𝐏_𝑚−1−(1−𝜒_𝑚)Φ]. (3.27) Thus, starting from the initial approximation, which is obtained from (3.15), higher-order approximations 𝑈_𝑚(𝜉) for 𝑚≥1 can be obtained through the recursive formula (3.27).

4. Successive Linearisation Method

In this section, we solve (2.6) using the successive linearisation method. The main assumptions underpinning the use of the successive linearisation method are the following. (i)The unknown function 𝐹(𝑦) maybe expanded as 𝐹(𝑦)=𝐹_𝑖(𝑦)+^𝑖−1_𝑚=0𝐹_𝑚(𝑦),𝑖=1,2,3,…, (4.1) where 𝐹_𝑖 are unknown functions and 𝐹_𝑚 (𝑚≥1) are approximations which are obtained by recursively solving the linear part of the equation that results from substituting (4.1) in the governing equation (2.6).(ii)𝐹_𝑖 becomes progressively smaller as 𝑖 becomes large, that is,lim_𝑖→∞𝐹_𝑖=0. (4.2) Substituting (4.1) in the governing equation gives𝐹_𝑖^(𝑖𝑣)+𝑎_1,𝑖−1𝐹_𝑖^′′′+𝑎_2,𝑖−1𝐹_𝑖^′′+𝑎_3,𝑖−1𝐹^′_𝑖+𝑎_4,𝑖−1𝐹_𝑖+Re(𝐹_𝑖𝐹_𝑖^′′′−𝐹^′_𝑖𝐹_𝑖^′′)=𝑟_𝑖−1, (4.3) where the coefficient parameters 𝑎_{𝑘,𝑖−1}, (𝑘=1,2,3,4), and 𝑟_𝑖−1 are defined as𝑎_1,𝑖−1=Re^𝑖−1_𝑚=0𝐹_𝑚+𝛼𝑦,𝑎_2,𝑖−1=−Re^𝑖−1_𝑚=0𝐹^′_𝑚𝑎+3𝛼,_3,𝑖−1=−Re^𝑖−1_𝑚=0𝐹_𝑚^′′,𝑎_4,𝑖−1=Re^𝑖−1_𝑚=0𝐹_𝑚^′′′,𝑟_𝑖−1=−(^𝑖−1_𝑚=0𝐹_𝑚^(𝑖𝑣)+𝛼𝑦^𝑖−1_𝑚=0𝐹_𝑚^′′′+3𝛼^𝑖−1_𝑚=0𝐹_𝑚^′′)−Re(^𝑖−1_𝑚=0𝐹_𝑚^𝑖−1_𝑚=0𝐹_𝑚^′′′−^𝑖−1_𝑚=0𝐹^′_𝑚^𝑖−1_𝑚=0𝐹_𝑚^′′). (4.4)

The SLM algorithm starts from the initial approximation𝐹₀1(𝑦)=21(3+𝛽)𝑦−2(1+3𝛽)𝑦³+𝛽𝑦⁴, (4.5) which is chosen to satisfy the boundary conditions (2.7)–(2.8). The parameter 𝛽 in (4.5) is an arbitrary constant which when varied results in multiple solutions. The subsequent solutions for 𝐹_𝑚, 𝑚≥1 are obtained by successively solving the linearized form of (4.3) and which is given as𝐹_𝑖^(𝑖𝑣)+𝑎_1,𝑖−1𝐹_𝑖^′′′+𝑎_2,𝑖−1𝐹_𝑖^′′+𝑎_3,𝑖−1𝐹^′_𝑖+𝑎_4,𝑖−1𝐹_𝑖=𝑟_𝑖−1, (4.6) subject to the boundary conditions𝐹_𝑖(0)=0,𝐹_𝑖^′′(0)=0,𝐹_𝑖(1)=1,𝐹^′_𝑖(1)=0. (4.7) Once each solution for 𝐹_𝑖, (𝑖≥1) has been found from iteratively solving (4.6) for each 𝑖, the approximate solution for 𝐹(𝑦) is obtained as𝐹(𝑦)≈^𝑀_𝑚=0𝐹_𝑚(𝑦), (4.8) where 𝑀 is the order of SLM approximation. Since the coefficient parameters and the right hand side of (4.6), for 𝑖=1,2,3,…, are known from previous iterations, (4.6) can easily be solved using analytical means or any numerical methods such as finite differences, finite elements, Runge-Kutta-based shooting methods, or collocation methods. In this paper, (4.6) is integrated using the Chebyshev spectral collocation method [32, 33] as described in the previous section.

Applying the Chebyshev spectral method to (4.6) leads to the matrix equation𝐀_𝑖−1𝐅_𝐢=𝐑_𝑖−1, (4.9) in which 𝐀 is an (𝑁+1)×(𝑁+1) square matrix and 𝐘 and 𝐑 are (𝑁+1)×1 column vectors defined by𝐀=𝐃⁴+𝐚_1,𝑖−1𝐃³+𝐚_2,𝑖−1𝐃²+𝐚_3,𝑖−1𝐃+𝐚_4,𝑖−1,𝐑_𝑖−1=𝐫_𝑖−1, (4.10) with𝐅_𝑖=[𝐹_𝑖(𝑥₀),𝐹_𝑖(𝑥₁),…,𝐹_𝑖(𝑥_𝑁−1),𝐹_𝑖(𝑥_𝑁)]^𝑇,𝐫_𝑖−1=[𝑟_𝑖−1(𝑥₀),𝑟_𝑖−1(𝑥₁),…,𝑟_𝑖−1(𝑥_𝑁−1),𝑟_𝑖−1(𝑥_𝑁)]^𝑇. (4.11) In the above definitions, 𝑁 is the number of collocation points, 𝑥=2𝑦−1, 𝐚_{𝑘,𝑖−1}, (𝑘=1,2,3,4) are diagonal matrices of size (𝑁+1)×(𝑁+1), and 𝐃=2𝒟. After modifying the matrix system (4.9) to incorporate boundary conditions, the solution is obtained as𝐘_𝑖=𝐀⁻¹_𝑖−1𝐑_𝑖−1. (4.12)

5. Results and Discussion

In this section, we compare the results obtained using the various methods: the SHAM, the SLM, and the numerical approximations with those obtained using the HAM in Dinarvand and Rashidi [30] and the homotopy-Páde method in Xu et al. [24]. The solution obtained using most numerical solutions depends on the initial approximation. Using different initial guesses can give rise to multiple solutions. Multiple solutions were obtained if the initial guess in (4.5) is used in the SHAM method in place of 𝐹₀(𝑦) in (3.1). In this paper, it was observed that using different values of 𝛽 results in multiple solutions. For the multiple solutions comparison was made against the HAM results of [24].

An optimal ℏ value can easily be sought that can considerably improve the convergence rate of the results. However, for comparison purposes we used ℏ=−1. It is however worth noting, as pointed out in Dinarvand et al. [2], that when ℏ=−1, the solution series obtained by the HAM is the same solution series obtained by means of the homotopy perturbation method.

In Table 1 we compare the values of 𝐹(𝑦) when 𝛼=−1 and Re=−2,0 and 2 with the numerical and the HAM results reported in Dinarvand and Rashidi [30]. In [30], convergence up to six decimal places was achieved at the sixth order of the HAM approximation for Re=0 and 2. In this study, the same level of convergence and accuracy was achieved at the first order approximation for the same values of Re. For Re=−2 the convergence of the homotopy analysis method series solution was achieved at the eighth order of approximation while with the spectral homotopy analysis method series solution gives the same level of convergence at the second order.

Table 1 

Comparison of the numerical results against the SHAM approximate solutions for 𝐹(𝑦) when 𝛼=−1 with 𝑁=60 and ℏ=−1.

In Table 2, we demonstrate the computational efficiency of the SHAM solution for large values of Re. As has been noted in the introduction, some semi-analytical methods fail to converge at large values of Re, for example, Dinarvand et al. [2] have shown that for |Re|>9.5 the HPM fails to converge. However, in using the SHAM convergence up seven decimal places is achieved at the third order of approximation for values of Re as large as Re=150. For 5≤Re<100 convergence up six decimal places is achieved at the second order. In Table 3, Re=2 is fixed and 𝐹^′′(1) evaluated for −2.5≤𝛼≤2.5. Convergence up to seven digits is achieved at the first order for 𝛼=0, at the second order for −2≤𝛼≤−0.5, and at the third-order approximation for 𝛼=−2.5. Clearly the SHAM gives faster convergence than the HAM under the same conditions.

Table 2 

Comparison of the numerical results against the SHAM approximate solutions for 𝐹^′′(1) when 𝛼=−1 with 𝑁=60 and ℏ=−1 for different values of Re.

Table 3 

Comparison of the numerical results against the SHAM approximate solutions for 𝐹^′′(1) when 𝛼=−1 with Re=10 and ℏ=−1 for different values of 𝑁.

Table 3 also gives a comparison of the SHAM approximate results with the numerical results generated using different values of 𝑁. It is evident however that for small values of 𝑁 (say less than 𝑁=20) the SHAM results are not very accurate. Accuracy however improves with an increase in 𝑁.

Table 4 gives a comparison of the numerical and the SHAM results of 𝐹^′′(1) when −2.5≤𝛼≤2.5. Convergence is generally achieved at either the third- or the fourth-order of the SHAM approximation. Figures 1 and 2 give a comparison of the numerical and the SHAM solutions for the characteristic mean-flow function 𝐹(𝑦)=−𝜈/𝑐 and 𝐹^′(𝑦)=𝜇𝑐/𝑥 at different Reynolds numbers and 𝛼. The mean-flow function 𝐹(𝑦), increases with increasing (positive) values of Re and 𝛼 while 𝐹^′(𝑦) decreases, which makes sense since 𝐹(𝑦) is inversely proportional to the injection or suction coefficient 𝑐=𝛼/Re while 𝐹^′(𝑦) is directly proportional. Figure 1 further illustrates the efficiency of the solution method with excellent agreement for Re as large as 200.

Table 4 

Comparison of the numerical results against the SHAM approximate solutions for 𝐹^′′(1) when Re=2 with 𝑁=60 and ℏ=−1 for different values of 𝛼.

(a)

(b)

Figure 1 

Comparison between numerical and SHAM approximate solution of $F(y)$ and $F^{\prime }(y)$ for different values of $\mathrm{Re}$ when $\alpha =-1$ when $\hslash =-1.14$ (for $\mathrm{Re}=-10$) and $\hslash =-1$ (for $\mathrm{Re}=\mathrm{0,200}$).

(a)

(b)

Figure 2 

Comparison between numerical and SHAM approximate solution of $F(y)$ and $F^{\prime }(y)$ for different values of $\alpha$ when $\mathrm{Re}=2$ when $\hslash =-1$ (for $\alpha =-5,-\mathrm{2.5,0},2.5$) and $\hslash =-0.94$ (for $\alpha =5$).

Using the initial approximation 𝐹₀(𝑦) in (4.5) with different values of 𝛽 in place of (3.1) in the SHAM implementation leads to multiple solutions. When 𝛽=0 and −5, the SHAM gives the multiple solutions observed in Xu et al. [24]. A comparison of the SHAM results against the HAM results reported in [24] is presented in Table 5. It is evident that the SHAM results converge much more rapidly than the HAM results of [24] for both branches of the solution.

Table 5 

Comparison between the multiple solutions of HAM results (see [24]) and the present SHAM results in the case of Re=−10 and 𝛼=4.

Tables 6 and 7 give, first, the analytical approximations of 𝐹^′(0)/Re and 𝐹^′′′(0)/Re for the two solutions obtained using the successive linearisation method. Secondly, the tables give a direct comparison of the convergence rates of the SLM and the [𝑚,𝑚] homotopy-Padé method used by Xu et al. [24]. The fourth-order SLM approximation gives the same level of accuracy as the twenty-fourth-order of the [𝑚,𝑚] homotopy-Páde approximation, which suggests that the successive linearisation method is much more computationally efficient and accurate compared to the [𝑚,𝑚] homotopy-Páde approximation (although it is not clear at this stage whether this could be attributed to the use of a more efficient initial guess).

Table 6 

Comparison of 𝐹^′(0)/Re and 𝐹^′′′(0)/Re obtained at different orders for the SLM, and the [𝑚,𝑚] homotopy-Páde approximation when Re=−10 and 𝛼=4. For the SLM first solution, 𝛽=−1, and 𝑁=25, and for the SLM second solution, 𝛽=−5, and 𝑁=25.

Table 7 

Comparison of 𝐹^′(0)/Re and 𝐹^′′′(0)/Re obtained at different orders for the SLM, and the [𝑚,𝑚] homotopy-Páde approximation when Re=−11 and 𝛼=3/4. For the SLM first solution, 𝛽=−5, and 𝑁=25, 𝛽=1, and 𝑁=25 for the SLM second solution, and for the SLM third solution 𝛽=3, and 𝑁=25.

6. Conclusion

In this paper, we have used the spectral homotopy analysis method (SHAM) and the successive linearisation method (SLM) to solve a fourth-order nonlinear boundary value problem that governs the two-dimensional Laminar flow between two moving porous walls. Multiple solutions recently reported in Xu et al. [24] have been obtained, depending on the initial approximation used. Comparison of the computational efficiency and accuracy of the results between the current methods and the previous homotopy analysis method results described in Dinarvand and Rashidi [30] and Xu et al. [24] has been made. Our simulations show that the convergence of the SHAM solution series to the numerical solution (up to six decimal place accuracy) is achieved at the second order (for Re=−2) and first order for 𝑅=0,2. In contrast to the standard homotopy analysis method (c.f. Dinarvand and Rashidi [30]) convergence was achieved at the eighth order (for Re=−2) and sixth order for Re=0,2. The SHAM is apparently more efficient because it offers more flexibility in choosing linear operators compared to the standard HAM. It is however important to note that if the same initial guess and linear operators were to be used, the two methods would give the same solution.

We have further shown that notwithstanding the acceleration of the convergence ratio of the homotopy series by means of the homotopy-Padé technique, the successive linearisation techniques is more computationally efficient (although this again could be due, in part, to the use of a different initial guess) and gives accurate results.

Acknowledgments

The authors wish to acknowledge financial support from the University of KwaZulu-Natal and the National Research Foundation (NRF). The authers also thank the anonymous reviewers whose helpful comments have contributed to the improvement of our work.