The BEM Solutions of MHD Flow and Heat Transfer in a Rectangular Duct with Temperature Dependent Viscosity

Elif Ebren Kaya* and M¨unevver Tezer-Sezgin

Department of Mathematics, Middle East Technical University, Ankara, Turkey
E-mail: ebren@metu.edu.tr
^*Corresponding Author

Received 05 March 2019; Accepted 17 July 2019;
Publication 05 August 2019

Abstract

The steady, laminar, fully developed magnetohydrodynamic (MHD) flow of an incompressible, electrically conducting fluid with temperature dependent viscosity is studied in a rectangular duct together with its heat transfer. Although the induced magnetic field is neglected due to the small Reynolds number, the Hall effect, viscous and Joule dissipations are taken into consideration. The momentum and the energy equations are solved iteratively. Firstly, the momentum equation is solved by using the boundary element method with a parametrix since the diffusion term contains variable viscosity parameter depending on the temperature exponentially. Next, the momentum equation is also solved by using the dual reciprocity boundary element method (DRBEM) for comparison. The energy equation is solved by using the DRBEM keeping all the terms containing the velocity as inhomogeneity. The temperature and the velocity behaviours are examined for several values of Hartmann number, dimensionless viscosity parameter, Brinkmann number and the Hall parameter. As Hartmann number is increasing, the velocity magnitude drops which is a well known property of the MHD duct flow. Increasing viscosity parameter reduces both the flow and the temperature magnitudes whereas the increase in the Hall parameter accelerates the flow and increases the fluid temperature.

Keywords: BEM, DRBEM, MHD duct flow, variable viscosity

1 Introduction

The problems of MHD flow and heat transfer via rectangular ducts have very important applications in areas such as designing heat transfer, cooling of electronic systems, chemical reactors, nuclear reactor, combustion systems and MHD generators. MHD generators of rectangular cross-section are widespread applications of MHD flow. They are made of tubes having insulating walls which are cut with two electrodes placed parallel to each other and perpendicular to insulating walls. This configuration resembles the MHD flow in rectangular duct [5]. Therefore, there have been many theoretical and experimental studies for the MHD flow in rectangular ducts. Some researchers and their methods are the followings. Türk and Tezer-Sezgin [13] have studied natural convection flow in square enclosures under magnetic field using the finite element method (FEM). In their study, the momentum equations include the magnetic effect and the induced magnetic field, due to the small magnetic Reynolds number, is neglected. They showed that, the FEM with quadratic elements enables one to solve MHD natural convection flow under the effect of a magnetic field for large values of Rayleigh and Hartmann numbers. Alsoy-Akgün and Tezer-Sezgin [2] have studied natural convection MHD flow equations in cavities by using the dual reciprocity boundary element method and the differential quadrature method (DQM). They have compared these two methods in terms of the flow and temperature behaviors concluding that the DRBEM and the DQM give almost the same accuracy. Although these studies differ from each other in their numerical implementations, they have a common feature which is the constant viscosity in the momentum equations considered in all of them. However, the momentum equations can contain a variable viscosity depending on space and time variables and even on the temperature. Sayed-Ahmed and Attia [10] have solved MHD flow and heat transfer with variable viscosity for Newtonian fluids in a rectangular duct with the Hall effect by using the finite difference method (FDM). Their results are confirmed by an experimental study of laminar flow and heat transfer in 2:1 rectangular duct filled with mineral oil, conducted by Xie and Hartnett [14] since the viscosity of the mineral oil changes dramatically. They obtained good agreement between the experimental data and the analytical solution. Shin et al. [11] also investigated the influence of variable viscosity of temperature-dependent fluids on the laminar heat transfer and friction factor in a 2:1 rectangular duct. The governing mass, momentum and energy equations were solved using a finite volume method (FVM). They get the excellent agreement with the experimental results conducted by Xie and Hartnett [14]. Evcin et al. [6] have also studied the same variable viscosity MHD flow and heat transfer problem by using mixed finite element method (FEM) together with optimal control techniques in order to control the system in desired velocity and temperature by the help of physically significant parameters of the system as control variables. Moreover, Prasad et al. [9] have presented FDM solution for the non-linear equations of MHD flow and heat transfer of an electrically conducting fluid over a stretching sheet with variable thickness. Cherlacola et al. [4] have examined the magneto hydrodynamic boundary layer flow with heat and mass transfer of Williamson nanofluid over a stretching sheet with variable thickness and variable thermal conductivity under the radiation effect. They have solved the differential equations subjected to the appropriate boundary conditions by using the spectral quasi-linearisation method (SQLM) and the sheet is non-flat in their study. Türk [12] also has proposed a Chebyshev spectral collocation method (CSCM) approximation to solve the thermally coupled MHD equations. In this study, the flow is subjected to a vertically applied magnetic field, and the presence of the induced magnetic field is also taken into consideration.

There are quite a number of fluid dynamics problems in which the fluid viscosity is varying with space, time and even with temperature. Among the viscosity depending on the temperature, exponential dependence is the most suitable according to the experimental results which has been indicated in [7]. Also, Capone and Gentile [3] have considered the exponentially varying temperature dependent viscosity and they have constructed the nonlinear stability analysis of convection for fluids with stress-free boundary conditions by using the Lyapunov direct method. In our case, the diffusion operator contains the viscosity function in its outer derivative terms. Thus, if the weighted residual statement is required for an equation containing a diffusion operator with a variable function, the corresponding fundamental solution or a weight function is required. For this purpose, Al-Jawary and Wrobel [1] derived such a weight function for variable coefficient diffusion equation. They have also considered the general partial differential equation with variable coefficients and they have used a parametrix (Levi function) which is usually available to solve the problem containing variable coefficients directly and accurately. In their study, with the help of the parametrix, the differential equation is reduced to a boundary-domain integral or integro-differential equation (BDIE or BDIDE).

In this paper, the steady, laminar, fully developed MHD flow of an incompressible, electrically conducting fluid with temperature dependent viscosity is solved in a rectangular duct together with its heat transfer by using both the boundary element method (BEM) with a parametrix (parametrix BEM) and the DRBEM with the fundamental solution of Laplace’s equation. Although, the induced magnetic field is neglected by taking small magnetic Reynolds number, the Hall effect, viscous and Joule dissipations are taken into consideration. The BEM procedure given in this paper uses the fundamental solution (parametrix) which treats the variable viscosity in the diffusion term directly [1]. That is, it is capable of addressing to the diffusion operator of the equations in its original form. Then, the momentum equation is also solved by using the dual reciprocity boundary element method using the fundamental solution of Laplace’s equation as is done in the energy equation by taking all the terms other than Laplacian as inhomogeneity. The velocity and the temperature profiles are obtained for several values of Hartmann number ($Ha$), Brinkmann number ($Br$), viscosity parameter ($B$) and the Hall parameter ($m$). The BEM implementations (both parametrix BEM and the DRBEM) discretizing only the boundary with constant elements capture the well known behavior of the MHD flow and the temperature of the fluid. Therefore, their computational costs are very small comparing with the other numerical methods, especially the DRBEM requires less CPU time since it involves no domain integral computation.

2 Mathematical Formulation of the Problem

