3D simulations of vertical displacement events in tokamaks: A benchmark of M3D-C¹, NIMROD, and JOREK

Vertical displacement events (VDEs) are axisymmetric instabilities that arise for elongated plasmas when the control of the vertical position is lost. In future large fusion devices such as ITER, major disruptions will unavoidably lead to the loss of plasma control¹ and a resulting VDE. Such events drive the plasma column into the wall, causing large thermal and electromagnetic loads to the plasma facing components and to the vacuum vessel.² During VDEs, additional 3D MHD instabilities can arise resulting in the localization of the loads and in net sideway forces. Moreover, it has been observed that these forces can rotate toroidally, leading to large forces on the vacuum vessel and its supporting structures when the rotation frequency resonates with the natural frequencies of the vessel.^3,4

In this respect, the development of validated 3D MHD codes including the effect of resistive walls is crucial to the study of these loads for future machines. Although axisymmetric codes such as DINA,⁵ TSC,⁶ or CarMa0NL⁷ allow studying many crucial aspects of disruptions, they cannot take into account these 3D effects. Moreover, for VDEs in which the vertical motion precedes the thermal energy loss (known as hot VDEs), the edge safety factor (q_a) is strongly reduced due to the slow decay of the plasma current during the vertical motion.⁸ In such cases, axisymmetric codes may impose $qa≥1$ by modifying the current profile, otherwise q_a drops below unity when the minor radius approximately drops by a factor of 2 (⁠$qa∼a2/Ip$⁠). When q_a falls below 2 or so, the configuration will no longer remain axisymmetric.^2,9 Non-axisymmetric modes will grow, producing horizontal forces and stochastic field lines. The minimum q_a and the level of current density are strongly influenced by MHD activity and parallel dynamics in the stochastic field line region. The dynamics affect the current density directly and by modifying the electrical resistivity through the fast temperature reduction.

In recent years, simulations of 3D VDEs have been performed by Strauss¹⁰ with the M3D code,¹¹ Sovinec and Bunkers¹² with NIMROD,¹³ Pfefferlé et al.¹⁴ with M3D-C¹ (Ref. 15), and Artola¹⁶ with JOREK.^17,18 However, no benchmark for 3D VDEs has been performed among these codes. Except for highly idealized cases, there are no analytical solutions for events that involve chaotic magnetic field lines and large vertical displacements. The work presented in this article is therefore an essential contribution to the verification and validation of these numerical codes.

In this work, we present such a benchmark for the codes JOREK, M3D-C¹, and NIMROD. These three codes are among the limited number of codes capable of simulating 3D MHD instabilities in tokamak geometry including resistive walls. The benchmark is addressing the full 3D dynamics and is partly based on the axisymmetric benchmark that was already performed between these three codes.¹⁹ The goal of this work is to demonstrate the consistency among these codes and to provide the fusion community with a useful benchmark for MHD simulations of disruptions.

This paper is organized as follows. In Sec. II, we revisit the models and the numerics used by the different codes. The setup and the parameters of the benchmark case are described in Sec. III. The comparison of the obtained results and the analysis of the 3D VDE simulations are presented in Sec. IV. Finally, we summarize our conclusions in Sec. V.

The single fluid MHD model used by the three codes for this benchmark is given by the following equations:

where E is the electric field, v is the mean flow velocity, B is the magnetic field, J is the current density vector, ρ is the ion density, p is the total pressure, and $T≡Te+Ti$ is the total temperature. The other parameters are the plasma resistivity (η), the stress tensor (⁠$τ¯$⁠), the particle diffusion tensor (⁠$D¯$⁠), the heat conductivity tensor (⁠$κ¯$⁠), the vacuum permeability (μ₀), and the ratio of specific heats (γ). The particle diffusion and heat conduction tensors are decomposed into parallel and perpendicular directions relative to the magnetic field (⁠$b=B/|B|$⁠),

where $I¯$ is the identity tensor. For the simulations presented in this paper, the Ohmic heating source, radiation losses and external particle, momentum and energy sources are not included to keep the setup simple. Except for the plasma resistivity, for which a Spitzer-like temperature dependence is used (⁠$η∝Te−3/2$⁠), the other coefficients (⁠$κ∥,κ⊥,D∥,D⊥$⁠) do not depend on any variable and are spatially constant. The exact values for these coefficients used in our computations can be found in Table II. The form of the stress tensor for each code is given in the code-specific Secs. II B–II D.

The MHD equations (1)–(6) are solved in a volume delimited by a resistive wall [see the blue line in Fig. 1(c)]. Inside the resistive wall, Ohm's law is simply

where η_w is the wall resistivity, which for this benchmark is spatially constant, and $Jw$ denotes the wall current density.

