Magnetohydrodynamic features and particle transport during thermal quench in HL-2A with density source

Disruptions primarily occur due to the rapid onset of global magnetohydrodynamic (MHD) instabilities on various operational limits within tokamaks, such as the Greenwald density limit, beta limit, and current limit.¹ Consequently, magnetic flux surfaces are abruptly destroyed, leading to the breakdown of plasma confinement. The consequences of tokamak disruptions and their related effects have undergone extensive study and evaluation.^1–5 However, in the event of occasional disruptions unexpectedly arising in future fusion reactors, mitigation measures must be implemented.

Presently, the estimated frequency of disruption events stands at 1 in 20 shots for ITER. During unmitigated disruption events in ITER, particularly in the course of full deuterium-tritium operations,³ as much as 350 MJ of thermal energy and up to 1 GJ of magnetic energy could be released and deposited onto the divertor and first wall (FW) surfaces within milliseconds. This could result in severe damage to the device. Consequently, the implementation of disruption mitigation strategies holds critical significance.

A variety of experimental and theoretical approaches have been explored for mitigating tokamak disruptions, with extensive investigation into the underlying physical mechanisms for effective mitigation.^1,4–6 Thermal load management represents one of the foremost concerns in disruption mitigation. During a disruption event, the majority of plasma thermal energy is lost during the thermal quench (TQ) phase through conduction, convection, and radiation, ultimately being deposited onto the divertor and FW surfaces. While impurity injection has been widely employed for TQ mitigation, introducing significant quantities of hydrogen isotopes is also a favorable option as it aids in suppressing runaway electrons. Uniform density increase is essential to prevent runaway electron growth in low-density regions. Moreover, the relatively mild cooling effect resulting from hydrogen isotope deposition is advantageous for reducing the presence of runaway electron seeds in the high-energy tail. In scenarios where mixed deuterium and impurity species injections are utilized, achieving a more even distribution of core density through deeper deposition will enhance the uniformity of radiation heat flux. Notably, significant pre-disruption thermal loss triggered by slow-growing MHD instabilities is good for disruption mitigation.¹ 

As previously mentioned, one of the pivotal aspects of disruptions is the generation of high-flux runaway electrons due to a large toroidal induction field when plasma current is instantaneously extinguished.^4,7 Recent tokamak experiments have elucidated the formation of runaway electrons and the conversion of runaway current subsequent to the onset of disruption in low-density plasmas.⁷ Elevating plasma density is recognized as a means to diminish or even prevent the generation of runaway electrons and the conversion of runaway current. Additionally, avoiding the large J × B forces resulting from in-vessel halo currents during vertical displacement events also involves controlling plasma density.^4,8 Thus, the evolution and control of plasma density during the TQ and current quench (CQ) phases of a major disruption event is of paramount importance in the context of mitigation. Ensuring the timely delivery of substantial numbers of electrons to the plasma core is crucial to avoiding or reducing the avalanche conversion from thermal current to runaway current, while also alleviating the thermal load on the divertor target. This delivery could potentially be achieved through the injection of substantial amounts of hydrogen isotopes.

Numerous disruption mitigation systems have employed various delivery methods, such as pellet injection (PI),⁹ massive gas injection (MGI),^10–14 shattered pellet injection (SPI),^6,15 dispersive shell pellet (DSP),^16,17 and liquid jets.¹⁸ Achieving high-efficiency mitigation relies on factors like deposition location and density mixing, both of which are significantly influenced by MHD activities during external density injection.

