Multi-machine benchmark of the self-consistent 1D scrape-off layer model DIV1D from stagnation point to target with SOLPS-ITER

To enable quantitative interactions between the SOL and external domains, geometric properties are defined on the grid of DIV1D along the principle x direction. As a basis for the geometric description we define the following invariant for heat flux conservation:

$$\begin{equation} \lambda_\mathrm{q}\left(x\right) \cdot \text{sin}\left(\theta\right)\left(x\right) \cdot B_\mathrm{trans}\left(x\right) = C \quad \mathrm{[m]}, \end{equation} \tag{ 6 }$$

where $\lambda_\mathrm{q}$ is the poloidal SOL width and $\text{sin}(\theta)$ is the inclination angle given by the poloidal over the toroidal magnetic field. The expansion of the heat flux channel described by DIV1D ($\varepsilon_\mathrm{f}$ governing B in [7]) is split into a contribution from the magnetic field $B_\mathrm{field}$ and one that is used to mimic cross-field transport $B_\mathrm{trans}$, see examples of cross-field transport in [17]. In this way, when $\text{sin}(\theta)$ remains unchanged the transport field $B_\mathrm{trans}$ can still widen the poloidal SOL width $\lambda_\mathrm{q}$. The area $A_{\mathrm{extern},i}$ of cells ${i \in [1,N]}$ in contact with external domains along the poloidal direction, but not at the target/stagnation plane, is approximated as follows:

$$\begin{equation} A_{\mathrm{extern},i} = \text{sin}\left(\theta\right)\left(x_{i}\right) \cdot \Delta x_{\mathrm{cb},i} \cdot 2 \pi R\left(x_{i}\right) \quad \mathrm{[m^2]}, \end{equation} \tag{ 7 }$$

where the length of cells $\Delta x_\mathrm{cb}$ parallel to the magnetic field is projected poloidally with the pitch angle $\text{sin}(\theta)$ and multiplied with the circumference of the plasma $2\pi R$ at the cell location x_i . The cell volumes V_i follow from multiplication:

$$\begin{equation} V_{i} = A_{\mathrm{extern},i} \cdot \lambda_\mathrm{q}\left(x_{i}\right) \quad \mathrm{[m^3]}. \end{equation} \tag{ 8 }$$

The extended DIV1D geometry requires additional information on the inclination angle $\text{sin}(\theta)$, major radius R, and cross-field widening of the heat flux channel (mimicked in the transport field $B_\mathrm{trans}$). The transport field $B_\mathrm{trans}$ is set through the transport expansion coefficient $\varepsilon_{\perp}$ starting from the X-point location L_X up to the target with connection length L:

$$\begin{equation} B_\mathrm{trans}\left(x\right) = \left[1+\frac{\left(\varepsilon_{\mathrm{\perp}}-1\right) \left(x-L_{X}\right)}{L-L_{X}}\right]^{-1}. \end{equation} \tag{ 9 }$$

The magnetic field strength $B_\mathrm{field}$ and $\text{sin}(\theta)$ together define the magnetic equilibrium of DIV1D, where the magnetic field $B_\mathrm{field}$ can be given as a vector input to DIV1D or set through the magnetic flux expansion factor ε_B , by substitution in equation (9). The field B as used in the equations of DIV1D is obtained by multiplication of magnetic and transport effects: $B = B_\mathrm{field} \cdot B_\mathrm{trans}$. The external area $A_\mathrm{extern}$ and volumes V are calculated on cell centers, other quantities $\lambda_\mathrm{q}$, $\text{sin}(\theta)$, B, and R are also calculated on cell faces.

The boundary conditions at the upstream point are considerably more simple. The upstream boundary is now given by a stagnation point where all fluxes are set to zero. i.e. the plasma heat $q_{\|}$, velocity $v_{\|}$ and neutral density derivative $\frac{\partial n_\mathrm{n}}{\partial x}$ are set to zero. The target boundary conditions remain unchanged with respect to those in [7]. Although not used in this work, DIV1D also supports simulation of the full SOL from target to target using sheath boundary conditions at both ends of the domain and core-flux source terms in the core-sol.

The upstream boundary condition no longer allows any particle, momentum, or heat fluxes to enter the DIV1D domain. Instead the fluxes established at the X-point originate from the core and vacuum external to the SOL. The flux of external neutrals into the SOL is determined as:

$$\begin{equation} \begin{aligned} \Gamma^{\vdash}_{nn} & = \left(n_\mathrm{b}-n_\mathrm{n}\right)A_\mathrm{extern} v_\mathrm{ex} \quad \mathrm{[s^{-1}]}, \end{aligned} \end{equation} \tag{ 10 }$$

where the cross-tube neutral flux $\Gamma_\mathrm{nn}^{\vdash}$ follows from the difference between external gas density $n_\mathrm{b}$ and neutral SOL density n_n , multiplied with the external area $A_\mathrm{extern}$ and an exchange velocity $v_\mathrm{ex}$. The external gas density can be defined for both the core and divertor common flux region (CFR), as well as for the divertor private flux region (PFR). We note that $v_\mathrm{ex}$ replaces the neutral exchange time $\tau_\mathrm{ex}$ in [7]. The external neutral gas density $n_\mathrm{b}$ for the divertor and core region can be set independently to mimic density differences caused by baffles. Similarly different values can be specified for the CFR and PFR.

The cross-sol neutral flux enters the div-sol from both the PFR and CFR. As such the div-sol neutral source is determined by:

$$\begin{equation} S^\textrm{div}_\mathrm{nn} = \left(\Gamma^{\vdash,\mathrm{PFR}}_\mathrm{nn} + \Gamma^{\vdash,\mathrm{CFR}}_\mathrm{nn}\right) \cdot V^{-1} \quad \mathrm{[s^{-1}m^{-3}]}, \end{equation} \tag{ 11 }$$

where fluxes are doubled (with for simplicity equal neutral density in PFR and CFR) and divided by the volumes V. For the neutral flux into the core-sol it is noted that a large fraction $f_\mathrm{core}^{\,\textrm{ion}}$ can propagate into the core before being ionized, reducing the neutral source in the core-sol as follows:

