Magnetohydrodynamic simulation study of impurity radiation-excited and driven tearing mode

1.   Introduction

High-density operation represents the fundamental operational mode for future large-scale fusion devices [1]. Experimental findings reveal a limit to enhancing plasma density, with the Greenwald density limit standing as the most widely recognized empirical scaling law [2, 3]. When plasma density approaches the Greenwald density limit, a variety of phenomena, such as MARFEs [4], current channel shrinkage [5], and MHD instability bursts [6] were observed. These phenomena can lead to a degradation of plasma confinement and can even result in a major disruption [7]. The physical mechanism behind the onset of the density limit remains an open problem. Both experimental [8] and simulation [9] studies have indicated that the tearing mode is the most dangerous mode for triggering density limit disruptions. Impurity radiation significantly influences the driving of tearing mode [10], with important applications in turbulent transport [11], MHD instability [12], divertor physics [13], and disruption mitigation [14]. Impurity modules have been developed and implemented in MHD codes such as JOREK [15], M3D-C1 [16], and NIMROD [17]. In this work, we incorporate impurity modules into the CLT code to study the excitation and driving mechanisms of tearing mode.

The stability of the tearing mode is dependent on the shape of the q profile [18]. Since the density limit is associated with the tearing mode, it seems that the density limit should be associated with the shape of the q profile. Although it is challenging to experimentally measure the shape of the q profile, numerous experiments have demonstrated that there is no apparent relationship between the density limit and the shape of the q profile [19–21]. It is generally believed that radiation collapse can lead to plasma contraction and create a current gradient on the resonance surface, thereby exciting the tearing mode [22, 23]. Further experimental and simulation results are necessary to support this viewpoint. In this work, the tearing mode excited by impurity radiation collapse is studied using MHD simulation.

Gate et al proposed a radiation-driven island model that links impurity radiative cooling with the growth of magnetic islands to explain the physical mechanism of density limit [24, 25]. Teng et al provided a quantitative prediction of the density limit by establishing plasma balance and internal power balance inside the magnetic island [26]. They semi-analytically calculated the growth of the radiation-driven tearing mode. Then they verified the semi-analytically model with 3D MHD simulations and found that cooling inside the magnetic island can significantly increase the width of the magnetic island [9]. In this work, the self-consistent process of impurity collapse and impurity radiation-driven magnetic island is studied by using 3D MHD simulations. This study provides a valuable complement to previous research.

This paper is organized as follows. In section 2, we present the KPRAD model with MHD equations used in CLT code. The simulation of impurity radiation-excited tearing mode is given in section 3. Simulations of impurity radiation-driven magnetic islands with the scanning of initial impurity parameters are presented in section 4. Section 5 gives the discussion and summary of this study.

2.   Numerical scheme

The three-dimensional toroidal resistive MHD equations used in the CLT code [27–30], coupled with the impurity radiation effect are given as follows:

where $\rho ,\;\boldsymbol{v},\;p,\;\boldsymbol{J},\;\boldsymbol{B},$ and E are the plasma density, the fluid velocity, the plasma pressure, the plasma current density, the magnetic field, and the electric field, respectively. ${D}_{\perp },\;{D}_{\parallel },\;{\kappa }_{\perp },\;{\kappa }_{\parallel },\;\nu ,$ and $\eta$ are the parallel and perpendicular plasma diffusion coefficient, the parallel and perpendicular thermal conductivities, the viscosity, and the resistivity, respectively. In this paper, the specific heat ratio is defined as $\mathrm{\Gamma }=5/3$. ${P}_{\mathrm{r}\mathrm{a}\mathrm{d}}$ represents the impurity radiative loss rate which can be computed through the KPRAD model [31]:

where Z, ${n}_{\mathrm{e}}$, ${n}_{Z}$, ${L}_{Z}$, ${I}_{Z}$, ${R}_{Z}$ and ${S} _{Z}$ are the impurity atomic number, the electron density, the impurity charge-state density, the radiation cooling rate, the ionization coefficient, the recombination coefficient, and the impurity source density, respectively. ${L}_{Z}$, ${I}_{Z}$, and ${R}_{Z}$ can be read from the ADPAK atomic database [32]. ${S} _{Z}$, as an impurity source term added to the plasma boundary region over time, can be manually adjusted.

All the variables in the CLT code are normalized as follows: $\rho /{\rho }_{00}\to \rho$, $\boldsymbol{v}/{v}_{\mathrm{A}}=\boldsymbol{v}$, $p/({B}_{00}^{2}/{\mu }_{0})\to p$, $\boldsymbol{B}/{B}_{00}\to \boldsymbol{B}$, $\boldsymbol{J}/({B}_{00}^{2}/{\mu }_{0}a)\to \boldsymbol{J}$, $\boldsymbol{E}/({v}_{\mathrm{A}}{B}_{00})\to \boldsymbol{E}$, $t/{t}_{\mathrm{A}}\to t$, ${D}_{\perp }/({a}^{2}/{t}_{\mathrm{A}}) \to {D}_{\perp }$, ${D}_{\parallel }/({a}^{2}/{t}_{\mathrm{A}})\to {D}_{\parallel }$, ${\kappa }_{\perp }/({a}^{2}/{t}_{\mathrm{A}})\to {\kappa }_{\perp }$, ${\kappa }_{\parallel }/({a}^{2}/{t}_{\mathrm{A}})\to {\kappa }_{\parallel }$, $\nu /({a}^{2}/{t}_{\mathrm{A}})\to \nu$, and $\eta /({a}^{2}/{t}_{\mathrm{A}})\to \eta$, respectively. $a$ is the minor radius, ${\rho }_{00}$ is the plasma density, ${B}_{00}$ is the magnetic field at the magnetic axis, ${t}_{\mathrm{A}}=a/{v}_{\mathrm{A}}$ is the Alfvén time, and ${v}_{\mathrm{A}}={B}_{00}/ \sqrt{{\mu }_{0}{\rho }_{00}}$ is the Alfvén speed, respectively.

In equation (2), the term $-\nabla \cdot \left[{\kappa }_{\perp }\nabla \left({p}_{0}\right)\right]$ represents a heat source, equivalent to the initial pressure diffusion. In equation (5), the term $-\eta {\boldsymbol{J}}_{0}$ represents the magnetic flux pump, the pump rate is equivalent to the resistive diffusion of the initial magnetic field. The presence of these heat and flux pumping terms prevents the initial equilibrium from diffusing, simplifying the analysis of instability development. To maintain $\nabla \cdot \boldsymbol{B}=0$, the CLT code independently solves Faraday’s law (equation (4)) and Ohm’s law (equation (6)).

In CLT, spatial derivatives are obtained using a fourth-order finite difference scheme, while time integration employs the fourth-order Runge–Kutta scheme. The cut-cell method is chosen for addressing boundary problems due to the physical boundary not coinciding with the grid points [33]. In this study, a fixed boundary condition is applied to all variables. The uniform grids with ($256\times 32\times 256$) in $(R,\phi ,Z)$ are used in all simulations in this paper.

3.   The excitation of the tearing mode with impurity radiative cooling

The initial equilibrium in this study is computed using the NOVA code [34]. The initial profiles of the safety factor, the plasma pressure, and the toroidal plasma current density used in this section are shown in figures 1(a)–(c). The initial impurity source distribution ${S} _{0}$ is shown in figure 1(d). In this study, we set ${S} _{Z}$ for different Z values all equal to ${S} _{0}$.

Figure  1.  The initial equilibrium profiles of (a) the safety factor, (b) the plasma pressure, and (c) the toroidal plasma current density, and (d) the initial impurity source distribution.

DownLoad: Full-Size Img PowerPoint

