2.1 Quasi-two-dimensional model

The core flow dynamics can be suitably approximated through the use of the quasi-two-dimensional MHD model proposed by Sommeria & Moreau (Reference Sommeria and Moreau1982) in the limit of high interaction parameter ( $N\equiv \mathit{Ha}^{2}/Re\gg 1$ ) and Hartman number ( $\mathit{Ha}\gg 1$ ). Electrically insulated side walls are assumed in the construction of the model, which leads to a two-dimensional Navier–Stokes equation augmented via a linear braking term representing friction due to the Hartmann layers. By employing the subscript $\bot$ to represent the horizontal components of the averaging of velocity and pressure in the $y$ -direction, in addition to the horizontal trace of operators, the quasi-two-dimensional MHD equations are given as

(2.6a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}_{\bot }}{\unicode[STIX]{x2202}t}=-(\boldsymbol{u}_{\bot }\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\bot })\boldsymbol{u}_{\bot }-\unicode[STIX]{x1D735}_{\bot }p_{\bot }+\frac{1}{Re}\unicode[STIX]{x1D6FB}_{\bot }^{2}\boldsymbol{u}_{\bot }-\frac{H}{Re}\boldsymbol{u}_{\bot }, & \displaystyle\end{eqnarray}$$

(2.6b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{\bot }=0, & \displaystyle\end{eqnarray}$$

$\boldsymbol{u}_{\bot }(x,z,t)$

$H\equiv n\mathit{Ha}(a/b)^{2}$

$n$

$n=2$

$b$

(2.7) $$\begin{eqnarray}\boldsymbol{u}_{\bot }=0\quad \text{at }x=\pm a.\end{eqnarray}$$

2.2 Base flow and linear perturbations

2.2.1 Three-dimensional flows

The three-dimensional MHD analysis employs a base flow that is steady, streamwise independent $\boldsymbol{u}_{0}(x,y)=w_{0}(x,y)$ and which conforms to the boundary conditions (2.4)–(2.5). The solution is obtained through time integration of the MHD governing equations (2.1) to a time-independent steady-state flow. The general solutions for the velocity, electric-potential and pressure fields can be represented by the summation of the base flow (denoted by the subscript $0$ ) and perturbation equations (denoted by the superscript $\prime$ ) such that

(2.8) $$\begin{eqnarray}[\boldsymbol{u},\unicode[STIX]{x1D719},p]=[\boldsymbol{u}_{0},\unicode[STIX]{x1D719}_{0},p_{0}]+[\boldsymbol{u}^{\prime },\unicode[STIX]{x1D719}^{\prime },p^{\prime }].\end{eqnarray}$$

The linearised governing MHD equations are constructed by substituting (2.8) into (2.1) and (2.3), whilst retaining only the lowest-order linear terms. The resulting equations are

(2.9a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}^{\prime }}{\unicode[STIX]{x2202}t}=-\text{DN}^{\prime }\boldsymbol{u}^{\prime }-\unicode[STIX]{x1D735}p^{\prime }+\frac{1}{Re}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}^{\prime }+\frac{Ha^{2}}{Re}[(-\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}^{\prime }+\boldsymbol{u}^{\prime }\times \boldsymbol{e}_{y})\times \boldsymbol{e}_{y}], & \displaystyle\end{eqnarray}$$

(2.9b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}^{\prime }=0, & \displaystyle\end{eqnarray}$$

(2.9c ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D719}^{\prime }=\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{u}^{\prime }\times \boldsymbol{e}_{y}), & \displaystyle\end{eqnarray}$$

where $\text{DN}^{\prime }\boldsymbol{u}^{\prime }=(\unicode[STIX]{x1D735}\boldsymbol{u}_{0})\boldsymbol{\cdot }\boldsymbol{u}^{\prime }+(\boldsymbol{u}_{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}^{\prime }$ is the linearised advection operator.

The proceeding transient growth analysis requires the construction of the adjoint form of the linearised governing MHD equations. Through an analogous derivation to that presented in Blackburn, Sherwin & Barkley (Reference Blackburn, Sherwin and Barkley2008), the adjoint reformulation of (2.9) is given as

(2.10a ) $$\begin{eqnarray}\displaystyle & \displaystyle -\frac{\unicode[STIX]{x2202}\boldsymbol{u}^{\ast }}{\unicode[STIX]{x2202}t}=-\text{DN}^{\ast }\boldsymbol{u}^{\ast }-\unicode[STIX]{x1D735}p^{\ast }+\frac{1}{Re}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}^{\ast }+\frac{Ha^{2}}{Re}[(-\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}^{\ast }+\boldsymbol{u}^{\ast }\times \boldsymbol{e}_{y})\times \boldsymbol{e}_{y}], & \displaystyle\end{eqnarray}$$

(2.10b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}^{\ast }=0, & \displaystyle\end{eqnarray}$$

(2.10c ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D719}^{\ast }=\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{u}^{\ast }\times \boldsymbol{e}_{y}), & \displaystyle\end{eqnarray}$$