In this study, we simulate the TQ processes during the pre-disruption phase to understand the impact of MHD-induced transport on disruption mitigation. The present model is idealized by presetting a pure deuterium (D) source with a constant deposition rate at different locations regardless of D injection process. This means that the simulations here are not directed against the specific delivery methods producing neutral density source. We aim to answer a simple question in this paper: how would different material deposition depths affect the efficiency of core density mixing during TQ core relaxation and associated MHD transport? Simulations are carried out using the 3D nonlinear reduced MHD code known as JOREK. By referencing equilibrium configurations of HL-2A plasmas, the simulations distinctly illustrate varying MHD responses and particle transport characteristics depending on the location of deposition (LoD), whether inside or outside the q = 2 rational surface. Initially, stochastic magnetic fields manifest in the region of the 2/1 magnetic island, expanding outward in the case when the D atoms deposited outside the 2/1 surface and, conversely, expanding inward in the case with the LoD inside the 2/1 surface. Notably, the latter case exhibits considerably higher core density than the others.

The remainder of this paper is organized as follows. The model of the JOREK code and detailed simulation settings are described in Sec. II. Then, simulation results and understanding of MHD activities and transport physics relevant to disruption mitigation are explained in Sec. III. Finally, the achievements of this study are summarized in Sec. IV.

The equilibrium configuration is established through Ohmic limiter discharges of HL-2A plasmas with a toroidal field of B_t = 1.4 T and a plasma current of I_p = 160 kA.²⁰ In the simulation, initial plasma profiles for electron density and temperature, along with ion temperature, are assumed as illustrated in Fig. 1, where T_i = T_e and central values of n₀ = 1.38 × 10¹⁹ m⁻³ and T₀ = 0.8 keV. The expressions for resistivity and parallel thermal conductivity are temperature-dependent, given by $η= η 0 T 0 / T 3 / 2$ and $κ ∥= κ ∥ 0 T 0 / T 5 / 2.$¹⁰ The scaling values of Spitzer resistivity $η 0=5.3× 10 − 7Ω m$ and Spitzer–Härm parallel heat conductivity $κ ∥ 0 = 3.4 × 10 29 m − 1 s − 1$ are chosen, approximately matching measured values of $η ex≈3.5× 10 − 7Ω m$ and $κ ∥ ex ≈ 2.0 × 10 29 m − 1 s − 1$⁠. Additionally, our simulations assume a constant deposition rate of 10²⁴ atoms per second. Note that this is an order of magnitude larger than Ar atoms in real experiments due to weaker electron cooling than impurity atoms.

The simulations commence with solving the Grad–Shafranov equation and obtaining the last closed flux surface (LCFS) and flux function ff ′ from Equilibrium Fitting (EFIT) reconstructed equilibria as inputs. The temperature and density profiles depicted in Fig. 1 replace the EFIT profiles before the D deposition phase. Subsequently, using these inputs, JOREK establishes equilibrium selected from the end of a sawtooth, which is in the steady-state stage, with flux coordinate grids in the poloidal plane, as illustrated in Fig. 2. The detailed simulation setup closely aligns with that described in Refs. 10 and 11. To investigate the influence of the LoD on disruption mitigation, simulations are conducted with varying deposition depths of the D density source. The D density source in the poloidal and toroidal directions follows a Gaussian distribution described by $f= exp − R − R S 2 + Z − Z S 2 Δ r S 2×exp(− ϕ − ϕ S Δ ϕ S 2)$⁠. Here, ( $R S, Z S, ϕ S)=( R S,0.01 m, 0.07 rad)$ signifies the deposition location for D density sources, while $Δ r S=0.05 m$ and $Δ ϕ S=2 rad$ denote the poloidal and toroidal extensions of these sources. Note that $Δϕ=2 rad$ is a wide toroidal deposition width, because the current 3D code is limited by toroidal resolution. This may have no significant effect on radial transport. Modifying the poloidal extension affects the timescale of MHD mode growth, while having limited impact on the maximum amplitude and subsequent transport. Five distinct cases are chosen with LoD $R j=$ 1.87, 1.91, 1.94, 1.97, and 2.02 (m), respectively, along the major radial direction, as shown in Fig. 2 with j = (a, b, c, d, e) corresponding to cases A, B, C, D, and E. The positions of the q = 1, 2, and 3 surfaces are indicated with red, orange, and green dotted curves in Fig. 2. Toroidal harmonics ranging from n = 0 to 5 are incorporated in our simulations. Additionally, it has been verified that case A is convergent, employing n = 0 to 10.

