The role of particle, energy and momentum losses in 1D simulations of divertor detachment

TL;DR

This paper introduces a new 1D divertor plasma simulation code, SD1D, to analyze the roles of recombination, radiation, and momentum exchange in plasma detachment, revealing that radiative power losses are crucial for detachment and identifying thresholds for plasma state transitions.

Contribution

The study develops a novel analytical framework and simulation approach to understand the complex physics of divertor detachment, emphasizing the importance of radiative power losses and momentum exchange.

Findings

Radiative power losses are essential for reducing target ion flux.

Recombination plays a minor role at flux rollover but becomes significant at very low temperatures.

Detachment thresholds depend on energy cost gradients and pressure-to-power flux ratios.

Abstract

A new 1D divertor plasma code, SD1D, has been used to examine the role of recombination, radiation, and momentum exchange in detachment. Neither momentum or power losses by themselves are found to be sufficient to produce a reduction in target ion flux in detachment (flux rollover); radiative power losses are required to a) limit and reduce the ionization source and b) access low-target temperature, T_target, conditions for volumetric momentum losses. Recombination is found to play little role at flux rollover, but as T_target drops to temperatures around 1eV, it becomes a strong ion sink. In the case where radiative losses are dominated by hydrogen, the detachment threshold is identified as a minimum gradient of the energy cost per ionisation with respect to T_target. This is also linked to thresholds in T_target and in the ratio of upstream pressure to power flux. A system of determining the detached condition is developed such that the divertor solution at a given T_target (or lack of one) is determined by the simultaneous solution of two equations for target ion current - one dependent on power losses and the other on momentum. Depending on the detailed momentum and power loss dependence on temperature there are regions of T_target where there is no solution and the plasma 'jumps' from high to low T_target states. The novel analysis methods developed here provide an intuitive way to understand complex detachment phenomena, and can potentially be used to predict how changes in the seeding impurity used or recycling aspects of the divertor can be utilised to modify the development of detachment.

Equations

$$\frac{\partial n}{\partial t}$$

$$\frac{\partial }{\partial t}\left(\frac{3}{2}p\right)$$

$$\frac{\partial }{\partial t}\left(m_{i}n{v}_{||}\right)$$

$${\mathbf{q}}_e$$

$$v_n=v_{||n}-{\left(\frac{B_{\varphi }}{B_{\theta }}\right)}^2\frac{{\partial }_{||}{p}_n}{\nu },$$

$$\Gamma \propto n_{target}\sqrt{T_{target}}\propto p_{target}/\sqrt{T_{target}}.$$

$$P_{in}=P_{imp}+\underset{P_{recl}}{\underbrace{P_{iz}+P_{target}}}=P_{imp}+\left(E_{iz}+\gamma T_{target}\right){\Gamma }^{pow},$$

$${\Gamma }^{pow}=\left(P_{in}-P_{imp}\right)/\left(E_{iz}+\gamma T_{target}\right).$$

$${\Gamma }^{mom}\propto n_{target}\sqrt{T_{target}}\propto p_{target}/\sqrt{T_{target}},$$

$$2p_{target}/p_{up}=1-f_{mom}=0.9{\left(1-\exp \left(-T_{target}/2.1\right)\right)}^{2.9}.$$

$$1-f_{mom}$$

$$\alpha$$

$${\Gamma }^{mom}=p_{up}\underset{F\left(T_{target}\right)}{\underbrace{\frac{\left(1-f_{mom}\right)}{\sqrt{8m_i{T}_{target}}}}}.$$

$$\frac{d{\Gamma }^{pow}}{dT}=-\frac{1}{E_{iz}+\gamma T_{target}}\frac{d{P}_{imp}}{d{T}_{target}}-\frac{P_{in}-P_{imp}}{{\left(E_{iz}+\gamma T_{target}\right)}^2}\left(\frac{d{E}_{iz}}{d{T}_{target}}+\gamma \right).$$

$${\frac{\partial E_{iz}}{\partial T_{target}}|}_{critical}<-\gamma .$$

$$p_{up}/P_{recl}$$

$$=$$

Full text

A system of determining the detached condition is developed such that the divertor solution at a given $T_{target}$ (or lack of one) is determined by the simultaneous solution of two equations for target ion current – one dependent on power losses and the other on momentum. Depending on the detailed momentum and power loss dependence on temperature there are regions of $T_{target}$ where there is no solution and the plasma ‘jumps’ from high to low $T_{target}$ states. The novel analysis methods developed here provide an intuitive way to understand complex detachment phenomena, and can potentially be used to predict how changes in the seeding impurity used or recycling aspects of the divertor can be utilised to modify the development of detachment.

††: Plasma Phys. Control. Fusion

Keywords: tokamak, detachment, SD1D

