This paper presents three-dimensional direct numerical simulations of liquid metal flow around a circular cylinder placed symmetrically in a rectangular duct, under a wide range of magnetic field intensities. Results are presented for values of the Hartmann number (based on the duct width) in the range of 0 ⩽ *Ha* ⩽ 1120, and the Reynolds number (based on the cylinder diameter and centerline velocity) in the range 0 ⩽ *Re*_{c} ⩽ 5000. The generated flow regimes and the associated critical values of parameters are investigated in detail through full three-dimensional simulations. The effect of the magnetic field on the wake structure is discussed in relation to the possible mechanisms for the generation or suppression of vortices, and to previous attempts to model magnetohydrodynamic flows using simplified two-dimensional models. Present results reveal a non-monotonic dependance of the critical Reynolds number for the onset of vortex shedding, with respect to the Hartmann number. For certain combinations of *Ha* and *Re* values, this work confirms the onset of a new flow regime, the existence of which has been recently suggested based on quasi-two-dimensional simulations. Unexpectedly, the spanwise distribution of the force coefficients along the cylinder is found to become more three-dimensional with increasing *Ha*. Furthermore, the three-dimensional nature of the present simulations reveals additional counter-intuitive features of the new regime that could not possibly had been captured by quasi-two-dimensional models. One such feature, shown here for the first time, is an increase in the flow unsteadiness with increasing intensity of the magnetic field.

## I. INTRODUCTION

The study of the flow of a conducting fluid over a confined cylinder placed in a rectangular duct, under the influence of an external homogeneous magnetic field,**B**, is of considerable practical interest. Liquid metal flows in confined arrangements under the presence of a magnetic field play a significant role in a wide range of applications in engineering and industrial processes, such as stirring of melts in the metallurgical industry and cooling of liquid metal blankets in fusion reactors. Introducing a cylindrical obstacle in the flow can induce vortices and enhance heat transfer rates in favor of these processes.^{1–4} Therefore, understanding the enhancement or suppression of flow instabilities, and the induction of secondary flows can provide valuable knowledge of practical importance.

Under these conditions, the wake flow characteristics are dominated by magnetohydrodynamic (MHD) phenomena. Thus, the flow character depends on two dimensionless parameters. One of these is the Reynolds number,*Re* = *UL*/ν, a purely hydrodynamic parameter giving a measure of the ratio of inertial to viscous forces. Here, *U* and *L* stand for a typical velocity and length scale of the flow, respectively, and ν for the kinematic viscosity of the fluid. The second is the interaction parameter, *N* = σ*LB*^{2}/ρ*U*, which represents the ratio of electromagnetic to inertial forces. Here σ and ρ stand for the electric conductivity and density of the fluid, respectively. The interaction parameter can be expressed in terms of the Reynolds number as *N* = *Ha*^{2}/*Re*, where

^{5,6}

The Lorentz force acts in two ways; first it modifies the velocity distribution of the primary flow and, second due to its dissipative nature, it dampens the fluctuations of the unsteady flow.^{5,6} In the case of a rectangular duct with insulating walls considered here, the electric currents induced in the bulk flow are closed within the fluid. At the walls perpendicular to the magnetic field dissipation takes place while at the center of the duct the Lorentz forces tend to adjust the velocity in the bulk to that close to the walls. As a result, the flow exhibits a large, uniform core, where the velocity becomes nearly constant, and relatively thin boundary layers, where the velocity falls to zero. The boundary layers that develop at the walls parallel to the magnetic field, are called Shercliff layers. At the walls perpendicular to the magnetic field, thinner boundary layers develop with sharp velocity gradients, called Hartmann layers.^{7} Depending on the intensity of the magnetic field, these layers can be extremely thin.

In such wall-bounded flows, if both the Hartmann number and the interaction parameter are very high, i.e., *Ha* ≫ 1 and *N* ≫ 1, any velocity variations parallel to the direction of the magnetic field are strongly suppressed by Joule dissipation, except in the viscous Hartmann layers.^{6,8} In these regions, because of the no-slip boundary condition, strong velocity gradients persist. Moreover, Lorentz forces diffuse momentum along the field lines and vortices whose axes are inclined to the magnetic field are strongly damped. On the other hand, vortices whose axes are parallel to the magnetic field are only influenced by friction in the Hartmann layer, which, however is significantly (*Ha* times) less efficient than Joule damping.^{6,8} These electromagnetic effects drive the flow toward a quasi-two-dimensional state, in which one can distinguish a core flow region, where both the gradient of the velocity and the velocity component in the direction of the magnetic field are negligible, and two thin Hartmann layers, where Lorentz forces are balanced by the viscous friction. Under these conditions, integrating the three-dimensional flow equations along the direction of the magnetic field, and adding a linear breaking term containing the net effect of the Lorentz forces, yields a model equation satisfying the two-dimensional hydrodynamic Navier-Stokes equation.^{9,10} An obvious advantage of this quasi-two-dimensional (Q2D) model is that it can reduce the computational cost of an initially three-dimensional problem by translating it to a two-dimensional one. Traditionally, most of the simulations of MHD flows in such geometrically complex domains, as the one considered in this study, rely on this model due to the high computational costs associated with full three-dimensional simulations.

While several studies in the literature deal with MHD flow passed a cylinder, only a few consider the case where the magnetic field is aligned with the axis of a confined cylinder. Thus, the only results previously reported in the literature for this case are from Refs. 3 and 4, and 11–13. Results for *N* ≫ 1, showed that the stability of the wake flow is significantly influenced by the magnetic field. Refs. 3,4,11,12 showed that for moderate and high values of the Hartmann number, *Ha* ⩾ 320, the onset of vortex shedding is delayed till higher Reynolds numbers. For a given blockage ratio β (cylinder diameter over duct height), a linear dependence of the corresponding critical Reynolds number on *Ha* has been observed.^{3,4,11} This behavior was attributed to the damping of the fluctuations of the unsteady flow by the magnetic field. Hussam *et al.*,^{4} carrying out a parametric study based on the duct blockage ratio, observed furthermore that with increasing β, the slope of the critical Reynolds number curve is also found to increase. This was associated with the interaction of the cylinder wake and the top and bottom walls.^{4,12} By increasing further the Reynolds number, experiments from Frank *et al.*^{11} and numerical simulations from Dousset and Pothérat^{3} have revealed the presence of an additional flow regime that follows the laminar periodic vortex shedding and does not have a counterpart in the purely hydrodynamic case. In this regime, the flow is characterized by irregular vortex shedding patterns, with secondary vortices being released from the top and bottom walls.^{3} Again, as the imposed magnetic field is increased, transition to this regime is delayed. Mück *et al.*^{13} showed that for interaction parameters in the range 2 ⩽ *N* ⩽ 10 transition from a time-dependent three-dimensional flow to a two-dimensional state occurred. After the transition to a quasi-two-dimensional flow, cigar-shape vortices (larger diameter at the center, smaller diameter near the Hartmann walls) were observed. They attributed the curvature of the vortices obtained near Hartmann walls not to the higher viscous damping in these regions, but rather to an inertial contribution to the electric potential in the core flow.

However, most of the aforementioned studies considering the case of a cylinder confined in a duct were not based on a fully three-dimensional analysis. Refs. 3,4,12 performed numerical simulations using a Q2D model. Even in the experimental study carried out by Frank *et al.*,^{11} a non intrusive measurement device was used to calculate the core flow quantities by measuring the electric potential only at one Hartmann wall. Only Mück *et al.*^{13} performed three-dimensional simulations. However, in that case an analytical wall boundary model was used, a fact that precluded the resolution of the thin Hartmann layers, therefore making impossible the detailed analysis of the interaction between the vortices and the walls. In addition, due to the limited available computer resources at that time, only two Reynolds numbers could be investigated (*Re* = 200 and 250), and the Hartmann number was varied in the range of 63 < *Ha* < 850.

The present work was motivated by the lack of detailed three-dimensional simulations addressing the case of confined MHD flows with flow-obstructions. In this study, we perform fully three-dimensional direct numerical simulations around a circular cylinder in a rectangular duct. The Hartmann and Schercliff layers developing along solid surfaces are fully resolved without using any wall model. Results are presented for values of the Hartmann number (based on the duct width) in the range of 0 ⩽ *Ha* ⩽ 1120, and the Reynolds number (based on the cylinder diameter and centerline velocity) in the range 0 ⩽ *Re*_{c} ⩽ 5000. Most of the studies in the literature have been based on the Q2D model, whose validity is restricted to high *Ha* numbers. Hence a main objective of this study is to fill in the gap in the area of low *Ha* for this type of flows. Fully three-dimensional direct numerical simulations are used to assess independently the performance and range of validity of simplified models, such as the Q2D model. Leveraging the information provided by the fully 3D simulations, we attempt to address physical effects, which, as already pointed out by Dousset and Pothérat,^{3} are not properly accounted for in the Q2D model due its inherent limitations. For example, we aim to shed light on the effect of the magnetic field intensity on the generated flow regimes. In addition, we describe the evolution of the flow characteristics (force coefficients and recirculation length) with *Ha* and *Re*, and discuss their spanwise distribution along the cylinder. Providing the full three-dimensional flow characteristics could also be beneficial to future modellers working towards the construction of advanced near-wall models for these type of flows.

The paper is organized as follows. In Sec. II, we define the test case considered and recall the governing equations, along with the proper boundary conditions. Section III is devoted to the numerical aspects of the paper; in particular, we describe the numerical algorithm and the development of the grid design, including its validation and verification. In Sec. IV, we present the evolution of the critical Reynolds number for the onset of vortex shedding with respect to *Ha*, while in Secs. VI and V results under steady and unsteady flow conditions are presented, respectively. Finally, Sec. VII presents the concluding remarks derived from the numerical experiments.

## II. PROBLEM STATEMENT AND FORMULATION

### A. Flow configuration and mathematical formulation

The geometry considered in this study is shown in Figure 1. The geometry consists of a circular cylinder of diameter *D* and length *W* placed symmetrically in a rectangular duct. The ratio of the cylinder length to its diameter defines the aspect ratio,

while the ratio of the cylinder diameter to the duct height *H* (distance between the duct walls perpendicular to the cylinder axis) defines the blockage ratio,

Results are reported for β = 1/4 and α = 4, allowing for a direct comparison with the numerical simulations of Dousset and Pothérat.^{3} The origin of the coordinates axes is placed at the cylinder's geometrical center, with the *z* direction being aligned with the cylinder axis. The duct inlet is placed at a distance of *L*_{i} = 12.5*D* upstream of the circular cylinder, while the outlet is located at *L*_{o} = 35.5*D* behind the body. This choice of parameters ensures minimal distortion of the flow structures due to the outlet boundary conditions, while maintaining a reasonable computational cost, as will be described further in Sec. III C.

The flow of an incompressible, electrically conducting Newtonian fluid, of density ρ, kinematic viscosity ν and electrical conductivity σ, is considered. An external homogeneous magnetic field of amplitude *B*_{0} is applied along the spanwise direction *z*. The walls of the cylinder and the duct are assumed to be electrically insulating. In most cases of liquid metal flow encountered in industrial applications, the magnetic Reynolds number is in general very small, and thus the induced magnetic field is found to be negligible when compared to the imposed magnetic field. Hence, in this study it is also assumed that the magnetic Reynolds number is much smaller than unity, *R*_{m} = μσ*UL* ≪ 1, where μ stands for the fluid magnetic permeability, therefore the *quasi-static* approximation is invoked.^{14} Under these assumptions, the non-dimensional magnetohydrodynamic equations governing the flow can be written as

The dimensionless flow variables appearing in Eqs. (3)–(6) are obtained from their dimensional counterpart by using the centerline inflow velocity, *U*_{c}, as the characteristic velocity scale and the cylinder's diameter, *D*, as the characteristic length scale. Therefore, the following variables, *t*, **u**, *p*, **j**, **B**, ϕ, which denote the dimensionless time, velocity, pressure, current density, imposed magnetic field, and electric potential are scaled with *D*/*U*_{c}, *U*_{c}, ρ*U*_{c}^{2}, σ*U*_{c}*B*_{0}, *B*_{0}, and *DU*_{c}*B*_{0}, respectively. Two dimensionless parameters appear in Eq. (4). The Reynolds number,*Re*_{c}, and the interaction parameter, *N*,

