Modelling of charge-exchange induced NBI losses in the COMPASS upgrade tokamak

Detailed integrated transport modelling with the METIS code [18] yields density and temperature profiles during the flat-top with 4 MW of NBI heating and 2 MW of ECRH heating. We consider a set of equilibria covering magnetic fields and current ranging from 2.5 T to 5 T and 0.8 MA to 2 MA, respectively. The pressure and current function (RB_φ ) radial profiles as obtained from the transport simulations in METIS are used as inputs for the free boundary code FIESTA [19]. The minor radius is approximately a ≃ 0.27 m.

Table 1 summarizes the relevant plasma parameters for the equilibrium considered in this study. The energy confinement time used is the one from [20] (ITERH98P(y,2)). Heat diffusion profile shapes are obtained from a Bohm gyro-Bohm scaling. The ratio of heat transport between electrons and ions is based on the ITG-TEM (ion temperature gradient-trapped electron mode) diagram. The value we retained for the temperature decay length is (based on [21]):

$$\begin{equation}{\lambda }_{T}=2.88{B}_{\mathrm{t}}^{-0.5}{q}_{\text{cyl}}^{0.974}\text{.}\end{equation} \tag{ 4 }$$

Here q_cyl = aB_t/(R₀ B_pol) where B_pol is the average value of the poloidal field at the separatrix. The density decay length in the SOL is taken as λ_n = 1.35λ_T .

In this publication we refer to two categories of scenarios: the baseline scan: A, B and C (magnetic field scan at low safety factor) and the density scan: D–F (scan in increasing line-average density at 5 T and lower current).

Typical profiles obtained during the main heating phase are represented for both categories of scenarios in figure 1. These were calculated for the most tangential injection of the NBI.

Zoom In Zoom Out Reset image size

Neutral gas arises in the boundary region of tokamak plasma either because of recycling from the walls or by external gas puffing, which is often used to control the plasma density. The neutral gas in this region can be either molecular or atomic in form. Our modelling is using simple results from the KN1D code [8] that covers a large range of physical processes but only gives solution for a profile in the mid-plane. We have extrapolated to the whole poloidal plane these profiles of neutral densities by arbitrarily making them a flux surface quantity. Furthermore, we have neglected local toroidal asymmetries induced by the beam deposition (known as halo neutrals [22]).

Due to geometric constrains and limited space in the openings around the torus, we are limited in tangency radius to about R_tan ≃ 0.65 m which does not allow for on-axis injection (R_axis ≃ 0.9 m). The NBIs will be located by pairs, on top of each other in each injector box (2 units of 1 MW of injection power per box). As such, they will have a vertical tilt. For simplicity, we have only considered the downward ≃5° inclination in our study. Results do not differ significantly for the paired ≃−5° upward inclination.

The neutral beam injector launches fast neutrals into the plasma. Most of them become fast ions as they collide with the plasma. The newly-developed ray-tracing code NUR calculates the initial ionization positions and velocity vector of the ions using data from [23]: its principle is given in appendix A. NUR also gives an estimate of the shine-through of the beam (see figure 2). As COMPASS upgrade is going to be operated with high plasma densities, the shine-through losses are low, as can be seen in figure 2: the fraction of power lost does not exceed 3.3% in the low density scenarios A and D. Therefore we are thereafter considering that all injected power by NBI is ionized without shine-through losses. In the remainder of the paper, we have considered that 1 MW was injected by the single NBI unit that we have modelled.

Zoom In Zoom Out Reset image size

Besides the tangency radius, the main beam geometry parameters used for the calculation of the deposition are the distance from the source to the tangency point (d_s ≃ 5.67 m), and the elevation of the source over the mid-plane (Z_s ≃ 0.425 m). We are considering a circular source grid of radius 10 cm as the source. According to the beam technical specifications, the beamlet source divergence is taken as 6.5 mrad. The beamlets are then converging at the focal point, located 5 m away from the source. This implies that the overall beam divergence will likely be higher than that of each individual beamlet. We also consider an aperture (at the narrowest part of the port) that can clamp the edges of the Gaussian distribution of the beam.