Dirichlet boundary conditions are applied for the fluid variables at the plasma–wall interface. The temperature and density B.C.s are simply $ρ=ρ(t=0)$ and $T=T(t=0)$⁠. These conditions do not capture the density and temperature variations along the plasma–wall interface that may take place in experiments. The extension of boundary conditions to more comprehensive modeling is the subject of ongoing work.^20,21 Although the three codes set a no-normal flow boundary condition at the plasma–wall interface (⁠$v·n=0,$ where n is the normal vector to the boundary), NIMROD and M3D-C¹ set all the velocity components to zero (⁠$v=0$⁠) whereas in JOREK, only the components normal to the wall and parallel to the magnetic field are set to zero (⁠$v·n=v∥=0$⁠). For all three codes, the modeling is analogous to the cold-wall fixed thermal-conduction computations of Sec. III C in Ref. 20. There, it is shown that rapid open-field heat transfer to cold walls keeps the halo region relatively thin, and the evolution is insensitive to boundary conditions on v. With the combination of $v·n=0, D≠0$⁠, and the fixed-ρ boundary condition, loss of mass occurs through particle diffusion at the plasma–wall interface and is decoupled from the magnetic field due to the temperature-dependent resistivity. The treatment of boundary conditions for the magnetic field is code dependent and is explained below.

The primary numerical properties of the codes are summarized in Table I. Sections II B–II D describe the particularities of the models used in JOREK, M3D-C¹, and NIMROD.

JOREK is a fully implicit non-linear extended MHD code for realistic tokamak geometries including open field-line regions.¹⁷ Although a full MHD model is available in JOREK,^22,23 the extension for resistive walls has not yet been implemented for that physics model. For this benchmark, we use a reduced MHD model that is based on the following ansatz for the magnetic field and the mean plasma velocity:

where ψ is the poloidal magnetic flux, $Φ$ is the electrostatic potential, and $F0=RBϕ$ is a constant representing the main reduced MHD assumption, which is that the toroidal field is fixed in time. However, poloidal currents are not fixed in time and they evolve according to the current conservation and momentum balance equations; only their contribution to the toroidal field is neglected. Note that the reduction of the equations is ansatz-based and does not result from geometrical ordering assumptions. The specific projection of Eqs. (1)–(6) together with the ansatz is explained in detail in Ref. 24. The form of the divergence of the stress tensor in the poloidal plane projection is $∇·τ¯|pol=R2μ∇ϕ×∇wϕ≈μ∇2vpol$⁠, where μ is the dynamic viscosity, $wϕ≡∇pol2Φ$ is the toroidal vorticity, and the subscript “$pol$” denotes the directions lying on the poloidal plane. For the parallel direction, the divergence of the stress tensor is chosen as $∇·τ¯|∥=μ∇2v∥$⁠.

The resistive wall and the boundary conditions for the magnetic field are included by coupling JOREK to the STARWALL code.^16,18,25 The coupling is performed by solving the full vacuum domain with a Green's function method, and therefore, the vacuum does not need to be meshed. STARWALL uses the thin wall assumption and discretizes the wall with linear triangular elements. Although the employed wall formalism is implemented for 3D thin walls including holes, the used setup is restricted to an axisymmetric wall. The no-normal flow boundary condition that is common to the three codes implies that the boundary conditions in JOREK must be $v∥=0$ and Φ = 0. Since in reduced MHD the poloidal electric field only depends on the electric potential (⁠$Epol=∇polΦ$⁠), the no-normal flow condition (Φ = 0) implies that the plasma–wall interface acts as a perfect conductor in the poloidal direction along the boundary (⁠$Epol×n=0$⁠). Note that this is not the case in M3D-C¹ and NIMROD, where at the boundary $Epol×n=ηJpol×n≠0$⁠. Finally, since the resulting momentum equation contains third-order derivatives on $Φ$⁠, an additional boundary condition is applied for the toroidal vorticity (⁠$wϕ=0$⁠).

JOREK discretizes the poloidal plane with quadrilateral bicubic Bezier elements using an isoparametric mapping.²⁶ The variation along the toroidal direction is represented by Fourier harmonics and for the time advance typically the Crank–Nicolson or the BDF5 Gears scheme is used.

M3D-C¹ is a versatile MHD code with a number of options. It can be run in the 2D (axisymmetric) mode, a linear-3D mode (for a single toroidal harmonic), or in a fully 3D nonlinear mode as was used for this paper. The magnetic vector potential and ion velocity are represented in terms of scalar quantities as $A=ψ∇ϕ+R2∇ϕ×∇f−F0 ln R∇Z$ and $v=R2∇U×∇ϕ+wR2∇ϕ+R−2∇polχ$⁠. Here, $(R,ϕ,Z)$ are cylindrical coordinates. The function f is constrained to be regular at the origin so that the constant F₀ is proportional to the total current in the toroidal field magnets. The code can be run in the two-variable reduced MHD option where only (⁠$ψ,U$⁠) are advanced in time, in the four-variable reduced MHD option where (⁠$ψ,f,U,w$⁠) are advanced in time, in the full single fluid MHD option (used in this paper), or in a two-fluid MHD option. Implicit hyper-resistivity terms can be included in the time advance. In this paper, hyper-resistivities with spatially constant values of $2.76×10−9$ and $1.10×10−9 Ω m3$ were used in the equations for ψ and f, respectively. The stress tensor in M3D-C¹ is $τ¯=μ[∇v+∇vT]$⁠. Although M3D-C¹ has the option to advance the pressure, Eq. (6), in these calculations we advanced the temperature instead, T, the same as described by Eq. (13) in Sec. II D.

