A Note on Improved Homotopy Analysis Method for Solving the Jeffery-Hamel Flow

Abstract

This paper presents the solution of the nonlinear equation that governs the flow of a viscous, incompressible fluid between two converging-diverging rigid walls using an improved homotopy analysis method. The results obtained by this new technique show that the improved homotopy analysis method converges much faster than both the homotopy analysis method and the optimal homotopy asymptotic method. This improved technique is observed to be much more accurate than these traditional homotopy methods.

1. Introduction

The mathematical study of the flow of a viscous incompressible two-dimensional fluid in a wedge-shaped channel with a sink or source at the vertex was pioneered by Jeffery [1] and Hamel [2]. The problem has since been studied extensively by, among others, Axford [3] who included the effects of an externally applied magnetic field and Rosenhead [4] who obtained a general solution containing elliptic functions.

Instability and bifurcation are other aspects of the Jeffery-Hamel problem that have attracted widespread interest; see, for example, Akulenko and Kumakshev [5, 6]. Three-dimensional extensions to and bifurcations of the Jeffery-Hamel flow have been made by Stow et al. [7] while McAlpine and Drazin [8] presented a normal mode analysis of two-dimensional perturbations of a viscous incompressible fluid driven between inclined plane walls by a line source at the intersection of the walls. Banks et al. [9] investigated various perturbations and the linear temporal stability of such flows and found evidence of a strong interaction between the disturbances up- and down-stream if the angle between the planes exceeds a certain Reynolds-number-dependent critical value. Makinde and Mhone [10] investigated the temporal stability of MHD Jeffery-Hamel flows. They showed that an increase in the magnetic field intensity has a strong stabilizing effect on both diverging and converging channel geometries. A review of the theory of instabilities and bifurcations in channels is given by Drazin [11].

As with most problems in science and engineering, the equations governing the Jeffery-Hamel problem are highly nonlinear and so generally do not have closed form analytical solutions. Nonlinear equations can, in principle, be solved by any one of a wide variety of numerical methods. However, numerical solutions rarely give intuitive insights into the effects of various parameters associated with a problem. Consequently, most recent studies of flows in diverging and converging channels have centred on the use of the Jeffery-Hamel flow equations as a testing and proving tool for the accuracy, reliability, and robustness of new techniques for solving nonlinear equations. Examples of such techniques include the summation series technique [12], the homotopy analysis method [13, 14], the decomposition method [15], the homotopy perturbation method [16], Hermite-Padé approximation [17], and the spectral-homotopy analysis method [18–21]. The study by Joneidi et al. [14] used three methods: the differential transform method (DTM), the homotopy analysis method (HAM), and the homotopy perturbation method (HPM) to solve the Jeffery-Hamel problem. The study confirmed that although both the DTM and the HPM give acceptable accuracy, the HAM is by far the superior method delivering faster convergence and better accuracy. Nonetheless, important improvements have been made to the HAM by Motsa et al. [18, 19]. The spectral modification of the HAM or SHAM proposed by Motsa et al. [18, 19] removes some of the prescriptive assumptions associated with the HAM and further accelerates the convergence rate of the method.

A recent study by Esmaeilpour and Ganji [22] reported on the solution of the Jeffery-Hamel problem using an optimal homotopy asymptotic method (OHAM). This method was developed by Marinca and Herişanu [23] for the approximate solution of nonlinear problems of thin film flow of a fourth-grade fluid. The OHAM was used by these authors and others to solve the equations for the steady flow of a fourth-grade fluid in a porous medium [24, 25] and has since been applied to many other nonlinear problems including the squeezing flow problem by Idrees et al. [26]. Studies by Islam et al. [27, 28] further suggest that the OHAM is more general than both the HAM and the HPM with the latter methods being special cases of the OHAM.

In this paper, we report on a new and improved method known as improved homotopy analysis method (IHAM) for solving general boundary value problems. The IHAM is an algorithm that seeks to improve the initial approximation that is later used in the HAM to solve the governing nonlinear equation resulting in significant improvement in the accuracy and convergence rate of the solutions. We seek to demonstrate the application of this method by solving the Jeffery-Hamel problem and to show its accuracy and rapid convergence by comparing the present solutions with those in the literature, including the HAM, SHAM, and the OHAM. The IHAM solutions are further compared with numerical solutions. Finally, we believe that this work will motivate further improvements to the homotopy-based semi-analytical methods.

2. Governing Equations

