2. Basic equations

special target plates [3]. Normally, it is done by using additional magnetic coils to form a

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

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

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).

divertor configuration (see Figure 1a).

4 Plasma Science and Technology - Basic Fundamentals and Modern Applications

vessel wall and erode it.

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 Efc of 3:5 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 Vx and V<sup>r</sup> of its velocity, where x and r are distances from the wall and opening axis (see Figure 1b).

Sfc <sup>¼</sup> <sup>n</sup> <sup>2</sup>km

ffiffiffi 2

∞

0

∞

0

Vfc � �, <sup>μ</sup><sup>1</sup> <sup>¼</sup> exp � Ux <sup>þ</sup> Uy

<sup>∂</sup><sup>r</sup> <sup>r</sup><sup>f</sup> cx � � <sup>¼</sup> Scx

ions of primary neutrals and cx atoms themselves, correspondingly, where

πV<sup>2</sup> th

cx is the total density of the source of cx atoms, with S<sup>0</sup>

exp � <sup>V</sup><sup>2</sup>

<sup>x</sup> <sup>þ</sup> <sup>V</sup><sup>2</sup> r V2

cxncx being the contributions due to charge exchange collisions with

the Heaviside function Θ ξð Þ¼ ≥ 0 1, Θ ξð Þ¼ < 0 0 and r<sup>0</sup> being the opening radius (see

�x∓ r ¼ const according to the continuity equation:

ð ∞

0

<sup>2</sup> <sup>Θ</sup> <sup>r</sup> <sup>þ</sup> <sup>x</sup> � <sup>r</sup><sup>0</sup> ð Þ<sup>ð</sup>

<sup>2</sup> <sup>Θ</sup> j j <sup>r</sup> � <sup>x</sup> � <sup>r</sup><sup>0</sup> ð Þ<sup>ð</sup>

<sup>μ</sup><sup>0</sup> <sup>¼</sup> exp � Ux � Uy � � � �

The velocity distribution function of cx atoms, f cx x; r; Vx; V<sup>r</sup>

V<sup>r</sup> r

Vx∂xf cx þ

nfcð Þ¼ x; r

þ Rfc

þ Rfc

2.3. Charge exchange cx atoms

cx <sup>þ</sup> <sup>S</sup><sup>1</sup>

cx <sup>¼</sup> nka

one can obtain

with

Figure 1b).

Here, Scx <sup>¼</sup> <sup>S</sup><sup>0</sup>

nbs <sup>þ</sup> nfc � �� and <sup>S</sup><sup>1</sup>

dis þ k m ion þ k m cx � �nm=<sup>4</sup>

and depends on r through nm. The particle densities nfc��, where the first subscript � corresponds to the sign of Vx and the second one – of Vr, change along the characteristics

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

<sup>p</sup> Vx∂lnfc�� <sup>¼</sup> SfcðÞ�<sup>l</sup> <sup>ν</sup>anfc��

where l is the length of the corresponding characteristics. The boundary conditions take into account that fc atoms are reflected with the probability Rfc from the wall, x ¼ 0, and are absent far from the wall, x ! ∞. For the total density of fc atoms, nfc ¼ nfcþþ þ nfcþ� þ nfc�þ þ nfc��,

> Sfcðy; j j r � x þ y Þ þ Sfcð Þ y; j j r þ x � y Vfc

> > Sfcð Þþ y; r þ x þ y Sfcð Þ y; j j r þ x � y Vfc

Sfcðy; j j r � x þ yÞ þ Sfcð Þ y; j j j j r � x � y Vfc

Vfc � �, Ux, <sup>y</sup> <sup>¼</sup>

μ0dyþ

μ1dyþ

http://dx.doi.org/10.5772/intechopen.76681

ð x, <sup>y</sup>

0

� �, is governed by the kinetic equation:

th ! � <sup>ν</sup>af cx, (3)

cx <sup>¼</sup> n km

cxnm þ k a cx �

μ1dy

νadz,

(2)

7

#### 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 Vbs,m and distinguish two groups of particles moving outward and toward the opening axis, with V<sup>r</sup> ¼ Vbs,m and V<sup>r</sup> ¼ �Vbs,m, correspondingly. The magnitudes of Vbs,m are predefined by the generation mechanism for the species in question. Within the plasma the point with certain coordinates x and r can be reached only by particles, coming from the wall points 0ð Þ ; r þ x for V<sup>r</sup> < 0 and 0ð Þ ; jr � xj for V<sup>r</sup> > 0, respectively. By taking into account the attenuation processes with neutrals in the plasma, one gets

$$n\_{bs,m}(\mathbf{x},\rho) = \frac{n\_{bs,m}(0,\rho+\mathbf{x}) + n\_{bs,m}(0,|\rho-\mathbf{x}|)}{2} \exp\left(-\int\_{V\_{bs,m}}^{\mathbf{x}} dy\right) \tag{1}$$

where <sup>ν</sup><sup>a</sup> <sup>¼</sup> n ka ion þ k a cx � � and <sup>ν</sup><sup>m</sup> <sup>¼</sup> n km dis <sup>þ</sup> <sup>k</sup><sup>m</sup> ion þ k m cx � � are the decay frequencies for atoms and molecules, correspondingly, with n being the plasma density assumed the same for electrons and ions, k a,m ion are the rate coefficients of the ionization by electrons, ka,m cx those for the charge exchange with ions, and k<sup>m</sup> dis 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<sup>r</sup> and of Vx. Henceforth, we distinguish four groups of these species with the velocity components Vx; V<sup>r</sup> � � ¼ �Vfc; �Vfc � �, where Vfc <sup>¼</sup> ffiffiffiffiffiffiffiffiffiffiffiffi Efc=m p and m is the atom mass. The source density of fc atoms of each group is as follows:

$$S\_{fc} = n \left( 2k\_{dis}^m + k\_{im}^m + k\_{cx}^m \right) n\_m / 4$$

and depends on r through nm. The particle densities nfc��, where the first subscript � corresponds to the sign of Vx and the second one – of Vr, change along the characteristics �x∓ r ¼ const according to the continuity equation:

$$\sqrt{2}V\_{\mathfrak{x}}\partial\_{l}\mathfrak{n}\_{\mathfrak{f}\mathfrak{c}\pm\pm} = S\_{\mathfrak{f}\mathfrak{c}}(l) - \mathfrak{v}\_{d}\mathfrak{n}\_{\mathfrak{f}\mathfrak{c}\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 Rfc from the wall, x ¼ 0, and are absent far from the wall, x ! ∞. For the total density of fc atoms, nfc ¼ nfcþþ þ nfcþ� þ nfc�þ þ nfc��, one can obtain

$$n\_{\hat{\boldsymbol{\mu}}}(\boldsymbol{\mu},\boldsymbol{\rho}) = \left[\frac{S\_{\hat{\boldsymbol{\mu}}}(\boldsymbol{y}, |\boldsymbol{\rho} - \boldsymbol{x} + \boldsymbol{y}|) + S\_{\hat{\boldsymbol{\mu}}}(\boldsymbol{y}, |\boldsymbol{\rho} + \boldsymbol{x} - \boldsymbol{y}|)}{V\_{\hat{\boldsymbol{\mu}}}} \mu\_0 d\boldsymbol{y} + \right] \tag{2}$$

$$\left[1 + \frac{R\_{\hat{\boldsymbol{\mu}}}}{2} \Theta(\boldsymbol{\rho} + \boldsymbol{x} - \boldsymbol{\rho}\_0)\right] \stackrel{\circ}{\underset{\boldsymbol{0}}{\operatorname{S}}} \frac{S\_{\hat{\boldsymbol{\mu}}}(\boldsymbol{y}, \boldsymbol{\rho} + \boldsymbol{x} + \boldsymbol{y}) + S\_{\hat{\boldsymbol{\mu}}}(\boldsymbol{y}, |\boldsymbol{\rho} + \boldsymbol{x} - \boldsymbol{y}|)}{V\_{\hat{\boldsymbol{\mu}}}} \mu\_1 d\boldsymbol{y} + \tag{2}$$

$$\left[1 + \frac{R\_{\hat{\boldsymbol{\mu}}}}{2} \Theta(|\boldsymbol{\rho} - \boldsymbol{x}| - \boldsymbol{\rho}\_0)\right] \stackrel{\circ}{\underset{\boldsymbol{0}}{\operatorname{S}}} \frac{S\_{\hat{\boldsymbol{\mu}}}(\boldsymbol{y}, |\boldsymbol{\rho} - \boldsymbol{x}| + \boldsymbol{y}) + S\_{\hat{\boldsymbol{\mu}}}(\boldsymbol{y}, ||\boldsymbol{\rho} - \boldsymbol{x}| - \boldsymbol{y}|)}{V\_{\hat{\boldsymbol{\mu}}}} \mu\_1 d\boldsymbol{y}$$

with

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 Efc of 3:5 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 Vx and V<sup>r</sup> of its velocity, where x and r are distances from the wall and opening axis (see Figure 1b).

6 Plasma Science and Technology - Basic Fundamentals and Modern Applications

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 Vbs,m and distinguish two groups of particles moving outward and toward the opening axis, with V<sup>r</sup> ¼ Vbs,m and V<sup>r</sup> ¼ �Vbs,m, correspondingly. The magnitudes of Vbs,m are predefined by the generation mechanism for the species in question. Within the plasma the point with certain coordinates x and r can be reached only by particles, coming from the wall points 0ð Þ ; r þ x for V<sup>r</sup> < 0 and 0ð Þ ; jr � xj for V<sup>r</sup> > 0, respectively. By taking into account the attenuation pro-

nbs,mð Þþ 0; r þ x nbs,mð Þ 0; j j r � x

dis <sup>þ</sup> <sup>k</sup><sup>m</sup>

ion are the rate coefficients of the ionization by electrons, ka,m

� �, where Vfc <sup>¼</sup> ffiffiffiffiffiffiffiffiffiffiffiffi

ion þ k m cx

molecules, correspondingly, with n being the plasma density assumed the same for electrons

dis that for the dissociation of molecules.

The fc atoms are born by the destruction of molecules and have positive and negative values both of V<sup>r</sup> and of Vx. Henceforth, we distinguish four groups of these species with the velocity

<sup>2</sup> exp �

ð x

0 @

� � are the decay frequencies for atoms and

νa,m Vbs,m dy

Efc=m p and m is the atom mass. The source

1

A (1)

cx those for the charge

0

2.1. Recycling molecules and atoms

cesses with neutrals in the plasma, one gets

nbs,mð Þ¼ x; r

� � ¼ �Vfc; �Vfc

density of fc atoms of each group is as follows:

ion þ k a cx � � and <sup>ν</sup><sup>m</sup> <sup>¼</sup> n km

where <sup>ν</sup><sup>a</sup> <sup>¼</sup> n ka

a,m

components Vx; V<sup>r</sup>

exchange with ions, and k<sup>m</sup>

2.2. Franck-Condon (fc) atoms

and ions, k

$$\mu\_0 = \exp\left(-\frac{\left|\mathcal{U}\_x - \mathcal{U}\_y\right|}{V\_{fc}}\right), \mu\_1 = \exp\left(-\frac{\mathcal{U}\_x + \mathcal{U}\_y}{V\_{fc}}\right), \mathcal{U}\_{x,y} = \int\_0^{x,y} \mathbf{v}\_x d\mathbf{z}\_x$$

the Heaviside function Θ ξð Þ¼ ≥ 0 1, Θ ξð Þ¼ < 0 0 and r<sup>0</sup> being the opening radius (see Figure 1b).

#### 2.3. Charge exchange cx atoms

The velocity distribution function of cx atoms, f cx x; r; Vx; V<sup>r</sup> � �, is governed by the kinetic equation:

$$V\_x \partial\_{\mathbf{x}} f\_{cx} + \frac{V\_{\rho}}{\rho} \partial\_{\rho} \left( \rho f\_{cx} \right) = \frac{S\_{cx}}{\pi V\_{th}^2} \exp\left( -\frac{V\_x^2 + V\_{\rho}^2}{V\_{th}^2} \right) - \mathbf{v}\_{\mathbf{a}} f\_{cx} \tag{3}$$

Here, Scx <sup>¼</sup> <sup>S</sup><sup>0</sup> cx <sup>þ</sup> <sup>S</sup><sup>1</sup> cx is the total density of the source of cx atoms, with S<sup>0</sup> cx <sup>¼</sup> n km cxnm þ k a cx � nbs <sup>þ</sup> nfc � �� and <sup>S</sup><sup>1</sup> cx <sup>¼</sup> nka cxncx being the contributions due to charge exchange collisions with ions of primary neutrals and cx atoms themselves, correspondingly, where

$$n\_{\infty} = \int\_{-\infty}^{\infty} \int\_{-\infty}^{\infty} dV\_{\text{x}} dV\_{\rho^2}$$

is the density of the latter; the velocity distribution of the source is assumed as the Maxwellian one with the ion temperature Ti, and Vth <sup>¼</sup> ffiffiffiffiffiffiffiffiffiffiffiffiffi 2Ti=m p is the ion thermal velocity.

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

To solve Eq. (3), <sup>V</sup><sup>r</sup> is discretized as �V<sup>l</sup> <sup>r</sup>, with V<sup>l</sup> <sup>r</sup> ¼ ΔV lð Þ � 1=2 , ΔV ¼ Vmax=lmax, and l ¼ 1, ⋯, lmax. The Vx distribution functions φ� <sup>l</sup> of particles with <sup>V</sup><sup>r</sup> in the range �V<sup>l</sup> <sup>r</sup> � ΔV=2 ≤ <sup>V</sup><sup>r</sup> <sup>≤</sup> � <sup>V</sup><sup>l</sup> <sup>r</sup> þ ΔV=2 are governed by the equations, following from the integration of Eq. (3) over these ranges:

$$V\_x \delta\_x \mathbf{q}\_l^{\pm} \pm \frac{V\_\rho^l}{\rho} \delta\_\rho (\rho \mathbf{q}\_l^{\pm}) = \frac{S\_{cx} \delta\_l}{2\sqrt{\pi} V\_{th}} \exp\left(-\frac{V\_x^2}{V\_{th}^2}\right) - \mathbf{v}\_d \mathbf{q}\_l^{\pm} \tag{4}$$

∂rð Þ rφ<sup>l</sup> ð Þy ≈ ∂rð Þ rφ<sup>l</sup> ð Þþ x ∂x∂rð Þ rφ<sup>l</sup> ð Þ� x ð Þ y � x , Uy ≈ Ux þ νað Þ� x ð Þ y � x ,

and (ii) by taking into account that the width of the region, occupied by cx atoms, is of several

γ<sup>l</sup> ≈ � ð Þ� Dl=r ∂rð Þ rφ<sup>l</sup> ,

Scxδ<sup>l</sup> ffiffiffi <sup>π</sup> <sup>p</sup> Vth

Equation (7) can be integrated as the first-order equation with respect to x, with the boundary conditions similar to those for γl, that is, φlð Þ¼ x ¼ 0 0 for Vx > 0 and φlð Þ! x ! ∞ 0 for

> ð ∞

Scxδ<sup>l</sup> ffiffiffi <sup>π</sup> <sup>p</sup> Vth

0

By neglecting the r derivative and by considering a single V<sup>r</sup> 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 � 1, δ<sup>l</sup> ¼ 1 and ncx ¼ η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

(ii) the characteristic dimension for the plasma parameter change is much larger than the typical MFPL for cx atoms Vth=νa. The same procedure can be performed in the case under

> ∂xð Þ ncxTi νam � �

the plasma with the velocity close to the ion thermal one, and (ii) neutral species are absent far

<sup>¼</sup> <sup>S</sup><sup>0</sup> cx � k a

th=ð Þ 2ν<sup>a</sup> . The boundary conditions to Eq. (9) imply as follows: (i) cx atoms leave

<sup>F</sup>αdu, with <sup>F</sup>αð Þ¼ <sup>u</sup> <sup>u</sup>�1exp �u<sup>2</sup> � <sup>α</sup>=<sup>u</sup> � � and <sup>α</sup> <sup>¼</sup> Uy � Ux

<sup>l</sup>¼<sup>1</sup> <sup>η</sup>l.

enough and the ionization rate is much smaller than that of the charge exchange, k

the consideration, providing the following 2D diffusion equation:

<sup>∂</sup><sup>r</sup> <sup>r</sup>∂rncx � � � <sup>∂</sup><sup>x</sup>

� Dth r

exp � <sup>V</sup><sup>2</sup> x V2

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

th ! � <sup>ν</sup>aφl: (7)

http://dx.doi.org/10.5772/intechopen.76681

<sup>r</sup> � ΔV=2 ≤ V<sup>r</sup> � � � � <sup>≤</sup> <sup>V</sup><sup>l</sup>

Iαdy � ηl, (8)

� � � <sup>r</sup> þ ΔV=2,

9

�=Vthð Þy . For the total

a ion ≪ k a cx, and

ionnncx, (9)

mean free path lengths (MFPL), that is, j j x � x<sup>0</sup> > j j Vx =ν<sup>a</sup> typically. As a result we get

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

with Dl <sup>¼</sup> <sup>V</sup><sup>l</sup>

<sup>η</sup><sup>l</sup> <sup>¼</sup> <sup>Ð</sup> ∞ �∞

where <sup>I</sup><sup>α</sup> <sup>¼</sup> <sup>Ð</sup>

with Dth <sup>¼</sup> <sup>V</sup><sup>2</sup>

∞ 0

2.3.2. Diffusion approximation

density of cx atoms, one has ncx <sup>¼</sup> <sup>P</sup><sup>l</sup>max

r � �<sup>2</sup>

=νa.

Vx∂xφ<sup>l</sup> � Dl

r ∂2 <sup>r</sup>ð Þ¼ rφ<sup>l</sup>

Vx < 0. Finally, for the density of cx atoms within the range V<sup>l</sup>

� Dl νar ∂2 <sup>r</sup> rη<sup>l</sup> � � <sup>¼</sup>

φldVx, one obtains the following integral-differential equation:

with

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

From the equation above, one gets the following ones for the variables φ<sup>l</sup> ¼ φ<sup>þ</sup> <sup>l</sup> þ φ� <sup>l</sup> and <sup>γ</sup><sup>l</sup> <sup>¼</sup> <sup>V</sup><sup>l</sup> <sup>r</sup> φ<sup>þ</sup> <sup>l</sup> � φ� l � �:

$$V\_x \delta\_x q \rho\_l + \frac{1}{\rho} \delta\_\rho (\rho \gamma\_l) = \frac{S\_{cx} \delta\_l}{\sqrt{\pi} V\_{th}} \exp\left(-\frac{V\_x^2}{V\_{th}^2}\right) - \mathbf{v}\_a q \rho\_l \tag{5}$$

$$
\partial\_{\mathbf{x}} \gamma\_l + \frac{\left(V\_{\rho}^{l}\right)^2}{V\_{\mathbf{x}} \rho} \partial\_{\rho} (\rho \mathbf{q} \mathbf{q}) = -\frac{\mathbf{v}\_a}{V\_{\mathbf{x}}} \gamma\_l \tag{6}
$$

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

$$\gamma\_l = -\frac{\left(V\_{\rho}^{l}\right)^2}{V\_x \rho} \int\_{\chi\_0}^{\chi} \partial\_{\rho}(\rho q\_l) \exp\left(\frac{\mathcal{U}\_y - \mathcal{U}\_x}{V\_x}\right) dy.$$

Here, x<sup>0</sup> is the position, where the boundary condition is imposed; since there is no influx of cx atoms from the wall, γlð Þ¼ x ¼ 0 0 for Vx > 0, neutral species are absent deep enough in the plasma and γlð Þ! x ! ∞ 0 for Vx < 0. The found expression for γ<sup>l</sup> can be made more convenient if (i) the x-variation of terms involved is approximated by linear Taylor series:

$$
\partial\_{\rho}(\rho q \rho\_l)(y) \approx \partial\_{\rho}(\rho q \rho\_l)(\mathbf{x}) + \partial\_{\mathbf{x}}\partial\_{\rho}(\rho q \rho\_l)(\mathbf{x}) \cdot (y - \mathbf{x}), \\
\mathcal{U}\_y \approx \mathcal{U}\_x + \mathbf{v}\_a(\mathbf{x}) \cdot (y - \mathbf{x}),
$$

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, j j x � x<sup>0</sup> > j j Vx =ν<sup>a</sup> typically. As a result we get

$$\gamma\_l \approx - \left( D\_l / \rho \right) \times \mathfrak{d}\_p (\rho q \varrho\_l)\_{\prime \prime}$$

with Dl <sup>¼</sup> <sup>V</sup><sup>l</sup> r � �<sup>2</sup> =νa.

ncx ¼

one with the ion temperature Ti, and Vth <sup>¼</sup> ffiffiffiffiffiffiffiffiffiffiffiffiffi

To solve Eq. (3), <sup>V</sup><sup>r</sup> is discretized as �V<sup>l</sup>

<sup>V</sup><sup>r</sup> <sup>≤</sup> � <sup>V</sup><sup>l</sup>

with

<sup>γ</sup><sup>l</sup> <sup>¼</sup> <sup>V</sup><sup>l</sup>

<sup>r</sup> φ<sup>þ</sup> <sup>l</sup> � φ� l � �:

of the flux density, one gets

these ranges:

lmax. The Vx distribution functions φ�

Vx∂xφ� <sup>l</sup> � Vl r r

<sup>δ</sup><sup>l</sup> <sup>¼</sup> erf � <sup>V</sup><sup>l</sup>

Vx∂xφ<sup>l</sup> þ

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

8 Plasma Science and Technology - Basic Fundamentals and Modern Applications

ð ∞ ð ∞

f cxdVxdVr;

<sup>r</sup> þ ΔV=2 are governed by the equations, following from the integration of Eq. (3) over

� erf � <sup>V</sup><sup>l</sup>

Scxδ<sup>l</sup> ffiffiffi <sup>π</sup> <sup>p</sup> Vth

The latter equation is straightforwardly integrated with respect to x, and for the r-component

∂rð Þ rφ<sup>l</sup> exp

Here, x<sup>0</sup> is the position, where the boundary condition is imposed; since there is no influx of cx atoms from the wall, γlð Þ¼ x ¼ 0 0 for Vx > 0, neutral species are absent deep enough in the plasma and γlð Þ! x ! ∞ 0 for Vx < 0. The found expression for γ<sup>l</sup> can be made more conve-

∂rð Þ¼� rφ<sup>l</sup>

( ) " # !

From the equation above, one gets the following ones for the variables φ<sup>l</sup> ¼ φ<sup>þ</sup>

∂rð Þ¼ rγ<sup>l</sup>

Vl r � �<sup>2</sup>

Vxr

ð x

x0

nient if (i) the x-variation of terms involved is approximated by linear Taylor series:

2Ti=m p is the ion thermal velocity.

<sup>l</sup> of particles with <sup>V</sup><sup>r</sup> in the range �V<sup>l</sup>

!

exp � <sup>V</sup><sup>2</sup> x V2 th

<sup>r</sup> � ΔV=2 Vth

exp � <sup>V</sup><sup>2</sup> x V2 th

!

νa Vx

Uy � Ux Vx � �

dy:

<sup>r</sup> ¼ ΔV lð Þ � 1=2 , ΔV ¼ Vmax=lmax, and l ¼ 1, ⋯,

� νaφ�

<sup>=</sup>erf <sup>V</sup>max Vth � � <sup>r</sup> � ΔV=2 ≤

<sup>l</sup> , (4)

:

� νaφl, (5)

γl, (6)

<sup>l</sup> þ φ�

<sup>l</sup> and

�∞

is the density of the latter; the velocity distribution of the source is assumed as the Maxwellian

�∞

<sup>r</sup>, with V<sup>l</sup>

∂<sup>r</sup> rφ� l � � <sup>¼</sup> Scxδ<sup>l</sup> 2 ffiffiffi <sup>π</sup> <sup>p</sup> Vth

<sup>r</sup> þ ΔV=2 Vth " # !

> 1 r

∂xγ<sup>l</sup> þ

Vl r � �<sup>2</sup>

Vxr

γ<sup>l</sup> ¼ �

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

$$V\_x \delta\_x q\_l - \frac{D\_l}{\rho} \delta\_\rho^2 (\rho q\_l) = \frac{S\_{cx} \delta\_l}{\sqrt{\pi} V\_{th}} \exp\left(-\frac{V\_x^2}{V\_{th}^2}\right) - \mathbf{v}\_a q\_l. \tag{7}$$

Equation (7) can be integrated as the first-order equation with respect to x, with the boundary conditions similar to those for γl, that is, φlð Þ¼ x ¼ 0 0 for Vx > 0 and φlð Þ! x ! ∞ 0 for Vx < 0. Finally, for the density of cx atoms within the range V<sup>l</sup> <sup>r</sup> � ΔV=2 ≤ V<sup>r</sup> � � � � <sup>≤</sup> <sup>V</sup><sup>l</sup> <sup>r</sup> þ ΔV=2, <sup>η</sup><sup>l</sup> <sup>¼</sup> <sup>Ð</sup> ∞ �∞ φldVx, one obtains the following integral-differential equation:

$$-\frac{D\_l}{\mathbf{v}\_a \rho} \partial\_\rho^2(\rho \eta\_l) = \int\_0^\bullet \frac{\mathbf{S}\_{c\mathbf{x}} \mathbf{\delta}\_l}{\sqrt{\pi} V\_{th}} I\_a d\mathbf{y} - \eta\_{l\nu} \tag{8}$$

where <sup>I</sup><sup>α</sup> <sup>¼</sup> <sup>Ð</sup> ∞ 0 <sup>F</sup>αdu, with <sup>F</sup>αð Þ¼ <sup>u</sup> <sup>u</sup>�1exp �u<sup>2</sup> � <sup>α</sup>=<sup>u</sup> � � and <sup>α</sup> <sup>¼</sup> Uy � Ux � � � �=Vthð Þy . For the total density of cx atoms, one has ncx <sup>¼</sup> <sup>P</sup><sup>l</sup>max <sup>l</sup>¼<sup>1</sup> <sup>η</sup>l.

#### 2.3.2. Diffusion approximation

By neglecting the r derivative and by considering a single V<sup>r</sup> 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 � 1, δ<sup>l</sup> ¼ 1 and ncx ¼ η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 a ion ≪ k a cx, and (ii) the characteristic dimension for the plasma parameter change is much larger than the typical MFPL for cx atoms Vth=νa. The same procedure can be performed in the case under the consideration, providing the following 2D diffusion equation:

$$-\frac{D\_{th}}{\rho}\partial\_{\rho}\left(\rho\partial\_{\rho}n\_{cx}\right)-\partial\_{x}\left[\frac{\partial\_{x}\left(n\_{cx}T\_{i}\right)}{\mathbf{v}\_{a}m}\right]=S^{0}\_{cx}-k^{\mu}\_{ion}n\mathbf{u}\_{cx}\tag{9}$$

with Dth <sup>¼</sup> <sup>V</sup><sup>2</sup> th=ð Þ 2ν<sup>a</sup> . 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, ncxð Þ¼ x ! ∞ 0. Additionally, ∂rη<sup>l</sup> ¼ 0 and ∂rncx ¼ 0 for r ¼ 0 and r ! ∞ for Eqs. (8) and (9), respectively.

g1ð Þþ u<sup>1</sup> u<sup>1</sup> þ u<sup>1</sup>

selected to reproduce Fαm;1;<sup>2</sup> ¼ F<sup>α</sup> um;1;<sup>2</sup>

<sup>I</sup><sup>α</sup> <sup>≈</sup> <sup>F</sup><sup>2</sup> α1 2F<sup>0</sup> α1 þ

þ F0 α1δ<sup>2</sup> <sup>1</sup> � F<sup>0</sup> α2δ<sup>2</sup> 2 <sup>10</sup> <sup>þ</sup> <sup>F</sup><sup>00</sup>

imate <sup>F</sup><sup>α</sup> by exp �u<sup>2</sup> � <sup>α</sup>=u<sup>2</sup>

to direct estimates.

2.3.4. Numerical procedure

one.

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi <sup>2</sup> <sup>þ</sup> <sup>6</sup>u<sup>2</sup> 1

� �, F<sup>0</sup>

αm δ3 <sup>1</sup> <sup>þ</sup> <sup>δ</sup><sup>3</sup> 2 120 þ

α ! 0, and it is not obvious that such a way is more time-saving than formula (10).

i. the source density Scxð Þ x; r is calculated, in the first approximation with ncx ¼ 0;

opening, the integral Iαð Þ x; y can be assessed once in advance. Then,

¼ 0, g1ð Þþ u<sup>2</sup> u<sup>2</sup> � u<sup>2</sup>

With properly selected initial approximations accurate enough, solutions of these equations

To assess the integrals Iα, we approximate function Fαð Þ u as a linear one at u ≤ u<sup>1</sup> and by polynomials of the fifth order at the intervals u<sup>1</sup> < u ≤ um and um < u ≤ u2. The coefficients are

ð Þ¼ <sup>u</sup>1,<sup>2</sup> 0. In addition, for <sup>u</sup> <sup>&</sup>gt; <sup>u</sup><sup>2</sup> the factor exp �u<sup>2</sup> � � leads to a very fast decay of <sup>F</sup><sup>α</sup> with increasing u. Therefore, for the aim to assess the corresponding contribution to Iα, we approx-

where δ1,<sup>2</sup> ¼ j j um � u1, <sup>2</sup> . Figure 2b shows I<sup>α</sup> versus α 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

Alternatively to the usage of formula (10), one can calculate the integral I<sup>α</sup> in advance and save the results in a table. Then, for any particular α, the integral I<sup>α</sup> is interpolated from the values in this table. However, some time is needed to find the proper α interval, and this time grows up with α ! 0 since I<sup>α</sup> ! ∞ as �lnα. Thus, the density of the data in the table has to be increased as

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

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

iii. the next approximation for ncxð Þ x; r and Scxð Þ x; r is computed, and the procedure is repeated till the relative error in, for example, ncxð Þ 0; 0 , becomes smaller than a desirable

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.

� �=u2. This results in the following expression:

Fα<sup>1</sup>δ<sup>1</sup> þ Fα<sup>m</sup>ð Þþ u<sup>2</sup> � u<sup>1</sup> Fα<sup>2</sup>δ<sup>2</sup>

<sup>α</sup>1,<sup>2</sup> ¼ F<sup>0</sup>

2 þ

ffiffiffi <sup>π</sup> <sup>p</sup> 2

<sup>α</sup>ð Þ u1, <sup>2</sup> , F<sup>00</sup>

1 � erfu<sup>2</sup> u2

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi <sup>2</sup> <sup>þ</sup> <sup>6</sup>u<sup>2</sup> 2

¼ 0:

http://dx.doi.org/10.5772/intechopen.76681

<sup>α</sup>ð Þ um , and F<sup>0</sup>

<sup>α</sup>ð Þ¼ um F<sup>00</sup>

α

11

(10)

q

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

<sup>α</sup><sup>m</sup> ¼ F<sup>00</sup>

exp � <sup>α</sup> u2 � �,

q

are found after three to five iterations with the Newton approach.

### 2.3.3. Assessment of the velocity space integral I<sup>α</sup>

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 ≤ y ≤ rw, and repeat this, in an iterative procedure, with respect to the whole profiles of ncxð Þx . In a 2D situation, the runs through the grid points in the r; V<sup>r</sup> � � phase subspace are involved additionally into calculations. These can be, however, efficiently parallelized. The integral I<sup>α</sup> 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αð Þ u is shown in Figure 2a for α ¼ 0:3, 1, and 3. One can distinguish four u intervals, where Fαð Þ u behaves principally differently: u ≤ u1, u<sup>1</sup> < u ≤ um, um < u ≤ u2, and u<sup>2</sup> < u. The characteristic points u1, <sup>2</sup> and um are defined by the conditions F<sup>0</sup> <sup>α</sup>ð Þ¼ um 0 and F<sup>00</sup> <sup>α</sup>ð Þ¼ u1,<sup>2</sup> 0 for the derivatives F<sup>0</sup> <sup>α</sup> ¼ �Fαg1=u<sup>2</sup> and <sup>F</sup><sup>00</sup> <sup>α</sup> <sup>¼</sup> <sup>F</sup>αg2=u<sup>4</sup> with

$$\mathcal{g}\_1(\mu) = 2\mu^3 + \mu - \mathfrak{a}, \; \mathcal{g}\_2(\mu) = \left(\mathfrak{g}\_1 + \mathfrak{u}\right)^2 - \left(2 + 6\mathfrak{u}^2\right)\mu^2.$$

For um one can use the Cardano formula for the real root of a cubic equation um <sup>¼</sup> ffiffiffiffiffiffiffiffiffiffi <sup>b</sup> <sup>þ</sup> <sup>a</sup> <sup>p</sup><sup>3</sup> � ffiffiffiffiffiffiffiffiffiffi <sup>b</sup> � <sup>a</sup> <sup>p</sup><sup>3</sup> with <sup>a</sup> <sup>¼</sup> <sup>α</sup>=4 and <sup>b</sup> <sup>¼</sup> ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi <sup>a</sup><sup>2</sup> <sup>þ</sup> <sup>1</sup>=<sup>216</sup> <sup>p</sup> . By searching for <sup>u</sup>1,2, we notice that for the roots of interest the equation g2ð Þ¼ u1,<sup>2</sup> 0 reduces to the following ones:

Figure 2. The functions <sup>F</sup>αð Þ¼ <sup>u</sup> <sup>u</sup>�1exp �u<sup>2</sup> � <sup>α</sup>=<sup>u</sup> � � for <sup>α</sup> <sup>¼</sup> 1 (solid curve), <sup>α</sup> <sup>¼</sup> <sup>0</sup>:3 (dashed curve), and <sup>α</sup> <sup>¼</sup> 3 (dasheddotted curve) (a); the integral <sup>I</sup><sup>α</sup> <sup>¼</sup> <sup>Ð</sup> ∞ 0 Fαð Þ u du estimated by the pass method (solid curve) and computed numerically (dashed curve) (b).

$$
\log\_1(u\_1) + u\_1 + u\_1\sqrt{2 + 6u\_1^2} = 0,\\
\mathcal{g}\_1(u\_2) + 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α, we approximate function Fαð Þ u as a linear one at u ≤ u<sup>1</sup> and by polynomials of the fifth order at the intervals u<sup>1</sup> < u ≤ um and um < u ≤ u2. The coefficients are selected to reproduce Fαm;1;<sup>2</sup> ¼ F<sup>α</sup> um;1;<sup>2</sup> � �, F<sup>0</sup> <sup>α</sup>1,<sup>2</sup> ¼ F<sup>0</sup> <sup>α</sup>ð Þ u1, <sup>2</sup> , F<sup>00</sup> <sup>α</sup><sup>m</sup> ¼ F<sup>00</sup> <sup>α</sup>ð Þ um , and F<sup>0</sup> <sup>α</sup>ð Þ¼ um F<sup>00</sup> α ð Þ¼ <sup>u</sup>1,<sup>2</sup> 0. In addition, for <sup>u</sup> <sup>&</sup>gt; <sup>u</sup><sup>2</sup> the factor exp �u<sup>2</sup> � � leads to a very fast decay of <sup>F</sup><sup>α</sup> with increasing u. Therefore, for the aim to assess the corresponding contribution to Iα, we approximate <sup>F</sup><sup>α</sup> by exp �u<sup>2</sup> � <sup>α</sup>=u<sup>2</sup> � �=u2. This results in the following expression:

$$\begin{split} I\_{\alpha} & \approx \frac{F\_{\alpha 1}^2}{2F\_{\alpha 1}'} + \frac{F\_{\alpha 1} \delta\_1 + F\_{am} (\mu\_2 - \mu\_1) + F\_{\alpha 2} \delta\_2}{2} + \\ & + \frac{F\_{\alpha 1}' \delta\_1^2 - F\_{\alpha 2}' \delta\_2^2}{10} + F\_{am}'' \frac{\delta\_1^3 + \delta\_2^3}{120} + \frac{\sqrt{\pi}}{2} \frac{1 - \operatorname{erf} \mu\_2}{\mu\_2} \exp\left(-\frac{\alpha}{\mu\_2}\right), \end{split} \tag{10}$$

where δ1,<sup>2</sup> ¼ j j um � u1, <sup>2</sup> . Figure 2b shows I<sup>α</sup> versus α 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<sup>α</sup> in advance and save the results in a table. Then, for any particular α, the integral I<sup>α</sup> is interpolated from the values in this table. However, some time is needed to find the proper α interval, and this time grows up with α ! 0 since I<sup>α</sup> ! ∞ as �lnα. Thus, the density of the data in the table has to be increased as α ! 0, and it is not obvious that such a way is more time-saving than formula (10).

#### 2.3.4. Numerical procedure

from the wall, that is, ncxð Þ¼ x ! ∞ 0. Additionally, ∂rη<sup>l</sup> ¼ 0 and ∂rncx ¼ 0 for r ¼ 0 and r ! ∞

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 ≤ y ≤ rw, and repeat this, in an iterative procedure, with respect to the whole profiles of ncxð Þx . In a 2D situation, the

calculations. These can be, however, efficiently parallelized. The integral I<sup>α</sup> 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αð Þ u is shown in Figure 2a for α ¼ 0:3, 1, and 3. One can distinguish four u intervals, where Fαð Þ u behaves principally differently: u ≤ u1, u<sup>1</sup> < u ≤ um, um < u ≤ u2, and u<sup>2</sup> < u. The characteristic points u1, <sup>2</sup> and um are defined by the

<sup>α</sup>ð Þ¼ u1,<sup>2</sup> 0 for the derivatives F<sup>0</sup>

<sup>g</sup>1ð Þ¼ <sup>u</sup> <sup>2</sup>u<sup>3</sup> <sup>þ</sup> <sup>u</sup> � <sup>α</sup>, g2ð Þ¼ <sup>u</sup> <sup>g</sup><sup>1</sup> <sup>þ</sup> <sup>u</sup> � �<sup>2</sup> � <sup>2</sup> <sup>þ</sup> <sup>6</sup>u<sup>2</sup> � �u<sup>2</sup>:

For um one can use the Cardano formula for the real root of a cubic equation um <sup>¼</sup> ffiffiffiffiffiffiffiffiffiffi

Figure 2. The functions <sup>F</sup>αð Þ¼ <sup>u</sup> <sup>u</sup>�1exp �u<sup>2</sup> � <sup>α</sup>=<sup>u</sup> � � for <sup>α</sup> <sup>¼</sup> 1 (solid curve), <sup>α</sup> <sup>¼</sup> <sup>0</sup>:3 (dashed curve), and <sup>α</sup> <sup>¼</sup> 3 (dashed-

Fαð Þ u du estimated by the pass method (solid curve) and computed numerically

� � phase subspace are involved additionally into

<sup>a</sup><sup>2</sup> <sup>þ</sup> <sup>1</sup>=<sup>216</sup> <sup>p</sup> . By searching for <sup>u</sup>1,2, we notice that for

<sup>α</sup> ¼ �Fαg1=u<sup>2</sup> and <sup>F</sup><sup>00</sup>

<sup>α</sup> <sup>¼</sup> <sup>F</sup>αg2=u<sup>4</sup>

for Eqs. (8) and (9), respectively.

conditions F<sup>0</sup>

<sup>b</sup> <sup>þ</sup> <sup>a</sup> <sup>p</sup><sup>3</sup> � ffiffiffiffiffiffiffiffiffiffi

dotted curve) (a); the integral <sup>I</sup><sup>α</sup> <sup>¼</sup> <sup>Ð</sup>

(dashed curve) (b).

with

2.3.3. Assessment of the velocity space integral I<sup>α</sup>

10 Plasma Science and Technology - Basic Fundamentals and Modern Applications

runs through the grid points in the r; V<sup>r</sup>

<sup>α</sup>ð Þ¼ um 0 and F<sup>00</sup>

<sup>b</sup> � <sup>a</sup> <sup>p</sup><sup>3</sup> with <sup>a</sup> <sup>¼</sup> <sup>α</sup>=4 and <sup>b</sup> <sup>¼</sup> ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

∞ 0

the roots of interest the equation g2ð Þ¼ u1,<sup>2</sup> 0 reduces to the following ones:

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αð Þ x; y can be assessed once in advance. Then,


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, ncxð Þ <sup>x</sup>; <sup>r</sup> , with a relative error below 10�<sup>6</sup> , 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 r directions, respectively, lmax =10 and Vmax ¼ 3Vth, 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αð Þ x; y 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 θ (see Figure 1a) and in Figure 4 x corresponds to the torus outboard, θ ¼ 0, where this distance is minimal. The profiles of the cx atom density found with these parameters and for the duct radius r<sup>0</sup> ¼ 0:05m, at the opening axis, r ¼ 0, and far from it, r ¼ 0:25m are shown in Figure 5. There is a noticeable difference between the profiles of ncx 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.

In the latter case, ncx 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 r0, at the distance h from its opening (see Figure 1b) is assessed. To do such assessment, consider an infinitesimally thin toroid within

Figure 5. The profiles of the cx atom densities far from the opening (solid curve) and at its axis (dashed curve) computed

Figure 4. The profiles of the plasma density (solid curve), the ion (dashed curve), and electron (dashed-dotted curve)

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

http://dx.doi.org/10.5772/intechopen.76681

13

temperatures versus the distance from the wall at the DEMO torus outboard.

in the diffusion approximation (a) and kinetically (b).

Noise-Free Rapid Approach to Solve Kinetic Equations for Hot Atoms in Fusion Plasmas http://dx.doi.org/10.5772/intechopen.76681 13

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.

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,

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

12 Plasma Science and Technology - Basic Fundamentals and Modern Applications

it is needed to solve the diffusion equation. Thus, for 500 and 30 grid points in the x and r directions, respectively, lmax =10 and Vmax ¼ 3Vth, 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αð Þ x; y has to be evaluated by each iteration. Even in this case kinetic computations for cx atoms are done faster by factor of 4.

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 θ (see Figure 1a) and in Figure 4 x corresponds to the torus outboard, θ ¼ 0, where this distance is minimal. The profiles of the cx atom density found with these parameters and for the duct radius r<sup>0</sup> ¼ 0:05m, at the opening axis, r ¼ 0, and far from it, r ¼ 0:25m are shown in Figure 5. There is a noticeable difference between the profiles of ncx 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

, kinetic calculations require much less CPU time than

ncxð Þ <sup>x</sup>; <sup>r</sup> , with a relative error below 10�<sup>6</sup>

3. Results of calculations

than for those computed kinetically.

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, ncx 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 r0, 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 dr, thickness dx, and the radius r, situated at the distance x from the wall. Cx atoms with the energies within the range E, E þ dE are generated in the toroid with the rate:

$$d\mathcal{R}\_{\rm cx}(\mathbf{x}, \rho, E) = \mathcal{S}\_{\rm cx}(\mathbf{x}, \rho) q\_{\rm i}(\mathbf{x}, E) 2\pi \rho d\rho d\mathbf{x} dE. \tag{11}$$

full-power year ð Þ nm=pfy , is hsp ¼ Γsp=nsp, with the particle density of molybdenum nsp. 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 hsp, found with ncx computed kinetically, does not unavoidably excides that obtained with ncx 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 Ti. This leads to noticeably smaller hsp (see Figure 7). In this case, oppositely to the situation at the outboard,

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

http://dx.doi.org/10.5772/intechopen.76681

15

Due to technical requirements, an acceptable mirror erosion rate in a future fusion reactor should not exceed 1nm=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 r<sup>0</sup> 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 r<sup>0</sup> ≪ 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

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).

hsp is larger if ncx is computed in the diffusion approximation.

where

$$q\_i(\mathbf{x}, E) = \frac{2\sqrt{E}}{\sqrt{\pi}T\_i^{1.5}(\mathbf{x})} \exp\left[-\frac{E}{T\_i(\mathbf{x})}\right]. \tag{12}$$

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 r ≤ rmax ¼ r0ð Þ 1 þ x=h 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

$$dj\_{\rm cx}(\mathbf{x}, \rho, E) = d\mathcal{R}\_{\rm cx} \frac{\mathbf{s} \cdot \exp(-\lambda/\mathbf{s})}{4\pi \left[ \left( \mathbf{x} + \mathbf{h} \right)^2 + \rho^2 \right]} \tag{13}$$

where <sup>s</sup> <sup>¼</sup> <sup>1</sup> <sup>þ</sup> <sup>r</sup><sup>2</sup>=ð Þ <sup>x</sup> <sup>þ</sup> <sup>h</sup> <sup>2</sup> h i�1=<sup>2</sup> is the cosine of the atom incidence angle with respect to the duct axis. The exponential factor in Eq. (13), with <sup>λ</sup>ð Þ¼ <sup>x</sup>; <sup>E</sup> Ux<sup>=</sup> ffiffiffiffiffiffiffiffiffiffiffiffi 2E=m p , 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<sup>1</sup> cx in the source density Scx.

The energy spectrum of cx atoms, hitting the mirror, is characterized by their flux density <sup>γ</sup> <sup>¼</sup> <sup>Ð</sup> djcx=dE in the energy range dE, where the integration is performed over r and x. By proceeding from <sup>r</sup> to <sup>s</sup>, according to the relation <sup>r</sup> <sup>¼</sup> ð Þ <sup>x</sup> <sup>þ</sup> <sup>h</sup> ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi <sup>1</sup>=s<sup>2</sup> � <sup>1</sup> <sup>p</sup> , one gets

$$\gamma(h,\rho\_0,E) = \int\_0^{r\_w} q\rho\_i(\mathbf{x},E)d\mathbf{x} \quad \int\_{h/\sqrt{h^2+\rho\_0^2}}^1 \frac{\mathbb{S}\_{cx}}{2} \exp\left(-\frac{\lambda}{s}\right)ds. \tag{14}$$

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

$$\Gamma\_{sp}(h,\rho\_0) = \int\_0^\infty dE \int \varphi\_i(\mathbf{x}, E) d\mathbf{x} \quad \int\_{h/\sqrt{h^2 + \rho\_0^2}}^1 \frac{S\_{c\mathbf{x}}}{2} \exp\left(-\frac{\lambda}{s}\right) Y\_{sp}(E, s) ds \tag{15}$$

Here, Ysp 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 nm per full-power year ð Þ nm=pfy , is hsp ¼ Γsp=nsp, with the particle density of molybdenum nsp. 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 hsp, found with ncx computed kinetically, does not unavoidably excides that obtained with ncx 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 Ti. This leads to noticeably smaller hsp (see Figure 7). In this case, oppositely to the situation at the outboard, hsp is larger if ncx is computed in the diffusion approximation.

the plasma, with a square cross section of the width dr, thickness dx, and the radius r, situated at the distance x from the wall. Cx atoms with the energies within the range E, E þ dE are generated

> E p ffiffiffi <sup>π</sup> <sup>p</sup> <sup>T</sup><sup>1</sup>:<sup>5</sup> <sup>i</sup> ð Þx

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 r ≤ rmax ¼ r0ð Þ 1 þ x=h can hit this mirror point. The contribution of the plasma volume in question to the density of the cx atom flux perpendicular

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 energy spectrum of cx atoms, hitting the mirror, is characterized by their flux density

djcx=dE in the energy range dE, where the integration is performed over r and x. By

ð 1

Scx

<sup>2</sup> exp � <sup>λ</sup>

<sup>h</sup><sup>=</sup> ffiffiffiffiffiffiffiffiffi <sup>h</sup>2þr<sup>2</sup> 0

p

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

ð 1

Scx

<sup>2</sup> exp � <sup>λ</sup>

s � �

<sup>h</sup><sup>=</sup> ffiffiffiffiffiffiffiffiffi h2 <sup>þ</sup>r<sup>2</sup> 0

p

Here, Ysp 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 nm per

<sup>φ</sup>ið Þ¼ <sup>x</sup>; <sup>E</sup> <sup>2</sup> ffiffiffi

djcxð Þ¼ x; r; E dRcx

duct axis. The exponential factor in Eq. (13), with <sup>λ</sup>ð Þ¼ <sup>x</sup>; <sup>E</sup> Ux<sup>=</sup> ffiffiffiffiffiffiffiffiffiffiffiffi

proceeding from <sup>r</sup> to <sup>s</sup>, according to the relation <sup>r</sup> <sup>¼</sup> ð Þ <sup>x</sup> <sup>þ</sup> <sup>h</sup> ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

r ðw

0

φið Þ x; E dx

φið Þ x; E dx

The latter process is taken into account by the contribution S<sup>1</sup>

γ h; r<sup>0</sup> ð Þ¼ ; E

ð ∞

0 dE r ðw

0

Γsp h; r<sup>0</sup> ð Þ¼

is the Maxwellian energy distribution function of ions.

14 Plasma Science and Technology - Basic Fundamentals and Modern Applications

dRcxð Þ¼ x; r; E Scxð Þ x; r φið Þ x; E 2πrdrdxdE: (11)

Tið Þx � �

: (12)

2E=m p , takes into account the

ds: (14)

Yspð Þ E;s ds (15)

cx in the source density Scx.

<sup>1</sup>=s<sup>2</sup> � <sup>1</sup> <sup>p</sup> , one gets

s � �

h i , (13)

exp � <sup>E</sup>

s � expð Þ �λ=s <sup>4</sup><sup>π</sup> ð Þ <sup>x</sup> <sup>þ</sup> <sup>h</sup> <sup>2</sup> <sup>þ</sup> <sup>r</sup><sup>2</sup>

is the cosine of the atom incidence angle with respect to the

in the toroid with the rate:

to the mirror surface is

where <sup>s</sup> <sup>¼</sup> <sup>1</sup> <sup>þ</sup> <sup>r</sup><sup>2</sup>=ð Þ <sup>x</sup> <sup>þ</sup> <sup>h</sup> <sup>2</sup> h i�1=<sup>2</sup>

where

<sup>γ</sup> <sup>¼</sup> <sup>Ð</sup>

follows:

Due to technical requirements, an acceptable mirror erosion rate in a future fusion reactor should not exceed 1nm=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 r<sup>0</sup> 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 r<sup>0</sup> ≪ 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

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 hsp scale is different than in Figure 6).

