Noise-Free Rapid Approach to Solve Kinetic Equations for Hot Atoms in Fusion Plasmas

Abstract

At the first wall of a fusion reactor, charged plasma particles are recombined into neutral molecules and atoms recycling back into the plasma volume where charge exchange (cx) with ions. As a result hot atoms with chaotically directed velocities are generated which can strike and erode the wall. An approach to solve the kinetic equation in integral form for cx atoms, being alternative to statistical Monte Carlo methods, has been speeded up by a factor of 50, by applying an approximate pass method to evaluate integrals, involving the ion velocity distribution function. It is applied to two-dimensional transfer of cx atoms near the entrance of a duct, guiding to the first mirror for optical observations. The energy spectrum of hot cx atoms, escaping into the duct, is calculated and the mirror erosion rate is assessed. Computations are done for a molybdenum first mirror under plasma conditions expected in the fusion reactor DEMO. Kinetic modeling results are compared with those found with a diffusion approximation valid in very cold and dense plasmas. For ducts at the torus outboard a more rigorous kinetic consideration predicts an erosion rate by a factor up to 2 larger than the diffusion approximation.

1. Introduction

In devices for thermonuclear fusion research, for example, of the Tokamak type, particles of hydrogen isotopes, deuterium and tritium, are in the form of a hot fully ionized plasma [1, 2]. To avoid a destruction of the machine wall, a special region, the so-called scrape-off layer (SOL), is arranged at the plasma edge, where particles stream along the magnetic field to the special target plates [3]. Normally, it is done by using additional magnetic coils to form a divertor configuration (see Figure 1a).

Figure 1.

The cross section of toroidally symmetric fusion device of the Tokamak type with the SOL region formed by the presence of a divertor (a) and the processes near the opening in the wall for a duct guiding to a first mirror (b).

By reaching the divertor target plates, plasma electrons and ions are recombined into neutral atoms and molecules which are finally exhausted from the device by pumps. However, in a future fusion reactor like DEMO [1, 2], only a minor fraction of 1% of neutrals generated at the targets will be pumped out. The rest of them is ionized again in the plasma near the targets. This “recycling” process significantly restrains the parallel plasma flow in the SOL [4]. Therefore, a considerable fraction of plasma particles lost from the plasma core will reach the vessel wall before they are exhausted into the divertor. Plasma fluxes to the wall saturate it with fuel particles in a time much shorter than the discharge duration, and a comparable amount of neutral species will recycle back from the wall into the plasma. Recycling neutrals are not confined by the magnetic field and penetrate at several centimeters into the SOL. Here, charge exchange (cx) collisions of them with ions generate atoms of energies much higher than that of primary recycling neutral particles. A noticeable fraction of such secondary cx atoms hit the vessel wall and erode it.

Statistical Monte Carlo methods [5] are normally used to model cx atoms at the edge of fusion devices. A crucial obstacle to apply these approaches for extensive parameter studies, for example, with the aim to optimize the duct geometry, has too long calculations needed to achieve reasonably small accident errors. This is, however, necessary, for example, to couple neutral parameters with the plasma calculations. In a one-dimensional geometry, an alternative approach, based on iteration procedure to solve the kinetic equation represented in an integral form, has been elaborated decades ago [6]. Being free from statistical noise and permitting the convergence of iterations to the error level defined by the machine accuracy, this method, nonetheless, is also time-consuming. The reason is the necessity to assess integrals in the velocity space from functions involving the ion velocity distribution function, additionally to integrations in the normal space. Recently [4], an approximate pass method has been applied to evaluate these integrals, and the acceleration of kinetic calculations by a factor of 50 has been achieved.

The amendments, outlined above, have allowed to perform calculations of the plasma parameters in the DEMO SOL with cx atoms described kinetically, by varying the input parameters, for example, the plasma transport characteristics, in a broad range [4]. In addition the results of these computations have been thoroughly compared with those obtained with cx species described in the so-called diffusion approximation. This approximation, often used in diverse edge modeling approaches to save CPU time (see, e.g., [7, 8]), is strictly valid under plasma conditions of low temperature and high density where the time between cx collisions of atoms with ions is much smaller than that till their ionization by electrons. Across the DEMO SOL, the plasma density and the temperatures of electrons and ions change, however, by orders of magnitude [4].

The usage of a diffusion approximation for cx atoms, generated from species recycling from the wall, becomes especially questionable by considering the situation near an opening in the vessel wall. Such openings will be made in a reactor for diverse purposes, for example, for ducts leading to the first mirrors, collecting light emitted by impurity species in the plasma (see Figure 1b). These installations are inaccessible for charged plasma particles, moving mostly along magnetic field lines and penetrating into the duct at a distance of 1 mm. However, hot cx atoms, unconfined by the magnetic field, can freely hit and erode these installations. The erosion rate is very sensitive to the energy spectrum of cx atoms which can be significantly dependent on the modeling approach. At the position of the duct opening, the inflow of recycling neutrals is actually absent, and the transfer of cx atoms has two- or even three-dimensional pattern. By calculating the density of cx atoms in the vicinity of a circular opening from a 2D diffusion equation, the erosion rate of a first mirror of Mo has been assessed in [9]. In the present paper, we extend the approaches, elaborated in [4] to model cx atoms kinetically in the 1D case, on a 2D geometry.

The results are compared with those of [9]. Although there cx atoms have been considered in the diffusion approximation to reduce CPU time, the new approach allows to perform more exact kinetic calculations even by orders of magnitude faster.

2. Basic equations

