5.7 Boundary conditions

5.7.1 Zero order term: the Fokker-Planck equation

Momentum dynamics

Internal boundaries Momentum dynamics is described in spherical coordinates, because collisions is the dominant physical process for current drive and matrices are therefore expected to be well conditionned. Consequently, internal boundaries must be specified at p = 0 and = 1 which correspond to the condition

(5.459)

Here Neumann type boundary conditions are used, and therefore only gradients must be specified at the internal boundaries⁴. One must thus determine at p = 0 and = 1 for the flux grid only, since the distribution and flux grids are interlaced by definition. Therefore, starting from grids definition given in Sec. 5.2,

(5.460)

as Δξ_0,1∕2 = Δξ_0,-1∕2, ξ_0,1∕2 = et ξ_0,0 = -1.

Much in the same way,

(5.461)

since Δξ_0,nξ+1∕2 = Δξ_0,nξ-1∕2. Furthermore, since by definition,

(5.462)

one obtains finally

(5.463)

A similar technique may be used for the momentum grid p. Hence

(5.464)

using Δp_-1∕2 = Δp_1∕2. Since Δp_1∕2 = p₁ - p₀, and p_1∕2 = one obtains finally,

(5.465)

because p₀ = 0.

It is important to note that internal boundary conditions are only needed in the evaluation of the cross-derivative terms, since in the discrete form of the Fokker-Planck equation using two grids, as shown, in Sec. 5.4.1, they are automatically fullfiled for other terms. Indeed, at fixed ψ_l+1∕2,

(k+1)

(5.471)

with ξ_0,0² = ξ_{0,nξ 0}² = 1

External boundaries The other boundaries are inserted into the problem in violation of the true physical picture, mainly because the momentum domain extends off to infinity, while only a subspace is considered in computations. The upper limit of the domain is therefore chosen so that the interesting physics may be accurately described, i.e. p_max ≫ p_th for studying the Maxwellian distribution function, p_max ≫ γ_max for the RF waves and p_max ≫ p_Dreicer for the runaway electrons. These conditions must ensure the conservative nature of the problem here addressed, i.e. that plasma cannot enter or leave the domain of integration. This corresponds to the local condition

(5.472)

where is the normal to the external boundary. Because of the symmetry of the collision operator, the subspace considered for computations is a sphere of radius p_max, = . When the condition (5.472) is fullfiled, external boundary conditions have almost a negligible influence, and the numerical solution in the subspace is usually close to the theoretical one. It is usually the case for most RF problems except for the Lower Hybrid wave in very hot plasmas, since the range of resonant interaction between the wave and the electrons is well localized in momentum space, far enough from the boundary so that no electron can leave the domain of integration.

The Fokker-Planck equation reduces to an hyperbolic equation in the vicinity of p_max (high velocity limit), and therefore no boundary condition must be specified. In this condition, the standard upstream differencing applies. Since for collisions, D_pp^C∕v ~ 1∕v⁴ and F_p^C ~ 1∕v², one can consider that D_pp^C ≈ 0 at p = p_max. Here, the pitch-angle scattering term D_ξξ^Crequires no special handling, because it causes diffusion in a direction which is parallel to the boundary.

Nevertheless, since cross-derivative terms potentially involves points outside the computational domain, the following rule is used, even if the exact value of Δp_max has no specific importance,

(5.473)

As by definition p_np+1∕2 = = , it results

(5.474)

or in equivalent manner,

(5.475)

Runaway electron problem If an electric field is present, then in the real problem, some electrons will runaway and leave the domain of integration. In that case, the condition S_p ⋅ = 0 is no more fullfiled, leading to an effective loss of electrons. Since the Fokker-planck differential equation is of hyperbolic type, no specific modification must be applied to the limit of the integration domain. However, one must ensure that the total number of electrons is kept constant. Several techniques may be used to avoid a decay of the number of electrons at the runaway rate Γ_R, as defined in Sec. 3.6.

