Thermoelectric magnetohydrodynamic flow in a liquid metal-infused trench

Abstract

We derive a mathematical model for steady, unidirectional, thermoelectric magnetohydrodynamic (TEMHD) flow of liquid lithium along a solid metal trench, subject to an imposed heat flux. We use a finite-element method implemented in COMSOL Multiphysics to solve the problem numerically, demonstrating how the fluid velocity, induced magnetic field and temperature change depending on the key physical and geometrical parameters. The observed flow structures are elucidated by using the method of matched asymptotic expansions to obtain approximate solutions in the limit where the Hartmann number is large and the trench walls are thin.

1. Introduction

The divertor is a vital component designed to absorb heat and waste products exhausted from a tokamak nuclear fusion reactor. The divertor must survive continuous extreme heat loads of around $10\ {\rm MW}\ {\rm m}^{-2}$, and a promising proposal is to coat it with a layer of flowing liquid lithium which can be constantly recycled. The liquid metal-infused trenches (LiMIT) concept has been devised to exploit the large heat flux experienced by the divertor to drive the lithium flow via thermoelectric effects. In this paper we derive and solve a model for the resulting flow along a single lithium-filled trench. Our aim is to determine how the flow properties depend on the applied magnetic field and heat flux, and how they can be beneficially influenced by varying the trench geometry.

The mathematical foundations of thermoelectric magnetohydrodynamic (TEMHD) duct flow were established by Shercliff (1979a,b). The proposal to use TEMHD effects to drive liquid lithium flow in fusion applications was pioneered at the University of Illinois at Urbana-Champaign (Jaworski 2009). The LiMIT concept (Ren et al. 2014; Xu, Curreli & Ruzic 2014; Fiflis et al. 2015) comprises of an array of parallel solid metal trenches that are filled with liquid lithium and cooled from below. As illustrated in figure 1, the combination of an applied magnetic field orthogonal to the trench walls and a temperature gradient parallel to the walls generates a thermoelectric current which, in turn, gives rise to a Lorentz force that drives the flow of lithium along the trench. The LiMIT concept has been tested experimentally in fusion devices (for example Ren et al. 2014; Fiflis et al. 2015; Xu et al. 2015; Andruczyk et al. 2020; Zuo et al. 2020; De Castro et al. 2021), albeit in the role of a ‘limiter’ rather than a divertor per se, and simulated numerically (Xu et al. 2014).

Figure 1. (a) Schematic of an array of trenches, indicating the orientation of the applied magnetic field. (b) Schematic of the heat flux, magnetic field and induced thermoelectric current in a single trench.

In this paper, we formulate a two-dimensional mathematical model for unidirectional liquid lithium flow in LiMIT which permits large parameter sweeps to be performed numerically and the trends to be understood using asymptotic analysis. In principle, our modelling approach can be adapted to describe other liquid metal plasma-facing components relying on TEMHD-driven flow, including free-surface ‘divertorlets’ (Fisher, Sun & Kolemen 2020; Saenz et al. 2022) and also liquid–metal concepts involving interfaces with porous media (Khodak & Maingi 2021). The divertorlet concept has some similarities to the LiMIT concept, except that the channels are aligned parallel to the toroidal magnetic field and an external current is applied to generate the required thermoelectric forcing.

Our model resembles classical models for magnetohydrodynamic (MHD) flow of liquid metal in ducts and channels, which have been of great theoretical interest ever since the 1950s (Shercliff 1953). Much of the early progress in MHD duct flow was based on asymptotic analysis in the limit of large Hartmann number ${\textit {Ha}}$ (see, for example, Walker, Ludford & Hunt (1971) and Hunt & Moreau (1976)). The basic structure in a rectangular duct consists of a core plug flow, with Hartmann layers of width $O({\textit {Ha}}^{-1})$ on walls orthogonal to the applied magnetic field and side layers of width $O({\textit {Ha}}^{-1/2})$ on walls parallel to the applied field (Temperley & Todd 1971). It is also well known that the conductivity of the duct walls can have a significant influence and give rise to velocity jets in the side layers (Hunt 1964).

More recently, advances in computational methods have made it possible to simulate MHD duct flow with novel geometries, including a fan-shaped insert (Kim, Lee & Lee 1997), sudden expansions (Mistrangelo & Bühler 2007) and other obstacles (Dousset 2014).

In § 2, we state the leading-order dimensionless governing equations and boundary conditions for steady, unidirectional flow of liquid lithium down a LiMIT-type trench. We also provide some preliminary numerical solutions for parameter regimes relevant to both theoretical and experimental set-ups. In § 3, we use asymptotic analysis in the limit of large Hartmann number to obtain approximate solutions to the problem when the trench walls are thin. In § 4, we perform several parameter sweeps to demonstrate how the flow properties depend on the geometrical and physical variables of interest. Finally, we provide a summary and discussion of our findings in § 5.

2. Mathematical model

2.1. Parameter values

Values and definitions of typical relevant thermophysical properties are given in table 1. The liquid metal inside the trench is taken to be lithium, and two candidate materials are considered for the trench itself: stainless steel 316 and tungsten.

Table 1. Mean values of relevant thermophysical properties for lithium, stainless steel and tungsten over temperatures in the range $[200\,^{\circ }{\rm C},400\,^{\circ }{\rm C}]$, aside from the Seebeck coefficient which has been averaged over $T\in [200\,^{\circ }{\rm C},277\,^{\circ }{\rm C}]$ due to lack of measurement data at higher temperatures (Kriessman 1953; Davison 1968; Choong 1975; Ho & Cho 1977; White & Collocott 1984; Itami, Shimoji & Shimokawa 1988; van der Marel et al. 1988; Fiflis et al. 2013; Tolias & EUROfusion MST1 Team 2017).

The set-up of a single periodic cell in an array of rectangular trenches is illustrated in figure 2. The geometry is characterised by the trench width $\mathcal {W}$, the heights ${t}_{B}$ and $\mathcal {H}$ of the base and the wide walls, the wall half-thickness ${t}_{S}$ and the height ${t}_{F}$ of the liquid layer above the trench. The dimensionless analogues of these quantities are scaled with $\mathcal {W}$, i.e. $(\mathcal {W},{t}_{B},\mathcal {H},{t}_{S},{t}_{F})=\mathcal {W}(2\omega,\varsigma,H,\tau,f)$. Following this scaling, we retain the parameter $\omega$ for the dimensionless trench half-width to allow ourselves the freedom to vary the width easily in the calculations below (thus taking $\mathcal {W}$ to be a typical but not necessarily exact value for the trench width).

Figure 2. Diagram summarising the coordinate system and the dimensional geometrical parameters associated with a rectangular trench (in red) with their dimensionless analogues (in black).

Typical trench geometrical parameters are shown in table 2, along with typical values of the applied magnetic field and heat flux, and the relevant dimensionless parameters. Three different cases are considered. The ‘default’ parameter regime has a typical geometrical set-up, with thermophysical properties relevant to a tungsten divertor, and is used as a base case in the computations performed below. The ‘SLiDE’ parameter regime is one relevant to the solid–liquid divertor experiment set-up, which involves a stainless steel divertor (Ruzic et al. 2011; Xu et al. 2013; Ren et al. 2014). The ‘ST40’ parameter regime is named after Tokamak Energy's spherical tokamak device, which uses a tungsten divertor, with impinging heat fluxes and applied magnetic fields relevant to a fusion reactor (McNamara et al. 2023). In all three cases the applied magnetic field is effectively horizontal, but the effects of the field being applied at a non-zero angle $\psi$ to the horizontal will also be examined in § 4.4.

Table 2. Default model parameter values, as well as typical values seen in the experimental set-up SLiDE, and proposed values for the ST40 tokamak (P.F. Buxton, private communication 2023). The superscripts $\mathfrak {L}$ and $\mathfrak {S}$ refer to values for the liquid lithium and the solid metal, respectively.

2.2. Governing equations