$$\begin{equation} S^\textrm{core}_\mathrm{nn} = \Gamma^\vdash_\mathrm{nn}\left(1-f_\mathrm{core}^{\,\textrm{ion}}\right)\cdot V^{-1} \quad \mathrm{[s^{-1}m^{-3}]}. \end{equation} \tag{ 12 }$$

The plasma sources that originate from the core are calculated using the total core heat and ion fluxes $\Gamma_\mathrm{core}~\mathrm{[ s^{-1}]}$ and $Q_\mathrm{core}~\mathrm{[J\,s^{-1}]}$, respectively. These fluxes are multiplied with a normalized spatial profile. This core source spatial profile is calculated as:

$$\begin{equation} d\left(x\right) = \left(1 - \left(x/L_\mathrm{core}\right)^2\right)^{\alpha} \cdot A_\mathrm{extern} \quad \mathrm{[m^2]}. \end{equation} \tag{ 13 }$$

This is normalized to $\bar{d}$ where the integral $\int_{x = 0}^{L_\mathrm{core}} \bar{d}(x) \mathrm{d}x$ equals unity. The profile shaping parameter α can be used to shape the distribution and create peaked or flat profiles, enabling varying flux distributions over the seperatrix due to different (neoclassical or anomalous) transport phenomena. Multiplication with the area $A_\mathrm{extern}$ in the normalized distribution allows direct division with volumes to obtain volumetric sources for heat and particles (respectively $S_\mathrm{ene}^\textrm{core}$ and $S_\mathrm{n}^\textrm{core}$) as follows:

$$\begin{equation} \begin{aligned} S_\mathrm{ene}^\textrm{core} = Q_\mathrm{core} \cdot \bar{d} \cdot V^{-1} \quad \mathrm{[s^{-1}m^{-3}]}\\ S_\mathrm{n}^\textrm{core} = \Gamma_\mathrm{core} \cdot \bar{d} \cdot V^{-1} \quad \mathrm{[s^{-1}m^{-3}]} \end{aligned} \end{equation} \tag{ 14 }$$

where the flux distributions $\bar{d}$ are equal for particles and heat. It is noted that ionized particles leave the core-sol as they flow into the far SOL (far-sol). This far-sol sink is omitted in the DIV1D equations as we enforce all ion fluxes to be redirected into the divertor. This is contrary for neutral cross-channel transport as set by the core-ionization fraction $f_\mathrm{core}^\textrm{ion}$ and exchange velocity $v_\mathrm{ex}$. The sources are considered in the main equations of DIV1D according to their subscripts, e.g. $S_\mathrm{ene}^\textrm{core}$ is added to the right hand side of the energy balance. The superscripts div and core denote if the source takes effect in the divertor or core-sol.

The goal for DIV1D is to produce realistic SOL solutions as function of core fluxes and external neutral densities. In order to systematically test DIV1D, the following section investigates the main SOL heat flux channel based on 2D solutions of SOLPS-ITER.

In this work, 1D SOL equilibria always represent the main heat flux channel and not single flux tubes. Due to cross-field transport, solutions of single flux tubes are unlikely to reflect macroscopic divertor plasma behavior [7]. This is illustrated in figure 2 where the full width at half the maximum (FWHM) of the heat flux [7] denotes the main heat-flux channel on the B2.5 grid coordinates (i.e. along grid cells x and y in figure 1). In figure 2, a dashed dotted line indicates the peak heat flux location. Along y toward the target, it can be seen that the channel widens in x and that the peak heat flux drifts radially outwards to larger x. Both effects are not captured by a single flux tube, while they are by the main heat flux channel. In the following paragraphs we discuss how the heat flux channel mapping is extended into the core-sol.

Zoom In Zoom Out Reset image size

In the core-sol the solutions are mainly determined by the perpendicular fluxes from the core. For example, the stagnation point is a result of the distribution of perpendicular core fluxes. Near the stagnation point parallel fluxes even become multi-directional, e.g. in the poloidal plane the bulk particle flow close to the seperatrix is down and further outward the bulk flow is up. To reduce the impact of multi-directional parallel fluxes and perpendicular flux distributions, we select a fixed number of B2.5 grid cells in the core-sol to represent the main heat flux channel (see also figure 2). The number of radial grid cells is chosen to align with the FWHM of the heat flux distribution just below the X-point, to circumvent 2D effects around the X-point [7]. Similar to [7], the main heat flux channel is used to map 2D SOLPS-ITER solutions of the SOL to 1D profiles along the leg. Along the x-coordinate in B2.5, the average, minimum and maximum values of quantities in the heat flux channel are used to construct 1D profiles along the leg (y-coordinate) where the extreme values represent lumped value intervals. It is noted that the mean values might not conserve physical quantities inside the selected regions of interest, e.g. internal energy. This 1D mapped SOLPS-ITER representation with fixed core-sol width and FWHM divertor-sol is used to inform DIV1D simulations and compare DIV1D with SOLPS-ITER in the remainder of this paper.

In this section we use the SOLPS-ITER solution of TCV simulation 150 683 to investigate the global heat and particle balances in the core-sol. Particularly interesting as boundary conditions for DIV1D are the fluxes that enter the core-sol over the seperatrix and leave into the far-sol. In literature the far-sol is typically denoted by changes in decay lengths and types of transport [21]. In this work the distinction between near SOL (near-sol)/core-sol and far-sol is not based on fundamental theory such as presented in [21], but simply by the domain that is averaged to obtain 1D profiles of plasma quantities from SOLPS-ITER.

