Highly accurate Coiﬂet wavelet-homotopy solution of Jeffery-Hamel problem at extreme parameters

The magnetised Jeffery-Hamel ﬂow due to a point sink or source in convergent and divergent channels is studied. The simpliﬁed governing equation ruled by the Reynolds number, the Hartmann number and the divergent-convergent angle with appropri-ate boundary conditions are solved by the newly proposed Coiﬂet wavelet-homotopy method. Highly accurate solutions are obtained, whose accuracy is rigidly checked. As compared with the traditional homotopy analysis method, our proposed technique has higher computational efﬁciency and larger applicable range of physical parameters. Results show that our proposed technique is very convenient to handle strong nonlinear problems without special treatment. It is expected that this technique can be further applied to study complex nonlinear problems in science and engineering involving into extreme physical parameters. Besides, the inﬂuence of physically important quantities on the ﬂow is discussed. It is found that wall stretching and shrinking exhibits totally different roles on the ﬂow development. The enhanced Lorenz force affects the ﬂow behaviours signiﬁcantly for both convergent and divergent cases.


Introduction
Incompressible viscous flow driven in convergent and divergent channels has numerous applications in the fields of fluid mechanics and civil, mechanical, bio-mechanical and environmental engineering. The famous Jeffery-Hamel problem [1,2] is usually used to describe such fluid motion, which was then extended by many researchers to different physical configurations such as radial flow between two inclined plane walls [3], flow in symmetrical channels with slightly curved walls [4], magnetohydrodynamic flow between convergent or divergent channels [5]. Like most classical problems in fluid mechanics, the class of Jeffery-Hamel problems does not admit analytical solutions. As a result, analytical and numerical approaches are common tools to give their solutions. In literature, many techniques including the spectral homotopy analysis method [6], the optimal homotopy asymptotic method [7], the Adomian decomposition method [8,9] and its modification [10], the perturbation method [11], the homotopy analysis method [12] and its inconsequential transformation [13,14], the predictor homotopy analysis method [15], the finite volume method [16] and the collocation method [17], have been suggested to attack suck kind of convergent or divergent channel flow problems.
Among those computational techniques, the homotopy analysis method [18,19] was regarded as a very accurate and efficient approach for handling nonlinear problems with strong nonlinearity. This is due to that it is not dependent upon small or large physical parameters, has freedom to express the solution based on its behaviours, and can adjust convergent rate via a convergence-control parameter.
Though having all above mentioned advantages, in our knowledge there is still a room to improve on this technique. Several obvious weaknesses of the homotopy analysis method are held. Firstly, it is very time-consuming for complicated nonlinear problems. Secondly, its convergent rate is rather slow when extreme large parameters are involved. Thirdly, it has limited solution expression forms. Yang and Liao [20,21] made efforts to overcome these issues by introducing the Coiflet wavelet functions [22][23][24] into the homotopy analysis method as the solution expressions. Their modified technique not only holds the advantages of strong nonlinear processing capability inherited by the homotopy analysis method and excellent local expression characteristics originated from the Coiflet wavelet, but also is very efficient and accurate for nonlinear equations with homogenous boundary conditions. Yu and Xu [25][26][27] established a generalized boundary modification approach based on the Coiflet wavelet which is suitable for both ordinary and partial differential equations subjected to non-homogenous boundary conditions. Their idea was further adopted by Wang et al. [28] for a electrohydro-dynamic flow problem and Chen [29] for a channel flow problem, respectively.
In this paper, we shall further develop the Coiflet wavelet-homotopy method for handling nonlinear problems with extreme physical parameters. A complete set of solution solving process will be established based on both the homotopy analysis method and the Coiflet wavelet based on the famous MHD Jeffery-Hamel problem. Note that in most above mentioned studies, few researchers have considered such extreme case. This makes our work novel and unique. Highly accurate solution for a large range of physical parameters will be given, whose validity and efficiency will be strictly checked. The outline of this is organized as follows: In Section 2, mathematical modeling and formulations of the MHD Jeffery Hamel flow is presented. In Section 3, the fundamentals of Coiflet wavelet-homotopy method is given.
In Section 4, results is verified and. In Section 5, the concluding remarks are drawn.