We consider the steady two-dimensional flow of a viscous incompressible fluid between two nonparallel rigid planes with angle 2𝛼. Assuming radial flow with velocity 𝐯=𝑢(𝑟,𝜃), the motion of such a fluid is described by the equations (see [14, 19, 22])𝜌𝑟𝜕()𝑢𝜕𝑟𝑟𝑢=0,𝜕𝑢1𝜕𝑟=−𝜌𝜕𝑝𝜕𝜕𝑟+𝜈²𝑢𝜕𝑟²+1𝑟𝜕𝑢+1𝜕𝑟𝑟²𝜕²𝑢𝜕𝜃²−𝑢𝑟²−1,𝜌𝑟𝜕𝑝+𝜕𝜃2𝜈𝑟²𝜕𝑢𝜕𝜃=0, (2.1)(2.2)(2.3) where 𝜌 is the fluid density, 𝑝 the pressure, and 𝜈 the kinematic viscosity. The continuity equation implies that ℎ(𝜃)=𝑟𝑢(𝑟,𝜃), (2.4) and by using the nondimensional parameters 𝑓𝑦=𝑓(𝜃)𝑓_max𝜃,𝑦=𝛼, (2.5) Equations (2.2)-(2.3) reduce to𝑓^′′′+2𝛼Re𝑓𝑓′+4𝛼²𝑓′=0, (2.6) where Re=𝛼𝑓_max𝜈=𝑈_max𝑟𝛼𝜈 (2.7) is the Reynolds number and 𝑈_max is the maximum velocity at the centre of the channel. The appropriate boundary conditions are𝑓(0)=1,𝑓^′(0)=0,𝑓(1)=0. (2.8)

3. Improved Homotopy Analysis Method (IHAM) Solution

In this section, we describe the use of the improved homotopy analysis method (IHAM) in the governing equation (2.6). To apply the IHAM, we assume that the solution 𝑓(𝑦) can be expanded as𝑓𝑦=𝑓_𝑖𝑦+^𝑖−1_𝑛=0𝑓_𝑛𝑦,𝑖=1,2,3,…, (3.1) where 𝑓_𝑖 are unknown functions whose solutions are obtained using the HAM approach at the 𝑖th iteration and 𝑓_𝑛, (1≤𝑛≤𝑖−1) are known from previous iterations. The algorithm starts with the initial approximation 𝑓₀(𝑦) which is chosen to satisfy the boundary conditions (2.8). An appropriate initial guess is𝑓₀𝑦=1−𝑦². (3.2) Substituting (3.1) in the governing equation (2.6)-(2.8) gives𝑓_𝑖^′′′+𝑎_1,𝑖−1𝑓^′_𝑖+𝑎_2,𝑖−1𝑓_𝑖+2𝛼Re𝑓_𝑖𝑓^′_𝑖=𝑟_𝑖−1, (3.3) subject to the boundary conditions𝑓_𝑖(0)=0,𝑓^′_𝑖(0)=0,𝑓_𝑖(1)=0, (3.4) where the coefficient parameters 𝑎_{𝑘,𝑖−1}, (𝑘=1,2) and 𝑟_𝑖−1 are defined as𝑎_1,𝑖−1=2𝛼Re^𝑖−1_𝑛=0𝑓_𝑛+4𝛼²,𝑎_2,𝑖−1=2𝛼Re^𝑖−1_𝑛=0𝑓^′_𝑛,𝑟_𝑖−1=−^𝑖−1_𝑛=0𝑓_𝑛^′′′+2𝛼Re^𝑖−1_𝑛=0𝑓^′_𝑛^𝑖−1_𝑛=0𝑓_𝑛+4𝛼^2𝑖−1_𝑛=0𝑓^′_𝑛. (3.5) Starting from the initial approximation (3.2), the subsequent solutions 𝑓_𝑖 (𝑖≥1) are obtained by recursively solving equation (3.3) using the HAM approach [29, 30]. To find the HAM solutions of (3.3), we begin by rewriting (3.3) as𝒩𝑓_𝑖𝑦=𝑟_𝑖−1, (3.6) where𝒩is a nonlinear operator defined by 𝒩𝑓_𝑖𝑦=𝑓_𝑖^′′′+𝑎_1,𝑖−1𝑓^′_𝑖+𝑎_2,𝑖−1𝑓_𝑖+2𝛼Re𝑓_𝑖𝑓^′_𝑖. (3.7) Let 𝑓_𝑖,0(𝑦) denote the initial guess for the unknown function 𝑓_𝑖(𝑦), and let ℏ≠0 be an auxiliary parameter. Using an embedding parameter 𝑞∈[0,1] we construct a homotopy ℋ𝐹_𝑖𝑦;𝑞;𝑓_𝑖,0𝑦,ℏ,𝑞=1−𝑞ℒ𝐹_𝑖𝑦;𝑞−𝑓_𝑖,0𝑦−𝑞ℏ𝒩𝐹_𝑖𝑦;𝑞−𝑟_𝑖−1, (3.8) where𝐿is an auxiliary linear operator with the property that ℒ[0]=0 and 𝐹_𝑖(𝑦;𝑞) is an unknown function. Upon equating ℋ to 0, we obtain the zero-order deformation equation 1−𝑞ℒ𝐹_𝑖𝑦;𝑞−𝑓_𝑖,0𝑦=𝑞ℏ𝒩𝐹_𝑖𝑦;𝑞−𝑟_𝑖−1. (3.9) When, 𝑞=0 (3.9) becomes ℒ𝐹_𝑖𝑦;0−𝑓_𝑖,0𝑦=0. (3.10) This equation holds provided 𝐹_𝑖𝑦;0=𝑓_𝑖,0𝑦 (3.11) as ℒ[0]=0. When𝑞=1, (3.9) is simplified to 𝒩𝐹_𝑖𝑦;1=𝑟_𝑖−1 (3.12) asℏ≠0. Equation (3.12) is the same as equation (3.6) provided 𝐹_𝑖𝑦;1=𝑓_𝑖𝑦. (3.13) It follows from (3.11) and (3.13) that as 𝑞 increases from 0 to 1, the unknown function 𝐹(𝑦;𝑞) varies continuously from initial guess 𝑓_𝑖,0(𝑦) to exact solution 𝑓_𝑖(𝑦) of (3.6).