Here we briefly state the basic equations governing liquid–metal TEMHD flow. In the non-relativistic limit where the displacement current is negligible, the magnetic flux density $\boldsymbol {B}$ satisfies Gauss’ and Ampère's Laws, namely

(2.1a,b)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{B}=0,\quad \boldsymbol{\nabla}\times\boldsymbol{B}=\mu\boldsymbol{J}, \end{equation}

where $\boldsymbol {J}$ is the electrical current density and $\mu$ is the magnetic permeability. Since the magnetic susceptibilities of lithium, SS 316 and tungsten are all very small (see table 1), we treat $\mu$ as a constant throughout. The incompressible Navier–Stokes equations, including the Lorentz body force, may thus be expressed in the form

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

(2.2b)$$\begin{gather}\rho\left(\frac{\partial\boldsymbol{u}}{\partial t}+(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{u}\right)=-\boldsymbol{\nabla}p +\rho\nu\nabla^2\boldsymbol{u}+\rho\boldsymbol{g}+\frac{1}{\mu}(\boldsymbol{\nabla}\times\boldsymbol{B})\times\boldsymbol{B}, \end{gather}$$

where $\boldsymbol {u}$ and $p$ denote the liquid velocity and pressure, while $\rho$, $\nu$ and $\boldsymbol {g}$ denote the density, kinematic viscosity and gravitational acceleration (all assumed constant). Finally, we assume that internal viscous and Ohmic heating are negligible compared with large externally applied heat flux $\mathcal {Q}$, so the temperature $T$ satisfies the heat equation

(2.3)\begin{equation} \rho c\left(\frac{\partial T}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} T\right)=-\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{q}, \end{equation}

where $c$ is the heat capacity and $\boldsymbol {q}$ the heat flux.

Thermoelectric effects enter through Ohm's and Fourier's constitutive laws (Shercliff 1979a), which are modified to

(2.4a,b) \begin{equation} \boldsymbol{J}=\sigma\left(\boldsymbol{E}+\boldsymbol{u}\times\boldsymbol{B}-S{\boldsymbol{\nabla}}T\right), \quad {\boldsymbol{q}}=-k{\boldsymbol{\nabla}}T+ST\boldsymbol{J}, \end{equation}

respectively, where $\boldsymbol {E}$ is the electric field, $\sigma$ is the electrical conductivity, $k$ is the thermal conductivity and $S$ is the Seebeck coefficient of the medium. The additional thermoelectric terms (proportional to $S$) can be derived from the Onsager reciprocal relations in thermodynamics (Callen 1948), and physically they arise as a result of the Seebeck, Peltier and Thomson effects (Shercliff 1979a).

The relative importance of the final term in (2.4b) is measured by the dimensionless grouping $\sigma S^2\mathcal {W}\mathcal {Q}/k^2$, which is always small (of order $10^{-2}$), so this term is neglected henceforth. Moreover, taking the curl of (2.4a) to eliminate the electric field also eliminates the Seebeck term, so the induction equation

(2.5)\begin{equation} \frac{\partial\boldsymbol{B}}{\partial t}=\boldsymbol{\nabla}\times(\boldsymbol{u}\times\boldsymbol{B})+\eta\nabla^2\boldsymbol{B} \end{equation}

is unaffected (where $\eta =1/(\sigma \mu )$ is the magnetic diffusivity). However, thermoelectric effects enter the model through the boundary conditions at the interface between the liquid lithium and the solid trench, where the Seebeck coefficient is discontinuous. Since the liquid velocity is zero at this interface, continuity of the tangential electric field leads to the boundary condition

(2.6)\begin{equation} \left[\boldsymbol{n}\times\left(\eta\boldsymbol{\nabla}\times\boldsymbol{B}+S\boldsymbol{\nabla}T\right)\right]_-^+=\boldsymbol{0}, \end{equation}

where $[{\cdot }]_-^+$ denotes the jump in a quantity across the solid–liquid interface, whose unit normal is $\boldsymbol {n}$. Thus, a magnetic field can be generated by a temperature gradient parallel to the interface.

2.3. Modelling assumptions

We proceed to set out a model of steady unidirectional flow of liquid lithium along a single trench within a periodic array. We assume that the geometry and the flow are all uniform along the trench, with velocity given by

(2.7)\begin{equation} \boldsymbol{u}(\boldsymbol{x})=w(x,y)\boldsymbol{k}, \end{equation}

in the coordinate system depicted in figure 2. We consider a uniform applied magnetic field $\mathcal {B}^{a}$ inclined at an angle $\psi$ to the $x$-axis, which induces a field $b$ along the trench, i.e.

(2.8)\begin{equation} \boldsymbol{B}(\boldsymbol{x})=\mathcal{B}^{{a}}\left(\boldsymbol{i}\cos\psi+\boldsymbol{j}\sin\psi\right)+b(x,y)\boldsymbol{k}. \end{equation}

As usual in unidirectional flow, the pressure $p$ in general varies linearly with distance along the trench. However, since we assume the free surface is uniform in the $z$-direction, it follows that the down-trench pressure gradient must be zero. The variations in pressure due to the induced magnetic field are measured by the modified Weber number ${We}=(\nu /\eta )(\rho \mathcal {W}\mathcal {U}^2/\gamma )$, where $\gamma$ is the surface tension and $\mathcal {U}$ is a typical scale for the velocity, which will be chosen below. Since this parameter is small in all of the cases considered here, we can take the pressure to be purely hydrostatic and, hence, the free surface to be flat (as depicted in figure 2).

The plasma above the trench is assumed to provide a uniform heat source $\mathcal {Q}$ at the free surface, while the cooling pipes below the trench are effectively treated as a perfect heatsink, so the base of the trench is held at a fixed temperature $T_0$. According to the references cited in table 1, the physical parameters $\sigma$, $k$ and $\nu$ could vary by up to 50 % over the temperature ranges of interest. Nevertheless, in the first instance we assume that all material properties are constants, i.e. that any temperature dependencies may be neglected.

2.4. Non-dimensionalisation

The problem is non-dimensionalised as follows, with dimensionless variables denoted with hats. We scale all lengths with the width of the trench $\mathcal{W}$, that is,

(2.9)\begin{equation} \left(x,y\right) =\mathcal{W}\left(\hat{x},\hat{y}\right). \end{equation}

As noted in § 2.1, although the dimensionless trench width may be set to unity without loss of generality, it is convenient to allow it to vary, and we denote it by $2\omega$. Since the pressure gradient along the trench must be zero, given our modelling assumptions, the only remaining unknowns are the temperature $T$ and the velocity $w$ and induced magnetic field $b$ along the trench. The temperature is scaled using the imposed heat flux $\mathcal{Q}$, while $w$ and $b$ are scaled using a typical velocity $\mathcal {U}$ which will be selected below. We thus set

(2.10a–c)\begin{equation} T =T_{0}+\left(\dfrac{\mathcal{Q}\mathcal{W}}{k^{\mathfrak{L}}}\right)\hat{T}, \quad w =\mathcal{U}\hat{w}, \quad b =\epsilon {Rm}\, \mathcal{B}^{{a}}\hat{b}, \end{equation}

where $\epsilon \ll 1$ and ${Rm}\ll 1$ are the reciprocal Hartmann number and magnetic Reynolds number, respectively, given by

(2.11a,b)\begin{equation} \epsilon =\dfrac{1}{{\rm Ha}}=\dfrac{1}{{\mathcal{B}}^{{a}}{\mathcal{W}}}\sqrt{\dfrac{\rho\nu}{\sigma^{\mathfrak{L}}}}, \quad {Rm} =\mathcal{UW}\mu\sigma^{\mathfrak{L}}. \end{equation}

The appropriate velocity scale is found to be proportional to the difference in Seebeck coefficient between liquid and solid, namely

(2.12)\begin{equation} \mathcal{U} =\dfrac{\mathcal{Q}(S^{\mathfrak{L}}-S^{\mathfrak{S}})}{k^{\mathfrak{L}}\mathcal{B}^{{a}}}. \end{equation}