Although the concept of magnetic fusion is based on the idea that charged particles can be infinitely long confined within closed magnetic surfaces, there are diverse physical mechanisms leading to the losses across these surfaces (see, e.g., [10]). Therefore, electrons and ions escape through the separatrix into the SOL and may reach the first wall of the machine vessel. Here, these species fast recombine into neutral atoms, with an energy comparable with that of ions. Partly, these atoms escape immediately back into the plasma, and these species are referred as backscattered (bs) atoms. The rest transmits the energy to the wall particles and may be bounded into molecules, desorbing back into the plasma as the wall, and is saturated with the working gas. Here, the densities of recycling backscattered atoms and desorbed molecules decay with the distance from the wall because diverse elementary processes such as ionization, charge exchange and dissociation in the latter case, caused by collisions with the plasma electrons and ions, lead to disintegration of neutrals [11]. By the dissociation and ionization of molecules, the so-called Franck-Condon fc atoms are generated of a characteristic energy $E_{\mathit{fc}}$ of $3.5\mathrm{eV}$ [3, 11]. By charge exchange collisions of the primary neutral species, listed above, with the plasma ions, cx atoms of higher energies and with chaotically oriented velocities are generated. Although, namely, the latter are of our particular interest, to describe them firmly, all the species introduced have to be considered. In the vicinity of a circular opening in the wall (see Figure 1b), any neutral particle is characterized by the components $V_x$ and $V_{\rho }$ of its velocity, where $x$ and $\rho$ are distances from the wall and opening axis (see Figure 1b).

2.1. Recycling molecules and atoms

Recycling species, lunched from the wall, comprise backscattered atoms and desorbed molecules. Henceforth, we assume that the $x$ components of their velocities have single values $V_{\mathit{bs},m}$ and distinguish two groups of particles moving outward and toward the opening axis, with $V_{\rho }=V_{\mathit{bs},m}$ and $V_{\rho }=-V_{\mathit{bs},m}$, correspondingly. The magnitudes of $V_{\mathit{bs},m}$ are predefined by the generation mechanism for the species in question. Within the plasma the point with certain coordinates $x$ and $\rho$ can be reached only by particles, coming from the wall points $\left(0,\rho +x\right)$ for $V_{\rho }<0$ and $\left(0,|\rho -x|\right)$ for $V_{\rho }>0$, respectively. By taking into account the attenuation processes with neutrals in the plasma, one gets

$$n_{\mathit{bs},m}\left(x,\rho \right)=\frac{n_{\mathit{bs},m}\left(0,\rho +x\right)+n_{\mathit{bs},m}\left(0,\left|\rho -x\right|\right)}{2}\text{exp}\left(-\underset{0}{\overset{x}{\int }}\frac{{\nu }_{a,m}}{V_{\mathit{bs},m}}\mathit{dy}\right)$$E1

where ${\nu }_a=n\left(k_{\mathit{ion}}^a+k_{\mathit{cx}}^a\right)$ and ${\nu }_m=n\left(k_{\mathit{dis}}^m+k_{\mathit{ion}}^m+k_{\mathit{cx}}^m\right)$ are the decay frequencies for atoms and molecules, correspondingly, with $n$ being the plasma density assumed the same for electrons and ions, $k_{\mathit{ion}}^{a,m}$ are the rate coefficients of the ionization by electrons, $k_{\mathit{cx}}^{a,m}$ those for the charge exchange with ions, and $k_{\mathit{dis}}^m$ that for the dissociation of molecules.

2.2. Franck-Condon (fc) atoms

The fc atoms are born by the destruction of molecules and have positive and negative values both of $V_{\rho }$ and of $V_x$. Henceforth, we distinguish four groups of these species with the velocity components $\left(V_x,V_{\rho }\right)=\left(\pm V_{\mathit{fc}},\pm V_{\mathit{fc}}\right)$, where $V_{\mathit{fc}}=\sqrt{E_{\mathit{fc}}/m}$ and $m$ is the atom mass. The source density of fc atoms of each group is as follows:

$$S_{\mathit{fc}}=n\left(2k_{\mathit{dis}}^m+k_{\mathit{ion}}^m+k_{\mathit{cx}}^m\right)n_m/4$$

and depends on $\rho$ through $n_m$. The particle densities $n_{\mathit{fc}\pm \pm }$, where the first subscript $\pm$ corresponds to the sign of $V_x$ and the second one – of $V_{\rho }$, change along the characteristics $\pm x\mp \rho =\mathit{\text{const}}$ according to the continuity equation:

$$\sqrt{2}{V}_x{\partial }_l{n}_{\mathit{fc}\pm \pm }=S_{\mathit{fc}}\left(l\right)-{\nu }_a{n}_{\mathit{fc}\pm \pm }$$

where $l$ is the length of the corresponding characteristics. The boundary conditions take into account that fc atoms are reflected with the probability $R_{\mathit{fc}}$ from the wall, $x=0$, and are absent far from the wall, $x\to \infty$. For the total density of fc atoms, $n_{\mathit{fc}}=n_{\mathit{fc}++}+n_{\mathit{fc}+-}+n_{\mathit{fc}-+}+n_{\mathit{fc}--}$, one can obtain