The NBI has injection energy fractions (expressed in number of deposited ions) of about 46% at 80 keV, 39% at 40 keV and 15% at 27 keV. Figure 2 illustrates the NBI deposition for the scenario with typical operational parameters (4400, scenario B). In that case, the plasma has a density of about $\left\langle {n}_{\mathrm{e}}\right\rangle \simeq$ 1.5 × 10²⁰ m⁻³: the shine-through is negligible and the deposition of the NBI generates a majority of passing particles with negligible first orbit losses (in an axisymmetric field).

The NUR code passes the distribution of newly born NBI ions as input to the orbit following code EBdyna in a format referred to as ion markers (corresponding for each ion to a given energy, velocity vector and position).

Coils are described as 3D wire filaments, i.e. without thickness. The current (at least during the first second of the flat-top) is considered to be mostly localized in the inner part of the coil, as anticipated by more sophisticated ANSYS modelling. A sketch of the 3D TF wire is shown in figure 3, next to the more realistic drawing as used in ANSYS modelling. The Biot–Savart modelling discretizes a wire as a combination of very small (∼1 mm) current elements. The Biot–Savart law yields:

$$\begin{equation}\mathbf{B}=\frac{{\mu }_{0}}{4\pi }{\int }_{\text{Coil}}\frac{\mathbf{I}{\times}\mathbf{r}}{{r}^{3}}\enspace \mathrm{d}l\end{equation} \tag{ 5 }$$

$$\begin{equation}\mathbf{A}=\frac{{\mu }_{0}}{4\pi }{\int }_{\text{Coil}}\frac{\mathbf{I}}{r}\enspace \mathrm{d}l\end{equation} \tag{ 6 }$$

with B = ∇ × A, I the current (oriented along the wire) in the coil and dl the current element length. The vector r is representing the position where the field is calculated. The average value of the vacuum perturbed field is subtracted and then the axisymmetric field from the equilibrium is added in order to get the total 3D magnetic field. The solution obtained from the Biot–Savart solver is compared with a more elaborate ANSYS calculation and yields similar results. The ripple δ_r = (B_max − B_min)/(B_max + B_min) is evaluated to be about 0.55% at the separatrix in the mid-plane, as illustrated in figure 3.

Zoom In Zoom Out Reset image size

Each NBI initial distribution is represented by 3 × 10⁵ markers as obtained by the NUR calculation. These markers are injected at regular time interval in the full-orbit simulation. The integrated evolution of the position of NBI D ions in the 3D TFR field is done using a grid in (R, Z, φ) repeated 16 times (the number of TF-coils). The algorithm used within the EBdyna code is given in appendix B.

An illustration of the evolution of the losses induced by the ripple is given on the right in figure 3. The initial distribution of the fast ions is the axi-symmetric NBI steady-state. Various parametrization of the simulation are considered. In particular, the TFR is switched off and on and the impact on the overall losses is clearly seen. When moving to collisional TFR simulation, one can notice, especially in the case of more tangential injection (R_tan = 45 cm), that the overall losses are increased due to refueling of the loss region in velocity space.

Radial transport of trapped fast ions shows a strong correlation to the ratio of the bounce (ω_b) and precession (ω_φ ) frequencies [24]. In particular, significant radial excursions are found at ω_b/ω_φ ≃ [16/1; 16/2; 16/3] ratios. On figure 4 the collisionless motion of the trapped particles injected with R_tan = 15 cm has been studied during 1 ms in terms of losses to the wall and also their rate of change of radial position δp_φ /δt (p_φ being the toroidal canonical angular momentum defined in expression (B.7)). The initial distribution of markers was the steady state fast ion distribution.

Zoom In Zoom Out Reset image size

The interaction of the NBI particles born from co-perpendicular injection geometry (e.g. R_tan = 15 cm) is complex: overlap of resonances seems to be causing enhanced stochastic losses. The increased losses at lower tangency radius are also explained by the larger amount of markers existing close the main resonances in ω_b/ω_φ ≃ [16/1; 16/2; 16/3], especially in the near-edge region (ρ_t > 0.75). The results of the TFR simulations are not changed significantly when we include the Coulomb collisions. The most noticeable effect is the enhancement of the losses due to the refilling of the loss region in velocity space (as was shown in figure 3, right).

A summary of the effect of the TFR for each scenario and injection geometry is given in the next section in table 2.

