Effects of symmetry on magnetohydrodynamic mixed convection flow in a vertical duct

The high boiling point and high thermal conductivity of liquid metals make them attractive as cooling fluids in technological applications with very strong heat fluxes. A well-known example of such an application is in the advanced nuclear fission reactors.¹ Our work is directly relevant to another area, namely, design of liquid metal blankets of magnetic-confinement nuclear fusion reactors. It is foreseen that such blankets serving the triple roles as heat exchangers, shields, and tritium breeders will be essential components of future reactors. The versions of the blanket consisting of complex systems of ducts and manifolds, in which a liquid metal (most likely a PbLi alloy) circulates, are believed to be the most promising candidates for this role.² 

The use of liquid metals in blankets leads to problems not typical for other types of coolants (gases, water, and salt melts): corrosion of walls and strong magnetohydrodynamic (MHD) resistance to the flow. It has been recently realized that yet another effect, namely, that of thermal convection occurring in the conditions of strong nonuniform heating and very strong imposed magnetic field, is an essential factor of blanket design. Experimental and numerical studies, such as Refs. 3–10, have shown that when conventional turbulence is suppressed by the magnetic field, the thermal convection affects the flow behavior in a profound, complex, and counter-intuitive way. In particular, convection instabilities often lead to quasi-regular high-amplitude fluctuations of temperature, which, if they appear in a fusion reactor, will completely change the blanket’s operation and, possibly, threaten its structural integrity via high-amplitude unsteady thermal stresses in its walls.

The magnetoconvection (this term is commonly used to describe the thermal convection with a strong effect of a magnetic field on the flow) brings on additional factors to be taken into account, such as the geometry of the channels, their orientation relative to gravity, magnetic field, and heat flux from the reactor chamber, and the structural features, such as the presence of auxiliary cooling circuits and electric conductivity of walls. It must also be stressed that a purely experimental study of the phenomenon is impossible. At present, there are no experimental facilities capable of recreating the extreme conditions of an operating reactor blanket. Furthermore, due to the opacity of liquid metals and limitations of measurement techniques (such as positioning and temporal and spatial resolution of the probes), it is impossible to obtain a complete picture of the flow. The only approach available to us is to conduct experiments at the parameters accessible in a laboratory and complement the experimental studies by numerical simulations designed to reveal the characteristics of flow and heat transfer and explain its physical mechanisms. This approach was adopted in some of the earlier studies, e.g., in Refs. 7 and 10. It is also followed in the present work.

We analyze the flow in a vertical rectangular duct under the influence of a horizontal magnetic field applied along the duct’s long side. The mean flow is directed downward. The constant-rate heating is applied to one or both walls parallel to the magnetic field. This configuration is present in many of the currently pursued designs of the liquid metal blanket modules.^11–15 Strong interest to it is also caused by the experimental results^5,16–18 recently obtained using the RK-2 facility. It has been found that at strong magnetic fields, the flows develop fluctuations of temperature with very high (several tens of degrees) amplitude. This has a significant influence on the overall heat transfer.⁶ We must also mention direct and potentially serious technological implications. As estimated in Ref. 19, the resulting cyclic thermal stresses can critically reduce the life-span of a heat exchange system, such as a blanket of a fusion reactor. In a real device, magnetic fields can be up to 10 T; therefore, we are focused on studying the strongest magnetic fields possible in experimental conditions.

It has been hypothesized in the experimental works that the fluctuations are caused by convection instabilities in which turbulence and, thus, turbulent mixing are suppressed by the magnetic field. The hypothesis has been confirmed and clarified in the numerical and theoretical studies.^10,20,21 In particular, good quantitative agreement between the experiments¹⁸ and simulations¹⁰ is found for the downward flow in a round pipe. The two-dimensional simulations conducted at very strong magnetic fields²¹ are also consistent with the experimental data. The key effect appears to be the unstable stratification, which develops in a duct with downward mean flow and constant-rate heating applied along its length. If the duct is infinite, the buoyancy force creates the so-called elevator modes, the exact solutions of the Boussinesq equations in the form of pairs of exponentially growing ascending/descending jets.²⁰ Such solutions exist in conventional hydrodynamic systems but are rarely realized due to their inevitable secondary instabilities and break down into turbulence. The three-dimensional instabilities and turbulent fluctuations are suppressed in systems with strong magnetic fields, which allows the jets to grow to very strong amplitudes. This also happens in ducts of finite size, where the classical elevator modes do not develop, but in a significant range of parameters, we observe the growth of pairs of strong and long ascending/descending jets and the formation of inflection points in the profile of streamwise velocity. Such flows are unstable to quasi-two-dimensional (often abbreviated to Q2D—the term for structures, which are nearly two-dimensional along the magnetic field lines outside of the boundary layers^4,9,21) rolls of Kelvin–Helmholtz type responsible for the high-amplitude fluctuations of temperature registered in the experiments.^10,21

It must be noted that the formation of inflection points and the instability to Q2D perturbation have also been found in vertical ducts with upward^4,22 or zero²³ mean flow, although in these cases the amplitude of the anomalous temperature fluctuations and the parameter range of their existence are substantially smaller than in the downward flow case.

The previous works have been predominantly focused on the configuration with a strongly asymmetric heating load. The heating has been applied to one wall or, in numerical models, internally, with the volumetric rate decreasing exponentially with the distance to one wall. At the same time, experiments, such as in Refs. 5 and 16, indicate that high-amplitude temperature fluctuations with distinct properties also exist when the heating is symmetric. Moreover, the structure of the applied heat flux has been found to strongly affect the pattern, amplitude, and typical frequencies of the measured temperature fluctuation signals.^5,16–18 The physical nature of this effect and, in general, the properties of the fluctuations caused by symmetric or mildly asymmetric heating have not been thoroughly studied and remain unclear.

The present study further explores the phenomenon of high-amplitude fluctuations in a downward flow with magnetoconvection with a particular emphasis on connection between the heating symmetries and the fluctuation properties. New experimental data obtained using the newly designed experimental facility RK-3²⁴ are presented. We also report, for the first time, the results of a high-resolution Direct Numerical Simulations (DNS) of this type of flow.