$$\begin{array}{l}{n}_{\mathit{fc}}\left(x,\rho \right)=\underset{0}{\overset{\infty }{\int }}\frac{S_{\mathit{fc}}\left(y,\left|\rho -x+y\right|\right)+S_{\mathit{fc}}\left(y,\left|\rho +x-y\right|\right)}{V_{\mathit{fc}}}{\mu }_0\mathit{dy}+\\ +\frac{R_{\mathit{fc}}}{2}\Theta \left(\rho +x-{\rho }_0\right)\underset{0}{\overset{\infty }{\int }}\frac{S_{\mathit{fc}}\left(y,\rho +x+y\right)+S_{\mathit{fc}}\left(y,\left|\rho +x-y\right|\right)}{V_{\mathit{fc}}}{\mu }_1\mathit{dy}+\\ +\frac{R_{\mathit{fc}}}{2}\Theta \left(\left|\rho -x\right|-{\rho }_0\right)\underset{0}{\overset{\infty }{\int }}\frac{S_{\mathit{fc}}\left(y,\left|\rho -x\right|+y\right)+S_{\mathit{fc}}\left(y,\left|\left|\rho -x\right|-y\right|\right)}{V_{\mathit{fc}}}{\mu }_1\mathit{dy}\end{array}$$E2

with

$${\mu }_0=\text{exp}\left(-\frac{\left|U_x-U_y\right|}{V_{\mathit{fc}}}\right),{\mu }_1=\text{exp}\left(-\frac{U_x+U_y}{V_{\mathit{fc}}}\right),U_{x,y}=\underset{0}{\overset{x,y}{\int }}{\nu }_a\mathit{dz},$$

the Heaviside function $\Theta \left(\xi \ge 0\right)=1,\Theta \left(\xi <0\right)=0$ and ${\rho }_0$ being the opening radius (see Figure 1b).

2.3. Charge exchange cx atoms

The velocity distribution function of cx atoms, $f_{\mathit{cx}}\left(x,\rho ,V_x,V_{\rho }\right)$, is governed by the kinetic equation:

$$V_x{\partial }_x{f}_{\mathit{cx}}+\frac{V_{\rho }}{\rho }{\partial }_{\rho }\left(\rho f_{\mathit{cx}}\right)=\frac{S_{\mathit{cx}}}{\pi V_{\mathit{th}}^2}\text{exp}\left(-\frac{V_x^2+V_{\rho }^2}{V_{\mathit{th}}^2}\right)-{\nu }_a{f}_{\mathit{cx}},$$E3

Here, $S_{\mathit{cx}}=S_{\mathit{cx}}^0+S_{\mathit{cx}}^1$ is the total density of the source of cx atoms, with $S_{\mathit{cx}}^0=n\left[k_{\mathit{cx}}^m{n}_m+k_{\mathit{cx}}^a\left(n_{\mathit{bs}}+n_{\mathit{fc}}\right)\right]$ and $S_{\mathit{cx}}^1={\mathit{nk}}_{\mathit{cx}}^a{n}_{\mathit{cx}}$ being the contributions due to charge exchange collisions with ions of primary neutrals and cx atoms themselves, correspondingly, where

$$n_{\mathit{cx}}=\underset{-\infty }{\overset{\infty }{\int }}\underset{-\infty }{\overset{\infty }{\int }}{f}_{\mathit{cx}}{\mathit{dV}}_x{\mathit{dV}}_{\rho };$$

is the density of the latter; the velocity distribution of the source is assumed as the Maxwellian one with the ion temperature $T_i$, and $V_{\mathit{th}}=\sqrt{2T_i/m}$ is the ion thermal velocity.

2.3.1. Integral-differential equations for the density of cx atoms

To solve Eq. (3), $V_{\rho }$ is discretized as $\pm V_{\rho }^l$, with $V_{\rho }^l=\mathrm{\Delta }V\left(l-1/2\right)$, $\mathrm{\Delta }V=V_{\max }/l_{\max }$, and $l=1,\cdots ,l_{\max }$. The $V_x$ distribution functions ${\phi }_l^{\pm }$ of particles with $V_{\rho }$ in the range $\pm V_{\rho }^l-\mathrm{\Delta }V/2\le V_{\rho }\le \pm V_{\rho }^l+\mathrm{\Delta }V/2$ are governed by the equations, following from the integration of Eq. (3) over these ranges:

$$V_x{\partial }_x{\phi }_l^{\pm }\pm \frac{V_{\rho }^l}{\rho }{\partial }_{\rho }\left({\rho \phi }_l^{\pm }\right)=\frac{S_{\mathit{cx}}{\delta }_l}{2\sqrt{\pi }{V}_{\mathit{th}}}\text{exp}\left(-\frac{V_x^2}{V_{\mathit{th}}^2}\right)-{\nu }_a{\phi }_l^{\pm },$$E4

with

$${\delta }_l=\left\{\text{erf}\left[-\left(\frac{V_{\rho }^l+\mathrm{\Delta }V/2}{V_{\mathit{th}}}\right)\right]-\text{erf}\left[-\left(\frac{V_{\rho }^l-\mathrm{\Delta }V/2}{V_{\mathit{th}}}\right)\right]\right\}/\text{erf}\left(\frac{V_{\max }}{V_{\mathit{th}}}\right).$$

From the equation above, one gets the following ones for the variables ${\phi }_l={\phi }_l^{+}+{\phi }_l^{-}$ and ${\gamma }_l=V_{\rho }^l\left({\phi }_l^{+}-{\phi }_l^{-}\right)$:

$$V_x{\partial }_x{\phi }_l+\frac{1}{\rho }{\partial }_{\rho }\left({\rho \gamma }_l\right)=\frac{S_{\mathit{cx}}{\delta }_l}{\sqrt{\pi }{V}_{\mathit{th}}}\text{exp}\left(-\frac{V_x^2}{V_{\mathit{th}}^2}\right)-{\nu }_a{\phi }_l,$$E5

