Magnetohydrodynamics Solver for a Two-Phase Free Surface Flow Developed in OpenFOAM

Abstract

A magnetohydrodynamics solver (“mhdCompressibleInterFoam”) has been developed for a compressible two-phase flow with a free surface by extending “compressibleInterFoam” solver within OpenFOAM suite. The primary goal is to develop a tool to simulate compression of magnetic fields in vacuum and simplified magnetized plasma targets by imploding rotating liquid metal liners in the context of a Magnetized Target Fusion (MTF) concept in pursuit by General Fusion Inc. At present, the solver is limited to axisymmetric problems and the magnetic field evolution is solved in terms of toroidal field component and poloidal flux functions. The solver has been validated and verified using a number of test cases for which analytical or other numerical solutions are provided. Those tests cases include: (i) compression of toroidal and poloidal magnetic fields in vacuum and cylindrical geometry, (ii) axisymmetric annular Hartmann flow, and (iii) compression of magnetized target initialized with a Grad–Shafranov equilibrium state in a cylindrical geometry. A methodology to incorporate conductive solid regions into simulation has also been developed. Capability of the code is demonstrated by simulating a complex case of compressing a magnetized target, which is injected during implosion of a rotating liquid metal liner with an initially soaked poloidal magnetic field. An application of the solver to simulate compression of a magnetized target in a geometry and parameters relevant to the Fusion Demonstration Plant (FDP) being developed by General Fusion Inc. is also demonstrated.

1. Introduction

Development of an MHD solver capable of handling compressible fluids with a free-surface has been prompted by simulation needs in General Fusion Inc. to advance development of the Fusion Demonstration Plant (FDP) [1]. The magnetized Target Fusion (MTF) scheme in pursuit by General Fusion Inc. relies on achieving nuclear fusion conditions during compression of magnetized plasma by imploding rotating liquid metal liner [1]. A schematic of the system is shown in Figure 1. Magnetized plasma is injected into a vacuum cavity inside a rotating liquid metal liner. Pistons, driven by high pressure gas, push on the outer surface of the rotating metal liner causing it to implode and compress a magnetized plasma target. In the vertical direction, the pistons are arranged into layers, and the trajectory of pistons can be varied from layer to layer. This allows to shape the liner as it implodes, starting up from a cylindrical one and ending up with more like a spherical one during the implosion. Simulating the entire system is challenging, as it requires combining the dynamics of magnetized plasma, the dynamics of the imploding liquid metal liner, and, to some extent, of the piston system. One of the main issues in simulating plasma and liner dynamics simultaneously is that they operate at different times scales; plasma simulations usually require a much smaller time step (in the order of $\mathrm{dt}=1\times {10}^{-11}\mathrm{s}$) compared with the hydrodynamics of the liner (in the order of $\mathrm{dt}=1\times {10}^{-8}\mathrm{s}$).

A great amount of work has been performed by simulating plasma evolution and hydrodynamics of liner implosion independently. Plasma evolution during formation and its stability during compression has been studied by dedicated plasma codes such as Versatile Advection Code (VAC) [2], NIMROD [3] and CORSICA [4] and DCON [5], whereas the shape and trajectory of the imploding liner have been provided as an input. Dynamics and stability of the rotating imploding liquid metal liner has been extensively simulated using the hydrodynamic “compressibleInterFoam” solver, which is part of the open source C++ libraries of OpenFOAM [6]. Simulating plasma evolution and liner hydrodynamics separately has its limitations, as some degree of coupling between magnetized plasma and liquid metal liner is expected during compression. This coupling is case dependent and it is difficult to predict in prior how significant it is in order to reliably simulate compression process. Current research focuses on investigating how significant the effect of magnetic field is on the liner dynamics during implosion. As the vacuum magnetic field (or magnetic target) gets compressed, the shape and trajectory of the liner may be affected by: (i) diffusion of magnetic field into the liner and the surrounding electrically conductive structure and (ii) the magnetic pressure pushing back from the compressed magnetic field. The effect of the magnetic field on the evolution of surface perturbations at the inner liner surface during implosion can be also studied by the including MHD effects into the simulations. More complex issues of how, for example, the stability of the plasma is affected by interaction of the magnetic field with the liquid metal liner is out of scope for this study.

Several concurrent approaches have been persuaded in General Fusion to achieve coupling between the magnetic field/plasma evolution and the liquid metal liner implosion dynamics. One of them is to simulate the plasma and liquid liner by separately dedicated codes and couple at the interface, i.e., VAC-OF coupling [7]. Another one, is to simulate the liner dynamics with a hydro code and then add part of the liner (along with velocity field) as an additional layer surrounding the plasma domain in the VAC plasma code [2]. With this approach, it is assumed that the liner trajectory is not affected by the plasma evolution, but the diffusion of the magnetic field into the liner is captured during compression.

The current approach is to extend the multiphase hydrodynamic solver, which we extensively use to simulate imploding liquid liner dynamics (“compressibleInterFoam” within OpenFOAM) to include the evolution of the magnetic fields, i.e., to obtain a two-phase MHD solver. “compressibleInterFoam” is a multiphase solver which uses VOF (Volume of Fluid) phase-fraction-based interface-capturing approach and is suitable for modeling two compressible immiscible fluids [8]. It has been successfully used in General Fusion to simulate imploding liners dynamics and study interface instabilities and means of their suppression [9,10,11]. Number of recent studies attempted to extend different hydrodynamic solvers within OpenFOAM to MHD solvers [12,13,14,15,16,17,18,19,20] or coupling OpenFOAM with another software such as Elmer [21,22] to tackle various problems requiring to take into account MHD effects. One of the main messages coming out of those works is that the approach of extending existing solvers within OpenFOAM toward MHD solvers is promising. We also hope that the new solver has a potential in the future to be applied and complement other theoretical and reduced-order studies which include magentic fields, such as in recent works [23,24], for example. The rest of the paper is organized as follows. The model equations for the newly developed solver (“mhdCompressibleInterFoam”) are summarized in Section 2. Results are presented and discussed in Section 3, which is subdivided into Section 3.1, Section 3.2 and Section 3.3 describing verification and validation cases, current capabilities of the code, and application to problems of interest, respectively. The verification and validation tests section is further subdivided into Section 3.1.1, Section 3.1.2 and Section 3.1.3 to discuss cylindrical compression of magnetic fields in vacuum, Hartmann annular flow, and compression of magnetic target in a cylindrical geometry, respectively. A brief summary and conclusions are provided in Section 4.

2. Governing Equations

The evolution of a magnetic field $\mathbf{B}$ of an electrically conductive fluid is described by an induction equation, Equation (1). The Induction Equation (1) is derived from Maxwell’s equations (Faraday’s and Ampere’s laws) combined with with Ohm’s law and eliminating electric field and the electric current from the equations. For details, see [25]

$$\frac{\partial \mathbf{B}}{\partial t}=\nabla \times \left(\mathbf{v}\times \mathbf{B}\right)-\nabla \times \left(\eta \nabla \times \mathbf{B}\right),\nabla ·\mathbf{B}=0,$$

(1)

where $\eta [{\mathrm{m}}^2/\mathrm{s}]$ is the magnetic diffusivity ($\eta ={\eta }^{\prime }/{\mu }_o$; ${\eta }^{\prime }-$ resistivity and ${\mu }_o-$ vacuum permeability). At present, the extension of the hydrosolver is limited to the axisymmetric problems. In this case, magnetic field can be presented in a form such that divergence-free condition is satisfied by construction. This approach is called poloidal–toroidal decomposition, and the details can be found in [26,27]. Solenoidal magnetic field $\mathbf{B}$ can be then written as Equation (2) in terms of toroidal function $F(r,z)=r{B}_{\varphi }$ ($B_{\varphi }=B_{\mathrm{tor}}$) and poloidal flux (per radian) function $\psi (r,z)$ (e.g., [28]):

$$\mathbf{B}=\underset{B_{pol}}{\underbrace{\nabla \psi \times \nabla \varphi }}+\underset{B_{tor}}{\underbrace{F\nabla \varphi }},$$

(2)

where $\varphi$ is a polar angle. Components of the magnetic field for axisymmetric problems are shown in Figure 2.

Induction Equation (1) is reduced to two scalar equations (Equations (3) and (4)) describing the evolutions of $F(r,z)$ and $\psi (r,z)$ (see [28] for details).

$$\frac{1}{r^2}\frac{\partial F}{\partial t}+\left[\nabla ·\left(\frac{\mathbf{v}}{r^2}F\right)-\nabla ·\left(\frac{\eta }{r^2}\nabla F\right)\right]-\nabla ·\left(\omega B_{pol}\right)=0,$$

(3)

where $\omega =V_{\varphi }/r$ is the fluid angular velocity.

$$\frac{\partial \psi }{\partial t}+\mathbf{v}·\nabla \psi -\eta \underset{{\mathsf{\Delta }}^{*}\psi }{\underbrace{\left[r^2\nabla ·\left(\frac{1}{r^2}\nabla \psi \right)\right]}}=0.$$

(4)

