Numerical Solutions of Heat Transfer for Magnetohydrodynamic Jeffery-Hamel Flow Using Spectral Homotopy Analysis Method

This paper studies heat transfer in a two-dimensional magnetohydrodynamic viscous incompressible flow in convergent/divergent channels. The temperature profile was obtained numerically for both cases of convergent/divergent channels. It was found that the temperature profile increases with an increase in Reynold number, Prandtl number, Nusselt number and angle of the wall but decreases with an increase in Hartmann number. A relatively new numerical method called the spectral homotopy analysis method (SHAM) was used to solve the governing non-linear differential equations. The SHAM 3rd order results matched with the DTM and shooting, showing that SHAM is feasible as a technique to be used.


Introduction
In recent years, the study of fluid flow and associated heat transfer phenomena has received increased attention from researchers and scientists due to many applications in industry and related areas. Viscous incompressible flow through a converging/diverging channel is commonly known as Jeffery-Hamel flow [1,2].This is an essential type of flow in the field of fluid mechanics. Applications of these types of flows include flow through rivers, different engineering processes and in the field of biology [2]. After the pioneering work of Jeffery and Hamel, many researchers have made contributions regarding the applications of this type of flow [3][4][5][6][7][8]. Magnetohydrodynamics (MHD) is the study of the interaction between a magnetic field and a moving conducting fluid. Applications of MHD flow include MHD generators, liquid metals, pumps, flow meters and several other engineering works [9][10][11]. Some of these applications involve Jeffery Hamel flows and thus motivate our study.
Motsa et al. [12] used the spectral homotopy analysis method (SHAM) to solve the MHD Jeffery Hamel problem, but their work did not take into consideration heat transfer. In this paper, we study heat transfer effects associated with MHD Jeffery-Hamel flow. A local similarity solution is used to convert the governing partial differential equations into ordinary differential equations before being solved using SHAM. The effects of various parameters on the temperature profile are discussed. Thus, the contributions of this paper endeavours to make are: (1) to reveal the behaviour of heat transfer and fluid flow towards both convergent and divergent channels; (2) to simulate the behaviour of the pertinent parameters such as Reynolds number, Eckert number, Prandtl number, Hartmann number, convergent/divergent angles that affect the fluid flow and temperature distributions, and to critically examine the graphical analysis of the effective parameters on the temperature distribution; and (3) to use SHAM incorporated with similarity transformations to synthesise chemical systems, understanding the behavior of the heat transfer and transport process. An accurate solution from the proposed method could be a stepping stone to establishing mathematical formulations to describe various microfluidic devices.

