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

(2.1)\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}

(2.2)\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0, \end{gather}

(2.3)\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

(2.4)\begin{equation} T(\boldsymbol{x},t) = \it {T_m}(x) + \theta(\boldsymbol{x},t) \end{equation}

of fluctuations $\theta$ and the mean-mixed temperature

(2.5)\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.

(2.6)\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

(2.7)\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

(2.8)\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$,

(2.9)\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

(2.10)\begin{equation} \boldsymbol{F }_b = Gr Re^{{-}2} \boldsymbol{e}_z T. \end{equation}

The Lorentz force is computed as

(2.11)\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

(2.12)\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,

(2.13)\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

(2.14)\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}

(2.15)\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0, \end{gather}

(2.16)\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

(2.17)\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

(2.18)\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

(2.19a,b)\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.