Self-consistent gyrokinetic modeling of turbulent and neoclassical tungsten transport in toroidally rotating plasmas Open Access

Tungsten (W), a high-Z element, is considered for use in the first-wall plasma-facing components in future tokamaks due to its high melting point, low erosion rate, and good thermal conductivity.¹ However, experiments have shown that the accumulation of W in the core region is frequently observed due to neoclassical inward convection, leading to a reduced plasma performance and energy loss through radiation.^2,3 In particular, even at low W concentration, radiative loss can cause complete extinction of the burning plasma.⁴ Consequently, it is crucial to mitigate the accumulation of W in tokamaks to improve the efficiency and performance of the reactor.

Over the past few decades, based on experimental observations of up–down radiative emission,^5,6 there has been a strong emphasis placed on the non-uniform distribution of impurity density across magnetic flux surfaces.^7–12 While the conventional neoclassical theory often assumes a uniform density distribution,^13,14 it is now widely accepted that an accurate prediction of neoclassical impurity transport requires a proper estimation of the poloidal asymmetry of impurity density, as it can enhance or reduce neoclassical impurity transport.^10,15–17

The toroidal rotation of plasma in tokamaks, whether driven intrinsically^18,19 or externally by neutral beam injection (NBI) heating, often results in a non-uniform distribution of W density. In particular, the centrifugal force generated by toroidal rotation drives impurities toward the outer wall, yielding a strong “in–out” poloidal asymmetry that leads to the enhancement of neoclassical transport.^10,20–22 Although the impact of toroidal rotation on turbulent transport is found to be much smaller than on neoclassical transport,²¹ the turbulent roto-diffusion generated by a radial gradient in the toroidal rotation²³ and the parallel velocity gradient (PVG) instabilities²⁴ can still play an important role in the case of a strongly rotating plasma.

The impact of toroidal rotation on both turbulent and neoclassical W transport was investigated in Ref. 21. The main approach used in the study was to compute turbulent and neoclassical transport separately. Turbulent transport was calculated using the gyrokinetic model in the toroidally co-moving frame, whereas neoclassical transport was evaluated using the drift-kinetic code NEO, but without considering frictional effects at the lowest order.^17,25

The aim of this paper, as an extension of the previous work on impurity transport from light to heavy impurities,²⁶ is to explore the effect of toroidal rotation on turbulent and neoclassical W transport using the full-F gyrokinetic code GYSELA.²⁷ By employing the multi-species collisional operator implemented in GYSELA,²⁸ we investigate the self-consistent treatment of both neoclassical and turbulent transport channels. The interplay between these channels can have a significant impact on impurity transport and its accumulation in the core of tokamaks.^29,30 To this end, we examine the impact of sonic (for the impurities) toroidal rotation on turbulent and neoclassical transport in the presence of poloidal asymmetry enhanced by centrifugal forces.

The remainder of the paper is organized as follows: First, in Sec. II, we introduce the GYSELA code and the source terms for studying impurity transport in a rotating plasma. Section III discusses the results of the nonlinear GYSELA simulations, with a focus on the poloidal asymmetry generated by the external momentum injection. In Sec. IV, the effect of sonic rotation on both turbulent and neoclassical W transport is explored, highlighting the inward neoclassical convection due to the poloidal asymmetry and the outward roto-diffusion due to the steepening of the velocity gradient. Finally, the paper concludes with a discussion of the main findings and conclusions in Sec. V.

The presence of the heat source $S heat$ in Eq. (1) is essential to allow flux-driven simulations to reach a statistical steady state. External torque is injected into the plasma via the momentum source $S mom$⁠. Each of these sources is designed to have a non-vanishing contribution to a single fluid moment equation, related to heat and parallel momentum in the present study.

It should be noted that, unlike gyrokinetic formulation in the co-moving frame that requires the inclusion of both the Coriolis drift and centrifugal effect in the derivation of guiding-center equations,^35–37 the equations of motion in GYSELA are solved in the laboratory frame. This implies that the effect of toroidal rotation is inherently present in the distribution function described in Eq. (8), without introducing any additional terms in the dynamics of the guiding center.

