Complete 3D MHD simulations of the current quench phase of ITER mitigated disruptions

Complete 3D simulations of the current quench phase of ITER disruptions are key to predict asymmetric forces acting into the ITER wall. We present for the first time such simulations for ITER mitigated disruptions at realistic Lundquist numbers. For these strongly mitigated disruptions, we find that the edge safety factor remains above 2 and the maximal integral horizontal forces remain below 1 MN. The maximal integral vertical force is found to be 13 MN and arises in a time scale given by the resistive wall time as expected from theoretical considerations. In this respect, the vertical force arises after the plasma current has completely decayed, showing the importance of continuing the simulations also in the absence of plasma current. We conclude that the horizontal wall force rotation is not a concern for these strongly mitigated disruptions in ITER, since when the wall forces form, there are no remaining sources of rotation.


Introduction
There is no general consensus on the predictions for the maximum integral wall forces acting on the vacuum vessel during ITER disruptions. Although it is broadly accepted that the maximum total vertical force (F z ) will be of the order of 80 MN [29,8], the prediction for the horizontal force (F h ) differs by orders of magnitude. Using the maximum values found in JET for F h (∼ 4 MN) [10] and extrapolating with Noll's formula [22], maximal horizontal forces of 40 MN are obtained. Dedicated studies using a source/sink model inspired by JET measurements [26] also find F h ∼ 40MN [7]. However, other models based on a 1/1 kink mode that interacts with the vacuum vessel only via eddy currents, produce forces an order of magnitude smaller [25]. The ratio between the vessel's current decay time (τ w ) and the current quench time (τ CQ ) was found to play a major role determining the horizontal force with 3D MHD simulations [28]. In that reference, the maximum force for ITER was F h ∼ 30 MN when τ CQ /τ w 1 and the minimum force was F h ∼ 4 MN when τ CQ /τ w 1. Two main effects could be the cause of the strong dependence of the vessel forces with τ CQ . In the following we will refer to the ITER vacuum vessel as "wall" for simplicity.
The first effect is related to the evolution of the edge safety factor (q a ). If q a remains above 2 during a dis-ruption, the 1/1 mode leading to large horizontal forces is not expected. Whether q a remains above 2 during a disruption depends on the competition between the plasma current decay and the vertical motion of the plasma column (i.e. plasma volume shrinkage). If the vertical motion takes place at constant plasma current (i.e. τ CQ /τ w 1), it is expected that q a decreases below unity [2] and that 1/1 modes become unstable. On the other hand, if the plasma current (I p ) decay is faster than the vertical motion, q a increases over time since q a ∝ 1/I p . However, even in the limit where τ CQ /τ w 1, elongated plasmas drift vertically in a time scale given by τ CQ [15,6]. Therefore, the evolution of q a will be ultimately determined by the function that describes the vertical position as a function of the plasma current (Z(I p )), which is dictated by the geometry of the plasma and the wall in this fast current quench limit.
The second effect arises from the penetration time of the magnetic field across the wall (τ w ). It was derived in [24], that the total wall force can be computed with a surface integral enclosing the wall and plasma volumes , J is the current density, B is the total magnetic field and n is a unit vector perpendicular to the toroidal surface. Therefore, for time scales much shorter than τ w , the vessel forces remain small (F ≈ 0) since B remains approximately unchanged outside the wall. However, as it is deduced from (1) and it is shown in this paper, the force can arise after the field penetrates the wall even in the absence of plasma currents (this was also observed in [31]). This fact has been largely ignored in the literature of 3D MHD disruptions and 3D simulations are typically not continued after the plasma current vanishes. In the present work we make sure to extend the simulation time beyond τ w to allow B to penetrate the vacuum vessel.
A major finding of our studies is that, in the τ CQ /τ w 1 limit, wall forces are static in the toroidal direction (i.e. non-rotating). This is due to the fact that when the wall forces are maximum (when the wall currents decay inside the vacuum vessel) the plasma current has already vanished and thus there is no electromagnetic source to drive the rotation. Therefore the issue of resonances between the force rotation and the natural frequencies of the vacuum vessel [27] is not relevant for this limit. To our knowledge this fact has not been identified in the literature so far.
The mitigation of electromagnetic loads in ITER [17] relies on the dissipation of magnetic energy from the plasma by electromagnetic radiation together with the experimental evidence that the wall forces decrease at smaller τ CQ [23]. Consequently, the ITER disruption mitigation system aims to reduce the plasma temperature during the CQ and thereby τ CQ by using massive material injection. In this paper, we demonstrate with 3D MHD simulations featuring realistic time scales, that the wall integral forces are indeed largely reduced for mitigated disruptive plasma conditions with respect to the maximal values extrapolated from simple models and empirical assumptions described above.
In section 2, we present the MHD model and in section 3, we present the simulation setup (initial conditions and parameters). Finally, in section 4, we present our results and we summarize our conclusions in section 5.