The laminar, steady flow of a viscous, incompressible, electrically conducting fluid is considered in a long pipe with the pipe-axis velocity $w=w\left(x,y\right)$ and the temperature $T=T\left(x,y\right)$ varying only in the cross-section (duct) $\Omega =\left[0,a\right]ü\left[0,b\right]$ of the pipe. The governing non-dimensional momentum and energy equations with viscosity coefficient $\overline{\mu }=e^{-BT}$ are given in the domain (duct) $\mathrm{\Omega <!--\; \Omega \; -->}$ as

$$\begin{array}{rl}& \frac{\mathrm{\partial <!--\; \partial \; -->}}{\mathrm{\partial <!--\; \partial \; -->}x}\left(\stackrel{¯<!--\; ¯\; -->}{\mathrm{\mu <!--\; \mu \; -->}}\frac{\mathrm{\partial <!--\; \partial \; -->}w}{\mathrm{\partial <!--\; \partial \; -->}x}\right)+\frac{\mathrm{\partial <!--\; \partial \; -->}}{\mathrm{\partial <!--\; \partial \; -->}y}\left(\stackrel{¯<!--\; ¯\; -->}{\mathrm{\mu <!--\; \mu \; -->}}\frac{\mathrm{\partial <!--\; \partial \; -->}w}{\mathrm{\partial <!--\; \partial \; -->}y}\right)=-<!--\; -\; -->1+\frac{H{a}^2}{1+m^2}w,(1)\\ & {\mathrm{\nabla <!--\; \nabla \; -->}}^{2}T+Br\stackrel{¯<!--\; ¯\; -->}{\mathrm{\mu <!--\; \mu \; -->}}[(\frac{\mathrm{\partial <!--\; \partial \; -->}w}{\mathrm{\partial <!--\; \partial \; -->}x}{)}^2+(\frac{\mathrm{\partial <!--\; \partial \; -->}w}{\mathrm{\partial <!--\; \partial \; -->}y}{)}^2]+\frac{H{a}^{2}Br}{1+m^2}{w}^2={\displaystyle \frac{w}{w_m}},(2)\end{array}$$

where viscosity parameter, $B$, Hall parameter, $m$, Hartmann number, $Ha$ and Brinkmann number, $Br$ are given as

Figure 1: Physical configuration of the problem.

$$\begin{array}{rl}B& =(bl\mathrm{\rho <!--\; \rho \; -->}{c}_p{w}_m{a}^2(d{T}_m/dz))/k,m=\mathrm{\sigma <!--\; \sigma \; -->}\mathrm{\beta <!--\; \beta \; -->}{B}_0,w_m={\int <!--\; \int \; -->}_0^b{\int <!--\; \int \; -->}_0^{a}wdxdy,(3)\\ Ha& =\frac{a{B}_0\sqrt{\mathrm{\sigma <!--\; \sigma \; -->}}}{\sqrt{{\mathrm{\mu <!--\; \mu \; -->}}_0}},Br=\frac{k{\mathrm{\mu <!--\; \mu \; -->}}_0^2(-<!--\; -\; -->dp/dz)}{\mathrm{\rho <!--\; \rho \; -->}{c}_p{w}_m(d{T}_m/dz)}.\end{array}$$

Here, $B_0$, $\sigma$ and $\beta$ are the applied magnetic field intensity, electric conductivity of the fluid and the Hall factor, respectively. $\rho$, $c_p$, $k$, ${\mu }_0$ and $l$, $a$, $b$ are the density, the specific heat capacity, the thermal conductivity, the viscosity at $T=T_m$ of the fluid and the characteristic lengths in the z- and x-, y-directions, respectively. $\frac{dp}{dz}$ and $\frac{d{T}_m}{dz}$ represent uniform pressure and temperature gradients in the pipe-axis direction, respectively. As is clearly seen from the definition of $\overline{\mu }$ above, the relation between the temperature and the viscosity is stated such that the increase of the temperature results in exponential decrease in the viscosity coefficient. This is due to the attractive interactions of the fluid molecules. Turning back to our concern, the boundary conditions $w=0$ and $T=0$ are imposed on the walls $\Gamma$ of the duct indicating no-slip velocity and the cold wall conditions.

3 The BEM Implementations

3.1 The Parametrix BEM

The velocity Equation (1) is solved by using the BEM with the fundamental solution of diffusion equation containing a variable coefficient which is called as a parametrix $P\left(x,x_i;y,y_i\right)=\frac{ln|r-r_i|}{2\pi \overline{\mu }\left(x_i,y_i\right)}$ given in [1], where $r=\left(x,y\right)$ and $r_i=\left(x_i,y_i\right)$ are the variable and fixed points in $\Gamma \cup \Omega$, respectively, $\Gamma$ defines the boundary of the duct $\Omega$. The boundary of the duct is discretized by $N$ number of constant elements. Multiplying the velocity Equation (1) by the parametrix $P\left(x,x_i;y,y_i\right)$ and integrating over the domain, the equation is obtained as

$$\begin{array}{rcll}& & {\int }_{\Omega }\left(\right.\frac{\partial }{\partial x}\left(\overline{\mu }\frac{\partial w}{\partial x}\right)+\frac{\partial }{\partial y}\left(\overline{\mu }\frac{\partial w}{\partial y}\right)\left)\right.P\left(x,x_i;y,y_i\right)d\Omega & \\ & & ={\int }_{\Omega }\left(\right.-1+\frac{H{a}^2}{1+m^2}w\left)\right.P\left(x,x_i;y,y_i\right)d\Omega ,i=1,\dots ,N& \text{(4)}\end{array}$$

where $i$ is the source point.

Applying Green’s second identity two times, the following discretized boundary integral equation is obtained

$$\begin{array}{rcll}& & c_i{w}_i-\frac{1}{2\pi }\sum _{j=1}^N{\int }_{{\Gamma }_j}w\left(r-r_i\right)\overrightarrow{n}|r-r_i{|}^2\overline{\mu }\left(x,y\right)\overline{\mu }\left(x_i,y_i\right)d{\Gamma }_j& \\ & & +\frac{1}{2\pi }\sum _{j=1}^N{\int }_{{\Gamma }_j}\overline{\mu }\left(x,y\right)\partial w\partial nln|r-r_i|1\overline{\mu }\left(x_i,y_i\right)d{\Gamma }_j& \\ & & +\frac{1}{2\pi \overline{\mu }\left(x_i,y_i\right)}\sum _{j=1}^{{\left(N∕4\right)}^2}{\int }_{{\Omega }_j}w\frac{\left(r-r_i\right)\vec{n}}{|r-r_i{|}^2}\partial \overline{\mu }\left(x,y\right)\partial nd{\Omega }_j& \\ & & =\frac{1}{2\pi \overline{\mu }\left(x_i,y_i\right)}\sum _{j=1}^{{\left(N∕4\right)}^2}{\int }_{{\Omega }_j}ln|r-r_i|\left(-1+\frac{H{a}^2}{1+m^2}w\right)d{\Omega }_j,& \\ & & i=1,\dots ,N& \text{(5)}\end{array}$$

where $i$ denotes the source point, ${}_j$ and ${}_j$ are constant boundary element and domain cell, respectively. The vector $\overrightarrow{n}$ is the outward normal vector to $\Gamma$.

The Equation (5) is rewritten in matrix-vector form,

$$\begin{array}{r}\mathbf{\text{H}}w-<!--\; -\; -->\mathbf{\text{G}}q=I^1(w)-<!--\; -\; -->I^2(w),\end{array}$$

(6)

where the vector $q=\overline{\mu }\left(x,y\right)\frac{\partial w}{\partial n}$, and the entries of the matrices H and G are given as