In the resistive wall region, only the magnetic vector potential variables are advanced in time, with Ohm's law: $E=ηwJ$⁠, where the wall resistivity may be a function of space even if it is set to a constant in this benchmark. In the vacuum region: $∇×∇×A=0$ is enforced. This region is directly meshed and solved by using perfect conductor boundary conditions [i.e., $B·n=B·n(t=0)$] at the outermost computational boundary containing the wall and the vacuum regions. The mesh and outermost boundary are given in the axisymmetric version of this benchmark (see purple curve of Fig. 1 in Ref. 19). The location of the outermost boundary is such that the VDE growth rates are close to the no-outer-boundary growth rates.¹⁹ 

The code uses an unstructured mesh with triangular elements in the (R, Z) plane and extruded in a structured manner in the toroidal direction. The mesh is normally adapted to the geometry and the problem, so that a finer mesh size is used in locations where large gradients are expected. All scalar variables are expanded in 3D finite elements that are a tensor product of the Bell²⁷ element in the (R, Z) plane and Hermite cubic elements in the toroidal angle $ϕ$⁠. This representation enforces continuity of each scalar and all of its first-derivatives across element boundaries. In the 3D nonlinear code, a split-implicit time advance is used where the velocity variables are advanced first, followed by the magnetic-field variables, the pressure(s), and then density.

The NIMROD code solves fluid-based models for magnetized plasma in non-reduced form. Unlike JOREK and M3D-C¹, it directly evolves the components of magnetic field B and flow-velocity v. A modified form of Faraday's law is used to advance B

where the second term on the right diffusively corrects any numerical divergence error. In addition, the pressure evolution equation (6) is replaced by an equivalent temperature equation using the ideal gas law p = nk_BT and Eq. (5)

Although many NIMROD computations use fixed thermal diffusivity coefficients $χ∼κ/n$⁠, the computations reported here and in Ref. 19 use fixed thermal conductivity coefficients $κ∥, κ⊥$ for the benchmarks with JOREK and M3D-C¹. NIMROD's stress tensor is formed with the traceless rate-of-strain tensor, $τ¯=μ[∇v+∇vT−(2/3)I¯∇·v]$⁠.

The VDE modeling uses the thin resistive wall approximation that is described in Ref. 12 to interface with a meshed external numerical vacuum response, where the magnetic representation is the same as in the internal plasma region. The meshed vacuum-field computation is similar to the approach used in M3D-C¹, and the shape of the outer region is chosen to approximate that used in the M3D-C¹ computations.¹⁹ However, NIMROD's thin-wall model interfaces the internal and external regions without resolving the wall's cross section. The equations conserve toroidal magnetic flux over the sum of the two regions, and the normal component of the equilibrium poloidal-B at the outer surface of external vacuum region is fixed for all time. As noted previously, the cold-wall fixed thermal-conductivity conditions makes extended-MHD computations insensitive to boundary conditions on v, and we set all components of $v=0$ at the resistive wall.

The system is advanced in time with an implicit leapfrog method that, in comparison with fully consistent implicit methods, requires smaller but less computationally intensive steps to achieve the same level of accuracy.²⁸ Like JOREK, NIMROD expands $ϕ$-dependent fields in toroidal Fourier harmonics. For the poloidal plane, the expansion is 2D C⁰ spectral elements with node spacing based on Legendre polynomials. The underlying polynomial representation is chosen at runtime, and except where noted below, all computations reported here use polynomials of degree 3. The representation does not satisfy the magnetic divergence constraint, identically, but the diffusive error correction term in Eq. (12) is effective at maintaining minimal error with NIMROD's high-order spatial representation.¹³ The NIMROD computations also use the stabilization method that is described in Ref. 29 so that the numerical representation reproduces MHD interchange with growth rates increasing to the converged value as resolution is increased.

As mentioned in Sec. I, the benchmark is a 3D version of the axisymmetric benchmark that was already published for these three codes.¹⁹ The chosen equilibrium is based on an NSTX experimental plasma (discharge No. 139536 at t = 309 ms) that was reconstructed with the EFIT code.³⁰ The equilibrium profiles, the separatrix, and the wall geometries and several scalar quantities describing the equilibrium are shown in Fig. 1. The temperature profile is calculated from the pressure profile (p) given in the geqdsk file with the expression $Te=Ti=946 eV(p/paxis)0.6$⁠. Note that this file is available in the supplementary material of Ref. 19 as well as the currents and geometry of the poloidal field coils that are needed to compute the free-boundary equilibrium. At the plasma–wall interface, the applied Dirichlet boundary conditions for the temperature and density are $Te,edge=14.6$ eV and $ne,edge=2.20×1018 m−3$⁠.