Differentiating (3.9) 𝑚 times with respect to 𝑞 and then setting 𝑞=0 and finally dividing the resulting equations by 𝑚! yields the 𝑚th-order deformation equationsℒ𝑓_𝑖,𝑚𝜂−𝜒_𝑚𝑓_{𝑖,𝑚−1}=ℏ𝑓^′′′_{𝑖,𝑚−1}+𝑎_1,𝑖−1𝑓^′_{𝑖,𝑚−1}+𝑎_2,𝑖−1𝑓_{𝑖,𝑚−1}+2𝛼Re^𝑚−1_𝑗=0𝑓_𝑖,𝑗𝑓^′′_{𝑖,𝑚−1−𝑗}−1−𝜒_𝑚𝑟_𝑖−1, (3.14) subject to the boundary conditions𝑓_𝑖,𝑚(0)=𝑓^′_𝑖,𝑚(0)=𝑓_𝑖,𝑚(1)=0, (3.15) where𝜒_𝑚=0,𝑚≤1,1,𝑚>1.. (3.16) If we define the linear operator ℒ byℒ𝐹_𝑖=𝜕𝑦;𝑞³𝐹_𝑖𝜕𝑦³, (3.17) then the initial approximation 𝑓_𝑖,0 that is used in the higher-order equations (3.14) is obtained by solving the differential equationℒ𝑓_𝑖,0=𝑓^′′′_𝑖,0=𝑟_𝑖−1, (3.18) subject to the boundary conditions𝑓_𝑖,0(0)=𝑓^′_𝑖,0(0)=𝑓_𝑖,0(1)=0. (3.19) Thus, starting from the initial approximation, which is obtained-from (3.18), higher order approximations 𝑓_𝑖,𝑚(𝑦) for 𝑚≥1 can be obtained through the recursive formula (3.14). We note that (3.14) forms a set of linear ordinary differential equations and can be easily solved analytically, especially by means of symbolic computation software such as Maple, Mathematica, Matlab, and others.

Expanding 𝐹_𝑖(𝑦;𝑞) in Taylor series about 𝑞=0 gives𝐹_𝑖𝑦;𝑞=𝐹_𝑖𝑦;0+^+∞_𝑘=1𝑞^𝑘𝜕𝑘!^𝑘𝐹_𝑖𝑦;𝑞𝜕𝑦^𝑘_𝑞=0. (3.20) If we set 𝑞=1 in (3.20) and define 𝑓_𝑖,𝑘𝑦1∶=𝜕𝑘!^𝑘𝐹_𝑖𝑦;𝑞𝜕𝑦^𝑘_𝑞=0 (3.21) for each 𝑘=1,2,…, then making use of (3.11) and (3.13) transforms (3.20) to 𝑓_𝑖𝑦=𝑓_𝑖,0𝑦+^+∞_𝑘=1𝑓_𝑖,𝑘𝑦. (3.22) Upon truncating the infinite series in (3.22), the solutions for 𝑓_𝑖 are then generated using the solutions for 𝑓_𝑖,𝑚 as follows:𝑓_𝑖=𝑓_𝑖,0+𝑓_𝑖,1+𝑓_𝑖,2+𝑓_𝑖,3+𝑓_𝑖,4+⋯+𝑓_𝑖,𝑚. (3.23) The [𝑖,𝑚] approximate solution for 𝑓(𝑦) is then obtained by substituting 𝑓_𝑖 (obtained from (3.23)) in (3.1).