$$H_{ij}=-<!--\; -\; -->{\displaystyle \frac{1}{2}}{\mathrm{\delta <!--\; \delta \; -->}}_{ij}+{\displaystyle \frac{1}{2\mathrm{\pi <!--\; \pi \; -->}}}{\displaystyle {\int <!--\; \int \; -->}_{{\mathrm{\Gamma <!--\; \Gamma \; -->}}_j}{\displaystyle \frac{(r-<!--\; -\; -->r_i)\stackrel{\to <!--\; \to \; -->}{n}}{|r-<!--\; -\; -->r_i{|}^2}}{\displaystyle \frac{\stackrel{¯<!--\; ¯\; -->}{\mathrm{\mu <!--\; \mu \; -->}}(x,y)}{\stackrel{¯<!--\; ¯\; -->}{\mathrm{\mu <!--\; \mu \; -->}}(x_i,y_i)}}d{\mathrm{\Gamma <!--\; \Gamma \; -->}}_j,H_{ii}=-<!--\; -\; -->\underset{j=1,j\ne <!--\; \ne \; -->i}{\overset{N}{\sum <!--\; \sum \; -->}}{H}_{ij},}$$

$$G_{ij}=12\pi {\int }_{{\Gamma }_j}ln|r-r_i|1\overline{\mu }\left(x_i,y_i\right)d{\Gamma }_j,G_{ii}=1\overline{\mu }\left(x_i,y_i\right)\frac{\xi }{2\pi }\left(ln\left(\frac{\xi }{2}\right)-1\right),$$

(7)

where $\xi$ represents the length of the constant boundary element.

The vectors $I^1\left(w\right)$ and $I^2\left(w\right)$ have the entries as

$$\begin{array}{rl}{I}_i^1(w)& =\frac{1}{2\mathrm{\pi <!--\; \pi \; -->}\stackrel{¯<!--\; ¯\; -->}{\mathrm{\mu <!--\; \mu \; -->}}(x_i,y_i)}\underset{j=1}{\overset{(N/4{)}^2}{\sum <!--\; \sum \; -->}}{\int <!--\; \int \; -->}_{{\mathrm{\Omega <!--\; \Omega \; -->}}_j}{w}_i\frac{(r-<!--\; -\; -->r_i)\stackrel{\to <!--\; \to \; -->}{n}}{|r-<!--\; -\; -->r_i{|}^2}{\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}\stackrel{¯<!--\; ¯\; -->}{\mathrm{\mu <!--\; \mu \; -->}}(x,y)}{\mathrm{\partial <!--\; \partial \; -->}n}}d{\mathrm{\Omega <!--\; \Omega \; -->}}_j(8)\\ I_i^2(w)& =\frac{1}{2\mathrm{\pi <!--\; \pi \; -->}\stackrel{¯<!--\; ¯\; -->}{\mathrm{\mu <!--\; \mu \; -->}}(x_i,y_i)}\underset{j=1}{\overset{(N/4{)}^2}{\sum <!--\; \sum \; -->}}{\int <!--\; \int \; -->}_{{\mathrm{\Omega <!--\; \Omega \; -->}}_j}ln|r-<!--\; -\; -->r_i|(-<!--\; -\; -->1+\frac{H{a}^2}{1+m^2}{w}_i)d{\mathrm{\Omega <!--\; \Omega \; -->}}_j,\end{array}$$

with $i=1,\dots ,N$ and $j=1,\dots ,{\left(N∕4\right)}^2$.

The unknown normal derivative vector $q$ on the boundary can be solved from the system (6) since $w=0$ everywhere on the boundary nodes by Gaussian elimination. Then, these boundary values are used in the Equation (5) by taking $c_i=1$ to compute the velocity values $w$ at each interior point $i$,

$$w_i=\sum _{j=1}^N{\overline{H}}_{ij}{w}_j-\sum _{j=1}^N{G}_{ij}{q}_j-\sum _{j=1}^{{\left(N∕4\right)}^2}{I}_i^1\left(w\right)+\sum _{j=1}^{{\left(N∕4\right)}^2}{I}_i^2\left(w\right),i=1,\dots ,L$$

(9)

where $L$ is the number of interior points and

$${\overline{H}}_{ij}=12\pi {\int }_{{\Gamma }_j}\left(r-r_i\right)\overrightarrow{n}|r-r_i{|}^2\overline{\mu }\left(x,y\right)\overline{\mu }\left(x_i,y_i\right)d{\Gamma }_j,i=1,\dots ,L,j=1,\dots ,N$$

3.2 The DRBEM Formulation

The velocity Equation (1) is rewritten as

$${\mathrm{\nabla <!--\; \nabla \; -->}}^{2}w={\displaystyle \frac{1}{\stackrel{¯<!--\; ¯\; -->}{\mathrm{\mu <!--\; \mu \; -->}}}}(-<!--\; -\; -->1+{\displaystyle \frac{H{a}^2}{1+m^2}}w-<!--\; -\; -->{\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}\stackrel{¯<!--\; ¯\; -->}{\mathrm{\mu <!--\; \mu \; -->}}}{\mathrm{\partial <!--\; \partial \; -->}x}}{\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}w}{\mathrm{\partial <!--\; \partial \; -->}x}}-<!--\; -\; -->{\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}\stackrel{¯<!--\; ¯\; -->}{\mathrm{\mu <!--\; \mu \; -->}}}{\mathrm{\partial <!--\; \partial \; -->}y}}{\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}w}{\mathrm{\partial <!--\; \partial \; -->}y}}),$$

(10)

to be able to use the fundamental solution of Laplace’s equation, $u^{*}=\frac{ln|r-r_i|}{2\pi }$ given in [8]. Multiplying by $u^{*}$ both sides of the equation and applying Green’s second identity two times the following boundary-domain integral equation is obtained

$$\begin{array}{rl}{c}_i{w}_i& -<!--\; -\; -->{\int <!--\; \int \; -->}_{\mathrm{\Gamma <!--\; \Gamma \; -->}}{q}^{\ast <!--\; \ast \; -->}wd\mathrm{\Gamma <!--\; \Gamma \; -->}+{\int <!--\; \int \; -->}_{\mathrm{\Gamma <!--\; \Gamma \; -->}}{u}^{\ast <!--\; \ast \; -->}\frac{\mathrm{\partial <!--\; \partial \; -->}w}{\mathrm{\partial <!--\; \partial \; -->}n}d\mathrm{\Gamma <!--\; \Gamma \; -->}(11)\\ & ={\int <!--\; \int \; -->}_{\mathrm{\Omega <!--\; \Omega \; -->}}{\displaystyle \frac{1}{\stackrel{¯<!--\; ¯\; -->}{\mathrm{\mu <!--\; \mu \; -->}}}}(-<!--\; -\; -->1+{\displaystyle \frac{H{a}^2}{1+m^2}}w-<!--\; -\; -->{\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}\stackrel{¯<!--\; ¯\; -->}{\mathrm{\mu <!--\; \mu \; -->}}}{\mathrm{\partial <!--\; \partial \; -->}x}}{\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}w}{\mathrm{\partial <!--\; \partial \; -->}x}}-<!--\; -\; -->{\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}\stackrel{¯<!--\; ¯\; -->}{\mathrm{\mu <!--\; \mu \; -->}}}{\mathrm{\partial <!--\; \partial \; -->}y}}{\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}w}{\mathrm{\partial <!--\; \partial \; -->}y}})u^{\ast <!--\; \ast \; -->}d\mathrm{\Omega <!--\; \Omega \; -->},\end{array}$$

