Numerical studies of plasma edge in W7-X with 3D FINDIF code

Modelling plasma edge transport in a complex three-dimensional (3D) geometry of stellarators poses significant challenges as it requires considerable computational resources and a complicated numerical setup. Several 3D fluid plasma edge codes have been developed to handle the complex geometry of stellarators, including EMC3 [2] and FINDIF [1].

EMC3-EIRENE applies a fluid model for the boundary plasma based on Braginskii’s equations with free model parameters for perpendicular transport. Neutral gas is described in the frame of the kinetic transport in coupled EIRENE code. Plasma fluid equations are solved using enhanced Monte Carlo (EMC) solver. The EMC3 code can use also a symmetry of an analyzed device to simplify calculations. On the other hand, the FINDIF code solves fluid Braginskii’s equations using pseudo-time evolution. FINDIF is also designed to perform calculations in full geometry independently from a symmetry of an analyzed fusion device. Thus, there would be a possibility in future studies to investigate the effect of non-symmetric toroidally impurity gas injection.

The present study aims to show the applicability of FINDIF code in a relevant range of physics parameters. To validate the code, the Wendelstein 7-X (W7-X) limiter configuration cases were simulated and compared with the numerical (EMC3-EIRENE) [3, 4] and experimental [4] results for similar plasmas. The scrape-off layer (SOL) magnetic topology is described with its impact on plasma edge transport. The radial decay lengths for density, electron temperature and parallel heat flux are calculated and related to the simple SOL model. There have also been performed a comparison of the obtained results with EMC3-EIRENE code in order to validate effects of the FINDIF code improvement against more well-established code. In undertaking development of the presented FINDIF code, the aim is to obtain an additional 3D edge modelling tool, next to existing EMC3 code, which in the process of further development will be able to simulate the complex transport of multi-impurity plasma as well as effect of plasma drifts.

In the next two sections, the description of the FINDIF code and simulations of W7-X edge plasmas are addressed. Magnetic field and heat load calculations and SOL decay lengths estimations are discussed in the another two sections, respectively. Finally, the last section presents some conclusions and perspectives pertaining to further FINDIF developments.

FINDIF is a finite difference, three-dimensional, multifluid plasma edge transport code [5]. In the model the parallel transport follows Braginskii formulas [6] in B2 form [7], while the cross-field transport is assumed to be predominantly anomalous and modelled by diffusive approximation [8]. In the current code version, mesh creation is largely automatized [9]. FINDIF had been used previously to analyze the energy transport in the 3D problem of TEXTOR-DED tokamak [1] as well as for rudimentary results for W7-X: investigation of the effects of ergodicity on energy transport by comparing a vacuum and a finite β case [10] and comparison of results for different magnetic configurations [11].

Equations describing plasma transport in the fluid approximation are cast in a generalized form of the convective-diffusion equation [12]. The equations are discretized in field-aligned, non-orthogonal local magnetic coordinate systems [12] to eliminate numerical diffusion resulting from strong parallel/perpendicular transport anisotropy. Currently, four equations are solved self-consistently and provide us with plasma ion density (n_i), parallel velocity (v_i), and ion and electron temperature (T_i and T_e).

The transport model must necessarily be supplied with boundary conditions (BC) relevant to the problem. The Dirichlet conditions are applied at the inner (‘core’) and outer (‘wall’) boundaries while the Bohm sheath boundary condition is used at plasma-facing components (‘plate’). Namely, temperatures and plasma density at the inner boundary model core plasma act as a particle and energy source for the edge. Outer boundary (artificial wall in our cases) is in the area where the plasma is generally rarefied and cool. The outer BC must be set in such a way that they do not breed unphysical energy/particle/momentum fluxes from outer space into the SOL. In the immediate plate vicinity, ion velocity distributions are absolutely non-Maxwellian and strong electric fields are present. We account for the aforementioned physics using Bohm BC. Within these BC, the plate is treated as an ideal density sink (electrons are absorbed and ions neutralized). Plasma velocity is fixed to the local sound speed. The heat conduction flux at the plate is set proportional to the convection flux (coefficients being encoded in ‘sheath transmission factors’).