Expressions for the Coulomb logarithm derived from [25] for the ions and [10] for the electrons are used:

$$\begin{equation}\begin{aligned}\hfill {{\Lambda}}_{\mathrm{D}\mathrm{i}}& =\frac{1}{2}\enspace \mathrm{ln}\left(\frac{{{\epsilon}}_{0}}{\mathrm{e}}\frac{{T}_{\mathrm{e}}}{{n}_{\mathrm{e}}}\right)-\mathrm{ln}\left(\frac{2{\mathrm{e}}^{2}}{4\pi {{\epsilon}}_{0}{m}_{\mathrm{D}}{v}^{2}}\right)\text{,}\hfill \\ \hfill {{\Lambda}}_{\mathrm{D}\mathrm{e}}& =15.2-\frac{1}{2}\enspace \mathrm{ln}\left(1{0}^{-20}{n}_{\mathrm{e}}\right)+\mathrm{ln}\left(1{0}^{-3}\enspace \mathrm{max}{\left({v}_{\text{the}},{v}_{\mathrm{f}}\right)}^{2}\right)\hfill \end{aligned}\end{equation} \tag{ 7 }$$

where T_e is the electron temperature in eV, v_the the corresponding thermal electron velocity and v_f the fast ion velocity. The saturation on velocity used on the Λ_De is a crude approximation that reflects the fact that the relative velocity of collisions remains constrained by the highest velocity of the reactants considered.

In order to completely describe the range of electron velocities available in the SOL, we use the formulas derived in [6] to express the collision rates on electrons (ν_{b e}) and on ions (ν_{b i}):

$$\begin{equation}\begin{gathered}{\nu }_{\text{be}}\simeq \left(\frac{{\mathrm{e}}^{4}{n}_{\mathrm{e}}}{4\pi {{\epsilon}}_{0}^{2}{m}_{\mathrm{D}}{m}_{\mathrm{e}}}\enspace \mathrm{ln}\enspace {{\Lambda}}_{\mathrm{b}\mathrm{e}}\right)\frac{\psi \left({x}_{\mathrm{e}}\right)}{{v}_{\mathrm{f}}^{3}}\\ {\nu }_{\text{bi}}\simeq \left(\frac{{\mathrm{e}}^{4}{n}_{\mathrm{i}}}{4\pi {{\epsilon}}_{0}^{2}{m}_{\mathrm{D}}^{2}}\enspace \mathrm{ln}\enspace {{\Lambda}}_{\mathrm{b}\mathrm{i}}\right)\frac{\psi \left({x}_{\mathrm{i}}\right)-{\psi }^{\prime }\left({x}_{\mathrm{i}}\right)}{{v}_{\mathrm{f}}^{3}}\end{gathered}\end{equation} \tag{ 8 }$$

where n_i indicates the background deuterium ion density, ${x}_{\mathrm{e}}={\left({v}_{\mathrm{f}}/{v}_{\text{the}}\right)}^{2}$, ${x}_{\mathrm{i}}={\left({v}_{\mathrm{f}}/{v}_{\text{thi}}\right)}^{2}$ (v_the and v_thi being the electrons and ions thermal velocities). The Maxwell integral ψ(x) and its derivative ψ'(x) are expressed as:

$$\begin{equation}\psi \left(x\right)=\frac{2}{\sqrt{\pi }}{\int }_{0}^{x}\sqrt{t}\enspace {\mathrm{e}}^{-t}\enspace {\mathrm{d}\mathrm{t}\quad \text{and}\quad \mathit{\psi }}^{\prime }\left(x\right)=\frac{2}{\sqrt{\pi }}\sqrt{x}\enspace {\mathrm{e}}^{-x}\end{equation} \tag{ 9 }$$

and for which relevant Taylor expansions can be used as given in [6]. Combined with the equation (8), this allows to express continuously the transition from v_f < v_the to v_f > v_the. The latter can occur in the SOL region.