A possibility is to perform the following substitution

(5.476)

in the Fokker-Planck equation, so that f₀ is independent of t. In that case, the momentum matrix coefficients M_{p,l+1∕2,i+1∕2,j+1∕2} must be modified according to the prescription

(5.477)

where Γ_R,l+1∕2 is the runaway rate calculated at ψ_l+1∕2 and time step k as defined in Sec. 5.2. As a consequence, the Fokker-Planck equation becomes slightly non-linear, since Γ_R,l+1∕2 is an integral over the distribution function f₀ determined at time step k. However, since the runaway population is usually very small as compared to the bulk, the non-linearity remains fairly weak. Nevertheless, this approach has two major drawbacks: the Fokker-Planck equation has no more an intrinsic conservative form, and in addition matrix M_p must be recalculated at each time step, since it is time dependent. All the advantages of the numerical implicit time scheme for a fast and stable rate of convergence are therefore lost. For this reason, this method is not considered in the code.

Another technique is to inject cold electrons at p = 0, in order to compensate runaway losses and keep the electron density constant at ψ_l+1∕2. Such an approach has the main advantage to have a clear physical meaning. Indeed, the loss of electrons will generate locally an electric field which in turn will force a flux of particles to compensate this local depletion. By definition, the number of electrons leaving the integration domain is simply Γ_R. In principle, without particle losses, p₀²S_{p,l+1∕2,0,j+1∕2} = 0, though S_{p,l+1∕2,0,j+1∕2} is infinite at p = 0, since, by definition, the Maxwellian distribution is an exact eigenfunction of the Fokker-Planck operator. Compensation of hot electron losses by cold ones leads to the relation

(5.478)

By definition, S_{p,l+1∕2,0,j+1∕2} is assumed uniform in pitch-angle, since no specific direction can be physicaly priviledged. Therefore,

(5.479)

and

(5.480)

for all j values.

Expression (5.467) is then modified according to the relation

Since Γ_R,l+1∕2^(k+1) is a weak function of f₀, it is possible to replace Γ_R,l+1∕2^(k+1) by its explicit form Γ_R,l+1∕2^(k), so that one may avoid to recalculate matrix M_{p,l+1∕2,i+1∕2,j+1∕2} at each iteration, an extremely time consuming procedure. Consequently, with this scheme, the zero order Fokker-Planck equation becomes simply

The condition to maintain density constant at each radial location ψ_l+1∕2 leads therefore to an additional term, that adds to the operator which describes damping of Maxwellian electrons on suprathermal ones. With this approach, the Fokker-Planck equation is slightly non-linear, provided the fraction of runaway electrons remain small as compared to the bulk one.

Furthermore, the electron distribution function is not a Maxwellian in the vicinity of p = 0, since an ad-hoc source term of particle with no velocity is added. This approach is therefore questionable for its validity, owing to the basic assumptions that are used to derive the linearized Fokker-Planck equation around a Maxwellian bulk. Therefore, in order to avoid this singularity, an alternative approach may be to add a source term normalized to Γ_R,l+1∕2^(k) but, whose distribution function is an exact Maxwellian corresponding to the electron temperature at ψ_l+1∕2.

In that case, in steady-state regime, the flux divergence ∇_p ⋅ S_p is no more null at each plasma location, but is given by the relation

(5.484)

where T_e = T_e,l+1∕2, which corresponds to the existence of an external Maxwellian source term. By definition,

	2π ∑ j=0nξ0-1 ∑ i=0np-1 exp
	× pi+1∕22Δp i+1∕2Δξ0,j+1∕2 = ΓR,l+1∕2(k)	(5.485)

and the Fokker-Planck equation becomes simply

Though this technique is in principle the most consistent with the underlying physics, it has also still some numerical drawbacks, because of the forward time differencing. It is well known that such terms put strong limitations on the time step value Δt, in order to preserve numerical stability. A coarse estimate of its level can be deduced from the condition

(5.487)