1 Introduction

It has long been recognised that divertor plasma detachment will be required in future fusion devices such as ITER and DEMO, in order to keep divertor target heat loads below technological limits (e.g. [1, 2, 3]). Modelling of divertor detachment in magnetic confinement fusion devices is often done using 2D models [4, 5, 6], but simplified analytic [7, 8] and 1D computational models [9, 10, 11, 12, 13] can provide insight into the underlying processes, and provide guidance for optimisation of future devices. Here the SD1D model is presented (section 2), which has been developed using BOUT++ [14, 15] to study detachment dynamics. It is a time-dependent code, which enables the study of the detachment process, feedback control of detachment, and the response to plasma transients such as ELMs, in addition to steady-state solutions.

Before applying this model to time-dependent problems, we first use SD1D to understand the roles of particle, power and momentum loss mechanisms involved in detached steady state solutions. The importance of power loss vs momentum loss to the detachment process has been debated in the literature [7, 8, 16] and studied experimentally [17, 18]. In addition, simplified models have been constructed which emphasise power losses to predict the detachment threshold [19, 20, 21, 22] and sensitivity of the detached region extent/location to external controls [21].

All simulations are carried out in MAST-Upgrade like geometry, with a parallel heat flux of $50$MW/m2 at the X-point, $30$m connection length comprising $10$m above the X-point and $20$m from X-point to target. These are typical of expected conditions in the first phase of MAST-Upgrade operation [23]. The effect of gradients in the total magnetic field (total flux expansion [21, 24]) is included, with an area expansion factor of 2 (ratio of the total field at the X-point to that at the target) between X-point and target in all cases shown here.

2 The SD1D model

A 1D time dependent fluid model [7, 11, 12] is solved for the plasma density $n$, parallel momentum $nv_{||}$ and static pressure $p=2enT$, assuming equal isotropic ion and electron temperatures $T=T_{e}=T_{i}$.

[TABLE]

where $\partial_{||}\equiv\mathbf{b}\cdot\nabla$. As discussed in [25], the use of isotropic pressure likely overestimates the magnetic mirror effect in low collisionality regimes. It is also known experimentally that in general $T_{i}>T_{e}$ [26, 27]. More sophisticated models retaining both parallel and perpendicular ion pressures have been developed [10, 28], together with separate ion and electron temperatures. Those models used a simpler neutral gas model than is employed here, discussed below. Here we focus mainly on high collisionality regimes, and leave removing these limitations to future work.

An external source of power $S_{E}$ injects energy at a constant rate into a volume above the X-point, in this case the first $10$m of the domain. The external source of particles $S_{n}$ is varied using a proportional-integral (PI) feedback controller to achieve a specified upstream plasma density. Coupling to neutrals occurs through particle sources and sinks (ionisation and recombination) represented by $S$; energy exchange $E$; radiation $R$ due to hydrogen excitation and impurity radiation; and friction forces $F$ due to ionisation, recombination and charge exchange. These hydrogenic rates are calculated using semi-analytic approximations [13, 29]. Unless otherwise stated, in these simulations a 1% fixed fraction carbon impurity model is used, based on coronal equilibrium and calculated using ADAS data [30]. In all results shown here where both are included, hydrogenic radiation exceeds the carbon impurity radiation.

A similar set of three equations is evolved for the neutral fluid density $n_{n}$, parallel momentum $n_{n}v_{||n}$ and pressure $p_{n}=en_{n}T_{n}$. The neutrals are not confined by the magnetic field, so transport of neutrals across the magnetic field can provide a way for neutrals to migrate upstream. To mimic this process in a 1D model, the effective parallel velocity is given by the sum of a parallel flow and parallel projection of a perpendicular diffusion:

[TABLE]

where the collision frequency $\nu$ includes charge exchange, ionisation, and neutral-neutral collisions. The factor $\left(\frac{B_{\phi}}{B_{\theta}}\right)^{2}$ is set to $10$ in all simulation results shown here.

One of the assumptions made in going to a one-dimensional model concerns the handling of momentum losses. In particular, fast charge-exchanged neutrals can leave the thin SOL, transferring their momentum directly to the walls of the device without interacting with the rest of the neutral gas. In most of the simulations shown here charged exchanged neutrals do not escape; momentum is conserved, so that the total pressure (plasma + neutrals) is constant.

3 Detachment development

The definition of detachment varies amongst publications. In this paper our focus is on the particle flux to the divertor target, $\Gamma$, which goes through at least two stages during detachment in the context of the simple equation for sheath ion flux:

[TABLE]

