Extended-MHD simulations of disruption mitigation via massive gas injection in SPARC

The initial conditions for the extended-MHD simulations are given by the SPARC PRD, which is a projected double-null diverted equilibrium with DT fueling in H-mode, maximizing fusion gain [24]. The geometry as well as pressure and safety factor profiles are shown in figure 1. In this paper $\psi_\mathrm N$ denotes the normalized poloidal flux. While this discharge is projected to operate in H-mode, the profiles do not have a pedestal. The plasma has a geometric major radius of $R_0 = {1.85}$ m with the magnetic axis being located at $R_{\mathrm m} = {1.89}$ m and the LCFS closely follows the inner wall. The central safety factor value is below one, and thus the plasma is prone to developing a sawtooth instability. The initial plasma current is $I_\mathrm p = {8.8}$ MA at an on-axis magnetic field of B₀ = 12.2 T, and the initial thermal and in-vessel poloidal magnetic energies are $W_\mathrm T = {20}$ MJ and $W_\mathrm M = {70}$ MJ, respectively.

Zoom In Zoom Out Reset image size

Nonlinear 3D simulations are performed with the initial value code M3D-C1, which implements an extended-MHD model derived from the Braginskii equations [25]. In particular, the extended-MHD model used in our simulations is given in terms of the equations

$$\begin{align} \frac{\partial n}{\partial t} + \nabla \cdot \left(n \mathbf{u}\right) & = \Sigma\,,\nonumber\\ \frac{\partial n_j}{\partial t} + \nabla \cdot \left(n_j \mathbf{u}\right) & = \sigma_j + D \nabla^2 n_j\,,\nonumber\\ n m_i \left(\frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} \right) & = \mathbf{J} \times \mathbf{B} - \nabla p - \nabla \cdot \Pi \nonumber\\ &\quad+ \mathbf{F} - m\mathbf{u} \Sigma\,,\nonumber\\ \frac{\partial p}{\partial t} + \mathbf{u} \cdot \nabla p + \Gamma p \nabla \cdot \mathbf{u} & = \left(\Gamma - 1\right) \Big[Q - \nabla \cdot \mathbf{q} + \eta J^2 \nonumber\\ &\quad- \mathbf{u} \cdot \mathbf{F} - \Pi : \nabla u - \frac{1}{2}m_i v^2 \Sigma\Big]\,,\nonumber\\ \mathbf{E} & = -\mathbf{u} \times \mathbf{B} + \eta \mathbf{J}\,,\nonumber\\ \mathbf{J} = \frac{1}{\mu_0} \nabla \times \mathbf{B},&\quad \frac{\partial \mathbf{B}}{\partial t} = - \nabla \times \mathbf{E}\,, \end{align} \tag{ 1 }$$

with $\Sigma : = \sigma + D \nabla^2 n$ and where, as usual, B denotes the magnetic field, p the pressure and n the main ion density. For an impurity with atomic number $Z\,n_j$ denotes the density for each charge state with $1 \unicode{x2A7D} j \unicode{x2A7D} Z$. Furthermore, u is the fluid velocity, J the current density, E the electric field, Γ the ratio of specific heats, m_i the main ion mass, η the resistivity, Π the viscous stress tensor, F are external forces, Q is a heat source term, σ an ion source term and D denotes particle diffusivity. The heat flux density is $\mathbf{q} = - \kappa \nabla T - \kappa_{\|} \frac{\mathbf{B} \mathbf{B}}{B^2} \nabla T$ with isotropic thermal conductivity κ, parallel thermal conductivity $\kappa_{\|}$ and temperature T. For further information on the M3D-C1 model equations see [15]. This system of equations is solved in 3D geometry using reduced quintic C¹ finite elements within the poloidal plane and cubic finite elements in toroidal direction [26, 27] with semi-implicit time stepping on an unstructured mesh with triangular elements in the poloidal plane and prism-shaped elements toroidally (see section 2.3 for more details on the mesh). Plasma resistivity is calculated according to the well known Spitzer resistivity model. Equilibrium plasma rotation can in principle be included, but is not considered in this study as it might have introduced further numerical challenges. In our model the stress tensor Π includes the effect of finite viscosity, which has a direct influence on plasma dynamics. The simulations are carried out fully three-dimensional as explained in subsection 2.3. The code uses a cylindrical $(R,\phi,Z)$ coordinate system, where R is the major radius, φ the toroidal angle and Z the vertical distance w.r.t. the midplane.