where $\text{DN}^{\ast }\boldsymbol{u}^{\ast }=(\unicode[STIX]{x1D735}\boldsymbol{u}_{0})^{\text{T}}\boldsymbol{\cdot }\boldsymbol{u}^{\ast }-(\boldsymbol{u}_{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}^{\ast }$ is the adjoint linear advection operator. $\boldsymbol{u}^{\ast }$ , $p^{\ast }$ and $\unicode[STIX]{x1D719}^{\ast }$ denote the velocity, pressure and electro-potential fields of the adjoint perturbation scenario, respectively.

Spatial homogeneity in the streamwise direction is assumed so that the time dependent three-dimensional linear perturbations for both the forward and adjoint systems are considered in the form of the respective decoupled orthogonal Fourier modes

(2.11) $$\begin{eqnarray}\displaystyle & \displaystyle [\boldsymbol{u}^{\prime },\unicode[STIX]{x1D719}^{\prime },p^{\prime }]=[\hat{u} ,\hat{v},{\hat{w}},\hat{\unicode[STIX]{x1D719}},\hat{p}](x,y,t)\boldsymbol{\cdot }\text{e}^{\text{i}kz}, & \displaystyle\end{eqnarray}$$

(2.12) $$\begin{eqnarray}\displaystyle & \displaystyle [\boldsymbol{u}^{\ast },\unicode[STIX]{x1D719}^{\ast },p^{\ast }]=[\hat{u} ^{\ast },\hat{v}^{\ast },{\hat{w}}^{\ast },\hat{\unicode[STIX]{x1D719}}^{\ast },\hat{p}^{\ast }](x,y,t)\boldsymbol{\cdot }\text{e}^{\text{i}kz}, & \displaystyle\end{eqnarray}$$

where $k=2\unicode[STIX]{x03C0}/\ell _{z}$ is the associated streamwise wavenumber in the $z$ -direction (a parameter to be varied in this study), and $\ell _{z}$ is the periodicity length of the domain, also in the $z$ -direction. The circumflex ( $\hat{\boldsymbol{\cdot }}$ ) represents the Fourier coefficient of the associated base variable.

2.2.2 Quasi-two-dimensional flows

Unlike the three-dimensional transient growth analysis, the method employing the SM82 model does not utilise Fourier modal expansion in the third dimension, but instead discretises the domain using spectral elements. For this case, the analytical solution for a one-dimensional steady Q2D base flow

(2.13) $$\begin{eqnarray}\boldsymbol{u}_{\bot ,0}(x)=\frac{\cosh \sqrt{H}}{\cosh \sqrt{H}-1}\left(1-\frac{\cosh (\sqrt{H}x)}{\cosh \sqrt{H}}\right),\end{eqnarray}$$

as defined in Pothérat (Reference Pothérat2007) is employed. Subjecting (2.13) to two-dimensional time-dependent linear perturbations of the form

(2.14) $$\begin{eqnarray}[\boldsymbol{u}_{\bot }^{\prime },p_{\bot }^{\prime }]=[u^{\prime },w^{\prime },p^{\prime }](x,z,t),\end{eqnarray}$$

the flow can then be described by the general solution

(2.15) $$\begin{eqnarray}[\boldsymbol{u}_{\bot },p_{\bot }]=[\boldsymbol{u}_{\bot ,0},p_{\bot ,0}]+[\boldsymbol{u}_{\bot }^{\prime },p_{\bot }^{\prime }].\end{eqnarray}$$

Hence, through substitution of (2.15) into (2.6), while keeping only the lowest-order linear terms, the linearised Q2D SM82 equations are given as

(2.16a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}_{\bot }^{\prime }}{\unicode[STIX]{x2202}t}=-\text{DN}\boldsymbol{u}_{\bot }^{\prime }-\unicode[STIX]{x1D735}_{\bot }p_{\bot }^{\prime }+\frac{1}{Re}\unicode[STIX]{x1D6FB}_{\bot }^{2}\boldsymbol{u}_{\bot }^{\prime }+\frac{H}{Re}\boldsymbol{u}_{\bot }^{\prime }, & \displaystyle\end{eqnarray}$$

(2.16b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}_{\bot }\boldsymbol{\cdot }\boldsymbol{u}_{\bot }^{\prime }=0, & \displaystyle\end{eqnarray}$$

where $\text{DN}^{\prime }\boldsymbol{u}_{\bot }^{\prime }=(\unicode[STIX]{x1D735}_{\bot }\boldsymbol{u}_{\bot ,0})\boldsymbol{\cdot }\boldsymbol{u}_{\bot }^{\prime }+(\boldsymbol{u}_{\bot ,0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\bot })\boldsymbol{u}_{\bot }^{\prime }$ is the Q2D linearised advection term. Similarly, the adjoint linearised SM82 equations are given by

(2.17a ) $$\begin{eqnarray}\displaystyle & \displaystyle -\frac{\unicode[STIX]{x2202}\boldsymbol{u}_{\bot }^{\ast }}{\unicode[STIX]{x2202}t}=-\text{DN}^{\ast }\boldsymbol{u}_{\bot }^{\ast }-\unicode[STIX]{x1D735}_{\bot }p_{\bot }^{\ast }+\frac{1}{Re}\unicode[STIX]{x1D6FB}_{\bot }^{2}\boldsymbol{u}_{\bot }^{\ast }+\frac{H}{Re}\boldsymbol{u}_{\bot }^{\ast }, & \displaystyle\end{eqnarray}$$

(2.17b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}_{\bot }\boldsymbol{\cdot }\boldsymbol{u}_{\bot }^{\ast }=0, & \displaystyle\end{eqnarray}$$

where $\text{DN}^{\ast }\boldsymbol{u}_{\bot }^{\ast }=(\unicode[STIX]{x1D735}_{\bot }\boldsymbol{u}_{\bot ,0})^{\text{T}}\boldsymbol{\cdot }\boldsymbol{u}_{\bot }^{\ast }-(\boldsymbol{u}_{\bot ,0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\bot })\boldsymbol{u}_{\bot }^{\ast }$ is the adjoint Q2D linearised advection term.

The streamwise periodicity length of the Q2D domain $\ell _{z,\bot }$ is determined from the three-dimensional transient growth optimal wavenumber $k_{opt}$ (as defined in § 2.3) at an equivalent $\mathit{Ha}$ such that $\ell _{z,\bot }=2\unicode[STIX]{x03C0}/k_{opt}$ .

2.3 Transient growth analysis

The present work is interested in the transient energy growth of linear perturbations over a finite time interval $\unicode[STIX]{x1D70F}$ (a parameter to be varied in this study). Outlined in this section is the methodology relating to the three-dimensional MHD case only. The formulation pertaining to the SM82 model is analogous to that presented next, and is therefore omitted for brevity.

The governing direct and adjoint three-dimensional MHD linearised equations, given respectively in (2.9) and (2.10), are solved subject to the boundary conditions (2.4)–(2.5) in combination with suitable initial conditions. The nature of linear transient growth is determined through a direct–adjoint eigenvalue system which is solved using a methodology consistent with Barkley, Blackburn & Sherwin (Reference Barkley, Blackburn and Sherwin2008). In this direct–adjoint system, an initial perturbation field providing a given gain in kinetic energy is found through the construction of an auxiliary eigenvalue problem

(2.18a,b ) $$\begin{eqnarray}\mathscr{A}^{\ast }(\unicode[STIX]{x1D70F})\mathscr{A}(\unicode[STIX]{x1D70F})\hat{\boldsymbol{u}}_{k}=\unicode[STIX]{x1D706}_{k}\hat{\boldsymbol{u}}_{k},\quad \Vert \hat{\boldsymbol{u}}_{k}\Vert =1,\end{eqnarray}$$