$${\partial }_x{\gamma }_l+\frac{{\left(V_{\rho }^l\right)}^2}{V_x\rho }{\partial }_{\rho }\left({\rho \phi }_l\right)=-\frac{{\nu }_a}{V_x}{\gamma }_l,$$E6

The latter equation is straightforwardly integrated with respect to $x$, and for the $\rho$-component of the flux density, one gets

$${\gamma }_l=-\frac{{\left(V_{\rho }^l\right)}^2}{V_x\rho }\underset{x_0}{\overset{x}{\int }}{\partial }_{\rho }\left({\rho \phi }_l\right)\text{exp}\left(\frac{U_y-U_x}{V_x}\right)\mathit{dy}.$$

Here, $x_0$ is the position, where the boundary condition is imposed; since there is no influx of cx atoms from the wall, ${\gamma }_l\left(x=0\right)=0$ for $V_x>0$, neutral species are absent deep enough in the plasma and ${\gamma }_l\left(x\to \infty \right)\to 0$ for $V_x<0$. The found expression for ${\gamma }_l$ can be made more convenient if (i) the $x$-variation of terms involved is approximated by linear Taylor series:

$${\partial }_{\rho }\left({\rho \phi }_l\right)\left(y\right)\approx {\partial }_{\rho }\left({\rho \phi }_l\right)\left(x\right)+{\partial }_x{\partial }_{\rho }\left({\rho \phi }_l\right)\left(x\right)\cdot \left(y-x\right),U_y\approx U_x+{\nu }_a\left(x\right)\cdot \left(y-x\right),$$

and (ii) by taking into account that the width of the region, occupied by cx atoms, is of several mean free path lengths (MFPL), that is, $\left|x-x_0\right|>\left|V_x\right|/{\nu }_a$ typically. As a result we get

$${\gamma }_l\approx -\left(D_l/\rho \right)\times {\partial }_{\rho }\left({\rho \phi }_l\right),$$

with $D_l={\left(V_{\rho }^l\right)}^2/{\nu }_a$.

By substituting the last expression into Eq. (5), this is reduced to the following one:

$$V_x{\partial }_x{\phi }_l-\frac{D_l}{\rho }{\partial }_{\rho }^2\left({\rho \phi }_l\right)=\frac{S_{\mathit{cx}}{\delta }_l}{\sqrt{\pi }{V}_{\mathit{th}}}\text{exp}\left(-\frac{V_x^2}{V_{\mathit{th}}^2}\right)-{\nu }_a{\phi }_l.$$E7

Equation (7) can be integrated as the first-order equation with respect to $x$, with the boundary conditions similar to those for ${\gamma }_l$, that is, ${\phi }_l\left(x=0\right)=0$ for $V_x>0$ and ${\phi }_l\left(x\to \infty \right)\to 0$ for $V_x<0$. Finally, for the density of cx atoms within the range $V_{\rho }^l-\mathrm{\Delta }V/2\le \left|V_{\rho }\right|\le V_{\rho }^l+\mathrm{\Delta }V/2$, ${\eta }_l=\underset{-\infty }{\overset{\infty }{\int }}{\phi }_l{\mathit{dV}}_x$, one obtains the following integral-differential equation:

$$-\frac{D_l}{{\nu }_a\rho }{\partial }_{\rho }^2\left({\rho \eta }_l\right)=\underset{0}{\overset{\infty }{\int }}\frac{S_{\mathit{cx}}{\delta }_l}{\sqrt{\pi }{V}_{\mathit{th}}}{I}_{\alpha }\mathit{dy}-{\eta }_l,$$E8

where $I_{\alpha }=\underset{0}{\overset{\infty }{\int }}{F}_{\alpha }\mathit{du}$, with $F_{\alpha }\left(u\right)=u^{-1}\text{exp}\left(-u^2-\alpha /u\right)$ and $\alpha =\left|U_y-U_x\right|/V_{\mathit{th}}\left(y\right)$. For the total density of cx atoms, one has $n_{\mathit{cx}}={\sum }_{l=1}^{l_{\max }}{\eta }_l$.

2.3.2. Diffusion approximation

By neglecting the $\rho$ derivative and by considering a single $V_{\rho }$ value, Eq. (8) is reduced to an integral one obtained in the 1D case (Eq. (16) in Ref. [6] or Eq. (11) in Ref. [4]). In this case $l\equiv 1,{\delta }_l=1$ and $n_{\mathit{cx}}={\eta }_1$. As it was demonstrated in Ref. [4], this integral equation can be reduced to an ordinary differential diffusion equation, if (i) the electron temperature is low enough and the ionization rate is much smaller than that of the charge exchange, $k_{\mathit{ion}}^a\ll k_{\mathit{cx}}^a$, and (ii) the characteristic dimension for the plasma parameter change is much larger than the typical MFPL for cx atoms $V_{\mathit{th}}/{\nu }_a$. The same procedure can be performed in the case under the consideration, providing the following 2D diffusion equation:

$$-\frac{D_{\mathit{th}}}{\rho }{\partial }_{\rho }\left(\rho {\partial }_{\rho }{n}_{\mathit{cx}}\right)-{\partial }_x\left[\frac{{\partial }_x\left(n_{\mathit{cx}}{T}_i\right)}{{\nu }_{a}m}\right]=S_{\mathit{cx}}^0-k_{\mathit{ion}}^a{\mathit{nn}}_{\mathit{cx}},$$E9

