Model for current drive induced crash cycles in W7-X

In the W7-X experiments with ECCD which we discuss in this article, the net-toroidal-current density j is composed of a contribution from the current drive j_CD and the ohmic plasma response j_screen. Transport calculations using the NTSS code [20] suggest that for these particular experiments the bootstrap current is negligible, and we therefore neglect pressure-driven currents. The wave-driven current j_CD is obtained by the TRAVIS ray-tracing code [19], which accounts for the magnetic geometry in the calculation of both the beam propagation and the current drive. The resulting j_CD profile is strongly peaked near the plasma core in the discharges we study.

Modeling the evolution of the ohmic current profile in complicated geometries (such as W7-X) is an important topic in fusion research which provides predictions for plasma parameters that are difficult to measure experimentally [21–24]. A derivation of the poloidal flux diffusion equation in a low-β stellarator is presented in appendix A.

Following [24], we write the evolution equation for the poloidal flux ψ as

$$\begin{equation}{\sigma }_{{\Vert}}\frac{\partial \psi }{\partial t}=\frac{{J}^{2}{R}_{0}}{{\mu }_{0}\rho }\frac{\partial }{\partial \rho }\frac{{S}_{11}}{J}\frac{\partial \psi }{\partial \rho }-\frac{{V}^{\prime }}{2\pi \rho }{j}_{\text{CD}},\end{equation} \tag{ 3.1 }$$

where V(ρ) is the volume enclosed inside the surface labelled by ρ and the prime indicates a derivative with respect to this coordinate. Furthermore, μ₀ is the permeability constant, R₀ the major plasma radius, σ_∥ denotes the parallel conductivity and S₁₁ is a term of a susceptance matrix defined through the equation (7) of [24]. S₁₁ is calculated using the appropriate VMEC equilibrium. Equation (3.1) is consistent with the ι diffusion equation (22) of [24], where we assumed stationary toroidal flux Φ, i.e. $\frac{\partial {\Phi}}{\partial t}=0$ and $\frac{\partial }{\partial t}\frac{{S}_{12}}{{S}_{11}}=0$. ψ relates to the toroidal current (I) as $\frac{\mathrm{d}\psi }{\mathrm{d}{\Phi}}=\frac{{\mu }_{0}I}{{S}_{11}{{\Phi}}^{\prime }}$, and $J=\frac{{\mu }_{0}F}{2\pi {R}_{0}{B}_{0}}$, where ${B}_{0}=\frac{{{\Phi}}_{\text{lcfs}}}{\pi {a}^{2}}$ and F is the poloidal current associated predominantly with the coils current. Φ_lcfs is the toroidal flux at the last closed flux surface. For volume average beta <0.5% J is within the range (0.996–1) and monotonic. Equation (3.2) represents Ohm's law for the net-toroidal current density

$$\begin{equation}j(r,t)={\sigma }_{{\Vert}}E(r,t)+{j}_{\text{CD}},\end{equation} \tag{ 3.2 }$$

where σ_∥ E = j_screen denotes the ohmic current density and j_CD that from current drive. The neoclassical conductivity σ_∥ is roughly two times smaller than the Spitzer value due to trapped-particle effects at different collision frequencies and is calculated using the NTSS code, as described in [25].

The boundary condition for equation (3.1) corresponds to a perfectly conducting wall located 35 cm from the last closed flux surface. In figure 1 the numerical solution of equation (3.1), j_CD and j_screen, and the sum of currents are plotted at different times.

There are several models for sawtooth crashes in tokamaks [26–30] which involve reconnection of the helical flux of a dominant mode, and we similarly rely on the redistribution of helical flux. The latter is defined by χ = qψ + ϕ, with ψ and ϕ being poloidal and toroidal fluxes, respectively, and q = m/n = 1/ι_res corresponding to the helicity of the dominant instability. We assume that this helical flux is suddenly redistributed when the rotational transform ι reaches a certain value, ι_relax . We shall refer to these events as 'relaxations'.