MHD model and assumptions
The basic model employed in this paper is a single temperature visco-resistive MHD model [12]  that uses the following reduced MHD ansatz for the plasma flow (v) and the magnetic field (B) where ψ is the poloidal magnetic flux and F 0 = RB φ is taken as a constant. Poloidal currents evolve according to current conservation and momentum balance, but we neglect their contribution to the toroidal field [12]. Successful axisymmetric and 3-dimensional benchmarks of VDEs have been performed with the full MHD codes NIMROD and M3D-C 1 to check the validity of these assumptions for the wall forces [16,5]. The quantities shown in equations (2)-(5) are the magnetic vector potential (A), the ion density (ρ), the total pressure (p),  Figure 2: Profiles at the start of the 3D simulations for the two cases considered in this paper. From top to bottom the profiles of the electron temperature, electron density, poloidally averaged toroidal current profile and safety factor profiles are shown. The radial coordinate is the normalized poloidal magnetic flux.
the total temperature (T ≡ T e + T i ), the electrostatic potential (Φ) and the current density (J). Other parameters in equations (2)-(5) are the plasma resistivity (η), the stress tensor(τ τ τ ), the thermal conductivity (κ κ κ) and the particle diffusion coefficients (D D D) and the ratio of specific heats (γ). The thermal conductivity coefficient tensor κ κ κ presents a high anisotropy (i.e. κ κ ⊥ ) while the particle diffusion coefficients are normally isotropic.
The equations are numerically solved with the fully implicit JOREK code [12,14]. Fourier harmonics are used to represent the toroidal direction and for the poloidal plane, third order quadrilateral Bezier finite elements are used (see figure 1). To take into account resistive wall effects, JOREK is implicitly coupled to the STARWALL code [19,13,6].
Our modeling for ITER mitigated disruption simulations does not consider either Ohmic heating or radiated power by impurities during the current quench. By neglecting them, we assume that their associated terms that should appear in equation (5) exactly cancel each other. In other words, we assume that the plasma magnetic energy is completely radiated. Such an assumption is realistic since radiation efficiencies > 90% are required to mitigate disruptions in ITER DT plasmas [17]. In our case, there are no sources of temperature and density, therefore they are determined by the initial conditions and their evolution through the convection and diffusion/conduction terms.
We model the ITER vacuum vessel composed by inner layer and outer axisymmetric layers using the thin wall approximation (see figure 1). The resistivity and thickness of each layer is given according to the vessel design specifications (η w = 0.8 µΩ and d w = 6 cm). 3D currents are allowed to flow in the inner layer, which is discretized with 75000 thin linear triangles. The outer layer is discretized with toroidal filaments and therefore only axisymmetric toroidal currents flow there. Since we simulate current quench times that are faster than the field penetration time of the inner vessel layer, 3D currents in the outer vessel do not play a role for the plasma dynamics in our case. The outer triangular support (OTS) and the divertor inboard rail (DIR) are also included since they are important passive components for vertical stability. This model for the ITER passive structures has already been benchmarked successfully with the code DINA for axisymmetric simulations [3,11]. Note that we are not taking into account the ITER blanket modules, which could potentially lead to additional stabilizing effects since dipolar currents with characteristic time scales of ∼10 ms can be induced in these modules. However, toroidal currents cannot circulate directly from module to module since they are insulated from each other and only electrically connected through the vacuum vessel. For these reasons, we do not expect that the presence of blanket modules changes the results presented in this work significantly (i.e. the maximum wall forces).
The computational polar grid (see figure 1) is composed of 100×200 radial and poloidal Bezier elements, respectively. The plasma computational boundary roughly matches the ITER first wall but does not represent the details of the divertor structures adequately. More refined boundaries including sharp edges will be included in future work.