with $D_{\mathit{th}}=V_{\mathit{th}}^2/\left(2{\nu }_a\right)$. The boundary conditions to Eq. (9) imply as follows: (i) cx atoms leave the plasma with the velocity close to the ion thermal one, and (ii) neutral species are absent far from the wall, that is, $n_{\mathit{cx}}\left(x\to \infty \right)=0$. Additionally, ${\partial }_{\rho }{\eta }_l=0$ and ${\partial }_{\rho }{n}_{\mathit{cx}}=0$ for $\rho =0$ and $\rho \to \infty$ for Eqs. (8) and (9), respectively.

2.3.3. Assessment of the velocity space integral $I_{\alpha }$

By inspecting Eq. (8), one can see the cause for large calculation time for kinetic computations for cx atoms. Even in the 1D case, for each grid point, one has to calculate enclosed double integrals over the ion velocity space and over the whole plasma volume, $0\le y\le r_w$, and repeat this, in an iterative procedure, with respect to the whole profiles of $n_{\mathit{cx}}\left(x\right)$. In a 2D situation, the runs through the grid points in the $\left(\rho ,V_{\rho }\right)$ phase subspace are involved additionally into calculations. These can be, however, efficiently parallelized. The integral $I_{\alpha }$ is appeared due to the generation of cx atoms with the ion Maxwellian velocity distribution and consequent attenuation by ionization and cx collisions. With a high accuracy, it can be assessed by an approximate pass method [12]. The function $F_{\alpha }\left(u\right)$ is shown in Figure 2a for $\alpha =0.3$, $1$, and $3$. One can distinguish four $u$ intervals, where $F_{\alpha }\left(u\right)$ behaves principally differently: $u\le u_1$, $u_1<u\le u_m$, $u_m<u\le u_2$, and $u_2<u$. The characteristic points $u_{1,2}$ and $u_m$ are defined by the conditions $F_{\alpha }^{\prime }\left(u_m\right)=0$ and $F_{\alpha }^{″}\left(u_{1,2}\right)=0$ for the derivatives $F_{\alpha }^{\prime }=-F_{\alpha }{g}_1/u^2$ and $F_{\alpha }^{″}=F_{\alpha }{g}_2/u^4$ with

$$g_1\left(u\right)=2u^3+u-\alpha ,g_2\left(u\right)={\left(g_1+u\right)}^2-\left(2+6u^2\right)u^2.$$

Figure 2.

The functions $F_{\alpha }\left(u\right)=u^{-1}\text{exp}\left(-u^2-\alpha /u\right)$ for $\alpha =1$ (solid curve), $\alpha =0.3$ (dashed curve), and $\alpha =3$ (dashed-dotted curve) (a); the integral $I_{\alpha }=\underset{0}{\overset{\infty }{\int }}{F}_{\alpha }\left(u\right)\mathit{du}$ estimated by the pass method (solid curve) and computed numerically (dashed curve) (b).

For $u_m$ one can use the Cardano formula for the real root of a cubic equation $u_m=\sqrt[3]{b+a}-\sqrt[3]{b-a}$ with $a=\alpha /4$ and $b=\sqrt{a^2+1/216}$. By searching for $u_{1,2}$, we notice that for the roots of interest the equation $g_2\left(u_{1,2}\right)=0$ reduces to the following ones:

$$g_1\left(u_1\right)+u_1+u_1\sqrt{2+6u_1^2}=0,g_1\left(u_2\right)+u_2-u_2\sqrt{2+6u_2^2}=0.$$

With properly selected initial approximations accurate enough, solutions of these equations are found after three to five iterations with the Newton approach.

To assess the integrals $I_{\alpha }$, we approximate function $F_{\alpha }\left(u\right)$ as a linear one at $u\le u_1$ and by polynomials of the fifth order at the intervals $u_1<u\le u_m$ and $u_m<u\le u_2$. The coefficients are selected to reproduce $F_{\mathit{\alpha m}\mathrm{,1,2}}=F_{\alpha }\left(u_{m\mathrm{,1,2}}\right)$, $F_{\alpha 1,2}^{\prime }=F_{\alpha }^{\prime }\left(u_{1,2}\right)$, $F_{\mathit{\alpha m}}^{″}=F_{\alpha }^{″}\left(u_m\right)$, and $F_{\alpha }^{\prime }\left(u_m\right)=F_{\alpha }^{″}\left(u_{1,2}\right)=0$. In addition, for $u>u_2$ the factor $\mathit{exp}\left(-u^2\right)$ leads to a very fast decay of $F_{\alpha }$ with increasing $u$. Therefore, for the aim to assess the corresponding contribution to $I_{\alpha }$, we approximate $F_{\alpha }$ by $\text{exp}\left(-u^2-\alpha /u_2\right)/u_2$. This results in the following expression:

$$\begin{array}{l}{I}_{\alpha }\approx \frac{F_{\alpha 1}^2}{2F_{\alpha 1}^{\prime }}+\frac{F_{\alpha 1}{\delta }_1+F_{\mathit{\alpha m}}\left(u_2-u_1\right)+F_{\alpha 2}{\delta }_2}{2}+\\ +\frac{F_{\alpha 1}^{\prime }{\delta }_1^2-F_{\alpha 2}^{\prime }{\delta }_2^2}{10}+F_{\mathit{\alpha m}}^{″}\frac{{\delta }_1^3+{\delta }_2^3}{120}+\frac{\sqrt{\pi }}{2}\frac{1-{\mathit{\text{erfu}}}_2}{u_2}\text{exp}\left(-\frac{\alpha }{u_2}\right),\end{array}$$E10

where ${\delta }_{1,2}=\left|u_m-u_{1,2}\right|$. Figure 2b shows $I_{\alpha }$ versus $\alpha$ found with the pass method, outlined above and evaluated numerically. By approximating the numerical results very accurately, the approximate pass method procedure allows to reduce the CPU time by a factor of 50 compared to direct estimates.