This model relies on the observation from resistive MHD calculations that, if the amplitude of the current density in the plasma core is high enough that the rotational ι exceeds unity, the corresponding linear eigenmode will be unstable. The helical magnetic field changes sign at the ι = 1 surface and, as the sawtooth models mentioned above propose, equal and opposite amounts of this flux then reconnect. The radial profiles of helical flux before and after such a reconnection event are shown in figure 6. Note that current sheets are formed during relaxation and results in a discontinuous slope. In addition to the n/m = 1 crashes triggered when ι exceeds unity by a sufficient margin, we also postulate n/m = 5/6 crashes when ι is small enough. Thus, whenever ι reaches either of two critical target values ι_relax (in our case: ι_relax > 1 or ι_relax < 5/6, see figure 1, top plot, marked in green), a reorganisation of the corresponding flux is triggered. The features of this model, in which the dominant (m, n) mode helical flux is conserved while the flux profile becomes monotonic, result in: (1) plasma volume conservation and (2) lowering of the magnetic energy. In addition to its effect on the rotational transform, each crash flattens temperature and density profiles over the radial range in which reconnection occurs.

Zoom In Zoom Out Reset image size

Between crashes, the temperature rises gradually in response to the heating, and the rotational transform evolves anew on the resistive time scale until ι again reaches a critical value and the next crash is triggered. In this way, a cycle of two different types of crashes ('small' and 'medium'-sized ones) arises with mode numbers corresponding to the two lowest-order natural resonances associated with the five-fold periodicity of W7-X [2]. This cycle seems to broadly recover the picture observed experimentally.

With this model, we first simulate discharge XP20171206.025 discussed in section 2. At the beginning of this experiment 14 kA ECCD was applied around r_eff ≈ 0.08 m, and temperature crashes started to appear shortly after the start of the ECCD. The ECCD current density was calculated with the TRAVIS ray-tracing code. The critical values of ι_relax were taken to be 1.1 (for crashes associated with the ι = 1 resonance surface) and 0.78 (for crashes associated with ι = 5/6) to match the behaviour of the experimentally observed current evolution (the time of the first crash, frequency of the crashes, the time at which the current exceeds saturation, and, approximately, the final maximum toroidal current). The temperature profile used in our simulation is shown in figure 7.

Zoom In Zoom Out Reset image size

The toroidal plasma current and magnetic energy calculated from the flux diffusion equation with and without relaxations are presented in figure 8. Note that the total toroidal current saturates at a lower value due to continuous magnetic energy dissipation via MHD crashes. Such an early saturation of the toroidal current is often observed in ECCD experiments in W7-X. Our model shows that the loss of magnetic energy via the MHD crashes can have a significant effect on the evolution of the total current. The released energy ends up in plasma thermal energy, but the amount of released energy is too small to have an observable effect in W7-X. Drops in the plasma energy associated with the toroidal plasma current, δW ≈ 0.9–1.2 kJ, corresponding to the medium-sized crashes can be seen in figure 8, (right plot). Small crashes (associated with the 5/6 mode) have no noticeable effect on the magnetic energy.

Zoom In Zoom Out Reset image size

From figure 3 it is evident that even though all the crashes lead to electron temperature flattening, the other crash parameters vary (frequency, amplitude). This important aspect is captured by the model.

A single crash cycle (zoom into figure 8, left plot) is shown in figure 9. In between two medium crashes (associated with the ι = 1 resonance) marked in red, there are 11 small crashes (the 5/6 resonance) marked in green. Vertical lines (red or green) mark the times at which crashes occur. This means that at these times ι reaches ι_relax and relaxation is triggered. After each crash, ι resumes its evolution on the resistive time-scale until the next relaxation occurs.

Zoom In Zoom Out Reset image size

It appears necessary to trigger a relaxation when ι < 1 in order to explain the occurrence of smaller crashes occur between medium ones, as observed experimentally (figure 3). Our numerical experiments show that these smaller crashes could also be due to other 'neighbouring' modes, i.e. for example associated with the ι = 4/5 resonance, which would break the periodicity of the device. Our simplified model is thus incapable to robustly distinguish between these possibilities.

Note that there is a relatively long time gap between the first medium-sized crash and the first small one in figure 9. This gap is also characteristic of the experimental observations 3 and arises because medium crashes flatten the ι-profile, pressing it to unity on the axis, making it take longer for the ι-profile to reach ι_relax = 0.78 again. The subsequent small crashes only relax the ι on axis to 5/6, so it takes less time to reach ι_relax = 0.78 afterwards. After several small crashes ι reaches ι_relax = 1.1 (off axis) and a medium-sized crash occurs again (see figure 10).

Zoom In Zoom Out Reset image size