Atomic processes such as ionization, recombination and radiation are included via the KPRAD module [22], which has been coupled with M3D-C1 [15]. This model does not assume a coronal equilibrium. Instead ionization and recombination determine the impurity charge state evolution. The radiated power, associated with charge state k is calculated as $P_{\mathrm{rad},k} = n_k n_\mathrm e L_k$, where n_k is the ion density of charge state k, $n_\mathrm e$ is the electron density and L_k is the radiation rate, which is taken from the ADPAK database [28]. Within this model loss of thermal energy can happen due to ionization, line radiation, bremsstrahlung and recombination radiation. Neutral impurities are injected at zero temperature and brought to the ion temperature $T_\mathrm i$ upon ionization. KPRAD uses an explicit time step with faster subcycling than a typical MHD time step such that multiple KPRAD time steps are carried out during a single advance of the MHD time step. Changes to the densities of individual ion charge states based on recombination and ionization rates are calculated within KPRAD and enter the M3D-C1 equation through the source/sink terms in the set of equations (1). Advection is calculated within M3D-C1 by solving a PDE during the MHD time step. More details on the coupling between M3D-C1 and KPRAD can be found in [15].

M3D-C1 uses an unstructured triangular mesh that can be divided into three distinct areas: (1) a plasma region, where the set of extended-MHD equations (1) is being solved. This region covers the plasma and extends beyond the last closed flux surface (LCFS) to the inner wall. In the region between the LCFS and the first wall, the plasma is considered to have very low pressure and temperature and thus very high resistivity. (2) A resistive wall region. (3) A vacuum region without plasma, but finite electromagnetic fields. Recent developments to the M3D-C1 code enable a more accurate representation of device structures. Unlike in previous studies of disruptions and their mitigation, the wall region is not modeled as a single thin layer, but consists of multiple shells comprised of spatially resolved conducting elements with realistic resistivity values. These three regions of the computational domain together with the device structures that are included in the M3D-C1 model are illustrated in figure 2. Including more conducting elements than only the inner wall allows us to better identify where eddy currents are induced inside the device and enables a more accurate calculation of the wall forces. SPARC has rows of ports at three poloidal locations on the outboard side. These ports are represented in the model in terms of axisymmetric regions of anisotropic resistivity, where the wall material's actual value of resistivity is applied to poloidal currents, but toroidal currents are suppressed by applying vacuum resistivity to the toroidal component of the current density within the port area. At the time when the simulations were carried out the SPARC design featured two passive plates, which are included in the wall model (see figure 2). However, these passive plates have since been removed and are no longer part of the SPARC design. The simulations presented here are fully three-dimensional and use prism-shaped elements in between toroidal planes. The number of toroidal planes used in a particular simulation depends on the number of toroidal injector locations and the toroidal spread of the gas plume, but ranges in between 12 and 24 for the cases presented here. M3D-C1 does not require these toroidal planes to be equidistant, and in order to better resolve the dynamics around the gas jets and to resolve the strong gradients, toroidal planes are packed around the injector locations as shown in figure 3. During the six injector simulation we increased the number of toroidal planes from 18 to 24 as the injected gas plume was moving toroidally. Two planes were added toroidally at each side of the injector sites to ensure that the strong gradients associated with the injected gas can still be resolved.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

One configuration that was considered for the SPARC DMS layout is six gas injector valves at three toroidal locations with an upper and lower gas valve. In the simulations we consider up-down symmetric configurations with six, four and two injectors, as well as up–down asymmetric configurations with one and five injectors. Within the poloidal plane the ends of the MGI barrels are modeled at R = 2.18 m and $Z = \pm{0.68}$ m. The toroidal location varies depending on the number of injectors. When one or two injectors are used the gas is injected only at a single toroidal position. With four injectors the toroidal spacing is 180^∘, with five and six injectors the toroidal spacing between injectors is 120^∘. The injected gas is modeled in terms of a gas plume with Gaussian shape and is toroidally localized. Smaller gas plumes provide a more realistic representation of the injected gas jets and thus allow for a more physical prediction. However, such small impurity plumes are more difficult to resolve in the simulations and require very small time steps in addition to using higher spatial resolution. In the process of finding plume sizes that would be as realistic as possible, but could still be resolved we adopted a half width of the Gaussian of 8 cm in poloidal direction and 80 cm in toroidal direction in the simulation with two injectors. In the six injector case we were able to reduce it to 4 cm and 24 cm, respectively. Figure 3 illustrates the injector layout in the simulations using the full six injector layout as an example. The gas plume has zero initial velocity, but is advected by the bulk plasma. At the time when these simulations were carried out the projection was for SPARC to use a mixture of neutral Ne and ionized D in a ratio of 1:10. It is expected that $\sim{10^{22}}$ Ne atoms must be delivered in a gas pulse with a duration of 2 ms for successful mitigation of the SPARC PRD. An injection rate of ${7.5 \times 10^{23}}$ particles of Ne per second is used which gives near the full delivery with six valves. In our simulations, the impurities are injected beginning at t = 0 at the nominal injection rate.

