Prediction of Better Flow Control Parameters in MHD Flows Using a High Accuracy Finite Difference Scheme

Show more

1. Introduction

The calculation and prediction of physical quantities to the desired accuracy is necessary in every engineering design because otherwise an under-performance together with an overshoot in the design cost. The study of magnetohydrodynamic flows with its heat transfer is important due to many practical applications in science and engineering. For instance, in many semiconducting single crystal growth processes, especially using the Czochralski technique, cylindrical shaped crystal is grown from the melt containing the material in the molten (liquid) state. In such systems, calculation of heat transfers is fairly complex due to the nonlinear coupling of the momentum with the heat diffusion and advection apart from radiative transfers. These semiconducting crystals are mostly used as substrates in electronic industry. However, during the crystal growth process, thermal stresses will be present if heat transfer is not well controlled [1] . Consequently, the resulting semiconducting single crystal substrates will have unavoidable dislocations which is an undesirable feature for applications. This means accurate calculation of the flow field, determination of temperature isotherms and hence precise estimates of heat transfer coefficients is of crucial importance in many practical applications. The accuracy of the final computed results from numerical simulation emerge from both the modeling (or approxi- mation methods) and also on the usage of a particular discretization method. In this paper, we address both the issues in order to have minimum approximation in modeling of physical system together with a highly accurate numerical solution by using a high order discretization procedure.

Many industrial applications include liquid metals in metallurgical processes such as stirring, pumping, casting, suppressing etc. [2] . The physical importance of the MHD flow is due to the coupling of Maxwell’s equations with momentum equation leading to mathematical complexity and fascinating physical phenomena such as flow acceleration or, suppression of flow, or modified heat transfer. The magnetic Reynolds number ${R}_{m}$ is an important physical parameter which varies over several orders of magnitude depending on the nature of the system under consideration. For example, ${R}_{m}$ is large in astrophysical and geophysical flows while in metallurgical flow problems or in crystal growth processes, it is very small. Numerous theoretical, numerical and experimental studies have been carried out on MHD flows and bluff body flow dynamics. When an electrically conducting fluid moves through a magnetic field the fluid will experience a Lorentz force, leading to Joule dissipation that eventually reduces the kinetic energy of the fluid. That is, the velocity distribution is regenerated by the Lorentz force and suppress the turbulence in the wake region. The Lorentz force can reduce the drag force effectively and suppress the flow separation in fluid flows as evidenced by the following data.

The solution methodology used for MHD flows can be classified in to three categories namely a) getting approximate analytical solutions by perturbation techniques or otherwise, b) getting the boundary layer solution and c) obtaining the full solution in the entire domain of interest. In the case of boundary layer solution approach, the solutions are obtained within the boundary layer of the flow, and the detailed flow velocities and temperature will not be available in the domain of fluid flow. Under the first category a), limited attempts are made using Laplace transform technique [3] and perturbation techniques [4] . Under the second category b), the boundary layer solutions for different types of MHD flows and associated heat transfer has been studied by many authors [3] [5] - [13] , wherein the governing nonlinear partial differential equations are transformed to a set of ordinary differential equations and they are solved using Runge Kutta methods. In addition, many authors have used lower order accurate finite difference methods for the study of MHD flows and in particular with heat transfer effects [14] . Various research groups target on different numerical schemes in order to capture the flow and wake dynamics. One important feature to capture in the analysis is the $Re$ at which separation starts. The orientation of magnetic field direction (parallel/transverse) to the free stream velocity have different contribution to stabilize the flow [15] [16] which is also reported experimentally [17] [18] . The effect of magnetic field on the stabilization of unsteady flows has been studied by a linear stability analysis and concluded that the magnetic field can have both damping and enhancing influence on the instability of the flow [19] [20] . In a recent experimental findings it is reported that the pressure drag coefficient can be increased with the help of magnetic field [21] . The numerical studies on the effect of magnetic field on the flow around cylinder has been reported and where the accuracy will not exceed second order [22] [23] [24] . The MHD flow past circular cylinder has been studied using an extension to immersed boundary method and concluded that the method is an efficient one to study the MHD problems [25] . Earth’s liquid core [26] , liquids involved in melting processes, crystal growth processing are some of the branches where the low Prandtl number study is necessary $\left(Pr\ll 1\right)$ . Fluids connected with everyday lives such as water, oils, alcohols have the $Pr~1$ while molten alloys like GaInSn will have a Prandtl number of 0.022. The forced convective heat transfer studies due to power-law fluid flow past cylinder is analyzed [27] and for the normal Newtonian fluid is analyzed [28] [29] . A Lax-Wendroff type matrix distribution method is used and the numerical solution to the MHD flow with heat transfer is studied [30] and concluded that the method is robust and accurate, where a combination of second order and third-order upwind schemes were used. In addition to change of drag coefficient due to magnetic field, experimental studies revealed that heat transfer rate can be reduced using applied magnetic field [31] [32] . Without magnetic field effects, the forced convection heat transfer from circular cylinder due to air and liquids has been analyzed experimentally and the data is correlated with $Pr$ and $Re$ [33] . It is also reported experimentally that an applied magnetic field will strongly reduce the intensity of velocity fluctuations in certain regions of the liquid metal flow [34] . The flow and heat transfer around a circular cylinder in presence of magnetic field has been analyzed experimentally [35] [36] and later numerically using a spectral method and observed that the applied magnetic field will reduce the oscillating amplitude of lift coefficients [37] . The study of liquid metal flow past a cylinder inside a duct provided means of understanding the effect of magnetic field on $Re$ and found that the imposed magnetic field will shift the appearance of flow instabilities to higher $Re$ [38] . There are a few reports on the use of high order numerical schemes in curvilinear systems. The problem of flow past impulsively started cylinder is analyzed using a semi- compact high order scheme and concluded that the scheme is capable of capturing the physical phenomena [39] . The classical problem of flow past cylinder (without external magnetic field) has been studied using a high order finite difference scheme wherein it is reported that the numerical scheme at a fairly coarser grid is sufficient to get accurate results [40] . In the present work we develop a new fourth order accurate discretization scheme for the full-MHD model in the streamfunction-vorticity formulation in cylindrical polar coordinate system and the effectiveness of the scheme is tested in terms of predicting better control parameters in the MHD flows. The computed velocities from $\psi $ and $\omega $ will be having higher numerical accuracy. Experimental measurement of velocities of the fluids of MHD flows with magnetic field can be found in the literature [41] .

2. Governing Equations for MHD Flows

In order to study the steady state flow properties and heat transfer in flows of electrically conducting fluids in presence of magnetic field, the governing equations would be the modified Navier-Stokes equation (in which additional body force terms due to magnetic field) which is coupled with the Maxwells equations of electrodynamics together with energy equation. The non- dimensionalized equations under consideration are

$\nabla \cdot q=0$ (1)

$\left(q\cdot \nabla \right)q=-\nabla p+\frac{2}{Re}{\nabla}^{2}q+\beta \left(j\times H\right)$ (2)

$\nabla \cdot H=0$ (3)

$\nabla \times H=j+\frac{\partial D}{\partial t}$ (4)

$\nabla \cdot j=0$ (5)

$j=\frac{{R}_{m}}{2}\left[E+\left(q\times H\right)\right]$ (6)

$\left(q\cdot \nabla \right)T=\frac{2}{RePr}{\nabla}^{2}T$ (7)

where the following definitions are used in the process of non-dimensionalization.