Another experimentally observed feature is that the frequency of the small events decreases between the medium crashes (figure 3). This is not noticeable in the model (figure 9). The reason for this discrepancy between the model and experiments is likely due to the fact that temperature flattening after a crash is not included in our calculations. One expects that right after a medium-sized crash the temperature in the core drops, which accelerates the resistive diffusion processes. Subsequently the temperature builds up, slowing down the frequency of the small crashes. In fact, the larger crashes seem to exhibit similar behaviour, presumably for the same reason (see, for example, figure 4). In sawtooth experiments on the CTH stellarator [31], it was observed that the sawtooth frequency increases with the fraction of the rotational transform due to external coils (rather than internal plasma current) [32]. This result appears to be consistent with the activity presented in figure 3. However, the toroidal current in W7-X is generally very low (<19 kA) so that the change of the edge rotational transform in the considered cases is insignificant (around 2%).

We also note that the shape of the total current evolution curve (see, for example, figure 9, blue curve) is affected by the resistive evolution of the current sheet arising after medium-sized crashes. An example of such a current sheet can be seen in figure 11.

Zoom In Zoom Out Reset image size

A final important observation is that the mixing area for medium-sized crashes gradually increases (see figure 11). The minor radius of the current sheet associated with these crashes increases in each cycle. This is consistent with the experimentally observed increase of plasma volume where the temperature profiles flatten. Also, the experimentally observed toroidal current evolution (figure 2, right plot) exhibits no wiggles or sawtooth-shaped events in the early stage of the discharge even though small and medium crashes are registered. They appear only later in the discharge. In line with the experimental data our model shows that only later crashes have an effect on the total toroidal current due to the increase of the mixing area.

To apply the theory to W7-X experiments, one needs first to calculate the pre-crash equilibrium parameters and, in particular, the invariants K, ϕ_e. For this purpose, we use the nonlinear ideal MHD equilibrium code VMEC (variational moments equilibrium code) [39], which allows us to produce a sequence of equilibria with varying toroidal current profiles.

The post-relaxation states corresponding to a range of different values of μ (see equation (4.3)) are, by assumption, the possible outcomes of large crashes in W7-X. To identify which state corresponds to a particular initial condition (calculated by VMEC), one finds the one with the same helicity. In the application of equation (4.6), one must ensure that ψ is defined to have the same value on the plasma boundary before and after the relaxation.

The stepped-pressure equilibrium code (SPEC) [40] is well suited for such calculations since it solves the Beltrami equation (4.3) in toroidal geometry as required [41]. SPEC calculates MHD equilibria as extrema of the multiregion, relaxed MHD (MRxMHD) energy functional described in [42]. In MRxMHD, the whole plasma volume Ω is partitioned into a number of sub-volumes, Ω_l , and each plasma volume Ω_l is bounded by the ideal interfaces I_l , which also act as heat transport barriers. The general form of the MRxMHD energy functional, F, is defined as

$$\begin{align}\hfill F& =\sum\limits _{l=1}^{{N}_{v}}\left\{{\int }_{{{\Omega}}_{l}}\left(\frac{p}{\gamma -1}+\frac{{B}^{2}}{2{\mu }_{0}}\right){\mathrm{d}}^{3}\tau \right.\hfill \\ \hfill & \left.\quad -\frac{{\mu }_{l}}{2}{\int }_{{{\Omega}}_{l}}(\mathbf{A}\cdot \mathbf{B}-{K}_{l,0}){\mathrm{d}}^{3}\tau \right\},\hfill \end{align} \tag{ 4.7 }$$

where p is the pressure in each region Ω_l , γ is the specific heat ratio, and K_l is the prescribed value of the magnetic helicity in Ω_l . In addition, the constraints on the enclosed toroidal and poloidal fluxes within each plasma region, Ω_l , and the ideal interfaces satisfying the boundary condition, $\mathbf{B}\cdot \hat{\mathbf{n}}=0$ for infinite conductivity, are considered.

Recently, the SPEC code was upgraded [43] to use Zernike radial basis functions, which makes it possible to calculate well-resolved W7-X equilibria much faster than previously.

While in ideal MHD the magnetic topology is continuously constrained, in MRxMHD the topology is only discretely constrained, thus allowing for partial relaxation. Using the SPEC code, one could solve equation (4.3) in the entire domain without the need for any interfaces between different sub-volumes or one can follow the analysis shown of references [44, 45], where plasma equilibria with internal interfaces were obtained. Here, we apply the first approach and let the plasma redistribute current throughout the whole plasma volume to predict the 'worst possible' outcome of large crashes in W7-X (section 4.3). The second approach may be applied to smaller crashes (as will be discussed in (section 4.4)) and also for the reconstruction and stability analysis of complicated pre-crash states (section 4.2).