During the attached phase $\Gamma$ rises $\propto T_{target}^{-1/2}$ as the target pressure $p_{target}$ is constant. During the first stage of detachment the rise of $\Gamma$ slows through a combination of ionization source loss and $p_{target}$ drop. We define the second stage, and the focus of most of this study, to be the point in any scan (e.g. upstream density, impurity seeding, SOL power) where $\Gamma$ ‘rolls over’; in other words the rate of change of $\Gamma$ becomes negative with respect to what is being varied. In this paper there are only density scans, although the carbon fraction in some cases is held constant leading to the carbon density rising.

The effects of a reference upstream density scan on target ion flux $\Gamma$ are shown in figure 1, and this scan will serve as a benchmark for other cases shown in this paper. We find that target ion flux rollover for the case labelled “Carbon + Hydrogen” occurs at an upstream density of around $1.89\times 10^{19}$m*-3* for fixed input power flux of $50$MW/m2. The “Carbon + Hydrogen” cases include radiative power losses from both hydrogen excitation radiation and 1% carbon impurity.

Also shown in figure 1 are results for a case where the cross-field neutral diffusion is doubled (equation 1b, “$\left(B_{\phi}/B_{\theta}\right)^{2}=20$”), and a case in which charge exchanged neutrals remove momentum from the simulation instead of transferring it to the neutral fluid (“CX neutrals escape”). Variation in the cross-field diffusion (due to changes in field-line pitch in equation 1b) has a modest effect on the results, but modification to the total (plasma + neutral) momentum loss has a large effect. This effect of neutrals escaping the SOL is examined further in section 5, where it is shown to be important in determining the total momentum loss fraction $f_{mom}$ (figure 4). This neutral loss process will depend on the device and configuration, so an improved model for neutral momentum loss is likely required for quantitative agreement with experiment.

In order to understand the behaviour of this model, and the underlying physical mechanisms, we here examine the plasma-neutral processes occuring during detachment. Of particular interest is the role of particle and energy loss (section 4), and the role of momentum loss (section 5) on the flux rollover at detachment.

SD1D has been benchmarked against the modified two-point model [31, 24], with good agreement for upstream densities below the flux rollover. As the momentum and power loss fractions ($f_{mom}$ and $f_{pow}$) increase, the result becomes sensitive to small errors in these quantities, so the agreement deteriorates: relative error is $3\times 10^{-5}$ at $n_{up}=1.8\times 10^{19}$m*-3* ($f_{mom}=0.22$, $f_{pow}=0.25$); $2$% after flux rollover at $n_{up}=1.9\times 10^{19}$m*-3* ($f_{mom}=0.73$, $f_{pow}=0.78$); $10$% at $n_{up}=3\times 10^{19}$m*-3* ($f_{mom}=0.87$, $f_{pow}=0.92$); $80$% error at $n_{up}=4\times 10^{19}$m*-3* ($f_{mom}=0.91$, $f_{pow}=0.98$). Some initial comparisons have also been made to the SOLPS-ITER code [6, 32], with agreement of the order of $20$%, but this is the subject of ongoing work, in particular to understand the validity of the fluid neutral model used here. Grid convergence tests have been performed for these simulations, and many of the operators used in SD1D have been verified using the Method of Manufactured Solutions (MMS) [15] though not yet the full model. Simulation results shown here used 800 grid cells, with convergence tests done at half and double resolution. Grid cells are packed close to the target so that the resolution parallel to the magnetic field is $3.8$mm at the target and $7$cm upstream.

4 The role of particle and energy loss

The parallel heat flux to the divertor target $q_{||}$ consists of plasma thermal energy, and the ionisation energy $E_{iz}$ released by surface recombination. Both of these heat fluxes are proportional to the particle flux $\Gamma$. Reducing target particle flux is therefore considered here as the important outcome of detachment for extending high power divertor lifetime.

Volume recombination of plasma ions at low temperatures (of the order of 1eV) provides a mechanism by which plasma ion flux to the target can be reduced, converting plasma flux to neutral flux before the plasma reaches the target. To examine the importance of recombination to these 1D results we have performed the same density scan but with recombination turned off.

Results are shown in figure 2, and indicate that recombination plays little role at the rollover in flux at an upstream density of $n_{up}\sim 2\times 10^{19}$m*-3*, but becomes important at higher upstream densities of around $4\times 10^{19}$m*-3*. This makes sense because target $T_{e}\simeq 1.2$eV at this higher upstream density, whereas at target flux rollover target $T_{e}\simeq 3.2$eV. This is in contrast to some previous results [33] in which flux rollover occurred at lower temperatures, but in agreement with recent experiments [34, 18] and simulations [35] of TCV. There are several ways in which target flux can be reduced, recombination being one of them, so discrepancies may be seen between studies in different regimes.

