Department of Applied Physics and Science Education

Modelling of shattered pellet injection experiments on the ASDEX Upgrade tokamak

by

Anshkumar Himanshu Patel

MSC THESIS

Assessment committee Member 1 (chair): prof.dr. R.J.E. Jaspers Member 2: dr. T.W. Morgan Member 3: dr. X.C Mi Advisory member 1: dr. G. Papp        Graduation Program: Science and Technology of Nuclear Fusion Capacity group: Science and technology of Nuclear Fusion Supervisor: prof.dr. R.J.E. Jaspers Date of defense: October 26, 2023 Student ID: 1671545 Study load (ECTS): 45 External supervisor(s): dr. G. Papp

The research of this thesis has been carried out in collaboration with Max Planck Institute for Plasma Physics, Garching.
This thesis is public and Open Access.

This thesis has been realized in accordance with the regulations as stated in the TU/e Code of Scientific Conduct.


Disclaimer: the Department of Applied Physics and Science Education of Eindhoven University of
Technology accepts no responsibility for the contents of MSc theses or practical training reports.

Abstract

A disruption mitigation system (DMS) is necessary for fusion-grade tokamaks like ITER in order to ensure the preservation of machine components throughout their designated operational lifespan. To address the intense heat and electromagnetic loads that occur during a disruption, a shattered pellet injection (SPI) system will be employed. This SPI system involves injecting material into the plasma through the means of a cryogenic pellet that is shattered on a bent tube before entering the plasma. The penetration and assimilation (ionized material that stays inside the plasma volume) of the injected material is influenced by various SPI parameters, including the fragment sizes, speeds, and composition of the shattered fragments. An SPI system was installed on the ASDEX Upgrade tokamak to study the effect of the aforementioned parameters. In this thesis, 1.5D simulations with the INDEX code have been utilised to conduct parametric scans, thus examining the influence of fragment sizes, velocities, and pellet composition on the efficacy of disruption mitigation.

When injecting only deuterium, I found material assimilation to be limited to the edge of the plasma with larger and faster fragments leading to higher assimilation. For mixed deuterium/neon injections, again, larger and faster fragments enabled higher assimilation. The amount of assimilated neon increased with increasing injected neon amounts but saturated for larger neon fraction pellets. I also carried out comparisons with previous experimental results of penetration, material assimilation and pre-TQ duration. Previous experimental results for pure deuterium injections indicated that larger and faster fragments exhibit greater penetration, aligning with findings from the simulations. Simulated material assimilation trends for pure deuterium injections were also found to be qualitatively similar to the experiments. Nonetheless, a major difference in the quantitative assimilation values was identified, likely associated with the experimental assimilation criterion. By assuming a semi-empirical TQ onset condition, pre-TQ duration dependence on fragment sizes and speeds for mixed deuterium/neon injections is also studied. I found that faster fragments lead to a shorter pre-TQ duration in the experiments, similar to the simulations. In the case of fragment size variation, larger fragments had a longer pre-TQ duration in the simulations, in contrast to the experiments.

Contents

1. 1 Introduction
2. 2 Theoretical background
   1. 2.1 Disruptions in tokamaks
   2. 2.2 Disruption mitigation
   3. 2.3 Progress on the disruption mitigation system
   4. 2.4 Shattered pellet injection at ASDEX Upgrade
   5. 2.5 Modelling of shattered pellet injections
   6. 2.6 Literature review
   7. 2.7 Objective of this thesis
3. 3 Methodology
   1. 3.1 Setup
      1. 3.1.1 INDEX
      2. 3.1.2 SPI in INDEX
      3. 3.1.3 AUG inputs for INDEX
   2. 3.2 Example: mixed deuterium/neon injection
   3. 3.3 Ablation dynamics for pure deuterium injections
   4. 3.4 Example: pure deuterium injection
4. 4 Simulation results
   1. 4.1 Effect of fragment sizes
      1. 4.1.1 Mixed Ne/D pellets
      2. 4.1.2 Pure deuterium pellets
   2. 4.2 Effect of fragment speeds
      1. 4.2.1 Mixed D/Ne pellets
      2. 4.2.2 Pure deuterium pellets
   3. 4.3 Effect of pellet composition
5. 5 Experimental comparisons
   1. 5.1 Penetration
   2. 5.2 Assimilation
   3. 5.3 Pre-TQ duration
6. 6 Summary and outlook

1 Introduction

To tackle the issue of growing energy demand [1] while mitigating the emission of greenhouse gases, nuclear fusion is a promising solution. The remainder of this section closely follows the discussion from ’Plasma Physics and Fusion Energy’ by Jeffrey P. Freidberg [2]. Fusion reactions involve merging two light nuclei, resulting in a combined mass that is smaller than the original nuclei. The excess mass is converted into usable energy for electricity generation. To obtain net electricity from fusion reactions, the fusing particles must be confined such that fusion power can overcome scattering losses [3]. Various confinement schemes have been developed and explored for fusion reactors. These schemes aim to achieve the extreme conditions of temperature and pressure necessary for initiating and sustaining fusion reactions. A widely researched form of confinement, magnetic confinement, utilizes magnetic fields to confine a mixture of hot charged particles in a plasma state. Magnetic confinement fusion reactors plan to confine a mixture of deuterium (D) and tritium (T).

$$\text{D}+\text{T}\rightarrow^{4}_{2}{}\text{H}_{\text{e}}(3.5\text{MeV})+\text{n}(14.1\text{MeV}).$$ (1)

The D-T mixture will be utilized as fuel due its highest reaction rate [4] among other fuels such as D-D, D-$\text{He}^{3}$. A single D-T fusion reaction emits an alpha particle and a neutron as shown in Equation 1. 80% of the total energy generated in fusion reactions is carried by neutrons while 20% is transferred to alpha particles that would be used to self-heat the plasma to sustain the fusion reaction. The energy of the emitted neutrons can be converted into heat to generate electricity through a turbine. Among various magnetic confinement approaches, tokamaks are the most extensively studied. The tokamak, such as ITER, utilizes magnetic coils arranged around a torus-shaped vessel, which generate a toroidal magnetic field to confine the plasma. A central solenoid induces current in the plasma which generates a secondary magnetic field along the poloidal direction. The two field components result in a helical magnetic field as shown in Figure 1. The helical magnetic field is required to confine the plasma without being affected by grad-B drifts that would otherwise push the plasma towards the plasma facing components. In an axisymmetric tokamak, the magnetic field lines form nested toroidal flux surfaces, where the magnetic flux through the surfaces is constant [5]. The plasma temperature and pressure on the flux surface is also constant as the energy and particle transport along the magnetic field lines is significantly faster than transport between the surfaces. The flux surfaces can also be denoted with a $q$ value which corresponds to the number of toroidal turns a field line on the flux surface makes to complete one poloidal turn.

Figure 1: Magnetic fields in a tokamak. The toroidal magnetic field is generated by toroidal field coils shown by grey color. The central solenoid induces current in the plasma which creates a poloidal magnetic field. The net magnetic field is helical. Image source: EUROfusion

Another magnetic confinement concept is the stellarator, which, unlike the tokamak, uses magnetic coils around the plasma vessel to generate the helical magnetic field directly without externally driven plasma current. All eyes are currently on ITER [6] which is designed to confine a D-T plasma in which $\alpha$ particle heating will dominate all other forms of plasma heating. It is foreseen to obtain inductive plasmas with an energy gain $Q$ = 10 beyond break-even where $Q$ is the ratio of the output fusion power to the input heating power provided to the plasma. It is also aimed to demonstrate the integrated operation of the technologies for a fusion power plant, test the components required as well as the concepts for a tritium breeding module for the tokamak. After that, DEMO [7] will be deployed at which stage nuclear fusion energy will be considered a viable energy source. However, the commercial realization of nuclear fusion power faces various engineering and technological challenges.

One challenge faced by tokamaks is that they are susceptible to disruptions [5], which are unwanted events that occur on a milliseconds timescale leading to a loss of thermal and magnetic confinement of the plasma. These events are a result of strong magnetohydrodynamic instabilities in the plasma. Disruptions in ITER can release substantial thermal and magnetic plasma energy, potentially damaging the machine components and reducing their lifespan [8]. Additionally, electrons with relativistic energies known as runaway electrons (RE) [9], are generated and accelerated in this phase and can also cause severe damage to the plasma-facing components (PFC) and even penetrate to deeper structures. To ensure the longevity of machine components, a disruption mitigation system (DMS) is essential as an investment protection system. In the case of ITER, a shattered pellet injection (SPI) system will be employed [10]. Cryogenic pellets composed of a hydrogen-neon mixture will be launched toward the tokamak and will shatter in a bent tube before entering the plasma. The resulting fragments will ablate, assimilate in the plasma which then radiates its energy uniformly, reduces heat loads on the tokamak components, increases plasma density and suppresses runaway electrons [10].

The penetration and assimilation of the injected material in the hot plasma depends on a few key parameters: the shattered fragment sizes, fragment speeds and the composition of the pellet. Here, the injected material that is deposited and stays within the plasma volume is considered to be assimilated. To study the effect of these parameters, the ASDEX Upgrade (AUG) tokamak recently installed a highly flexible SPI system. In this thesis, I carry out parametric scan simulations of fragment sizes, speeds and pellet composition and study the plasma response using the INDEX disruption simulation code. Furthermore, I also present comparisons of the simulation results with AUG experimental SPI campaign results. The rest of thesis is structured as follows: In chapter 2, the theoretical background and available literature regarding disruptions and SPI is summarized and a research question is defined. Chapter 3 describes the INDEX code, relevant AUG inputs for INDEX and introduces two example SPI simulations in AUG. With this knowledge, simulation results of the parametric scans are presented in chapter 4. Chapter 5 compares the simulation trends of penetration and assimilation with the experimental measurements and chapter 6 summarizes the main results in the thesis.

2 Theoretical background

2.1 Disruptions in tokamaks

The ITER tokamak, like other tokamaks, will be prone to disruptions. These events can result from e.g. a control failure when the magnetohydrodynamic stability limits are exceeded. These instabilities create perturbations to the magnetic field structure. If these perturbations become large enough, they can overlap with each other, causing reconnection of the magnetic field lines and driving significant heat and particle transport across the plasma. During the reconnection event, the plasma current increases slightly (around $10\%$) due to the redistribution of the plasma current profile [11]. The perturbation can lead to an influx of impurities into the plasma that can radiate energy through impurity line radiation. A disruption is usually characterized by two phases, the thermal quench (TQ) and the current Quench (CQ) phase. These different phases of an unmitigated disruption are shown in a schematic diagram in Figure 2. In the TQ phase, the large-scale perturbations drive heat transport across the plasma, cooling it down rapidly (of sub-milliseconds to a few milliseconds) from several keV (kiloelectronvolts) down to a few eVs or tens of eVs. At the end of the TQ phase, the magnetic flux surfaces may start to re-form and cross-field transport losses may be reduced.

Post-TQ plasma temperature is determined by the power balance between ohmic heating by the plasma current and impurity radiation loss [8]. The drop in plasma temperature increases the resistivity of the plasma, causing the plasma current to start to decay. The sudden increase in resistivity also gives rise to a high in-plasma electric field [11]. The decay of the plasma current depends on the resistivity and inductance of the plasma and the surrounding vessel elements. In unmitigated disruptions, the radiated power can be less due to the lower amount of impurities entering the plasma. As a result, the plasma current decays slowly, and the plasma column is displaced vertically following the loss of vertical stability control, resulting in the generation of halo currents [8] and causing high electromagnetic (EM) loads on the vacuum vessel. Alternatively, if the radiated power is too high, the plasma current can decay faster, generating eddy currents in the surrounding structures and intense EM loads between components [8]. The induced electric field is large and can enable electrons to overcome the Coulomb collisional drag force in the plasma, causing them to gain relativistic energies (MeV range), which are then termed runaway electrons (REs) [9]. Depending on the induced electric field and the background plasma parameters, the generated REs can carry a significant fraction of the plasma current. The high energy RE beam can melt plasma facing components and even penetrate to deeper vessel structures.

Figure 2: Schematic diagram of plasma current ($I_{p}$) and plasma thermal energy ($E_{th}$) during a disruption event in ITER and the different phases of disruptions.

Disruptions are already of concern in present-day large tokamaks such as JET and will be more detrimental in ITER because the stored plasma energy will be much higher than that of present tokamaks [8]. For reference, in JET, the thermal energy ($\leq 12$ MJ) and the poloidal magnetic energy ($\sim 10$ MJ) are dissipated in unmitigated TQ and CQ phases in the form of heat to the plasma-facing components (PFCs) on time scales of $\sim 1$ ms and $>20$ ms, respectively [12]. In ITER, a significant fraction of the thermal energy ($\leq 350$ MJ) and the poloidal magnetic energy ($\leq 1300$ MJ) can be deposited on the PFCs [13] during unmitigated disruptions. Taking into account the increase in the wetted divertor area that will receive the heat loads and also the linear increase in the time scale of the TQ energy deposition, the expected peak power density will still be higher in ITER than in JET by a factor of 4 during the TQ [14]. Considering that JET already runs close to the material limits, the concern for PFC lifetime is even higher for ITER due to the shorter deposition timescales and higher asymmetries in deposition than JET.