The conservation of helicity implies that the pre-crash helicity needs to be calculated in order to use it as a constraint for the calculation of the post-crash state. Rigorous identification of the pre-crash 'equilibrium' is a difficult problem that could, for instance, be done through nonlinear resistive MHD modeling. We take a step in this direction, improving our rather arbitrary ι_relax crash condition (discussed in section 3.3), by performing linear stability analysis within a sequence of MRxMHD equilibria calculated using SPEC for the evolving current profile. The stability of MRxMHD equilibria in SPEC can be assessed via the Hessian matrix $\mathcal{H}$, which describes the second-order variation of the MRxMHD energy functional. By determining the eigenvalues and eigenvectors of the Hessian matrix, $\mathcal{H}$, which satisfy

$$\begin{equation}\sum\limits _{i=1}^{N}{\mathbf{v}}_{i}^{\text{T}}\cdot \mathcal{H}\cdot {\mathbf{v}}_{i}=\sum\limits _{i=1}^{N}{\lambda }_{i}{v}_{i}^{2},\end{equation} \tag{ 4.8 }$$

where v_i and λ_i are the eigenvectors and eigenvalues respectively, one can predict the stability of an equilibrium. Negative eigenvalues imply instability. Linear stability analysis of MRxMHD can reproduce both ideal and resistive MHD stability results [46–48]. General derivations of the SPEC-Hessian matrix in toroidal geometry are discussed in [49].

Figure 13 shows a plot of the largest negative eigenvalue (λ) for the 5/5-mode against the peak ι value, which grows between crashes as a result of the evolution of the current profile. The equilibrium in these calculations has been represented with M_max = 6, N_max = 2 modes per period in the poloidal and toroidal directions, respectively, and by 4 volumes, which allows for a reasonable approximation of the pre-crash ι profile in the SPEC code. This is motivated by the corresponding periodicity of the W7-X plasma but disallows toroidal mode numbers other than 5. Such modes could, in principle, be destabilised by ECCD. This limitation is due to numerical difficulties associated with strongly shaped plasmas with multiple volumes in the SPEC code.

Zoom In Zoom Out Reset image size

The corresponding resonance at ι = (N N_fp)/M = 1 corresponds to a 5/5-island chain. The pre-crash ι-profile presented in the next subsection, figure 14, purple curve, represents an earliest unstable equilibrium, which we will use as an initial pre-crash state.

Zoom In Zoom Out Reset image size

To obtain a realistic initial pre-crash equilibrium with significant toroidal current, we produce a sequence of free-boundary VMEC equilibria with varying toroidal current and toroidal current profiles. The goal is to find an equilibrium with an ι-profile which touches unity (consistent with the stability analysis discussed above) and corresponds to a substantial total toroidal current. The current diffusion process was simulated using the NTSS code [20], and the ECCD profile was taken to match the experimentally observed toroidal current within the model discussed in section 3 (varying amplitude and position of the current drive). Since the pressure gradient in these experiments is small, almost no contribution from the bootstrap current is expected, and this current was thus omitted in the calculations. A result of such calculations is presented in figure 14 purple curve. This case corresponds to an experimental equilibrium with 9.8 kA of net toroidal current.

To cross-validate the results obtained by VMEC and SPEC, we first perform a benchmark of a vacuum (zero total current) W7-X case, figure 14. The agreement between these curves is very good. Note that this is a peculiar case where the single-volume SPEC calculation agrees with VMEC results. The reason for this agreement is that the W7-X vacuum state is a Taylor relaxed state and possesses nearly perfect magnetic surfaces. Generally, one would need a considerably larger amount of volumes in SPEC to avoid formation of magnetic islands (which are not modelled by VMEC).

Conserving the pre-crash helicity K, we run SPEC to find a possible outcome of the non-linear event happening in the plasma. The relaxed ι, figure 14, green curve, is quite flat, monotonically increasing towards the edge, and has a similar shape to the vacuum ι-profile. We remark that the plasma boundaries for vacuum and finite current SPEC calculations are different: we used a sequence of free-boundary VMEC calculations starting from the vacuum ι-profile to describe a pre-crash state with a large toroidal current. Please note, that the SPEC solution for the post-crash state depends on the VMEC reconstruction of the pre-crash state because the plasma boundary, the helicity and the toroidal and poloidal fluxes from this reconstruction are required as input to SPEC.

Also, there is an increase of value of the rotational transform at the edge, ι_edge, after the crash, which indicates an increased total toroidal current. A positive plasma current makes the edge magnetic islands move into the plasma (see figure 15). This result is consistent with the experimental observation of a movement of the divertor strike line in W7-X discharges [50] and can be seen in the step-wise increase of the toroidal current after a crash in figure 12.