This scaling reflects the way the flow is driven by the heat flux $\mathcal{Q}$ applied to the free surface, through the Seebeck effect, while increasing the applied magnetic field too much slows down the flow, as observed qualitatively by Xu et al. (2014).

The dimensionless versions of the governing equations and boundary conditions, evaluated in each domain and on each boundary, are summarised in figure 3. For simplicity, we restrict attention to trenches with reflectional symmetry about $\hat {x}=0$. The problem contains several more dimensionless parameters, namely

(2.13a–c)\begin{equation} \varGamma = \frac{\epsilon \mathcal{W}^2g\sin\theta}{\nu\mathcal{U}}, \quad \varSigma_{\mathfrak{SL}}=\dfrac{\sigma_{\mathfrak S}}{\sigma_{\mathfrak L}}, \quad K_{\mathfrak{SL}}=\dfrac{k_{\mathfrak S}}{k_{\mathfrak L}}. \end{equation}

The solid-to-liquid electrical conductivity ratio is denoted by $\varSigma _{\mathfrak {SL}}$, so the regimes ${\varSigma _{\mathfrak {SL}}\ll 1}$ and $\varSigma _{\mathfrak {SL}}\gg 1$ correspond to the limits of perfectly electrically insulating or conducting trench walls, respectively. A similar analogy holds for the thermal conductivity ratio $K_{\mathfrak {SL}}$. Both of these constants are $O(1)$ in practice. The dimensionless constant $\varGamma$ measures the importance of gravity in driving the flow along the trench. It is generally small (see table 2) but not completely negligible, so we retain this term for the time being.

Figure 3. A summary of the steady, unidirectional TEMHD trench problem, assuming symmetry about $\hat {x} = 0$.

Under the assumptions set out in § 2.3, the induction equation (2.5) and momentum equation (2.2b) in the liquid reduce to (2.18a) and (2.18b) (in figure 3), respectively. Because we have neglected the thermoelectric heat flux and viscous dissipation, the thermal problem decouples from the rest of the equations. The temperature $T$ satisfies Laplace's equation in the lithium and the trench wall, with standard continuity conditions at the interface between them, and is driven by a dimensionless unit heat flux at the free surface. Indeed, the whole system is driven by this heat flux, which is the only inhomogeneous term in the problem. The resulting temperature gradient induces a magnetic field through the Seebeck boundary condition at the solid–liquid interface, and the induced field in turn generates the flow along the trench.

As in previous MHD modelling of liquid–metal flow in fusion applications (e.g. Lunz & Howell 2019), the boundary condition (2.15a) (in figure 3) on the induced magnetic field assumes that the region above the trench is effectively a vacuum. From the scalings (2.10), we estimate the thermoelectric current density generated in the lithium by

(2.21)\begin{equation} |\boldsymbol{J}|\sim\frac{\epsilon\mathcal{Q}\sigma^{\mathfrak{L}}\left(S^{\mathfrak{L}}-S^{\mathfrak{S}}\right)}{k^{\mathfrak{L}}} \sim 200\ \mathrm{kA}\ \mathrm{m}^{-2}. \end{equation}

In comparison, Jaworski, Khodak & Kaita (2013) suggests a typical steady state current density of order $10\ \mathrm {kA}\ \mathrm {m}^{-2}$ entering from the plasma in a fusion device, with transient currents of up to $100\ \mathrm {kA}\ \mathrm {m}^{-2}$ possible during disturbances such as so-called edge-localised modes. It is therefore reasonable to neglect such effects in our steady state model, though they could be included in future refinements, especially if we want to consider time dependence and linear stability (Fiflis et al. 2016).

Due to the ease of design and manufacture, we focus on rectangular trenches in this paper. The geometry of such a trench can be characterised by several geometrical parameters, which are shown visually in figure 2 and whose typical values are listed in table 2.

2.5. Numerical solutions

The problem shown in figure 3 was solved numerically using finite element methods implemented in COMSOL Multiphysics simulation software. The numerical simulations were tested for convergence with respect to mesh refinement and robustness with respect to small changes in the physical parameters. The solutions were found to be well converged provided the mesh was sufficiently refined in boundary layers and near sharp corners in the fluid and solid domains. With the leading-order problem shown in figure 3, it was possible to complete a simulation on a desktop computer obtaining the temperature, induced magnetic field and the fluid velocity in approximately 5–10 s.

Plots of the temperature, velocity and magnetic field in the three different parameter regimes described in § 2.1 are shown in figure 4. In all three cases, the temperature appears to be approximately a linear function of $\hat {y}$, reaching its peak at the free surface where the external heat flux is applied, although with some variation in the $\hat {x}$-direction near the walls due to the change in thermal conductivity. Meanwhile, the velocity plots (figure 4b,e,h) show an obvious boundary layer structure which arises because the MHD problem (shown in figure 3) becomes singular as $\epsilon \to 0$. There are Hartmann boundary layers of width $O(\epsilon )$ at the vertical walls, as well as a weaker ‘side’ layer of width $O(\epsilon ^{1/2})$ near the horizontal boundaries at the base and the tops of the walls (see, for example, Temperley & Todd (1971) for explanations of these scalings with respect to the Hartmann number). Outside these layers, it appears that the velocity inside the trench approaches a constant as $\epsilon \to 0$. The boundary layers can also be seen in the induced magnetic field plots (figure 4c,f,i), which show a magnetic field dipole. The boundary layer structure will be further elucidated in § 3.

Figure 4. Plots of the (a,d,g) temperatures, (b,e,h) velocities and (c,f,i) induced magnetic fields for the ‘default’, ‘SLiDE’ and ‘ST40’ parameter cases.

Another important feature particularly evident in figure 4(b,h) is the appearance of a conducting-base velocity jet near the bottom of the trench. This jet arises due to the fact that the electrically conducting base is in direct contact with a vacuum, which forces the electrical current loops to close partially in the fluid. The velocity jet appears near the wall that is parallel to the applied magnetic field, as observed by Hunt (1964).

Regarding the SLiDE set-up, we note that the velocities shown in figure 4(e) are much faster than observed in practice (Xu et al. 2013) and that the maximum temperature shown in figure 4(d) is close to the boiling point of lithium (1615 K). In practice a temperature limit of ${\approx }450\,^\circ \text {C}$ (${\approx }720\ \text {K}$) is often imposed to prevent excessive lithium evaporation (e.g. Smolentsev 2021). This behaviour occurs because the magnetic field in this case is relatively weak (0.07 T) while the impinging heat flux is fusion-relevant ($10\ {{\rm MW}\ {\rm m}}^{-2}$), both of which increase the velocity scaling in (2.12). However, we note that the heat flux in the SLiDE experiments is confined to a narrow electron beam rather than being applied uniformly along the trench, as assumed in our model.

2.6. Thin wall approximation

In general, analytical progress can be made with the steady problem summarised in figure 3 only when the trench sidewalls and base are thin. In the limit where the wall thickness $\tau$ and base thickness $\varsigma$ tend to zero, the governing equations need to be solved only in the liquid domain, with the walls and base represented by effective boundary conditions, as illustrated in figure 5. Here we introduce several further dimensionless parameters, namely

(2.22a–d)\begin{equation} c_{{S}}=\dfrac{\epsilon}{\varSigma_{\mathfrak{SL}}^{{S}}\tau}, \quad c_{{B}}=\dfrac{\epsilon^{1/2}}{\varSigma_{\mathfrak{SL}}^{{B}}\varsigma}, \quad \kappa_{{S}}=\dfrac{1}{K_{\mathfrak{SL}}^{{S}}\tau}, \quad \kappa_{{B}}=\dfrac{K_{\mathfrak{SL}}^{{B}}}{\varsigma}, \end{equation}

in which the superscripts ‘$S$’ and ‘$B$’ denote the parameter values on the sidewalls and the base, respectively. At present, the sidewalls and base of the trench are made of the same material; however, the most interesting distinguished limit occurs when the above constants are all $O(1)$ while $\epsilon$, $\tau$ and $\varsigma$ tend to zero.