The safety factor value is ${q}_{0}=1.2$ on the magnetic axis and ${q}_{\mathrm{e}}=4$ at the edge. The plasma pressure on the magnetic axis is ${p}_{0}=8\times {10}^{-4}$. This initial equilibrium has a broad uniform pressure distribution in the central region. The stability parameter of the tearing mode [35] is ${\mathrm{\Delta }}'_{0}=-3.08$, indicating that the $m/n=2/1$ tearing mode is stable in this equilibrium, where m and n are the poloidal and toroidal mode numbers. The parameters in this simulation are ${D}_{\perp }=1\times {10}^{-6}$, ${D}_{\parallel }=1\times {10}^{-3}$, ${\kappa }_{\perp }=1\times {10}^{-10}$, ${\kappa }_{\parallel }=1$, $\nu = 1\times {10}^{-6}$, and $\eta \propto {Z}_{\mathrm{e}\mathrm{f}\mathrm{f}}{T}_{\mathrm{e}}^{-3/2}$. Here, ${Z}_{\mathrm{e}\mathrm{f}\mathrm{f}}$ represents the effective charge number of impurity, and ${T}_{\mathrm{e}}$ is the electron temperature, both of which can be calculated using the KPRAD model [31]. In this simulation, carbon ($Z=6$) is chosen as the impurity. The $m/n=2/1$ perturbation is applied at the beginning of this simulation.

For this equilibrium, it will be stable for all n modes without impurity. We chose carbon impurity to be applied to the plasma boundary in this case, and figure 2 illustrates the evolutions of the kinetic energy and the growth rates for different toroidal modes. The high-n modes are excited by impurity radiative cooling. The duration of these modes in the linear stage is the same. Although the growth rate of the $n=1$ mode is smaller than those of the other modes, the amplitude of the $n=1$ mode is still much larger than those of the other modes because the $n=1$ mode is initially imposed. At the nonlinear stage, the 2/1 tearing mode becomes destabilized and boosted by the higher-n modes due to the nonlinear mode coupling.

Figure  2.  Evolutions of (a) the kinetic energy and (b) the growth rates for different toroidal modes.

DownLoad: Full-Size Img PowerPoint

The contour plots of the perturbed temperature at four typical moments are shown in figure 3, while figure 4 illustrates the toroidal electric field of different modes at the same moments. At $t=2576{t}_{\mathrm{A}}$, the plasma temperature is cooling due to the impurity radiation, and the temperature gradient starts to build up gradually just outside the $q=2$ surface, as shown in figure 3(a). In the early nonlinear stage, the $m/n=5/2,\ 6/2,\ 9/3$ modes are exited while the $m/n=2/1$ mode remains almost unchanged, as shown in figures 4(a) and (b). At $t=2790{t}_{\mathrm{A}}$, it marks the onset of the nonlinear stage (figure 2(b)). The temperature gradient between the $q=2$ and $q=3$ surface becomes very large (figure 3(b)), which causes the unstable higher-n modes quickly develop to become the dominant mode (figure 4(b)). With the further buildups of the temperature gradient inside the $q=3$ surface, the higher-n mode ($m/n=12/4$) becomes a dominant mode at $t=3005{t}_{\mathrm{A}}$, as shown in figures 3(c) and 4(c). At the same time, the $m/n=2/1$ tearing modes also become destabilized and developed on the $q=2$ surface. However, the amplitude of the 2/1 tearing mode is considerably smaller than that of the 6/3 mode (figure 4(c)). Therefore, it is difficult to identify the 2/1 tearing mode at this stage as shown in figure 3(c). As the impurity penetrates into the magnetic island and reaches a certain level as shown in figure 5, the impurity radiation cooling results in the quick growth of the 2/1 magnetic island. At $t=3220{t}_{\mathrm{A}}$, the amplitude of the 2/1 tearing mode becomes similar to that of the 6/3 mode (figure 4(d)), and the structure of the 2/1 tearing mode is discernible (figure 3(d)). Due to the cooling between the $q=2$ and $q=3$ surfaces, a temperature gradient forms at the $q=2$ surface and inside the $q=3$ surface, leading to the emergence of high-n modes ($m/n=12/4,\ 9/3$) inside the $q=3$ surface, as shown in figures 4(c) and (d).

