JINTRAC integrated simulations of ITER scenarios including fuelling and divertor power flux control for H, He and DT plasmas

The ITER Research Plan [1] foresees a progression of plasma operational campaigns in hydrogen/helium, deuterium, and ultimately deuterium–tritium mixtures for the demonstration fusion power production. In the non-active operational phase (pre-fusion plasma operation) commissioning of the tokamak systems such as diagnostics, fuelling, heating and current drive systems, in vessel coils and plasma control systems will take place. Early operational experience and experimental results in hydrogen (H), helium (He), and then later in deuterium (D) will in addition give insights on approaches to optimise deuterium–tritium (DT) plasma scenarios. This is important to secure ITER's success in its primary objective of realising high fusion performance. The plasma density profile and its build-up are recognised to be very important in establishing the best route to achieve the target fusion gain, Q = 10, in ITER H-modes utilising the existing additional heating [1, 2].

In present tokamaks, gas puffing and recycling can effectively fuel the core plasma since neutral particles can pass through the relatively cool, sparse edge plasma. In contrast, the hotter and denser ITER edge plasma will be practically impenetrable to hydrogenic neutrals, which will not reach the separatrix the but will be ionised already in the far scrape-off-layer (SOL) [3–6]. Thus, pellet injection, although its deposition is rather peripheral in ITER, will be required to fuel the core plasma in ITER.

So far ITER integrated plasma simulation studies of fuelling and particle transport in ITER scenarios have been performed by introducing the sources from the scrape-off layer through boundary conditions/SOLPS parameterisations in core plasma modelling suites (e.g. [7–12]). Steady-state Q = 10 coupled core-edge transport calculations have also been carried out [13]. There has also been work focussed mainly on the magnetic evolution of the plasma where plasma transport was rudimentarily modelled without any coupling to the SOL/divertor plasma [14–16]. Here we here report for the first time on simulations including coupled core-edge transport calculations that follow the entire ITER plasma discharge evolution past the X-point formation for a wide range of ITER scenarios. This modelling approach allows us to study and optimize the fuelling of the plasma and to adjust the additional heating waveforms in H, He and DT plasmas to achieve the required stationary plasma conditions while avoiding operational limits, e.g. excessive neutral beam (NB) shine-through, excessive loads on the first wall panels and divertor power fluxes. In ITER the steady-state divertor power handing is limited to 10 MW m⁻², although fast transients below 20 MW m⁻² may be tolerated if their duration is short not to consume a sizeable lifetime of the ITER divertor due to recrystallisation causing mechanical fracture [17].

For non-active H and He plasmas we carried out L-mode simulations to identify robust scenarios for which the density is sufficiently high to commission and operate the NB at the highest beam voltages with acceptable first wall loads; namely a line-average density of 4.5 × 10¹⁹ m⁻³ in H and 3.0 × 10¹⁹ m⁻³ in He [18]. It is important to note that for fuelling of H plasmas both H gas fuelling and pellet fuelling are available. For He plasmas only gas fuelling can be used to fuel He; it is possible to inject H pellets, but this can lead to dilution and increase the density required to inject NB with acceptable shine-through loads on the first wall. In addition, the L–H power threshold also increases with the H dilution in He plasmas [19–22] and thus injecting H pellets into He plasmas can increase the required power to access H-mode beyond the capabilities in the initial phase of ITER operation, where H-mode plasmas are expected to be first explored.

In our simulations of 15 MA/5.3 T ITER baseline scenario plasmas with 50/50 DT composition, the main emphasis was the optimisation of the fuelling and impurity seeding to achieve high Q operational conditions with acceptable power loads to the first wall and divertor. This requires a sufficient density to inject the full power NB during the L-mode phase and an optimised fuelling scheme after the plasma enters the H-mode phase.

Previous studies have shown that an excessive density rise in the initial H-mode phase can decrease the core plasma temperature and prevent the fusion power build-up thus causing the return of the plasma to L-mode [23]. On the other hand, a too low plasma density in this phase can lead to the divertor power flux exceeding the design value of 10 MW m⁻². Therefore our H-mode access simulations have optimised the fuelling rate and impurity seeding rates to reach a robust H-mode in the minimal amount of time with high fusion gain (Q ∼ 5) and this is then followed by a higher fuelling rate in which the plasma density increases to ∼85%–90% of the Greenwald limit at which Q = 10 is expected to be achieved in ITER. Extensive fuelling optimisation studies for the stationary phases of ITER DT H-modes over a range of plasma currents and toroidal fields are described in [24].

While the control of the density and divertor power fluxes through divertor radiation by impurity seeding in the access to high Q operation requires some optimisation, this is much more challenging for the exit phase of high Q operation. In this phase the plasma density decreases, which can lead to a significant increase of the divertor power fluxes with values in excess of 10 MW m⁻². Thus, careful adjustment of both plasma fuelling (pellet and gas fuelling) and impurity (Ne) seeding is required in the exit phase to prevent excessive divertor power fluxes as well as the radiative collapse of the plasma solution due to uncurbed divertor radiation. This radiative collapse is reflected in our simulations by a numerical instability that corresponds to the experimental situation in which full plasma detachment takes place, possibly followed by a thermal instability such as a MARFE. Although in this paper we use the experimental terminology for this numerical instability, to speed up the simulations we have turned off part of the required physics (e.g. neutral–neutral collisions and some hydrogenic molecular processes) to model these physical phenomena in detail for ITER.

To model the above ITER plasma conditions and scenarios comprehensively and self-consistently, we have employed the integrated suite of core and SOL/divertor transport codes JINTRAC [25] developed at JET. This suite connects JETTO/SANCO [26, 27], a 1.5D core transport solver that includes impurities, with EDGE2D/EIRENE [28, 29], a 2D SOL/edge multi-fluid solver that comprises plasma interactions with the ITER Be wall and W divertor including atomic processes [30].

For simplicity, the only impurities that we have modelled are Be and Ne. At the time of writing, we are carrying out simulations including W, which will be presented in a future publication. We did not consider seeding N instead of Ne, even though the former seems to be a better divertor radiator in present experiments, since it can lead to tritiated ammonia formation, contributing to in-vessel tritium retention [31]. Indeed recent simulations show that N and Ne provide similar capabilities regarding divertor radiation in ITER due to the high plasma edge temperatures compared to present devices [32] and consequently Ne remains the preferred choice for ITER and in our simulations.

In the initial phase of our ITER JINTRAC simulation studies, we conducted a considerable benchmark between EDGE2D and SOLPS 4.3 [33], which showed that the SOL/divertor physics description was in agreement between the two codes [34] (see appendix B). Moreover, we have also significantly modified EDGE2D to model plasmas with He as the majority species. This has so far been validated against SOLPS 4.3 and there is a recent activity to validate it against JET He plasma experimental data. In this paper we present what to our knowledge are the first fully integrated core-edge simulations of He plasmas, applied to ITER in this specific case. Another important outcome of this code development work is that EDGE2D is now capable of handling plasmas with two main hydrogenic species, which is necessary for the evaluation of DT mixtures for different gas and pellet fuelling schemes, and recycling for DT plasmas.