In this section, we describe a set of nonlinear GYSELA simulations to explore the effect of toroidal rotation on impurity transport. To this end, we consider W as the main impurity due to its large mass A_W = 184 and high charge state Z_W = 74 when fully ionized. However, we restrict our analysis to the partially ionized state Z_W = 40 as W is not fully stripped even at the core temperature of tokamaks. We conduct the simulations in a circular geometry with a numerical resolution of (⁠ $N r × N θ × N φ × N v ∥ × N μ ) = ( 256 × 512 × 32 × 128 × 64$⁠), using a time step of $Δ t = 16 ω c i − 1$⁠. The velocity space parameters are set to $0 < μ < 20$ and $− 7 < v ∥ s < 7$⁠, to account for both pitch angle scattering and energy diffusion of the collision operator, as well as supersonic impurities.²⁸ The nonlinear gyrokinetic Vlasov equation defined in Eq. (1) is solved for both main ions and W, while electrons are simplified to follow a Boltzmann response within a flux surface. The initial radial profiles of density and temperature are established at $R / L n = 2.2$ and $R / L T = 6$ for both species, triggering ITG turbulence. Other parameters, such as normalized gyro-radius $ρ * = ρ i / a = 1 / 190$⁠, inverse aspect ratio $ϵ = a / R 0 ≃ 0.2$⁠, and safety factor $q 95 ≃ 3$⁠, are kept constant throughout the present study. The radial profile of collisionality is determined by the normalized collisionality $ν s * = ν i q R 0 / ( v T s ϵ 3 / 2 )$⁠. For the analysis, we set the main ions in the banana regime with $ν D * = 0.1$ and W in the Pfirsch–Schlüter (PS) regime with $ν W * ≃ 23$ near the edge. The W concentration is restricted to the trace limit (⁠ $C W = 10 − 6$⁠), ensuring that the background turbulence remains unaffected by its presence.³⁸ 

In this study, by increasing the value of S_mom, higher values of Mach number are achieved. The radial profile of the momentum source remains unchanged with $ρ s , mom = 0.1$ and $L s , mom = 0.1$⁠. In addition, the amplitude and radial profile of the heat flux $S heat$ are kept constant throughout the analysis, with $S E = 0.017 , ρ s , heat = 0.012$⁠, and $L s , heat = 0.1$⁠. Heat and momentum are injected only for deuterium, indicating that the transfer of heat and kinetic energy occurs via inter-species collisional effects. Figure 1 shows the radial profile of the different source terms for the case with $S mom = 10 − 2$⁠. No inner boundary condition is applied to cover the magnetic axis $( ρ = 0 )$⁠, while Dirichlet boundary conditions are chosen at the separatrix $( ρ = 1 )$ to satisfy $ϕ sep = 0$⁠. A radial diffusion term, $D ( F ¯ s ) = 1 / ( r B ∥ s * ) ∂ r ( D r B ∥ s * ∂ r F ¯ s )$⁠, with a radial diffusion coefficient $D r = 0.03$⁠, is applied in the buffer region to mitigate numerical oscillations in the edge region.

To investigate its effect on density distribution and impurity transport, which is typically proportional to $∝ ( 1 + M 2 )$⁠,⁸ W ions are accelerated to the supersonic regime. Although the Mach number of main ions in fusion plasma is known to be low, typically $M D < 0.4$⁠, the Mach number of W can exceed the unity $( M W > 1 )$ due to the mass ratio $M W = A W / A D M D$⁠.

The numerical computation time required to reach a steady state is proportional to the gyro-Bohm scaling $∝ ρ * − 3$⁠,³⁹ which makes it infeasible to attain a steady state for multi-species simulations. To overcome this limitation, we first perform a simulation without W impurities for a sufficiently long time $( ∼ 10 5 ω c i − 1 )$ to develop the background turbulence and to saturate the toroidal velocity of main ions. Then, we carry out the second part of the simulation with W impurities for more than $10 5 ω c i − 1$ to get convergence on the toroidal rotation of W through collisional effects. Once the simulations are close to the quasi-steady state, all quantities are analyzed with a time window of $∼ 3000 ω c i − 1$ and are averaged over flux surfaces.

The radial profiles of the Mach number for both species are illustrated in Fig. 2 with different amplitudes of the momentum source. In the absence of the momentum injection S_mom = 0, the toroidal velocity, produced by turbulent intrinsic rotation,^18,19 exhibits low and flat profile for both species in the countercurrent direction, i.e., $M D ∼ − 0.02$ and $M W ∼ − 0.6$⁠. Once the momentum source is activated, the toroidal rotation increases and converges to a similar shape as the source profile illustrated in Fig. 1. W is propelled into the supersonic regime for amplitudes stronger than $S mom = 4 × 10 − 3$⁠. In the core region, the value of M_W reaches up to $M W ≃ 4$ with $S mom = 10 − 2$⁠. Despite the injection of momentum, the Mach number and its gradient for deuterium remain at low levels, with values of $M D < 0.4$ and $| d v φ , D / d r | ≃ 0.4$⁠, respectively. According to the criterion for the PVG instability to occur, which is $| d v φ , D / d r | > R / ( L n τ + 1 )$ with $τ = T e / T i$⁠,²⁴ the present nonlinear simulations are not subject to PVG instabilities. Similar Mach numbers have been observed in JET experiments with NBI injection,⁴⁰ where the Mach number for deuterium was found to be as high as $M D ≃ 0.7$⁠, yielding $M W ≃ A W / A D M D ≃ 4.5$⁠.