where *Ha* is the Hartmann number, defined as

The drag and lift coefficients are both determined by considering the viscous and the pressure forces on the cylinder surface,

Here *F*_{D} and *F*_{L} are the drag and lift forces per unit length of the cylinder, defined as the streamwise and crossflow components of

where *n*_{j} denotes the component of the unit normal vector of an element on the cylinder surface pointing in the direction of *x*_{j}. The Strouhal number, characterizing the vortex shedding phenomenon, is based on the dominant frequency of the lift coefficient

The base pressure coefficient is defined as

where *p*_{b} is the spatial average, with respect to the spanwise coordinate *z*, of pressure at the rear stagnation point of the cylinder (180° from the front), and *p*_{∞} is the free-stream pressure at the inlet boundary (*x* = −12.5*D*, *y* = 0). Finally, the recirculation length, *L*_{r}, is defined as the average distance, along the center-plane (*y* = 0), from the rear stagnation point of the cylinder to the point where the streamwise velocity is zero.

### B. Boundary conditions

Boundary conditions are given at the domain's inlet by prescribing the analytical solution for an MHD flow in a rectangular duct with insulating walls, in the absence of a cylinder, as described by Müller and Bühler.^{5} A no-slip boundary condition is imposed on the cylinder surface and the duct walls,

At the outlet, a convective boundary condition is applied, which serves to minimize reflective effects and avoid the distortion of the flow structures leaving the domain, while ensuring mass conservation,

Here *U*_{conv} ≡ ∫*u*_{j}*n*_{j}*ds*/∫*ds*, with the integrals being calculated over the outlet surface.

For the electric potential, the condition ∂ϕ/∂*n* = 0 is imposed at the duct walls to ensure that no current crosses the insulating walls. At the inlet and outlet of the duct, the following boundary condition is used:

ensuring that no electric currents enter or leave the domain.^{13} This boundary condition reduces to the requirement of zero normal component of the current at the inlet and outlet.

## III. NUMERICAL ASPECTS

### A. Numerical method

The current computations have been performed using the unstructured collocated nodal-based finite-volume code, named CDP, developed at the Center of Turbulence Research (Stanford/NASA Ames). A second-order accurate centered-difference scheme, with skewness corrections, is used for the discretization of the diffusive and nonlinear terms. The mass and momentum equations (3) and (4) are coupled using a fractional-step method. The time integration of the flow equations is done using a Crank-Nicholson scheme for the diffusive terms. The nonlinear terms are treated semi-implicitly. Details of the code have been described extensively in Refs. 15,16.

For this study, we have complemented this code with a module to compute the Lorentz force and include it in the momentum balance. An important feature of the MHD module is the use of a consistent and conservative scheme for the computation of the Lorentz force.^{17} A very detailed description of this scheme and other numerical techniques used in the MHD module is reported in Albets-Chico *et al.*^{18} The implementation and integration of this module in CDP has been validated for several MHD flows and in various geometries, including MHD channel flow^{19} and MHD pipe flow.^{20} In Subsection III B, we consider the case of fully-developed MHD flow in a rectangular duct, in the absence of a cylinder, in order to further validate the code and ensure a proper grid design and density.

### B. 2D fully-developed MHD duct flow

In order to provide the appropriate grid for the different flow regimes analyzed in this work, and in particular to establish proper resolution criteria regarding the Hartmann and Shercliff layers, a series of preliminary numerical computations are performed in a rectangular duct without the presence of the cylinder. Additionally, this step further validates the numerical code and its ability to deal with basic MHD flow dynamics, as the numerical results can be assessed against the known semi-analytical solution for the duct flow without a cylinder.^{5}

For this reason, the configuration of Figure 1 is first considered without the presence of a cylinder. The flow is assumed to be fully developed, hence a two-dimensional case is addressed. Periodic boundary conditions are specified for the velocity in the streamwise direction, while a fixed pressure drop is used in the whole domain as a forcing term and computed iteratively so as to enforce a constant mass flow. The Reynolds number is fixed at *Re*_{c} = 100, while Hartmann numbers in the range *Ha* = 0–1120 are considered.

Four different grid configurations are tested, differing in the spatial resolution at the Hartmann and Shercliff walls, as summarized in Table I. In the direction normal to these walls, linear grid stretching is applied. The non-dimensional velocity profiles obtained using the four different grids are shown in Figure 2. In this figure, the velocity is scaled by *W*^{2}ν^{−1}ρ^{−1}( − ∂*p*/∂*x*), where ∂*p*/∂*x* is the pressure gradient driving the flow. This scaling relates the velocity profile with the magnitude of the pressure drop required to produce it. This is a rather strict but necessary way of presenting the results, because otherwise one could obtain very similar-looking velocity profiles, even with a very coarse mesh, just by fixing the mass flow rate. However, in such a case, the viscous stresses at the walls, and hence the pressure drop, would be over-predicted. Therefore, it is necessary not only to obtain an accurate velocity profile, but also to predict accurately the pressure drop needed to drive the fluid through the duct.

Grid . | N_{y} × N_{z}
. | Δy_{min}
. | Δz_{min}
. |
---|---|---|---|

G1 | 45 × 45 | 0.02 | 0.02 |

G2 | 65 × 65 | 0.01 | 0.01 |

G3 | 85 × 85 | 0.005 | 0.005 |

