Modelling of NSTX hot vertical displacement events using $M3D-C1$

Of the many reasons a plasma discharge disrupts,¹ Vertical Displacement Events (VDEs) lead to the most severe forces and stresses on the vacuum vessel and Plasma Facing Components (PFCs), especially in the presence of toroidal asymmetries and rotation.² After loss of positional control, the plasma column drifts across the vacuum vessel and comes in contact with the first wall, at which point the stored magnetic and thermal energy is abruptly released. Hot VDEs, where a high fraction of the pre-disruption energy is released into the wall, are more damaging compared with cold VDEs where the thermal and current quench precedes the vertical instability (often causing it) and a large fraction of the plasma energy dissipates prior to hitting the vessel.

Vessel forces have been extensively modelled in 2D but, with the constraint of axisymmetry, the fundamental 3D effects that lead to toroidal peaking, sideways forces, field-line stochastisation, and halo current rotation have been vastly overlooked. Predictive modelling capabilities^3–10 and engineering tools^11–15 are under development. The highly nonlinear multi-physics of disruptions is tackled via increasingly demanding numerical computations. The modelling of VDEs is particularly challenging because the induction and advection terms are as important as conduction and diffusion terms in the solved equations. There is no general method for solving these systems for arbitrary parameters; the range of applicability of all algorithm is limited. The most common approach to VDE simulations is to impose the evolution of fields such as plasma current or plasma position^16–20 and solve for the remaining fields. Given the stiffness of the underlying equations, there is a risk that forcing the dynamics spoils the predictive capability of those numerical tools and misleads physical interpretation of their results. One would hope to be able to model VDEs without interfering with their natural evolution.

In this work, we present the main results of an intense VDE modelling activity using the implicit 3D extended MHD code M3D-C1 and share our experience with the multi-domain and highly nonlinear physics encountered. At the culmination of code development by the M3D-C1 group over the last decade, highlighted by the inclusion of a finite-thickness resistive vacuum vessel within the computational domain,²¹ a series of 3D nonlinear simulations are performed using realistic transport coefficients.

The modelling is based on the NSTX shot #139536, a well-diagnosed frozen VDE,²² that has previously been simulated with the $M3D$ code.²³ The sequence of events during this hot VDE can be summarised from experimental traces as follows. At 300 ms into the discharge, the feedback control system is partially deactivated in order to trigger a vertical instability. The current centroid drifts towards the wall in about 50 ms, displaying a largely exponential time trace. The total toroidal current remains constant during the vertical drift phase but, as soon as the plasma comes in contact with the wall, the 600 kA plasma current decays in about 5 ms. The toroidal wall current rises to almost 400 kA and decays over a 10 ms period soon after. The NSTX shunt tiles²⁴ measure, during the current quench only, an axisymmetric normal wall current of 40 kA and an n = 1 component of 20 kA. These experimental figures of merit serve as modelling targets.

The $M3D-C1$ simulations are initialised with the $EFIT$²⁵ data from NSTX shot #139536 at 300 ms. The experimental PF currents are used as initial input for $M3D-C1$'s Grad-Shafranov solver; the code slightly adjusts the PF currents through Picard iterations to obtain exact force balance on the finite element mesh. The equilibrium obtained, and thus the vertical instability, depend weakly on the mesh. The external field produced by the PF coils remains static throughout the simulation. The possibility of including a fixed loop voltage has not been used because no significant effect was observed in preliminary studies. Ideally, one would want to reproduce the feedback control system; however, the goal in this series of simulations is to highlight the dominant effects, the timings and sensitivity to various parameter changes without the complexity of a time-dependent drive.

The MHD fields are discretised in the toroidal direction with Hermite cubic splines on 24 equidistant planes. The unstructured triangular mesh, which is the same on each poloidal plane, provides support for the C¹ finite element solution in R and Z coordinates. Mesh points can be packed where fine structures and sharp gradients are expected to form. In the case of downward VDEs, the triangulation is made denser at the bottom of the vessel yielding an anisotropic mesh, as depicted on Fig. 1.

