Experiments and non-linear MHD simulations of hot vertical displacement events in ASDEX-Upgrade

An overview of the evolution of the thermal energy, the plasma current, the vertical position, the edge safety factor q₉₅ and the halo current magnitude I_h is given in figure 1. The initial displacement is triggered by the control coils as described above, but the rate of the displacement differs slightly between the cases. However, the initial growth rate is the same and determined by the resistive decay time of the currents in the structures. Later repetitions of shot # 39655 showed that the time evolution of the axis position can be reproduced with a high accuracy for this setup.

Zoom In Zoom Out Reset image size

The safety factor at the edge (q₉₅) decreases with the displacement as the edge plasma is scraped off at constant current [36]. When q₉₅ reaches a value of around −2.5, −3, −2.5 and −4, respectively (the order in the following is according to table 1), fast MHD instabilities trigger the TQ that leads to a small spike in I_p due to the current flattening. The thermal energy cannot be reconstructed with certainty because it relies on equilibrium reconstruction and is not reliable during disruptions [37]. Instead, the start of the TQ is determined by magnetic signals. The values for the axis position, q₉₅ and thermal energy are obtained from function parametrization [38] based on pre-calculated equilibria, and thus are not fully reliable during disruptions when the plasma is significantly displaced and halo currents are present. The CQ sets in after the TQ with a linear decay time ($\tau_{\mathrm {CQ}} = \frac{t_{20}-t_{80}}{0.6}$) of 2.83, 4.17, 3.17 and 4.0 ms, respectively. An overview of the main characteristics is given in table 2.

Table 2. The current quench time is shorter for the 0.8 MA discharge due to the larger radiation fraction. The halo current fractions (HF) are comparable, while the halo current asymmetries vary more strongly.

Case	$\tau_{\mathrm {CQ}}$ (ms)	HF (%)	TPF	$f_{\mathrm {rad}}$
# 39655	2.82	27.6	1.50	28%
# 39718	4.17	29.4	1.62	18%
# 39720	3.17	34.0	1.68	21%
# 39722	4.00	31.4	1.98	23%

In order to prevent arcing events during disruptions, resistors are installed in the passive stabilizing loop (PSL) [39], which allow for a net current to form in the upper and lower loops at large toroidal electric fields as they occur during the CQ. The halo current reaches values of up to 250 kA and will be characterized more in detail in the following.

The diagnostics for halo currents have already been described in [8] and all derived quantities like the halo fraction (HF) and toroidal peaking factor (TPF) are calculated as defined there. The former quantity is the ratio between the maximum poloidal halo current and the initial plasma current, while the latter is the ratio of the largest halo current in one sector to the average halo current. The lower divertor geometry has changed since the paper was published and the updated geometry is shown in figure 2. Shunt measurements are installed in the upper and lower divertors, the heat shield and in the tiles in front of the PSL. As the plasma is forced to move downward, only the lower divertor measurements are relevant. There are four sectors (4, 10, 12 and 14) with complete poloidal coverage with shunts that were used to calculate all figures of merit unless stated otherwise. In sector 4, the measurement is possible at two tiles to assess toroidal variations over a short distance. Seven of the 64 tiles in the lower divertor have been excluded from the analysis due to malfunctioning. Instead, their values are replaced by the data from the closest tile in toroidal direction. At two poloidal positions in the outer (DUAm) and inner divertor (DUIm), toroidal coverage is good with measurements in 10 out of 16 sectors.

Zoom In Zoom Out Reset image size

2.2.1. Magnitude and toroidal distribution.

In the hot VDE discharges considered here, halo fractions of 28%–34% have been reached, meaning that around one-third of the initial plasma current flows through the wall structures. Small halo current magnitudes are measured already during the late phase of the displacement, but the largest values occur during the CQ phase due to the loss of poloidal flux inside the core [4] as shown in figure 3, where the vertical lines mark the times when I_p has decayed to 90% and 10%, respectively. A toroidal peaking factor (TPF) at the time of maximum halo current of 1.5–2 is observed, which is, as well as the HF, in the typical range for AUG VDEs [8, figure 12]. Observations in various machines have found that TPF × HF stays below a limit of 0.74 [40], which is also the case for these experiments. While the HF is similar for the discharges, the TPF varies by 30% in between them as shown in figure 4(a). Spikes in the halo current occur simultaneously with large values of the toroidal peaking factor (see figure 3) during the CQ phase, which is a sign of MHD activity. The discharges # 39720 and # 39722 show the largest TPF over time and simultaneously a large relative spike in the halo currents.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