$$\begin{equation}\begin{aligned}\hfill {\nu }_{\text{be}}& \simeq \left(\frac{{\mathrm{e}}^{4}{n}_{\mathrm{e}}}{\pi {{\epsilon}}_{0}^{2}{m}_{\mathrm{D}}{m}_{\mathrm{e}}}\enspace \mathrm{ln}\enspace {{\Lambda}}_{\mathrm{b}\mathrm{e}}\right)\frac{1}{{v}_{\text{the}}^{3}}\frac{1}{3\sqrt{\pi }}\hfill \\ \hfill & \quad {\times}\left(1-\frac{3}{5}{\left(\frac{{v}_{\mathrm{f}}}{{v}_{\text{the}}}\right)}^{2}+\frac{3}{14}{\left(\frac{{v}_{\mathrm{f}}}{{v}_{\text{the}}}\right)}^{4}\right)\text{,}\enspace {v}_{\mathrm{f}}{< }{v}_{\text{the}}\hfill \\ \hfill {\nu }_{\text{be}}& \simeq \left(\frac{{\mathrm{e}}^{4}{n}_{\mathrm{e}}}{4\pi {{\epsilon}}_{0}^{2}{m}_{\mathrm{D}}{m}_{\mathrm{e}}}\enspace \mathrm{ln}\enspace {{\Lambda}}_{\mathrm{b}\mathrm{e}}\right)\frac{1}{{v}_{\mathrm{f}}^{3}}\left(1-\left(\frac{2}{\sqrt{\pi }}\sqrt{{x}_{\mathrm{e}}}\enspace {\mathrm{e}}^{-{x}_{\mathrm{e}}}\right)\right.\hfill \\ \hfill & \quad {\times}\left.\left(1+\frac{2}{2{x}_{\mathrm{e}}}-\frac{1}{4{x}_{\mathrm{e}}^{2}}\right)\right)\text{,}\enspace {v}_{\mathrm{f}}{ >}{v}_{\text{the}}\text{.}\hfill \end{aligned}\end{equation} \tag{ 10 }$$

Of course, the usual expression for the most common case regarding ν_{b e} is that of electrons verifying v_f < v_the.

The average momentum loss is simply expressed as:

$$\begin{equation}\frac{\mathrm{d}{v}_{\mathrm{f}}}{\mathrm{d}t}=\left({\nu }_{\text{be}}+{\nu }_{\text{bi}}\right){v}_{\mathrm{f}}\simeq \left(\delta {v}_{\mathrm{e}}+\delta {v}_{\mathrm{i}}\right)/\delta t\enspace \text{.}\end{equation} \tag{ 11 }$$

where δv_e = ν_be v_f and δv_i = ν_bi v_f refer to the norm of velocity lost during the small time δt to the electrons and the ions respectively. For completeness—as derived for the NUBEAM code in [26]—the standard deviation related to the variation of the norm of velocity during the small time δt is included as:

$$\begin{equation}\frac{\left\langle {\Delta}{v}_{\mathrm{f}}^{2}\right\rangle }{\delta t}=2{\nu }_{\mathrm{b}\mathrm{e}}\left(\frac{e}{{m}_{\mathrm{f}}}\left({T}_{\mathrm{e}}+{\left(\frac{{v}_{\mathrm{f}}}{{v}_{\mathrm{c}}}\right)}^{3}\right)\right)\end{equation} \tag{ 12 }$$

with

$$\begin{equation}{v}_{\mathrm{c}}=\frac{\sqrt{\pi }}{4}{\left(\frac{{m}_{\mathrm{e}}}{{m}_{\mathrm{D}}}\right)}^{1/3}{v}_{{\text{th}}_{\mathrm{e}}}\end{equation} \tag{ 13 }$$

and we arbitrarily separate this spread on electrons and ions (for power deposition calculations) according to:

$$\begin{equation*}\delta {v}_{\sigma \mathrm{e}}=\sqrt{\left(\delta t{\nu }_{\mathrm{b}\mathrm{e}}\right)\frac{2\enspace \mathrm{e}}{{m}_{\mathrm{f}}}{T}_{\mathrm{e}}}\end{equation*}$$

and

$$\begin{equation}\delta {v}_{\sigma \mathrm{i}}=\sqrt{\left(\delta t{\nu }_{\mathrm{b}\mathrm{e}}\right)\frac{2\enspace \mathrm{e}}{{m}_{\mathrm{f}}}{T}_{\mathrm{i}}{\left(\frac{{v}_{\mathrm{f}}}{{v}_{\mathrm{c}}}\right)}^{3}}\text{.}\end{equation} \tag{ 14 }$$