From particle conservation the target flux is given by the sum of the flux from the SOL into the divertor and the ionisation flux, less the ions lost to recombination. Other mechanisms such as cross-field transport may also contribute, but are not included in this 1D model, and will typically average to zero over the divertor volume. In these simulations the recycling coefficient is set to $f_{recycle}=0.99$ so the upstream flux is small, and in the absence of significant recombination the drop in target flux must be due to a drop in ionisation. The input power is fixed, so this drop in ionisation is due to either a reduction in power available for ionisation (due to increased impurity radiation losses), or an increase in the effective energy cost per ionisation $E_{iz}$, or both. This effective ionisation cost is calculated as the ratio of the total power lost to hydrogen ionisation and excitation, divided by the rate of ionisations. $E_{iz}$ increases as $T_{e}$ falls, due to excitation radiation, in this simulation from $\sim 31$eV per ionisation at $n_{up}=1.4\times 10^{19}$m*-3* to $\sim 80$eV per ionisation at rollover (averaged over the domain; shown in more detail in figure 6).

One way to quantify this power limitation picture of detachment is in terms of the fraction of power available which goes into hydrogen excitation and ionisation, with the remainder being transmitted through the sheath [18]. In these simulations the radiation and ionisation regions overlap, complicating the analysis. Nevertheless, at flux rollover the fraction of the power into the divertor going to ionisation increases from $\sim 40$% to $\sim 75$% during the strong drop in $T_{target}$, supporting the association between power limitation and flux rollover [17, 36, 8, 33, 18].

To test the importance of power dissipation, figure 3 data labelled “No radiation ($E_{iz}=13.6$eV)” shows the effect of removing all power radiation mechanisms, so that the energy cost per ionisation is only the ionisation potential $E_{iz}=13.6$eV, and there is no impurity radiation. In this case rollover during detachment does not occur in the range of upstream density studied, though the target temperature falls to $3.3$eV at an upstream density of $7.5\times 10^{19}$m*-3*, comparable to the target temperature at flux rollover in the case with radiation included, at an upstream density of $1.89\times 10^{19}$m*-3*. The fraction of available power going to ionisation only reaches $38$% in this case, which though significant is less than the original case discussed above.

A similar result is found when a fixed ionisation cost of $E_{iz}=30$eV is used, to account for some hydrogen excitation radiation, also shown in figure 3. In this case the fraction of available power going to ionisation reaches $59$%. Assuming all of the divertor input power, $P_{in}$, is used for ionisation, and taking into account the $3.5$eV Franck-Condon energy with which dissociated neutral atoms are recycled, the upstream input power would lead to a target flux of $P_{in}/\left(E_{iz}-3.5\mathrm{eV}\right)=5.9\times 10^{24}$m*-2s-1* (in the limit of zero target temperature). The highest fluxes reached in the $E_{iz}=30$eV scan, $\sim 4\times 10^{24}$m*-2s-1*, imply that around $68$% of the input power is going to ionisation. Though the rate of change of target flux with upstream density is reducing, no rollover ($d\Gamma/dn_{up}<0$) in the target flux is seen in this case. A scan with $E_{iz}=60$eV still does not produce a flux rollover.

This lack of flux rollover in cases with fixed ionisation cost can be understood by considering the power balance and temperature dependence of power losses. In the absence of significant recombination the target particle flux is given by the ionisation rate. The input power into the simulation, $P_{in}$, therefore goes either to impurity radiation ($P_{imp}$), ionisation ($P_{iz}$), or to target heat flux ($P_{target}$) as given in equation 1d

[TABLE]

where $\Gamma$ is the target particle flux, labelled “pow” here since this expression is derived from power balance. The sheath heat transmission coefficient $\gamma$ is set to $6$ in these simulations, and $T_{target}$ the temperature at the target (assumed to be the same as at the sheath entrance). $P_{recl}$ is the power available for ionisation, some of which is used for ionisation $P_{iz}=E_{iz}\Gamma^{pow}$, and the remainder going to the target in the form of kinetic energy $P_{target}=\gamma T_{target}\Gamma^{pow}$. $E_{iz}$ is the energy cost of each ionisation, consisting of the ionisation potential (13.6eV) and the radiation due to excitations preceding ionisation. From equation 1d we can write down the power balance constraint for the target ion flux

[TABLE]

If $E_{iz}$ and $P_{in}$ are fixed, and impurity radiation is negligible, then with $\Gamma$ cannot be reduced without increasing $T_{target}$, clearly not the desired outcome. The detached solution with a reduction of both target temperature and particle flux is therefore not available. This is discussed in more detail in section 5.