In figures 3 and 4 the respective heat and ion particle flux over the seperatrix (SEP inflow) into the core-sol are depicted together with the outflow into the far-sol (far-SOL outflow). The distance to the target covers the core-sol from stagnation point to X-point, from 18 m to 8 m respectively. It can be seen that a large fraction of the fluxes that flow into the near-sol also flow out to the far-sol. The fraction of ions that flow out to the far-sol is larger than the fraction of heat flowing to the the far-sol. This is likely because the decay length for heat is smaller than for density [21]. Additionally, large fluxes can be observed near the X-point. Finally it is noted that, despite the large outflow into the far-sol, the majority of heat flows into divertor through parallel transport. The following subsection uses the heat flux channel mapping to determine inputs/settings for DIV1D.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The goal for DIV1D is to obtain realistic target and upstream SOL plasma quantities as function of core-fluxes (e.g. from core codes) and external neutral reservoir densities. In this section we extract settings for DIV1D from 1D mapped SOLPS-ITER solutions given that the main heat flux channel has a fixed number of B2.5 grid cells that define the width of the core-sol.

• The heat and particle sources from the core—Q_core and $\Gamma_{\textrm{core}}$—are obtained from the integral of cross-channel fluxes into the main heat flux channel as displayed in figures 3 and 4. The fluxes into the far-sol are not considered because a majority of these particles is still transported into the divertor. Also the majority of the heat is found to be transported into the divertor.
• The geometric values for $\text{sin}(\theta)$ and R can be specified separately in the div-sol and core-sol. The values are obtained by the average of the mapped 1D SOLPS-ITER profiles in the respective domains. Alternatively single values are used or full arrays are interpolated on the DIV1D grid.
• The transport field $B_\mathrm{trans}$ is extracted from SOLPS-ITER solutions by comparing the area of the heat flux channel in front of the target and just below the X-point while correcting for the contribution of magnetic flux expansion.
• The poloidal SOL width $\lambda_\mathrm{q}$ is constructed by setting a value around the outer midplane and augmenting this with the transport field $B_\mathrm{trans}$ and magnetic pitch angle $\text{sin}(\theta)$ (i.e. following equation (6)).
• The core ionization fraction $f_\mathrm{core}^\textrm{ion}$ is used as a fitting parameter, and expected to be related to the ratio between the ionization mean free path and poloidal SOL width $\lambda_\mathrm{q}$.
• The neutral exchange velocity $v_\mathrm{ex}$ is a fit parameter that is expected to have the order of magnitude of the neutral thermal velocity.
• The external neutral reservoir density $n^\textrm{div, core}_\mathrm{b}$ is obtained by averaging over the outer B2.5 cells for the core-sol and div-sol. This is an important difference with respect to [7], where the median value of the neutral particle density inside the SOL was selected. Another important remark is that in this work an effective neutral density is taken as only the atomic density.

A detailed list of (expected) DIV1D settings extracted from SOLPS-ITER for this paper is available in appendix A. Important to note are the three main parameters that are used to fit DIV1D solutions: the SOL width λ_q ; the core ionization fraction $f_\mathrm{core}^\textrm{ion}$; and the neutral exchange velocity $v_\mathrm{ex}$. The following section utilizes DIV1D settings as derived from mapped SOLPS-ITER solutions to simulate the SOL in 1D on multiple devices.

For the TCV the DIV1D model with core-sol is evaluated on a set of SOLPS-ITER simulations from [6] that represent a density ramp (see figure 5 in [7] where the same set of simulations was used). These SOLPS-ITER simulations are extensively compared to experiments in [6]. The inputs for DIV1D that follow from SOLPS-ITER simulations are as follows: core-sol length $L_\mathrm{core} = 10\,\textrm{m}$; div-sol length $L_\mathrm{div} = 8\,\textrm{m}$; transport expansion $\varepsilon_{\perp} = 2.2$; a homogeneous carbon concentration $\xi_\mathrm{C} = 0.03$, and heat flux from the core $Q = 270\,\textrm{kW}$ distributed with α = 1 following equation (13). Geometric descriptions of pitch angle $\text{sin}(\theta)$, major radius R, magnetic field strength $B_\mathrm{field}$ are arrays taken from mapped SOLPS-ITER profiles along the SOL. Neutral densities in reservoirs adjacent to the core and divertor (PFR and CFR) depend only on the atomic density in the edge of the B2.5 grid in SOLPS-ITER simulations. Detailed values are listed in table A1.

A number of parameters are chosen to match DIV1D with 1D mapped SOLPS-ITER heat flux channel profiles. The core-sol width is chosen at $\lambda_\mathrm{q} = 0.0153\,\textrm{m}$ to match the accumulated parallel heat flux at the X-point. The core ionization fraction $f_\mathrm{core}^\textrm{ion} = 0.8$ was used to match the upstream density and the neutral exchange velocity $v_\mathrm{ex} = 5\,\textrm{km s}^{-1}$ was used to obtain agreement in the div-sol profiles. Target recycling is set as $f_\mathrm{rec} = 0.5$ to reduce the target neutral density. Even though the same SOLPS-ITER solutions were used as in [7], the input parameters for DIV1D were changed in two ways in this work: (1) the neutral exchange velocity was changed; (2) the neutral reservoir density is here taken at the boundary of the B2.5 grid where in [7] one takes the median of the atomic density inside the SOL.

Results on three simulations in a density ramp for high-recycling conditions are presented in figure 5, where DIV1D is depicted together with SOLPS-ITER. It can be seen that the profiles for heat flux $q_\|$, temperature T and electron density $n_\mathrm{e}$ of SOLPS-ITER and DIV1D are similar with DIV1D often laying inside the value intervals of the 1D mapped SOLPS-ITER solutions. The neutral atom density $n_\mathrm{D0}$ of DIV1D lies close to the value intervals from the X-point up to a few meters in front of the target, but is very different in the core-sol. Another large discrepancy is observed in the parallel velocity where DIV1D overestimates the velocity by up to 10 km s⁻¹. Also notable is that the temperature of ions and electrons in the core-sol is not the same in SOLPS-ITER and that DIV1D is consistently close to the electron temperature.

Zoom In Zoom Out Reset image size

