# Evaluation of ideal MHD mode stability of CFETR baseline scenario

###### Abstract

The CFETR baseline scenario is based on a H-mode equilibrium with high pedestal and highly peaked edge bootstrap current, along with strong reverse shear in safety factor profile. The stability of ideal MHD modes for the CFETR baseline scenario has been evaluated using NIMROD and AEGIS codes. The toroidal mode numbers (n=1-10) are considered in this analysis for different positions of perfectly conducting wall in order to estimate the ideal wall effect on the stability of ideal MHD modes for physics and engineering designs of CFETR. Although, the modes (n=1-10) are found to be unstable in ideal MHD, the structure of all modes is edge localized. Growth rates of all modes are found to be increasing initially with wall position before they reach ideal wall saturation limit (no wall limit). No global core modes are found to be dominantly unstable in our analysis. The design of and strong reverse shear in profile is expected to prevent the excitation of global modes. Therefore, this baseline scenario is considered to be suitable for supporting long time steady state discharge in context of ideal MHD physics, if ELMs could be controlled.

## 1 Introduction

Besides being a partner in ITER [1], China has recently proposed to design and potentially build China Fusion Engineering Test Reactor (CFETR) [2]. The goal is to address the physics and engineering issues essential for bridging the gap between ITER and DEMO, and to promote the advancement towards fusion reactor. These issues include efficient breeding of tritium after capturing high energy neutrons into lithium blanket and exploring option for DEMO blanket and divertor solutions. CFETR is expected to achieve high annual duty factor of and demonstrate tritium self-sufficiency with target tritium breeding rate greater than [3, 4].

A conceptual engineering design of CFETR including different coils and remote maintenance systems was prepared in the beginning [5]. The initial parameters of CFETR was set up through running a 0-D system code [2], and later these are optimized involving different other system codes (GASC and TESC) [3]. The preliminary design of snowflake divertor for CFETR has been made, and simulation work is carried out to evaluate heat flux onto the divertor [6]. To fulfill different physics goals, the CFETR has been designed for two steady-state scenarios - baseline and advanced scenarios. Baseline CFETR scenario is designed to achieve moderate fusion power (200 MW) applying a fully non-inductive current drive, giving more importance towards challenging annual duty factor . So, the idea is to achieve these targets with a conservative stable physics scenario first, before finally moving to the advanced scenario. Advanced design is aimed at higher fusion power and gain close to fusion reactor with challenging fraction of non-inductive bootstrap current drive. A detailed comparison between these two scenarios using different system codes analysis has been reported in a recent article [4]

Due to the goal of achieving high and high fraction of non-inductive bootstrap current in CFETR, both pressure and current driven instabilities are likely to threaten steady state operation. To confirm the viability of long duration steady state operation in CFETR scenarios, a thorough evaluation of the stability of all ideal MHD modes is essential, so that a stable parameter space could be determined. The strong reverse shear in safety factor profile and the optimized design of are expected to stabilize different devastating global core modes, such as and internal kink modes. The requirement of moderate to high fusion power gain in CFETR, would require higher pedestal top pressure value resulting in a steeper gradient in edge pressure profile. The aim for fully non-inductive operation, has proposed requirement of and of bootstrap current fraction to baseline and advanced scenarios respectively, whereas the ITER steady state is designed to be (see Table-1 of [7]). These requirements lead to high pedestal and peaked edge current, which are expected to drive the excitation of edge localized modes (ELMs). The repetitive expulsion of stored plasma energy and particles due to ELMs, would degrade plasma confinement and damage divertor and first wall components. For reactor scale machines, the sizes of ELMs are projected to be larger than those in current tokamaks [8, 9]. Thus, stability analysis of ELMs is essential for further evaluating and optimizing the design of CFETR baseline scenario.

The present article reports the stability analysis of the ideal MHD modes for CFETR baseline scenarios using the initial value extended-MHD code NIMROD [10] and the eigen-value code AEGIS [11]. In the ideal MHD model, the stability of modes are evaluated using both NIMROD and AEGIS codes, and the growth rates are compared. Also, in another calculation, we have used Spitzer resistivity to represent realistic resistive regimes of CFETR scenarios. The effect of conducting wall on the growth rates of modes has been studied with different positions of CFETR wall. The objective is to find the no wall limit of ideal mode growth rates and to provide physics base for the engineering design on the optimal choice of wall position of CFETR.

