⇦Chinese Physics Letters, 2021, Vol. 38, No. 5, Article code 055201 PTC: Full and Drift Particle Orbit Tracing Code for $\alpha$ Particles in Tokamak Plasmas Feng Wang (王丰), Rui Zhao (赵睿), Zheng-Xiong Wang (王正汹)*, Yue Zhang (张悦), Zhan-Hong Lin (林展宏), Shi-Jie Liu (刘仕洁), and CFETR Team Affiliations Key Laboratory of Materials Modification by Laser, Ion and Electron Beams (Ministry of Education), School of Physics, Dalian University of Technology, Dalian 116024, China Received 3 December 2020; accepted 5 March 2021; published online 2 May 2021 Supported by the National Natural Science Foundation of China (Grant Nos. 11975068 and 11925501), and the National Key Research and Development Program of China (Grant No. 2017YFE0300501).
*Corresponding author. Email: zxwang@dlut.edu.cn Citation Text: Wang F, Zhao R, Wang Z X, Zhang Y, and Lin Z H et al. 2021 Chin. Phys. Lett. 38 055201    Abstract Fusion born $\alpha$ particle confinement is one of the most important issues in burning plasmas, such as ITER and CFETR. However, it is extremely complex due to the nonequilibrium characteristics, and multiple temporal and spatial scales coupling with background plasma. A numerical code using particle orbit tracing method (PTC) has been developed to study energetic particle confinement in tokamak plasmas. Both full orbit and drift orbit solvers are implemented to analyze the Larmor radius effects on $\alpha$ particle confinement. The elastic collisions between alpha particles and thermal plasma are calculated by a Monte Carlo method. A triangle mesh in poloidal section is generated for electromagnetic fields expression. Benchmark between PTC and ORBIT has been accomplished for verification. For CFETR burning plasmas, PTC code is used for $\alpha$ particle source and slowing down process calculation in 2D equilibrium. In future work, 3D field like toroidal field ripples, Alfvén and magnetohydrodynamics instabilities perturbation inducing $\alpha$ particle transport will be analyzed. DOI:10.1088/0256-307X/38/5/055201 © 2021 Chinese Physics Society Article Text Energetic particles with energy from several tens to hundreds of keV to MeV are common in plasma physics, and they are extremely important in magnetic confinement fusion (MCF) plasmas.[1] In recent years, Chinese scientists have conducted many outstanding works in theory,[2] simulation[3–6] and experiment.[7,8] After several decades of research and development, the ITER project, which is a magnetic confinement device with power plant scale, has begun to assemble and plans to do deuterium-tritium operation in the next two decades. The primary target of the ITER project is to produce $500$ MW of fusion power from $50$ MW of power injected.[9] $\alpha$ particles from fusion reaction carry about $100$ MW power and this power is supposed to heat the plasmas to maintain the fusion reaction, which means the self-heating power is twice the injection power, and the heating efficiency from $\alpha$ particle is critical for high-$Q$ operation. Meanwhile, compared with the ionized fuel particles ($D^+$, $T^+$), elastic collisions for $\alpha$ particle are weak, which implies that the $\alpha$ particle keeps an unperturbed orbit trajectory for a long time and shows strong kinetic effects. In recent years, the Chinese fusion community has put forward the Chinese Fusion Engineering Testing Reactor (CFETR) project,[10,11] and the integrated modeling for core plasma is performed preferentially.[12] For CFETR burning plasmas with $200$ MW $\alpha$ particle power, it is important to make it clear whether the heating flux on the divertor or the first wall induced by $\alpha$ particle is destructive. Hao[13] and Zhao[14] have already calculated the heat flux induced by $\alpha$ particle TF ripple loss in two different operation scenarios. Particle simulation is a common tool in plasma physics. The particle-in-cell (PIC) method[15] is popular, which describes a mass number of charged particles in their self-consistent electric and magnetic fields. Gyrokinetic PIC is an improved version for magnetically confined plasma.[16] The PIC method is a powerful tool for instability problems in plasmas. However, for a long time scale physics, such as $\alpha$-particle slowing down process, the computational expense is unacceptable. Another option to describe the kinetic equation for particles is the test particle method, where the particle orbit is followed in the same way of PIC method, while the electromagnetic fields are not solved self-consistently. It is an effective numerical method to analyze the energetic particle phase space characteristic and confinement property in tokamak plasmas. Since the early 1980s, a series of numerical codes using test particle method were developed, such as NUBEAM,[17] ORBIT[18] and ASCOT.[19] These codes were mainly used to study ions phase space dynamics with moderate energies (tens of keV to hundreds of keV). For higher energetic particles (above $1$ MeV), Gyro-orbit effects could be important. In the last few years, test particle codes with full orbit option have been raised, such as SPIRAL,[20] orBIT (CUEBIT),[21] VENUS-LEVIS.[22] For CFETR-like devices, the analysis of $\alpha$ particle confinement is necessary for the design. Under this background, the particle orbit tracing code (PTC) has been developed.
cpl-38-5-055201-fig1.png
Fig. 1. Diagram of PTC code.
The diagram of the PTC code is shown in Fig. 1. The major component of the PTC code is to follow particle orbits in 3D electromagnetic field, and the PTC code has both the full orbit and drift orbit models for particle tracing. For the full orbit, the equations of motion can be written as $$\begin{aligned} \frac{d \boldsymbol{v}}{d t} &=\frac{q}{m}(\boldsymbol{v} \times \boldsymbol{B} +\boldsymbol{E}) + \dot {\boldsymbol{v}}_{\rm c}, \\ \frac{d \boldsymbol{x}} {d t} &=\boldsymbol{v}, \end{aligned}~~ \tag {1} $$ where $\dot{\boldsymbol{v}}_{\rm c}$ is the velocity change induced by elastic collisions. For the drift orbit model, the equations of gyro-center motion[23] can be written as $$\begin{aligned} \frac{d \boldsymbol{X}}{dt} &= \frac{\boldsymbol{B}^*}{B_\parallel ^*} v_\parallel + \frac{\mu}{m\varOmega B_\parallel^*}\boldsymbol{B} \times \nabla B + \frac{1}{B B_\parallel^*} \boldsymbol{E} \times \boldsymbol{B},\\ \frac{dv_\parallel}{dt} &= -\frac{\mu}{m} \frac{\boldsymbol{B}^*}{B_\parallel ^*} \cdot \nabla B + \frac{q}{m} \frac{\boldsymbol{B}^*}{B_\parallel ^*} \cdot \boldsymbol{E} +\dot{v}_{\parallel {\rm c}},\\ \frac{d \mu}{d t}&= \dot{\mu}_{\rm c}, \end{aligned}~~ \tag {2} $$ where $$\begin{aligned} \boldsymbol{B}^*&=\boldsymbol{B} +B \frac{v_\parallel}{\varOmega} \nabla \times \boldsymbol{b},\\ B_\parallel ^* &\equiv \boldsymbol{b} \cdot \boldsymbol{B}^* =B \Big(1+\frac{v_\parallel}{\varOmega} \boldsymbol{b} \cdot \nabla \times \boldsymbol{b}\Big), \end{aligned}~~ \tag {3} $$ and $\boldsymbol{b}=\boldsymbol{B}/B$, $\varOmega=qB/m$, $\mu=v_\perp ^2/2B$. The $\dot{v_\parallel}_{\rm c}$ and $\dot{\mu}_{\rm c}$ are the elastic collision inducing phase space changes. In principle, the full orbit model includes both cyclotron motion and guiding center drift motion, and the particle orbit can be described more accurately. However, for $\alpha$ particles in the CFETR-like tokamak, the gyro frequency of $\alpha$ particles is around $50$ MHz, and the gyroradius is around 0.5–5 cm, depending on the energy and position of the particle. This implies that the calculation of full orbit motion is time consuming, while the gyroradius effects can be either crucial or negligible for $\alpha$ particles in different physics processes. To balance the computation time and the physics effects, an option for automatic transition between the full orbit model and the drift orbit model is developed by monitoring the particle energy. The basic idea is that a critical characteristic size is set up for automatic transition between these two models: when the gyro radius is much larger than this characteristic size, then the full orbit model is used, otherwise the drift orbit model is selected. For the particles near the first wall, the gyro radius can also have some effects on particle interactions with the wall. Therefore, when particles are close to the wall, the full orbit model is used. The equations of motion are solved by the classical 4th order Runge–Kutta method, and for the full orbit model, the Boris method is an alternative option for more accuracy in energy conservation. The elastic collisions operator between test particles and background plasmas is solved by the Monte Carlo method, which is described below. Since the PTC is a test particle code, the electromagnetic field is defined externally. The $g$-eqdsk file[24] is used to define the equilibrium magnetic field. The poloidal magnetic flux $\psi$ in 2D uniform ($R$, $Z$) mesh is defined in the $g$-eqdsk file, and the poloidal current profile $f=RB_{\phi}$ is defined in 1D uniform $\psi$ mesh. The bilinear interpolation is used to represent the magnetic field in continuous space, which makes the divergence free condition satisfy well. For the full orbit model, the magnetic field is expressed as $$ B_R=\frac{1}{R} \frac{\partial \psi}{\partial Z},~~~ B_Z=-\frac{1}{R} \frac{\partial \psi} {\partial R},~~~ B_\phi= \frac{f}{R}.~~ \tag {4} $$ For drift orbit model, by definition, the gradient and curl of magnetic field can be expressed as $$\begin{aligned} \nabla B ={}&\nabla \sqrt{\Big(\frac{1}{R}\frac{\partial\psi}{\partial Z}\Big)^{2}+\Big(\frac{1}{R}\frac{\partial\psi}{\partial R}\Big)^{2}+\Big(\frac{f}{R}\Big)^{2}}, \\ \nabla \times \mathbf b ={}&\Big(\frac{1}{R} \frac{\partial b_Z}{\partial \phi}-\frac{\partial b_\phi}{\partial Z}\Big)\hat{\boldsymbol e}_R \\ &+\Big(\frac{\partial b_\phi}{\partial R}-\frac{1}{R}\frac{\partial b_R}{\partial\phi})\hat{\boldsymbol e}_Z \\ &+\Big(\frac{\partial b_R}{\partial Z}-\frac{\partial b_Z}{\partial R}\Big) \hat{\boldsymbol e}_\phi. \end{aligned}~~ \tag {5} $$ For tokamak plasmas, the equilibrium is symmetric in toroidal direction, which means $(\nabla B)_\phi =0$ and $(\nabla \times \boldsymbol{b})_\phi=0$. Figure 2 shows the orbit trajectory comparison of an $\alpha$ particle in the equilibrium magnetic field of the CFETR 2019 steady-state scenario. The contour, red line, and black line represent the flux surface, last closed flux surface (LCFS), and limiter, respectively. The blue-solid line and black-dashed line plot the full orbit and drift orbit calculated by the PTC code. Compared with the benchmark drift orbit (the red dot line) calculated by ORBIT code, the PTC drift orbit fits quite well with the trajectory of ORBIT code, especially in the region near the tip of the banana orbit, as shown on the right. It is proved that the input equilibrium module and particle push module in the PTC are correct.
cpl-38-5-055201-fig2.png
Fig. 2. Orbit comparison for an $\alpha$ particle in the CFETR 2019 steady-state scenario.
The perturbation field can be defined in an analytical or numerical model. For toroidal field (TF) ripple, an analytical model is used in the code. Following the Yushmanov's work,[25] the TF ripple can be approximated by $$\begin{aligned} \delta \boldsymbol{B}& = B_0 R_0 \Big[\frac{1}{N} {\sin}(N\phi)\nabla \delta + \delta {\cos}(N\phi)\nabla \phi\Big], \\ \nabla \delta &= \frac{\partial \delta}{ \partial R} \hat{\boldsymbol e}_R +\frac{\partial \delta}{\partial Z}\hat{\boldsymbol e}_Z, \end{aligned}~~ \tag {6} $$ where $\delta = \delta_0 {\cosh}(\alpha Z) J_N(\alpha R)$, and $\alpha$ is the profile for analytical magnetic perturbation, $J_N$ is the Bessel function of order $N$. In this format, the perturbation of poloidal field is included, and the divergence free condition for TF ripple is satisfied. For general 3D perturbation, analytically, the perturbation has the form in flux coordinate, $$ \delta \psi =\sum_{m,n}\delta\psi_{m,n}=\sum_{m,n}\hat{\psi}_{m,n}(r) \exp[-i(\omega t+n\phi+m\theta)].~~ \tag {7} $$ Since the code mainly adopts the cylindrical coordinates, the numerical Jacobian of flux coordinate in a triangle mesh for a poloidal section is implemented. As shown in Fig. 3(a), the computational domain is divided into two regions: core region and scrape-off-layer (SOL) region. In the core region, the vertexes are set on flux surfaces with an equidistant grid from magnetic axis to the last closed flux surface in low field side $R$ space. On each flux surface, the vertexes are placed with equidistant arc length in the poloidal section. The number of vertexes increases with minor radius to ensure that the vertexes are evenly spaced in general. For the SOL region, the inner boundary is the last closed surface, and the vertexes for the last closed surface are shared with the core region. The outer boundary is from $g$-eqdsk file. The meshes are generated by a constrained Delaunay triangulation method. With the triangle mesh, the geometry of tokamak poloidal section can be described precisely, and the flux coordinate can be interpolated conveniently.
cpl-38-5-055201-fig3.png
Fig. 3. (a) The PTC triangle mesh and (b) $\alpha$ particle source for the CFETR 2019 steady-state scenario.
The primary target of PTC code is to model the $\alpha$ particles from D-T fusion reactions. For particles with a velocity distribution, the fusion rate per volume $dR/dV$ is proportional to the fusion reactivity $\langle\sigma v\rangle$. The $\alpha$ particle source is defined by fusion rate per volume,[26,27] $$ \frac{dR}{dV}=\frac{n_i n_j}{1+\delta_{ij}}\langle\sigma v\rangle,~~ \tag {8} $$ where $n_{i}$ and $n_{j}$ are the fusion fuel's particle densities, and $\delta_{ij}$ is the Kronecker symbol. For particles with distributions $f(\boldsymbol{v}_i)$ and $f(\boldsymbol{v}_j)$, and the relative velocity $\boldsymbol{v}_r=\boldsymbol{v}_i-\boldsymbol{v}_j$, the fusion reactivity is defined by $$ \langle\sigma v\rangle =\int \int f({{\boldsymbol v}_i}) f({{\boldsymbol v}_j})\sigma (|\boldsymbol{v}_r|)|\boldsymbol{v}_r| d\boldsymbol{v}_i d\boldsymbol{v}_j.~~ \tag {9} $$ The nuclear reactions' cross section is fitted by $$ \sigma =\frac{S(E)}{E\exp(B_{\rm G}/\sqrt{E})},~~ \tag {10} $$ where $B_{\rm G}=\pi\alpha Z_{1}Z_{2}\sqrt{{2m_{r}c^{2}}}$ is the Gamov constant, $\alpha=e^{2}/\hbar c=1/137.03604$. The $S$ values are fitted by $$ S(E)=\frac{A_1+E[A_2+E[A_3+E(A_4+EA_5)]]}{1+E[B_1+E[B_2+E(B_3+EB_4)]]}.~~ \tag {11} $$ For D-T reaction, $B_{\rm G}=34.3827~(\sqrt{\rm{keV}})$, in energy range $E=[0.5, ~550]$ keV, $A_1=6.927 \times 10^4$, $A_2=7.454\times 10^8$, $A_3=2.05\times 10^6$, $A_4=5.2002\times 10^4$, $B_1=6.38\times 10^1$, $B_2=-9.95\times 10^{-1}$, and $B_3=6.981\times 10^{-5}$, $B_4=1.728\times 10^{-4}$. Using the center-of-mass reacting nuclei system, the energy of the system is $$ \epsilon =\frac{1}{2}m_r {v_r^2}.~~ \tag {12} $$ Here $v_r=|\boldsymbol{v}_r|$, and $m_r =\frac{m_i m_j}{m_i +m_j}$ is the reduced mass of the system. For Maxwellian velocity distributions, we have $$ f_i(v_i)=\Big(\frac{m_i}{2\pi k_{_{\scriptstyle \rm B}} T}\Big)^{3/2} \exp \Big(-\frac{m_i v_i ^2}{2k_{\rm B} T}\Big).~~ \tag {13} $$ Writing the volume element in velocity space as $d\boldsymbol{v}=4\pi v^2 dv$ and using the center-of-mass energy, we can finally reach $$\begin{alignat}{1} \langle \sigma v \rangle =\frac{4\pi}{\sqrt{ 2\pi m_r (k_{_{\scriptstyle \rm B}} T)^3} } \int_0 ^{\infty} \sigma(\epsilon) \epsilon \exp(-\epsilon/k_{_{\scriptstyle \rm B}} T) d \epsilon.~~ \tag {14} \end{alignat} $$ An example of D-T reaction fusion rate is shown in Fig. 3(b), which is based on the CFETR 2019 steady-state scenario. In tokamak core plasmas, the classical relaxation process for $\alpha$ particles is mainly dominated by elastic collisions. A Monte Carlo method is used for coulomb collisions[28] in full orbit phase space. With a local rectangular coordinate system $(\hat{e}_1, \hat{e}_{2},\hat{e}_3)$, where the particle initial velocity ($\boldsymbol{v}_0$) is parallel to $\hat{e}_1$, the collision operator can be approximated with finite time step $\Delta t$: $$ \boldsymbol{v}_1 =\boldsymbol{v}_0 +\Delta {\boldsymbol v}_\parallel \hat{e}_1 +\Delta {\boldsymbol v}_{\perp 1}\hat{e}_2 +\Delta {\boldsymbol v}_{\perp 2}\hat{e}_3,~~ \tag {15} $$ where $\Delta {\boldsymbol v}_\parallel = -\nu_{\rm sd}\Delta t_{\rm ec}{\boldsymbol v}_0 +\sqrt{\nu_\parallel \Delta t_{\rm ec} }{\boldsymbol v}_0 N_1$, $\Delta {\boldsymbol v}_{\perp1,2}=\sqrt{\frac{\nu_\perp \Delta t_{\rm ec}}{2}}{\boldsymbol v}_0 N_{2,3}$. $N_{i}$ are Gaussian random numbers with $\langle N_i\rangle =0$, $\langle N_i^2\rangle =1$. The coefficients $\nu_{\rm sd}$, $\nu_\parallel$ and $\nu_\perp$ correspond to slowing down, parallel diffusion and transverse diffusion, respectively. $$\begin{align} \nu_{\rm sd} &= \sum_{\beta} \frac{A_D^\beta}{2{\boldsymbol v}^3} \Big(1+\frac{m_\alpha}{m_\beta}\Big) \Big(\phi(\chi_\beta) -\chi_\beta \phi'(\chi_\beta)\Big), \\ \nu_\parallel &=\sum_{\beta} \frac{A_D^\beta}{{\boldsymbol v}^3}G(\chi_\beta), \\ \nu_\perp &= \sum_{\beta}\frac{A_D^\beta}{{\boldsymbol v}^3} [\phi(\chi_\beta) -G(\chi_\beta)].~~ \tag {16} \end{align} $$ Here $\chi_{\beta} ={\boldsymbol v}/{\boldsymbol v}_{\rm s,\beta}$ with ${\boldsymbol v}_{\rm s,\beta}=\sqrt{2T_\beta/m_\beta}$, $\phi(\chi)$ is the error function. The subscript $\alpha$ is corresponding to the test particle and $\beta$ is corresponding to the back-ground particle. $$ G(\chi) = \frac{\phi(\chi) -\chi \phi'(\chi)}{2\chi^2},~~ \tag {17} $$ $$ A_D^\beta = \frac{8\pi q_\alpha ^2 q_\beta ^2}{m_\alpha^2} n_\beta \ln \varLambda_\beta.~~ \tag {18} $$ The Coulomb logarithm is defined as $$ \ln \varLambda_\beta = \ln \frac{b_{\max}}{b_{{\min}}}.~~ \tag {19} $$ Here $b_{\max}~(b_{{\min}})$ is a typical maximum (minimum) impact parameter for the collision. For an ion moving through a plasma, $\bar{{\boldsymbol v}}_\beta ^2 =\bar{{\boldsymbol v}}_{\rm th,\beta} ^2 +{\boldsymbol v}^2$, and $$ b_{\max} = \Big[ \sum_\beta\frac{\omega_\beta^2}{\bar{{\boldsymbol v}}_\beta ^2}\Big]^{-1/2},~~ \tag {20} $$ where $\omega_p =\sqrt{4\pi n_\beta q_\beta^2/m_\beta}$ and ${\boldsymbol v}_{\rm th,\beta}=\sqrt{T_\beta/m_\beta}$. The lower limit $b_{{\min}}$ refers to the minimum impact parameter in the collision, and it is the larger one of $$ b_{{\min},\beta}^{\rm CL}=\frac{q_\alpha q_\beta}{m_{\alpha \beta} \bar{{\boldsymbol v}}_\beta^2}~~ \tag {21} $$ and $$ b_{{\min}}^{\rm QM}=\frac{\hbar}{2m_{\alpha \beta} \bar{{\boldsymbol v}}_\beta}.~~ \tag {22} $$
cpl-38-5-055201-fig4.png
Fig. 4. The $\alpha$ particle slowing down test for mean energy versus time (a), and the pitch angle scattering test for mean pitch angle versus time (b) with back-ground electrons $n_{\rm e}=1\times 10^{20}$ m$^{-3}$, $T_{\rm e}=20$ keV.
For $\alpha$ particles in typical tokamak plasma that $\nu$ is in the range $[1\times 10^{-3}, 1]$, the time step for collisions can be larger than the time step for particle orbit push, and the orbit average values of collision coefficients are used to increase the accuracy of the method. For the drift orbit model, the drift orbit phase space is transited to the full orbit phase space by inducing a random gyro phase $\theta$. After applying the collision operator, the phase space is transited back to the drift orbit phase space. A test run including 1000 $\alpha$ particles with initial energy $3.5$ MeV is set up to validate the collision operator. As shown in Fig. 4(a), after about 3 s, the average energy of $\alpha$ particles is slowing down to $30$ keV, which is about $1.5 T_{\rm e}$. As shown in Fig. 4(b), the pitch angle scattering effect is also tested. In this test, the particles are initially loaded with velocity parallel to $Z$ direction, then the mean pitch angle $\varLambda_Z = \frac{\boldsymbol{v}}{v}\cdot \boldsymbol{e}_Z =1$ at $t=0$, and with pitch scattering effect, the mean pitch angle converges to $0$. In summary, $\alpha$ particle physics is one of the most important key issues for magnetic confinement fusion. To better understand the burning plasmas in ITER or CFETR-like devices, a new numerical code (PTC) has been developed based on test particle method. The PTC code includes both the full orbit and drift orbit models to track charged particles in tokamak with 3D field. With unstructured triangle meshes, the tokamak geometry is described accurately including core plasma and SOL region. The $\alpha$ particle source is calculated using the fusion cross sections, and the elastic collision is modeled by a Monte Carlo method, which implies that this code can be used to calculate $\alpha$ particle heating for the background plasmas, and $\alpha$ particle neoclassical transport. For $\alpha$ particle transport in the SOL region, atomic reaction can be important and will be considered in the near future. References Energetic Particles in Magnetic Confinement Fusion PlasmasGyrokinetic theory of the nonlinear saturation of a toroidal Alfvén eigenmodeToroidal Alfvén eigenmode triggered by trapped anisotropic energetic particles in a toroidal resistive plasma with free boundaryLow-frequency fishbone driven by passing fast ions in tokamak plasmasEffect of Magnetohydrodynamic Perturbations on the Orbit Loss of Alpha Particles in Tokamak PlasmaVerification of Energetic-Particle-Induced Geodesic Acoustic Mode in Gyrokinetic Particle SimulationsReview of the experiments for energetic particle physics on HL-2AEnergetic Particle Physics on the HL-2A Tokamak: A ReviewEnergetic ions in ITER plasmasOverview of the present progress and activities on the CFETRProgress of the CFETR designSelf-consistent modeling of CFETR baseline scenarios for steady-state operationFusion α loss due to toroidal field ripple perturbation in CFETRAlpha particle ripple loss in CFETR steady-state scenarioParticle simulation of plasmasGyrokinetic approach in particle simulationThe tokamak Monte Carlo fast ion module NUBEAM in the National Transport Code Collaboration libraryHamiltonian guiding center drift orbit calculation for plasmas of arbitrary cross sectionASCOT: Solving the kinetic equation of minority particle species in tokamak plasmasA description of the full-particle-orbit-following SPIRAL code for simulating fast-ion experiments in tokamaksFull orbit computations of ripple-induced fusion α-particle losses from burning tokamak plasmasVENUS-LEVIS and its spline-Fourier interpolation of 3D toroidal magnetic field representation for guiding-centre and full-orbit simulations of charged energetic particlesProperties of energetic-particle continuum modes destabilized by energetic ions with beam-like velocity distributionsThe CHEASE code for toroidal MHD equilibriaEffects of magnetic ripple on 3D equilibrium and alpha particle confinement in the European DEMOImproved formulas for fusion cross-sections and thermal reactivitiesFOCUS: A full-orbit CUDA solver for particle simulations in magnetized plasmas
[1] Chen W and Wang Z X 2020 Chin. Phys. Lett. 37 125001
[2] Qiu Z Y, Chen L, and Zonca F 2019 Nucl. Fusion 59 066024
[3] Yang S X et al. 2018 Nucl. Fusion 58 046016
[4] Yu L M et al. 2019 Nucl. Fusion 59 086016
[5] Wu L N and Yu G Y 2002 Chin. Phys. Lett. 19 1312
[6] Chen Y et al. 2020 Chin. Phys. Lett. 37 095201
[7] Ding X T and Chen W 2018 Plasma Sci. Technol. 20 094008
[8] Shi P W et al. 2021 Chin. Phys. Lett. 38 035202
[9] Pinches S D et al. 2015 Phys. Plasmas 22 021807
[10] Wan Y X et al. 2017 Nucl. Fusion 57 102009
[11] Zhuang G et al. 2019 Nucl. Fusion 59 112010
[12] Chen J L et al. 2017 Plasma Phys. Control. Fusion 59 075005
[13] Hao B L et al. 2020 Sci. Sin.-Phys. Mech. Astron. 50 065201
[14] Zhao R et al. 2020 Plasma Phys. Control. Fusion 62 115001
[15] Dawson J M 1983 Rev. Mod. Phys. 55 403
[16] Lee W W 1983 Phys. Fluids 26 556
[17] Pankin A et al. 2004 Comput. Phys. Commun. 159 157
[18] White R B and Chance M S 1984 Phys. Fluids 27 2455
[19] Hirvijoki E et al. 2014 Comput. Phys. Commun. 185 1310
[20] Kramer G J et al. 2013 Plasma Phys. Control. Fusion 55 025013
[21] McClements K G 2005 Phys. Plasmas 12 072510
[22] Pfefferlé D et al. 2014 Comput. Phys. Commun. 185 3127
[23] Todo Y 2006 Phys. Plasmas 13 082503
[24] Lütjens H et al. 1996 Comput. Phys. Commun. 97 219
[25] Pfefferlé D et al. 2016 Nucl. Fusion 56 112002
[26] Bosch H S and Hale G M 1992 Nucl. Fusion 32 611
[27]Atzeni S and Meyer-ter-Vehn J 2004 The Physics of Inertial Fusion: Beam Plasma Interaction Hydrodynamics, Hot Dense Matter (Oxford: Clarendon Press)
[28] Clauser C F, Farengo R, and Ferrari H E 2019 Comput. Phys. Commun. 234 126