Figure  3.  Contour plots of perturbed temperature at (a) $t=2576{t}_{\mathrm{A}}$, (b) $t=2790{t}_{\mathrm{A}}$, (c) $t=3005{t}_{\mathrm{A}}$, and (d) $t=3220{t}_{\mathrm{A}}$.

DownLoad: Full-Size Img PowerPoint

Figure  4.  The toroidal electric fields of different modes at (a) $t=2576{t}_{\mathrm{A}}$, (b) $t=2790{t}_{\mathrm{A}}$, (c) $t=3005{t}_{\mathrm{A}}$, and (d) $t=3220{t}_{\mathrm{A}}$.

DownLoad: Full-Size Img PowerPoint

To illustrate the driving effect of impurities on the tearing mode, figure 5 presents the temporal evolutions of total impurity density at O-point and X-point, as well as the width of the 2/1 magnetic island. During the formation of the magnetic island, impurities enter the island through the X-point, accumulate at the O-point, and generate radiative cooling, driving the growth of the island. Consequently, the impurity density at the O-point is higher than that at the X-point. Before $t=2000{t}_{\mathrm{A}}$, the impurity density within the island remains relatively low. Therefore, the driving effect of the impurity contraction and radiation cooling can be ignorable, which is why the magnetic island given by the initially 2/1 perturbation remains almost unchanged. After $t=2000{t}_{\mathrm{A}}$, the impurity density at the O-point and the magnetic island width both increase simultaneously with time, while the impurity density at the X-point does not show a corresponding growth pattern relative to the island width. This observation suggests that the growth of the tearing mode is primarily driven by impurities within the magnetic island. The magnetic island triggered by the impurity radiation could serve as a seed island for the neoclassical tearing mode (NTM), which is a subject for further studies.

Figure  5.  The temporal evolutions of the impurity density at O-point and X-point of the 2/1 magnetic island (blue line), as well as the width of the 2/1 magnetic island (red ‘x’ symbol).

DownLoad: Full-Size Img PowerPoint

To study the formation of the 2/1 tearing mode, figure 6 presents the contour plot of the perturbed toroidal electric field, the Poincaré plot of magnetic field lines, the profile of the toroidal current density, and the q profile at $t=3005{t}_{\mathrm{A}}$, respectively. The 2/1 tearing mode is difficult to distinguish from the mode structure in figure 6(a), but it can be identified through the Poincaré plot of magnetic field lines in figure 6(b). It is evident that the 2/1 magnetic island is formed on the $q=2$ surface. Impurity radiative cooling can lead to temperature and plasma contraction. Plasma contraction results in the gradient of the plasma current, as shown in figure 6(c). When the gradient of the toroidal current density reaches the $q=2$ surface, the 2/1 tearing mode can be excited. At this point, the q profile flattens out on the $q=2$ surface, as shown in figure 6(d).

Figure  6.  (a) Contour plot of the perturbed toroidal electric field, (b) Poincaré plot of magnetic field lines, (c) profile of the toroidal current density and (d) the q profile at $t=3005{t}_{\mathrm{A}}$.

DownLoad: Full-Size Img PowerPoint

4.   The driving of the tearing mode by scanning impurity parameters

The initial profiles of the safety factor and the plasma pressure used in this section are shown in figure 7. The stability parameter of the tearing mode is ${\mathrm{\Delta }}'_{0} > 0$, which means the $m/n=2/1$ tearing mode is unstable in this equilibrium. We scan the impurity source concentration ${S} _{0}$ and the impurity atomic number Z to study the relationship between the impurity cooling effect and the growth rate of the tearing mode. The parameters in the simulations are ${D}_{\perp }=1\times {10}^{-6}$, ${D}_{\parallel }= 1\times {10}^{-3}$, ${\kappa }_{\perp }=1\times {10}^{-6}$, ${\kappa }_{\parallel }=1$, $\nu =1\times {10}^{-6}$, and $\eta \propto {Z}_{\mathrm{e}\mathrm{f}\mathrm{f}}{T}_{\mathrm{e}}^{-3/2}$. Random perturbations are applied at the beginning of these simulations.