The outline of the paper is as follows: in section 2 we present the JINTRAC modelling suite and our modelling assumptions. Our ITER simulation results for the non-active phase are described in section 3 and the simulations of DT plasmas in section 4 with a summary and conclusions in section 5.

The simulations that we present here follow most of the time evolution of the ITER DT and He H-mode reference plasma scenarios from the initial L-mode at X-point formation at low current through the build-up to the stationary ELMy H-mode phase, including the L–H transition, the initial ELM-free H-mode and, finally, the stationary ELMy H-mode phase with controlled ELMs. For DT plasmas we have also studied the H–L transitions, with particular emphasis on keeping the divertor power load and divertor plasma temperature under control. Overall, the main emphasis of the simulations in He and DT has been on the fuelling requirements to achieve an appropriate plasma density and to control its evolution for: (1) the unrestricted use of the NB in L-mode plasmas to trigger the L–H transition; (2) reaching and sustaining ELMy H-mode conditions; and (3) for 15 MA DT plasmas, reaching a Greenwald density fraction of 85%, which is required to obtain Q = 10, the ultimate goal of ITER.

For H plasmas a set of dedicated studies was performed to evaluate the possibility of achieving the minimum density required for the unrestricted use of the ITER NB at highest injection energy (870 keV for H neutrals), in L-mode plasmas. This is important for the use of the NB system in the initial phases of ITER operation; the required density is 4.3 × 10¹⁹ m⁻³ for H plasmas, which is substantially larger than for DT (2.5 × 10¹⁹ m⁻³). Hence, these simulations have been carried out at constant values of field and current rather than for the whole plasma scenario as was done for the He and DT H-mode plasmas.

Contrary to existing tokamaks, the high recycling regime in ITER (i.e. a much larger divertor flux than the outflux from the core plasma) persists even at low separatrix densities and high divertor temperatures [35], due to the large divertor plasma dimensions which lead to a very efficient neutral ionisation. Hence, for most plasma conditions in ITER a clear majority of the recycled neutrals are ionised in the divertor plasma and, typically less than 1% of the total neutral recycling flux actually reaches, and is ionised in the core plasma [5]. This renders gas fuelling and the ensuing increased recycling very ineffective to fuel the core plasma and to increase the core plasma density. Indeed, core plasma fuelling is foreseen to be performed by the injection of pellets in ITER. This is modelled with the modules NGPS [36] in L-mode and HPI2 [37] in H-mode. Both models are integrated in JETTO and describe the ablation and deposition of the pellet mass. HPI2 also takes into account the post-ablation drift and the fast transport phenomena due to the ∇B drift of the ablation cloud. As shown in figure 1 the EDGE2D/EIRENE model includes a pump and stationary conditions are established when the pumped neutral flux is equal to the gas + pellet fuelling flux and NB particle source. The 2D grid used in EDGE2D was calculated by GRID2D using a ITER reference equilibrium [38] covering the whole plasma volume. This reference equilibrium was developed within an iterative approach of exchanging engineering and physics information obtained from the magnetic and kinetic plasma scenarios developed by using DINA [39] and CORSICA [40].

Zoom In Zoom Out Reset image size

It should be noted that in the present work we have assumed that recycled neutrals can only cool the plasma and therefore do not contribute significantly to the recycling of energy in the SOL/divertor. This is a sound assumption for recycling on carbon or beryllium plasma facing materials, where fast reflections (from an ion returning as an atom) are less probable than thermal reflections (where the ion returns as part of a molecule). In addition, with a carbon wall the recycled power fluxes are also reduced by the formation of carbohydrate molecules and by the relatively high wall retention of hydrogenic species, which increases the power absorption of the wall. The lower retention in W walls and the small hydrogen to W mass ratio imply that fast reflections become more probable from a W wall, which increases the recycled energy flux. In appendix A we present an analysis of the impact of such effects by using a model that is more sophisticated in its description of the recycled energy flux from W than the one used in most of the studies in the paper. The assessment in appendix A shows that there is a minor effect on the core and SOL plasma of the increased energy recycling coefficient whereas in the divertor this effect is more significant and makes the divertor less detached for given upstream plasma conditions in the SOL. Hence, the presented simulations are more conservative (provide a lower limit) with respect to the value of the separatrix density at which full divertor detachment (leading to thermal instability) sets in, which could be slightly higher than our estimates due to W recycling effects. This is supported by a similar detailed analysis of material recycling effects on ITER divertor performance with SOLPS-ITER [41].

In addition to the modelling of the fuelling of ITER plasmas and the control of the density evolution, our simulations have also addressed the control of divertor power loads. In ITER it is required to maintain the divertor peak power fluxes below their design limit of 10 MW m⁻². For a range of plasma scenarios conditions (typically high power/high current plasmas) the achievement of this power flux level has necessitated the introduction of Ne into the plasma, which has been modelled by a gas puff in the outer divertor target region. The resulting Ne radiation in the SOL/divertor can be very effective in reducing the power flux at the target and decrease the divertor plasma temperature, which is important to limit W sputtering. The latter is not included in our simulations but an ion temperature of ∼5 eV at the location of maximum ion flux is taken as a proxy for significant W sputtering to take place.

For the lower power phases of the plasma scenarios (e.g. L-mode phases) we have frequently observed that injection of Ne was not required to achieve power fluxes under the ITER design limit; this could be achieved by increasing the edge plasma density by gas fuelling alone. On the other hand, in the phases where large SOL power fluxes occur together with a fast evolution of the edge plasma parameters (e.g. H-mode termination phases), it is extremely complex to adjust gas fuelling and Ne injection to maintain the divertor power flux under the design limit while avoiding a strong detachment of the divertor. This is due to the strong coupling between core and edge plasma (e.g. through the edge power flux that determines the level of confinement which, in turn, is affected by the increased core plasma radiation from the injected Ne). In this paper we present the first results of the control of the divertor power fluxes in such phases, which has been achieved by manually tuning the gas fuelling and Ne injection level by trial and error. A more systematic study of power flux control in H-mode termination phases in ITER scenarios requires appropriate feedback loops for gas puffing and Ne injection to be implemented in JINTRAC; this work is in progress and will be used for future ITER simulations.

The simulations presented here consider a single isotopic species in the SOL and divertor to represent D and T. At the separatrix the resulting ion density and neutral influx from EDGE2D is transferred to the core plasma model in JETTO as being composed of 50% D and 50% T. On the other hand, the total DT plasma particle outflux crossing the separatrix from JETTO is transformed into a single hydrogenic isotopic specie when transferred to EDGE2D. Successful core + SOL/divertor test simulations of JET DT plasmas have been carried out with a more recent version of EDGE2D, which can handle plasmas with two main hydrogenic species [42]. The results of this version will be benchmarked against JET hydrogen–deuterium experimental plasmas before it is used predictively for ITER.