where $q^{*}=\frac{\partial u^{*}}{\partial n}$, $i=1,\dots ,N$.

The integrand on the right hand side of the Equation eqvel is considered as inhomogeneity $b_1\left(x,y\right)$. This inhomogeneity term can be expanded by a series of a radial basis functions as

$$b^1(x,y)={\displaystyle \frac{1}{\stackrel{¯<!--\; ¯\; -->}{\mathrm{\mu <!--\; \mu \; -->}}}}(-<!--\; -\; -->1+{\displaystyle \frac{H{a}^2}{1+m^2}}w-<!--\; -\; -->{\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}\stackrel{¯<!--\; ¯\; -->}{\mathrm{\mu <!--\; \mu \; -->}}}{\mathrm{\partial <!--\; \partial \; -->}x}}{\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}w}{\mathrm{\partial <!--\; \partial \; -->}x}}-<!--\; -\; -->{\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}\stackrel{¯<!--\; ¯\; -->}{\mathrm{\mu <!--\; \mu \; -->}}}{\mathrm{\partial <!--\; \partial \; -->}y}}{\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}w}{\mathrm{\partial <!--\; \partial \; -->}y}})=\underset{j=1}{\overset{N+L}{\sum <!--\; \sum \; -->}}{\mathrm{\alpha <!--\; \alpha \; -->}}_j{f}_j(x,y),$$

(12)

where $L$ denotes the number of interior points, $f_j\left(x,y\right)$ are the radial basis functions which are connected to the particular solutions ${û}_j$ with ${\nabla }^2{û}_j=f_j$, the coefficients ${\alpha }_j$ are undetermined constants.

Substituting $f_j={\nabla }^2{û}_j$ in b1 and applying Green’s second identity two times again for the right hand side of the Equation eqvel, the following boundary integral equation is obtained.

$$c_i{w}_i-{\int }_{\Gamma }{q}^{*}wd\Gamma +{\int }_{\Gamma }{u}^{*}\frac{\partial w}{\partial n}d\Gamma =\sum _{j=1}^{N+L}{\alpha }_j\left(\right.c_i{û}_{ij}-{\int }_{\Gamma }{q}^{*}{û}_{j}d\Gamma +{\int }_{\Gamma }{u}^{*}{\widehat{q}}_{j}d\Gamma \left)\right.,$$

(13)

where ${\widehat{q}}_j=\frac{\partial {û}_j}{\partial n}=\frac{\partial {û}_j}{\partial x}\frac{\partial x}{\partial n}+\frac{\partial {û}_j}{\partial y}\frac{\partial y}{\partial n}$, $i=1,\dots ,N$.

Constructing the matrices $Û$, $\widehat{Q}$ and coordinate matrix $F$ by taking the vectors $\widehat{u_j}$, $\widehat{q_j}$ and $f_{ij}=1+r_{ij}$, $r_{ij}$ being the distance from the point $i$ to the point $j$, as columns respectively, and evaluating the values of $b^1$ at $N+L$ points, a set of linear equations as $b^1=F\alpha$ is obtained. The resulting $\mathbf{\text{H}}$ and $\mathbf{\text{G}}$ matrices are extended to $\left(N+L\right)\times \left(N+L\right)$ matrices by adding zero and identity matrices [8]. Thus, the Equation (13) can be written now as a $\left(N+L\right)+\left(N+L\right)$ system.

$$\mathbf{\text{H}}w-<!--\; -\; -->\mathbf{\text{G}}{\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}w}{\mathrm{\partial <!--\; \partial \; -->}n}}=(\mathbf{\text{H}}\stackrel{^<!--\; ^\; -->}{\mathbf{\text{U}}}-<!--\; -\; -->\mathbf{\text{G}}\stackrel{^<!--\; ^\; -->}{\mathbf{\text{Q}}}){\mathbf{\text{F}}}^{-<!--\; -\; -->1}{b}^1$$

(14)

where the vector $b^1$ has entries as

$$\begin{array}{cc}{b}_i^1& =\{{\displaystyle \frac{1}{\stackrel{¯<!--\; ¯\; -->}{\mathrm{\mu <!--\; \mu \; -->}}}}[-<!--\; -\; -->1+{\displaystyle \frac{H{a}^2}{1+m^2}}{w}_i-<!--\; -\; -->({\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}\stackrel{¯<!--\; ¯\; -->}{\mathrm{\mu <!--\; \mu \; -->}}}{\mathrm{\partial <!--\; \partial \; -->}x}}{)}_i({\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}w}{\mathrm{\partial <!--\; \partial \; -->}x}}{)}_i-<!--\; -\; -->({\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}\stackrel{¯<!--\; ¯\; -->}{\mathrm{\mu <!--\; \mu \; -->}}}{\mathrm{\partial <!--\; \partial \; -->}y}}{)}_i({\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}w}{\mathrm{\partial <!--\; \partial \; -->}y}}{)}_i]\},\\ & i=1,\dots <!--\; \dots \; -->,N+L.\end{array}$$

(15)

and the vectors in this equation are multiplied entrywise.

The space derivatives of the velocity $w$ and viscosity coefficient $\overline{\mu }$ in the vector $b^1$ are computed by using the coordinate matrix $\mathbf{\text{F}}$ as

$$\begin{array}{rl}{\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}w}{\mathrm{\partial <!--\; \partial \; -->}x}}& ={\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}\mathbf{\text{F}}}{\mathrm{\partial <!--\; \partial \; -->}x}}{\mathbf{\text{F}}}^{-<!--\; -\; -->1}w,{\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}w}{\mathrm{\partial <!--\; \partial \; -->}y}}={\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}\mathbf{\text{F}}}{\mathrm{\partial <!--\; \partial \; -->}y}}{\mathbf{\text{F}}}^{-<!--\; -\; -->1}w,{\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}\stackrel{¯<!--\; ¯\; -->}{\mathrm{\mu <!--\; \mu \; -->}}}{\mathrm{\partial <!--\; \partial \; -->}x}}={\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}\mathbf{\text{F}}}{\mathrm{\partial <!--\; \partial \; -->}x}}{\mathbf{\text{F}}}^{-<!--\; -\; -->1}\stackrel{¯<!--\; ¯\; -->}{\mathrm{\mu <!--\; \mu \; -->}}\text{and}\\ & {\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}\stackrel{¯<!--\; ¯\; -->}{\mathrm{\mu <!--\; \mu \; -->}}}{\mathrm{\partial <!--\; \partial \; -->}y}}={\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}\mathbf{\text{F}}}{\mathrm{\partial <!--\; \partial \; -->}y}}{\mathbf{\text{F}}}^{-<!--\; -\; -->1}\stackrel{¯<!--\; ¯\; -->}{\mathrm{\mu <!--\; \mu \; -->}}.\end{array}$$

(16)

Then, the Equation (14) becomes