where $\unicode[STIX]{x1D706}_{k}$ and $\hat{\boldsymbol{u}}_{k}$ denote eigenvalues and normalised eigenvectors, respectively. The state-transition operator $\mathscr{A}(\unicode[STIX]{x1D70F})$ represents the terms on the right-hand side of (2.9a ). Therefore, the operator describes the time evolution of an arbitrary initial perturbation $\boldsymbol{u}^{\prime }(0)$ to $t=\unicode[STIX]{x1D70F}$ such that $\boldsymbol{u}^{\prime }(\unicode[STIX]{x1D70F})=\mathscr{A}(\unicode[STIX]{x1D70F})\boldsymbol{u}^{\prime }(0)$ . Similarly, $\mathscr{A}^{\ast }(\unicode[STIX]{x1D70F})$ is the equivalent adjoint evolution operator of $\mathscr{A}(\unicode[STIX]{x1D70F})$ that evolves an equivalent adjoint variable $\boldsymbol{u}^{\ast }(\unicode[STIX]{x1D70F})$ , as solved via (2.10a ), backwards in time from $t=\unicode[STIX]{x1D70F}$ to $t=0$ . The eigenvectors $\hat{\boldsymbol{u}}_{k}$ , which represent the initial perturbation field $\boldsymbol{u}^{\prime }(0)$ that generates growth $\unicode[STIX]{x1D706}_{k}$ over time $\unicode[STIX]{x1D70F}$ , are obtained by determining the equivalent real and orthonormal right singular vectors produced by a singular value decomposition of $\mathscr{A}^{\ast }(\unicode[STIX]{x1D70F})\mathscr{A}(\unicode[STIX]{x1D70F})$ . The singular value decomposition is determined using an implicitly restarted Arnoldi iterative method, the particulars of which can be found in Lehoucq, Sorensen & Yang (Reference Lehoucq, Sorensen and Yang1998). It is important to note that the operator $\mathscr{A}^{\ast }(\unicode[STIX]{x1D70F})\mathscr{A}(\unicode[STIX]{x1D70F})$ is not explicitly constructed: instead, the action of $\mathscr{A}^{\ast }(\unicode[STIX]{x1D70F})\mathscr{A}(\unicode[STIX]{x1D70F})$ on the field $\boldsymbol{u}^{\prime }$ is determined through forward time integration of an initial condition over a time interval $\unicode[STIX]{x1D70F}$ via (2.9), and then subsequently backwards from $\unicode[STIX]{x1D70F}$ to the initial time under the adjoint system in (2.10).

For the present transient growth analysis, it is the eigenvector producing the maximal growth for all $k$ and $\unicode[STIX]{x1D70F}$ that is of particular interest. This eigenvector is commonly referred to as the optimal mode and this terminology is thus used hereafter. The optimal energy gain $G_{max}$ occurs at an optimal time interval $\unicode[STIX]{x1D70F}_{opt}$ having optimal streamwise wavenumber $k_{opt}$ .

3.1 Spatial discretisation

A high-order spectral element method following that described in Karniadakis, Israeli & Orszag (Reference Karniadakis, Israeli and Orszag1991) is employed for spatial discretisation of the governing equations. The domain is meshed using quadrilateral macro-elements with internally applied Lagrangian polynomial functions, the order of which, $N_{p}$ , is varied to control spatial resolution in the spanwise $x$ – $y$ cross-plane for three-dimensional investigations, and in the streamwise $x$ – $z$ plane for the Q2D case. Linear grading of macro-elements towards all solid surfaces is employed to help resolve regions that experience the largest flow gradients. The spacing of spectral elements is scaled with respect to the $\mathit{Ha}$ -dependent boundary layer thickness, ensuring that a minimum of 8 spectral elements span its height. For the Q2D analysis, a minimum of 20 macro-elements span the streamwise and spanwise directions. An example of the macro-element distribution in the $x$ – $y$ cross-plane for three-dimensional investigations at $\mathit{Ha}=30$ is shown in figure 1.

3.2 Temporal discretisation

A third-order backwards differentiation scheme involving an operator-splitting method is employed for temporal advancement of both three-dimensional and quasi-two-dimensional governing equations. This is analogous to the method described in Karniadakis et al. (Reference Karniadakis, Israeli and Orszag1991). To the best of the authors’ knowledge, this temporal discretisation scheme has not been previously implemented for three-dimensional MHD flows. The implicit solution of the electric-potential field at each time step based on the high-order projection of the velocity field, circumvents common issues pertaining to charge conservation encountered in finite-volume-based schemes (Ni et al. Reference Ni, Munipalli, Morley, Huang and Abdou2007). As such, presented here is a full outline for the nonlinear three-dimensional governing MHD equations in (2.1). The three-dimensional linearised equations (2.9) and (2.10), as well as the corresponding SM82 equations in (2.16) and (2.17) are solved analogously.

The MHD governing equations (2.1) are evaluated at future time step ( $r+1$ ), with the time derivative term approximated via backwards differentiation. The semi-discrete form of equation (2.1) is

(3.1) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\displaystyle \unicode[STIX]{x1D6FC}_{0}\boldsymbol{u}^{(r+1)}-\mathop{\sum }_{q=1}^{J}\unicode[STIX]{x1D6FC}_{q}\boldsymbol{u}^{(r-q+1)}}{\unicode[STIX]{x0394}t}\nonumber\\ \displaystyle & & \displaystyle \quad =\left[-(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}-\unicode[STIX]{x1D735}p+\frac{1}{Re}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}+\frac{Ha^{2}}{Re}[(-\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}+\boldsymbol{u}\times \boldsymbol{e}_{y})\times \boldsymbol{e}_{y}]\right]^{(r+1)},\end{eqnarray}$$

where $J$ denotes the order of the scheme (e.g.  $J=3$ in the present formulation) and $\unicode[STIX]{x1D6FC}_{q}$ are the corresponding coefficients (i.e.  $\unicode[STIX]{x1D6FC}_{0}=11/6$ , $\unicode[STIX]{x1D6FC}_{1}=-3$ , $\unicode[STIX]{x1D6FC}_{2}=3/2$ and $\unicode[STIX]{x1D6FC}_{3}=-1/3$ ).

An operator-splitting approach is used to split the resolution of equation (3.1) into three sub-steps. As a precursor, an explicit projection of the velocity field to the future ( $r+1$ ) time is evaluated,

(3.2) $$\begin{eqnarray}\widetilde{\boldsymbol{u}}=\mathop{\sum }_{q=0}^{J-1}\unicode[STIX]{x1D6FE}_{q}\boldsymbol{u}^{(r-q)},\end{eqnarray}$$

where for $J=3$ the $\unicode[STIX]{x1D6FE}_{q}$ coefficients are $\unicode[STIX]{x1D6FE}_{0}=3$ , $\unicode[STIX]{x1D6FE}_{1}=-3$ , $\unicode[STIX]{x1D6FE}_{2}=1$ . Introducing intermediate velocity fields $\boldsymbol{u}^{\dagger }$ and $\boldsymbol{u}^{\ddagger }$ , the semi-discrete projection for the advection, pressure and viscous diffusion terms are then respectively given by

(3.3a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\displaystyle \boldsymbol{u}^{\dagger }-\mathop{\sum }_{q=1}^{J}\unicode[STIX]{x1D6FC}_{q}\boldsymbol{u}^{(r-q+1)}}{\unicode[STIX]{x0394}t}=-(\widetilde{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\widetilde{\boldsymbol{u}}+\frac{Ha^{2}}{Re}[(-\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}^{(r+1)}+\widetilde{\boldsymbol{u}}\times \boldsymbol{e}_{y})\times \boldsymbol{e}_{y}], & \displaystyle\end{eqnarray}$$

(3.3b ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\boldsymbol{u}^{\ddagger }-\boldsymbol{u}^{\dagger }}{\unicode[STIX]{x0394}t}=-\unicode[STIX]{x1D735}p^{(r+1)}, & \displaystyle\end{eqnarray}$$

(3.3c ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x1D6FC}_{0}\boldsymbol{u}^{(r+1)}-\boldsymbol{u}^{\ddagger }}{\unicode[STIX]{x0394}t}=\frac{1}{Re}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}^{(r+1)}. & \displaystyle\end{eqnarray}$$