The simulations concerning TQ physics commence with equilibria as established above. As the D density source is activated, MHD fluctuations exhibit exponential growth until the nonlinear stage is attained, during which the n = 1 mode is assumed to play a dominant role. Magnetic islands form and gradually widen, coinciding with an increase in magnetic fluctuations. Subsequently, these islands merge, leading to stochastic behavior. Consequently, magnetic flux surfaces undergo disruption, causing a rapid decline in the core temperature—an event termed TQ onset.

During the later stages of TQ, a distinct surge in magnetic fluctuation energy (⁠ $E ̃ M$⁠) is observable just before the abrupt increase in plasma current (I_p) occurs. Such characteristics are displayed in Fig. 3, where a comparison among different cases involving five distinct LoDs is provided. Notably, apart from a minor disruption with the LoD extending beyond the q = 3 surface (case E), two discernibly different types of disruption behaviors emerge, contingent upon LoD outside or inside the q = 2 surface. Notably, the simulations with LoD at R_c and R_d (cases C and D), positioned between the q = 2 and q = 3 surfaces, exhibit highly sharp magnetic perturbation spikes and subsequent plasma current surges. In contrast, cases A and B with LoDs inside the q = 2 surface display relatively smoother $E ̃ M$ spike and I_p surges. While the amplitudes of spikes in cases D and B might appear similar, the distinct sharp spike in magnetic fluctuation energy (⁠ $E ̃ M$⁠) observed in case D leads to a greater maximum value of the I_p spike. These findings suggest a potential critical role of the q = 2 surface in disruption mitigation. The ensuing section delves into exploring the underlying physical processes and mechanisms contributing to the contrast between these two categories of disruption behaviors. This exploration involves an analysis and comparison of MHD responses, including field line stochasticity, island dynamics, and core particle transport.

Focusing on the physical processes within the pre-TQ and TQ phases, the simulation results are scrutinized, employing case C (LoD outside q = 2) and case A (LoD inside q = 2) for comparative assessment. These two representative cases encapsulate the aforementioned divergence in disruption behaviors—one with LoDs positioned outside the q = 2 surface and the other within it. Subsequently, we outline the TQ progression in these two scenarios, categorizing the temporal evolution into four distinct phases: initial, pre-TQ, main TQ, and late TQ phases, corresponding to the columns in Figs. 4 and 5, spanning from left to right.

Figures 4 and 5 present contour plots illustrating the electron temperature T_e and density n_e, alongside Poincaré plots illustrating magnetic field lines (pointing downward) across a series of time instances (from left to right columns) within each simulation run. In the initial state (t = 0.17 ms), the first column shows distributions as the D density source is just introduced at either R_c or R_a.

Transitioning to the pre-TQ phase (the second column in Figs. 4 and 5), as D atoms undergo ionization, the maximum electron temperature remains largely unaltered as shown in Fig. 6 through $T e(t≤0.26 ms)$ profiles on the midplane. However, the outer region experiences cooling and flattening, leading to lower electron temperature in proximity to outer region.

Simultaneously, the plasma density surrounding the q = 2 surface exhibited substantial augmentation on the low field side (LFS) in both instances. However, their depositions occurred relatively outward in case C and inward in case A, respectively. During this phase, the localized cooling induced by deposition leads to increased resistivity and ensuing current perturbation.^10,21 Consequently, n = 1 magnetohydrodynamic (MHD) fluctuations are activated by resistive tearing and/or kink instabilities, as shown in Fig. 3(a). It is challenging to differentiate between resistive tearing and kink instabilities, and both mechanisms may be at play. Notably, the m/n = 2/1 mode materializes around the q = 2 surface in both cases. The O-point of the magnetic island aligns with the LoD region, aligning with observations from prior studies.^10,11 Furthermore, the m/n = 3/1 magnetic island is also stimulated and subsequently merges with the 2/1 island in case C. In contrast, a more pronounced excitation of the m/n = 1/1 magnetic island is evident in case A, while the flux surfaces remain coherent in the q > 2 region, as distinctly illustrated by the fluctuating T_e contour plots in the initial column of Fig. 8. Of note, field line stochasticity originates primarily in the m/n = 2/1 island region, extending outside the q = 2 surface region in case C or within the q = 2 region in case A.

