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

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.


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 downstream 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 Heris ¸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.

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 v u r, θ , the motion of such a fluid is described by the equations see 14, 19, 22 1 r ∂u ∂r where ρ is the fluid density, p the pressure, and ν the kinematic viscosity.The continuity equation implies that and by using the nondimensional parameters where is the Reynolds number and U max is the maximum velocity at the centre of the channel.The appropriate boundary conditions are

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 f y can be expanded as where f i are unknown functions whose solutions are obtained using the HAM approach at the ith iteration and f n , 1 ≤ n ≤ i − 1 are known from previous iterations.The algorithm starts with the initial approximation f 0 y which is chosen to satisfy the boundary conditions 2.8 .An appropriate initial guess is Substituting 3.1 in the governing equation 2.6 -2.8 gives subject to the boundary conditions where the coefficient parameters a k,i−1 , k 1, 2 and r i−1 are defined as

3.5
Starting from the initial approximation 3.2 , the subsequent solutions f i i ≥ 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 where N is a nonlinear operator defined by Let f i,0 y denote the initial guess for the unknown function f i y , and let / 0 be an auxiliary parameter.Using an embedding parameter q ∈ 0, 1 we construct a homotopy where L is an auxiliary linear operator with the property that L 0 0 and F i y; q is an unknown function.Upon equating H to 0, we obtain the zero-order deformation equation 1 − q L F i y; q − f i,0 y q N F i y; q − r i−1 . 3.9 When, q 0 3.9 becomes This equation holds provided as L 0 0. When q 1, 3.9 is simplified to as / 0. Equation 3.12 is the same as equation 3.6 provided F i y; 1 f i y .

3.13
It follows from 3.11 and 3.13 that as q increases from 0 to 1, the unknown function F y; q varies continuously from initial guess f i,0 y to exact solution f i y of 3.6 .Differentiating 3.9 m times with respect to q and then setting q 0 and finally dividing the resulting equations by m! yields the mth-order deformation equations subject to the boundary conditions where .

3.16
If we define the linear operator L by L F i y; q ∂ 3 F i ∂y 3 , 3.17 then the initial approximation f i,0 that is used in the higher-order equations 3.14 is obtained by solving the differential equation subject to the boundary conditions Thus, starting from the initial approximation, which is obtained-from 3.18 , higher order approximations f i,m y for m ≥ 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 F i y; q in Taylor series about q 0 gives

3.20
If we set q 1 in 3.20 and define for each k 1, 2, . .., then making use of 3.11 and 3.13 transforms 3.20 to

3.22
Upon truncating the infinite series in 3.22 , the solutions for f i are then generated using the solutions for f i,m as follows:

3.23
The i, m approximate solution for f y is then obtained by substituting f i obtained from 3.23 in 3.1 .

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 y for different orders m, n 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.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 y, the absolute error using the OHAM is much larger than that obtained using the IHAM.The largest absolute error at any angle y using the IHAM is 5.6 × 10 −7 while, by comparison, the largest absolute error when using the OHAM is 6.958 × 10 −3 .
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 −9 while the largest absolute error incurred by using the sixteenth-order HAM at any angle y is approximately a hundred times larger at 4.7 × 10 −7 .
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 .Figure 2 shows the effect of the Reynolds number on f 0 , which is related to wall shear stress, for different values of α.It can be seen from the figure that f 0 decreases monotonically for α > 0. In the case of α < 0, f 0 increases for a range of Re until it reaches a peak then decreases.

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 Mathematical Problems in Engineering 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.

Table 1 :
Comparison between the numerical results and the order m, n IHAM approximate results using −1 for f y against the OHAM results reported in 22 when Re 50, α 5.

Table 2 :
Comparison between the errors of the order m, n IHAM results against the OHAM results reported in 22 for Re 50 and α 5.

Table 3 :
Comparison between the numerical results and the order m, n IHAM approximate results using −1 for f y against the HAM results reported in 14 when Re 110, α 3.

Table 4 :
Comparison between the errors of the order m, n IHAM results against the HAM results reported in 14 for Re 110 and α 3.