(3.4) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D719}^{(r+1)}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\widetilde{\boldsymbol{u}}\times \boldsymbol{e}_{y}).\end{eqnarray}$$

Taking the divergence of equation (3.3b ) and enforcing the divergence-free constraint on the second intermediate velocity field (i.e.  $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}^{\ddagger }=0$ ) yields a Poisson equation for the pressure,

(3.5) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}^{2}p^{(r+1)}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left(\frac{\boldsymbol{u}^{\dagger }}{\unicode[STIX]{x0394}t}\right).\end{eqnarray}$$

Equation (3.3c ) is recast as a Helmholtz equation,

(3.6) $$\begin{eqnarray}\left(\unicode[STIX]{x1D6FB}^{2}-\unicode[STIX]{x1D6FC}_{0}\frac{Re}{\unicode[STIX]{x0394}t}\right)\boldsymbol{u}^{(r+1)}=-\frac{Re}{\unicode[STIX]{x0394}t}\boldsymbol{u}^{\ddagger },\end{eqnarray}$$

for solution of the final velocity field $\boldsymbol{u}^{(r+1)}$ .

Dirichlet and Neumann boundary conditions on $\unicode[STIX]{x1D719}$ , $p$ and $\boldsymbol{u}$ are respectively imposed through equations (3.4)–(3.6). In addition to the boundary conditions defined in (2.4)–(2.5), the above equations are solved subject to high-order pressure field Neumann-type boundary conditions enforced on all domain perimeters. This allows for third-order time accuracy to be preserved (Karniadakis et al. Reference Karniadakis, Israeli and Orszag1991).

The present numerical scheme is built upon an existing solver that has been rigorously validated for two- and three-dimensional incompressible Navier–Stokes flows and Q2D MHD flows (Hussam, Thompson & Sheard Reference Hussam, Thompson and Sheard2012a ,Reference Hussam, Thompson and Sheard b ; Hussam & Sheard Reference Hussam and Sheard2013; Vo et al. Reference Vo, Montabone, Read and Sheard2015; Cassells et al. Reference Cassells, Hussam and Sheard2016; Hamid, Hussam & Sheard Reference Hamid, Hussam and Sheard2016; Leigh, Tsai & Sheard Reference Leigh, Tsai and Sheard2016; Ng et al. Reference Ng, Vo, Hussam and Sheard2016; Sheard, Hussam & Tsai Reference Sheard, Hussam and Tsai2016). Nonetheless, a careful eye was kept on the conservation properties of the numerical algorithm, specifically, the discretised divergence of the velocity and current fields. The domain-integrated norms for both values did not exceed $10^{-9}$ in any of the simulations presented herein.

3.3 Numerical validation

To ensure an accurate computation of the three-dimensional and quasi-two-dimensional steady MHD base flows, the velocity profiles were validated against analytical solutions as developed by Hunt & Stewartson (Reference Hunt and Stewartson1965) and Pothérat (Reference Pothérat2007), respectively. The integrated error between the numerical and analytical velocity field solutions was found to be below 0.001 % of the peak velocity for all cases presented here. Furthermore, convergence of wall shear stresses $\unicode[STIX]{x1D70F}_{zx}$ and $\unicode[STIX]{x1D70F}_{zy}$ , the standard ${\mathcal{L}}^{2}$ norm and the current magnitude $|\,j|$ of better than 0.03 % (where 0 % convergence is exact) was satisfied for polynomial orders of $N_{p}=8$ , which is thus employed hereafter. It is worth noting here that adequate convergence was obtained for domain-integrated measures, such as ${\mathcal{L}}^{2}$ and $|\,j|$ , at much lower polynomial orders ( $N_{p}\leqslant 5$ ). However, accurate resolution of the Hartmann and Shercliff boundary layer dynamics, as measured through wall shear stresses, becomes the most stringent obstacle to numerical precision due to the large velocity gradients of these layers. Adequate grid independence was also seen for the quasi-two-dimensional model at $N_{p}=8$ , and this component of the numerical scheme has also been previously validated in works by Hamid et al. (Reference Hamid, Hussam, Pothérat and Sheard2015) and Cassells et al. (Reference Cassells, Hussam and Sheard2016).

Table 1.

Grid sensitivity data for three-dimensional MHD base flow computations at $Re=15\,000$ for two Hartmann number regimes: $\mathit{Ha}=30$ and $\mathit{Ha}=800$ (shown in the parentheses). The parameter $N_{p}$ indicates the polynomial order of the macro-element shape function. The wall shear stress inside a single Hartmann or Shercliff layer is indicated respectively by $\unicode[STIX]{x1D70F}_{zx}$ and $\unicode[STIX]{x1D70F}_{zy}$ . The integrals of velocity magnitude ${\mathcal{L}}^{2}$ and current density $|\,j|$ over the flow domain are also given.

The accuracy of eigenmode predictions are dependent on the precision of the Krylov subspace iteration employed in the Arnoldi method, in addition to the resolution of the numerical methods. For the former, the iterative precision of the eigenvalue $\unicode[STIX]{x1D706}$ and eigenvector $\hat{\boldsymbol{u}}$ computations was measured by the residual

(3.7) $$\begin{eqnarray}r=\Vert \mathscr{A}^{\ast }\mathscr{A}\hat{\boldsymbol{u}}-\unicode[STIX]{x1D706}\hat{\boldsymbol{u}}\Vert ,\end{eqnarray}$$

where $\Vert \boldsymbol{\cdot }\Vert$ is the standard vector norm. Eigenvalue convergence was deemed sufficient when $r<10^{-7}$ . For the latter, the grid independence of eigenvalue computations were ensured through a sufficiently high-order polynomial degree shape function $N_{p}$ . For Hartmann numbers $\mathit{Ha}=30$ and $800$ at $Re=15\,000$ , the optimal eigenvalue predictions as a function of $N_{p}$ for both three-dimensional and quasi-two-dimensional models are reported in table 2. For both Hartmann numbers, eigenvalue convergence of better than $0.002\,\%$ was achieved at $N_{p}=8$ , which is employed hereafter for all transient growth analyses.

Table 2.

Grid sensitivity data for optimal eigenvalue computations using both three-dimensional and quasi-two-dimensional (shown in parentheses) MHD models at $Re=15\,000$ for two Hartmann number regimes: $\mathit{Ha}=30$ and $\mathit{Ha}=800$ . The parameter $N_{p}$ indicates the polynomial order of the macro-element shape function.

