Modelling of the Ignitor scrape-off layer including neutrals

Ignitor is a tokamak project aimed at achieving ignition. In the reference scenario, plasma-surface interactions are controlled by a Mo first-wall/limiter, which constitutes a simple engineering solution but, at the same time, a special challenge for edge plasma modelling. Here the ASPOEL plasma fluid code, already applied to Ignitor in the recent past, is coupled with the neutral Monte Carlo code EIRENE. We study the effects of the neutrals on the plasma density and temperature profiles in the Ignitor scrape-off layer, and compute the particle and heat loads onto the Ignitor first-wall limiter.


Introduction
The main goal of the Ignitor experiment [1] is that of establishing the "reactor physics" (i.e. the physics of power producing reactors) in regimes close to ignition, where the "thermonuclear instability" can set in with all its associated non linear effects. The machine (see Table 1) is characterized by an optimal combination of high magnetic field (B T < 13 T), compact dimensions (R 0 ≅ 1.32 m), relatively low aspect ratio (R 0 /a ≅ 2.8), and considerable elongation and triangularity of the plasma cross section. The optimal central density of the fusing nuclei to achieve ignition is estimated to be about 10 21 m -3 . The corresponding line-averaged density is well below the known density limit (related to the average plasma current density) for magnetically confined plasmas, as the plasma current I p allowed by the machine design can reach 11 MA. Ignition can be achieved by ohmic heating alone, shortly after the end of the plasma current rise. The peak temperature at ignition is expected to be about T e0 ≅ T i0 ≅11 keV for an energy confinement time τ E ≅ 0.6 sec; the relatively low beta poloidal, i.e. the ratio of the plasma energy density (pressure) to the poloidal magnetic field energy density, is consistent with favourable conditions for macroscopic plasma stability. An ICRH system operating between 80 and 120 MHz is also available, to accelerate the attainment of ignition by applying modest amounts of RF power during the current rise, or to facilitate access to the H-mode regime in double X-point configurations at lower plasma currents.
According to an independent assessment [2], "ignition in Ignitor would be reached under nonsteady state conditions, as in fact ignition is intrinsically a time evolving event. The ignited regime would last a few confinement times (2-5× E , E 0.5 s) and many alpha particle slowing down times (50-100 sd ), which is adequate from the physics point of view". Here we shall refer to the operating scenario described in figure 1, which was studied in [3].    From the point of view of plasma-surface interactions (PSI), Ignitor constitutes a significant and somewhat special challenge to modelers, as it adopts a first wall/limiter (FWL), which is rather different from the now standard divertor configuration of many present and future tokamaks. While 0D-1D models of PSI in Ignitor have been developed long time ago with emphasis on high-Z impurity screening in high-density limiter tokamaks [4], and also validated against FTU data [5], the issues related to 2D modeling accounting for the singularity introduced by the tangency line between FWL and last closed flux surface (LCFS) at the inboard equatorial plane of the machine have been only addressed in recent work: a new code, ASPOEL, was developed for that purpose by our group at Politecnico di Torino [6], [7], and also applied to and validated against divertor far scrape-off layer (SOL) data [8].
Here we extend the work presented in [7] by coupling ASPOEL with the Monte-Carlo code for neutrals EIRENE [9]. The coupling algorithm is described and the main effects introduced by the presence of the neutrals on both plasma profiles and FWL loads will be highlighted. The main purpose of the paper is of computational nature, namely to demonstrate how the coupled tools can work in a realistic FWL geometry, while no claim is made as to the direct physics relevance for Ignitor of the computational test conditions considered here.

Approach to edge modeling in FWL configurations
The ASPOEL code adopts a control-volume finite-element (CVFE) hybrid approach [10]; it uses dual CV (polygonal) and FE (triangular) grids, where the FE give the natural interpolation scheme for the approximation of the fluxes on the CV faces. In view of the strong anisotropy of the transport coefficients introduced by the magnetic field, it is well known that it is essential to align the triangles with the local field direction [11]; therefore, the mesh generation in ASPOEL is done such as to guarantee perfect alignment of one side of each triangle. Approaches using less stringent alignment mesh generation criteria [12] are also being pursued by other authors, but will not be discussed here.
The basic strategy of mesh generation adopted here starts from an equilibrium computed with the MAXFEA code [13], see figure 2; we then use the CARRE code [14] to build a first grid which ignores the FWL geometry; this grid is finally adapted to the FWL by adding magnetic (poloidal) surfaces at the intersection between the FWL and the given radial surfaces (other approaches to the same problem have been proposed more recently [15]). ASPOEL solves the standard set of two-fluid equations for electrons and a single fuel ion species, consisting of: continuity, quasi-neutrality, diffusive Ansatz for the (anomalous) fluxes in radial direction, parallel (to the magnetic field) momentum balance, electron and ion energy balance. For the sake of demonstration purposes we assume a pure plasma with deuterium as the only ion species: the model is at this time unable to treat the two fuel ion species (D, T) simultaneously, or any impurity species. Steady state conditions in the edge plasma are assumed, considering that the time scale of relevance in the SOL (~ ms) is much faster than that of the evolution of the main plasma (~ 0.1 -1 s).
Particle, momentum and energy sources are provided by EIRENE. ASPOEL computes the (background) plasma conditions and the plasma (D + ) flux onto it; EIRENE follows the neutral particles recycling from the FWL (e.g. D, D 2 ) as a result of the plasma bombardment, until they suffer a collision. The main atomic physics processes included in the model for the case at hand are: electron impact (including ionization), charge-exchange, elastic collisions and volume recombination, for both atomic and molecular neutral hydrogenic particles. About 30000 neutral particles were launched from the FWL and followed to perform the statistics + a comparable number arising from volume recombination. According to now standard practice in the B2-EIRENE coupling [16], EIRENE was called every ASPOEL iteration.