4. Conclusions

double dotted line).

Nomenclature

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 func-

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: ng ¼ 0 (solid line), 3 � 1017m�<sup>3</sup> (dashed line), 3 � <sup>10</sup>18m�<sup>3</sup> (dotted line), 1019m�<sup>3</sup> (dashed-dotted line), and 3 � 1019m�<sup>3</sup> (dashed-

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

http://dx.doi.org/10.5772/intechopen.76681

17

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

Dl,th Diffusivity of cx atoms in kinetic and diffusion

descriptions, correspondingly

tions, involving the Maxwellian velocity distribution of plasma ions.

the fraction of cx atoms which can hit and erode the mirror.

that it cannot contribute to mirror erosion. The reduction of the flux j <sup>E</sup> of incoming cx atoms with the energy E by collisions with the gas in the duct is governed by the equation:

$$d\mathbf{j}\_E/dl = -\sigma\_{el}(E)n\_\\$\mathbf{j}\_{E'} $$

where <sup>σ</sup>el <sup>≈</sup> <sup>3</sup>:<sup>8</sup> � <sup>10</sup>�<sup>19</sup>E�0:<sup>14</sup> <sup>m</sup><sup>2</sup>;eV � � 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\_{sp}(h,\rho\_0) = \int dE \int \rho\_i d\mathbf{x} \int \limits\_{h/\sqrt{h^2 + \rho\_0^2}}^{r\_w} \frac{S\_{ct}}{2} \exp\left(-\frac{\lambda + \sigma\_{cl} n\_\circ h}{s}\right) Y\_{sp}(E,s)ds \tag{16}$$

In Figure 8, the calculated dependences of the erosion rate hsp on h and r<sup>0</sup> are shown for several magnitudes of the density ng of the working deuterium gas in the mirror duct. One can see that the enhancement of ng above a level of 2 � 1019m�<sup>3</sup> should lead to the reduction of hsp below the target level of 1nm=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 eV. Such cold dense plasma cloud can affect the transfer of cx atoms in the plasma near the opening in the wall.

Noise-Free Rapid Approach to Solve Kinetic Equations for Hot Atoms in Fusion Plasmas http://dx.doi.org/10.5772/intechopen.76681 17

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: ng ¼ 0 (solid line), 3 � 1017m�<sup>3</sup> (dashed line), 3 � <sup>10</sup>18m�<sup>3</sup> (dotted line), 1019m�<sup>3</sup> (dashed-dotted line), and 3 � 1019m�<sup>3</sup> (dasheddouble dotted line).