For the simulations presented here we have developed a special version [5, 24] of the standard JETTO Bohm-gyro-Bohm transport model tuned to GLF23 [43], by including a collisionality dependent particle pinch, v_pinch(ν*) = f(ν*) ⋅ 0.5Dr/a², where D is the cross-field particle diffusivity, ν* is the normalized electron collisionality and $f\left({\nu }^{\ast }\right)=\mathrm{min}\left(1,\mathrm{exp}(1-{\nu }^{\ast }/{v}_{\text{th}})\right.$.The constant v_th = 0.04 was obtained by matching GLF23 density predictions to BgB density predictions of ITER plasmas at different current and field. This pinch reproduces the density peaking observed in GLF23 associated with anomalous transport fluxes at low collisionalities. The transport of the impurities, Be sputtered from the main chamber wall and Ne from seeding to keep the divertor power fluxes acceptable are described by SANCO for the core plasma. JINTRAC assumes that the impurity transport can be described by the linear combination of neoclassical transport of the impurities plus anomalous transport which is assumed to be similar to that of the main ion species. We have used NCLASS [44] to calculate the neoclassical transport for all particle species (main ions and impurities).

As there is currently no simple model that accurately can describe the cross-field transport in the SOL, we have assumed in EDGE2D that the transport has a tanh decay from the separatrix values given by the core transport model assumptions, to the fixed values, D = 0.3 m² s⁻¹ and χ_e = χ_i = 1 m² s⁻¹ in the far SOL, frequently used in ITER SOL simulations. Additional effects from plasma drifts in the SOL/divertor have not been taking into account. If an edge transport barrier is present, the transport coefficients at the separatrix are extended 5 mm into the SOL. Moreover, simulations of JET L-mode plasmas in the high-recycling regime show that transport in this range, specifically flat profiles with χ_i,e = 0.5 m² s⁻¹, D_⊥ = 0.15 m² s⁻¹, is needed in order to match the density and temperature at the outer midplane separatrix, and at the outer divertor target [45, 46]. The cross-field transport we have used is conservative with respect separatrix densities and onset of detachment, as simulations with lower transport, still within experimental observations in tokamaks, obtained higher core densities for the same gas rate and heating power, and were less prone to detachment [6].

Even if we have used ESCO, the internal fixed-boundary equilibrium solver in JETTO, to advance the core equilibrium in our modelling, we had to fix the poloidal shape of the EDGE2D grid for technical reasons. Consequently, in our current ramp-up and ramp-down simulations the plasma poloidal shape is assumed to remain unchanged and we have altered the pitch of the magnetic field in the SOL to reproduce the change in connection length associated with the changes to the poloidal field as the plasma current varies during these phases.

The auxiliary heating schemes in the ITER baseline are radio frequency (RF) heating by ECRH (electron cyclotron resonant heating) and ICRF (ion cyclotron radio frequency heating) and NB (neutral beam) injection, all of which are included in the modelling. We have applied PENCIL [47] to calculate the NB power deposition profiles in DT plasmas and ASCOT [48] for He plasmas to get accurate power deposition and beam shine-through.

ITER considers a range of ICRF heating schemes which range from He³ minority and second harmonic T at a toroidal field value of 5.3 T to H minority at 2.65 T for He plasmas. Although the module PION [49], included in JINTRAC, is, in principle, able to calculate the ICRF heating profiles, in the simulations presented here we have, for simplicity, adopted an ad-hoc Gaussian ICRF deposition profile centred on-axis with the power equally allocated to the ions and electrons; this is a good approximation for DT plasmas at 5.3 T while it overestimates slightly the core ion heating for H minority of He plasmas.

For ECRH heating we have also applied an ad-hoc Gaussian power deposition but this time centred at a normalised minor radius 0.2. The recent installation of GRAY [50] in JINTRAC will allow to calculate the ECRH heating and current drive self-consistently in future ITER studies. However, as the main aim of this study is to determine fuelling requirements for ITER scenarios and to investigate schemes for divertor power load control, the use of ad-hoc RF deposition profiles is an acceptable simplification.

In DT plasmas we took the α-particle heating into account by applying the Mako model [51, 52] since the additional heating it produces is essential for the 15 MA/5.3 T DT plasma to reach and sustain H-mode. Moreover, the He ash was modelled self-consistently in the core and SOL/divertor plasma and pumped out in the divertor.

Since at this time, no first-principle model exists that can correctly predict confinement transitions, we have relied on experimental observations and scaling laws to determine the power level at which the plasma accesses or exits the H-mode confinement regime in ITER scenarios. One such scaling law is the Martin L–H threshold scaling [53]

$$\begin{equation}{P}_{\mathrm{L}\mathrm{H}}=2.15{e}^{\pm 0.107}{n}_{\mathrm{e}20}^{0.782\pm 0.035}{B}_{\mathrm{T}}^{0.803\pm 0.032}{S}^{0.941\pm 0.019}\end{equation} \tag{ 1 }$$

where P_LH is the threshold power in MW, n_e20 the line-average density in 10²⁰ m⁻³, B_T the magnetic field in T, and S the plasma surface area in m². An L–H confinement transition is assumed to occur in our simulations when

$$\begin{equation}{P}_{\text{Loss}}-{P}_{\mathrm{L}\mathrm{H}}\geqslant 0,\end{equation} \tag{ 2 }$$

where P_Loss, the loss power crossing the separatrix is

$$\begin{equation}{P}_{\text{Loss}}={P}_{\text{OHM}}+{P}_{\text{abs}}-\mathrm{d}W/\mathrm{d}t-{P}_{\text{rad}}-{P}_{\text{Floss}}\end{equation} \tag{ 3 }$$

and P_OHM is the ohmic power, P_abs is the absorbed power, W is the total plasma energy, P_rad is the total radiated power and P_Floss is the fast ion losses. Conversely an H–L transition occurs in the simulations when P_Loss − P_LH ⩽ 0.

The Martin scaling was derived for deuterium plasmas. Experimental observations in AUG, DIII-D, and JET [20–22] indicate that for He plasmas the L–H transition threshold is higher by a factor of 1.44 than in D plasmas, a factor which we have applied in our simulations of L–H transitions in He.

The edge transport barrier (ETB) width has been determined with the EPED1 model [54, 55] and we have assumed this value to remain fixed throughout the L–H and H–L transitions, i.e. whenever P_Loss − P_LH ⩾ 0. In addition, we allow the level of remaining anomalous transport within the ETB to depend on the surplus of heat flux through the separatrix, P_Loss, compared to power threshold P_LH so that when

$$\begin{equation}{P}_{\text{loss}}-{P}_{\mathrm{L}\mathrm{H}} > 0\end{equation} \tag{ 4 }$$

$$\begin{align}\hfill \chi ,{D}_{\text{ETB}}& =\chi ,{D}_{\text{anomolous}}{e}^{-\left({P}_{\text{loss}}-{P}_{\mathrm{L}\mathrm{H}}\right)/\left({P}_{\mathrm{L}\mathrm{H}}{{\Delta}}_{\mathrm{L}\mathrm{H}}\right)}\hfill \\ \hfill & \quad +\chi ,{D}_{\text{Neoclassical}},\hfill \end{align} \tag{ 5 }$$

where Δ_LH defines the reduction of the anomalous transport with the edge power flow margin over the power threshold. To be consistent with JET experiments we have set Δ_LH = 0.1 [56], which also allows for smoother transitions with less dithering between L- and H-mode.

Following the L–H transition the pressure in the ETB rises and the ballooning stability parameter $\alpha =-(2{\mu }_{0}R{q}^{2}/{B}^{2})\left.\mathrm{d}p/\mathrm{d}r\right)$ rises. When it gets to a specified value, α_crit, evaluated on the basis of edge ideal MHD stability with the codes MISHKA [57] and ELITE [58], we trigger the continuous ELM model [59]. This consists of an ad hoc elevation of transport in the ETB to maintain α = α_crit throughout the stationary phase of the ELMy H-mode. This is an approximation of the expected plasma behaviour of ITER ELMy H-modes where ELM energy losses need to be controlled to very small values through ELM suppression [60].