The current density $\mathbf{J}$ (Ampere’s law), Lorentz force, and Ohmic heating (Joule heating) $P_{ohm}$ can be then expressed in terms of F and $\psi$ as in Equations (5)–(7), respectively (see [28] for details).

$$\mathbf{J}=\frac{1}{{\mu }_o}\nabla \times \mathbf{B}=\frac{1}{{\mu }_o}\left(-{\mathsf{\Delta }}^{*}\psi \nabla \varphi +\nabla F\times \nabla \varphi \right),$$

(5)

$$\begin{array}{cc}\hfill \mathrm{Lorentz}\mathrm{force}=\mathbf{J}\times \mathbf{B}=& -\frac{\nabla \left(F^2\right)}{2{\mu }_o{r}^2}-\frac{1}{{\mu }_o{r}^2}{\mathsf{\Delta }}^{*}\psi \nabla \psi +\frac{1}{{\mu }_o}\left(B_{pol}·\nabla F\right)\nabla \varphi ,\hfill \\ \hfill & {\mathsf{\Delta }}^{*}\psi =r^2\nabla ·\left(\frac{1}{r^2}\nabla \psi \right).\hfill \end{array}$$

(6)

$$P_{ohm}=\frac{\eta }{{\mu }_o}\left(\frac{{|\nabla F|}^2}{r^2}+\frac{{\left({\mathsf{\Delta }}^{*}\psi \right)}^2}{r^2}\right).$$

(7)

Equations (3) and (4) were added to the set of the governing equations of the “compressibleInterFoam” solver, which itself is an extension (to account for compressibility) of the extensively validated solver ’interFoam’ [6]. Governing equations (conservation laws) implemented in the “compressibleInterFoam” solver are listed below. Momentum and energy equations have been modified to include Lorentz force (Equation (6)) and Ohmic heating (Equation (7)), respectively. Those terms can be activated separately in gas and liquid regions. For description of conservation laws (Equations (8)–(10)) and details of their implementation in OpenFOAM, see the recent book by Greenshields and Weller [29].

The mass continuity equation is:

$$\frac{\partial \rho }{\partial t}+\nabla ·\left(\rho \mathbf{v}\right)=0.$$

(8)

The momentum equation is:

$$\begin{array}{cc}\hfill \frac{\partial \left(\rho \mathbf{v}\right)}{\partial t}+\nabla ·\left(\rho \mathbf{v}\mathbf{v}\right)=& -\nabla \left(p+\frac{2}{3}\mu \nabla ·\mathbf{v}\right)+\left[\nabla ·\left(\mu \nabla \mathbf{v}\right)+\nabla \mathbf{v}·\nabla \mu \right]+\rho \mathbf{g}\hfill \\ \hfill & +\mathrm{surface}\mathrm{tension}+\underset{\mathrm{Lorentz}\mathrm{force}(\mathrm{Equation}\left(6\right))}{\underbrace{\mathbf{J}\times \mathbf{B}}}\hfill \end{array}$$

(9)

The energy equation is:

$$\frac{\partial \left(\rho C_{p}T\right)}{\partial t}+\nabla ·\left(\rho \mathbf{v}{C}_{p}T\right)=\nabla ·\left(k\nabla T\right)+\underset{\mathrm{source}\mathrm{terms}}{\underbrace{S_T}}+\underset{\mathrm{Ohmic}\mathrm{heating}(\mathrm{Equation} \left(7\right))}{\underbrace{P_{ohm}}}$$

(10)

In Equations (8)–(10), $\rho$, $\mu$, and $\mathbf{g}$ are density, viscosity, and gravity acceleration, respectively; $C_p$ and k are the heat capacity and thermal conductivity.

In VOF-like methods, in addition to the conservation equations for mass, momentum, and energy (Equarions (8)–(10)), an equation for the phase fraction (Equation (11)) is also to be solved [30]. Additional details of implementation of the volume-of-fluid (VOF) method in OpenFOAM can be found in the Ph.D. Thesis of Rushe [31]. A detailed description of the interFoam solver can be found in Deshande et al. [32].

The phase fraction continuity equation is:

$$\frac{\left(\rho {\alpha }_i\right)}{\partial t}+\nabla ·\left(\rho {\alpha }_i\mathbf{v}\right)=0,$$

(11)

where subscript i denotes the phase and $\alpha$ is the phase fraction.

There are several variants of the VOF approach; see, for example, Ferziger and Pericć [30]. In OpenFOAM, both fluids are treated as a single fluid “mixture” whose properties vary in space according to the volume fraction of each phase (Equation (12)). With this approach, the interface is not treated as a boundary and no boundary conditions are prescribed to it; instead, it is treated as a discontinuity in fluid properties [30].

$$\rho ={\alpha }_1{\rho }_1+{\alpha }_2{\rho }_2,\mu ={\alpha }_1{\mu }_1+{\alpha }_2{\mu }_2,k={\alpha }_1{k}_1+{\alpha }_2{k}_2,\eta ={\alpha }_1{\eta }_1+{\alpha }_2{\eta }_2,$$

(12)

where subscripts 1 and 2 denote the two fluids (e.g., liquid and gas) and ${\alpha }_1+{\alpha }_2=1$.

It is worth noting, that “mixture” is an important methodological approach that is particularly useful for applications. One of the areas where the “mixture” approach, assisted by extended irreversible thermodynamic frameworks, leads to the promising results with significantly lower computational loads is that of magnetorheological fluids [33,34].

Finally, the equation of state (Equation (13)) for each phase has to be provided to close the system:

$${\rho }_i={\rho }_i(p,T),$$

(13)

where several options are available for each phase depending on the application of interest.

Equations (3)–(13) are the governing equations for the new “mhdCompressibleInterFoam” solver. At present, constant fluid properties (viscosity, thermal conductivity, magnetic diffusivity) are used in simulations.

3. Results

3.1. Verification Tests

3.1.1. Compression of Toroidal and Uniform Vertical Poloidal Vacuum Magnetic Fields

In this section, two verification tests for cylindrical compression of vacuum magnetic fields by imploding liquid liner are considered. In the first test, a compression of a toroidal magnetic field introduced into a vacuum (gas) region at the start of simulation is discussed, while in the second test, the uniform vertical (poloidal) magnetic field is added in addition to the toroidal one. One should note that with the way VOF approach is implemented, it is not possible to simulate real vacuum (or low density plasma), as density ratio ${\rho }_{\mathrm{liquid}}/{\rho }_{\mathrm{gas}}$ should not be too extreme for VOF to work properly. Using low density values in the gas also leads to a decrease in the time step, which, in turn, increases computational cost. Therefore, to mimic the vacuum from electromagnetic point of view, the following is done: (i) the initial hydrodynamic pressure in the gas region is set to a sufficiently low value such that its effect on liner trajectory is negligible, and (ii) the magnetic field in the gas region is prescribed by an external procedure or by using a high value of magnetic diffusivity (the magnetic diffusivity of vacuum is infinite).

In the first verification test, a toroidal field in the vacuum (gas) region is instantaneously produced by the electrical current $I_o$ (driven by the external power supply) that runs through the conducting shaft placed in the center of the cylinder. The external power source is then disconnected. With such a setup, a total toroidal flux in the system (vacuum and liquid liner) remains constant. At the start of the simulation, a toroidal magnetic field is present only in the vacuum (gas) region, starting then to diffuse into the liquid liner. The exact amount of the initial toroidal flux that diffuses into the liner during implosion depends on the compression time and trajectory. During the implosion, the toroidal field also gets compressed, leading to the increase in the magnetic pressure acting on the liner, which in turn alters the compression trajectory.

A cross-section of the computational setup used for the cylindrical compression of vacuum magnetic fields test cases is shown in Figure 3. To simulate 2-dimensional axisymmetric cases in OpenFOAM, the geometry is specified as a wedge of small angle ($5^{\circ }$ in this work) and 1 cell thick running along the plane of symmetry. Front and back wedge planes must be specified as separate patches of wedge type (see OpenFOAM documentation [6] and reference [29] for more details). Position of the central shaft and direction of the electrical current $I_{shaft}=I_o\overrightarrow{z}$ are shown in the figure (Figure 3a) along with a direction of generated by the shaft current toroidal magnetic field $B_{tor}$. Direction of the uniform vertical magnetic field $B_{pol}=B_o\overrightarrow{z}$ (for the poloidal field compression test case discussed later) is also shown in Figure 3a). The initial position and geometry of the liquid liner and vacuum (gas) cavity are shown by red and blue colors, respectively, in Figure 3a. Inner boundary of the computational domain marked by number “3”, corresponds to the central shaft with radius $R_{shaft}$. No-slip wall boundary condition is imposed at this boundary. At top and bottom boundaries, marked by “1” in Figure 3a, a slip wall boundary condition is used. At the outer boundary (marked by number “2”), a pressurized pushing gas causing liner implosion is prescribed as pressure profile $P_{out}\left(t\right)$, that is constant for this test case. Zero-gradient boundary condition is specified for toroidal function F at all boundaries. Physically, this corresponds to boundaries with perfect electric conductivity and results in conservation of total toroidal magnetic flux inside computational domain.The instant when the toroidal magnetic field is highly compressed and the liquid liner is near its turnaround point is shown in Figure 3b. As mentioned earlier, a low pressure gas is also initially present in the cavity (in addition to the toroidal magnetic field), but the parameters are chosen in such a way that the effect of the magnetic field on the liner trajectory is a dominant one.

