1. 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 Reference Ozoe2005; Davidson Reference Davidson2016). 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. Reference Abdou, Morley, Smolentsev, Ying, Malang, Rowcliffe and Ulrickson2015). 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 (quasi-2-D) 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 Reference Branover1978; Sommeria & Moreau Reference Sommeria and Moreau1982; Zikanov et al. Reference Zikanov, Krasnov, Boeck, Thess and Rossi2014; Davidson Reference Davidson2016). 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. Reference Genin, Zhilin, Ivochkin, Razuvanov, Belyaev, Listratov and Sviridov2011; Vetcha et al. Reference Vetcha, Smolentsev, Abdou and Moreau2013; Zikanov, Listratov & Sviridov Reference Zikanov, Listratov and Sviridov2013; Belyaev et al. Reference Belyaev, Sardov, Melnikov and Frick2021). The term magneto-convective fluctuations (MCFs) proposed for the phenomenon by Belyaev et al. (Reference Belyaev, Sardov, Melnikov and Frick2021) will be used in this paper. As discussed in detail in the review of Zikanov et al. (Reference Zikanov, Listratov, Razuvanov, Belyaev, Frick and Sviridov2021) 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. (Reference Zikanov, Listratov and Sviridov2013) and Zhang & Zikanov (Reference Zhang and Zikanov2014). 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 quasi-two dimensional (see, e.g. Smolentsev Reference Smolentsev2021; Zikanov et al. Reference Zikanov, Listratov, Razuvanov, Belyaev, Frick and Sviridov2021, 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 50 K 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. Reference Zikanov, Listratov, Razuvanov, Belyaev, Frick and Sviridov2021). 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. Reference Belyaev, Poddubnyi, Razuvanov and Sviridov2018) possibly under the condition of significantly reduced strength of the wall material (Kolmakov et al. Reference Kolmakov, Terent'ev, Prosvirnin, Chernov and Leont'eva-Smirnova2016). 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 an 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 Reference Kirillov and Muraviev1997). It is also important as an archetypal system, in which the MCFs were first identified (in Genin et al. Reference Genin, Zhilin, Ivochkin, Razuvanov, Belyaev, Listratov and Sviridov2011 and Zikanov et al. Reference Zikanov, Listratov and Sviridov2013, where they were named anomalous fluctuations) and explained.
Figure 1. Flow geometry and coordinate system. The arrows marked by letters 
$\boldsymbol {g}$, 
$\boldsymbol {B}$ and 
$q$ denote, respectively, the orientations of the gravity acceleration, magnetic field and wall heating.
Similar systems for either ducts or round pipes have been studied experimentally (Genin et al. Reference Genin, Zhilin, Ivochkin, Razuvanov, Belyaev, Listratov and Sviridov2011; Belyaev et al. Reference Belyaev, Ivochkin, Listratov, Razuvanov and Sviridov2015; Sahu et al. Reference Sahu, Courtessole, Ranjan, Bhattacharyay, Sketchley and Smolentsev2020) and numerically (Zikanov et al. Reference Zikanov, Listratov and Sviridov2013; Zhang & Zikanov Reference Zhang and Zikanov2014; Vo, Pothérat & Sheard Reference Vo, Pothérat and Sheard2017; Listratov et al. Reference Listratov, Ognerubov, Zikanov and Sviridov2018). The flow is controlled by four dimensionless parameters: the Reynolds, Prandtl, Grashof and Hartmann numbers,
\begin{equation} \textit{Re} = \frac{Ud}{\nu},\quad \textit{Pr} = \frac{\nu}{\chi}, \quad \mathit{Gr} = \frac{g\beta q d^4}{\nu^2\kappa}, \quad \mathit{Ha} = Bd \sqrt{{\frac{\sigma}{\rho\nu}}}, \end{equation}
with the duct half-width 
$d$, the mean streamwise velocity 
$U$, the kinematic viscosity 
$\nu$, the temperature diffusivity 
$\chi$, the acceleration due to gravity 
$g$, the coefficient of thermal expansion 
$\beta$, the heat flux of constant rate 
$q$, the thermal conductivity 
$\kappa$, the electrical conductivity 
$\sigma$ and the mass density 
$\rho$. Rectangular duct geometry adds the aspect ratio 
$\varGamma = {2d}/{h}$ as a parameter, where 
$h$ is the height of the duct.
In the linear stability analysis of the Poiseulle–Rayleigh–Bénard duct flow, interesting results were obtained with a transverse magnetic field performed by Vo et al. (Reference Vo, Pothérat and Sheard2017). Two-dimensional (2-D) approximation valid in the limit of strong magnetic field presented later in this paper was used. One important result of Vo et al. (Reference Vo, Pothérat and Sheard2017) 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 (
${\sim }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. Reference Genin, Zhilin, Ivochkin, Razuvanov, Belyaev, Listratov and Sviridov2011; Belyaev et al. Reference Belyaev, Ivochkin, Listratov, Razuvanov and Sviridov2015) and explained in the linear stability analysis and direct numerical simulations (DNS) by Zikanov et al. (Reference Zikanov, Listratov and Sviridov2013). Flows of mercury with 
$\textit {Pr} \approx 0.022$, 
$\textit {Re}$ up to 
$10^5$, 
$\mathit {Gr}$ up to 
$10^8$ and 
$\mathit {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 quasi-2-D rolls. Transport of the rolls by the mean flow generated the MCFs.
The analysis was extended in the numerical simulations of Zhang & Zikanov (Reference Zhang and Zikanov2014). Flows in a horizontal duct of aspect ratio 
$\varGamma = 1$ with bottom heating and a transverse magnetic field at 
$\textit {Pr} = 0.0321$, 
$\textit {Re} = 5000$, 
$50 \le \mathit {Ha} \le 800$ and 
$10^5 \le \mathit {Gr} \le 10^9$ were investigated. The instability leading to the formation of quasi-2-D rolls similar to those found in the pipe flow was detected at sufficiently high 
$Gr$ and 
$Ha$.
Investigations of Zhang & Zikanov (Reference Zhang and Zikanov2014) 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 
$\mathit {Gr}$ type characterized by quasi-2-D distributions of velocity and temperature dominated by spanwise rolls appeared at 
$Gr$ below a certain 
$Gr^{\ast }(Ha)$. At higher 
$\mathit {Gr}$, stronger convection resulted in three-dimensional (3-D) 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 \sim 10^{10} {-} 10^{12}$) and magnetic fields (
$Ha \sim 10^{4}$) (see, e.g. Smolentsev, Moreau & Abdou Reference Smolentsev, Moreau and Abdou2008; Smolentsev et al. Reference Smolentsev, Moreau, Bühler and Mistrangelo2010). Such extreme parameters present serious obstacles to analysis, because neither laboratory experiments nor 3-D 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 (Reference Zhang and Zikanov2014). The extrapolation predicted the existence of MCFs at the typical blanket parameters. It also predicted that the flow would likely be of the low 
$\mathit {Gr}$ type at 
$\mathit {Gr} \le 10^{10}$ and of the high 
$\mathit {Gr}$ type at higher 
$\mathit {Gr}$. The experiments in the pipe flow (see the review of recent results in Zikanov et al. Reference Zikanov, Listratov, Razuvanov, Belyaev, Frick and Sviridov2021), on the contrary, indicate that MCFs may disappear at high 
$\mathit {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 analyse the MCF effect in this range. Linear stability analysis and DNS of flows in a horizontal duct with 
$\textit {Pr} = 0.025$, 
$Re = 5000$, 
$10^8 \le \mathit {Gr} \le 10^{10}$ and 
$10^3 \le \mathit {Ha} \le 10^4$ are performed. The study follows the work of Zhang & Zikanov (Reference Zhang and Zikanov2014) but differs by much larger values of 
$Gr$ and 
$Ha$ and the aspect ratio 
$\varGamma = 3.5$ selected to match the new experimental facility (see, e.g. Belyaev et al. Reference Belyaev, Sviridov, Batenin, Biryukov, Nikitina, Manchkha, Pyatnitskaya, Razuvanov and Sviridov2017), on which the same configuration is to be explored at 
$\mathit {Ha} \lesssim 10^3$ and 
$\mathit {Gr} \lesssim 10^8$ in the near future. Another essential difference between our work and the work by Zhang & Zikanov (Reference Zhang and Zikanov2014) is that we carry out an in-depth analysis of the accuracy of the 2-D approximation applied to quasi-2-D flows at such high 
$Ha$.
2. 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 
$\varGamma = 3.5$ (see figure 1). A spatially uniform and time-independent magnetic field 
$\boldsymbol {B}=B\boldsymbol {e}_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.
2.1. 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 Reference Davidson2016). The approximation implies that the imposed magnetic field 
$\boldsymbol {B}$ is much stronger than the perturbations of the magnetic field 
$\boldsymbol {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/\kappa$ as the temperature scale, 
$B$ as the scale of the magnetic field strength and 
$dUB$ as the scale of electric potential. The equations can be written as
\begin{gather} \frac{\partial \boldsymbol{u}}{\partial t} + (\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla})\boldsymbol{u} ={-}\boldsymbol{\nabla} p -\boldsymbol{\nabla} \hat{p} -\boldsymbol{\nabla} \tilde{p} + \frac{1}{Re} \nabla^2 \boldsymbol{u} + \boldsymbol{F}_b + \boldsymbol{F}_L, \end{gather}
\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0, \end{gather}
\begin{gather}\frac{\partial \theta}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \theta = {\frac{1}{\textit{Re}\textit{Pr}}} \nabla^2 \theta - u_x \frac{\text{d}T_m}{\text{d}\,x}, \end{gather}
where 
$\boldsymbol {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, Garandet & Msoreau Reference Alboussière, Garandet and Msoreau1993; Lyubimova et al. Reference Lyubimova, Lyubimov, Morozov, Scuridin, Hadid and Henry2009; Zikanov et al. Reference Zikanov, Listratov and Sviridov2013; Zhang & Zikanov Reference Zhang and Zikanov2014; Zikanov et al. Reference Zikanov, Listratov, Razuvanov, Belyaev, Frick and Sviridov2021) 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
\begin{equation} T(\boldsymbol{x},t) = \it {T_m}(x) + \theta(\boldsymbol{x},t) \end{equation}
of fluctuations 
$\theta$ and the mean-mixed temperature
\begin{equation} T_m(x) = \frac{\displaystyle\int_{A} u_xT \,\text{d}A}{\displaystyle\int_{A} u_x \,\text{d}A}= A^{{-}1} \int_{A} u_xT \,\text{d}A , \end{equation}
where 
$A = 2 h/d$ is the cross-section area of the duct. One can also use the decomposition into fluctuations and simple mean temperature 
$\bar {T}(x)=A^{-1}\int _{A} T \,\text {d}A$. Applying the energy balance between the wall heating and the streamwise convection heat transfer, we find that 
${T_m}(x)$ and 
$\bar {T}(x)$ are linear functions with the same derivative, i.e.
\begin{equation} \frac{\text{d}T_m}{\text{d}\,x} = \frac{\text{d}\bar{T}}{\text{d}\,x} = \frac{\varPi}{A RePr} = \frac{\varGamma}{2RePr} , \end{equation}
where 
$\varPi = 2$ is the perimeter of the heated portion of the wall.
The total pressure 
$P$ is presented in (2.1) as
\begin{equation} P = \hat{p}(x) + \tilde{p}(x,z) + p(\boldsymbol{x},t), \end{equation}
where 
$p(\boldsymbol {x},t)$ is the field of pressure fluctuations statistically homogeneous in the streamwise direction, and 
$\hat {p}$ is a linear function of 
$x$ corresponding to the spatially uniform streamwise gradient 
${\text {d}{\hat {p}}}/{\text {d}x}$ 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
\begin{equation} \tilde{p}(x,z) = \frac{\text{d}T_m}{\text{d}\,x}\frac{Gr}{Re^2} xz = \frac{\varPi}{A\textit{Re}\textit{Pr}}\frac{Gr}{Re^2} xz \end{equation}
arises due to the buoyancy force caused by the mean-mixed temperature 
$T_m$,
\begin{equation} \boldsymbol{F}_{b, m} = Gr Re^{{-}2} \boldsymbol{e}_z T_m, \end{equation}
where 
$\boldsymbol {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 field 
$\tilde {p}$, such that its vertical gradient balances 
$\boldsymbol {F}_{b, m}$. The pressure field is a 2-D 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 top–bottom 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. Reference Zikanov, Listratov and Sviridov2013; Zhang & Zikanov Reference Zhang and Zikanov2014, Reference Zhang and Zikanov2017; Zikanov et al. Reference Zikanov, Listratov, Razuvanov, Belyaev, Frick and Sviridov2021).
The buoyancy force in (2.1) is
\begin{equation} \boldsymbol{F }_b = Gr Re^{{-}2} \boldsymbol{e}_z T. \end{equation}The Lorentz force is computed as
\begin{equation} \boldsymbol{F}_L = Ha^2 Re^{{-}1} \boldsymbol{j}\times \boldsymbol{e}_y, \end{equation}
where 
$\boldsymbol {e}_y$ is the unit vector along the imposed magnetic field (see figure 1). The electric current 
$\boldsymbol {j}$ is determined by the Ohm's law
\begin{equation} \boldsymbol{j} ={-}\boldsymbol{\nabla} \phi + (\boldsymbol{u} \times \boldsymbol{e}_y), \end{equation}
where the electric potential 
$\phi$ is a solution of the Poisson equation expressing the instantaneous electric neutrality of the fluid,
\begin{equation} \nabla^2\phi = \boldsymbol{\nabla} \boldsymbol{\cdot} (\boldsymbol{u} \times \boldsymbol{e}_y). \end{equation} The inlet-exit conditions are those of periodicity of the velocity 
$\boldsymbol {u}$, temperature fluctuations 
$\theta$, pressure fluctuations 
$p$ and potential 
$\phi$.
2.2. Two-dimensional approximation
Flows with a very strong imposed magnetic field are considered, so the Hartmann number and the Stuart number satisfy 
${Ha} \gg 1$ and 
$N \equiv Ha^2/Re\gg 1$, respectively. The flows are anticipated to have quasi-2-D form with nearly zero gradients along the magnetic field lines except in the thin Hartmann layers at the walls perpendicular to the field. The 2-D approximation proposed by Sommeria & Moreau (Reference Sommeria and Moreau1982) can be applied in this asymptotic limit. The problem can be expressed in terms of the variables integrated wall-to-wall along the direction of the magnetic field, leading to 2-D dynamics for 
$y$-averaged quantities. The approximation has been verified and examined by Pothérat, Sommeria & Moreau (Reference Pothérat, Sommeria and Moreau2000, Reference Pothérat, Sommeria and Moreau2005), and utilized in numerical studies of liquid metal flows in rectangular ducts (see, e.g. Pothérat Reference Pothérat2007; Smolentsev, Vetcha & Moreau Reference Smolentsev, Vetcha and Moreau2012; Vetcha et al. Reference Vetcha, Smolentsev, Abdou and Moreau2013; Vo et al. Reference Vo, Pothérat and Sheard2017; Zhang & Zikanov Reference Zhang and Zikanov2018).
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 2-D solutions, while the full solutions obtained numerically without resorting to the model will be designated as 3-D solutions.
The SM82 model is derived for flows with 
${Ha} \gg 1$ and 
$N \gg 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 quasi-2-D 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. (Reference Smolentsev, Moreau and Abdou2008). As demonstrated in this and in the following studies (see, e.g. Gelfgat & Molokov Reference Gelfgat and Molokov2011; Vetcha et al. Reference Vetcha, Smolentsev, Abdou and Moreau2013; Zhang & Zikanov Reference Zhang and Zikanov2018), the model can be extended to a 2-D approximation of temperature if the imposed heat flux is perpendicular to the magnetic field.
The SM82 version of (2.1)–(2.3) is
\begin{gather} \frac{\partial \boldsymbol{u}}{\partial t} + (\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla})\boldsymbol{u} ={-}\boldsymbol{\nabla} p -\boldsymbol{\nabla} \hat{p} -\boldsymbol{\nabla} \tilde{p} + \frac{1}{Re} \nabla^2 \boldsymbol{u} - \frac{Ha}{Re}\boldsymbol{u} + \frac{Gr}{Re^2}T\hat{\boldsymbol{e}}_z, \end{gather}
\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0, \end{gather}
\begin{gather}\frac{\partial \theta}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \theta = {\frac{1}{\textit{Re}\textit{Pr}}} \nabla^2 \theta - u_x \frac{\text{d}T_m}{\text{d}\,x}, \end{gather}
where all the flow variables are now 2-D fields obtained by wall-to-wall averaging. The term 
$({Ha}/{Re})\boldsymbol {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 3-D model.
2.3. 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, Zikanov & Boeck (Reference Krasnov, Zikanov and Boeck2011), and later developed, tested and applied to high 
${Ha}$ flows and flows with thermal convection in numerous works including those by Krasnov, Zikanov & Boeck (Reference Krasnov, Zikanov and Boeck2012), Zhao & Zikanov (Reference Zhao and Zikanov2012), Zikanov et al. (Reference Zikanov, Listratov and Sviridov2013), Zhang & Zikanov (Reference Zhang and Zikanov2014) and Gelfgat & Zikanov (Reference Gelfgat and Zikanov2018). 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. Reference Ni, Munipalli, Huang, Morley and Abdou2007; Krasnov et al. Reference Krasnov, Zikanov and Boeck2011). The computational grid is clustered towards the walls according to the coordinate transformation in the horizontal direction 
$y = \tanh (A_y \eta )/\tanh (A_y)$ and in the vertical direction 
$z = \tanh (A_z \xi )/\tanh (A_z)$. Here 
$\eta$ and 
$\xi$ 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 Reference Zikanov2019). 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 2-D equations for Fourier components conducted on the transformed grid (see Krasnov et al. Reference Krasnov, Zikanov and Boeck2011).
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.
2.4. 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 3-D structure with the rolls aligned with the magnetic field and, thus, 
$x$-dependent velocity and temperature at high 
$Ha$. Following Zikanov et al. (Reference Zikanov, Listratov and Sviridov2013) and Zhang & Zikanov (Reference Zhang and Zikanov2014), we treat the problem as that of the instability of the laminar steady-state streamwise-uniform base flow 
$\boldsymbol {U}(y,z)$, 
$\varTheta (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 a 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 § 3.1.
The linear stability analysis is conducted using a modified version of the numerical model described in § 2.3. We follow evolution of perturbations – solutions of the equations linearized around the base flow 
$\boldsymbol {U}(y,z)$, 
$\varTheta (y,z)$, 
$P(y,z)$. Individual Fourier modes determined by their streamwise wavelength 
$\lambda$ are computed. This is practically achieved by setting the length of the computational domain to 
$\lambda$ 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 
$\lambda$. 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
\begin{equation} \gamma = \frac{1}{2E'}\frac{\textrm{d}E'}{\textrm{d}t}, \end{equation}
where 
$E'=\langle f^2\rangle$, 
$\langle \cdots \rangle$ 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 §§ 3.3 and 3.4.
2.5. 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 § 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 
$\delta _{Ha} \sim \mathit {Ha}^{-1}$ at the vertical walls and the Shercliff layers of thickness 
$\delta _{Sh} \sim \mathit {Ha}^{-1/2}$ at the top and bottom walls, is critically important for accurate representation of the flow behaviour.
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
\begin{equation} \frac{\textrm{d}\hat{p}}{\textrm{d} x} = A^{{-}1}(\tau_{Ha}+\tau_{Sh}), \end{equation}
where 
$\tau _{Ha}$ and 
$\tau _{Sh}$ are the computed values of the integrated friction forces at the Hartmann and Shercliff walls of the duct, respectively, expressed as
\begin{equation} \tau_{Ha}=\tau_y = \frac{1}{Re}\sum_{y={\pm} 1}\int_{{-}1/\varGamma}^{1/\varGamma}\frac{\partial U_x}{\partial y}\,\textrm{d}z,\quad \tau_{Sh}=\tau_z = \frac{1}{Re}\sum_{z={\pm} 1/\varGamma}\int_{{-}1}^1\frac{\partial U_x}{\partial z}{\textrm{d} y}. \end{equation}
Values of 
$\tau _{Ha}$, 
$\tau _{Sh}$ and the error 
$\epsilon$, 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 \times N_z = 192 \times 96$, 
$A_y = 4.0$ and 
$A_z = 2.0$ is sufficient. The maximum and minimum grid steps of such a grid are 
$\Delta y_{min}\approx 0.0001$, 
$\Delta y_{max}\approx 0.042$, 
$\Delta z_{min} \approx 0.0009$, 
$\Delta z_{max} \approx 0.012$. The Hartmann and Shercliff layers are resolved by, respectively, 
$9$ and 
$32$ grid points.
Table 1. Grid sensitivity study conducted for 
$Ha = 1200$ and 
$Gr = 10^{8}$; 
$\tau _{Ha}$ and 
$\tau _{Sh}$ are the wall friction forces, 
$\epsilon$ is the absolute error of the balance (2.18). The number of grid points inside the Hartmann and Shercliff boundary layers are 
$N_{Ha}$ and 
$N_{Sh}$, respectively.
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 § 3.1 for a discussion). The summary of the grids used in the simulations is presented in table 2.
Table 2. Parameters of the computational grids used in the simulations for 
$Gr = 10^8 {-} 10^{10}$. The number of grid points inside the Hartmann and Shercliff boundary layers are 
$N_{Ha}$ and 
$N_{Sh}$, respectively.
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 
$\lambda \lesssim 2$, which needs to be increased to 
$N_x = 64$ for greater 
$\lambda$.
Computational domains of length 
$L_x = 4{\rm \pi}$ or 
$L_x = 2{\rm \pi}$ 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 \times 10^{-3}$, are used in linear stability and DNS simulations.
3. Results
3.1. Base flow
The structure of the base flow, as it is defined in § 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 flow for a particular regime (quasi-2-D or 3-D flows to be discussed shortly), the maximum and minimum values of 
$U_x$, integral quantities, such as the wall friction force 
$\textrm {d}\hat {p}/{\textrm {d} x}$, the volume-averaged kinetic energies of streamwise and transverse velocities
\begin{equation} E_x = A^{{-}1}\int_{{-}1/\varGamma}^{1/\varGamma}\int_{{-}1}^1 U_x^2\,\textrm{d}y\,\textrm{d}z,\quad E_t = A^{{-}1}\int_{{-}1/\varGamma}^{1/\varGamma}\int_{{-}1}^1 (U_y^2+U_z^2)\,\textrm{d}y\,\textrm{d}z, \end{equation}and the mean square of temperature perturbations
\begin{equation} E_{\theta} = A^{{-}1}\int_{{-}1/\varGamma}^{1/\varGamma}\int_{{-}1}^1 \varTheta^2 \,\textrm{d}y\,\textrm{d}z. \end{equation}
Figure 2. Base flow at 
$Ha = 1000$, 
$Gr = 10^8$ (
$a$), 
$Ha = 1000$, 
$Gr = 10^9$ (
$b$), 
$Ha = 1000$, 
$Gr = 10^{10}$ (
$c$), 
$Ha = 2000$, 
$Gr = 10^{10}$ (
$d$) and 
$Ha = 3000$, 
$Gr = 10^{10}$ (
$e$). Vector fields and streamlines of transverse circulation (
$u_y$, 
$u_z$) are shown in the left column (not in (
$a$), since the velocity's amplitude is virtually zero in this case). The middle and right columns show distributions of temperature 
$\varTheta$ and streamwise velocity 
$U_x$, respectively. Solid and dashed isolines in the middle column indicate positive and negative values, respectively. The wall heating is at 
$z = -0.2857$, and the magnetic field is in the 
$y$-direction.
Figure 3. Quasi-2-D base flow at 
$Ha = 5000$, 
$Gr = 10^9$ (
$a$), 
$Ha = 5000$, 
$Gr = 10^{10}$ (
$b$), 
$Ha = 10\,000$, 
$Gr = 10^9$ (
$c$) and 
$Ha = 10\,000$, 
$Gr = 10^{10}$ (
$d$). Distribution of streamwise velocity 
$u_x$ is shown. The wall heating is at 
$z = -0.2857$, and the magnetic field is in the 
$y$-direction.
Table 3. Integral characteristics and type of computed base flow states.
Similar data for flows with lower values of 
$Gr$ and 
$Ha$ in a duct with 
$\varGamma = 1$ can be found in Zhang & Zikanov (Reference Zhang and Zikanov2014).
The flow structure is predominantly determined by the effect of magnetoconvection. Similarly to the findings of Zhang & Zikanov (Reference Zhang and Zikanov2014), we observe two regimes of the base flow depending on whether 
$Gr$ is smaller or larger than a certain threshold 
$Gr^{\ast }(Ha)$. The quasi-2-D regime observed at 
$Gr < Gr^{\ast }$ is characterized by the transverse convection-induced 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 (
$\varTheta \approx \varTheta (z)$ and 
$U_x \approx 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 figures 2(
$a$) and 3.
The 3-D regime observed at 
$Gr > Gr^{\ast }(Ha)$, when the strength of the magnetic field is insufficient to suppress convection circulation, is characterized by significant transverse flow and fully 2-D variations of temperature. (see figure 2b–e).
To avoid confusion, it is pertinent to repeat the terminology here. The terms ‘3-D’ and ‘quasi-2-D’ 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, ‘2-D’ and ‘quasi-one-dimensional’.
We find the quasi-2-D 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 3
$b$,
$d$). The cause of the asymmetry has been explained in § 2.1.
The 3-D regimes are only found in a limited range of moderate values of 
$Ha$ at 
$Gr = 10^9$ and 
$10^{10}$ (see table 3). The transverse circulation consist of a single roll (see figure 2
$b$,
$e$) or two symmetric rolls (see figure 2
$c$,
$d$). The single roll has no preferred circulation direction and may appear in the solution either as shown in the figure 2(
$b$) or as a symmetric reflection with respect to the vertical midplane (see figure 2
$e$ for an example). The circulation causes visible 2-D distributions of 
$\varTheta$ and 
$U_x$, and, at the same values of 
$Gr$, a decrease of 
$E_\theta$ as a result of mixing (see table 3). We observe stronger top–bottom asymmetries or even formation of reverse flow in the top portion of the duct as the strength of convection increases at fixed 
$Ha$ (see figure 2
$a$–
$c$ for an example). As an illustration of the asymmetry, the minimum and maximum values of 
$u_x$ are shown in table 3.
The classification of the flow regimes into quasi-2-D and 3-D regimes can also be described in terms of the values of the average kinetic energy of transverse circulation 
$E_t$. Our data shown in table 3 and figure 4 are in a good qualitative agreement with the results of Zhang & Zikanov (Reference Zhang and Zikanov2014). The observed differences can be attributed to the substantially different studied ranges of 
$Ha$ and 
$Gr$ and different aspect ratio.
Figure 4. The average kinetic energy of transverse circulation in the base flow 
$E_t$ as a function of 
$Gr/Ha^2$. Circles indicate the numerical results of Zhang & Zikanov (Reference Zhang and Zikanov2014) at 
$\varGamma = 1.0$. Stars and crosses indicate, respectively, quasi-2-D and 3-D regimes found in this work for the flow at 
$\varGamma = 3.5$. Values of 
$E_t$ and 
$Gr/Ha^2$ for each computed flow can be found in table 3.
We see in figure 4 that, at a fixed Reynolds number considered in this work, the intensity of the transverse circulation is well approximated by a function of the single control parameter – the combination 
$Gr/Ha^2$. Analysing the flow structures at various values of 
$E_t$ we find a clear demarcation between 3-D and quasi-2-D regimes. Here 
$E_t$ is greater than, approximately, 
$10^{-4}$ in 3-D regimes. Quasi-2-D flows all have values of 
$E_t$ less than 
$10^{-6}$.
Our interest in this study is primarily in the quasi-2-D regimes. The 3-D regimes are not considered in the rest of the paper.
3.2. Applicability of the SM82 model
In this section we investigate the applicability of the SM82 model to analysis of magnetoconvection instability.
3.2.1. 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
\begin{gather} {U_x}(z) =\frac{Re}{Ha} \left( Cc_1 + \frac{A}{\varGamma}s_1 - Az - C \right), \end{gather}
\begin{gather}\varTheta(z) = \frac{Re}{2Ha^2} ( C\varGamma c_1 + A s_1 ) - \frac{Re}{Ha} \frac{\varGamma}{4} \left( \frac{A}{3} z^3 + C z^2 \right) + a_1z, \end{gather}where
\begin{gather} a_1 ={-} \frac{Re}{Ha^{3/2}} \frac{\varGamma}{2} ( C t_1 + A t_1^{{-}1} ) + \frac{Re}{2Ha} \left( \frac{A}{2\varGamma} + C \right), \end{gather}
\begin{gather}c_1 = \frac{ \cosh(\sqrt{Ha}z) } { \cosh(\sqrt{Ha}/\varGamma) } ,\quad s_1 = \frac{ \sinh(\sqrt{Ha}z)} { \sinh(\sqrt{Ha}/\varGamma) }, \quad t_1 = \tanh( \sqrt{Ha}/\varGamma ), \end{gather}
\begin{gather}A = \frac{\varGamma}{2RePr} \frac{Gr}{Re^2}, \quad C = \frac{Ha}{Re\varGamma} \frac{1}{t_1/\sqrt{Ha} - 1/\varGamma}. \end{gather} The profiles (3.3) and (3.4) are shown for 
$Gr = 10^8$, 
$10^9$, 
$10^{10}$ and several values of 
$Ha$ in figure 5. For comparison, solutions of the 3-D equations 
$U_x$ and 
$\varTheta$ obtained for the computed base flow solutions described in § 3.1 are shown for the midplane 
$y=0$. The terms 2-D and 3-D correspond to, respectively, approximate SM82 and computed solutions.
Figure 5. Base flow profiles at 
$Gr = 10^8$ (
$a$,
$d$), 
$Gr = 10^9$ (
$b$,
$e$) and 
$Gr = 10^{10}$ (
$c$,
$f$). The SM82 model solutions (3.3), (3.4) and the distributions along the midplane 
$y=0$ obtained in the full numerical solutions of § 3.1 are denoted, respectively, as 2-D and 3-D. The curves for 
$Ha = 10^3$ are only shown at 
$Gr = 10^8$ (
$a$,
$d$), since at higher 
$Gr$ the base flow is not quasi-two dimensional (see § 3.1). The top and bottom rows show the profiles of streamwise velocity and temperature, respectively. The insets in (
$d$) and (
$e$) show zoomed-in illustrations of the temperature profiles near the top wall.
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=10\,000$.
The main reason for the discrepancy between the results of the 2-D model and 3-D calculations at such a high 
$Ha$ is the geometry of the flow. The duct has a large aspect ratio 
$\varGamma$ 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 2-D and 3-D solutions are almost indistinguishable from each other at 
$\varGamma \le 1$.
3.2.2. Linear stability analysis
In order to verify applicability of the SM82 model to linear stability analysis of quasi-2-D flows, instability of several high 
$Ha$ flows was evaluated twice: once using the full 3-D model of the base flow and perturbations and once entirely in the framework of the 2-D SM82 model. The results are presented in figure 6 and table 4. We see good agreement between predictions of 2-D and 3-D 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 
$\gamma$ for the two models is 
$33\,\%$ at 
$Ha = 2000$, 
$Gr = 10^{9}$, 
$11\,\%$ at 
$Ha = 3000$, 
$Gr = 10^{9}$ and 
$3\,\%$ at 
$Ha = 10\,000$, 
$Gr = 10^{9}$ . The situation is less clear for flows at 
$Gr = 10^{10}$ (see figure 6
$c$ and the last six lines of table 4). Here we only see a qualitative agreement. The shape of the 
$\gamma (\lambda )$ curves, the wavelength of the most unstable mode, and the effect of 
$Ha$ on stability are similar in the 3-D and 2-D solution. The quantitative agreement is, however, poor, with the difference between the values of 
$\gamma$ found for the two models being about 
$50\,\%$.
Figure 6. Rates of exponential growth 
$\gamma$ shown as functions of the axial wavelength 
$\lambda$ at 
$Ha = 1000, 2000, 10\,000$, 
$Gr = 10^8$ (
$a$), 
$Ha = 2000, 3000, 10\,000$, 
$Gr = 10^9$ (
$b$) and 
$Ha = 4000, 5000, 10\,000$, 
$Gr = 10^{10}$ (
$c$). The results of 3-D and 2-D (SM82) models are denoted as filled and empty circles, respectively.
Table 4. Results of the linear stability analysis of 3-D and 2-D flow solutions. Rates of exponential growth 
$\gamma$ are shown as functions of the axial wavelength 
$\lambda$. The results of the SM82 model are underlined. The growth rates are determined as in (2.17).
The quantitative disagreement between the base flow profiles and, as an evident consequence, stability properties found in the 3-D and 2-D models is difficult to interpret. The velocity and temperature distributions computed in the framework of the 3-D model clearly show that the base flow is quasi-two dimensional and nearly perfectly unidirectional at 
$Ha$ higher than approximately 
$4000$ (see figure 3
$b$,
$d$ and values of 
$E_t$ in table 3). As illustrated in § 3.3, fields of growing perturbations also remain quasi-two dimensional 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 = 10\,000$ 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 quasi-two dimensional. The profiles 
$U_x(y, z = \textrm {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.
Figure 7. Base flow at 
$Gr = 10^8$ (
$a,b$), 
$Gr = 10^9$ (
$c$,
$d$) and 
$Gr = 10^{10}$ (
$e$,
$f$) for 
$Ha = 10\,000$. The left and right columns show, respectively, distributions of streamwise velocity 
$U_x$ and the streamwise component of the Lorentz force 
$F_{Lx}$ near the Hartmann wall. The red dashed line shows the boundary of the Hartmann layer of thickness 
$\delta _{Ha} = \mathit {Ha}^{-1}$.
3.2.3. Nonlinear flows
The results presented so far in this section indicate that the SM82 2-D 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 2-D and 3-D models at 
$Gr = 10^{10}$, 
$Ha = 10\,000$ is illustrated in figures 8 and 9 and discussed below. The procedure of computing nonlinear flows is described in § 3.4. Here we only mention that the same numerical resolution is used in 2-D and 3-D models. The shorter wavelength domain length 
$L_x = 2{\rm \pi}$ is used in the 3-D model. This rather small length has no significant effect on the flow evolution as it has been confirmed in the additional 2-D simulations conducted with 
$L_x = 4{\rm \pi}$ (shown in figure 8) and 
$L_x = 2{\rm \pi}$.
Figure 8. Flow structure in nonlinear regime at 
$Gr = 10^{10}$, 
$Ha = 10\,000$. The instantaneous distributions of 
$u_z$ (
$a$,
$d$), 
$u_x$ (
$b$,
$e$) and 
$\theta$ (
$c$,
$f$) obtained in 2-D and 3-D (plotted in the midplane 
$y = 0$) models are shown in (
$a$–
$c$) and (
$d$–
$f$), respectively. The profiles of 
$u_x$ and 
$\theta$ obtained by averaging over 
$x$ and time (with the base flow profiles from figure 5
$c$,
$f$) are shown, respectively, in (
$g$) and (
$h$).
Figure 9. Time signals of temperature (
$a$,
$c$) and streamwise velocity (
$b$,
$d$) at 
$y=0$ (in the 3-D flow) in the nonlinear regime at 
$Ha = 10^4$, 
$Gr = 10^{10}$ are shown for 2-D (
$a$,
$b$) and 3-D (
$c$,
$d$) models.
The simulations show that the 3-D flow remains quasi-two dimensional at these values of 
$Ha$ and 
$Gr$. At the same time, its reproduction by the 2-D SM82 model is inaccurate in some aspects. The instantaneous distributions of velocity components and temperature shown in figure 8 clearly illustrate the difference. Two-dimensional approximations (see figure 8
$a$–
$c$) are similar to 3-D flows (see figure 8
$d$–
$f$) in terms of the largest typical streamwise wavelength (about 
$1.5$) but demonstrate a 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 9(
$b$) and 9(
$d$), we also see that the 2-D approximation substantially underestimates the typical amplitude of velocity fluctuations. As an example, the standard deviation for the signal of 
$u_x$ in the middle of the duct (
$z = 0$) is 
$0.33$ for the 2-D model and 
$0.71$ for the 3-D model. Computed values of volume-average kinetic energy of the fluctuations (not shown) confirm this conclusion. The values of energy found in the 2-D approximation are about an order of magnitude lower than in the actual flow computed in the framework of the 3-D model. Interestingly, the misrepresentation of the structure of velocity and temperature fluctuations by the 2-D 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 figure 8(
$g$,
$h$) show a strong change in comparison with the base flow, but only moderate differences between the 2-D and 3-D results.
3.2.4. Applicability of the SM82 model: summary
We conclude that the SM82 approximation accurately represents quasi-2-D flows at moderately large 
$Gr$ (
$10^8$ and 
$10^9$ in our system). The accuracy deteriorates at higher 
$Gr$ even though the flow remains quasi-two dimensional. An example of this is observed at 
$Ha \gtrsim 4000$ and 
$Gr = 10^{10}$. The base flow profiles are clearly different between the 2-D and 3-D models (see figure 5
$c$,
$f$). The 2-D 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 6
$c$). The values of the growth rate 
$\gamma$ are, however, substantially underestimated by the 2-D model (see figure 6
$c$ and table 4). The nonlinear flow states resulting from the instability are predicted incorrectly by the 2-D 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 2-D and 3-D models.
We do not have a satisfactory explanation for this effect and leave its further exploration for future studies. It should be mentioned that the quality of the 2-D approximation improves with increasing 
$Ha$. As an example, the values of the linear instability growth rate 
$\gamma$ shown in table 4 are underpredicted by the 2-D model by about 
$35\,\%$ at 
$Ha = 4000$ and by about only 
$18\,\%$ at 
$Ha = 10\,000$.
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 3-D structures are present in the flow for which the diffusion length is shorter than the size of the domain (Pothérat & Klein Reference Pothérat and Klein2014, Reference Pothérat and Klein2017).
In the remaining part of this paper, the discussion of the instability and nonlinear states is based on the SM82 2-D approximation for flows at 
$Gr = 10^8$ and 
$10^9$ and on the full 3-D solutions for flows at 
$10^{10}$.
The numerical solution for the 2-D 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.7a,b) for the base flow.
3.3. 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 
$\lambda$ is varied with step 
$0.1$ in the simulations. The computed values of the exponential growth rate 
$\gamma$ as a function of the wavelength 
$\lambda$ for various combinations of 
$Ha$ and 
$Gr$ are shown in figure 10. Two trends of the linear stability behaviour were proposed by Zhang & Zikanov (Reference Zhang and Zikanov2014): (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 (Reference Zhang and Zikanov2014) to modification of the base flow, namely to suppression of the transverse circulation resulting in weaker mixing and stronger unstable temperature stratification.
Figure 10. Rates of exponential growth 
$\gamma$ shown as functions of axial wavelength 
$\lambda$ at 
$Gr = 10^8$ (
$a$), 
$Gr = 10^9$ (
$b$) and 
$Gr = 10^{10}$ (
$c$) for various values of 
$Ha$. Results of 2-D SM82 approximation are shown for 
$Gr = 10^8$ and 
$10^9$. Results of 3-D computational analysis are shown for 
$Gr = 10^{10}$. Additional plots in (
$b$) and (
$c$) show zoomed-in areas around (
$\lambda _{max}, \gamma _{max}$).
Table 5. Results of linear stability analysis. Wavelengths 
$\lambda _{max}$ and exponential growth rates 
$\gamma _{max}$ of the fastest growing modes are shown as functions of 
$Ha$ and 
$Gr$. Only the data for flow regimes identified as quasi-two dimensional in the analysis of the base flow are shown.
In order to investigate these trends for our system and parameter range, we present the exponential growth rate 
$\gamma _{max}$ of the fastest growing modes and the corresponding wavelengths 
$\lambda _{max}$ in figure 11 and table 5. Our results are clearly consistent with the first trend, 
$\gamma _{max}$ increases and 
$\lambda _{max}$ decreases with growing 
$Gr$. However, we find a different behaviour in regard of the second trend. The increase of 
$Ha$ leads to a noticeable or slight decrease of the exponential growth rates at 
$Gr = 10^8$ or 
$Gr = 10^9$, respectively (see figure 10
$c$). For 
$Gr = 10^{10}$, 
$\gamma$ is nearly insensitive to the values of 
$Ha$.
Figure 11. The wavelength 
$\lambda _{max}$ (
$a$) and the exponential growth rate 
$\gamma _{max}$ (
$b$) of the fastest growing perturbations as a function of 
$Ha$. The data are taken from table 5. The results of the current work and that of Zhang & Zikanov (Reference Zhang and Zikanov2014) are denoted as filled and unfilled elements, respectively. Results of 2-D SM82 approximation are shown for 
$Gr = 10^8$ and 
$10^9$. Results of 3-D computational analysis are shown for 
$Gr = 10^{10}$.
We conclude that the counterintuitive behaviour of stronger instability at higher 
$Ha$ is not observed in quasi-2-D flows at high 
$Ha$ considered in our study. It is, nevertheless, noteworthy that 
$\gamma _{max}$ does not decrease with 
$Ha$ at 
$Gr = 10^{10}$. It increases slightly at 
$Ha \lesssim 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 5(
$c$). Here 
$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. (Reference Zikanov, Listratov and Sviridov2013); Zhang & Zikanov (Reference Zhang and Zikanov2014). 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. (Reference Zikanov, Listratov and Sviridov2013) and Zhang & Zikanov (Reference Zhang and Zikanov2014), 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 \geq 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 quasi-2-D instabilities by Zhang & Zikanov (Reference Zhang and Zikanov2014). 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$.
Figure 12. Spatial structure of the fastest growing instability modes during the stage of exponential growth at 
$Gr = 10^8$ (
$a$), 
$Gr = 10^9$ (
$b$) and 
$Gr = 10^{10}$ (
$c$,
$d$) for 
$Ha = 10\,000$. Results of 2-D SM82 approximation are shown for 
$Gr = 10^8$ and 
$10^9$. Results of 3-D computational analysis are shown for 
$Gr = 10^{10}$. Perturbations of temperature and vector fields of velocity perturbations (
$u'$, 
$w'$) in the vertical midplane (
$y=0$) are shown in (
$a$–
$c$). Perturbations of vertical velocity in the horizontal midplane section 
$z=0$ are shown in (
$d$). Solid and dashed isolines indicate positive and negative values, respectively.
3.4. 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. A DNS approach based on direct solution of the nonlinear governing equations is utilized. The 2-D SM82 model (2.14)–(2.16) and the computational domain of length 
$L_x = 4{\rm \pi}$ are used for flows at 
$Gr=10^8$ and 
$10^9$. Full 3-D equations (2.1)–(2.3) are solved and the domain is reduced to 
$L_x = 2{\rm \pi}$ at 
$Gr=10^{10}$. Other parameters of the computational model are described in § 2.5. Each simulation starts with the streamwise-independent base flow (see § 3.2) computed at the same 
$Gr$ and 
$Ha$, to which random small-amplitude (
${\sim }10^{-3}$) random perturbations of velocity and temperature are added.
Figure 13. Time signals of the kinetic energy of the streamwise velocity obtained in the DNS of flows at 
$Gr = 10^8$ (
$a$), 
$Gr = 10^9$ (
$b$) and 
$Gr = 10^{10}$ (
$c$) for 
$Ha = 5000$ and 
$Ha = 10\,000$. Results of 2-D SM82 approximation are shown for 
$Gr = 10^8$ and 
$10^9$. Results of 3-D computational analysis are shown for 
$Gr = 10^{10}$.
Figure 14. Flow structure in nonlinear regime at 
$Gr = 10^8$ (
$a$–
$c$) and 
$Gr = 10^9$ (
$d$–
$f$) for 
$Ha = 10\,000$. The instantaneous distributions of 
$u_z$ (
$a$,
$d$), 
$u_x$ (
$b$,
$e$) and 
$\theta$ (
$c$,
$f$) obtained in the 2-D model are shown. The profiles of 
$u_x$ and 
$\theta$ obtained by averaging over 
$x$ and time are shown, respectively, in (
$g$) and (
$h$).
Figure 15. Time signals of temperature measured at top and bottom walls and in the middle of the duct in fully developed flows at 
$Gr = 10^8$, 
$10^9$, 
$10^{10}$ are shown for 
$Ha = 5 \times 10^3$ in (
$a$,
$c$,
$e$) and for 
$Ha = 10^4$ in (
$b$,
$d$,
$f$). Results of 2-D SM82 approximation are shown for 
$Gr = 10^8$ and 
$10^9$ (
$a$–
$d$). Results of 3-D computational analysis are shown for 
$Gr = 10^{10}$ (
$e$–
$f$).
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 2-D model at 
$Gr = 10^8, 10^9$ and at least 
$100$ time units for the 3-D 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}$.
The 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 8
$d$ and 14
$a$,
$d$). The rolls cause variations of temperature (see figures 8
$f$, 14
$c$ and 14
$f$). Transport of the rolls by the mean flow is a known reason of MCFs in horizontal channels (Zikanov et al. Reference Zikanov, Listratov and Sviridov2013, Reference Zikanov, Listratov, Razuvanov, Belyaev, Frick and Sviridov2021; Zhang & Zikanov Reference Zhang and Zikanov2014).
Comparison of the flow structures in figures 8(
$d$–
$f$) and 14 reveal the effect of the value of 
$Gr$ on convection rolls. As anticipated, an increase of 
$Gr$ leads to a higher non-dimensional amplitude of the velocity fluctuations. This results in stronger vertical mixing as illustrated by the streamwise- and time-averaged profiles in figure 14(
$g$,
$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 8(
$f$), 14(
$c$) and 14(
$f$), 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. (Reference Zikanov, Listratov, Razuvanov, Belyaev, Frick and Sviridov2021), 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 \times 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 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 the 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 that the duct half-width 
$d = 5$ cm and using the physical properties of PbLi at 
$573$ K (Zikanov et al. Reference Zikanov, Listratov, Razuvanov, Belyaev, Frick and Sviridov2021). 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/\kappa$, 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. Reference Belyaev, Sviridov, Batenin, Biryukov, Nikitina, Manchkha, Pyatnitskaya, Razuvanov and Sviridov2017), 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. Reference Zikanov, Listratov, Razuvanov, Belyaev, Frick and Sviridov2021). 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$.
4. Concluding remarks
We have analysed mixed convection in a liquid metal flow in a duct with bottom heating and a 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 (Reference Zhang and Zikanov2014).
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 2-D approximation by Sommeria & Moreau (Reference Sommeria and Moreau1982) 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 2-D 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 they are obtained for a single configuration of a horizontal duct flow with bottom heating and a transverse magnetic field. At the same time, there are multiple indications that similar behaviours can be observed in other configurations related to the existing designs of liquid metal blankets of fusion reactors. This will need to be explored in future studies.
Further study of MCFs at high Hartmann and Grashof numbers is warranted by their practical importance and theoretical significance. Many interesting possible directions of future work can be suggested. We mention two of them. One is the exploration of the phenomenon for other geometries, where strong MCFs are known to exist, for example, for downward flow in a vertical duct. Another particularly interesting direction is the analysis of the effects of finite thermal and electrical conductivities of the walls.
Acknowledgements
The authors are thankful to D. Krasnov for continuing assistance with the numerical model and S. Molokov for interesting and useful discussions.
Funding
The work of R.A. and O.Z. is supported by the US NSF (grant CBET 1803730 ‘Extreme magnetoconvection’). The work of Y.L. is supported by the Ministry of Science and Education of the Russian Federation (grant 14.Z50.31.0042) and by the Russian Foundation for Basic Research (grant NNIO 18-508-12005).
Declaration of interests
The authors report no conflict of interest.
Open access