The radial profile of the flux-averaged density and temperature for both the main ions and W are displayed in Fig. 3. For deuterium, the density and temperature profiles are not significantly affected by the external momentum injection with the considered amplitudes of $2 × 10 − 3 ≤ S mom ≤ 10 − 2$⁠. The slight increase in T_D is likely due to the balance between the heat source and the outward heat flux, rather than the effect of toroidal rotation. Unlike deuterium, we observe significant changes in W density with increasing the amplitude of S_mom, while the temperature profiles remain unaffected and saturate at the same level as the main ions, indicating energy transfer through collisional effects. The overall density profile implies a central accumulation in cases of strong rotation.

Figure 4 describes the radial profile of the root mean square (rms) of the perturbed electric potential for different values of momentum injection. In the absence of momentum injection, the turbulence intensity remains below 4%. However, the injection of momentum tends to increase fluctuation levels in the region where the source is applied. Fluctuation levels in the region where $ρ > 0.5$ remain unchanged across different amplitudes, while an enhancement in fluctuation levels is observed in the inner region where $( ρ < 0.5 )$⁠. This indicates that the injected momentum does not stabilize the turbulence; instead, it results in a slight increase in turbulent levels in the inner region. The background ion temperature gradient (ITG) turbulence consistently remains throughout this study.

Previous works with GYSELA have demonstrated the stabilizing effect on turbulence via the formation of a transport barrier generated by a vorticity source term.^41,42 Furthermore, recent research with the JET-ILW (ITER-like wall) tokamak reports that steep temperature gradients at the plasma periphery can prevent W accumulation in the core region by enhancing outward neoclassical convection.⁴³ The impact of transport barrier on impurity transport using GYSELA is under investigation.

Non-uniform distributions of W obtained from the GYSELA simulations of rotating plasma are shown in Fig. 5. Three different cases are compared to elucidate the effect of toroidal rotation in driving an enhanced in-out asymmetry. In the absence of momentum injection S_mom = 0, an up–down asymmetry of W is observed mainly due to the background turbulence.⁴⁴ The degree of asymmetry is approximately 20%. As the amplitude of the source term increases, the geometry of asymmetry changes from up–down to in–out due to the centrifugal effects. In particular, the case with $S mom = 10 − 2$ yields an excess of approximately 40% of W localized in the outboard region.

The first dimensionless parameter, defined in Eq. (13), corresponds to a modified Mach number. The geometrical factors defined in Eqs. (11) and (12) scale quadratically with respect to the value of M. The second dimensionless parameter, denoted as g, characterizes the steepness of the density and temperature gradient profiles for the main ions. In the standard neoclassical theory, this parameter is assumed to be small, indicating a weak frictional force. However, its impact becomes dominant when impurities localize in the outboard region, particularly for high-Z impurities (⁠ $g ∝ Z 2$⁠) or in the pedestal region.

Figure 6 represents the non-uniform distribution of W in a toroidally rotating plasma, obtained from the GYSELA simulation on the left in Fig. 6, and its comparison with analytical expressions from Eqs. (9)–(12). The other parameters in the analytical expressions are evaluated from the same GYSELA simulation with $S mom = 10 − 2$⁠. Both theoretical approaches, namely, the Hamiltonian in a rotating frame and the drift-kinetic equation, are in qualitative agreement with the GYSELA simulations regarding the enhanced in-out asymmetry due to the centrifugal effects.

The expression for the density distribution derived in Eq. (9) shows a level of poloidal asymmetry comparable to that of GYSELA, whereas a factor of 1/2 is observed when using the expression derived from the drift-kinetic equation. The main difference between the two expressions is the inclusion of the poloidally varying electric potential $ϕ 1 ( ψ , θ )$⁠, which arises from toroidal rotation and turbulence. The discrepancies with neoclassical predictions can be attributed to the fact that the expressions in Eqs. (10)–(12) rely on neoclassical estimates of the poloidal ions velocity.^46–49 These estimates are not always consistent with experimental observations, suggesting that other mechanisms, such as intrinsic rotation induced by turbulence, may be important in this process.