Additionally, we compare the upstream electron density $n_\mathrm{e,omp}$ and temperature $T_\mathrm{omp}$ of DIV1D to that of SOLPS-ITER. Figure 6 depicts the upstream density as function of neutral density in the divertor CFR gas reservoir $n^\textrm{div}_\mathrm{b,CFR}$. It can be seen that the trend of DIV1D is similar to that of SOLPS-ITER although DIV1D underestimates the upstream density $n_\mathrm{e,omp}$ by roughly 20% for the neutral densities $n^\textrm{div}_\mathrm{b,{CFR}}$ above $9\cdot 10^{17}\,\mathrm{m^{-3}}$. In figure 7 the outer midplane (upstream) temperature $T_\mathrm{omp}$ versus upstream electron density $n_\mathrm{e,omp}$ is compared for 1D mapped SOLPS-ITER and DIV1D heat flux channel solutions. It can be seen that DIV1D temperatures are inside the value interval of electrons in 1D mapped SOLPS-ITER.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the AUG tokamak we calibrate the DIV1D model on three SOLPS-ITER simulations reported in [19]. These simulations represent divertor plasmas from the onset of detachment to complete detachment and are validated on experimental discharges #27 100 and #34 821. The inputs for DIV1D again follow from mapping SOLPS-ITER solutions, with detailed values given in table A2. The most important settings of DIV1D include a power of roughly $750\,\textrm{kW}$, no impurities, small flux expansion $\varepsilon_{B} = 1.1$ and very small transport based expansion $\varepsilon_{\perp} = 1.06$. Geometrical parameters are set as: connection length $L = 35\,\textrm{m}$, div-sol length $L_\mathrm{div} = 10\,\textrm{m}$, $\text{sin}(\theta) = 0.035$, and major radius $R = 1.8\,\textrm{m}$. Interesting to note is that both heat and particle fluxes over the seperatrix change between simulations and the core-sol widens. These varying parameters are listed in table 2 along with changes in the external neutral gas reservoir density in the PFR as extracted from SOLPS-ITER solutions.

Table 2. Important settings of DIV1D for the comparison for AUG in figure 8

MDS	$\Gamma_\mathrm{core}$	$Q_\mathrm{core}$	$n^\textrm{div,PFR}_\mathrm{b}$	$\lambda_\mathrm{q}$
182 202	0.6	703	1.1	0.05
182 203	0.9	719	4.6	0.08
182 204	1.3	835	7.3	0.15
Unit	$(10^{22}\,\textrm{s}^{-1})$	(kJ s−1)	$(10^{19}\,\textrm{m}^{-3})$	(m)

To calibrate DIV1D, we select the SOL width $\lambda_\mathrm{q}$ to match the X-point heat flux. To match the particle balance we chose the neutral exchange velocity $v_\mathrm{ex} = 0.56\,\textrm{km s}^{-1}$; target recycling $f_\mathrm{rec} = 0.5$; and core ionizing fraction $f_\mathrm{core}^\textrm{ion} = 0.6$. The varying inputs for the DIV1D simulations are listed in table 2. Note here that we vary the SOL width $\lambda_\mathrm{q}$ beyond the values found from 1D mapped SOLPS-ITER. Equally important is that we list the change in PFR neutral density for these SOLPS-ITER simulations as the neutral density in the PFR exceeds that in the CFR as a result of the vertical targets in these AUG experiments. In SOLPS-ITER this neutral balance is tuned to roughly match experiments by introducing artificial baffling and pumping to mimic the changing neutral transport characteristics during a discharge, e.g. because the wall transitions from a sink to a source [22].

Varying the setting for the midplane SOL width, results of three DIV1D simulations for AUG are depicted in figure 8 ranging from attached to detached divertor plasmas. The attached state can be recognized in simulation $\#$182 202 by the high target temperature (${\approx}20\,\textrm{eV}$) that relates to a high target velocity and low recycling. This low-recycling relates to the drop in electron density near the target. For details relating to the detachment state the reader is referred to [11, 23] and references therein. There are two distinguishable dashed blue lines for $n_\mathrm{b}$ in the divertor below the X-point where the high and low value respectively represent the PFR and CFR atomic neutral density as used in the DIV1D simulation. Compared to the mapped SOLPS-ITER solutions, it can be seen that DIV1D matches the electron temperature. Whereas the ion temperature is consistently higher in the core-sol. The parallel heat flux of DIV1D also closely follows SOLPS-ITER. The connection length of DIV1D was set constant across simulations while the connection length of 1D mapped heat flux channel profiles varied from 36 to 32 m. Consequently one can observe local minima in the heat flux around 33 m for #182 203 and #182 204. Finally, it can be seen that DIV1D overestimates the velocity around the X-point and that DIV1D overestimates the upstream electron density in #182 202. This could be due to drift effects, see e.g. [24], not considered in DIV1D [7].

Zoom In Zoom Out Reset image size

For the MAST-U tokamak, DIV1D is calibrated on SOLPS-ITER simulation #67 989 from [20]. This simulation represents an L-mode high power edge plasma for a double-null configuration with a super-X tightly baffled divertor geometry. Unfortunately it has not been compared to experimental data yet in [20]. Again, most inputs for DIV1D are set using the 1D mapped SOLPS-ITER solution of the main heat flux channel.

For the MAST-U simulation we use the following base settings: a power of $Q_\mathrm{core} = 1.250\,\textrm{kW}$, carbon impurity concentration $\textrm{f}_\mathrm{C} = 0.02$, flux expansion $\varepsilon_{B} = 2.6$ and transport expansion $\varepsilon_{\perp} = 3.0$. The geometry is defined by a connection length $L = 28.9\,\textrm{m}$, div-sol length $L_\mathrm{div} = 23.1\,\textrm{m}$, $\text{sin}(\theta) = 0.05$, and major radius $R = 0.97\,\textrm{m}$. The large change in major radius toward the target is expressed in the flux expansion $\varepsilon_{B} = 2.6$ and causes significant changes in the value of the field B in the DIV1D equations. It is also noted that the very detached equilibrium of MAST-U divertor plasmas obstruct the derivation of transport expansion $\varepsilon_{\perp}$ from 1D mapped SOLPS-ITER solutions, consequently this value is relatively uncertain compared to those found for TCV and AUG.