For simplicity, a rectangular vacuum vessel is chosen instead of the complicated geometry of the NSTX vacuum vessel [see Fig. 1(c)]. The rectangular shape is defined by the four vertices $(R=0.24 m,Z=±1.4 m)$ and $(R=1.6 m, Z=±1.4 m)$⁠. Note that the wall corners are not exactly sharp in the simulations because Fourier harmonics are used to represent the wall contour. The thickness of the resistive wall is $Δw=0.015$ m and the wall resistivity is $ηw=3×10−5Ω$ m. Where a thin-wall approximation is used, the effective thin wall resistivity $ηw/Δw$ is used.

The simulation is divided into two phases: an axisymmetric (2D) run and a 3D run. The values for different parameters common to both phases and specific to each phase are specified in Table II. In the axisymmetric phase, the X-point drifts downwards until the plasma becomes limited by the lower part of the wall. The end of this phase and the start of the 3D phase are therefore determined by the change of geometry of the last closed flux surface (LCFS) (i.e., it becomes defined by a limiter point instead of the lower X-point). The 2D phase is identical to the initial phase for the non-linear benchmark presented in Ref. 19 but with the wall resistivity, plasma resistivity and particle diffusion and heat conduction coefficients increased by a factor of 10. Note that during the full simulation, the plasma resistivity is a factor 10 larger than the Spitzer's value. Such an increase was required to reduce the computational cost of the simulations. Although the main goal of this work is not an experimental comparison, the chosen parameters lead to vertical displacements of 1 m in a timescale of 10 ms, which is similar to the VDEs observed in experiments³¹ where such displacements take place within 5–20 ms. In the 2D phase, the heat conduction and particle diffusion coefficients are such that the thermal energy does not decay during the vertical motion. Since the parallel density transport is typically governed by parallel convection, for simplicity the parallel diffusion coefficient is set to be equal to the perpendicular particle diffusion coefficient. In the 3D phase, we significantly increase the diffusion/conduction coefficients to smooth the sharp pressure gradients that arise due to the fast shrinking of the LCFS. Although for that phase the diffusion/conduction coefficients do affect the thermal energy decay, it will be shown in Sec. IV A that the final collapse of plasma energy is governed by the 3D MHD activity once the magnetic field topology becomes chaotic.

In JOREK, a polar grid is used with increased resolution in the region where the plasma becomes limited by the wall. The number of Bézier elements used in the plasma region is 22 000, and the number of linear triangular elements to mesh the thin wall is 48 000. For the toroidal direction, 11 Fourier modes were used with $n∈[0,10]$⁠. A resolution scan in the number of Fourier modes is shown in Sec. IV B. Time steps of the order of 5–10 Alfvén times were used during the 2D phase and time steps of 0.1–0.2 Alfvén times were used for the 3D phase.

The M3D-C¹ calculation used the unstructured poloidal plane mesh shown in Ref. 19, which has 17 424 vertices on each plane, each with 12 degrees of freedom for each scalar field in 3D (6 in 2D). The 3D phase used 16 toroidal Hermite cubic finite elements for each of the triangular vertices. A time step of $Δt=0.5τA$ was used, except for a period of 0.048 ms, starting at time 10.185 ms (1.233 ms after the start of the 3D) when it was halved to $Δt=0.25τA$⁠, to avoid numerical instability. After this time, an “upwind” second-order toroidal diffusion term was added to the scalar convection terms for the pressure, density, and magnetic field, equal to $0.125×δx×vϕ$⁠, where $δx=0.4$ m was the approximate distance between toroidal planes.

The NIMROD computations have been run with 39 000 quadrilateral bicubic elements for the region inside the resistive wall. The meshing is largely rectangular, but a packed layer of smaller elements is used near the resistive wall and conforms to the curved corners. On average, the bicubic C⁰ elements have nine distinct nodes when taking into account the sharing of element-border nodes. Thus, the inner region has 1 complex degree of freedom per Fourier harmonic for n, T, and each component of v and B on 351 000 nodes. Like the JOREK computation, toroidal variations are represented with Fourier harmonics $n∈[0,10]$⁠. Regarding the time-advance, although NIMROD's implicit leapfrog uses implicit advection, the temporal staggering tends to lose accuracy for dynamics if the time step does not satisfy the Courant–Friedrichs–Lewy (CFL) condition for flow-velocity. Thus, the time step is adjusted, dynamically, to satisfy this condition. For the most dynamic part of the 3D computations $Δt>0.015τA$⁠, which is hundreds of times larger than an explicit stability limit. In other parts of the 3D evolution, time steps are more than an order of magnitude larger.

In this section, the simulated case is presented, the results are compared between the three codes, and the physics of the case is analyzed. In Sec. IV A, the results for the different phases of the simulation are shown together with a discussion on MHD stability. In Sec. IV B, a resolution scan is performed to demonstrate that the case is sufficiently converged in toroidal resolution. In Sec. IV C, the origin of the wall forces is explored as well as the influence of the parallel heat conductivity, the viscosity and the initial asymmetric perturbations on the results.