The SPARC PRD equilibrium exhibits an on-axis safety factor below unity with a rational q = 1 surface at $\psi_N = 0.25$, as seen in figure 1(b). This leads to the onset of a sawtooth instability early in the simulations. At the same time, the presence of the gas plumes due to MGI induces a low-n MHD perturbation, where n is equal to the number of toroidal injection locations, i.e. n = 1 for two injectors and n = 3 for six injectors. The magnetic energy associated with different toroidal modes in these two cases is shown in figure 4. With only two injectors the n = 1 mode is dominant as expected from both, the sawtooth, and also the n = 1 perturbation introduced by the gas injection in a single toroidal plane. Mode growth is slow during the pre-TQ and amplitudes peak during the TQ. In the six injector case the n = 3 mode due to gas injection grows first and shortly after it saturates the n = 1 mode becomes dominant. Here, the pre-TQ is defined as the period between the start of injection and the time when the thermal energy $W_\mathrm {th}$ starts to drop significantly. During the TQ a strong n = 2 perturbation develops together with a broad band of weaker modes at higher n. Figure 5 shows the toroidal current density plotted in the injection plane at $\phi = {60}^{\circ}$ in the simulation with six gas injectors. In the pre-TQ and early TQ edge perturbations are seen that arise due to the presence and toroidal periodicity of the injected impurities. In these phases toroidal eddy currents are formed in the inner vacuum vessel in between the port regions, where such currents are suppressed due to anisotropic resistivity. In the late TQ the plasma becomes stochastic, which coincides with the onset of a broad spectrum of toroidal modes seen in figure 4(b).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The initial stored energy of this SPARC PRD equilibrium is $W_\mathrm {th} = {20}$ MJ. In the simulations this energy is lost via radiation due to atomic processes as modeled with KPRAD, but also through thermal conduction through the first wall (red line in figure 2). The amount of energy that is lost via conduction strongly depends on the value of thermal conductivity κ that is imposed on the plasma. Within the extended-MHD model κ is not calculated self-consistently, but provided based on empirically determined values. Throughout most of the simulations and where not otherwise stated we use a value of $\kappa = {1.54\times 10^{21}}\,{\mathrm{m}^{-1}\, \mathrm{s}^{-1}}$ for isotropic thermal conductivity, while parallel thermal conductivity is set to a value of $\kappa_{\|} = {1.54\times 10^{26}}\,{\mathrm{m}^{-1}\, \mathrm{s}^{-1}}$. Figure 6 shows the temporal evolution of the total radiated power in all considered injector scenarios ranging from one to the maximum of six injectors⁵. Small sudden bursts of radiation are seen in the pre-TQ throughout all cases. While it is not entirely clear what causes these bursts, one explanation could be the cold injected gas reaching rational flux surfaces where it is rapidly transported along the field lines and thus able to radiate promptly. Nevertheless, as can be seen from figure 6(b) these radiation spikes do not lead to a noticeable loss of thermal energy, and thus do not matter for the overall energy balance. In the pre-TQ the loss of thermal energy is seen to scale with the number of injectors and thus also the total rate of injected impurities. The onset of the sawtooth starting at $\lt {1}$ ms also does not correlate with an increase in radiated power. This can be explained by a lack of impurity mixing, as the impurities were not transported to the core and remain at the plasma edge at this early point in time. It can be seen in figure 6(b) that the one and five injector cases loose their thermal energy more rapidly than the other cases. This is a result of numerical difficulties in these simulations that required the isotropic thermal conductivity to be raised to a value of $\kappa = {1.39\times 10^{22}}\, {\mathrm{m}^{-1}\, \mathrm{s}^{-1}}$ during the early TQ to get converged results. This increased the loss of energy by means of conduction through the first wall. These difficulties were not seen in the two and six injector simulations, where a radiative TQ is observed. In the two cases at the full nominal injection rate, i.e. two injectors and six injectors the initial stored thermal energy is lost within $\approx {4}$ ms after injection started. With six injectors, however, a small amount of energy is lost through conduction in the pre-TQ as a result of imposing a temporarily larger κ in this phase. Despite starting the actual TQ from a lower $W_\mathrm {th}$ with six injectors the TQ finishes about 0.2 ms later than in the two injector case. With 6 injectors it takes the plasma slightly longer to loose all of its thermal energy compared with 2 injectors and $P_\mathrm {rad}$ peaks at a later time, as seen from figure 6. At the end of the TQ radiated power is still large, which is due to a conversion of Ohmic heating to radiation. In both cases more than 95% of the energy are radiated, while only about 2.1 MJ are conducted through the wall.