$$\begin{array}{rl}\mathbf{\text{H}}w& -<!--\; -\; -->\mathbf{\text{G}}{\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}w}{\mathrm{\partial <!--\; \partial \; -->}n}}=(\mathbf{\text{H}}\stackrel{^<!--\; ^\; -->}{\mathbf{\text{U}}}-<!--\; -\; -->\mathbf{\text{G}}\stackrel{^<!--\; ^\; -->}{\mathbf{\text{Q}}}){\mathbf{\text{F}}}^{-<!--\; -\; -->1}\{{\displaystyle \frac{1}{\stackrel{¯<!--\; ¯\; -->}{\mathrm{\mu <!--\; \mu \; -->}}}}[-<!--\; -\; -->1+{\displaystyle \frac{H{a}^2}{1+m^2}}w-<!--\; -\; -->{\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}\mathbf{\text{F}}}{\mathrm{\partial <!--\; \partial \; -->}x}}{\mathbf{\text{F}}}^{-<!--\; -\; -->1}\stackrel{¯<!--\; ¯\; -->}{\mathrm{\mu <!--\; \mu \; -->}}{\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}\mathbf{\text{F}}}{\mathrm{\partial <!--\; \partial \; -->}x}}{\mathbf{\text{F}}}^{-<!--\; -\; -->1}w\\ & -<!--\; -\; -->{\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}\mathbf{\text{F}}}{\mathrm{\partial <!--\; \partial \; -->}y}}{\mathbf{\text{F}}}^{-<!--\; -\; -->1}\stackrel{¯<!--\; ¯\; -->}{\mathrm{\mu <!--\; \mu \; -->}}{\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}\mathbf{\text{F}}}{\mathrm{\partial <!--\; \partial \; -->}y}}{\mathbf{\text{F}}}^{-<!--\; -\; -->1}w]\}\end{array}$$

(17)

The components of $\mathbf{\text{H}}$ and $\mathbf{\text{G}}$ matrices for a constant element in DRBEM are the same as the matrices $\mathbf{\text{H}}$ and $\mathbf{\text{G}}$ given in the Equation (7) if $\overline{\mu }\left(x,y\right)=1$ is taken in parametrix BEM. Thus, the Equation (14) is rewritten in matrix-vector form

$$\begin{array}{r}\stackrel{¯<!--\; ¯\; -->}{\mathbf{\text{H}}}w-<!--\; -\; -->\mathbf{\text{G}}q=(\mathbf{\text{H}}\stackrel{^<!--\; ^\; -->}{\mathbf{\text{U}}}-<!--\; -\; -->\mathbf{\text{G}}\stackrel{^<!--\; ^\; -->}{\mathbf{\text{Q}}}){\mathbf{\text{F}}}^{-<!--\; -\; -->1}{\displaystyle \frac{1}{\stackrel{¯<!--\; ¯\; -->}{\mathrm{\mu <!--\; \mu \; -->}}}}(-<!--\; -\; -->1),\end{array}$$

(18)

where the entries of the matrix $\overline{\mathbf{\text{H}}}$ are given as

$$\begin{array}{rl}{\stackrel{¯<!--\; ¯\; -->}{H}}_{ij}& =H_{ij}-<!--\; -\; -->R_{ij},\text{where}\\ \mathbf{\text{R}}& =(\mathbf{\text{H}}\stackrel{^<!--\; ^\; -->}{\mathbf{\text{U}}}-<!--\; -\; -->\mathbf{\text{G}}\stackrel{^<!--\; ^\; -->}{\mathbf{\text{Q}}}){\mathbf{\text{F}}}^{-<!--\; -\; -->1}{\displaystyle \frac{1}{\stackrel{¯<!--\; ¯\; -->}{\mathrm{\mu <!--\; \mu \; -->}}}}({\displaystyle \frac{H{a}^2}{1+m^2}}-<!--\; -\; -->{\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}\mathbf{\text{F}}}{\mathrm{\partial <!--\; \partial \; -->}x}}{\mathbf{\text{F}}}^{-<!--\; -\; -->1}\stackrel{¯<!--\; ¯\; -->}{\mathrm{\mu <!--\; \mu \; -->}}{\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}\mathbf{\text{F}}}{\mathrm{\partial <!--\; \partial \; -->}x}}{\mathbf{\text{F}}}^{-<!--\; -\; -->1}-<!--\; -\; -->{\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}\mathbf{\text{F}}}{\mathrm{\partial <!--\; \partial \; -->}y}}{\mathbf{\text{F}}}^{-<!--\; -\; -->1}\stackrel{¯<!--\; ¯\; -->}{\mathrm{\mu <!--\; \mu \; -->}}{\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}\mathbf{\text{F}}}{\mathrm{\partial <!--\; \partial \; -->}y}}{\mathbf{\text{F}}}^{-<!--\; -\; -->1}),\\ {\stackrel{¯<!--\; ¯\; -->}{H}}_{ii}& =-<!--\; -\; -->\underset{j=1,j\ne <!--\; \ne \; -->i}{\overset{N}{\sum <!--\; \sum \; -->}}{\stackrel{¯<!--\; ¯\; -->}{H}}_{ij}.\end{array}$$

(19)

Then, the energy Equation (2) is solved by using the dual reciprocity boundary element method. Similar to the flow Equation laplace1, all the terms other than ${\nabla }^{2}T$ are treated as inhomogeneities. Multiplying by $u^{*}$ both sides of the equation and applying Green’s second identity two times on the left hand side of the equation, the following boundary-domain integral equation is obtained

$$\begin{array}{rl}{c}_i{T}_i& -<!--\; -\; -->{\int <!--\; \int \; -->}_{\mathrm{\Gamma <!--\; \Gamma \; -->}}{q}^{\ast <!--\; \ast \; -->}Td\mathrm{\Gamma <!--\; \Gamma \; -->}+{\int <!--\; \int \; -->}_{\mathrm{\Gamma <!--\; \Gamma \; -->}}{u}^{\ast <!--\; \ast \; -->}\frac{\mathrm{\partial <!--\; \partial \; -->}T}{\mathrm{\partial <!--\; \partial \; -->}n}d\mathrm{\Gamma <!--\; \Gamma \; -->}\\ & ={\int <!--\; \int \; -->}_{\mathrm{\Omega <!--\; \Omega \; -->}}(-<!--\; -\; -->Br\stackrel{¯<!--\; ¯\; -->}{\mathrm{\mu <!--\; \mu \; -->}}[(\frac{\mathrm{\partial <!--\; \partial \; -->}w}{\mathrm{\partial <!--\; \partial \; -->}x}{)}^2+(\frac{\mathrm{\partial <!--\; \partial \; -->}w}{\mathrm{\partial <!--\; \partial \; -->}y}{)}^2]-<!--\; -\; -->\frac{H{a}^{2}Br}{1+m^2}{w}^2+\frac{w}{w_m})u^{\ast <!--\; \ast \; -->}d\mathrm{\Omega <!--\; \Omega \; -->}\\ & ={\int <!--\; \int \; -->}_{\mathrm{\Omega <!--\; \Omega \; -->}}{b}^2(x,y)u^{\ast <!--\; \ast \; -->}d\mathrm{\Omega <!--\; \Omega \; -->},\end{array}$$

(20)

where

$$b^2(x,y)=-<!--\; -\; -->Br\stackrel{¯<!--\; ¯\; -->}{\mathrm{\mu <!--\; \mu \; -->}}[(\frac{\mathrm{\partial <!--\; \partial \; -->}w}{\mathrm{\partial <!--\; \partial \; -->}x}{)}^2+(\frac{\mathrm{\partial <!--\; \partial \; -->}w}{\mathrm{\partial <!--\; \partial \; -->}y}{)}^2]-<!--\; -\; -->\frac{H{a}^{2}Br}{1+m^2}{w}^2+\frac{w}{w_m}=\underset{j=1}{\overset{N+L}{\sum <!--\; \sum \; -->}}{\mathrm{\beta <!--\; \beta \; -->}}_j{f}_j(x,y)$$