We compare the vertical position of the magnetic axis, the total toroidal current (including the halo current) and the total thermal energy for the three codes in Fig. 2. The total toroidal halo current is also shown for JOREK and M3D-C¹. NIMROD's post-processing tools to compute the toroidal halo current have not been developed yet. As expected for the 2D phase (⁠$t′≡t−t3D,start<0$⁠), the level of agreement is similar to the axisymmetric VDE benchmark reported in Ref. 19. The time axis has been shifted so $t′=0$ when the plasma becomes limited in the three codes. As was shown in Ref. 19, such synchronization is necessary because the different numerical perturbations initiating the VDE can lead to different timescales even if the growth rates are in good agreement.

In this benchmark case, the complete thermal quench is induced by the vertical displacement (see Sec. IV A 2) and takes place when the plasma volume has decreased by more than a factor of two (⁠$Vp/Vp0=0.26$ at $t′=1.16$ ms in JOREK). Due to this sequence of events, the case can be classified as a hot VDE, which features the largest wall forces observed in experiments. During the 2D phase, the plasma current and the thermal energy do not decay due to the high plasma temperature and the small diffusion/conduction coefficients. The chosen increase in the diffusion/conduction coefficients at the start of the 3D phase leads to a loss of a large fraction of the thermal energy. Note, however, that once 3D MHD instabilities set in, the decay rate of the thermal energy is dominated by the 3D effects and not by the choice of diffusion/conduction coefficients (see β_N time trace in Fig. 3). Such an increase in the diffusion/conduction coefficients was necessary to prevent the formation of large edge pressure gradients that arise when the plasma volume decreases. Simulations without the increase in the diffusion/conduction coefficients showed that high-n edge localized modes become unstable due to the large gradients. To account for the high toroidal mode numbers adequately, the numerical resolution would have to have been significantly larger and additional terms (e.g., diamagnetic flows terms) would have to have been included in the equations, greatly increasing the complexity of the benchmark. As this work is a benchmark exercise which primarily aims to compare the 3D wall forces (which are observed to be given by low-n mode numbers), the high-n edge localized modes are avoided with the choice of the diffusion/conduction coefficients.

Linear stability calculations were performed at different vertical positions during the 2D phase. All three codes find that the initial equilibrium is unstable to resistive edge instabilities localized at the q = 3 and q = 4 rational surfaces, but as the equilibrium profiles evolve due to current diffusion and the vertical motion, these modes are stabilized. For example, when performing the stability calculations at $Zaxis=−0.15$ m, corresponding to $t′=−1.4$ ms in Fig. 3, no instability was found.

In Fig. 3, the magnetic energies of the dominant toroidal harmonics are shown as a function of time. The magnetic energy of a toroidal harmonic n is defined as $Wmagn≡(∫|Bpoln|2dV)/(2μ0)$⁠, where $Bpoln$ is the poloidal magnetic field contribution given by the harmonic n. Figure 3 shows that when q₉₅ has dropped to a value around two, several toroidal harmonics become unstable. Before this time, a weakly growing n = 1 mode is observed, which is a mix of an external $m/n=3/1$ kink mode and a resistive 2/1 mode. Once the q = 2 surface moves into the open field line region, low-n external kink modes become unstable. In Fig. 4, the mode structures of n = 1, 2, 3 are shown, and it can be seen that these instabilities are associated with the q = 2 surface since the 2/1, 4/2 and 6/3 structures are dominant. Note that this is consistent with external kink modes³² that become unstable when

where q_a denotes the cylindrical edge safety factor; in our case $qa≈1.8<m/n=2$⁠. Higher-n modes are initially excited by non-linear coupling, and higher n modes remain sub-dominant over most of the evolution (see Fig. 5).

All three codes show that the plasma becomes very unstable between 0.85 and 1.10 ms after the plasma becomes limited (see Figs. 3 and 5). Also, there is agreement on the fact that n = 1 is the dominant mode for all but a short period of time. At the beginning of the saturation phase, the n = 2 mode is important and moreover, the n = 3 mode energy can exceed the n = 1 mode energy during short transient phases. Although the mode behavior is similar, differences among the codes of about a factor of 2 are found on the energy saturation level of the dominant n = 1 mode (see also Fig. 7). The growth of the n = 1 mode can also be significantly different as can be inferred when comparing n = 1 time traces of M3D-C¹ with the other two codes in Fig. 3.

The evolution of the thermal pressure during the thermal quench is shown in Fig. 6. The pressure evolution indicates that the external kink modes grow from the plasma edge and start penetrating into the plasma core. The latter process can be observed in the Poincaré plots also shown in Fig. 6. The external modes deform the plasma boundary and create an ergodic layer from about $t′=1.00$ ms onward. Since these modes grow toward the plasma core, several core plasma regions become areas with open field lines. This can be observed in the zones with low density of points in the Poincaré plots, which imply that the field lines typically intersect the wall before completing a full toroidal turn. The regions with open field lines are the ones losing pressure at a faster rate and therefore they have a lower pressure. This is because the thermal energy is quickly lost along the magnetic field lines due to the fast parallel heat conduction. The fast loss of thermal energy is sensitive to the magnetic topology evolution, and in the simulations, it takes 0.14 ms in JOREK, 0.24 ms in M3D-C¹ and 0.18 ms in NIMROD. Therefore, maximum differences of a factor of about 1.7 are found on the obtained duration of the thermal quench. This range of times is consistent with NSTX disruptions where the loss of core confinement occurs in a timescale of $∼0.2$ ms.³³ 