Transitioning to the main TQ phase, as shown in the third column of Figs. 4 and 5, magnetic energy fluctuations are predominantly governed by the n = 1 component and its nonlinear saturation. Simultaneously, there is a pronounced drop in electron temperature within the core region, as shown in Figs. 6 and 7, while electron density continues to rise near the deposition region of the source. These MHD activities and transport responses are characteristic of the TQ process. During this stage, MHD activities are primarily influenced by the m/n = 2/1 magnetic island, although in case A, the m/n = 1/1 mode assumes greater significance, as evident from the mode structure (see Fig. 8). The region of stochastic magnetic fields rapidly expands inward in case C and outward in case A, encompassing the entire poloidal cross section of the plasma. This divergence is attributed to the earlier onset of m/n = 3/1 instabilities in case C and the larger size of the m/n = 1/1 island in case A.

In the late TQ phase, the TQ process concludes, marked by a substantial reduction in electron temperatures and the emergence of entirely stochastic field lines, as portrayed in the final column of Figs. 4, 5, and 7. Magnetic fluctuation energy experiences a sharp ascent, followed by a sudden decline, resulting in a surge in plasma current (I_p), as demonstrated in Fig. 3. Notably, the current growth time is considerably shorter, and the magnitude of the surge is more pronounced in case C. Additionally, it is worth noting that during this phase, the core plasma density significantly rises in case A, distinguishing it from case C. The ensuing subsection delves into an investigation of the underlying mechanisms.

To summarize, the TQ process, mitigated through externally injected D density sources at different deposition locations, exhibits common characteristics. These include MHD responses predominantly influenced by the 2/1 magnetic island with stochastic magnetic field behavior, complete flattening of electron temperature, the rise in plasma density around the 2/1 surface on the LFS, rapid release of fluctuating magnetic energy, and a surge in plasma current in the later stages. Such shared behaviors align with typical features of the TQ disruption process. Notably, a compelling contrast between case C and case A arises due to the distinct LoDs—outside or inside the q = 2 surface. It is evident that case A experiences significant excitation of the m/n = 1/1 magnetic island, and the magnetic flux surfaces in the q > 2 region remain intact for a notably longer period compared to case C. These differences imply a more intense major disruption in case C relative to case A, suggesting that the latter presents a more effective scenario for disruption mitigation due to improved core density mixing and relatively milder MHD. The initial appearance of field line stochasticity occurs in the outer region encircling the q = 2 surface in case C, in contrast to its inner occurrence in case A.

The underlying mechanisms behind the distinct MHD activity responses and particle transport to different LoDs of the D density source are further elucidated as follows.

First, why does the expansion of magnetic stochasticity in the 2/1 island region occur outward in case C but inward in case A?

The ionization of the D density source leads to an increase in electron density alongside a reduction in electron temperature, ultimately elevating resistivity and introducing variations in plasma current. Consequently, the dissipative m/n = 2/1 mode (likely a visco-resistive tearing mode) becomes unstable in both scenarios, resulting in the formation of a magnetic island. Its O-point region nearly coincides with the LoD on the low magnetic field side. Subsequently, the localized plasma current tends to flatten within the island region, as shown in Fig. 9 at t = 0.22 ms. During this juncture, a 1/1 mode may emerge in both cases, whereas the 3/1 mode manifests solely in case C due to the influence of deposition depth on MHD responses, as evidenced by Fig. 8 and Poincaré plots in Figs. 4 and 5.

