MHD Analysis on the Physics Design of CFETR Baseline Scenarios

Abstract

The China Fusion Engineering Test Reactor (CFETR), currently under intensive physics and engineering designs in China, is a major project representative of the low-density steady-state pathway to the controlled fusion energy. One of the primary tasks of the physics design for CFETR is the assessment and analysis of the magnetohydrodynamic (MHD) stability of the proposed design schemes. Comprehensive MHD stability assessment of the CFETR baseline scenarios have led to preliminary progress that may further benefit engineering designs. For CFETR, the electron cyclotron current drive (ECCD) power and current required for the full stabilization of neoclassical tearing mode (NTM) have been predicted in this work, as well as the corresponding controlled magnetic island width. A thorough investigation on resistive wall mode (RWM) stability for CFETR is performed. For 80% of the steady state operation scenarios, active control methods may be required for RWM stabilization. The process of disruption mitigation with massive neon injection on CFETR is simulated. The time scale of and consequences of plasma disruption on CFETR are estimated, which are found equivalent to International Thermonuclear Experimental Reactor (ITER). Major MHD instabilities such as NTM and RWM remain challenge to steady state tokamak operation. On this basis, next steps on CFETR MHD study are planned on NTM, RWM, and shattered pellet injection (SPI) disruption mitigation.

Appendices

Appendix 1 MARS Code

MARS-F code [21] is based on the single fluid, linearized resistive MHD model,

$$\begin{aligned}&(\gamma +in\Omega )\varvec{\xi }={\varvec{v}}+(\varvec{\xi }\cdot \varvec{\nabla })R^2\nabla \phi , \end{aligned}$$

(6)

$$\begin{aligned}\rho (\gamma +in\Omega ){\varvec{v}}&= -\varvec{\nabla }p+{\varvec{j}}\times {\varvec{B}}+{\varvec{J}}\times {\varvec{b}} \nonumber \\&\quad -\rho [2\Omega {\hat{Z}}\times {\varvec{v}}+({\varvec{v}}\cdot \varvec{\nabla }\Omega )R^2\varvec{\nabla }\phi ]\nonumber \\&\quad -\varvec{\nabla }\cdot \varvec{\Pi }, \end{aligned}$$

(7)

$$\begin{aligned}&(\gamma +in\Omega ){\varvec{b}}=\varvec{\nabla }\times ({\varvec{v}}\times {\varvec{B}}-\eta {\varvec{j}})+({\varvec{j}}\cdot \varvec{\nabla }\Omega )R^2\varvec{\nabla \phi }, \end{aligned}$$

(8)

$$\begin{aligned}&(\gamma +in\Omega ){\varvec{p}}=-{\varvec{v}}\cdot \varvec{\nabla }P-\Gamma {P}\varvec{\nabla }\cdot {\varvec{v}}, \end{aligned}$$

(9)

$$\begin{aligned}&{\varvec{j}}=\varvec{\nabla }\times {\varvec{b}}, \end{aligned}$$

(10)

where R and $\phi$ are the plasma major radius and geometric toroidal angle, and ${\hat{Z}}$ is unit vectors along the vertical direction in the poloidal plane, respectively. The variables $(\varvec{\xi },{\varvec{v}},{\varvec{j}},{\varvec{b}},\rho ,p)$ represent the plasma perturbed displacement, velocity, current, magnetic field, density and pressure, respectively. The corresponding equilibrium quantities are denoted by $({\varvec{J}},{\varvec{B}},P)$. $\Omega$ is the angular frequency of the plasma flow along the toroidal angle, and n is the toroidal harmonic number. $\varvec{\Pi }$ is a viscous stress tensor, which is associated with the viscous force damping, such as the parallel sound wave damping.

Although the growth rate of the RWM is very slow, it eventually sets the upper limit on plasma pressure for the long pulse or steady-state advanced tokamak operations. MARS code is designed to compute the growth rate of the RWM and how to control it by the passive (the plasma rotation and drift kinetic resonances) and active (the feedback control system) methods.

MARS code has been benchmarked and extensively applied to model RWM and compare with the experimental observation. For examples, a kinetic version of MARS found low-rotation threshold when applied to model a DIII-D discharge with balanced beam injection, agreeing with experimental observations [46]. MARS-F and its coupling to CARIDDI (CarMa) found quantitative agreement between the computed RWM growth rate and the experiments in RFX [47]. MARS-F has also modeled the resonant field amplification for a series of JET plasmas which agrees with experimental measurements [48]. Reference [49] has compared the unstable RWM regime obtained using MARS-K with that in DIII-D experiments, revealing the impact of energetic particle losses and toroidal rotation drop in destabilizing the mode. Finally, the MARS-K modeling of stable RWM induced resonant field amplification quantitatively agrees with DIII-D experiments [50].

Appendix 2 AEGIS Code