When taking only the two poloidal positions with almost full toroidal coverage (DUIm and DUAm) into account, the TPF reaches a value of 3. As only a small poloidal area is taken into account here, this value is not relevant for force measurements but allows to better judge the asymmetry. Figure 4(b) shows the toroidal angle of the sector with the largest current for # 39655. From this, a smooth transition of the peak halo current between the sectors can be observed, which is a sign of a mode rotating on the time scale of the CQ. The direction of the rotation varies between the discharges. However, the spikes in the TPF in # 39720 and # 39722 are mostly due to toroidally localized bursts in MHD activity and not the sign of a rotating mode.

2.2.2. Poloidal distribution.

Figure 5 shows the poloidal distribution of the halo currents over time. At the divertor baffle, close to the contact point, the sign of the halo current reverses over time, indicating a shift of the plasma toward the high field side (HFS). As the measurement only shows the integrated values over the tiles, fluctuations and cancellations can occur. The tiles DUIu, DUMi and DUAu lie in the divertor legs (see figure 2), and thus are shadowed from the plasma, leading to less halo currents in this area.

Zoom In Zoom Out Reset image size

As each of the tiles covers a relatively broad area, it is not straightforward to define a halo width. In fact, a sharp drop in the halo current density occurs at the gap in between the heat shield and the divertor as shown in figure 2. The heat shield practically does not carry any halo current in accordance with the previous observation, that over 90% of the halo currents flow into the divertor [8]. In all of the shots, the outermost tiles still have relevant magnitudes of halo current, while there is no halo current in the heat shield. Therefore, the halo current width does not seem to vary strongly with the plasma current, contrary to what was observed before in hot VDE experiments in COMPASS [19].

ASDEX Upgrade is equipped with foil bolometers and diodes in several sectors [41, 42] to measure the radiated energy of the plasma as well as radiation asymmetries. Figure 6 shows the total radiated energy obtained from the foil bolometers at three different toroidal locations during the disruption. It can be seen that significant radiation asymmetries can occur during hot VDEs. A fast camera in sector 16 confirms a melting event, which is responsible for the toroidally localized radiation. The radiated $E_\textrm{rad}$ energy averaged over the three sectors is 0.38 MJ, 0.14 MJ, 0.16 MJ and 0.08 MJ for the four discharges, respectively. This corresponds of a radiation fraction $f_\textrm{rad}$ varying from 18% to 28% of the total energy, which is lower than previously observed for VDEs in AUG [43]. However, those previous experiments had a larger carbon fraction due to the divertor material. The radiation fraction is defined as follows:

$$\begin{align} f_\textrm{rad} &= \frac{E_\textrm{rad}}{E_\textrm{th}+E_\textrm{mag}} \end{align} \tag{ 1 }$$

$$\begin{align} E_\textrm{mag}&=\frac{1}{2} L I_p^2 = \frac{1}{2} \mu_0 R \left(\mathrm{log} \frac{8R}{a}-2+\frac{l_i}{2}\right) I_p^2, \end{align} \tag{ 2 }$$

where the thermal energy was taken from the diamagnetic diagnostics, R and a are the major and minor radius, respectively, and l_i is the internal inductance. The energies at the time of the VDE start are considered. When the time scale of the heating is comparable to the one of the disruption, it has to be taken into account. However, this is negligible in this Ohmic L-mode case.

Zoom In Zoom Out Reset image size

The total energy of the plasma ($E_\textrm{tot}$) is transferred into radiation ($E_\textrm{rad}$), coupled into the conductive structures ($E_\textrm{coupled}$) and conducted into the PFCs ($E_\textrm{cond}$) [44]:

$$\begin{equation} E_\textrm{tot}=E_\textrm{th}+E_\textrm{mag} = E_\textrm{coupled} + E_\textrm{rad} + E_\textrm{cond}. \end{equation} \tag{ 3 }$$