The process of calibrating DIV1D on SOLPS-ITER simulations of MAST-U is similar to that on TCV and AUG. Firstly, changing the SOL width $\lambda_{q} = 0.077$ to match the X-point parallel heat flux $q_{\|}$. Secondly, adapting the neutral exchange velocity $v_\mathrm{ex} = 1\,\textrm{km s}^{-1}$ and core ionization fraction $f_\mathrm{core}^\textrm{ion} = 0.9$ to adjust the particle balance. Target recycling is set to $f_\mathrm{rec} = 0.9$, but solutions are not very sensitive to this parameter in the considered simulation.

The resulting calibrated DIV1D simulation for MAST-U is depicted together with the corresponding 1D mapped SOLPS-ITER heat flux channel solution in figure 9. It can be seen that the temperature T and parallel heat flux $q_{\|}$ lie almost within the value intervals of the SOLPS-ITER solutions. On the other hand, the plasma density $n_\mathrm{e}$ is over-predicted between 26 and 20 m to the target while the parallel velocity $v_{\|}$ is under-estimated by DIV1D.

Zoom In Zoom Out Reset image size

A highlight for the MAST-U super-X divertor is the volumetric recombination region, which can be seen from 14 m to the target as a peak and decay in the plasma density $n_\mathrm{e}$. As the temperature T falls below 0.5 eV, recombination becomes a dominant sink in DIV1D causing the plasma density $n_\mathrm{e}$ to reduce two orders of magnitude—ending below $10^{18}\,\mathrm{m^{-3}}$ whereas the solutions of SOLPS-ITER ion densities remain above $1\cdot 10^{19}\,\mathrm{m^{-3}}$.

In previous work on validating DIV1D in [7] the domain was constrained to a region below the X-point. The reason was that the X-point plasma quantities in 1D mapped SOLPS-ITER profiles show large variations. The extension of DIV1D with the core-sol enables it to self-consistently determine the plasma density and parallel heat flux at the X-point. Important in this regard is that the representation of the core-sol (e.g. total volume, length, width) allows translation of the core fluxes and external neutral reservoir density into correct X-point plasma quantities (i.e. parallel heat flux and plasma density) and upstream plasma quantities for temperature and density. As such, when calibrating DIV1D on SOLPS-ITER we are required to alter the SOL width $\lambda_\mathrm{q}$ and core ionization fraction $f_\mathrm{core}^\textrm{ion}$ for a match in the X-point heat flux $q_\|$ and density $n_\mathrm{e}$, respectively. The translation of core fluxes into matching upstream parameters by DIV1D is demonstrated for TCV in figures 6 and 7. Similarly a match is achieved for the X-point heat flux $q_{\|}$ on TCV, AUG, and MAST-U whereas the X-point density $n_\mathrm{e}$ is not consistently matched (see figures 5, 8 and 9).

In calibrating DIV1D on SOLPS-ITER we use the external atomic densities while it would also be possible to use molecular densities. The molecular densities outside the SOL typically exceed atomic densities in the used SOLPS-ITER simulations. When the external molecule densities exceed the neutral density inside the SOL over the entire domain, this complicates calibrating DIV1D. Especially near the target, high external neutral densities inhibit neutrals in DIV1D from leaving the SOL to recirculate. As such we find high external densities can cause unrealistic neutral compression near the target. The addition of these molecular or high external density effects—to 1D models such as DIV1D—is not addressed in this paper.

After calibration on 1D mapped SOLPS-ITER solutions, the 1D simulations presented in this work represent realistic SOL solutions. Realistic since the calibrated DIV1D simulations on TCV, AUG, and MAST-U match the heat flux, temperature, and plasma density of 1D mapped SOLPS-ITER profiles within roughly 50% for plasma regions with temperatures above 1 eV. However, on MAST-U there are still large discrepancies for divertor plasmas at temperatures below 1 eV, possibly due to molecular effects which can have a large role in highly dissipative divertor plasmas [25]. Also the velocity $v_{\|}$ shows quite some discrepancies, overestimating it for TCV (see figure 5), underestimating it in MAST-U (see figure 9), whereas on AUG the velocity seems to match from X-point to target but not in the core-sol (see figure 8). Discrepancies in the parallel velocity might be addressed by investigating the particle and momentum balance in more detail.

For the calibration of DIV1D on SOLPS-ITER simulations for AUG there are two additional interesting observations. Firstly, it seems infeasible to use a single core-sol width in DIV1D when calibrating DIV1D on several AUG SOLPS-ITER simulations. This is also follows from the SOLPS-ITER simulations that feature varying fall-off lengths at the outer midplane. Extensive studies on fall-off lengths [26–29], including those of neutrals [30], might provide a way to constrain the core-sol width of DIV1D with a physics basis in the future. Secondly, for AUG it was required to consider the neutral density in the PFR to calibrate DIV1D on SOLPS-ITER. In the PFR the neutral density is significantly lower than in the CFR. The large difference in neutral density between CFR and PFR is caused by the divertor geometry and the placement of valves and pumps, something that has received attention in the SOLPS-ITER modeling [22]. These effects are not modeled self-consistently in the present work but are enforced through the external neutral densities, where molecules are explicitly omitted.

In extending the DIV1D model to the core-sol, several simplifying assumptions are made. In particular considering only a single SOL width in the core-sol for all quantities is not in-line with SOL literature, see e.g. [21]. As a consequence, the DIV1D model must be calibrated on existing solutions of the SOL to obtain matching core-sol solutions and therefore to obtain relevant upstream conditions for the div-sol solution. In this work the main calibration parameters are the core-sol width $\lambda_\mathrm{q}$, neutral exchange velocity $v_\mathrm{ex}$ and fraction of neutrals that ionize in the core $f_\mathrm{core}^\textrm{ion}$. A very clear limitation for now is that one fit of DIV1D only aligns with SOLPS-ITER in a finite region, not ranging from fully attached to fully detached states.