Lastly, we have applied the Kadomtsev model [61] to model the sawtooth reconnection. A sawtooth crash is prompted when the safety factor drops below unity at the magnetic axis. Although, a fixed minimum time (0.1 s) must have passed since the previous crash to avoid too frequent crashes destabilising the code. However, in most of our simulations sawtooth crashes normally occur much more infrequently than this.

We have concentrated our non-active fuelling modelling studies on the achievement of appropriate L-mode densities in H (4.3 × 10¹⁹ m⁻³) and He (2.7 × 10¹⁹ m⁻³) to permit unrestricted operation and commissioning of the NB over the whole range of injected hydrogen ion energy (500–870 keV) up to the highest values, and on the achievement of ELMy H-mode plasmas in He. Experimental evidence from several tokamaks for hydrogenic plasmas, such as JET [19], show that P_LH ∝ 1/A_eff where A_eff is the effective atomic mass and, hence, the threshold power for hydrogen plasmas is expected to be significantly higher than for DT. For ITER this implies that H H-modes are improbable at heating powers below 100 MW at full current and field. Therefore, H-mode studies in the non-active phase are expected to take place at reduced currents and fields and most likely in He to take advantage of the reduced power threshold compared to H plasmas (${P}_{\mathrm{L}\mathrm{H}}^{\mathrm{H}}=1.4{P}_{\mathrm{L}\mathrm{H}}^{\mathrm{H}\mathrm{e}}=2{P}_{\mathrm{L}\mathrm{H}}^{\mathrm{D}}$). Predictions for 7.5 MA/2.65 T H-modes in ITER show that the ELMs could generate energy loads above the melting limits for the W divertor monoblock target and, hence, such plasmas provide highly relevant conditions for the understanding of H-mode plasma behaviour in ITER and for the commissioning and refinement of the ELM control system to demonstrate that ITER can achieve the required level of ELM control before the start of nuclear operations [1].

3.1. Hydrogen L-mode simulations of 15 MA/5.3 T plasmas

Our simulations of 10 MA and 15 MA hydrogen L-mode plasmas at 5.3 T foresee that fuelling with gas only will be insufficient to obtain the density of 4.3 × 10¹⁹ m⁻³ required for unrestricted NBI operation, even when the plasmas are heated with the maximum baseline level of 40 MW RF power. Indeed, for both currents we only reached line-average densities of up to ∼2.8–3.1 × 10¹⁹ m⁻³ for a gas rate between 2.5–3.0 × 10²² s⁻¹ (figure 2). For higher gas fuelling rates the divertor strongly detaches, which could lead to the triggering of radiation instabilities such as a MARFE. Moreover, the separatrix density saturates and consequently the core density does not increase with gas fuelling beyond these rates.

Zoom In Zoom Out Reset image size

To increase the density further, we have switched to dominant pellet fuelling by reducing the gas rate to a very low value of 5.0 × 10²⁰ s⁻¹ and injecting pellets (90 mm³/5.3 × 10²¹ particles) into the plasma. By applying an initial pellet frequency of 3.3 Hz and then increase it to 4 Hz at 182.5 s, a 15 MA/5.3 T hydrogen L-mode plasma can reach a line-average density of ∼4.6 × 10¹⁹ m⁻³ even when it is heated by only 20 MW of ECRH (figure 2). To ensure that this plasma remains thermally stable the gas fuelling needs to be kept at very low values. As a rule of thumb, target temperatures less than 1 eV indicate a high level of detachment in our simulations and a likelihood of high recombination rates, the ionisation front creeping close to the x-point and thermal instability of the solution. For these pellet-fuelled hydrogen plasmas, the inner target temperature remains above 1 eV although the outer target temperature can momentarily decrease below this value due to the intermittent particle outflux from the plasma caused by the injected pellets. It is important to note that these simulations do not contain extrinsic impurities such as Ne seeding because the divertor power loads remained well under 10 MW m⁻²; the Be content was also very low in these simulations with Z_eff ≈ 1 throughout the core and edge plasma.

3.2. Simulation of He H-mode scenarios at 7.5 MA/2.65 T

We have performed simulations of ITER helium H-mode plasma 7.5 MA/2.65 T scenarios including the plasma current ramp-up right after X-point formation and the access to stationary H-mode. The starting conditions correspond to a diverted He L-mode plasma at 3 MA/2.65 T with 20 MW ECRH heating and a volume-average electron density of 1.0 × 10¹⁹ m⁻³; the modelled ramp-rate is 200 kA s⁻¹ leading to the 7.5 MA flat top current being reached 22.5 s after X-point formation. To model the current carrying equilibrium we rescaled the poloidal flux ψ with increasing current and adjusted correspondingly the edge field line pitch in the SOL while maintaining a constant plasma shape. The latter is also used as boundary condition for ESCO, which regularly updates the core equilibrium. We also considered different levels of gas puff to investigate the resulting increase in density. From the results of the simulations it is clear that the density ramp rate during the current ramp-up is exceptionally sensitive to the prescribed gas rate. If the rate is too high (>6 × 10²¹ s⁻¹), it can cause early detachment of the divertor and thermal instability of the numerical solution, whereas if it is too low (<1 × 10²¹ 1 s⁻¹), the density cannot be increased together with the plasma current as required for the later access to H-mode in this scenario. To ensure the achievement of the required density we have implemented a numerical feedback scheme in the JINTRAC suite to change the level gas puff to provide a density ramp at a given Greenwald fraction, which we have chosen to be ∼50%.