Alternatively to the usage of formula (10), one can calculate the integral $I_{\alpha }$ in advance and save the results in a table. Then, for any particular $\alpha$, the integral $I_{\alpha }$ is interpolated from the values in this table. However, some time is needed to find the proper $\alpha$ interval, and this time grows up with $\alpha \to 0$ since $I_{\alpha }\to \infty$ as $-\text{ln}\alpha$. Thus, the density of the data in the table has to be increased as $\alpha \to 0$, and it is not obvious that such a way is more time-saving than formula (10).

2.3.4. Numerical procedure

The procedure to solve Eq. (9) numerically is organized as follows: For the problem in question, with the plasma profiles assumed unaffected by the processes in the vicinity of the duct opening, the integral $I_{\alpha }\left(x,y\right)$ can be assessed once in advance. Then,

1. the source density $S_{\mathit{cx}}\left(x,\rho \right)$ is calculated, in the first approximation with $n_{\mathit{cx}}=0$;

2. for all x the first term on the right-hand side of the equation is computed as a function of $\rho$, and the second-order ordinary Eq. (8) is solved by the method outlined in [10] for all $V_l$;

3. the next approximation for $n_{\mathit{cx}}\left(x,\rho \right)$ and $S_{\mathit{cx}}\left(x,\rho \right)$ is computed, and the procedure is repeated till the relative error in, for example, $n_{\mathit{cx}}\left(0,0\right)$, becomes smaller than a desirable one.

A typical behavior of the error with the iteration number for calculations presented in the next section is shown in Figure 3. One can see that with consequent iterations the calculation error reduces steadily without any noise and the machine accuracy can be actually achieved.

Figure 3.

The reduction of the relative error in calculation with the iteration number.

Normally, a diffusion approximation is used to describe cx atoms to reduce CPU. The present approach to solve kinetic equation makes this unnecessary. To get the density of cx atoms, $n_{\mathit{cx}}\left(x,\rho \right)$, with a relative error below ${10}^{-6}$, kinetic calculations require much less CPU time than it is needed to solve the diffusion equation. Thus, for 500 and 30 grid points in the $x$ and $\rho$ directions, respectively, $l_{\max }$ =10 and $V_{\max }=3V_{\mathit{th}}$, the CPU times for a 1 GHz processor are 30 s and 112 min, correspondingly. Of course, the time for kinetic calculations increases significantly if neutrals affect the plasma background, and $I_{\alpha }\left(x,y\right)$ has to be evaluated by each iteration. Even in this case kinetic computations for cx atoms are done faster by factor of 4.

3. Results of calculations

Figure 4 shows the profiles of the plasma parameters, assumed homogeneous on the flux surfaces, across the surfaces, computed in [13] with the input parameters from the European DEMO project [1, 2]. The distance $x$ from the first wall to a certain flux surface depends on the poloidal angle $\theta$ (see Figure 1a) and in Figure 4 $x$ corresponds to the torus outboard, $\theta =0$, where this distance is minimal. The profiles of the cx atom density found with these parameters and for the duct radius ${\rho }_0=0.05\mathrm{m}$, at the opening axis, $\rho =0$, and far from it, $\rho =0.25\mathrm{m}$ are shown in Figure 5. There is a noticeable difference between the profiles of $n_{\mathit{cx}}$ calculated in the diffusion approximation (5a) and kinetically (5b). In particular, the difference between the density profiles at these two positions is significantly larger in the diffusion approximation than for those computed kinetically.

Figure 4.

The profiles of the plasma density (solid curve), the ion (dashed curve), and electron (dashed-dotted curve) temperatures versus the distance from the wall at the DEMO torus outboard.

Figure 5.

The profiles of the cx atom densities far from the opening (solid curve) and at its axis (dashed curve) computed in the diffusion approximation (a) and kinetically (b).

In the latter case, $n_{\mathit{cx}}$ is still noticeable in much deeper and hotter plasma regions. This circumstance is of importance for the erosion of installations in ducts due to physical sputtering because this is very sensitive to the energy of impinging particles. To demonstrate this the erosion rate of a Mo mirror, imposed into the duct of three different radii ${\rho }_0$, at the distance $h$ from its opening (see Figure 1b) is assessed. To do such assessment, consider an infinitesimally thin toroid within the plasma, with a square cross section of the width $\mathit{d\rho }$, thickness $\mathit{dx}$, and the radius $\rho$, situated at the distance $x$ from the wall. Cx atoms with the energies within the range $E,E+\mathit{dE}$ are generated in the toroid with the rate:

$${\mathit{dR}}_{\mathit{cx}}\left(x,\rho ,E\right)=S_{\mathit{cx}}\left(x,\rho \right){\phi }_i\left(x,E\right)2\mathit{\text{πρdρdxdE}}.$$E11

where

$${\phi }_i\left(x,E\right)=\frac{2\sqrt{E}}{\sqrt{\pi }{T}_i^{1.5}\left(x\right)}\text{exp}\left[-\frac{E}{T_i\left(x\right)}\right].$$E12

is the Maxwellian energy distribution function of ions.

Since cx atoms move in all directions, some of them hit the mirror inside the duct (see Figure 1b). Henceforth, we analyze the surface position at the duct axis. One can see that cx atoms generated only in toroids with $\rho \le {\rho }_{\max }={\rho }_0\left(1+x/h\right)$ can hit this mirror point. The contribution of the plasma volume in question to the density of the cx atom flux perpendicular to the mirror surface is