Another concern is the avalanche multiplication of runaway electrons that occurs when an existing seed of REs generates more REs in the presence of the sufficiently high electric field after the TQ. The avalanche multiplication rate increases roughly by a factor of 10 with an increase in 1 MA of plasma current [9]. Therefore, the multiplication rate for ITER which will have a maximum 15 MA plasma current will be higher than that of JET, where the highest currents are $\sim 4$ MA, by a factor of $10^{11}$. This can lead to the formation of runaway beams that carry up to 10 MA current with a beam energy of 20 MJ [8, 11, 15, 16]. The disruption mitigation system (DMS) at ITER is responsible for addressing threats posed by the different phases of disruption that could affect the operation of tokamak components until the end of their intended lifespan. The basic mechanism of the DMS involves injecting impurities in the plasma to radiate energy uniformly and to increase the plasma density for RE avoidance. To summarize, the main concerns to be addressed by the DMS are the following:

• In the TQ phase: Loss of thermal energy of the plasma on milliseconds timescale that can cause surface melting of PFCs due to high heat loads and avoid generation of runaway electrons after TQ,

• In the CQ phase: Intense electromagnetic forces arising from halo and eddy currents in the vacuum vessel and the blanket module, respectively,

• If a RE beam is formed: Generation of runaway electron beams that can carry a substantial amount of plasma current, which deposits energy on the PFCs that is highly localized and can cause localized melting and even damage to deeper components such as cooling channels [8].

2.2 Disruption mitigation

The DMS in ITER is a final resort option to be utilized when disruption avoidance through plasma control is no longer possible. The requirements for reducing heat and electromagnetic loads to prevent PFC melting, prevent RE formation, and suppress a runaway electron beam, if it is formed in the CQ phase, are described in Lehnen et. al. [8] and Hollmann et. al. [17]. The DMS must reliably dissipate 350 MJ of thermal energy and up to 1 GJ of magnetic energy, while maintaining a CQ duration between 50 and 150 ms to avoid extreme heat and electromagnetic loads. Additionally, the RE current must not exceed 150 kA when it strikes the ITER wall [18].

While there are multiple methods of disruption mitigation [17], the more researched methods used for disruption mitigation in present tokamaks are Massive Gas injection (MGI) and Shattered Pellet Injection (SPI). Both methods involve injecting large amounts of material in the plasma to dissipate thermal and magnetic energy through line radiation, increasing the plasma density to prevent RE formation, and dissipating RE energy by collisional dissipation and synchrotron radiation if a RE beam is formed [8]. The ITER DMS concept relies on SPI technology, which was chosen due to the deeper and faster material assimilation in the plasma core compared to MGI before the TQ, as discussed in Lehnen et. al. [10].

The SPI system involves the formation and acceleration of cryogenic pellets that are shattered into smaller fragments against a tilted surface before entering the plasma. The shattered fragments ablate into neutral gas particles because of the hot background plasma, which then ionize and assimilate within the plasma. Here, assimilation refers to the ablated and ionized particles that are deposited within the plasma volume. The efficiency of the SPI system for disruption mitigation depends on multiple parameters of the SPI system. These parameters include elements like pellet injection schemes (discussed later), optimum fragment sizes and velocities, injection location, pellet compositions and required quantities. These elements are being experimentally evaluated on different tokamaks [19, 20, 21] and are also being studied using simulations [22, 23, 24]. A schematic of the SPI injection system with its various stages is shown in Figure 3. There are two different injection schemes envisioned for the DMS before the TQ [25]: single and staggered injections. A schematic diagram of the plasma parameters during mitigated disruption for the two injection schemes is shown in Figure 4.

Figure 3: Schematic diagram of the shattered pellet injection system and its different stages. $\alpha$ denotes the shattering angle.

(a) Single mixed Ne/H injection.

(b) Staggered injection using pure H injection followed by mixed Ne/H injection.

Figure 4: Schematic diagram of plasma parameters during a mitigated disruption event in ITER and for two different schemes: [Top] Single Ne/H injection, [Bottom] Staggered injection. $I_{p}$ is plasma current, $E_{th}$ is plasma thermal energy, $T_{e}$ is plasma electron temperature and $n_{e}$ is the plasma electron density.

Single injections involve injecting a mixed H/Ne pellet. ITER will utilize hydrogen rather than deuterium, which is utilized in present devices [26, 21] due to the classification of deuterium as a nuclear material [27]. The neon will radiate the plasma energy through line radiation i.e. radiative cooling, while hydrogen assimilates in the plasma to increase the plasma density before the disruption for RE avoidance as shown in 3(a). However, mixed injection of H/Ne pellets has been shown to lead to shorter delay to the TQ that can restrict material assimilation [28]. It is possible that multiple pellets of the same composition will be injected using different injectors, depending on the material assimilation requirements. The second injection scheme, known as the staggered injection scheme, is shown in 3(b). In this scheme, a nearly or completely pure H pellet is injected first that reduces the plasma temperature and increases plasma density by diluting the plasma at constant pressure i.e. dilution cooling, which suppresses RE generation without a radiative collapse [29, 30]. For the second injection, mixed H/Ne pellet(s) can be injected before the TQ to mitigate heat and electromagnetic loads [25]. This scheme has the advantage of allowing more time for the hydrogen to assimilate compared to mixed H/Ne pellets. However, it has to be assessed whether the second mixed pellet can ablate and assimilate in a cooled plasma after the first injection. The second injection then radiates plasma thermal and magnetic energy before the disruption to avoid excessive loads. For both injection schemes, a pre-TQ phase with duration $\Delta\text{t}_{\text{pre-TQ}}$ is defined as the phase from when the first pellet fragments appear in the plasma until the TQ as indicated in Figure 4.

The pellets generated by the ITER SPI system will have a diameter of 28.5 mm, a length-to-diameter ratio of 2 which contains $\sim 2\times 10^{24}$ hydrogen atoms per pellet [31]. There will be 27 injectors at different toroidal and poloidal locations in six different ports. 24 of these injectors in three large equatorial ports will serve the function of mitigating thermal loads, increasing plasma density for RE avoidance, and dissipating RE energy in the pre-TQ phase. The remaining 3 injectors will be housed in three upper ports that are $120^{\circ}$ apart and serve the purpose of mitigation of the CQ in the event that the DMS has not been triggered before the TQ. The upper injectors can inject smaller fragments and even gaseous material compared to the equatorial launchers to ensure sufficient assimilation in the colder post-thermal plasma. The design of the system allows to implement shattering angles in the range of $12.5^{\circ}\text{ to }30^{\circ}$. The pellet velocity will be determined by the desired fragment size and velocity distribution and is expected to be between 200 and 600 m/s [31].

2.3 Progress on the disruption mitigation system

Apart from engineering studies relevant for development, testing and commissioning of the DMS elements, there are various physics relevant studies required to assess the following aspects of the DMS:

1. Adequacy of injection locations,

2. Efficiency of multiple simultaneous injections,

3. Efficacy of runaway electron avoidance and mitigation schemes,

4. Optimum fragment sizes and velocities for highest assimilation,

5. Required quantities, pellet composition, injection sequences.

Experiments are being carried out on present tokamaks with SPI systems which include KSTAR, JET, DIII-D, J-TEXT and ASDEX Upgrade (AUG) [27, 26, 21]. Combined experimental results from these tokamaks and modelling studies will be utilized to assess the aforementioned aspects of the DMS. While ITER has multiple injections at different locations described in the previous section, the injection location(s) can affect the final assimilation of the material and also affect the asymmetries in radiation associated with SPI which arise because of spatially concentrated line radiation from the fragments near the PFCs. For the sake of completeness, the information that is expected to be obtained from different tokamaks is briefly described.

KSTAR’s SPI system has two identical injectors which are 180 degrees toroidally separated [27]. JET’s SPI system has two injectors at the same injection point. DIII-D has two injectors which are 120 degrees toroidally separated [26]. These experiments with different injection locations and number of injectors are being used to address SPI aspects such as the adequacy of injection location to avoid first wall melting due to localised line radiation. The effectiveness of the mitigation process on the difference in the arrival timings of multiple pellets can also be tested. Additionally, owing to the different plasma size and energy content of these devices, an extrapolation to ITER SPI can be made in terms of material assimilation, 0D radiation fraction (fraction of heat dissipated before disruptions by SPI) and 3-D radiation asymmetries. Another important set of parameters of the SPI system are the fragment sizes and speeds. ASDEX Upgrade plays an important role here due to its SPI system being capable of adjusting fragment sizes and speed distributions independently. This is possible due to its unique combination of shattering tubes with different shattering angles [32] (discussed in Section 3). Along with the adequacy of injection locations and the number of injections, optimal fragment sizes and speeds must be determined for maximum penetration, assimilation, radiation fraction and distribution of radiation load.

Apart from fragment sizes and speeds, the composition of the pellet also affects the penetration and assimilation of the injected material. The pellets at ITER will be formed by a mixture of Z=1 hydrogen (H) and Z=10 neon (Ne) depending on the injection scheme [33]. It is well known from pellet fuelling experiments that pellets injected from the low field side (LFS) of the toroidal magnetic field, as will be the case in ITER, experience a drift that causes a shift of the deposition profile relative to the ablation profile towards the LFS [34, 35]. This drift is known as the plasmoid drift. Ideally, pellet assimilation should happen in the core of the plasma to radiate maximum energy, however, the plasmoid drift can prevent this from happening [36]. It was recently shown [37] that adding a small fraction of neon to the pellet reduces outward drift and facilitates deeper deposition than pure H pellets. The drift effect is minimized due to the line radiation from neon and is discussed in further detail in Section 3.3. However, larger fractions of neon can also disrupt the plasma before the material can assimilate and also cause a very fast CQ leading to intolerable electromagnetic loads [23]. The shorter pre-TQ duration for a mixed neon injection is due to the formation of a radiative cold front in the edge plasma by the injected impurities that destabilise tearing modes at low-order rational $q$ surfaces. When the radiative cold front reaches the q=2 surface, the modes can grow near the edge, causing enhanced radial transport and the collapse of the hot core plasma, similar to MGI-triggered disruptions [38, 39, 40, 41, 42, 43].

2.4 Shattered pellet injection at ASDEX Upgrade

ASDEX Upgrade (AUG) is a medium-sized tokamak with a divertor configuration and a full tungsten wall located at the Max Planck Institute for Plasma Physics, Garching, Germany [44]. AUG installed a highly flexible SPI system in late 2021 to provide inputs for the design and optimisation of the ITER DMS [32]. The AUG SPI system is a triple-barrel system where each independent guide tube can be equipped with different shatter heads. The pellets in AUG were made with pure deuterium, pure neon and mixed deuterium/neon injections with varying compositions. After extensive laboratory testing [45], the three shatter heads (shown in Figure 7) selected included a 46 mm long, circular cross-section with shattering angle $25^{\circ}$ which generate a fragment plume with an increased spatial spread. The remaining two shatter heads had a similar 78 mm long, rectangular cross section but with two different shattering angles $12.5^{\circ}$ and $25^{\circ}$ respectively. The fragment size distribution was studied to be strongly dependent on the impact velocity $v_{\perp}$ of the pellet with respect to the shatter head where $v_{\perp}=v_{\text{pellet}}\cdot\sin\alpha$. Lab experiments[45] showed that the mean fragment size $\langle d_{\text{frag}}\rangle$ decreases for increasing normal velocity. The mean fragment velocity was found to be close to the pellet velocity and initial analysis shows that it lies between the pellet velocity and the parallel component of pellet velocity with respect to the shattering angle i.e. $\langle v_{\text{frag}}\rangle\in[v_{\parallel},v_{\text{pellet}}]$ where $v_{\parallel}=v_{\text{pellet}}\cdot\cos\alpha$. Injecting pellets with similar pre-shatter speed, one can generate varying fragment size distributions by using shatter heads with different shattering angles. Similarly, by comparing injections with similar impact velocity from different shattering angles, different fragment velocity distributions can be created. An experimental campaign was conducted in 2022 where $\sim$ 240 tokamak discharges were allocated to the SPI experiments. Some of the primary machine, plasma, and SPI parameters are shown in Table 1.

Table 1: Machine and plasma parameters of AUG for discharges in this report.

  Plasma major radius $R_{0}$ 1.65 m Plasma minor radius $a$ 0.5 m Toroidal magnetic field strength $B_{T}$ $1.8$ T Plasma current $I_{p}$ $\leq 800$ kA Heating power $P_{\text{total}}$ 27 MW H-mode plasma energy content $W_{th}$ $500\text{~{}kJ}<W_{th}<800$ kJ Pellet length $L_{\text{pellet}}$ 4.5 mm $\leq$ $L_{\text{pellet}}\leq$ 10.5 mm Pellet diameter $d_{\text{pellet}}$ 4 mm / 8 mm