The change in norm of velocity vector after the small time δt is thus:

$$\begin{equation}\delta {v}_{\mathrm{f}}=\left(\delta t\left({\nu }_{\mathrm{b}\mathrm{e}}+{\nu }_{\mathrm{b}\mathrm{i}}\right)\right){v}_{\mathrm{f}}+\mathcal{N}\left(1,1\right)\delta {v}_{\sigma \mathrm{e}}+{\mathcal{N}}^{\prime }\left(1,1\right)\delta {v}_{\sigma \mathrm{i}}\enspace \text{,}\end{equation} \tag{ 15 }$$

where $\mathcal{N}\left(1,1\right)$ and ${\mathcal{N}}^{\prime }\left(1,1\right)$ refer to random numbers following the normal distribution of mean 0, variance 1 and standard deviation 1. The energy given to the electrons and the ions are then respectively:

$$\begin{equation*}\delta {\mathcal{E}}_{\mathrm{e}}=\frac{{m}_{\mathrm{D}}}{2}\left({v}_{\mathrm{f}}\left(\delta {v}_{\mathrm{e}}\right)-{\left(\delta {v}_{\mathrm{e}}\right)}^{2}\right)\end{equation*}$$

and

$$\begin{equation}\delta {\mathcal{E}}_{\mathrm{i}}=\frac{{m}_{\mathrm{D}}}{2}\left({v}_{\mathrm{f}}\left(\delta {v}_{\mathrm{i}}\right)-{\left(\delta {v}_{\mathrm{i}}\right)}^{2}\right)\end{equation} \tag{ 16 }$$

where the fluctuations related to the standard deviation have been omitted for conciseness.

In addition to energy loss, the fast ions experience perpendicular diffusion in velocity space during their thermalization (pitch angle diffusion). We can neglect the change of direction induced by the electrons (lower by a factor m_e/m_D). The deviation of the velocity vector due to the collisions with ions is based on the characteristic time for a π/2 deviation given in [6] and is obtained by applying at each time step the change in solid angle:

$$\begin{equation*}\delta {\Omega}=\mathcal{N}\left(1,1\right)\frac{\pi }{2}\sqrt{\left(\delta t\right){\nu }_{\perp \mathrm{i}}}\end{equation*}$$

with

$$\begin{equation}{\nu }_{\perp \mathrm{i}}={\nu }_{\mathrm{i}}\left(\psi \left({x}_{\mathrm{i}}\right)+{\psi }^{\prime }\left({x}_{\mathrm{i}}\right)-\frac{\psi }{2{x}_{\mathrm{i}}}\right)\end{equation} \tag{ 17 }$$

with the remaining degree of freedom being taken randomly in a uniform distribution between 0 and 2π.

The slowing down process due to Coulomb collisions is applied to the injected markers until they reach the energy of 1.5T_i when they are then considered 'thermalized' and removed from the simulation. Their remaining energy is however still considered 'given' to the plasma ions in the total power balance. In order to be computationally efficient, we only apply the collision operators every 1 × 10⁻⁸ s which is found to be enough to reflect the local values of density and temperature along the particle orbit. A typical simulation of the 3 × 10⁵ markers lasts a few hours on a hundred processors. Once a steady-state of the beam distribution is obtained after the time τ_s, one can average a few subsequent outputs of the code in order to obtain a detailed map of the velocity space and density of particles. This average is typically taken over an extra ∼1.5 ms of simulation.

An illustration of the EBdyna output in the axisymmetric case is given in terms of deposited power and velocity space distribution and compared with the output from the NUBEAM code in figure 5 for scenario B. The velocity space distributions have very similar patterns in this specific example as well as in all other scenarios. The NBI density and power deposition radial profiles are in excellent agreement except in the pedestal region where it is found that the full orbit code EBdyna predicts a larger heating contribution to electrons than NUBEAM.

Zoom In Zoom Out Reset image size

The charge exchange processes we have considered in our study are:

$$\begin{equation}\begin{gathered}{\mathrm{D}}^{+}+\mathrm{D}\to \mathrm{D}+{\mathrm{D}}^{+}\\ {\mathrm{D}}^{+}+{\mathrm{D}}_{2}\to \mathrm{D}+{\mathrm{D}}_{2}^{+}\end{gathered}\end{equation} \tag{ 18 }$$