Assuming f_{0,l+1∕2,i+1∕2,j+1∕2}^(k) ≃ f_{0M,l+1∕2,i+1∕2,j+1∕2}, it turns out that

(5.488)

where n_e,l+1∕2 = n_e is the local electron density. In normalized units, as defined in Sec.6.3, time step of the order of Δt ~ 10⁺³ - 10⁺⁴ are used to reach in few iterations the steady-state solution. Therefore, Γ_R,l+1∕2^(k)∕n_e,l+1∕2 must be lower than 10^-5 - 10^-6 in order to avoid onset of numerical instabilities, which corresponds usually to normalized Ohmic electric field less than 0.05.

Finally, the ultimate and most simple approach to deal with the loss of runaway electrons is to enforce the electron distribution function at each time step. Defining a numerical density n_e,l+1∕2^(k) at time step k ,

(5.489)

this procedure corresponds to the following replacement

(5.490)

This procedure is equivalent to reinject electrons in the plasma so that local density is steadily maintained. The main advantage of this method is that it preserves the implicit or backward numerical time advancing. Furthermore, it avoids flux calculations in order to determine the runaway rate at ech time step, a procedure which slows down the rate of convergence. Finally, the normalization procedure is general and may be applied for any process that leads to local electron losses, like radial transport, magnetic ripple losses. For this reason, this method is chosen, if needed, for the code. However, even if this technique is numericaly very powerful, and physicaly justified, it may hinder some artificial numerical electron losses resulting from an improper numerical flux balance in each cell. Consequently, it is crucial to benchmark the Fokker-Planck code without this option, in order to verify that the conservative scheme is well satisfied and that electron losses remain always small.

Lower hybrid wave and very high temperature plasmas For the Lower Hybrid current drive problem, difficulties similar to the runaway problem may take place in the core region of the plasma, when the electron temperature is very high. For this quasi-electrostatic wave, the resonant interaction with the plasma takes place along the magnetic field line, as discussed in Sec. 4.3.2, and therefore, a tail of very energetic electrons is pulled out from the bulk in the parallel direction p_∥ as for the Ohmic electric field. Unlike the Ohmic case, the domain of interaction for the Lower Hybrid wave is bounded in momentum space, and is given by the well known wave-particle resonance condition

(5.491)

where the relativistic factor γ = , β _th = v_th∕c and p² = p _∥² + p_⊥². On the axis p_⊥ = 0, the resonance domain is given by the relation

(5.492)

where

(5.493)

and

(5.494)

When p_⊥≠0, the lower bound of the resonance condition becomes

(5.495)

while the upper limit is

(5.496)

For very cold plasmas, i.e. when β_th²p_max² ≪ 1 since p_⊥≤ p_max,

(5.497)

and resonance domain boundaries are straight lines parallel to p_∥ = 0 axis. In that case, electron losses are naturally very weak, since at the intersection of the resonance and integration domains, the wave-induced flux is nearly tangent to the boundary of the integration domain. Since p_max ranges usualy between 20 - 30, β_th must be much less than 0.01, which corresponds to electron temperatures less than 0.06keV . Such a condition is never encountered in tokamak plasmas, except close to the very edge.

Because of relativistic corrections, lower and upper boundaries of the Lower Hybrid resonance domain are hyperbolic functions of p_⊥, which indicates that the wave induced particle flux always crosses the integration domain. The curvature of the resonance domain boundaries as a function of p_⊥ depends directly of β_th as shown in Fig. 5.3 and therefore, one may expect that particle flux induced by the wave may be nearly perpendicular to the integration domain boundary. This is the worst situation, for which important losses of electrons may take place if the quasilinear diffusion is very large between p_{∥ min}^LH and p_{∥ max}^LH, as for the runaway problem. In such a regime, a tail of very energetic electrons is pulled out from the bulk up to the integration domain, and the accuracy of the current driven by the wave as calculated by the code is highly questionable, when a significant fraction of this population leaves the integration domain.