In this picture there are three ways in which particle flux $\Gamma$ can be reduced: recombination; an increasing $E_{iz}$ as $T_{target}$ drops; or impurity radiation $P_{imp}$, which reduces power available for ionisation. In the original “Carbon + Hydrogen” case shown in figure 3 the dominant mechanism is an increase in $E_{iz}$ (hydrogen radiation) as $T_{target}$ drops. In the case labelled “Carbon $E_{iz}=30$eV” with a $1$% fixed fraction carbon impurity and $30$eV fixed ionisation cost, particle flux rollover is also achieved, with the fraction of available power going to ionisation ($P_{iz}/P_{recl}$) around $52$% at rollover ($n_{up}=3.5\times 10^{19}$m*-3*). High power fractions going to ionisation appears to be associated with flux rollover, consistent with the power limitation picture, but is not by itself sufficient to produce flux rollover. To understand this further we next examine the role of momentum losses and the description of target ion flux as in equation 1e but based on momentum ($\Gamma^{mom}$).

5 The role of momentum loss

As briefly mentioned in section 3, one of the robust experimental characteristics of detachment is loss of target plasma pressure, implying a loss of momentum to neutrals or other sinks. Simple arguments indicate that the target particle flux varies as in equation 1c, reproduced here as equation 1f and labelled $\Gamma^{mom}$ to indicate that the origin of this relationship is the momentum balance:

[TABLE]

where the pressure $p_{target}$ and temperature $T_{target}$ are evaluated at the divertor target. In the absence of momentum losses, a sonic boundary condition implies that the static target plasma pressure is proportional to the upstream pressure $p_{up}=2p_{target}$. As mentioned, a fall in the target pressure, faster than the drop in $\sqrt{T_{target}}$, is therefore necessary to reduce target particle flux, overcoming the increase in flux associated with a drop in target temperature. We note here that the drop in $p_{target}$ can be due either to volumetric momentum losses or a drop in the upstream pressure, or some combination of both [18]. For a steady-state solution to exist, the plasma must self-organise into a configuration which satisfies both this momentum constraint as well as particle flux and power balance constraints discussed in section 4.

A measure of the pressure loss from X-point to target is the ratio of target pressure $p_{target}$ to upstream pressure $p_{up}$, shown in figure 4. As has been shown previously [37, 38, 16], this pressure ratio is found to be mainly a function of target temperature, and here is shown to be independent of the power dissipation model (Carbon, Hydrogen radiation, or fixed $E_{iz}$). In addition, a case with no flux expansion (Constant B field), and a case where 50% of recycled neutrals are distributed evenly along the divertor leg (Redistributed neutrals) follow the same trend. The one case which follows a different trend in figure 4 is labelled “CX neutrals escape”, in which charge-exchanged neutrals escape the SOL, lose their momentum, and return at a random location along the divertor leg. This results in a loss of total (plasma + neutral) momentum, modifying $f_{mom}$. As discussed in section 2, the treatment of neutral momentum is likely to be important in quantatively matching experiment.

As in [7, 16], the pressure ratio and hence the fraction of total momentum lost, $f_{mom}$, fits a function of the form in equation 1g and shown in figure 4:

[TABLE]

The fit coefficients found here are comparable to those from experimental data e.g. C-MOD [37] $1-f_{mom}=1.2\left(1-\exp\left(-T_{target}/2.3\right)\right)^{4}$ and SOLPS modelling of AUG [16] H-mode $1-f_{mom}=0.8\left(1-\exp\left(-T_{target}/2\right)\right)^{1.2}$. Possible sources of discrepancy include radial momentum transport, and variations in how plasma momentum is transferred to neutrals and the walls of the device discussed further below. As discussed elsewhere [38, 16], this strong dependence on target temperature is consistent with a simple model in which neutrals ionise in a narrow region close to the target, so that the neutral density and electron temperature can be considered homogenous over the interaction region. In this Self-Ewald model [39], given in equation 1ha, $f_{mom}$ is only a function of ionisation and charge exchange rate coefficients, which to first approximation depend only on the temperature:

[TABLE]

where $\left<\sigma v\right>_{iz}$ and $\left<\sigma v\right>_{cx}$ are the ionisation and charge-exchange rate coefficients. This Self-Ewald solution overestimates the momentum loss, as shown in figure 4, since it assumes that all charge exchange events remove momentum from the system, in addition to occuring in an isothermal environment.

The drop in pressure from upstream to target shown in figure 4 is commonly associated with detachment, but is not sufficient to lead to rollover: Here all scans with conserved momentum follow the same pressure ratio curve, but the scans with fixed energy cost per ionisation do not have a roll over in target flux, as shown in figure 3.

To further understand the lack of rollover in target flux in cases with fixed energy cost per ionisation $E_{iz}$, figure 5 shows the target flux against target temperature for three scans with no impurity radiation and different fixed values for $E_{iz}$.