where the neutral reactants are in the ground (not excited) state. We denote σ_D and ${\sigma }_{{\mathrm{D}}_{2}}$ respectively the cross sections corresponding to these reactions : the equations found in [7, 27] are used to model the cross section as a function of reactants relative energy is given in appendix C. The probability for an ion not to undergo the charge exchange reaction after the time t is then by definition:

$$\begin{equation}\mathcal{P}\left(t\right)=\prod\limits _{j}\enspace \mathrm{exp}\left(-\frac{\delta t}{{\tau }_{{\text{cx}}_{j}}}\right)\quad \text{with}\quad {\tau }_{{\text{cx}}_{j}}=\frac{1}{\left\langle {\sigma }_{{\mathrm{D}}_{j}}v\right\rangle {n}_{{0}_{j}}}\end{equation} \tag{ 19 }$$

where ${n}_{{0}_{j}}$ is the relevant background neutral density that are evaluated at successive instant t = jδt and v is the velocity of the ion at time t. $\left\langle {\sigma }_{{\mathrm{D}}_{j}}v\right\rangle$ refers to the reaction rate obtained by averaging the cross section values against a Maxwell–Boltzman distribution in velocities of the background neutrals or ions (see appendix C for details). For molecular deuterium, we only need the cold background approximation and the expression is simplified.

Conversely, the probability for a neutral to become ionized because of charge exchange is obtained by replacing ${n}_{{0}_{j}}$ with the local ion density. Furthermore, a fast neutral can be re-ionized because of the ion impact ionization process:

$$\begin{equation}\mathrm{D}+{\mathrm{D}}^{+}\to {\mathrm{D}}^{+}+{\mathrm{D}}^{+}{\mathrm{e}}^{-}\end{equation} \tag{ 20 }$$

although this process is less significant than the CX one, as is discussed in appendix C. Figure 6 illustrates the distribution of—eventually multiple—CX neutralizations and subsequent re-ionizations. The CX neutralizations are strongly localized on the LFS in the pedestal and near SOL region whilst the re-ionizations are much more scattered all over the poloidal plane. This example, based on scenario A, has a power loss fraction of 15.3% due to CX neutrals lost. This can be seen on the right figure where the lost power density at the LFS mid-plane is the difference between the blue and red curves. The total fraction of the NBI lost to the wall (neutrals and ions) is 24.4% in this configuration.

Zoom In Zoom Out Reset image size

Using the steady state distribution of the NBI, we can track the resulting markers on their orbit for about 1 ms in order to discriminate their orbit categories. Such a classification is illustrated in figure 7 for the scenario B (average parameters) and the tangential injection R_tan = 65 cm and R_tan = 30 cm. The resulting classification gives us in particular the fraction of trapped particles for all configurations (scenarios and scan in R_tan). This fraction is plotted in figure 7 where we also indicate the average orbit width δ_r of the trapped population and their average Larmor radii.

Zoom In Zoom Out Reset image size

It is important to note that—in the axi-symmetric case—in considered scenarios, promptly ejected particles numbers are negligible when R_tan ⩾ 45 cm when CX processes are ignored. The losses can become however important for 0 ⩽ R_tan ⩽ 30 cm (co-perpendicular injection geometry). Table 2 gives the percentage of collisional steady-state losses for the 6 scenarios in the axi-symmetric case and compares it with the enhanced value obtained when the TFR is applied.

Before we move on to our detailed study of power deposition and losses when CX processes are introduced, we evaluate the NBI deposition when profiles are only the results of the Coulomb collision process when the TFR is applied. The set of profiles obtained for the total power deposited to plasma (excluding the final thermalization power) is given in figure 8 for all scenarios at R_tan = 0.65 m. We also indicate there the total amount of power deposited in the core (volume integration for positions inside ρ_t ≃ 0.45) and corresponding total power to the ions for all scenarios and co-injection tangency radii. We identify an optimum geometry at R_tan = 30 cm for scenarios C and E, R_tan = 45 cm for scenarios B, D and F (higher density) and R_tan = 65 cm for scenario A. The total power given to ions overall increases in all cases when we move to more tangential geometries, as core deposition is favorable to ion heating.