4. Results and Discussion

The accuracy and reliability of the IHAM is determined by comparing the current results with the optimal homotopy analysis method (OHAM) of Esmaeilpour and Ganji [22], the homotopy analysis results (HAM) of Joneidi et al. [14], and the numerical results obtained using the MATLAB bvp4c routine which is a boundary value solver based on the adaptive Lobatto quadrature scheme [31, 32]. In Table 1 we show a comparison of the current velocity results at various dimensionless angles 𝑦 for different orders [𝑚,𝑛] of the IHAM. The Reynolds number is fixed at Re=50, and we consider a diverging channel with angle 𝛼=5^∘. Convergence of the current results to the numerical results to six decimal places is achieved at order [2,2] while convergence to ten decimal places is achieved at order [2,3]. On the other hand, convergence of the OHAM is evidently slow, and up to the twentieth order of approximation, the method mostly fails to converge to the numerical results.

Table 1 

Comparison between the numerical results and the order [𝑚,𝑛] IHAM approximate results (using ℏ=−1) for 𝑓(𝑦) against the OHAM results reported in [22] when Re=50, 𝛼=5.

The slow convergence of the OHAM is confirmed in Table 2 where the absolute errors of the two methods in relation to the numerical solution are given. At any angle 𝑦, the absolute error using the OHAM is much larger than that obtained using the IHAM. The largest absolute error at any angle 𝑦 using the IHAM is 5.6×10⁻⁷ while, by comparison, the largest absolute error when using the OHAM is 6.958×10⁻³.

Table 2 

Comparison between the errors of the order [𝑚,𝑛] IHAM results against the OHAM results reported in [22] for Re=50 and 𝛼=5.

Table 3 gives a comparison of the convergence rate of the current method with the HAM [14]. Convergence of the IHAM to the numerical results is found to be rapid, with agreement to ten decimal places being archived at order [3,3] of the IHAM algorithm. The sixteenth-order HAM on the other hand fails to achieve the accuracy of both the IHAM and the numerical results. The poor performance of the HAM is further confirmed in Table 4 where the largest absolute error incurred by the IHAM at order [3,2] is 7.2×10⁻⁹ while the largest absolute error incurred by using the sixteenth-order HAM at any angle 𝑦 is approximately a hundred times larger at 4.7×10⁻⁷.

Table 3 

Comparison between the numerical results and the order [𝑚,𝑛] IHAM approximate results (using ℏ=−1) for 𝑓(𝑦) against the HAM results reported in [14] when Re=110, 𝛼=3.

Table 4 

Comparison between the errors of the order [𝑚,𝑛] IHAM results against the HAM results reported in [14] for Re=110 and 𝛼=3.

Figure 1 shows the effect of the Reynolds number on the fluid velocity for different values of 𝛼. It can be seen from the figure that the fluid velocity increases with Reynolds numbers in the case of convergent channels (𝛼<0) but decreases with Re in the case of divergent channels (𝛼>0). This observation is consistent with the observations made in [13, 14, 19].

(a)

(b)

Figure 1 

Velocity profile 𝑓(𝑦) for different values of Re when 𝛼=5 and 𝛼=−5, respectively.

Figure 2 shows the effect of the Reynolds number on 𝑓^′′(0), which is related to wall shear stress, for different values of 𝛼. It can be seen from the figure that 𝑓^′′(0) decreases monotonically for 𝛼>0. In the case of 𝛼<0, 𝑓^′′(0) increases for a range of Re until it reaches a peak then decreases.

(a)

(b)

Figure 2 

Variation of 𝑓(0) against the Reynolds number when 𝛼=5 and 𝛼=−5, respectively.

5. Conclusion

In this brief note we have proposed an improved homotopy analysis method (IHAM) for the solution of general nonlinear differential equations. We have compared the performance of the new algorithm against the optimal homotopy analysis method (OHAM), the standard homotopy analysis method (HAM), and numerical approximations by solving the Jeffery-Hamel problem for Newtonian flow in converging/diverging channels. Numerical computations show that the IHAM is accurate and converges to the numerical approximations at lower orders compared to both the HAM and the OHAM. However, at this juncture, we cannot say with certainty that this method is better than these other existing methods. We recommend that we need to use this improved homotopy analysis method (IHAM) to solve more nonlinear differential equations.

Acknowledgments

The authors wish to acknowledge financial support from the University of Swaziland, University of KwaZulu-Natal, University of Venda, and the National Research Foundation (NRF).