2.5 Modelling of shattered pellet injections

In support of the experimental activities carried out at various tokamaks, additional simulation activities are also required for the physics-based extrapolation from the experimental results to ITER. The primary focus of the simulation activities include three-dimensional (3D) MHD simulations performed using codes JOREK [46], M3D-C1 and NIMROD [24] to understand the plasma dynamics during mitigated disruptions and development of models to improve understanding of generation and mitigation of runaway electrons. These activities are ongoing and have been successful in interpreting present experimental results and extrapolating to future tokamaks [47, 24, 48]. JOREK simulations for AUG are also being carried out presently [49]. However, as 3D modelling is computationally expensive, only a limited number of simulations to explore the design parameters can be carried out before the design completion of the DMS. On the opposite end of the spectrum in terms of simplicity, zero dimensional (0D) modelling using particle and energy balance can be used to assess key trends of material assimilation and radiation loss [50]. While the simplicity allows for assessing key trends of the SPI assimilation, it does not take into account the plasma size information which is relevant for assessing material penetration, radial distribution of assimilated material and extrapolating to other present devices and larger devices such as ITER. To complement the aforementioned modelling approaches, SPI simulations based on one-dimensional (1D) transport simulations can also provide inputs for parameter spaces that can be studied in more detail with more complex simulations. The simplified modelling permits extensive parametric scans taking into account plasma size and SPI configuration information. This thesis focuses on such simulations using the disruption simulation code INDEX [23] to study and match the experimental trends of penetration, and assimilation of the injected material for different SPI parameters, these being different fragment sizes, speeds, and pellet composition.

2.6 Literature review

Previous modelling and experimental activities provide some indications about the effect of fragment sizes, speeds and pellet composition on the efficiency of the DMS which are discussed in this section. Previous modelling activities reported the following:

• ITER simulations with the INDEX code have shown that larger fragment sizes, higher velocity, and higher velocity dispersion lead to deeper penetration and higher assimilation of the material [23]. Additionally, benchmarking of INDEX with axisymmetric JOREK simulations for ITER SPI scenario was carried out and agreement between total ablation rate, radiated power and plasma profile evolution was found when a similar SPI setup (of fragment sizes, speeds, and injection timing) is used.

• Simulations with the 3D MHD code NIMROD indicate that a mixed D/Ne injection might be more benign than pure deuterium injection in terms of conducted heat loads on the plasma facing components after pellet injection [24].

• 3D JOREK simulations [46] studied the MHD response for ITER SPI. An important ’size-effect’ was pointed out which indicates that even for similar plasma thermal energy density, electron temperature, and ratio of injected atoms to plasma volume, devices with larger major radius and stronger magnetic field would cool down slower after SPI.

Analysis of the SPI experiments at AUG is currently ongoing, which includes studies on the effect of fragment size, speed and composition on pre-TQ duration, material assimilation, toroidal radiation asymmetry after injection, and radiated energy fraction [51, 52, 53]. Further studies of the laboratory shattering experiments [45, 54] are also being carried out to determine fragment size and speed distribution for different pellet velocities and shattering angles and also different shatter heads. The results of the experimental data analysis relevant for this thesis is discussed in Section 5. Additional modelling activities involving JOREK [49] and DREAM [55] simulations for AUG SPI discharges are also being carried out.

2.7 Objective of this thesis

The research question that this report aims to answer is the following: To what extent can the experimental trends observed in penetration, assimilation of shattered fragments in AUG be corroborated with disruption simulations?

This thesis focuses on validating the relationship between injected material penetration and assimilation in the ASDEX Upgrade (AUG) tokamak with three key SPI parameters: fragment sizes, speeds, and pellet composition. The simulations carried out in the thesis will focus on the pre-TQ phase when the majority of material penetration and assimilation takes place [23]. In this thesis, assimilation is linked to the ablated and ionized material that remains within the plasma volume.

The INDEX code presents a viable approach for conducting extensive parametric scans of SPI parameters, taking into account material penetration and plasma size information, while remaining computationally efficient compared to 3D codes. As ITER simulations with INDEX progress [23], it is crucial to validate the SPI models within INDEX using existing devices.

By successfully comparing INDEX simulations with experimental trends, the simulation results can be used to identify ideal pellet and fragment parameters for deepest penetration, highest and uniform assimilation. Successful matching can be followed by identifying relevant quantities that should be used to extrapolate to ITER SPI in future work. Unsatisfactory matching with experimental results can also provide insights on missing aspects of experiments not being included in the simulations which can significantly affect plasma dynamics that require new or improved models.

3 Methodology

3.1 Setup

3.1.1 INDEX

To simulate mitigated disruptions, I have used the INDEX (Integrated Numerical Disruption EXperiment) [37] code. INDEX assumes that the plasma temperature and pressure depend on only a single spatial coordinate, the toroidal magnetic flux ($\rho$) which is a proxy for the plasma radius. The plasma properties are assumed to be axisymmetric and hence INDEX is referred to as a 1.5D code.

The INDEX code solves a magnetic diffusion equation on 1D flux coordinates and Grad-Shafranov (G-S) equilibrium on the 2D space with the boundary conditions given by solving a series of circuit equations of toroidally continuous filaments [56].The magnetic diffusion equation (refer to Equation 6) describes the time evolution of the magnetic fields in the plasma while the Grad-Shafranov equation describes the magnetic field and pressure configuration in the poloidal plane of the tokamak. The plasma equilibrium is evolved on a resistive timescale. This timescale corresponds to the temporal evolution of the equilibrium on the order of the time associated with the diffusion and dissipation of magnetic fields due to finite plasma resistivity. In essence, the G-S equation indirectly accounts for the ideal Magnetohydrodynamics (MHD) effects required for plasma stability.

Hence, the INDEX code can be used to simulate the pre-TQ phase all the way to towards the end of the CQ phase, usually till fast vertical displacement events [57]. The transport module models the plasma ion species, and additional impurity ion species for different charge states based on OpenADAS ionisation and recombination rates. The particle and energy transport of the ion species $\alpha$ with the $i$th charge state and energy transport for electrons that are used are:

Particle balance

$$\frac{1}{V^{\prime}}\frac{\partial}{\partial t}\left(n_{\alpha}^{i+}V^{\prime}\right)=\frac{1}{V^{\prime}}\frac{\partial}{\partial\rho}V^{\prime}\left\langle|\nabla\rho|^{2}\right\rangle D_{\alpha}\frac{\partial n_{\alpha}^{i+}}{\partial\rho}+S_{\alpha}^{i+},$$ (2)

Ion energy balance

$$\frac{3}{2}\left(V^{\prime}\right)^{-5/3}\frac{\partial}{\partial t}\left(p_{\alpha}V^{\prime 5/3}\right)=\frac{1}{V^{\prime}}\frac{\partial}{\partial\rho}V^{\prime}\left\langle|\nabla\rho|^{2}\right\rangle\times\left[\kappa_{\alpha}\frac{\partial p_{\alpha}}{\partial\rho}+T_{\alpha}\left(\frac{3}{2}D_{\alpha}-\kappa_{\alpha}\right)\frac{\partial n_{\alpha}}{\partial\rho}\right]\quad+P_{\mathrm{ex}}^{\alpha e}+\sum_{\beta}P_{\mathrm{ex}}^{\alpha\beta},$$ (3)

Electron energy balance

$$\frac{3}{2}\left(V^{\prime}\right)^{-5/3}\frac{\partial}{\partial t}\left(p_{e}V^{\prime 5/3}\right)\\ =\frac{1}{V^{\prime}}\frac{\partial}{\partial\rho}V^{\prime}\left\langle|\nabla\rho|^{2}\right\rangle\times\left[\kappa_{e}\frac{\partial p_{e}}{\partial\rho}+T_{e}\left(\frac{3}{2}D_{e}-\kappa_{e}\right)\frac{\partial n_{e}}{\partial\rho}\right]\\ \quad+P_{\text{ohm }}-P_{\text{rad }}-P_{\text{ion }}+\sum_{\alpha}P_{\mathrm{ex}}^{e\alpha},$$ (4)

where $\rho$ is the magnetic surface label, $V(\rho)$ denotes the plasma volume enclosed by the magnetic surface, $V^{\prime}=\frac{dV}{d\rho}$, and the bracket $\langle\bullet\rangle$ denotes the surface averaging process. $n_{\alpha}^{i+}$ is the the density of the ion species $\alpha$ charge state with charge state $i$. $p_{e}$ is the electron pressure. A single scalar pressure for all the ion species is considered $p_{\alpha}=\left(\Sigma_{i=1}^{Z}n_{\alpha}^{i+}\right)T_{\alpha}$ is considered where $T_{\alpha}$ is the temperature of each ion species $\alpha$. In INDEX, quasineutrality is assumed and hence the electron density is calculated at each time step from the ion density. $P_{\text{ohm}}$ is the ohmic heating power and $P_{ex}^{\alpha\beta}$ is the energy transfer between different ion species $\alpha$ and $\beta$. The ionisation and recombination rates used to calculate the particle source and sink terms $S_{\alpha}^{i+}$ and the radiative cooling rates used to calculate the radiation and ionisation losses, $P_{\text{rad}}$ and $P_{\text{ion}}$, are extracted from the OpenADAS database [58]. $D_{\alpha}$ and $\kappa_{\alpha}$ are the particle diffusion and heat diffusivity coefficients. In addition to the 1D transport equations above, a neutral particle source equation is also solved:

$$\frac{1}{V^{\prime}}\frac{\partial}{\partial t}\left(n_{\alpha}^{0}V^{\prime}\right)=\frac{1}{V^{\prime}}\frac{\partial}{\partial\rho}V^{\prime}\left\langle|\nabla\rho|^{2}\right\rangle D_{\alpha}^{0}\frac{\partial n_{\alpha}^{0}}{\partial\rho}+S_{\alpha}^{0}+S_{\alpha}^{\mathrm{SPI}}$$ (5)

where $S_{\alpha}^{0}$ and $S_{\alpha}^{SPI}$ are the neutral particle source/sink terms associated with ionisation/recombination and the particle source due to the SPI respectively. Similar to recent pre-TQ SPI simulations with INDEX [23], neutral particle sources introduced in the pre-disruption hot plasma ionize instantaneously and do not influence the plasma significantly in the pre-TQ phase. Hence, the neutral particle transport $D_{\alpha}^{0}$ is set to be equal to the ion diffusion coefficient of the same species. The 1D coordinate $\rho$ is related to the toroidal magnetic flux $\chi=2\pi\oint\mathbf{B}\cdot d\mathbf{S}_{\varphi}$ as $\rho=\sqrt{\frac{\chi}{\chi_{s}}}$ where the subscript s denotes the plasma surface and $\phi$ is the toroidal angle. The current density profile is coupled with the plasma density and temperature through the magnetic diffusion equation which is expressed in terms of the rotational transform $\nu=\frac{d\psi}{d\chi}$ as

$$\left.\frac{\partial\nu}{\partial t}\right|_{\chi}=\frac{\partial}{\partial\chi}\left[\frac{\eta}{\mu_{0}}\frac{(d\chi/dV)^{2}}{A^{2}}\frac{\partial}{\partial\chi}(AK\nu)-\frac{\eta}{d\chi/dV}\left\langle\mathbf{j}_{\mathrm{ex}}\cdot\mathbf{B}\right\rangle\right]$$ (6)

where $K(\psi)=\left<\frac{|\nabla V|^{2}}{R^{2}}\right>$ and $A(\psi)=\left<\frac{1}{R^{2}}\right>$ are the surface averaged metric coefficients and $\eta(\psi)$ denotes neoclassical conductivity [59]. The last term in Equation 6 is used to model externally or runaway driven current densities. While runaway electrons can play a significant role in the CQ phase, they are inconsequential in the pre-TQ phase and hence are not considered in the simulations carried out in this thesis. By coupling Equation 2, Equation 3, Equation 4, Equation 5 with Equation 6 through $\eta(\psi)$, the evolution of the current density profile can be self-consistently solved taking into account changes in the plasma temperature and density due to the assimilation of the injected material. The INDEX free-boundary equilibrium solver is similar to that of the DINA code [60, 61].

3.1.2 SPI in INDEX

The INDEX code has a SPI module which has been applied to various tokamaks such as DIII-D [62], KSTAR, ITER [23] and now AUG. To generate a fragment plume, first INDEX generates an initial fragment size distribution ($r_{p}(t=0)$) calculated using a Bessel function based on a statistical fragmentation model [63] with the probability distribution:

$$P(r_{p})=\frac{r_{p}K_{0}(\kappa_{p}r_{p})}{I},\quad I=\int_{0}^{\infty}r_{p}K_{0}(\kappa_{p}r_{p})dr=\kappa_{p}^{-2}$$ (7)