Zoom In Zoom Out Reset image size

It is found that, when we apply the CX reactions (e.g. in figure 11), the core deposited power is almost unchanged in all scenarios, all beam geometries and all realistic values of background neutral densities. As such, the right top panel in figure 8 is a realistic indication of the heating performance of each configuration.

Comparing the characteristic time of the CX process τ_CX with the thermalization time τ_s of the ions can give us an idea of how deeply the NBI deposition will be affected by the neutrals. This is illustrated in figure 9 where we can see that we can expect all ions with average radial position beyond ρ_t ∼ 0.85 to be neutralized at least once in scenarios C and F, beyond ρ_t ∼ 0.8 in scenarios B and E and beyond ρ_t ∼ 0.75 in scenarios A and D. This position thus tends to move inward with lower density, reflecting the depth of penetration of neutrals (n₀ ≳ 10¹⁵ m⁻³). When the ratio τ_s/τ_CX is beyond 2, only few NBI fast ions can exist in the considered radial region. Comparing the scenarios, in particular C (Ip = 2.0 MA) with F (Ip = 0.8 MA), we can see the beneficial influence of higher plasma current and thus increased edge plasma density: the amount of neutrals in the edge region is reduced so that τ_CX is locally increased.

Zoom In Zoom Out Reset image size

It is of interest to study the influence of the magnitude of the background neutral density on our results since the value obtained from equation (3) is empirical. Furthermore, the magnitude of the neutral background in SOL has been observed to be correlated to modifications of the SOL transport by changes in the ExB flows (see [15]): it is thus interesting to study what are the consequences of additional content of neutrals in the SOL.

In figure 9, we selected scenario B and the tangency radii of 30 cm and 65 cm in order to do a scan in boundary value of molecular deuterium. The CX induced losses is strongly correlated to the magnitude of the neutral density. However, beyond some value, the dependence becomes weaker due to the disappearance of all fast ions in the edge region and the impossibility for the neutrals to get too deep in the core plasma. We also quantify the change in the ion losses. The tangential injection (R_tan = 65 cm) is found to be more sensitive to a further increase in neutral density: this is partly explained by the larger amount of markers initially deposited in the edge region. Furthermore, re-ionized markers tend to end up on loss orbits: the ion losses also increase with neutral density unlike in the R_tan = 30 cm case.

The CX losses are analysed for all scenario and tangency radii in the range R_tan = 15 cm to 65 cm. The results are plotted in figure 10, on the left. Two trends of CX losses can be separated based on the amount of power deposited respectively in the pedestal and mid-edge region (see figure 8, right):

• Scenarios A, C and F have low values of τ_s/τ_CX w.r.t. 1 at ρ_t ≃ 0.8 (figure 9): their CX losses are dominated by the deposited power in the pedestal
• Scenarios B, D and E have large values of τ_s/τ_CX at ρ_t ≃ 0.8 of 1 or above: their CX losses are dominated by the deposited power in the mid-edge region

Zoom In Zoom Out Reset image size

The influence of density is seen when one considers the trend within the density scan (scenarios D–F) at the most tangential co-injection geometries (R_tan ⩾ 45 cm) for which there is a low amount of power deposited in the pedestal: the drop of temperature due to the increase of density lowers τ_s and the lower penetration of the neutrals increases τ_CX: this contributes to the overall decreasing CX losses with increasing density (τ_s/τ_CX at ρ_t ≃ 0.8 is lowered from 2.2 (scenario D) to 0.1 (scenario F) as seen in figure 9).

On the right of figure 10, we consider a scan in beam geometry. When TFR losses are important, the overall ion orbit losses are found to be reduced by the application of the CX processes (comparing black and red line in the right figure in 10 for scenario B). This is due to the large number of neutralizations in the edge region, preventing them to be lost on longer time scales because of the TFR. The figure indicates an achievable candidate for counter injection at about R_tan = −30 cm (≃38% losses). It is of interest that balanced injection can likely be achieved in COMPASS upgrade. Indeed, the promising QH-mode observed in low momentum plasmas in DIII-D could provide a viable solution to sustaining a pedestal in a fusion reactor [28] (provided its applicability at higher heating power can be extended).

