3D non-linear MHD simulation of the MHD response and density increase as a result of shattered pellet injection

PAPER 3D non-linear MHD simulation of the MHD response and density increase as a result of shattered pellet injection To cite this article: D. Hu et al 2018 Nucl. Fusion 58 126025 View the article online for updates and enhancements. Recent citations Interpretive MHD modeling of dispersive shell pellet injection for rapid shutdown in tokamaks V.A. Izzo - Shattered pellet penetration in low and high energy plasmas on DIII-D R. Raman et al - Physics of runaway electrons in tokamaks Boris N. Breizman et al - This content was downloaded from IP address 87.218.223.151 on 13/07/2020 at 12:11 1 © 2018, ITER Organization Printed in the UK 1. Introduction Shattered pellet injection (SPI) is the baseline concept for the ITER disruption mitigation system, the aim of which is to deplete the large thermal and magnetic energy stored within the plasma homogeneously by radiative losses so as to prevent localized energy deposition on the device (e.g. localized heat flux to plasma facing components by plasma deposition or runaway electrons beam strike). Specifically, the ITER SPI system has the capability to inject up to 1025 atoms into the plasma [1]. For the thermal quench (TQ) mitigation, the objective is to deplete the thermal energy by radiation as well as modify the conductive heat flux by dilution, thus mitigate the heat flux to the plasma facing components (PFC), and also to raise the core electron density to be high enough to prevent the hot-tail runaway seeds generation. For the current quench (CQ) mitigation, the objective is to reduce the heat flux through the plasma halo and to reduce electro-magnetic loads, both through appropriate levels of the radiated power. Nuclear Fusion 3D non-linear MHD simulation of the MHD response and density increase as a result of shattered pellet injection D.Hu1, E.Nardon2, M.Lehnen1, G.T.A.Huijsmans2,3, D.C.van Vugt3,a andJET Contributors4,b 1 ITER Organization, Route de Vinon sur Verdon, CS 90 046,13067 Saint Paul-lez-Durance, Cedex, France 2 CEA, IRFM, F-13108 Saint-Paul-Lez-Durance, France 3 Eindhoven University of Technology, De Rondom 70 5612 AP Eindhoven, Netherlands 4 EUROfusion Consortium, JET, Culham Science Centre, Abingdon, OX14 3DB, United Kingdom of Great Britain and Northern Ireland E-mail: [email protected] Received 19 July 2018, revised 2 September 2018 Accepted for publication 25 September 2018 Published 22 October 2018 Abstract The MHD response and the penetration of a deuterium shattered pellet into a JET plasma is investigated via the non-linear reduced MHD code JOREK with the neutral gas shielding (NGS) ablation model. The dominant MHD destabilizing mechanism by the injection is identified as the local helical cooling at each rational surface, as opposed to the global current profile contraction. Thus the injected fragments destabilize each rational surface as they pass through them. The injection penetration is found to be much better compared to MGI, with the convective transport caused by core MHD instabilities (e.g. 1/1 kink) contributing significantly to the core penetration. Moreover, the injection with realistic JET SPI system configurations is simulated in order to provide some insights into future operations, and the impact on the total assimilation and penetration depth of varying injection parameters such as the injection velocity or fineness of shattering is assessed. Further, the effect of changing the target equilibrium temperature or q profile on the assimilation and penetration is also investigated. Such analysis will form the basis of further investigation into a desirable configuration for the future SPI system in ITER. Keywords: tokamak, MHD instability, shattered pellet injection, JOREK, reduced MHD, simulation, disruption mitigation (Some figuresmay appear in colour only in the online journal) D. Hu etal 3D non-linear MHD simulation of the MHD response and density increase as a result of shattered pellet injection Printed in the UK 126025 NUFUAU © 2018, ITER Organization 58 Nucl. Fusion NF 10.1088/1741-4326/aae614 Paper 12 Nuclear Fusion IOP International Atomic Energy Agency a While visiting at ITER Organization, France. b See the author list of [29]. 2018 1741-4326 1741-4326/18/126025+18$33.00 https://doi.org/10.1088/1741-4326/aae614 Nucl. Fusion 58 (2018) 126025 (18pp) D. Hu etal 2 Should a runaway beam form during the current quench, despite the effort to suppress the seed generation, the ITER SPI system can also inject large quantities of argon to dissipate the runaway energy through collisions and line radiation [2]. TQ mitigation will largely determine also the CQ properties, namely the electron temperature and density, and thus the efficiency of CQ mitigation [3]. To achieve both objectives, it is desirable to deliver the material right into the plasma core, since this would result in both a more uniform radiative heat flux to the wall, and a higher core electron density to prevent runaway electron formation. This injection penetration in turn is related to the magneto-hydrodynamic (MHD) response from the plasma, as the MHD modes can play significant roles in the inward transport via the plasma convection and the destruction of flux surfaces. Thus it is of great interest to ensure high experimental availability of the ITER tokamak for experiments to acquire a better understanding about the MHD dynamics and associated injection penetration during SPI. This investigation will also form the basis of future self-consistent prediction of the heat flux onto the plasma facing components in a SPI-mitigated disruption. In this paper, the aforementioned MHD perturbation and density increase are studied by modelling a deuterium SPI into JET target plasmas. JET will be equipped in 2018 with a SPI system that will serve as an important demonstration and extrapolation tool for the future ITER SPI design. Although the most likely injection material for thermal quench mitigation would be a mix of deuterium and neon, here we only look into the injection of pure deuterium to serve as a first look into the interaction between the injection and the ensuing MHD modes. The system of interest is described by the reduced non-linear MHD equationscombined with a diffusive neutral species, solved by the 3D code JOREK [4, 5]. The ablation of the fragments is modelled by the neutral gas shielding (NGS) model [6, 7]. To better demonstrate the principle of the MHD destabilization mechanism caused by the SPI, simulations with equatorial injection will be carried out first. The dominant mechanism will be shown to be the local helical cooling, thus helical current perturbation, at each low order rational resonant surface. Later on, the realistic JET SPI configuration with injection from an upper vertical port will be used to provide insight into the upcoming JET experiments, as well as to demonstrate the impact of various injection parameters on the assimilation and penetration of the injection. The rest of the paper will be arranged as follows. In section2, our system of interest will be described and the governing reduced MHD equationswill be introduced, as well as the NGS model describing the ablation process. In section3, the MHD modes excited by the injection will be investigated. With the understanding gained regarding the MHD instabilities, we proceed to compare the difference between SPI and massive gas injection (MGI) behaviors in section4. We then explore the impact of varying injection parameters with the real JET SPI configuration in section5. Discussion and conclusion regarding the MHD behavior and the implications for future SPI operation will be presented in section6. 2. The system of interest In this section, we introduce our governing equations and assumptions as well as the standard target equilibrium and the injection configurations. 2.1. The governing equationsand the assumptions We model the system by considering the reduced MHD equationscombined with diffusive neutral species [4]. In the tokamak coordinates (R,Z,φ) , the magnetic field and velocity field can be expressed as follows B=F0∇φ+∇ψ×∇φ, (1) v =v B− R 2∇ u ×∇ φ . (2) Here, F0/R is the toroidal magnetic field and F0 is approximately seen as constant in our study, while ψ is the poloidal magnetic flux. Further, u is the flow potential for the E×B flow, v is the parallel velocity scaled by the magnetic field. The governing equationsare then: ∂ψ ∂t =η(Te)∆ ∗ψ − R { u,ψ }− F0 ∂u ∂φ, (3) j=∆ ∗ψ,jφ=−j/R, (4) R∇·  R2ρ∇pol ∂ u ∂t  = 1 2  R2|∇polu|2,R2ρ  +  R4ρω,u  +{ψ,j } −F0 R ∂j ∂φ +ρT,R2+Rµ(Te)∇2ω −∇· ρρ n S ion (T e ) − ρ2α rec (T e )  R2 ∇pol u  , (5) ω= 1 R ∂ ∂R R∂ u ∂R +∂ 2u ∂Z 2 , (6) ∂ρ ∂t=−∇ · (ρv)+∇·  D⊥∇⊥ρ+D∇ρ  +ρρnSion (Te)−ρ 2 αrec (Te), (7) ∂(ρ T ) ∂t=−v·∇(ρT)−γρT∇·v+∇·  κ⊥∇⊥T+κ∇T  +2 3R 2η(Te)j2 − ξionρρnSion (Te) − ρρnPL(Te) − ρ2PB(Te) , (8) ρ B2∂ v  ∂t=−ρ F 0 2R2 ∂ ∂φ  B2v2   −ρ 2R  B2v2 ,ψ  − F 0 R2 ∂(ρ T ) ∂φ + 1 R{ψ,ρT } +B2µ  (T e ) ∇ 2 pol v  +  ρ2α rec (T e ) − ρρ n S ion (T e )  B2v  , (9) ∂ρ n ∂t = ∇· (Dn ·∇ ρn) − ρρnSion (Te)+ρ2αrec (Te)+Sn . (10) In the above equations, equation(3) is the induction equation, equation(4) is the result of Ampère’s law with the permeability absorbed into the current density, equation(5) is the vorticity equation, while equation(6) is the definition of vorticity. Moreover, equation(7) is the continuity equation, equation(8) Nucl. Fusion 58 (2018) 126025 D. Hu etal 3 is the pressure equation, equation(9) is the parallel momentum equation, and finally equation(10) is the diffusive neutral species density equation. Here, we have defined T≡Te+Ti and Te=Ti , we also used ∆∗ ψ ≡R2∇·R−2∇ ψ  , and the dissipative coefficients η∝T − 3 / 2 e and µ∝T − 3 / 2 e are the resistivity and viscosity. Here, due to numerical reasons, we have used an artificially large resistivity which is ten times larger than the Spitzer resistivity. Furthermore κ⊥ and κ∝ T 5 / 2 e are the perpendicular heat conductivity and parallel Braginskii heat conductivity respectively [8]. The parallel and perpendicular plasma diffusion coefficients and the neutral diffusion coefficient are D , D⊥ and Dn respectively, but we used D = 0m 2s−1 in our study, so that the parallel density relaxation is purely carried out by convective flows as the strong convective flow dominate over the diffusion process. As for the coefficients governing the interaction between the plasma and neutrals, Sion(T) is the ionization rate and αrec(T) is the recombination coefficient, the detailed form of which are described in [5]. Further, ξion is the normalized deuterium ionization energy, PL(Te) is the neutral line radiation coefficient and PB(Te) the bremsstrahlung radiation [5]. We assume that the newly ionized deuterium thermalize immediately so that the plasma always remains Maxwellian. The {f,g} in the above equationsdenotes the Poisson bracket with {f,g}≡R(∇f×∇g)·∇φ . Thus equations (3)–(10) form our governing equations. To close the equations, we still have to specify the neutral source term Sn caused by the ablation of fragments. We are not concerned with the shattering process itself and will treat the fragments as they are already generated. To this end, we consider the strongly shielded NGS model in a Maxwellian plasma [6, 7]. The principle of this model is simple, that is, we consider a given heat flux coming down along the field line towards the ablating fragment, the ablation rate of the fragment must be such that it maintains a certain line integrated neutral density along the field line to deplete the incoming heat flux so that the actual flux arriving at the fragment surface is negligible [7]. Hence for a given background electron temperature Te and density ne, the ablation rate for a spherical deuterium fragment with radius rp is ∂ t N(s− 1 )=4.12 × 10 16 r 4 / 3 p (m)n 1 / 3 e (m− 3 )T 1.64 e (eV) . (11) Here, ∂tN is the number of ablated atoms per second. In our model, we deposit those ablated neutrals around the fragment with the following gaussian shape S n∝exp  −(R−Rf) 2 +(Z−Zf) 2 ∆r2 NG  ×exp  −  φ−φf ∆φNG 2. (12) Here, Rf, Zf and φf are the spatial position of fragments, while we choose the neutral cloud parameter ∆rNG =2 cm and ∆φNG =0.6 . This deposition is artificially elongated along the toroidal direction due to limited resolution in toroidal harmonics, but since the most relevant plasma response are those of lower mode numbers, this lack of higher order harmonics is not expected to have a major impact in our investigation. The evolution of fragment size is then governed by the conservation of mass, that is n p4πr2 p ∂ ∂t rp= ∂ ∂t N , (13) with np≃5.958 ×1028 m−3 being the atom density of the deuterium fragments. We treat every single fragment separately, with the ambient density and temperature taken from the surrounding of the fragments. Two non-Maxwellian fluid effects may affect the ablation rate. One is the self-limiting effect arising from the finite amount of hot electrons within the flux tube which results in a truncated tail in the distribution function as ablation goes on, thus decreasing the ablation rate [9]. The other is the slow thermalization of the hot ambient electron population with the cold electrons generated by ionization, which tends to create a long tail distribution, resulting in enhanced ablation compared with the NGS prediction. The two competing effects will begin to manifest themselves once the thermal equilibration time exceeds the timescale of pellet fragments crossing the flux tubes. The comparison of those two timescales will be discussed in section5.5. Furthermore, the fragments are treated as travelling through the plasma without drag. This is a reasonable assumption since the density difference between the fragments and the plasma is very large even after the injection: np / ne > O106 , so that for the drag force to manifest itself within 1 ms, the fragment radius must be approximately smaller than 10−6 m, by the time of which it becomes irrelevant to the plasma conditions and time evolution that we are considering here. 2.2. The target equilibrium and the injection configurations We use the JET pulse No. 86887 as a template for the so called ‘standard equilibrium’, with q0 = 0.935 and q95 = 2.9. The toroidal magnetic field Bt≃2 T, and the total plasma current is Ip≃2 MA. The plasma is in L-mode before injection, with central electron temperature Te(0)≃1.25 keV, and central plasma density ne ( 0 ) ≃2.9 ×1019 m−3 . We have chosen a relatively low core temperature equilibrium since it allows us to use relatively more realistic resistivity and also to better compare with previous MGI results using the same equilibrium [4]. It also serve to mimic the confinement degradation before the disruption [3]. As a reference, the electron temperature, density as well as stored thermal energy of our standard equilibrium are compared with that of a typical JET H-mode with NBI power PNBI ≃18 MW in table1 [10]. No background impurity radiation is assumed. This particular equilibrium is stable to large scale tearing modes ( m⩾2 ), and the toroidal coupling to the m = 2 mode renders the 1/1 ideal internal kink stable for such a low β plasma [11, 12]. The 1/1 resistive kink, although always unstable, is numerically observed to have a natural growth rate γ−1∼O ( 5 ms) , the inverse timescale of which is already longer than the whole injection time we are concerned with here. Thus it is treated as practically stable in the absence of SPI. Nucl. Fusion 58 (2018) 126025 D. Hu etal 4 Midplane cuts of the electron temperature profile, the electron density profile, the pressure profile and the toroidal current density profile are shown in figure1. The ne and Te profiles are generated by using the Thomson scattering data and the equilibrium is constructed by EFIT data as described in [5]. The chained and dashed red lines represent the major radius at the midplane for the q = 2 and q = 1 surfaces, respectively. As mentioned in section1, we consider an equatorial injection to demonstrate the principal mechanism of MHD destabilization by SPI, and a realistic JET SPI configuration to demonstrate the impact of the injection parameters on the injection penetration and assimilation. Sketches of the two injection configurations are shown in figure2, with the red lines outlining the spread cone of the trajectories of the fragments. As a further note, the grid size used in our investigation is 101 in ψ direction and 128 in θ direction for the closed field line region, while the scrape-off layer has a grid size of 4 only. This small number of grid elements covering the open field line region is justified by our emphasis on the core MHD activity. The reference parameters of the above two injection configurations are as follows. For the equatorial injection, the injection is carried out from the low field side (LFS) pointing purely along major radius as shown in figure2(a). This equatorial configuration is only meant to demonstrate the interaction between the injection and the MHD modes then compare with that of the MGI case, as it is easy to tell when do the fragments arrive on a given rational surface. The more realistic investigation is left for the JET-like configuration described later. The total injection amount is 5×1022 particles, equally shattered into 100 fragments each with radius 1.26 mm. The injection speed is 500 ±100 m s−1 with a flat distribution function and a spread vertex angle 40 degrees. It should be noted that the spread angle is unrealistically large in this equatorial case, but numerical investigations with different spread angles and the speed spread shows that, for the values chosen, which are reasonable evaluations of those to be achieved in experiments, they make little difference to the MHD destabilization and consequentially the injection penetration. As for the JET-like injection, the injection is carried from upper LFS and pointing downwards [15] as shown in figure 2(b), and the reference injection direction (the axis of the velocity spreading cone) is within the (R,Z) plane. The total injection quantity is set to be 3.6 ×1022 deuterium atoms, corresponding to the medium sized injection as per the JET SPI design [15, 16]. Moreover, the injected quantity is shattered into 100 fragments with the following size distribution [17] P(rp)= r p K 0(κp r p) I ,I ≡∞ 0 rpK0(κprp)dr=κ−2 p , (14) where K0 is the modified Bessel function of the second kind, and κp is the inverse of the characteristic fragment size which is determined by requiring npNp ∞ 0 P(rp) 4 3πr3 pdrp=N , (15) with Np being the total number of fragments, N being the total injected particles. Thus κ p=  N 6π2n p N p−1 / 3 . (16) Moreover, the fragment velocity is set to be 200 ±40 m s−1, with a vertex angle of 20 degrees. The reference speed and vertex angle are chosen according to JET SPI system design [15, 16], while the distribution of velocity is chosen ad hoc in want of deeper theoretical understanding or experimental observation at present. As a further note, although the characteristic fragment size is related with the injection velocity [15], here we have chosen the two separately due to the lack of quantitative understanding regarding their relationship. Later on in section5, we will deviate from the reference injection parameters for the JET-like configuration to see the impact of varying injection quantities, injection speed and shattering fineness on the MHD activity and the injection penetration. Moreover, equilibria with a different q profile and different electron temperature are also investigated to see the impact on the assimilation rate. 3. MHD response caused by the SPI Macroscopic current driven modes are the major players in the post-injection MHD response due to their global mode structure. Those large scale modes are destabilized by the current density displacement as a result of the drastic electron temperature change after SPI or MGI, which occurs on the local resistive timescale τη∼l2/η , with l being the length scale of said displacement and η∝T−3 / 2 e is the Spitzer resistivity. The community has long established that the propagation of the cooling front along the minor radius, and consequentially the global current contraction contribute greatly to large scale MHD excitations [18, 19], though numerical investigation of deuterium MGI also pointed out the importance of local helical cooling to the growth of corresponding helical modes [4]. Indeed, as will be found in this paper, as long as the fragment travelling timescale is smaller or comparable with the current contraction timescale, the local current perturbation will dominate over the global current contraction as the main MHD destabilizing mechanism during SPI. This is due to the fact that the local perturbation length scale is much smaller than that of the global contraction. Table 1. The comparison of temperature, density, as well as total stored thermal energy before the injection between our standard equilibrium and a typical JET H-mode with NBI power PNBI ≃18 MW [10]. Parameters Standard equilibrium Typical JET H-mode Central Te (keV) 1.25 5.5 Central ne ( 1019 m−3 )2.9 5 Pedestal Te (keV) N/A 1.8 Pedestal ne ( 1019 m−3 )N/A 2.8 Thermal energy (MJ) 0.6 4.5 Nucl. Fusion 58 (2018) 126025 D. Hu etal 5 Due to the mild radiation coefficient from hydrogen isotopes [20], the dominant cooling mechanism in our investigation is the plasma dilution caused by the fragments ablation. The post-injection electron temperature is still on the order of 100 eV during the pre-TQ phase, as can be seen in figure3 where the mid-plane cut of the electron temperature profile evolution for the equatorial SPI is shown. Here by ‘preTQ’ we mean before the final collapse of the core temperature within the q = 1 surface. From this it can be estimated that the timescale for the global current profile to have a 10 cm (a) (b) (c) (d) Figure 1. Profiles at the midplane of the standard target equilibrium for (a) the electron temperature profile, (b) the electron density profile, (c) the pressure profile and (d) the toroidal current density profile. The red chained and dashed lines denote the major radius of the q = 2 and q = 1 surfaces at the midplane, respectively. Figure 2. The sketch of SPI configurations for (a) the equatorial injection and (b) the JET-like injection. Nucl. Fusion 58 (2018) 126025 D. Hu etal 6 contraction along the minor radius is about 10 ms. This is longer than the travelling timescale of the fragments assuming a velocity of 200 m s−1, thus the global current contraction is unlikely to have a major contribution to the MHD excitation. The above statement is supported by the n = 0 current density profile evolution after the equatorial SPI, as is shown in figure4. It can be seen that the mean current density profile does not exhibit strong contraction even just before the onset of the thermal quench, and significant profile variation only occurs after the core current is flattened by non-linear v×B induced hyper-resistivity after the thermal quench [13, 14]. It should be noted, however, that ‘jagged’ features have developed along the minor radius which correspond to the positions of low order rational surfaces, as indicated on figure4 by vertical lines. Those are the result of the local helical cooling as will be shown later in this section. The aforementioned local helical cooling is essentially caused by the geometry of the magnetic field. As the fragments enter the plasma and begin to ablate, they induce a rapid cooling along field lines by parallel heat conduction. The typical timescale of such cooling can be estimated by considering the Braginskii heat conduction [8] of a plasma with 300 eV electron temperature and 1020 m−3 density, and a connection length Lc≡2πRq ≃50 m. The resulting parallel cooling time is then τ∼10−5 s. On irrational surfaces, this will ultimately result in a more or less uniform cooling of the whole flux surface, as the field lines will not connect with themselves. Near rational surfaces, however, field lines connect with themselves after several toroidal turns, and the fragments induce a helical cooling structure instead. This structure will decay on the perpend icular transport timescale, which is about hundreds of microseconds. The resulting n=0 temperature perturbation is shown in figure5, where the dominant m = 2, m = 5 and m = 3 components correspond to the 2/1, 5/3 and 3/2 helical harmonics, respectively. There is also a faint trace of m = 1 component in the plasma core, which is caused by the 1/1 component of plasma displacement. This 1/1 component is likely to be the result of mode beating, such as the beating of the 3/2 and the 2/1 mode. Such a helical cooling structure will induce a corresponding negative helical current perturbation, which is greatly destabilizing for resonant modes. Such a destabilizing effect is due to the above described helical perturbation modifying the local mode structure near the resonant surface in such a way that it lowers ψ s− /ψ s and increases ψ s + /ψ s , resulting in a increase of the stability criterion ∆  ≡ ψ s ψs  + − , driving the tearing instability [21]. Here, ψ s− , ψ s+ and ψs are the radial gradient of the perturbed magnetic flux at both sides of the resonant surface and the value of the perturbed flux at the resonant surface, respectively. This mechanism is essentially the same as the cooling island mechanism proposed by White, Gates and others to explain the sudden growth of islands in density limit disruptions [22, 23], where the radiation within islands results in similar helical cooling structures, causing the destabilization of the islands. Furthermore, the additional bootstrap current profile modification within the island as a result of the pressure profile change can also play a role in the destabilization. This local cooling mechanism implies that the fragments will destabilize successive rational surfaces as they travel across the plasma, generating a broad spectrum of magnetic perturbations. If those surfaces are packed densely enough, the resulting overlapping islands will cause large transport along field lines and thus a significant decrease of plasma confinement. This can be seen from the Poincaré plots of magnetic field lines shown in figure6. The black cross in the figures represents the approximate location of the fragment cloud ‘vanguard’, although there exists some spread both within the poloidal plane and along the toroidal direction. Nonetheless, it can be seen that islands open up as the fragments pass by, and stochasticity follows the vanguard of the fragment cloud closely as they dive into the plasma core. A similar effect is previously reported for the fuelling/triggering pellet triggering of medium-n modes, although there the triggering mechanism is attributed to the local pressure increase rather than the helical cooling [24]. This difference in mechanism is due to the fact that the fuelling/triggering pellet is simply too small to cool the plasma down drastically, thus the helical cooling effect is minimal. The thermal quench is ultimately triggered when the fragments enter the q = 1 surface as shown in figure6(d) and excite the 1/1 kink, which destroys the core confinement completely. 4. MHD modes and injection penetration compared to MGI With the above understanding of the MHD destabilization mechanism of SPI, we can proceed to investigate the MHD spectrum as a result of the injection. In this section, we compare the SPI result with that of a similar quantity MGI. The MGI configuration is the same as the one described in [4], with a small injection quantity of 4.8 ×1021 deuterium atoms. The MGI neutral source is considered as a stationary source at the plasma edge, the spatial shape of which is set to best match with the interferometer data and the temporal shape is set to conform the vacuum solution of gas flow. As for SPI, we Figure 3. The mid-plane electron temperature profile at the beginning of the simulation, just before the thermal quench and just after the thermal quench. The black chained and dashed lines represent the q = 2 and q = 1 surfaces respectively. Nucl. Fusion 58 (2018) 126025 D. Hu etal 7 use the equatorial injection configuration, as described in section2.2, but the total injection quantity is 6.25 ×1021 atoms. Note that there is still some difference between the total injection quantity of the MGI and SPI cases due to historical reasons, but the injected quantities are similar and thus we expect that the differences between SPI and MGI reported here are due to the different material injection schemes and not the the slightly different amounts of material injected. Thus, this is unlikely to alter the following comparison significantly. The magnetic and kinetic perturbation energy of n = 1 to n = 5 harmonics for both cases are shown in figure7. As a reference, the n = 0 mean magnetic energy, which is not shown on the figure, has the order O(1) . For the SPI case, the spectrum of MHD perturbations before the thermal quench is broad, as can be seen from figure7(a) where there is hardly one order of magnitude difference between the magnetic perturbation energy of different harmonics. The thermal quench is triggered at t = 0.9 ms when the fragments penetrate into the q = 1 surface. This penetration time is somewhat longer (0.9 ms compared to 0.7 ms) than that is obtained in section3 and this is due to the reduced amount considered in these simulations, since the ‘vanguard’ fragments are burnt up before they can reach the q = 1 surface. The above behavior is in contrast to that of the MGI case as shown in figure7(b) where the 2/1 mode is dominant, and the thermal quench is triggered at t = 5.2 ms after the island growth is large enough to destabilize the 3/2 mode which then overlaps with it, destroying the flux surfaces [4]. This difference in the MHD response is due to the SPI penetrating much deeper before triggering the thermal quench compared to the MGI case, for which the penetration is limited to the q = 2 surface. Thus, in the MGI case only the 2/1 mode is being destabilized, until it grows to a substantial amplitude to nonlinearly destabilize other modes. This difference in penetration can be readily seen by looking at the density profile evolution and the average density increase within each flux surface, as shown in figure8, where the comparison of the density profile, as well as the average density increase before and after the thermal quench is presented for both cases. Here, the thermal quench occurs from t = 0.844 ms to t = 1.26 ms for the SPI case, while for the MGI case, it happens from t = 5.37 ms to t = 6.07 ms. From figure 8(a), it can be seen Figure 4. The mid-plane current density profile multiplied by the major radius R at the beginning of the simulation, just before the TQ (t = 0.7 ms) and about 1 ms after the TQ. The global current contraction is very limited before the TQ, while after the TQ the core current distribution is flattened by the hyper-resistivity. Figure 5. The n=0 relative electron temperature perturbation after the injection, the black point represents the position of the fragment cloud. Resonant helical cooling is evident at the corresponding rational surfaces, with the most dominant three components being the 2/1, 5/3 and 3/2, and a faint sign of 1/1 component at the core, which is caused by the plasma displacement. Nucl. Fusion 58 (2018) 126025 D. Hu etal 8 (a) (b) (c) (d) Figure 6. The Poincaré plot of magnetic field lines at (a) t = 0 ms, (b) t = 0.245 ms, (c) t = 0.567 ms, (d) t = 0.669 ms. The black cross represents the approximate position of fragments cloud vanguards. (a) (b) Figure 7. The magnetic and kinetic perturbation energy of n = 1–5 modes for (a) the SPI and (b) the MGI. The black chained lines in both cases indicate the time of the onset of the TQ. Nucl. Fusion 58 (2018) 126025 D. Hu etal 15 the other hand, τeq increases with T3 / 2 e , so that this timescale becomes comparable with that of the fragments crossing as the electron temperature approaches 10 keV. This could have strong impact on the ablation since the incoming heat flux ‘seen’ by the fragments is carried by the still hot background electrons with longer mean free path, rather than the thermalized ones. Meanwhile, such a long tail is self-limiting due to the finite amount of hot electrons within a flux tube, as the neutral cloud will absorb those incoming hot electrons and thus tends to truncate the hot tail [9]. The exact result of above competition should be subject to further examination should we investigate the pellet ablation in a Te⩾10 keV plasma. In this section, we look at an equilibrium with twice the electron temperature compared to the standard case we introduced in section 2.2. The toroidal current profile and the density profile remain unchanged. The core temperature is Te(0)=2.5 keV, the central safety factor q(0)=0.943 and the edge safety factor q95 = 2.697. Naively, if the plasma temperature was constant, which would correspond to the no thermalization limit, the ablation rate would have the simple scaling ∂tN∝Te ( t = 0 ) 1.64 according to equation(11). Once thermalization happens, however, this power dependence will not hold. Moreover, once the plasma is cooled down to Te∼O(100 eV) , the mean free path of the plasma electrons is reduced to λe≃O(0.1 m) , so that the core fragments in a fragments’ cloud are effectively shielded from the background plasma by the periphery fragments, reducing any further ablation. The comparison of total assimilation between the high temperature equilibrium and the standard equilibrium for the same SPI configuration is shown in figure17. It can be seen that the increased temperature is indeed increasing the total assimilation, although the power dependence of the latter on the former is less than unity. 5.6. JET-like injection into a rotating plasma Intuitively, the MHD destabilization mechanism detailed in section3 would be compromised in a fast toroidally rotating plasma, as the toroidal rotation will tend to average out the helical cooling effect and result in a more uniformly cooled plasma. In principle, such an effect would begin to manifest when the toroidal rotation frequency is comparable with the inverse timescale of fragments crossing the flux tubes, which can be approximated by τc≃rng/vp , with rng being the radius of the neutral cloud and vp being the fragment velocity. For vp≃200 m s−1 and rng ≃2 cm , the crossing timescale would be τc≃10−4 s, corresponding to a 10 kHz rotation frequency. On the other hand, however, once the islands opened up, ablation of the fragments will more easily increase the density within the island than outside of the island, thus the X-point cooling is small compared to the in-island cooling, causing a net helical cooling structure. This picture is similar to that of the unmodulated ECCD island stabilization, where the unmodulated current drive in the O-point dominates that in the X-point, leading to a net stabilizing effect [27]. Thus, we would expect the toroidal rotation to suppress mode excitation before island formation, but to lose the stabilizing effect in the presence of a significant island. Furthermore, as the fragments ablate, the toroidal rotation will decrease significantly due to momentum conservation, thus removing this stabilization mechanism once substantial density increase occurs. To show this, we investigate SPI cases for which a rigid body plasma rotation is mimicked by artificially moving the fragments along the toroidal direction with rotation frequency f = 10 kHz and f = 1 kHz respectively. As a comparison, we also look into a case where the resistivity and neutral source are artificially set to be toroidally symmetric so that the n = 0 current contraction is the only MHD destabilizing mechanism. The other parameters are the same as for the standard JET-like injection, and the impact of toroidal rotation on the equilibrium has been neglected. The comparison of the MHD response with the non-rotating case is shown in figure18. It can be seen that for the f = 10 kHz rotating case, before the onset of the TQ, the rotational average has significant impact for the n > 1 modes as their amplitude is significantly lower compare with the non-rotation case. Beginning from t = 2.4 ms however, significant growth of n > 1 modes occurs which ultimately leads to the first TQ at t = 3.0 ms. This is then followed by the final collapse of the core temperature at t = 3.50 ms. Since the two TQ are very close to each other, they are nearly indistinguishable on figure18(b). For the slowly rotating case, the MHD spectrum looks more or less the same compared with the non-rotating one, and the time at which the first and second TQ occur is also similar, suggesting the toroidal rotation is not important once it’s slow enough. For the n = 0 current contraction only case, the TQ is triggered after a significant amount of current is contracted into the q = 1 surface at t = 3.50 ms. Comparing figures18(b) and (d), it is evident that the second TQ of the fast rotating case happens at the same time as the TQ in the n = 0 current contraction case, suggesting the current contraction plays an important role in the final collapse of core temperature. Meanwhile, the helical effect is suppressed during the early part of the injection and only manifests itself very close to the final TQ, when the islands begin to form. Thus, Figure 17. The comparison of the total assimilated particles for JET-like SPI into the standard and high Te equilibria. The higher temperature is seen to be beneficial for the assimilation, though the power dependence of total assimilation on the initial temperature is less than one. Nucl. Fusion 58 (2018) 126025 D. Hu etal 16 it can be said that once the toroidal rotation is faster than the inverse time scale of fragments crossing the flux tube, then the rotation will indeed suppress the helical destabilizing effect, but as figure18(c) has shown, once the rotation has slowed down, it is not important anymore. As a further note, here we are considering injection into an initially MHD free rotating plasma, injection into a rotating plasma with existing modes requires further investigation in the future. 6. Summary and conclusion The MHD instabilities and the density response during deuterium SPI into a JET L-mode plasma has been investigated in this work. The main focus has been on the macroscopic current driven modes that are responsible for triggering the thermal quench and for the convective density mixing. The evolution of the plasma following injection of deuterium by SPI and MGI has been investigated using the 3-D non-linear reduced MHD code JOREK, combined with the strongly shielded NGS model to describe the fragments ablation. It is found that the MHD destabilization by the deuterium SPI is dominated by the the local helical cooling and current perturbation instead of the n = 0 current contraction. This helical effect is driven by the almost adiabatic local decrease of the electron temperature near rational surfaces due to the pellet ablation/ionization. Hence, the SPI fragments destabilize successive layers of rational surfaces as they fly into the plasma core, result in broad-spectrum MHD perturbation and widespread field line stochasticity in the wake of the fragments. It is also found that the SPI shows superior penetration after the TQ compared with the MGI even when the fragment trajectories do not cross the magnetic axis, as the former deposit the injected material right into the q = 1 surface before and during the onset of the TQ, where the 1/1 kink convection is most efficient. This results in sufficient core density mixing. As a consequence, the SPI enjoys a shorter time difference between the core temperature and density relaxation compared with its MGI counterpart, although some localized regions with cold, hollow density and good flux surfaces still exist, which may be vulnerable to runaway electron formation. Further investigations reveals several important injection parameters which impact on the penetration and assimilation of injected deuterium. The shattering fineness of the SPI is found to increase the total assimilation of injected atoms, (a) (b) (c) (d) Figure 18. The comparison of the perturbed magnetic and kinetic energy of n = 1–5 modes for the JET-like SPI in (a) a stationary plasma, (b) a rigid body rotating plasma with rotation frequency 10 kHz, (c) a rigid body rotating plasma with rotation frequency 1 kHz and (d) a plasma where the n = 0 current contraction is the only MHD destabilizing mechanism. The red and blue chained lines indicate the onset time of the first and the second TQ. Nucl. Fusion 58 (2018) 126025 D. Hu etal 17 so long as the fragments can still be considered drag-less. Increased injection velocity can both improve the heat load mitigation by prevent premature triggering of the TQ, and shorten the time difference between the core density and temper ature relaxation due to more efficient density convection by the 1/1 kink. Varying the spreading angle or the velocity dispersion has no significant impact. Last, but not least, the property of the target equilibrium is shown to be crucial for achieving better injection efficiency. For instance, higher target plasma temperature results in increased assimilation. Meanwhile, off-axis deuterium SPI into an equilibrium with q(0)>1 shows no strong global density convection, thus poor core penetration. This, however, does not imply off-axis SPI would be ineffective for the advanced scenarios, as their q profile and temperature profile is different from the cases we investigated here. The detailed study for both deuterium and impurity SPI into such a scenario is left for future work. Finally, fast toroidal plasma rotation shows suppression of helical cooling before islands formation, though such stabilizing effect becomes insignificant once the plasma has slowed down. Although the above studies are carried out for a JET plasma, they also contribute to the future ITER DMS design, since the basic MHD processes are expected to remain the same. Disruption mitigation in ITER will rely on the injection of radiating impurities such as neon and argon. Inclusion of such impurity species in JOREK as well as SPI simulation into high temperature H-mode are therefore of high priority and will be done in the near future. Acknowledgments The authors thank L. Baylor, D. Pfefferlé, A. Boozer, B. Breizman, E. Joffrin, J. Snipes and P. Parks for fruitful discussion. ITER is the Nuclear Facility INB no. 174. The views and opinions expressed herein do not necessarily reflect those of the ITER Organization or the European Commission. This publication is provided for scientific purposes only. Its contents should not be considered as commitments from the ITER Organization as a nuclear operator in the frame of the licensing process. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 under grant agreement No 633053. This work is carried out partly on the supercomputer MARCONI operated by Cineca, and also on CURIE operated by CEA. ORCID iDs D.C. van Vugt https://orcid.org/0000-0002-1108-3927 References [1] BaylorL.R. etal 2015 Disruption mitigation system developments and design for ITER Fusion Sci. Technol. 68211–5 [2] LehnenM. etal 2015 Disruptions in ITER and strategies for their control and mitigation J. Nucl. Mater. 46339–48 [3] HenderT.C. etal 2007 ITER progress chapter 3: MHD stability, operational limits and disruptions Nucl. Fusion 47S128–202 [4] NardonE. etal 2017 Progress in understanding disruptions triggered by massive gas injection via 3D non-linear MHD modelling with JOREK Plasma Phys. Control. Fusion 59014006 [5] FilA. etal 2015 Three-dimensional non-linear magnetohydrodynamic modeling of massive gas injection triggered disruptions in JET Plasma Phys. Plasmas 22062509 [6] ParksP.B. and TurnbullR.J. 1978 Effect of transonic flow in the ablation cloud on the lifetime of a solid hydrogen pellet in a plasma Phys. Fluids 211735 [7] GálK. etal 2008 Role of shielding in modelling cryogenic deuterium pellet ablation Nucl. Fusion 48085005 [8] BraginskiiS.I. 1965 Reviews of Plasma Physics vol 1, ed M.A. Leontovich (New York: Consultants Bureau) p 205 (Russian transl.) [9] HoulbergW.A., MiloraS.L. and AttenbergerS.E. 1988 Neutral and plasma shielding model for pellet ablation Nucl. Fusion 28595 [10] LehnenM. etal 2011 Disruption mitigation by massive gas injection in JET Nucl. Fusion 51123010 [11] BussacM.N. etal 1975 Internal kink modes in toroidal plasmas with circular cross sections Phys. Rev. Lett. 351638 [12] WhiteR.B. 2006 The Theory of Toroidally Confined Plasmas 2nd edn (London: Imperial College Press) p 119 [13] StraussH.R. 1986 Hyperresistivity produced by tearing mode turbulence Phys. Fluids 293668 [14] CraddockG.G. 1991 Hyperresistivity due to densely packed tearing mode turbulence Phys. Fluids B 3316 [15] BaylorL. 2017 Developments in shattered pellet technology and implementation on JET and ITER PPPL TSD Workshop Report (Princeton, NJ, USA) (http://tsdw.pppl.gov/ Talks/2017/Lexar/Monday%20Session%201/Baylor.pdf) [16] LehnenM. 2016 ‘SPI pellet sizes for JET’, Memorandum for Project 16008 (internal report) [17] ParksP. 2016 Modeling dynamic fracture of cryogenic pellets GA Report GA-A28352 General Atomics (https://doi. org/10.2172/1344852) [18] HollmannE.M. etal 2007 Observation of q-profile dependence in noble gas injection radiative shutdown times in DIII-D Phys. Plasmas 14012502 [19] ThorntonA.J. etal 2012 Plasma profile evolution during disruption mitigation via massive gas injection on MAST Nucl. Fusion 52063018 [20] PostD.E. 1995 A review of recent developments in atomic processes for divertors and edge plasmas J. Nucl. Mater. 220–2143–57 [21] HuD. and ZakharovL.E. 2015 Quasilinear perturbed equilibria of resistive unstable current carrying plasma J.Plasma Phys. 81515810602 [22] GatesD.A. and Delgato-AparicioL. 2012 Origin of tokamak density limit scaling Phys. Rev. Lett. 108165004 Nucl. Fusion 58 (2018) 126025 D. Hu etal 18 [23] WhiteR.B., GatesD.A. and BrennanD.P. 2015 Thermal island destabilization and Greenwald density limit Phys. Plasma 22022514 [24] FutataniS. etal 2014 Non-linear MHD modelling of ELM triggering by pellet injection in DIII-D and implications for ITER Nucl. Fusion 54073008 [25] SmithH. etal 2006 Runaway electrons and the evolution of the plasma current in tokamak disruptions Phys. Plasmas 13102502 [26] SergeevV.Yu. etal 2006 Studies of the impurity pellet ablation in the high-temperature plasma of magnetic confinement devices Plasmas Phys. Rep. 32 363–77 [27] La HayeR.J. 2006 Neoclassical tearing modes and their control Phys. Plasmas 13055501 [28] MehlhornT.A. 2011 NRL Plasma Formulary (Washington, DC: Naval Research Laboratory) p 34 [29] LitaudonX. etal 2017 Nucl. Fusion 57102001 Nucl. Fusion 58 (2018) 126025