The coupled energy can be calculated from the change of currents in the vacuum vessel [45], where the energy transferred to the PF coil system has been neglected here. A rough estimate of the latter by considering the change of PF coil currents is an order of magnitude smaller than the other energy sinks, without, however, considering the power supply system. The mutual inductance between the plasma and the conductive structures required for $E_\textrm{coupling}$ has been calculated with STARWALL by replacing the plasma core with a set of axisymmetric filaments. The calculation is performed for an already significantly displaced plasma to take into account the change of distances to the structures during the VDE, but neglecting the further change in mutual inductance during the TQ and CQ phases. The radiated energy has been shown in the previous section. As the measurements of the heat flux on the divertor were not available, the conducted heat is estimated from the missing energy, which makes up to 60% of the total energy for the hot VDE discharges as shown in figure 7. However, the active coil system and the change of the self- and mutual inductance of the plasma have been neglected in this estimation.

Zoom In Zoom Out Reset image size

The vacuum vessel is suspended by eight rods equipped with strain gauges to measure their elongation [46]. This value can be related to a force by Hooke's law given by:

$$\begin{align} F = \varepsilon E A \end{align} \tag{ 4 }$$

with ε (µm m⁻¹) being the elongation, E = 200 GPa the elasticity module and A the cross-sectional area of 552 mm² per rod. However, as the peak of the elongation is delayed with respect to the maximum of the halo and vessel currents, inertia effects play a role. This value can therefore only give an estimation for the real force and a detailed mechanical model of the vessel system would be required for a more realistic result. Additionally, there are measurements of the horizontal and vertical distance between the vacuum vessel at the port location and the supporting structures of the TF coils. The oscillation of the VV, in particular the ports, can thus influence the measurements.

Figure 8(a) shows the vertical force as calculated above for shot # 39722. The vertical force starts with the displacement of the plasma (2.97 s) as eddy currents are induced in the vessel and interact with the magnetic field. In the first steeper phase (2.97 s–2.99 s), the force is due to the variation in the control coils. The force then steadily increases and reaches its peak after the disruption, which is in contrast to the simulation results that will be shown in section 5.4, where the maximum of the force occurs with only a small delay to the peak of the halo currents, indicating that inertia is relevant. Afterward, the sign reverses due to oscillations of the vacuum vessel. The shape of the oscillations is the same for all the discharges, but the amplitude of the maximum scales with the square of the plasma current as shown in figure 9.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The horizontal movement is smaller and the measurements are influenced by vessel and port oscillations. In opposite sectors, the sign of the displacement is reversed, indicating a small global displacement in the range of 0.1 mm. In the discharges # 39655 and # 39720, the data are noisy, whereas in # 39718 and # 39722 (see figure 8(b)), the signal indicates a measurable horizontal shift. In figure 9, the maximum horizontal displacement is shown ($\sqrt{(d_x^2 + d_y^2}$ with d_x and d_y being the maximum of two opposing sectors at each time instant). Thus, the sideways shift increases linearly with I_p , showing that the TPF alone can only be an indicator at constant plasma current. Based on the calculations in [8], the impulse can be obtained from the maximum elongation by assuming that the period of the horizontal force is much smaller than the vessel damping time. A displacement of 0.11, 0.08, 0.06 and 0.05 mm would correspond to an impulse of 126, 90, 72 and 55 Ns, respectively, and scale thus with $I_p B_t$ according to Noll's formula [47].