Figure 5. Summary of the thin-wall problem for a rectangular trench.

We note that the temperature in the thin-wall problem is still decoupled from the MHD problem. In the regime where $K_\mathfrak {SL}^{S}$ is large, the solution to the temperature problem in figure 5 is given by

(2.23)\begin{equation} \hat{T}^{\mathfrak{L}}(\hat{x},\hat{y})=\hat{y}+\dfrac{1}{\kappa_{{B}}}, \end{equation}

which is a linear function of $\hat {y}$, in approximate agreement with the full numerical solutions shown in figure 4. In general, however, the temperature must still be determined numerically, so in the subsequent analysis we solve the remaining MHD problem, treating the liquid temperature as an in principle known function.

3. Asymptotic analysis

3.1. Outer regions and boundary layers

We use matched asymptotic expansions to analyse the thin-wall two-dimensional TEMHD trench problem set out in figure 5 in the limit as $\epsilon \rightarrow 0$. To simplify matters, we take $\psi =0$, so that the applied magnetic field is strictly in the toroidal direction in a tokamak configuration. The bulk flow inside the trench ($0<\hat {y}< H$) is primarily driven by thermoelectric effects in the boundary layer of thickness $O(\epsilon )$ at the trench wall $\hat {x}=\omega$. The velocity in the outer region above the trench ($\hat {y}>H$) is driven by gravity, and is an order of magnitude larger. These two regions are matched through a crossover region near $\hat {y}=H$ of thickness $O(\epsilon ^{1/2})$, and there is also a ‘side’ layer near the base $\hat {y}=0$ of thickness $O(\epsilon ^{1/2})$. We will state the leading-order problems satisfied in these two layers, although they are not in general amenable to analytical solution. We do not consider the further inner regions near the corners of the domain, first noted by Todd (1967), as they do not affect the leading-order outer solution.

3.2. Outer problem above the trench

Above the trench walls, by integrating the Navier–Stokes equation in figure 5 from $\hat{x}=0$ to $\hat{x}=\omega$, we find

(3.1)\begin{equation} \dfrac{\rm d}{{\rm d}\hat{y}}\int_{0}^{\omega}\hat{w}(\hat{x},\hat{y})\,{\rm d}\hat{x}=\dfrac{\varGamma\omega}{\epsilon}(H+f-\hat{y}) \quad\text{in }H < \hat{y} < H+f. \end{equation}

From this solvability condition, we deduce that $\hat {{w}}$ must be order $1/\epsilon$ above the trench. To match with an $O(1)$ velocity field $\hat {w}$ inside the trench, we must have

(3.2)\begin{equation} \int_{0}^{\omega}\hat{w}(\hat{x},\hat{y})\,{\rm d}\hat{x}= \dfrac{\varGamma}{2\epsilon}(\hat{y}-H)(H+2f-\hat{y})+C \quad\text{in }H<\hat{y}< H+f, \end{equation}

where $C$ is an order-one constant.

We also have that $\hat {w}$ is independent of $\hat {x}$ to leading order, and thus

(3.3)\begin{equation} \hat{w}(\hat{x},\hat{y})\sim \dfrac{\varGamma}{2\epsilon}(\hat{y}-H)(H+2f-\hat{y})+O(1) \quad\text{in }H<\hat{y}< H+f. \end{equation}

By seeking further terms in an asymptotic expansion, one can show that the induced magnetic field $\hat {b}^{\mathfrak{L}}$ above the trench is zero to all algebraic orders in $\epsilon$, whilst the corrections to the velocity are all constants.

3.3. Outer problem in the trench

The governing equations in the bulk of the trench are shown inside the shadowed box in figure 5. At leading order, solutions to this problem inside the trench are given by

(3.4)\begin{equation} \hat{w}(\hat{x},\hat{y})\sim\hat{w}_0(\hat{y}), \quad \hat{b}(\hat{x},\hat{y})\sim-\varGamma\hat{x} \quad\text{in }0<\hat{y}< H. \end{equation}

The outer problem on its own does not completely determine the form of $\hat {w}_{0}(\hat{y})$, which must therefore be found by matching with the Hartmann layer at the trench wall.

3.4. Inner problem near the right wall

Letting $\hat{x}=\omega -\epsilon \hat {X}_{R}$, with ${\hat {w}(\hat{x},\hat{y})}={\hat {w}_{R}(\hat{X}_{R},\hat{y})}$ and ${\hat {b}^{\mathfrak {L}}(\hat{x},\hat{y})}={\hat {b}_{R}^{\mathfrak {L}}(\hat{X}_{R},\hat{y})}$, one obtains the leading-order problem

(3.5a)$$\begin{gather} \dfrac{\partial^{2}{\hat{b}^{\mathfrak{L}}_{{R},0}}}{\partial\hat{X}_{R}^{2}}-\dfrac{\partial{\hat{w}_{{R},0}}}{\partial\hat{X}_{R}} =0 \quad \text{in }\hat{X}_{R}>0,\quad0<\hat{y}< H, \end{gather}$$

(3.5b)$$\begin{gather}\dfrac{\partial^{2}{{\hat{w}_{{R},0}}}}{\partial\hat{X}_{R}^{2}}-\frac{\partial{\hat{b}^{\mathfrak{L}}_{{R},0}}}{\partial\hat{X}_{R}} =0 \quad \text{in }\hat{X}_{R}>0,\quad0<\hat{y}< H, \end{gather}$$

(3.5c)$$\begin{gather}{\hat{w}_{{R},0}} =0{\quad} \text{on }\hat{X}_{R}=0,\quad0<\hat{y}< H, \end{gather}$$

(3.5d)$$\begin{gather}-\dfrac{\partial{\hat{b}^{\mathfrak{L}}_{{R},0}}}{\partial\hat{X}_{R}}+c_{{S}}{\hat{b}^{\mathfrak{L}}_{{R},0}}=\dfrac{\partial{\hat{T}^{\mathfrak{L}}}}{\partial\hat{y}}(\omega,\hat{y}) {\quad} \text{on }\hat{X}_{R}=0,\quad0<\hat{y}< H, \end{gather}$$

along with leading-order matching conditions

(3.5e)$$\begin{gather} {\hat{w}_{{R},0}(\hat{X}_{R},\hat{y})} \to\hat{w}_{0}(\hat{y})\text{ as }\hat{X}_{R}\to\infty, \end{gather}$$

(3.5f)$$\begin{gather}{{\hat{b}_{{R},0}^{\mathfrak{L}}(\hat{X}_{R},\hat{y})}} \to-\varGamma\omega\text{ as }\hat{X}_{R}\to\infty. \end{gather}$$

The problem (3.5) has solution

(3.6a)$$\begin{gather} {\hat{b}^{\mathfrak{L}}_{{R},0}}(\hat{X}_{R},\hat{y}) =\dfrac{1}{1+c_{{S}}}\left(\dfrac{\partial{\hat{T}^{\mathfrak{L}}}}{\partial\hat{y}}\left(\omega,0\right)+\omega\varGamma c_{{S}}\right){\rm e}^{-\hat{X}_{R}}-\omega\varGamma, \end{gather}$$

(3.6b)$$\begin{gather}{\hat{w}_{{R},0}}(\hat{X}_{R},\hat{y}) =\dfrac{1}{1+c_{{S}}}\left(\dfrac{\partial{\hat{T}^{\mathfrak{L}}}}{\partial\hat{y}}\left(\omega,0\right)+\omega\varGamma c_{{S}}\right)\left(1-{\rm e}^{-\hat{X}_{R}}\right), \end{gather}$$

which also determines the leading-order outer velocity in the trench, namely

(3.6c)\begin{equation} \hat{w}_{0}(\hat{y})=\dfrac{1}{1+c_{{S}}}\left(\dfrac{\partial{\hat{T}^{\mathfrak{L}}}}{\partial\hat{y}}\left(\omega,0\right)+\omega\varGamma c_{{S}}\right). \end{equation}