Figure  7.  The initial equilibrium profiles of (a) safety factor and (b) plasma pressure.

DownLoad: Full-Size Img PowerPoint

We choose a typical case from these simulations to illustrate the behaviors of the impurity. In this case, we have ${S} _{0}=5\times {10}^{-5}$ and $Z=6$. The evolutions of the kinetic energy and the growth rates for different toroidal modes are shown in figure 8. In this case, the linear growth rate for $m/n=2/1$ is ${\gamma }_{}=0.0075$ that is slightly larger than ${\gamma }_{}=0.0061$ for the case without the impurity. This is because the temperature perturbation caused by impurity radiation contracts near the $q=2$ surface during the linear stage (figure 9(c)). Although the peak of the temperature perturbation is not exactly at the $q=2$ surface, it still drives the growth of the 2/1 mode. In this case, only the 2/1 tearing mode is linearly unstable. While the $n > 1$ modes exhibit larger growth rates during the linear stage, they are identified as beat modes since ${\gamma }_{n=2}={2\gamma }_{n=1}$, ${\gamma }_{n=3}={3\gamma }_{n=1}$, and ${\gamma }_{n=4}={4\gamma }_{n=1}$. Impurities can flow into the interior of the magnetic island, causing the impurity radiative cooling outside the $q=2$ surface to remain weak, which prevents the excitation of the ballooning mode.

Figure  8.  Evolutions of (a) the kinetic energy and (b) the growth rates for different toroidal modes.

DownLoad: Full-Size Img PowerPoint

The contour plots of the total impurity density, the perturbed temperature, the perturbed toroidal electric field and impurity radiative power at two typical moments are shown in figure 9. At $t=1932{t}_{\mathrm{A}}$, the region of the maximum impurity is close to the $q=2$ surface where a large temperature gradient is built-up due to impurity radiative cooling as shown in figures 9(a) and (c). At this moment, the modes with $n > 1$ just starts to beat-up, so figure 9(e) primarily shows the evident 2/1 tearing mode. At $t=3650{t}_{\mathrm{A}}$, during the formation of the magnetic island, impurities are drawn into the island. The impurity density within the magnetic island is influenced by the poloidal distribution of impurities. As the impurity penetrates into the magnetic islands, the temperature cooling inside the magnetic island can increase the saturation width of the magnetic island. In this case, the saturation width of the 2/1 magnetic island is $w/a=0.165$ that is larger than $w/a=0.157$ without impurity.

Figure  9.  Contour plots of (a) the total impurity density, (c) the perturbed temperature, and (e) perturbed toroidal electric field at $t=1932{t}_{\mathrm{A}}$, and contour plots of (b) the total impurity density, (d) the perturbed temperature, and (f) the perturbed toroidal electric field at $t=3650{t}_{\mathrm{A}}$.

DownLoad: Full-Size Img PowerPoint

The linear growth rate and the saturation width of 2/1 magnetic island with scanned Z and ${S} _{0}$ are shown in figure 10. Although different Z have different radiation cooling curves, they show some commonalities in these cases. The parameter ${S} _{0}{\cdot Z}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{2}$ can be regarded as the impurity radiative loss rate in equation (7). Since the perturbed temperature is directly influenced by impurity cooling, ${S} _{0}{\cdot Z}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{2}$ and min(T) are linearly dependent under certain conditions, where min(T) is the minimum perturbed temperature. The relationship between the growth rate and min(T) is also linear (figure 10(b)), but the relationship between the width of the magnetic island and min(T) is not linear (figure 10(d)) due to the different time durations required for the magnetic island to saturate. There is a certain relationship between the saturation width of the magnetic island and the impurity parameters, as shown in figure 10(c). It is implied that it is possible to predict the saturation width of the magnetic island based on the initial values of Z and ${S} _{0}$. The saturation width of the magnetic island tends to converge as the minimum perturbed temperature decreases, as shown in figure 10(d). Generally, it is possible to consider $w/a=0.2$ as the criterion for disruption [26]. However, the saturation width of the magnetic island in these cases is significantly less than $w/a=0.2$. In these cases, the driving effect of impurity radiation on the magnetic island is weak.

