⇦Chinese Physics Letters, 2021, Vol. 38, No. 4, Article code 045203 Energetic Particle Transport Prediction for CFETR Steady State Scenario Based on Critical Gradient Model Yunpeng Zou (邹云鹏)1, V. S. Chan (陈锡熊)2,3, Wei Chen (陈伟)1*, Yongqin Wang (王雍钦)1, Yumei Hou (侯玉梅)1, and Yiren Zhu (朱毅仁)1 Affiliations 1Southwestern Institute of Physics, Chengdu 610041, China 2General Atomics, San Diego, CA 92186, USA 3School of Nuclear Science and Technology, University of Science and Technology of China, Hefei 230026, China Received 1 December 2020; accepted 2 February 2021; published online 6 April 2021 Supported by the National Key R&D Program of China (Grant No. 2019TFE03020000), the National Natural Science Foundation of China (Grant Nos. 11875021, 12005054, and 12005055), the Sichuan Science and Technology Program (Grant No. 2020jqqn0070).
*Corresponding author. Email: chenw@swip.ac.cn Citation Text: Zou Y P, Chen X X, Chen W, Wang Y Q, and Hou Y M et al. 2021 Chin. Phys. Lett. 38 045203    Abstract The critical gradient mode (CGM) is employed to predict the energetic particle (EP) transport induced by the Alfvén eigenmode (AE). To improve the model, the normalized critical density gradient is set as an inverse proportional function of energetic particle density; consequently, the threshold evolves during EP transport. Moreover, in order to consider the EP orbit loss mechanism in CGM, ORBIT code is employed to calculate the EP loss cone in phase space. With these improvements, the AE enhances EPs radial transport, pushing the particles into the loss cone. The combination of the two mechanisms raises the lost fraction to 6.6%, which is higher than the linear superposition of the two mechanisms. However, the loss is still far lower than that observed in current experiments. Avoiding significant overlap between the AE unstable region and the loss cone is a key factor in minimizing EP loss. DOI:10.1088/0256-307X/38/4/045203 © 2021 Chinese Physics Society Article Text The China fusion engineering test reactor (CFETR)[1] is designed to be the next-generation fusion facility; one of its aims is self-sustainable burning state with a fusion gain of $Q\approx20$–30.[2] The behavior of energetic particles (EPs), particularly $\alpha$ particles generated by the fusion reaction, influences heating power, bootstrap current, helium ash exhausting, etc.[3] Unfortunately, the integrated simulation based on the OMFIT framework[4] does not include an EP model for EP instability and transport estimation. Instead, in this framework, the ONETWO code[5] is employed to roughly predict the $\alpha$ particle's radial profile. A number of experiments[6–8] have proved that there is an obvious distinction between the classical slowing down model and realistic distribution. The primary reason for this is that Alfvén eigenmodes (AEs) driven by EP could enhance EP transport conversely, in a manner resembling a feedback system. Currently, various simulation methods, such as MEGA,[9,10] are being developed to accurately calculate EP distribution. However, due to the huge computational effort and long simulation time required, these methods are inconvenient for CFETR design. Thus, an efficient method, critical gradient model (CGM),[11] is employed to estimate $\alpha$ particle redistribution due to AE instabilities. According to the current analytical theory,[12] AE is destabilized by the EP pressure gradient. Since the contribution of temperature gradient is slight[13] on threshold calculation, neglecting the effect of the temperature gradient, there should be a critical density gradient for marginal stable AE. The CGM therefore assumes that if the EP density gradient exceeds the threshold, the EP transport would be enhanced, and the diffusivity would increase significantly. The existence of the critical gradient and EP stiff transport is confirmed by the simulation[11,14] and has been experimentally verified.[15] GYRO code[16] was employed to calculate the local critical gradient, which was then input into ALPHA,[17] a 1D EP density transport code. Diffusivity was assumed to be a linear function of the distinction between EP gradient and threshold. This combination of GYRO and ALPHA was then used to estimate the EP redistribution.[18] Owing to the significant computational demands of GYRO, a parallelized wrapper, TGLFEP,[13] has been employed in its place for the purpose of efficient threshold scanning. In addition, the ALPHA code has been rewritten to EPtran[19] by expending EP distribution into $(\psi, E, v_{\parallel}/v)$ space. Thus, this work employs a combination of TGLFEP and EPtran to predict EP transport. The key issue for the CGM is critical gradient calculation. Since both GYRO and TGLFEP function as local stability analysis programs, the frequency and growth rate are calculated using local plasma parameters. For each flux surface, we select AEs from multiple numerical output based on mode structure and frequency, choosing the AE with maximum growth rate to calculate the critical gradient. In previous versions of the CGM, a parameter for density gradient was used to estimate AE marginal stability. We have reduced the density gradient until the AE growth rate is equal to zero, i.e., the critical density gradient is calculated to satisfy the condition of $\gamma_{\rm AE}=0$. During the scanning process, the characteristic length of EP density $[L_{n_{\rm EP}}=-n_{\rm EP}/(dn_{\rm EP}/dr)]$ is fixed. Clearly, the assumption of fixed characteristic length means that the profile shape could not change, which is inconsistent with the reality. Although a correction factor was added,[13] the profile was modified only slightly, with an additional artificial parameter required for the correction. Instead, we select the normalized density gradient for marginal stability estimation, and the normalized critical density gradient is set as a function of EP density. To relax this assumption, we resort to an analytic formula of AE growth rate,[20] $$\begin{alignat}{1} \frac{\omega_{\rm i}}{k_{\parallel}v_{\rm A}}={}&-q_0^2\Big[\frac{\beta_{\rm c}}{2}(G_{\rm mi}^{\rm T}+G_{\rm me}^{\rm T})\\ &+\beta_\alpha(G_{\rm s\alpha}^{\rm T}-nq_0\delta_{\alpha}H_{\rm sa}^{\rm T})\Big].~~ \tag {1} \end{alignat} $$ Here, only $\beta_\alpha$ and $\delta_\alpha$ [$\propto dp_\alpha/(p_{\alpha}dr)$] will be modified during EP transport. Since the drive term includes both pressure and pressure gradient, $L_{n_{\rm EP}}$ should be allowed to vary during EP profile evolution. For marginal stability, Eq. (1) can be rewritten as $$ \frac{a}{L_{n_{\rm EP}}^{\rm th}}=\frac{k_1}{n_{\rm EP}}+k_2,~~ \tag {2} $$ where $k_1$ represents the Landau damping of bulk plasma, $k_2$ is the Landau damping of the EP itself, and the temperature gradient drive of the EP. The TGLFEP code confirmed the inversely proportional relationship via parameter scanning (Fig. 1). The threshold calculated by TGLFEP are depicted as points, fitted by the dashed curve. Consequently, the coefficients $k_1$ and $k_2$ are obtained by the fitting curve, and then input into the transport code, EPtran. Note that the coefficients $k_1$ and $k_2$ require to be calculated for each surface. Meanwhile, the diffusivity is rewritten based on Eq. (9) in Ref. [19] as follows: $$ D_{\rm EP}^{\rm rr}=D_{\rm AE}\Big[\frac{a}{L_{n_{\rm EP}}}-\frac{a}{L_{n_{\rm EP}}^{\rm th}}\Big]_{>0}+D_{\rm ITG/TEM},~~ \tag {3} $$ where $[x]_{>0}=0$ if $x < 0. D_{\rm AE}$ is an artificial transport coefficient, set to be large enough to drive EP transport generally. The value affects the transported profile only slightly. Here ${a}/{L_{n_{\rm EP}}}$ represents the normalized density gradient during transport, and ${a}/{L_{n_{\rm EP}}^{\rm th}}$ is the normalized critical density gradient. $D_{\rm ITG/TEM}$ denotes the micro-turbulent transport of EP obtained via the Angioni model.[21] By virtue of this improved method of threshold calculation, the threshold of each surface can be modified during EP transport. Furthermore, the previous CGM only considered EP loss induced by diffusion at the last closed flux surface (LCFS); as a result, most of the EPs are deposited near the edge rather than being lost. Indeed, the non-uniformity of the magnetic field leads to drifts in the guiding center. Both passing particles and trapped particles may escape from the LCFS due to the wide orbit width. The behavior of the EP depends on location, energy, and pitch angle. To estimate whether or not a particle is lost, ORBIT code[22] is therefore employed to calculate the loss cone in phase space.
cpl-38-4-045203-fig1.png
Fig. 1. Normalized critical gradient for $\rho=0.5$ and $\rho=0.6$ as an inverse proportional function of EP density. The TGLFEP simulation results are depicted as points, fitted by the dashed curves.
We have now established a workflow to calculate EP distribution via multiple programming methods. EFIT code[23] is employed to construct equilibrium, which is then input into the other programs. TGLFEP code calculates the threshold for AE marginal stability. The mode structure is obtained via the MHD-kinetic code, NOVA-K,[24] and amplitude is calculated using MEGA[9] for AE saturation. The structure of the perturbation is then input into ORBIT for loss cone calculation. Finally, both threshold and loss cone are read using EPtran to obtain an EP redistribution prediction. With these optimizations of the CGM, $\alpha$ particle transport and redistribution is predicted for a CFETR steady-state scenario. CFETR is designed in a large size with a major radius ($R$) of 7.2 m, and minor radius ($a$) of 2.2 m. The magnetic field is $B=6.5$ T, and the plasma current is $I=12$ MA. The temperature and density on the axes are $T_{\rm i}(0)=23.5$ keV, $T_{\rm e}(0)=24.5$ keV, and $n_{\rm e}(0)=12.5\times10^{19}$ m$^{-3}$, respectively. Reverse magnetic shear is designed [see the red dashed curve in Fig. 2(a)] so that $q(0)=6.36$, $q_{{\min}}=3.55$, and $q(a)=9.33$. The isotope fraction is assumed to be $D:T=50:50$, so that the initial distribution of $\alpha$ particles is calculated by means of the fusion reaction, via the classical slowing down model:[25] $$ F_{\rm s}(v)=\frac{S_0\tau_{\rm s}}{4\pi}\frac{H(v_\alpha-v)}{v_{\rm c}^3+v^3},~~ \tag {4} $$ where $S_0$ is the source intensity, $\tau_{\rm s}$ is the slowing-down time, $v_{\rm c}$ is the crossover velocity, and $v_\alpha$ is the $\alpha$ particle birth speed (3.5 MeV). The stability of TAE is analyzed using NOVA-K. The unstable toroidal mode number is wide, due to large size of the device, and the $n=6$–10 TAEs are dominant. Owing to relatively flat $q$ profile within $\sqrt{\psi}=0.2$–0.8, multiple modes can be driven by a single toroidal mode number at different radial locations. The $n=5$ continuum and mode location are depicted in Fig. 2(a), for example. Two black solid curves depict the continuum, with the TAE gap in the middle of the two curves. The red dashed curve depicts the $q$ profile. Three $n=5$ TAEs are distinguished by different colors, and the mode structures are shown in the insets. The two inside modes (yellow and blue) have higher growth rates than the outside mode (orange). So many unstable modes could induce significant EP transport; however, we only consider the single toroidal mode number, which is limited by the simulation method. In this way, the dominant toroidal mode number is selected for each flux surface. The stability of TAE dominates the profile of the threshold. The initial threshold profile (no transport) is depicted by the red dashed curve in Fig. 2(b), and the normalized density gradient of classical slowing down distribution is depicted by the black solid curve. In the range from $\sqrt{\psi}=0.4$–0.8, exceeding the threshold indicates that the TAE is unstable. During the transport, the threshold increases in the core [$\sqrt{\psi} < 0.5$] and decreases near the edge [$\sqrt{\psi}>0.5$].
cpl-38-4-045203-fig2.png
Fig. 2. (a) Shear Alfvén continuum and mode location for $n=5$ in the TAE gap. The red dashed curve depicts the $q$ profile. Three TAEs are depicted using different colors, and the corresponding structures are shown in the insets. (b) Threshold for $\alpha$ particle in CFETR steady state scenario (red dashed curve). The classical slowing-down distribution is denoted by black solid curve for comparison.
The boundary of the loss cone is depicted in Fig. 3 as a function of $\sqrt{\psi}$ and pitch angle, and the EP energy is denoted by different colors. Here, the initial EP distribution is selected to be uniform, i.e., $(\psi, E, v_{\parallel}/v)$, at the outer mid-plane. The simulation indicates that particles with a negative pitch angle are more easily lost, and the loss boundary peaks at $v_{\parallel}/v =-0.6$. The prediction of loss cone in CFETR is smaller than recent experimental values[26] due to the high magnetic field, but is slightly larger than that given by ITER,[26] since the orbital width is proportional to $q$. In contrast to the monotonically increasing $q$ profile of the ITER scenario, a flat $q$ profile with high $q_{\min}$ induces serious orbit loss. In experimental terms, the loss mechanism could be enhanced by AE perturbation, and consequently the loss cone expands. However, in the CFETR scenario, the unstable AE region and loss cone are almost separated, so that the influence of AE on the loss cone is negligible. Here, the collision between $\alpha$ particle and background plasma is neglected, which modifies the energy and pitch angle of the $\alpha$ particle in motion.
cpl-38-4-045203-fig3.png
Fig. 3. The loss cone as a function of $(\psi, E, v_{\parallel}/v)$ for the CFETR steady-state scenario. EP energy is depicted using different colors. The points are obtained by ORBIT for loss boundary, and corresponding curves indicate linear fitting. The region on the right-hand side of the curve relates to loss, and the left-hand side refers to confinement.
The $\alpha$ particle deposition is denoted by a red curve in Fig. 4(a), taking into consideration both transport induced by AE, and orbit loss. For comparison, the classical slowing-down distribution is shown as a black curve, the prediction of the previous CGM is represented as a blue curve, and prediction of ONETWO is given as a purple curve, which is employed in the current OMFIT framework. The ONETWO code uses a more exact slowing down model, which includes impurity effects, the fusion reaction between beam deuterium and background tritium, and in which Coulomb collision is related to local electron density and temperature. The profiles where $\rho>0.5$ are enlarged in the inset. The slowing down distribution is noticeably higher than the two transported profiles in the core, indicating the obvious transport of $\alpha$ particles subsequent to generation. Evidently, the previous CGM without FOW effect induces $\alpha$ particle deposition near the edge, similarly to the ITER prediction.[27] Compared with the previous CGM, the improved CGM optimizes the transport profile in the unstable AE region and edge. Moreover, in Fig. 4(b), the contribution of AE (green dash curve) and FOW (yellow dash curve) effects can be shown separately, as well as being evaluated in combination. The classical slowing down and transported profile are depicted using the same colors as those in Fig. 4(a). From the insets, we note that a high magnetic field favors a reduction in the EP loss cone, i.e., good EP confinement, so that FOW-induced loss is only  $1\%$. AEs primarily change the shape of the EP profile, so that the loss fraction is $3.8\%$. In combination, unstable AEs push particles outwards into the loss cone region to enhance prompt loss, raising the loss fraction to $6.6\%$. Fortunately, the loss depicted by the CFETR scenario is much lower than the DIII–D discharge ($50\%$). We attribute this to the well separation between the AE unstable region and the loss cone region.
cpl-38-4-045203-fig4.png
Fig. 4. (a) Behavior of $\alpha$ particle deposition (red curve), considering unstable AE and orbit loss. For comparison, the classical slowing-down distribution is shown as a black curve, the prediction by the previous CGM is represented by a blue curve, and the purple curve depicts the ONETWO result, which is used in current integrated simulations, based on the OMFIT framework. (b) The contribution of AE (green dashed curve) and FOW (yellow dashed curve) mechanisms are depicted separately. Profiles with $\sqrt{\psi}> 0.5$ are depicted in the insets.
So far, we have only considered the spatial transport of $\alpha$ particles. Without energy diffusion, the number of particles with high energy in the current simulation is larger than that in reality. Since particles with higher energy are more easily lost due to their larger orbital width, the number of lost particles is overestimated. In addition, helium ash is not included in the model, i.e., all $\alpha$ particles are regarded as energetic particles. During the Coulomb collision between EP and background plasma, the EP continuously loses energy, and finally transfers to helium ash. The diffusion coefficients of EP and ash are quite different. The energy-dependent mechanism can be neglected in current experimental devices, since the injected energy of NBI is only $O(10)$ keV. However, it will become significant for 3.5 MeV $\alpha$ particles in future fusion devices. In summary, the $\alpha$ particle redistribution of the CFETR steady-state scenario is predicted by our improved CGM, which takes into account the fact that threshold evolution relieves the abnormally flat profile in the unstable AE region. The FOW effect adds an additional sink for EP loss. The simulation indicates that unstable AE enhances EP radial transport, pushing the particles outwards, so that the particles move into the loss cone and escape. Separate mechanisms only result in slight EP loss, but the combination of AE and FOW raises the loss to $6.6\%$. The lost fraction is still much lower than that measured in current experiments, benefitting from a well-separated unstable AE and loss cone region. References Overview of the present progress and activities on the CFETRProgress of the CFETR designEnergetic Particles in Magnetic Confinement Fusion PlasmasIntegrated modeling applications for tokamak experiments with OMFITVariational moment solutions to the Grad–Shafranov equationFast-ion transport by Alfvén eigenmodes above a critical gradient thresholdTransport of energetic ions due to sawteeth, Alfvén eigenmodes and microturbulenceNonlinear magnetohydrodynamic effects on Alfvén eigenmode evolution and zonal flow generationMulti-phase simulation of fast ion profile flattening due to Alfvén eigenmodes in a DIII-D experimentGyrokinetic simulations of mesoscale energetic particle-driven Alfvénic turbulent transport embedded in microturbulenceExcitation of the toroidicity‐induced shear Alfvén eigenmode by fusion alpha particles in an ignited tokamakAlfvén eigenmode stability and critical gradient energetic particle transport using the Trapped-Gyro-Landau-Fluid modelNonlinear verification of a linear critical gradient model for energetic particle transport by Alfvén eigenmodesPhase-space dependent critical gradient behavior of fast-ion transport due to Alfvén eigenmodesAn Eulerian gyrokinetic-Maxwell solverPrediction of the fusion alpha density profile in ITER from local marginal stability to Alfvén eigenmodesDevelopment and validation of a critical gradient energetic particle driven Alfven eigenmode transport model for DIII-D tilted neutral beam experimentsKinetic transport simulation of energetic particlesStability of Alfvén gap modes in burning plasmasGyrokinetic calculations of diffusive and convective transport of α particles with a slowing-down distribution functionHamiltonian guiding center drift orbit calculation for plasmas of arbitrary cross sectionReconstruction of current profile parameters and plasma shapes in tokamaksNOVA: A nonvariational code for solving the MHD stability of axisymmetric toroidal plasmasTurbulent transport of alpha particles in reactor plasmasAlfvén eigenmode stability and fast ion loss in DIII-D and ITER reversed magnetic shear plasmasPrediction of Alfvén eigenmode energetic particle transport in ITER scenarios with a critical gradient model
[1] Wan Y et al. 2017 Nucl. Fusion 57 102009
[2] Zhuang G et al. 2019 Nucl. Fusion 59 112010
[3] Chen W and Wang Z 2020 Chin. Phys. Lett. 37 125001
[4] Meneghini O et al. 2015 Nucl. Fusion 55 083008
[5] Lao L L et al. 1981 Phys. Fluids 24 1431
[6]Van Zeeland M A et al. 2011 Phys. Plasmas 5 056114
[7] Heidbrink W W et al. 2017 Phys. Plasmas 24 056109
[8] Pace D C et al. 2011 Nucl. Fusion 51 043012
[9] Todo Y, Berk H L and Breizman B N 2010 Nucl. Fusion 50 084016
[10] Todo Y et al. 2014 Nucl. Fusion 54 104012
[11] Bass E M and Waltz R E 2010 Phys. Plasmas 17 112319
[12] Fu G and Van Dam J W 1989 Phys. Fluids B: Plasma Phys. 1 1949
[13] He S, Waltz R E and Staebler G M 2017 Phys. Plasmas 24 072305
[14] Bass E M and Waltz R E 2017 Phys. Plasmas 24 122302
[15] Collins C S et al. 2017 Nucl. Fusion 57 086005
[16] Candy J and Waltz R E 2003 J. Comput. Phys. 186 545
[17] Waltz R E and Bass E M 2014 Nucl. Fusion 54 104006
[18] Waltz R E et al. 2015 Nucl. Fusion 55 123012
[19] He S and Waltz R E 2016 Nucl. Fusion 56 056004
[20] Betti R and Freidberg J P 1992 Phys. Fluids B: Plasma Phys. 4 1465
[21] Angioni C and Peeters A G 2008 Phys. Plasmas 15 052307
[22] White R B and Chance M S 1984 Phys. Fluids 27 2455
[23] Lao L L et al. 1985 Nucl. Fusion 25 1611
[24] Cheng C Z and Chance M S 1987 J. Comput. Phys. 71 124
[25] Estrada-Mila C, Candy J and Waltz R E 2006 Phys. Plasmas 13 112303
[26] Van Zeeland M A et al. 2012 Nucl. Fusion 52 094023
[27] Bass E M and Waltz R E 2020 Nucl. Fusion 60 016032