It has been found that depending on the parameters and the structure of the applied heat flux, the fluctuation signals have different forms, which implies different structures of the flow. At the same time, the connection between the heating configurations and the nature of the fluctuating flow states has not been thoroughly analyzed so far. The behavior, spatio-temporal properties, and range of existence of this phenomenon largely remain unclear.

This paper is organized as follows: the experimental and computational procedures are described in Sec. II; the results are presented in Sec. III; and finally, summarizing thoughts are outlined in Sec. IV.

Experimental studies are carried out using the RK-3 facility (HELMEF),²⁴ which consists of an electromagnet and several mercury loops of closed type. The electromagnet DEM-1 with cooled copper windings creates a magnetic field with a nearly uniform main (transverse) component in the gap of length about 1 m and width 80 mm between the magnet poles (see Fig. 1). The maximum induction of the magnetic field is 1.8 T.

A sealed centrifugal pump is used to drive the flow along the loop. Mercury is moved from the storage tank into the calming chamber, where pressure fluctuations from the pump are suppressed by the gas cushion located above the liquid surface. Next, the liquid flows through test section positioned in the gap of the magnet. After exiting from the test section, the liquid passes through two pipe-in-pipe heat exchangers, which are necessary to re-establish the working temperature of the fluid after it has been heated in the test section, and an electromagnetic flowmeter. Then, the liquid flows through the flowmeter valve block. The amplitude of pressure pulsations at the entrance into the test section is monitored by using a pressure sensor. Platinum resistance thermometers are utilized to monitor the liquid’s temperature at the entrance into the test section and the exit from the mixing chamber.

A schematic representation of the test section is shown in Fig. 1. Its main components are the input chamber, the honeycomb, the main part with plate heaters and sensors, the exit chamber, and the immersion measuring probe. The main part has a rectangular cross section of the inner dimensions 16 × 56 mm² with stainless steel (AISI 316) walls of 2 mm thickness (see Fig. 2).

The technique of scanning temperature measurements described in detail in Ref. 25 is used to study the flow. The probe is attached to the coordinate and moving mechanism and inserted through the exit part of the test section (see Fig. 1). It allows us to measure temperature at an arbitrary point within a given cross section. The sensor is a K-type thermocouple with a 0.25 mm hot junction, which allows us to measure temperature with 0.2 °C absolute accuracy and 0.01 °C sensitivity. The data acquisition system used with the thermocouple has a frequency of 100 Hz.

For this work, measurements are performed in the cross section indicated by the dashed-dotted line in Fig. 1 on a grid covering the entire cross section, as described in Sec. III A.

The duct is heated by two independently controlled flat metal heaters (see Fig. 2). When assembling the test section, a special heat-conducting paste is applied between the duct wall and the heater to ensure low contact resistance. Fiberglass wrapping is applied on the outer side of the heated section of the duct. In addition, two semi-cylindrical fluoroplastic clamps are laid on top of each heater to ensure better clamping to the walls. Heat flux sensors and thermocouples are mounted on the inner surfaces of the heating plates to control the heat loss and temperature of the heaters.

Due to the specifics of the manufacturing technology, there is a longitudinal weld on the inner surface on one of the wide sides of the duct. In order to prevent extraneous disturbances in the flow, the weld was lapped to a uniformly flat surface during the manufacturing of the duct.

The test section includes a honeycomb, which can be inserted inside the upstream portion of the test section, as illustrated in Fig. 1. The honeycomb consists of 105 stainless steel tubes with wall thickness 0.5 mm, inner diameter 2 mm, and length 80 mm. The purpose of the honeycomb is to analyze the effect of inlet conditions. It is known (see, e.g., the classical experiment^26,27 and the recent computational study²⁸) that in isothermal flows, the evolution of MHD turbulence along a duct is strongly affected by the state of the flow when it first experiences the significant impact of the magnetic field. In particular, if turbulence is created by a honeycomb upstream of the zone of the strong magnetic field and has time to develop before such an impact occurs (the configuration A in Fig. 1), monotonic suppression of velocity fluctuations along the duct is observed. If, however, the honeycomb is located within the strong magnetic field zone (the configuration B in Fig. 1), low-frequency velocity fluctuations of very high amplitude are registered downstream. As visualized in DNS,²⁸ the fluctuations are a manifestation of the transformation of the honeycomb-created jets and shear layers into quasi-two dimensional flow structures (vortex sheets and then large-scale vortices) caused by the magnetic field. Introducing a movable honeycomb into our experiment allows us to see whether an analogous effect of the initial state is present in flows with thermal convection.

Our last comment concerns the necessary steps of installation of the test section. After installing the measuring probe with a coordinate and movement mechanisms in the test section, the coordinates of all four duct walls are determined by endoscope inspection. Next, the honeycomb and the inlet section are installed. At the end, the assembled test section is mounted between the poles of the electromagnet and connected to the main RK-3 (HELMEF) loop.

The research conducted at the RK-3 facility is a continuation of the experimental program previously conducted at the RK-2 facility.⁶ Table I compares the main parameters of the two experiments. The values of the physical properties of mercury, such as kinematic viscosity ν, electric and thermal conductivities σ and λ, thermal diffusivity χ, and the coefficient of thermal expansion β are taken from Ref. 29 at mean temperature between inlet and outlet (about 30 °C). The other physical parameters used in this table are the mean flow velocity U and the thickness d_$w$ and electrical conductivity σ_$w$ of the wall.

We note that the average heat flux in the form q = (q₁ + q₂)/2 is used in the definitions of the Grashof number, Nusselt number, and other heating-related dimensionless parameters. The shorter half-width of the duct a is used as the typical length scale. Each of these choices is, evidently, not unique. Our selection is made to assure consistency with earlier experimental works conducted for various ducts and heating configurations.^5,6,17–19

We see in Table I that the main difference between the RK-3 and RK-2 facilities is in the strength and configuration of the imposed magnetic field. The magnet of the new facility allows us to reach almost twice larger values of the Hartmann number. It also produces a significantly longer segment of a nearly uniform magnetic field.