Figure  10.  The linear growth rates of the 2/1 tearing mode with (a) ${S} _{0}{\cdot Z}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{2}$ and (b) min(T) under different Z, and the saturation widths of the 2/1 magnetic island with (c) ${S} _{0}{\cdot Z}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{2}$ and (d) min(T) under different Z. Here ${Z}_{\mathrm{e}\mathrm{f}\mathrm{f}}$ is the effective charges number when the magnetic island is saturated and min(T) is the minimum perturbed temperature.

DownLoad: Full-Size Img PowerPoint

5.   Discussion and summary

The excitation and driving of the 2/1 tearing mode by impurity radiative cooling are studied through three-dimensional toroidal MHD simulations in this study. The density limit is commonly believed to be associated with the growth of the tearing mode. Therefore, we investigate the role of impurity in the density limit by examining their effect on the growth of the tearing mode.

We first conduct MHD simulations with the impurity applied at the boundary region under an equilibrium condition where the tearing mode is stable. It is observed that impurity radiation-induced plasma contraction results in a current gradient that excites the tearing mode at the resonance surface. Therefore, the process of the impurity-excited tearing mode is verified by MHD simulation. This result indicates that the impurity-excited tearing mode is not constrained by the shape of the q profile. That could be the reason why the density limit is not dependent on the q profile. During the formation of the magnetic island, impurities enter the island from the X-point and accumulate at the O-point, where they undergo radiative cooling that drives the growth of the magnetic island.

Impurity radiative cooling can drive both the linear growth rate of the tearing mode and the increase in magnetic island width. During the linear stage, the temperature cooling contracts towards the resonance surface, leading to increase the plasma pressure gradient and current and consequently enhance the linear growth rate of the tearing mode. When the magnetic island forms and it entrains nearby impurities into the island, the impurity radiation cooling further enhances the growth and the saturation level of the magnetic island. There is a correlation between different initial impurity parameters (Z, ${S} _{0}$) and the saturation width of the magnetic island. Since there are varying radiative cooling curves for different Z values, this correlation may only apply in the cases of low beta. As the impurity radiation increases, the saturation width of the magnetic island tends to converge. In these cases, the maximum saturation width of the magnetic island is considerably less than the criterion for disruption. We note that the island growth rate can be calculated using the modified Rutherford equation [36]

where ${\Delta }'_{\mathrm{c}\mathrm{l}\mathrm{a}\mathrm{s}\mathrm{s}\mathrm{i}\mathrm{c}}\left(w\right)$ is the classical term determined by the shear of q, ${\Delta }'_{\mathrm{r}\mathrm{a}\mathrm{d}}\left(w\right)$ is the current perturbation term determined by radiation, and ${\Delta }'_{\mathrm{A}}\left(w\right)$ is the island asymmetry term. Since the initial equilibrium in these cases is unstable for the tearing mode, there is a strong gradient in the current profile. The effect of high-Z impurities accumulating towards the core is not considered in our MHD model, resulting in relatively weak impurity radiative cooling inside the magnetic island. The q profile and the current profile undergo only minor changes in these MHD simulations, so the saturation width of the magnetic island does not change significantly. These results imply that the saturation width of the magnetic island is mainly determined by classical effects under the equilibrium of the unstable tearing mode. The influence of boundary impurity radiation contraction on magnetic island growth is weak in this model.