Equations for each variable are solved in an iterative, self-consistent manner using quasi-time evolution. The Dirichlet conditions are applied at the inner and outer boundaries while the Bohm sheath boundary condition is used at plasma-facing components. The numerical scheme is implicit along the lines and explicit across the field. The implicit scheme for fast transport offers better numerical stability and allows for larger time steps. The explicit treatment of cross-field terms reduces the computer memory demands.

Using non-orthogonal local magnetic coordinate systems directly impacts the form of a computational mesh, which consists of points along magnetic field lines. The resulting meshes are unstructured, non-orthogonal and non-uniform.

The W7-X stellarator is an optimized HELIAS-type stellarator with five-fold symmetry. The design aimed at the successful neoclassical transport reduction and efficient island divertor operation, in which the divertors are placed at the edge magnetic islands that originate from small radial field perturbations. However, as divertor graphite components had not yet been installed, to test the machine operation, W7-X was started in a limiter configuration in its initial plasma operation phase (OP1.1) [13–15]. The five poloidal, uncooled graphite limiters were placed on the inboard side in the bean-shape planes and a dedicated magnetic configuration without big islands in the SOL was chosen. It resulted in SOL consisting of three adjacent regions of flux tubes with different connection lengths [3]. The magnetic topology of the configuration is discussed in section “SOL decay lengths”. The mesh created for this configuration consists of approximately 3 million points, located at 120 toroidal cuts, which corresponds to a toroidal resolution of 3°, and building up approximately 12 000 magnetic field lines. It spans over the whole torus. Therefore, the new FINDIF tools properly capture the 3D SOL.

The core BC are imposed at inner simulation boundary (ISB), which coincides with a flux surface and is located inside last closed flux surface (LCFS), a few centimetres into the core. This approach is similar to the method employed concerning tokamak edge codes [7, 16]. We set X_ISB where X = n_i, v_i, T_i, T_e, common for all points at ISB. Due to the perturbing effect of the near-by SOL, the profiles at LCFS are far from flat, despite a strong smoothing effect of parallel transport. Deeper into the core, the distributions do get flatter. That is why it is acceptable to assume flat profiles on ISB. The values ‘at LCFS’ we refer to later are defined as flux surface averages.

For simplicity, we assume T_e,ISB = T_i,ISB = T_ISB and χ_e = χ_i = χ in the calculations. The simulations consisted of three types of scans:

– S1 – density scan with fixed T_ISB = 150 eV, D = 0.5 m²·s^–1, χ = 1.5 m²·s^–1 and n_ISB from 6 x 10¹⁸ m^–3 to 1.5 x 10²⁰ m^–3;

The S4 scan was performed at a constant D = 0.5 m²·s^–1, and χ = 1.5 m²·s^–1, whereas the n_ISB was varied in the range from 2.0 x 10¹⁹ m^–3 to 9.0 x 10¹⁹ m^–3 and the T_ISB was adjusted to keep fixed input power P = 3.2 MW.

While the choice of ISB is somewhat arbitrary, LCFS is unambiguously defined by the magnetic configuration. Any code–code or code–experiment comparisons must be done in terms of (average) temperature and density on the LCFS. Figure 1 presents all simulations labeled by (n_LCFS, T_LCFS). And indeed, fixed-T_ISB scan (S1) to a very good approximation can be thought of as a scan with fixed T_LCFS (clearly, T_LCFS < T_ISB). As for S2 scan, nLCFS visibly, but slowly, falls with boundary temperature. In case of S3 scan, varying anomalous transport coefficients relatively strongly, but monotonically, affects (n_IBS, T_IBS) to (n_LCFS, T_LCFS) mapping. As for S4, total power crossing LCFS is equal to P_inp, because neither external power sinks nor sources in the core region are implemented in FINDIF.

Fig. 1.