In the He current ramp-up simulation at about 236 s, the density has already almost risen to the target value that is modelled for steady-state L-mode plasmas at 7.5 MA/2.65 T (figure 3). Because the plasma current had gotten to no more than 5.5 MA, the simulation should be prolonged at a constant density with varying current to obtain the flat top plasma current value of 7.5 MA. However, at this core plasma density, current and ECRH power the inner divertor was nearly detached, so we could not increase the current and finish the simulation due to numerical instabilities. Future work will look at how to fine-tune the He fuelling rate to achieve a viable current ramp-up. Note that the steady-state He L-mode simulations at 7.5 MA/2.65 T with 20 MW of ECRH could only reach line-average densities of ∼2.1 × 10¹⁹ m⁻³ with a He gas fuelling rate of 7.0 × 10²¹ s⁻¹ before the He radiation front in the divertor started to move to the X-point and render the divertor conditions unstable.

Zoom In Zoom Out Reset image size

This density is lower than the target feedback value for the current ramp at 7.5 MA; 50% of the Greenwald limit corresponds to ∼3.0 × 10¹⁹ m⁻³, which is similar to the required value for unrestricted NB operation (at full power of 33 MW) in He plasmas, which in turn is essential for H-mode access in these plasma conditions.

In order to determine if shine-through losses would be a problem in the H-mode access phase, we have modelled the L–H transition density rise after 33 MW NB power is injected in a 7.5 MA/2.65 T He plasma with 20 MW of ECRH heating and a line-average density of ∼2.1 × 10¹⁹ m⁻³. The NB deposition and shine-through losses (${P}_{\mathrm{N}\mathrm{B}}^{\mathrm{S}\mathrm{T}}$) have been evaluated with ASCOT and, as shown in figure 4, it is found that the plasma density increases rather quickly after the H-mode transition so that shine-through losses reduce from 20% to 3%–6% in about 5 s. We therefore expect that the resulting power fluxes on the first wall remain acceptable in this phase [18].

Zoom In Zoom Out Reset image size

To test the sensitivity of H-mode access in He to threshold assumptions, we have performed simulations for two values of the L–H threshold in He; one optimistic assuming the same threshold for He and D plasmas as observed in ASDEX-Upgrade, and another more realistic, assuming that the He H-mode threshold is higher than that of D by a factor 1.4 as determined through experiments on JET, DIII-D, Alcator C-Mod, etc. As anticipated, the lesser L–H threshold in the optimistic case permits a swifter transition to ELMy H-mode (figure 4) with the edge stability ballooning parameter reaching α = α_crit = 2.1 in ∼3.5 s after switching on the NB. On the other hand, after 5 s the case with the a less optimistic threshold still remains in the ELM-free phase of the transition (α ≈ 1.4 < α_crit = 2.1) as the remnant anomalous transport remains too high to achieve significantly improved confinement. In this case a faster access to ELMy H-mode would require the use of the available 20 MW ICRH besides the 33 MW of NB and 20 MW of ECRH used in these simulations.

The current ramp-up and ramp-down phases are important for the achievement of ITER fusion production goals. We have therefore performed full simulations of the Q = 10 baseline scenario from X-point formation in the current ramp-up to X-point to limiter transition in the current ramp-down. A full account of these simulations can be found in [62]. During the current ramp-up the plasma current is increased towards the target value for the high Q burning plasma phase. In the ramp-down phase the plasma must exit burning conditions and needs to be terminated in a controlled way while avoiding disruptions both when the termination can be performed over long timescales and when the timescales are shorter following, for instance, malfunctions of one or more ancillary systems (e.g. heating and current drive systems).

It has been observed that the SOL width varies strongly with plasma current in H-mode and likely in L-mode as well. In H-mode, our simulations capture this effect as the continuous ELM model determines the H-mode edge transport, which is tuned to stay at the MHD limit once the pressure (α_crit) is high enough. Therefore, the effective edge transport coefficients that extend 5 mm into the SOL in H-mode, scale inversely with I_p and thus so will the edge e-folding lengths [24]. These will be similarly affected in L-mode, as in this case the edge transport is dominated by the Bohm contribution, which is inversely proportional to the local magnetic field and thus also to the enclosed plasma current.

Our simulations of the baseline scenarios have been performed from early X-point formation at 3 MA/5.3 T at 10 s, which allows the application of auxiliary heating to lower the flux consumption in the ramp-up phase. The simulated current ramp-up rate of 200 kA s⁻¹ leads to 15 MA being reached at 70 s (figure 5, left), which is well within the capabilities of the ITER poloidal field system. At this stage the additional heating is increased to trigger the H-mode transition which finally results in flat-top 15 MA/5.3 T ELMy H-mode plasma condition with a Greenwald density fraction of 85%, and a fusion gain of Q = 10. Our simulations have demonstrated that such a plasma scenario is indeed feasible within the physics basis and models included in the JINTRAC modelling suite (which has been compared to experiment) and without exceeding physics or design limits of the ITER components (e.g. divertor power fluxes). This has been demonstrated for both flat-top conditions corresponding to specific instants of the Q = 10 scenario as well as for the scenario as a whole, as illustrated in (figure 5).

Zoom In Zoom Out Reset image size

The studies of the achievable density in ITER DT L-modes show that for a given input power, the core density builds up with the gas fuelling rate and then saturates [5, 6, 63], following the same qualitative trends as for H and He plasmas in section 3. If the gas rate is further increased, the divertor plasma enters a deeply detached state that leads to the thermal instability of the solution and the termination of the simulation. Our JINTRAC modelling suggests that in the current ramp-up of the Q = 10, a Greenwald density fraction of less than 35% can be anticipated with gas fuelling and merely 20 MW of RF heating [6, 63]. With more RF heating higher fractions can be achieved with gas fuelling only, especially at lower currents. Nevertheless, it is more effective to use pellet fuelling if higher densities are required during the ramp-up, in qualitative agreement with the H L-mode simulations in section 3. It is important to note, however, that during the DT phase the ITER NB will be operated in D at 1 MeV (for the full power of 33 MW) and that the density required for unrestricted NB operation in DT plasmas (∼2.2 × 10¹⁹ m⁻³) is much lower than for He and H so that at ∼35% of n_GW, unrestricted NB operation is possible for currents above 8 MA. At present NB injection at such low levels of plasma currents is not considered in the Q = 10 scenario development. Still, if necessary, full NB power at lower currents would be feasible by dominant pellet fuelling as demonstrated in section 3.

Aligned with prior JINTRAC core-only predictions [23], our fully integrated core + SOL/divertor simulations indicate that the build-up of the fusion power, and corresponding α heating, is essential (on top of additional heating) to provide sufficient margin above the L–H threshold power in order to access quasi-stationary ELMy H-mode burning plasma conditions [63]. This requires a careful optimisation of the plasma fuelling in the H-mode access phase so that the core ion temperature can increase to values (>10 keV) at which DT fusion reaction rates are significant and can contribute to the total heating power. This leads to further improvement of the H-mode performance through the concomitant reduction of the transport in the ETB. This is illustrated in figure 5 (left) where two pellet fuelling rates are applied after the L–H transition in the ITER Q = 10 scenario. For the highest rate the plasma returns to L-mode after a brief, ELM-free phase following the H-mode transition because of the insufficient α heating and edge power flow.