where $K_{0}$ is the modified Bessel function of the second kind, $\kappa_{p}$ is the inverse of the characteristic fragment size and is related to the total number of molecules in the pre-shattered pellet and to the number of fragments generated after shattering [47]. INDEX pseudo-randomly samples $N_{f}$ number of different fragment sizes $r_{i}(i=1,2...N_{f})$ from the probability distribution in Equation 7. Assuming each fragment has a spherical shape, it then calculates the number of molecules in each fragment $N_{i}=\frac{4}{3}n_{m}\pi r_{i}^{3}$ where $n_{m}$ is the molecular density of the pellet material. Since only a finite number of fragments are sampled, the sum of the molecules in the fragments do not match with $N_{p}$. To compensate for this, the fragment size distribution is renormalized to match the pellet material amount resulting in the number of molecules in the new fragments being $N_{i}^{\prime}*\frac{N_{p}}{\Sigma_{i}N_{i}}N_{i}$ and the new sizes being $r_{i}^{\prime}=(3N_{i}^{\prime}/4\pi n_{m})^{1/3}$. Since the fragments sizes are sampled using pseudo-random numbers, different realisations of the same probability distribution can be generated. Experimentally, the pellet shattering process is also a stochastic one so the shattered fragments will have slightly different sizes even if the pellet injection parameters are the same. The effect of different realisations of fragment size distributions from the same probability distribution is only briefly discussed in this thesis and is being further investigated in DREAM simulations [55]. However, to ensure that trends in penetration and assimilation are due to the changes in the shattering parameters and not due to a particular outlier realization of the fragment size and speed probability distribution, 3-5 different realizations are used for simulation sets of parametric scans. However, to systematically eliminate the effect of outlier realizations, substantially more realizations should be utilized.

The initial velocity distribution of the fragments is modelled as a Gaussian distribution with a user-defined standard deviation and the pellet speed as the mean. The standard deviation of the fragment speed distribution for all the simulations in this thesis was set to 40% based on JET SPI experiments [64]. It is expected that different pellet speeds and shattering angles might modify the standard deviation [45], however additional experimental inputs are required for taking this into account. After the generation of fragments with corresponding sizes and speeds, the fragments travel towards the plasma with a toroidal and poloidal spread which was set to $\pm 20^{\circ}$ for all the simulations in this thesis [45]. INDEX maps the 3D location of the fragments at each time step to the 1D flux coordinates used in the particle and energy transport equations taking into account the toroidal geometry. The ablation of shattered fragments in the plasma is modelled using the Neutral Gas Shielding (NGS) model which has been developed and used previously for fuelling pellets [65]. Experimentally, the dependence of the penetration depth of non-shattered pellets for fuelling the core plasma [66, 67] follows a scaling close to that obtained from the NGS model. The recession rate of the fragment radii used in INDEX has been proposed in Parks et. al. [68] and can be calculated for neon-deuterium mixed pellets with molar ratio $X$ of Deuterium as

$$\frac{dN_{\text{mix }}}{dt} =\frac{C\lambda(X)}{f_{W}(1-X)+X}n_{e}^{1/3}T_{e}^{5/3}r_{p}^{4/3}\text{ particles s }^{-1},$$ (8)

$$\lambda(X) =A+\tan(BX)$$

where $f_{W}$ is the mass ratio between the neon atom and the deuterium molecule, $n_{e}$ and $T_{e}$ are the electron density (in $m^{-3}$) and electron temperature (in eV), respectively at the location of the fragment with radius $r_{p}$ (in m). Using Equation 8, the number of ablated deuterium and neon atoms are deposited in the plasma are estimated using:

$$\frac{dN_{D}}{dt}=2X\frac{dN_{mix}}{dt};\quad\frac{dN_{Ne}}{dt}=(1-X)\frac{dN_{mix}}{dt}.$$ (9)

The numerical coefficients in Equation 8 are [23, 68] $f_{W}\equiv W_{Ne}/W_{D_{2}}=20.183/4.0282$, $A=27.0837$, $B=1.48709$, $C=4.062\times 10^{14}$. During a simulation, INDEX traces the motion of multiple pellet markers with given fragment velocities, calculates the recession in fragment sizes using the NGS model and assumes the ablated particles are deposited as a flux-surface-averaged neutral particle source. The particle source is convoluted into 1D radial flux coordinates using a Gaussian shape centred at the particle source location and with a width of a few centimetres in real space based on previous measurements of fuelling pellet ablation [69].

3.1.3 AUG inputs for INDEX

Before running INDEX simulations for AUG SPI, I set up various input parameters and profiles. These inputs include geometric and electrical properties of various ohmic, poloidal field coils, and vessel elements that influence the plasma equilibrium convergence in INDEX. A set of initial conditions of electron, temperature and current density radial profiles are also required as input. A corresponding plasma equilibrium for the above inputs is also required. Coordinates of plasma facing components (PFCs) act as boundaries which when crossed by the plasma cause the simulation to stop. For the SPI, coordinates of the beginning and the end of the shattering tubes is also required. I implemented the geometrical and electrical properties of various AUG coil and vessel elements in the INDEX code, which are shown in Figure 5. The rectangular elements outside the vacuum vessel in Figure 5 are the ohmic and the poloidal field coils. The coils with prefix ”OH” refer to ohmic heating coils, prefix ”V” are shaping coils, prefix ”CoI” are position control coils and prefix ”PSL” are passive stabilising loop coils [70]. The plasma equilibrium surfaces of constant poloidal flux are shown in colours ranging from blue to yellow. Blue scatter points on the last closed flux surface show a few of the points provided as input for the plasma equilibrium in INDEX. Black scatter points on the plasma facing components are the boundary points provided to terminate the simulation if crossed. The orange arrow line shows the injection vector for one of the $25^{\circ}$ injection tubes in AUG.

Figure 5: AUG geometrical components implemented in the INDEX code shown in a poloidal cross section.

The 1D profiles of electron temperature, density and current density used for initial equilibrium convergence in INDEX are shown in Figure 6. These profiles were taken from #40355, which was a typical H-mode discharge, at 2.25s. The same profiles are used as inputs for JOREK simulations of AUG SPI [49].

(a) Electron temperature (left axis) and electron density (right axis) profiles

(b) Current density (left axis) and resultant safety factor q (right axis)

Figure 6: Radial plasma profiles used for initial convergence of plasma equilibrium in INDEX.

As mentioned before, AUG has 3 different shatter tubes whose start and end coordinates (X,Y,Z in machine coordinates) are provided in Table 2. An image of the shattering tubes inside the machine is also shown in Figure 7. For all the simulations carried out in this thesis, the particle diffusion coefficient was set to $D_{\alpha}=2$ m^2/s and the thermal diffusivity for ions and electrons was set to $\kappa_{\alpha}=4.5$ m^2/s. It should be noted that the diffusion coefficients were assumed to be radially uniform. In this thesis, the simulations were conducted without considering any inherent machine impurities (like tungsten) to prioritize the exploration of parametric variations in the SPI parameters.

Figure 7: The three shatter tubes inside the vessel. Figure taken from Heinrich et. al. [71].

Table 2: Start and end coordinates of the shattering tubes at AUG (in machine XYZ coordinates with the geometric centre of the machine as the origin)

Shattering tube (shattering angle, cross-section shape)	Shattering tube start coordinate (X,Y,Z) m	Shattering tube end coordinate (X,Y,Z) m
GT1 $(25^{\circ}$, rect.)	(2.340, -0.510, 0.347)	(2.227, -0.489, 0.308)
GT2 $(25^{\circ}$, circ.)	(2.233, -0.458, 0.343)	(2.295, -0.452, 0.324)
GT3 $(12.5^{\circ}$, rect.)	(2.341, -0.490, 0.385)	(2.271, -0.469, 0.361)

For comparison of assimilation measurements, Thomson scattering [72, 73] and CO₂ laser interferometer [74, 75] diagnostics were used in AUG. Initial assessment of the validity of the diagnostic signals and their limitations during SPI was reported in my ITER internship report [51]. To compare the INDEX simulations to the experiments, I implemented synthetic diagnostics of these two diagnostics in the INDEX code. The line of sights of these diagnostics along with some other relevant diagnostics are shown in the poloidal cross section in Figure 8. The COO interferometers have a vertical line of sight shown by the black solid lines. The Thomson scattering diagnostic has the same poloidal location as the core channel of the interferometer.

Figure 8: Relevant line of sight diagnostics for SPI in the AUG poloidal cross section. Black dashed arrow shows the injection vector for the GT1 shatter tube.

3.2 Example: mixed deuterium/neon injection

In order to introduce the simulation procedure, an example pre-TQ simulation with an injection of 10% neon / 90% deuterium (molecular number composition) is discussed. It starts at $t_{\text{sim}}=0$ ms when the fragments are assumed to be at the shattering point. In the present simulation, 200 fragments were sampled for a pellet of length of 10 mm and diameter of 8mm. The fragment velocities were sampled from a Gaussian distribution centred at $\langle v_{\text{frag}}\rangle=230$ m/s with $\Delta v/v=40\%$. The resultant fragment size and speed distributions are shown in Figure 9. The average size from the fragment size distribution was 1.14 mm. In the remainder of this thesis, fragment size and speed distributions for multiple distributions and realizations will be shown as continuous distributions for plotting purposes (example: Figure 18) however it should be kept in mind that the distributions have been discretized.

Figure 9: Binned fragment size and speed distributions for the example simulation.

The time trace of plasma current evolution is shown in Figure 11. An initial decay in the plasma current is caused by two reasons. First, the current profile is modified as the simulation progresses due to the absence of a bootstrap current model in INDEX which plays a significant role in H-mode discharges. Second, while the loop voltage (which drives the current) is taken into account in the circuit equations (refer to Section 3.1.1), the present simulation does not include ”realistic” AUG coil currents actions to maintain constant $I_{p}$. To compare the pre-TQ duration in the experiments to the simulation, a condition similar to the experiment was considered. Experimentally, the arrival of the fragments in the plasma is expected when a significant rise in an edge AXUV channel (Figure 8) is observed. Similarly, in the simulation, the arrival of the fragments in the plasma is considered when the first fragment crosses a circle centred at the plasma magnetic axis and a radius equal to the distance between the plasma centre and the point of intersection of the edge AXUV channel and the injection vector, which is $\sim$ 50 cm. The time of fragments entering the plasma is marked in Figure 11 at 0.8 ms with a black dashed line. Several vertical lines are marked as well which correspond to time stamps whose plasma profiles are plotted in Figure 11. The radial profiles are plotted against the normalized poloidal flux label $\psi_{N}=(\psi-\psi_{a})/(\psi_{s}-\psi_{a})$ where $s$ and $a$ denote the magnetic axis and the last-closed flux surface, respectively. Hence $\psi(0)$ corresponds to the magnetic axis while $\psi(1)$ corresponds to the last closed flux surface. The choice of this label was motivated by previous INDEX simulations [23] that were used for benchmarking with JOREK results.

As discussed in Section 2.2, the pre-TQ duration is an important parameter that can limit the amount of material assimilation in the plasma before the onset of the disruption. Based on the discussion in Section 2.3, the TQ onset is expected when the plasma temperature becomes less than 10 eV inside the q=2 surface. Such a semi-empirical TQ onset condition is required as the 1D transport simulations do not include any description of linear or nonlinear MHD instabilities. However, it is expected that the axisymmetric current contraction that can lead to the TQ can be captured by simulating the particle energy balance and the resulting current density profiles. The onset of the TQ condition in this example is marked in Figure 11 at $\sim$ 2.15 ms.

Figure 10: Time trace of plasma current evolution for a mixed deuterium/neon injection.

Figure 10: Time trace of plasma current evolution for a mixed deuterium/neon injection.

Figure 11: Profiles for various SPI relevant parameters against normalized poloidal flux: (A) electron temperature, (B) electron density, (C) neon density, (D) neutral source function, (E) radiated power, (F) Current density. Different colors indicate different time points of the simulation.

In Figure 11, the plasma response to the SPI is shown. The physical events are given in chronological order:

1. As the fragments penetrate the plasma, initial rise in the neutral particle source function at the edge is observed in Figure 11 (D) in the purple profile corresponding to 1.2 ms in the simulation. An increase in the radiated power is also observed in Figure 11 (E) due to line radiation associated with neon.

2. At 1.8 ms, corresponding to red colored profiles, a significant neutral particle source is deposited towards the core as shown in Figure 11 (D). Following behind the neutral particle source in space, an increase in the electron and neon density can be observed in Figure 11 (B), (C) respectively. A further rise in intensity in the radiated power can be observed in Figure 11 (E) which only lags slightly behind the ablation profile in space. It should be noted that I have observed similar movement of the ablation front and radiation front in space in preliminary analysis of experimental SPI discharges [51]. Because of the radiated power, a drop in the plasma temperature can also be observed in Figure 11 (A). However, the formation of a cold front ($<$ 10 eV) does not occur yet.