The simple geometry of DIV1D, also makes it fairly simple to qualitatively evaluate different divertor configurations. This is illustrated by the possibility to calibrate on plasmas in the super-X divertor of MAST-U. Similarly future work could investigate an X-point target configuration that features significant poloidal flux expansion [3]. From the DIV1D equations it can be seen that poloidal expansion increases the SOL width and decreases the opportunity of recycled target neutrals to escape the SOL, see also [22]. The effects of poloidal flux expansion may therefore be attributed to more than the effect on the equations of the charged species as typically used in analytic modeling for alternative divertor configurations (e.g. in [31]). It is noted, however, that the geometry of DIV1D is strictly limited to configurations where there are clearly distinct main heat flux channels.

The extension of DIV1D omitted separate treatment of ion temperatures. The difference in upstream ion and electron temperatures is typical for low collision edge plasmas and expected to be present in reactor scale devices [32]. Interesting for reactors could then also be the impact of kinetic effects on SOL solutions [33]. However, we expect the omission of ion temperatures to have limited impact on the foreseen use of DIV1D for optimization and control. In such applications it is expected to be sufficient to assume a certain ratio in upstream temperature between ions and electrons.

The range of devices and scenario's that DIV1D reproduced in this work provides a basis to toward mimicking 1D mapped SOLPS-ITER heat flux channel solutions on JET [34], DEMO [35], ITER [36], STEP [37], and SPARC [38], as well as more elaborate analysis on TCV, AUG, and MAST-U. The computational speed of DIV1D (or another 1D code matching SOLPS-ITER) make it attractive for studies in optimization and preliminary assessment of tokamak exhaust designs. In this workflow, one is required to first provide an equilibrium with SOLPS-ITER and then explore its surroundings with DIV1D. We do stress that most engineering guarantees are still provided by SOLPS-ITER.

In terms of general 1D modeling efforts, many of the recommendations in [7] still hold, including drifts [24, 39], molecules [10], and better impurity emission rates [40]. Additionally, one could estimate peak fluxes by superimposing skewed Gaussians or distributions from EMC3-Lite [41]. Also, the addition of an ion energy equation could be explored to improve the description in the core-sol. Finally, one could investigate the connection of 1D models (such as DIV1D) with surrounding domains. To a great extent that has been the goal of this work, to extend DIV1D with the core-sol for a connection with the core and to investigate cross-field neutral transport in order to connect DIV1D with external neutral reservoirs. For the connection to the core one should consider that the SOL width may no longer be a free parameter because the power fall-off length is likely governed by core plasma conditions [27].

Finally, as DIV1D is being developed specifically to aid exhaust control efforts on fusion devices, it should reproduce exhaust dynamics. To investigate these dynamics, DIV1D can be perturbed around stationary reference solutions in a similar way that system identification experiments are performed [42]. Considering these dynamics, the equilibration timescales of the SOL observed in DIV1D are on the order of a few milliseconds [43]. On the other hand, experimental observations of the exhaust dynamics in response to gas valves are on the order of confinement timescales: in AUG timescales for reattachment are on the order of 100 ms [44] whereas for TCV typical exhaust controllers operate below 10 Hz [2] (i.e. around timescales of 100 ms). As such, the relatively fast timescale of DIV1D on its own is unlikely to explain the relative slow timescales observed in AUG and acted upon in TCV. We believe the slow experimental timescales to involve the core and external neutral reservoirs, which can be coupled to DIV1D given the extension with core SOL and geometry as presented in this paper. Such dynamic investigations are part of future research.

The equations remain discretized on a nonequidistant grid as employed also in the SD1D code [14]. For N grid cells the boundaries $x_{\textrm{cb},i}$ counting from i = 0 at the X-point on the left to i = N at the target on the right are given by [14]

$$\begin{equation} x_{\textrm{cb},i} = L \left( {\left(2 - \delta_{x,\textrm{min}}\right) i \over N} - {\left(1 - \delta_{x,\textrm{min}}\right) i^2 \over N^2} \right) \end{equation} \tag{ B.1 }$$

where L is the total distance from X-point to the target and $\delta_{x,\textrm{min}}$ is a parameter that sets the ratio between the smallest grid cell at the target and the average grid cell size. The cell centers are defined as

$$\begin{equation} x_{i} = {x_{\mathrm{cb},i} + x_{\mathrm{cb},i-1} \over 2}, \quad\quad\hbox{for}~{i = 1 \cdots N}. \end{equation} \tag{ B.2 }$$

The widths of the grid cells are given by

$$\begin{equation} \Delta x_{\mathrm{cb},i} = x_{\mathrm{cb},i} - x_{\mathrm{cb},i-1}, \quad\quad\hbox{for}~{i = 1,\ldots ,N}. \end{equation} \tag{ B.3 }$$

Similarly, the distance between cell centers defines

$$\begin{equation} \Delta x_{i} = x_{i+1} - x_{i}, \quad\quad\hbox{for}~{i = 1,\ldots ,N-1}. \end{equation} \tag{ B.4 }$$

Following the cell-centered finite volume method, variables are calculated on the cell centers and fluxes are calculated on the cell boundaries.

The primary variables that are evolved in the code are the plasma density n, the plasma momentum $P \equiv n m v_\parallel$, the total internal energy $E \equiv 3 n k T$, and the neutral density $n_\mathrm{n}$. For evaluation of the ODE solver the normalized unknown variables are stacked in a vector Y of length $4{N}$ as

$$\begin{align} Y_{i} \,&\equiv {\displaystyle n_i \over n_\mathrm{norm} } &\quad\quad\quad\hbox{for}~{i = 1,\ldots ,N}, \nonumber \\ Y_{N + i} \,&\equiv {\displaystyle P_{i} \over n_\mathrm{norm} m c_\mathrm{norm} } &\quad\quad\hbox{for}~{i = 1 \cdots N}, \nonumber\\ Y_{2N + i} \,&\equiv {\displaystyle E_{i} \over 3 n_\mathrm{norm} k T_\mathrm{norm} } &\quad\quad\hbox{for}~{i = 1 \cdots N}, \nonumber\\ Y_{3N + i} \,&\equiv {\displaystyle n_\mathrm{n} \over n_\mathrm{norm} } &\quad\quad\hbox{for}~{i = 1 \cdots N}, \end{align} \tag{ B.5 }$$