In JOREK and NIMROD, the total wall force is calculated with the expression given by³⁴ 

where n is the normal vector to a closed toroidal surface enclosing the wall. The force in M3D-C¹ is calculated directly with the integral over the vessel elements (⁠$∫J×BdV$⁠). In Fig. 7, the total vertical force, the total horizontal force and the toroidal phase of the horizontal force in the three codes are shown as a function of time (including the magnetic energy of the n = 1 mode for reference). The vertical force (F_z) is already significant before the onset of the non-axisymmetric instabilities, indicating that it is dominated by 2D effects. The magnitude and the time evolution of this force are in fair agreement among the codes. The maximum horizontal force $Fh=Fx2+Fy2$ among the codes varies from 1.3 kN to 3.5 kN and the time of the peak value varies by 0.2 ms. Although the horizontal force F_h appears together with the rise of the n = 1 magnetic energy, their correlation is not obvious. For example, JOREK shows a horizontal force a factor approximately 2.7 larger than NIMROD although the maximum n = 1 magnetic energy is similar. The toroidal phase of the horizontal force indicates that F_h is not able to complete a full toroidal turn in any of the codes after the thermal quench, showing that the force is rotating only very slowly during the phase when asymmetric forces develop.

One of the main features that these codes offer is the capability to calculate self-consistent 3D halo currents. In Fig. 8, the current density normal to the wall is shown as a function of the toroidal angle (⁠$ϕ$⁠) and the R coordinate on the bottom side of the rectangular wall (⁠$Z=−1.4$ m). At time t = 1.04 ms filamentary structures appear implying that the n > 1 harmonics are important to determine the 3D halo current. In particular, a strong n = 3 mode structure is present in the normal current density as this mode has the largest energy at this time in JOREK (see Fig. 5). Note that across the limiter point (⁠$R∼0.8$ m), the current filaments are strong enough to change the sign of the n = 0 current density. Later, when the horizontal force reaches a maximum (t = 1.12 ms), the halo currents show a dominant n = 1 structure. Although NIMROD and M3D-C¹ plots are not displayed, they qualitatively show the same behavior.

Asymmetric VDEs are often characterized by the toroidal asymmetries of the plasma current (I_p). Experimental evidence shows that when measuring I_p at different toroidal locations³⁵ (e.g., using Rogowski coils at different toroidal sectors), VDEs can cause asymmetries on the plasma current of the order of 10%–20% (⁠$ΔIp/Ip$⁠). In Fig. 9, the plasma current is computed at ten equidistant poloidal planes (⁠$ϕ=2πk/10$ with $k=0,1,. . . ,9$⁠) in the three codes. In NIMROD, the maximum I_p asymmetries are of the order of 1.5%–2%, in M3D-C¹ they are on the order of 1%–1.4%, and in JOREK they are smaller than 0.01%.

Although it does not seem to affect key quantities such as the wall forces, the fact that JOREK shows no I_p asymmetries has been investigated as part of this work. The coupling with the resistive wall code¹⁸ does not include the effect of halo currents. This implies that there is a jump of the normal current density across the plasma–wall interface. Although the normal current density arriving to the interface is finite (⁠$J·n|in≠0$⁠) as shown in Fig. 8, the normal current density leaving the interface is zero (⁠$J·n|out=0$⁠). Due to $∇·J=0$⁠, the fact that $J·n|out=0$ automatically implies that $∂ϕIp=0$⁠. The jump on the normal current density is compensated by tangential surface currents flowing on the JOREK's plasma–wall interface that are hidden in the formalism. These currents are not force-free, and therefore, they produce an important contribution to the wall forces (i.e., the halo current contribution).

The JOREK simulation was repeated with a larger number of toroidal harmonics, $n∈[0,20]$ instead of $n∈[0,10]$⁠. The results are shown in Fig. 10, which indicate that high-n mode numbers do not significantly change the evolution of the important figures of merit (such as the plasma current, the thermal energy and the wall forces). NIMROD simulations with $n∈[0,10]$ and $n∈[0,21]$ are also shown in the figure. They also indicate convergence with respect to toroidal representation, but we note that this comparison had been completed with resistivity depending only on the n = 0 component of T and with parallel viscosity, that is, a factor of 100 larger than the scalar viscosity value. While another NIMROD computation with $n∈[0,21]$ and having all the parameters of the three-code comparison shows a similar level of convergence for the NIMROD computations well into the thermal-quench phase, it has not been run sufficiently far to display. We also observe that although the wall force converges with toroidal resolution within each code, the force does not converge between codes.

A check on poloidal resolution is most easily accomplished with NIMROD, where the degree of the spectral-element polynomials can be changed while keeping the same mesh of elements. A lower-resolution computation of the 3D phase with biquadratic elements (156 000 nodes) and the parameters of Sec. IV A produces a peak n = 1 energy that is 37% larger but a peak horizontal force that is only 2.2% larger than the bicubic results.