The localized flattening of plasma current could possibly signify the linear and/or nonlinear excitation of a resonant mode. To ascertain the role of nonlinear mode coupling in the TQ process, two-dimensional (2D) Fourier transforms are conducted for the (m, n) components of the poloidal flux.

As observed in Figs. 10 and 11, multiple primary mode structures within case C enter a nonlinear phase, encompassing instances (a) and (b) during the main TQ, as well as (c) and (d) nearing the late TQ stage, with n > 1 mode structures also visible in (b) and (d), respectively. It is well established that the 1/1, 2/1, and 3/1 modes arise due to MHD activity. Subsequently, the 2/1 mode is coupled both inwardly (from the q = 2/1 surface to the magnetic axis) with the 1/1 mode and outwardly (from the q = 2/1 surface to the edge) with the 3/1 mode. This interaction results in the generation of the 3/2, 5/2, and n > 2 modes during a nonlinear phase. In the initial stage (t ≤ 0.27 ms) of the main TQ event, the amplitude of the 3/1 mode surpasses that of the 1/1 mode, leading to a more pronounced magnitude of the 5/2 mode compared to the 3/2 mode. Subsequently, as the n > 1 modes situated between the 2/1 and 3/1 modes become more significant than those existing between the 1/1 and 2/1 modes, an overlap of magnetic islands transpires within this region. Consequently, the manifestation of magnetic field line stochasticity originates in this area within case C. Notably, the plasma current profiles in Fig. 9(a) exhibit a collapse mainly outside the core region. This highlights the contribution of the 5/2 mode to the edge collapse.

Furthermore, it becomes evident that in case C, the current profiles depicted in Fig. 9(a) demonstrate that the flattened j_φ,n = 0 region undergoes a radial inward collapse. This, in turn, prompts the successive rapid expansion of the 3/2, 4/3, 5/4, and 6/5 magnetic islands. This progression indicates that initially, the linear excitation of a resonant mode predominates, gradually giving way to the enhanced nonlinear excitation of a resonant mode, particularly evident near the late TQ stage, where the 3/2 mode experiences substantial growth, as demonstrated in Fig. 11(d). As these n > 1 harmonic modes overlap with the 2/1 island, flux surfaces are rendered stochastic beyond the core. Ultimately, the amplitude of the 3/2 mode becomes sufficiently robust to disrupt the isolating magnetic surfaces. Consequently, the region confined within the q = 1 surface rapidly experiences stochastic behavior, thus culminating in the end TQ phase. This suggests that the onset of TQ can be triggered by an avalanche in the current profile, and transport induced by the 1/1 mode. The term “current spike” arises from the global flattening of the current profile during the core relaxation process of the TQ, while the concept of “current profile avalanche” refers to the redistribution of current that propels the region of strong current gradient inward, ultimately leading to a more peaked current profile just prior to the TQ's initiation. It is important to note that this framework for TQ initiation has been proposed in earlier studies.^10,12,22 In light of these simulation findings, the 3/2 mode emerges as a significant contributor to core collapse.

As observed in Figs. 6(b) and 10(b), a 1/1 mode crash might transpire at t = 0.28 ms. This event is heralded by the swift expansion of n > 1 harmonic modes. At the point of the crash, all harmonic modes attain a similar amplitude, aligning with the typical behavior of the nonlinear phase of the internal kink mode. Within this nonlinear phase, four time points are selected, leading to the demonstration of nonlinear mode structures for the primary (m, n) components of the poloidal flux within case A, as shown in Fig. 12. Panels (a)–(d) correspond to different moments at t = 0.31, 0.35, 0.47, and 0.57 ms, respectively. During this progression, the 1/1 mode consistently grows, primarily driven by variations on the 1/1 surface, while the 3/1 mode gains prominence at t = 0.35 ms. Subsequently, the 3/2 mode grows due to nonlinear interactions between the 1/1 and 2/1 modes. Simultaneously, higher-order modes, including n > 2 modes such as 4/3, 5/3, 5/4, and 6/5, witness amplification, with their amplitudes even surpassing that of the 3/2 mode.