The mesh is decomposed in three regions:²¹ (i) the inside of the vacuum vessel (in red on Fig. 1) where the extended MHD equations below are solved, (ii) the finite-thickness resistive wall (2 cm in all cases here, in blue in Fig. 1) where only the magnetic field is evolved according to $∂tB+∇×(ηwallj)=0$ and $∇×B=μ0j$⁠, and (iii) the exterior vacuum region (in green in Fig. 1) where $j=0$ is enforced. By construction of the mesh and because the wall resistivity ηwall is constant, the wall is axisymmetric. The capability of modelling ports, breaks, or passive conductors with a spatially dependent wall resistivity exists and will be used in future simulations.

The VDE simulations are evolved according to single-fluid extended MHD equations. $M3D-C1$ allows for two-fluid effects and more sophisticated closures, which might be important to model the plasma–wall interaction and scrape-off layer physics. Here, the equations solved comprise the continuity equation for the fluid density n

the momentum equation for the fluid mean flow $v$

and the energy equation for the fluid isotropic pressure p

coupled to Faraday, Ampère, and Ohm's laws for the magnetic field $B$

In the equations above, m is the ion mass, the stress tensor is closed by $Π=μ(∇v+∇vT)+λ(∇·v)I$ and the heat flux by $q=−κ⊥∇T−κ||bb·∇T$ where $b=B/B$ and the temperature is $T=p/n$⁠. The adiabatic constant is $Γ=5/3$⁠.

The boundary of the computational domain (in magenta in Fig. 1) acts as a perfect conductor for the magnetic field, where the normal component is held fixed to its initial value. Section III B demonstrates that the boundary is sufficiently far away not to affect the evolution of the VDE. Dirichlet boundary conditions are applied on the fluid variables at the resistive wall: density is set to $nedge$⁠, pressure to $pedge$ and the flow is assumed to vanish (no slip).

Inside the vacuum vessel (red region in Fig. 1), the resistivity is a spatially varying function through the temperature according to the modified Spitzer expression

The offset $Toffset$ effectively lowers the temperature just for computing the plasma resistivity. While its effect is negligible in the core, this bias provides control over the resistivity in the so-called halo region beyond the last-closed flux-surface. The halo region is initialised as a uniform cold and low-density plasma through the choice of edge pressure, $pedge∼8$ Pa, and edge density, $nedge=0.01×1020 m−3$⁠, so that numerical instabilities and negative pressure overshoots are avoided during the vertical drifting phase (advection). The halo temperature is thus set to $Tedge=pedge/nedge∼25 eV$ and the halo resistivity is approximately $ηhalo=1.4×10−5 Ω m$⁠. Dividing by the surface area of the halo region, the resistance of the halo region competes with that of the wall, such that the halo plasma has a damping effect on instabilities, as seen in Secs. III A and IV. In reality, the halo is thin, sparse, and about three times colder, which is the reason for using $Toffset$⁠.

The density diffusion D, heat conductivity (parallel $κ||$⁠, perpendicular $κ⊥$⁠), viscosity coefficients $μ,λ$⁠, and the ratio to Spitzer resistivity η₀ are input parameters. Diffusion prevents sharp structures and negative overshoots from accumulating. For numerical stability, these transport coefficients are usually higher than realistic values. It is argued, however, that, as long as the correct timings are respected in the simulations, the VDE evolution is only weakly affected by the adjustment of transport coefficients.

The equations listed above are solved using implicit time-stepping,²⁶ a feature that is necessary to perform stable simulations over long timescales and accurately resolve the advection-driven dynamics as well as the wide separation (and gradients) between Alfvén and resistive dynamics. Indeed, characteristic Lundquist numbers attained in the simulations are $Sp∼106−107$ in the plasma core, $Shalo∼100−1000$ in the halo region and $Swall∼104$ in the resistive wall.

Experimental traces indicate that the plasma remains stable (axisymmetric) during the vertical drifting phase and preserves its energy content, flux-surfaces, and profiles. This observation is exploited to save computational resources and deploy the heavier 3D simulations only during the phase when the plasma contacts the wall. Our simulations are thus first carried out in 2D up until the linear stability analysis, performed in parallel at frequent intervals, predicts the fast growth of $n≠0$ modes. This roughly coincides with when the X-point reaches the wall and the plasma becomes limited. The nonlinear simulations are resumed in 3D, a few time steps before the time of contact, and followed until the plasma current has fully quenched. While 2D nonlinear VDE simulations run to completion within a week on a few dozen processors, the 3D nonlinear segments, that are launched on thousands of processors, require months of computation.