The liquid metal flow in the test section of the experiment is modeled as a flow of an incompressible electrically conducting Newtonian fluid in a duct of rectangular cross section (see Fig. 2). The imposed magnetic field is steady, non-uniform, and has two components (B_x, B_y), approximated using the model of Ref. 30, which provides a curl- and divergence-free magnetic field. The main component B_y has the same profile as in the experimental setup (see Fig. 1). The second component B_x is much weaker, especially inside the homogeneous region; however, it provides a proper closure for the realistic entry.

For practically all laboratory and technological flows of liquid metals, including the flows in the experiments presented in this paper and in blankets of fusion reactors, the magnetic Reynolds and Prandtl numbers are small: Re_m = UD/η ≪ 1 and Pr_m = ν/η ≪ 1. Here, σ and μ₀ are the electric conductivity of the fluid and magnetic permittivity of vacuum. At small Re_m and Pr_m, the analysis can be carried out in the framework of the quasi-static approximation of magnetohydrodynamic interactions in which the additional magnetic field induced by the electric currents flowing in the fluid is assumed to be small in comparison to the imposed field and, therefore, neglected in the expressions for Ohm’s law and Lorentz force (see, e.g., Ref. 31).

The non-dimensional governing equations are

The buoyancy and Lorentz forces are

and the electric current is given by Ohm’s law with the electric potential determined from the conservation of charge,

Here, u, p, T, j, and ϕ are the fields of velocity, pressure, temperature, electric currents, and electric potential; e_g is the unit vector in the direction of gravity; and b is the non-dimensionalized imposed magnetic field. As the typical scales, we use the mean velocity U for velocity, the shorter half-width D of the duct for length, D/U for time, ρU² for pressure, qD/λ for temperature, the maximum strength B₀ of the transverse component B_y for the magnetic field, and UB₀D for the electric potential. Further information, in particular the definition of q, is provided in Table I.

The non-dimensional parameters of the problem are Re, Pr, Gr, and Ha defined in Table I. Another commonly used MHD parameter is the Stuart number N = Ha²Re⁻¹, which estimates the ratio between the Lorentz and inertia forces.

The boundary conditions at the walls are no-slip for velocity and perfect electric insulation for the electric potential,

The approximation of perfectly insulating walls is justified by very small values of the wall-to-fluid conductance ratio in the experiment, as shown in Table I. The Neumann boundary conditions are applied for temperature, $∂T∂n|y=±b=0$ at the adiabatic side walls (shorter side), and $∂T∂n|z=−a=−q2*$ and $∂T∂n|z=a=−q1*$ at the heated walls (longer side). Here, $q1*$ and $q2*$ are the non-dimensionalized heat fluxes.

At the inlet section, we prescribe a pattern of streamwise velocity imitating the pattern of round jets exiting the honeycomb (see Fig. 3). A similar procedure was used in our earlier studies.^28,32 The number of jets and their diameter correspond to those of the experimental setup. Since, in the experiment, the flow in a single tube has Re ≈ 1200, the jet is represented by a laminar parabolic velocity profile.

The temperature distribution at the inlet is assumed to be uniform and unperturbed (which is a reasonably good approximation of the experimental conditions) and is set as T = 0. In addition, to better reproduce the process of mixing of the honeycomb jets and make it less grid-dependent, random perturbations of amplitude 10⁻³ are added to the inlet velocity at every time step.

The commonly used choices of the velocity and temperature conditions at the exit are ∂f/∂x = 0 and ∂f/∂t + U∂f/∂x = 0 (convective condition). Here, U is the mean velocity, x is the streamwise coordinate normal to the exit cross section, and f stands for a velocity component or temperature. In our earlier studies of similar geometries,^28,32 the condition of zero streamwise gradient was applied. In the present simulations of MHD flows with strong convection effects, we have found that this condition leads to a strong reduction of the time step Δt needed to avoid numerical instability when high-amplitude perturbations appear in the exit section. The convective outflow condition has been found to be far less sensitive to perturbations and, thus, has been used for both velocity and temperature in all the simulations reported in this paper.

The Neumann boundary condition ∂ϕ/∂x = 0 for the electric potential is applied at the inlet and exit of the duct. This implies that electric current j can exit the domain and, thus, is assumed to form virtually closed loops outside.

The in-house DNS code developed for incompressible HD and MHD flows in rectangular geometries is used for this study. The solver has been successfully applied in a number of prior studies for simulations of turbulent and transitional MHD and convection flows at high Re and Ha (see, e.g., Refs. 8, 9, 28, and 32–34). A detailed description of the method is provided in Ref. 35, with further details concerning its extension to flows in domains with inlet–exit conditions explained in Ref. 34. Only the main features of the discretization procedures and the new features implemented in the course of this work are discussed below.

The solver is based on a finite difference scheme. The time discretization of the momentum equation (1) and the convective outflow condition uses the fully explicit Adams–Bashforth/backward-differentiation method of second order.³⁶ The temperature equation (3) is solved with a semi-implicit approach. The non-linear terms and the streamwise part of the conduction term are integrated explicitly using the same scheme as for the momentum equation. The rest of the Laplace operator is integrated implicitly. This helps us to maintain a reasonably large time step Δt under the conditions of low Pr and small wall-normal grid sizes while retaining the flexibility of treating mixed types of inlet/exit conditions. The incompressibility condition is satisfied by applying the standard projection method,³⁷ which requires a solution of an elliptic problem for pressure at every time step.

For the spatial discretization we use a highly conservative scheme proposed for incompressible flows³⁸ and later extended to the case of low-Re_m MHD.^39,40 The scheme is of the second order of approximation. The discretization is conducted directly on a collocated non-uniform rectangular grid discussed below. The solution variables $v$⁠, p, j, T, and ϕ are all stored, and the governing equations (1)–(5) are all approximated at the same grid points. Correct coupling between the vector and scalar fields (velocity–pressure and current–potential) is achieved via the use of the flux variables for velocity and electric currents obtained by interpolation to the staggered half-integer points. Solvers for the 3D elliptic problems for pressure and electric potential use a combination of the fast cosine transform in the x-direction and the 2D cyclic reduction method⁴¹ in the (y, z)-plane. The cosine expansion matches the boundary conditions on ∂f/∂x imposed on pressure and electric potential at the inlet and exit. For temperature, we solve a 2D elliptic problem in each (y, z)-plane, except for the inlet and exit sections.