Zoom In Zoom Out Reset image size

Towards the end of the TQ an increase in plasma current is observed between 3 ms and 4 ms just before the start of the CQ. The toroidal current as a function of time is shown in figure 6(c). In both scenarios the CQ is very fast with a duration of only about 2.2 ms in the case of two injectors and 2 ms with six injectors, respectively. The observed CQ is faster than the prediction of the ITPA scaling [29], which yields a minimum CQ duration of 3 ms for the given SPARC equilibrium. However, we find that $T_\mathrm e$ drops rapidly with the onset of the CQ and during the second half of the CQ $T_\mathrm e$ reaches values below 6 eV, which is smaller than the values of $T_\mathrm e$ that were seen in experiments for example on EAST [30]. The rapid drop of $T_\mathrm e$ observed in the simulation together with the low value of $T_\mathrm e$ towards the end of the CQ would lead to increased plasma resistivity according to the Spitzer resistivity model and could explain the fast CQ. M3D-C1 modeling of KSTAR found shorter CQ durations than seen empirically [31]. Further validation of the predicted CQ durations in M3D-C1 is a topic of future works. At 6 ms, when the toroidal current is close to zero, the KPRAD module was switched off and at the same time the plasma viscosity was increased drastically in M3D-C1 to terminate the plasma dynamics with the aim to speed up the simulation and calculate the wall forces in the post-CQ.

There is no substantial vertical movement of the plasma during the rapid shutdown, i.e. we observe neither a hot or cold vertical displacement event (VDE). This is shown in figure 7 for the simulations with two and six injectors. The quantity shown here is the vertical position of the current centroid, which can be challenging to determine depending on the magnetic geometry. With six injectors this code diagnostic indicates a maximum vertical displacement of about 10 cm, which is smaller than the vertical displacement of 24 cm with only two injectors. However, Poincare plots show an even smaller vertical displacement of only a few cm in both cases. The absence of a VDE is beneficial for the forces on the plasma vessel as eddy currents as well as halo currents associated with the vertical plasma motion are absent.

Zoom In Zoom Out Reset image size

Overall, the forces on the vessel structure are moderate as seen in figure 8, which shows the horizontal and vertical axisymmetric force, as well as the non-axisymmetric sideways force integrated over the wall region. The strongest force is the axisymmetric radial force, which peaks at 40 MN during and after the CQ and points inwards in both injector cases. After the peak, the force is seen to slowly decay, which is consistent with the projected wall time of ${\sim}{50}$ ms on SPARC. The axisymmetric vertical force is two orders of magnitude smaller than the radial force, and thus differences between the two injector and six injector cases can be ignored. The non-axisymmetric sideways force (in horizontal direction) differs in the two considered injector scenarios. With only two injectors the force reaches peak values of above 2 MN, whereas with the full six injector layout the force is only about half of this value. However, the sideways force is one order of magnitude smaller than the axisymmetric radial force and thus in terms of the wall forces the six injector case performs only marginally better.

Zoom In Zoom Out Reset image size