In this Q = 10 scenario the divertor power loads remains under the 10 MW m⁻² limit without utilising extrinsic impurity seeding (Ne) in the ELM-free phase up to Q ∼ 5. For higher Q phases and quasi-stationary conditions, Ne seeding is essential to maintain acceptable divertor power fluxes. Undeniably, the need to maintain sufficient divertor radiation to lower the divertor power flux to tolerable values introduces specific complications when considered together with core plasma fuelling with pellets and the associated particle flux transients.

It should be noted that even if we only show the total Ne concentration in the SOL and divertor in figure 5, the local Ne levels are determined by the SOL/divertor transport and, in general, increases towards the separatrix and higher densities is also seen at the inner divertor leg. In addition our simulations show that although Ne penetrates to the magnetic axis, the total core penetration of Ne remains low (<0.24%) in the current ramp-up simulations and therefore its associated radiation (<5 MW) does not adversely affect the quality of the H-mode or produced fusion power.

Our simulations of pellet fuelling of 15 MA/5.3 T ITER ELMy H-mode Q = 10 plasmas with acceptable divertor power fluxes, which requires a separatrix density ∼5.5 × 10¹⁹ m⁻³ and Ne seeding, show that the pellet size must be carefully optimised to achieve a fully integrated scenario (including fuelling, particle and power exhaust); pellets with a particle content larger than 5.5 × 10²¹ atoms lead to deep transient divertor detachment after pellet injection and the thermal collapse of the plasma solution [24, 63].

The current ramp-down is typically expected to be twice as long as the ramp-up phase. This is to allow for a slow exit from H-mode and ramp-down of the current, which ensures good vertical stability control (not included in simulations). At the beginning or during the current ramp-down phase the plasma exits the burning conditions. Throughout the ramp-down itself, the plasma current is reduced carefully to guarantee a controlled termination that avoids vertical stabilities and disruptions [64]. This is a complex phase as the control of vertical stability is linked to plasma shape and radial plasma position control. Consequently fast losses of plasma energy (such as during H–L transitions) that can trigger both large divertor power fluxes and radial inwards shifts of the plasma, might lead to the loss of vertical stability, if the limits of the vertical stabilisation systems in ITER are exceeded.

To study plasma behaviour in the ramp-down phase we have considered a current ramp-down evolution from 15 MA/5.3 T Q = 10 burning plasma to 3 MA/5.3 T in ∼90–120 s, depending on the plasma current at which the H–L transition takes place (figure 5, right). Our model of the ETB transport modification with a margin above the H-mode threshold (equation (5)) results in a period of ∼5 s in which the plasma energy decreases from H-mode values to L-mode values in the exit phase. This timescale agrees with present experiments [23, 65] and allows the sustainment of good radial position control in this phase.

On the other hand, power fluxes to the divertor can be very substantial in the initial phase of the H-mode termination. To avoid the divertor power fluxes to exceed the limit of 10 MW m⁻² and W sputtering to significantly increase (W is not modelled in our simulations but an ion temperature of ∼5 eV at the location of maximum ion flux is taken as a proxy for significant W sputtering to take place), it is necessary to use Ne seeding.

In later, low power phases of the ramp down, where the edge power flow decreased, we had to carefully tune down or even turn off Ne seeding to ensure that the Ne level in the plasma is reduced by the pumping system sufficiently fast to avoid deep detachment of the divertor. In our modelling we did the tuning of plasma fuelling, Ne seeding, etc, by manual adjustment. Figure 5 (right) illustrates the results for two examples: a complete switch-off of the additional heating at the end of the 15 MA flat top; and for an H–L transition in the current ramp-down itself (at around 10 MA).

In conclusion, these simulations show that integrated scenarios for burn termination with stable plasma behaviour and acceptable divertor power fluxes can be obtained with the ITER heating & current drive and fuelling/impurity seeding systems. However, the ambition to progress further in the optimisation of this complex phase has prompted the development of more sophisticated divertor power load feedback control loops in the JINTRAC modelling suite. This will avoid the computationally intensive trial and error approach that we followed to achieve the results in figure 5.

In this paper, we have presented integrated core-SOL/divertor transport modelling with JINTRAC of the non-active and nuclear phases of ITER plasma scenarios, as well as the first core-SOL/divertor simulations of He plasmas and of the full DT Q = 10 scenario including the current ramp-up/down diverted phases. In these studies of ITER plasmas, avoiding deep detachment was a primary consideration and indeed often the limiting factor in the our efforts to reach the desired densities.

We have demonstrated that pellet fuelling is the most efficient way to increase plasma density in L-mode plasmas and that it will be required to provide unrestricted NB operation in 15 MA/5.3 T hydrogen L-mode plasmas in order to avoid full divertor detachment. For He plasmas gas fuelling alone is sufficient to provide access to H-mode conditions at 7.5 MA/2.65 T by reducing the shine-through losses to operationally acceptable values through the associated density increase in the L–H transitions. Nevertheless, as the higher density also raises the L–H power threshold, reaching ELMy H-mode with only 33 MW NB and 20 MW ECRH is uncertain and will depend on the actual He L–H power threshold in ITER. For DT plasmas the use of pellet fuelling is less essential for unrestricted NB operation due to the lower NB shine-through losses in such conditions.

Only when using pellet fuelling, our results predict that, commissioning and operation of the ITER NB system to full power is possible in 15 MA/5.3 T L-mode hydrogen plasmas with 20 MW of ECRH. Similarly, He H-mode plasmas at 7.5 MA/2.65 T with an additional heating level of 53–73 MW can be achieved, which allows the characterisation of H-mode plasmas and the demonstration of ELM control schemes in the non-active phase before ITER DT operation.

Our integrated simulations demonstrate that access to high Q conditions in ITER DT plasmas requires an optimised build-up of the α heating after the L–H transition because the margin above the L–H threshold provided by the additional heating is moderate for high current 5.3 T DT plasmas in ITER. This optimised α heating build-up can be achieved by tuning the density rise through pellet fuelling with an initial low fuelling rate phase to allow the ion temperature to build up in the plasma to, typically Q ∼ 5, followed by a larger rate to provide the density rise to the quasi-stationary Q = 10 ELMy H-mode conditions with 85% of n_GW.

For the high Q (typically Q > 5) transient and quasi-stationary phases to reach the Q = 10 scenario, Ne seeding is required to maintain the divertor power fluxes under the 10 MW m⁻² limit. This introduces a complex interaction between power flux control, particle control and fuelling when pellets are used to fuel the plasma due to the intrinsic particle fluxes that they cause. Pellets beyond a given size can trigger transient deep divertor detachment leading to the thermal collapse of the solution. The physics is subtle, and fully integrated simulations are therefore required when considering fuelling optimisation of ITER high Q plasmas; individual optimisation of separate aspects of particle and power control may be misleading.