The solver is parallelized using the hybrid MPI-OpenMP approach.

The simulations have been conducted in the domain with dimensions L_x × L_y × L_z = 120 × 7 × 2, which reproduces the experimental test section from the exit of the honeycomb to a position sufficiently far downstream of the measurement cross section (see Fig. 1). The computational grid is orthogonal, uniform in the x-direction and clustered toward the duct’s walls in the (y, z)-plane by a coordinate transformation. The grid must be strongly clustered in flows with strong transverse magnetic fields in order to resolve the Hartmann and sidewall boundary layers. At the same time, accurate reproduction of mixing, shear instabilities, and development of turbulence in the inlet zone adjacent to the honeycomb requires fine resolution in the bulk region.

We have conducted a series of preliminary simulations aimed to check different types of grid stretching (e.g., hyperbolic tangent and Gauss–Lobatto) and their impact on the mixing of honeycomb jets and formation of turbulence. It has been found that the best results, both in terms of the shortest length of mixing and uniformity of the turbulent velocity profile, are yielded by the modified Gauss–Lobatto grid clustering,

where −1 ⩽ η ⩽ 1 and −1 ⩽ ζ ⩽ 1 are the coordinates in which the grid is uniform. This type of grid clustering and the grid with N_x × N_y × N_z = 7200 × 512 × 128 points were used for all the simulations presented below. The choice of the grid parameters also relies on our prior study of MHD duct flow with a honeycomb,²⁸ where a computational grid with similar properties was used.

Results were obtained independently using the experimental and numerical approaches. Experiments were designed to recreate and further investigate the phenomena reported in Refs. 5 and 16. The simulations reproduced the conditions of the experiment and provided a detailed insight into the nature of the experimentally observed phenomena.

Temperature measurements were made in the duct’s cross section indicated by the dashed-dotted line in Fig. 1. The movable single-point probe described in Sec. II A of this paper was used. At every location, the temperature signal was recorded for not less than 30 s of the evolution of fully developed flow. The parameter range shown in Table I was covered.

Several flow regimes were studied in detail. Temperature was measured sequentially at points of a grid covering the entire duct’s cross section. We present the results of the experiments at Re = 5000 and Gr ≈ 1.25 × 10⁶. It must be noted that Gr can be maintained in the experiments with a certain accuracy estimated in the case of the results presented here as 1.15 × 10⁶ < Gr < 1.35 × 10⁶. This amounts to a reasonable uncertainty of ±8% from the mean value Gr = 1.25 × 10⁶. The latter will be used in the following as a reference value for both experiments and numerical simulations.

As examples of a detailed study, Figs. 4 and 5 present the results obtained at Re = 5000, Ha = 0, 200, and 375, Gr ≈ 1.25 × 10⁶, and with the honeycomb configuration B (see Fig. 1). The grid points are indicated by dots. The shown two-dimensional distributions of the mean value, standard deviation, and skewness of temperature are the results of the post-processing performed by the Kriging gridding method in the Surfer software package. Examples of the temperature signals measured at one location are given in Fig. 6.

In the flows with zero magnetic field [see Figs. 4(a)–4(c), 5(a)–5(c), and 6(a)], the measured data are consistent with the conventional picture of a turbulent flow. The temperature fluctuations have moderate amplitudes [see Figs. 4(b), 5(b), and 6(a)]. The skewness [see Figs. 4(c) and 5(c)] demonstrates moderately strong excursions of cold fluid toward the heated walls caused by turbulent fluctuations. An important feature of the mean temperature fields [see Figs. 4(a) and 5(a)] is the slight asymmetry of the cold central zone in the y-direction. It is especially pronounced in the flow with two-sided heating [see Fig. 4(a)] and can be attributed to the natural asymmetry of the flow field averaged on a relatively short (30 s) time period.

An imposed magnetic field completely changes the picture. Both systems respond to a moderately strong field (Ha = 200) by developing temperature fluctuations with the standard deviation and skewness substantially higher than in the flows at Ha = 0. The nature of the fluctuations is different in the two heating configurations. In the case of one-sided heating, the signal consists of long periods of temperature growth interrupted by abrupt drops [see Fig. 6(c)]. This is similar to the typical behavior found earlier in studies of systems with one-sided heating.^{5,6,10,16–18,21} The fluctuations in the duct with two-sided heating have noticeably smaller amplitude and a nearly harmonic pattern [see Fig. 6(b)].

At Ha = 375, the temperature fluctuations in the flow with two-sided heating are fully suppressed [see Figs. 4(h) and 4(i)]. An exception is a small area in the corner, where small-amplitude perturbations appear to remain. The suppressed mixing leads to stronger gradients of the mean temperature [see Fig. 4(g)]. In the case of one-sided heating, the temperature fluctuations exist at Ha = 375. Moreover, their amplitude is slightly higher than at Ha = 200 [see Fig. 6(d)].

In addition to the detailed analysis of flow regimes at Ha = 0, 200, and 375, we conducted a parametric study of the effect of Ha. Measurements in the bulk (y ≈ 0, z ≈ 0) and at two locations closer to the sidewall layers (y ≈ 0, z ≈ ±0.6) were used. Temperature fluctuations are characterized by the non-dimensional standard deviation $σT*$ in Fig. 7. The data for Re = 5000 with one-sided and two-sided heating as well as two honeycomb configurations A and B are shown. For comparison, we also present data of the previous experiments⁶ performed at the RK-2 facility. It must be noted that the experiments⁶ were carried out at the same flow and heating rates as the new experiments, but the non-dimensional parameters were slightly different due to the higher average temperature of mercury in the test section. In particular, the Grashof number was evaluated as Gr = (1.6 ± 0.1) × 10⁶.