The Adaptive Eigenfunction Independent Solution shooting (AEGIS) code employs the adaptive shooting method in the radial direction and Fourier decomposition in the poloidal direction [22]. Therefore, the AEGIS code has high resolution near the singular surfaces for the study of MHD instabilities. The AEGIS code has been used to study the linear behaviors of RWMs in ITER and the earlier smaller-sized design of CFETR scenarios [51, 52].

The following perpendicular MHD equation was solved in AEGIS,

$$\begin{aligned} -\rho _m(\omega +n\Omega +i\gamma _p)^2 \varvec{\xi _{\perp }} = \delta {\varvec{J}}\times {\varvec{B}}+\delta {\varvec{B}}\times {\varvec{J}}-\varvec{\nabla }\delta P, \end{aligned}$$

where $\rho _m$ is the total apparent mass density, $\omega$ the mode frequency, n the toroidal mode number, $\Omega$ the toroidal rotation frequency, $\gamma _p$ is a small parameter used to heal the numerical singularity while calculating the Alfvén damping, $\xi$ is the fluid displacement, with subscript $\perp$ denoting the perpendicular component to the magnetic field, and ${\varvec{J}}$, ${\varvec{B}}$, and P are the equilibrium current density, magnetic field, and plasma pressure, respectively.

Appendix 3 MD Code

The reduced MHD model implemented in the MD code is given as follows [53]

$$\begin{aligned} \frac{\partial \psi }{\partial t} &= {} [\psi ,\phi ]-\partial _z\phi -S^{-1}_A (j-j_b-j_d)+E_{z0}, \end{aligned}$$

(11)

$$\begin{aligned} \frac{\partial u}{\partial t}&= {} [u,\phi ]+[j,\psi ]+\partial _z j+R^{-1}\nabla ^2_\perp u, \end{aligned}$$

(12)

$$\begin{aligned} \frac{\partial p}{\partial t}&= {} [p,\phi ]+\chi _\parallel \nabla ^2_\parallel p + \chi _\perp \nabla ^2_\perp p +S_0 , \end{aligned}$$

(13)

where $\psi$ and $\phi$ are the magnetic flux and electrostatic potential, $j=-\nabla ^2_\perp \psi$ and $u=\nabla ^2_\perp \phi$ are the current density and vorticity in the axial direction, respectively. The bootstrap current density is proportional to the pressure gradient as in $j_b=-f\frac{\sqrt{\varepsilon }}{B_\theta }\frac{\partial p}{\partial r}$, with f measuring the strength of bootstrap current fraction, which is defined as $f_b=\int _0^a j_brdr/\int _0^a j_zrdr$. $S_A=\tau _\eta /\tau _A$ and $R=\tau _\nu /\tau _A$ are the magnetic Reynolds number and kinematic Reynolds number, respectively, where $\tau _\eta =a^2\mu _0/\eta$, $\tau _\nu =a^2/\nu$ and $\tau _A=\sqrt{\mu _0\rho }a/B_0$ are the resistive diffusion time, the viscous diffusion time, and the Alfvén time, respectively. $\chi _\parallel$ and $\chi _\perp$ are the parallel and perpendicular transport coefficients. The source terms $E_{z0}=S_A^{-1}(j_0-j_{b0})$ and $S_0=-\chi _\perp \nabla ^2_\perp p_0$ in Eqs. (11) and (13) are chosen to balance the diffusion of equilibrium Ohm current and pressure, respectively. The length, time and velocity are normalized by the plasma minor radius a, Alfvén time $\tau _A$ and Alfvén velocity $V_A=B_0/\sqrt{\mu _0\rho }$ respectively. The Poisson bracket is defined as $[f,g]={\hat{z}}\cdot \nabla f\times \nabla g$.

Appendix 4 TM8 Code

The TM8 code has been used to model the physics of the ECCD stabilization of NTM [54], the drift-tearing modes [55], the double tearing modes [56], the mode coupling [57], the stochastic field [58], the resonant magnetic perturbation [59], and the error field [60]. The corresponding simulation results compare well with experiments.

The reduced MHD model implemented in the TM8 code includes the Ohm’s law, the plasma vorticity equation, and the plasma pressure evolution equation [61]

$$\begin{aligned}&\frac{\partial \psi }{\partial {t}}+{\varvec{v}}\cdot \varvec{\nabla }\psi =E-\eta (j_p-j_b-j_d), \end{aligned}$$

(14)

$$\begin{aligned}&\rho \left( \frac{\partial }{\partial t}+{\varvec{v}}\cdot \varvec{\nabla }\right) \nabla ^2\Theta ={\varvec{e}}_t\cdot (\varvec{\nabla }\psi \times \varvec{\nabla }j_p)+\rho \mu \nabla ^4\Theta , \end{aligned}$$

(15)