Zoom In Zoom Out Reset image size

Using our new post-crash equilibrium, we calculate the total toroidal current in the system:

$$\begin{equation}{\int }_{0}^{2\pi }\mathbf{B}\cdot \frac{\partial \mathbf{r}}{\partial \theta }\enspace \mathrm{d}\theta ={\int }_{0}^{2\pi }{B}_{\theta }\enspace \mathrm{d}\theta ={\mu }_{0}{I}_{\text{tor}},\end{equation} \tag{ 4.9 }$$

where B_θ is the covariant poloidal component of the magnetic field. This gives I_tor ≈ 10.4 kA resulting in a current jump ΔI ≈ 600 A, which is in relatively good agreement with the experiment (figure 12).

In the previous subsection, we considered the very largest crashes in W7-X, which encompass the entire plasma volume and can therefore plausibly be modelled by complete Taylor relaxation. We now return to smaller crashes, previously described in section 3, and explore the question whether they, too, could be described by Taylor relaxation instead of the Kadomtsev-type model proposed earlier. These crashes do not encompass the entire plasma and we therefore consider the Taylor relaxation to happen only in a subset of the full plasma volume. Such models have been proposed in the tokamak literature in order to explain puzzling observations why the safety factor q = 1/ι remains less than unity in the plasma centre throughout the sawtooth cycle [51].

Unlike in section 3.3, where we relied on the free parameter ι_relax to identify the pre-crash state, we now use stability analysis of the equilibrium as a more rigorous identification of the pre-crash state. Whereas this state will thus be identified from linear stability as computed in section 4.2, there is no obvious way to choose the mixing volume, which will simply be taken to be equal to that used in section 3, in order to compare the results of the two models. In both models, the resistive evolution between crashes is calculated as described in section 3.

The purple curve in figure 16 shows the ι-profile after a medium crash, as obtained by the SPEC code using 4 volumes. Please note again, that the SPEC solution depends on the VMEC reconstruction of the pre-crash state, in that the plasma boundary, the helicities and toroidal and poloidal fluxes from this reconstruction are required as input to SPEC. The gray curve shows the corresponding profile, obtained by the VMEC code, at the marginal stability point according to the stability calculation described in (4.2). Note that ι has a local maximum very close to unity. The pre-crash equilibrium has a total toroidal current of 7.4 kA and the minor radius of the mixing region is chosen to be r_mix = 0.23 m. Conserving the total magnetic helicity in this region and keeping ι fixed outside (by prescribing the same K, ψ_e and ϕ_e), we find a post-crash ι-profile which has a discontinuity due to a current sheet, containing a total current δI_sheet = 2180 A, at r = r_mix.

Zoom In Zoom Out Reset image size

In figure 17 we compare the current evolution obtained with the two models for last 4 medium-sized crashes happening around 4.5, 5.0, 5.5 and 5.9 seconds. The blue curve corresponds to the model described in 3.3 with ι_relax = 1.07. This value is different from the one chosen in section 3.3 since here we discuss another discharge. The ECCD position and power were taken to match the experimentally observed toroidal current evolution. The red curve shows the corresponding outcome when partial Taylor relaxation is instead postulated when the plasma becomes unstable. In this case, both the pre-crash stability analysis and the post-crash ι and current profiles are obtained from the SPEC code. Between crashes, the flux diffusion is modelled using the same simplified 1D model until ι-profile becomes unstable again (as discussed in section 4.2).

Zoom In Zoom Out Reset image size

The time traces for the total current in figure 17 are almost the same, but the predictions for the post-crash state are very different. Figure 18 shows a comparison between post-crash ι (left plot) and current j (right plot) profiles (for the crash happening around 5 seconds) obtained with the two different models. Even though both models predict similar current sheet formation (figure 18, right plot) at the boundary of affect region, the shapes of the ι-profiles are quite different (figure 18, left plot). The model based on the reorganisation of helical flux has a strict rule for the magnetic reconnection: flux surfaces with identical value of helical flux (but different direction of the helical field) reconnect. Such reconnection only happens if the pre-crash ι- profile significantly exceeds unity, and results in a post-crash profile that exactly equals unity on the magnetic axis (figure 18, left plot, blue curve). On the other hand, the Taylor-relaxed state is free from these constraints (figure 18, left plot, orange curve) and leads to a flat current profile in the relaxation region. Fully non-linear MHD simulations could perhaps indicate which one of the two scenarios is more realistic.

Zoom In Zoom Out Reset image size