The rest of this paper is organized as follows. In the second section, the equilibrium profiles of baseline scenario are introduced. In the third section, the resistive single-fluid MHD model in the NIMROD code is described. In the first subsection of the fourth section, the ideal MHD results from NIMROD is described with the benchmarking between NIMROD and AEGIS codes. In the second subsection of the fourth section, the influence of Spitzer profile on the stability of ideal MHD modes are shown. Finally, the main points are summarized and the conclusion is drawn.

## 2 Equilibrium of CFETR Baseline Scenario

We consider the equilibrium of CFETR baseline scenario in our calculation. The necessary physics and engineering parameters of this scenario was first set up through 0-D system code analysis. Then, this equilibrium has self-consistently been generated through multi-dimensional integrated modeling in OMFIT framework using the auxiliary heating source in a combination of electron cyclotron wave and neutral beam injection [7]. The plasma size is slightly smaller than ITER, with a major radius of m and a minor radius of m. The toroidal magnetic field (T) and the plasma current (10 MA) at magnetic axis are listed in Table 1 of reference [7], among others. Since the baseline case is not designed for demonstrating high fusion gain, the normalized is set to be , well below the no-wall limit where is the plasma inductance. This is expected to help this equilibrium to lie within stability limits of ideal MHD global modes. The plasma profiles of electron number density, ion temperature, safety factor and current density are shown as functions of square root of the normalized poloidal magnetic flux. Both density (Fig. 1a) and temperature (Fig. 1b) profiles show an edge pedestal region inside LCFS. Safety factor (q) profile has strong reverse shear region (Fig. 1c) and with low core current in order to avoid sawtooth crash. The current density profile has highly peaked edge current due to high fraction of bootstrap current (Fig. 1d).

## 3 Magneto-hydrodynamic (MHD) Model in NIMROD and AEGIS

The MHD equations used in our NIMROD calculations are:

(1) | |||

(2) | |||

(3) | |||

(4) | |||

(5) |

where u is the center-of-mass flow velocity with particle density and ion mass , is the combined pressure of electron () and ion (), represents resistivity, and is the ion viscous stress tensor. The initial value NIMROD code has been consistently used in studying different macroscopic phenomena in both fusion and space plasmas [12, 13, 14].

The AEGIS code solves ideal MHD eigen-value equation employing adaptive shooting method along radial direction and Fourier decomposition in poloidal and toroidal direction. This code has been efficiently used before in evaluating stability of low- modes in presence of both conducting and resistive walls [11, 15]. In AEGIS, ideal MHD formalism has been used to evaluate linear stability of toroidal modes , where the plasma region within separatrix is modeled to have zero resistivity, and the vacuum region extended from separatrix to conducting wall, does not contain any plasma or current, amounting to infinite resistivity. On the contrary, NIMROD uses the resistive MHD model for both the hot core plasma within separatrix and the low density, low temperature plasma of vacuum-like halo region between separatrix and conducting wall. So, for the purpose of comparison with the ideal MHD results, a hyperbolic tangent resistiy profile is adopted in NIMROD to represent the lowly resistive core plasma and highly resistive vacuum region. Employing this resistivity model, a comparison in growth rates is drawn between NIMROD and AEGIS results for the modes.

## 4 Results of Ideal MHD Stability Analysis

### 4.1 Step function profile of resistivity

#### 4.1.1 Ideal MHD stability analysis in NIMROD

The dimensionless parameter to model the perfectly conducting ideal core plasma and infinitely resistive vacuum-like region is the Lundquist number defined as , where resistive diffusion time with being the permeability of free space, resistivity, minor radius and Alfven time with being the radius of the magnetic axis, and the values of magnetic field and mass density at magnetic axis respectively. The profile of Lundquist number (inverse of resistivity) is specified as a function of the normalized poloidal flux with step-like hyperbolic tangent form shown in Fig. 2, where the Lundquist numbers in plasma and vacuum regions are denoted as and respectively. Following the same procedure described in earlier references [12, 16], was scanned to find its value in the ideal MHD regime. The value of is set to be and then growth rates of toroidal modes are calculated with conducting wall located at , where is the minor radius of plasma. As shown in Fig. 3a, the growth rate (normalized in Alfven time) of toroidal mode increases with mode number , with the fastest growing one being .