Experimental observations^50–52 have revealed that a strong poloidal asymmetry, which is enhanced by the centrifugal force, can lead to a significant accumulation of W in the plasma core, primarily due to the inward neoclassical convection. Given these observations, Sec. IV will focus on exploring both neoclassical and turbulent transport of W in toroidally rotating plasma.

In general, neoclassical impurity transport depends on various plasma parameters for the main ions and impurities, i.e., $Γ W neo ∝ f ( N i , N W , T i , T W , ν z i )$⁠. However, unlike the standard neoclassical theory, a non-uniform distribution of impurity density has to be included to accurately predict the neoclassical flux. In Fig. 8, the effect of inhomogeneous W density on neoclassical W flux is represented by performing a scan of different values in δ and Δ from $N ̃ W = δ cos θ + Δ sin θ$⁠. The reconstructed 2D asymmetry of W distribution is then inserted in Eq. (19) to estimate the neoclassical flux, while other terms, measured from the GYSELA simulations, are kept constant. The poloidal asymmetry scan clearly indicates an increase in inward neoclassical flux as the asymmetry grows. In the case of up–down asymmetry [Fig. 8(b)], the trend is found to be symmetric, while the in–out asymmetry study [Fig. 8(a)] reveals a strong impact of impurity localization in the outboard region on neoclassical transport, even causing a change in the order of magnitude.

In the present study, the impact of toroidal rotation on turbulent flux is primarily mediated through the roto-diffusion term. The roto-diffusion term, proportional to the radial gradient of toroidal velocity, is known to be significant for W as the magnitude of the coefficient scales as $C u ∝ A W / Z W$⁠.

Figure 9 illustrates the cumulative turbulent and neoclassical fluxes obtained from GYSELA simulations at various levels of momentum injection. The fluxes are integrated both in time and volume within the range of $0.1 < ρ < 0.9$⁠. Consistently with the previous findings in Fig. 8, the neoclassical inward convection resulting from an enhanced in-out asymmetry of W density is observed to be proportional to the magnitude of toroidal rotation. A recent analysis of W accumulation in KSTAR experiments using the drift-kinetic code NEO⁶⁰ revealed that neoclassical inward convection is maximized at a specific toroidal rotation value, beyond which the outward convection increases.⁵¹ The non-monotonic relationship between the toroidal rotation and neoclassical inward convection can be attributed to the different dependencies of the convection terms on the main ion density and temperature profiles, as well as collisionality. An analytical investigation of neoclassical transport shows that the temperature screening of impurities is enhanced at sufficiently high Mach number and low collisionality, while this effect decreases as collisionality increases.⁶¹ In the present study, however, the impact of toroidal rotation on the main profiles is found to be negligible, as depicted in Fig. 3, and only a high collisional case with W is considered. This explains the linear relationship observed between neoclassical inward convection and the magnitude of toroidal rotation in Fig. 9.

Although the toroidal rotation has a direct impact on neoclassical transport owing to the presence of strong density inhomogeneity, its influence on turbulent transport in Fig. 9 (right) is observed to be relatively small. The enhancement of outward turbulent flux resulting from the increase in the toroidal rotation, as depicted in Fig. 9, can be attributed to the increased turbulence intensity observed in Fig. 4 and the turbulent pinch terms defined in Eq. (24), i.e., roto-diffusion and thermo-diffusion terms. Although both terms are known to drive outward convection in the case of ITG turbulence, the roto-diffusion predominates over the thermo-diffusion term, as the former is proportional to $A W / Z W = 4.6$⁠, while the latter is proportional to $1 / Z W = 0.025$⁠. In fact, for heavy impurities with strong toroidal rotation, the magnitude of roto-diffusion can be as significant as the diffusive part.²³ 

The impact of the roto-diffusion term on the W turbulent flux becomes apparent when we compare the time evolution of the radial gradient of toroidal velocity and the turbulent flux for W, as depicted in Fig. 10. The radial gradient of toroidal rotation gradually increases over time in proportion to the injected momentum. For W, the values of $| d v φ , W / d t |$ are observed to exceed unity for amplitudes above $S mom = 4 × 10 − 3$ due to its heavy mass, whereas the gradient remains low at $| d v φ , D / d r | ≃ 0.4$ for the main ions. As shown in Fig. 9, the outward turbulent flux emerges at $t = 10 5 ω c i − 1$ for amplitudes greater than $S mom = 6 × 10 − 3$⁠. The results demonstrate that the influence of the roto-diffusion term becomes more significant compared to other terms in Eq. (24) as the absolute value of $| d v φ , W / d r |$ gradually increases. Additionally, it is noteworthy that centrifugal effects also have a non-negligible impact on other convection terms, particularly in the case of heavy impurities.⁵⁵ Consequently, the sign and magnitude of the turbulent flux depend on the interplay between different convection terms, which typically results in an outward flux in the case of ITG turbulence.