This difficulty becomes rapidly important when β_th increases as a consequence of the forward peaking of the electron dynamic because of the relativistic corrections. In that case, the total kinetic energy comes into play, and therefore electrons with a high p_⊥ component will interact resonantly with the wave at larger p values because of their increasing mass.

It is difficult to introduce an unambiguous criterion which defines limits beyond which electron losses become important, since there is an interplay between the shape of the resonant domain, and also the absolute level of the quasilinear diffusion rate. As discussed in Sec. 4.3.2, its value is given by the relation

As shown in Fig. 5.3, the only possibility for particle losses to remain at an acceptable small level whatever the quasilinear diffusion level, is that the upper boundary of the Lower Hybrid resonance domain crosses the integration domain at p_⊥∕p_max > 1∕. This purely geometrical effect put very strong limitations on the lowest level of the parallel refractive index n_{∥ min}^LH that can be studied. Indeed, , p_{∥ max}^LH,0 is therefore given by the relation,

(5.499)

or

(5.500)

In term of parallel refractive wave index, the previous condition is equivalent to the simple relation

(5.501)

where the usual resonance relation

(5.502)

is used. In the limit β_th²p_max² ≫ 1, v_{∥ max}^LH becomes independent of p_max

(5.503)

which equivalent to

(5.504)

This result confirms that increasing p_max is useless in order to reduce particle losses, if the quasilinear diffusion coefficient D_b,n^LH(0) remains very large over all the quasilinear resonance domain, except if the parallel refractive index is small. In most calculations, p_max ≃ 20 - 30, and since in the core of the plasma β_th ranges between 0.14 and 0.17, β_th²p_max² ≃ 7.8 - 26, which fully justify this approximation.

As discussed for the runaway problem, particle losses are compensated by the normalization of the electron distribution function at each time step. However, when n_{∥ min}^LH exceeds , a warning is send by the code, in order to indicate that strong particle losses may take place if the quasilinear diffusion coefficient is large at the boundary of the integration domain.

Treatment of the trapped region Since characteristic times for collision, RF quasilinear diffusion and Ohmic electric field acceleration are all much larger than the bounce time in the weak collision “banana” regime, f₀ = f₀ when ≤ ξ_0T . Therefore, symmetrically place points around ξ₀ = 0 are equivalent in the trapped region, a physical property that is described by the term _T in the bounce averaging procedure in Sec. 2.2.1. This effect introduces a new internal boundary corresponding to the trapped-passing transition, and an external one at ξ₀ = 0, since by symmetry,

(5.505)

The Fokker-Planck equation is then solved in a reduced domain of integration, where the momentum phase space corresponding to -ξ_0T ≤ ξ₀ ≤ 0 is removed. Once determined in this reduced domain, the distribution function is recovered on the whole space, by enforcing the condition f₀ = f₀ in the interval ≤ ξ_0T .

A specific treatment of the new internal boundary must be considered, so that flux continuity is correctly described between co- and counter-passing regions, despite the reduction of the phase space domain. Let define grid points

(5.506)

where j_T0 = ∕2 = ∕2, as indicated in Fig.5.4. Here j_T,l+1∕2^- is the index that corresponds to the trapped/passing boundary ξ₀ = -ξ_0T , while j_T,l+1∕2⁺ corresponds to the symmetric boundary with respect to the axis ξ₀ = 0, i.e. ξ₀ = ξ_0T .

By definition, grid points 1 and 3 are in the trapped region, just placed below and above co- and counter-passing boundaries respectively at opposite sides of the axis ξ₀ = 0 and with the same distance to this axis, while grid points 2 and 4 have a similar arrangement but are placed in the co- and counter-passing regions respectively right after the trapped-passing boundary. The grid point 5 is placed on the right side of the axis ξ₀ = 0, in the trapped sub-domain 0 ≤ ξ₀ ≤ ξ_0T .

