Magnetoconvection in a horizontal duct flow at very high Hartmann and Grashof numbers

Direct numerical simulations and linear stability analysis are carried out to study mixed convection in a horizontal duct with constant-rate heating applied at the bottom and imposed transverse horizontal magnetic field. A two-dimensional approximation corresponding to the asymptotic limit of very strong magnetic field effect is validated and applied, together with full three-dimensional analysis, to investigate the flow's behaviour in the previously unexplored range of control parameters corresponding to typical conditions of a liquid metal blanket of a nuclear fusion reactor (Hartmann numbers up to $10^4$ and Grashof numbers up to $10^{10}$). It is found that the instability to quasi-two-dimensional rolls parallel to the magnetic field discovered at smaller Hartmann and Grashof numbers in earlier studies also occurs in this parameter range. Transport of the rolls by the mean flow leads to magnetoconvective temperature fluctuations of exceptionally high amplitudes. It is also demonstrated that quasi-two-dimensional structure of flows at very high Hartmann numbers does not guarantee accuracy of the classical two-dimensional approximation. The accuracy deteriorates at the highest Grashof numbers considered in the study.


Introduction
Combined convection and magnetohydrodynamic (MHD) effects dramatically change the nature of flows of electrically conducting fluids. The combination appears in many technological applications such as metallurgy, liquid metal batteries, and growth of semiconductor crystals (Ozoe 2005;Davidson 2016). Another prominent example is the liquid metal blankets of nuclear fusion reactors where an electrically conducting fluid (e.g., a PbLi alloy) serves as a coolant, radiation shield and tritium breeder (Abdou et al. 2015). A distinctive feature of this system is that the convection and magnetic filed effects are both exceptionally strong.
Many aspects of the transformation of flows of electrically conducting fluid under the influence of a strong magnetic field, such as suppression of turbulent fluctuations, anisotropic or quasi-two-dimensional (Q2D) states with zero or weak velocity gradients along the field lines, formation of MHD boundary layers, and delay of laminar-turbulent transition, are relatively well understood (see, e.g., Branover 1978;Davidson 2016;Sommeria & Moreau 1982;Zikanov et al. 2014). This paper addresses a recently discovered and still poorly understood phenomenon -the high amplitude fluctuations in flows in ducts and pipes (see, e.g., Genin et al. 2011;Vetcha et al. 2013;Zikanov et al. 2013;Belyaev et al. 2021). The term magneto-convective fluctuations (MCFs) proposed for the 2 R. Akhmedagaev, O. Zikanov, and Y. Listratov phenomenon by Belyaev et al. (2021) will be used in this paper. As discussed in detail in the review of Zikanov et al. (2021) and references therein, the fluctuations have been detected in experimental and computational studies of a large variety of systems: pipes and ducts of various orientations with respect to gravity, various heating arrangements, and various configurations of the magnetic field.
The fluctuations were called anomalous in some earlier works, e.g., by Zikanov et al. (2013) and Zhang & Zikanov (2014). This term now appears imprecise and somewhat misleading since it has been understood that the fluctuations are rather common. They occur in a wide variety of magnetoconvection flows. It must also be mentioned that, in a broader context, the magnetoconvective fluctuations are a part of the general phenomenon of large-amplitude fluctuations commonly found in flows, where turbulence is suppressed by a strong magnetic field and flow fields are strongly anisotropic or Q2D (see, e.g., Smolentsev 2021;Zikanov et al. 2021, for discussion and references).
The nature of the magnetoconvective fluctuations can be briefly described as follows. They appear in the conditions of a very strong magnetic field effect, i.e. in the range of Hartmann numbers, where turbulence is fully suppressed by magnetic damping. In experiments, the MCFs are manifested by oscillations of temperature with very high amplitude (up to 50K in some cases) and typical frequencies much lower than the frequencies of turbulence-induced fluctuations. Specific properties of the MCFs vary with the flow's configuration and values of the control parameters (Zikanov et al. 2021). The effect has potentially serious consequences for design and operation of liquid metal blankets of future fusion reactors. Should the fluctuations appear in an actual blanket, they may lead to strong and unsteady thermal stresses in the walls (see, e.g., Belyaev et al. 2018) possibly under the condition of significantly reduced strength of the wall material (Kolmakov et al. 2016). Due to their possibly very large amplitude, the stresses will threaten the structural integrity of a fusion reactor system. Significant effects on heat transfer, transport of tritium, and wall corrosion are also anticipated. As we discuss later in this section, it is yet impossible to say how realistic these expectations are, since no experiments or computations at very high Ha and Gr typical for reactor conditions have been conducted so far.
Flows in a rectangular duct with heating applied at the bottom and imposed transverse horizontal magnetic field (see figure 1) are considered in this paper. The configuration is not found in currently developed specific designs of liquid metal blankets of fusion reactors, although it may occur in future designs of an upper divertor and top blanket modules (Kirillov & Muraviev 1997). It is also important as an archetypal system, in which the MCFs were first identified (in Genin et al. (2011) and Zikanov et al. (2013), where they were named anomalous fluctuations) and explained.
Similar systems for either ducts or round pipes have been studied experimentally (Genin et al. 2011;Belyaev et al. 2015;Sahu et al. 2020) and numerically (Zikanov et al. 2013;Zhang & Zikanov 2014;Vo et al. 2017;Listratov et al. 2018). The flow is controlled by four dimensionless parameters: the Reynolds, Prandtl, Grashof and Hartmann numbers, Magnetoconvection in a horizontal duct flow Figure 1. Flow geometry and coordinate system. The arrows marked by letters g, B and q denote, respectively, the orientations of the gravity acceleration, magnetic field and wall heating.
Interesting results were obtained in the linear stability analysis of the Poiseulle-Rayleigh-Bérnard duct flow with a transverse magnetic field performed by Vo et al. (2017). Two-dimensional approximation valid in the limit of strong magnetic field presented later in this paper was used. One important result of Vo et al. (2017) is relevant to our work even though different boundary conditions were used. It was demonstrated that the convection instability occurs at moderate and high Grashof number (approximately above 10 6 ) at the Hartmann numbers (∼ 10 4 ) typical for reactor blanket conditions. The presence of MCFs in a horizontal round pipe with a lower half of the wall heated was detected in experiments (Genin et al. 2011;Belyaev et al. 2015) and explained in the linear stability analysis and direct numerical simulations (DNS) by Zikanov et al. (2013). Flows of mercury with Pr ≈ 0.022, Re up to 10 5 , Gr up to 10 8 , and Ha up to 500 were investigated. It was shown that at a strong magnetic field the suppression of flow structures having large gradients along the field lines resulted in the most unstable modes in the form of convection rolls with axes aligned with the field. The instability led to development of convection structures in the form of Q2D rolls. Transport of the rolls by the mean flow generated the MCFs.
The analysis was extended in the numerical simulations of Zhang & Zikanov (2014). Flows in a horizontal duct of aspect ratio Γ = 1 with bottom heating and transverse magnetic field at Pr = 0.0321, Re = 5000, 50 Ha 800 and 10 5 Gr 10 9 were investigated. The instability leading to the formation of Q2D rolls similar to those found in the pipe flow was detected at sufficiently high Gr and Ha.
Investigations of Zhang & Zikanov (2014) conducted in the broader range of parameters than for the pipe flow demonstrated existence of two distinct secondary flow regimes. The realization of the regimes depended on the relative strength of the convection and MHD effects. The low-Gr type characterized by Q2D distributions of velocity and temperature dominated by spanwise rolls appeared at Gr below a certain Gr * (Ha). At higher Gr , stronger convection resulted in three-dimensional (3D) flow states combining the spanwise rolls with streamwise ones (the geometrically preferred convection structure in pipes and ducts with bottom heating).
Flows of liquid metals in fusion reactor blankets and divertors are subject to very strong effects of convection (Gr ∼ 10 10 − 10 12 ) and magnetic fields (Ha ∼ 10 4 ) (see, e.g., Smolentsev et al. 2008Smolentsev et al. , 2010. Such extreme parameters present serious obstacles to analysis, because neither laboratory experiments nor 3D simulations of unsteady flow regimes in realistic blanket or divertor geometries can, at this moment, achieve such values. In an attempt to reach the typical blanket flow conditions, the data on two types of the secondary flow regime in a horizontal duct were extrapolated to high Gr and Ha by Zhang & Zikanov (2014). The extrapolation predicted existence of MCFs at the typical blanket parameters. It also predicted that the flow would likely be of the low-Gr type at Gr 10 10 and of the high-Gr type at higher Gr . The experiments in the pipe flow (see the review of recent results in Zikanov et al. 2021), on the contrary, indicate that MCFs may disappear at high Ha, so the extrapolation can be wrong. The nature of the convection flow at the parameters corresponding to ducts in blankets and divertors of an operating fusion reactor remains unknown, setting up the motivation for the present study.
The focus of our investigation is on the magnetoconvection in the range of very high Gr and Ha including the values typical for a reactor blanket and divertor. To the best of our knowledge, this study is the first to analyze the MCF effect in this range. Linear stability analysis and DNS of flows in a horizontal duct with Pr = 0.025, Re = 5000, 10 8 Gr 10 10 and 10 3 Ha 10 4 are performed. The study follows the work of Zhang & Zikanov (2014) but differs by much larger values of Gr and Ha and the aspect ratio Γ = 3.5 selected to match the new experimental facility (see, e.g., Belyaev et al. 2017), on which the same configuration is to be explored at Ha 10 3 and Gr 10 8 in the near future. Another essential difference between our work and the work by Zhang & Zikanov (2014) is that we carry out an in-depth analysis of the accuracy of the twodimensional approximation applied to Q2D flows at such high Ha.

Presentation of the problem
The flow of an incompressible, Newtonian, viscous, electrically conducting fluid (a liquid metal) with constant physical properties is considered. The fluid moves through a horizontal duct of aspect ratio Γ = 3.5 (see figure 1). Spatially uniform and timeindependent magnetic field B = Be y is imposed in the horizontal transverse direction. All walls are perfectly electrically insulated. The top and side walls are perfectly thermally insulated. The bottom wall is subject to uniform heating with the heat flux of constant rate q. The no-slip boundary conditions for velocity are applied at the walls.

Physical model
The Boussinesq and quasi-static approximations are applied. The quasi-static approximation is valid at small Reynolds and Prandtl numbers and usually utilized in numerical and theoretical studies of MHD flows of liquid metals (Davidson 2016). The approximation implies that the imposed magnetic field B is much stronger than the perturbations of the magnetic field b induced by the electric currents caused by the fluid motion. The induced magnetic field can be neglected in the expressions of the Lorentz force and Ohm's law. Furthermore, the induced field is assumed to adjust instantaneously to changes of velocity field.
The governing equations are rendered non-dimensional using the duct half-width in the magnetic field direction d as the length scale, mean streamwise velocity U as the velocity scale, wall heating-based group qd/κ as the temperature scale, B as the scale of the magnetic field strength, and dU B as the scale of electric potential. The equations can be written as: Magnetoconvection in a horizontal duct flow where u is the velocity field. The decompositions of the temperature and pressure fields commonly used in studies of mixed convection in ducts and pipes (see, e.g., Alboussière et al. 1993;Lyubimova et al. 2009;Zikanov et al. 2013;Zhang & Zikanov 2014;Zikanov et al. 2021) are applied. The decompositions are convenient, since they allow one to recast the problem in terms of the fluctuation fields, which are statistically uniform in the streamwise direction and, thus, study the flow in a relatively short segment of the channel with periodic inlet-exit conditions. The temperature field is written as a sum of fluctuations θ and the mean-mixed temperature where A = 2h/d is the cross-section area of the duct. One can also use the decomposition into fluctuations and simple mean temperatureT (x) = A −1 A T dA. Applying the energy balance between the wall heating and the streamwise convection heat transfer, we find that T m (x ) andT (x) are linear functions with the same derivative: where Π = 2 is the perimeter of the heated portion of the wall. The total pressure P is presented in (2.1) as where p(x, t) is the field of pressure fluctuations statistically homogeneous in the streamwise direction, andp is a linear function of x corresponding to the spatially uniform streamwise gradient dp/dx applied as a flow-driving mechanism. In the simulations discussed in this paper, the gradient is adjusted at every time step to maintain constant mean velocity. The second term of the decomposition becomes necessary in numerical models of mixed convection in non-vertical channels with periodic inlet-exit conditions. The component arises due to the buoyancy force caused by the mean-mixed temperature T m : where e z is the unit vector opposite to the direction of gravity (see figure 1). The force has a non-zero curl and, therefore, modifies the velocity field. Its action on the flow can be described by introducing the pressure fieldp, such that its vertical gradient balances F b,m . The pressure field is a two-dimensional function increasing with the streamwise coordinate x and vertical coordinate z. Its z-dependent x-gradient, which appears in the respective momentum equation, generates a flow in the positive x-direction in the lower part of the channel and in the negative x-direction in the upper part. The result is a topbottom asymmetry of the streamwise velocity profile and of the associated convection heat flux, which can dramatically change the structure of the flow at high Gr and Ha (see Zikanov et al. 2013;Zhang & Zikanov 2014Zikanov et al. 2021). The buoyancy force in (2.1) is 6 R. Akhmedagaev, O. Zikanov, and Y. Listratov (2.10) The Lorentz force is computed as (2.11) where e y is the unit vector along the imposed magnetic field (see figure 1). The electric current j is determined by the Ohm's law where the electric potential φ is a solution of the Poisson equation expressing the instantaneous electric neutrality of the fluid: (2.13) The inlet-exit conditions are those of periodicity of the velocity u, temperature fluctuations θ, pressure fluctuations p, and potential φ.

Two-dimensional approximation
Flows with a very strong imposed magnetic field are considered, so the Hartmann number and the Stuart number satisfy Ha 1 and N ≡ Ha 2 /Re 1, respectively. The flows are anticipated to have Q2D form with nearly zero gradients along the magnetic field lines except in the thin Hartmann layers at the walls perpendicular to the field. The 2D approximation proposed by Sommeria & Moreau (1982) can be applied in this asymptotic limit. The problem can be expressed in terms of the variables integrated wallto-wall along the direction of the magnetic field, leading to 2D dynamics for y-averaged quantities. The approximation has been verified and examined by Pothérat et al. (2000Pothérat et al. ( , 2005, and utilized in numerical studies of liquid metal flows in rectangular ducts (see, e.g., Pothérat 2007;Smolentsev et al. 2012;Vetcha et al. 2013;Vo et al. 2017;Zhang & Zikanov 2018).
The often applied abbreviation SM82 will be used for the model in the following. The y−independent solutions obtained in the framework of the model will be referred to as 2D solutions, while the full solutions obtained numerically without resorting to the model will be designated as 3D.
The SM82 model is derived for flows with Ha 1 and N 1, in domains with electrically insulating walls and constant wall-to-wall distance in the field direction. It utilizes the fact that the Lorentz force becomes nearly zero in the bulk region of Q2D flows in such geometries, and that the effect of the magnetic field on the flow is largely reduced to thin Hartmann layers and can be accurately modelled by the linear friction term in the momentum equation.
It must be noted that the original SM82 model was developed for isothermal flows. Its extension to flows with heat transfer and temperature variations was, to our best knowledge, first proposed by Smolentsev et al. (2008). As demonstrated in this and in the following studies (see, e.g., Gelfgat & Molokov 2011;Vetcha et al. 2013;Zhang & Zikanov 2018), the model can be extended to 2D approximation of temperature if the imposed heat flux is perpendicular to the magnetic field.
The SM82 version of (2.1)−(2.3) is where all the flow variables are now 2D fields obtained by wall-to-wall averaging. The term (Ha/Re)u represents the effect of friction in the Hartmann layers. The same notation as in (2.1)−(2.3) is used. The boundary conditions on velocity and temperature on the remaining two wall are the same as in the 3D model.

Numerical method
The governing equations (2.1)−(2.3) and (2.14)−(2.16) are solved numerically using the finite difference scheme introduced by Krasnov et al. (2011), and later developed, tested and applied to high-Ha flows and flows with thermal convection in numerous works including those by Krasnov et al. (2012), Zhao & Zikanov (2012), Zikanov et al. (2013), Zhang & Zikanov (2014), and Gelfgat & Zikanov (2018). The spatial discretization is of the second order and nearly fully conservative with regards to the mass, momentum, electric charge, kinetic energy and thermal energy conservation principles (Ni et al. 2007;Krasnov et al. 2011). The computational grid is clustered towards the walls according to the coordinate transformation in the horizontal direction y = tanh(A y η)/ tanh(A y ) and in the vertical direction z = tanh(A z ξ)/ tanh(A z ). Here η and ξ are the transformed coordinates, in which the grid is uniform, and A y and A z are the coefficients determining the degrees of clustering. The time discretization is implicit for the conduction and viscosity terms and based on the Adams-Bashforth/backward-differentiation method of the second order and the standard projection algorithm (see, e.g., Zikanov 2019). The nonlinear convection and body force terms are treated explicitly. The elliptic equations for potential, pressure, temperature and velocity components are solved using the Fourier decomposition in the streamwise coordinate and the direct cyclic reduction solution of the 2D equations for Fourier components conducted on the transformed grid (see Krasnov et al. 2011).
The algorithm is parallelized using the hybrid MPI-OpenMP approach. The MPI memory distribution is along the y−coordinate in the physical space and along the streamwise wavenumber in the Fourier space.

Approach to linear stability analysis
The base flow needs to be selected before conducting the linear stability analysis. We note that an archetypal structure of a laminar flow with convection in a horizontal channel heated from below is a superposition of the streamwise flow u x (y, z) and one or several streamwise-uniform convection rolls (u y (y, z), u z (y, z)). At high Ha, the structure can be modified by the magnetic field and replaced by a 3D structure with the rolls aligned with the magnetic field and, thus, x-dependent velocity and temperature at high Ha. Following Zikanov et al. (2013) and Zhang & Zikanov (2014), we treat the problem as that of the instability of the laminar steady-state streamwise-uniform base flow U (y, z), Θ(y, z), P (y, z) to x-dependent perturbations.
The base flow is calculated by artificially imposing uniformity in the streamwise direction, i.e. by applying x-averaging after every time step. In order to assure that a fully developed state of the base flow is reached, each solution is computed for sufficiently long time. Long evolution, in some regimes up to 1000 time units, is typically required in order to arrive at this state. No unsteady base flow solutions have been detected in the studied range of parameters. The steady-state solutions are discussed in section 3.1.
The linear stability analysis is conducted using a modified version of the numerical model described in section 2.3. We follow evolution of perturbations -solutions of the equations linearized around the base flow U (y, z), Θ(y, z), P (y, z). Individual Fourier modes determined by their streamwise wavelength λ are computed. This is practically achieved by setting the length of the computational domain to λ and filtering out all the Fourier modes except the zero mode corresponding to the base flow and the first mode corresponding to the perturbations of wavelength λ. All simulations start with random noise distributions of velocity and temperature.
The linear instability is identified by the exponential growth of the perturbations with the growth rate determined as where E = f 2 , . . . stands for volume averaging, and f stands for perturbations of a velocity component or temperature. The growth rate coefficient is recorded after its values computed for all three velocity components and temperature coincide with each other and remain constant within the third digit after the decimal point for at least 100 time units. The results of the linear stability analysis are presented in sections 3.3 and 3.4.

Grid sensitivity study
The grid sensitivity study has been conducted for the base flow. A detailed description of the various flow regimes is provided in section 3.1. For the present discussion, it is sufficient to say that accurate resolution of the internal flow structure, along with two boundary layers: the Hartmann layers of thickness δ Ha ∼ Ha −1 at the vertical walls and the Shercliff layers of thickness δ Sh ∼ Ha −1/2 at the top and bottom walls, is critically important for accurate representation of the flow behavior.
As an example, the results obtained at Ha = 1200, Gr = 10 8 are presented in table 1. In a fully developed steady-state flow, the integrated Lorentz and buoyancy forces are zero. The wall friction must be balanced by the driving pressure gradient according to where τ Ha and τ Sh are the computed values of the integrated friction forces at the Hartmann and Shercliff walls of the duct, respectively, expressed as Values of τ Ha , τ Sh and the error , with which the computed solution satisfies (2.18), found on various grids are compared in table 1. On the basis of these data, we conclude that the grid with N y × N z = 192 × 96, A y = 4.0 and A z = 2.0 is sufficient. The maximum and minimum grid steps of such a grid are ∆y min ≈ 0.0001, ∆y max ≈ 0.042, ∆z min ≈ 0.0009, ∆z max ≈ 0.012. The Hartmann and Shercliff layers are resolved by, respectively, 9 and 32 grid points.
The parameters of the grids in the entire studied parameter range of Ha and Gr have been determined in the same way. It has been found that the value of Gr does not affect the selection at fixed Ha. This effect can be explained by the presence of very strong magnetic fields which fully suppress transverse circulation (see section 3.1 for a discussion). The summary of the grids used in the simulations is presented in table 2.
A grid sensitivity study has also been conducted to determine the minimum number of grid points N x required in the linear stability analysis. It has been found that the growth of linear unstable modes is accurately reproduced at N x = 32 for modes with λ 2, which needs to be increased to N x = 64 for greater λ.
Computational domains of length L x = 4π or L x = 2π are used in DNS. As we will see below, these lengths are substantially larger than the streamwise wavelength of the fastest growing instability modes. The flow structures have been accurately resolved with, respectively, 384 or 192 grid points in the x−direction.
The time steps adjusted to secure numerical instability and, thus, varying with Ha and Gr, but never exceeding 2.5 × 10 −3 , are used in linear stability and DNS simulations.

Base flow
The structure of the base flow, as it is defined in section 2.4, for several typical cases is illustrated in figures 2 and 3. The results for all the completed simulations are summarized in table 3. The table shows the type of the flow for a particular regime (Q2D or 3D flows to be discussed shortly), the maximum and minimum values of U x , integral quantities, such as the wall friction force dp/dx, the volume-averaged kinetic energies of streamwise and transverse velocities and the mean square of temperature perturbations Similar data for flows with lower values of Gr and Ha in a duct with Γ = 1 can be found in (Zhang & Zikanov 2014).
The flow structure is predominantly determined by the effect of magnetoconvection. Similarly to the findings of Zhang & Zikanov (2014), we observe two regimes of the base flow depending on whether Gr is smaller or larger than a certain threshold Gr * (Ha). The Q2D regime observed at Gr < Gr * is characterized by the transverse convectioninduced circulation entirely suppressed by the strong magnetic field. The distributions of the temperature and streamwise velocity are nearly one-dimensional outside of the Hartmann boundary layers (Θ ≈ Θ(z) and U x ≈ U x (z)). In the absence of transverse circulation, the distributions of temperature are determined by the balance between the heat conduction and the heat convection by U x . Examples of this regime are shown in figure 2a and figure 3.
The 3D regime observed at Gr > Gr * (Ha), when the strength of the magnetic field is insufficient to suppress convection circulation, is characterized by significant transverse flow and fully 2D variations of temperature. (see figure 2b − e).
To avoid confusion, it is pertinent to repeat the terminology here. The terms 3D and Q2D are used in this paper to describe the general flow transformation caused by the magnetic field, i.e. suppression of velocity and temperature gradients along the magnetic field lines in the core of the duct and formation of thin Hartmann boundary layers. The base flow, in which streamwise uniformity is also imposed, becomes, respectively, 2D and Q1D.
We find the Q2D regime in the larger part of the explored range of Ha and Gr including the most interesting cases of large Ha (see the rightmost column in table 3). The total friction force increases at stronger magnetic fields. The visible effect of convection is the asymmetry of the velocity profile, with U x larger in the bottom than in the top half (see, e.g., figure 3b, d). The cause of the asymmetry has been explained in section 2.1.
Our interest in this study is primarily in the Q2D regimes. The 3D regimes are not considered in the rest of the paper.

Applicability of the SM82 model
In this section, we investigate the applicability of the SM82 model to analysis of magnetoconvection instability.

Base flow
For a streamwise-uniform, unidirectional, steady-state base flow, the SM82 model equations (2.14)−(2.16) are reduced to a system of linear ordinary differential equations. The solution satisfying the boundary conditions is: The symbols 2D and 3D correspond to, respectively, approximate SM82 and computed solutions.
Good agreement between the computed solutions and the solutions of the SM82 model is evident at Gr = 10 8 and 10 9 . There are some deviations between the computed and model profiles of U x , but they are small and decrease with increasing Ha. The agreement is significantly worse in the case of the flows at Gr = 10 10 . Here we observe large deviations between the computed and model curves at Ha = 5000 and smaller, yet still significant deviations at Ha = 10000.
The main reason for the discrepancy between the results of 2D model and 3D calculations at such a high Ha is the geometry of the flow. The duct has a large aspect ratio Γ and the magnetic field oriented along the long side. Deviations from two-dimensionality become more pronounced in such geometries. To verify this explanation, we performed additional simulations, which revealed that velocity and temperature profiles in 2D and 3D solutions are almost indistinguishable from each other at Γ 1.

Linear stability analysis
In order to verify applicability of the SM82 model to linear stability analysis of Q2D flows, instability of several high−Ha flows was evaluated twice: once using full 3D model of the base flow and perturbations and once entirely in the framework of the 2D SM82 model. The results are presented in figure 6 and table 4. We see good agreement between predictions of 2D and 3D models at Gr = 10 8 and 10 9 . The accuracy improves with growing Ha. As an example, the average relative difference between the values of the growth rate γ for the two models is 33% at Ha = 2000, Gr = 10 9 , 11% at Ha = 3000, Gr = 10 9 and 3% at Ha = 10000, Gr = 10 9 . The situation is less clear for flows at Gr = 10 10 (see figure 6c and the last six lines of table 4). Here we only see a qualitative agreement. The shape of the γ(λ) curves, the wavelength of the most unstable mode, and the effect of Ha on stability are similar in the 3D and 2D solution. The quantitative agreement is, however, poor, with the difference between the values of γ found for the two models being about 50%.
The quantitative disagreement between the base flow profiles and, as an evident consequence, stability properties found in the 3D and 2D models is difficult to interpret. The velocity and temperature distributions computed in the framework of the 3D model clearly show that the base flow is Q2D and nearly perfectly unidirectional at Ha higher than approximately 4000 (see figure 3b, d and values of E t in table 3). As illustrated in section 3.3, fields of growing perturbations also remain Q2D at such high Ha. Additional calculations performed with larger grids and longer times of flow evolution did not lead to significant changes. The deviations from quasi-two-dimensionality and inaccuracy of the numerical model are, therefore, excluded as possible reasons.
A further useful, albeit not fully explaining illustration is provided in figure 7. Computed base flow distributions of U x and the streamwise component of the Lorentz force F Lx are shown for Ha = 10000 and Gr = 10 8 , 10 9 and 10 10 within and near the Hartmann boundary layer. We see that the strong vertical variation of U x existing at Gr = 10 10 extends toward the Hartmann wall and causes a respective variation of the Lorentz force. It must be noted that this picture does not contradict to the identification of the flow as Q2D. The profiles U x (y, z = const), if taken outside the sidewall layers at the horizontal walls and scaled by the respective maximum values of U x , collapse into one curve with a flat core and Hartmann boundary layers.

Nonlinear flows
The results presented so far in this section indicate that the SM82 2D model may also inaccurately describe the nonlinear flow regimes developing as a result of the instability at Gr = 10 10 . As a test of this possibility, comparison between the results of 2D and 3D models at Gr = 10 10 , Ha = 10000 is illustrated in figures 8 and 9 and discussed below. The procedure of computing nonlinear flows is described in section 3.4. Here we only mention that the same numerical resolution is used in 2D and 3D models. The shorter wavelength domain length L x = 2π is used in the 3D model. This rather small length has no significant effect on the flow evolution as it has been confirmed in the additional 2D simulations conducted with L x = 4π (shown in figure 8) and L x = 2π.
The simulations show that the 3D flow remains Q2D at these values of Ha and Gr. At the same time, its reproduction by the 2D SM82 model is inaccurate in some aspects. The instantaneous distributions of velocity components and temperature shown in figure 8 clearly illustrate the difference. 2D approximations (see figure 8a-c) are similar to 3D flows (see figure 8d-f ) in terms of the largest typical streamwise wavelength (about 1.5) but demonstrate noticeably less regular pattern and higher energy in shorter wavelengths. The difference is reflected by the point signals of velocity and temperature shown in figure  9. It is also observed in the power spectrum density graphs of velocity and temperature (not shown). Comparing figures 9b and 9d, we also see that the 2D approximation substantially underestimates the typical amplitude of velocity fluctuations. As an example, the standard deviation for the signal of u   the 2D approximation are about an order of magnitude lower than in the actual flow computed in the framework of the 3D model. Interestingly, the misrepresentation of the structure of velocity and temperature fluctuations by the 2D model does not lead to a similarly strong inaccuracy in the prediction of the effect of mixing by fluctuations. Profiles of average streamwise velocity and temperature in figures 8(g, h) show strong change in comparison with the base flow, but only moderate differences between the 2D and 3D results.

Applicability of the SM82 model: summary
We conclude that the SM82 approximation accurately represents Q2D flows at moderately large Gr (10 8 and 10 9 in our system). The accuracy deteriorates at higher Gr even though the flow remains Q2D. An example of this is observed at Ha 4000 and Gr = 10 10 . The base flow profiles are clearly different between the 2D and 3D models (see figures 5c, f ). The 2D linear stability analysis is qualitatively correct in the sense that it correctly predicts the principal type of the unstable perturbations and the streamwise wavelength of the most unstable mode (see figure 6c). The values of the growth rate γ are, however, substantially underestimated by the 2D model (see figure 6c and table 4). The nonlinear flow states resulting from the instability are predicted incorrectly by the 2D model, which adds artificial irregularity and short-wave fluctuations and underestimates the amplitude of velocity fluctuations (see figures 8 and 9). It should be noted that the discrepancy is not due to irregularities of the model or our computational procedure. Calculations carried out at lower Gr and high Ha reveal a perfectly good agreement in the bulk region between the 2D and 3D models.
We do not have a satisfactory explanation to this effect and leave its further exploration for future studies. It should be mentioned that the quality of the 2D approximation It must be mentioned that this is not the only example of the model's breakdown. The model is known to break down when any 3D structures are present in the flow for which the diffusion length is shorter than the size of the domain (Pothérat & Klein 2014.
In the remaining part of this paper, the discussion of the instability and nonlinear states is based on the SM82 2D approximation for flows at Gr = 10 8 and 10 9 and on the full 3D solutions for flows at 10 10 .
The numerical solution for the 2D approximation is obtained using a modified version of the code which solves (2.14)−(2.16). The code has been verified through comparison of its results with the analytical solution of (3.3)−(3.7) for the base flow.

Results of linear stability analysis
The results of the linear stability analysis are summarized in figures 10−12 and table 5. We need to mention that the wavelength λ is varied with step 0.1 in the simulations.
The computed values of the exponential growth rate γ as a function of the wavelength λ for various combinations of Ha and Gr are shown in figure 10. Two trends of the linear stability behaviour were proposed by Zhang & Zikanov (2014): (1) a higher growth rate and shorter wavelength appear at higher Gr; (2) an increase of Ha leads to a higher growth rate. The second, apparently counterintuitive effect was attributed by Zhang &  Zikanov (2014) to modification of the base flow, namely to suppression of the transverse circulation resulting in weaker mixing and stronger unstable temperature stratification.
In order to investigate these trends for our system and parameter range, we present the exponential growth rate γ max of the fastest growing modes and the corresponding wavelengths λ max in figure 11 and table 5. Our results are clearly consistent with the first trend, γ max increases and λ max decreases with growing Gr. However, we find a different behaviour in regard of the second trend. The increase of Ha leads to noticeable or slight decrease of the exponential growth rates at Gr = 10 8 or Gr = 10 9 , respectively (see figure 10c). For Gr = 10 10 , γ is nearly insensitive to the values of Ha.
We conclude that the counterintuitive behaviour of stronger instability at higher Ha is not observed in Q2D flows at high Ha considered in our study. It is, nevertheless, noteworthy that γ max does not decrease with Ha at Gr = 10 10 . It increases slightly at Ha 6000 and remains practically constant at higher Ha. The enhanced friction in the Hartmann boundary layers is compensated by another effect, the only plausible candidate for which is the strong modification of the streamwise velocity profile visible in figure 5c. U x strongly grows with Ha in the bottom half of the duct, i.e., its part with the strongest unstable temperature gradient.
The most unstable mode is oscillatory. This was also observed in the earlier works of Zikanov et al. (2013); Zhang & Zikanov (2014). Point signals of temperature and velocity oscillate in time with constant frequency, This is caused by the transport of the rolls by mean flow. We have computed the phase velocity as the ratio of the axial wavelength to the period of oscillations of a signal at a given point to illustrate this effect. We found that, similarly to findings of Zikanov et al. (2013) and Zhang & Zikanov (2014), it varies little with Ha and Gr for the most unstable modes, and has the value close to the mean velocity value 1.
The findings have critical implications for design and operation of the fusion reactor systems, since they indicate that strength of the convection instability is not diminished by strong magnetic fields at Gr 10 10 typical for reactor blankets.
The spatial structure of the unstable modes is illustrated in figure 12 for Gr = 10 8 , 10 9 , 10 10 and Ha = 10 4 . The structures are qualitatively similar to those found for Q2D instabilities by Zhang & Zikanov (2014). Consistent with the first trend mentioned above and with the base flow modification illustrated in figure 5 is the fact that the energy of growing perturbations becomes contained to the lower part of the duct at higher values of Gr.

Results of DNS of nonlinear flows
The results concerning the nonlinear flow regimes are illustrated in figures 8(d, e, f ), 13, 14, and 15. DNS approach based on direct solution of the nonlinear governing equations is utilized. The 2D SM82 model (2.14)−(2.16) and the computational domain of length L x = 4π are used for flows at Gr = 10 8 and 10 9 . Full 3D equations (2.1)−(2.3) are solved and the domain is reduced to L x = 2π at Gr = 10 10 . Other parameters of computational model are described in section 2.5. Each simulation starts with the streamwiseindependent base flow (see section 3.2) computed at the same Gr and Ha, to which random small-amplitude (∼ 10 −3 ) random perturbations of velocity and temperature are added.
The typical flow evolution is illustrated by the curves of average kinetic energy shown in figure 13. The flow reaches a fully developed state after the instability and initial development. The evolution of the fully developed flow is computed for at least 500 time units for the 2D model at Gr = 10 8 , 10 9 and at least 100 time units for the 3D model at Gr = 10 10 . At this stage, the integral parameters fluctuate around steady means (at Gr = 10 9 and 10 10 ) or remain steady (at Gr = 10 8 ). The amplitudes of the fluctuations are small at Gr = 10 9 and large, but still moderate at Gr = 10 10 .
Structure of fully developed flows is illustrated in figures 8 and 14. The velocity field shows finite-amplitude roll-like structures (hereafter referred to as rolls) resulting from the instability, which are superimposed on a streamwise-independent mean flow (see figures 8d and 14a, d). The rolls cause variations of temperature (see figures 8f , 14c and 14f ). Transport of the rolls by the mean flow is a known reason of MCFs in horizontal channels (Zikanov et al. 2013;Zhang & Zikanov 2014;Zikanov et al. 2021).
Comparison of the flow structures in figures 8d-f and 14 reveal the effect of the value of Gr on convection rolls. As anticipated, increase of Gr leads to higher nondimensional amplitude of the velocity fluctuations. This results in stronger vertical mixing as illustrated by the streamwise-and time-averaged profiles in figure 14g, h. In particular, a nearly uniform vertical distribution of average temperature with a thin (but still much thicker than the Shercliff layer) boundary layer at the bottom is observed at Gr = 10 10 .
As we discussed earlier, MCFs caused by the instability have potentially critical implications for design and operation of liquid-metal components of nuclear fusion reactors. The DNS results allow us, for the first time, to evaluate the properties of the MCFs at the high values of Gr and Ha corresponding to the actual reactor conditions.
In addition to the instantaneous temperature distributions in figures 8f , 14c and 14f , the discussion will be based on the point-signals of temperature measured at the top and bottom walls and in the middle of the duct (see figure 15). As discussed, e.g., by Zikanov et al. (2021), measuring such signals is the most reliable and commonly used tool for studying MCFs in experiments.
The evident conclusion from the DNS data is that MCFs are fully present in flows with Gr = 10 8 , 10 9 , 10 10 and the highest values Ha = 5 × 10 3 and 10 4 considered in this study. The fluctuations are observed in the entire duct. The temperature signals are regular and dominated by one or several low frequencies (the typical period is 2-3 non-dimensional 24 R. Akhmedagaev, O. Zikanov, and Y. Listratov  time units at Gr = 10 8 and 10 9 ). The signal is less regular and characterized by higher typical frequencies at Gr = 10 10 . Interestingly, the non-dimensional amplitude of the temperature fluctuations decreases noticeably with growing Gr. Comparison of the signals in the two columns of figure  15 demonstrates that value of Ha has practically no effect on the MCFs. This can be attributed to the effect of nonlinearity, which distributes energy of the fluctuations to a range of streamwise modes.
Considering the practical implications, it is interesting to evaluate the parameters of the MCFs in dimensional units. We will do that for the temperature signals at the bottom of the duct (z = −0.2857) assuming the duct half-width d = 5 cm and using the physical properties of PbLi at 573 K (Zikanov et al. 2021). The wall heat rate is q = 10.56 kW m −2 at Gr = 10 8 , q = 105.6 kW m −2 at Gr = 10 9 and q = 1056 kW m −2 at Gr = 10 10 . We find, by applying the temperature scale qd/κ, that the largest amplitude of fluctuations of wall temperature is in the range 5 − 6 K at Gr = 10 8 , 44 − 62 K at Gr = 10 9 , and somewhat unrealistic 180 − 300 K at Gr = 10 10 . The typical time period of the fluctuations is 5.64 s at Gr = 10 8 , 6.45 s at Gr = 10 9 , and 4.51 s at Gr = 10 10 for Ha = 10 4 . Similar evaluations have been done for the future experiments on the recently built experimental facility (see, e.g., Belyaev et al. 2017), in which liquid mercury flows in the duct with the half-width of d = 2.8 cm. The physical properties of Hg are taken at 303 K (Zikanov et al. 2021). The wall heat rate is q = 9.59 kW m −2 at Gr = 10 8 and q = 95.9 kW m −2 at Gr = 10 9 . The results of nonlinear simulations allow us to predict the largest amplitude of fluctuations of temperature in the middle of the duct. The amplitudes are in the range of 2 − 3 K at Gr = 10 8 and 12 − 18 K at Gr = 10 9 . The typical time period of the fluctuations is 1.1 s at Gr = 10 8 and 1.58 s at Gr = 10 9 for Ha = 1000.

Concluding remarks
We have analysed mixed convection in a liquid metal flow in a duct with bottom heating and transverse magnetic field. The analysis is extended to much higher values of Ha and Gr than the previous analysis of similar effects by Zhang & Zikanov (2014).
The main conclusion of our work is that magnetoconvective fluctuations appear at the parameters anticipated for operational regimes of blankets and divertors of future fusion rectors. The fluctuations are not suppressed or even significantly reduced in amplitude by the very strong magnetic field. The amplitude remains high, reaching tens or hundreds degrees K (depending on the value of Gr) in a typical duct geometry. This has significant far-reaching implications for mixing, heat and mass transfer, and structural integrity of reactor components. The most dangerous modes of the instability have the form of rolls localized in the lower half of the duct and having the streamwise wavelength measured in horizontal half-widths of the duct, approximately, between 0.8 and 1.4 at Gr = 10 8 , 0.6 at Gr = 10 9 , and 0.4 at Gr = 10 10 .
Another conclusion concerns applicability of the two-dimensional approximation by Sommeria & Moreau (1982) to flows with thermal convection. We have found that the approximation may become inaccurate at high values of Gr even though the flow remains quasi-two-dimensional. Full reasons of this phenomenon remain to be understood. One of the reasons is, clearly, the geometry of the flow. The 2D model tends to be less accurate if applied to flows in ducts with larger aspect ratios and the magnetic field parallel to the long side. In general, the conclusion is important as a warning against application of the model without a proper verification.
It is pertinent to stress that the conclusions must be considered as preliminary because