The parametric study confirms the principal difference between the flows with one-sided and two-sided heating. One aspect is the significant (3–4 times) difference in the amplitude of the temperature fluctuations. Comparison of the spatial distributions of the standard deviation in Figs. 4(e) and 5(e) shows that the amplitude difference is not specific to a point of measurement, but a general property of the entire flow domain.

We also see that in the case of two-sided heating, the high-amplitude fluctuations exist only in a limited range of Ha, approximately between Ha = 100 and Ha = 300. No such high-Ha cutoff is found when heating is one-sided.

The data in Figs. 6 and 7 show that the location of the honeycomb does not have a significant impact on temperature fluctuations. The presence (configuration A) or absence (configuration B) of small-scale turbulence at the entrance into the magnetic field appears not to influence the convection instability in a noticeable way. The only exception, for which we do not have an explanation at the moment, is the slight deviation between the data for configurations A and B observed in the duct with one-sided heating at Ha = 100.

Figure 7 shows good agreement between the new experiments and the experiments conducted on the RK-2 facility.⁶ The fact that very close amplitudes of fluctuations are found at close but still different (≈1.25 × 10⁶ vs ≈1.6 × 10⁶) values of Gr on two facilities with no common parts is viewed by us as (i) a validation of the experimental model and (ii) confirmation that development of high-amplitude temperature fluctuations in strong magnetic fields is a robust physical phenomenon not limited to a narrow parameter range and not related to the features of a particular experimental setup.

One significant difference between the data obtained at RK-2 and RK-3 is, however, clearly visible in Fig. 7 and has to be explained. The fluctuation amplitude in the duct with one-sided heating was observed to sharply decrease when Ha approached the highest value Ha = 200 achievable at the RK-2.⁶ It was hypothesized on the basis of those experiments that the fluctuations were limited to a range of moderate values of Ha and would not exist at the order-of-magnitude larger values typical for a fusion reactor blanket. We see in Fig. 7 that no such decrease is observed at the new facility. The fluctuation amplitude continues to increase to the highest tested value Ha = 375. We conclude that the hypothesis is not confirmed. The drop of the amplitude at the RK-2 facility was probably caused by magnet features, most likely a shorter zone of the uniform magnetic field, which meant the reduced length of MHD flow development.

The experimental data^{5,6,16–18,20} suggest that the temperature fluctuations detected at strong magnetic fields are caused by large-scale instabilities in a flow modified by natural convection effects in the conditions of suppressed turbulence. The structure and dynamics of the instability modes cannot be fully revealed in experiments. Besides that the experimental findings end up with other unanswered questions, such as where and how the instabilities appear, whether they spread over the entire domain or remain localized, and what is the nature of the two clearly distinct families of unstable solutions obtained in the experiments with one-sided and two-sided heating.

The results of the high-resolution numerical simulations directed toward answering these question are presented below. It must be noted that our work is the first attempt of three-dimensional DNS analysis of this type of flow. The previous DNS studies addressed flows in round pipes¹⁰ or two-dimensional flows in the framework of an asymptotic high-Ha approximation.²¹ 

Due to the substantial computational cost of simulations, the analysis is focused on the effect of Gr and the type of wall heating. Fixed values Re = 5000 and Ha = 200 are used. Furthermore, since the experiments do not reveal any significant differences between the flows with two honeycomb configurations, only one (configuration A) is chosen for analysis. To facilitate a direct comparison with the experiment on the RK-3 facility, we present below the numerical results for the reference value Gr = 1.25 × 10⁶ (see Figs. 4–7).

As typical examples, structures of fully developed flows with one-sided and two-sided heating at Gr = 1.25 × 10⁶ are shown in Figs. 8–10. We see in Fig. 8 that turbulent fluctuations are suppressed in the zone of the magnetic field. Further downstream, velocity fluctuations of distinctly Q2D form develop. The discussion of their structure will be based on the two-dimensional illustrations of flow fields in the plane y = 0 (see Figs. 9 and 10).

Figure 9 shows the instantaneous distributions of velocity components u_x and u_z and temperature T. To better illustrate the correlation between the fluctuations of u_z and T, their plots are combined in the left panels of Figs. 9(a) and 9(b). A closer look at the downstream part of the duct is provided by Fig. 10, which combines distributions of the velocity component u_z with vector fields of velocity perturbations u′ with respect to the time-averaged mean.

We see that in both flow configurations, the profile of u_x is deformed by the buoyancy force, with velocity strongly reduced and even reversed near the heated walls. Soon after the deformation becomes pronounced, approximately at x = 40 in the case of one-sided heating and at x = 100 in the case of two-sided heating, the flow develops instability in the form of large-scale vortices. As argued in earlier works^{4,10,20,21,23,42} and discussed later in this paper, the most plausible explanation of the instability is the development of inflection points of the streamwise velocity profile.

As the vortices move along the duct, they cause the fluctuations of temperature registered in the experiment. To validate this statement, the computed single-point temperature signal is compared with the experimental data in Fig. 11. We see a reasonably good agreement for both heating configurations. Qualitatively, the shapes of the experimental and computational curves are similar. The quantitative parameters are also close. The typical fluctuation amplitudes are about 20 °C in the case of one-sided heating and about 5 °C in the case of two-sided heating in both DNS and experiment. The dominant frequencies are also in agreement.

Further comparison between DNS and experiment is presented in Figs. 12 and 13. The DNS fields are used to compute the distributions of the same statistical properties (mean value, standard deviation, and skewness) of temperature at the same cross section x = 76 as in the experiment (see Figs. 4 and 5). Time averaging over 320 non-dimensional units (about 40 s in the dimensional units of the experiment) is used in the DNS, which is similar to the averaging over not less than 30 s in the experiment. The temperature signals in Figs. 6 and 11 indicate that such averaging periods are sufficient in the case of two-sided heating, but only marginally so when the heating is one-sided.

