Three-dimensional numerical simulations of magnetohydrodynamic flow around a confined circular cylinder under low, moderate, and strong magnetic fields

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²/ρ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²/Re, where

$Ha=LBσ/ρν$ is the Hartmann number, the square of which characterizes the ratio of electromagnetic to viscous forces. Any movement of a conducting fluid under the influence of a magnetic field, causes an electric field to develop and as a result electric currents are generated. In turn, the interaction of these currents with the imposed magnetic field gives rise to an MHD force, known as the Lorentz force.^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.⁷ 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.,⁴ 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.¹¹ and numerical simulations from Dousset and Pothérat³ 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.³ Again, as the imposed magnetic field is increased, transition to this regime is delayed. Mück et al.¹³ 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.,¹¹ 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.¹³ 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,³ 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.

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.³ 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.5D upstream of the circular cylinder, while the outlet is located at L_o = 35.5D 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₀ 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.¹⁴ 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², σU_cB₀, B₀, and DU_cB₀, 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.5D, 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.

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.⁵ A no-slip boundary condition is imposed on the cylinder surface and the duct walls,

$uwall=0$⁠.

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_jn_jds/∫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.¹³ This boundary condition reduces to the requirement of zero normal component of the current at the inlet and outlet.

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.¹⁷ A very detailed description of this scheme and other numerical techniques used in the MHD module is reported in Albets-Chico et al.¹⁸ The implementation and integration of this module in CDP has been validated for several MHD flows and in various geometries, including MHD channel flow¹⁹ and MHD pipe flow.²⁰ 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.

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.⁵ 

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²ν⁻¹ρ⁻¹( − ∂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.

Results are compared with the analytical solution in terms of a series expansion provided by Müller and Bühler.⁵ 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.

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,²¹ 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.005D and 0.002D, 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 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

$C¯D$⁠, recirculation length

$\bar{L}_r$

$L¯r$⁠, base pressure coefficient

$\bar{C}_{pb}^{\prime }$

$C¯pb′$⁠, and the Strouhal number 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,

$C_{pb}^{\prime }$

$Cpb′$ in which the reference pressure is at the front stagnation point of the cylinder (0° from the front). Results show that a domain independent solution can be achieved with the domains considered. For example, replacing the shortest domain, GL-S, with the longest one, GL-L, resulted in only 2.3% change in

$\bar{L}_r$

$L¯r$⁠, and 0.3% in

$\bar{C}_{pb}^{\prime }$

$C¯pb′$⁠, whereas the corresponding changes were 1.9% and 0.2% when replacing the intermediate domain GL-M with the longest domain. In the case of

$\bar{C}_D$

$C¯D$ and 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 2D from the outlet. Based on these results, one can conclude that the intermediate domain GL-M is sufficiently long to resolve the flow.

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.2D to 0.1D. 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⁶ and 10.5 × 10⁶ 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.5D, y = 0.5D, 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.²² 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.

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 numberRe_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³ 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 numberRe_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.52Ha. This relation is in very good agreement with the results of Dousset and Pothérat,³ who report a relation of Re_cr ≃ 0.43Ha. A similar linear relationship between Re and Ha was also observed by Frank et al.¹¹ (Re_cr ≃ 0.47Ha) and Hussam et al.⁴ (Re_cr ≃ 0.5Ha) 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,

$Reb,crmin≈90$⁠, in the vicinity of 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²³ 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_{cr}^{min}<Re<Re_{cr}(Ha)$

$Recrmin<Re<Recr(Ha)$ [shaded area in Figure 6(b)], in which for a fixed value of 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² and the case of a duct flow.¹⁰ 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.⁵ 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.¹⁰ 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¹² and the critical Reynolds number becomes proportional to Ha, as already pointed out in our results.

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³ 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.

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,³ 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³⁴ 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.³⁵ 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.

Figure 10, shows the evolution of the base pressure coefficient, C_pb, with respect to Ha and Re_c. Dousset and Pothérat³ 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.

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.

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.,¹¹ and Dousset and Pothérat,³ 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.³⁰ 

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,³ which indicates when the newly identified flow regime occurs. Therefore, a total of three simulations are carried out, namely cases U1, U2, and U3, using grids GL-M, GF-M, and GX-M, respectively. Case U1 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 U2 corresponds to moderate (Ha = 320) and case U3 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 U2 and U3.

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 U3, results between the two studies are in good agreement, however for case U1 large discrepancies can be observed. These can be attributed to the low interaction parameter associated with case U1, 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.³ 

Figure 16 shows snapshots of contour plots of spanwise vorticity, ω_z, in the center-plane z = 0, for the case U1 (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.

In case U2, inertial effects are amplified with respect to case U1, 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 < 5D, 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 λ₂ criterion. The λ₂ is defined as the second eigenvalue of S² + X², where S and X denote the symmetric and antisymmetric parts of the velocity gradient tensor, respectively.³⁸ Plotting iso-surfaces of λ₂ 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 λ₂ normalized by its absolute minimum, λ_{2, min}, for the case U2. 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.¹³ This can be verified from Figure 18, if we focus our attention further downstream. Indeed, from approximately x ≃ 5D 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.¹³ 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 U2, 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³ 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.

Case U3 corresponds to the same Reynolds number (Re_c = 5000) as case U2, 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 (U2) at a lower Hartmann number (Ha = 320, Figure 17). In the near wake of the cylinder, i.e., x < 10D, 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 > 10D, 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₁, KV₂). The induced secondaryvortices 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 λ₂, 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 U2, 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 U3 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.

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 U1, the perfectly time-periodic nature of the flow is displayed in Figure 21(a). The irregular vortex shedding observed in cases U2 and U3, 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.

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.5D, y = 0.5D), at the center downstream (x = 10D, y = 0) and further downstream (x = 20D, y = 0). Power spectral densities have been obtained using Welch's averaged periodogram method,⁴² while a Hamming window was applied to each overlapping segment of data. Results are shown for all cases examined in this subsection, namely cases U1, U2, and U3.

For the case U1 (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⁴³ for hydrodynamic flow over an unconfined cylinder in the laminar flow regime.

In the case U2 (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.5D, 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⁻³, 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 = 10D and x = 20D, 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 U3 (Ha = 1120, Re_c = 5000) the spectra exhibit different features. Close to the cylinder, at x = 1.5D, 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 U2, 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 = 10D and x = 20D. 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 U2 and U3 (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 U3), 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.

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,³ 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.

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.