where $c_\mathrm{norm} = \sqrt{2 k T_\mathrm{norm}/m}$ is the sound speed at the normalizing temperature $T_\mathrm{norm}$. The normalized density $n_\mathrm{norm}$ can be a single value or an array based on the density n_i .

The code now features second order accurate interpolation to the cell faces/boundaries and first order accurate extrapolation at the target to replace mirror cells. This is detailed in equations (B.6), where quantities u represent velocity, temperature, and density. Values on the cell boundaries are obtained from interpolation using the adjacent cell centers. For the cell boundaries at the edge of the domain, the velocity, temperature, and densities are linearly extrapolated for simplicity. Note that there are N + 1 cell boundaries, starting from zero.

$$\begin{align} u_{i} = \,& \frac{1}{2} \left(u_{i}\Delta x_{\mathrm{cb},i+1} + u_{i+1}\Delta x_{\mathrm{cb},i} \right) / \left(2\Delta x_{i}\right) \nonumber \\ \,& \hbox{for}~{i = 1 \cdots N-1}, \nonumber\\ u_0 = \,& u_1 - \left(u_2 - u_1\right) \frac{\Delta x_1}{0.5 \Delta x_\mathrm{cb,1}} \nonumber \\ u_{N} = \,& u_{N} - \left(u_{N} - u_{N-1}\right) \frac{\Delta x_{N-1}}{0.5 \Delta x_{\mathrm{cb},N}}. \end{align} \tag{ B.6 }$$

The advected part of the fluxes is calculated similar to the numerical scheme as used in SD1D [14] with a combination of slope and flux limiters with a lax flux damping discontinuities.

B.3.1. Slope limiter.

The slope limiter determines a limited slope of the solution inside a cell based on slopes over cell boundaries to its neighbors: $s_{i} = \frac{u_{i+1}-u_{i}}{\Delta x_{i}}$. The limited slopes σ_i of the cells are obtained, compliant to the non-equidistant grid as follows [47]:

$$\begin{align} \sigma_{i}& = \frac{u_{i+1}-u_{i}}{\Delta x_{i}} \phi\left(A_{i},B_{i},\theta_{i}\right) , \end{align} \tag{ B.7 }$$

$$\begin{align} A_{i} &= \frac{\Delta x_{\mathrm{cb},i-1} + \Delta x_{\mathrm{cb},i}}{\Delta x_{\mathrm{cb},i}+ \Delta x_{\mathrm{cb},i+1}}, \end{align} \tag{ B.8 }$$

$$\begin{align} B_{i} &= \frac{2 \Delta x_{\mathrm{cb},i}}{\Delta x_{\mathrm{cb},i}+ \Delta x_{\mathrm{cb},i+1}}, \end{align} \tag{ B.9 }$$

$$\begin{align} \theta_{i} & = \frac{u_{i}-u_{i-1}}{u_{i+1}-u_{i}} \nonumber\\ & \quad \hbox{for}~{i = 2 \cdots N-1}. \end{align} \tag{ B.10 }$$

Where φ is the following enhanced MinMod slope limiter [47]:

$$\begin{align} \phi\left(A_{i},B_{i},\theta_{i}\right) = B_{i}/A_{i} \cdot \text{min}\left(\theta_{i}, A_{i}\right). \end{align} \tag{ B.11 }$$

For clarification, figure B1 compares behavior of normal slope limiters and enhanced slope limiters on differences and slopes for a non-equidistant grid. It is noted that numerical schemes with both normal and enhanced slope limiters are total variation diminishing (TVD) stable when applied to non-equidistant grids—e.g. in [7, 14]—as both satisfy the following sufficient condition [47, 48]:

$$\begin{equation} 0 \unicode{x2A7D} \phi_{i}\left(\theta\right) \unicode{x2A7D} 2, \quad 0 \unicode{x2A7D} \frac{\phi_{i}\left(\theta\right)}{\theta} \unicode{x2A7D} 2, \quad \forall {i}, \theta. \end{equation} \tag{ B.12 }$$

The boundary cells lack a neighbor such that the standard slope limiter routine cannot be applied. For simplicity, we use the slopes to neighbors $\sigma_1 = s_1, \sigma_{N} = s_{N-1}$. Alternatively one could use extrapolated quantities for cells mirrored behind the boundaries. Using the enhanced slope limiter φ, the advected quantities f are extrapolated to left (L) and right (R) boundary as follows:

$$\begin{equation} f_{i}^{\,\textrm{L,R}} = f_{i} \pm \sigma_{i} \cdot \frac{1}{2}\Delta x_{\mathrm{cb},i}. \end{equation} \tag{ B.13 }$$

Zoom In Zoom Out Reset image size

B.3.2. Flux limiter and Lax flux.

The flux limiter and Lax flux remain the same as in Dudson et al [14]. The advected fluxes on the cell boundaries $F_{\mathrm{cb},i}$ are calculated as:

$$\begin{align} F_{\mathrm{cb},i}^{\,\textrm{H}} =\,& \frac{v_{\mathrm{cb},i}}{2} \left(f_{i}^{\,\textrm{R}} + f_{i+1}^{\,\textrm{L}} \right) + \underbrace{\frac{c_{\mathrm{s,cb},i}}{2} \left(f_{i}^{\,\textrm{R}} - f_{i+1}^{\,\textrm{L}} \right)}_\mathrm{Lax ~flux} \end{align} \tag{ B.14 }$$

$$\begin{align} F_{\mathrm{cb},i}^{\textrm{L}} = \,& v_{\mathrm{cb},i\,}f_{i}^{\,\textrm{R}}\quad \quad \hbox{for}~v_{{\mathrm{cb},i}}> c_{\textrm{s,cb,i}} \end{align} \tag{ B.15 }$$