In this section, we study the origin of the wall forces with JOREK. As it was indicated in Sec. IV A 6, normal electric currents do not enter the resistive wall (⁠$J·n|wall=0$⁠) and, instead, flow along the plasma–wall interface in JOREK. In this sense, the plasma domain's boundary is not force-free. This feature allows one to clearly separate between halo and eddy currents. Since only eddy currents flow on the JOREK-STARWALL wall, the wall force caused by eddy currents can be computed with

where $f≡[(B·n)B−nB2/2]/μ0$ is the magnetic stress tensor projected in the normal direction to the wall. We have used the previous formula to calculate the different contributions to the total wall force as shown in Fig. 11. Although the eddy currents have a significant contribution, the horizontal force (mainly given by F_y) is governed by halo currents. Moreover, the eddy current force is in the opposite direction with respect to the halo current force, and thus it reduces the total horizontal force. In the case of the vertical force, both halo and eddy currents have a similar contribution before the thermal quench. After the thermal quench, the plasma resistivity is significantly increased due to the drop in the plasma temperature. As a consequence, the core plasma current decays and it induces halo currents in the open field line region, which has appreciable temperature due to the imposed B.C. $Te,edge=14.6$ eV. As the induction of halo current and its associated wall force increases at a higher rate than the total wall force (which varies on the resistive wall time), the eddy current contribution to the total force must decrease as observed in Ref. 36. As a general conclusion for both forces, the JOREK simulation shows that the maximum forces are due to halo currents.

In this section, a case with increased parallel thermal conductivity and a case with decreased viscosity are presented. The results shown in Fig. 12 indicate that the viscosity plays a minor role for the evolution of I_p and of the thermal energy. Although not large, reducing the viscosity produces perceptible effects on the wall force. Increasing the parallel conductivity produces a more dramatic effect on the results. The dynamics become faster and the loss of thermal energy takes 0.043 ms instead of 0.14 ms. Although it is expected that the increase in $κ∥$ leads to a faster thermal quench (the heat is conducted faster along the open field lines), it is not clear that the thermal quench time has a linear dependence with the parallel heat conduction coefficient due to the effect of the thermal pressure on the field line topology. In this case when $κ∥$ was increased by a factor of 10 the thermal quench time was reduced by a factor of 3. The magnitude of the horizontal force is significantly larger during the thermal quench although the maximum forces measured are similar (3.5–4.2 kN).

The three codes use different asymmetric perturbations at the start of the 3D phase to excite $n≥1$ activity. JOREK and NIMROD use spatially smooth perturbations, whereas M3D-C¹ uses mesh-scale random noise. Thus, even with perturbations of roughly the same amplitude, the projections onto any growing modes differ. Moreover, initial asymmetric perturbations of steady symmetric fields are typically considered arbitrary in MHD computations, as long as they are sufficiently small. In this application, however, the VDE continually changes the symmetric state, hence the response of the asymmetries. To investigate how the perturbations affect the outcome of our simulations, the NIMROD and JOREK computations have been repeated with amplitudes of the perturbations increased, relative to the computations described in Secs. IV A, IV B, IV C 1, and IV C 2.

Figure 13 compares the evolution of the plasma current, thermal energy, n = 1 and 2 magnetic energies, and the magnitude of horizontal force from NIMROD computations with perturbations at the start of the 3D phase that differ in amplitude by 1000. The small-perturbation case is the one described throughout Sec. IV A. The comparisons show both qualitative and quantitative differences in the results. With the large perturbation, the (2,1) mode reaches large amplitude before other toroidal harmonics become unstable. Therefore, the MHD activity is dominated by the (2,1) mode and its helical harmonics, unlike the small-perturbation case, which shows multi-helicity activity. This is evident in Fig. 14, which shows pressure in the $ϕ=0$ plane, at 1.01 ms of the large-perturbation case and at 1.05 ms of the small-perturbation case. JOREK results also show the dominance of the (2,1) mode structure instead of the multi-helicity activity when starting with large $n≥1$ perturbation amplitudes, as confirmed by the comparison in Fig. 14. Apart from the timing of the thermal quench, the quantitative impacts of starting from the large perturbation in the NIMROD computations include peak n = 1 magnetic fluctuation energy being seven times larger, the thermal quench taking 0.12 ms instead of 0.18 ms, and the peak magnitude of the horizontal force being 80% larger.

We note that both sets of perturbations at the start of the 3D phase are small enough to avoid immediate nonlinear effects. The early-time growth rates of $γτA≈6×10−3$ for the n = 1 mode differ by $≈0.1%$⁠. Using this growth rate with $τA=1.2 μ$ s, we estimate that in the absence of the VDE evolution, even the large perturbation case would not reach saturation amplitude until at least 2.5 ms after the start of the 3D phase. The increasing growth rates with the VDE evolution and the large destabilization that occurs when q₉₅ decreases below two after 0.85 ms make the timing of the disruption somewhat insensitive to the initial perturbation. Nonetheless, our results show that the treatment of asymmetric perturbations can be an important consideration for asymmetric VDE computations. In physics studies, 3D perturbations would normally be present in a self-consistent way either arising from preexisting MHD like neoclassical tearing modes and/or the injected material from a disruption mitigation system such that a dependency on arbitrary initial conditions is not present.