3. At 2.2 ms (dark yellow color profiles), the penetration of material continues which is ablated in the still sufficiently hot plasma. This can be observed in the ablation profile in Figure 11 (D). More material is deposited which can be observed in Figure 11 (B,C). Further radiative cooling of the plasma leads to the formation of a cold front towards the edge of the plasma. A small bump can be observed in the radiated power profile in Figure 11 (E) which occurs at the forefront of the cold front. The small bump in the current density profile in Figure 11 (F), arises due to the contraction of the current profile from the edge due to the movement of the cold front. As the TQ onset condition is satisfied at $\sim$ 2.15 ms, strong MHD activity can be expected at this time which would enhance transport significantly. As discussed in Section 2, the onset of the TQ also leads to a small spike in the plasma current (also seen in experiments, refer Figure 13) which is not observed in the present simulations. This discrepancy arises as the present simulations do not utilize a hyper-resistivity model [76] that can emulate the enhanced plasma transport.

With the fragments entering at 0.8 ms, the pre-TQ duration for this example simulation is 1.35 ms. This method of estimating the pre-TQ duration is used later in this thesis for comparing trends of simulated pre-TQ duration(s) with experimental pre-TQ measurements. The simulation stops at $\sim$ 2.6 ms due to difficulties in plasma equilibrium convergence.

As discussed above, a common feature of the mixed deuterium/neon injections is that the increase in ablation and radiation follows the movement of the fastest fragments however, the cold front lags behind in space. While this can be seen in the plasma profiles in Figure 11 (A,D), this distinction can be made more clear by studying the electron temperature and density profiles in space and time as shown in Figure 12. One of the black dotted lines represent the trajectory of one of the fastest fragments in the distribution which appears at the plasma edge at $\sim$ 1.1 ms and the other black dotted line appearing at $\sim$ 1.8 ms shows the trajectory of the fragment whose speed was at the 25^th percentile within the distribution indirectly indicating the position of the slowest few fragments. While the plasma density starts increasing immediately after the entry of the fastest fragments, plasma cooling happens on a longer timescale. The edge and intermediate plasma regions have a temperature on the order of 100 eV after the first fragments and eventually the formation of a cold front occurs at $\sim$ 2 ms as marked by the black dashed line. When the cold front reaches the $q=2$ surface which corresponds to $\psi_{N}\approx 0.86$, the onset of the TQ is expected.

(a) Electron temperature space and time evolution

(b) Electron density space and time evolution

Figure 12: Density and temperature profile evolution in space and time. Black dotted lines mark the trajectories of the fastest and slowest fragments. Black dashed line marks the movement of the cold front. Vertical black dotted line marks the crossing of the cold front beyond the $q=2$ surface hereby satisfying the TQ onset condition.

3.3 Ablation dynamics for pure deuterium injections

Pure deuterium pellet injections experience a plasmoid drift effect when injected from the magnetic LFS [37]. In this section, I will briefly describe the plasmoid drift, which was first observed with fueling pellets [77]. The injected fragments start to ablate when they enter the plasma which leads to the formation of a cigar shaped plasmoid along the field lines. The plasmoid expansion rate is slower than the rate of energy input into the fragments leading to over-pressurised plasmoids. The plasmoids become vertically polarised due to charge separation in a non-uniform magnetic field and experience a $\vec{E}\times\vec{B}$ drift [78] in the direction of the major radius. This drift can reduce the assimilation (ionized material within the plasma volume) of the injected material, even limiting the assimilation to the edge of the plasma [37]. This effect is undesirable as increasing core density is important to avoid runaway electron formation. This drift is not observed for mixed D/Ne injections. Modelling [79] results indicate that the line radiation from neon for a mixed neon pellet has a significant impact on limiting the temperature and pressure of the ionized plasmoid. As a result, the ionized plasmoid is expected to homogenize along the magnetic field lines where the ablated material is released, without causing substantial cross-field drift motion.

In INDEX, the drift estimation is not based on first principles but instead on a back-averaging model [80]. A user-defined fraction of the ablated atoms at radius $\rho^{\prime}=r/a$ is deposited uniformly over the exterior plasma volume where $\rho>\rho^{\prime}$. The back-averaged density increase caused by single shards is then calculated as

$$\Delta n_{\text{dep}}=\Delta N/\left[(1+\beta)V_{p}(1)-V_{p}(\rho^{\prime})\right]$$ (10)

where $\Delta N$ is the number of ablated atoms and $V_{p}(\rho^{\prime})$ is the interior plasma volume with $\rho<\rho^{\prime}$. The extent of the radial drift displacement is determined by the factor $\beta$. Physically, $\beta$ corresponds to the ratio of the volume outside the plasma where the particles can be deposited to the volume inside the plasma. Hence, it controls the spatial extent of the particle deposition. Due to the lack of a self-consistent drift model in INDEX, an interpretative approach is required to determine the value of $\beta$. In the present setup, the same $\beta$ value is used for all the fragments. A user-defined fraction of shifted atoms can be set to any value between 0 and 1 where 0 corresponds to all the deposited atoms being shifted from the ablation location (the extent of the shift depending on $\beta$) and 1 corresponds to the atoms ablating at the fragment location. For pure deuterium injections presented later in this report, the fraction of locally deposited particles was set to 0.

Figure 13: Time traces of plasma and diagnostic signals for #40743: (a) Plasma current in MA, (b) AXUV bolometer signal whose LOS passes very close to the injection location, (c) Soft X-Ray (SXR) amplitude from a channel passing through the core of the plasma. Black dashed line marks the time stamp of Thomson scattering measurement.

I carried out a set of simulations for this thesis with values of $\beta=2,3,4$ (refer to Equation 10) to match to the Thomson scattering measurements of #40743. Plasma and diagnostic parmeters of #40743 are shown in Figure 13. Relevant input parameters for the simulation are shown in Table 3. The fragments first enter the plasma in the simulation at $\sim$ 0.6 ms ($t_{\text{sim}}=0$ is the shattering time). The corresponding comparison of synthetic and experimental Thomson scattering profiles is plotted in Figure 14. For reference, the pre-injection plasma profile is also plotted and compared with the experimental pre-injection profiles shown by dark blue points and dark blue line profiles. The plasma profiles after most of the material assimilation has finished are shown with light blue colors which happens around $\sim 3$ ms, 2.4 ms after the fragments enter the plasma. It can be noted that a lower $\beta=2$ value leads to higher assimilation in the plasma compared to $\beta=4$. However, the initial assimilation is limited to the edge in all 3 cases. Following the maximum material assimilation towards the edge, the density profile evolution is mainly linked to the particle diffusion. The density profiles at $t_{sim}=8$ ms are also shown with red lines in the same figure and the corresponding experimental Thomson scattering profile is plotted with red scatter points. At $t_{\text{sim}}=8$ ms, the figure shows that the core assimilation is best matched for $\beta=4$ with the edge density being slightly underestimated. However, a precise match of the simulated density will also be affected by the value of the particle diffusion coefficient and the recycling from outside the last closed flux surface. This would require additional inputs from experiments or full MHD simulations.

Comparison of the assimilation measurements for the core and edge interferometer channels in the experiments and the synthetic diagnostic are shown in Figure 15. First, in the experimental interferometer signals, a peak density rise can be observed for the core and edge interferometer channels. This density rise is presumed to be associated with the characteristic of the line of sight measurements during SPI. As the fragments ablate, they form plasmoids and the interferometer line of sight passing through a dense plasmoid can lead to overestimation of the line integrated measurement. Eventually, as the plasmoid completely assimilates in the background plasma, the interferometer signals decay to an intermediate value and the final assimilated material can be obtained. This behaviour will be confirmed with 3D JOREK simulations [49]. The INDEX synthetic diagnostic signals are closer to the final density assimilation measurements rather than the peak density measurements. The implications of this observation for comparison of experimental and simulated assimilation measurements are discussed in Section 5.2.

For the optimal determination of the $\beta$ parameter, the synthetic Thomson scattering diagnostic was preferred to the matching of the interferometer signal due to its localized measurement characteristic. Hence, the optimal value of $\beta$ was fixed to be $\beta\approx 4$. The same $\beta$ value was also used to match Thomson scattering profiles after SPI in DIII-D experiments [62]. A different value of $\beta$ might be required for different fragment sizes and fragment speeds, however, due to the lack of Thomson scattering measurements to corroborate, a fixed value is utilized.

Table 3: Input parameters for the $\beta$ scan.

  Mean number of fragments 400 Mean fragment size 0.6781 mm Mean fragment velocity 270 m/s Pellet composition(% Ne/D) 0/100 Velocity dispersion ($\Delta v/v$) 40% Local deposition fraction 0 Pellet length 4.67 mm Pellet diameter 8 mm