Before performing the reduction of the domain of integration, the symmetry inside the trapped region must be enforced according to the relation

(5.507)

where j_T0 = ∕2, with j_T,l+1∕2^- < j^- < j_T0 or an in equivalent manner j_T0 < j⁺ < j_T,l+1∕2⁺. This procedure allows to correctly account for some processes that are not symmetric in p_∥, like wave-particle interaction.

The reduction of the domain of integration in the momentum phase space leads to introduce a new matrix, _p whose coefficients are defined as follows

(5.508)

with j^′ = j, when 0 ≤ j < j_T,l+1∕2^-, and j^′ = j - n_ξ0T,l+1∕2, when j_T0 ≤ j ≤ n_ξ0 - 1. Here, n_ξ0T,l+1∕2 = ∕2 corresponds to the number of grid points removed in the pitch-angle direction, and with this definition, it can be easily cross-checked that j_T0 - n_ξ0T,l+1∕2 = j_T,l+1∕2^-.

Furthermore, since grid point 3 is no more considered in the calculations while grid points 3 and 1 are equivalent, flux relations between neighboring grid points 4 and 3 must then be replaced by new relations wich link grid points 4 and 1. Much in the same way, flux relations between grid points 2 and 3 no more exist, and therefore, the flux relation between grid points 2 and 1 must be modified accordingly. This procedure leads to add six new set of coefficients to the matrix _p, namely

(5.509)

and

(5.510)

while matrix coefficients

(5.511)

and

(5.512)

where j_T,l+1∕2^′- = j_T,l+1∕2^- and j_T,l+1∕2^′+ = j_T,l+1∕2⁺ - n_ξ0T,l+1∕2 must be modified accordingly.

As shown in a graphic form in Figs. 5.5 and 5.6,

(5.513)

since the distribution function at grid 4 is no more connected with the one at grid point 3 or 5, but to the grid point 1 only. Furthermore, since grid points 1 and 3 are equivalent, their clockwise flux links with grid point 2 are only half in order to avoid counting them twice. Therefore,

Moreover, since counterclockwise flux that links grid points 3 and 4 is equivalent by symmetry with clockwise flux that links grid points 1 and 2 , one obtains,

Finally, for the external boundary at ξ₀ = 0, all matrix coefficients at grid point j_T0 + 1∕2 are forced to zero, except two,

Magnetic ripple losses Super-trapped electrons with p_∥∕p_⊥≪ 1 are very sensitive to all the details of the magnetic configuration, and may be lossed in the local magnetic wells because of the finite number of toroidal magnetic field coils. This effect requires in principle a full 3 -D spatial calculation, that is beyond the present theoretical frame here described, as for the stellarator magnetic configuration. Nevertheless, some interesting results may be obtained, by considering that super-trapped electrons whose banana tips fall in the magnetic well drift vertically in an irreversible manner and leave therefore the plasma. This process takes place provided their kinetic energy is large enough so that the detrapping probability by collisions remain negligible.

Describing this effect is consequently equivalent to introduce new external boundaries in the trapped region, in order to describe in an ad-hoc manner this physical effect. This approach is justified because the drifting time across the plasma is short enough for the locally trapped electrons in the magnetic ripple to neglect their contributions to the overall momentum dynamics.

For this purpose, an Krook term is introduced in the Fokker-Planck equation

(5.518)

where ν_{MR,l+1∕2,i+1∕2,j^′+1∕2} is an effective loss frequency, whose value is a function of the super-trapped domain as shown in Fig. 5.7. This loss frequency varies therefore with the radial location ψ_l+1∕2, but also with momentum p_i+1∕2 and pitch-angle ξ_{0,j^′+1∕2}. Outside this domain, ν_{MR,l+1∕2,i+1∕2,j^′+1∕2} is set to a zero so that losses are always negligible on the time scale for the steady-state distribution function to build-up. Conversely, inside the domain characterized by the pitch-angle boundary 0 ≤ ξ_{0,j^′+1∕2} ≤ ξ_{0,jMR,l+1∕2^′++1∕2} and p_i+1∕2 ≥ p_MR,l+1∕2, where j_MR,l+1∕2^′+ is the index below which particle are lossed in the magnetic ripple, and p_MR,l+1∕2 is the momentum detrapping threshold by collisions, ν_{MR,l+1∕2,i+1∕2,j^′+1∕2} is usually taken constant and much larger than τ_t^-1 which is the smallest time scale here considered, namely the transit time. The matrix is then simply modified according to the following rule