For the three-dimensional eigenvalue computations, streamwise wavenumbers were investigated in the range $0\leqslant k\leqslant 80$ , with the local maxima resolved to an accuracy of at least $0.1\,\%$ of the peak value. To ensure that the optimal growths were captured sufficiently, and a monotonic decay in amplifications were achieved at higher $\unicode[STIX]{x1D70F}$ , the analysis was conducted over time intervals extending to $\unicode[STIX]{x1D70F}=5\unicode[STIX]{x1D70F}_{opt}$ for both three-dimensional and quasi-two-dimensional models. The optimal time intervals were resolved to at least 0.2 time units. The chosen numerical methodology and resolution were able to closely capture the same optimal growth rates predicted for hydrodynamic flow ( $\mathit{Ha}=0$ ), as well as the scaling in the range of $10\leqslant \mathit{Ha}\leqslant 50$ at $Re=5000$ obtained by Krasnov et al. (Reference Krasnov, Zikanov, Rossi and Boeck2010). The linearised component of this solver has also been previously verified and implemented in works such as Hussam et al. (Reference Hussam, Thompson and Sheard2012b ), Tsai et al. (Reference Tsai, Hussam, Fouras and Sheard2016) and Sapardi et al. (Reference Sapardi, Hussam, Pothérat and Sheard2017).

4.1 Optimal transient growth

Transient growth is observed for all investigated cases. For three-dimensional analyses, the energy gain appears to be split into three broad regimes governed by distinctly different dynamics. For weak magnetic fields where the base velocity profile remains of a form consistent with Poiseuille duct flow (i.e.  $\mathit{Ha}\leqslant 1$ ), there is negligible change in the attainable growth $G_{max}$ , which still scales as ${\sim}Re^{2}$ and is largely independent of $\mathit{Ha}$ . Throughout this weak-MHD regime the time interval for maximum growth $\unicode[STIX]{x1D70F}_{opt}$ is similarly $\mathit{Ha}$ -independent with scaling ${\sim}Re$ . Conversely, a small, but finite, streamwise wavenumber dependence develops for non-zero field strengths, which remains mainly independent of $Re$ for all $\mathit{Ha}$ cases investigated. Both of these respective behaviours are correspondingly depicted in figure 3.

Figure 2.

Optimal transient growth as a function of Hartmann number for both three-dimensional and quasi-two-dimensional models at $Re=5000$ (——) and $15\,000$ (– – –). Also provided are the optimal transient growths predicted by Krasnov et al. (Reference Krasnov, Zikanov, Rossi and Boeck2010) using a three-dimensional analysis, and by Pothérat (Reference Pothérat2007) using the SM82 model. The horizontal axis is logarithmically and linearly scaled for $\mathit{Ha}\geqslant 2.2\times 10^{-2}$ and $\mathit{Ha}<2.2\times 10^{-2}$ , respectively.

Figure 3.

(a) Optimal time interval and (b) streamwise wavenumber as a function of Hartmann number for $Re=5000$ (——) and $15\,000$ (– – –). The black, grey and white symbols in (b) represent the variation of the eigenvector periodicity in the magnetic field direction as a function of Hartmann number. The periodicity in this direction is quantified through the wavenumber $\unicode[STIX]{x1D6FD}=2\unicode[STIX]{x03C0}/l_{y}$ , where $l_{y}$ is the wavelength of the eigenvector in the vertical $y$ -direction. The horizontal axis scaled as per figure 2. The vertical axis on (b) is logarithmically and linearly scaled for $k_{opt}\geqslant 2.2\times 10^{-7}$ and $k_{opt}<2.2\times 10^{-7}$ , respectively.

Figure 4.

Optimal transient growth plotted against the Reynolds based on (a) the Hartmann layer and (b) Shercliff layer thickness $R_{H}$ and $R_{S}$ , respectively, for both three-dimensional and quasi-two-dimensional models at $Re=5000$ (——) and $15\,000$ (– – –).

For $\mathit{Ha}>1$ , as the base flow topology morphs away from the hydrodynamic solution, the growth is strongly suppressed and correlated with a sharp decrease in the time horizon and streamwise wavelength for maximum amplification. Here the optimal time interval develops a strong dependence on $\mathit{Ha}$ and is now relatively independent of $Re$ . It is within this moderate-field-strength regime (hereafter moderate-MHD regime) in which close agreement with Krasnov et al. (Reference Krasnov, Zikanov, Rossi and Boeck2010) in terms of the optimal growth $G_{max}$ and scaling of $G_{max}\propto \mathit{Ha}^{-1.5}$ is seen for $10\leqslant \mathit{Ha}\leqslant 50$ at $Re=5000$ . However, given the narrow Hartmann number range over which this relation holds, it is unlikely that this represents a true power law scaling. The streamwise wavenumber of optimal disturbances is found to depend only on the Hartmann number (similarly for $\mathit{Ha}\lesssim 1$ ), indicating that eigenvector modulation is due primarily to electro-magnetic, rather than inertially driven, phenomena.

A final asymptotic regime is reached when the Reynolds number at the scale of the Hartmann layer thickness $R_{H}=Re/\mathit{Ha}$ reaches a transitional value of no lower than $R_{H}\approx 33.\dot{3}$ (this finding is supported with supplementary investigations at $Re=2000$ and $10\,000$ which are not published here). Below this value, $G_{max}$ exclusively follows a relationship governed by the Reynolds number built upon the Shercliff layer thickness $R_{S}=Re/\mathit{Ha}^{1/2}$ . The respective regimes of governance of these parameters are readily illustrated in figure 4(a,b) in which $G_{max}$ is plotted against $R_{H}$ and $R_{S}$ for a subset of the Hartmann numbers investigated. It is within this high Hartmann number regime where the flow is close to quasi-two-dimensional (hereafter Q2D-MHD regime) and a remarkable agreement with SM82 model predictions are found. Here, not only is the scaling of energy growth $G_{max}\propto \mathit{Ha}^{-0.44\,\pm \,0.07}$ closely captured by SM82, but also, strikingly, the amplification magnitudes are almost identical. The scaling of energy amplification is consistent with that known to dictate Q2D MHD critical linear and energy stability dynamics governed by Shercliff layer thickness $\unicode[STIX]{x1D6FF}_{S}=O(\mathit{Ha}^{-1/2})$ (Pothérat Reference Pothérat2007; Vo, Pothérat & Sheard Reference Vo, Pothérat and Sheard2017). The tendency of transient growth dynamics towards two-dimensionality is further illuminated by the transition of the optimal time interval and wavelength to a scaling of $\mathit{Ha}^{-1/2}$ and $\mathit{Ha}^{1/2}$ , respectively. The mechanisms leading to these effects are discussed in more detail throughout §§ 4.2 and 4.3.

4.2 Optimal mode topology

In this section the optimal eigenvectors are analysed in terms of their spatial form and initial componental energy distribution. Attention is focused on the three-dimensional optimal modes and how their characteristics change over transition from weak-MHD to Q2D-MHD regimes. In addition, a comparison between three-dimensional and SM82 optimal mode topology upon evolution to $\unicode[STIX]{x1D70F}_{opt}$ in the Q2D-MHD regime is presented. The initial disturbance topologies leading to optimal transient growth are visualised as iso-surfaces of the vertical component of vorticity $\unicode[STIX]{x1D714}_{y}$ in figure 5 for $0\leqslant \mathit{Ha}\leqslant 300$ at $Re=5000$ and figure 6 for $50\leqslant \mathit{Ha}\leqslant 800$ at $Re=15\,000$ . Displayed in figure 7 is the distribution of the total initial kinetic energy across the horizontal $\mathscr{E}_{x}$ , vertical $\mathscr{E}_{y}$ and streamwise $\mathscr{E}_{z}$ velocity components (i.e.  $u,v$ and $w$ ).