Our integrated JINTRAC simulations of the exit from burn in the Q = 10 scenario show that the flexibility of the ITER systems can provide conditions in which plasma position control can be maintained and power fluxes to the divertor remain under the acceptable limit. This involves careful optimisation of fuelling and impurity (Ne) seeding to ensure it provides the necessary divertor radiation in the phases with high edge power flow (initial phase of the H-mode termination) and that Ne content is sufficiently low in later phases of the termination when the edge power fluxes are lower (end of H-mode and L-mode phases).

While these initial simulations indicate that the flexibility of the ITER fuelling, heating and impurity seeding systems is appropriate to ensure a robust access to/exit from and sustainment of Q = 10 burning plasmas, the plan for a detailed optimisation of heating and fuelling for a wide range of ITER high Q scenarios have prompted new feedback loops to be implemented in JINTRAC to improve the productivity of the code by avoiding much of the trial and error that has been required to perform the present work.

The main author wants to thank Xavier Bonnin at ITER for valuable insights and remarks on the manuscript. JINTRAC was used under license agreement between Euratom and CCFE, Ref. Ares(2014)3576010—28/10/2014. The full SOLPS simulation data mentioned appendix B are obtainable from the public ITER IMAS database. This work was funded jointly by the RCUK Energy Programme (Grant No. EP/I501045) and by ITER Task Agreement C19TD51FE implemented by Fusion for Energy under Grant GRT-502. To obtain further information on the data and models underlying this paper, whose release may be subject to commercial restrictions, please contact . The views and opinions expressed do not necessarily reflect those of Fusion for Energy which is not liable for any use that may be made of the information contained herein. The views and opinions expressed herein do not necessarily reflect those of the ITER Organization or Fusion for Energy.

A.1. Introduction

A.2. L-mode comparisons for ITER H plasmas

In this case we study the effect of the different modelling assumption of energy recycling as described in appendix A.1 on a pellet-fuelled 10 MA/5.3 T pure hydrogen plasma with 20 MW of ECRH. To ensure acceptable divertor power fluxes, gas is fuelled at a rate of 5.0 × 10²⁰ s⁻¹ besides the fuelling by pellets (90 mm³/5.3 × 10²¹ particles), which we inject at a frequency of 4 Hz.

Figure A.1 shows that the volume-averaged core ion density and on-axis ion temperature are not affected by the divertor energy recycling assumptions but that there is a very minor effect on the outer-mid-plane (OMP) electron and ion densities and temperatures. This is supported by the nearly identical OMP profiles at the end of the simulation (figure A.2). Only at the divertor plates can we observe the impact of the different energy recycling assumptions. At the outer divertor (figure A.3) the densities and temperatures are somewhat higher for the correct W energy recycling assumptions, whereas at the inner divertor (figure A.4) we only see an increase in the target plasma temperatures. We have also obtained very similar results for pellet-fuelled (90 mm³/5.3 × 10²¹ particles, injection frequency 3.3–4 Hz) 15 MA/5.3 T pure hydrogen plasma simulations with 20 MW of ECRH and with a gas rate of 5.0 × 10²⁰ s⁻¹. The main differences between the 10 and 15 MA cases are: (1) a slight increase of a few eV in the OMP ion temperature close to the plasma boundary for 15 MA; (2) a ∼0.1 × 10¹⁹ m⁻³ increase of both ion and electron OMP separatrix densities for 15 MA; and (3) a more pronounced increase in the inner target temperatures at 15 MA.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A.3. H-mode comparisons for ITER DT plasmas

In this section, we examine the impact of the energy-recycling modelling assumptions on a 15 MA/5.3 T DT H-mode with 33 MW of NB and 20 MW of ECRH heating. On top of a gas rate of 1.0 × 10²⁰ s⁻¹ we simulate pellet fuelling (50/50 DT, 74 mm³/4.6 × 10²¹ particles) to achieve a target density of ∼1.0 × 10²⁰ m⁻³, as required for Q = 10 operation in ITER. JINTRAC automatically adjusted the frequency of the pellets to reach and maintain this target through appropriate feedback. These simulations also included impurities, Be sputtered from the main vessel wall and Ne from seeding (1.0 × 10²⁰ s⁻¹) to reduce the divertor power fluxes, although the maximum power flux fluctuated between 5–30 MW m⁻² at the outer target and between 5–15 MW m⁻² at the inner target. Even if this case uses the latest assumptions for the ion energy source model, it should be noted that the power loads values quoted above are overestimated compared to the real values, since JINTRAC does not correct for the neutral reflected energy and this effect will decrease the power loads.

For the two different energy-recycling assumptions, the core time traces of the ion volume-average density and on-axis temperature are quite similar (figure A.5). This holds true for the time traces of the OMP separatrix temperatures as well. However, the OMP separatrix density with W-like energy recycling assumptions evolves with time whereas it remains rather stationary with Be-like energy recycling assumptions. This is likely due to a slight difference in the particle flux crossing the separatrix from the core into the SOL due to the higher energy of the recycled neutrals and increased core ionisation, which allows for a build-up of the OMP density in the W energy recycling assumptions case. This is also the cause of the different timings of the pellet injections as a pellet is not triggered unless the density goes below its target value by the feedback loop. Besides these small differences, the OMP SOL profiles are nearly identical (figure A.6). There are small variations (⩽10% at peak values) in the OT divertor profiles (figure A.7), while we can observe a more significant difference in the IT divertor ion temperature (∼30% at peak value) even if the other IT profiles differ by about 10% at peak values (figure A.8). It is important to remember that the plasmas modelled correspond to well attached conditions with minimal ionisation and recombination sources at the strike points (not shown). These plasma conditions are expected to show the largest differences in plasma behaviour due to the two different energy recycling assumptions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Contrary to the L-mode case in appendix A.2, where the divertor target temperatures increased when switching from a Be to a W divertor, in H-mode the temperatures dropped slightly. The difference occurs as for the same upstream conditions, the divertor temperatures tend to increase with more energetic recycling of neutrals, as illustrated in the L-mode case. However, in the H-mode case shown here, this is not so apparently obvious, as the time history (figure A.5) shows a more marked impact on the upstream conditions. This arises as the heating of the divertor region triggers higher impurity seeding to control the target heat load, and this cooling effect trumps the effect of the higher recycling energy.

A.4. Conclusions

B.1. Introduction

B.2. Grid generation and equilibrium

For the generation of the grid we used a 15 MA reference case, for which the equilibrium was provided by the ITER Organisation via a GEQDSK-File with a Ψ-map (poloidal magnetic flux) resolution of 65 × 129.

We did several improvements to the grid generation algorithm of the GRID2D code to provide an optimal starting point for the actual comparison between EDGE2D and SOLPS4.3. For instance, the algorithm in GRID2D to transform the metric coefficient H_ρ [28] from an orthogonal grid to a non-orthogonal grid, which can represent the tilted target geometry, had some inaccuracies. These inaccuracies lead to oscillations in H_ρ . We have now corrected the algorithm in GRID2D, which removes the oscillations in H_ρ and thus from the simulated pressure profiles.