The fine-tuning of input parameters is necessary to ensure that the appropriate physics regime is reached and that experimental targets are deliberately met. Sensitivity scans are efficiently performed with 2D nonlinear simulations. Figure 2 shows a series of axisymmetric runs where input parameters are scanned to alter the timing of events.

From an MHD point of view, the plasma must be in excellent force balance with Eddy currents to evolve on timescales that are about three orders of magnitude slower than Alfvénic (This statement extended to 3D explains the slow evolution of external kink modes in devices like JET.^27,28). The plasma column thus moves as a sequence of equilibria, drifting due to the relaxation of both wall currents²⁹ and plasma current.³⁰ The latter mechanism is subdominant in hot VDEs, but may become an important aspect of VDE mitigation strategies in e.g., ITER for which the wall time is likely to be longer than the plasma current decay time.

Comparing the orange, green, and red curves on the top plot of Fig. 2, the VDE is generally slower when the wall resistivity is decreased. As a counter-example, the vertical motion is faster in red cases than blue, even though the wall resistivity is three times lower in the former. The plasma resistivity is multiplied by a factor 30 in all cases except the blue, for which the exact Spitzer value is used. A faster plasma current decay contributes to accelerate the vertical motion. The blue cases are more consistent with the fact that the drifting phase lasts for approximately 50 ms while the plasma current remains almost flat. The red cases better reproduce the fact that the toroidal wall current reaches more than half of the initial plasma current during the fast current quench phase.

A higher halo temperature is seen to stabilise the VDE, as seen by comparing the solid curves with the dashed curves, where $Thalo=25$ eV (⁠$Toffset=0$⁠) and $Thalo=9.9$ eV (⁠$Toffset=15$ eV), respectively. The effect is stronger in the blue series, where the plasma current remains closer to its initial value before touching the wall. A more thorough investigation of the effect of the halo temperature is presented in Sec. III A.

The characteristic VDE time-scale is, to lesser extent, sensitive to transport coefficients, as seen by comparing the red dashed and dash-dotted curves whose particle diffusion and heat conductivity each differ by a factor 10 (combined difference of a factor 100). The gradual loss of thermal energy in all cases is an unavoidable consequence of the relatively high perpendicular heat conductivity required for numerical stability (possibly aggravated by the mesh resolution). The decrease in temperature, which causes the plasma resistivity to increase, may lead to an acceleration of the vertical motion.

Viscosity has a weak stabilising effect on the VDE, as noticed when comparing the dotted and dashed blue curves; this coefficient is about 100 times smaller in the former case, yielding a slightly faster VDE. Viscosity helps evacuate strongly sheared flows and is often increased to overcome the accumulation of strong gradients. The influence of viscosity on the onset of instabilities (in particular 3D) remains to be determined.

As qualitatively shown on Fig. 2, the characteristic timescale of the vertical drifting phase depends on the halo temperature. Indeed, the open field-line region is filled with a low temperature plasma. The latter is in direct thermal contact with the wall through high parallel heat conductivity $κ||$ and tends to remain cold and resistive. However, the cross-sectional area of the halo region is quite large such that its resistance is comparable to that of the wall. The decay rate of currents in the halo region is then of the same order as the growth rate of some instabilities. In this case, the halo is stabilising since the currents it carries have time to oppose to flux changes (positive or negative depending on the drive). In reality, the impact of the halo region on the vertical instability is minor since the scrape-off-layer is thin, sparse, and cold.

The instantaneous growth rate, $γVDE(Z)=Ż/(Z−Zref)$⁠, gives an estimate of the VDE timescale. Computed at a fixed vertical position, the scaling of the instantaneous growth rate is expected to be linear as a function of the wall resistivity,²⁹ at least asymptotically when the wall tends to a perfect conductor. Figure 3 shows the dependency on the wall resistivity of the instantaneous VDE growth rate at $Z=0.1 (m)$ for three different halo temperatures. The blue curve on Fig. 3, which corresponds to $Toffset=0$⁠, deviates significantly from linear behaviour (depicted by the dashed black curve on Fig. 3), because the halo region is able to oppose the downward motion of the plasma. The effect is strongest when the VDE is fastest and the wall resistivity exceeds the effective halo resistivity (see vertical dotted lines on Fig. 3). By adjusting the offset temperature [Eq. (6)] to decrease the effective halo temperature, linear scaling is partially recovered at low wall resistivity, as seen on the red and orange curves on Fig. 3. The behaviour of the VDE growth rate in 2D nonlinear simulations is consistent with the previous linear stability analysis.²¹ 