Mathematical Description
We consider a magnetohydrodynamic flow of an incompressible viscous conducting fluid due to either a source or a sink at the intersection of two stretchable or shrinkable nonparallel plane walls with an included angle 2β. The cases β < 0 and β > 0, respectively, denote the convergent and the divergent channels. The physical sketch is shown in Fig.1. The polar coordinates (r, θ, z) are employed, where r, θ and z denote the polar radius, polar angle and polar axis respectively. The flow is assumed to be laminar and symmetrical. The variable magnetic field B is assumed to be perpendicular to the radial direction, whose distribution is [B 0 /r, 0, 0]. The fluid parameters are assumed to remain unchanged along z direction, as a result the flow is only dependent on r and θ. These assumptions indicate that the velocity field of the flow has the form V = [u r , 0, 0], where u r = u(r, θ) is a function of r and θ.
The governing equations for the flow are given, as suggested by Jeffery [1], as where t denotes the time, V, ρ, µ, p and F are the velocity vector, the fluid density, the dynamic viscosity, the pressure and the Lorenz force (total electromagnetic force) respectively. Here the Laplacian operator is given by Based on the Ohm's low, the Lorenz force is expressed by where J is the electrical current density, defined by in which σ is the electrical conductivity, E is the applied electrical field that is enforced to zero (refer to Xu et al. [31]).
On the other hand, the walls are assumed to be radially stretched or shrunk in accordance with where S r is the stretching/shrinking rate. On the centreline, it holds where u c is a reference velocity on the centreline.
For the steady case, the aforementioned equations are expanded, after a series of physical assumptions such as low magnetic Reynolds number and symmetrical flow, as 1 ρr where ν is the kinematic viscosity. Note that here the strongly radial flow assumption is applied so that it holds u θ = 0. Then the scaling parameters are defined as By means of these transformations, the continuity equation (8) is automatically satisfied, and the set of momentum equations (9) and (10) is reduced to subject to the boundary conditions where the prime denotes the derivative to η, C w , R e and H a are the stretching or shrinking parameter, the Reynolds number and the Hartmann number, which are defined by Note that C w > 0 and C w < 0 respectively represent the stretching and shrinking wall cases.
Physically important quantity for this problem is the local skin friction coefficient, defined by where τ w is the wall shear stress, governed by Substituting Eq.(16) into Eq. (17), with considering the scaling transformation (11), we obtain It is clear that f ′ (1) is the core part of the local skin friction.

Generalized Coiflets wavelet expression
Here the Coiflet wavelet modification technique suggested by Yu et al. [27] is employed to solve Eq. (12).
By this approach, the function g(x) (x ∈ [0, 1]) is expressed, for a prescribed resolution level j, by where k denotes the kth interpolation point, and g(x k ) is the value of g(x) at x k . α 0,s and α 1,s , ϖ 0,s (x) and ϖ 1,s (x) are coefficients associated with non-homogenous boundary conditions. For known boundary conditions at certain derivative orders a and/or b and in which with p 0,s,k and p 1,s,k being the matrix elements of the following two matrixes Here the relations P 0 = 2 i j p 0,s,k and P 1 = 2 i j p 1,s,k are held for s = 0, 1, 2, 3.
From the definition of g(x), it is known that its derivatives are only dependent on the weight function ϕ k (x). Therefore we are able to proceed the nth-order derivatives of g(x) by directly differentiating the weight function ϕ k (x) n times with respect to x.

Linearization
In the framework of the homotopy analysis method, the solution of Eq. (12) can be decomposed into The mth-order deformation equation is constructed, for m ≥ 1, as subject to the boundary conditions where c 0 is a homotopy convergence-control parameter, L is the linear operator, chosen as and The Coiflet wavelet projections are expressed by where P j denotes the Coiflet wavelet projection, h k (η) is the weight function defined by The mth-order deformation equation is then expressed via the Coiflet wavelet as Eq. (29) is solved by means of the Galerkin method. Multiplying both sides of Eq.(29) by h k (η) and integrating in the domain [0,1], we obtain wheref in which, a k,l and b k,l are defined as with the connection coefficients being given by