(21)

is also approximated using the radial basis functions $f_j\left(x,y\right)$. ${\beta }_j$ are undetermined coefficients to be determined from the system $F\beta =b^2$ when $b^2$ is collocated at $N+L$ points.

Substituting $f_j={\nabla }^2{û}_j$ in b2 and therefore in eqtemp, and applying Green’s second identity two times again for the right hand side of the Equation eqtemp, the following boundary integral equation is obtained

$$\begin{array}{rl}{c}_i{T}_i& -<!--\; -\; -->{\int <!--\; \int \; -->}_{\mathrm{\Gamma <!--\; \Gamma \; -->}}{q}^{\ast <!--\; \ast \; -->}Td\mathrm{\Gamma <!--\; \Gamma \; -->}+{\int <!--\; \int \; -->}_{\mathrm{\Gamma <!--\; \Gamma \; -->}}{u}^{\ast <!--\; \ast \; -->}\frac{\mathrm{\partial <!--\; \partial \; -->}T}{\mathrm{\partial <!--\; \partial \; -->}n}d\mathrm{\Gamma <!--\; \Gamma \; -->}\\ & =\underset{j=1}{\overset{N+L}{\sum <!--\; \sum \; -->}}{\mathrm{\beta <!--\; \beta \; -->}}_j(c_i{\stackrel{^<!--\; ^\; -->}{u}}_{ij}-<!--\; -\; -->{\int <!--\; \int \; -->}_{\mathrm{\Gamma <!--\; \Gamma \; -->}}{q}^{\ast <!--\; \ast \; -->}{\stackrel{^<!--\; ^\; -->}{u}}_{j}d\mathrm{\Gamma <!--\; \Gamma \; -->}+{\int <!--\; \int \; -->}_{\mathrm{\Gamma <!--\; \Gamma \; -->}}{u}^{\ast <!--\; \ast \; -->}{\stackrel{^<!--\; ^\; -->}{q}}_{j}d\mathrm{\Gamma <!--\; \Gamma \; -->}).\end{array}$$

(22)

Thus, the equation (22) results in the matrix-vector system

$$HT-<!--\; -\; -->G{\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}T}{\mathrm{\partial <!--\; \partial \; -->}n}}=(H\stackrel{^<!--\; ^\; -->}{U}-<!--\; -\; -->G\stackrel{^<!--\; ^\; -->}{Q})F^{-<!--\; -\; -->1}{b}^2$$

(23)

where the vector $b^2$ has entries as

$$b_i^2=-Br\overline{\mu }\left[\right.\left(\right.\frac{\partial w}{\partial x}{\left)\right.}_i^2+\left(\right.\frac{\partial w}{\partial y}{\left)\right.}_i^2\left]\right.-\frac{H{a}^{2}Br}{1+m^2}{w}_i^2+\frac{w_i}{w_m},i=1,\dots ,N+L.$$

(24)

The squares of the vectors ${\left(\frac{\partial w}{\partial x}\right)}_i$, ${\left(\frac{\partial w}{\partial y}\right)}_i$ are computed by multiplying corresponding entries.

The squares of the first order derivatives of $w$ in $b^2$ vector in the equation b2 are also computed with the help of coordinate matrix $\mathbf{\text{F}}$ as

$$({\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}w}{\mathrm{\partial <!--\; \partial \; -->}x}}{)}^2=({\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}\mathbf{\text{F}}}{\mathrm{\partial <!--\; \partial \; -->}x}}{\mathbf{\text{F}}}^{-<!--\; -\; -->1}w{)}^2,({\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}w}{\mathrm{\partial <!--\; \partial \; -->}y}}{)}^2=({\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}\mathbf{\text{F}}}{\mathrm{\partial <!--\; \partial \; -->}y}}{\mathbf{\text{F}}}^{-<!--\; -\; -->1}w{)}^2,$$

(25)

Then, the Equation (23) becomes

$$\begin{array}{rl}\mathbf{\text{H}}T-<!--\; -\; -->\mathbf{\text{G}}{\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}T}{\mathrm{\partial <!--\; \partial \; -->}n}}& =(\mathbf{\text{H}}\stackrel{^<!--\; ^\; -->}{\mathbf{\text{U}}}-<!--\; -\; -->\mathbf{\text{G}}\stackrel{^<!--\; ^\; -->}{\mathbf{\text{Q}}}){\mathbf{\text{F}}}^{-<!--\; -\; -->1}(-<!--\; -\; -->Br\stackrel{¯<!--\; ¯\; -->}{\mathrm{\mu <!--\; \mu \; -->}}[({\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}\mathbf{\text{F}}}{\mathrm{\partial <!--\; \partial \; -->}x}}{\mathbf{\text{F}}}^{-<!--\; -\; -->1}w{)}^2+({\displaystyle \frac{\mathrm{\partial <!--\; \partial \; -->}\mathbf{\text{F}}}{\mathrm{\partial <!--\; \partial \; -->}y}}{\mathbf{\text{F}}}^{-<!--\; -\; -->1}w{)}^2]\\ & -<!--\; -\; -->\frac{H{a}^{2}Br}{1+m^2}{w}^2+\frac{w}{w_m}).\end{array}$$

(26)

Now, first the parametrix BEM discretized flow Equation (6) and the DRBEM discretized temperature Equation (26) are going to be solved iteratively giving the flow and the temperature of the fluid. Then, also the DRBEM discretized Equations (18) and (26) are going to be solved iteratively to determine the velocity and the temperature of the MHD flow to be able to compare with the parametrix BEM results.

4 Numerical Results and Discussions

The matrix-vector Equations (6) and (18) for the flow are solved by using the BEM using a fundamental solution which is a parametrix and by using the fundamental solution of Laplace’s equation in the DRBEM, respectively. Then, the matrix-vector Equation (26) for the temperature is solved by using the DRBEM. The equations are solved iteratively by taking $w=0$ and $T=0$ initially. The domain of the problem, the cross-section of the pipe is taken as a square $\left[0,1\right]\times \left[0,1\right]$ and discretized by using $N=100$ (25 on each side), and $N=180$ (45 on each side) constant boundary elements and $L=625$, and $L=2025$ interior nodes (interior domain mesh is obtained from the lines passing through the end points of the constant boundary elements) to obtain the solution using the BEM with a parametrix, and by using $N=180$ constant boundary elements and $L=2025$ interior nodes to obtain the solution using the DRBEM. The velocity and temperature isolines on the cross-section $\left[0,1\right]\times \left[0,1\right]$ of the pipe are simulated for several values of the Brinkman number as $Br=0,1$, the Hartmann number range $0\le Ha\le 10$, the viscosity parameter $B=0,1,2$ and the Hall parameter $m=0,3,8$. The velocity of the flow is relaxed as $w^{n+1}=\alpha w^{n+1}+\left(1-\alpha \right)w^n$ with a relaxation parameter $0<\alpha <1$ to slow down the sharp decrease of velocity magnitude in using the parametrix BEM. The domain integrals in the parametrix BEM formulation are computed by using the composite trapezoidal rule.