Figure 5.

Vertical component of vorticity $\unicode[STIX]{x1D714}_{y}$ iso-surfaces of the optimal eigenvector fields producing maximum transient amplification $G_{max}$ for (a–f) $\mathit{Ha}=0$ , $5$ , $10$ , $50$ , $150$ and $300$ at $Re=5000$ . Blue and yellow contours represent positive and negative vorticity, respectively. Contour levels are adjusted to approximately 90 % of the maximum magnitude of $\unicode[STIX]{x1D714}_{y}$ . The flow is from left to right in the positive $z$ -direction, with the magnetic field oriented vertically in the positive $y$ -direction. For clarity, the iso-surfaces are only plotted for $0\leqslant x\leqslant 1$ .

Figure 6.

Vertical component of vorticity $\unicode[STIX]{x1D714}_{y}$ iso-surfaces of the optimal eigenvector fields producing maximum transient amplification $G_{max}$ for (a–d) $\mathit{Ha}=50$ , $150$ , $300$ and $800$ at $Re=15\,000$ . Iso-surfaces specified as per figure 5.

Figure 7.

Distribution of initial eigenvector energy across horizontal $\mathscr{E}_{x}$ , vertical $\mathscr{E}_{y}$ and streamwise $\mathscr{E}_{z}$ velocity components (i.e.  $u,v$ and $w$ ), as a function of $\mathit{Ha}$ for (a) $Re=5000$ and (b) $Re=15\,000$ . Energy magnitudes are normalised using the total perturbation kinetic energy.

In the absence of a magnetic field, the optimal eigenvectors present as streamwise independent and aligned counter-rotating vortices, as shown in figure 5(a). For weak magnetic field strengths $0<\mathit{Ha}\lesssim 1$ , the only qualitative change in topology is that of slight wavelength dependence (i.e.  $k_{opt}\rightarrow 0$ as $\mathit{Ha}\rightarrow 0$ ). Thus, the initial energy in the streamwise velocity component is substantially smaller compared to the spanwise components, and the modes retain close similarity to the algebraic instabilities producing maximum transient amplification in wall-bounded hydrodynamic shear flows (Schmid & Henningson Reference Schmid and Henningson2001; Biau, Soueid & Bottaro Reference Biau, Soueid and Bottaro2008). In contrast, optimal mode topology in the moderate-MHD regime appears to be governed by $R_{H}$ . Upon transition to this regime, the optimal modes form obliquely aligned streamwise rolls overlapping within the side wall boundary layers and slanted out from the wall in the upstream direction, which is readily highlighted in figure 5(b,c). Streamwise independence is no longer present, and as seen in figure 7, a substantial increase in $\mathscr{E}_{z}$ relative to $\mathscr{E}_{x}$ and $\mathscr{E}_{y}$ is observed. Interestingly for both $Re$ at $\mathit{Ha}\gtrsim 2$ , $\mathscr{E}_{x}$ and $\mathscr{E}_{y}$ start to deviate away from one another due to a relative increase and decrease in magnitude, respectively. By $\mathit{Ha}\approx 30$ , the energy contained in the streamwise component is larger than the horizontal constituent for both Reynolds numbers.

It is within the moderate MHD-regime (i.e.  $350\gtrsim R_{H}\gtrsim 33.\dot{3}$ ) where $\mathscr{E}_{z}$ becomes globally the largest component, albeit only at relatively large field strengths of $\mathit{Ha}\approx 100$ and $\mathit{Ha}\approx 300$ for $Re=5000$ and $15\,000$ , respectively. The eigenvectors also present again as streamwise aligned oblique modes in this regime, yet now with increased vertical periodicity. This behaviour is depicted in figures 5(d) and 6(b,c), along with a quantification presented in figure 3(b). This change in vertical wavelength represents a counter-intuitive modification to the largest field-parallel length scale of the disturbance structures $\ell _{y}$ . The corresponding decrease in $\ell _{y}$ with increasing $\mathit{Ha}$ in this parameter space is in stark contrast to the natural progression of eddy anisotropy, which is often expected as field strengths are increased. Indeed anisotropy is caused by the Lorentz force increasing momentum diffusion of structures perpendicular to the field and the corresponding diffusivity scales as $\boldsymbol{B}_{0}^{2}$ (Sommeria & Moreau Reference Sommeria and Moreau1982; Pothérat, Sommeria & Moreau Reference Pothérat, Sommeria and Moreau2000; Pothérat & Klein Reference Pothérat and Klein2014).

However, upon transition to the final asymptotic regime $R_{H}\lesssim 33.\dot{3}$ , a breakdown in the periodicity of vortices in the vertical direction occurs, and a recovery of the parallel length scale tending towards that of the duct height $\ell _{y}\sim 2a$ is seen. As illustrated through figures 5(e,f) and 6(d), when the relative field strength is increased, $\ell _{y}$ correspondingly grows so that it exceeds that of the duct height (i.e.  $\ell _{y}\gtrsim 2a$ ); resulting in the eigenvector becoming nearly two-dimensional. Quantification of this behaviour is seen in the sharp decline of $\mathscr{E}_{y}$ in this regime due to the Lorentz force diffusing velocity differentials perpendicular to the field direction. Further to this, transition to this regime sees a change in the hierarchy of spanwise energy components such that $\mathscr{E}_{x}$ now dominates. The dominance of both the initial streamwise and horizontal energies over that of $\mathscr{E}_{y}$ is an additional signal to the transition to two-dimensionality of the optimal modes.

Figure 8.

Vertical component of vorticity $\unicode[STIX]{x1D714}_{y}$ iso-surfaces of the perturbation fields at time of maximum transient amplification $\unicode[STIX]{x1D70F}_{opt}$ for $\mathit{Ha}=800$ at (a) $Re=5000$ and (b)  $Re=15\,000$ . The leftmost figures in each of (a) and (b) are the iso-surfaces for three-dimensional analysis, the rightmost figures are a uniform extrusion of the $x$ – $z$ planar Q2D optimal disturbance fields in the $y$ -direction. Iso-surface colours as per figure 5. For clarity, the iso-surfaces are only shown for $0.75\leqslant x\leqslant 1$ .