$q=\frac{{q}^{\prime}}{{U}_{\infty}},\text{\hspace{1em}}p=\frac{p\text{'}}{\rho {U}_{\infty}^{2}},\text{\hspace{1em}}r=\frac{{r}^{\prime}}{a},\text{\hspace{1em}}T=\frac{{T}^{\prime}-{T}_{\infty}}{{T}_{s}-{T}_{\infty}},$

$H=\frac{{H}^{\prime}}{{H}_{\infty}},\text{\hspace{1em}}E=\frac{{E}^{\prime}}{{E}_{\infty}},\text{\hspace{1em}}j=\frac{{j}^{\prime}}{{j}_{\infty}}$

in which the primed variables denote the respective quantities with dimensions. The Alfvén number $\beta $ is a dimensionless quantity characterizing the flow in presence of magnetic field and it is the ratio of the speed of the Alfvén wave to the speed of the main stream fluid. The last term in Equation (2) is the body force term which is nonlinearly coupled with other equations in the list. Equation (1) is due to incompressibility condition while Equations ((3) and (4)) are the Maxwell’s equations to be satisfied for the applied magnetic field $H$ . In the Amperes law (4), the second term is due to electric displacement vector $D$ . Since displacement current density is negligible in fluid flows the second term in the right hand side of Equation (4) can be dropped. For steady state conditions, the electrodynamic continuity equation is (5) and the conduction current density $j$ has to satisfy the Ohms law which is Equation (6). First it is noted that Equations ((1) to (6)) are coupled. Now, if heat transfer analysis is to be carried out, then the energy transport equation to be solved is given by Equation (7). If we solve (7) by treating the velocity $q$ as coupled with Equations ((1) to (6)) then we can get the natural convection heat transfer properties. On the other hand, if we first solve the set of Equations ((1) to (6)) and provide this solution to Equation (7) as its input, then we can study the forced convection properties. In any case, first we need a discretization scheme to solve the governing equations. In the following, we propose to device a solution method based on streamfunction-vorticity approach and which is suitable for two- dimensional flow simulations. In cylindrical polar system, the velocity and applied magnetic field are

$q=\left({q}_{r}\mathrm{,}{q}_{\theta}\mathrm{,0}\right)$ (8)

$H=\left({h}_{r}\mathrm{,}{h}_{\theta}\mathrm{,0}\right)$ (9)

To satisfy the incompressibility condition, the velocity components can be expressed in terms of a scalar streamfunction $\psi $ given by

${q}_{r}\left(r\mathrm{,}\theta \right)=\frac{1}{r}\cdot \frac{\partial \psi \left(r\mathrm{,}\theta \right)}{\partial \theta}$ (10)

${q}_{\theta}\left(r\mathrm{,}\theta \right)=-\text{\hspace{0.05em}}\frac{\partial \psi \left(r\mathrm{,}\theta \right)}{\partial r}$ (11)

Since Equations ((1) and (3)) exhibit solenoidal nature, a definition similar to (10) and (11) can be used by making use a scalar called magnetic streamfunction, denoted by $A$ as below.

${h}_{r}\left(r\mathrm{,}\theta \right)=\frac{1}{r}\cdot \frac{\partial A\left(r\mathrm{,}\theta \right)}{\partial \theta}$ (12)

${h}_{\theta}\left(r\mathrm{,}\theta \right)=-\frac{\partial A\left(r\mathrm{,}\theta \right)}{\partial r}$ (13)

By taking curl of Equation (8), the vorticity is obtained as

$\nabla \times q=\omega $ (14)

which, upon using (10) and (11), we get the differential equation for velocity as below.

$\frac{{q}_{\theta}}{r}+\frac{\partial {q}_{\theta}}{\partial r}-\frac{1}{r}\left(\frac{\partial {q}_{r}}{\partial \theta}\right)=\omega $ (15)

and where $\omega $ has components only in the z-direction and is a function of $r$ and $\theta $ of the cylindrical polar coordinate system. We make use of Equation (14) in (2) by applying the curl operator to (2) which results in

$-\nabla \times \left(q\times \omega \right)=-\frac{2}{Re}\cdot \nabla \times \left(\nabla \times \omega \right)+\beta \cdot \left[\nabla \times \left\{\nabla \times H\right\}\times H\right]$ (16)

where we have also used (4) after neglecting the displacement vector. In the present work, we consider the case of applying only an external magnetic field but no electric field is applied. In addition, for the case of two-dimensional flows, it can be shown that $E=0$ if only magnetic field is applied to the electrically conducting flow. Then, from (4) and (6), we have

$\nabla \times H=\frac{{R}_{m}}{2}\cdot \left(q\times H\right)$ (17)

Using (17) in (16) and then expanding the curl operator in the cylindrical system and by making use of (9) we obtain the differential equation for vorticity as given below.

$\begin{array}{l}\left[{q}_{r}\left(\frac{\partial \omega}{\partial r}\right)+\frac{{q}_{\theta}}{r}\left(\frac{\partial \omega}{\partial \theta}\right)\right]-\frac{2}{Re}\left[\frac{{\partial}^{2}\omega}{\partial {r}^{2}}+\frac{1}{r}\left(\frac{\partial \omega}{\partial r}\right)+\frac{1}{{r}^{2}}\left(\frac{{\partial}^{2}\omega}{\partial \theta}\right)\right]\\ =\frac{\beta {R}_{m}}{2}\cdot [-\frac{{h}_{r}^{2}{q}_{\theta}}{r}-2{h}_{r}{q}_{\theta}\left(\frac{\partial {h}_{r}}{\partial r}\right)-{h}_{r}^{2}\left(\frac{\partial {q}_{\theta}}{\partial r}\right)+\frac{{h}_{r}{h}_{\theta}{q}_{r}}{r}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+{h}_{\theta}{q}_{r}\left(\frac{\partial {h}_{r}}{\partial r}\right)+{h}_{r}{q}_{r}\left(\frac{\partial {h}_{\theta}}{\partial r}\right)+{h}_{r}{h}_{\theta}\left(\frac{\partial {q}_{r}}{\partial r}\right)+\frac{2{h}_{\theta}{q}_{r}}{r}\left(\frac{\partial {h}_{\theta}}{\partial \theta}\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{{h}_{\theta}^{2}}{r}\left(\frac{\partial {q}_{r}}{\partial \theta}\right)-\frac{{h}_{r}{h}_{\theta}}{r}\left(\frac{\partial {q}_{\theta}}{\partial \theta}\right)-\frac{{h}_{\theta}}{r}\left(\frac{\partial {h}_{r}}{\partial \theta}\right)-\frac{{q}_{\theta}{h}_{r}}{r}\left(\frac{\partial {h}_{\theta}}{\partial \theta}\right)]\end{array}$ (18)

Now, if we simplify the Equation (17) using (8) and (9) we get the differential equation for magnetic field as follows.

$\frac{{h}_{\theta}}{r}+\frac{\partial {h}_{\theta}}{\partial r}-\frac{1}{r}\cdot \left(\frac{\partial {h}_{r}}{\partial \theta}\right)=\frac{{R}_{m}\text{\hspace{0.05em}}\left({q}_{r}{h}_{\theta}-{q}_{\theta}{h}_{r}\right)}{2}$ (19)

Similarly, Equation (7) is simplified using (10) and (11) to yield the differential equation for temperature as given below.

${q}_{r}\left(\frac{\partial T}{\partial r}\right)+{q}_{\theta}\left(\frac{\partial T}{\partial \theta}\right)=\frac{2}{Re\cdot Pr}\left[\frac{{\partial}^{2}T}{\partial {r}^{2}}+\frac{1}{r}\left(\frac{\partial T}{\partial r}\right)+\frac{1}{{r}^{2}}\left(\frac{{\partial}^{2}T}{\partial {\theta}^{2}}\right)\right]$ (20)

Finally we have a set of four coupled partial differential Equations ((15), (18)-(20)) having several nonlinear terms in the vorticity differential equation. The components of $q$ and $H$ appearing in these four equations are given by (10) to (13). The above set of equations are generally valid for MHD flows in cylindrical geometry.

3. Discretization Procedure

In this section we outline the discretization procedure for the discretization of governing equations. The implementation of fourth order finite differences to the governing equations for the present problem requires considerable exercise because of curvilinear coordinate system. From the Taylor series, for the x and h-direction, if we denote h and k as the grid spacing along radial and angular direction then, the fourth order accurate finite differences for the first and second derivatives can be written as

$\frac{\partial \Psi}{\partial \xi}={D}_{\xi}\Psi -\frac{{h}^{2}}{6}\frac{{\partial}^{3}\Psi}{\partial {\xi}^{3}}+\mathcal{O}\left({h}^{4}\right)$ (21)

$\frac{{\partial}^{2}\Psi}{\partial {\xi}^{2}}={D}_{\xi}^{2}\Psi -\frac{{h}^{2}}{12}\frac{{\partial}^{4}\Psi}{\partial {\xi}^{4}}+\mathcal{O}\left({h}^{4}\right)$ (22)

$\frac{\partial \Psi}{\partial \eta}={D}_{\eta}\Psi -\frac{{k}^{2}}{6}\frac{{\partial}^{3}\Psi}{\partial {\eta}^{3}}+\mathcal{O}\left({k}^{4}\right)$ (23)

$\frac{{\partial}^{2}\Psi}{\partial {\eta}^{2}}={D}_{\eta}^{2}\Psi -\frac{{h}^{2}}{12}\frac{{\partial}^{4}\Psi}{\partial {\eta}^{4}}+\mathcal{O}\left({k}^{4}\right)$ (24)

where ${D}_{\xi}$ , ${D}_{\xi}^{2}$ , ${D}_{\eta}$ and ${D}_{\eta}^{2}$ denote the standard second order central difference operators which are given by

${D}_{\xi}{\Psi}_{i,\text{\hspace{0.05em}}j}=\frac{{\Psi}_{i+1,\text{\hspace{0.05em}}j}-{\Psi}_{i-1,\text{\hspace{0.05em}}j}}{2h}$ (25)

${D}_{\xi}^{2}{\Psi}_{i,\text{\hspace{0.05em}}j}=\frac{{\Psi}_{i+1,\text{\hspace{0.05em}}j}-2{\Psi}_{i,\text{\hspace{0.05em}}j}+{\Psi}_{i-1,\text{\hspace{0.05em}}j}}{{h}^{2}}$ (26)

${D}_{\eta}{\Psi}_{i,\text{\hspace{0.05em}}j}=\frac{{\Psi}_{i,\text{\hspace{0.05em}}j+1}-{\Psi}_{i,\text{\hspace{0.05em}}j-1}}{2k}$ (27)

${D}_{\eta}^{2}{\Psi}_{i,\text{\hspace{0.05em}}j}=\frac{{\Psi}_{i,\text{\hspace{0.05em}}j+1}-2{\Psi}_{i,\text{\hspace{0.05em}}j}+{\Psi}_{i,\text{\hspace{0.05em}}j-1}}{{k}^{2}}$ (28)

The second order central difference operators for cross-derivatives can be obtained from Taylor series, and are given by

${D}_{\xi}{D}_{\eta}{\Psi}_{i\mathrm{,}\text{\hspace{0.05em}}j}=\frac{1}{4hk}\left[{\Psi}_{i+\mathrm{1,}\text{\hspace{0.05em}}j+1}-{\Psi}_{i+\mathrm{1,}\text{\hspace{0.05em}}j-1}-{\Psi}_{i-\mathrm{1,}\text{\hspace{0.05em}}j+1}+{\Psi}_{i-\mathrm{1,}\text{\hspace{0.05em}}j-1}\right]$ (29)

${D}_{\xi}^{2}{D}_{\eta}{\Psi}_{i,\text{\hspace{0.05em}}j}=\frac{1}{2{h}^{2}k}\left[{\Psi}_{i+1,\text{\hspace{0.05em}}j+1}-{\Psi}_{i+1,\text{\hspace{0.05em}}j-1}+{\Psi}_{i-1,\text{\hspace{0.05em}}j+1}-{\Psi}_{i-1,\text{\hspace{0.05em}}j-1}-2\left({\Psi}_{i,j+1}-{\Psi}_{i,j-1}\right)\right]$ (30)

${D}_{\xi}{D}_{\eta}^{2}{\Psi}_{i,\text{\hspace{0.05em}}j}=\frac{1}{2h{k}^{2}}\left[{\Psi}_{i+1,\text{\hspace{0.05em}}j+1}+{\Psi}_{i+1,\text{\hspace{0.05em}}j-1}-{\Psi}_{i-1,\text{\hspace{0.05em}}j+1}-{\Psi}_{i-1,\text{\hspace{0.05em}}j-1}-2\left({\Psi}_{i+1,\text{\hspace{0.05em}}j}-{\Psi}_{i-1,\text{\hspace{0.05em}}j}\right)\right]$ (31)

$\begin{array}{c}{D}_{\xi}^{2}{D}_{\eta}^{2}{\Psi}_{i\mathrm{,}\text{\hspace{0.05em}}j}=\frac{1}{{h}^{2}{k}^{2}}[{\Psi}_{i+\mathrm{1,}\text{\hspace{0.05em}}j+1}+{\Psi}_{i+\mathrm{1,}\text{\hspace{0.05em}}j-1}+{\Psi}_{i-\mathrm{1,}\text{\hspace{0.05em}}j+1}+{\Psi}_{i-\mathrm{1,}\text{\hspace{0.05em}}j-1}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-2\left({\Psi}_{i+\mathrm{1,}\text{\hspace{0.05em}}j}+{\Psi}_{i-\mathrm{1,}\text{\hspace{0.05em}}j}+{\Psi}_{i\mathrm{,}\text{\hspace{0.05em}}j+1}+{\Psi}_{i\mathrm{,}\text{\hspace{0.05em}}j-1}\right)+4{\Psi}_{i\mathrm{,}\text{\hspace{0.05em}}j}]\end{array}$ (32)

3.1. Discretization of Streamfunction Equation

Equation (15) along with definitions (10) and (11) gives

$\frac{1}{r}\left[\frac{\partial \psi}{\partial r}\right]+\frac{\partial}{\partial r}\left[\frac{\partial \psi}{\partial r}\right]-\frac{1}{r}\left[\frac{\partial}{\partial \theta}\left\{\frac{1}{r}\frac{\partial \psi}{\partial \theta}\right\}\right]=\omega $ (33)

Now the transformation $r={\text{e}}^{\text{\pi}\xi}$ and $\theta =\text{\pi}\eta $ are implemented followed y the usage of D-operator from Equations ((21) to (28)) so that the above equation becomes.

${D}_{\xi}^{2}{\psi}_{i\mathrm{,}j}+{D}_{\eta}^{2}{\psi}_{i\mathrm{,}j}+{s}_{i\mathrm{,}j}={\left(\frac{{h}^{2}}{12}\frac{{\partial}^{4}\psi}{\partial {\xi}^{4}}+\frac{{k}^{2}}{12}\frac{{\partial}^{4}\psi}{\partial {\eta}^{4}}\right)}_{i\mathrm{,}j}+\mathcal{O}\left({h}^{4}\mathrm{,}{k}^{4}\right)$ (34)

where the RHS of the above equation is the truncation error term and $\mathcal{O}\left({h}^{4}\mathrm{,}{k}^{4}\right)$ indicates that the above equation is of the order of ${h}^{4}$ in $\xi $ direction and of the order of ${k}^{4}$ in $\eta $ direction. and

${s}_{i,j}={\text{\pi}}^{2}{\text{e}}^{2\text{\pi}\xi}{\omega}_{i,j}$ (35)

Since the Equation (34) inherently contain higher order derivatives of the stream function, they are computed as

$\frac{{\partial}^{3}\psi}{\partial {\xi}^{3}}=-\frac{{\partial}^{3}\psi}{\partial \xi \partial {\eta}^{2}}-\frac{\partial s}{\partial \xi}=-{D}_{\xi}{D}_{\eta}^{2}\psi -{D}_{\xi}s$ (36)

$\frac{{\partial}^{4}\psi}{\partial {\xi}^{4}}=-\frac{{\partial}^{4}\psi}{\partial {\xi}^{2}\partial {\eta}^{2}}-\frac{{\partial}^{2}s}{\partial {\xi}^{2}}=-{D}_{\xi}^{2}{D}_{\eta}^{2}\psi -{D}_{\xi}^{2}s$ (37)

$\frac{{\partial}^{3}\psi}{\partial {\eta}^{3}}=-\frac{{\partial}^{3}\psi}{\partial {\xi}^{2}\partial \eta}-\frac{\partial s}{\partial \eta}=-{D}_{\xi}^{2}{D}_{\eta}\psi -{D}_{\eta}s$ (38)

$\frac{{\partial}^{4}\psi}{\partial {\eta}^{4}}=-\frac{{\partial}^{4}\psi}{\partial {\xi}^{2}\partial {\eta}^{2}}-\frac{{\partial}^{2}s}{\partial {\eta}^{2}}=-{D}_{\xi}^{2}{D}_{\eta}^{2}\psi -{D}_{\eta}^{2}s$ (39)

Substituting Equations (36)-(39) in the discretized streamfunction Equation (34) and applying the D-operator from Equations (25)-(32), we get the fourth order accurate discretized representation of Equation (33) as given below.

$\begin{array}{l}2\left(2z-{h}^{2}-{k}^{2}\right){\psi}_{0}+\left({k}^{2}-2z\right)\left({\psi}_{1}+{\psi}_{3}\right)+\left({h}^{2}-2z\right)\left({\psi}_{2}+{\psi}_{4}\right)+z{\displaystyle \underset{i=5}{\overset{8}{\sum}}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{\psi}_{i}\\ =-{h}^{2}{k}^{2}{\text{\pi}}^{2}{\text{e}}^{2\text{\pi}\xi}\left\{\left(1+4{\text{\pi}}^{2}x\right){\omega}_{0}+4\text{\pi}x\text{\hspace{0.05em}}{D}_{\xi}\left[\omega \right]+\left(x{D}_{\xi}^{2}+y{D}_{\eta}^{2}\right)\left[\omega \right]\right\}\end{array}$ (40)

where $x={h}^{2}/12$ , $y={k}^{2}/12$ and $z=x+y$ , $\underset{i=0}{\overset{8}{\sum}}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{\psi}_{i$ represents the 8 nearest neighboring points and a centre point in the computational domain.

3.2. Discretization of Vorticity Equation

Equation (18) is simplified using Equations (10)-(13) and then $r={\text{e}}^{\text{\pi}\xi}$ with $\theta =\text{\pi}\eta $ is applied. Then, we can arrive at the following equation in terms of $\xi $ and $\eta $ .

$\frac{{\partial}^{2}\omega}{\partial {\xi}^{2}}+\frac{{\partial}^{2}\omega}{\partial {\eta}^{2}}+c\frac{\partial \omega}{\partial \xi}+d\frac{\partial \omega}{\partial \eta}=F$ (41)

where the coefficients ${c}_{i\mathrm{,}j}$ , ${d}_{i\mathrm{,}j}$ and ${F}_{i\mathrm{,}j}$ are

$c=-\frac{Re}{2}\left(\frac{\partial \psi}{\partial \eta}\right)$ (42)

$d=\frac{Re}{2}\left(\frac{\partial \psi}{\partial \xi}\right)$ (43)

and

$\begin{array}{c}F=\frac{-{\text{e}}^{-2\text{\pi}\xi}}{4{\text{\pi}}^{2}}Re\cdot {R}_{m}\cdot \beta [-\frac{\partial \psi}{\partial \eta}\frac{\partial A}{\partial \eta}\frac{{\partial}^{2}A}{\partial {\xi}^{2}}-\frac{\partial \psi}{\partial \xi}\frac{\partial A}{\partial \xi}\frac{{\partial}^{2}A}{\partial {\eta}^{2}}+2\pi \frac{\partial \psi}{\partial \eta}\frac{\partial A}{\partial \eta}\frac{\partial A}{\partial \xi}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{\partial A}{\partial \xi}\frac{\partial \psi}{\partial \eta}\frac{{\partial}^{2}A}{\partial \xi \partial \eta}+\frac{\partial \psi}{\partial \xi}\frac{\partial A}{\partial \eta}\frac{{\partial}^{2}A}{\partial \xi \partial \eta}-2\frac{\partial A}{\partial \xi}\frac{\partial A}{\partial \eta}\frac{{\partial}^{2}\psi}{\partial \xi \partial \eta}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-2\text{\pi}\frac{\partial \psi}{\partial \xi}{\left(\frac{\partial A}{\partial \eta}\right)}^{2}+{\left(\frac{\partial A}{\partial \eta}\right)}^{2}\frac{{\partial}^{2}\psi}{\partial {\xi}^{2}}+{\left(\frac{\partial A}{\partial \xi}\right)}^{2}\frac{{\partial}^{2}\psi}{\partial {\eta}^{2}}]\end{array}$ (44)

On substitution of Equations (21)-(24) in the vorticity Equation (41) we get the discretized form as below.

${D}_{\xi}^{2}{\omega}_{i,j}+{D}_{\eta}^{2}{\omega}_{i,j}+{c}_{i,j}{D}_{\xi}{\omega}_{i,j}+{d}_{i,j}{D}_{\eta}{\omega}_{i,j}+{\tau}_{i,j}={F}_{i,j}$ (45)

and the truncation error $\tau $ in (45) is

${\tau}_{i\mathrm{,}j}=-{\left(c\frac{{h}^{2}}{6}\frac{{\partial}^{3}\omega}{\partial {\xi}^{3}}+\frac{{h}^{2}}{12}\frac{{\partial}^{4}\omega}{\partial {\xi}^{4}}+d\frac{{k}^{2}}{6}\frac{{\partial}^{3}\omega}{\partial {\eta}^{3}}+\frac{{k}^{2}}{12}\frac{{\partial}^{4}\omega}{\partial {\eta}^{4}}\right)}_{i\mathrm{,}j}+\mathcal{O}\left({h}^{4}\mathrm{,}{k}^{4}\right)$ (46)

Now, to eliminate the third and fourth derivatives of $\omega $ that are present in the truncation error term, we differentiate (41) once and twice with respect to $\xi $ and $\eta $ respectively, to yield the following.

$\frac{{\partial}^{3}\omega}{\partial {\xi}^{3}}=-\frac{{\partial}^{3}\omega}{\partial \xi \partial {\eta}^{2}}-c\frac{{\partial}^{2}\omega}{\partial {\xi}^{2}}-d\frac{{\partial}^{2}\omega}{\partial \xi \partial \eta}-\frac{\partial c}{\partial \xi}\frac{\partial \omega}{\partial \xi}-\frac{\partial d}{\partial \xi}\frac{\partial \omega}{\partial \eta}-\frac{\partial F}{\partial \xi}$ (47)

$\begin{array}{c}\frac{{\partial}^{4}\omega}{\partial {\xi}^{4}}=-\frac{{\partial}^{4}\omega}{\partial {\xi}^{2}\partial {\eta}^{2}}+c\frac{{\partial}^{3}\omega}{\partial \xi \partial {\eta}^{2}}-d\frac{{\partial}^{3}\omega}{\partial {\xi}^{2}\partial \eta}+\left({c}^{2}-2\frac{\partial c}{\partial \xi}\right)\frac{{\partial}^{2}\omega}{\partial {\xi}^{2}}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\left(cd-2\frac{\partial d}{\partial \xi}\right)\frac{{\partial}^{2}\omega}{\partial \xi \partial \eta}+\left(c\frac{\partial c}{\partial \xi}-\frac{{\partial}^{2}c}{\partial {\xi}^{2}}\right)\frac{\partial \omega}{\partial \xi}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\left(c\frac{\partial d}{\partial \xi}-\frac{{\partial}^{2}d}{\partial {\xi}^{2}}\right)\frac{\partial \omega}{\partial \eta}+c\frac{\partial F}{\partial \xi}-\frac{{\partial}^{2}F}{\partial {\xi}^{2}}\end{array}$ (48)

$\frac{{\partial}^{3}\omega}{\partial {\eta}^{3}}=-\frac{{\partial}^{3}\omega}{\partial {\xi}^{2}\partial \eta}-c\frac{{\partial}^{2}\omega}{\partial \xi \partial \eta}-d\frac{{\partial}^{2}\omega}{\partial {\eta}^{2}}-\frac{\partial c}{\partial \eta}\frac{\partial \omega}{\partial \xi}-\frac{\partial d}{\partial \eta}\frac{\partial \omega}{\partial \eta}-\frac{\partial F}{\partial \eta}$ (49)

$\begin{array}{c}\frac{{\partial}^{4}\omega}{\partial {\eta}^{4}}=-\frac{{\partial}^{4}\omega}{\partial {\xi}^{2}\partial {\eta}^{2}}-c\frac{{\partial}^{3}\omega}{\partial \xi \partial {\eta}^{2}}+d\frac{{\partial}^{3}\omega}{\partial {\xi}^{2}\partial \eta}+\left(cd-2\frac{\partial c}{\partial \eta}\right)\frac{{\partial}^{2}\omega}{\partial \xi \partial \eta}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\left({d}^{2}-2\frac{\partial d}{\partial \eta}\right)\frac{{\partial}^{2}\omega}{\partial {\eta}^{2}}+\left(d\frac{\partial c}{\partial \eta}-\frac{{\partial}^{2}c}{\partial {\eta}^{2}}\right)\frac{\partial \omega}{\partial \xi}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\left(d\frac{\partial d}{\partial \eta}-\frac{{\partial}^{2}d}{\partial {\eta}^{2}}\right)\frac{\partial \omega}{\partial \eta}+d\frac{\partial F}{\partial \eta}-\frac{{\partial}^{2}F}{\partial {\eta}^{2}}\end{array}$ (50)

Now we substitute (42)-(50) in (45) and then apply the D-operators from (25)-(32). Then, the resulting equation is simplified to arrive at the fourth order accurate representation of vorticity Equation (41) as below.

$\begin{array}{l}-8\left(e{k}^{2}+f{h}^{2}-2z\right){\omega}_{0}+\left(4{k}^{2}e+2h{k}^{2}g-8z-4hzc\right){\omega}_{1}\\ +\left(4{h}^{2}f+2{h}^{2}ko-8z-4kzd\right){\omega}_{2}+\left(4{k}^{2}e-2h{k}^{2}g-8z+4hzc\right){\omega}_{3}\\ +\left(4{h}^{2}f-2{h}^{2}ko-8z+4kzd\right){\omega}_{4}+\left(4z+2hcz+2kdz+hkl\right){\omega}_{5}\\ +\left(4z-2hcz+2kdz-hkl\right){\omega}_{6}+\left(4z-2hcz-2kdz+hkl\right){\omega}_{7}\\ +\left(4z+2hcz-2kdz-hkl\right){\omega}_{8}\\ =\frac{-{h}^{2}{k}^{2}Re\text{\hspace{0.05em}}{R}_{m}\text{\hspace{0.05em}}\beta {\text{e}}^{-2\text{\pi}\xi}}{{\text{\pi}}^{2}}\left\{G+\left(x\text{\hspace{0.05em}}c\text{\hspace{0.05em}}{D}_{\xi}+y\text{\hspace{0.05em}}d\text{\hspace{0.05em}}{D}_{\eta}\right)\mathrm{[}G\mathrm{]}+\left(x{D}_{\xi}^{2}+y{D}_{\eta}^{2}\right)\left[G\right]\right\}\end{array}$ (51)

where

$\begin{array}{c}G=-{D}_{\eta}\left[\psi \right]{D}_{\eta}{D}_{\xi}^{2}\left[A\right]-{D}_{\xi}\left[\psi \right]{D}_{\xi}{D}_{\eta}^{2}\left[A\right]+2\text{\pi}{D}_{\eta}\left[\psi \right]{D}_{\eta}{D}_{\xi}\left[A\right]\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+{D}_{\xi}\left[A\right]{D}_{\eta}\left[\psi \right]{D}_{\xi}{D}_{\eta}\left[A\right]+{D}_{\xi}\left[\psi \right]{D}_{\eta}\left[A\right]{D}_{\xi}{D}_{\eta}\left[A\right]\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-2{D}_{\xi}\left[A\right]{D}_{\eta}\left[A\right]{D}_{\xi}{D}_{\eta}\left[\psi \right]-2\text{\pi}{D}_{\xi}\left[\psi \right]{D}_{\eta}\left[A\right]{D}_{\eta}\left[A\right]\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+{D}_{\eta}\left[A\right]\text{\hspace{0.05em}}{D}_{\eta}\left[A\right]{D}_{\xi}^{2}\left[\psi \right]+{D}_{\xi}\left[A\right]{D}_{\xi}\left[A\right]{D}_{\eta}^{2}\left[\psi \right]\end{array}$ (52)

where the values of $\psi $ around ${\psi}_{i\mathrm{,}j}$ are indexed as shown in the 2D compact stencil (Figure 1). Further, the coefficients required in the above equations are as follows.

$e=1+x\left\{{c}^{2}-Re{D}_{\xi}{D}_{\eta}\left[\psi \right]\right\}$ (53)

$f=1+y\left\{{d}^{2}+Re{D}_{\xi}{D}_{\eta}\left[\psi \right]\right\}$ (54)

Figure 1. Two dimensional 9-point compact stencil used for discretization.

$g=\frac{-Re}{2}\left[\left\{{D}_{\eta}+y\text{\hspace{0.05em}}d\text{\hspace{0.05em}}{D}_{\eta}^{2}+x\text{\hspace{0.05em}}c\text{\hspace{0.05em}}{D}_{\xi}{D}_{\eta}+z\text{\hspace{0.05em}}{D}_{\xi}^{2}{D}_{\eta}\right\}\left[\psi \right]+y\text{\hspace{0.05em}}{\text{\pi}}^{2}{\text{e}}^{2\text{\pi}\xi}{D}_{\eta}\left[\omega \right]\right]$ (55)

$o=\frac{Re}{2}\left[\left\{{D}_{\xi}+x\text{\hspace{0.05em}}c\text{\hspace{0.05em}}{D}_{\xi}^{2}+y\text{\hspace{0.05em}}d\text{\hspace{0.05em}}{D}_{\xi}{D}_{\eta}+z\text{\hspace{0.05em}}{D}_{\xi}{D}_{\eta}^{2}\right\}\left[\psi \right]+x\text{\hspace{0.05em}}{\text{\pi}}^{2}{\text{e}}^{2\text{\pi}\xi}{D}_{\xi}\left[\omega \right]\right]$ (56)

$l=z\text{\hspace{0.05em}}c\text{\hspace{0.05em}}d+Re\left\{x\text{\hspace{0.05em}}{D}_{\xi}^{2}-y\text{\hspace{0.05em}}{D}_{\eta}^{2}\right\}\left[\psi \right]$ (57)

In addition to the above, we ill represent Equations ((42) and (43)) in fourth order accurate discrete form as below

$c=-Re\left\{\frac{1}{2}{D}_{\eta}\left[\psi \right]+y\text{\hspace{0.05em}}{D}_{\xi}^{2}{D}_{\eta}\left[\psi \right]+y\text{\hspace{0.05em}}{\text{\pi}}^{2}{\text{e}}^{2\text{\pi}\xi}{D}_{\eta}\left[\omega \right]\right\}$ (58)

$d=Re\left\{\frac{1}{2}{D}_{\xi}\left[\psi \right]+x\text{\hspace{0.05em}}{D}_{\xi}{D}_{\eta}^{2}\left[\psi \right]+x\text{\hspace{0.05em}}{\text{\pi}}^{2}{\text{e}}^{2\text{\pi}\xi}{D}_{\xi}\left[\omega \right]\right\}$ (59)

3.3. Discretization of Magnetic Streamfunction

Similar to the case of streamfunction equation, first the equation for magnetic field (19) is expressed in terms of magnetic stream function using Equations ((12) and (13)) and the velocities are replaced using (10) and (11), so that we get the magnetic streamfunction equation as given below

$\frac{{\partial}^{2}A}{\partial {\xi}^{2}}+\frac{{\partial}^{2}A}{\partial {\eta}^{2}}+C\frac{\partial A}{\partial \xi}+D\frac{\partial A}{\partial \eta}=0$ (60)

where the coefficients C and D are

$C=-\frac{{R}_{m}}{2}\left(\frac{\partial \psi}{\partial \eta}\right)$ (61)

$D=\frac{{R}_{m}}{2}\left(\frac{\partial \psi}{\partial \xi}\right)$ (62)

Substituting (21)-(24) in the magnetic streamfunction Equation (60), we get

${D}_{\xi}^{2}{A}_{i\mathrm{,}j}+{D}_{\eta}^{2}{A}_{i\mathrm{,}j}+{C}_{i\mathrm{,}j}{D}_{\xi}{A}_{i\mathrm{,}j}+{D}_{i\mathrm{,}j}{D}_{\eta}{A}_{i\mathrm{,}j}+{\zeta}_{i\mathrm{,}j}=0$ (63)

and the truncation error in the Equation (63) is

${\zeta}_{i\mathrm{,}j}=-{\left(C\frac{{h}^{2}}{6}\frac{{\partial}^{3}A}{\partial {\xi}^{3}}+\frac{{h}^{2}}{12}\frac{{\partial}^{4}A}{\partial {\xi}^{4}}+D\frac{{k}^{2}}{6}\frac{{\partial}^{3}A}{\partial {\eta}^{3}}+\frac{{k}^{2}}{12}\frac{{\partial}^{4}A}{\partial {\eta}^{4}}\right)}_{i\mathrm{,}j}+\mathcal{O}\left({h}^{4}\mathrm{,}{k}^{4}\right)$ (64)

It may be noted that the coefficients in the magnetic streamfunction equations are denoted by C and D while those in vorticity equation are denoted by c and d. Also, the coefficient D in the above equation may not confuse with the operators like ${D}_{\xi}$ because operators will always have a suffix $\xi $ or $\eta $ . Elimination of higher order derivatives of A in (64) is done by differentiating the magnetic streamfunction Equation (60) once and twice with respect to $\xi $ and $\eta $ respectively, to yield the following.

$\frac{{\partial}^{3}A}{\partial {\xi}^{3}}=-\frac{{\partial}^{3}A}{\partial \xi \partial {\eta}^{2}}-C\frac{{\partial}^{2}A}{\partial {\xi}^{2}}-D\frac{{\partial}^{2}A}{\partial \xi \partial \eta}-\frac{\partial C}{\partial \xi}\frac{\partial A}{\partial \xi}-\frac{\partial D}{\partial \xi}\frac{\partial A}{\partial \eta}$ (65)

$\begin{array}{c}\frac{{\partial}^{4}A}{\partial {\xi}^{4}}\mathrm{=}-\frac{{\partial}^{4}A}{\partial {\xi}^{2}\partial {\eta}^{2}}+C\frac{{\partial}^{3}A}{\partial \xi \partial {\eta}^{2}}-D\frac{{\partial}^{3}A}{\partial {\xi}^{2}\partial \eta}+\left({C}^{2}-2\frac{\partial C}{\partial \xi}\right)\frac{{\partial}^{2}A}{\partial {\xi}^{2}}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\left(CD-2\frac{\partial D}{\partial \xi}\right)\frac{{\partial}^{2}A}{\partial \xi \partial \eta}+\left(C\frac{\partial C}{\partial \xi}-\frac{{\partial}^{2}C}{\partial {\xi}^{2}}\right)\frac{\partial A}{\partial \xi}+\left(C\frac{\partial D}{\partial \xi}-\frac{{\partial}^{2}D}{\partial {\xi}^{2}}\right)\frac{\partial A}{\partial \eta}\end{array}$ (66)

$\frac{{\partial}^{3}A}{\partial {\eta}^{3}}=-\frac{{\partial}^{3}A}{\partial {\xi}^{2}\partial \eta}-C\frac{{\partial}^{2}A}{\partial \xi \partial \eta}-D\frac{{\partial}^{2}A}{\partial {\eta}^{2}}-\frac{\partial C}{\partial \eta}\frac{\partial A}{\partial \xi}-\frac{\partial D}{\partial \eta}\frac{\partial A}{\partial \eta}$ (67)

$\begin{array}{c}\frac{{\partial}^{4}A}{\partial {\eta}^{4}}\mathrm{=}-\frac{{\partial}^{4}A}{\partial {\xi}^{2}\partial {\eta}^{2}}-C\frac{{\partial}^{3}A}{\partial \xi \partial {\eta}^{2}}+D\frac{{\partial}^{3}A}{\partial {\xi}^{2}\partial \eta}+\left(CD-2\frac{\partial C}{\partial \eta}\right)\frac{{\partial}^{2}A}{\partial \xi \partial \eta}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\left({D}^{2}-2\frac{\partial D}{\partial \eta}\right)\frac{{\partial}^{2}A}{\partial {\eta}^{2}}+\left(D\frac{\partial C}{\partial \eta}-\frac{{\partial}^{2}C}{\partial {\eta}^{2}}\right)\frac{\partial A}{\partial \xi}+\left(D\frac{\partial D}{\partial \eta}-\frac{{\partial}^{2}D}{\partial {\eta}^{2}}\right)\frac{\partial A}{\partial \eta}\end{array}$ (68)

Now we first substitute (61)-(68) in Equation (63) and then apply D-operators from Equations (25)-(32) to arrive at the fourth order accurate discretized representation of the magnetic streamfunction Equation (60) as given below.

$\begin{array}{l}\left[-\alpha {D}_{\xi}^{2}-{\beta}^{\prime}{D}_{\eta}^{2}+\gamma {D}_{\xi}+\lambda {D}_{\eta}+\mu {D}_{\xi}{D}_{\eta}\right]A\\ -\left(x+y\right)\left[{D}_{\xi}^{2}{D}_{\eta}^{2}-C\text{\hspace{0.05em}}{D}_{\xi}{D}_{\eta}^{2}-D{D}_{\xi}^{2}{D}_{\eta}\right]A=0\end{array}$ (69)

where the coefficients $\alpha $ , ${\beta}^{\prime}$ , $\gamma $ , $\lambda $ , $\mu $ , $C$ and $D$ appearing above are given by

$\alpha =\left(1+x\text{\hspace{0.05em}}{C}^{2}\right)-\left[{R}_{m}x\text{\hspace{0.05em}}{D}_{\xi}{D}_{\eta}\right]\psi $ (70)

${\beta}^{\prime}=\left(1+y{D}^{2}\right)-\left[{R}_{m}y{D}_{\eta}^{2}\right]\psi $ (71)

$\begin{array}{l}\gamma =C+\frac{{R}_{m}}{2}x\left[{D}_{\xi}^{2}{D}_{\eta}-C\text{\hspace{0.05em}}{D}_{\xi}{D}_{\eta}\right]\psi -\frac{{R}_{m}}{2}y\left[{D}_{\xi}^{2}{D}_{\eta}-D\text{\hspace{0.05em}}{D}_{\eta}^{2}\right]\psi \\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}-\left[\frac{{R}_{m}}{2}y{\text{\pi}}^{2}{\text{e}}^{2\text{\pi}\xi}{d}_{\eta}\right]\omega \end{array}$ (72)

$\begin{array}{c}\lambda =D+\frac{{R}_{m}}{2}x\left[{D}_{\xi}{D}_{\eta}^{2}+C\text{\hspace{0.05em}}{D}_{\xi}^{2}\right]\psi -\frac{{R}_{m}}{2}y\left[{D}_{\xi}{D}_{\eta}^{2}+D\text{\hspace{0.05em}}{D}_{\xi}{D}_{\eta}\right]\psi \\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\left[\frac{{R}_{m}}{2}x{\text{\pi}}^{2}{\text{e}}^{2\text{\pi}\xi}\left(2\text{\pi}+{D}_{\xi}\right)\right]\omega \end{array}$ (73)

$\mu ={R}_{m}\left[-x{D}_{\xi}^{2}+y{D}_{\eta}^{2}\right]\psi -\left(x+y\right)C\text{\hspace{0.05em}}D$ (74)

In addition to the above, we ill represent Equations ((61) and (62)) in fourth order accurate discrete form as below

$C=\frac{{R}_{m}}{2}\left[{D}_{\eta}+2y{D}_{\xi}^{2}{D}_{\eta}\right]\psi +\left[{R}_{m}y{\text{\pi}}^{2}{\text{e}}^{2\text{\pi}\xi}{D}_{\eta}\right]\omega $ (75)

$D=-\frac{{R}_{m}}{2}\left[{D}_{\xi}+2x{D}_{\xi}{D}_{\eta}^{2}\right]\psi -\left[{R}_{m}x{\text{\pi}}^{2}{\text{e}}^{2\text{\pi}\xi}\left(2\pi +{D}_{\xi}\right)\right]\omega $ (76)

3.4. Discretization of Energy Equation

In Equation (20) velocities are eliminated using (10) and (11) then, $r={\text{e}}^{\text{\pi}\xi}$ and $\theta =\text{\pi}\eta $ transformations are employed. This leads to the following equation for temperature.

$\frac{{\partial}^{2}T}{\partial {\xi}^{2}}+\frac{{\partial}^{2}T}{\partial {\eta}^{2}}+{c}^{\prime}\frac{\partial T}{\partial \xi}+{d}^{\prime}\frac{\partial T}{\partial \eta}=0$ (77)

where the coefficients, ${c}^{\prime}$ and ${d}^{\prime}$ are given by

${c}^{\prime}=-\frac{RePr}{2}\left(\frac{\partial \psi}{\partial \eta}\right)$ (78)

${d}^{\prime}=\frac{RePr}{2}\left(\frac{\partial \psi}{\partial \xi}\right)$ (79)

Equations (21)-(24) are applied on the above equation so that we have

${D}_{\xi}^{2}{T}_{i\mathrm{,}j}+{D}_{\eta}^{2}{T}_{i\mathrm{,}j}+{{c}^{\prime}}_{i\mathrm{,}j}{D}_{\xi}{T}_{i\mathrm{,}j}+{{d}^{\prime}}_{i\mathrm{,}j}{D}_{\eta}{T}_{i\mathrm{,}j}+{\gamma}_{i\mathrm{,}j}=0$ (80)

The truncation error in the Equation (80) is

${\gamma}_{i\mathrm{,}j}=-\left({c}^{\prime}\frac{{h}^{2}}{6}\frac{{\partial}^{3}T}{\partial {\xi}^{3}}+\frac{{h}^{2}}{12}\frac{{\partial}^{4}T}{\partial {\xi}^{4}}+{d}^{\prime}\frac{{k}^{2}}{6}\frac{{\partial}^{3}T}{\partial {\eta}^{3}}+\frac{{k}^{2}}{12}\frac{{\partial}^{4}T}{\partial {\eta}^{4}}\right)+\mathcal{O}\left({h}^{4}\mathrm{,}{k}^{4}\right)$ (81)

In order to higher derivatives of T present in the truncation error terms, we differentiate the energy Equation (77) once and twice with respect to $\xi $ and $\eta $ respectively and we get the following

$\frac{{\partial}^{3}T}{\partial {\xi}^{3}}=-\frac{{\partial}^{3}T}{\partial \xi \partial {\eta}^{2}}-{c}^{\prime}\frac{{\partial}^{2}T}{\partial {\xi}^{2}}-{d}^{\prime}\frac{{\partial}^{2}T}{\partial \xi \partial \eta}-\frac{\partial {c}^{\prime}}{\partial \xi}\frac{\partial T}{\partial \xi}-\frac{\partial {d}^{\prime}}{\partial \xi}\frac{\partial T}{\partial \eta}$ (82)

$\begin{array}{c}\frac{{\partial}^{4}T}{\partial {\xi}^{4}}=-\frac{{\partial}^{4}T}{\partial {\xi}^{2}\partial {\eta}^{2}}+{c}^{\prime}\frac{{\partial}^{3}T}{\partial \xi \partial {\eta}^{2}}-{d}^{\prime}\frac{{\partial}^{3}T}{\partial {\xi}^{2}\partial \eta}+\left({{c}^{\prime}}^{2}-2\frac{\partial {c}^{\prime}}{\partial \xi}\right)\frac{{\partial}^{2}T}{\partial {\xi}^{2}}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\left({c}^{\prime}{d}^{\prime}-2\frac{\partial {d}^{\prime}}{\partial \xi}\right)\frac{{\partial}^{2}T}{\partial \xi \partial \eta}+\left({c}^{\prime}\frac{\partial {c}^{\prime}}{\partial \xi}-\frac{{\partial}^{2}{c}^{\prime}}{\partial {\xi}^{2}}\right)\frac{\partial T}{\partial \xi}+\left({c}^{\prime}\frac{\partial {d}^{\prime}}{\partial \xi}-\frac{{\partial}^{2}{d}^{\prime}}{\partial {\xi}^{2}}\right)\frac{\partial T}{\partial \eta}\end{array}$ (83)

$\frac{{\partial}^{3}T}{\partial {\eta}^{3}}=-\frac{{\partial}^{3}T}{\partial {\xi}^{2}\partial \eta}-{c}^{\prime}\frac{{\partial}^{2}T}{\partial \xi \partial \eta}-{d}^{\prime}\frac{{\partial}^{2}T}{\partial {\eta}^{2}}-\frac{\partial {c}^{\prime}}{\partial \eta}\frac{\partial T}{\partial \xi}-\frac{\partial {d}^{\prime}}{\partial \eta}\frac{\partial T}{\partial \eta}$ (84)

$\begin{array}{c}\frac{{\partial}^{4}T}{\partial {\eta}^{4}}=-\frac{{\partial}^{4}T}{\partial {\xi}^{2}\partial {\eta}^{2}}-{c}^{\prime}\frac{{\partial}^{3}T}{\partial \xi \partial {\eta}^{2}}+{d}^{\prime}\frac{{\partial}^{3}T}{\partial {\xi}^{2}\partial \eta}+\left({{d}^{\prime}}^{2}-2\frac{\partial {d}^{\prime}}{\partial \eta}\right)\frac{{\partial}^{2}T}{\partial {\eta}^{2}}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\left({c}^{\prime}{d}^{\prime}-2\frac{\partial {c}^{\prime}}{\partial \eta}\right)\frac{{\partial}^{2}T}{\partial \xi \partial \eta}+\left({d}^{\prime}\frac{\partial {c}^{\prime}}{\partial \eta}-\frac{{\partial}^{2}{c}^{\prime}}{\partial {\eta}^{2}}\right)\frac{\partial T}{\partial \xi}+\left({d}^{\prime}\frac{\partial {d}^{\prime}}{\partial \eta}-\frac{{\partial}^{2}{d}^{\prime}}{\partial {\eta}^{2}}\right)\frac{\partial T}{\partial \eta}\end{array}$ (85)

Now substituting Equations (78)-(85) in the truncation error terms of Equation (77), and applying operator D to the derivatives, we will get the fourth order discretized form of energy equation, as given below

$\begin{array}{l}\left[-\varsigma {D}_{\xi}^{2}-\rho {D}_{\eta}^{2}+\varpi {D}_{\xi}+\vartheta {D}_{\eta}+\epsilon {D}_{\xi}{D}_{\eta}\right]T\\ -\left(x+y\right)\left[{D}_{\xi}^{2}{D}_{\eta}^{2}-C\text{\hspace{0.05em}}{D}_{\xi}{D}_{\eta}^{2}-D{D}_{\xi}^{2}{D}_{\eta}\right]T=0\end{array}$ (86)

where the coefficients $\varsigma $ , $\rho $ , $\varpi $ , $\vartheta $ , $\epsilon $ , ${c}^{\prime}$ and ${d}^{\prime}$ appearing above are given by

$\varsigma =\left(1+x\text{\hspace{0.05em}}{{c}^{\prime}}^{2}\right)-\left[RePr\text{\hspace{0.05em}}x\text{\hspace{0.05em}}{D}_{\xi}{D}_{\eta}\right]\psi $ (87)

$\rho =\left(1+y{{d}^{\prime}}^{2}\right)-\left[RePr\text{\hspace{0.05em}}y\text{\hspace{0.05em}}{D}_{\eta}^{2}\right]\psi $ (88)

$\begin{array}{c}\varpi ={c}^{\prime}+\frac{RePr}{2}x\left[{D}_{\xi}^{2}{D}_{\eta}-{c}^{\prime}{D}_{\xi}{D}_{\eta}\right]\psi -\frac{RePr}{2}y\left[{D}_{\xi}^{2}{D}_{\eta}-{d}^{\prime}{D}_{\eta}^{2}\right]\psi \\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-\left[\frac{RePr}{2}y{\text{\pi}}^{2}{\text{e}}^{2\text{\pi}\xi}{d}_{\eta}\right]\omega \end{array}$ (89)

$\begin{array}{c}\vartheta ={d}^{\prime}+\frac{RePr}{2}x\left[{D}_{\xi}{D}_{\eta}^{2}+{c}^{\prime}{D}_{\xi}^{2}\right]\psi -\frac{RePr}{2}y\left[{D}_{\xi}{D}_{\eta}^{2}+{d}^{\prime}{D}_{\xi}{D}_{\eta}\right]\psi \\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\left[\frac{RePr}{2}x{\text{\pi}}^{2}{\text{e}}^{2\text{\pi}\xi}\left(2\text{\pi}+{D}_{\xi}\right)\right]\omega \end{array}$ (90)

$\epsilon =RePr\left[-x{D}_{\xi}^{2}+y{D}_{\eta}^{2}\right]\psi -\left(x+y\right){c}^{\prime}{d}^{\prime}$ (91)

In addition to the above, we will represent Equations ((78) and (79)) in fourth order accurate discrete form as below

${c}^{\prime}=\frac{RePr}{2}\left[{D}_{\eta}+2y{D}_{\xi}^{2}{D}_{\eta}\right]\psi +\left[RePr\text{\hspace{0.05em}}y{\text{\pi}}^{2}{\text{e}}^{2\text{\pi}\xi}{D}_{\eta}\right]\omega $ (92)

${d}^{\prime}=-\frac{RePr}{2}\left[{D}_{\xi}+2x{D}_{\xi}{D}_{\eta}^{2}\right]\psi -\left[RePr\text{\hspace{0.05em}}x{\text{\pi}}^{2}{\text{e}}^{2\text{\pi}\xi}\left(2\text{\pi}+{D}_{\xi}\right)\right]\omega $ (93)

3.5. Boundary Conditions and Implementation

We evaluate the performance of the above discretization scheme by running numerical experiments for the MHD flow around a cylinder, since it is one of the classic problems in MHD flow where the flow features have been studied by many authors. Hence we can compare our high accuracy results and discuss on the performance of proposed discretization scheme. In order to cluster more grid points near the surface of the cylinder, an exponential transformation is employed. It is also assumed that the magnetic field will not penetrate across the boundary of the cylinder. The boundary conditions used are the following. On the surface of the cylinder, $\psi =0$ , $\partial \psi /\partial r=0$ , $\omega ={\partial}^{2}\psi /\partial {r}^{2}$ and $A=0$ . Far away from the cylinder which is denoted by $\xi \to \infty $ , the imposed conditions are $\psi ~rsin\left(\theta \right)$ , $\omega \to 0$ , $A\to rsin\left(\theta \right)$ . The axis of symmetry is given by $\eta =0$ and $\eta =1$ , on which lines, the boundary conditions imposed are $\psi =0$ , $\omega =0$ , and $A=0$ . For solving the energy equation, we impose a uniform temperature $T=1$ on the cylinder surface and $T=0$ at far-field distances. On the axis of symmetry, $\partial T/\partial \eta =0$ is imposed. The derivatives appearing in the boundary conditions for streamfunction vorticity and temperature are

$\frac{\partial \psi}{\partial \xi}=0$ (94)

$\frac{{\partial}^{2}\psi}{\partial {\xi}^{2}}=-{\text{\pi}}^{2}{\omega}^{2}$ (95)

$\frac{\partial T}{\partial \eta}=0$ (96)

The above derivatives are replaced with suitable one-sided fourth order accurate finite differences as given below

${\psi}_{1,j}=\frac{\left(48{\psi}_{2,j}-36{\psi}_{3,j}+16{\psi}_{4,j}-3{\psi}_{5,j}\right)}{25}$ (97)

${\psi}_{n+1,j}=\frac{-\left(48{\psi}_{n,j}-36{\psi}_{n-1,j}+16{\psi}_{n-2,j}-3{\psi}_{n-3,j}\right)}{25}$ (98)

${T}_{i,1}=\frac{\left(48{T}_{i,2}-36{T}_{i,3}+16{T}_{i,4}-3{T}_{i,5}\right)}{25}$ (99)

${T}_{i,m+1}=\frac{-\left(48{T}_{i,m}-36{T}_{i,m-1}+16{T}_{i,m-2}-3{T}_{i,m-3}\right)}{25}$ (100)

The vorticity boundary condition used on the surface of the cylinder is the Briley’s formula

${\omega}_{1,j}=\frac{-\left(108{\psi}_{2,j}-27{\psi}_{3,j}+4{\psi}_{4,j}\right)}{18{h}^{2}}$ (101)

The algebraic system of Equations ((40), (51), (69) and (86)) exhibits diagonal dominance and hence the Gauss-Seidel iterative procedure is used to solve the discretized coupled equations. The convergence criterion is set as the norm of dynamic residuals and computations are terminated if the norm is less than 10^{−8}. An additional difficulty faced due to the high accuracy numerical scheme is that, at every iteration, before the Equation (51) can be solved, the values of
$G$ ,
${D}_{\xi}\left[G\right]$ ,
${D}_{\eta}\left[G\right]$ ,
${D}_{\xi}^{2}\left[G\right]$ ,
${D}_{\eta}^{2}\left[G\right]$ and all coefficients therein, Equations (53)-(58) should be evaluated at the required grid points. This same is true while solving for magnetic stream function, A and the temperature T equations. It is also observed that the presence of several nonlinearities on the right side of (51) which makes the numerical solution time consuming and parameter restricted. This means a higher computational cost in terms of memory is required, especially at high resolution grids like
$512\times 512$ . In addition, it takes more time for convergence when compared to the use of traditional unwinding schemes. The lift and drag are the two forces acting on the cylinder due to fluid flow. Because of the two counter rotating symmetric vortices behind the circular cylinder, the lift forces get canceled. The drag coefficients (pressure, viscous and total) are calculated using the relations

${C}_{P}=\frac{4\text{\pi}}{Re}{\displaystyle {\int}_{0}^{1}}{\left(\frac{\partial \omega}{\partial \xi}\right)}_{\xi =0}sin\left(\text{\pi}\eta \right)\text{d}\eta $ (102)

${C}_{V}=-\frac{4\text{\pi}}{Re}{\displaystyle {\int}_{0}^{1}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{\omega}_{\xi =0}sin\left(\text{\pi}\eta \right)\text{d}\eta $ (103)

${C}_{D}={C}_{P}+{C}_{V}$ (104)

The heat flux $q\text{\hspace{0.05em}}\left(\theta \right)$ , local Nusselt number $Nu$ and mean Nusselt number $Nm$ are calculated using

$q\text{\hspace{0.05em}}\left(\theta \right)=\frac{k\left({T}_{s}-{T}_{\infty}\right)}{a}{\left(\frac{\partial T}{\partial \xi}\right)}_{\xi =0}$ (105)

$Nu\text{\hspace{0.05em}}\left(\theta \right)=\frac{2aq\left(\theta \right)}{k\left({T}_{s}-{T}_{\infty}\right)}=\frac{-2}{\text{\pi}}{\left(\frac{\partial T}{\partial \xi}\right)}_{\xi =0}$ (106)

$Nm={\displaystyle {\int}_{0}^{1}}{\left(\frac{\partial T}{\partial \xi}\right)}_{\xi =0}\text{d}\eta $ (107)

In the above, all derivatives are approximated with fourth order accurate central differences. However, for points near the boundary, appropriate forward or backward fourth order finite difference expressions are used. Since fine grid solutions are considered, the integral in the above is carried out using the Simpson’s rule.

4. Code Validation

We have computed the fourth order accurate results for $Re=40$ with magnetic Reynolds number ranging $0\le {R}_{m}\le 2$ for different values of the magnetic field strengths $\beta $ . The values of Prandtl number used in the present study are $Pr=0.065$ , 0.73, 5 and 7 which will typically represent liquids such as liquid lithium, air, aqueous KOH solution or some refrigerants and saline water. The numerical results are presented using the one obtained from the finest grid $256\times 256$ , while the coarser grids employed are $128\times 128$ and $64\times 64$ . The convergence criterion used in this study is that the norm of the dynamic residuals $\epsilon <{10}^{-8}$ . The far-field distance or the computational domain size is chosen as 41 times the radius of cylinder, which is a fairly a large space to acquire accurate results. Table 1 lists the values of drag coefficient for different

Table 1. Grid independence test: Drag coefficients obtained in 64 × 64, 128 × 128 and 256 × 256 grids. The results show that the computation in 128 × 128 grid is optimum.

Table 2. Comparison of flow parameters for Re = 40 with available data in literature for zero magnetic field case, that is for $\beta =0.$ Here, ${l}_{S}$ is the separation length measured from the center of the circular cylinder.

$Re$ and for three different grids in presence of magnetic field. From the table, it is noted that $128\times 128$ grid is found to be optimum. Further refining of grids, that is employing $512\times 512$ , etc. will not give more accuracy because the accuracy of the numerical scheme is decided by method of discretization and not by the fineness of grid. Table 2 shows the comparison of flow parameters obtained using our numerical scheme with the data available in the literature [24] [26] [27] [36] [39] [40] [42] - [60] . For the zero magnetic field $\left(\beta =0\right)$ case, good agreement of results can be seen wherein. However, we observe a consistent and slightly higher value for the drag coefficient obtained in the

Table 3. Comparison of results for Mean Nusselt number (Nm) in the absence of external magnetic field for Re = 40 and Pr = 0.73.

present work which is attributed to the higher order numerical scheme employed in contrast to traditional first or second order accurate data in the literature. Further the data for mean Nusselt number ${N}_{m}$ is also compared against other reported data [28] [58] - [72] and it is tabulated in Table 3. It is observed that out of the data shown in table, 71% of researchers have reported the value of mean Nusselt number close to our value of 3.29, whereas others have reported a higher value close to 3.48. The only value which departs from the rest is that reported in [68] .

5. Results and Discussion

The response of separation length $\left({l}_{s}\right)$ and the separation angle $\left({\theta}_{s}\right)$ of the recirculation bubble with the Alfvén number $\left(\beta \right)$ is presented in Figure 2. The separation length reduces strongly with small magnetic fields given by $\beta <1$ . Further increase of $\beta $ suppresses the separation only to a little extent. The ratio of inertial to Lorentz force is considered as friction parameter given by $Re/Ha$ is also equal to $\sqrt{Re/N}$ . The Hartmann layer is damped by the inverse friction parameter [38] and hence wake suppression is observed. From the plot of separation angle versus $\beta $ (Figure 2) it is noted that the flow of fluid with higher conductivity can be easily controlled with a weak magnetic field meaning that if the conductivity of the fluid is low, then a relatively strong magnetic field

Figure 2. Comparison of the separation length ${l}_{s}$ and separation angle ${\theta}_{s}$ computed using the present fourth order scheme with second order scheme.

is required for the same level of wake suppression. The Lorentz force will directly change the fluid velocity and according to Equation (2) and hence the vorticity of the fluid changes due to applied magnetic field. Together with a gradual decrease in the length of the standing vortex near the cylinder, a growing region of positive vortices are seen when the magnetic field strength is increased (not shown) and this is also predicted theoretically [73] . The positive contours of $\omega $ are not present when $N=0$ . The region of positive vorticity gets directed towards the center and extends in the upstream with increasing magnetic field or with increasing the conductivity of the fluid. This pushes the negative vorticity region be form at far distance [17] . Since the induced magnetic field is taken into consideration in the formulation of the problem, the magnetic stream function A is computed according to Equation (60) during the solution procedure, which is substituted in Equations ((12) and (13)) to get the effective or total magnetic field components. The contours of the lines are presented in Figure 3 for ${R}_{m}=1$ with varying $\beta $ . From this contour plot, we observer that the magnetic field is no longer uniform in the entire flow region. This is due to the current density which got generated in the fluid by virtue of its motion and, in turn, the applied magnetic field is modified. In regions near $\theta ={90}^{\circ}$ , the induced magnetic field is added to the applied field, while in the other upstream and in the wake regions, the induced magnetic field has opposed the applied field resulting in a lower magnetic field. Also, as $\beta $ increases, we observe more changes to the resultant magnetic field.

The dependence of pressure drag coefficient ${C}_{P}$ , viscous drag coefficient ${C}_{V}$ , and total drag coefficient ${C}_{D}$ on $\beta $ and ${R}_{m}$ are shown in Figure 4. The bottom-most plot of this graph show results from the present fourth order scheme. The top and middle plots show the comparison of results obtained from the present work with that obtained using a second order accurate scheme. Table 4 shows a detailed comparison of drag coefficients and flow separation parameters. For small values of magnetic Reynolds number ${R}_{m}$ the total drag coefficient decreases when $\beta \approx 0.5$ and then the drag coefficient increases for higher $\beta $ . This is in contrast to the result obtained using second order method where it is predicted that the drag coefficient minimum takes place when $\beta \simeq 1$ , which are noted from the red color data points of top and middle plots. The increase of drag coefficient following the ${\left(N/Re\right)}^{1/2}$ power-law with N is predicted theoretically [17] and the existence of a minimum drag coefficient for certain magnetic field is also reported [74] .

Figure 3. Contours of resultant magnetic field for $Re=40$ , ${R}_{m}=1$ . Here $\beta =1$ , 2 and 3 (top to bottom).

Figure 4. Comparison of second order and fourth order accurate calculation of total drag coefficient ${C}_{D}$ (top and middle). It is noted that for higher values of ${R}_{m}$ , the fourth order results show a larger deviation from the second order results. The bottom one shows the effect of magnetic Reynolds number on ${C}_{D}$ .

The forced convective heat transfer in the fluid flow is studied for $Re=40$ with $0\le {R}_{m}\le 2$ , $0\le \beta \le 50$ for Prandtl numbers 0.065, 0.73, 5, 7 which may represent typical fluids like liquid metals, air, KOH and water respectively. The surface Nusselt number on the surface of the cylinder and its effect on the magnetic field is shown in Figure 5 (bottom plot). In the absence of external magnetic field $Nu$ decreases on the surface of the cylinder until the separation point (green circles in the plot). Beyond this angle, an enhancement of heat transfer is noted which is due to the recirculation region. When the magnetic field is increased the minimum of the Nu in the plot is shifted towards $\theta =0$

Table 4. Drag coefficients and flow separation parameters for Re = 40 and its comparison with second order results.

Figure 5. Angular variation of Nusselt number and its influence on magnetic Reynolds number and the strength of magnetic field. Here, $Re=40$ and $Pr=0.73$ (top) and $Pr=7$ (bottom).

which is the rear stagnation point (red, black and blue lines). That is, along with suppression of flow separation, the enhancement of heat transfer is also suppressed. These effects proportionally increase for fluid with higher Prandtl numbers. The application of external magnetic field causes the thickness of the viscous boundary layer to increase. The maximum heat transfer region lies in the front stagnation point. At this point the magnetic field casus a decrease in Nu. The effect of electrical conductivity on the Nu is shown in Figure 5 (top plot). By considering the low magnetic case $\left(\beta =0.5\right)$ , in the upstream region, the heat transfer is higher for a fluid with lower conductivity. However, the reverse in true in the downstream region. This kind of non-monotonic behavior can be nullified, if the same experiment is repeated for a higher magnetic field, say $\beta =1$ . These effects are due to the change in radial and angular velocities of the fluid by the application of external magnetic field. Figure 6 shows the variations of mean Nusselt number with $\beta $ . For weak magnetic field strengths, the mean Nusselt number reduces when compared to the zero magnetic field case. This happens until a critical field value of $\beta =0.5$ [Figure 6 (bottom)]. For $\beta =0.65$ , the mean Nusselt number converges to a single value ${N}_{m}=2.6$ , independent of electrical conductivity of the fluid. The critical field divides the plot into two regions: For $\beta <0.4$ , the mean Nu increases with conductivity of the fluid. For $\beta >0.65$ the higher conductivity of the fluid aids to increase the heat transfer. Thus the mean Nusselt number ${N}_{m}$ depends non-monotonically with the electrical conductivity ${R}_{m}$ of the fluid. In the region $0.4\le \beta \le 0.65$

Figure 6. Effect of magnetic field and magnetic Reynolds number on the mean Nusselt number for $Re=40$ and ${R}_{m}=1$ .

Figure 7. Variation of temperature with radial distance, $\xi $ along $\theta =\text{\pi}/4$ line and its influence on magnetic field strength (top) and Prandtl numbers (bottom). Here $Re=40$ and ${R}_{m}=1$ .

non-monotonic behavior takes place.

The variation of temperature with radial distance from the surface of the cylinder is plotted in the Figure 7 along a specific direction $\theta =\text{\pi}/4$ . This direction is chosen because along this line the fluid flow will be mixing with the recirculation zone. The top plot shows its dependence with magnetic field. The slope of the plot at any point will be the temperature gradient along the radial direction. It is noted that the temperature gradient significantly increases with the increase of magnetic field only near the cylinder body. The effect of $Pr$ on the temperature and hence the thermal gradient is shown in Figure 7 (bottom). Steep thermal gradients occur in fluids with higher $Pr$ leading to higher heat transfer. The thermal boundary layer thickness for the present problem is taken as the thickness of the fluid over which 99% of the cylinder body temperature is lost. It is found that for $Pr=0.73$ , the magnetic field reduces the thermal boundary layer thickness from 5.03 to 3.6. It may be noted that the x-axis of the plot $\xi $ denotes the distance from the center of cylinder in multiples of its radius. For a given magnetic field, say $\beta =2$ , the thermal boundary layer thickness reduces drastically with increase of $Pr$ . Figure 8 shows the effect of Prandtl number on temperature isotherms in the absence of external magnetic field. Since a fluid with higher $Pr$ has the tendency to take away more heat, plume formation takes place around the separation point and maximum heat transfer takes in that region. In addition, the higher the $Pr$ of the fluid, the

Figure 8. Effect of $Pr$ on the thermal field in the absence of magnetic field. Here, $Re=40$ and $Pr=0.065$ , 0.73, 5 and 7 respectively (top to bottom).

larger is the density of isotherm lines in the upstream region, which means that the thermal boundary layer thickness is reduced. This is not the case when the magnetic field strength is varied for a fixed $Pr$ , where the thermal boundary layer thickness does not significantly change.

6. Demonstration of Numerical Accuracy

In order to evaluate the accuracy of our numerical results, we have taken the

Figure 9. A plot of divided difference of drag coefficient $\text{d}\left({C}_{D}\right)/\text{d}h$ versus the grid step size h.

drag coefficient values ${C}_{D}$ from different grids of grid step size h. Then, its divided difference with respect to h is calculated as

$\frac{\Delta \left({C}_{D}\right)}{dh}=\frac{{\left({C}_{D}\right)}_{{h}_{1}}-{\left({C}_{D}\right)}_{{h}_{2}}}{{h}_{1}-{h}_{2}}=y$ (108)

A plot between y and the grid step size h is made on a log-log scale which is shown in Figure 9. The three points a, b and c in the figure correspond to the values obtained in four different grids. It is found that the result follow a straight line behavior with a slope equal to 3. This is true for zero magnetic field case as well as for flow with any ${R}_{m}$ and $\beta $ considered in this study. In particular, the plot shows the $N=0$ and $N=0.5$ cases along with the one available in the literature [40] . Since the quantity in y-axis itself is equivalent to one slope (because of divided differences) the overall slope of the line is equal to four, indicating that the numerical results obtained in the present study have the desired fourth order accuracy.

7. Conclusion

In this paper, we have studied the problem of electrically conducting flow over a circular cylinder in the presence a steady magnetic field using a high accuracy finite difference scheme. When compared to the results obtained using second order mixed with upwind schemes, it is found that relatively lower magnetic field is sufficient for suppressing the flow separation. In the present study, it is found that the suppression of flow separation depends on both conductivity of the fluid and the strength of magnetic field, while the study using the low ${R}_{m}$ approximation reveals that the separation is nearly independent of conductivity of the fluid. For small values of magnetic Reynolds number ${R}_{m}$ the total drag coefficient decreases when $\beta \approx 0.5$ and then the drag coefficient increases for higher values of $\beta $ . In the absence of external magnetic field, the minimum heat transfer takes place at the separation point. When magnetic field is incre- ased, the minimum of the Nusselt number $Nu$ is shifted towards the rear stagnation point. The mean Nusselt number ( ${N}_{m}$ ) depends non-monotonically with the electrical conductivity ( ${R}_{m}$ ) of the fluid in the region $0.4\le \beta \le 0.65$ . The degradation of heat transfer by the application of external magnetic field up to the values of $\beta \approx 0.5$ , further there is an increase of heat transfer is observed. The temperature gradient significantly increases with the increase of magnetic field only near the cylinder body. High thermal gradients occur in fluids with higher $Pr$ leading to higher heat transfer. The overall accuracy of the derived numerical scheme is estimated to be four using divided difference principle.

Acknowledgements

One of the authors, R.Sivakumar thank UGC for supporting this work vide major research project grant letter F. No. 37-312/2009 (SR) dated January 12, 2010, and also DST for their FIST funding vide order SR/FST/PSII-021/2009 dated August 13, 2010 towards establishing a basic computing cluster facility. Finally Abin Rejeesh would like to thank DST for providing INSPIRE fellowship vide letter No. DST/INSPIRE Fellowship/2011/[361] dated January 16, 2012.

Nomenclature

$A$ = Magnetic stream function.

$a\mathrm{,}d$ = Radius and diameter of cylinder.

$c$ = Specific heat capacity of fluid.

${C}_{P}\mathrm{,}{C}_{V}\mathrm{,}{C}_{D}$ = Pressure, viscous and total drag coefficient.

$E$ = Electric field.

${\rm H}\mathrm{,}{H}_{\infty}$ = Induced Magnetic field and applied far field magnetic field.

${h}_{r}\mathrm{,}{h}_{\theta}$ = Radial and transverse magnetic field components.

$J$ = Electric current density.

${l}_{s}$ = Separation length of the wake.

$N$ = Interaction parameter or Stuart number = $a\sigma {H}_{\infty}^{2}/\left(\rho {U}_{\infty}\right)$ .

$Nu$ = Local Nusselt number:

${N}_{m}$ = Average Nusselt number

$p$ = Dimensionless pressure.

$Pr$ = Prandtl number = $\nu /\kappa $ .

$q$ = heat flux.

${q}_{r}\mathrm{,}{q}_{\theta}$ = Radial and transverse velocity components.

$Re$ = Reynolds number = $2a{U}_{\infty}/\nu $

${R}_{m}$ = Magnetic Reynolds number = $2a{U}_{\infty}{\mu}_{m}\sigma $

$T$ = Dimensionless temperature.

${T}_{s}$ = Surface temperature of the cylinder.

${T}_{\infty}$ = Far field free stream temperature.

${U}_{\infty}$ = Uniform free stream velocity in the far field:

$\alpha $ = Thermal expansion coefficient of the fluid:

$\beta $ = Alfven number.

${\delta}_{T}$ = Thermal boundary layer thickness:

$\kappa $ = Thermal diffusivity of the fluid:

${\mu}_{m}$ = Magnetic permeability of the fluid.

$\mu $ = Dynamic viscosity of the fluid:

$\nu $ = Kinematic viscosity of the fluid = $\mu /\rho $

$\rho $ = Density of the fluid:

$\sigma $ = Electrical conductivity of the fluid:

${\theta}_{s}$ = Separation angle of the wake.

$\left(r\mathrm{,}\theta \right)$ = 2D cylindrical polar coordinates:

$\left(\xi \mathrm{,}\eta \right)$ = Coordinates after $r={\text{e}}^{\text{\pi}\xi}$ , and $\theta =\text{\pi}\eta $ transformation.

$\left(\psi \mathrm{,}\omega \right)$ = Streamfunction and vorticity.

References

[1] Dupret, F., Nicodeme, P., Ryckmans, Y., Wouters, P. and Crochet, M. (1990) Global Modelling of Heat Transfer in Crystal Growth Furnaces. International Journal of Heat and Mass Transfer, 33, 1849-1871.

https://doi.org/10.1016/0017-9310(90)90218-J

[2] Hussam, W.K., Thompson, M.C. and Sheard, G.J. (2011) Dynamics and Heat Transfer in a Quasi-Two-Dimensional MHD Flow Past a Circular Cylinder in a Duct at High Hartmann Number. International Journal of Heat and Mass Transfer, 54, 1091-1100.

https://doi.org/10.1016/j.ijheatmasstransfer.2010.11.013

[3] Mbeledogu, I.U. and Ogulu, A. (2007) Heat and Mass Transfer of an Unsteady MHD Natural Convection Flow of a Rotating Fluid Past a Vertical Porous Flat Plate in the Presence of Radiative Heat Transfer. International Journal of Heat and Mass Transfer, 50, 1902-1908.

https://doi.org/10.1016/j.ijheatmasstransfer.2006.10.016

[4] Das, S.S., Satapathy, A., Das, J.K. and Panda, J.P. (2009) Mass Transfer Effects on MHD Flow and Heat Transfer Past a Vertical Porous Plate Through a Porous Medium under Oscillatory Suction and Heat Source. International Journal of Heat and Mass Transfer, 52, 5962-5969.

https://doi.org/10.1016/j.ijheatmasstransfer.2009.04.038

[5] Sparrow, E.M. and Cess, R.D. (1961) The Effect of a Magnetic Field on Free Convection Heat Transfer. International Journal of Heat and Mass Transfer, 3, 267-274.

https://doi.org/10.1016/0017-9310(61)90042-4

[6] Hossain, M.A. (1992) Viscous and Joule Heating Effects on MHD-Free Convection Flow with Variable Plate Temperature. International Journal of Heat and Mass Transfer, 35, 3485-3487.

https://doi.org/10.1016/0017-9310(92)90234-J

[7] Seddeek, M.A. (2002) Effects of Radiation and Variable Viscosity on a MHD Free Convection Flow Past a Semi-Infinite Flat Plate with an Aligned Magnetic Field in the Case of Unsteady Flow. International Journal of Heat and Mass Transfer, 45, 931-935.

https://doi.org/10.1016/S0017-9310(01)00189-2

[8] Tsai, R. and Huang, J.S. (2009) Heat and Mass Transfer for Soret and Dufour’s Effects on Hiemenz Flow through Porous Medium onto a Stretching Surface. International Journal of Heat and Mass Transfer, 52, 2399-2406.

https://doi.org/10.1016/j.ijheatmasstransfer.2008.10.017

[9] Khan, W.A. and Pop, I. (2010) Boundary-Layer Flow of a Nanofluid past a Stretching Sheet. International Journal of Heat and Mass Transfer, 53, 2477-2483.

https://doi.org/10.1016/j.ijheatmasstransfer.2010.01.032

[10] Noor, N.F.M., Abbasbandy, S. and Hashim, I. (2010) Heat and Mass Transfer of Thermophoresis MHD Flow Over an Inclined Radiative Isothermal Permeable Surface in the Presence of Heat Source/Sink. International Journal of Heat and Mass Transfer, 55, 2122-2128.

https://doi.org/10.1016/j.ijheatmasstransfer.2011.12.015

[11] Makinde, O.D., Khan, W.A. and Khan, Z.H. (2013) Buoyancy Effects on MHD Stagnation Point Flow and Heat Transfer of a Nanofluid Past a Convectively Heated Stretching/Shrinking Sheet. International Journal of Heat and Mass Transfer, 62, 526-533.

https://doi.org/10.1016/j.ijheatmasstransfer.2013.03.049

[12] Rashidi, M.M., Abelman, S. and Mehr, N.F. (2013) Entropy Generation in Steady MHD Flow due to a Rotating Porous Disk in a Nanofluid. International Journal of Heat and Mass Transfer, 62, 515-525.

https://doi.org/10.1016/j.ijheatmasstransfer.2013.03.004

[13] Sheikholeslami, M., Abelman, S. and Ganji, D.D. (2014) Numerical Simulation of MHD Nanofluid Flow and Heat Transfer Considering Viscous Dissipation. International Journal of Heat and Mass Transfer, 79, 212-222.

https://doi.org/10.1016/j.ijheatmasstransfer.2014.08.004

[14] Ganesan, P. and Palani, G. (2004) Finite Difference Analysis of Unsteady Natural Convection MHD Flow past an Inclined Plate with Variable Surface Heat and Mass Flux. International Journal of Heat and Mass Transfer, 47, 4449-4457.

https://doi.org/10.1016/j.ijheatmasstransfer.2004.04.034

[15] Josserand, J., Marty, P. and Alemany, A. (1993) Pressure and Drag Measurements on a Cylinder in a Liquid Metal Flow with an Aligned Magnetic Field. Fluid Dynamics Research, 11, 107-117.

https://doi.org/10.1016/0169-5983(93)90010-8

[16] Mutschke, G., Shatrov, V. and Gerbeth, G. (1998) Cylinder Wake Control by Magnetic Fields in Liquid Metal Flows. Experimental Thermal and Fluid Science, 16, 92-99.

https://doi.org/10.1016/S0894-1777(97)10007-3

[17] Lahjomri, J., Caeran, and Alemani, A. (1993) The Cylinder Wake in a Magnetic Field aligned with the Velocity. Journal of Fluid Mechanics, 253, 421-448.

https://doi.org/10.1017/S0022112093001855

[18] Henoch, C. and Stace, J. (1995) Experimental Investigation of a Salt Water Turbulent Boundary Layer Modified by an Applied Streamwise Magnetohydrodynamic body Force. Physics of Fluids, 7, 1371-1383.

https://doi.org/10.1063/1.868525

[19] Mutschke, G., Gerbeth, G., Shatrov, V. and Tomboulides, A. (1997) Two and Three-Dimensional Instabilities of the Cylinder Wake in an Aligned Magnetic Field. Physics of Fluids, 9, 3114-3116.

https://doi.org/10.1063/1.869429

[20] Mutschke, G., Gerbeth, G., Shatrov, V. and Tomboulides, A. (2001) The Scenario of Three-Dimensional Instabilities of the Cylinder Wake in an External Magnetic Field: A Linear Stability Analysis. Physics of Fluids, 13, 723-734.

https://doi.org/10.1063/1.1344895

[21] Golubev, A.E., Karsilnikov, E.Y. and Lushchik, V.G. (2009) An Experimental Investigation of Liquid Metal Flow Past a Cylinder in Longitudinal Magnetic Field. High Temperature, 47, 432-437.

https://doi.org/10.1134/S0018151X09030171

[22] Sekhar, T.V.S., Sivakumar, R. and Ravikumar, T.V.R. (2006) Flow around a Circular Cylinder in an External Magnetic Field at High Reynolds Number. International Journal of Numerical Methods for Heat and Fluid Flow, 16, 740-759.

https://doi.org/10.1108/09615530610679084

[23] Sekhar, T.V.S., Sivakumar, R., Kumar, H. and Ravikumar, T.V.R. (2007) Effect of Aligned Magnetic Field on the Steady Viscous Flow past a Circular Cylinder. Applied Mathematical Modelling, 31, 130-139.

https://doi.org/10.1016/j.apm.2005.08.011

[24] Sekhar, T.V.S., Sivakumar, R. and Ravikumar, T.V.R. (2009) Effect of Magnetic Reynolds Number on the Two-Dimensional Hydromagnetic Flow around a Cylinder. International Journal for Numerical Methods is Fluids, 59, 1351-1368.

https://doi.org/10.1002/fld.1870

[25] Grigoriadis, D.G.E., Sarris, I.E. and Kassinos, S.C. (2010) MHD Flow past a Circular Cylinder Using the Immersed Boundary Method. Computers and Fluids, 39, 345-358.

https://doi.org/10.1016/j.compfluid.2009.09.012

[26] Tuann, S.Y. and Olson, M.D. (1978) Numerical Studies of the Flow around a Circular Cylinder by Finite Element Method. Computers and Fluids, 6, 219-240.

https://doi.org/10.1016/0045-7930(78)90015-4

[27] Bharti, R.P., Chhabra, R.P. and Eswaran, V. (2006) Steady Flow of Power Law Fluids Across a Circular Cylinder. The Canadian Journal of Chemical Engineering, 84, 406-421.

https://doi.org/10.1002/cjce.5450840402

[28] Bharti, R.P., Chhabra, R.P. and Eswaran, V. (2007) Steady Forced Convection Heat Transfer from a Heated Circular Cylinder to Power-Law Fluids. International Journal of Heat and Mass Transfer, 50, 977-990.

https://doi.org/10.1016/j.ijheatmasstransfer.2006.08.008

[29] Bharti, R.P., Chhabra, R.P. and Eswaran, V. (2007) A Numerical Study of the Steady Forced Convection Heat Transfer from an Unconfined Circular Cylinder. Heat Mass Transfer, 43, 639-648.

https://doi.org/10.1007/s00231-006-0155-1

[30] Senturk, K., Tessarotto, M. and Aslan, N. (2009) Numerical Solutions of Liquid Metal Flows by Incompressible Magneto-Hydrodynamics with Heat Transfer. Numerical Methods in Fluids, 60, 1200-1221.

https://doi.org/10.1002/fld.1928

[31] Uda, N., Miyazawa, A., Inoue, S., Yamaoka, N., Horike, H. and Miyazaki, K. (2001) Forced Convection Heat Transfer and Temperature Fluctuations of Lithium under Transverse Magnetic Field. Journal of Nuclear Science and Technology, 38, 936-943.

https://doi.org/10.1080/18811248.2001.9715120

[32] Genin, L., Dorofeev, D.I., Zhilin, V.G., Ivochkin, Y.P., Listratov, Y.I., Razuvanov, N.G. and Sviridov, V.G. (2007) An Experimental Investigation into the Heat Transfer Developed Along a Uniformly Heated Tube Carrying a Flow of Liquid Metal and Placed in Transverse Magnetic Field. Thermal Engineering, 54, 223-230.

https://doi.org/10.1134/S0040601507030093

[33] Sanitjai, S. and Goldstein, R.J. (2004) Forced Convection Heat Transfer from a Circular Cylinder in Crossflow to Air and Liquids. International Journal of Heat and Mass Transfer, 47, 4795-4805.

https://doi.org/10.1016/j.ijheatmasstransfer.2004.05.012

[34] Andreev, O., Kolesnikov, Y. and Thess, A. (2006) Experimental Study of Liquid Metal Channel Flow Under the Influence of a Nonuniform Magnetic Field. Physics of Fluids, 18, 65108.

https://doi.org/10.1063/1.2213639

[35] Kim, S.J. and Lee, C.M. (2000) Investigation of the Flow around a Circular Cylinder under the Influence of an Electromagnetic Force. Experiments in Fluids, 28, 252-260.

https://doi.org/10.1007/s003480050385

[36] Kim, S.J. and Lee, C.M. (2001) Control of Flows around a Circular Cylinder: Suppression of Oscillatory Lift Force. Fluid Dynamics Research, 29, 47-63.

https://doi.org/10.1016/S0169-5983(01)00019-3

[37] Yoon, H.S., Chun, H.H., Ha, M.Y. and Lee, H.G. (2004) A Numerical Study on the Fluid Flow and Heat Transfer Around a Circular Cylinder in an Aligned Magnetic Field. International Journal of Heat and Mass Transfer, 47, 4075-4087.

https://doi.org/10.1016/j.ijheatmasstransfer.2004.05.015

[38] Dousset, V. and Potheart, A. (2008) Numerical Simulations of a Cylinder Wake under a Strong Axial Magnetic Field. Physics of Fluids, 20, 17104.

https://doi.org/10.1063/1.2831153

[39] Sanyasiraju, Y.V.S.S. and Manjula, V. (2005) Flow past an Impulsively Started Circular Cylinder Using a Higher-Order Semicompact Scheme. Physical Review E, 72, 016709.

https://doi.org/10.1103/PhysRevE.72.016709

[40] Raju, B.H.S. and Sekhar, T.V.S. (2012) Higher Order Compact Scheme Combined with Multigrid Method for Momentum, Pressure and Energy Equations in Cylindrical Geometry. The Open Numerical Methods Journal, 4, 46-58.

https://doi.org/10.2174/1876389801204010046

[41] Ricou, R. and Vives, C. (1982) Local Velocity and Mass Transfer Measurements in Molten Metals Using an Incorporated Magnet Probe. International Journal of Heat and Mass Transfer, 25, 1579-1588.

https://doi.org/10.1016/0017-9310(82)90036-9

[42] Takami, H. and Keller, H.B. (1969) Steady Two-Dimensional Viscous Flow of an Incompressible Fluid Past a Circular Cylinder. Physics of Fluids, 12, 51-56.

https://doi.org/10.1063/1.1692469

[43] Dennis, S.C.R. and Chang, G.Z. (1970) Numerical Solutions for Steady Flow Past a Circular Cylinder at Reynolds Number up to 100. Journal of Fluid Mechanics, 42, 471-489.

https://doi.org/10.1017/S0022112070001428

[44] Collins, W.M. and Dennis, S.C.R. (1973) Flow past an Impulsively Started Circular Cylinder. Journal of Fluid Mechanics, 60, 105-127.

https://doi.org/10.1017/S0022112073000078

[45] Nieuwstadt, F. and Keller, H.B. (1973) Viscous Flow Past Circular Cylinders. Computers and Fluids, 1, 59-71.

https://doi.org/10.1016/0045-7930(73)90026-1

[46] Coutanceau, M. and Bouard, R. (1977) Experimental Determination of the Main Features of the Viscous Flow in the Wake of a Circular Cylinder in Uniform Translation. Journal of Fluid Mechanics, 79, 231-256.

https://doi.org/10.1017/S0022112077000135

[47] Fornberg, B. (1980) A Numerical Study of Steady Viscous Flow Past a Circular Cylinder. Journal of Fluid Mechanics, 98, 819-855.

https://doi.org/10.1017/S0022112080000419

[48] Rogers, S.E. and Kwak, D. (1988) An Upwind-Differencing Scheme for the Incompressible Navier-Stokes Equations. Applied Numerical Mathematics, 8, 43-64.

https://doi.org/10.1016/0168-9274(91)90097-J

[49] He, X. and Doolen, G. (1997) Lattice Boltzmann Method on Curvilinear Coordinate System. Journal of Computational Physics, 34, 306-315.

https://doi.org/10.1006/jcph.1997.5709

[50] Park, J., Kwon, K. and Choi, H. (1998) Numerical Solutions of Flow Past a Circular Cylinder at Reynolds Numbers up to 160. KSME International Journal, 12, 1200-1205.

https://doi.org/10.1007/BF02942594

[51] Ye, T., Mittal, R., Udhayakumar, H.S. and Shyy, W. (1999) An Accurate Cartesian Grid Method for Viscous Incompressible Flows with Complex Immersed Boundaries. Journal of Computational Physics, 156, 209-240.

https://doi.org/10.1006/jcph.1999.6356

[52] Niu, X.D., Chew, Y.T. and Shu, C. (2003) Simulation of Flows Around an Impulsively Started Circular Cylinder by Taylor Series Expansion- and Least-Square Based Lattice Boltzmann Method. Journal of Computational Physics, 188, 176-193.

https://doi.org/10.1016/S0021-9991(03)00161-X

[53] Silva, L.E.A.L.F., Silveira-Neto, A. and Damasceno, J.J.R. (2003) Numerical Simulation of Two-Dimensional Flows over a Circular Cylinder Using the Immersed Boundary Method. Journal of Computational Physics, 189, 351-370.

https://doi.org/10.1016/S0021-9991(03)00214-6

[54] Tseng, Y.H. and Ferziger, J.H. (2003) A Ghost-Cell Immersed Boundary Method for Flow in Complex Geometry. Journal of Computational Physics, 192, 593-623.

https://doi.org/10.1016/j.jcp.2003.07.024

[55] Chakraborty, J., Verma, N. and Chhabra, R.P. (2004) Wall Effects in the Flow Past a Circular Cylinder in a Plane Channel: A Numerical Study. Chemical Engineering and Processing: Process Intensification, 43, 1529-1537.

https://doi.org/10.1016/j.cep.2004.02.004

[56] Ding, H., Shu, C., Yeo, K.S. and Xu, D. (2004) Simulation of Incompressible Viscous Flow Past a Circular Cylinder by Hybrid FD Scheme and Meshless Least Square Based Finite Difference Method. Computer Methods in Applied Mechanics and Engineering, 193, 727-744.

https://doi.org/10.1016/j.cma.2003.11.002

[57] Shu, C., Liu, N. and Chew, Y.T. (2007) A Novel Immersed Boundary Velocity Correction-Lattice Boltzmann Method and Its Application to Simulate Flow Past a Circular Cylinder. Journal of Computational Physics, 226, 1607-1622.

https://doi.org/10.1016/j.jcp.2007.06.002

[58] Udhayakumar, S., Sekhar, T.V.S. and Sivakumar, R. (2015) Numerical Experiments on the Study of Mixed Convection Flow in Cylindrical Geometry. Numerical Heat Transfer, Part A: Applications, 68, 870-886.

https://doi.org/10.1080/10407782.2015.1023129

[59] Udhayakumar, S., Abin Rejeesh, A.D., Sekhar, T.V.S. and Sivakumar, R. (2016) Numerical Investigation of Magnetohydrodynamic Mixed Convection Over an Isothermal Circular Cylinder in Presence of an Aligned Magnetic Field. International Journal of Heat and Mass Transfer, 95, 379-392.

https://doi.org/10.1016/j.ijheatmasstransfer.2015.11.041

[60] Udhayakumar, S., Abin Rejeesh, A.D., Sekhar, T.V.S. and Sivakumar, R. (2016) Study of Directional Control of Heat Transfer and Flow Control in the Magnetohydrodynamic Flow in Cylindrical Geometry. International Journal of Heat and Fluid Flow, 61, 482-498.

https://doi.org/10.1016/j.ijheatfluidflow.2016.06.011

[61] Dennis, S.C.R., Hudson, J.D. and Smith, N. (1968) Steady Laminar Forced Convection from a Circular Cylinder at Low Reynolds Numbers. Physics of Fluids, 11, 933-940.

https://doi.org/10.1063/1.1692061

[62] Apelt, C.J. and Ledwich, M.A. (1979) Heat Transfer in Transient and Unsteady Flows Past a Heated Circular Cylinder in the Range . Journal of Fluid Mechanics, 95, 761-777.

https://doi.org/10.1017/S0022112079001683

[63] Badr, H.M. (1983) A Theoretical Study of Laminar Mixed Convection form a Horizontal Cylinder in a Cross Stream. International Journal of Heat and Mass Transfer, 26, 639-653.

https://doi.org/10.1016/0017-9310(83)90014-5

[64] Jafroudi, H. and Yang, H.T. (1986) Steady Laminar Forced Convection from a Circular Cylinder. Journal of Computational Physics, 65, 46-56.

https://doi.org/10.1016/0021-9991(86)90003-3

[65] Chen, C.K., Yang, Y.T. and Wu, S.R. (1994) Laminar Mixed Convection from a Circular Cylinder Using a Body-fitted Coordinate System. Journal of Thermophysics and Heat Transfer, 8, 695-701.

https://doi.org/10.2514/3.600

[66] Aldoss, T.K., Ali, Y.D. and Al-Nimr, M.A. (1996) MHD Mixed Convection from a Horizontal Circular Cylinder. Numerical Heat Transfer, Part A: Applications, 30, 379-396.

https://doi.org/10.1080/10407789608913846

[67] Lange, C.F., Durst, F. and Breuer, M. (1998) Momentum and Heat Transfer from Cylinders in Laminar Crossflow at . International Journal of Heat and Mass Transfer, 41, 3409-3430.

https://doi.org/10.1016/S0017-9310(98)00077-5

[68] Sparrow, E.M., Abraham, J.P. and Tong, J.C.K. (2004) Archival Correlations for Average Heat Transfer Coefficients for Non-Circular and Circular Cylinders and for Spheres in Crossflow. International Journal of Heat and Mass Transfer, 47, 5285-5296.

https://doi.org/10.1016/j.ijheatmasstransfer.2004.06.024

[69] Soares, A.A., Ferreira, J.M. and Chhabra, R.P. (2005) Flow and Forced Convection Heat Transfer in Crossflow of Non-Newtonian Fluids over a Circular Cylinder. Industrial and Engineering Chemistry Research, 44, 5815-5827.

https://doi.org/10.1021/ie0500669

[70] Biswas, G. and Sarkar, S. (2009) Effect of Thermal Buoyancy on Vortex Shedding Past a Circular Cylinder in Cross-Flow at Low Reynolds Number. International Journal of Heat and Mass Transfer, 52, 1897-1912.

https://doi.org/10.1016/j.ijheatmasstransfer.2008.08.034

[71] Soares, A.A., Couto, N.D., Naia, D., Goncalves, N.J. and Rouboa, A. (2012) Numerical Investigation of Effects of Buoyancy Around a Heated Circular Cylinder in Parallel and Contra Flow. Journal of Mechanical Science and Technology, 26, 1501-1513.

https://doi.org/10.1007/s12206-012-0310-1

[72] Vimala, S., Damodaran, S., Sivakumar, R. and Sekhar, T.V.S. (2016) The Role of Magnetic Reynolds Number in MHD Forced Convection Heat Transfer. Applied Mathematical Modelling, 40, 6737-6753.

https://doi.org/10.1016/j.apm.2016.02.019

[73] Leibovich, S. (1967) Magnetohydrodynamic Flow at a Rear Stagnation Point. Journal of Fluid Mechanics, 29, 401-413.

https://doi.org/10.1017/S0022112067000916

[74] Goldsworthy, F.A. (1961) Magnetohydrodynamic Flows of a Perfectly Conducting Viscous Fluid. Journal of Fluid Mechanics, 11, 519-528.

https://doi.org/10.1017/S0022112061000706