The VDE is initiated by the virtual position controller in JOREK [49], displacing the plasma by 2 cm in 20 ms with the fast control coils like in the experiment. The resulting coil currents are close to the experimental ones. Afterward, all coil currents are frozen except for the OH coil, which still varies like in the experiment. Only the results of shot # 39655 are presented in this section, but the axisymmetric evolution of the other discharges showed a similar behavior when comparing the results with the experimental observations. The vertical growth rate depends on the initial equilibrium, in particular on the axis position and current profile, as well as on the geometry and resistance of the passive conductors and the SOL temperature. The profile of the latter can be influenced by the parallel and perpendicular heat conduction coefficients, and additionally by the temperature on the boundary of the computational domain. The resulting growth rates of a scan in these quantities for # 39655 are shown in figure 10. A higher parallel or perpendicular heat conduction reduces the temperature in the SOL, so that halo currents are reduced and the plasma moves faster, as it has been studied before [27, 31]. For appropriate parameters, the vertical position evolves consistently to the experiment and also the separatrix shape is consistent with equilibrium reconstructions at different points in time (see figure 11). During the first drift phase, also the PSL current evolution shown in figure 12 evolves with the displacement as the experimental one, showing that the conductive structures are described correctly.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to induce TQ artificially in the axisymmteric simulations, the perpendicular heat and particle diffusion coefficients are increased to a constant value that is two orders of magnitude larger than the initial one. Furthermore, the hyper-resistivity [51] is increased to mimic the redistribution of current as it happens during the fast reconnection phase of the TQ. This leads to a spike in the plasma current, due to the conservation of magnetic helicity [52] as shown in the overview of the simulation results in figure 13. After the thermal energy has decreased to around 10%, the perpendicular heat diffusion is decreased to a value of 2 m² s⁻¹ throughout the domain and the value of the hyper-resistivity is reduced again to a value that does not influence the physical evolution. At the final phase of the CQ, before the last closed flux surface is lost, the latter parameter is increased for numerical reasons, as in reality, MHD activity would lead to a flattening of the current profile, when q₉₅ decreases below a threshold. This will also be shown in section 5.

Zoom In Zoom Out Reset image size

After the TQ, Ohmic heating leads to reheating due to the enhanced resistivity by the reduced temperature and the large current density, which leads to the increase in thermal energy shown in figure 13(a). The post-TQ temperature is established by a balance of Ohmic heating, heat conduction and radiation [53], which requires a knowledge of pre-existing impurities, sputtering and realistic boundary conditions. This is, however, beyond the scope of this work. For VDEs, also the vertical growth rate is important which determines the rate of flux surface opening [44] and is linked with the SOL temperature by the halo currents. A colder SOL does not allow for large halo currents and thus accelerates the vertical motion.

Using circuit equations, one can calculate the temperature that is consistent with the observed CQ times:

$$\begin{equation} L \frac{\mathrm d I}{\mathrm d t} = -R I, \end{equation} \tag{ 5 }$$

where L is a matrix of the self- and mutual inductance of the system, I the currents of the plasma and resistive structures and R the resistance matrix of the conductors and the plasma assuming Spitzer resistivity. An estimate of the CQ time $\tau_{CQ}$ using (5), neglecting the effect of the shrinking plasma cross section, shows that a a temperature in the range of 5–8 eV is expected in order to obtain the same $\tau_{CQ}$ as observed in the experiment. A CQ in JOREK with a predefined uniform temperature confirms this result. Neither impurity radiation nor a realistic heat flux is included due to the missing convection and a reduced parallel conduction as a result of the Dirichlet boundary conditions fixing the boundary temperature at a few electron volts, so that in total, the energy loss is underestimated. Consequently, the CQ is ten times longer than in the experiment.

4.2.1. Heat flux dependence.

To estimate the influence of the heat flux, simulations with a minimum value for the parallel conduction have been carried out. $\chi_\parallel$ is set to the value corresponding to a threshold temperature $T_{\mathrm{min,\chi}}$ if the electron temperature falls below this value:

$$\begin{equation} \chi_{\parallel\mathrm{min}} = \chi_\parallel (T_{\mathrm{min,\chi}}) \ \textrm{if}\ T_e < T_{\mathrm{min,\chi}}. \end{equation} \tag{ 6 }$$

Here, $T_{\mathrm{min,\chi}}$ has been varied in the range of 10 eV–50 eV. Figure 14 shows how $\tau_{CQ}$ reduces when increasing the heat flux to the boundary. At the same time, the movement of the current centroid is faster, which leads to a faster reduction of the region of closed flux surfaces and decrease in q₉₅. The integrated heat flux over time varies from 0.44 MJ to 0.53 MJ, which is close to the simple estimation of section 2.4 of 0.68 MJ. Also, the Poynting flux through the computational boundary of 0.38 MJ–0.24 MJ is in the same order as the estimated coupling energy of 0.28 MJ. For a minimum conduction of 50 eV, $\tau_{CQ}$ is still a factor of 2 larger. When adding a uniform background radiation as another energy sink (crosses in figure 14(b)), $\tau_{CQ}$ can be reduced further. This shows that a realistic heat flux is crucial for modeling the correct SOL temperature and CQ dynamics. Figure 13 shows the thermal energy, the plasma current, the magnetic axis evolution and q₉₅ in comparison to the experiment for different minimum parallel conduction and without background impurity radiation. The halo current magnitude is similar between simulation and experimental measurement and does not vary strongly with the CQ time. Altogether, the evolution of the axis position, q₉₅ and halo current magnitude is consistent with the observations.

Zoom In Zoom Out Reset image size

In other simulations, the resistivity, thus also the Ohmic heating, is increased by a factor of 5 to obtain the same effect, which corresponds to a reduction of the temperature by a factor of 2.9. With the right factor, also the CQ time can be reproduced and the influence on the halo current magnitude or the maximum force is low. The halo fraction varies from 20% to 30% depending on the case. For the case closest to the experimental results, the HF is 26% and the CQ time is 3.6 ms compared to a value of 2.82 ms in the experiment with a HF of 27.6%.

4.2.2. Force calculation.

The vertical force is calculated from the EM stress tensor as introduced in [54, equation (11)]. The order of magnitude is similar compared to the experiment as shown in figure 15, where a force of around 175 kN is reached for # 39655. The forces will be analyzed more in detail in section 5. A variation of the boundary temperature, the perpendicular diffusion profiles and the resistivity lead to a range of CQ times. The shorter the CQ, the more the vertical force is reduced as shown in figure 15 because the vertical moment ($M_\textrm{IZ} = \int j_\phi Z \mathrm{d}R\mathrm{d}Z$) of the current decreases.

Zoom In Zoom Out Reset image size

With the 3D effects now present in the simulations, MHD instabilities arising due to the vertical motion can trigger the TQ self-consistently. The mechanism of the thermal energy loss is described in the following. Like in the 2D simulations, the outer layers of the plasma core are scraped off at almost constant plasma current, resulting in a reduction of the edge safety factor. Furthermore, the current profile steepens at the edge and leads to a triggering of MHD modes. In particular, a 2/1 mode is triggered as well as a smaller, more localized 3/1 mode which can be seen in figure 16. Different to previous observations [28, 31], the mode is not only localized in the edge but reaches into the core. Figure 17 shows how the magnetic topology becomes stochastic starting from the edge, leading to a fast loss of temperature in this region. This becomes evident in the temperature profile in figure 18(a), where the first three time slices correspond to the first three Poincaré plots in the previous figure. As the movement continues, q₉₅ reduces further while the 2/1 mode grows in size (figures 17(b) and (c)) until the core becomes fully stochastic at around 100.64 ms. This increases the rate of thermal energy loss by parallel heat transport along magnetic field lines, which now connect the boundary to the core. This can be seen from the steepening of the slope of the thermal energy loss in figure 18(b).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The TQ takes place on a time scale of $\tau_{{\textrm {TQ}},80-20} = 50\,\mu\text{s}$. The temperature flattens to a final temperature of 250 eV as shown in figure 18 as the parallel heat conduction becomes inefficient at lower temperatures due to the $T_e^{5/2}$ dependence of $\chi_\parallel$. After the TQ, the closed flux surfaces reform gradually in the core (figure 17(e)) and a reheating takes place.

As described in section 4.2, the post-TQ temperature results from a balance of Ohmic heating, impurity radiation and the heat flux at the boundary. The halo current decay is determined by the SOL temperature and contributes to the rate of flux surface opening by the scraping off mechanism. Thus, a small change of the temperature can have a large influence on the CQ time and secondary MHD activity after the TQ. To see the effect of impurities, a scan in tungsten content has been performed varying from 10⁻⁵ to 10⁻³ of the original particle content, where tungsten is included as a uniform background impurity density as the divertor consists of tungsten tiles. In the simplified impurity model used here, the radiation $L_\textrm{rad}$ is given by:

$$\begin{equation} L_{\textrm{rad}} = n_e n_{\textrm{imp}} f_{\textrm{imp}} (T_e, n_e), \end{equation} \tag{ 7 }$$

with the uniform impurity density $n_\textrm{imp}$ and, the temperature T_e and density $n_e$ dependent, linear radiation function $f_\textrm{imp}$ of the impurity based on ADAS data. The effects of the impurities on n_e , $Z_\textrm{eff}$ or η are not taken into account. Furthermore, density has been lost during the VDE and TQ phase due to the boundary conditions, while a constant or increasing electron density is expected due to ionization of impurities. Thus, the radiation is underestimated. Simulations with constant particle content have been repeated, where no large difference in the global dynamics is observed.

Apart from this, the heat conduction is underestimated in these simulations with Dirichlet boundary conditions. Using the Spitzer resistivity and not including radiation or realistic heat fluxes, the CQ time is a factor of ten larger. Thus, a self-consistent inclusion of sputtering and impurity radiation as well as a realistic heat flux is important to perform realistic MHD simulations.

Therefore, a scaling of the resistivity by a factor of 5 is applied as described in section 4, which corresponds to a reduction of the global temperature by a factor of 2.9. Figure 19 shows the evolution of the plasma current and poloidal halo currents in the simulation for two cases with a background impurity density n_W of 3 × 10¹⁴ m⁻³, which is too low to affect the CQ time. In the second case, the resistivity is scaled by a factor of 5, resulting in a faster CQ. It can be seen that the maximum of the halo currents does not vary with the CQ time.

Zoom In Zoom Out Reset image size

5.2.1. MHD events during the CQ phase.

In the experiments, spikes of the halo current during the CQ phase (see figure 3) are a sign of MHD activity. This is also observed in the simulations. Two cases with different amount of background impurity density are compared in figure 20 as an example. In the case with a tungsten density n_W of 3 × 10¹⁴ m⁻³, the edge safety factor decreases to a value below 2 and stays, however, above unity. While in the other case with n_W   =  3 × 10¹⁶ m⁻³, a crash occurs when q₉₅ reaches 2 and q on axis drops below 1, flattening the current profile inside the q = 2 surface and reducing the temperature. As a consequence of the resistivity increase due to the temperature reduction and enhanced by the scaling of η, the temperature recovers fast and a second crash occurs when q₉₅ crosses 2 again. The associated TPFs (figure 23) also behave accordingly. In the case with larger MHD activity, the TPF stays at a higher value for a longer time. On the other hand, in the case with less tungsten, the halo currents become axisymmetric again at the time of largest halo current. While the impurity density is not sufficient to reduce the CQ time, it can influence the MHD dynamics by localized cooling.

Zoom In Zoom Out Reset image size

5.2.2. Simulations at low q₉₅.

The poloidal halo currents in the reduced MHD models are reconstructed from the force balance assumption ($\,j\times B = \nabla p$) and agree well with full MHD models [28]. Note, however, that the calculation of the halo current magnitude differs between experiment and simulation. In the experiment, the average of each divertor side is taken, whereas the axisymmetric halo current magnitude in JOREK is calculated in the following way:

$$\begin{equation} I_h = \frac{1}{2}\iint \mid \mathbf{j}\cdot \mathbf{n} \mid \mathrm{d}S, \end{equation} \tag{ 8 }$$

where the absolute of the current normal to the computational domain ${\mathbf{j}} \cdot {\mathbf{n}}$ is integrated on the boundary. Consequently, the halo current magnitude also includes oscillations on small scales, which are smoothed out in the experimental data. This difference can be seen in figure 19, where the TQ leads to small localized current structures that are not visible in the shunt current signals, but cause a peak of the halo currents in the simulation. However, this does not influence the magnitude during the CQ as the current has a dominant sign on either side of the contact point. Figure 21 shows that the contact point of the plasma with the structures is at the same location between the simulation and experiment, whereas the halo current width in the simulation is larger. These simulations do not contain the detailed geometry of the divertor, which has an influence on the spread of the current. Furthermore, sheath boundary conditions are necessary for a more realistic evolution of the SOL parameters. In particular, the limit on the ion saturation current has an influence on the halo current width [26], but is not accounted for in the present work as this is numerically very challenging in 3D simulations. As the halo current has a similar magnitude as in the experiment, this is not supposed to have a large effect on the global vertical force that only scales with the effective poloidal path l_eff of the halo currents described by [8]:

$$\begin{equation} F_z = I_h B_\phi l_{\textrm{eff}}.\end{equation} \tag{ 9 }$$

However, the halo current width determines local heat loads and forces. Also, the toroidal distribution of the halo currents is of particular interest for the local forces and heat loads.

Zoom In Zoom Out Reset image size

Figure 22 shows the toroidal distribution of the halo currents for the case with a tungsten background impurity density n_W of 3 × 10¹⁶ m⁻³ at different times during the TQ and CQ. During the reconnection phase of the TQ, the current distribution shows localized structures, which are not observable by shunt currents that only provide spatially integrated halo currents. In the local measurements by Langmuir probes, however, a spike at the time of the TQ is measured, showing that those currents also occur in the experiment. The TQ is followed by the CQ phase, where a stationary n = 1 pattern develops in the simulated halo currents that is also present in the eddy currents in the wall (not shown here). The capturing of the halo current rotation requires more realistic boundary conditions [31]. The TPF between the experiment and JOREK is comparable, but its evolution depends on the MHD activity during the CQ phase and decays faster (see figure 23).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is possible to decompose the force in a halo current and passive structure component [28] by changing the surface on which the stress tensor is computed. Figure 24 shows the evolution of the different components, where the eddy current component also includes the PSL contribution. Here, the halo current force is calculated by placing the integration surface just outside the boundary of the computational domain, while the total force is calculated by placing the boundary outside the vacuum vessel. The contribution from the wall currents and PSL is then obtained by substracting the halo contribution from the total one. The vertical force in the 3D simulation is similar to the axisymmetric simulations due to the fact that it depends mostly on the n = 0 evolution of the plasma. As the inertial effects of the vessel are not included in the model, the peak of the force (figure 24) is not delayed with respect to the maximum of the halo currents and also the vessel oscillations, responsible for the negative peak, are not part of the modeling. During the initial displacement phase starting at 2.97 s, only the eddy currents in the PSL and the VV produce a vertical force. At the time of the CQ, the halo currents are the main contributor to the vertical force because the induced wall currents by the I_p spike during the TQ reduce the eddy current contribution. After the plasma current has decayed, the total vertical force is given only by the wall current contribution which now mainly come from currents induced in the CQ.

Zoom In Zoom Out Reset image size

The horizontal force ($F_h = \sqrt{F_x^{\,2} + F_y^{\,2}}$) is two orders of magnitude smaller than the vertical one and not included in the figure. In contrast, JET routinely shows sideway forces in the same order of magnitude as the vertical one [7]. The finding of a small horizontal force in the simulation is consistent with the observation that the horizontal displacement is small in AUG and that the edge safety factor stays above 2 during the CQ phase of the simulation. The impulse of the vessel has been estimated in section 2.5 and is around 126 Ns for the discharge # 39655. The impulse calculated from the 3D simulation is, however, an order of magnitude lower with around 13 Ns, which can be the result of the asymmetries that seem to decay faster than in the experiment and the lack of a kink mode due to the large safety factor throughout the CQ phase. Furthermore, the wall model in STARWALL only considers an axisymmetric wall with isotropic resistivity, so that the wall currents are allowed to flow freely, whereas in reality, the vessel resistivity is not uniform. When estimating a force, the discrepancy is similar. For a statically applied force, a displacement $\Delta r$ of 0.1 mm corresponds to a force of 7.1 kN. The total horizontal force in JOREK is in the order of kN and thus is also smaller.

In addition, the occurrence of the horizontal force in JET is correlated with a toroidal asymmetry of the plasma current [7], which is at the moment not fully reproducible with the current JOREK model as I_p is conserved from toroidal plane to plane as it was shown previously in [28]. The presence of an I_p asymmetry can also not be confirmed in AUG due to the saturation of the poloidal field probes during the disruption.