A primary aim of the present work is to determine the validity of SM82 in reproducing three-dimensional transient growth predictions. As such, the disturbance topology in the Q2D-MHD regime using both MHD approaches at time of maximum energy gain $\unicode[STIX]{x1D70F}_{opt}$ is presented for comparison in figure 8 for $\mathit{Ha}=800$ at (a) $Re=5000$ and (b) $Re=15\,000$ . Here a uniform vertical extrusion of the SM82 disturbance fields is provided to help indicate the resultant averaging of the three-dimensional solution in the $y$ -direction as per (2.16). Consistent with the excellent agreement of optimal energy growth between both MHD models, here the structure of the disturbance fields are again similar: both feature overlapping alternating-sign vorticity structures adjacent to the side wall, and a series of opposite-signed vertically aligned vortices adjacent to the near-wall vorticity structures. Some differences are observed: the outer regions of the side wall layers display a qualitative discrepancy between models in the creation of cylindrical field-aligned vortices. For the full-MHD model, a quadratic three-dimensionality persists in the vortices to resemble the barrel-like turbulent structures theorised in Pothérat et al. (Reference Pothérat, Sommeria and Moreau2000) and reported in Mück et al. (Reference Mück, Günther, Müller and Bühler2000) and Pothérat (Reference Pothérat2012). In addition, the difference in structures across the two sub-layers within the Shercliff layer are reminiscent of that found in turbulent MHD flows by Krasnov, Zikanov & Boeck (Reference Krasnov, Zikanov and Boeck2012), albeit with slightly different topology dimensionality. Nonetheless, the closeness of the eigenvalue predictions highlights an interesting facet of the SM82 formulation. Although the exact nature of the flow may present small but noticeable differences, the averaging of the flow dynamics along the field direction results in the mean integrated quantities being suitably approximated. As evidenced in figure 8, this is even more so for higher field strengths, where at low $R_{S}$ (applicable to fusion applications) the difference in topology is suggested here to be negligible.

4.3 Transient growth mechanisms

This section quantifies the mechanisms leading to transient amplifications in an attempt to address whether instability in these duct flows can result exclusively from two-dimensional dynamics. An analysis of the energy growth in orthogonal velocity components is provided in § 4.3.1. To understand the fundamental flow dynamics which modulates mode topology and alters attainable energy gains, an analysis of the physical time scales for inertial and electro-magnetic phenomena is given in § 4.3.2. However, before proceeding with this analysis, the fundamental mechanisms that lead to transient growth wall-bounded flows are revisited.

For low- $\mathit{Ha}$ MHD and hydrodynamic flows, transient amplifications can exist due to the growth of vertically aligned and streamwise aligned modes, which are respectively analogous to the well-known Orr–Sommerfeld and Squire modes arising from analysis of the linear stability of parallel shear flows (Orr Reference Orr1907; Landahl Reference Landahl1980). Each of these has its own fundamental mechanisms for transient energy growth. The two-dimensional (here magnetic field-aligned) Orr–Sommerfeld modes undergo modest transient amplification due to the Orr mechanism; the wall-normal velocity disturbance leans upstream against the shear flow and rises up into Tollmien–Schlichting wave packets which are subsequently propagated downstream (Orr Reference Orr1907). In contrast, the streamwise invariant Squire modes achieve significant transient amplification via the lift-up mechanism: low streamwise velocity fluid closer to the duct wall is lifted up by the streamwise vortices, with higher velocity fluid subsequently replacing it (Landahl Reference Landahl1980). Beyond these fundamental mechanisms, a coupling exists between modes in which the transient growth is further amplified by the forcing of streamwise aligned vorticity by wall-normal velocity (Schmid & Henningson Reference Schmid and Henningson2001).

4.3.1 Componental energy analysis

To characterise the dominance and coupling of specific transient growth mechanisms at varying $\mathit{Ha}$ over evolution to $t=\unicode[STIX]{x1D70F}_{opt}$ , the componental kinetic energy growth for the three-dimensional optimal disturbance is analysed. The energy gain in the horizontal $\mathscr{G}_{x}$ , vertical $\mathscr{G}_{y}$ and streamwise $\mathscr{G}_{z}$ velocity components is shown as a function of $\mathit{Ha}$ for $Re=5000$ and $15\,000$ in figures 9(a) and 9(b), respectively.

Figure 9.

Amplification of perturbation kinetic energy $\mathscr{G}_{x}$ , $\mathscr{G}_{y}$ and $\mathscr{G}_{z}$ in the horizontal, vertical and streamwise velocity components (i.e.  $u,v$ and $w$ ), respectively, as a function of $\mathit{Ha}$ for (a)  $Re=5000$ and (b)  $Re=15\,000$ .

For low field strengths, transient growth is confined solely to the streamwise velocity components $\mathscr{G}_{z}$ . This is strongly indicative of lift-up mechanisms dominating transient amplifications, and is consistent with the presiding hydrodynamic response. However, for $\mathit{Ha}>1$ , a significant reduction in $\mathscr{G}_{z}$ correlated with an increase in spanwise components $\mathscr{G}_{x}$ and $\mathscr{G}_{y}$ signals impediment of this mechanism and a growing influence of the Orr mechanism forcing side wall-normal velocity. Regardless, the largest transient gains remain driven by lift-up responses, and only for $500\gtrsim R_{H}\gtrsim 33.\dot{3}$ and $3000\gtrsim R_{H}\gtrsim 50$ is significant growth found in all three velocity components at $Re=5000$ and $Re=15\,000$ , respectively. At transition to the Q2D-MHD regime, Orr related mechanisms become the largest component of transient growth for both Reynold numbers. Hence, for $R_{H}\lesssim 33.\dot{3}$ the lift-up mechanisms are largely suppressed, and the dominant two-dimensional Orr mechanisms can be sufficiently captured by the SM82 formulation. Transient growth is produced by the non-modal interaction of linear eigenmodes: revisiting figure 2, the substantially higher optimal growths found in the three-dimensional duct at smaller $\mathit{Ha}$ compared to the SM82 predictions can be understood in terms of the SM82 model having only a subset (the Orr modes) of the eigenmodes available to the three-dimensional model to interact for transient growth.

Figure 10.

Vertical component of vorticity $\unicode[STIX]{x1D714}_{y}$ iso-surfaces of the perturbation fields at time of maximum transient amplification $\unicode[STIX]{x1D70F}_{opt}$ for (a–d) $\mathit{Ha}=10$ , $50$ , $150$ and $300$ at $Re=5000$ . Iso-surfaces specified as per figure 5.

To provide a qualitative picture of these effects, the energy contained in orthogonal velocity components is qualitatively captured through comparison of the optimal eigenvectors for $Re=5000$ in figure 5(c–f), and their respective topology after time evolution to $\unicode[STIX]{x1D70F}_{opt}$ as seen in figure 10. For moderate-MHD regimes (figure 10 a,b), the perturbation field persists as streamwise aligned vortices consistent with lift-up mechanisms, yet now with their orientation slanted downstream due to bulk shear caused by a mild Orr related response. For higher field strengths pertaining to the Q2D-MHD regime (figure 10 c,d), modulation of the optimal modes is again split into two regions inside the Shercliff boundary layers with varying degrees of two-dimensionality. However, the dominant disturbance structures remain aligned with the magnetic field, and their mean radial growth is strongly reminiscent of Orr mechanisms providing moderate transient growth. Here suppression of the lift-up response through the lack of any streamwise orientation is qualitatively seen. The correlation of this behaviour with the preclusion of streamwise energy growth $\mathscr{G}_{z}$ further explains why the SM82 model is an excellent predictor of transient growth mechanisms towards large $\mathit{Ha}$ .