The insertion of limiters into the magnetic field causes the SOL to consist of three dominant types of lines as distinguished by target-to-target connection lengths. The connection length pattern on cuts φ_tor = 0° and φ_tor = 12° is shown in Fig. 2. The shortest field lines with L_C = 36 m (blue) start and end on the same limiter; the lines with L_C = 43 m (green) perform one toroidal turn and hit the neighboring limiter; and the longest field lines, L_C = 79 m (red), hit the neighboring limiter after two toroidal turns. On a cut, the regions for the first and second toroidal turns of the longest field lines are separated by a region of medium-length lines.

Fig. 2.

The magnetic topology strongly impacts the plasma parameters’ distribution. In Fig. 3, plasma pressure and Mach number on the cut are shown. The heterogeneous patterns of the pressure and plasma flows follow the complex magnetic topology of the SOL. These results are in agreement with the EMC3-EIRENE numerical findings [3], as the same pressure and Mach number variations were observed that correlate with the connection lengths.

Fig. 3.

Heat load on plasma-facing components, both its total and its distribution, is a crucial issue for designing and operating a fusion device. Therefore, it is of utmost importance that FINDIF properly captures the heat flux magnitude and distribution correctly. In FINDIF, the parallel heat flux reaching the limiter is calculated as 1qV=ncS(γeTe+γiTi)$${q_V} = n{c_S}\left( {{\gamma _{\rm{e}}}{T_{\rm{e}}} + {\gamma _{\rm{i}}}{T_{\rm{i}}}} \right)$$ where cs=(Te+Ti)/m$${c_s} = \sqrt {\left( {{T_e} + {T_{\rm{i}}}} \right)/m} $$ is the sound speed, and γ_e = 4.5 and γ_i = 3.5 are the sheath heat transmission coefficients for electrons and ions, respectively. To calculate the deposited heat flux, the magnetic field incident angle on the limiter surface is included as 2qdepo =q∥sinα$${q_{{\rm{depo\;}}}} = {q_\parallel }\sin \alpha $$

The heat flux distribution on the limiter surface calculated for the baseline is shown in Fig. 4. The similar pattern is obtained for all of the performed calculations with peak values in the range of 1.4–40 MW·m^–2. The vertical stripes present with two maxima, on the higher left and lower right regions. This characteristic pattern results from a few different factors. Firstly, the stripe of zero deposited heat in the middle appears because the magnetic field lines are parallel to the limiter watershed there. Secondly, the outside parts, being radially further away from the LCFS, also get reduced heat load. Finally, the heat maxima correspond to the regions of long magnetic field lines, which get more cross-field fuelling of heat and particles than the shorter ones. These outcomes are in agreement with the results described in the studies of Effenberg et al. [3] and Niemann et al. [4] qualitatively, matching magnetic field and heat load patterns.

Fig. 4.

The peak heat load increases with the corresponding variable parameter for S1, S2 and S3 scans (Fig. 5). To justify this correlation, P has been estimated for all simulations, showing that it also increases in the scans. For the S4 scan with fixed input power, the peak heat load decreases with increasing n_ISB. This outcome is consistent with the results described in the study of Effenberg et al. [3], where it was shown that for the same P, the higher the density at the separatrix, the lower the peak-deposited heat flux as the profiles across the limiter surface flatten.

Fig. 5.

The physics modelled by FINDIF corresponds well to the simple SOL model given in the study of Stangeby [17], which is justified for limiter cases, as the recycling hydrogen is more likely to be ionized inside the LCFS due to the direct vicinity of the target surface and the main plasma. The simple SOL model assumes deep plasma ionization and therefore the only particle source in SOL is due to cross-field transport. There is also no volumetric recombination, neutral friction, radiative and charge-exchange cooling included. In this model, the characteristic decay length of density is estimated using particle balance as 3λn=(DLCcS)1/2$${\lambda _{\rm{n}}} = {\left( {D{L_C}{c_S}} \right)^{1/2}}$$

In case q rises as nT^3/2 (compare Eq. (1)) the relation between plasma density, temperature and parallel heat flux decay lengths is given by 41λq=1λn+32λT$${1 \over {{\lambda _q}}} = {1 \over {{\lambda _{\rm{n}}}}} + {3 \over {2{\lambda _T}}}$$