In the simulations it is seen that the injected impurities stay localized around the sites where the gas was injected. Figure 9 shows electron density $n_\mathrm e$ and radiated power $P_\mathrm {rad}$ at the plane of injection at different times during the simulation. While the size of the impurity cloud increases somewhat in size over time, there is no considerable radial, toroidal nor poloidal transport of the impurities throughout the TQ. As mentioned before, impurities are injected at zero temperature and then heated by the bulk plasma. Since the simulations use a single-fluid model, there is no ambipolar electric field resulting from electron pressure gradients, which might drive propagation of impurities to the surrounding plasma [32]. The instantaneous temperature equilibration causes a local decrease in temperature and pressure at the injection sites thus keeping the injected impurity atoms from spreading. Furthermore, toroidal equilibrium rotation is not considered in these simulations as simulation time was limited and the incomplete physics basis on rotation braking following MGI challenged the development of a reduced model. If finite equilibrium rotation was included, the impurities would be expected to better distribute along the flux surfaces in toroidal direction leading to less accumulation and also lower toroidal peaking as is discussed in the next section.

Zoom In Zoom Out Reset image size

To estimate the toroidal spread of radiation and to get an estimate on peak radiative heat loads on the first wall we calculate the TPF for radiated power. At a given time t the TPF can be defined as the peak value of the radiation field integrated in the poloidal plane at any toroidal location divided by the toroidal average

$$\begin{align} TPF = \frac{\left[ \iint P_\mathrm {rad}\left(R,\phi,Z\right) \mathrm{d} R \mathrm{d} Z \right]_\mathrm {peak}}{\frac{1}{K} \sum_k \iint P_\mathrm {rad}\left(R,\phi,Z\right) \mathrm{d} R \mathrm{d} Z}\,,\end{align} \tag{ 2 }$$

where $P_\mathrm {rad}$ is the radiated power, K is the number of equidistant toroidal planes used to evaluate the TPF and k is an index that labels these planes. Note that the toroidal planes in the above equation are different from the (non-equidistant) toroidal planes in the M3D-C1 simulations. The finite elements used in M3D-C1 allow the radiation field to be evaluated at any toroidal angle and the calculation of the TPF is performed using equidistant toroidal angles with K = 16 and K = 24 for the two and six injector cases, respectively. These values are chosen to ensure that the toroidal angles where the radiation field is evaluated coincide with the locations of impurity injection. To evaluate radiation distribution in the different injector cases, we will calculate the TPF as a function of time to determine its overall maximum value, but also its value at the time of peak radiated power. Note that this method does not determine the magnitude of heat loads across the surface area of the first wall, which can be calculated via ray tracing and is beyond the scope of this work. However, the calculated radiation fields can be used as input to ray tracing codes to map out the radiation loads at each point of the wall. Calculating the TPF is nevertheless an appropriate measure to evaluate the effect of different injector configurations on the toroidal peaking of radiation.

Figure 10(a) shows the TPF as a function of time in simulations with two and six injectors at the same nominal injection rate. At the beginning of the injection at t = 0 the TPF is very large, as the impurities are most localized (and so is radiation) at this moment. Since the total amount of impurities in the plasma is very small at this early stage, $P_\mathrm {rad}$ is low enough for the TPF to not be of concern. As can be seen by comparing figures 4 and 10 the TPF is not correlated with MHD activity. This can be explained by poor impurity mixing into the plasma core. The TPF drops rapidly and remains mostly at values of $2-4$ during the pre-TQ and early TQ. After a short increase to a value of 4 during the TQ in the six injector case, the TPF reduces to values of about 1.5–2.5 at the time when radiated power peaks and radiated power would be the most damaging. When only two injectors are used the TPF is considerably larger at a value of 7, which is reached when the total radiated power peaks. Hence, using the six injector configuration leads to more equally distributed heat loads than using only two injectors. Heat loads are most dangerous to plasma-facing components when they are peaked at a moment when the radiated power is large. Thus, the product of these two quantities provides a good estimate of how dangerous the heat loads are in the different injector configurations. This is illustrated in figure 10(b), and it can be seen that with only two injectors the heat loads are more concentrated and six injectors perform more favorably. We recall that in these simulations the impurities are advected by bulk plasma flow only. Since the impurities are injected at zero temperature there is no local increase in pressure and with the absence of an ambipolar electric field the TPF is overestimated and will likely be lower. This effect could be included in two-fluid simulations, which are reserved for future work. In addition, equilibrium toroidal rotation is zero. Another contributing factor to a high TPF is the zero initial velocity of the injected impurities, which keeps the impurities from further spreading radially into the plasma core. All of these effects would lead to a stronger motion of impurity particles away from the injector sites and a further reduction of the TPF could be expected when these effects were included. This should be subject of future simulations.

Zoom In Zoom Out Reset image size