Figures 12 and 13 reveal certain differences between the DNS and experimental data. All of them are attributed to the imperfections of the experimental procedure and numerical model. One already mentioned is the shortness of the averaging time in the case of one-sided heating, which leads to statistical uncertainty, especially in the calculation of skewness. We should also mention that in DNS, the averaging is done simultaneously for all spatial points, while consecutive time periods, one for each spatial point, are used in the experiment. Another evident imperfection is the relative roughness of the grid of measurement points in the experiment.

Other less obvious imperfections are related to the effects of walls. The DNS model uses the idealization in which a wall is either perfectly insulated or heated at a perfectly uniform and steady rate. In the experiment, the walls have finite thermal conductivity and impose some thermal damping on temperature fluctuations. Furthermore, it was discovered after completion of the work presented here that the imposed wall heat flux in the experiment was slightly higher near the corners than in the middle of the walls. Finally, the accuracy of temperature measurements is known to be negatively affected in the vicinity of the walls.

Considering all these imperfections, the experimental and DNS pictures of the flow appear to be in a reasonably good qualitative and quantitative agreement with each other. This demonstrates the robustness of the physical mechanisms of convection instabilities explored here.

One feature of the computed flow in the duct with two-sided heating is the low amplitude of temperature fluctuations in the stripe around the centerline z = 0 [see Fig. 12(d)]. The respective flow fields in Fig. 9 show that the fluctuations of the temperature field, indeed, do not penetrate this stripe at x = 76 (this, however, happens further downstream). Interestingly, the experiment shows strong fluctuations in this zone at the same flow parameters [see Fig. 12(c)]. As we discuss in detail below, this can be attributed to the imperfect symmetry of the experimental flow.

The differences caused by different heating configurations are clearly visible in Figs. 8–13. In the case of two-sided heating, the deformation of the streamwise velocity profile and the pattern of vortices are symmetric with respect to the central plane z = 0. No such symmetry exist in the case of one-sided heating. The profile is deformed, and the vortices first appear near the heated wall. The vortices move into the central part of the duct and toward the cold wall as they are transported downstream by the mean flow.

The fluctuation properties are noticeably different between the two heating configurations. One obvious aspect is that the amplitude of temperature fluctuations at the same downstream location is about four times higher in the case of one-sided heating. Another aspect is the pattern of Q2D vortices and corresponding fluctuations of the temperature field. It is regular in the case of two-sided heating (see Figs. 9–11). On the contrary, an irregular pattern is found in the duct with one-sided heating. Furthermore, in that case, the vortices form clusters separated by zones of unperturbed streamwise flow [visible approximately between x = 80 and 100 in the snapshot in Fig. 9(a)]. These features are consistent with the findings of the earlier studies of similar systems.^10,21

Additional simulations were performed to analyze the effect of Gr and flow symmetry on the instability. The results are illustrated in Figs. 14–18. The fully developed flow states were computed in a wide range of Gr. A flow obtained in an earlier simulation or a laminar isothermal flow was used as an initial condition. It should be mentioned that the choice of the initial state affects the first phase of the flow evolution, but not its fully developed form. In particular, no hysteresis is found in two series of simulations with gradually increasing or decreasing Gr and the flow field from one simulation used to start the next.

Flows with both types of heating are stable at small Gr [see Fig. 15(a) as an example]. After the initial development, such a flow attains a state in which all the velocity and temperature fluctuations are fully suppressed in the entire zone of the strong magnetic field (approximately at x > 30). The stability of such states was verified in the simulations by computing them during more than 1000 non-dimensional time units.

As an example of transition to the instability, Fig. 14 shows temperature signals for the parametric study of both one-sided and two-sided heating configurations. For the one-sided heating case [Fig. 14(a)], we see that the flow is stable at Gr = 0.3 × 10⁶ and unstable at Gr = 0.4 × 10⁶. Further increase in Gr leads to (i) earlier appearance of temperature fluctuations after the initial phase, (ii) higher fluctuation amplitude, and (iii) shorter typical time periods between the abrupt temperature drops. All these features can be attributed to the faster growth of the ascending/descending jets caused by strong buoyancy forces. Interestingly, the dependency of the fully developed state on the value of Gr weakens at Gr above, approximately, 1.0 × 10⁶. As an example, the signal for Gr = 1.0 × 10⁶ in Fig. 14 is very similar to the signal at Gr = 1.25 × 10⁶ (see Fig. 11) and higher Gr (not shown) in terms of the fluctuation amplitude and dominant frequencies.

The similarity with the flows at Gr = 1.25 × 10⁶ and higher Gr is also found for the spatial structure of the velocity and temperature fields (not shown). Furthermore, the instability shows convincing similarity with the instabilities in systems with one sided-heating, but different geometries^10,21 and much higher Ha.²¹ We conclude that further increase in Gr within the computationally feasible range is unlikely to significantly change the flow’s behavior and leave the matter for future studies.

In the case of two-sided heating [see Fig. 14(b)], the first traces of instability with an extremely weak amplitude appear at Gr = 1.2 × 10⁶, albeit the flow remains generally stable. At a higher Gr = 1.3 × 10⁶, one can already see the onset of a fully developed unstable regime in the form of an almost ideal harmonic time-signal. With the further growth of Gr, the oscillations exhibit deviation from the harmonic form, visibly evolving additional frequencies, as shown in Fig. 14(b) for Gr = 1.5 × 10⁶. This can be attributed to the onset and further development of longitudinal oscillations between the neighboring pairs of Q2D vortices. In fact, the smallest value of Gr at which the instability is found in simulations is Gr = 1.25 × 10⁶ (see Figs. 8–10). We note that this result is not fully consistent with the experiment, where the temperature signal corresponding to instability was detected at Gr = 10⁶ [see Fig. 6(b)].