In the following calculations, exponential decay of density, electron temperature and parallel heat flux with effective radius has been assumed. Effective radius is computed using the following expression: 5reff =∫0VdV′S(V′)$${r_{{\rm{eff\;}}}} = \smallint _0^V{{dV'} \over {S\left( {V'} \right)}}$$

The calculated density decay lengths for FINDIF simulations are plotted in Fig. 6 as points with colors denoting the connection lengths L_C = 36 m (blue), L_C = 43 m (green) and L_C = 79 m (red). Predictions of Eq. (3), scaled by a factor of 1.4 and using temperatures at a separatrix for calculation, are shown as solid lines with the same color coding. The density decay length is found to not change in the S1 scan (except for n_ISB below 10¹⁹ m^–3), decrease in the S2 scan and increase in both the S3 and S4 scans with the variable parameters n_ISB, T_ISB, D and n_ISB, respectively. These changes of λ_n correspond well to the trends predicted by the simple SOL model. The impact of connection length is overestimated in the SOL model related to values from the simulations, which is likely due to the averaging effect of transport between neighboring regions of different L_C.

Fig. 6.

The calculated electron temperature and parallel heat flux decay lengths for FINDIF simulations are plotted in Figs. 7 and 8. Both λ_Te and λ_q follow the same trends as λ_n except for n < 10¹⁹ m^–3, where electron temperature decay length decreases for lower densities, while parallel heat flux decay length is constant. The latter results from opposite low-density trends of λ_n and assuming the relation of Eq. (4) (for the simulations performed the approximate relation Te ≈ T_i holds).

Fig. 7.

Fig. 8.

In the study of Niemann et al. [4], several combinations of experimentally available parameters have been used to calculate a scaling of parallel heat flux decay length in the form of the power laws. For every set of used parameters, a connection length factor was included and a consistent value of ~0.22 for the respective exponent was determined. The scaling with n and T_e,edge showed little to no impact of these quantities. The same reference also states the scaling calculated using EMC3-EIRENE. Similarly to these experimental and numerical results, FINDIF also shows the impact of L_C on λ_q. In the study of Effenberg et al. [3], a λ_q proxy defined as decay length along the limiter surface, λ_sf, has been used. For EMC3-EIRENE simulations, in which a constant P has been taken as a boundary condition, it has been determined that λ_sf increases with n_LCFS, which is a trend observed also for the results of this work.

To compute the power input into the simulation domain, we perform numerical integration of the scalar product of heat flux density vector with a reversed (that is, pointing inwards) versor normal to ISB.

The results confirm the capability of FINDIF to simulate in a wide range of physical parameters. The dedicated tools for field line tracing and mesh creation are validated by comparison with previously published results, as the pattern of magnetic field lines on the limiter and toroidal cross-sections agree with cited results. The transport code in FINDIF, despite being a simpler physical model without neutrals and impurities, shows outcomes of physical quantities consistent with the references. The heat flux pattern on the limiter is impacted by the connection length, the exponential fall-off of heat flux and the magnetic field lines’ incident angle.

The peak heat flux density decreases with plasma density for a fixed power simulation scan. The pressure profile shows modulation with connection length due to stronger particle fuelling of longer magnetic field lines. The Mach number variation in the toroidal cross-section is also caused by the magnetic topology, as the midpoints, which closely match stagnation points, are shifted toroidally for regions of different connection lengths and of connecting different limiters. The results have shown that FINDIF qualitatively captures the heat flux distribution.

The density, temperature and parallel heat flux radial decay lengths trends are feasibly explained by comparison with the simple SOL model, as its assumptions closely match the model used in the FINDIF simulations. The decay lengths values are largely impacted by the connection length, which is a common result for the experimental and simulation references as well as the analytical model.

The model used currently in FINDIF, corresponding to the simple SOL model as shown, is suitable for the limiter case. The most important next step for FINDIF is to include neutral particles and impurity radiation to facilitate a better description of particle and energy transport in the boundary plasma.