Validation
To ensure the accuracy of our proposed approach, the quadratic Riccati differential equation [30] is used as an illustrative example, which is written as subject to the initial conditions The above nonlinear system admits the exact solution In the framework of the homotopy analysis method, the initial guess and the linear operator are chosen as By means of our above-mentioned approach, under the definitions we are able to obtain the following Mth-order wavelet solutioñ To check its accuracy and efficiency, the averaged square error Err s,M and the relative error Err r,M are introduced, which are written by where t k = k/2 j ,Ỹ M andỸ E represent the Mth-order wavelet approximation and the exact solution respectively. Table 1 presents the averaged square error at different computational order M and resolution level j.
All results are obtained by using the same desktop computer, with its configuration being the processor Intel(R) Core(TM) i3-5015U CPU @ 2.10GHz, and installed memory(RAM) 4.00GB. It is seen from the table that both the computational order M and the resolution level j play key roles on error reduction.
Further to check the accuracy of the Coiflet wavelet-homotopy solutions, we compare the 120th-order  Table 2. This confirms the validity of our proposed method.

Results
In this section, the accuracy of our Coiflet solution is verified. The check procedure is divided into two steps. We first examine the residual error of the solution itself. It has known that Eq.(12) does not admit an analytical solution, we therefore introduce the relative error as its convergent indictor. In doing so, we define where f M (η k ) and f M +1 (η k ) are the Mth-order and the (M+1)th-order Coiflet wavelet-homotopy solutions.
For the sake of simplicity, we consider the particular case C w = 0. Other mathematical and physical parameters are chosen as c 0 = 0.1, Ha = 1000, Re = 50 and β = 5, respectively. As shown in Fig.2, the relative residual error decreases rapidly as the computational order evolves for all considered resolution levels. The maximum relative error is less than 1 × 10 −7 at 120th order computation.
We then compare our results with the published ones in literature. It is seen from Table 3 and Table   4 that our solution matches the published ones for both convergent and divergent channels. This further confirms the validity and accuracy of our proposed approach. Note that here the range of c 0 is chosen from the range [-0.5, -0.08] and Ha = 0.  In addition, we discuss the efficiency of our proposed technique. It has known that the traditional Coiflet wavelet approach [22] is equivalent to the first order Coiflet wavelet-homotopy computation at c 0 = −1. As a result we only need to compare our proposed method with the traditional optimized homotopy analysis method. The comparison is listed in Table 5 and Table 6, in which we find that our proposed approach is more than 100 times faster than the traditional homotopy analysis method when the values of f ′ (1) obtained by both techniques are same for four decimal places. Note that in these cases the convergence control parameter is chosen from the range [-0.5, -0.2] for both approaches.  Table 7 and Table 8. For small and medium Ha, all computational approaches give very excellent results. As Ha with varying integration distances and step sizes remain same, they are deemed as stable. The Newtonian iteration is executed with the convergent criteria of 1.0 × 10 −8 for all considered cases.
The physical aspects of the problem are discussed in the following part. The flow velocity behaviors under the influence of different physical parameters are portrayed Fig.3-Fig.8. In Fig.3 and Fig.4, the effects of the stretching/shrinking parameter C w on the flow velocity f (η) driven in convergent channel 50 and β = −5.
Clearly the case C w = 0 represents flow driven through a channel with stationary walls, which has been discussed by Esmaeilpour et al. [7]. And the case C w = 1 physically indicates that the stretching velocity of the channel wall is identical to the flow velocity. For the case C w > 1, the velocity of fluid increases faster than the centreline velocity u c . It means that the fluid particles intensifies very close to the channel wall owing to the influence of the viscosity and inertia of the fluid, just like the polymers concentrate close to the convergent channel rather than the centreline of the channel. For the case C w < 0, the wall stretching direction is opposite to the fluid flow direction, the flow is retarded owing to the resistance of the reverse force generated by the shrinking walls. It is found that in both cased the velocity profiles 50 and β = 5.
increase as C w enlarges. In addition, the effects of stretching/shringking wall parameters on the wall skin friction (wall shear stress) for both convergent and divergent channels are tabulated in the Table 9 and    solutions for extreme wide ranges of physical parameters. Besides, the physical mechanism of the flow problem has been presented and discussed. Some major physical findings of the current work are as follows: • For both convergent and divergent channel cases, wall stretching drags more fluid from the channel wall, while wall shrinking pulls the fluid to the wall and causes the back flow.
• The wall stretching of a divergent channel gives similar flow pattern as the wall shrinking of a convergent channel, or vice-versa.
• The Hartmann number Ha plays a key role on the flow velocity distribution. While opposite trends are obtained for stretching convergent channel case and shrinking divergent channel case.