Problem definition and boundary conditions
We consider the region between the FWL and the magnetic surface located 2 mm inside the LCFS at the inboard equatorial plane (as customary, see e.g. [16], not only the SOL but also a portion of the main plasma is included in the computational domain, in order to relax the effect of conditions imposed at the inner boundary -indeed, even a larger portion should possibly be profitable in this case, since the LCFS is always very close to the FWL in Ignitor). Up-down symmetry is assumed, so that only the portion of the domain above the equatorial plane will be considered.
For the purpose of the computational test to be presented in this paper, boundary conditions (BC) have been imposed as follows: at the inner boundary, Dirichlet conditions n = 10 20 m -3 , T e = T i = 15 eV V || = 0 m/s; at the FWL, the standard Bohm-Chodura criterion V || = C S , the Teilhaber-Birdsall criterion [17] V r = 0.01C S , the generalized conditions on energy transmission through the Debye sheath Γ En e = 5.5 * Γ n T e and Γ En i = 3.5 * Γ n T i . Concerning these BC it should be noted that somewhat higher n and T at the inner boundary could be more physically relevant for Ignitor; also, as the BC at the FWL are rather uncertain, in view of the continuous transition between tangency and finite incidence of the magnetic field line, it should possibly be useful to perform a sensitivity study of their effect on the solution.
In order to guarantee a minimum of consistency with the main plasma conditions, it is important to consider the main plasma power balance (first principle of thermodynamics) For the Ohmic ignition scenario in FWL configuration under consideration, we can neglect P ICRH , while alpha + ohmic heating -variation of the plasma internal energy have been estimated to sum up to a total of ~ 20 MW [3]. As the power advected and conducted from the main plasma to the SOL is computed by ASPOEL starting from the above-mentioned Dirichlet BC, we can obtain from (1) the power radiated in the main plasma, which will contribute the background of the FWL heat load distribution, see below. Radial (anomalous) diffusivity values of 0.3, 0.3, 4.5 and 0.5 m 2 /s were assumed here for particles, momentum, electron and ion energy, respectively.

Results
We present here a selection of the results obtained from the coupled ASPOEL-EIRENE simulation. The runs are performed on a fixed grid with ~ 10000 nodes. Reduction of the number of nodes by a factor of 2 in each direction has a weak effect (e.g., a reduction of ~ 10 % of the peak heat flux on the FWL), so that the results below can be considered reasonably grid independent.
The plasma density is significantly affected by the presence of neutrals, with an increase by almost an order of magnitude with respect to the case without neutrals, especially in the region close to the top (and bottom) of the plasma chamber, and peak density more than a factor of 10 above the inner boundary BC, see figure 3. This is due to the strong recycling in that region, related in turn to the relatively high plasma flux, see figure 5 below, reaching the FWL in the only portion of it which significantly departs from tangency to the field lines: angles between poloidal field and the tangent to the wall attain an absolute maximum of almost 10 degrees in this region, see figure 6 below. The distribution of the electron temperature in front of the FWL, see figure 4, is also significantly affected by the neutrals, with ionizations leading to strong cooling of the plasma in the top (and bottom) regions of the machine. The ion temperature (not shown) behaves similarly. As, however, the temperature reduces by a factor of four, while the density increases by a factor of almost 10, the particle flux to the FWL, which scales as nT 1/2 , is also peaked near the top (bottom) of the machine, see figure 5.
For the same reasons, impurity production from sputtering at the FWL should be fairly limited in these conditions, as in the region of highest plasma outflux plasma ions reach the FWL with energies close to the Mo sputtering threshold.
The distribution of the heat load on the FWL, due to plasma advection and conduction from the SOL, as well as radiation from the main, is peaked close to the inboard equatorial plane, see figure 7. The load distribution is qualitatively similar to that computed in [7], the main quantitative difference being due to the inclusion here of the neutrals, as well as of the background radiation from the main plasma. About 2/3 of the power reaching the FWL in the plasma channel are due to the electrons. Since impurities in the SOL (and the associated radiation) are not included in the model, the computed peak load is a rather pessimistic estimate.

Conclusions
Neutrals have been included in the 2D computational model of the Ignitor SOL by coupling the edge plasma fluid code ASPOEL and the neutral Monte Carlo code EIRENE. The case of a firstwall/limiter configuration was considered, which constitutes a special challenge for edge plasma modeling. The computed peak heat flux on the FWL, close to the inboard equatorial plane of the machine, increases when neutrals are accounted for, but including impurities in the SOL model should reduce the peak value by a non-negligible amount.
A high density / low temperature region is created by ionization of recycling neutrals near the top / bottom of the machine. This temperature reduction could strongly limit the production of impurities by sputtering, but this feature will have to be checked by a self-consistent calculation.