Initial conditions and used parameters
The chosen initial conditions correspond to those of a strongly mitigated disruptive plasma in which the ther-  Table 1: Parameters used during the current quench phase. For the temperature dependent parameters the reference temperature is Te0 = 300 eV. Note that except for η and κ the coefficients are spatially constant. * Due to numerical problems, µ had to be increased at some points of the simulation by a factor of 10 (see section 4.4).
mal quench (TQ) has already taken place. A 15 MA / 5.3 T L-mode plasma was considered and the initial conditions were constructed with the following procedure.
In the first place, the pre-TQ L-mode reference equilibrium was computed. Secondly, in order to achieve a large electron density, the density profile was re-scaled by a factor of 20 while the plasma temperature was rescaled by a factor of 1/20. This procedure allows to keep the pressure profile unchanged and to obtain an identical Grad-Shafranov equilibrium as for the reference Lmode but with higher electron densities and lower temperatures. Finally, the perpendicular thermal transport coefficient (κ ⊥ ) was increased in order to simulate an artificial thermal quench, leading to the profiles shown in figure 2 which are used as starting points for the 3D simulations. By adjusting κ ⊥ , the two different cases shown in figure 2 were constructed. These two cases are used to study the sensitivity of the current quench dynamics to differences in the temperature and current density profiles. The resulting temperatures and densities can be expected after the thermal quench following injection of Neon via shattered pellet injection (SPI). Note that the chosen temperatures are close to the Neon radiation peak (∼ 30 eV). Once the profiles shown in figure 2 are established, the thermal diffusion coefficients are reduced to avoid a further temperature decay and the 3D simulation of the CQ starts. During the CQ, 11 toroidal Fourier harmonics are considered ranging from 0 to 10 (n ∈ [0, 10]). The plasma parameters used during the current quench phase are listed in table 1.
The safety factor (q) profiles are not completely flattened inside the q = 2 surface, which are typically observed to be the relaxed MHD states after the thermal quench [21]. These MHD unstable profiles are chosen intentionally to attain a self-consistent relaxation of the q and J φ profiles during the first phase of the CQ. Rather than starting from arbitrarily stable current density profiles, we consider that in order to obtain a J φ profile similar to the one that would be obtained after a TQ, it is more convenient to start with unstable profiles and allow for a self-consistent relaxation.
We apply Dirichlet boundary conditions for the fluid variables at the plasma-wall interface (ρ, T and v). The temperature and density B.C.s are simply n e (t) = 10 20 m −3 and T e (t) = 1 eV. In addition we set a nonormal flow boundary condition at the plasma-wall interface (v·n = 0) and no parallel velocity (v = 0). Note that in our formulation, the v · n = 0 condition is equivalent of that of an ideal wall in the poloidal direction (Φ = 0). It also implies that energy and particles cannot lost through the boundary by convection.

ITER current quench simulations
In this section we describe and compare the two simulated cases. The main difference between the two cases are the different profiles shown in figure 2, which took t = 5 ms and t = 2.1 ms to establish for case #1 and #2 by running JOREK axisymmetrically, respectively. Otherwise the employed parameters are the identical (see table 1).