4.3.2 Time scale analysis

This section provides an analysis of the specific role electromagnetic and inertial phenomena play in determining transient growth and mode topology. To this end, the following non-dimensional time scale definitions are introduced:

1. (i) Vertical mean-flow energy transfer

   (4.1) $$\begin{eqnarray}\unicode[STIX]{x1D70F}_{z}=\frac{\ell _{z}}{U_{0}},\end{eqnarray}$$
   where $\ell _{z}$ is the streamwise wavelength of the optimal modes.

2. (ii) Streamwise mean-flow energy transfer

   (4.2) $$\begin{eqnarray}\unicode[STIX]{x1D70F}_{y}=\frac{\ell _{y}}{U_{0}},\end{eqnarray}$$
   where $\ell _{y}$ is the maximum vertical radius of a single vortex present in the optimal mode $\unicode[STIX]{x1D714}_{y}$ iso-surface.

3. (iii) Viscous dissipation

   (4.3) $$\begin{eqnarray}\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D708}}=Re[(\ell _{x}^{2})^{-1}+(\ell _{y}^{2})^{-1}+(\ell _{z}^{2})^{-1}]^{-1},\end{eqnarray}$$
   where $\ell _{x}$ is a characteristic length taken orthogonal to both $\boldsymbol{e}_{y}$ and $\boldsymbol{u}_{0}$ , here taken to be the Shercliff layer thickness $\unicode[STIX]{x1D6FF}_{S}$ .

4. (iv) Two-dimensionalisation by the Lorentz force

   (4.4) $$\begin{eqnarray}\unicode[STIX]{x1D70F}_{L}=\frac{1}{N}\left[1+\left(\frac{\ell _{y}}{\ell _{x}}\right)^{2}+\left(\frac{\ell _{y}}{\ell _{z}}\right)^{2}\right].\end{eqnarray}$$

5. (v) Hartmann friction

   (4.5) $$\begin{eqnarray}\unicode[STIX]{x1D70F}_{H}=R_{H}.\end{eqnarray}$$

6. (vi) Total bulk dissipation

   (4.6) $$\begin{eqnarray}\unicode[STIX]{x1D70F}_{d}=[\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D708}}^{-1}+\unicode[STIX]{x1D70F}_{L}^{-1}]^{-1}.\end{eqnarray}$$

All time scale quantities here are calculated at the optimal time period $\unicode[STIX]{x1D70F}_{opt}$ . All equations (4.1)–(4.6) are considered in non-dimensional form with dimensional scales consistent with those defined in § 2.

Figure 11.

Inverse time scales for (a)  $Re=5000$ and (b)  $15\,000$ are respectively plotted as a function of $\mathit{Ha}$ .

Equations (4.1)–(4.6) largely fall into one of two categories; those relating purely to inertial phenomena, and those influenced by electromagnetic forces. The former consists of the time scale for significant transference of energy from the mean base flow to the perturbation field in the vertical and streamwise directions $\unicode[STIX]{x1D70F}_{z}$ and $\unicode[STIX]{x1D70F}_{y}$ , respectively (analogous to the eddy turnover time, see Frisch (Reference Frisch1995)), in addition to the time scale for the dissipative effects of viscosity to act on the perturbation fields $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D708}}$ . The latter category comprises the time scale for the Lorentz force to suppress velocity differentials between transverse planes $\unicode[STIX]{x1D70F}_{L}$ , as well as that of energy dissipation due to friction in the Hartmann boundary layers $\unicode[STIX]{x1D70F}_{H}$ (derived from the SM82 model, see Pothérat et al. (Reference Pothérat, Sommeria and Moreau2000)). The total dissipative effects in the bulk flow due to both viscous and Lorentz forces are characterised by $\unicode[STIX]{x1D70F}_{d}$ . The inverse of each time scale defined in (4.1)–(4.6) is presented for $0.1\leqslant \mathit{Ha}\leqslant 800$ at $Re=5000$ and $15\,000$ in figure 11(a,b), respectively. Larger magnitudes of $\unicode[STIX]{x1D70F}^{-1}$ indicate a given physical mechanism modulating the dynamics of the linearised flow at a faster rate.

For all regime parameters investigated, the fastest acting mechanisms are one of the two inertial transfer times scales. For low- and moderate-MHD regimes, the lift-up related $\unicode[STIX]{x1D70F}_{z}$ acts most quickly, whereas the transition to the Q2D-MHD regime sees a switch to the Orr mechanism related streamwise equivalent $\unicode[STIX]{x1D70F}_{y}$ . This highlights an interesting point, despite the Lorentz force acting to bias the eigenspectrum towards Q2D eigenmodes with increasing Hartmann number, the growth mechanisms stem largely from inertial effects. In line with that presented for hydrodynamic flows by Barkley et al. (Reference Barkley, Blackburn and Sherwin2008), this suggests that the non-orthogonality driving transient growth emanates from the asymmetry of the advection operator describing inertial forces, rather than the operator relating to pressure (and here, Lorentz forces). For $Ha\gtrsim 5$ , the hierarchy of time scales changes in two significant ways. Firstly, Joule dissipation now acts more quickly than viscous dissipation. Secondly, the vertical inertial transfer time becomes faster than all other dissipative mechanisms at play. This change marks the transition to a regime where the transient growth is dominated by electromagnetic effects in the bulk. From this point onwards, the maximum amplification results from both inertial transfer terms, which are now also predominantly impeded by Joule dissipation.

Importantly, at $R_{H}\approx 33.\dot{3}$ the Joule dissipation in the bulk suddenly falls below the level of the dissipation in the Hartmann layer derived from the SM82 model. This signals a transition to two-dimensional dynamics where bulk velocity gradients along the magnetic field have been smoothed out by the action of the Lorentz force, and the dissipation concentrates where the strongest gradients remain – in the Hartmann layer. It is noteworthy that significant gradients still exist in the bulk of the optimal modes at the cross-over point between the two regimes. Nevertheless, the transition also corresponds to the point where transient growth becomes controlled by $R_{S}$ , which expresses the ratio between two-dimensional inertia and the Lorentz force, rather than $R_{H}$ in the three-dimensional MHD regime. Thus, although the point where bulk dissipation drops below the level of Hartmann layer dissipation is an unmistakable signature of a transition to quasi-two-dimensional MHD dynamics, the optimal mode retains noticeable three-dimensional features at values of $\mathit{Ha}$ well beyond this transition (Pothérat & Klein Reference Pothérat and Klein2014). Though somewhat counter-intuitive, the occurrence of two-dimensional dynamics has been recently observed in three-dimensional MHD turbulence by Baker et al. (Reference Baker, Pothérat, Davoust and Debray2018), and underlines that the link between topological and dynamical dimensionality is anything but obvious.