As discussed in section 4, power balance dictates that the target flux and target temperature follows equation 1d, shown as solid black lines for each of the three values of $E_{iz}$. In steady state this solution must also be consistent with momentum balance (equation 1f), shown as dashed red curves for three different upstream densities. The intersection of the solid black (power) and dashed red (momentum) curves is therefore the consistent steady-state solution.

The expression for flux from momentum balance (equation 1f) can be written as the upstream pressure multiplied by a function $F\left(T_{target}\right)$ which only depends on target temperature (equation 1hi):

[TABLE]

In all SD1D cases discussed here the variation in upstream temperature is small, $57-61$eV, so $p_{up}\propto n_{up}$. The dashed red curves in figure 5 therefore scale up and down in proportion to the upstream pressure ($\sim$ density), whilst retaining a constant shape $F\left(T_{target}\right)$. The black curves scale in proportion to the input power, whilst also retaining a constant shape. As these inputs are changed, the intersection of the two curves moves. Since for large $T_{target}$ the flux from momentum balance (equation 1f) goes like $\Gamma^{mom}\sim 1/\sqrt{T}$, whilst the flux from power balance (equation 1d) goes like $\Gamma^{pow}\sim 1/T$, there is always a solution for these scans. It can be seen in figure 5 that as the upstream density is increased the red curve ($\Gamma^{mom}$) moves upwards; this causes the intersection point to move to the left, corresponding to a lowering of target temperature and an increase in target flux, which asymptotically approaches the limit $\Gamma_{max}=P_{in}/E_{iz}$. Note that recombination could produce a rollover in these scans, by removing flux before it reaches the target, but is not sufficient in the cases studied here.

Having understood why the cases with fixed $E_{iz}$ do not have a rollover in target flux, we now turn to the more experimentally relevant cases which do have a flux rollover. From figure 5 it can be seen that in order to have a rollover in target flux, the black curve derived from power balance (equation 1d), must change slope so that $\frac{d\Gamma^{pow}}{dT}>0$. This is in addition to the requirement for momentum loss which determines the shape of the red curve.

Assuming a fixed input (upstream) power, $P_{in}$, the power to the recycling region is given by $P_{recl}=P_{in}-P_{imp}$ (equation 1d), where $P_{imp}$ is the impurity radiation [18]. The gradient of the flux with respect to target temperature is given in equation 1hj.

[TABLE]

The gradient of $\Gamma^{pow}$ can therefore be changed by either impurity power radiation which depends on target temperature, or through an $E_{iz}$ which depends sufficiently strongly on target temperature.

Impurity radiation can drive a change in sign of $\frac{d\Gamma^{pow}}{dT}$, but is complicated by the dependence of $P_{imp}$ on plasma density. Here we consider the case without impurity radiation, $P_{imp}=0$. This is relevant to the scan labelled “Carbon + Hydrogen” since in those cases the hydrogen radiation dominates the power loss. The criterion for a change in gradient of $\Gamma^{pow}$, and so a rollover of target flux, is therefore given by equation 1hk:

[TABLE]

The effective value of $E_{iz}$ from the SD1D simulation with varying $E_{iz}$ (labelled “Carbon + Hydrogen”) is shown in figure 6, averaged over the field line for each case corresponding to the solution $T_{target}$.

Since $E_{iz}$ is an increasingly steep function of temperature as the target temperature falls, equation 1hk can be related to a threshold for rollover at $T_{target}=7.84$eV, $E_{iz}=60.8$eV in the red curve of figure 6 (point-wise value, not averaged over the domain). By combining equations 1hi and 1d, this threshold in the gradient of $E_{iz}$ can also be related to a critical ratio of upstream pressure to recycling power $p_{up}/P_{recl}$, which has been proposed as a detachment threshold [36, 33, 18], and which is a function only of target temperature as shown in equation 1hla:

[TABLE]

where equation 1g has been used for $f_{mom}$. This value is smaller than the $\sim 17$N/MW threshold found in [33] for SOLPS4.3 simulations of pure Deuterium DIII-D-like equilibria, but larger than that derived and measured [18]. This difference may be due to different effective $E_{iz}$, which here uses a simplified semi-analytic model, to differences in $\gamma$, or to a difference in $f_{mom}$, which varies somewhat between devices and simulations [16]. The effective $E_{iz}$ averaged over the domain (data points in figure 6) deviates from $E_{iz}$ at the target temperature (curve in figure 6), which also shifts equation 1hla to higher $p_{up}/P_{recl}$.