Figure 14: Comparison of experimental (#40743) and synthetic Thomson scattering measurements. Blue and red circular scatter points show the experimental plasma density measurements before and after SPI respectively. Plotted lines show the simulated plasma density profiles at different times.

Figure 15: Synthetic interferometer signals (green solid) for different $\beta$ values compared with experimental signals (black dotted) from the core channel (top) and edge channel (bottom) of #40743.

3.4 Example: pure deuterium injection

With the back-averaging parameter, pure deuterium injection simulations can be carried out. Since the plasma response to pure deuterium injections is significantly different to mixed deuterium-neon injections, an example simulation is described in this section. The same input SPI parameters were used as in Section 3.3 with $\beta=4$. Particle recycling from outside the last closed flux surface (LCFS) is turned off for all pure deuterium simulations. The justification for this choice is that the back-averaging model artificially deposits the ablated particles outside the LCFS and in the simulations, they can diffuse inside after deposition leading to a constant fuelling of particles and increase in plasma density. While experimentally, no significant rise in the plasma density is observed after the ablation of the fragments has finished within 3-5 ms [51]. The time trace of the simulated plasma current is shown in Figure 16. Profiles at the time stamps marked in Figure 16 with different colour vertical lines are shown in Figure 17. Again, the physical events are described below in chronological order.

1. The fragments enter the plasma at $\sim$ 0.6 ms as marked in Figure 16. As the fragments enter from the plasma edge, a rise in the neutral source function can be observed in Figure 17 (c) in the yellow curve.

2. At a later time step at 2 ms, the neutral source function penetrates to a deeper location as marked by the green curve in Figure 17 (c). Since the back-averaging model displaces the ablated material radially outwards, the penetration of the source function is limited. Green curves in Figure 17 (a), (b) show a rise in plasma density and a drop in plasma temperature at the edge of the plasma as a result of dilution cooling.

3. Further along in the simulation at 3 ms, the neutral particle source function in turquoise is strongly reduced indicating that the majority of the ablation has finished. The ablation is limited to the edge plasma resulting from the last arriving fragments. Subsequently, all the injected material is deposited however only a limited fraction is deposited inside the plasma (determined by $\beta$). Afterwards, the changes in the plasma density and temperature profiles are mainly caused by diffusion.

4. Finally, at 8 ms (corresponding to magenta colour), no source function can be observed in Figure 17(c). As the density profile evolves, the electron density diffuses inwards and the electron temperature drops towards the core. Due to the nature of dilution cooling, an inversion in the profiles can be observed i.e. as the material diffuses inwards leading to a density rise and temperature drop in the core, the edge density drops and edge temperature increases slightly.

Based on the observations above, the main difference in the plasma response for pure deuterium injections, compared to mixed deuterium/neon injections, is the absence of a radiative cold front as no impurities are injected. With negligible radiation, the current profile also remains intact. Furthermore, due to the back-averaging model, the material assimilation is limited to the edge and core density increases by diffusion. An excluded element in the simulations that might affect the onset of the TQ in pure deuterium simulations is the inclusion of intrinsic impurities.

Figure 16: Time trace of plasma current evolution for pure deuterium injection.

Figure 17: Profiles for various SPI relevant parameters against normalised poloidal flux for a pure deuterium injection.

4 Simulation results

Having established the transport and SPI models of INDEX and various AUG inputs, the main goal of this thesis, the effect of SPI parameters on penetration and assimilation can be assessed. The results shown in this thesis mainly focus on the effect of SPI on the pre-TQ phase. The following subsections address the effect of fragment sizes, speeds and pellet composition on penetration and assimilation in the plasma. The effects of different fragment sizes and speeds are categorised by the pellet composition and studied separately for mixed neon/deuterium and pure deuterium pellets due to the difference in plasma response.

4.1 Effect of fragment sizes

4.1.1 Mixed Ne/D pellets

For mixed Ne/D injection, the effect of fragment sizes is first reported in this section. A set of simulations with injections of 10% Ne and 90% D was carried out. The fragment sizes were varied indirectly by sampling a different number of fragments. For constant pellet volume, sampling more fragments leads to a smaller average fragment size. Binned distributions of fragment sizes and speeds are shown in Figure 18. 5 different realizations of fragment size and speed distributions were generated for all the cases. The mean fragment velocity was set to be 230 m/s for all the cases. The relation of the size distributions for varying number of fragments and the mean fragment diameter is shown in Figure 19. The relevant parameters for this set of simulations are summarised in Table 4.

Figure 18: Different fragment size distributions (top) and fragment speed distributions (bottom) for varying number of sampled fragments. Different distributions with the same color indicate different realizations of the probability distribution for a fixed number of sampled fragments.

Figure 19: Relation between number of fragments and the mean fragment diameter of the sampled fragments.

Table 4: Parameters for mixed deuterium/neon injections fragment size scan.

  Mean fragment velocity 230 m/s Pellet composition(% Ne/D) 10%/90% Velocity dispersion ($\Delta v/v$) 40% Local deposition fraction 1 Pellet length 10 mm Pellet diameter 8 mm

As mentioned before, the role of the DMS before disruption is to radiate away the plasma energy and increase the plasma density to avoid runaway electron formation. To study the assimilation of material for different fragment sizes in the core of the plasma, the time traces of volume averaged electron and neon density (sum of all charge states) are shown in Figure 20. It can be observed that smaller fragments start to assimilate faster in the plasma as can be observed with the earlier increase in free electron and neon density inside the $q=2$ surface. This behaviour can be attributed to smaller fragments having a relatively higher surface area that leads to faster ablation (also refer to Equation 8). As a result, the thermal energy of the plasma also starts to drop faster for smaller fragments as shown in 20(a). It should be noted that the electron and neon density time traces in Figure 20, for the case of large fragments (N=20,30,40,50,70), can cross over each other at a given point in simulation time. Due to larger statistical variations associated with sampling a small number of fragments, such behaviour can be expected depending on the distribution of fragment sizes and their corresponding speeds. Regardless, a general trend of faster assimilation for smaller fragments still holds when the fragment sizes vary drastically. The circular markers in Figure 20 indicate the time when the TQ onset condition (defined in Section 3) is satisfied. The pre-TQ duration(s) are separately plotted against the number of fragments in 20(b).

Figure 20: Average electron density (a) and neon density (b) inside the q=2 surface. Faster fragments start to assimilate quicker than larger fragments however maximum assimilation before the TQ onset (circular markers) is lower for smaller fragments.

(a)

(b)

Figure 21: Time traces of plasma thermal energy (top) and pre-TQ duration for different fragment sizes (bottom). Smaller fragments lead to a shorter pre-TQ duration.

In 20(b), smaller fragments lead to an earlier TQ onset. This behaviour can be understood by the faster deposition of neon leading to quicker radiative cooling of the edge plasma. As smaller fragments cool down the plasma edge faster than larger fragments, the penetration of material before the TQ onset is decreased leading to limited core density rise for smaller fragments. The radial profiles of temperature and density at the onset of the TQ condition are plotted against the normalised poloidal flux in Figure 22. Larger fragments are better for deeper penetration and higher assimilation before the onset of the TQ.

Figure 22: Electron temperature profiles (top) and density profiles (bottom) at TQ onset condition. Larger fragments can penetrate deeper and radiate away energy from deeper in the plasma leading to lower core temperatures and higher core density at TQ onset. q=2 surface is at $\psi_{N}\approx 0.87$.

4.1.2 Pure deuterium pellets

Pure deuterium injections would be utilised in the staggered injection scheme as discussed in Section 2.2. Since pure deuterium fragments experience a plasmoid drift, their simulation was carried out with a back-averaging model which is discussed in Section 3.3. Similar to the mixed D/Ne injections discussed in the Section 4.1.1, a scan of different fragment sizes was carried out by sampling a different number of fragments for a fixed pellet volume. The mean fragment speed was also set to 230 m/s. The main input parameters of the simulation are summarized in Table 5 and the fragment size and speed distributions are shown in Figure 23. Since these simulations do not involve injecting neon, no radiative cooling front is developed at the plasma edge. As a result, the TQ onset condition discussed in Section 3.2 is not satisfied in any of the simulations. Experimentally, it was found that pure deuterium injections with large 8 mm diameter pellets mostly lead to a disruption. A possible reason for disruptions in pure deuterium injections could be due to density limits [81], however, this has not been concluded yet. Hence, for the comparisons of the material assimilation, it should be kept in mind that experimental discharges might disrupt before the maximum simulated material assimilation is achieved. Another element in the simulations that might affect the onset of the TQ in pure deuterium simulations is the presence of intrinsic impurities which is not considered in the present simulations.

Table 5: Static parameters for fragment size scan for pure deuterium injections.

  Number of fragments [20, 50, 100, 300, 600, 800, 1000, 2000] Mean fragment diameters [2.77, 1.94, 1.50, 1.00, 0.77, 0.69, 0.64, 0.51] mm Pellet composition(% Ne/D) 10%/90% Mean fragment velocity 230 m/s Velocity dispersion ($\Delta v/v$) 40% Local deposition fraction 0 Back averaging parameter $\beta$ 4 Pellet length 10 mm Pellet diameter 8 mm

Figure 23: Different fragment size distributions (top) and fragment speed distributions (bottom) for varying number of sampled fragments.

The average electron density inside the q=2 surface and inside the plasma volume is shown in Figure 24. A trend of higher assimilation for larger fragments (corresponding to fewer number of fragments) can be observed. Additionally, as discussed in Section 3.4, material deposition in these simulations is limited to the plasma edge and core density only rises later due to diffusion. This phenomenon is indicated in Figure 24 as the volume averaged density is initially higher than the density inside q=2 surface for a given simulation. The maximum volume averaged density is reached when the the majority of ablation has taken place. After this stage, the electron density diffuses inwards and at one point the density inside q=2 surface overcomes the volume averaged density. The crossover happens earlier for larger fragments.

Figure 24: Time traces of electron density inside the plasma volume (solid lines) and inside q=2 surface (dashed lines). Larger fragments lead to a higher maximum density rise.

While the smaller fragments reach the maximum assimilated density prior to larger fragments, the magnitude of maximum assimilated material is lower compared to larger fragments. This is because smaller fragments ablate away quicker than faster fragments and penetration is shallower compared to larger fragments. To show this, the electron temperature and density profiles at time of maximum volume averaged density for each simulation is shown in Figure 25. It can be noticed that the smaller fragments lead to only edge assimilation while larger fragments can lead to comparatively deeper penetration when majority of the ablation has finished. Consequently, the plasma temperature is slightly lower towards the core for larger fragments as a result of dilution cooling of the plasma. After the maximum volume averaged density is reached, the plasma density profile is governed mainly by diffusion.

Figure 25: (a) Electron temperature profiles, (b) electron density profiles at the time of maximum volume averaged density for each case of different fragment sizes. Larger fragments can penetrate deeper and lead a higher core density rise.

4.2 Effect of fragment speeds

4.2.1 Mixed D/Ne pellets

Apart from fragment sizes, I also studied the effect of fragment speeds by carrying out a set of simulations with varying mean fragment speeds but the same fragment sizes for 10% neon and 90% deuterium composition. 200 fragments were sampled for each simulation corresponding to a mean fragment diameter of 1.137 mm. The resulting fragment size and speed distributions are shown in Figure 26. A fixed velocity dispersion of $\Delta v/v=40\%$ is used for all the fragment speed distributions in 25(b) so larger mean fragment speeds also have a larger standard deviation. This can be observed by comparing the speed distributions of $\langle v_{\text{{frag}}}\rangle=300,700$ m/s in blue and orange. The main input parameters of the simulation are summarized in Table 6.

(a)

(b)

Figure 26: Different fragment size distributions (top) and fragment speed distributions (bottom) for varying number of sampled fragment speeds. Different distributions with the same color indicate different realizations of the probability distribution for a fixed number of sampled fragments.

Table 6: Parameters for fragment speed scan for mixed deuterium/neon injections.

  Number of fragments 200 Mean fragment diameters $\sim 1.11$ mm Mean fragment velocity [300, 400, 500, 600, 700] m/s Pellet composition(% Ne/D) 10%/90% Velocity dispersion ($\Delta v/v$) 40% Local deposition fraction 1 Pellet length 10 mm Pellet diameter 8 mm

To study the assimilation, the electron and neon density is plotted inside the q=2 surface in Figure 27. Observing the time traces for different mean fragment speeds, it can be noticed that faster fragments lead to higher core assimilation and the assimilation happens quicker as well. This trend can be understood by the movement of faster fragments penetrating the plasma quicker and being able to penetrate deeper before the deposited neon cools down the plasma hereby reducing ablation. Similar to Figure 20, the circular markers indicate the time points where the TQ condition is satisfied. The dependence of pre-TQ duration on the mean fragment velocity is shown in Figure 28. A decrease in the pre-TQ duration can be observed for faster fragments.

Figure 27: Average electron density (top) and neon density (bottom) inside the q=2 surface. Faster fragments start to assimilate quicker and lead to a higher material assimilation at TQ onset (circular markers).

Figure 28: Pre-TQ duration as a function of mean fragment velocity. Faster fragments lead to a shorter pre-TQ duration.

The dependence of radiated thermal energy on the mean fragment velocity is noted in Figure 29. Faster fragments lead to a quicker drop in the plasma thermal energy due to the quicker neon accumulation, which can also be noticed in Figure 27 and hence also starts to decrease faster due to radiative cooling. The dropping rate of thermal energy is slightly faster for faster mean fragments. Similar to previous analysis, the circular (scatter) points indicate the TQ onset markers. By plotting the remaining thermal energy at TQ onset against the mean fragment speed in 28(b), faster fragments are noted to radiate away more energy before TQ onset. The higher fraction of radiated energy before the TQ for faster fragments is understood by the radiation of energy from the core because of the deeper core neon deposition by faster fragments. This is evident in the plasma profiles at the TQ onset condition, shown in Figure 30 where faster fragments lead to a larger reduction in the plasma temperature towards the core compared to slower fragments. Hence, even though faster fragments lead to an earlier TQ onset, more material can be assimilated and more of the plasma thermal energy can be radiated before the TQ onset as compared to slower fragments.

(a) Plasma thermal energy time traces for different mean fragment velocities.

(b) Remaining thermal energy at TQ onset for different mean fragment speeds. Black dashed line shows the initial thermal energy. Faster fragments radiate away more energy before TQ onset.

Figure 29:

Figure 30: Electron temperature profiles (top) and density profiles (bottom) at TQ onset condition. Faster fragments can penetrate deeper and radiate away energy from deeper in the plasma leading to lower core temperatures and higher core density at TQ onset. q=2 surface is at $\psi_{N}\approx 0.87$.

4.2.2 Pure deuterium pellets

Again using the back-averaging model, pure deuterium injection simulations were carried out by sampling different fragment speed distributions with the same fragment size distribution. The input fragment size and speed distributions for this set of simulations is shown in Figure 31 . Other relevant parameters of the simulation are mentioned in Table 7. Similar to the deuterium injections in Section 4.1.2, there is no formation of a cold front due to the lack of injected and intrinsic impurities.

Figure 31: Fragment size (top) and speed (bottom) distributions for varying mean fragment velocity.

Table 7: Simulation parameters for fragment speed scan of pure deuterium injections.

  Number of fragments 400 Mean fragment size 0.874 mm Mean fragment velocity [200, 300, 400, 500, 600, 700] m/s Pellet composition(% Ne/D) 0%/100% Velocity dispersion ($\Delta v/v$) 40% Local deposition fraction 0 $\beta$ 4 Pellet length 10 mm Pellet diameter 8 mm

To study the dependence of assimilation of the injected deuterium, the time traces of the average electron density inside the plasma volume and inside the q=2 surface for simulations with different mean fragment velocities are plotted in Figure 32. Faster fragments start to assimilate faster as a result of entering and penetrating the plasma quicker than slower fragments. The maximum material assimilation inside the plasma volume is also higher for faster fragments as compared to slower ones. In line with previous deuterium simulations discussed in Section 4.1.2, initial assimilation of the injected material is limited to the plasma edge which can be indirectly observed in Figure 32 with the volume-averaged density being higher than the average density inside the q=2 surface. In conformity with the simulations discussed in Section 4.1.2, the maximum volume averaged density is reached when majority of the ablation has finished after which the plasma density evolves only due to diffusion. As the assimilated material diffuses inwards, the average density inside the q=2 surface overcomes the volume averaged density, this crossover happening quicker for faster fragments.

Figure 32: Time traces of electron density inside the plasma volume (solid lines) and inside q=2 surface (dashed lines). Faster fragments lead to higher maximum assimilation compared to slower fragments.

The role of fragment speed is very straightforward; faster fragments can penetrate deeper (and also quicker) in the plasma before completely ablating or drifting out. Because of the higher plasma temperature towards the core, more material is ablated and assimilated compared to slower fragments. This effect can be directly observed by comparing the plasma density profiles at the point of maximum volume averaged density for each simulation as shown in Figure 33. A higher electron density can be observed for faster fragments as compared to slower fragments consistently across the plasma.

Figure 33: Electron temperature profiles (top) and electron density profile (bottom) at TQ onset for mean fragment velocities.

4.3 Effect of pellet composition

Finally, I carried out a scan of the the pellet composition by varying the neon fraction in the pellet from 0.085% to 100% with multiple simulations with trace amounts of neon ($\leq 1\%$). The lower limit in the neon composition was decided based on the pellet composition in the 2022 AUG experimental campaign. The main parameters of the simulation are shown in Table 8. Five different realisations of fragment size and speed distributions were tested for each case. It should be kept in mind that while the addition of small amounts of neon has been shown to suppress the plasmoid drift, it isn’t clear yet what fraction of neon would be sufficient to suppress the drift completely and research is ongoing. In the present simulations, it has been assumed that mixed D/Ne fragments with even trace neon fractions will have no plasmoid drift. The effect of this assumption on the simulation results is discussed later in this section.

Table 8: Parameters for pellet composition scan.

  Number of fragments 200 Mean fragment diameter $\sim$ 1.11 mm Mean fragment velocity 230 m/s Velocity dispersion ($\Delta v/v$) 40% Local deposition fraction 1 Pellet composition(% Ne) [0.085, 0.1, 0.25, 0.5, 1, 5, 10, 25, 50] %

In line with the previous comparisons, first, the time traces of electron density and neon density inside the q=2 surface for injections with different neon fractions are shown in Figure 34. Similar to Figure 20, the circular markers indicate the time point when the TQ onset condition is satisfied. The pre-TQ duration is plotted separately against the neon fraction in Figure 35. It can be observed that the electron assimilation is higher for injected fragments with less neon and the neon assimilation is higher for injected fragments with more neon. To study the assimilation of neon in the core of the plasma, 35(a) shows the the assimilated neon atoms in the plasma inside the q=2 surface as a function of the injected neon atoms. For trace neon fractions, the assimilated neon atoms increases linearly with increasing neon fraction. However, at larger neon fractions, the assimilated neon quantity increases at a much slower rate and is almost saturated. This trend can be explained by swift cooling of the plasma for high amount of injected neon. Beyond a certain amount of neon, 5% 8 mm pellets or $>10^{21}$ atoms from 35(a), the plasma cools down very drastically (also refer to Figure 37) and hence only a limited amount of material can be ablated leading to a self-regulated saturation. To study the dependence of the ablation fraction i.e. ratio of ablated neutral atoms to the number of injected atoms on the injected neon, time traces of ablation fraction for varying neon fractions are plotted in 35(b). The reduced ablation fraction for larger neon fractions is very apparent. Correspondingly, decreasing pre-TQ duration for larger neon fraction can also be observed in Figure 35.

Figure 34: Average electron density (top) and neon density (bottom) inside the q=2 surface. Electron assimilation is higher for injections with less neon while neon assimilation is higher for injections with more neon.

Figure 35: Pre-TQ duration as a function of neon fraction. Higher neon fraction injections lead to a shorter pre-TQ duration.

(a)

(b)

Figure 36: (a) Neon atoms inside the q=2 surface at the TQ onset against neon fraction in the injected pellet. The assimilated neon increases linearly with injection neon amount but saturates for higher injected neon quantities; (b) Ablation fraction time traces for varying neon fractions.

Since large fraction mixed neon pellets will have to radiate away plasma energy before the TQ onset, the thermal energy time evolution was studied for different neon fractions and is shown in Figure 37. While the initial fragments enter the plasma at $\sim$ 0.7 - 0.8 ms, the decay of the thermal energy due to radiative cooling starts at $\sim$ 1.1 ms. A faster rate of thermal energy loss due to radiative cooling can be observed for increasing neon fractions. A behaviour similar to the assimilation can be seen i.e. for trace neon quantities ($<1\%$), large differences are observed in the rate at which thermal energy decreases. However, for larger neon fractions, the thermal energy is radiated at very similar rates although the trend of higher decay rate for higher neon fraction persists. Looking at the TQ onset markers in Figure 37, a non-monotonic trend of remaining thermal energy can be observed when going from small to large neon fractions which has a minima at $0.5\%$ neon cases. This can be more clearly seen by relating the remaining thermal energy at TQ onset to the neon fraction as plotted in Figure 38.

Figure 37: Time traces of plasma thermal energy evolution for varying neon fraction.

Figure 38: Remaining thermal plasma energy at TQ onset condition. Black solid line shows initial thermal energy. A non-monotonic behaviour can be observed with minimum values for trace neon injections.

The cause of the non-monotonic behaviour can be understood by looking at the plasma temperature and radial profiles at the TQ onset which are plotted in Figure 39. For the trace neon injections (in blue), Figure 39 shows an inside out temperature collapse where the core temperature drops lower than 10 eV. While for large neon fractions, an outside-in temperature collapse occurs as a result of the radiative cold front moving in. The inside-out temperature collapse for trace neon injections can be understood as follows. As the fragments penetrate the plasma edge, they start to ablate, however, the ablation rate remains low (refer to Equation 8), leading to a slow deposition of neon atoms for trace neon fragments. Further, as the fragments travel deeper in the core, the ablation rate increases due to the hot plasma leading to a higher core deposition of neon atoms as compared to the edge. As a result, the neon starts radiating strongly from the core while the edge still remains relatively hot compared to the higher neon fraction injections. As a result of the inside-out temperature collapse, a larger fraction of the thermal energy can be radiated before the TQ onset condition is satisfied. As mentioned before, the plasmoid drift in the simulations has been enabled for trace neon injections. This assumption can affect the trend of the stronger core deposition for trace neon injections as discussed above. The trends that would likely be affected would be the non-monotonic trends in remaining thermal energy at TQ onset (Figure 38) and the inside-out temperature collapse for trace neon injection (Figure 39).

Figure 39: Electron temperature profiles (top) and electron density profile (bottom) at TQ onset for different neon fractions. An inside-out temperature collapse can be observed for trace neon fraction injections in contrast with higher neon fraction injections where the cold front moves in from the edge.

5 Experimental comparisons

A preliminary comparison with experimental data is reported in this section. I carried out some of the experimental analysis regarding penetration and assimilation and the findings are available in my ITER internship report [51]. However, I also rely on other experimental results carried out by other members of the AUG and ITER teams. Comparisons between the simulations and experiments are carried out for dependence of fragment sizes, speeds and composition on penetration, material assimilation. Additionally, the simulated pre-TQ duration is also compared with experimental pre-TQ duration.

5.1 Penetration

Starting with penetration, fast visible cameras [21, 32, 51] measured $D_{\alpha}$ radiation in the experimental campaign and were used to gauge the penetration of the fragments for pure deuterium injections. I found that larger fragments (with similar mean fragment velocities) showed stronger ablation towards the core for pure deuterium injections. In comparing fragment speeds, higher fragment speeds penetrated much deeper and also showed stronger ablation towards the core as compared to slower fragments. Representative experimental penetration plots are showed in Figure 40 and Figure 41. For both these figures, the penetration was studied by mapping the fast camera and AXUV bolometer line of sights to a radial location (y-axis in Figure 40, Figure 41) based on their intersection with the injection vector. More details of the mapping can be found in my ITER internship report [51]. Figure 40 compares the penetration for different fragments sizes. Higher radiation was observed from the intermediate plasma region for larger fragments as compared to smaller fragments. Similarly, faster fragments penetrated much deeper in the plasma compared to slower fragments as shown in Figure 41.

(a) #40655: 2D AXUV mesh plot, larger fragments ($D_{2}$ SPI, $12.5^{\circ}$ rect., 442.9 m/s)

(b) #40656: 2D AXUV mesh plot, smaller fragments ($D_{2}$ SPI, $25^{\circ}$ rect., 471 m/s)

Figure 40: 2D AXUV plots for #40655 [top] and for #40656 [bottom]. Larger fragments (#40655) lead to deeper penetration as compared to smaller fragments (#40656). In both plots, vertical blue dashed and red dashed lines mark the $I_{p}$ dip and spike respectively. Horizontal lines mark different radial locations of interest; Black dotted: Separatrix, Red dotted: q=3, yellow dotted: q=2, blue dotted: q=1.5 and cyan dotted: Magnetic axis $R_{0}$.

(a) 41012 2D camera plot, slower fragments ($D_{2}$ SPI, $25^{\circ}$ rect., 342 m/s, with $D_{\alpha}$ filter)

(b) 40772 2D camera plot, faster fragments ($D_{2}$ SPI, $12.5^{\circ}$ rect., 700 m/s, without $D_{\alpha}$ filter)

Figure 41: 2D camera plots for #41012 [top] and for #40772 [bottom]. Faster fragments (#40772) lead to deeper penetration as compared to slower fragments (#41012). In both plots, vertical blue dashed and red dashed lines mark the $I_{p}$ dip and spike respectively. Horizontal lines mark different radial locations of interest; Black dotted: Separatrix, Red dotted: q=3, yellow dotted: q=2, blue dotted: q=1.5 and cyan dotted: Magnetic axis $R_{0}$. Note the the two plots have different time duration.

Comparing the pure deuterium simulation results in Figure 25 and Figure 33, larger and faster fragments penetrate and deposit material deeper in the plasma at the time when most of the ablation has finished. A more direct indication of fragment penetration can be obtained from the centre of mass of the fragments shown in Figure 42. Again, larger and faster fragments lead to higher penetration as can be seen by the final location of the centre of mass along the injection vector. Quantitative matches of material penetration are challenging to perform due to limitations in experimental measurements of penetration of individual fragments and limited viewing angles.

(a) Variations in fragment sizes (small sizes for a higher number of fragments).

(b) Variation in fragment speeds.

Figure 42: Simulated centre of mass of the injected fragments for variations in fragment sizes and speeds plotted in the poloidal cross section. Different contour levels are surfaces of constant toroidal flux. Black arrow shows the injection vector. Deeper penetration of the injected material can be observed for larger and faster fragments.

5.2 Assimilation

Assimilation measurements in the AUG experimental campaign for the pure deuterium injections were made using the COO core and edge interferometer channels shown in Figure 8. Trends in absolute assimilation for different fragment sizes and speeds in pure deuterium injections were studied in Patel, 2023 and Jachmich et. al., 2023 [53, 51] and are shown in Figure 43. It should be noted that the assimilation estimates are provided with the peak density rise from the interferometer (refer to Figure 15). With increasing impact velocity, i.e. decreasing mean fragments size, an overall decrease in assimilation can be observed. At the same perpendicular velocity, the $12.5^{\circ}$ sq. shatter tube indicated by diamond data points would have a higher mean fragment velocity indicated by a higher injection speed.

Figure 43: Peak rise in electron density from core interferometer channel (refer to Figure 8) normalized to injected deuterium atoms against impact velocity. Figure taken from S. Jachmich et. al. [53]. With increasing impact velocity ($\sim$ decreasing fragment sizes) the absolute assimilation decreases. Comparing data points with higher injection velocities, marked by ’$V_{\text{inj}}$’, higher assimilation is found for faster injections.

For assimilation comparison with the simulations, the absolute assimilation for the set of pure deuterium injection simulations carried out in Section 4.1.2 and Section 4.2.2 can be referred to. While the trends shown in the aforementioned sections indicate larger and faster fragments being better for assimilation, a quantitative comparison of the assimilation fraction can be carried out. Hence, the assimilation fraction for the size and speed scans from the synthetic interferometer diagnostic is shown in Figure 44.

(a) Fragment sizes.

(b) Fragment speeds.

Figure 44: Simulated absolute electron assimilation fraction for varying fragment sizes and speeds.

Assimilation measurements from the synthetic interferometer signals match the experimental trends of larger and faster fragments being better for higher assimilation. However, a striking difference in the absolute assimilation fraction between the simulations and the experiments is that the maximum experimental assimilation fractions lie between 25% and 75% while the simulated assimilation fraction lie between 6% and 12%. The main reason behind this significant discrepancy is that the experimental assimilation measurements take into account the peak density rise in the interferometer signal which has been preliminary observed to $\sim$ 2-8 times larger than the value after the peak rise (refer to Section 3.3) although a more rigorous experimental analysis is being carried out by the AUG and ITER teams. The extent of density rise for fragments with different sizes or speeds might also be affected by this procedure. Regardless, the initial experimental analysis suggests larger and faster fragments being better for increased assimilation as also observed in the simulations.

For mixed D/Ne injections, comparison with qualitative assessments of neon assimilation can be carried out. S. Jachmich et. al. [53] used the plasma current decay rate in the CQ phase to study the assimilated neon in the plasma. Higher injected neon amounts were reported to lead to faster current decay indicating a higher amount of assimilated neon which was also observed in simulations (see also 35(a)). First results on the radiated energy fraction ($f_{\text{rad}}$) for SPI induced disruptions had been reported for AUG [52]. The radiated energy fraction was defined as

$$\mathrm{f}_{\text{rad }}=\frac{\mathrm{W}_{\text{rad }}}{\mathrm{W}_{\text{mag }}+\mathrm{W}_{\text{th }}-\mathrm{W}_{\text{coupled }}+\mathrm{W}_{\text{heating }}}$$ (11)

and should ideally be close to 1 for an efficient mitigation scheme. Here, $W_{\text{rad}}$ is the radiated energy, $W_{\text{mag}}$ is the magnetic energy, $W_{\text{th}}$ is the thermal energy, $W_{\text{coupled}}$ is the magnetic energy coupled with the surrounding structures and $W_{heating}$ accounts for external heating sources. The results from the 2022 AUG experimental campaign are shown in Figure 45 and indicate that the radiation fraction increases with increasing amount of injected neon, however, for $n_{\text{atoms,Ne}}>10^{21}$, most of the plasma thermal energy is radiated and the curve saturates [52]. A similar trend of simulated neon assimilation in the pre-TQ phase can be observed in 35(a). It should be noted that the plasma dynamics during the TQ and CQ phase might influence the neon assimilation and hence the radiated energy fraction. By carrying out full disruption simulations in INDEX, a more direct comparison can be made.

Figure 45: Radiated energy fraction ($f_{\text{rad}}$) as a function of the injected neon atoms. Different colors indicate the pellet diameter and different markers indicate the shattering tube used for the injection. $f_{\text{rad}}$ increases with increasing neon concentration. Figure taken from P. Heinrich et. al. [52].

5.3 Pre-TQ duration

The pre-TQ duration is an important parameter for the DMS as a long pre-TQ duration can enable multiple injections to meet the material assimilation requirements at ITER. Hence, the pre-TQ duration for AUG SPI discharges was compared to the simulated pre-TQ duration for different fragment sizes and speeds. For the simulated discharges, the pre-TQ duration is defined in Section 3.2 and is related to the cold front reaching the $q=2$ surface. The formation of a cold front in simulations is only found in simulated mixed D/Ne injections and therefore, the pre-TQ duration can only compared for these injections. While qualitative trends in pre-TQ duration are compared, it should be noted that quantitative comparison might require MHD simulations to take into account the delay from MHD mode onset (refer to Section 2) to grow enough to lead to a TQ.

For the mixed Ne/D injections, I have plotted the pre-TQ duration in Figure 46 for different fragment speeds and similar fragment sizes that were available for 10% and 40% mixed D/Ne injections. At similar impact velocity, higher pre-shattering pellet speeds ($\sim$ higher fragment speeds) lead to a shorter pre-TQ duration. A shorter pre-TQ duration for faster fragments has been observed in simulations in Section 4.1.1, particularly in Figure 28. The absolute value of the pre-TQ duration in the experiments differ from those found in the simulations however an exact comparison is difficult because of (a) the specific TQ onset condition used in the simulations, (b) differences in the modelled and experimental fragment size and speed distributions [45] that modify the time taken by the cold front to reach the $q=2$ surface, (c) time taken by MHD activity growth for the TQ onset. Nevertheless, the parametric scans in this thesis can reproduce the shorter pre-TQ trend for faster fragments.

Figure 46: Experimental pre-TQ duration for varying perpendicular velocity for 10% neon and 40% neon pellets. $v_{pellet}$ is the pre-shattered pellet speed. Larger pellet speed ($\sim$ mean fragment speed) leads to a shorter pre-TQ duration.

To assess the effect of fragment sizes, I have plotted the pre-TQ duration against the mean fragment velocity in Figure 47. Comparable injections with similar fragment sizes but different fragment speeds were available for 10% neon mixed deuterium/neon injections. The mean fragment velocity is assumed to be the mean of the pellet velocity and the parallel component of the pellet velocity with respect to the shattering tube as per the discussion in Section 2.4. For the same mean fragment velocity, a higher impact velocity $v_{\perp}$ i.e. smaller fragments lead to a longer pre-TQ duration. This is in contrast to the simulated pre-TQ duration for varying fragment sizes as shown in 20(b) where the pre-TQ duration is longer for larger fragments. Two possible reasons can be attributed to the contrasting trends of pre-TQ duration. First, shattering parameters leading to large fragment sizes in experiments (lower $v_{\perp}$) were often accompanied with a significantly higher number of small fragments as compared to the theoretical fragmentation model [45] used in present simulations. In the presence of smaller fragments, which ablate quicker, the edge cooling might be accelerated leading to modification in the pre-TQ duration. Hence, the fragmentation model itself might have limited validity in the case of low impact velocities. Second, larger fragments might perturb the plasma more significantly due to more localized cooling, leading to stronger MHD activity and causing a disruption quicker in the experiments [46]. To assess the impact of both these parameters, additional inputs are required from pellet shattering experimental analysis and 3D MHD modelling activities.

The comparisons above were carried out with the experimental data points at 200 m/s due to the simulations carried out in Section 4.2 having similar mean fragment velocities. However, at $\sim$ 450 m/s in Figure 47, a set of data points (green circle, black square) indicate that larger fragments lead to a longer pre-TQ, in line with the simulations. The possibility of these points being valid data points (and not statistical outliers) have to be considered motivating the requirement of more experimental data points.

Figure 47: Experimental pre-TQ duration for varying mean fragment velocity for 10% neon pellets. Larger impact velocity ($\sim$ smaller fragments) lead to a longer pre-TQ duration.

6 Summary and outlook

I employed the 1.5D INDEX code in this thesis to model the plasma response before the thermal quench in the ASDEX Upgrade tokamak during shattered pellet injections. I studied the influence of three critical SPI parameters: fragment sizes, velocities, and pellet composition, with a particular focus on their impact on the penetration and assimilation of the injected material within the plasma. To simulate the effects of plasmoid drift during pure deuterium injections, a back-averaging model was used. I carried out interpretive simulations to determine the back-averaging parameter by matching assimilation measurements from experimental and synthetic Thomson scattering measurements.

In the context of mixed deuterium/neon injections, it is observed that larger and faster fragments result in increased penetration depth and higher material assimilation for an equal quantity of injected material. By applying a semi-empirical onset condition for thermal quench (TQ), the pre-TQ duration can be estimated. For varying fragment sizes, the pre-TQ duration is longer with increasing fragment sizes. Conversely, when investigating varying fragment speeds, faster fragments correspond to a shorter pre-TQ duration. The phenomenon of larger fragments leading to a longer pre-TQ phase can be attributed to their greater surface-area-to-volume ratio, enabling them to traverse deeper into the plasma and deposit more material before the plasma cools down, resulting in higher assimilation and an extended pre-TQ duration. On the other hand, faster fragments can achieve deeper penetration and higher assimilation compared to their slower counterparts, even for a shorter pre-TQ duration. While experimental data of the pre-TQ duration aligns with simulation results for variations in fragment speeds, there is an inconsistency for variations in fragment sizes. A possible reason for this mismatch may stem from disparities between the experimental and modelled fragment size distributions or MHD effects. Quantitative comparisons of the pre-TQ duration might require MHD modelling to take into account the time taken by destabilized MHD modes to grow and lead to the TQ.

For pure deuterium injections, I carried out simulations assuming the same back-averaging parameter for all fragment sizes and speeds. In line with the mixed deuterium/neon injections, larger and faster fragments lead to higher penetration and assimilation although the material assimilation is initially limited to the plasma edge. The core plasma density increases later on a diffusive time scale after most of the material has assimilated. Experimental comparisons of the material penetration in pure deuterium injections are in line with the simulated trends. Trends of experimental material assimilation also align qualitatively, with larger and faster fragments leading to higher assimilation. However, the absolute experimental assimilation values appear roughly 3-6 times higher than the simulated assimilation. The main source of this discrepancy can arise from the limitation in the experimental assimilation values. More accurate matching of assimilation trends requires further analysis of the experimental data.

I also carried out simulations with variations in the pellet composition. Neon assimilation increased for increasing amounts of injected neon content. However, a self-regulating saturation of the neon assimilation was observed for neon content larger than $10^{21}$ atoms. The underlying cause of this effect is higher amounts of neon leading to accelerated cooling of the plasma, which, in turn, imposes limitations on further neon assimilation. Comparisons between simulated neon assimilation and indirect experimental measurements of neon assimilation were carried out. Experimental CQ decay rates increased with increasing neon content indicating higher neon assimilation as seen in the simulations. I noticed a similar trend for the experimental radiated energy fraction, which also increased with higher neon content. The pre-TQ duration monotonically decreased with increasing neon concentration. Trace neon injections ($<10^{21}$ atoms) lead to an inside out temperature collapse in contrast with higher neon concentration injections where the plasma edge cooled down first. As a result, the highest amount of plasma energy was radiated before the TQ onset for trace neon injections. An important limitation for the trace neon injections is that the plasmoid drift can limit the material assimilation to the edge, however it wasn’t included for these simulations and can affect the inside-out temperature collapse phenomenon. This limitation also motivates the use of improved plasmoid drift models.

As mentioned throughout the thesis, there are various other additional activities that will improve the simulations and matching with experiments to quantify the effect of different fragment sizes and speeds. These effects are summarised and further elaborated upon below:

• Experimental vs modelled fragment size distribution
  As pointed out in Section 5, disparities between the fragment sizes in experimental and modelled distributions have been found [45]. In the case of pellets with low impact velocities, the modelled fragment size distribution underestimates the amount of small fragments. As smaller fragments ablate faster as discussed in Section 4.1.1, an underestimation of smaller sized fragments can affect the experiments differently than simulations leading to changes in assimilation and pre-TQ duration.

• Gas production
  During a SPI, there are two possible avenues of gas entering the plasma. If the injector uses propellant gas to accelerate the pellet, as is the case in AUG, it is possible that the propellant gas can enter the vessel before the pellet. Additionally, gas can also be produced during the shattering process. It has already been shown in DIII-D SPI experiments [82] that the gas entering the plasma can reduce the plasma energy content and increase the plasma density before the fragments arrive in the plasma. The effect of gas production can be taken into account in the simulations to further improve matching with experimental pre-TQ duration.

• Plasmoid drift
  INDEX uses the back-averaging model to simulate the plasmoid drift affecting pure deuterium injections. A fixed back-averaging parameter $\beta$ has been used for all fragment sizes and speeds in this work. Plasma density measurements during pellet assimilation or inputs from 3D MHD simulations can assist in determining $\beta$ for different fragment sizes and speeds. Improved plasmoid drift models can also lead to more accurate estimates of assimilation for pure deuterium and trace neon injections.

• Assimilation measurements
  As discussed in Section 3.3, the present assimilation estimates have been obtained using the peak interferometer signal values. If the peak value corresponds to the line of sight measurement passing through a dense plasmoid which eventually dissipates in the background plasma, then the final assimilation measurements might be more reflective of the assimilated material. To quantitatively validate the assimilation measurements of INDEX, the final experimental density should be estimated.

• Current Quench simulations
  To quantify the assimilated neon content in the experiments, CQ simulations of the AUG discharges can be carried out. These simulations can be initialized with a post-TQ plasma with assimilated deuterium and neon. As the plasma current decay rate in the CQ phase is strongly dependent on the impurity radiation rate, the amount of assimilated neon determines the current decay rate. By carrying out a scan of assimilated neon, the experimental current decay rate can be matched with the simulated current decay rate. This process can be carried out for injections with different amounts of assimilated neon to get an estimate of the assimilated neon amount. While some preliminary simulation setup work was carried out in the duration of this thesis in collaboration with Dr. Akinobu Matsuyama, additional setup of plasma control in the CQ needs to be implemented.

The present simulations suggest that larger and faster fragments are beneficial for deeper penetration and higher assimilation for both pure deuterium and mixed deuterium/neon injections. However, for better validation of the INDEX code, quantitative comparisons of material assimilation for different fragment sizes and speeds should be carried out. This will require improved experimental assimilation measurements. Furthermore, CQ simulations can help in assessing the quantitative neon assimilation. Further simulations and experiments are also required to confirm the effect of different fragment sizes in determining the trends of pre-TQ duration for mixed deuterium/neon injections. After validating the INDEX code using AUG results on the effect of different fragment sizes, speeds and pellet composition, further simulations of different SPI parameters can be carried out for ITER plasma discharge scenarios. However, quantitative extrapolations of the aforementioned parameters to ITER depends on various other parameters, some of which are mentioned below:

• Thermal energy content
  The assimilation of material in the plasma depends strongly on the thermal energy content of the plasma due to the dependence of the ablation rate on the plasma temperature. It has been reported in AUG SPI discharges as well [53] that a higher thermal energy content leads to an increased assimilation. Hence, optimal fragment fragment parameters might depend on the target plasma. For this purpose, the INDEX code should be validated on different tokamaks with different thermal energy content.

• Pre-existing MHD activity
  In the simulations carried out in this thesis, SPI was applied to an H-mode ”healthy plasma” without any MHD activity due to a higher availability of experimental data to compare to. However, in practice, the DMS might be triggered when the plasma already has MHD activity and/or might have transitioned to L-mode [11]. Wieschollek et. al. [83] studied the role of exiting MHD activity on plasma dynamics during SPI in AUG using JOREK and showed that pre-existing large magnetic islands might affect the pre-TQ duration.