For this first verification case, the toroidal magnetic field (represented by the toroidal function $F(r,z)$) is calculated by an external procedure using conservation of the total toroidal flux. Hence, the toroidal function F in the gas region is overwritten at each time step, while Equation (3) is solved in the liquid. The test case parameters are also chosen to ensure that the magnetic field is not diffused through the entire thickness of the liquid liner during the implosion such that no magnetic flux exiting through the outer liner surface.

The procedure for calculating the toroidal function $F_{gas}\left(t\right)$ in the gas cavity is detailed below. The initial toroidal flux ${\mathsf{\Phi }}_0$ is given by:

$${\mathsf{\Phi }}_o=I_o·L_o;L_o=\frac{{\mu }_o}{{\left(2\pi \right)}^2}{\int }_{V_{\mathrm{gas}}}\frac{1}{r^2}dV,$$

(14)

where $I_o$ is the initial shaft current, $L_o$ is the initial inductance of the gas cavity, and ${\mu }_o=4\pi \times {10}^{-7}\mathrm{H}/\mathrm{m}$ is the vacuum permeability. The relation between the initial shaft current and the initial value of toroidal function in the gas cavity is given by (according to the Biot–Savart law):

$$F_{\mathrm{gas}}(t=0)={\mu }_o/\left(2\pi \right)·I_o.$$

(15)

Using Equations (14) and (15), the initial toroidal flux ${\mathsf{\Phi }}_o$ can be calculated as:

$${\mathsf{\Phi }}_o=\frac{1}{2\pi }·F_{\mathrm{gas}}(t=0)·{\int }_{V_{\mathrm{gas}}}\frac{1}{r^2}dV.$$

(16)

Once the magnetic field starts diffusing into a liquid liner, the total conserved toroidal flux consists of that in the gas cavity and that in the liquid liner:

$${\mathsf{\Phi }}_o={\mathsf{\Phi }}_{\mathrm{total}}\left(t\right)=\underset{\mathrm{flux}\mathrm{in}\mathrm{gas}}{\underbrace{\frac{1}{2\pi }{F}_{\mathrm{gas}}\left(t\right)·{\int }_{V_{\mathrm{gas}}}\frac{1}{r^2}dV}}+\underset{\mathrm{flux}\mathrm{in}\mathrm{liquid}}{\underbrace{\frac{1}{2\pi }{\int }_{V_{\mathrm{liquid}}}\frac{F\left(t\right)}{r^2}dV}},$$

(17)

where the first integral is calculated over the volume occupied by the gas cavity and the second one is over the volume of the liquid liner. Thus, the value of $F_{\mathrm{gas}}$ at any given time t can be calculated from Equation (17) as:

$$F_{\mathrm{gas}}\left(t\right)=\left({\mathsf{\Phi }}_o-\frac{1}{2\pi }{\int }_{V_{\mathrm{liquid}}}\frac{F\left(t\right)}{r^2}dV\right)/\frac{1}{2\pi }{\int }_{V_{\mathrm{gas}}}\frac{1}{r^2}dV.$$

(18)

The simulation procedure can be summarized as follows:

• Step 1: For a given initial shaft current $I_o$ and size of the gas cavity, calculate the initial value of the toroidal function in the gas $F_{\mathrm{gas}}(t=0)$ by Equation (15) and the initial toroidal flux ${\mathsf{\Phi }}_o$ by Equation (16).
• Step 2: Perform one iteration with “mhdCompressibleInterFoam”.
• Step 3: Calculate an updated value of the toroidal function in the gas $F_{\mathrm{gas}}\left(t\right)$ by Equation (18).
• Step 4: Overwrite the toroidal function F for all the grid cells inside the gas cavity with the value calculated in Step 3.
• Step 5: Repeat steps 2 to 4 for the duration of the simulation.

Examples of inner liner surface trajectories obtained with the above prescribed procedure are given in Figure 4. The following set of parameters has been used for those simulations: initial radius of the inner liner surface $R_{\mathrm{in}}=0.9$m, initial radius of the outer liner surface $R_{\mathrm{out}}=1.3$m, cylinder height $H=2\mathrm{m}$, central shaft radius $R_{\mathrm{shaft}}=0.2$m, constant pushing pressure $P_{\mathrm{out}}=1\times {10}^6$Pa, initial gas cavity pressure $P_{\mathrm{in}}=1\times {10}^4$Pa, and constant liner density of ${\rho }_{\mathrm{liquid}}=1000\mathrm{kg}/{\mathrm{m}}^3$. Structured orthogonal mesh has been used with a uniform $1$mm grid resolution in the radial direction. In the azimuthal direction, $1$mm grid resolution has been imposed at the radius of the initial inner liner surface. Isothermal calculation has been run, i.e., all liquid phase properties (including magnetic diffusivity) have been held constant throughout the simulation. For the gas phase, this implies isothermal compression. Here, it is worth reiterating that for this test case, the dynamics of the gas is not of an interest, provided that gas has a negligible effect on the liner trajectory. What is of interest is the effect of the toroidal magnetic field on the liner trajectory and the amount of the magnetic flux diffused into the liner during implosion.

In Figure 4a,b, results are shown for the initial shaft currents of $I_o=2$MA and $3$MA, respectively. Two sets of curves (dashed and dash dot) corresponding to the two values of magnetic diffusivity (${\eta }_{\mathrm{liquid}}=0.4{\mathrm{m}}^2/\mathrm{s}$ and ${\eta }_{\mathrm{liquid}}=0$) are plotted for each value of the initial shaft current. Also plotted by the solid line is a trajectory obtained in a pure hydrodynamic simulation, i.e., for a zero value of the shaft current ($I_o=0$). Results obtained with the newly developed solver “mhdCompressibleInterFoam” are shown by the red lines and the ones obtained with in-house 1D code are shown by the black lines. (Description of the in-house 1D code is provided in Appendix A. Equations describe only dynamics of the liner with pushing pressure and cavity pressure are used as boundary conditions. Gas cavity pressure is recalculated based on the inner liner surface position assuming adiabatic (or isothermal) compression. Toroidal function in the gas $F_{\mathrm{gas}}$ is calculated from total toroidal flux conservation the same way is in “mhdCompressibleInterFoam”).

One can see that with no magnetic pressure acting on the liner (solid lines in Figure 4), it implodes faster, goes deeper into compression, and is characterized by a sharp rebound occurring close to the central shaft. Results obtained with “mhdCompressibleInterFoam” (that reduces to the original “compressibleInterFoam” solver for the case of $I_o=0$) are in good agreement with those produced with our in-house 1D code during the entire implosion phase. As in this case, the liner comes close to the shaft, capturing rebound using OpenFOAM is challenging, as it requires fine resolution and small times steps. Hence, the OpenFOAM simulation is limited to the imploding stage only.

By examining the trajectories obtained with toroidal magnetic field in the gas cavity, one can observe the following: (i) with an increase in the magnetic field strength, the liner turn around radius increases for the same value of magnetic diffusivity, (ii) for a higher value of magnetic diffusivity, the liner goes deeper into a compression as some portion of the magnetic flux diffuses into the liner reducing overall magnetic push-back, (iii) there is a good agreement between “mhdCompressibleInterFoam” results and the in-house 1D code until the turnaround point for all cases. It can also be noticed that for the case of zero magnetic diffusivity (dashed lines), there is a slight deviation between “mhdCompressibleInterFoam” and 1D code. This happens because even when the magnetic diffusivity is set to zero, some magnetic flux getting diffused into the liner due to the numerical diffusion. Finally, we should emphasize, that with the current numerical setup, it is not always possible to run the simulation much longer beyond the turn around point. This happens because during the rebound stage, the magnetic field in the liner close to the inner surface might become higher than in the gas cavity. This might put liquid under tension and cause delamination of the inner liner surface [35]. Mechanism of the liner surface delamination is not captured correctly at present, as it requires a cavitation model and addition of the third phase (cavitated liquid) into simulation. This is outside the scope of the current work, whose primary goal is to simulate liner implosion up-to the turn around point.

The second verification case is an extension of the first one, where a uniform vertical magnetic field is also present in the gas cavity at the start of the simulation. The computational setup is the same as in the previous case and is shown in Figure 3. A vertical uniform poloidal field is initially present only in the gas cavity, it then diffuses into a liquid liner and also gets compressed during the liner implosion similar to the case with a toroidal magnetic field. For a uniform vertical poloidal magnetic field $B_{\mathrm{pol}}=B_o\overrightarrow{z}$, a poloidal flux function $\psi$ is of a form $\psi =c_1{r}^2+c_2$, where $c_1$ and $c_2$ are constants. The initial distribution of the poloidal flux function in the gas cavity ${\psi }_{\mathrm{gas}}(t=0)$ is then prescribed as:

$${\psi }_{\mathrm{gas}}(t=0)={\psi }_{\mathrm{shaft}}·\left(R_{\mathrm{in}}^2-r^2\right)/\left(R_{\mathrm{in}}^2-R_{\mathrm{shaft}}^2\right);{\psi }_{\mathrm{shaft}}=0.5·B_o·\left(R_{\mathrm{in}}^2-R_{\mathrm{shaft}}^2\right),$$

(19)

where ${\psi }_{\mathrm{shaft}}$ is the value of the poloidal flux function at the shaft, $B_o$ is the strength of the uniform vertical poloidal field, $R_{\mathrm{in}}$ and $R_{\mathrm{shaft}}$ are the initial radii of the inner liner surface and the central shaft, respectively. In this verification test, the equations for the evolution of magnetic field (Equations (3) and (4)) have been solved in both gas and liquid regions. Here a high value of magnetic diffusivity inside the gas is used to mimic vacuum. This is instead of using an external procedure to calculate magnetic field in the gas region, as it has been done in the previous verification test. At the inner boundary $r=R_{\mathrm{shaft}}$ (marked by “3” at Figure 3a), a zero gradient boundary condition is imposed for the toroidal function F (superconductor shaft) and a fixed value of $\psi ={\psi }_{\mathrm{shaft}}$ for the poloidal flux function. As in the previous test case, it is assumed that the magnetic field does not reach the back side of the liner, hence $\psi =0$ and a zero gradient is prescribed for F at the outer $r=R_{\mathrm{out}}$ boundary. Zero-gradient boundary condition is imposed at top and bottom boundaries (marked by “1” at Figure 3a) for both toroidal function F and poloidal flux function $\psi$. Several tests with different values of magnetic diffusivity in the gas have been carried out and the results have converged to those obtained with using an external procedure to calculate the magnetic field in the cavity for the values of magnetic diffusivity ${\eta }_{\mathrm{gas}}>20$.

Trajectories of the inner liner surface during the compression of poloidal, toroidal and combination of both are shown in Figure 5. Results obtained with “mhdCompressibleInterFoam” are shown by the symbols, whilst those obtained with the in-house 1D code (see Appendix A for the details) are plotted by the solid lines. Case parameters are similar to the first verification test; $R_{\mathrm{in}}=0.9$m, $R_{\mathrm{out}}=1.3$m, $R_{\mathrm{shaft}}=0.2$m, $H=2\mathrm{m}$, $P_{\mathrm{out}}=1\times {10}^6$Pa, $P_{\mathrm{in}}=1\times {10}^4$Pa, ${\rho }_{\mathrm{liquid}}=1000\mathrm{kg}/{\mathrm{m}}^3$, ${\eta }_{\mathrm{gas}}=100{\mathrm{m}}^2/\mathrm{s}$, ${\eta }_{\mathrm{liquid}}=0.4{\mathrm{m}}^2/\mathrm{s}$. Black, red and blue curves correspond respectively to (i) $B_o=0.4$T, $I_o=0$, (ii) $B_o=0$, $I_o=2$MA and (ii) $B_o=0.4$T, $I_o=2$MA cases. First of all, one can see that for all three cases, there is good agreement between the results obtained with “mhdCompressibleInterFoam” and those of 1D code during the implosion phase. For a chosen set of parameters, effects of poloidal only (black) or toroidal only (red) vacuum fields on the liner trajectory are similar. When combined (blue), the liner’s turnaround radius increases, the time it takes to reach a turnaround point also increases and the turnaround is more gradual.

Radial profiles of the poloidal flux function $\psi$ and toroidal function F at several times during the compression are shown in Figure 6 and Figure 7 for the case where the poloidal and toroidal vacuum magnetic fields are initially present in the cavity (case is shown by the blue line in Figure 5; $B_o=0.4$T, $I_o=2$MA). Results obtained with “mhdCompressibleInterFoam” solver are shown by the red lines and those obtained with 1D code by the blue ones. MHD equations are solved in both liquid liner and gas cavity region with “mhdCompressibleInterFoam” solver, whereas the 1D code solved only the dynamics of the liquid liner. “Ends” of the blue lines correspond to the radial position of the inner surface of the liner at different times; it can be seen that they correlate with the points on the blue curve in Figure 5 at those times. The latest time for which profiles are shown, is about the liner turnaround point. Diffusion of both toroidal and poloidal magnetic fields into the liquid liner is clearly seen, and there is good agreement between “mhdCompressibleInterFoam” and in-house 1D code for the magnetic field profiles inside the liner. In the gas cavity region, a parabolic profile of the poloidal flux function $\psi$ corresponding to the uniform vertical field, can be also seen at all times presented. A constant value of the toroidal function F in the gas cavity, which increases as the toroidal field gets compressed, is clearly seen in Figure 7. It is also worth reiterating, that for a given set of parameters, the magnetic field has not diffused through the entire thickness of the liquid liner, such that no magnetic flux leaves through the outer liner surface during the implosion.

Finally, we would like to mention, that a number of grid convergence tests have been run for the presented test cases. Grid cell size was varied between $0.5\mathrm{mm}$ and $5\mathrm{mm}$. The effect of the time step has been also studied. Main outcome of these studies is that pretty high resolution grid is required to obtain a converged solution for the diffusion of magnetic field into the liner, especially for the small values of magnetic diffusivity. Even finer resolution is required to accurately resolve thermal layer (due to Ohmic heating) in the liquid liner if it is of interest. To achieve grid independent results, a higher grid resolution is required compared to the pure hydrodynamic simulations of imploding liner. A grid with $1\mathrm{mm}$ square cells has been found appropriate for our problems of interest. Running simulations on structured meshes with square cells (as much as possible) and keeping a time step small and nearly constant during the simulation was unsurprisingly found to lead the best results.

3.1.2. Hartmann Annular Flow

In this section the “mhdCompressibleInterFoam” solver is validated against an analytical solution for a Hartmann Annular Flow. Hartmann annular flow is an extension of a classical MHD Hartmann flow benchmark problem [36] to the case of an annular channel with a rectangular cross-section [37], see a schematic of a computational setup in Figure 8. Geometry of the test case is defined by the inner and outer radii of the annular channel $R_{\mathrm{in}}$ and $R_{\mathrm{out}}$, respectively, and by the height of the channel H. A channel is filled with an electrically conducting fluid which is initially at rest. A non-slip boundary condition is applied at all walls of the channel. A vertical uniform magnetic field $B_o$ acting on a fluid is imposed through the boundary condition on the poloidal flux function (per radian) $\psi$; $\psi (R_{\mathrm{in}},z)=0.5B_o{R}_{\mathrm{in}}^2$—at the inner wall, $\psi (R_{\mathrm{out}},z)=0.5B_o{R}_{\mathrm{out}}^2$—at the outer wall, $\psi (r,0)=\psi (r,H)=0.5B_o{r}^2$ — at the bottom and top walls. The initial poloidal field can be either imposed at the start of the simulation (to speed up reaching the steady-state) or let it to diffuse from the boundaries as simulation progresses. To induce a poloidal current, a fixed value boundary condition for a toroidal function F is applied at the top and bottom walls as: $F(r,0)=-F_o$ and $F(r,H)=F_o$. A zero gradient boundary condition for F is applied at the inner and outer radii walls (superconductor).

With such a setup, an electrically conductive fluid starts to swirl inside the channel due to the Lorentz force (in particular due to the third term in Equation (6)) and eventually reaches a steady state when balanced by the viscous forces. An example of swirling velocity contours at the channel cross-section is shown in Figure 8, where magnitude of the swirling velocity varies from the maximum (dark blue) inside the channel to zero (red) at the wall. In the middle section of the annular channel, away from the inner and outer radii walls (schematically indicated by the green rectangular in Figure 8), this flow problem can be analytically solved for a steady state. In this central region the laminar steady state flow is described by the following equations (see [37] for details):