When we consider the impact of the CX losses on the plasma, one might want to know how the significant percentage of power lost (figure 11) will impact the plasma in terms of heating and performances. We have thus implemented a 'toy model' in the transport solver METIS that reproduces the edge loss and core increase pattern of the ratio of deposited power with CX-losses and without (P_dep,CX/P_dep, noCX). We find similarities between our results and those obtained in other simulations (for DIII-D) using the SPIRAL code [29]. As can be noted on figure 11, whilst the power deposited in the edge drops significantly, the power deposited in the core of the plasma is slightly increasing. That is because some of the neutrals originating from the CX processes are re-ionized in a more central location.

Zoom In Zoom Out Reset image size

We consider the calculation for two NBI geometries in scenario B: one with tangential injection (R_tan = 65 cm) and one with co perpendicular injection (R_tan = 30 cm). The results of this simple simulations are shown in figure 11 in the right panel. We conclude that the overall drop in plasma stored energy and temperature is relatively small, even when a very large neutral density is considered. The simple transport model predicts a uniformly radially distributed drop in ion temperature of about 5%. The losses do not seem to impact significantly the edge temperature: the drop of heat source in the edge where the gradient is constrained by the energy stored in the pedestal is not causing a drop of edge temperature but increases the cooling rate of the edge and thus decreases the core temperature. Both core electron and ion temperatures are observed to drop when increasing the neutral density (which corresponds to an increase of the edge cooling rate). The slight increase in core heating does not compensate for the lost power. As was shown previously in figure 9, the impact of a very large background neutral density (n_0D2 = 4n₀ ≃ 8 × 10¹⁸ m⁻³ case) is more problematic in the tangential injection case whilst no further drop in terms of core T_i is obtained for the co-perpendicular (R_tan = 30 cm).

An estimate of the density of power deposited on the wall was calculated using a simplified 2D geometry of a preliminary design of the wall. The results are found in figure 12 where we have separated the incident fast neutrals and ions. Unlike the ions, the CX-born neutrals will mostly not be blocked by the outer wall limiter and will actually spread on a much larger area, that includes all of the openings in the vacuum vessel. A more detailed study of the actual power loss will be done in the future for design purposes using a complete 3D model of the components in the tokamak chamber. This is beyond the scope of this publication. We can however quantify an upper estimate for the heat flux on the limiter generated by 1 MW of NBI heating in scenario B: about 50 kW m⁻² (co-injection) and 200 kW m⁻² (counter-injection). These values remain well within the acceptable range of heat loads for tungsten based components (less than 1 MW m⁻² in total for 4 MW of NBI heating). These heat loads and their energy distribution could be used for modelling the sputtering occurring at the surface of the limiter (and other in-vessel components) in order to evaluate the beam-induced impurity source. The simulations will also support the design of fast ion loss detectors.

Zoom In Zoom Out Reset image size

Shoulder formation is related to a change of filament transport in the SOL from sheath-limited to inertial regime, which is characterized by a substantial increase in density (and to a lesser extent temperature) decay length in the far SOL [30] (meaning by this denomination the region of the SOL away from separatrix by more than one conventional decay length). Increased density decay length leads to increased transport to the first wall, which can in principle lead to higher heat fluxes on the plasma-facing components. It is on the other hand possibly favourable w.r.t. The power deposition length on the divertor [31].

One of the interest of the COMPASS upgrade tokamak will be its ability to study configurations that involve a far-SOL shoulder. We have constructed an approximate shoulder configuration by multiplying the density decay length by a factor 5 in the far-SOL (for SOL positions beyond one conventional decay length as obtained from 4).

We have applied a shoulder to scenario F, which has a large Greenwald fraction of ∼0.5. The resulting profiles are given in figure 13 where we compare them with the configuration that does not have a shoulder. We have applied the KN1D modelling to both configurations and obtained different neutral background profiles.

Zoom In Zoom Out Reset image size

We can then use the complete full-orbit modelling of EBdyna to study the difference in the CX processes because of the changes in the SOL. When a shoulder is present, the increased plasma density results in an outer ionization front away from the separatrix. This significantly reduces the neutral density in the pedestal. As a result, the edge CX losses are reduced by the presence of the shoulder.