The target flux for this scan with varying $E_{iz}$ (labelled “Carbon

• Hydrogen” in figures 1 and 3) is shown in figure 7. This is similar to figure 5, but now the black and grey curves which represent the power balance constraint (equation 1d) are not monotonic and so a flux rollover is possible. The grey curve shows the power balance obtained by using the target temperature in the $E_{iz}$ rates used in SD1D. This fits the simulation results well at high target temperatures, but at an upstream density $n_{up}\simeq 1.88\times 10^{19}$m*-3* a transition is seen in the target temperature. Figure 7 shows that the reason for this transition is that at the corresponding target temperature $T_{target}=11.3$eV the gradients of the grey and red curves are such that $\frac{\partial\Gamma^{pow}}{\partial T}>\frac{\partial\Gamma^{mom}}{\partial T}$ so that these curves no longer intersect as the upstream density is increased. This loss of steady-state solution results in a rapid change in density and temperature profiles, approximately doubling the peak density. This moves the peak in ionisation away from the target, modifying the power dissipation curve (from the grey to the black line in figure 7) with an average $E_{iz}$ which deviates from the value at $T_{target}$, as shown in figure 6. The solution moves to a different point on the same red curve at a similar upstream density ($1.88\times 10^{19}$m*-3* to $1.8825\times 10^{19}$m*-3*), but now with a target temperature of $3.8$eV rather than $11.3$eV.

Rapid transitions into detachment are common (e.g. [40, 41, 42]) but are undesirable in a tokamak divertor, both for design and control. The above discussion implies that power dissipation mechanisms which dissipate strongly at high temperatures, particularly at target temperatures above the peak in $\Gamma^{mom}$ (equation 1hi), have the potential to result in sudden transitions. This is because at these high temperatures $\frac{\partial\Gamma^{mom}}{\partial T}<0$ so that sufficiently strong power dissipation can produce a situation where $\frac{\partial\Gamma^{pow}}{\partial T}>\frac{\partial\Gamma^{mom}}{\partial T}$, the black and red curves in figure 7 no longer intersect, and a smooth change is no longer possible. For these simulations the peak in $\Gamma^{mom}$ occurs at $6.0$eV when the fit in equation 1g is used. This analysis points to a possible solution: mechanisms which remove momentum at higher target temperatures, perhaps including radial transport, which modify $f_{mom}\left(T_{target}\right)$ and increase $\frac{\partial\Gamma^{mom}}{\partial T}$, would reduce the likelihood of sudden transitions.

We conclude that when recombination is not significant, momentum loss without sufficient power loss results in a lack of target flux rollover (figure 5); power loss without sufficient pressure loss can result in rapid transitions between states (figure 7). More specifically, a change in gradient of $\Gamma^{pow}$ w.r.t $T_{target}$ is needed for target flux rollover. If this change in gradient is driven by hydrogen excitation radiation, then there is a threshold in the gradient of $E_{iz}\left(T_{target}\right)$ (equation 1hk) which can be related to a threshold in $T_{target}$ (figure 6) and to a threshold in $p_{up}/P_{recl}\simeq 12.6$ N/MW (equation 1hla). If this threshold is reached at a higher temperature than the rollover in $\Gamma^{mom}$ at $\sim 6$eV then a rapid transition may occur (figure 7), since a smooth variation in plasma state is no longer possible.

5.1 Removing charge exchange

In previous sections the effect of removing power dissipation mechanisms has been studied, identifying the potential for transitions to occur when insufficient momentum is removed. We now test the effect of removing the dominant plasma momentum loss mechanism in these simulations, charge exchange. At rollover the pressure loss due to charge exchange is $16$ times larger than the sum of all other effects; the next most significant effect is ionisation ($15.4$ times smaller) which acts in the opposite direction. Given the discussion above of equation 1f, it might be expected that turning off charge exchange would raise target flux and inhibit the roll-over during detachment. What is seen in figure 8 is more nuanced: At low densities the target flux is indeed increased relative to the case with charge exchange, but the sharp drop in target temperature and roll-over of target flux associated with detachment is seen to occur at a lower upstream density.

The plasma profiles are qualitatively different in the cases with and without charge exchange: When charge exchange is included, conservation of momentum between plasma and neutral fluids results in a static solution shown in figure 9a, in which the plasma fluid pressure is balanced by a neutral cushion in front of the target: since neutrals do not escape, the total plasma + neutral pressure is conserved. When charge exchange is turned off, this static solution is not possible because the required momentum exchange cannot occur. Instead in the solution shown in figure 9b a part of the plasma momentum is lost to neutral pressure through ionisation ($102$ Pa) and some through the magnetic mirror effect [25] ($20$ Pa), but the majority of the upstream static plasma pressure ($313$ Pa) is balanced by plasma dynamic pressure $m_{i}nv_{||}^{2}$ ($172$ Pa). Power is radiated in a quasi-neutral thermal front (gradient) region, and acceleration of the plasma into the low pressure region behind the thermal front converts internal energy to kinetic energy. In the absence of momentum loss this results in a cold supersonic flow of plasma to the target, reminiscent of a de Laval rocket nozzle [43]. Note that in this case anisotropic pressure (viscous) effects are likely to be important [25], but are not included here. The model is quasi-neutral by construction, so if space charge effects play a role in the thermal front then these would not be observed here.