This simple formula for the bulk velocity demonstrates the importance of the temperature gradient in driving TEMHD flow. As will be demonstrated below, the expression (3.6c) shows good agreement with numerical results for the bulk velocity, even when the sidewall thickness $\tau$ is not especially small. Equation (3.6c) also resembles that obtained in Xu et al. (2015) in the limit of large ${\textit {Ha}}$.

3.5. Composite approximations

The solutions we have obtained thus far allow us to create composite expressions which provide good approximations to the velocity, magnetic field and temperature profiles inside the trench, even when the walls are no longer thin, provided we are not close to any horizontal boundary layers near the top or bottom of the trench. In $0<\hat {y}< H$ and $0<\hat {x}<\omega -\tau$ we have

(3.7a)$$\begin{gather} {{{\hat{w}}}}(\hat{x},\hat{y}) \approx\dfrac{1}{1+c_{{S}}}\left(\dfrac{\partial{\hat{T}^{\mathfrak{L}}}}{\partial\hat{y}}(\omega-\tau,\hat{y})+(\omega-\tau)\varGamma c_{{S}}\right)\left(1-2{\rm e}^{-\left(\omega-\tau\right)/\epsilon}\cosh\left(\dfrac{\hat{x}}{\epsilon}\right)\right), \end{gather}$$

(3.7b)$$\begin{gather}{\hat{b}^{\mathfrak{L}}}(\hat{x},\hat{y}) \approx\dfrac{2}{1+c_{{S}}}\left(\dfrac{\partial\hat{T}^{\mathfrak{L}}}{\partial\hat{y}}(\omega-\tau,\hat{y})+(\omega-\tau)\varGamma c_{{S}}\right){\rm e}^{-(\omega-\tau)/\epsilon}\sinh\left(\dfrac{\hat{x}}{\epsilon}\right)-\varGamma\hat{x}. \end{gather}$$

From (3.7b), one can construct an approximate expression for the overall induced magnetic field, including that inside the solid, assuming that it is approximately linear inside the trench wall, i.e.

(3.8)\begin{equation} {\hat{b}}(\hat{x},\hat{y})\approx\begin{cases} {\hat{b}^{\mathfrak{L}}}(\hat{x},\hat{y}) & \text{for }0\le\hat{x}\le\omega-\tau,\\ {\hat{b}^{\mathfrak{L}}}(\omega-\tau,\hat{y})\left(\dfrac{\omega-\hat{x}}{\tau}\right) & \text{for }\omega-\tau\le\hat{x}\le\omega. \end{cases}\end{equation}

Furthermore, we can use the solution for the temperature, given by (2.23), to construct an approximate vertical temperature profile in the middle of the trench. This approximation is given by

(3.9)\begin{equation} T(0,y)\approx\begin{cases}\displaystyle \frac{\hat{y}}{K_{\mathfrak{SL}}^{{B}}} & \text{if }0<\hat{y}\le\varsigma,\\ \displaystyle \hat{y} + \frac{\varsigma}{K_{\mathfrak{SL}}^{{B}}}-\varsigma & \text{if }\varsigma<\hat{y}\le H+f,\end{cases} \end{equation}

and applies when $K_{\mathfrak {SL}}^{{B}}$ is small.

3.6. Inner problem near the base

Upon rescaling $\hat{y}=\epsilon ^{1/2}\hat{Y}_{B}$, with ${\hat {w}(\hat{x},\hat{y})}={\hat {w}_{B}(\hat {x},\hat{Y}_{B})}$ and ${\hat {b}^{\mathfrak {L}}(\hat{x},\hat{y})}={\hat {b}_{B}^{\mathfrak {L}}(\hat{x},\hat{Y}_{B})}$, we obtain leading-order problem

(3.10a)$$\begin{gather} \dfrac{\partial^{2}{\hat{b}^{\mathfrak{L}}_{{B},0}}}{\partial\hat{Y}_{B}^{2}}+\dfrac{\partial{\hat{w}_{{B},0}}}{\partial\hat{x}} =0{\quad} \text{in }\hat{Y}_{B}>0,\quad0<\hat{x}<\omega, \end{gather}$$

(3.10b)$$\begin{gather}\dfrac{\partial^{2}{{\hat{w}_{{B},0}}}}{\partial\hat{Y}_{B}^{2}}+\frac{\partial{\hat{b}^{\mathfrak{L}}_{{B},0}}}{\partial\hat{x}}+\varGamma =0{\quad} \text{in }\hat{Y}_{B}>0,\quad0<\hat{x}<\omega, \end{gather}$$

with boundary conditions

(3.10c)$$\begin{gather} \hat{w}_{{B},0}= \dfrac{\partial{\hat{b}^{\mathfrak{L}}_{{B},0}}}{\partial\hat{Y}_{B}}-c_{{B}}{\hat{b}^{\mathfrak{L}}_{{B},0}} =0,\quad \text{on}\ \hat{Y}_{B}=0, \end{gather}$$

(3.10d)$$\begin{gather}\dfrac{\partial{\hat{w}_{{B},0}}} {\partial\hat{x}}= {\hat{b}^{\mathfrak{L}}_{{B},0}} =0, \quad \text{on}\ \hat{x}=0, \end{gather}$$

and far-field conditions

(3.10e)\begin{equation} {\hat{w}_{{B},0}} \to\dfrac{1}{1+c_{{S}}}\left(\dfrac{\partial{\hat{T}^{\mathfrak{L}}}}{\partial\hat{y}}\left(\omega,0\right)+\omega\varGamma c_{{S}}\right), \quad {\hat{b}^{\mathfrak{L}}_{{B},0}} \to-\varGamma\hat{x}\quad \text{as}\ \hat{Y}_{B}\to\infty, \end{equation}

as well as the matching condition

(3.10f)\begin{equation} \left(1+c_{{S}}\right){\hat{w}_{{B},0}}+c_{{S}}{{\hat{b}^{\mathfrak{L}}_{{B},0}}}=\dfrac{\partial{\hat{T}^{\mathfrak{L}}}}{\partial\hat{y}}\left(\omega,0\right)\quad\text{at }\hat{x}=\omega. \end{equation}

The effective boundary condition (3.10f) is derived by analysing the ‘inner bottom right-hand corner’ problem below in § 3.7.

The problem (3.10) appears to be analytically intractable in general, but we show in the Appendix how it may be solved in some limiting cases, and that the resulting solutions display typical physical phenomena, including wall jets and the possibility of return flow.

3.7. Inner problem in bottom right-hand corner

After rescaling both $\hat{x}=\omega -\epsilon \hat{X}_{R}$ and $\hat{y}=\epsilon ^{1/2}\hat{Y}_{B}$, with ${\hat {w}(\hat{x},\hat{y})}={\hat {w}_{BR}(\hat {X}_{R},\hat{Y}_{B})}$ and ${\hat {b}^{\mathfrak {L}}(\hat{x},\hat{y})}={\hat {b}_{BR}^{\mathfrak {L}}(\hat{X}_{R},\hat{Y}_{B})}$, we obtain the leading-order problem

(3.11a)\begin{equation} \dfrac{\partial^{2}{\hat{b}^{\mathfrak{L}}_{{BR},0}}}{\partial\hat{X}_{R}^{2}}-\dfrac{\partial{\hat{w}_{{BR},0}}}{\partial\hat{X}_{R}} = \dfrac{\partial^{2}{\hat{w}_{{BR},0}}}{\partial\hat{X}_{R}^{2}}-\frac{\partial{\hat{b}^{\mathfrak{L}}_{{BR},0}}}{\partial\hat{X}_{R}} =0,{\quad} \text{in }\hat{X}_{R}>0,\quad\hat{Y}_{B}>0,\end{equation}

along with boundary and matching conditions