Mathematical Model
Let us consider incompressible steady two-dimensional viscous fluid flow between non-parallel plates making angle 2α, as shown in Figure 1. The plates are deemed to be divergent if α > 0 and convergent if α < 0. Suppose that the flow is along the radial direction and is a function of r and θ only, that is v = (u (r θ,),0). Heat transfer effects on the flow are considered with T as temperature assumed to be a function of r and θ. The continuity equation and the Navier-Stokes equations, when written in polar coordinates, are of the form [12] ρ r ∂ ∂r (ru(r, θ)) = 0, and the energy equation is subject to the boundary conditions [27] where ρ is the density, v be the coefficient of kinematic viscosity, B 0 be the electromagnetic induction, σ be the conductivity of the fluid, k is the thermal conductivity, P be the pressure and c p is the specific heat at constant pressure, U max is the maximum velocity at the center of the channel (r = 0). Equations (1)-(4) constitute the governing partial differential equations for the problem under consideration and (5) are the associated boundary condition. Physically, when u θ = 0, the expression of the velocity in radial flow is given as [27] u(r, where f (θ) is an arbitrary function of θ. Mathematically, Equation (6) can be obtained from Equation (1) by integrating both sides with respect to r.    After using Equation (6), with boundary conditions Now Equations (7)-(9) are converted into non-dimensional form after using the following dimensionless variables [7], where T w is the wall temperature, F(y) is the dimensionless velocity parameter which can be obtained by dividing f (θ) to its maximum value f max = r U max . G(y) is the dimensionless temperature. Substituting Equation (10) into Equations (7)-(9) and after simplification, we have finally subject to boundary conditions where Re = αU max v and H 2 = σB 0 2 ρv be the Reynolds number and square of the Hartmann number respectively. Also, the Eckert number Ec and Prandtl numbers Pr are defined as [32]

Solution by Using Spectral Homotopy Analysis Method
In the next section, the nonlinear boundary value problem describes in Equations (11)-(14) will be solved by using the spectral homotopy analysis method [12]. In order to apply the spectral homotopy analysis method to the problem, we suppose the initial guess satisfying the boundary conditions Equations (13) and (14) for the above problem is The problem domain is first transformed from [0, 1] to [−1, 1] using the following mapping from y to x as To convert non-homogeneous boundary conditions into homogeneous form, define the following transformations After substituting Equation (17), in Equations (11)- (14) the boundary value problem reduces to the form subject to boundary conditions where primes denote the derivatives with respect to x and For initial approximation, solution of the non-homogeneous linear part of the governing Equation (18) is obtained by using the Chebyshev pseudospectral method, that is subject to boundary conditions Here U 0 (x) is approximated as a truncated series of Chebyshev polynomial as where T k is the k th Chebyshev polynomial, U k are coefficients and x 0 , x 1 , . . . , x N are Gauss-Lobatto collocation points [5] defined as and s order derivative of U 0 (x) at this collocation points are Here D is the Chebyshev spectral differentiation matrix [5]. Substituting Equations (25)-(27) into Equations (23) and (24), given subject to boundary conditions Here where T represents the transpose and diag is a diagonal matrix of order (N + 1) × (N + 1). Now, we utilise the boundary conditions (29), by deleting the first and last rows and columns of matrix A and deleting the first and last rows of U 0 and ϕ. On the resulting last row of modified matrix A, the boundary condition (30) is imposed and replacing the resulting last row of the modified matrix ϕ by zero. Hence [U 0 (x 1 ), U 0 (x 2 ), . . . , U 0 (x N−1 )], it can be calculated by which is considered as an initial approximation for the solution of Equation (18), using the spectral homotopy analysis method. Now, define a linear operator as where q is the embedding parameter, which is defined by the homotopy analysis method and U(x; q) is an unknown function. So, the zeroth-order deformation equation is given by Makukula et al. [33] ( where h is convergence controlling parameter and N[U(x; q)] is a nonlinear operator defined as The m th order deformation equations are subject to the boundary conditions where By applying the Chebyshev pseudospectral transformation [5] on Equations (36)-(38), we have and subject to the boundary conditions where Now for boundary conditions (41), the first and last rows of P m−1 and ϕ, and the first and last rows and first and last columns of A in (39) are deleted. Also, to incorporate the boundary condition (42), the first row of Ψ and the first row and column B are deleted. Boundary condition (44) is imposed on the last row of the modified matrix B, and the last row of the modified matrix Ψ is replaced by zero. Therefore, the recursive formulas for m ≥ 1, is given by, Hence by using the initial approximation, obtained by Equation (32), we can find the higher-order approximation U m (x) and V m (x) for m ≥ 1, using recursive formula (Equation (49)) and (Equation (50)). Thus, once U m is obtained, we can calculate V m from Equation (50). The techniques used in this section are based on the techniques of spectral homotopy analysis method as described in Motsa et al. [12] and Motsa et al. [34].

Results and Discussion
This variation of temperature profile with viscous dissipation under the effects of various parameters like; Eckert number Ec, Reynold number Re, angle α and Prandtl number Pr are analyzed. The spectral homotopy analysis method is used to obtain the results for temperature profile G(y) with different values of parameters as shown graphically in Figures 2-11. Figures 2-6 shows the effects of parameters on temperature profile by considering diverging channel flow. As shown in Figure 2, which maps the effects of Reynold number Re on temperature profile, it is evident that temperature increases with the increase in Reynold number. There is a smooth increase in the temperature profile with varying Prandtl number Pr and Eckert number Ec, as in Figures 3 and 4. Figure 5 shows that temperature profile increases when the angle of wall increases. Now, the case is quite the opposite for Hartmann number in the diverging channel, that is, with the increase in Hartmann number H, temperature profile decreases. This is clear from Figure 6. Figures 7-11 show the effects of this parameter on temperature profile when converging channel flow is under consideration. It is worth noting in Figure 7, with the increase in Reynold number Re, temperature profile decreases but does not show similar behaviour to the case of diverging channel flow. However, Figure 8, shows the behaviour of Prandtl number Pr on temperature profile is very similar to the case of diverging channel flow. Temperature profile also increases in the case of increasing Eckert number Ec shown in Figure 9, as in the case of the diverging channel flow profile. In the case of converging channel flow temperature profile is quite opposite to that of the converging channel with the increasing magnitude of wall angle α depicted in Figure 10. Also, Figure 11 shows the effects of the Hartmann number H on the temperature profile. With the increment of the Hartmann number, there is a slight decrease in the temperature profile. Tables 1 and 2 show the effects of magnetic field on skin friction and Nusselt number taken for different values of Hartmann number H. Numerical values obtained by SHAM are also calculated using an approximate analytical technique known as the differential transform method (DTM) and a traditional numerical technique, known as the shooting method.

Conclusions
In this study, heat transfer effects were analyzed in MHD converging/diverging channel flows. These effects are shown graphically in Figures (2)- (11) and numerically in Tables 1 and 2. A newly developed spectral homotopy analysis method (SHAM) by Motsa et al. [12] was used to calculate the results. The numerical values for skin friction and Nusselt number for varying Hartmann number are depicted in Tables 1 and 2. They were also calculated using the DTM and the shooting method. SHAM third-order results matched with the DTM and the shooting method, showing that SHAM is feasible as a technique to be used. The comparison of SHAM results up to third-order SHAM with both of the techniques showed excellent agreement up to six decimal places.

Conclusions
In this study, heat transfer effects were analyzed in MHD converging/diverging channel flows. These effects are shown graphically in Figures 2-11 and numerically in Tables 1 and 2. A newly developed spectral homotopy analysis method (SHAM) by Motsa et al. [12] was used to calculate the results. The numerical values for skin friction and Nusselt number for varying Hartmann number are depicted in Tables 1 and 2. They were also calculated using the DTM and the shooting method. SHAM third-order results matched with the DTM and the shooting method, showing that SHAM is feasible as a technique to be used.

Future Work
The present study has ignored the effect of nanofluid and porous medium, which are also relevant to various chemical dynamical systems. Future studies will aim to examine hybrid nanofluid over a thermoelastic porous medium [31,35,36].