The edge of the computational domain acts as a perfect conductor located at a given distance beyond the external PF coils, enforcing Dirichlet boundary conditions on the magnetic field in the outer vacuum region. The computational boundary may have a stabilising effect on the vertically displacing plasma if positioned too close, due to artificial mirror currents that freeze the normal component of the magnetic field. Figure 4 compares two scans of the instantaneous VDE growth rate varying the wall resistivity at $Z=0.1(m)$ for different computational boundaries. The red curve on Fig. 4 relates to the mesh shown on Fig. 1. The black curve is for a mesh whose boundary is expanded outwards by 150%. The comparison demonstrates that the location of the computational boundary has almost no effect on the evolution of the VDE in the slower (more realistic) cases.

The resistive wall, which is 2 cm thick in our simulations to be consistent with drawings of the NSTX vacuum vessel (Ref. 22, Fig. 1), may develop skin currents. It is easily verified that the mesh is sufficiently dense to resolve such sharp structures. The skin depth due to the vertical instability is estimated to be $δ=2ηwall/γVDEμ0≳10$ cm, while the typical dimension of a triangle within the unstructured mesh is $∼1$ cm. Figure 5 confirms that the VDE growth rate is identical when the mesh is denser by a factor 4.

In this section, the evolution of two 3D simulations, whose settings and histories slightly differ, is discussed. The first case is launched with higher values of halo temperature $Thalo=25$ eV (⁠$Toffset=0$⁠) and wall resistivity $ηwall=4.9×10−5 Ω m$ corresponding to the green curve on Fig. 2. The second case is run with $Thalo=9.9$ eV and $ηwall=1.9×10−6$⁠, identical to the dashed-red curve on Fig. 2.

The moment the column comes in contact with the wall coincides with a change of physics regime. During this phase, the plasma current decays abruptly and the plasma is prone to non-axisymmetric deformation. The onset of 3D instabilities occurs because (i) the plasma violently becoming kink unstable as flux-surfaces are peeled off and the edge q drops below a certain threshold (⁠$qedge≲2$⁠), (ii) the compression of flux near the wall leading to current density clumps and edge modes (peeling-ballooning type). In both scenarios, flux-surfaces are destroyed and thermal energy is efficiently released through parallel heat transport along stochastic field-lines, as exemplified on Fig. 6 showing a sequence of Poincaré sections for the green 3D case. In both simulations, toroidal modes are triggered near the edge of the plasma, breaking the flux-surfaces from the outside-in. The resulting downfall of temperature makes the plasma resistivity soar, thereby precipitating the current quench. In response to the rapid decay of toroidal current in the core, current is driven in the open field-line region (halo) connecting to the wall.

Figure 7 provides a detailed comparison between the 2D and 3D evolution for the two cases. The time $tmaxIw$⁠, when the toroidal wall current reaches its maximum value in the 2D simulations (Fig. 2), is used as a reference to normalise the x-axis and overlay the green and red cases of Fig. 2 in a meaningful way.

In the green case, the plasma remains axisymmetric until it thermalises with the wall temperature. Toroidal modes appear only at the end of the slow thermal quench, when the plasma has considerably shrunk. The halo region, whose temperature is relatively high, acts as a stabilising shell by broadening the contact area with the wall. The effect of non-axisymmetric toroidal modes on the course of events is weak since the plasma and wall currents return to their axisymmetric value after the brief appearance of 3D modes. Harmonics rise and decay in proportion with a dominant n = 2 component.

In the red case, the onset of toroidal modes comes at an early stage. Their presence is far more disruptive, as seen by the fact that both the plasma thermal energy and plasma current are abruptly released. Accordingly, the toroidal wall current shoots up in order to preserve total flux. The time-window for toroidal modes is wider than in the green case. Toroidal harmonics arrive in waves, cascading from high to low mode numbers. The n = 4 mode is the first to reach its peak value, followed by a strong n = 2 mode which dominates the spectrum until an n = 1 ultimately takes over.