(3.11b)$$\begin{gather} {\hat{w}_{{BR},0}} =0, \quad -\dfrac{\partial{\hat{b}^{\mathfrak{L}}_{{BR},0}}}{\partial\hat{X}_{R}}+c_{{S}}{\hat{b}^{\mathfrak{L}}_{{BR},0}}=\dfrac{\partial{\hat{T}^{\mathfrak{L}}}}{\partial\hat{y}}\left(\omega,0\right),\quad \text{on }\hat{X}_{R}=0, \end{gather}$$

(3.11c)$$\begin{gather}{\hat{w}_{{BR},0}}(\hat{X}_{R},\hat{Y}_{B}) \to{\hat{w}_{{B},0}}(\omega,\hat{Y}_{B}), \quad {{\hat{b}^{\mathfrak{L}}_{{BR},0}}}(\hat{X}_{R},\hat{Y}_{B}) \to{\hat{b}^{\mathfrak{L}}_{{BR},0}}(\omega,\hat{Y}_{B}) \quad \text{as }\hat{X}_{R}\to\infty \end{gather}$$

and

(3.12a)$$\begin{gather} {\hat{w}_{{BR},0}} =0,\quad\text{on }\hat{Y}_{B}=0, \end{gather}$$

(3.12b)$$\begin{gather}\dfrac{\partial{\hat{b}^{\mathfrak{L}}_{{BR},0}}}{\partial\hat{Y}_{B}}+c_{{B}}{\hat{b}^{\mathfrak{L}}_{{BR},0}} =0,{\quad} \text{on }\hat{Y}_{B}=0. \end{gather}$$

The solution to problem (3.11) is given by

(3.13a)$$\begin{gather} {\hat{w}_{{BR},0}}(\hat{X}_{R},\hat{Y}_{B}) ={\hat{w}_{{B},0}}(\omega,\hat{Y}_{B})\left(1-{\rm e}^{-\hat{X}_{R}}\right), \end{gather}$$

(3.13b)$$\begin{gather}{{\hat{b}^{\mathfrak{L}}_{{BR},0}}}(\hat{X}_{R},\hat{Y}_{B}) ={\hat{b}^{\mathfrak{L}}_{{B},0}}(\omega,\hat{Y}_{B})+{\hat{w}_{{B},0}}(\omega,\hat{Y}_{B}){\rm e}^{-\hat{X}_{R}}, \end{gather}$$

where the boundary condition (3.11b) is satisfied provided

(3.14)\begin{equation} \left(1+c_{{S}}\right){\hat{w}_{{B},0}}(\omega,\hat{Y}_{B})+c_{{S}}{{\hat{b}^{\mathfrak{L}}_{{B},0}}}(\omega,\hat{Y}_{B})=\dfrac{\partial{\hat{T}^{\mathfrak{L}}}}{\partial\hat{y}}\left(\omega,0\right). \end{equation}

Thus, as promised, we obtain the matching condition (3.10f) required to solve the inner problem near the base.

The solution (3.13) depends only parametrically on $\hat {Y}_{B}$, and therefore cannot in general satisfy the final boundary conditions (3.12). This observation indicates the presence of a further inner layer inside the corner region. We do not present the governing equations in this inner–inner region, as they do not appear to be analytically tractable and do not affect the outer solution, at leading order at least (Cook, Ludford & Walker 1972).

4. Parameter sweeps

4.1. Parameter values and geometrical ratios

We now carry out some parameter sweeps on the full steady trench problem shown in figure 3, keeping all parameter values fixed at the ‘default’ values in table 2, unless specified otherwise. The numerical solutions shown (solid) include comparison with predictions (3.7)–(3.9) provided by the asymptotics (dashed). To facilitate varying the trench geometry while avoiding unphysical parameter regimes, we introduce various physically relevant geometrical ratios: the aspect ratio $r_{{A}} = {H}/{\omega }$; the overfilling ratio $r_{{F}} = {\omega f}/{(\omega -\tau )(H-\varsigma )}$; the sidewall fraction $r_{{S}}={\tau }/{\omega }$; the base fraction $r_{{B}}={\varsigma }/{H}$.

In performing parameter sweeps, we show the effect of changing each parameter on the horizontal and vertical velocity profiles, as well as the horizontal magnetic field profiles. In this context, by ‘horizontal’ and ‘vertical’ profiles, we mean one-dimensional slices through the bisecting lines $\hat {y}=\tfrac {1}{2}(H+\varsigma )$ and $\hat {x}=0$, respectively. The vertical magnetic field profiles are omitted because $\hat {b}(0,\hat {y})$ is identically zero. The temperature is plotted only in cases where it is affected by the parameter being swept over.

4.2. Changing the horizontal dimensions $\tau$ and $\omega$

The effects of changing the sidewall thickness $\tau$ while keeping the total trench width $\omega$ fixed are shown in figure 6(a–c). The horizontal velocity profiles (figure 6a) demonstrate the anticipated flow structure, with the velocity approximately uniform outside Hartmann layers at the trench walls Increasing $\tau$ has the effect of both increasing the net wall conductivity and squeezing the fluid into a narrower channel, and thus decreases the maximum velocity of the fluid enclosed between the walls. For $\tau \gtrsim 0.35$, the Hartmann layers appear to merge and create an approximately parabolic velocity profile. The dashed curves show that the numerical results are in excellent agreement with the asymptotic prediction (3.7a).

Figure 6. Numerical plots of (a,d) horizontal velocity profiles halfway up the trench wall, (b,e) vertical velocity profiles along the $\hat {y}$-axis and (c,f) horizontal magnetic field profiles halfway up the trench wall, for dimensionless sidewall thickness $\tau \in \{ 0.005, 0.05, 0.1,\ldots,0.45, 0.495 \}$ (in (a–c)) and dimensionless trench width $\omega \in \{0.4,0.6,\ldots,1.8,2 \}$ (in (d–f)).

The horizontal magnetic field profiles in figure 6(c) show approximate piecewise linear behaviour, which becomes steeper and more pronounced when the walls become thinner, resulting in a stronger electrical current density passing through them, since their electrical conductivity remains constant. Outside the Hartmann layers, the induced field is very small as expected, given the small value of $\varGamma$ here, namely 0.0299.

Figure 6(b) shows that the boundary layer at the bottom of the trench becomes thinner and less curved as the wall thickness increases, suppressing the velocity jet. The asymptotic prediction (3.7a) works well inside the trench and away from the bottom boundary layer. Above the trench, the velocity seems to approach a constant value somewhat lower than that inside the trench. Making the walls thicker causes the speed of the flow above the walls to decrease because of the increased drag due to the tops of the trench walls. Also, the ‘crossover region’ near the tops of the walls becomes narrower and less curved as the sidewalls are brought closer together.

Figure 6(d–f) shows the effects of varying the dimensionless trench width $\omega$ whilst keeping the sidewall thickness $\tau$ fixed. In figure 6(d), we see that the approximately uniform core velocity predicted by the asymptotics (dashed curves) stays almost constant, as now only the trench width is varying, but not the effective wall conductivity. However, the numerical solutions (solid curves) show that the velocity inside the trench continues to increase slightly with increasing $\omega$, with a central bulge in the velocity profile. Figure 6(e) shows the cause of this discrepancy: as the trench becomes wider, the conducting-wall jet near the base becomes elongated, thus being captured by the central horizontal line $\hat {y} = \tfrac {1}{2}(H+\varsigma )$. The intensity of the magnetic field dipole shown in figure 6(f) does not vary as significantly as in figure 6(c), because the effective conductivity of the sidewalls does not change as they move apart. However, the magnetic field gradient in the bulk changes sign as $\omega$ increases, again because of infiltration by the bottom boundary layer.

4.3. Changing the vertical dimensions $\varsigma$ and $H$