Figure 2(a, b) show the velocity and temperature behaviors as Hartmann number is increasing in which the velocity is obtained from the parametrix BEM and the DRBEM, respectively. It is observed that as Hartmann number increases, the velocity magnitude drops, due to the damping effect of the magnetic field with increasing intensity, this is a well-known flattening tendency of the MHD duct flow. It is observed that as Hartmann number increases, the temperature magnitude also drops. The magnitudes of the velocity obtained from the parametrix BEM are slightly less than the velocity magnitudes obtained with the DRBEM. This may be due to the computations of the domain integrals.

From Figure 3(a, b) and also from Figure 4(a, b), it is also observed that, as the viscosity parameter is increasing, the magnitudes of the velocity and the temperature drop due to the increment in the viscosity for different Hall parameter and Brinkmann number values. The same is observed as, the increase in the velocity magnitudes in the parametrix BEM is less than the increase in the velocity magnitudes obtained from the DRBEM.

From Figure 5(a, b) and also from Figure 6(a, b) it is seen that when the Hall parameter increases, the damping effect of the magnetic force decreases due to the term $\sigma ∕\left(1+m^2\right)$. That is, the velocity magnitude increases. It is also observed that as the Hall parameter increases, the effective conductivity decreases, which reduces the Joule dissipation. Therefore, the magnitude of the temperature also increases.

The effects of the Hall parameter and Hartmann number on the Nusselt number, $Nu$, are given on Tables 1–3 where $Nu=-\frac{1}{4T_m}$ and $T_m=\frac{1}{w_m}{\int }_0^1{\int }_0^{1}wTdxdy$ respectively, for the parametrix BEM and the DRBEM results.

It is observed that as Hartmann number increases, the values of Nusselt number are increasing as shown in Tables 1–3 since the values of the temperature decrease. However, as the Hall parameter is increasing, the values of Nusselt number are decreasing since the values of the temperature increase. But, it can be seen that, the effect of Hall parameter on Nusselt number may be neglected for small values of Hartmann number ($Ha<2$).

Figure 2: Equavelocity and isolines for $m=0$, $Br=0$ and $B=1$. (a) pBEm, (b) DRBEM.

Figure 3: Equavelocity and isolines for $Ha=3$, $m=0$ and $Br=0$. (a) pBEm, (b) DRBEM.

Figure 4: Equavelocity and isolines for $Ha=3$, $m=3$ and $Br=1$. (a) pBEm, (b) DRBEM.

Figure 5: Equavelocity and isolines for $Ha=3$, $Br=0$ and $B=1$. (a) pBEm, (b) DRBEM.

Figure 6: Equavelocity and isolines for $Ha=3$, $Br=1$ and $B=1$. (a) pBEm, (b) DRBEM.

The solutions obtained by the parametrix BEM and the DRBEM in terms of $w_{max}$ and $T_{min}$ slightly differ (still with almost $10^{-2}$ accuracy) in the values, they all catch the same behaviour of the velocity and the temperature of the MHD duct flow. When the number of boundary elements is increased from 100 to 180, Nusselt numbers are computed more accurately than DRBEM when FEM results [6] are taken as a basis. As the number of boundary elements is increased form 100 to 180 in parametrix BEM, the Nusselt numbers are close with almost $10^{-2}$ accuracy as shown in Tables 1–3 and in agreement with FEM results given in [6]. However, $m=0$ case needs relaxation parameter in using parametrix BEM and parametrix BEM has to compute domain integrals. This is why it takes more CPU time compared to DRBEM. Thus, the Table 4 shows that solving this MHD flow problem with temperature dependent viscosity by using the DRBEM is more time saving than solving with the parametrix BEM formulation due to the domain integrals computations. For this reason the DRBEM is preferable for solving this MHD flow problem.

Table 1: The effect of the Hall parameter, $m$, and Hartmann number, $Ha$, on Nusselt number, $Nu$, by using the parametrix BEM with 100 boundary elements (Br $=$ 0, B $=$ 1)

m	Ha	0.0	1.0	2.0	3.0	4.0	5.0
0.0	3.5662	3.6113	3.7533	3.7772	3.8586	3.9340
3.0	3.5662	3.5711	3.5896	3.6082	3.6247	3.6774
5.0	3.5662	3.5680	3.5738	3.5873	3.5981	3.6101
8.0	3.5662	3.5669	3.5692	3.5730	3.5830	3.5889

Table 2: The effect of the Hall parameter, $m$, and Hartmann number, $Ha$, on Nusselt number, $Nu$, by using the parametrix BEM with 180 boundary elements (Br $=$ 0, B $=$ 1)

m	Ha	0.0	1.0	2.0	3.0	4.0	5.0
0.0		3.6125	3.6595	3.8099	3.8342	3.9202	4.0000
3.0		3.6125	3.6176	3.6368	3.6563	3.6735	3.7294
5.0		3.6125	3.6144	3.6204	3.6344	3.6458	3.6583
8.0		3.6125	3.6132	3.6156	3.6196	3.6299	3.6361

Table 3: The effect of the Hall parameter, $m$, and Hartmann number, $Ha$, on Nusselt number, $Nu$, by using the DRBEM with 180 boundary elements (Br $=$ 0, B $=$ 1)

m	Ha	0.0	1.0	2.0	3.0	4.0	5.0
0.0		3.6276	3.6417	3.6825	3.7467	3.8294	3.9250
3.0		3.6276	3.6290	3.6333	3.6403	3.6500	3.6623
5.0		3.6276	3.6282	3.6298	3.6325	3.6363	3.6411
8.0		3.6276	3.6278	3.6285	3.6296	3.6311	3.6330

Table 4: CPU Times in pBEM and DRBEM with 180 boundary elements (Br $=$ 0, B $=$ 1)

Figure 7: Midline velocity profiles at $y=0.5$ using DRBEM.

Figure 8: Midline velocity profiles at $y=0.5$ using DRBEM.

Figures 7 and 8 show the velocity profiles along the midline of the rectangular duct obtained by using the DRBEM for increasing Hartmann number, viscosity parameter and Hall parameter. It can be seen from the Figure 7(a), as Hartmann number increases, the velocity magnitude drops for fixed values of the Hall parameter, Brinkman number and viscosity parameter. From the Figure 7(b) it is observed that, as the viscosity parameter is increasing, the magnitude of the velocity drops for fixed values of Hartmann number, Hall parameter and Brinkman number. Also, it can be seen from the Figure 8(a, b) that, as the Hall parameter increases, the flow magnitude also increases for Hartmann number, viscosity parameter and Brinkman number values as in the refence papers [10] and [6].

Figure 9: Midline temperature profiles at $y=0.5$ using DRBEM.

Figure 9 shows the temperature profiles along the midline of the rectangular duct obtained by using the DRBEM for increasing Hartmann number and viscosity parameter. It can be seen from the Figure 9(a), as Hartmann number increases the temperature magnitude drops for fixed values of the Hall parameter, Brinkman number and viscosity parameter. This is why the terms which are multiplied by these parameters form the force term for the diffusion equation for the temperature. Figure 9(b) the effect of viscosity parameter is almost negligible in the magnitude of the temperature for fixed values of Hartmann number, Hall parameter and Brinkman number.