G4 | 85 × 104 | 0.005 | 0.002 |

Grid . | N_{y} × N_{z}
. | Δy_{min}
. | Δz_{min}
. |
---|---|---|---|

G1 | 45 × 45 | 0.02 | 0.02 |

G2 | 65 × 65 | 0.01 | 0.01 |

G3 | 85 × 85 | 0.005 | 0.005 |

G4 | 85 × 104 | 0.005 | 0.002 |

Results are compared with the analytical solution in terms of a series expansion provided by Müller and Bühler.^{5} The solution was obtained for the case of an MHD flow in a rectangular duct with insulating walls and a strong transverse magnetic field. For brevity, only results for *Ha* = 320 and 1120 are shown. For lower Hartmann numbers even better agreement with the semi-analytical solution is obtained, as the electrodynamic layers are less challenging from a numerical point of view. For *Ha* = 320, the maximum local discrepancy between the velocity profiles predicted numerically on grid G3 with respect to the ones obtained analytically is less than 2%. This percentage difference is further reduced to less than 0.4% when the analytical solution is compared to the results obtained using grid G4. Hence, one can conclude that the intermediate grid G3 is sufficiently fine to resolve the flow for *Ha* ⩽ 320, as lower Hartmann numbers present thicker electrodynamic layers. On the other hand, for higher values of *Ha* (*Ha* = 1120), the maximum local discrepancy between the analytical velocity profiles with respect to the ones obtained numerically on grid G3 is around 17%, which is unacceptably high. This fact makes imperative the use of the finer grid G4 for this range of Hartmann numbers, so as to reduce the numerical error to less than 3%. In general, when using grid G3 for *Ha* ⩽ 320 and grid G4 for 320 < *Ha* ⩽ 1120 the agreement is very good for all the examined cases, confirming the ability of our numerical code, together with the appropriate grids, to accurately resolve the Hartmann and Shercliff layers.

### C. 3D mesh details and computational domain

In this subsection, we present the mesh design for the configuration of interest (see Figure 1) that is, in the presence of a circular cylinder, together with a domain dependence analysis.