The effects of changing the base thickness $\varsigma$ and wall height $H$ are demonstrated in figure 7, assuming that the free surface thickness $f$ is kept constant. The horizontal velocity plots in figure 6(a,d) show that the bulk velocity approaches a limiting value as either $\varsigma \to 0$ or $H\to \infty$, and in either case the asymptotic and numerical predictions are in generally good agreement. However, as either $\varsigma \to H$ or $H\to \varsigma$, the velocity profile features local maxima near the trench walls, caused by the conducting-base jet intersecting the central horizontal line $\hat {y}=\tfrac {1}{2}(H+\varsigma )$, thus revealing a velocity profile characteristic of conducting-wall MHD duct flow (Hunt & Leibovich 1967). The asymptotic approximation (3.7a) breaks down in either of these limits where the base becomes close to the top of the walls, and the velocity ultimately collapses to zero as the TEMHD effect vanishes along with the sidewalls. Moreover, the boundary layers near the bottom and top of the trench merge and annihilate one another.

Figure 7. Numerical plots of (a,d) horizontal velocity profiles halfway up the trench wall, (b,e) vertical velocity profiles along the $\hat {y}$-axis and (c,f) horizontal magnetic field profiles halfway up the trench wall, for dimensionless base thickness $\varsigma \in \{ 0.01, 0.1, 0.2,\ldots,0.9,0.99 \}$ (in (a–c)) and dimensionless wall height $H \in \{0.2,0.4,\ldots,1.8,2 \}$ (in (d–f)).

In figure 6(b), we see that making the base thicker shifts the velocity profile upwards, slightly strengthening the conducting-base jet and thus increasing the average flow speed, provided it is less than approximately 90 % of the wall height. Increasing $H$ has a similar shifting effect on the velocity profile near the top as shown in figure 6(e). In both cases, the bulk velocity profile demonstrates a uniform shape away from the jet and the top of the wall, and in this region the asymptotics and numerics show good agreement.

In figure 6(c,f), we see that the numerical solutions agree well with the asymptotic prediction that the induced magnetic field is close to zero outside the Hartmann layers, except when the base and top of the walls are brought very close together and the corresponding boundary layers start to affect the outer flow.

4.4. Geometrical parameter space

We now examine the effects on the global flow properties of varying the trench geometry. We consider three averaged quantities that describe the global flow characteristics, namely the average flow rates above and below the trench walls, defined by

(4.1a)$$\begin{gather} \hat{\bar{w}}_{{A}} =\dfrac{1}{2\omega f}\int_{H}^{H+f}\int_{-\omega}^{\omega}\hat{w}(\hat{x},\hat{y})\,{\rm d}\kern0.06em \hat{x}\,{\rm d}\hat{y}, \end{gather}$$

(4.1b)$$\begin{gather}\hat{\bar{w}}_{{B}} =\dfrac{1}{2(\omega-\tau)(H-\varsigma)}\int_{\varsigma}^{H}\int_{-\omega+\tau}^{\omega-\tau}\hat{w}(\hat{x},\hat{y})\,{\rm d}\kern0.06em \hat{x}\,{\rm d}\hat{y}, \end{gather}$$

respectively, and the average surface temperature, defined by

(4.1c)\begin{equation} \hat{\bar{T}}_{{F}}=\int_{-\omega}^{\omega}\hat{T}^{\mathfrak{L}}(\hat{x},H+f)\,{\rm d}\kern0.06em \hat{x}.\end{equation}

For the LiMIT concept to be successful, reasonably high flow rates both above and within the trench should be maintained, whilst the surface temperature should not get too high.

We start by varying the horizontal wall dimensions $\omega$ and $r_{S}=\tau /\omega$ in figure 8(a–c). Figure 8(a) demonstrates that the average velocity above the walls is largely unaffected by the trench width, but increasing the trench wall thickness, and thus the effective conductivity, significantly reduces the flow rate. The velocity inside the trench is maximised for a relatively thin trench with a sidewall fraction of around 10 %, as shown in figure 8(b). When the sidewall fraction is less than approximately 5 %, the average velocities above and below the trench walls take similar values. Otherwise, however, the two average velocities can differ markedly, so that the flow speed inside the trench cannot readily be inferred by observation of the free surface. As figure 8(c) shows, making the trench wider does not greatly affect the average surface temperature. However, making the walls thinner does result in a hotter free surface, because there is a smaller volume of solid wall through which heat can conduct (recall that the solid thermal conductivity is roughly three times that of the liquid). On the other hand, making the trench walls too thick results in very little flow either above or below the walls.

Figure 8. The dependence of the average velocity above and inside the trench ((a,d) and (b,c), respectively) and the average free-surface temperature (c,f) on the trench width $\omega$ and sidewall fraction $r_{S}$ (a–c) and on the trench height $H$ and base fraction $r_{B}$ (d–f), for $\omega \in [0.4,2]$, $H\in \{0.2,2\}$ and $r_{S},r_{B}\in [0.01,0.99]$. Schematics of the trench geometry are shown in the corners of each plot.

Next, we turn our attention to figure 8c(d–f) in figure 8, where the effects of varying the wall height $H$ and the base fraction $r_{B}$ are examined. The velocities above and below the walls appear to show similar distributions. However, figure 8(d) reveals a curved band of optimal values for the flow above the trench, where there is only relatively moderate flow inside the trench. On the left-hand side of this curved region, we see that very little flow is generated by a trench with a very high base and walls resembling short stumps. On the other hand, an optimal value of the flow below the walls in figure 8(e) is achieved when the trench is tall, but with a relatively thin base.

This latter regime is also where the free surface is hottest, as shown in figure 8(f), due to there being less conducting wall for heat flux to pass through. In practice, we suggest that trench dimensions along the optimal band in figure 8(d) represent a good trade-off between increasing the average velocities above and below the walls, and reducing the temperature of the free surface. This region of parameter space corresponds to a configuration with a thick divertor plate and only slightly higher trench walls.

4.5. Changing the film thickness $f$ and the trench angle $\theta$

The effects of varying the film thickness $f$ and the trench slope angle $\theta$ are shown in figure 9. Recall that the trench slope angle $\theta$ enters the problem through the gravitational parameter $\varGamma$ defined by (2.13).

Figure 9. Plots showing the effect of varying (a) the film thickness $f$, (b) the trench slope angle $\theta$, on the vertical velocity profile.

In figure 9(a) we see that overfilling the trench has a marked impact on the nature of the flow above the walls without significantly affecting the flow inside the trench. The velocity above the trench may be much larger than that inside, and approximately quadratic in $\hat {y}$, with its maximum at the free surface, in qualitative agreement with the asymptotic prediction (3.3). Recalling that $\theta = 0.1$ in these solutions, we infer that even moderate trench slope angles can result in large flow rates above the walls. In practice, the intention is to use thermoelectric effects as the principal driver of the flow rather than relying on gravity, not least because similar divertor plates may be used at the top of the tokamak, where gravity acts in the opposite direction. In any case, we hypothesise that a situation where the liquid above the trench flows several times faster than that inside the trench may be susceptible to free-surface instability and thus should be avoided in practice.

As figure 9(b) shows, the velocity profile between the walls is mostly unaffected by the slope angle $\theta$, because TEMHD is the dominant driving mechanism there. On the other hand, increasing $\theta$ increases the flow speed above the trench walls, and in the case where $f = 0.1$ as plotted, while making $\theta < 0$ slows down the flow above the trench and can lead to backflows. The fluid above the walls experience both the drag from the TEMHD-driven flow inside the trench and the gravitational forcing, and when $\theta <0$ these effects act in opposite directions. From the numerical simulations and asymptotic estimates, one can in principle determine how much the trench should be filled for a given slope angle, so that velocities above and below the trench remain approximately equal.

4.6. Varying the applied field $\mathcal {B}^{a}$

In figure 10, we explore the effects of varying the reciprocal Hartmann number $\epsilon$. Recall that $1/\epsilon$ is proportional to the applied magnetic field strength ${\mathcal {B}}^{a}$, while the gravitational parameter $\varGamma$ is independent of ${\mathcal {B}}^{a}$. Provided the material and geometrical properties are kept constant, varying $\epsilon$ is thus mathematically equivalent to varying ${\mathcal {B}}^{a}$. However, when considering the practical implications of the results, we must recall that the velocity scaling (2.12) contains a factor of $1/{\mathcal {B}}^{a}$. We therefore also plot the results using dimensional variables to determine what would be observed in practice as the magnetic field strength is increased.

Figure 10. Plots illustrating the behaviour of dimensionless solutions in (a–c) for $\epsilon \in \{0.5,0.2,0.1,0.05, 0.02,0.01,0.005,0.002,0.001\}$, and the dimensional solutions in (d–f) for ${\mathcal {B}}^{a}\in \{0.1,0.2,0.5,1,2, 5,10\}$ T.

In figure 10(a), we clearly observe the singular behaviour of the horizontal velocity profiles as $\epsilon \to 0$, where the velocity approaches a constant value, determined by the vertical temperature gradient, outside Hartmann layers whose width is of $O(\epsilon )$. In figure 10(b), we observe that the conducting-wall jet at the base becomes narrower and stronger as $\epsilon$ decreases. Meanwhile, the velocity above the trench walls decreases with decreasing $\epsilon$. We also see in figure 10(c) that, as $\epsilon$ decreases, the induced magnetic field approaches a piecewise linear function of $\hat {x}$, negligible inside the fluid but with a finite gradient inside the walls. In all cases, the asymptotic predictions (dashed curves) become increasingly accurate with decreasing $\epsilon$, as expected.

However, when considering dimensional quantities, increasing the magnetic field strength ${\mathcal {B}}^{a}$ has a retarding effect not apparent in figure 10(a–c). In figure 10(d), we see that increasing the magnetic field strength does create boundary layers near the walls, but simultaneously decreases the maximum velocity. This effect is also captured in figure 10(e) where, even though the velocity jet near the trench base still appears, it is much less prominent. Figure 10(f) shows that increasing the applied magnetic field also has a suppressing influence on the induced magnetic field. This perhaps counterintuitive behaviour is consistent with our scaling (2.10c) for the induced magnetic field, which is proportional to $1/{\mathcal {B}}^{a}$.

4.7. Changing $\psi$, the applied field inclination

Although the applied magnetic field is likely to be horizontal in practice, due to the toroidal symmetry of a tokamak, for completeness we illustrate the effect of the applied field angle $\psi$ on the velocity and the induced magnetic field in figures 11 and 12, respectively. The most obvious observation from figures 11(f) and 12(f) is that the flow and magnetic field are reversed, as expected, when $\psi ={\rm \pi}$ and the applied magnetic field is pointing from right to left, rather than left to right.

Figure 11. Surface plots illustrating the effect of changing the applied magnetic field angle $\psi$ on the dimensionless velocity field: (a) $\psi =\frac {1}{6}{\rm \pi}$; (b) $\psi =\frac {1}{3}{\rm \pi}$; (c) $\psi =\frac {1}{2}{\rm \pi}$; (d) $\psi =\frac {2}{3}{\rm \pi}$; (e) $\psi =\frac {5}{6}{\rm \pi}$; (f) $\psi ={\rm \pi}$.

Figure 12. Surface plots illustrating the effect of changing the applied magnetic field angle $\psi$ on the dimensionless induced magnetic field: (a) $\psi =\frac {1}{6}{\rm \pi}$; (b) $\psi =\frac {1}{3}{\rm \pi}$; (c) $\psi =\frac {1}{2}{\rm \pi}$; (d) $\psi =\frac {2}{3}{\rm \pi}$; (e) $\psi =\frac {5}{6}{\rm \pi}$; (f) $\psi ={\rm \pi}$.

In figure 11(a) with $\psi ={\rm \pi} /6$, the typical Hunt velocity jet shifts over to the bottom left-hand corner, and the thickness of the boundary layer near the base decreases. We notice that this conducting jet appears to be stretched out in the direction of the magnetic field, and a small stagnant region also emanates from the top right-hand corner. It appears that the flow profiles are being distorted along dividing subcharacteristics, since the leading-order TEMHD trench problem (see figure 3) becomes hyperbolic in the limit $\epsilon \to 0$. This phenomenon, where the velocity gradient appears to become discontinuous across dividing characteristics passing through the corners, was described by Alty (1971) and Morley & Roberts (1996).

In figure 11(b), with increasing $\psi$, the velocity becomes close to zero near the subcharacteristic passing through the top right-hand corner, in a neighbourhood of which the flow begins to reverse. In figure 11(c), with $\psi ={\rm \pi} /2$, there is a perfectly symmetrical velocity dipole with the liquid in the right-hand half of the trench flowing at an equal and opposite speed to that in the left-hand half, and the average velocity is equal to zero. The same behaviour as in figure 11(a,b) is demonstrated in figure 11(d,e), albeit with reversed direction. In any case, it appears that the net velocity in the trench is maximised when the applied magnetic field is orthogonal to the temperature gradient, as expected.

Turning our attention to figure 12, we notice that the induced field is close to zero outside a triangular region near the bottom right-hand corner in figure 12(a,b). The dividing line follows a subcharacteristic inclined at an angle $\psi$ to the horizontal, which encloses a region of high induced field strength, counteracting the dipole formed by the rest of the magnetic field profile. In figure 12(c), the applied field is parallel to the trench walls, where the induced magnetic field exhibits side layers of thickness of $O(\epsilon ^{1/2})$. These layers have a stronger electrical current through the walls, and hence a stronger Lorentz force which drives the fluid at a higher speed in both the left- and right-hand sides of the trench, although in opposite directions.

5. Conclusions

In this paper, we performed a hybrid numerical/asymptotic analysis of steady unidirectional TEMHD flow inside a trench akin to those seen in LiMIT. We obtained asymptotic solutions for the velocity, magnetic field and temperature in the limit of large Hartmann number and thin trench walls, demonstrating how temperature gradients drive fluid flow. We obtained numerical solutions to the problem in different parameter cases corresponding to both laboratory and fusion-relevant applications, and these demonstrated very good agreement with the asymptotic predictions.

Several interesting future avenues for research have arisen as a result of the modelling work we have presented. One could attempt to model the cooling tubes beneath the trench in more detail, e.g. by replacing the Dirichlet boundary condition on the temperature on the base with a Robin boundary condition. We considered only steady solutions in this paper, but unsteady effects may be important in practice, for example during start-up of a tokamak, when the applied magnetic field is ramped up from zero to a several teslas over a number of minutes. Also of practical interest is the effect on the flow in the trench of plasma instabilities such as edge-localised modes, which can cause large localised transient spikes in heat flux and/or applied field, as well as significant current into the free surface from the plasma. The practical applicability of our work could be improved by including temperature dependence in the material properties of lithium which, as noted in § 2.3, can vary appreciably (though not dramatically) over relevant temperature ranges.

We have focused on rectangular trenches in this modelling work, because they are the most straightforward to design and manufacture. However, three-dimensional printing capabilities can in principle produce novel trench shapes, including cross-sections that vary along the trench. It would be interesting to pose a shape optimisation problem, for example, to maximise the average flow speed inside the trench, given engineering constraints on the wall and base properties.

We have assumed that the physically imposed properties, including the impinging heat flux and applied magnetic field, are all uniform in the flow direction along the trench. In practice, the impinging heat flux is likely to be concentrated around a strip near the centre of the trench (described by a Gaussian in Xu et al. (2014)), and the applied magnetic field is also likely to vary slightly along the trench in inverse proportion to distance from the symmetry axis. The next step is to relax the assumption of unidirectional flow and take into account variations along the trench. Since the trench can be assumed to be much longer than it is wide, a slowly varying version of the model can be derived by considering a sequence of two-dimensional slices, each of which resembles the purely two-dimensional problem analysed in this paper. The aim would be to obtain a shallow-water type theory (as in, e.g. Fiflis et al. 2016; Lunz & Howell 2019; Lunz 2020), in which the applied field, heat flux and film height all vary parametrically with distance $z$ along the trench.