$$\begin{align} F_{\mathrm{cb},i}^{\,\textrm{L}} = \,& v_{\mathrm{cb},i\,}f_{i+1}^{\,\textrm{L}}\quad \quad \hbox{for}~v_{{\mathrm{cb},i}}< -c_{\textrm{s,cb},i} \end{align} \tag{ B.16 }$$

where in addition to the slope limiter φ, a flux limiter $\varphi_{\mathrm{cb},i}$ is applied in the form of:

$$\begin{equation} F_{\mathrm{cb},i} = F_{\mathrm{cb},i}^\textrm{L}+\varphi_{\mathrm{cb},i}\left(F_{\mathrm{cb},i}^\textrm{H}-F_{\mathrm{cb},i}^\textrm{L}\right) \end{equation} \tag{ B.17 }$$

which switches to a lower (superscript L) order upwind scheme for supersonic velocities, i.e. $|v_{\mathrm{cb},i}| \gt c_{\mathrm{s,cb},i}$. For subsonic flow, the higher (superscript H) order central scheme is used and a Lax flux is added to damp discontinuities allowed in the slope of the solution by the slope limiters. Note thus that there are two types of limiters, a slope limiter φ_i different for each cell due to the non-equidistant grid and flux limiter $\varphi_{\mathrm{cb},i}$ switching between integration methods. Pressure gradients are discretized using a centered scheme for non-equidistant grids. For more details see [14, 47] and references therein.

In terms of equations this entails the following finite differences for the density, momentum, energy and neutral equations with i denoting the cell with cb, i and cb, i − 1, the upper and the lower cell boundary respectively:

$$\begin{align} \dot{Y}_{i} \equiv& - B_\mathrm{i} \left( \frac{\Gamma_{n,{i}}}{B_{\mathrm{cb},i}} -\frac{\Gamma_{n,{i-1}}}{B_{\mathrm{cb},i-1}} \right)/ \Delta x_{\mathrm{cb},i} + S_{\mathrm{n},i}, \end{align} \tag{ B.18 }$$

$$\begin{align} \dot{Y}_{N + i} \equiv& - B_{i} \left( \frac{\Gamma_{\mathrm{mom},i}}{B_{\mathrm{cb},i}} -\frac{\Gamma_{\mathrm{mom},i-1}}{B_{\mathrm{cb},i-1}} \right)/\Delta x_{\mathrm{cb},i} \nonumber \\ & - \frac{\left(p_i-p_{i-1}\right)}{2\Delta x_{i-1}} - \frac{\left(p_{i+1}-p_{i}\right)}{2\Delta x_{i}} + S_{\textrm{mom},i}, \end{align} \tag{ B.19 }$$

$$\begin{align} \dot{Y}_\mathrm{2N + i} \equiv& - B_{i} \left( \frac{q_{\|,{i}}}{B_{\mathrm{cb},i}} -\frac{q_{\|,{i}-1}}{B_{\mathrm{cb},i-1}} \right)/\Delta x_{\mathrm{cb},i}, \nonumber \\ &+ \frac{v_{\|,i}\left(p_{i}-p_{i-1}\right)}{2\Delta x_{i-1}} + \frac{v_{\|,i}\left(p_{i+1}-p_{i}\right)}{2\Delta x_{i}} + S_{\textrm{ene},i}, \end{align} \tag{ B.20 }$$

$$\begin{align} \dot{Y}_{3N + i} \equiv& \left( D\left(n_{\mathrm{cb},i},T_{\mathrm{cb},i}\right)\frac{n_{\mathrm{n},i+1}-n_{\mathrm{n},i}}{\Delta x_{i}\Delta x_{\mathrm{cb},i}} \right. \nonumber \\ &\left. - D\left(n_{\mathrm{cb},i-1},T_{\mathrm{cb},i-1}\right)\frac{n_{\mathrm{n},i}-n_{\mathrm{n}, i-1}}{\Delta x_{i-1}\Delta x_{\mathrm{cb},i}} \right) + S_{\mathrm{neutral},i}. \end{align} \tag{ B.21 }$$

At the boundaries the pressure gradients are calculated by forward or backward differences whereas the fluxes are determined by the boundary conditions.

At the target(s), the fluxes are calculated using the Bohm condition $v_\mathrm{target} \unicode{x2A7E} c_\mathrm{s,target }$ with $c_\mathrm{s} = \sqrt{2kT/m}$, i.e. $v_\mathrm{Bohm} = \text{max}(|v_\mathrm{target}|, c_\mathrm{s,target }) \text{sign}(v_\mathrm{target })$ for isothermal flow toward the target. The target velocity and temperature are extrapolated from previous cell centers using equation (B.6). Accordingly, fluxes at the target(s) are:

$$\begin{align} \Gamma_{n,\textrm{t}} & = v_\mathrm{Bohm}\cdot n_\mathrm{t}/B_\mathrm{t} \end{align} \tag{ B.22 }$$

$$\begin{align} \Gamma_\mathrm{mom,t } & = v_\mathrm{Bohm}\cdot n_\mathrm{t}m v_\mathrm{Bohm} / B_\mathrm{t} \end{align} \tag{ B.23 }$$

$$\begin{align} q_\mathrm{t} & =\text{sign}\left(v_\mathrm{t }\right)c_\mathrm{s,t}\cdot \gamma n_\mathrm{t} e T_\mathrm{t } / B_\mathrm{t}. \end{align} \tag{ B.24 }$$

The fluxes at the stagnation point are set to zero for a zero gradient boundary for the solved quantities.

In the implementation discussed so far there are two issues to be noted. Firstly, the Bohm criterion $v_\mathrm{target}\unicode{x2A7E} c_\mathrm{s,target}$ is not always satisfied because the velocity in the last cell in front of the target drops for an unknown reason. Secondly, for high resolution grids, numerical oscillations appear. To combat these issues, a numerical viscosity is implemented to penalize second order derivatives in the velocity profile and deviations from the Bohm criterion.