CQ phase of case #1
For case #1, the axisymmetric run lasted about 5 ms to establish the profiles shown in figure 2. At t = 5 ms the 10 non-axisymmetric modes were initialized to "noise" level and the 3D simulation started. The case is unstable to a variety of tearing modes due to the large plasma resistivity and the initial current profile. These modes are linearly unstable, have similar growth rates (∼ 2 − 3 ms −1 ) and appear at different rational surfaces (q = m/n) described by the toroidal (n) and poloidal periodicities (m) as shown in figure 3. The most unstable modes have low−n toroidal mode numbers and the higher−n modes remain subdominant. The evolution of the poloidal magnetic energy for a selection of dominant mode numbers (n = 1, 2, 3) is shown in figure 4.
At t = 10 ms, the modes energies start to saturate and reach their maximum value at t = 13 ms. Note that at the time of mode saturation, the total plasma current has already decayed to 11-12 MA due to the low plasma temperature. Larger "noise" levels for the initial perturbations could cause the saturation phase to occur at a higher I p . In this respect, more realistic initial conditions for the non-axisymmetric modes could be achieved by also simulating the thermal quench, however this is out of the scope of this work. Later on, we will compare this case to case #2, in which, the modes saturate at a higher plasma current (14 MA).
The mode activity causes a slow flattening of the current profile as it can be observed in the evolution of the internal inductance (l i (3)) in figure 4. The MHD modes lead to the destruction of a significant fraction of the magnetic flux surfaces as shown in the Poincare plots of figure 5. Almost complete ergodization of the field line topology is found at t = 13 ms while at t = 23 ms fluxsurfaces reappear at the plasma core. The confinement of runaway electrons in the 3D perturbed fields produced with this simulation is studied in [30], showing that they will be quickly depleted at t = 13 ms but that they could re-appear in the re-formed core flux surfaces (t = 23 ms).
As I p decays, the plasma moves vertically upwards and toroidal current in the halo region is induced (see I p , Z axis and I halo,φ traces in figure 4). The CQ time defined as τ CQ ≡ (t 20% − t 80% )/0.6 [9] is 47 ms, which is close to the envisaged minimum CQ time for mitigated disruptions (50 ms) [29]. Such evolution leads to an edge safety factor (q 95 trace) that remains in a range of 3-4. The effects of the halo currents in the evolution of q 95 will be discussed in section 4.4.
Since q 95 stays above the value of 3, modes that could result in large horizontal wall forces are not observed (i.e. m/n = 1/1 modes). Accordingly, the horizontal force remains at the noise level (see F 2 x + F 2 y trace in figure  4). The maximum toroidal peaking factor of the poloidal halo currents (TPF trace in figure 4) is 1.2 and arises when the halo current fraction (HF) is still very small (2%). We define the halo fraction as the total poloidal halo current (I halo,θ ) normalized by the pre-disruptive toroial plasma current. The TPF is defined as TPF ≡ max 1 2 φ |J · n|Rdl / (I halo,θ /2π) (8) where J·n is the normal current density into the wall and dl is the differential poloidal length of the wall's contour at a given φ. The maximum value of the HF×TPF product is 0.09, which is significantly lower than the largest products of 0.75 observed in current experiments [9] as it is expected for mitigated disruptions. A value of the HF×TPF product below 0.15 corresponds to category I (frequent occurrence) electromagnetic transients considered for the ITER design [17].