$$\left\{\begin{array}{c}\frac{B_o}{{\mu }_o}\frac{\partial F}{\partial z}+\rho \nu \frac{{\partial }^{2}L}{\partial z^2}=0\hfill \\ B_o\frac{\partial L}{\partial z}+\eta \frac{{\partial }^{2}F}{\partial z^2}=0,\hfill \end{array}\right.$$

(20)

where $L\equiv r{V}_{\varphi }(r,z)$ is an angular momentum, F(z) is the toroidal function, $B_o$ is the strength of the uniform vertical magnetic field, $\eta$ is magnetic diffusivity, ${\mu }_o$ is vacuum permeability, $\rho$ and $\nu$ are density and kinematic viscosity of the fluid. For the case of an annular channel entirely filled with liquid, the above system is solved with the following set of boundary conditions:

$$\begin{array}{c}F(z=0)=-F_o\hfill \\ F(z=H)=F_o\hfill \\ L(z=0)=0\hfill \\ L(z=H)=0.\hfill \end{array}$$

(21)

The shapes of the angular momentum L and toroidal function F profiles depend on the Hartmann number that is defined as:

$$Ha=B_{o}H/\sqrt{{\mu }_0\rho \eta \nu }.$$

(22)

The Hartmann number represents the ratio of electromagnetic and viscous forces, hence, depending on the value of $Ha$ number, a flow regime is dominated by the viscous forces (low $Ha$) or by electromagnetic forces (high $Ha$).

Validation tests have been carried out for the following set of parameters: $R_{\mathrm{in}}=0.1$m, $R_{\mathrm{out}}=0.15$m, $H=0.01$m, $\rho =1000\mathrm{kg}/{\mathrm{m}}^3$, $\nu =1{\mathrm{m}}^2/\mathrm{s}$, $\eta =1{\mathrm{m}}^2/\mathrm{s}$, $F_o=0.01\mathrm{T}·\mathrm{m}$ and three different strengths of the vertical uniform magnetic field $B_o=1\mathrm{T}$, $40\mathrm{T}$ and $80\mathrm{T}$. For those parameters, the corresponding test cases Hartmann numbers are $Ha=0.28$, 11 and 22. Steady state laminar flow simulations have been run on a computational grid with square cells at different resolutions. Grid convergence has been confirmed and the results obtained for the grid with $1\mathrm{mm}$ cell size are presented below. Vertical angular momentum and toroidal function profiles at several radial positions (away from the inner and outer walls) along with analytical solution obtained by solving Equation (20) with boundary conditions given by Equation (21) are plotted in Figure 9 for different Hartmann numbers. Analytical solutions have been obtained for the same set of parameters used in numerical simulations, i.e., units of measurement for individual physical quantities are identical to those used in simulations.

First of all, one can see, that there is an excellent agreement between the results obtained with the “mhdCompressibleInterFoam” solver and the analytical solution for all the test cases. A slight deviation from the analytical solution is observed for the angular momentum profile plotted at the radial position of $r=0.14\mathrm{m}$ at Hartmann number $Ha=0.28$ (part a). This is likely because viscous forces dominate this regime, and the influence of the outer radius wall extends to that radial position. As it can be expected, at a low Hartmann number, the velocity exhibits a parabolic shape profile and the poloidal electrical current (proportional to $\nabla F$) is uniformly distributed through out the channel. With the increase in Hartmann number, development of the Hartmann layers near the top and bottom walls is clearly observed and the velocity profile changes from the parabolic to the flat top shape. Hartmann layers become thinner with a further increase in Hartmann number such that electrical currents are concentrated inside those thin layers near the top and bottom walls, whereas the central part of the channel is electrical current free (area of constant toroidal function F).

3.1.3. Compression of Magnetic Target in a Cylindrical Geometry

In this section a more complex verification test for the compression of a simplified magnetic target is discussed. Results of an axisymmetric compression simulation of a simplified magnetic target by an imploding liquid metal liner driven by a pressurized gas in a cylindrical geometry are presented and compared to those obtained with VAC code. A schematic of the initial simulation setup is shown in Figure 10a, while a highly compressed target is shown in Figure 10b. Parameters for this test case have been chosen within the range relevant to General Fusion’s FDP. Geometry of the test case is defined by the radii of inner $R_{\mathrm{in}}$ and outer $R_{\mathrm{out}}$ liner surfaces, radius of the central shaft $R_{\mathrm{shaft}}$ and height of the cylinder H.

Magnetic target to be compressed is specified as an initial condition in terms of the poloidal flux function per radian $\psi (r,z)$ and toroidal function $F(r,z)$ and is calculated as a solution of Grad-Shafranov (GS) equilibrium equation [39,40] given by Equation (23). For a rectangular cross-section and pressure and toroidal function profiles given by Equation (24), Grad-Shafranov equation can be solved analytically (as detailed below).

$$\frac{{\partial }^2\psi }{\partial r^2}-\frac{1}{r}\frac{\partial \psi }{\partial r}+\frac{{\partial }^2\psi }{\partial z^2}+\frac{1}{2}\frac{d{F}^2}{d\psi }+{\mu }_o{r}^2\frac{dp}{d\psi }=0,$$

(23)

where ${\mu }_o$ is the magnetic permeability, $p\left(\psi \right)$ is the pressure and $F\left(\psi \right)$ is a toroidal function. The nature of the equilibrium is mainly determined by the choices of the two functions $F\left(\psi \right)$ and $p\left(\psi \right)$ as well as the boundary conditions. Here it is assumed that:

$$F^2=F_o^2+{\lambda }^2{\psi }^2,p=p_o,$$

(24)

where $F_o,\lambda$ and $p_o$ are some constants. Here $F_o$ corresponds to the value of the initial central shaft current $I_{\mathrm{shaft}}$ and pressure is assumed to be constant and equal to $p_o$. Substituting Equation (24) into Equation (23) leads to a specific form of Grad-Shafranov equation:

$$\frac{{\partial }^2\psi }{\partial r^2}-\frac{1}{r}\frac{\partial \psi }{\partial r}+\frac{{\partial }^2\psi }{\partial z^2}+{\lambda }^2\psi =0,$$

(25)

with boundary condition of $\psi =0$ at all boundaries of the domain. A rectangular cross-section initially occupied by the magnetic target, is constrained by the central shaft radius $r=R_{\mathrm{shaft}}$, initial radius of the inner liner surface $r=R_{\mathrm{in}}$, bottom wall $z=0$ and top wall $z=H$ (see Figure 10a). General solution to Equation (25) is:

$$\psi =\left((ar{J}_1\left(kr\right)+br{Y}_1\left(kr\right)\right)sin\left(\frac{\pi z}{H}\right),$$

(26)

where $J_1$ and $Y_1$ are Bessel functions. Radial wave-number k and amplitude constants a and b are determined from boundary and initial conditions. For given radii $R_{\mathrm{shaft}}=r_1$ and $R_{\mathrm{in}}=r_2$, wave number k is determined from boundary condition ($\psi =0$ at $r=r_1$ and $r=r_2$) using:

$$J_1\left(k{r}_1\right)Y_1\left(k{r}_2\right)-J_1\left(k{r}_2\right)Y_1\left(k{r}_1\right)=0.$$

(27)

Then full solution is:

$$\begin{array}{cc}\hfill & \psi =c\left(r{J}_1\left(k{r}_1\right)Y_1\left(kr\right)-r{Y}_1\left(k{r}_1\right)J_1\left(kr\right)\right)sin\left(\frac{\pi z}{H}\right),\hfill \\ & {\lambda }^2=k^2+\frac{{\pi }^2}{H^2},\hfill \\ & F=\sqrt{F_o^2+{\lambda }^2{\psi }^2}.\hfill \end{array}$$

(28)

where c in an amplitude constant determining the maximum poloidal flux.

Parameters for the verification case are given below in

Table 1. Test case parameters for comparison between “mhdCompressibleInterFoam” solver and “layered” VAC code.

	Geometry
Shaft radius	$R_{shaft}$	$0.15\mathrm{m}$
Inner liner radius	$R_{in}$	$1.5\mathrm{m}$
Outer liner radius	$R_{out}$	$2.11\mathrm{m}$
Height of the liner	H	$2\mathrm{m}$
	Liquid liner (lithium)
Density	${\rho }_{\mathrm{lithium}}$	$500{\mathrm{kg}/\mathrm{m}}^3$
Viscosity	${\mu }_{\mathrm{lithium}}$	$3.645\times {10}^{-4}\mathrm{Pa}$
Magnetic diffusivity	${\eta }_{\mathrm{lithium}}$	$0.2{\mathrm{m}}^2/\mathrm{s}$
EoS	constant density
Pushing pressure	$p_{\mathrm{out}}$	$7\times {10}^6\mathrm{Pa}$
	Magnetic target
Initial pressure (uniform)	$p_{\mathrm{target}}$	$1\times {10}^4\mathrm{Pa}$
Initial density (uniform)	${\rho }_{\mathrm{target}}$	$0.11\mathrm{kg}/{\mathrm{m}}^3$
Viscosity	${\mu }_{\mathrm{target}}$	$0.35\mathrm{Pa}\mathrm{s}$
Magnetic diffusivity (constant)	${\eta }_{\mathrm{target}}$	$1{\mathrm{m}}^2/\mathrm{s}$
Initial magnetic configuration
Shaft current	$I_{\mathrm{shaft}}$	$2\mathrm{MA}$
Maximum poloidal flux (per radian)	${\psi }_{max}$	$0.1056{\mathrm{Tm}}^2$
Parameters in Equations (27) and (28)
Wave number	k	$2.62731/\mathrm{m}$
Eigenvalue	$\lambda$	$3.0611/\mathrm{m}$
Amplitude constant	c	$0.12\mathrm{Tm}$
Shaft current constant	$F_o$	$0.4\mathrm{Tm}$
	Simulation type
Laminar, iso-thermal, high target viscosity, constant magnetic diffusivity

. Following assumptions have been made: laminar flow, incompressible liquid liner, constant magnetic diffusivity (different values for gas and liquid regions) and, for simplicity, iso-thermal. Viscosity in the gas region has been set to an artificially high value in order to damp waves bouncing inside the cavity during compression. Simulation has been run on an uniform computational grid with $1$mm cell size. Simulation has been run for the set of boundary conditions (see Figure 3 for notation): non-slip wall at the inner boundary ($r=R_{shaft}$), slip wall boundary condition at top and bottom walls ($z=0$ and $z=H$), prescribed pressure at the outer boundary ($r=R_{outer}$, fixed value of $\psi =0$ and zero-gradient for F at all boundaries).

Evolution of the magnetic field during the compression obtained with “mhdCompressibleInterFoam” solver is compared to that obtained with “layered VAC” code for the same set of parameters and physical models. “Layered VAC” code is in-house modified version of VAC code, which allows to simulate a layer of liquid metal in addition to the plasma region. Current limitation of compression simulations with VAC code is that the trajectory and shape of the liner should be known in prior (usually from hydrodynamic simulation), which then used as a moving boundary for VAC. In the “layered VAC”, the velocity field in the liquid liner layer is also provided as an input to the simulation. This allows to examine diffusion of the magnetic field into the liner during compression. With such a numerical setup we can compare evolution of the magnetic target and diffusion of the magnetic field into the liner between the two codes. What it is not possible to compare, is the effect of the magnetic field on the trajectory of the liner, as it has to be prescribed in VAC simulation.

Simulation with “layered VAC” code has been run on a structured non-uniform grid; grid resolution at the plasma-liquid metal interface is $1$mm (the same as for “mhdCompressibleInterFoam” solver) and is much coarser for the rest of the domain ($1$cm). Both codes have been run for identical set of parameters and physical models but viscosity in the gas. Inviscid fluid model is used in VAC simulation, whereas in “mhdCompressibleInterFoam”, high viscosity is set in the gas region. Initial condition is shown in Figure 11a, where magnetic target is plotted by the contours of the poloidal function $\psi$. Initial radial profiles of the poloidal function $\psi$ and toroidal function F at the equator ($z=1\mathrm{m}$) (Equation (28) with parameters in Table 1) are shown in Figure 11b.

Comparison between “mhdCompressibleInterFoam” solver and “layered VAC” for the evolution of the magnetic field during compression is presented in Figure 12. Radial profiles of toroidal function F and poloidal flux function $\psi$ at the equator at several instances during compression are given in the left and right columns of the Figure 12, respectively. First of all, one can see that there is good agreement between the results obtained with “mhdCompressibleInterFoam” solver and “layered VAC” for the evolution of the magnetic field in both gas and liquid liner regions, which is encouraging. It it also important to note that viscosity seems to have a little effect on the evolution of the magnetic field. This indicates that using a high viscosity value in “mhdCompressibleInterFoam” to damp the waves and high velocities that might occur in the vicinity of the central shaft, does not significantly influence magnetic field evolution but it does help to stabilize computation and avoid undesirable decrease in the time step.

Comparison of the toroidal function F and poloidal flux function $\psi$ at the plasma-liquid metal interface at equator ($z=1$m) during compression is plotted in Figure 13. Results are shown for two different values of volume of fraction $\alpha$, given by solid ($\alpha =0.5$) and broken ($\alpha =0.95$) red lines. Corresponding data from “layered VAC” simulation is shown by hollow circular symbols. It can be seen that there is good agreement between the results obtained with both codes. Following the toroidal function F value at the interface (Figure 13a), the next can be observed: (i) it starts at $0.4\mathrm{T}\mathrm{m}$ which corresponds to the initial shaft current of $I_{\mathrm{shaft}}=2\mathrm{MA}$; (ii) it gradually goes up for the majority of the compression and then rapidly increases towards the end. This behavior goes in line with compression trajectory that is characterized by rapid acceleration of the liner during the last millisecond of the compression. The poloidal flux function $\psi$ value at the interface (Figure 13b) can be interpreted as the amount of poloidal flux being diffused into the liner at a certain time during compression. Initially, there is no poloidal flux in the liquid liner (${\psi }_{\mathrm{interface}}(t=0)=0$). Then the poloidal flux starts to diffuse into the liner, and for the parameters and trajectory of the test case, about one third of the initial flux ends up in the liner at the time of maximum compression time (${\psi }_{\mathrm{interface}}\approx 0.035\mathrm{Wb}/\mathrm{rad}$ at $t=8.9\mathrm{ms}$ compared to ${\psi }_{\max }\approx 0.1\mathrm{Wb}/\mathrm{rad}$ of the initial target).

3.2. Incorporating Electrically Conductive “Solid” Regions into Simulation

To simulate compression of the magnetic targets in setups relevant to the experimental apparatus, it might be necessary to partly incorporate the surrounding electrically conductive structure into simulation. An example of such a computational setup, that is an extension of the one was considered in the previous Section 3.1.3 (Figure 10), is shown in Figure 14. It is worth mentioning, that in the compression scheme under design in General Fusion, magnetized plasma is injected into the vacuum cavity inside the liquid metal liner at some point during implosion. This is because a certain volumetric compression of the plasma is to be achieved in a certain time in order to be on a path to achieve fusion conditions. Hence, plasma is injected when the liner is already moving and accelerated to a certain velocity. A schematic of such a system is shown in Figure 14a. It contains the same parts as the test case in Figure 10, i.e., vacuum (gas) region surrounded by the liquid liner (area in-between the two yellow lines marked as baffles in the figure). A central shaft region is added on the left, bottom plate marked as “structure” at the bottom (separated by the yellow line from the vacuum/liner regions), top plate similar to the bottom one at the top, and finally, an opening to inject the plasma into the cavity (blue region above the upper yellow line). A simplified compression scenario can be summarized as follows: (i) evacuated cylindrical cavity inside the rotating liquid liner—as a starting point, though liner rotation is not considered at this stage, (ii) pressurized gas pushes at the outer liquid liner surface causing the liner to implode, (iii) when the inner liner surface reaches outer radius of the plasma injection system opening (Figure 14b), magnetized plasma is injected into the cavity, (iv) liner continues to implode compressing the plasma. It is important to note, that in the real system, imploding liner will need to bridge the injection opening, and some portion of the liquid metal will inevitably end up in the opening etc. Here, for simplicity, it is assumed that once magnetic target is injected, the opening is closed instantaneously, such that no liquid is allowed to enter injection area and implosion carries on as it would be in a pure cylindrical geometry. Including parts of the electrically conductive solids into simulation allows us to study diffusion of the magnetic field into the structure during compression. Electrical currents generated in the structure may lead to a strong heating of the material etc. Also, some magnetic field configuration can be present in the structure, liner and evacuated cavity prior to the implosion. Ability to simulate evolution of the magnetic field in such configurations is important for overall understanding of MHD effects during compression process. A simplified example of such a case is a fully soaked uniform vertical magnetic field, that is discussed later in this section.

For our purposes, considering rigid solid regions is sufficient, so in the regions designated as “solids”, we only solve magnetic field evolution equations (and potentially temperature), whereas phase fraction continuity and momentum equation are not solved. Solid regions are set up as a liquid phase, with a magnetic diffusivity corresponding to that of a solid material. Initial velocity field is set to zero and with momentum equation not being solved, it remains zero in those regions for the entire run. In order to specify boundary conditions at the interface between “solid” and gas/liquid regions the following procedure has been implemented: (i) Internal zero thickness baffles have been defined between solid and gas/liquid regions. Example of such baffles is shown in Figure 10 by the yellow lines. (ii) Periodic boundary conditions have been applied at those baffles for the magnetic field components, therefore, baffles are not seen by the magnetic field. For the hydrodynamic variables “wall” boundary conditions are used and, depending on the case, “slip” or “no slip” wall boundary conditions are applied. (iii) Momentum equation is not solved for the grid cells inside “solid” regions to ensure zero velocity. For the problem considered, the internal baffle defined at $z=2\mathrm{m}$ (dashed yellow line in Figure 14b) also ensures that no liquid enters the injection opening as liner bridges the gap.

Another aspect of those kind of simulations is calculation of the initial magnetic target. It has been accomplished by developing a stand alone solver in OpenFOAM that finds Grad-Shafranov equilibrium (by solving Equation (23)) for an arbitrary geometry of the gas cavity. As an alternative, one can use any of already existing Grad-Shafranov solvers and then interpolate solution on the OpenFOAM mesh. To run a simulation where magnetic target is injected at some time though the implosion, the following procedure has been adapted:

• Start the simulation with low density and pressure gas inside the cavity (to minimize effect of the gas on the liner trajectory). Increase viscosity of the gas to damp waves bouncing in the cavity to obtain a smooth velocity field.
• Stop simulation at the time of the magnetic target injection. Calculate magnetic field for the target by finding Grad-Shafranov equilibrium using the “GradShafranovFoam” solver.
• Restart the simulation and follow magnetic target compression.

It is worth noting, that in reality, there is no gas in the cavity prior to plasma injection, and the gas velocity field is an artifact. Resetting gas velocity to zero at the time of magnetic target injection, caused waves in the cavity because of discontinuity in the velocity field. The best results have been achieved by restarting the simulation with the velocity field as it has been obtained prior to injection of the target. The assumption here is that the velocity field does not have a strong effect on the evolution of the magnetic field. The good agreement between “mhdCompressibleInterFoam” solver (run with high viscosity of the gas) and “layered VAC” (run as inviscid gas) provides some support for this assumption.

A simulation for a numerical setup shown in Figure 14 has been carried out on a uniform rectangular mesh with $1\mathrm{mm}$ size cells with simulation parameters listed in

Table 2. Case parameters for the setup including solid regions (see Figure 14).

	Geometry
Shaft radius	$R_{shaft}$	$0.15\mathrm{m}$
Inner liner radius	$R_{in}$	$1.5\mathrm{m}$
Outer liner radius	$R_{out}$	$2.11\mathrm{m}$
Height of the liner	H	$2\mathrm{m}$
Solid regions	$r<0.9\mathrm{m}$ , $r>1.3\mathrm{m}$, $2\mathrm{m}<z\le 3\mathrm{m}$, $-0.2\mathrm{m}<z\le 0\mathrm{m}$
Injection opening	$0.9\mathrm{m}\le r\le 1.3\mathrm{m}$, $2\mathrm{m}<z\le 3\mathrm{m}$
	Liquid liner (lithium)
Density	${\rho }_{\mathrm{lithium}}$	$500{\mathrm{kg}/\mathrm{m}}^3$
Viscosity	${\mu }_{\mathrm{lithium}}$	$3.645\times {10}^{-4}\mathrm{Pa}$
Magnetic diffusivity	${\eta }_{\mathrm{lithium}}$	$0.2{\mathrm{m}}^2/\mathrm{s}$
EoS	constant density
Pushing pressure	$p_{\mathrm{out}}$	$7\times {10}^6\mathrm{Pa}$
	Magnetic target
Initial pressure (uniform)	$p_{\mathrm{target}}$	$1\times {10}^4\mathrm{Pa}$
Initial density (uniform)	${\rho }_{\mathrm{target}}$	$0.11\mathrm{kg}/{\mathrm{m}}^3$
Viscosity	${\mu }_{\mathrm{target}}$	$0.35\mathrm{Pa}\mathrm{s}$
Magnetic diffusivity (constant)	${\eta }_{\mathrm{target}}$	$1{\mathrm{m}}^2/\mathrm{s}$
Initial magnetic configuration—Grad-Shafranov equilibrium (injected in $t=4.1\mathrm{ms}$)
Shaft current	$I_{\mathrm{shaft}}$	$2\mathrm{MA}$
Maximum poloidal flux (per radian)	${\psi }_{max}$	$0.1{\mathrm{Tm}}^2$
Shaft current constant	$F_o$	$0.4\mathrm{Tm}$
	Simulation type
Laminar, iso-thermal, high target viscosity, constant magnetic diffusivity

. Geometry and parameters of the liquid liner and magnetic target are similar to those in Section 3.1.3 (comparison between “mhdCompressibleInterFoam” and “layered” VAC), albeit magnetic target is injected during the compression. Fixed value of $\psi =0$ and a zero gradient for F boundary conditions have been applied at all boundaries. Slip wall boundary condition were set at internal baffles (yellow lines in Figure 14) and no slip wall boundary condition was imposed at the solid shaft-gas boundary. Results for evolution of the magnetic field during compression are shown in Figure 15 starting from the moment of magnetic target injection at $t=4.1\mathrm{ms}$. Contours of the poloidal function $\psi$ are shown in the first column within the range $0\le \psi \le 0.1\mathrm{Wb}/\mathrm{rad}$. Evolution of the toroidal function F is shown in the second column for the range of $0\le F\le 1.5\mathrm{Tm}$. By following the evolution of the $\psi$ contours, it can be seen, that poloidal field diffuses into the liner and central shaft as magnetic target gets compressed. From the plasma point of view, once a contour of $\psi$ touches the liner/wall, it is not considered to be part of plasma any longer, and one of the ways to identify the plasma, is to consider the area limited by the last closed flux contour inside the gas. Maximum value of $\psi$ at gas-liquid/solid interfaces is an indicator of the amount of flux lost by plasma during compression.

By following the evolution of F contours, diffusion of the toroidal magnetic field into the liquid liner and solid regions is observed. Electrical poloidal current is flowing in the liner and solid regions where toroidal magnetic field is diffused. It is interesting to note, that the toroidal magnetic field diffused into a bottom/top solid regions at the early stages of the compression remains there, while the liquid liner continues to implode (e.g., Figure 15f), such that the electrical current continues to flow in the solid and in the thin layer at vicinity of the liquid liner-solid interface. Amplification of the toroidal flux as it gets compressed is also seen in the figure. By calculating amount of toroidal flux in the liner and solid regions, portion of the flux being diffused during compression can be estimated.

A setup shown in Figure 15, can be further extended by imposing some magnetic field configuration in the liner and solid regions prior to plasma injection. As an example, a vertical uniform magnetic field that is fully soaked into the liner and solid regions is now a starting point for the simulation. Initial distribution of the poloidal function $\psi$ corresponding to the uniform vertical field is then given by:

$$\psi (t=0,r)=0.5B_o{r}^2,$$

(29)

where $B_o$ is strength of the vertical magnetic field. As in the previous case, simulation is first run until the inner surface of the liner reaches the outer radius of the injection opening. Magnetic target is then calculated by solving the Grad–Shafranov equation for the gas cavity domain. Boundary condition for the poloidal flux function is specified using poloidal flux function field obtained in “mhdCompressibleFoam” simulation at the time of target injection. Example of such a simulation is shown in Figure 16 starting from the moment of plasma target injection. Geometry and parameters for the simulation are kept the same as in the previous case (see Table 2) and the strength of initially imposed vertical magnetic field is $|B_o|=0.1\mathrm{T}$. The poloidal flux function $\psi$ corresponding to the uniform vertical magnetic field and that of magnetic target should be of opposite signs, i.e., for a positive ${\psi }_{\max }$ of the target there is a negative $\psi$ (opposite direction) of the vertical field.

In Figure 16 the contours of the poloidal flux function and corresponding toroidal function are shown in the left and right columns, similar to the Figure 15. Magnetic target configuration at the moment of injection is shown in Figure 16a for the range of $-0.23\le \psi \le 0.1\mathrm{Wb}/\mathrm{rad}$. Negative values of $\psi$ correspond to the uniform vertical field and $\psi =0.1\mathrm{Wb}/\mathrm{rad}$ is equal to the initial ${\psi }_{\max }$ of the magnetic target. The magnetic plasma target is limited by the last closed $\psi$ contour, such that the vertical magnetic field provides a buffer between the target and inner surface of the liner. One can also see that the vertical field remains frozen in the liner, such that magnetic field lines bend as liner implodes, and, depending on the strength of the field, may deform the liner surface. For the current set of parameters, a target remains separated from the liner surface because of a buffer created by the vertical magnetic field throughout the entire compression. Evolution of the toroidal field is similar to the case without vertical field, e.g., compare right columns in Figure 15 and Figure 16. In the future, different initial configurations of the poloidal field obtained in the experiments can be used as a starting point for simulations to investigate their effect on the liner dynamics and plasma compression characteristics.

3.3. Including Rotation of the Liquid Liner into Simulations

In all test cases considered in the previous sections, rotation of the liquid metal liner has not been included into simulations. In the real compression system, however, the liquid liner is rotating at a certain speed prior to implosion. In fact, the initial cylindrical cavity is formed due to the rotation of the liquid metal. In addition, sufficient rotation of the liner helps to stabilize its inner surface against Rayleigh-Taylor instability that might potentially develop deep into compression [11]. To extend the previous test case to encounter for the rotation of the liner, it is important to correctly specify starting conditions for the simulation based on what is expected to happen in the real system. Going back to the numerical setup in Figure 14, but with an initially fully soaked vertical uniform magnetic field and a rotating liner, we now have stationary and rotating electrically conductive parts in the system. Therefore, if those parts are in electric contact, toroidal field and accompanied induced currents are to be generated in the system prior to the implosion (due to the last term in Equation (3)). Example of such currents is shown in Figure 17. Here it is assumed that a rotating liquid liner has electric contacts with top and bottom stationary regions prior to the liner implosion. As such, there is a discontinuity in rotational speed $\mathsf{\Omega }$ at the interfaces. In the case of soaked poloidal field this results in generation of the toroidal field in the vicinity of the interfaces. Developed toroidal field configuration also depends on the boundary conditions. In the example shown in Figure 17, it was assumed that the value of toroidal function is zero at the solid region boundaries, i.e., it is an insulator. To obtain the result shown in Figure 17, a steady state simulation has been run with imposed constant rotation speed of the liner. In general, this type of simulations is required to estimate torque etc. due to the induced currents and to generate initial magnetic field configuration in the system. It is worth noting, that the generated toroidal field is of the opposite signs in vicinity of the upper and lower interfaces between rotating liner and stationary solids. This means that effect of the induced currents on the liner will be different near the top and bottom of the liner.

Example of compression simulation that includes both rotating liner and stationary solid regions is shown in the Figure 18. In this simulation, an initially fully soaked vertical magnetic field with strength $|B_o|=0.05\mathrm{T}$ is compressed by the liquid liner initially rotating at $\mathsf{\Omega }=40\mathrm{rad}/\mathrm{s}$. Here only compression of a vacuum magnetic field is considered, i.e., no magnetic target is injected at any stage of the compression. To represent vacuum, magnetic diffusivity in the cavity (and injection path) has been set to ${\eta }_{\mathrm{vacuum}}=10{\mathrm{m}}^2/\mathrm{s}$. All other parameters and geometry remain the same as in the previously considered non-rotating liner (see Table 2). Two different scenarios are considered. In the first one, left column of the figure (Figure 18a,c,e,g), the rotating liner and top and bottom stationary regions are assumed to be in electrical contact prior to the implosion (situation previously discussed and shown in Figure 17), whereas, in the second one, right column of the figure (Figure 18b,d,f,h), rotating and stationary parts are not in an electrical contact initially, but they come into contact, once the liner starts to implode. Contours of both poloidal flux function $\psi$ and toroidal function F are shown in Figure 18. Starting conditions for both simulations are shown correspondingly in Figure 18a,b. Uniform vertical field is shown by the contours of the poloidal function $\psi$ (straight vertical lines) and is identical in both simulations. As previously discussed, a toroidal magnetic field (induced currents) is initially present in the first scenario (Figure 18a), whereas, in the second scenario no induced currents are initially present (Figure 18b).

First of all, it is worth noting that for the simulations presented in Figure 18, implosion is notably slower than in the previously considered cases (e.g., Figure 15 and Figure 16), despite a similar set of parameters. This difference is due to the initial rotation of the liner, as some part of the pushing pressure (that is kept identical in all cases) goes to counter the rotation, and, for the rotation speed of $\mathsf{\Omega }=40\mathrm{rad}/\mathrm{s}$, this part is not negligible. Several observations can be made inspecting Figure 18: (i) Both toroidal and poloidal components of the magnetic field initially present in the liner, move with the liner, leading to high gradients in the vicinity of the interfaces between moving and stationary parts. Those gradients, in turn, are directly related to the Lorentz force acting on the liner. This might have a noticeable effect on the rotational speed of the liner near the interfaces. If to zoom in onto top corner of the inner liner surface for both simulations (compare for example Figure 18e,f), it can be seen that the position of this corner is not exactly the same. This can only be due to the different distribution of induced currents, as everything else is identical between the two simulations. (ii) Depending on the assumptions made with regard to the electrical contact between rotating liner and stationary parts, liner surface shape during the implosion might be noticeably different, and the effect is more pronounced with an increase in rotation speed and strength of the initial poloidal field. (iii) Shape of the liner deviates from initial cylindrical shape deep into compression (e.g., Figure 18g), therefore if the exact shape of the liner is of importance, initial magnetic configuration needs to be treated with care.

3.4. Example of Magnetic Target Compression in General Fusion Inc. Fusion Demonstration Plant (FDP) like Geometry

In this section compression simulation of the magnetic target in General Fusion FDP like geometry is briefly presented. The methodology is similar to that discussed in the previous sections; the liner starts to implode (with or without initially soaked magnetic field) and when its inner surface reaches the designated position, magnetic target is injected into the cavity and subsequently compressed. In FDP type compression, an initially cylindrical cavity evolves towards a spherical shape during the implosion. Almost all of the liquid metal required to form the imploding liner is initially stored in the bores arranged into a number of rows in the vertical direction (see Figure 1, [1]). The entire vessel is rotating prior to the implosion. Highly pressurized gas pushes on the liquid in the bores. This forces liquid metal to flow out of the bores into the chamber forming the rotating liquid liner. Shaping of the liner is achieved by varying pressurized gas pushing profiles applied to each row of the bores in such a manner, that the initially cylindrical inner surface of the liner transforms towards a spherical one as liner implodes.

A simplified example of such a compression is shown in Figure 19, starting from the time of target injection. Solid regions and injection opening are omitted for simplicity, and the target is injected at the moment when the upper corner of the inner surface of the liner reaches radial position of $r=1.3\mathrm{m}$ (initial radius of a cylindrical inner surface for this case is $r=1.5\mathrm{m}$). Simulations have been performed for both without and with initially fully soaked vertical field. Initial conditions for the target compression simulations are shown in Figure 19a for $|B_o|=0$ and Figure 19b for $|B_o|=0.05\mathrm{T}$. Here it is assumed that a uniform vertical field is imposed at the moment of target injection, therefore, vertical field contours are the straight vertical lines to begin with. It can be seen that the inner surface of the liner deviates from the cylindrical shape at the moment of injection because of different pressure pulses applied at different bores, i.e., higher pressure near the top and bottom and lower pressure with a slight time delay at the center. Magnetic target configuration has been calculated by solving Grad-Shafranov equation for the curved shape of the cavity. As in the previous tests, magnetic configuration has been defined by Equation (24), with $F_o=0.4\mathrm{T}\mathrm{m}$ ($I_{\mathrm{shaft}}=2\mathrm{MA}$) and ${\psi }_{max}=0.1\mathrm{Wb}/\mathrm{rad}$. Injected magnetic target are shown by $\psi$ contours with the last closed $\psi$ surface emphasized with the green color. Volume inside the green line is what is considered to be plasma. Initial rotational speed of the liner has been set to $\mathsf{\Omega }=34\mathrm{rad}/\mathrm{s}$ with liner properties set to that of lithium. Generation of toroidal field in the liner due to contact between rotating liner and stationary regions, however, has been omitted in those simulations. Constant values of magnetic diffusivity have been used for gas and liquid metal phases (${\eta }_{\mathrm{gas}}=1{\mathrm{m}}^2/\mathrm{s}$, ${\eta }_{\mathrm{lithium}}=0.2{\mathrm{m}}^2/\mathrm{s}$). Simulations have been run using uniform $2\mathrm{mm}$ resolution mesh. Boundary conditions are the same as in the previous cases. Following target compression in FDP geometry presented in Figure 19, one can observe shaping of the liner’s inner surface towards a more spherical one during implosion, trapping the plasma target in the center of the cavity. Following the last closed flux surface denoted by the green line, it can be seen that shape of the plasma target can be manipulated to some extent by imposing uniform vertical field prior to the compression. Wiggles developed on the contours of the $\psi$ associated with uniform vertical magnetic field in the liner (middle column), are due to boundary conditions at the bores walls. Initially vertical lines are deformed as $\psi$ field values remain frozen at the bore walls, while liquid is pushed along the bore. The next step will be to extend the current setup to include part of surrounding electrically conducting structure and simulate more realistic configurations.

4. Summary and Conclusions

An “mhdCompressibleInterFoam” solver for simulating electrically conducting compressible and immiscible fluids using the Volume of Fluid (VOF) has been developed based on the “compressibleInterFoam” solver within OpenFOAM. The new code has been verified using several reference test cases for which analytical/reduced order solutions exist and also by comparing results with those obtained with VAC (Versatile Advection Code) for a more complex problem of compressing magnetic target by an imploding liquid metal liner in a cylindrical geometry. Methodology to include stationary regions to mimic parts of conducting solid structure surrounding liquid metal liners has been proposed and applied to a simplified setup relevant to the MTF scheme for FDP developed by General Fusion Inc. So far the new solver appears to be robust and the results are encouraging.

The next step is to apply “mhdCompressibleInterFoam” to a wide range of relevant problems which include conducting fluids and solids to better understand solver capabilities and limitations. The main purpose of this solver is to bridge the gap between pure hydrodynamic simulations of liner dynamics and detailed simulations of plasma evolution which are conducted using dedicated plasma codes. One of the crucial assumptions here, that evolution of magnetic fields is not overly sensitive to the fine details of plasma evolution. Preliminary data, coming from comparison between “mhdCompressibleInterFoam” results and those produced with VAC code using a much more sophisticated plasma model, indicates that it is likely to be the case (at least for some range of parameters). If the assumption holds, then to study interaction of magnetic fields with conducting liquid and solid using a simplified plasma model is justified.

Stability of the inner surface of the liner during implosion is also of great importance for the overall success of MTF scheme in pursuit by General Fusion. Hydrodynamic interface instabilities and means of their suppression have been extensively studied in General Fusion by combination of numerical simulations and linear stability theory. We are optimistic that with our new “mhdCompressibleInterFoam” solver some of the effects of the magnetic field on the stability of the interface can be studied (for the case of $n=0$ mode). Unfortunately, $n=1$ and higher modes instabilities of the liner surface that might develop because of $3D$ plasma instabilities are out of reach for the current version of the code.

The main limitations of the current solver are: (i) solver can only be applied to the axi-symmetric problems. Therefore, investigation of any three-dimensional instabilities (in plasma or liner surface) is not possible for the time being; (ii) only simplified magnetic target model has been implemented and tested. As such, detailed plasma evolution characteristics (such as temperature) can not be simulated at present; (iii) VOF approach, as it is currently implemented in OpenFOAM, is numerically difficult for high density ratios between fluids. This puts a limitation on the kind of plasma that can potentially be simulated. For example, simulations of low density plasma are likely to be out of reach even if more complex plasma models are implemented.

With respect to the future solver development, an implementation of more sophisticated plasma models is something we will be investigating. Extension of the solver to full $3D$ is also being considered.