Another upgrade we did to the GRID2D algorithm is that it now produces a smooth transition between the compression in the private flux region and the SOL region. This avoids the formation of strongly distorted cells during the grid generation process for ITER, where GRID2D needs a stronger grid compression in the SOL region, to comply with the shoulder region of the inner and outer targets, than the maximum compression applicable for the private flux region.

Moreover, as part of the comparison we also implemented the grid used by SOLPS4.3 with an increased grid resolution around the first few cells at the inner and outer target.

Figure B.1 shows an example of the impact of the new grid generation algorithm on the simulation results obtained by EDGE2D using the Mach number profile around the outer target. With the old grid (figure B.1, left) there is a strong distortion of the Mach number profiles around the shoulder of the outer divertor target, which is not present with the new grid (figure B.1, right).

Zoom In Zoom Out Reset image size

B.3. Vacuum vessel

B.4. EIRENE neutral model

B.5. Development on EDGE2D and important settings

B.6. Development on EDGE2D for He plasmas and important settings

B.7. Simulation results for 20 MW H density scan

For the following simulations and those in appendices B.8 and B.11, we use the settings described in appendix B.5 together with the reduced NIMBUS-like neutral model in EIRENE shown in appendix B.4. For the density scan with 20 MW input power we use the outer midplane (OMP) separatrix density obtained from the SOLPS4.3 simulations (table B.1) to control the EDGE2D OMP separatrix density to the same value, using a control puff.

Table B.1. Outer midplane separatrix density and reference run numbers for the SOLPS and EDGE2D simulations of the 20 MW H density scan.

nOMP,sep (m−3)	SOLPS simulation	EDGE2D simulation
0.573 × 1019	IDS #102292	10140
0.923 × 1019	IDS #102294	10141
1.234 × 1019	IDS #102296	10142
1.326 × 1019	IDS #102297	10143
1.767 × 1019	IDS #102299	10144
2.065 × 1019	IDS #102308	10145

For the 20 MW H density scans (figures B.2–B.5) show that both codes produce similar results apart from the two highest densities (n_OMP,sep = 1.767 × 10¹⁹ m⁻³ and n_OMP,sep = 2.065 × 10¹⁹ m⁻³). In the lower density range, EDGE2D and SOLPS4.3 tend to give comparable inner and outer target densities (figure B.3) and temperatures for both electrons and ions, with the biggest deviation for the inner target ion temperature (figure B.4). Moreover, the static target pressure is quite similar in both codes. For the power load on the targets (figure B.5; contribution from plasma species only), EDGE2D tends to give higher power loads than SOLPS4.3, with a better agreement on the outer target.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The disagreement between EDGE2D and SOLPS4.3 for the two highest densities was expected due to the differences in the neutral models used for EIRENE. As the two high density cases are already at and beyond the rollover, the additional processes such as molecular assisted recombination (MAR), which are only included in the full EIRENE (Kotov2008) neutral model, become significant.

B.8. Simulation results for 40 MW H density scan

The simulation cases for the H density scan at 40 MW are summarised in table B.2.

Table B.2. Outer midplane separatrix density and reference run numbers for the SOLPS and EDGE2D simulations of the 40 MW H density scan.

nOMP,sep (m−3)	SOLPS simulation	EDGE2D simulation
0.989 × 1019	IDS #102281	10146
1.808 × 1019	IDS #102277	10147
2.237 × 1019	IDS #102279	10148
2.799 × 1019	IDS #102301	10149

For the 40 MW H density scan the outer midplane profiles show very good agreement between EDGE2D and SOLPS4.3, apart from the outer midplane density outside of the separatrix that drops slightly faster in the EDGE2D simulations than in SOLPS4.3 (figure B.6).

Zoom In Zoom Out Reset image size

At the inner and outer targets, the profiles (figures B.7 and B.8) are fairly similar between EDGE2D and SOLPS4.3, though the peak of the density profile is shifted further away from the strike point in EDGE2D. There is a reasonable agreement between the codes for the power loads at the targets. However, at the inner target, EDGE2D still predicts a noticeable higher power load than SOLPS4.3 for the lowest separatrix density (figure B.9).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

B.9. Divertor neutral screening

In order to compare the divertor neutral screening, i.e. the fraction of the neutrals which escape from the divertor region, between the EDGE2D and SOLPS4.3 simulations, we compare 0-dimensional key parameters of the simulations, such as ionisation in different simulation regions (main SOL, inner and outer divertor), neutral influx to the core and recycling at targets and main chamber walls, for the 20 MW and 40 MW density scans.

Figure B.10 shows the recycling flux (particle outflux at the target returned as a neutral influx) at the inner (IT) and outer target (OT) versus the outer midplane separatrix density for the 20 MW density scan of EDGE2D and SOLPS4.3. Figure B.11 shows the same data for the 40 MW density scan is. Also here both codes are in the same recycling regime before the roll-over.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure B.12 shows the neutral influx over the separatrix into the core for the 20 MW and 40 MW H density scans. For the 40 MW density scan both codes predict nearly identical neutral core influxes. For the 20 MW density scan, the neutral core influxes agree in both codes below recycling flux roll-over. After the roll-over, EDGE2D gives a much higher neutral core influx, which is probably related to the different neutral model used in EIRENE by the EDGE2D code.

Zoom In Zoom Out Reset image size

For the 40 MW density scan, the divertor screening is very similar in EDGE2D and SOLPS4.3 (figure B.13). Only at the highest density, the divertor screening starts to differ between EDGE2D and SOLPS4.3. In the 20 MW density scan, the divertor screening is also similar in both codes before the recycling flux reaches the roll-over. After the roll-over the results from the codes diverge (note that the divertor screening (0.48 (SOLPS4.3)/0.17 (EDGE2D)) for the highest density is far outside of the plot).

Zoom In Zoom Out Reset image size

B.10. Simulation results for 40 MW at a single He OMP density

The simulation cases for a single OMP He separatrix density at 40 MW are summarised in table B.3.

Table B.3. Outer midplane separatrix density and reference run numbers for the SOLPS and EDGE2D simulations for a single He OMP density at 40 MW.

nOMP,sep (m−3)	SOLPS simulation	EDGE2D simulation
0.567 × 1019	IDS #112358	10656

A good match between EDGE2D and SOLPS4.3 for the OMP profiles is obtained for the electron and He⁺⁺ (figure B.14) densities and electron and ion temperatures. Compared to SOLPS4.3, EDGE2D predicts a lower He⁺ density outside the separatrix. The He⁺ is negligible inside the separatrix and then increases rapidly reaching a maximum 4 cm outside the separatrix of ∼5.4 × 10¹⁶ m⁻³ in SOLPS4.3 and ∼2.2 × 10¹⁶ m⁻³ in EDGE2D before falling off again towards the wall.

Zoom In Zoom Out Reset image size

At the outer target, He and electron densities (figure B.15), and electron temperature are very similar between the codes (figure B.16), with a minor discrepancy between the codes seen in the ion temperature (figure B.17).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

B.11. Conclusions of EDGE2D/SOLPS4.3 comparison