(5.519)

With this method, it is not necessary to introduce additional specific conditions and moreover the implicit time scheme is fully preserved. By definition, the Fokker-Planck equation is no more conservative, but particle losses are compensated by the normalization technique at each time step, as for the runaway problem. This is an acceptable technique, provided level of losses due to magnetic ripple is small as compared to the electron bulk density.

Spatial dynamics

Internal boundary For the dynamics in configuration space, since the problem has only one dimension, the internal boundary must be only specified at ψ = ψ₀ = 0, which corresponds to the condition

(5.520)

Using the same approach as for the momentum dynamics, Δψ must be defined at this location. Hence

(5.521)

using Δψ_-1∕2 = Δψ_1∕2. Since Δψ_1∕2 = ψ₁ - ψ₀, and ψ_1∕2 = one obtains finally,

(5.522)

because ψ₀ = 0.

It is important to note that internal boundary conditions are also only needed for the evaluation of cross-derivative terms, since in the discrete form of the Fokker-Planck equation using two grids, as shown, in Sec. 5.4.2, it is automatically fullfiled for other terms. Indeed, at fixed p_i+1∕2,ξ_0,j+1∕2,

(k+1)

External boundary As discussed for the momentum dynamics, a primary condition of the code is that electron density profile is preserved. As for the particle losses in momentum space, this can be performed by injecting an equivalent number of particles at p = 0, in order to compensate the particle flux leaving the integration domain because of diffusion across magnetic flux surfaces. It corresponds to a reaction of the bulk, because of the occurence of a small electric field, generated by the electron transport so that the overall electron density profile is kept constant. By this method, a steady state regime may be achieved, and it is not necessary to specify which type of mechanism is at play, diffusion or convection.

Let define the local variation of the electron density Δn_e,l+1∕2 at time step k because of spatial transport by the simple relation

(5.525)

It can be expressed in term of a loss rate

by introducing the integration time step Δt and the flux balance gives on the discrete mesh is

(5.526)

assuming like for losses in momentum space that the flux has no pitch-angle dependence at p₀. The zero order Fokker-Planck equation becomes simply

Finally, for cross-derivative terms, one have to determine Δψ_nψ at the boundary of the integration domain. It is readily given by relation,

(5.529)

as done for p, since ψ_nψ = ψ_a.

Treatment of the trapped region The radial transport equation is greatly complicated by the presence of trapped electrons, since the diffusion across magnetic field lines may contribute to trap or detrap particles who are in the close vicinity of the trapped/passing boundary, because of the radial dependence of this boundary, as shown in Fig. 5.8. By definition, this process takes place intrinsically at fixed magnetic moment, a fundamental assumption of the flux conservative form of the bounce averaged Fokker-Planck equation, as shown in Sec.3.5.1.

The matrix build-up is therefore quite complex, since one has to take into account that counter-passing electrons may be trapped when they diffuse towards lower magnetic regions of the plasma, while conversely, trapped electrons may be detrapped by crossing magnetic flux surfaces in the high magnetic field side direction.

Let define three neighboring grid points at radial locations ψ_l-1∕2, ψ_l+1∕2 and ψ_l+3∕2 and corresponding trapped/passing boundaries ξ_{0T,jT,l-1∕2^±,} ξ_{0T,jT,l+1∕2^±,} and ξ_{0T,jT,l+3∕2^±}, since non-uniform pitch-angle grid mesh is designed so that each boundary is exactly placed on the flux grid. By definition, since