$${\mathit{dj}}_{\mathit{cx}}\left(x,\rho ,E\right)={\mathit{dR}}_{\mathit{cx}}\frac{s\cdot \text{exp}\left(-\lambda /s\right)}{4\pi \left[{\left(x+h\right)}^2+{\rho }^2\right]},$$E13

where $s={\left[1+{\rho }^2/{\left(x+h\right)}^2\right]}^{-1/2}$ is the cosine of the atom incidence angle with respect to the duct axis. The exponential factor in Eq. (13), with $\lambda \left(x,E\right)=U_x/\sqrt{2E/m}$, takes into account the destruction of cx atoms in ionization and charge exchange collisions on their way through the plasma. By ionization the atom disappears. By charge exchange a new cx atom is generated. The latter process is taken into account by the contribution $S_{\mathit{cx}}^1$ in the source density $S_{\mathit{cx}}$.

The energy spectrum of cx atoms, hitting the mirror, is characterized by their flux density $\gamma =\int {\mathit{dj}}_{\mathit{cx}}/\mathit{dE}$ in the energy range $\mathit{dE}$, where the integration is performed over $\rho$ and $x$. By proceeding from $\rho$ to $s$, according to the relation $\rho =\left(x+h\right)\sqrt{1/s^2-1}$, one gets

$$\gamma \left(h,{\rho }_0,E\right)=\underset{0}{\overset{r_w}{\int }}{\phi }_i\left(x,E\right)\mathit{dx}\underset{h/\sqrt{h^2+{\rho }_0^2}}{\overset{1}{\int }}\frac{S_{\mathit{cx}}}{2}\text{exp}\left(-\frac{\lambda }{s}\right)\mathit{ds}.$$E14

The density of the outflow of mirror particles eroded by physical sputtering with cx atoms is as follows:

$${\Gamma }_{\mathit{sp}}\left(h,{\rho }_0\right)=\underset{0}{\overset{\infty }{\int }}\mathit{dE}\underset{0}{\overset{r_w}{\int }}{\phi }_i\left(x,E\right)\mathit{dx}\underset{h/\sqrt{h^2+{\rho }_0^2}}{\overset{1}{\int }}\frac{S_{\mathit{cx}}}{2}\text{exp}\left(-\frac{\lambda }{s}\right)Y_{\mathit{sp}}\left(E,s\right)\mathit{ds}$$E15

Here, $Y_{\mathit{sp}}$ is the sputtering yield, whose dependence on $E$ and $s$ is calculated by applying semiempirical formulas from Ref. [14]. The erosion rate, measured henceforth in $\mathit{nm}$ per full-power year $\left(\mathit{nm}/\mathit{pfy}\right)$, is $h_{\mathit{sp}}={\Gamma }_{\mathit{sp}}/n_{\mathit{sp}}$, with the particle density of molybdenum $n_{\mathit{sp}}$. It is shown in Figure 6 versus the distance $h$ between the wall position and the mirror point in question. One can see that in the former case the erosion rate is by a factor of 2 larger than in the latter one. It is of importance to notice that $h_{\mathit{sp}}$, found with $n_{\mathit{cx}}$ computed kinetically, does not unavoidably excides that obtained with $n_{\mathit{cx}}$ calculated in the diffusion approximation. The results above have been gotten for a mirror positioned in a duct at the plasma outboard, that is, at the largest major radius $R$ (see Figure 1a). Here, the local gradients of the plasma parameters are the largest because the distance between two particular flux surfaces has the minimum. Therefore, cx atoms penetrate, before they are ionized, into plasma regions with higher ion temperatures. The situation is different at the torus top where the corresponding distance has maximum and neutral species are attenuated already at significantly low $n$ and $T_i$. This leads to noticeably smaller $h_{\mathit{sp}}$ (see Figure 7). In this case, oppositely to the situation at the outboard, $h_{\mathit{sp}}$ is larger if $n_{\mathit{cx}}$ is computed in the diffusion approximation.

Figure 6.

The erosion rate of a Mo mirror versus the distance from the wall position to the mirror surface for ducts positioned at the torus outboard computed with the density of cx atoms found in the diffusion approximation (a) and kinetically (b).

Figure 7.

The erosion rate of a Mo mirror versus the distance from the wall position to the mirror surface for ducts positioned at the torus top computed with the density of cx atoms found in the diffusion approximation (a) and kinetically (b). (Note that the $h_{\mathit{sp}}$ scale is different than in Figure 6).

Due to technical requirements, an acceptable mirror erosion rate in a future fusion reactor should not exceed $1\mathit{nm}/\mathit{pfy}$. As one can see in the case of a duct at the torus outboard, this target level is exceeded significantly for all $h$ and ${\rho }_0$ under consideration. Seeding of the working gas into the duct is considered as a possible way to diminish the erosion rate below the maximum level allowed. The energy of cx atoms, coming into the duct, is reduced through elastic collisions with gas molecules, before cx atoms hit the mirror. Elastic collisions between neutral species lead to scattering on large angles. Consider a cx atom which, without gas in the duct, can hit the mirror directly. If the gas is seeded, in a narrow duct in question, with ${\rho }_0\ll h$, any collision of the cx atom with gas molecules leads to such a change of the atom velocity that with an overwhelming probability the atom will many times strike the duct wall before it gets the mirror. Through the collisions with the wall, the cx atom loses its energy so dramatically that it cannot contribute to mirror erosion. The reduction of the flux $j_E$ of incoming cx atoms with the energy $E$ by collisions with the gas in the duct is governed by the equation:

$${\mathit{dj}}_E/\mathit{dl}=-{\sigma }_{\mathit{el}}\left(E\right)n_g{j}_E,$$

where ${\sigma }_{\mathit{el}}\approx 3.8\cdot {10}^{-19}{E}^{-0.14}\left[m^2,\mathit{eV}\right]$ is the cross section for elastic collisions between cx atoms and molecules of hydrogen isotopes [15]. Consequently, the density of the outflow of mirror particles eroded by physical sputtering with cx atoms is modified, compared with Eq. (15), as follows:

$${\Gamma }_{\mathit{sp}}\left(h,{\rho }_0\right)=\underset{0}{\overset{\infty }{\int }}\mathit{dE}\underset{0}{\overset{r_w}{\int }}{\phi }_i\mathit{dx}\underset{h/\sqrt{h^2+{\rho }_0^2}}{\overset{1}{\int }}\frac{S_{\mathit{cx}}}{2}\text{exp}\left(-\frac{\lambda +{\sigma }_{\mathit{el}}{n}_{g}h}{s}\right)Y_{\mathit{sp}}\left(E,s\right)\mathit{ds}$$E16

In Figure 8, the calculated dependences of the erosion rate $h_{\mathit{sp}}$ on $h$ and ${\rho }_0$ are shown for several magnitudes of the density $n_g$ of the working deuterium gas in the mirror duct. One can see that the enhancement of $n_g$ above a level of $2\cdot {10}^{19}{\mathrm{m}}^{-3}$ should lead to the reduction of $h_{\mathit{sp}}$ below the target level of $1\mathit{nm}/\mathit{pfy}$. The question, to what extent the local plasma parameters may be changed by the outflow of the gas from the duct, has to be investigated in the future on the basis of approaches developed in [16]. There, it has been demonstrated that the ionization of the gas outflowing into the SOL can lead to dramatic growth of the local density and cooling of the plasma to a temperature of $1\mathrm{eV}$. Such cold dense plasma cloud can affect the transfer of cx atoms in the plasma near the opening in the wall.

Figure 8.

The erosion rate of a Mo mirror versus the distance from the wall position to the mirror surface for ducts positioned at the torus outboard computed with cx atoms treated kinetically for different gas densities in the duct: $n_g=0$ (solid line), $3\cdot {10}^{17}{\mathrm{m}}^{-3}$ (dashed line), $3\cdot {10}^{18}{\mathrm{m}}^{-3}$ (dotted line), ${10}^{19}{\mathrm{m}}^{-3}$ (dashed-dotted line), and $3\cdot {10}^{19}{\mathrm{m}}^{-3}$ (dashed-double dotted line).

4. Conclusions

The iteration approach to solve 1D kinetic equation for cx atoms, proposed decades ago [6], has been elaborated further to describe the transport of these species in a 2D geometry, in the vicinity of a circular opening in the wall of a fusion reactor. Unlike the Monte Carlo methods, this approach does not generate statistical noise so that calculation errors can be reduced to the level restricted by the machine accuracy. In order to perform calculations for a broad range of input parameters and do a thorough comparison with the results, obtained in the diffusion approximation for cx atoms [13], the solving procedure has been accelerated by a factor of 50, by applying an approximate pass method to assess integrals in the velocity space from functions, involving the Maxwellian velocity distribution of plasma ions.

The found possibility to speed up kinetic calculations is of importance, in particular, to perform firm assessments of the erosion rate of the first mirrors in future fusion reactors like DEMO. For a mirror located at the torus outboard, more accurate kinetic calculations predict by a factor of 2 higher erosion rate than the approximate diffusion approach. The erosion rate can be reduced very strongly either by putting the mirror duct at the torus top or by seeding the working gas into the duct. In the latter case, the elastic collisions with molecules in the gas reduce significantly the fraction of cx atoms which can hit and erode the mirror.

Nomenclature

Dl,th	Diffusivity of cx atoms in kinetic and diffusion descriptions, correspondingly
E	Energies of ions and atoms
Fαu≡1uexp−u2−αu,Iα≡∫0∞Fαudu
h	Distance from duct opening to the mirror
hsp	Mirror erosion rate in nm per full-power year
γhρ0E	Flux density of cx atoms onto the mirror
Γhρ0	Outflow density of sputtered mirror particles
kcxa,m	Charge exchange rate coefficients
kdism	Molecule dissociation rate coefficient
kiona,m	Ionization rate coefficients
λxE=Ux/2E/m
m	Mass of atoms
n	Plasma density
nm,bs,fc,cx	Densities of different neutral species—desorbed molecules, backscattered atoms, Franck-Condon, and charge exchanged atoms, respectively
νm,a	Decay frequencies for molecules and atoms
RcxxρE	The generation rate of cx atoms in the energy range E,E+dE within an infinitesimally thin plasma toroid of the width dρ, thickness dx, the radius ρ, and situated at the distance x from the wall
ρ	Distance from the opening axis
ρ0	The opening radius
s	Cosine of the cx atom incidence angle with respect to the duct axis
Sfc,cx	Source densities of fc and cx atom species
Scx0,1	cx atom source density contributions due to charge exchange of primary neutral species and cx atoms
σel	Cross section for elastic collisions between cx atoms and molecules of hydrogen isotopes
Te,i	Electron and ion temperatures
Ux,y≡∫0x,yνadz
Vm,bs,fc	Absolute velocities of neutral species
Vx,ρ	Velocity components perpendicular to the wall and to the opening axis, respectively
Vth=2Ti/m	Thermal velocity of ions
φixE	Energy distribution function of ions
x	Distance from the wall