Then growth rates of all these modes are calculated after varying the position of conducting wall in a wide range starting from close to LCFS to a wall distance from magnetic axis at . The growth rates of reach the no-wall limit at close proximity of LCFS, which indicates that the ideal MHD growth rate for this baseline case does not depend much on the conducting wall position (Fig. 3b). Perturbed pressure and radial magnetic field for are edge localized at the edge pedestal near separatrix as shown in Fig. 4 (the location of separatrix is indicated by black lines of poloidal flux contour). From Fig. 4b, the poloidal mode structure has poloidal mode number for , and thus the rational surface can be identified as which locates at the pedestal. An apparent difference in mode structure between the mode and the mode is noticeable from Fig. 4, where the mode has broader radial structure than .

#### 4.1.2 Ideal MHD stability analysis in AEGIS and comparison with NIMROD

A comparison between NIMROD and AEGIS results is performed for modes . The ideal MHD growth rates of modes have been evaluated using AEGIS code for the same equilibrium discussed in Section 2. The comparison of normalized growth rates is shown in the Fig. 5 for two different wall locations at and . It is clear that has good agreement in growth rate between NIMROD and AEGIS. Modes have slight differences in growth rates between these codes. Perturbed radial displacements of mode calculated in AEGIS are plotted in Fig. 6 for different poloidal harmonics. Both real and imaginary part of these eigenfunctions have one harmonic to be external kink mode and others may have internal mode structures peaked around rational surfaces.

### 4.2 Spitzer model profile of resistivity

#### 4.2.1 Stabilizing role of resistivity profile

The stability of modes has been re-calculated after considering the Spitzer resisitivity profile that is , where , , denote the electron temperature, resistivity at magnetic axis, and the electron temperature profile respectively. Now, our equilibrium configuration has a resisitivity profile covering whole simulation domain depending on the radial profile of electron temperature. The inclusion of resistivity profile is expected to make the numerical modeling more accurate for predicting the stability of CFETR baseline design. In a recent article, resistivity has been reported to have both stabilizing and destabilizing effects on ideal MHD edge localized modes [17]. The linear calculation in NIMROD has checked the stability of modes after placing conducting wall at and found modes to be unstable, where the mode is always stable (blue curve in Fig. 7b). Presence of Spitzer resistivity profile leads to lower the growth rates of modes and the stabilization of the mode as compared to the results shown in Fig. 3a using the hyperbolic tangent profile. This result is consistent with earlier studies using NIMROD in the context of other tokamaks equilibria such as NSTX and JT-60U [13, 16].

#### 4.2.2 Growth rate variation with wall position and shape

The effect of conducting wall position on growth rate of all modes has been evaluated after considering self-similar wall configuration and Spitzer resistivity profile. The wall position has been varied in the calculation until the no wall limit of growth rate is reached. The normalized growth rates of modes are plotted in the Fig. 7a with wall position changing from close to separatrix to . The stable position of conducting wall is found to be at , a little away from plasma boundary. No mode is found unstable inside this position of wall. The growth rates of all modes vary with wall position in a similar way. Initially, it increases rapidly until the wall position is reached. Afterwards, it gradually approaches the no-wall limit value. While the wall positions for all modes transitioning to no-wall limit are basically same, the no wall limit growth rate increases monotonically with mode numbers from to .

The results in previous paragraph are calculated for self-similar wall. The growth rate calculations of different modes have also been carried out after considering recently proposed real shaped wall configuration of CFETR. The present wall position is near to the wall location of , but shape is different from regular self-similar wall. A clear stabilizing effect of real shape of wall is found as compared to the self-similar wall at (Fig. 7b). The growth rates of for two different wall positions with self-similar wall shape are plotted together with proposed wall shape with using same Spitzer resistivity profile for all three cases. High- growth rates are close to those with self-similar wall at , whereas low- rates are similar to the self-similar wall at .

#### 4.2.3 Density profile vs. uniform density

The influence of non-uniform pedestal density profile on the stability of edge modes has been studied. Density pedestal has driven the edge localized modes more unstable, as overall growth rate of all modes increases higher than the uniform density case (Fig. 8). The growth rates of modes are nearly same for both density cases but more different for modes. The higher the toroidal mode number is, the stronger is the influence of density pedestal on growth rate. Here, level of uniform density is kept same as the value of density profile at magnetic axis, therefore the normalizing Alfven time scale () is same for both density cases.