(5.530)

the following relation holds

(5.531)

Here j_T0 is the index of the grid point corresponding to the axis ξ₀ = 0 in the trapped region at ψ_l+1∕2. This defines naturaly eight different regions in momentum space at ψ_l+1∕2, where spatial fluxes must be considered. Since matrix is reduced in momentum space, because the region j_l+1∕2^- < j_l-1∕2^- < j_T0,l+1∕2 is not considered in the calculations, one has therefore to consider effectively only six regions, namely corresponding to

(5.532)

and

(5.533)

or using relations j_T,l+1∕2^′- = j_T,l+1∕2^- and j_T,l+1∕2^′+ = j_T,l+1∕2⁺ - n_ξ0T,l+1∕2

(5.534)

As for the momentum dynamics, the reduction of the domain of integration in the momentum phase space leads to introduce a new matrix, _ψ. Obviously, the main diagonal is fully preserved, and therefore

(5.535)

provided j^′ < j_T,l+1∕2^′-, while

(5.536)

when for j^′ > j_T,l+1∕2^′+. Only off-diagonal matrix coefficients have to be modified, because of the presence of trapped and circulating electrons. Indeed, these coefficients reflect the particle flux transfer between different location.

If j^′ < j_T,l+3∕2^′-, matrix coefficients are given by

(5.537)

and for j_T,l+3∕2^′- < j^′ < j_T,l+1∕2^′-, the simple relation still holds for one matrix coefficient only

(5.538)

since the trapped region is narrower at ψ_l-1∕2 as compared to ψ_l+1∕2. But since counter-circulating electrons become trapped by diffusion from ψ_l+1∕2 → ψ_l+3∕2 and electron dynamics in the domain j_T,l+3∕2^- < j < j_T0 is removed and replaced by its counterpart in the positive pitch-angle region j_T0 < j < j_T,l+3∕2⁺, it turns out that

(5.539)

and the flux link must be replaced by

(5.540)

taking into account of the mirror indexing around the axis ξ₀ = 0, and the fact that matrix blocks corresponding to different radii have different sizes.

In a general way, when j^′ > j_T,l+1∕2^′-, the following rule applies

(5.541)

all corresponding matrix coefficients being replaced by new ones according to the relations

(5.542)

wether electrons remains trapped or co-passing, or become trapped by the outward radial transport or detrapped by the inward cross-field diffusion. Since n_ξ0T,l+1∕2 - n_ξ0T,l-1∕2 > 0 and n_ξ0T,l+1∕2 - n_ξ0T,l+3∕2 < 0,

(5.543)

which indicates that diagonals shrink when trapped/passing boundary is crossed, i.e. j^′ > j_T,l+1∕2^′-. Therefore, matrix _ψ contains much more diagonals than M_ψ which has only three ones.

At ψ = 0, there is no need to specify the spatial flux S_ψ, with the two grids technique. Since by construction R_min,0B_{P min,0}S_{ψ0,i+1∕2,j+1∕2}^(k+1) = 0, no particle accumulation occurs naturaly at ψ , and at this singular point f_{0,1∕2,i+1∕2,j+1∕2}^(k+1) = f_{0,-1∕2,i+1∕2,j+1∕2}^(k+1) is well fullfiled.

At the plasma edge ψ_a, i.e. at the boundary of the integration domain, there is also no need to specify any condition on the flux of fast particles are leaving the plasma definitively, as for the runaway problem. Consequently, even without any external perturbation, the fact that radial transport of fast transport leads to strong departure from the Maxwellian momentum dependence at all radius.

As for the dynamics in momentum space, Δψ must be specified at the boundaries of the domain of integration. One finds that

(5.544)

since ψ₀ = 0, and

(5.545)

noting that ψ_nψ = ψ_a.

5.7.2 Up to first order term: the Drift Kinetic equation