CQ phase of case #2
For case #2, the axisymmetric run lasted about 2.1 ms to establish the profiles shown in figure 2 and the nonaxisymmetric modes were initialized to noise level at that point in time. Similar to case #1, several tearing modes are unstable as shown in figure 3 but here a stronger dominance of the 2/1 mode is observed. The evolution of the magnetic energy for a selection of dominant mode numbers (n = 1, 2, 3) is shown in figure 6.
Contrary to case #1, the saturation of the modes energy takes place at a larger fraction of the initial plasma current (I p ≈ 14 MA at t = 7 ms). The non-linear growth of non-axisymmetric modes causes a large and quick flattening of the current profile as indicated by the l i (3) trace in figure 6. It is noteworthy that the large drop in inductance (∆l i (3) = 0.26 in 3.3 ms) is followed by an increase of the current in the halo region and a spike in the total current. This suggests that the current flattening caused by the MHD activity takes place also beyond the LCFS, and that scrape-off layer currents may be playing a important role when describing the typical I p -spikes that are routinely observed in disruptive plasmas. The distributions of the current and temperature before and after the flattening caused by the MHD activity are presented in figures 7 and 8 respectively.
The resulting current quench time for this simulation is τ CQ = 45.5 ms, which is also close to the minimum allowed by ITER specifications (50 ms). In contrast with case #1, the plasma drifts vertically downwards as I p decays (note the minus sign in front of the Z axis trace in figure 6). We attribute this change of direction to the stronger drop in l i (3) for case #2. It was found in [18] that large drops of the internal inductance during major disruptions in ITER lead to downward displacements. Similarly for the ASDEX-Upgrade tokamak [20], it was found that drops of l i in lower single-null plasmas cause a "vertical dragging effect" that moves the plasma downwards.
Similar to case #1, q 95 remains in the range of 3 − 4 during the CQ and the horizontal force remains at noise level (below 1 MN). The maximum TPF is 1.7 and arises when the halo fraction is 5.7%. The maximum HF×TPF product is 0.17 which is two times larger than in case #1. This is still far from the largest experimental products (0.75), but slightly above the level for category I electromagnetic transients in ITER (0.15).

Wall forces after the current quench
We continue the simulation after the plasma current has completely decayed. During that phase, induced wall currents during the CQ relax and decay, leading to the penetration of the magnetic field across the wall and to the rise of net wall forces. As seen in figure 9, the vertical force arises in a time scale given by the L/R time of the vacuum vessel (∼ 500 ms). Note that the time traces of the vertical forces are very similar to the ones obtained in [1]  distribution of the current density in the inner vacuum vessel is roughly identical (see figure 10). In addition, we checked that the current distribution in the other passive components are also similar for both cases. Since the wall force and the final current distribution is independent of the vertical motion direction, we conclude that the vertical force arises due to the diffusion of the net toroidal wall current along the wall contour and its interaction with the magnetic field produced by the poloidal field coils. Such force would in principle not arise in the case of a fully up-down symmetric plasma centered in an up-down symmetric wall before the CQ starts, which points to the advantage of having up/down symmetric plasmas and reactor designs to minimize disruption forces. However as depicted in figure 1, the plasma and wall structures are not up-down symmetric in ITER and therefore, it is natural that these asymmetries lead to net wall forces due to a net toroidal current. In any case, the magnitude of the maximum F z is a factor 6 smaller than the largest  estimates [29]. Finally the horizontal force completely decays at the moment the plasma current vanishes, confirming that the found horizontal force is noise arising from a finite employed resolution.