Further increase in Gr leads to larger fluctuation amplitudes [see Figs. 16(b) and 16(c)], although they remain significantly (about three times) smaller than in the case of one-sided heating. This result is consistent with the experimental findings. Comparing the plots for Gr = 2.5 × 10⁶ [Fig. 15(b)] and Gr = 1.25 × 10⁶ [Figs. 9(b) and 10(b)], we see that the flow structure does not change in a principal way with the growth of Gr. The instability creates the same regular pattern of large vortices occupying the central part of the duct. The only difference is that in the flow at Gr = 2.5 × 10⁶, the vortices appear much higher upstream and evolve to full amplitude by the time they reach the cross section x = 76, where temperature is measured in the experiment.

Additional simulations have been performed to reveal how the flow’s behavior in the duct with two-sided heating is affected by small deviations from perfect symmetry. The results give another perspective for the asymmetry effect. Furthermore, they allow us to test the possible reasons for the inconsistency between the simulations and experiment in the case of two-sided heating, in particular, for the fact that the computed flow is stable at Gr = 10⁶, while the experiment shows instability [cf. Figs. 15(a) and 6(b)]. Finally, the explored asymmetry effects are directly relevant to the true orientation and heat fluxes in the ducts of a fusion reactor blanket.

Two possibilities, which are relatively easy to mimic in the simulations, are considered: a slight (10°) deviation from the vertical orientation and a slight asymmetry of the wall heating (⁠$q1*=1.1$ and $q2*=0.9$ instead of $q1*=q2*=1$⁠). In both cases, the deviations from the ideal symmetry are by ∼10%. Flows at Gr = 10⁶, Re = 5000, and Ha = 200 are computed. A 10° deflection from the vertical orientation does not produce any significant effect. The computed flow remains stable and similar to that in Fig. 15(a). An instability is, however, triggered by the asymmetric wall heating. The pattern of the Q2D vortices and the duct area occupied by them are very similar to those in the case of perfectly symmetric heating at Gr = 1.25 × 10⁶ [see Fig. 9(b)], yet the arrangement of the vortices is somewhat less regular. A similar effect is shown by the temperature fluctuation signal in Fig. 16(d). Because of irregularities, the signal does not seem truly harmonic and exhibits closer similarity to the experimental signal shown in Fig. 6(b).

As a further illustration of the properties of unstable flow regimes, Figs. 17 and 18 show the time signals of the streamwise and transverse velocity components u_x and u_z. The component u_y in the direction of the magnetic field is very weak in Q2D flows (with amplitude about two orders of magnitude weaker than of u_z in our simulations) and, therefore, not shown. We see in Figs. 17 and 18 that the behavior of velocity perturbations is similar to that of temperature perturbations. At the same parameters, the flow with one-sided heating shows much higher fluctuation amplitude than the flow with two-sided heating. The signals are regular and nearly harmonic at two-sided heating [see Figs. 17(b), 17(c), 18(b), and 18(c)]. In the case of one-sided heating [see Figs. 17(a) and 18(a)], we see irregular signals with bursts of u_z accompanying abrupt changes of u_x. A 10% deviation from symmetry in the two-sided heating case results in the fluctuation signals reminding more the signal at one-sided than at two-sided heating [see Figs. 17(d) and 18(d)].

We start by summarizing the main experimental and numerical results presented so far. Flows with both one-sided and two-sided heating show instabilities leading to the development of strong Q2D vortices in the downstream portion of the duct. The transport of the vortices by the mean flow is manifested by high-amplitude fluctuations of temperature and components u_x and u_z of velocity. The cases of one-sided and two-sided heating are principally different in several aspects. At the same values of Re and Ha, a duct with one-sided heating becomes unstable at lower Gr comparing to two-sided heating. At the same Re, Ha, and Gr, the unstable flows have different fluctuation amplitudes: 3–4 times higher at one-sided than at two-sided heating. The time signals registered in unstable flows are of different types: irregular signals with bursts of u_z and abrupt changes of u_x and T in the case of one-sided heating and regular, nearly harmonic, smooth signals in the case of two-sided heating. We have also found that even a small deviation from the perfect symmetry of a two-sided heating leads to a behavior with clear features of the one-sided heating case.

Before presenting the convection instability mechanism, we have to exclude the possibility that the observed phenomena are caused by other reasons, specifically by the hydrodynamic instability of an MHD duct flow or by the effects of entry into the magnetic field.

The first of these possibilities is easily excluded by looking at the non-dimensional parameters. It is known^33,43–45 that loss of stability of a laminar base state and transition to turbulence in an MHD flow first occur in the range of the Hartmann-thickness based Reynolds number R ≡ Re/Ha between approximately 200 and 500, depending on the strength and shape of introduced perturbations. In this work, we observe instability at much smaller R, e.g., at R = 25 when Re = 5000 and Ha = 200 or at R = 13.33 when Re = 5000 and Ha = 375.

The second possibility refers to the distortion of the streamwise velocity profile that occurs when a flow enters a zone of the strong transverse magnetic field. The shear layers developing near the duct’s walls parallel to the magnetic field may give rise to Q2D vortices as shown, e.g., in Ref. 28. These vortices, however, rapidly decay due to viscous and Joule dissipations, and gradual weakening of the shear layers, so they inevitably disappear a short distance downstream. In our simulations, vortices first appear far downstream of the entrance into the magnetic field and show no tendency to decay.

We conclude that the instability phenomena presented in this paper cannot be explained by isothermal MHD effects. It is virtually beyond doubt that they are caused by the thermal convection mechanism.

The mechanism is schematically illustrated in Fig. 19. It was described in earlier studies for flows with one-sided wall heating or strongly asymmetric internal heating.^4,5,10,18,21 Its modification for the case of perfectly or nearly symmetric heating is, for the first time, presented below.

A heat flux applied to walls or internally in a downward duct flow causes buoyancy forces, which change the profile of streamwise velocity. In an infinite duct, strong buoyancy leads to exponentially growing elevator mode solutions stabilized and, thus, realizable in a strong magnetic field.²⁰ Such solutions are impossible in realistically finite ducts, but the profile is still modified significantly, with strongly reduced streamwise velocity or even reversed flow in the zones of intense heating. This is illustrated schematically in Fig. 19 and, using our computed flow fields, in Fig. 20. As suggested earlier,^42,46,47 the modified Q2D flows are prone to instabilities of two types: the Kelvin–Helmholtz instability associated with the inflection points that appear in the streamwise velocity profile or the boundary layer instability. Our study, in agreement with earlier works,^4,5,10,18,21 indicates that in the downward flow with convection and in the parametric range considered here, only the Kelvin–Helmholtz instability mechanism is active. Q2D vortices are produced at inflection points and then grow in size and migrate in the transverse direction as they are transported downstream.