Since the first order bounce-averaged drift kinetic equation may be expressed in the same conservative form as for zero order Fokker-Planck bounce-averaged equation, most of the boundary conditions apply in a natural way without additional specifications, in particular at p = 0 and = 1, as a consequence of the two grids technique, one for the fluxes and the other for the distribution function.

In matrix form, the general expression of the first order bounce averaged drift kinetic equation

(5.546)

becomes

Here, is the iteration index number, while C and are the electron-electron self-collision operators that make the equation non-linear.

Internal boundaries for function g

In the weak collision “banana” regime, it is shown in Sec.3.4 that g = 0 for trapped electrons. Therefore, the drift kinetic equation must be solved in a reduced domain of integration, as for f₀, but with different and more simple internal boundary conditions. Indeed, the only constraint is to enforce the condition g = 0 for ≤ ξ_0T .

Let define grid points

(5.550)

as indicated in Fig.5.9. Here j_T,l+1∕2^- is the index that corresponds to the trapped/passing boundary ξ₀ = -ξ_0T , while j_T,l+1∕2⁺ corresponds to the symmetric boundary with respect to the axis ξ₀ = 0, i.e. ξ₀ = ξ_0T .

Since g is prescribed in the trapped domain, all the region may be removed from the calculations, and the reduction of the domain of integration in the momentum phase space leads to introduce a new matrix _p whose coefficients are defined as follows

(5.551)

with j^′ = j, when 0 ≤ j < j_T,l+1∕2^-, and j^′ = j - m_ξ0T,l+1∕2, when j_T,l+1∕2⁺ < j ≤ n_ξ0. Here, m_ξ0T,l+1∕2 = j_T,l+1∕2⁺ - j_T,l+1∕2^- corresponds to the number of grid points removed in the pitch-angle direction, and with this definition, it is easy to cross-check that j_T,l+1∕2^′+ -j_T,l+1∕2^′- = j_T,l+1∕2⁺ -j_T,l+1∕2^--m_ξ0T,l+1∕2 = 0. Since grid points 2 and 4 are not connected, the following condition applies

(5.552)

and conversely

(5.553)

where j_T,l+1∕2^′- = j_T,l+1∕2^′+ = j_T,l+1∕2^′. As internal boundaries match already existing diagonals of the matrix M_{p,l+1∕2,i+1∕2,j+1∕2}, _{p,l+1∕2,i+1∕2,j^′+1∕2} remains a simple nine diagonal matrix, and is therefore less complex than the 15 diagonals matrix _{p,l+1∕2,i+1∕2,j^′+1∕2} corresponding to the Fokker-Planck equation.

Internal boundaries for

Since no inversion is required for calculating the vector _{l+1∕2,i+1∕2,j+1∕2}, it is consequently evaluated at all grid points, including those corresponding to trapped electron domain. Consequently, the reduction along the pitch-angle direction that results from the condition g = 0 is only performed after the calculation of _{l+1∕2,i+1∕2,j+1∕2} over the full ξ₀ grid. Once done, the prescription on the index number j^′ that holds for the matrix conversion

(5.554)

may be then applied in a similar way for _{l+1∕2,i+1∕2,j+1∕2}, which leads to define a new vector _{p,l+1∕2,i+1∕2,j^′+1∕2} whose size is reduced

(5.555)

and finally the drift kinetic equation on the reduced grid is simply given by the relation

(5.556)

where all the trapped electron domain is removed.

External boundaries

The problem of external boundaries discussed in Sec. 5.7.1 for the zero-order Fokker-Planck equation may be applied in a straightforward manner to the first-order drift kinetic ones, since it may be expressed in a similar conservative form. Consequently, if the zero-order Fokker-Planck equation conserves the total number of particles without specific conditions, the solution g may be found without additional constraints at the boundaries of the integration domain. Conversely, if Dirichlet conditions are considered for the Fokker-Planck equation, like enforcing the Maxwellian solution at p = 0, similar and consistent conditions must be also applied for the drift kinetic equation.