$$\begin{aligned}&\frac{3}{2}\left( \frac{\partial }{\partial t}+{\varvec{v}}\cdot \varvec{\nabla }\right) p=\varvec{\nabla }\cdot (\chi _{\parallel }\nabla _{\parallel }p)+\varvec{\nabla }\cdot (\chi _{\perp }\nabla {\perp }{p})+Q, \end{aligned}$$

(16)

where ${\varvec{v}}=\varvec{\nabla }\Theta \times {\varvec{e}}_t$, $\Theta$ is the stream function, ${\varvec{e}}_t$ the unit vector in the toroidal direction, and $j_p=-\nabla ^2\psi -2nB_{0t}/(mR)$, $j_b=-c_b\frac{\sqrt{\varepsilon }}{B_{\theta }}\frac{\partial p}{\partial r}$, and $j_d$ are the plasma current density, the bootstrap current density, and the current density driven by ECW in the ${\varvec{e}}_t$ direction, respectively.

Appendix 5 NIMROD/KPRAD Code

In the NIMROD code [62], the 3D extended MHD model is coupled with an atomic and radiation physics model from the KPRAD code [41, 42, 63], and the implemented equations for the coupled impurity-MHD model are as follows:

$$\begin{aligned}&\rho \frac{d\mathbf {V}}{dt} = - \nabla p + \mathbf {J} \times \mathbf {B} + \nabla \cdot (\rho \nu \nabla \mathbf {V}), \end{aligned}$$

(17)

$$\begin{aligned}&\frac{d n_e}{dt} + n_e \nabla \cdot \mathbf {V} = \nabla \cdot (D \nabla n_e) + S_{ion/rec}, \end{aligned}$$

(18)

$$\begin{aligned}&\frac{d n_i}{dt} + n_i \nabla \cdot \mathbf {V} = \nabla \cdot (D \nabla n_i) + S_{ion/3-body} , \end{aligned}$$

(19)

$$\begin{aligned}&\frac{d n_Z}{dt} + n_Z \nabla \cdot \mathbf {V} = \nabla \cdot (D \nabla n_Z) + S_{ion/rec} , \end{aligned}$$

(20)

$$\begin{aligned}&n_e \frac{d T_e}{dt} = (\gamma - 1)[n_e T_e \nabla \cdot \mathbf {V} + \nabla \cdot \mathbf {q_e} - Q_{loss}], \end{aligned}$$

(21)

$$\begin{aligned}&\mathbf {q}_e = -n_e[\kappa _{\parallel } {\hat{b}} {\hat{b}} + \kappa _{\perp } ({\mathscr {I}} - {\hat{b}} {\hat{b}})] \cdot \nabla T_e, \end{aligned}$$

(22)

$$\begin{aligned}&\mathbf {E} + \mathbf {V} \times \mathbf {B} = \eta \mathbf {j}. \end{aligned}$$

(23)

Here, $n_i$, $n_e$, and $n_Z$ are the main ion, electron, and impurity ion number density respectively, $\rho$, $\mathbf {V}$, $\mathbf {J}$, and p the plasma mass density, velocity, current density, and pressure respectively, $T_e$ and $\mathbf {q}_e$ the electron temperature and heat flux respectively, D, $\nu$, $\eta$, and $\kappa _{\parallel } (\kappa _{\perp })$ the plasma diffusivity, kinematic viscosity, resistivity, and parallel (perpendicular) thermal conductivity respectively, $\gamma$ the adiabatic index, $S_{ion/rec}$ the density source from ionization and recombination, $S_{ion/3{\text {-}}body}$ also includes contribution from 3-body recombination, $Q_{loss}$ the energy loss, $\mathbf {E} (\mathbf {B})$ the electric (magnetic) field, ${\hat{b}}=\mathbf {B}/B$, and ${\mathscr {I}}$ the unit dyadic tensor.

All particle species share a single temperature $T=T_e$ and fluid velocity $\mathbf {V}$, which assumes instant thermal equilibration among the main ions, the impurity ions, and the electrons. Pressure p and mass density $\rho$ in momentum equation (17) include contributions from the impurity species. Each charge state of impurity ion density is tracked in the KPRAD module and used to update the source/sink terms in the continuity equations due to ionization and recombination. Both convection and diffusion terms are included in each continuity equations where all the diffusivities are assumed same. Quasi-neutrality is maintained through the condition $n_e=n_i+\sum Z n_{z}$, where Z is the charge of impurity ion. The energy loss term $Q_{loss}$ in Eq. (21) is calculated from KPRAD module based on a coronal model, which includes contributions from bremsstrahlung, line radiation, ionization, recombination, Ohmic heating, and intrinsic impurity radiation. Anisotropic thermal conductivities are temperature dependent, i.e. $\kappa _{\parallel } \propto T^{5/2}$ and $\kappa _{\perp } \propto T^{-1/2}$. Similarly, the temperature dependence of resistivity $\eta$ is included through the Spitzer model.