Figure 8 presents a time sequence of the current density on the divertor plate of the device, where the toroidal component is depicted in black arrows and the normal in red/blue colours. The toroidal angle has been straightened out in the y-axis while the x-axis represents the radial length along the width of the divertor plate. Both 3D simulations exhibit rich patterns and similar features. Only the red case is reported on Fig. 8. The contact line, identified as the separation between red and blue colour, stays roughly at the same radial position. At first, the toroidal current points in the opposite direction from the plasma current, in response to the vertical plasma motion. As the toroidal current transitions towards the co-current direction, high-n patterns are visible in the normal component of the current density near the contact line. These footprints cascade into an n = 2 perturbation, present also in the toroidal component. The non-axisymmetric patches ultimately decay and the current density on the divertor plate flows in the co-current direction as the plasma vanishes into the wall. In comparison, shunt tiles, located at discrete poloidal and toroidal locations at the bottom of NSTX vessel, routinely measure the induced currents in the open-field line region connecting with the wall. The toroidal resolution of these shunt tiles is enough to resolve the n = 1 component only. While the amplitude of normal currents qualitatively match between the experimental traces and the simulation diagnostics, the mode number, duration, timing, and rotation of halo current substantially differ. Almost no global (n = 1) rotation of the halo currents is observed in our simulations. At a fixed poloidal location, the patterns rotate only due to a shearing effect which arises from scraping-off of outer flux-surfaces and reaching lower values of the q-profile. On the other side of the contact line at equal distance, the pattern rotates in the opposite direction. The lack of rotation in the simulations is a consequence of solving single fluid MHD model and the presence of a strictly axisymmetric wall. Including more physics (two-fluid, sheath, 3D wall geometry) at the plasma-wall boundary would provide the drive (or locking mechanism) for n = 1 rotation. Only then can our simulations be compared with experimental halo rotation databases. The numerical implementation of the realistic plasma–wall interaction is, however, far from trivial, especially in fully 3D nonlinear simulations.

The modelling of NSTX hot VDE shot #139536 was performed using the $M3D-C1$ code, recently upgraded to include a finite-thickness resistive wall within the computational domain. The sequence of events observed in the experiment was matched by scanning the input parameters until realistic conditions were established. The drifting phase, during which the plasma is immune to non-axisymmetric modes, was approached using nonlinear 2D simulations up until the X-point reached the wall. Beyond this point, simulations were resumed to include 3D effects for the remainder of the VDE, a task about an order of magnitude more computationally intensive.

The series of 2D simulations revealed an acute sensitivity of the timing of events and the drive of the vertical instability to changes in input parameters, especially the choice of wall resistivity and halo temperature. The instantaneous growth rate was shown to deviate from an expected linear relation with wall resistivity. The halo region was found to stabilise the vertical motion of the plasma column. The use of a small offset to lower the temperature in the Spitzer expression conveniently increased the halo resistivity to realistic values such that linear scaling of the VDE growth rate was recovered at low wall resistivity.

The 3D evolution of two cases was analysed from the point of contact until the vanishing of the plasma current. The halo region was shown to be responsible for a late onset of toroidal modes in the case where the halo temperature was unbiased. In the case where the effective halo temperature was brought to realistic values by an appropriate temperature offset, instabilities lead to an abrupt quench of thermal and magnetic energy. Poincaré plots revealed that the flux-surfaces were progressively destroyed by inwards-penetrating edge modes. The stochastisation of field-lines led to parallel heat transport, thereby rapidly cooling the plasma. The fast decay of plasma current that follows from the increased resistivity induces currents in the open-field line region. The patterns created by the normal component of the current density to the divertor plate was seen to cascade from high-n to low-n mode numbers.

With an experimental estimate of the width of the halo region, the correct halo temperature could be inferred for producing more accurate simulations. Determining the halo width requires an increase in the poloidal resolution of shunt tiles in present-day machines. The topic of halo current rotation was vastly unaddressed in this series of simulations. Limiting factors include boundary conditions, the physics of plasma–wall interaction and the axisymmetry of the resistive wall. A proper benchmark with experiments of the complex rotation and shearing patterns of halo currents, however, requires an increase in the toroidal resolution of shunt tiles in today's devices.

See supplementary material for movies of the VDE evolution. Videos include for the red and green case discussed in Sec. IV (i) a poloidal cross-section of the time-varying toroidal current density, temperature, and particle density, (ii) the time-varying patterns from halo currents on the divertor plate.

The authors would like to acknowledge stimulating discussions with C. E. Myers, A. Boozer, H. Strauss, C. Sovinec, and J. Bialek. The authors were supported by the Department of Energy Contract No. DE-AC02-09CH11466 during the course of this work.