However, while magnetic surfaces within the q = 1 surface remain intact, magnetic stochasticity originates within the region situated between the 1/1 and 2/1 surfaces in case A. In contrast, the plasma current profile is flattened inside the 1/1 surface and rearranged within the region between the 1/1 and 2/1 surfaces, maintaining its integrity akin to the q > 2 region in Fig. 9(b) (t = 0.31 ms).

Furthermore, in case A, it is apparent from Fig. 9(a) that a steeper current gradient develops between the 2/1 and 1/1 islands. This raises the question: how is the main TQ phase triggered in this instance? Figures 12(b) and 12(c) show the growth of the 3/1 and 4/1 modes owing to the broadened current profile, subsequently leading to the generation of the 5/2 perturbation through nonlinear coupling of the 3/1 and 2/1 modes. Consequently, these modes (5/2, 3/1, and 4/1) collaborate with the 2/1 mode, fostering the extension of field line stochasticity beyond the q = 2 surface. It is also evident that the current profile experiences a collapse beyond the 2/1 island region at t = 0.47 ms. However, due to the notably swifter growth of the 1/1 mode in this instance, coupled with the delayed excitation of the 3/1 mode, the current profile between the 1/1 and 2/1 islands remains almost constant for $t≤$ 0.47 ms. This contributes to an elongated TQ process (see Figs. 3, 6, and 7), potentially representing a more efficient disruption mitigation scenario.

In the late stages, the 1/1 mode expands further, rapidly followed by significant growth in the stronger 3/2 perturbation due to nonlinear coupling of the 1/1 and 2/1 modes, as shown in Fig. 12(d). Consequently, the magnetic confinement within the q = 1 region is disrupted by the mounting perturbations of the 1/1, 3/2, and 2/1 modes. This prompts field lines in the region to rapidly become stochastic, leading to the expulsion of hot plasma from the X-point and the influx of cold plasma through the O-point.^13,23 This cumulative effect results in the stochastic behavior encompassing the entire region. Consequently, the main TQ phase transpires through a fusion of stochasticity and interactions among the 1/1, 3/2, and 2/1 modes. Parallel to case C, the 3/2 mode is again a pivotal contributor to core collapse, while the 5/2 mode proves crucial for edge collapse. However, the nonlinear excitation of a resonant mode takes on a more dominant role during the main TQ phase in case A.

In the preceding subsection (Sec. III A), we delved into five distinct cases characterized by varying LoDs, as shown in Fig. 2, with a specific emphasis on cases A and C, which correspond to being inside and outside the q = 2 surface, respectively. The evolution of magnetic perturbation energy (⁠ $E ̃ M$⁠) and the corresponding plasma current (I_p) are shown in Figs. 3(a) and 3(b), respectively. Notably, the blue line illustrates a rapid surge in the plasma current (I_p) in case C. This abrupt increase is closely tied to a substantial release of magnetic energy, a characteristic that aligns well with the classical explanation of TQ onset. Conversely, milder I_p surges are observed for other cases, such as case A, owing to weaker MHD activity. Further scrutiny of simulation outcomes underscores the pivotal role of the 3/2 mode in inducing core collapse, while the 5/2 mode holds paramount significance for edge collapse.

In order to grasp the underlying physics of the TQ process, this section takes a closer look at particle transport attributed to D density deposition. Evidently, the rise in core density is heavily contingent on the LoD, a threshold which is evident in Fig. 13(a). Clearly, the core density in case A far exceeds that of the other cases. Consequently, the threshold should fall within the confines of the q = 2 surface. Upon crossing this threshold, the core density experiences a swift escalation akin to case A.