The proposed mechanism explains the differences between the flows with symmetric and asymmetric heating observed in our study. At the same Gr, the actual heating rate applied to a wall in the one-sided heating case is twice as strong as in the case of two-sided heating. The distortion of the streamwise velocity profile and, thus, the strength of the inflection point are substantially higher (see Fig. 20). The unstable growth is faster and leads to vortices, which are stronger and appear at smaller Gr.

In a subtler way, the proposed mechanism also explains why the perturbation signal is nearly harmonic in the case of two-sided heating, but irregular, with abrupt changes in the case of one-sided heating. As schematically illustrated in Fig. 19 and visible in plots of computed flow fields, most clearly in Fig. 10, large coherent vortices form principally different patterns in the two cases. In the case of one-sided heating, each vortex developing near the inflection point at a certain streamwise location creates transverse mixing thus disrupting the inflection point profile. The instability at this location is prevented for some time needed for the profile to be restored by the buoyancy force. The result of this dynamics is a pattern of large vortices, all with the same sense of rotation, separated from each other by some distances [see Fig. 19(a)]. Their nonlinear dynamics leads to the formation of clusters [see Fig. 10(a)] and irregular fluctuation signals [see Figs. 6, 11(a), 16(a), 17(a), and 18(a)]. Each abrupt change of the perturbation variables, for example, an abrupt drop of temperature in Figs. 6, 11(a), and 16(a) or a spike of transverse velocity in Fig. 18(a), corresponds to the time moment, when a well-mixed zone of a large vortex arrives at the measurement location.

The situation is principally different in the case of two-sided heating. Now, the streamwise velocity profile has two perfectly symmetric inflection points, each creating vortices of the same strength but opposite sense of rotation. In a manner similar to that in other flows with perfect reflection symmetry, e.g., in the von Karman flow past a cylinder, the vortices form a regular staggering pattern. Motion of the pattern in the streamwise direction generates the nearly harmonic signals we see in Figs. 16(b), 16(c), 17(b), 17(c), 18(b), and 18(c).

Interestingly, even a weak asymmetry of heating leads to the predominant growth of vortices at one, stronger inflection point and, thus, to the situation more closely reminding one-sided rather than two-sided heating case. We have found that the results with a 10% asymmetry of heating are more consistent with the experiment than the results obtained with a perfect symmetry. This clearly suggests that the ideal symmetry is inaccessible in a real experiment or a real application. Computational results obtained for a “symmetrically” heated duct should be viewed as a special, degenerative case.

The results of this work furnish strong support to the view that the appearance of high-amplitude low-frequency temperature fluctuations is a common and robust feature of magnetoconvection in downward flows in vertical ducts. The flow instabilities leading to the fluctuations have been reproduced with three independent case-studies: two experimental facilities and DNS. The studies are in a reasonably good (considering the imperfections of the experimental procedure and the idealized character of the numerical model) quantitative agreement with each other in terms of time-averaged temperature field, and amplitude and typical frequencies of the fluctuations. Different inlet conditions did not change the regions of existence or structure of the fluctuations.

The results of the high-resolution DNS confirm that in the considered range of Stuart numbers N ∼ 10, the unstable flow is Q2D. Simulations performed in the framework of the asymptotic 2D approximation, as it was done, e.g., in Ref. 4 or Ref. 21, should yield reasonably good results valid at least as preliminary estimations of the potential flow behavior. The advantage of such simulations would be that they can reach much higher values of Ha and Gr than in 3D simulations or laboratory experiments, possibly in the range anticipated for blankets of future fusion reactors.

The numerical and experimental studies presented in this paper highlight the main properties of large-scale flow fluctuations. In the case of one-sided heating, temperature and velocity at a given point show long periods of gradual variation interrupted, from time to time, by abrupt changes. The picture is completely different in the case of two-sided heating, where the temperature and velocity signals are nearly harmonic. This difference between the two systems is due to the different mechanisms of generation of secondary vortices: at a single or at two symmetric inflection points of the streamwise velocity profile.

The potential impact of the observed phenomena on the operation of a fusion reactor blanket is twofold. On the one hand, the observed regular Q2D structures can produce additional flow mixing and significantly enhance heat transfer, possibly up to the levels observed in turbulent flows. On the other hand, the system’s dynamics is very sensitive to its symmetry properties. The irregular bursts, observed in the cases of asymmetric heating even when the deflection from perfect symmetry is weak, can create unexpectedly high thermal loads or hot spots at the walls, which would pose potentially serious problems for the structural integrity of the blanket. Furthermore, our results show that the asymmetry caused by the duct’s orientation has a much weaker effect than asymmetry of heating. This may have direct implications for design of fusion blanket elements and deserves future study.

Finally, we note a rather good agreement between the experiment and DNS. Numerical simulations offer the possibility of a more detailed and complete analysis of the flow, often inaccessible in the experiment, as the information can be extracted not only locally, but in the entire domain. From this point of view, DNS, in particular, and numerical simulations, in general, must be seen as an essential tool of planning and extending the experimental program.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The authors are grateful for the invaluable help in setting up the work provided by recently deceased Professor V. G. Sviridov.

This work was supported by the RFBR NNIOa 18-508-12005 and DFG KR 4445/2-1 research projects within a joint Russian–German collaboration program. O.Z. was supported by the U.S. NSF (Grant No. CBET 1803730). Experiments are performed using the RF USF “Mercury MHD testing facility.” The simulations are performed at the Computer Center of TU Ilmenau and at the SuperMUC cluster at the Leibniz Rechenzentrum Garching (Large Scale Project Nos. pr62se and pn68ni of the Gauss Centre for Supercomputing).