Influence of halo currents and other parameters
In our model the temperature at the plasma-wall interface is initialized with an electron temperature of 1 eV. As the plasma is scraped-off due to the vertical motion and as the ergodic field lines connect the plasma core with the halo region, the SOL temperature increases. As shown in figure 8, the temperature can reach values of 10-15 eV in this region during the CQ. Such temperatures are consistent with experimental temperature measurements at the plasma-wall interface during disruptions [4]. At these temperatures, the Spitzer-Harm parallel conductivity (κ ∝ T 2.5 e ) remains moderate and allows significant temperatures / temperature gradients in the SOL to carry the corresponding heat-flux. Estimates for ITER field line lengths of 100 m at T e = 10 eV and n e = 10 20 m −3 give characteristic parallel conduction times of the order of 70 ms.
Parallel convection was found not to be sufficient to decrease the SOL temperature to lower values, although we could not run the case at lower parallel viscosities due to the appearance of numerical instabilities in the v variable. In addition, the parallel viscosity had to be increased by an extra factor of 10 at some points of the simulations to deal with such numerical problems (specifically at t = 14.3 ms for case #1 and at t = 19.8 ms for case #2). However, when comparing re-run fragments of the simulation with and without the increase of µ we did not observe a significant difference in the evolution of Figure 10: Current density distribution in the inner vessel at φ = 0 when the vertical force is maximum. The blue and red arrows represent the current density magnitude of case #1 and case #2 respectively. q 95 , the magnetic energies or the wall forces. Note here that convection can transport thermal energy towards the plasma-wall interface, but this energy can only be lost there by conduction since we set a no-normal flow condition (v · n = 0).
We note that stronger energy sinks in the halo region could decrease the temperature and consequently the total halo currents. In principle, the total halo currents can change the evolution of q 95 . When the halo current flow is limited to very low values due to very low temperatures in the halo, the current that otherwise would be induced in the halo can be partially re-induced in the core, thus slowing down the I p decay with respect to the vertical motion and leading to smaller q 95 values that can potentially drive larger asymmetries in the halo currents. The maximal toroidal halo currents found here (4.1 MA for case #1 and 8.4 MA for case #2) and the maximal poloidal halo current fractions (8.7% for case #1 and 16.6% for case #2) are within (or very close for case #2) to the ITER specifications for category I electromagnetic transients.
More advanced SOL models including Ohmic heating, impurity transport and radiation, neutral particles and sheath boundary conditions are needed to reliably predict the temperatures and densities in the halo region and thus the correct evolution of q 95 . In future works we will pursue such simulations.

Conclusions
We have presented complete 3D MHD simulations of the current quench phase of ITER mitigated plasmas in a comprehensive manner. The chosen parameters include realistic Lundquist numbers (or Spitzer resistivities) and parallel heat conductivity. Due to the high computational cost of such simulations (of the order of 5 million core.hours per simulation) we present two cases featuring an upwards and a downwards vertical displacement. Additional scans will be performed in further follow-up studies.
Both cases show that the integral wall forces will be highly reduced for mitigated plasmas with respect to the largest expected values. We find horizontal forces below 1 MN, with the largest extrapolations from JET being ∼ 40 MN. We attribute the absence of the horizontal forces to the lack of 1/1 modes as q 95 remains above the value of 3 during the current quench. Such beneficial evolution of q 95 originates from the fast current quench time (∼ 50 ms) compared to the wall current decay time (∼ 500 ms) and the induced toroidal halo currents. Further analysis including more advanced models of the halo region will be pursued in future studies. Apart from more refined models of the SOL, we will also pursue studies including the toroidal asymmetries of the ITER vacuum vessel that have not been taken into account in this work. The vessel asymmetries need further assessment since they can potentially impact the current flow and thus the asymmetries of forces during ITER disruptions.
A maximum vertical force of 11-14 MN is found regardless of the direction of the plasma vertical motion. Such force is attributed to the diffusion of the net wall current in an up-down asymmetric vessel and its interaction with the poloidal field coils. The simulated mitigated disruptions show a reduction of the vertical force by a factor of 6 with respect to the largest expected values (80 MN), which confirms the effectiveness of the ITER mitigation strategy and that halo currents are within the category I boundary. For disruptions with CQ times much shorter than the wall's resistive time, our simulations also point to the fact that disruption forces can be minimized if the pre-TQ plasma and vessel structures are up/down symmetric, which has implications for DEMO reactor designs.
We also highlight the potential importance of considering the SOL to describe the I p -spike phenomenon that is typically observed in tokamak disruptions. As it has been shown in this work, the flattening of the toroidal current and temperature profiles due to the MHD activity takes place beyond the Last Closed Flux Surface.
Finally we conclude that the horizontal wall force rotation is not a concern for these highly mitigated plasma disruptions in ITER, since when the wall forces form, there are no remaining sources of rotation (i.e. the plasma current has decayed already).