A benchmark case for a three-dimensional vertical displacement event in an NSTX plasma has been presented. The case has been run with the MHD codes JOREK, M3D-C¹, and NIMROD for a three-code comparison. The full run is divided into two phases: an axisymmetric run (2D) and a 3D run. Good agreement was found during the 2D phase for several figures of merit such as the plasma current, the thermal energy and the vertical position (as it was already checked in Ref. 19). The 3D phase was initiated when the plasma became limited by the wall instead of the lower X-point.

In spite of pronounced differences between physics models and numerical methods, a wide range of 3D features predicted by the three codes are in qualitative agreement. For example, the three codes predict that the plasma becomes unstable to low-n external kink modes 0.85–1.1 ms after the plasma becomes limited by the wall. This happens when the q value at the last closed flux surface falls below 2. The growth of these modes lead to the stochastization of the magnetic field lines causing a loss of thermal energy on a timescale of 0.14 ms in JOREK, 0.24 ms in M3D-C¹ and 0.18 ms in NIMROD. During the thermal quench, similar filamentary structures can be observed in the pressure plots for the three codes (see Fig. 6). The total vertical force on the wall is in excellent agreement between the codes, and they all find the predicted maximum 3D horizontal force to be an order of magnitude smaller, in the range of 1.3–3.5 kN. Moreover, the horizontal force is only slowly rotating (less than one toroidal turn) after the thermal quench. Additional JOREK studies reveal that when the forces reach their maximum, they originate from halo currents.

The halo currents show mid-n filamentary patterns with large enough amplitude to reverse the sign of the normal current (⁠$J·n$⁠) given by the n = 0 component (see Fig. 8). At the moment when the force reaches its maximum, the n = 1 component becomes dominant. A scan in the number of Fourier harmonics used indicated that the wall force is sufficiently converged in toroidal resolution $n∈[0,10]$ to provide an adequate description of the modeled physics processes. Additional simulations show that the wall forces are weakly dependent on viscosity and that the choice of the parallel heat conductivity strongly influences the duration of the thermal quench but not the magnitude of the forces.

The nature of VDEs implies that the entire profile and its stability properties evolve over the course of the event. We have found that this leads to both qualitative and quantitative sensitivities to the asymmetric perturbations that are applied to excite $n≥1$ activity. In this case, large perturbations at the start of the 3D phase lead to (2,1)-dominated saturation before other modes are strongly destabilized, whereas small perturbations leave the (2,1) sufficiently small that saturation is a multi-helicity process. This type of sensitivity differs from other nonlinear tokamak stability problems, where the background profile is held relatively constant. The NIMROD and JOREK results that test sensitivity to initial perturbation show the same trends, but we infer that such details at the start of the 3D phase may contribute to the quantitative discrepancies in the peak fluctuation energy and horizontal force that occur much later in the VDE evolution. For physics studies outside this benchmark scenario, initial 3D perturbations would be given by preexisting MHD activity or disruption mitigation such that the results do not depend on somewhat arbitrary initial conditions.

Finally, important differences were observed for the toroidal asymmetry of I_p. The present plasma–wall coupling²⁴ of the JOREK's model is not able to reproduce any I_p asymmetries, which are of the order of a few per cent in NIMROD and M3D-C¹. Nevertheless, the ansatz-based reduced MHD model is able to capture many of the 3D features correctly, even for the large β spherical tokamak plasma considered here.

The consistent results among the three codes bring confidence for their use in disruption studies. Moreover, important post-processing diagnostics were developed and validated during this work. Future efforts will focus on benchmarks of more complex simulations including Ohmic heating, radiation, realistic Spitzer resistivity, impurity injection, neutral particles, and more advanced boundary conditions. We note that in order to have predictive capabilities and to validate the model predictions against the experiment, so that they can be applied with confidence to ITER, these additional effects need to be included (especially radiation and Ohmic heating). Nonetheless, the present work is an important step toward realistic simulations.

This work was supported by the ITER Monaco Fellowship Programme. The authors acknowledge access to the EUROfusion High Performance Computer (Marconi-Fusion) through EUROfusion funding to perform the presented simulations.

ITER is the Nuclear Facility INB No. 174. This paper explores physics processes during the plasma operation of the tokamak when disruptions take place; nevertheless, the nuclear operator is not constrained by the results presented here. The views and opinions expressed herein do not necessarily reflect those of the ITER Organization.

This work was supported by U.S. DOE Contract Nos. DE-AC02–09CH11466 and DE-SC0018001, and the SciDAC Center for Tokamak Transient Simulations. Some of the computations presented in this report used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02–05CH11231. The support from the EUROfusion Researcher Fellowship programme under the task Agreement No. WP19–20-ERG-DIFFER/Krebs is gratefully acknowledged. Part of this work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014–2018 and 2019–2020 under Grant Agreement No. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

The authors would like to thank Guido Huijsmans and Michael Lehnen for fruitful discussions about the contents of this article.

The data that support the findings of this study are available from the corresponding author upon reasonable request.