The rarefaction front at which plasma static pressure drops moves upstream as the upstream density is increased. This results in a larger drop in the upstream temperature: Between upstream densities $n_{up}=1.7\rightarrow 3\times 10^{19}$m*-3* the upstream temperature falls from $58.8\rightarrow 58.2$eV if charge exchange is included, and $58\rightarrow 43$eV if charge exchange is excluded. Over the same range of upstream densities the front moves from $L_{||}\simeq 29$m to $L_{||}\simeq 13$m, close to the X-point at $10$m. The drop in upstream temperature is consistent with a simple 2-point scaling [7] with parallel connection length $T^{upstream}\sim L_{||}^{2/7}$, where $L_{||}$ is the location of the density peak. Such a use of the 2-point scaling treats the density peak as a ’virtual target’ [44, 45].

The transition at lower upstream density in the case without charge exchange, compared to the case with charge exchange, can be understood by considering the power and momentum balance as done in figures 5 and 7.

This is shown in figure 10: Removing charge exchange momentum losses reduces $f_{mom}$ and so increases $\Gamma^{mom}$ (equation 1hi). This moves the location where the $\Gamma^{mom}$ curve intersects the $\Gamma^{pow}$ curve, so that the same target temperature occurs at a lower upstream density. This lowers the upstream density at which a transition is observed, relative to the case with charge exchange (figure 7). Setting $f_{mom}=0$, the momentum balance curves (red dashed lines) in figure 10 cease to insersect the power balance curve at an upstream density of around $n_{up}=1.7\times 10^{-3}$m*-3*, at which point a transition is observed in the SD1D simulation.

It seems unlikely that this solution with supersonic flow could ever occur in tokamak experiments, though the possibility of transitions to supersonic flow in divertor plasmas has been shown previously [25]. The dramatic change in plasma solutions when momentum exchange is modified demonstrates the importance of momentum exchange to determining the divertor detachment thresholds and dynamics. It also shows the utility of the analysis in figure 10 to explaining the observed simulation results.

6 Conclusions

We have undertaken steady-state, 1D, simulations of divertor plasma detachment, specifically target ion flux rollover, utilizing a new computational tool, SD1D. The importance of recombination, radiative power loss, and momentum exchange in the detachment process with emphasis on the behavior of the target ion current rollover has been evaluated.

We find that for MAST-Upgrade like simulation parameters recombination does not play a significant role at flux rollover (target temperature $T_{target}\sim 3-5$eV), though is significant when $T_{target}$ drops to $\sim 1$eV.

It is shown that in these simulations momentum loss, as characterised by a drop in total pressure along the magnetic field, is insufficient, by itself, to produce flux rollover during detachment. Impurity and/or hydrogenic radiative losses are also required to increase at low temperatures. In the particular case studied where excitation power losses dominate over impurity radiation, this corresponds to $dE_{iz}/dT_{target}<-\gamma$. We conclude that when recombination is not significant, momentum loss without sufficient power loss results in a lack of target flux rollover.

The precise dependence of momentum and power losses on target temperature has implications for the availability of steady state solutions. We have found that when the target flux is characterized in terms of power balance, $\Gamma^{pow}\left(T_{target}\right)$, and momentum balance, $\Gamma^{mom}\left(p_{up},T_{target}\right)$, the allowed intersections of the two curves readily predict available steady state solutions. In some cases (e.g. power dissipation occurring without sufficient momentum removal) we find that there are regions of $T_{target}$ where there are no solutions – leading to rapid drops in temperature for essentially no increase in upstream density.

One implication of this work is that modification of the shape of $\Gamma^{pow}\left(T_{target}\right)$ and $\Gamma^{mom}\left(p_{up},T_{target}\right)$ can be used to influence detachment behavior. One can envision causing such shape modifications by changing the escape probability of CX neutrals carrying momentum (e.g. changing divertor target geometry) or changing the temperature dependence of radiative losses (e.g. switching from one seeded impurity to another).

The SD1D code is available at https://github.com/boutproject/SD1D. All inputs, data and processing scripts are available at DOI 10.5281/zenodo.1410281.

Acknowledgements

This work has received funding from the EPSRC under grant EP/N023846/1. B. Lipschultz was funded in part by the Wolfson Foundation and UK Royal Society through a Royal Society Wolfson Research Merit Award as well as by the RCUK Energy Programme (EPSRC grant number EP/I501045).

This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission