Effect of insulator end cap thickness on time-dependent Hartmann flow in a rotating mirror

Abstract

We present a framework for analysing plasma flow in a rotating mirror. By making a series of physical assumptions, we reduce the magnetohydrodynamic (MHD) equations in a three-dimensional cylindrical system to a one-dimensional system in a shallow, cuboidal channel within a transverse magnetic field, similar to the Hartmann flow in ducts. We then solve the system both numerically and analytically for a range of values of the Hartmann number and calculate the dependence of the plasma flow speed on the thickness of the insulating end cap. We observe that the mean flow overshoots and decelerates before achieving a steady-state value, a phenomenon that the analytical model cannot capture. This overshoot is directly proportional to the thickness of the insulating end cap and the external electric field, with a weak dependence on the external magnetic field. Our simplified model can act as a benchmark for future simulations of the supersonic mirror device CMFX (centrifugal magnetic fusion experiment), which will employ more sophisticated physics and realistic magnetic field geometries.

1. Introduction

Rotating mirrors present a novel research direction in magnetic confinement fusion (Post 1987). This approach typically employs a large background magnetic field combined with an externally supplied current, generating plasma rotation at supersonic speeds. Compared with more complex magnetic configurations, this method offers a potentially simpler design for fusion reactors. Previous experiments have demonstrated the feasibility of supersonic rotating mirrors for nuclear fusion, as reported in the literature (Reid et al. 2014; Ellis et al. 2005).

The CMFX (centrifugal magnetic fusion experiment), funded by the ARPA-E BETHE (Breakthroughs Enabling THermonuclear-fusion Energy) program (Romero-Talamás et al. 2021), aims to achieve a significant fusion triple product — with the product of plasma density, temperature and confinement time greater than $10^{20}\ ({\rm keV}\ {\rm s})\ {\rm m}^{-3}$. The success of this endeavour would mark a substantial advance in the development of practical fusion energy.

The effectiveness of this approach may depend on the choice of materials for the plasma end caps. These disk-shaped components are critical for maintaining axial plasma confinement and enabling magnetic field rotation. Therefore, the material characteristics of the end caps are important for the efficiency and viability of the fusion process.

In this study, we examine how the electrical properties and thickness of the end-cap materials affect plasma flow. To simplify the analysis, we assume a model that transforms the problem from laminar flow in a cylindrical geometry to laminar flow in a shallow, cuboidal channel within a transverse magnetic field, similar to Hartmann flow in ducts (Hartmann & Lazarus 1937a,b).

Previous research, including the work by Huang (2004), has studied the impact of the Hartmann boundary layer on plasma flow and magnetohydrodynamic (MHD) stability. Hassam & Huang (2019) have also explored the effect of perfectly conducting walls and the emergence of small-scale physical oscillations in a one-dimensional (1-D) MHD channel flow. Recent work by Zhang et al. (2022) has analysed flow regimes in rectangular channels over a wide range of magnetic and fluid Reynolds numbers. However, a comprehensive solution that encompasses end-cap fields and their impact on flow remains elusive. Our simplified model aims to establish an analytical relationship, providing a benchmark for future simulations employing more sophisticated physics and realistic magnetic field geometries. The principles derived here are adaptable to broader scenarios using complex numerical or analytical models. Figure 1 illustrates the set-up of the supersonic mirror device.

Figure 1. A simplified version of the supersonic mirror in a cylindrical coordinate system $(r, \theta, z)$. The background field $B_0 \hat{\boldsymbol{z}}$ is generated by external magnets (not shown). The device comprises an inner electrode, which is a solid conducting rod of radius $R_0$, and an outer conducting shell of radius $R_1$. On both ends, the grey part of the end cap denotes an insulator, whereas the white part denotes an imperfect conductor. Plasma remains in the annular region between the two electrodes. Due to the external potential difference between the electrodes, the radial current $j_{\mathrm {ext}} \hat{\boldsymbol{r}}$ flows through the plasma; coupled with the background field $\boldsymbol{B}_0 = B_0 \hat{\boldsymbol{z}}$ causes the plasma to rotate in the azimuthal ($\hat{\boldsymbol{\theta }}$) direction.

To describe the evolution of the flow under an external azimuthal force, we use the MHD model. In § 2, we explain the three-dimensional (3-D) MHD model and the assumptions used to simplify it to a 1-D model. We then solve the 1-D model in all three media: imperfect conductor, insulator and plasma. We present the evolution of the time-dependent solution at different values of the Hartmann number. In § 3, we demonstrate how the mean plasma flow overshoots and decelerates, before reaching a steady-state value, and how the thickness of the insulator and the external fields affect the mean flow. In § 4, we present our conclusions.

2. Solving the MHD model: theory and analysis

In this section, we present a theoretical framework for analysing plasma behaviour in a rotating mirror. We adopt the incompressible MHD model, which assumes a constant plasma density, to describe plasma dynamics. To facilitate our analysis, we initially simplify the model to a 1-D representation based on a set of well-defined assumptions appropriate to the CMFX device. The resulting simplified 1-D model is then solved numerically, and further simplified and solved analytically to provide insights into the plasma behaviour under the influence of the rotating magnetic field in the mirror set-up.

For plasma, we start with the incompressible MHD (Freidberg 2014) model

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

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

(2.3)\begin{gather}\left( \frac{\partial}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} \right) \frac{p}{\rho^{\gamma}} = 0, \end{gather}

(2.4)\begin{gather}-\frac{\partial \boldsymbol{B}}{\partial t} = \boldsymbol{\nabla} \times \boldsymbol{E}, \end{gather}

(2.5)\begin{gather}\mu_{p} \boldsymbol{j} = \boldsymbol{\nabla} \times \boldsymbol{B} - \frac{1}{c^2}\frac{\partial \boldsymbol{E}}{\partial t}, \end{gather}

(2.6)\begin{gather}\boldsymbol{E} ={-}\boldsymbol{u}\times \boldsymbol{B} +\eta \boldsymbol{j}, \end{gather}

(2.7)\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{B} = 0, \end{gather}

where $\rho$ is the plasma density, $\boldsymbol{u}$ is the plasma flow velocity, $\boldsymbol{B} = \boldsymbol{B}_0 + \boldsymbol{B}_1$ is the total magnetic field, $\boldsymbol{B}_0$ is the static background field generated by the external coils, $\boldsymbol{B}_1$ is the time-dependent field generated by the plasma, $\boldsymbol{j}$ is plasma-generated current density, $p$ is the plasma pressure and $\boldsymbol{F} = \boldsymbol{j}_{\mathrm {ext}} \times \boldsymbol{B}_0$ is the external force on the plasma generated by the current $\boldsymbol{j}_{\mathrm {ext}}$ generated due to the externally applied potential difference. The constants $\mu$, $\eta$ and $\mu _{p}$ are the dynamic viscosity, magnetic field diffusivity and magnetic permeability of the plasma, respectively. In equilibrium, we assume a uniformly magnetized static plasma with $p = p_0, \rho = \rho _0, \boldsymbol{u} = 0, \boldsymbol{B} = \boldsymbol{B}_0$.

Similarly, inside the plasma wall materials described in figure 1, we solve Maxwell's equations

(2.8)\begin{gather} -\frac{\partial \boldsymbol{B}}{\partial t} = \boldsymbol{\nabla} \times \boldsymbol{E}, \end{gather}

(2.9)\begin{gather}\mu_{c} \boldsymbol{j} = \boldsymbol{\nabla} \times \boldsymbol{B} - \frac{1}{c^2}\frac{\partial \boldsymbol{E}}{\partial t}, \end{gather}

where $\mu _{c}$ is the magnetic field permeability of the imperfect conductor. By definition, inside an insulator, the current

(2.10)\begin{equation} \boldsymbol{j} = 0,\end{equation}

and for an imperfect conductor that satisfies Ohm's law,

(2.11)\begin{equation} \boldsymbol{j} = \sigma_{c} \boldsymbol{E},\end{equation}

where $\sigma _{c}$ denotes the conductivity of the imperfect conductor. The last two equations are the result of the electrical properties of the end-cap materials. The electric and magnetic fields in all three regions: imperfect conductor, insulator and plasma, must be continuous, since there is no charge accumulation or surface current on the interfaces. We begin by normalizing the model equations and describing the various assumptions we make to simplify the 3-D MHD model.

2.1. Normalization and cylinder-to-slab transformation

In this section, we simplify the model by normalizing it as follows:

(2.12a–g)\begin{align} \tilde{\boldsymbol{B}} = \frac{\boldsymbol{B}}{B_0}, \quad \tilde{\boldsymbol{u}} = \frac{\boldsymbol{u}}{v_{A}}, \quad v_{A} = \frac{B_0}{\sqrt{\mu_{p} \rho_0}}, \quad \tilde{\boldsymbol{E}} = \frac{\boldsymbol{E}}{v_{A} B_0}, \quad \tilde{t} = \frac{t v_{A}}{L}, \quad \tilde{z} = \frac{z}{L}, \quad \tilde{r} = \frac{r}{R_0}, \end{align}

where $B_0$ is the magnetic field due to the coils, $v_A$ is the Alfvén speed, $L$ is axial length and $R_0$ is the radial size of the device. Solving (2.1)–(2.11) analytically in a 3-D cylindrical geometry is not possible. To reduce the complexity of this model, we adopt an asymptotic ordering that differentiates between multiple scales and various quantities by employing the small parameter $\epsilon$:

(2.13a–c)\begin{gather} \tilde{u}_r, \tilde{u}_z \sim \epsilon \tilde{u}_{\theta} \sim \epsilon^2, \quad \tilde{B}_r, \tilde{B}_z \sim \epsilon \tilde{B}_{\theta} \sim \epsilon \tilde{B}_{{v}0} \sim \epsilon^2, \quad \tilde{E}_{r} \sim \epsilon \tilde{E}_{{v}0} \sim \epsilon^2, \end{gather}

(2.14a–d)\begin{gather}\hat{\boldsymbol{z}}\boldsymbol{\cdot} \boldsymbol{\nabla} \sim 1/L \sim (\hat{\boldsymbol{r}}, \hat{\boldsymbol{\theta}})\boldsymbol{\cdot} \boldsymbol{\nabla} \sim 1/R_0, \quad (R_1-R_0)/R_0 \sim \epsilon,\quad \beta \sim \epsilon^2,\quad \lvert \boldsymbol{F} \rvert \sim \epsilon. \end{gather}

In these orderings, the subscript letter $r, \theta$ or $z$ denotes the component of a quantity, while the subscript ${v}$ denotes the external fields in the vacuum region, and $\beta = 2 \mu _o p/B_0^2$ is the ratio of plasma pressure to magnetic pressure. The assumptions $\hat{r}\boldsymbol{\cdot } \boldsymbol{\nabla } \sim 1/R$ and $(R_1-R_0)/R_0 \sim \epsilon$ imply a shallow device with small spatial gradient that enables us to eliminate any radial variation in our annular domain. By making assumptions that remove radial variation and impose azimuthal symmetry, the problem in three dimensions can be effectively transformed into a configuration that resembles a 1-D slab or rectangular channel.

We also assume that the azimuthal flow dominates the other components of the flows and that the plasma-generated azimuthal magnetic field is greater than the rest of the components of the field. Such behaviour of the plasma leads to a set of equations where all the nonlinear terms involving interaction between the flow and fields completely vanish, transforming our model into a set of linear equations. Another important assumption to be used is that the power supply generates electric and magnetic fields $E_{v}$ and $B_{v}$ that are of the same order as the plasma-generated fields. These fields will be an important part of our conclusion in the last section.

The conceptualization and simplification of the 3-D domain is further elucidated in figure 2. In the slab limit, we substitute the subscripts $r, \theta$ and $z$ with $x, y$ and $z$, respectively. Note that such equations only involve the radial component of the electric field $\tilde {E}_r(z)$ (or $\tilde {E}_x(z)$ in a slab), the azimuthal components of the magnetic field $\tilde {B}_{\theta }(z)$ (or $\tilde {B}_y(z)$ in a slab) and plasma flow $\tilde {u}_{\theta }(z)$ (or $\tilde {u}_y(z)$ in a slab) throughout this paper. To further simplify the notation, we omit $\,\tilde{}\,$ in normalized quantities. In the following section, we begin by solving Maxwell's equation in an imperfectly conducting end cap.

Figure 2. (a) Rotating mirror set-up of a large aspect ratio mirror device in the form of an annular cylinder. (b) A simplified slab (rectangular channel) model showing a cut section of the mirror in panel (a). The region outside the end caps is treated as a vacuum with external vacuum fields $E_{v}, B_{v}$. The equilibrium field $\boldsymbol{B}_0$ points in the $z$-direction whereas the plasma flow $u_y$ and plasma-generated magnetic field $B_y$ are in the azimuthal $y$-direction.

Throughout the paper, we impose even parity on the flow and electric field, and odd parity on the magnetic field about the $z=0$ plane,

(2.15a–c)\begin{equation} u_y(z) = u_y({-}z), \quad E_x(z) = E_x({-}z), \quad B_y(z) ={-}B_y({-}z). \end{equation}

2.2. Time-dependent fields in the vacuum region

We start by assuming the fields outside the device as

(2.16)\begin{equation} \left.\begin{array}{c} E_{v} = E_{{v}0} [1 - \exp({-}c_0 t)],\\ B_{v} = B_{{v}0} \operatorname{sgn}(z) [1 - \exp({-}c_0 t)], \end{array}\right\}\end{equation}

where $B_{{v}0}, E_{{v}0}$ are the steady-state vacuum electric and magnetic fields due to external sources such as the power supply, and we can use the outside fields as boundary conditions to calculate the coefficients in the rest of the media.

The external electric and magnetic fields arise from the voltage difference and the current flowing through the power lines. The time-dependent form of the fields is similar to that of the current and potential of a resistor circuit. The parameter $c_0$ corresponds to the time taken for a typical experimental discharge, typically hundreds of Alfvén time periods.

These fields will be used as boundary conditions for the time-dependent fields inside the imperfect conductor in the next section.

2.3. Time-dependent solution inside the imperfect conductor

Our analysis begins by solving Maxwell's equations within the imperfect conductor. We employ the magnetic and electric fields in the vacuum region as boundary conditions to fully determine the field characteristics up to the interface between the imperfect conductor and the insulator. We use the induction equation,

(2.17)\begin{equation} \frac{\partial B_{y}}{\partial t} ={-}\frac{\partial E_x}{\partial z},\end{equation}

and Ampere's law,

(2.18)\begin{equation} E_x = \frac{1}{{S}_{c}} j_{x} ={-}\frac{1}{{S}_{c}}\partial_z B_{y}, \quad {S}_{c} = v_{A} L \sigma_{c} \mu_{c}, \end{equation}

where $\mu _{c}$ is the magnetic field permeability of the imperfect conductor, $\sigma _{c}$ is the electric conductivity and ${S}_{c}$ is the Lundquist number which is a measure of magnetic diffusion inside a material. We have neglected the displacement current term because it scales as $(v_{A}/c)^2 \ll 1$ which is small compared with the Alfvénic time scale $t \sim 1$ of the problem.

The general fields inside the imperfect conductor can be written as

(2.19)\begin{equation} \left.\begin{array}{c} \displaystyle E_x ={-}\dfrac{D_3}{{S}_{c}} + \sqrt{\dfrac{c_0}{{S}_{c}}} \exp({-}c_0 t)[ D_5 \cos(\sqrt{{S}_{c} c_0}z) + D_6 \operatorname{sgn}(z) \sin(\sqrt{{S}_{c} c_0}z) ],\\ B_y = D_3 z + D_4 - \exp({-}c_0 t)[ D_5 \sin(\sqrt{{S}_{c} c_0} z) - D_6 \operatorname{sgn}(z) \cos(\sqrt{{S}_{c} c_0} z) ]. \end{array}\right\} \end{equation}

Matching the tangential fields $E_x$ and $B_y$ with the outside fields $E_{v}$ and $B_{v}$, respectively, we obtain and simplify the fields inside the conductor as

(2.20)\begin{align} E_x & = E_{{v}0} - \exp({-}c_0 t) \left[E_{{v}0} \cos(\sqrt{{S}_{c}c_0}(z-\operatorname{sgn}(z)(0.5+d_{i} +d_{c})))\vphantom{\frac{c_0}{{S}_{c}}}\right.\nonumber\\ & \quad \left.+\sqrt{\frac{c_0}{{S}_{c}}} B_{{v}0} \operatorname{sgn}(z) \sin(\sqrt{{S}_{c}c_0}(z-\operatorname{sgn}(z)(0.5+d_{i}+d_{c})))\right], \end{align}

(2.21)\begin{align} B_y & ={-}{S}_{c} E_{{v}0} (z - \operatorname{sgn}(z)(0.5 +d_{i} +d_{c})) + B_{{v}0} \operatorname{sgn}(z) \nonumber\\ & \quad +\exp({-}c_0 t)\left[\sqrt{\frac{{S}_{c}}{c_0}}E_{{v}0} \sin(\sqrt{{S}_{c}c_0}(z-\operatorname{sgn}(z)(0.5+d_{i}+d_{c})))\right.\nonumber\\ & \quad - \left. \vphantom{\sqrt{\frac{{S}_{c}}{c_0}}}B_{{v}0} \operatorname{sgn}(z) \cos(\sqrt{{S}_{c}c_0}(z-\operatorname{sgn}(z)(0.5+d_{i}+d_{c})))\right]. \end{align}

Equations (2.20) and (2.21) represent the electromagnetic fields within the imperfect conductor. We note that the expressions include a $\operatorname {sgn}$ function that is discontinuous at $z = 0$. However, since the domain of the imperfect conductor does not include the point $z = 0$, the fields remain continuous inside the conductor and the use of a $\operatorname {sgn}$ function is valid. Note that these solutions are only valid for finite values of $S_{c}$. These results will serve as boundary conditions for the subsequent analysis of the fields within the insulator, as explained in the following section.

2.4. Time-dependent solution inside the insulator

After solving for the field inside the imperfect conductor, we solve for the electromagnetic fields inside the insulator end cap. The fields inside the perfect insulator must satisfy the induction equation

(2.22)\begin{equation} \frac{\partial B_{y}}{\partial t} ={-}\frac{\partial E_x}{\partial z}, \end{equation}

and Ampere's law

(2.23)\begin{equation} j_{x} ={-}\partial_z B_{y} = 0. \end{equation}

In the insulator, the magnetic field is spatially uniform, whereas the electric field varies in response to the time-dependent magnetic field. The tangential components of the magnetic $B_{y}$ and electric $E_{x}$ fields are

(2.24)\begin{align} E_{x} & = E_{{v}0} + \exp({-}c_0 t) \left\{\left[{-}E_{{v}0} \cos(\sqrt{{S}_{c} c_0}d_{c}) + \sqrt{\frac{c_0}{{S}_{c}}} B_{{v}0} \sin(\sqrt{{S}_{c} c_0}d_{c}) \right]\right.\nonumber\\ & \quad - c_0 (\lvert z \rvert - 0.5 - d_{i})\left.\left[\sqrt{\frac{{S}_{c}}{c_0}} E_{{v}0} \sin(\sqrt{{S}_{c} c_0}\, d_{c}) + B_{{v}0} \cos(\sqrt{{S}_{c} c_0}\, d_{c}) \right] \right\}, \end{align}

(2.25)\begin{align} B_y & = \operatorname{sgn}(z)({S}_{c} E_{{v}0} d_{c} + B_{{v}0}) \nonumber\\ & \quad - \operatorname{sgn}(z) \exp({-}c_0 t)\left[\sqrt{\frac{{S}_{c}}{c_0}} E_{{v}0} \sin(\sqrt{{S}_{c} c_0} d_{c}) + B_{{v}0} \cos(\sqrt{{S}_{c} c_0} d_{c}) \right]. \end{align}

The magnetic field stays the same as the boundary between the imperfect conductor and the insulator, but the electric field changes linearly with the insulator thickness due to the induction equation.

Similar to the previous section, we have neglected the displacement current term because it scales as $(v_{A}/c)^2 \ll 1$ which is small compared with the Alfvénic time scale $t \sim 1$ of the problem. However, it is important to note that close to $t=0$, the fields $E_{v} = B_{v} = 0$, whereas the fields inside the imperfect conductor and insulator have finite values which violates causality. This violation of causality arises from the omission of the displacement current term. In Appendix A, we show that incorporating the displacement current resolves this issue, and neglecting the displacement current does not influence the long-time solution or any outcomes of this study.

The effect of these external fields on the plasma dynamics is explored in § 3.

2.5. Time-dependent solution inside the plasma

After solving the fields in the imperfect conductor and the insulator, we solve the equation inside the plasma. This will completely define the solution in all three media.

The lowest order equations correspond to the equilibrium: $p = p_0, \rho = \rho _0, \boldsymbol{B} = \boldsymbol{B}_0$ and $\boldsymbol{u} = 0$, and time-dependent quantities arise only at first order. The plasma satisfies (2.1)–(2.7) which, using the orderings defined in (2.13a–c) and (2.14a–d) for a shallow-channel case with dominant azimuthal plasma flow and magnetic fields, can be reduced to a set of coupled 1-D partial differential equations:

(2.26)\begin{gather} \partial_t u_{y} = \partial_z B_{y}+ \frac{1}{Re} \partial^2_{z} u_{y} + F_0(1 - \exp({-}c_0 t)), \quad {Re} = \frac{\rho_0 v_{A} L}{\mu}, \end{gather}

(2.27)\begin{gather}\partial_t B_{y} = \partial_{z} u_{y} + \frac{1}{Rm} \partial^2_{z}B_{y}, \quad {Rm} = \frac{v_{A} L \mu_{p}}{\eta}, \end{gather}

(2.28)\begin{gather}E_{x} ={-}u_{y} - \frac{1}{Rm} \partial_{z} B_{y}, \end{gather}

subject to time-dependent boundary conditions on the insulator–plasma interface. The flow $u_{y} = u_{y}(z)$ and field $B_{y} = B_{y}(z)$, and the external forcing term $F_y = F_0 (1 - \exp (-c_0 t))$. Note that the form of the forcing function is chosen to be similar to that of the start-up stage of a rotating mirror. The forcing and external fields $E_{v}, B_{v}$ are generated by discharging a voltage source, typically an array of capacitors, and therefore have the same time evolution, i.e. $(1-\exp (-c_0t))$. For simplicity, we assume that the electrical properties of the plasma do not change during operation.

Coupled equations (2.26) and (2.27) describe the evolution of the azimuthal plasma flow and magnetic field, while (2.28) determines the electric field. We also introduce a new dimensionless parameter, the Hartmann number, defined as ${Ha} \equiv \sqrt {{Re}{Rm}}$, which will be used in the subsequent analysis. To further understand this model, we numerically solve it and provide analysis of the solution in the next section.

3. Results from the simplified 1-D Hartmann model

In this section, we solve the simplified 1-D Hartmann flow model and present a detailed analysis of our results. We will then explain how plasma flow depends on the insulator wall thickness.

We numerically solve (2.26) and (2.27) using the spectral numerical code Dedalus (Burns et al. 2020) implementing time-dependent boundary conditions on electric and magnetic fields obtained from the solution in the insulating end caps, given by (2.25) and (2.24) at $z = \pm 0.5$. We apply a no-slip boundary condition to the flow. We use a Chebyshev grid with $n_z = 128$ grid points to resolve rapidly varying spatial features, such as the boundary layer. The Chebyshev polynomial is well suited for this problem, as the grid points are mostly clustered at the edge of the domain, which efficiently resolves the boundary layer. For the time-stepper, we use a second-order backward difference (SBDF2) time-stepping routine. An implicit method such as SBDF2 helps us avoid the unphysical oscillations associated with a stiff system such as this. We solve the model for two values of the Hartmann number and present the results in figure 3.

Figure 3. Plasma flow and magnetic field profiles as functions of cylindrical distance $z$ at three different time values with (a) ${Re} = 50, {Ha} = 50$ and (b) ${Re}= 50, {Ha} = 500$. For this solution, we have chosen $F_0 = 0.5, c_0 = 10^{-2}, {S}_{c} = 10^{4}, d_{c} = 5 \times 10^{-3}, d_{i} = 10^{-2}, E_{{v}0} = -1, B_{{v}0} = 1$. Due to non-ideal effects, the plasma forms a sharp boundary layer near the insulating end caps. Note that the fields and flows have been scaled by $1/\epsilon$ to avoid adding factors of $\epsilon$ to all quantities on the $y$ axis.

The plasma behaviour is governed by ideal MHD in the core, with non-ideal effects dominating the boundary layer. The forcing term is balanced by the curvature of the magnetic field in the core region (around $z=0$), which is frozen in the fluid, flowing with a nearly uniform speed. However, in the boundary layer, the plasma can move relative to the magnetic field lines, which allows it to slip over the boundary. The thickness of the boundary layer is proportional to $1/\sqrt {Ha}$ – a system with a high ${Ha}$ has a thinner boundary layer. Additionally, the core flow speed appears to be unaffected by the Hartmann number.

Since ${Ha} \gg 1$ for supersonic rotating plasmas, the presence of a boundary layer introduces a new length scale that we can use to solve the problem analytically (Bender & Orszag 2013). This process, the analytical solution, and its comparison with the numerical solution are presented in Appendix B.

An important feature that we observe is the overshooting of the mean core plasma flow

(3.1)\begin{equation} \bar{u}_y = 2 \int_{{-}0.25}^{0.25} \, {\rm d} z u_y, \end{equation}

beyond its steady-state value. We demonstrate this phenomenon in figure 4. As the boundary conditions are time dependent, the system adjusts to the new boundary values through the generation of Alfvén waves. However, time-dependent boundary conditions lead to an overshooting of the mean flow before it decelerates to a steady-state value. We observe this phenomenon in our model over a wide range of input parameters, and overshooting continues to occur in systems with high Hartmann numbers, regardless of whether the mean flow is calculated over the entire domain or just the core. Therefore, for the CMFX device, it is crucial to ensure a sub-Alfvénic mean flow, as approaching Alfvénic speeds can induce various instabilities and reduce the confinement time (Teodorescu et al. 2010).

Figure 4. Mean core plasma flow speed $\bar {u}_y = 2\int _{-0.25}^{0.25} \, \textrm{d} z u_y$ as a function of time $t$ for (a) ${Ha} = 50$ and (b) ${Ha} = 500$, and compare the numerical and the simplified analytical solutions. The inset shows the initial part of the numerical and analytical solutions. The solutions agree well, but only close to the steady state. However, the analytical model cannot capture the Alfvénic dynamics in the beginning, or the overshooting and subsequent deceleration (shaded region) of the flow. The parameters used for these figures are the same as those used in figure 3.

Finally, for the same parameter values used in figure 3, we also plot the dependence of the mean flow normalized by the external electric field $\bar {u}_y/E_{{v}0}$ as a function of the thickness of the insulator $d_{i}$ in figure 5 for two different values of the external magnetic field $B_{{v}0}$. We find that the magnitude of the plasma flow reduces with increasing insulator end cap thickness and is sensitive to the values of the external electric and magnetic fields, especially during the ramp-up phase of the device. However, the dependence of the mean flow on the insulator thickness is strongly affected by the external vacuum electric field $E_{{v}0}$. Therefore, to avoid overshooting of the plasma flow from the target value, one must carefully choose the insulator end cap thickness. As the system approaches steady state, the dependence of the mean flow on the insulator goes down and at steady state, they are completely independent.

Figure 5. Dependence of the mean flow $\bar {u}_y$ normalized by the vacuum electric field $E_{{v}0}$ against the insulator end cap thickness $d_i$ at different times and for different values of the external magnetic field. (a) ${Ha} = 500, B_{{v}0} = 1$ and (b) ${Ha} = 500, B_{{v}0} = 10$. The flow velocity has a strong dependence on the thickness of the insulator at the beginning. We see that the flow transitions from overshooting to undershooting the steady-state value with increasing insulator thickness, around $d_i = 2$ in panel (a) and $d_i = 2.5$ in panel (b). Hence, to avoid overshooting the flow velocity, one must choose a thicker insulator end cap. Note that all plasma-generated flows and fields have been scaled up by $1/\epsilon$ because of their small size compared with the respective background quantities.

To accurately determine the flow magnitude, we must accurately and self-consistently calculate the fields $E_{{v}0}$ and $B_{{v}0}$. However, in general, the electric field $E_{{v}0} = E_{{v}0}(V_0, R_0, R_1)$ is a complex function of $V_0$, the potential difference across the electrodes, and $R_0$, $R_1$, the radial positions of the electrodes. The current $I_{{v}0} = I_{{v}0}(V_0, \sigma _{p}, \mu _{p})$ in the external wires is an unknown function of $V_0$ and the electrical properties of the plasma. The magnetic field between the circular electrodes will be $B_{{v}0} = B_{{v}0}(I_{{v}0}, R_0, R_1, \sigma _{p}, \mu _{p})$. Hence, the values of the fields depend on the electrical properties of the plasma and how they change over time. To better understand it, one would require a nonlinear, multi-physics solver.

4. Conclusions

We simplified the 3-D MHD model by transforming the equations describing a supersonic rotating plasma in cylindrical geometry into a 1-D slab model by making a set of assumptions appropriate for the CMFX. We then numerically and analytically solved our simplified model and calculated the dependence of the plasma flow speed on the thickness of the insulating end cap. Our analysis shows that the mean flow speed of the plasma is linearly reduced with the thickness of the insulator wall. This reduction in mean flow is stronger for large external electric fields. Moreover, to avoid the plasma flow from overshooting beyond its steady-state value, one must carefully choose the insulator end cap thickness. Therefore, the design of the CMFX device must carefully take into account the electrical properties of the end caps and the external electric and magnetic fields.

This work opens many avenues for future research. A key step forward is to extend our technique to calculate the effect of the insulator thickness on the plasma flow by solving the 3-D equations in all the different materials simultaneously using a multiphysics solver. Our simplified model can act as a benchmarking tool for these sophisticated solvers. The analysis could also be repeated with different boundary conditions, realistic magnetic field geometry, and with temperature dependence, gyroviscosity, Hall effects and kinetic effects.

Furthermore, our analysis has potential implications for ongoing liquid lithium-metal experiments, as detailed by Saenz et al. (2022), which are crucial for liquid divertor concepts in future tokamak fusion power plants. The similarity of these experiments with our 1-D slab model solution further underscores the relevance of our findings in practical fusion applications.