Finally, the Figures 10–13 show velocity and temperature profiles obtained from both the parametrix BEM and the DRBEM for increasing Hartmann number, viscosity parameter and Hall parameter along the midline of the rectangular duct. Although the velocity magnitudes obtained by the parametrix BEM is less than the velocity magnitudes obtained by the DRBEM, the velocity profiles show the same flattening tendency behavior of the MHD duct flow. From Figure 11 it is seen that, as the viscosity parameter increases, the difference between the velocity magnitudes obtained from the parametrix BEM and the DRBEM increases due to the increment in the viscosity term $\overline{\mu }=e^{-BT}$. We can see from Figure 11 that, as Hartmann number increases, the difference between the velocity magnitudes obtained from the parametrix BEM and the DRBEM increases due to the increment in the nonlinearity term $\frac{H{a}^2}{1+m^2}w$, and from Figure 12, when Hall parameter is increasing, the effect of the nonlinearity term decreases and so the difference between the velocity magnitudes obtained from the parametrix BEM and the DRBEM decreases.

Figure 10: Midline velocity profiles at $y=0.5$, and $Ha=3$, $m=5$ and $Br=1$.

Figure 11: Midline velocity profiles at $y=0.5$, and $Br=1$, $B=1$, $m=3$.

Figure 12: Midline velocity profiles at $y=0.5$, and $Br=1$, $B=1$, $Ha=3$.

Figure 13: Midline temperature profiles at $y=0.5$, and $Br=1$, $B=1$, $m=3$

However, the difference in the velocity magnitudes from both methods results from the extra domain integral computations in the parametrix BEM procedure. The computations of the domain integrals in parametrix BEM must be numerical, and the composite trapezoidal rule has been used. They cause at least computational and round-off errors. This may be responsible for the discrepancies in parametrix BEM and DRBEM velocity values near the maximum. The increase of the number of constant elements, $N$, from 100 to 180 does not further reduce this discrepancy in the velocity values. However, from Figure 13, it can be seen that, the coincidence of the two methods, the parametrix BEM and the DRBEM is very well in terms of midline temperature values.

5 Conclusion

The steady, laminar, fully developed MHD flow of an incompressible, electrically conducting fluid with temperature dependent viscosity is studied in a rectangular duct together with its heat transfer. The flow and the energy equations of the problem are solved iteratively. Since the diffusion term contains variable viscosity parameter depending on the temperature exponentially, the velocity equation is first solved by using the boundary element method with a parametrix which treats the variable viscosity in the diffusion term directly. Thus, the use of this fundamental solution makes it possible to discretize the diffusion operator of the equation in its original form. The velocity and the energy equations are also solved by using the DRBEM with the fundamental solution of Laplace’s equation. Although the obtained results from parametrix BEM and the DRBEM differ in $w_{max}$ and $T_{min}$ values, both method capture the well known flattening property of the MHD duct flow as Hartmann number is increasing. As the value of Hartmann number is increasing, the velocity magnitude drops, however, increasing the Hall parameter, increases the velocity and temperature magnitudes. In addition, increase in the viscosity parameter reduces the magnitudes of the velocity and the temperature due to the exponential increment in the viscosity inside the diffusion term. Parametrix BEM with the same number of elements requires more CPU time than DRBEM. In that sense, the DRBEM is still preffered in solving MHD flow and heat transfer in rectangular pipes.

Acknowledgements

This work was supported by [The National Ph.D Fellowship Programme of The Scientific and Technological Research Council of Turkey (TÜBÍTAK)] under the Grant [2211].

References

[1] M. A. Al-Jawary, L. C. Wrobel, Numerical solution of two-dimensional mixed problems with variable coefficients by the boundary-domain integral and integro-differential equation methods, Engineering Analysis with Boundary Elements, 35 (2011), pp. 1279–1287.

[2] N. Alsoy-Akgün, M. Tezer-Sezgin, DRBEM and DQM solutions of natural convection flow in a cavity under a magnetic field, Progress in Computational Fluid Dynamics, 13 (2013), pp. 270–284.

[3] F. Capone, M. Gentile, Nonlinear stability analysis of convection for fluids with exponentially temperature-dependent viscosity, Acta Mechanica, 107 (1994), pp. 53–64.

[4] S. R. Cherlacola, K. Naikoti, M. M. Rashidi, MHD flow and heat transfer characteristics of Williamson nanofluid over a streching sheet with variable thickness and variable thermal conductivity, Transactions of A. Razmadze Mathematical Institute, 171 (2017), pp. 195–211.

[5] L. DRAGOŞ, Magnetofluid Dynamics, Abacus Press: England, 1975.

[6] C. Evcin, Ö. Uğur, M. Tezer-Sezgin, Determining the optimal parameters for the MHD flow and heat transfer with variable viscosity and Hall effect, Computers and Mathematics with Applications, (2018).

[7] D. R. Jenkins, Rolls versus squares in thermal convection of fluids with temperature-dependent viscosity, Journal of Fluid Mechanics, 178 (1987), pp. 491–506.

[8] P. W. Partridge, C. A. Brebbia, L. C. Wrobel, The Dual Reciprocity Boundary Element Method, Computational Mechanics Publications, Southampton, Boston, 1992.

[9] K. V. Prasad, K. Vajravelu, and H. Vaidya, Hall effect on MHD flow and heat transfer over a stretching sheet with variable thickness, International Journal for Computational Methods in Engineering Science and Mechanics, 17 (2006), pp. 288–297.

[10] M. E. Sayed-Ahmed, H. A. Attia, MHD Flow and Heat Transfer in a Rectangular Duct with Temperature Dependent Viscosity and Hall Effect, International Communication of Heat Transfer, 27 (2000), pp. 1177–1187.

[11] S. Shin, Y. I. Cho, W. K. Gingrich, W. Shyy, Numerical study of laminar heat transfer with temperature dependent fluid viscosity in a 2:1 rectangular duct, International Journal of Heat and Mass Transfer, 36 (1993), pp. 4365–4373.

[12] Ö. Türk, Chebyshev Spectral Collocation Method Approimation to Thermaly Coupled MHD equations, Süleyman Demirel University Journal of Natural and Applied Sciences, 22 (2018), pp. 355–366.

[13] Ö. Türk, M. Tezer-Sezgin, FEM solution of natural convection flow in square enclosures under magnetic field, International Journal of Numerical Methods for Heat & Fluid Flow, 23 (2013), pp. 844–866.

[14] C. Xie, J. P. Hartnett, Influence of variable viscosity of mineral oil on laminar heat transfer in a 2:1 rectangular duct, International Journal of Heat and Mass Transfer, 35 (1992), pp. 641–648.

Biographies

Elif Ebren Kaya received her BSc in Mathematics from the Middle East Technical University, Turkey in 2013. She has been a PhD student and a Research Asistant in the Department of Mathematics, Middle East Technical University, Turkey since 2013.

Münevver Tezer-Sezgin received her BSc in Mathematics from the Middle East Technical University, Turkey in 1974. She received her MSc and PhD degrees in Applied Mathematics in 1978 and 1980 from the University of Saskatchewan, Canada and the University of Calgary, Canada, respectively. She is retired and currently working as a part time Professor in the Department of Mathematics, Middle East Technical University, Turkey. She also has held one year visiting position at the University of Victoria, Canada, and short term visiting position at the Technical University of Darmstadt, Germany. She is the holder of Mustafa N. Parlar 1990 Research and 2014 Science Awards in applied mathematics and engineering.

European Journal of Computational Mechanics, Vol. 28 1&2, 97-122.
doi: 10.13052/ejcm1958-5829.28125
© 2019 River Publishers

Abstract

1 Introduction

2 Mathematical Formulation of the Problem

3 The BEM Implementations

3.1 The Parametrix BEM

3.2 The DRBEM Formulation

4 Numerical Results and Discussions

5 Conclusion

Acknowledgements

References

Biographies