The influence of the grid resolution on the computed physical characteristics of the purely hydrodynamic flow was examined in detail in a previous paper of ours,^{21} where we considered the case of a circular cylinder confined in a channel with a slightly smaller blockage ratio. Grid development for the current configuration was guided by that study, with a few modifications to compensate for the different blockage ratio and the additional confinement imposed by the walls at the cylinders ends. Furthermore, results related to the grid resolution of the Shercliff and Hartmann layers found in Subsection III B were considered when generating the grid. In view of the exponential decay of the thickness of the Hartmann layers with increasing *Ha* (δ_{H} ∼ 1/*Ha*), and in order to accommodate the wide range of *Ha* numbers investigated herein, two grids, differing only in the resolution of the Hartmann layers have been used. Hence, based on the results obtained for the fully-developed case in Subsection III B, two grids were generated, namely GL-M and GH-M, with a minimum grid spacing, along the direction of the magnetic field, of Δ*z*_{min} = 0.005*D* and 0.002*D*, respectively.

The grid naming convention adopted here is as follows. The first segment of the grid name (two letters preceding the hyphen) indicate the relative grid resolution, i.e., GL- and GH- stand for lower and higher grid resolution, respectively. The second segment of the grid name (letter following the hyphen) indicates the domain size, i.e., S, M, and L stand for short, medium and long domains, respectively. The short (-S) and long (-L) domains were used only during the domain size independence study, otherwise the medium domain (-M) was used in all other runs. The parameters of the different grid configurations tested are summarized in Table II.

Grid . | Cases . | L_{i}
. | L_{o}
. | Δn_{H}
. | Δn_{S, C}
. | N_{z}
. | Total nodes . | |
---|---|---|---|---|---|---|---|---|

GL-S | Ha = 320 | Re = 2000 | 8.5 | 22.5 | 0.005 | 0.005 | 85 | 4 535 090 |

GL-M | Ha ⩽ 320 | Re ⩽ 5000 | 12.5 | 35.5 | 0.005 | 0.005 | 85 | 5 568 265 |

GL-L | Ha = 320 | Re = 2000 | 16.5 | 45.5 | 0.005 | 0.005 | 85 | 6 419 115 |

GH-M | Ha > 320 | Re ⩽ 5000 | 12.5 | 35.5 | 0.002 | 0.005 | 104 | 6 812 936 |

GF-M | Ha = 320 | Re = 5000 | 12.5 | 35.5 | 0.005 | 0.005 | 85 | 8 600 000 |

GX-M | Ha = 1120 | Re = 5000 | 12.5 | 35.5 | 0.002 | 0.005 | 104 | 10 507 120 |

GR-T | Ha = 320 | Re = 5000 | 7 | 14 | 0.0015 | 0.0015 | 129 | 19 762 978 |

Grid . | Cases . | L_{i}
. | L_{o}
. | Δn_{H}
. | Δn_{S, C}
. | N_{z}
. | Total nodes . | |
---|---|---|---|---|---|---|---|---|

GL-S | Ha = 320 | Re = 2000 | 8.5 | 22.5 | 0.005 | 0.005 | 85 | 4 535 090 |

GL-M | Ha ⩽ 320 | Re ⩽ 5000 | 12.5 | 35.5 | 0.005 | 0.005 | 85 | 5 568 265 |

GL-L | Ha = 320 | Re = 2000 | 16.5 | 45.5 | 0.005 | 0.005 | 85 | 6 419 115 |

GH-M | Ha > 320 | Re ⩽ 5000 | 12.5 | 35.5 | 0.002 | 0.005 | 104 | 6 812 936 |

GF-M | Ha = 320 | Re = 5000 | 12.5 | 35.5 | 0.005 | 0.005 | 85 | 8 600 000 |

GX-M | Ha = 1120 | Re = 5000 | 12.5 | 35.5 | 0.002 | 0.005 | 104 | 10 507 120 |

GR-T | Ha = 320 | Re = 5000 | 7 | 14 | 0.0015 | 0.0015 | 129 | 19 762 978 |

Grid GL-M has been used for the cases where *Ha* ⩽ 320, and grid GH-M for the last two cases with higher Hartmann numbers, *Ha* = 640 and *Ha* = 1120. *A posteriori*, we determined that at least 6 mesh nodes were present inside the Hartmann layers and 18 nodes in the Shercliff layers for all the analyzed cases. For both grids, a hyperbolic tangent-based function was used to stretch the cell sizes in a clustered region close to the cylinder, and a linear grid stretching was applied in the direction normal to the duct walls, as shown in Figure 3. Further upstream and downstream of the cylinder, a uniform grid was used in the streamwise direction.

In order to examine the influence of the computational domain on the predicted flow fields, a domain dependence study was also performed. For this reason, two more grids were tested, namely GL-S and GL-L, with a shorter and longer computational domain, respectively. These grids are based on the spacing of grid GL-M. Simulations were performed at *Ha* = 320, and *Re*_{c} = 2000, over periods of at least 180 dimensionless time units, corresponding to about 45 vortex shedding cycles, with time averaging of results performed over the last 25 shedding cycles, once the flow had reached a statistically steady state.

The discrepancy between the time-averaged values of drag coefficient

*St*resulting from the use of different domain sizes is shown in Table III. Because,

*C*

_{pb}[see Eq. (12)] depends on the distance between the points used to compute the pressure drop, we define a different base pressure coefficient,

*St*almost no variations can be seen. Furthermore, Figure 4 displays the mean streamwise velocities at different one-dimensional slices, for domains GL-S, GL-M, and GL-L. From this comparison it can be concluded that the effect of the convective outlet is very small, and is isolated within a distance of 2

*D*from the outlet. Based on these results, one can conclude that the intermediate domain GL-M is sufficiently long to resolve the flow.

Grid . | Ha
. | Re_{c}
. | $\bar{L}_r$ $L\xafr$
. | $\bar{C}_{pb}^{\prime }$ $C\xafpb\u2032$
. | $\bar{C}_D$ $C\xafD$
. | St
. |
---|---|---|---|---|---|---|

GL-S | 320 | 2000 | 0.490 | 3.149 | 1.715 | 0.260 |

GL-M | 320 | 2000 | 0.488 | 3.152 | 1.714 | 0.260 |

GL-L | 320 | 2000 | 0.479 | 3.159 | 1.714 | 0.260 |

GL-M | 320 | 5000 | 0.524 | 3.028 | 1.627 | 0.241 |

GF-M | 320 | 5000 | 0.514 | 3.054 | 1.644 | 0.248 |

GR-T | 320 | 5000 | 0.513 | 3.063 | 1.661 | 0.249 |

GH-M | 1120 | 5000 | … | 3.727 | 2.148 | 0.265 |

GX-M | 1120 | 5000 | … | 3.728 | 2.152 | 0.265 |

Grid . | Ha
. | Re_{c}
. | $\bar{L}_r$ $L\xafr$
. | $\bar{C}_{pb}^{\prime }$ $C\xafpb\u2032$
. | $\bar{C}_D$ $C\xafD$
. | St
. |
---|---|---|---|---|---|---|

GL-S | 320 | 2000 | 0.490 | 3.149 | 1.715 | 0.260 |

GL-M | 320 | 2000 | 0.488 | 3.152 | 1.714 | 0.260 |

GL-L | 320 | 2000 | 0.479 | 3.159 | 1.714 | 0.260 |

GL-M | 320 | 5000 | 0.524 | 3.028 | 1.627 | 0.241 |

GF-M | 320 | 5000 | 0.514 | 3.054 | 1.644 | 0.248 |

GR-T | 320 | 5000 | 0.513 | 3.063 | 1.661 | 0.249 |

GH-M | 1120 | 5000 | … | 3.727 | 2.148 | 0.265 |

GX-M | 1120 | 5000 | … | 3.728 | 2.152 | 0.265 |

As explained below, it was necessary to design additional grids that would enable us to investigate the presence of the fourth flow regime, which has been reported to come in succession after regular vortex shedding.^{3,11} For this purpose, simulations have been performed for an even higher Reynolds number,*Re*_{c} = 5000, and for *Ha* = 320 and 1120. The high value of the Reynolds number in this regime was found to promote the presence of finer structures, especially at *Ha* = 320 (see Figure 18). To capture these fine structures, grids GL-M and GH-M that were used in cases where *Ha* = 320 and *Ha* = 1120, respectively, were refined even further. Grid refinement was performed only in the streamwise direction, in the vicinity of the cylinder and further downstream, decreasing the maximum value of Δ*x*_{max} from 0.2*D* to 0.1*D*. As summarized in Table II, this yielded two finer grids, namely GF-M (Fine resolution, Medium domain) and GX-M (eXtra fine resolution, Medium domain), with a total of 8.6 × 10^{6} and 10.5 × 10^{6} nodes, respectively. Given the irregular character of the flow in this regime, special care was needed to ensure the validity of statistical quantities. For this purpose, simulations were carried out for longer periods of at least 400 dimensionless time units (100 vortex shedding cycles), while the time-averaging was performed over a minimum period of 300 time units (75 shedding cycles). As we shall discuss in Sec. VI, in the case of *Ha* = 320, *Re* = 5000, the presence of small-scale structures in the near-wake is an inherent and important feature of the flow. In an effort to ensure proper resolution of the smallest scales of motion on grid GF-M, an even finer grid was generated in a shorter domain, namely GR-T (Reference resolution, Truncated domain). The latter grid is up to three times finer in the near-wake region than grid GF-M (see Table II). Because of the extremely high computational cost associated with the use of this grid, only one simulation was performed to provide a reference case at *Ha* = 320, and *Re*_{c} = 5000.

Overall, the agreement between the results obtained with the different grids is very satisfactory, with the discrepancies being below 2.2% at *Ha* = 320, and 0.3% at *Ha* = 1120, as shown in Table III. Figure 5 shows the power spectral density of the crossflow velocity component *u*_{y}, obtained at *x* = 1.5*D*, *y* = 0.5*D*, *z* = 0 using grids GF-M and GR-T. The observed spectra, as we shall see in Sec. VI, are consistent with the instantaneous visualizations of the flow, where the presence of small scales is clearly observed. While the two spectra are in good agreement with each other up to the end of the inertial range, the sharp drop of the energy content at the high-frequency end of the spectrum happens at higher frequencies in the GR-T spectra than in the GF-M spectra. To clarify this issue, we have used the estimates of the size of the smallest scales in MHD turbulence suggested by Pothérat and Dymkou.^{22} Using these estimates, we believe that grid GF-M provides marginal resolution of the smallest scales along the direction of the magnetic field in the turbulence-like regions of the flow. While the effect on the global features of the flow discussed herein is expected to be unimportant, this observation should be kept in mind when interpreting visualizations of turbulence-like structures in the flow that are shown later.

## IV. CRITICAL *Re*, ONSET OF VORTEX SHEDDING

We have first investigated the effect of *Ha* on the laminar flow regime and the transition from steady to unsteady flow. Figure 6(a) shows the critical Reynolds number*Re*_{c,cr} (based on the centerline velocity) for the onset of vortex shedding as a function of the Hartmann number. In order to determine *Re*_{c,cr} for a given value of *Ha*, a series of independent simulations were performed for different values of *Re* starting from a zero velocity field in all cases. Each critical Reynolds number is defined in an interval represented by an error bar. The lower end of the bar corresponds to the simulation with the highest *Re* that leads to a steady flow-regime, while the higher end corresponds to the simulation with the lowest *Re* that leads to unsteady flow. Results are compared with those obtained by Dousset and Pothérat^{3} using a Q2D model. However, it should be noted that in their study, *Re*_{cr} was obtained by performing two series of numerical computations, where *Re* was sequentially increased in small steps till unsteadiness was reached, and then gradually decreased again till steadiness was recovered. Figure 6(b) shows an expanded view in the area of low values of Hartmann number, this time showing the critical Reynolds number*Re*_{b,cr} based on the bulk velocity. Keeping *Re*_{b} fixed while varying *Ha* corresponds to the same mass flow rate, which is important when presenting the transition from steady to unsteady flow.

As seen from these figures, for *Ha* > 80 the critical Reynolds number becomes proportional to the Hartmann number with the linear relation *Re*_{cr} ≃ 0.52*Ha*. This relation is in very good agreement with the results of Dousset and Pothérat,^{3} who report a relation of *Re*_{cr} ≃ 0.43*Ha*. A similar linear relationship between *Re* and *Ha* was also observed by Frank *et al.*^{11} (*Re*_{cr} ≃ 0.47*Ha*) and Hussam *et al.*^{4} (*Re*_{cr} ≃ 0.5*Ha*) for configurations with different blockage and aspect ratio, i.e. β = 0.1, α = 5, and β = 0.2, α = 2.5, respectively. The relative differences between the results obtained by Ref. 3 and the current study (see Table IV), may be at least partly attributed to the different methods used for obtaining the critical Reynolds number. Nevertheless, the agreement between the quasi-two-dimensional model and current three-dimensional results gradually improves with increasing *Ha*. This is expected, due to the enlargement of the quasi-inviscid core flow, which improves the validity of the Q2D model.

Surprisingly, as one decreases the Hartmann number below *Ha* = 80, the critical Reynolds number decreases slightly until it reaches a minimum value,

*Ha*

_{cr}≈ 35. From that point on, it starts to increase up to

*Re*

_{b,cr}≈ 150 as

*Ha*is further decreased towards the hydrodynamic case,

*Ha*= 0. The value of

*Re*

_{cr}in the purely hydrodynamic limit is in reasonable agreement with the experimental investigation of Nishioka and Sato

^{23}who found a critical value of

*Re*

_{b,cr}≈ 140 for a similar aspect ratio, α = 4.5. This counter-intuitive behavior indicates that there is a group of parameters in the range of

*Ha*<

*Ha*

_{cr}, and

*Re*and with increasing

*Ha*, the flow will go from a steady to an unsteady state. To verify the unexpected effect of generating unsteadiness by increasing

*Ha*, we have performed a series of additional numerical simulations for a few selected

*Re*values (within the shaded area), as opposed to the original simulations that were done for fixed

*Ha*values. For each

*Re*value, the new simulations were carried out by increasing

*Ha*in small steps, each time allowing the flow to reach a statistically steady state. By doing so, we have verified that when crossing the original stability boundary, as given in Figure 6(b), unsteadiness was obtained. Of course, when

*Ha*is increased significantly steadiness eventually prevails again.

To the best of our knowledge, this is the first time that such a behavior at low Hartmann numbers is observed, at least for the case of a wake flow with the magnetic field aligned with the cylinder axis. Other studies in the literature of similar flow configurations approached the problem by making use of the Q2D model,^{3,4,11} and as a result, were limited to relatively high Hartmann numbers (*Ha* ⩾ 320), where *Re*_{cr} is only found to increase with increasing *Ha*. Nonetheless, it is important to note that similar non-monotonic scenarios have been demonstrated in the past for other configurations, such as the case of a free cylinder with the magnetic field being parallel to the oncoming flow^{2} and the case of a duct flow.^{10} For example, Ref. 10, while investigating the case of a shear driven flow in a duct with non-uniform conductance of the Hartmann walls, observed a similar behaviour at moderate values of *Ha* and attributed that effect to the velocity deficit across the shear layers. Despite the fact that once again a Q2D model was used in his study, results were expected to apply qualitatively. Even though that study was addressing a different configuration, our results are qualitatively similar and follow similar physical explanations, as we shall see in the following paragraphs.

The overall behavior can be explained taking into account two physical aspects of the flow: how laminar vortex shedding is triggered in the wake of a circular cylinder and how the magnetic field affects the stability of the flow.

The first aspect has been addressed several times in the past and it is well-known: The onset of the von Kármán vortex flow results from a global Hopf bifurcation of the steady flow.^{24,25} The driving mechanism for this self-excited instability is associated with the inception of small disturbances located symmetrically at approximately the end of the separation bubble. In turn, these disturbances generate sinusoidal waves traveling against the basic flow on the sides of the recirculation region.^{26–28} If the blockage ratio is small and the aspect ratio is large, the presence of walls, either at the ends of the cylinder or parallel to the cylinder axis, does not appear to change the bifurcation properties in a significant manner. However, as the blockage ratio increases or the aspect ratio becomes smaller, the walls confine the streamlines near the cylinder, and the flow becomes more stable as the propagation of infinitesimal disturbances becomes suppressed.^{29–31} In addition to this effect, the local acceleration of the flow in the vicinity of the cylinder and along the centerline of the duct, presumably stabilizes the flow, leading to a delay of the onset of vortex shedding.^{12,31–33} From all the aforementioned, it is evident that the rectangular duct aspect and blockage ratios become critical parameters when considering the stability of the flow.

Concerning the second aspect, the magnetic field acts in two ways. On one hand, it increases the damping of the fluctuations of the unstable flow through Joule dissipation, as pointed out by Müller and Bühler.^{5} On the other hand, it modifies the velocity distribution of the primary flow, hence also acting on the flow stability through inertia.

For small Hartmann numbers, where *N* ∼ 1, viscous and Joule dissipation are of the same order of magnitude and overall dissipation is stronger than in the hydrodynamic case. This would suggest that the flow stability should increase with increasing *Ha*, however the current results show the opposite for *Ha* < 35. The unexpected destabilization of the flow with increasing *Ha* could be explained by the second action of the Lorentz force, viz. the modification of the velocity distribution before the cylinder, as shown in Figure 7. As already mentioned, the instability mechanism leading to the onset of vortex shedding is affected by local acceleration in the vicinity of the cylinder. Therefore, it is useful to look at the effect of *Ha* on the fluid velocity on the upstream side of the cylinder. As shown in Figure 6(b) (dotted line), with increasing Hartmann number, the centerline velocity *U*_{c} decreases rapidly at first, i.e., for low values of *Ha*, and approaches asymptotically the bulk flow velocity for larger values of *Ha*. The rapid decrease of the centerline velocity *U*_{c} at low *Ha* can be translated to a deceleration of the flow in the neighborhood of the cylinder, which effectively enhances instabilities and destabilizes the flow. This is in line with the explanation given by Bühler.^{10} In other words, with increasing *Ha*, the action of Lorentz forces leads to a more uniform velocity distribution along the cross-section of the duct (as shown in Figure 7), resembling a flow configuration with a higher aspect ratio, and smaller blockage ratio. This action mimics the effect of moving the walls of the duct further away from the cylinder. Consequently, having already explained the significance of aspect and blockage ratio on the flow stability, increasing *Ha* enhances the propagation of disturbances, which eventually destabilizes the flow leading to lower values for the critical Reynolds numbers. This effect becomes clearer when comparing the reduction of the flow stability at fixed bulk flow, and this explains the motivation for using *Re*_{b,cr} instead of *Re*_{c,cr} in Figure 6(b).

In the case of high Hartmann numbers, where *N* ≫ 1, the Lorentz forces become dominant, and the main effect of the magnetic field is the stabilization of the flow, by opposing vortical motions in planes lateral to the direction of the magnetic field. In this case, increasing *Ha* leads to a significant reduction of the energy amplification of disturbances due to Hartmann damping^{12} and the critical Reynolds number becomes proportional to *Ha*, as already pointed out in our results.

## V. GLOBAL FLOW FEATURES IN THE STEADY REGIME

In this section we present and discuss results from numerical simulations with time-independent flow, i.e. under parameters that produce steady mean flow. Results from the present three-dimensional simulations are compared against the numerical results of Dousset and Pothérat^{3} obtained using a quasi-two-dimensional model. It should be noted that in their study, they are reporting only results for high *Ha*, which is necessary for their model to hold, so a direct comparison is only possible for the cases where *Ha* ⩾ 320. An overview of the relative differences between the results obtained by the two studies (with respect to the current study) is given in Table V.

### A. Effect of *Ha* on the recirculation length *L*_{r}

The variation of the spanwise-averaged recirculation length, *L*_{r}, as a function of the Reynolds number,*Re*_{c}, and Hartmann number, *Ha*, is shown in Figure 8. In Figure 8(a), an almost linear evolution of *L*_{r} with *Re*_{c}, for each *Ha*, is observed. In the case of low Hartmann numbers, up to *Ha* = 80, the magnetic field seems to have little effect on *L*_{r}, and all data collapse onto a single curve that depends only on *Re*_{c}. Increasing *Ha* further causes the data to depart from that curve, while the slope slightly decreases. As shown in Figure 8(a), an increase of *Ha* at a fixed value of *Re* causes a decrease in *L*_{r}. This is attributed to the action of the Lorentz forces as we shall discuss in more detail in the following paragraph. Interestingly, if we plot *L*_{r} as a function of *Re*_{c}/*Ha*^{0.8}, as shown in Figure 8(b), a universal scaling law appears, this time for high Hartmann numbers, *Ha* ⩾ 320. This collapse for high *Ha*, was first demonstrated by Dousset and Pothérat,^{3} with whom the current results agree very well. As seen from Table V, the relative differences between the two studies are around 3%.

In addition to the spanwise-averaged recirculation length, *L*_{r}, it is also useful to look at the distribution of the recirculation length along the span, *L*_{r}(*z*), as shown in Figure 9. This allows us to analyze in greater detail the evolution of the recirculation length with increasing *Ha*. In this plot, we compare *L*_{r}(*z*) for different *Ha*, while keeping the flow rate constant, i.e. *Re*_{b} constant. In the case of low *Ha* values [see Figure 9(a)], we can see that, as we move away from the centerline and closer to the Hartmann walls, there is a large discrepancy between the values of *L*_{r}(*z*). Contrary to the previous observation, where the magnetic field was revealed to have little effect on the *average* value of the recirculation length, here we can see that an increasing magnetic field induces significant differences in the distribution of the recirculation length. Particularly at *Ha* = 20, besides the peak reached by *L*_{r}(*z*) at midspan, two local peaks exist near the Hartmann walls. As recently explained by Dousset and Pothérat^{34} while considering the case of a *truncated* rectangular cylinder, these two peaks are due to the presence of a pair of counter-rotating streamwise vortices, referred to as base vortices, located in the vicinity of each Hartmann wall downstream of the cylinder. The strength of these vortices is diminished with decreasing boundary layer thickness upstream of the cylinder.^{35} Therefore, with the increase of *Ha*, as the Hartmann boundary layer becomes rapidly thinner (δ_{H} ∼ 1/*Ha*), the base vortices are weakened, which translates into the damping of the peaks near the walls of the duct. In addition to the weakening of the base vortices, Lorentz forces act in the direction opposite to the flow in the core region of the duct, while diffusing vorticity along the magnetic field. These effects explain why increasing the intensity of the magnetic field for a given mass flow rate leads towards lower values of *L*_{r} and smaller variations along the span. At higher Hartmann numbers [see Figure 9(b)], Lorentz forces are enhanced, and as a result, the recirculation region displays increasingly a relatively flat distribution, i.e. the flow is driven towards a two-dimensional flow.

### B. Evolution of the base pressure coefficient *C*_{pb}

Figure 10, shows the evolution of the base pressure coefficient, *C*_{pb}, with respect to *Ha* and *Re*_{c}. Dousset and Pothérat^{3} demonstrated that for high Hartmann numbers, *Re*/*Ha* is a governing parameter for which all data collapse on a single curve. In this study, we observe that the *Re*/*Ha* universality holds also for all the low values of *Ha* investigated here for the first time. Even though at low Hartmann numbers a shift towards higher values of *Re*/*Ha* is observed, the slope remains the same and *Re*/*Ha* can be considered as a governing parameter of *C*_{pb}, for the entire range of parameters studied here. As shown in Table V, when the current results are compared with the results obtained with the Q2D model, once again, a remarkable agreement between the two studies can be observed. The relative difference between the two studies ranges between 4% and 5% for *Ha* = 320 and 640, and drops to 1% for *Ha* = 1120. This demonstrates that the Q2D model produces a very accurate averaging of the physical phenomena along the spanwise direction, provided that the magnetic field is sufficiently strong for the model to hold.

The spanwise distribution of the base pressure coefficient, *C*_{pb}(*z*), normalized with its spanwise average value can be seen in Figure 11. In contrast to the spanwise distribution of the recirculation length, where a stronger magnetic field leads to a flatter distribution of *L*_{r}(*z*) (see Figure 9), the base pressure coefficient is found to become less uniform for higher values of *Ha*. As the Hartmann number increases, the pressure progressively drops at the mid span and increases close the side walls. This unexpected trend is linked to the three-dimensional distribution of the Lorentz forces along the surface of the cylinder. For example, Figure 12 shows the current density streamlines, passing close to the rear part of the cylinder surface at *Ha* = 1120. The surface of the cylinder is colored by the crossflow component of the current density, *J*_{y}. As clearly seen, due to the electrically insulating flow insert, the current density paths are deflected, generating a three-dimensional distribution of currents around the cylinder. Along the rear side of the insert, this current pattern becomes more complex, with a region of negative values of the cross-flow component of current density *J*_{y} appearing between regions of positive values. In turn, this non-uniform distribution of the current density field at the rear part, induces Lorentz forces which act in the positive streamwise direction at both ends of the cylinder and towards the upstream direction at the core region of the flow. These non-uniform forces effectively redistribute the pressure field at the rear part of the cylinder leading to the “three-dimensionality” of the base pressure coefficients as shown in Figure 12. As the Hartmann number is increased, Lorentz forces become stronger, and the variation of the base pressure coefficient along the spanwise direction is enhanced. Nevertheless, it is important to note that the spanwise variation of *C*_{pb}(*z*) is relatively small. Even at the highest Hartmann number considered in this investigation, *Ha* = 1120, the relative difference between the minimum and maximum values of *C*_{pb}(*z*) is found to be around 10% with respect to the minimum value.

### C. Influence of *Ha* on the drag coefficient *C*_{D}

Figure 13 shows the evolution of the drag coefficient, *C*_{D}, as a function of *Re*_{c}, and *Ha*. Focusing at the range of high *Ha* values in this figure (*Ha* ⩾ 320), a universal scaling law is found to hold for the drag coefficient as a function of *Re*/*Ha*^{0.8}, as was the case with *L*_{r}. While the same slope is valid for values of *Ha* < 320, results in this case are shifted relative to the curve followed at high *Ha*. Here again, the agreement between the current three-dimensional simulations and the Q2D model, is very good, as shown in Table V. Nevertheless, the discrepancy between the two studies is found to be larger for *C*_{D} values than for *C*_{pb}, and this difference is explained in the following paragraphs.

The distribution along the cylinder span, of the drag coefficient, *C*_{D}(z), normalized with its spanwise average is shown in Figure 14, where different Hartmann numbers are compared at *Re*_{b} = 50. Similarly to the base pressure coefficient, an increasing variation of the drag coefficient along the spanwise direction is observed with increasing *Ha*. Indeed, *C*_{D}(*z*) displays a rather flat distribution at low Hartmann numbers and gradually becomes highly non-uniform as *Ha* is increased. This evidence shows again how the budget of forces becomes more “three-dimensional” with increasing *Ha* due to the stronger action of the Lorentz forces. The physical mechanism underlying such a non-uniformity can be better understood by examining the spanwise variations of the pressure and viscous components of the total drag, as shown in Figure 15 for selected values of Hartmann numbers. In all cases, the magnitude of the total drag is dominated by the pressure forces around the surface of the cylinder, which on average correspond to about 80% of the total. At *Ha* = 20, the influence of Lorentz forces on the pressure drag is small with the latter displaying a relatively flat profile, therefore the spanwise variation of *C*_{D}(*z*) is due almost entirely to the effect of viscous stresses, that are obviously affected by the wall-boundary layer at both ends of the cylinder. At *Ha* = 1120, owing to the stronger magnetic field, the Lorentz forces affect considerably the distribution of the pressure forces, as previously explained. As a result, the contribution of the pressure drag to the spanwise variation of the total drag, outweighs significantly the viscous forces. Interestingly, unlike *C*_{pb}(*z*), the relative difference between the minimum and maximum value of *C*_{D}(*z*) for *Ha* = 1120, is almost 200% with respect to its minimum value. This is attributed to the fact that *C*_{pb}(*z*) is simply a point measurement, while *C*_{D}(*z*) is an integral quantity of both the viscous and pressure forces over the surface of the cylinder and therefore more indicative of the Lorentz forces effect. From an engineering point of view, the increased non-uniformity of the spanwise distribution of *C*_{D} is important since it can lead to uneven stress distributions. Yet, this unexpected trend could not have been inferred by examining the velocity distribution, which on the contrary becomes more two-dimensional with increasing *Ha*.

In view of this strikingly “three-dimensional” distribution of the drag coefficient, one would have expected that the Q2D-calculated value for the average *C*_{D} would have differed significantly from the result produced by this fully three-dimensional simulation. However, as shown in Table V, the discrepancies between the two studies are quite small, and not significantly larger than those for the *C*_{pb}. This fact suggests that for steady-flow configurations, the averaging of the three-dimensional governing equations – as done by the Q2D model – can capture the integral effects of the MHD flow, thus producing sensible and reasonably accurate values for the averaged parameters, such as *C*_{D} and *C*_{pb}.

## VI. FLOW CHARACTERISTICS OF THE UNSTEADY REGIME

The third and final part of the current study concerns the unsteady flow regime. Within this regime very rich phenomena appear, which are of interest partially due to the antecedents in the literature. As already mentioned in the introduction, the results of Frank *et al.*,^{11} and Dousset and Pothérat,^{3} revealed the presence of an additional flow regime that comes in succession after the laminar periodic vortex shedding. This regime is characterized by irregular vortex patterns and by the shedding of secondary vortices from the top and bottom walls, and it does not have a counterpart in the purely hydrodynamic case at such moderate blockage ratio.^{30}

The aim of this section is to verify the presence of this regime and to describe additional flow features that can only be observed using a full three-dimensional analysis. Furthermore, the analysis includes a discussion of unsteady regimes at low interaction parameters, *N* ∼ 1, never addressed before.

Due to the increased computational time needed for such three-dimensional flow simulations, and especially the need to integrate to a sufficiently long physical time to achieve a statistically steady-state, only a limited number of simulations are carried out. The choice of simulations is guided by the current results from Sec. IV, which show when regular vortex shedding starts to take place, and by the parametric study of Dousset and Pothérat,^{3} which indicates when the newly identified flow regime occurs. Therefore, a total of three simulations are carried out, namely cases *U*1, *U*2, and *U*3, using grids GL-M, GF-M, and GX-M, respectively. Case *U*1 is performed, at *Ha* = 320, *Re*_{c} = 2000, *N* ≃ 3, matching the von Kármán vortex flow. The other two simulations are performed for a fixed higher Reynolds number (*Re*_{c} = 5000) while varying the Hartmann numbers. Thus, case *U*2 corresponds to moderate (*Ha* = 320) and case *U*3 to strong magnetic field intensity (*Ha* = 1120). These sets of parameters correspond to the forth flow regime according to Ref. 3 and yield interaction parameters of *N* ≃ 1 and *N* ≃ 16, respectively, for cases *U*2 and *U*3.

The results obtained in the unsteady flow regime are summarized in Table III, while in Table VI they are compared to the data reported in Ref. 3. For case *U*3, results between the two studies are in good agreement, however for case *U*1 large discrepancies can be observed. These can be attributed to the low interaction parameter associated with case *U*1, as the Q2D model used in Ref. 3 is likely to produce invalid results at this low interaction parameter, as pointed out by the authors.^{3}

### A. Case *U*1

Figure 16 shows snapshots of contour plots of spanwise vorticity, ω_{z}, in the center-plane *z* = 0, for the case *U*1 (*Ha* = 320, *Re*_{c} = 2000). In this figure, a typical von Kármán vortex street is observed. Vortices are shed in a regular manner from alternate sides of the cylinder, forming two rows that move alternately clockwise and anticlockwise. Nevertheless, when the vortex trajectories are contrasted to the corresponding confined case of flow over a cylinder in a *channel*,^{21,33,36} a fundamental difference between the two cases can clearly be observed. Unlike the channel flow case, where the trajectories of the two rows of vortices intersect and their position with respect to the symmetry line is inverted,^{21,33,36} in this case no inversion of the shed vortices is observed. As reported in previous studies,^{5,11,37} under the action of the magnetic field the wake of the cylinder is narrowed and the distance between the vortices is reduced. Consequently, shed vortices reach the vicinity of the walls only further downstream. As the inversion of the von Kármán vortices is caused by interactions with the top and bottom walls soon after they are shed,^{33,36} no inversion of vortices is observed in this case.

### B. Case *U*2

In case *U*2, inertial effects are amplified with respect to case *U*1, by increasing the Reynolds number to *Re*_{c} = 5000, while keeping the Hartmann number constant at *Ha* = 320. As seen in Figures 17(a) and 17(b), vortex shedding remains periodic, but appears irregular. Looking at the different instances of the flow, the wake still consists of an arrangement of vortices shed from the cylinder, however the wake is further narrowed and loses its symmetry. Shed vortices maintain their vorticity for longer distances, follow an irregular path and their paths intersect crossing each other. It is interesting to note that in the near wake region, *x* < 5*D*, generated vortices are fragmented and small-scale structures are observed. This can be seen more clearly in Figure 17(c), where a closer look into the near wake region of Figure 17(a) is presented. Here, the primary vortices appear heavily distorted, and small structures appear in the separated shear layers. Moving further downstream, these smaller structures are eventually damped, and one finds that the structure of the wake returns to a more coherent state.

The visualization of flow structures in this rather complicated three-dimensional wall-bounded shear flow can be facilitated using the λ_{2} criterion. The λ_{2} is defined as the second eigenvalue of *S*^{2} + *X*^{2}, where *S* and *X* denote the symmetric and antisymmetric parts of the velocity gradient tensor, respectively.^{38} Plotting iso-surfaces of λ_{2} provides an effective way to visualize the primary and induced vortices of the cylinder wake, while excluding the wall shear region.

Figure 18 shows snapshots of iso-surfaces of λ_{2} normalized by its absolute minimum, λ_{2, min}, for the case *U*2. Iso-surfaces are colored by the spanwise vorticity component, ω_{z}, in order to reveal the spanwise rotation direction of each vortical structure. The spanwise rollers essentially identify the primary vortex cores. At these values of *Ha* = 320 and *Re*_{c} = 5000, the corresponding interaction parameter is *N* ≃ 1, and the flow is expected to approach a two-dimensional state according to Mück *et al.*^{13} This can be verified from Figure 18, if we focus our attention further downstream. Indeed, from approximately *x* ≃ 5*D* and further downstream, the primary vortex cores are aligned in the spanwise direction, parallel to the magnetic field direction, although they display a slight curvature in the bulk of the flow. Only in the very thin viscous Hartmann layers close to the side walls, three-dimensional effects are visible. This tendency towards a two-dimensional flow is in very good agreement with the observations of Mück *et al.*^{13} However, in the near wake region the flow is strikingly more complex and highly three-dimensional. The flow in the separated shear layers is highly fragmented and hosts a plethora of small-scale vortices. Inevitably, as these small turbulent-like structures travel downstream, they become suppressed by the magnetic field and eventually are completely damped.

The structure of the near wake can be attributed to the presence of the circular cylinder, which acts as a turbulence promoter, thereby facilitating the generation of three-dimensional flow structures (especially at this high Reynolds number). Looking at a collection of various snapshots (video 1), the small-scale vortices in the near wake seem to arise from instabilities originating at the shear layers, before mixing with the primary Kármán vortices. This development of three-dimensional structures on the scale of the shear layer vortices bears resemblance to flow characteristics during the initial stages of the shear-layer transition regime, as described in several purely hydrodynamic studies for the case of a free cylinder at similar Reynolds numbers.^{39–41} In hydrodynamic scenarios, the shear layers separating from the sides of the cylinder eventually become turbulent. Of course, the full progression of this phenomenon cannot be seen here due to the damping action of the magnetic field.

One could also argue that the irregular character of vortex shedding resembles the new flow regime reported in previous studies for similar values of *Ha* and *Re*_{c}.^{3,11} For this range of parameters, Ref. 3 reported irregular shedding of primary vortices that interacted strongly with the boundary layers at the walls leading to shedding of secondary vortices from the top and bottom walls. However, this is not supported by the present results. Despite the irregular shedding of primary vortices seen in case *U*2, the wake remains narrow and the interaction between the primary vortices and the Shercliff layers at the top and bottom walls of the duct is weak. As a consequence, the flow fails to produce boundary layer separation at the duct walls and, eventually, no shedding of secondary vortices from the top and bottom walls is observed. Secondary vortices shed from the top and bottom walls, represent a primary feature of the fourth flow regime reported in Ref. 3, and therefore a main difference between the two analyses. It is nevertheless important to underline the fact that Dousset and Pothérat^{3} suggested that for such a low interaction parameter, i.e., *N* ≃ 1, three-dimensional effects are likely to appear, therefore rendering the validity of the Q2D model questionable, which in turn raises questions about their description of the flow regime. Indeed, as shown in Figure 18, the current three-dimensional computations show that, for these set of parameters, three-dimensional structures are a salient feature of the flow. Current results reveal that the Q2D model is not able to capture these specific three-dimensional effects, as they are caused by small and turbulent-like events that cannot be accurately captured by the averaging of the equations. A more detailed discussion concerning these differences and the respective physical reasons behind them will be addressed in Subsection VI E.

### C. Case *U*3

Case *U*3 corresponds to the same Reynolds number (*Re*_{c} = 5000) as case *U*2, but at a higher Hartmann number value of *Ha* = 1120. Intuitively, under these conditions of stronger electrodynamic damping, one would expect the shed vortices to be more suppressed and to diffuse more rapidly as they convect downstream, thus reducing the irregularity of the overall vortex shedding pattern. Surprisingly, as seen in Figure 19, for this set of parameters the shedding pattern of the vortices is actually even more irregular than the one obtained in the previous case (*U*2) at a lower Hartmann number (*Ha* = 320, Figure 17). In the near wake of the cylinder, i.e., *x* < 10*D*, the lateral spacing between shed vortices is increased, with the vortices following an irregular oblique trajectory closer to the top and bottom walls. As they are progressively convected further downstream, i.e., *x* > 10*D*, they interact with the duct walls inducing secondary vortices from the top and bottom walls (see vortex core *WV* in Figure 19) and inevitably are reflected back in the opposite direction (see vortex cores *KV*_{1}, *KV*_{2}). The induced *secondary*vortices are, however, weak and are quickly dissipated soon after they are shed. Nevertheless, their influence on the overall phenomenology is significant. Another interesting feature of the flow is the appearance of secondary vortices (see vortex core *SV*) alongside with the main vortices. At times, a primary vortex core, soon after it is shed, deforms, and eventually is torn apart, giving birth to a smaller secondary vortex. This intermittent secondary vortex, interacts with its surrounding vortices and eventually disrupts significantly the trajectories of the primary vortex cores. This effect occurs spontaneously from either side of the cylinder in a consistent and systematic way. The formation of the secondary vortex can be clearly seen in Figure 19(e), where a close-up view of Figure 19(a) in the vicinity of the cylinder is presented. Additionally, in this figure one can identify shed vortices that appear less fragmented and more coherent at the shear layers.

This is also evident in Figure 20, where iso-surfaces of λ_{2}, normalized by its absolute minimum, λ_{2,min}, are presented for the same case. In this figure, the stronger magnetic damping (as compared to the previous case *U*2, shown in Figure 18) is obvious. Close to the cylinder for example, the shear layers are devoid of small three-dimensional structures and more coherent vortices are formed. Only a weak undulation of the primary vortices is observed. Further downstream, vortex rollers are almost perfectly aligned along the spanwise direction. These effects can clearly be attributed to the flattening of the velocity profiles along the spanwise direction as shown before, and as expected due to the high interaction parameter *N* ≃ 16. It is interesting to note that for this Hartmann number the vortex alignment highly resembles that of a flow configuration without side-walls or even of unconfined cases. The current flow phenomenology for case *U*3 is in better agreement with the results of Ref. 3, although the flow is still clearly less irregular. In order to understand the differences and the physical reasons behind them, a detailed spectral analysis is incorporated next in the discussion.

### D. Effect on lift and drag coefficient

The regular or irregular nature of the flow is reflected by the time history of the lift, *C*_{L}, and drag, *C*_{D}, coefficients as shown in Figure 21. For case *U*1, the perfectly time-periodic nature of the flow is displayed in Figure 21(a). The irregular vortex shedding observed in cases *U*2 and *U*3, is clearly depicted on the time history of the lift and drag coefficient, as seen in Figures 21(b) and 21(c). Both signals keep a periodic character, however, lose their coherence.

### E. Spectral analysis

Another powerful way to identify the occurrence of vortex structures, and investigate the presence of irregular flow phenomena, is through a spectral analysis of velocity signals gathered from probes placed at different streamwise positions. Figure 22 displays the power spectral density of the crossflow velocity, *u*_{y}, obtained from probes positioned at the mid-plane (*z* = 0) at several (*x*, *y*) locations: in the vicinity of the cylinder's shear layer (*x* = 1.5*D*, *y* = 0.5*D*), at the center downstream (*x* = 10*D*, *y* = 0) and further downstream (*x* = 20*D*, *y* = 0). Power spectral densities have been obtained using Welch's averaged periodogram method,^{42} while a Hamming window was applied to each overlapping segment of data. Results are shown for all cases examined in this subsection, namely cases *U*1, *U*2, and *U*3.

For the case *U*1 (*Ha* = 320, *Re*_{c} = 2000), the power spectrum exhibits the expected strong peak at the fundamental frequency of the flow (Strouhal frequency, *St*), corresponding directly to the passage of a laminar periodic vortex array through the measuring probe. The rest of the peaks observed represent harmonic frequencies (in the centerline only odd harmonics are observed), in line with the observations of Kovasznay^{43} for hydrodynamic flow over an unconfined cylinder in the laminar flow regime.

In the case *U*2 (*Ha* = 320, *Re*_{c} = 5000), the strong presence of the Strouhal frequency in the spectra can still be observed, suggesting that a vortex street continues to be present. However, the spectrum close to the cylinder, at *x* = 1.5*D*, is broadened indicating the presence of smaller scales of motion. The scaling of the energy spectrum in the inertial region is of the form *f*^{−3}, which is characteristic of MHD turbulence.^{9,44} This broad range of frequencies illustrates the development of finer-scale three-dimensional structures in the shear layers, in line with our previous observations (see Figure 17). Looking at the spectra obtained from the probes positioned further downstream, at *x* = 10*D* and *x* = 20*D*, it is evident that the small-scale fluctuations are slowly suppressed, and the spectra are dominated by a primary frequency of oscillation and its harmonics. The weakly turbulent behaviour of the flow gives place to larger, regular, and smoother vortices that correspond to the primary vortices. Nevertheless, non-harmonic components are also present due to an irregular periodic signal, supporting the idea of the irregular vortex shedding process described before.

Finally, in the case *U*3 (*Ha* = 1120, *Re*_{c} = 5000) the spectra exhibit different features. Close to the cylinder, at *x* = 1.5*D*, the spectrum is characterized by fewer high frequencies and harmonics of the primary shedding frequency can clearly be identified. This spectral feature support the visualizations of the flow shown in Figure 20, where the cylinder wake was found to be dominated by coherent vortices with shear layers devoid of small-scale three-dimensional structures. More importantly, when compared to case *U*2, the energy content of low frequencies is strongly enhanced (by almost an order of magnitude) indicating the presence of more energetic primary vortices in the flow. This is associated with the irregular character of the flow at this set of parameters, as will be explained in the following paragraph. Irregular vortex shedding is also supported, by looking at the spectra further downstream, at *x* = 10*D* and *x* = 20*D*. Although a strong peak is present indicating periodic vortex shedding, the absence of clearly identifiable harmonic frequencies, together with the broadening of the power spectrum is linked with the presence of irregular vortex flow patterns. Furthermore, the power density contained in the higher frequencies has been greatly reduced, indicating that no turbulent-like phenomena are present. This complete suppression of turbulence is of course not surprising, given the stronger magnetic field.

The presented spectral analysis, in combination with the flow visualizations of cases *U*2 and *U*3 (shown in Subsections VI B and VI C, respectively), allow us to finally propose a mechanism explaining the surprising enhancement of flow irregularity with increasing Hartmann number and the differences observed between the current results and quasi two-dimensional studies. This explanation is based on the MHD modification of small-scale three-dimensional motions. At low interaction parameters, *N* ≃ 1, the combination of strong inertial forces due to the high Reynolds number (*Re*_{c} = 5000) and of a flow obstruction promotes transitional or even turbulent characteristics near the cylinder. The presence of small scale three-dimensional structures causes very strong Joule dissipation, which eventually extracts energy from the large-scale primary vortex cores. Through this process the primary vortices become less energetic, leading to a more regular vortex shedding pattern. For a given Reynolds number, increasing the Hartmann number suppresses the intensity of the smallest three-dimensional flow structures due to increased Hartmann damping. Consequently, increasing *Ha* results in more energetic primary vortices, which retain most of their inertia since less energy is transferred to smaller scales of motion. This could explain why for the higher *Ha* number considered (case *U*3), the vortex cores start interacting with the top and bottom walls, thus inducing secondary vorticity and leading to a more irregular flow pattern. The Q2D model fails to capture the enhancement of the Joule dissipation due to the presence of highly three-dimensional structures, and as a result, over-predicts the inertia retained in the primary vortex cores, especially at these high Reynolds numbers. This phenomenology can explain why the Q2D model yields more irregular flow characteristics than those observed in the current three-dimensional analysis for the same set of parameters.

Nevertheless, additional work would be required in order to clarify the trends at *Ha* values even higher than those considered here. A further increase of the Hartmann number, could damp even more the three-dimensional effects, yielding a more irregular flow regime, closer to the one identified by the Q2D approach. On the other hand, the increased magnetohydrodynamic damping could be sufficiently strong to suppress vortex shedding, stabilize the flow, and eventually lead to a steady flow regime. These considerations require further extensive simulations that would help understand the complex interaction between the inertial three-dimensional effects, the mechanisms of magnetic dumping, and the way they affect the generated flow patterns.

## VII. CONCLUSIONS

Direct numerical simulations of three-dimensional flow of a liquid metal around a circular cylinder placed in a rectangular duct, under a wide range of magnetic field intensities have been performed. The blockage and aspect ratio were kept constant at β = 1/4, α = 4. The Reynolds number was varied between 50 ⩽ *Re*_{c} ⩽ 5000 for Hartmann numbers in the range of 0 ⩽ *Ha* ⩽ 1120.

A non-monotonous behavior of the critical Reynolds number for the onset of vortex shedding with respect to the Hartmann number has been shown. For *Ha* > 80, a linear dependence of the critical *Re* is observed with respect to *Ha*. However, when decreasing the Hartmann number below *Ha* = 80, a non-monotonous relation is observed, where initially *Re*_{cr} decreases slightly until it reaches a minimum value around *Ha* ≃ 35. As *Ha* is decreased even further, *Re*_{cr} starts to increase again. This unexpected phenomenon is explained by the effect of the Lorentz forces on the velocity profile, which becomes flatter for increasing Hartmann numbers. A flatter velocity profile implies a reduction of the velocity deficit for a given bulk flow, which essentially mimics the effect of moving the walls further away, therefore, promoting the destabilization of the flow. This effect could be of practical importance in various engineering applications where enhancement of heat and mass transfer or mixing is desired.

Within the steady flow regime, the evolution of the recirculation length and of the drag and base pressure coefficients with respect to *Ha* and *Re* was presented. The predictions of the Q2D model for these global flow characteristics are in good agreement with the current three-dimensional numerical results when averaged over the cylinder length. However, the current three-dimensional computations display an increase in the three-dimensionality of the spanwise distribution of forces along the cylinder with increasing *Ha*, an effect that cannot possibly be captured in the Q2D analysis. Hence, depending on the range of parameters, an accurate stress analysis of the cylinder, and presumably of other flow inserts, could require a fully three-dimensional treatment.

We have also showed that for a high Reynolds number,*Re*_{c} = 5000, even for relatively high values of the interaction parameter (*N* up to 16), inertia effects in the near wake play a significant role in the formation and transport of vortices. For *Ha* = 320 (*N* ≃ 1), inertial forces lead to turbulence-like instabilities, which redistribute the energy towards small scales, eventually reducing the energy of the shedding vortex cores. When the Hartmann number is increased to *Ha* = 1120 (*N* ≃ 16), the presence of a stronger magnetic field acts like a turbulent suppressor, hence primary vortices retain most of the inertia and remain more energetic as they move downstream. While the irregular character of the flow at *Ha* = 1120 has been described by Dousset and Pothérat,^{3} the direct comparison of the flow at *Ha* = 1120 and *Ha* = 320 leads to the surprising new observation of enhanced flow irregularity with increasing *Ha*. This effect also explains the differences observed between the current three-dimensional results and the Q2D model predictions. The quasi-two-dimensional model cannot capture the enhancement of the Joule dissipation due to the presence of highly three-dimensional structures in the flow. As consequence, it overestimates the energy retained in the primary vortex cores and predicts a much more irregular vortex shedding for a given set of parameter values than what has been observed in the present work. Nevertheless, it is reasonable to expect this behavior to be non-monotonic, as for even higher magnetic fields the flow should become steady again, since the energy of the large scales would eventually be also damped by stronger magnetic fields.

## ACKNOWLEDGMENTS

The authors would like to thank Professor A. Pothérat for kindly making available Q2D results from Ref. 3. This work has been performed under the UCY-CompSci project, a Marie Curie Transfer of Knowledge (TOK-DEV) grant (Contract No. MTKD-CT-2004-014199), funded by the CEC under the 6th Framework Program, and under the contract of association (Contract No. ERB 5005 CT 99 0100) between the European Atomic Energy Community and the Hellenic Republic. Partial support has also been received through a Center of Excellence grant from the Norwegian Research Council to the Center for Biomedical Computing.

## REFERENCES

_{D}= 3900