Figure 11 presents the W accumulation obtained from the total flux, $Γ W tot = Γ W neo + Γ W turb$⁠, within the radial range of $ρ = 0.1$ and $ρ = 0.3$ for different values of the injected momentum. No significant changes in the core W density are observed up to the time $t = 0.5 × 10 5 ω c i − 1$⁠, which corresponds to the momentum transfer time from the main ions to W. Following this transfer time, a more significant accumulation of W is observed, approximately 20% of increase in the $S mom = 10 − 2$ scenario. The central accumulation of W in GYSELA simulations is commonly observed even in the absence of external toroidal rotation.²⁶ The numerical results obtained in the present study are consistent with experimental evidence from different tokamaks^51,53,54 where central NBI heating yields a central deposition of W. To prevent the accumulation of W within the core region, the proper use of ICRH and ECRH systems can be helpful in controlling the W within the core region.^{9,33,52,62–66} This is achieved by flattening the density profile of main ions and enhancing the temperature gradient of the bulk plasma, which facilitates the outward turbulent diffusion of W and an increase in the temperature screening effect. In addition, the temperature anisotropy that arises from ICRH is known to reduce poloidal asymmetry, thereby decreasing inward neoclassical convection.^9,67

In the present paper, the effects of toroidal rotation on both turbulent and neoclassical W transport are addressed using the nonlinear, global, and flux-driven 5D gyrokinetic code GYSELA. The injection of toroidal momentum via the external source term drives W to the supersonic regime, while the toroidal velocity of the main ions remains low. By leveraging the multi-species collisional operator and external source terms implemented in GYSELA, both turbulent and neoclassical W transport in a toroidally rotating plasma are treated self-consistently. Unlike the conventional approach, which treats neoclassical and turbulent transport separately, the self-consistent treatment used in this study captures the interplay between these two channels of W transport.

Through a series of nonlinear simulations involving different levels of momentum injection, it has been observed that centrifugal forces cause W to accumulate in the outboard region, creating an in-out poloidal asymmetry. This asymmetry has been found to significantly enhance the neoclassical inward convection, leading to a central accumulation of W in the case of strong plasma rotation. The core W accumulation accelerated by the centrifugal force and the resulting poloidal asymmetry, cannot be offset by the turbulent flux. However, as the momentum injection continues, the roto-diffusion term becomes as significant as the diffusive part once the radial gradient of the rotation velocity exceeds a certain threshold, driving outward turbulent flux in the case of ITG turbulence.

The numerical results obtained in the present study are in qualitative agreement with theoretical predictions. However, for a more comprehensive modeling of W transport in GYSELA, it could be interesting to include the time evolution of the density and temperature profiles of main ions, taking into account their response to other heating systems. These additional heating systems are known to mitigate core W accumulation by reducing poloidal W asymmetry and enhancing the temperature screening effect.

This work has been carried out within the framework of the EUROfusion Consortium, partially funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No. 101052200–EUROfusion). The Swiss contribution to this work has been funded by the Swiss State Secretariat for Education, Research and Innovation (SERI). Views and opinions expressed are, however, those of the author(s) only and do not necessarily reflect those of the European Union, the European Commission or SERI. Neither the European Union nor the European Commission nor SERI can be held responsible for them. The simulations presented herein were carried out on the CINECA Marconi supercomputer under FUA35 GSNTITE project.

The authors have no conflicts to disclose.

Kyungtak LIM: Conceptualization (equal); Formal analysis (equal); Software (equal); Validation (lead); Visualization (lead); Writing – original draft (lead). Xavier Garbet: Conceptualization (equal); Formal analysis (equal); Supervision (equal). Y. Sarazin: Conceptualization (equal); Formal analysis (equal); Supervision (equal). Etienne Gravier: Conceptualization (equal); Formal analysis (equal); Supervision (equal). Maxime Lesur: Conceptualization (equal); Formal analysis (equal); Supervision (equal). Guillaume LO-CASCIO: Conceptualization (equal). Timothé Rouyer: Conceptualization (supporting).

The data that support the findings of this study are available from the corresponding author upon reasonable request.