#### 4.2.4 Mode structures

The detailed structure of modes are shown in contour plots of Fig. 9a-d for self-similar wall located at , and also these are shown in Fig. 10a-d for the proposed wall geometry of CFETR. The perturbed pressure and radial component of magnetic field quantities are plotted in 2D (R-Z) plane for modes . All unstable modes have radial structure only localized at the edge pedestal region which predicts them to be of peeling-ballooning types. The location of all modes is close to the inside of separatrix which is indicated by black lines of poloidal flux contour. The positions and shapes of mode structure of these two different wall configurations remain same in (Figs. ). The spatial structure of the mode is more radially localized than that of the mode.

#### 4.2.5 Convergence test

A thorough convergence has been checked for radial and poloidal grid numbers, time step () and polynomial degree of finite element basis used in NIMROD calculation. The growth rates of modes remain almost same for poloidal grid number range (Fig. 10a) and radial grid number range (Fig. 10b). From time step to (Fig. 10c) the variation in growth rate remains within . Although there is moderate difference in growth between polynomial degree and for mode , but polynomial degrees and have almost same growth rates (Fig. 10d). These results show a good numerical convergence in our calculation.

## 5 Summary and Discussions

In summary, our analysis on the linear stability of CFETR baseline equilibrium, finds the excitation of edge localized modes at the pedestal region but no global modes are found to be dominantly unstable. Two different resisitivity models have been employed in the calculation, namely the hyperbolic tangent profile and the Spitzer resistivity profile. The growth rates of have been separately calculated and compared for these resisitivity models. In the ideal MHD model using hyperbolic tangent resistivity profile, modes are found to be unstable with edge localized mode structure.

The effect of conducting wall position on the stability of ideal MHD modes have been evaluated. A noticeable difference is found between the results from two resistivity profiles. In Spitzer resistivity profile case, all modes become stabilized before wall position but for hyperbolic tangent profile, all modes remain unstable even if the wall is placed at plasma boundary.

On basis of our analysis, the baseline scenario of CFETR equilibrium is not expected to dominantly unstable to global ideal MHD modes. This might help to avoid disruption event caused by such ideal MHD instabilities. But, due to steep pedestal gradient and peaked edge current, this scenario can be susceptible to the medium to large size ELMs.

This present calculation draws an overall picture of unstable linear ideal MHD modes with perfectly conducting wall in the CFETR baseline scenario, which are dominantly edge-localized modes in nature. To achieve long duration steady state operation maintaining fixed , efficient methods need to be investigated for controlling ELMs. Further characteristics of ELMs need to be determined from nonlinear simulation. In addition, the effect of toroidal flow on ELMs in this CFETR baseline scenario is planned to be examined as another potential element for changing ELM characteristics.

## 6 References

## References

- [1] Aymer R et al. 2002 Plasma Phys. Control. Fusion 44 519
- [2] Wan B et al. 2014 IEEE Trans. Plasma Sci. 42 495
- [3] Chan V S et al. 2015 Nucl. Fusion 55 023017
- [4] Shi N et al. 2016 Fusion Eng. Des. 112 47
- [5] Song Y T et al. 2014 IEEE Trans. Plasma Sci. 42 503
- [6] Mao S F et al. 2015 J. Nucl. Mat. 463 1233
- [7] Jian X et al. 2017 Nucl. Fusion 57 046012
- [8] Loarte A et al. 2007 Progress in the ITER Physics Basis Chapter 4: Power and particle control Nucl. Fusion 47 S203
- [9] Leonard A W 2014 Phys. Plasmas 21 090501
- [10] Sovinec C R, Glasser A and the NIMROD Team 2004 J. Comp. Phys. 195 355
- [11] Zheng L J et al. 2006 J. Comput. Phys. 211 748
- [12] Burke B J et al. 2010 Phys. Plasmas 17 032103
- [13] King J R et al. 2016 Phys. Plasmas 23 062123
- [14] Zhu P and Raeder J 2013 Phys. Rev. Lett. 110 235005
- [15] Zheng L J et al. 2005 Phys. Rev. Lett. 95 255003
- [16] Banerjee D et al. 2017 Nucl. Fusion 57 076005
- [17] Banerjee D et al. 2017 Phys. Plasmas 24 054501