The rationale behind the substantial core density increase in case A is elucidated as follows. To assess the effects of the rational surface, we manipulate the initial equilibrium by shifting the q = 1 flux surface outward. Case F (LoD at R_f = 1.935 m) and case G (LoD at R_g = 1.975 m) correspond, respectively, to the relatively same q = 1 surface positions as case A and case B. It is equally apparent that the increment in core density is contingent on the LoD, relative to the q = 2 surface position, rather than the radial position in meters [Fig. 13(b)]. Consequently, while the volumetric impact does contribute due to the varying LoDs, it should not claim dominance.

These distinct MHD responses can induce disparate mode structures, consequently influencing the efficacy of core density profile enhancement. We can map out the radial mode structure (magnitude distribution) of the poloidal flux and velocity calculated through mφ_im/r (where φ_im represents the imaginary part of the potential, m signifies the poloidal mode number, and r denotes the minor radius) for 1/1 and 2/1 modes (refer to Figs. 11, 12, and 14). Upon comparing the mode structures of case A and case C, even though the flux amplitude of the 1/1 mode is lower than that of the 2/1 mode, the 1/1 mode exhibits a higher velocity than the 2/1 mode in both cases. Notably, the radial width and velocity amplitude of the 1/1 mode are substantially greater in case A than in case C, a characteristic that facilitates the transport of particles to the core. This observation underscores the importance of the velocity of the 1/1 mode in elevating core density.

For a more comprehensive examination of particle transport, the net radial particle flux, denoted as ⟨nV⟩, is presented in Fig. 15, where ⟨…⟩ signifies the flux-surface average, $V=− R 2∇u×∇ϕ+ v | | B$ represents the plasma velocity, and $n≈ n e$ is the plasma density. The first term in the plasma velocity arises from the E × B drift, while the second denotes the parallel component along the magnetic field. Subsequently, we compute and illustrate the two terms separately across three cases to demonstrate that the net particle flux is primarily a result of E × B convection, predominantly exhibiting an inward pinch that exceeds the motion induced by the LoD source toward the core. The average particle flux generated by diffusion and/or convection is further analyzed within the region bounded by the q = 1 flux surface. In this specific region, two cases corresponding to LoDs either inside or outside the q = 2 flux surface are selected. The diffusion coefficient of the stochastic magnetic field is estimated by D_M = D_stc_s, with D_st representing the stochastic magnetic diffusion coefficient previously established in prior work^24,25 and c_s symbolizing the ion sound velocity. The results yield D_M = 9.6 m²/s in case C and D_M = 3.4 m²/s in case A. The heightened D_M in case C compared to case A can be attributed to more pronounced magnetic perturbations in the former. However, the density gradient within the core region is lower in case C than in case A [as seen in Fig. 13(a)]. Consequently, the average particle flux generated by diffusion approximates to be the same in both cases (refer to Table I). $Q n , ( E × B ) r$ and $Q n , ( v | | B ) r$ correspond to the averaged particle flux resulting from the plasma velocity of $− R 2 ∇ →u× ∇ →ϕ$ and $v | | B$ within a minor radius of 0.15 m, as shown in Fig. 15. The analysis suggests that E × B convection attributed to the 1/1 mode takes precedence, as illustrated in Table I.

The outcomes of the five simulation runs, barring the instance of a minor disruption, unveil two distinct TQ process types stemming from varying LoDs. The principal conclusions drawn from this study can be concisely summarized as follows:

1. Certain shared characteristics characterize the TQ process, including MHD responses primarily governed by the 2/1 mode, field line stochasticity arising from island overlap, a complete flattening of the electron temperature profile, a rise in plasma density around the LoD, a rapid release of magnetic fluctuation energy, followed by a sharp surge in plasma current. These findings underline the pivotal role of nonlinear mode coupling in shaping the development of these features, particularly in the ultimate collapse of TQ. Notably, the 5/2 mode is instrumental in edge collapse, while the 3/2 mode plays a vital role in core collapse.

2. A comparative analysis between case C (LoD outside the q = 2) and case A (LoD inside the q = 2) highlights distinct disparities. In case C, magnetic field line stochasticity manifests first in the region exterior to the q = 2 flux surface, whereas in case A, this phenomenon is initially observed in the region within the q = 2 surface, attributable to the LoD-induced cooling effect. Consequently, the current profile progressively collapses inward along the radial direction in case C, while in case A, a steeper current gradient forms between the 2/1 and 1/1 islands during the interval $t≤0.47 ms$⁠. The latter scenario results in an extended TQ process, potentially indicating a more effective disruption mitigation scenario.

3. Furthermore, a discernible trend emerges where the rise in core density is contingent upon the LoD, hinging on the placement inside or outside of the q = 2 surface. Upon the density source crossing this LoD depth threshold, a rapid escalation in core density ensues. While the magnetic perturbation energy release in case C is notably potent, the favorable mode structure observed in case A may hold greater significance for fostering core density augmentation.

4. Distinct mode structures can arise from varying deposition positions, substantially influencing particle transport. Consequently, despite the 1/1 mode exhibiting a lower amplitude than the 2/1 mode, the particle velocity engendered by the 1/1 mode surpasses that resulting from the 2/1 mode in both cases. The amplitude and radial extent of the 1/1 mode in case A significantly surpass those in case C, thereby promoting inward particle transport toward the core. The analysis of net radial particle flux underscores the prevalence of E × B convection in governing particle transport.

The temporal duration of TQ and the enhancement in core density are inherently tied to the LoD, holding significant implications for HL-2A/HL-3 disruption mitigation experiments. These could potentially be appropriated for other fusion devices. The primary focus of this study has been on the LoD. The central objective is to ascertain a more optimal distribution of source deposition locations. Subsequently, the forthcoming endeavor will involve determining how delivery methods can effectively achieve such a distribution.

In future works with more realistic modeling of the injection dynamics, various factors can influence the deepest achievable material deposition. For instance, plasmoids generated following pellet across the discharge may undergo substantial radial drift, commonly referred to as the E × B drift. On the one hand, this drift could be potentially mitigated under certain impurity ratios.²⁶ On the other hand, this drift toward the primary rational surface might be truncated due to the shortened connection length.²⁷ Consequently, as long as the shards penetrate the q = 2 surface, the deposition should ideally be in proximity to or within the q = 2 surface. Additionally, a conceivable influence is the so-called rocket effect,²⁸ where non-uniform ablation of the pellet's front and back caused by the aforementioned plasmoid drift imparts net momentum on the pellet shards, leading to a trajectory shift. These effects could potentially result in relatively shallow penetration. Furthermore, the timescale of the simulated TQ can be qualitatively compared with experimental findings. While both fall within the order of a few tenths of a millisecond, it is important to acknowledge the simulation's smaller scale. Moreover, the simulated TQ's timescale is marginally shorter than the experimental counterpart, likely due to the omission of injection dynamics in the simulation process, or perhaps the influence of $∇ ∥B$ on mode dynamics in a reduced MHD model.²⁹ These aspects warrant further discussion and exploration.

This work was supported by the National Key R & D Program of China under Grant Nos. 2019YFE03010003, 2018YFE0303102, 2022YFE03100004, and 2022YFE03050004. A part of this work was supported by the National Natural Science Foundation of China under Grant Nos. 11905004, 12275071, 12375210, and U1967206 and Sichuan Science and Technology Program under Grant No. 2022JDRC0014. This work was partly carried out on Tianhe-3F operated by NSCC-TJ.

The authors have no conflicts to disclose.

Shilin Hu: Conceptualization (equal); Data curation (equal); Funding acquisition (equal). Di Hu: Investigation (equal); Writing – review & editing (equal). Jiquan Li: Investigation (equal); Writing – review & editing (equal). G. Z. Hao: Investigation (equal); Writing – review & editing (equal). yunbo Dong: Writing – review & editing (equal). Guido T.A. Huijsmans: Writing – review & editing (equal).

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