SPLEND1D, a reduced one-dimensional model to investigate the physics of plasma detachment

Février, O.; Gorno, S.; Theiler, C.; Carpita, M.; Durr-Legoupil-Nicoud, G.; von Allmen, M.

doi:10.1063/5.0202986

Studying the process of divertor detachment and the associated complex interplay of plasma dynamics and atomic physics processes is of utmost importance for future fusion reactors. While simplified analytical models exist to interpret the general features of detachment, they are limited in their predictive power, and complex two-dimensional (2D) or even three-dimensional (3D) codes are generally required to provide a self-consistent picture of the divertor. As an intermediate step, one-dimensional (1D) models of the scrape-off layer (SOL) can be particularly insightful as the dynamics are greatly simplified, while still self-consistently including various source and sink terms at play, as well as additional important effects such as flows. These codes can be used to shed light on the physics at play, to perform fast parameter scans, or to interpret experiments. In this paper, we introduce the SPLEND1D (Simulator of PLasma ENabling Detachment in 1D) code: a fast and versatile 1D SOL model. We present in detail the model that is implemented in SPLEND1D. We then employ the code to explore various elements of detachment physics for parameters typical of the Tokamak à Configuration Variable, including the atomic physics and other processes behind power and momentum losses, and explore the various hypotheses and free parameters of the model.

I. INTRODUCTION

Plasma power and particle exhaust is a crucial issue for future fusion reactors. If unmitigated, the target heat fluxes in ITER and DEMO are expected to greatly exceed the $10 {MW / m}^{2}$ limit that is considered necessary for tolerable steady-state conditions. It will thus be necessary to operate in an, at least partially, detached regime.^1–3 In such a regime, the total plasma pressure develops strong parallel gradients along the field lines in the scrape-off layer (SOL), driven by volumetric and perpendicular momentum losses,^4,5 and a significant fraction of the plasma power is dissipated by volumetric power sinks. This results in a reduction of the target heat and particle flux densities, as well as target temperature, which is important to reduce target erosion.^6,7 The detachment process typically sets in at a target temperature below 5 eV,^5,8,9 where plasma-neutral interactions are further enhanced, further reducing target heat flux, temperature, and ion flux.

While simplified analytical models exist to interpret the general features of detachment, they are limited in their predictive power, and complex 2D or even 3D codes are generally required to provide a self-consistent picture of the divertor. However, the results of these codes can be difficult to interpret. Furthermore, the numerical cost of running such simulations generally prohibits large parameter scans that could be useful to shed light on the role of the various terms entering the equations. Therefore, as an intermediate step, one-dimensional (1D) models of the scrape-off layer (SOL) can be particularly insightful as the dynamics are greatly simplified, while still self-consistently including various source and sink terms at play, as well as additional important effects such as flows. Typically based on a projection of the Braginskii's equations¹⁰ along a field line, these models simplify the perpendicular dynamics either by ignoring it, or by simplifying it with some simple ansatz. Such models have been typically employed in the context of a specific study.^11–13 More recently, efforts have focused on implementing dedicated 1D codes that can be shared, reused and coupled with other codes or integrated into workflows for integrated modeling. Examples of such codes are SOLF1D,¹⁴ SD1D,¹⁵ and DIV1D.¹⁶ In this, paper, we introduce the Simulator of PLasma ENabling Detachment in 1D (SPLEND1D) code, a fast and versatile 1D model used to explore the complex interplay of atomic physics and plasma dynamics underlying the detachment process. The purpose of this paper is twofold. First, to present in great detail the model implemented in the SPLEND1D code and demonstrate the speed and robustness of the code. Second, to apply the code for the study of a base case scenario for parameters typical for the Tokamak à Configuration Variable (TCV).¹⁷ The neutral particle source is increased in this scenario to reach the onset of detachment, similarly to the experimental onset of detachment through increased plasma fueling. The different processes are subsequently investigated, thus providing a first illustration of the SPLEND1D capabilities and possible applications. This paper is organized as follows. In Sec. II, the derivation of the models for the charged species (plasma) and the neutral species is presented, as well as the numerical methods, boundary conditions, and implemented source terms. Section III presents a reference base case that is used to demonstrate the capabilities of the SPLEND1D code in terms of ease of use and result interpretation. We also present some measurements of the code accuracy and runtime. In Sec. IV, we simulate detachment of the base case through upstream density ramps, and highlight the mechanisms at play in the model that enable the onset of momentum and energy losses. Section V investigates the role of some free parameters of the model, such as the impurity concentration, the neutral confinement time, the heat-flux limiters, and the choice of boundary condition for the parallel velocity. Finally, in Sec. VI, we present some advanced studies that have been enabled by SPLEND1D, investigating the choice of the neutral model, the effect of separating ion and electron energy equations, and the results of time-dependent simulations. These advanced studies demonstrate SPLEND1D's aptness in interpreting current TCV experiments, for example to study the role of connection length on divertor detachment.¹⁸ Comparisons to other 1D codes in the community are discussed throughout the text.

II. MODEL

In this first section, we describe the model implemented in SPLEND1D. The equations solved by SPLEND1D are based on the Braginskii equations,¹⁰ which are typically used to describe the evolution of the plasma in the SOL in 1D, 2D or 3D codes, using a fluid approximation. In particular, we present in detail the assumptions used to derive the SPLEND1D code.

A. Geometry

We consider a one-dimensional model for the SOL. The geometry of this 1D SOL is axisymmetric, i.e., uniform along

φ

in Fig. 1, but otherwise arbitrary. In particular, both the magnitude (B) of the magnetic field (B) as well as its pitch angle are allowed to vary along a field line. We introduce the curvilinear coordinate s along the magnetic field, such that

\vec{d s}

is parallel to the unit vector

b = \frac{B}{B}

⁠, Fig. 1. Here, s is defined to increase toward the target, Fig. 1. We denote

α (s)

the local pitch angle. It can be related to the components of the magnetic field by

tan α = \frac{B_{θ}}{B_{φ}},

(1)

where

B_{θ} (s)

and

B_{φ} (s)

are the poloidal and toroidal components of the magnetic field, Fig. 1.

FIG. 1.

View large Download slide

Definition of the geometry used in SPLEND1D. b is the unit vector parallel to the magnetic field. α corresponds to the pitch angle. s is a curvilinear coordinate along the flux tube, oriented along ds. θ and $φ$ are two spatial coordinates that will be relevant to the 2D neutral model that is presented in Sec. II C.

B. Plasma model

1. Equations

In the following, we assume toroidal symmetry within the system, such that

\partial / \partial φ = 0

for all considered quantities. The plasma fluid velocity is assumed purely parallel to the magnetic field, and to be the same for ions and electrons, such that there is no parallel current. The velocity vector V is written as

V = u_{| |} b

⁠. The conductive and viscous heat fluxes are also assumed to be parallel to the magnetic field. We denote n the plasma density, assuming quasi-neutrality (⁠

n_{i} = n_{e} = n

with Z = 1), and m_e (respectively m_i) the electron (respectively ion) mass. The neutral mass m_n is taken equal to the ion mass, m_n = m_i. T_i, and T_e are the ion and electron temperatures. We assume high enough collisionality, such that, for both ions and electrons, the pressure is isotropic. With these assumptions, we project the Braginskii equations along b, which results in the following continuity equation, total parallel plasma momentum equation (obtained by summing up the electron and ion momentum equations assuming

m_{e} ≪ m_{i}

⁠), and electron and ion energy equations:

\frac{\partial n}{\partial t} + B \frac{\partial}{\partial s} (\frac{n u_{| |}}{B}) = S_{p}^{n} + H_{P},

(2)

\frac{\partial}{\partial t} (m_{i} n u_{| |}) + B \frac{\partial}{\partial s} (\frac{m_{i} n u_{| |}^{2}}{B}) = - \frac{\partial}{\partial s} (p_{e} + p_{i}) + S_{p, | |}^{u} - b \cdot \nabla \cdot Π,

(3)

\frac{\partial}{\partial t} (\frac{3}{2} n T_{e}) + B \frac{\partial}{\partial s} (\frac{\frac{5}{2} n T_{e} u_{| |} + q_{| |, e}^{cond}}{B}) = u_{| |} \frac{\partial p_{e}}{\partial s} + Q_{e} + S_{e}^{E} + S_{imp}^{E} + H_{e},

(4)

\frac{\partial}{\partial t} (\frac{3}{2} n T_{i} + \frac{1}{2} m_{i} n u_{| |}^{2}) + B \frac{\partial}{\partial s} (\frac{\frac{5}{2} n T_{i} u_{| |} + \frac{1}{2} m_{i} n u_{| |}^{3} + q_{| |, i}^{cond}}{B}) = - u_{| |} \frac{\partial p_{e}}{\partial s} + Q_{i} + S_{i}^{E} + H_{i} - \nabla \cdot (u_{| |} Π_{| |}) .

(5)

Here, quantities are defined in SI units. As for the temperatures, T stands for

k_{B} T

⁠, where

k_{B}

is the Boltzmann constant.

p_{e} = n T_{e}, p_{i} = n T_{i}

are the static pressures of the electrons and ions, respectively.

S_{p}^{n}, S_{p, | |}^{u}, S_{e}^{E}

and

S_{i}^{E}

are particle, momentum, electron energy and ion energy source terms resulting from ionization, recombination, charge-exchange and excitation reactions. They will be described in Sec. II D. The term

S_{imp}^{E}

in Eq. corresponds to the energy loss due to impurity radiation, also described in Sec. II D.

q_{| |, e}^{cond}

(respectively

q_{| |, i}^{cond}

⁠) is the electron (respectively ion) parallel conductive heat flux, whose expressions will be given in Sec. II B 2.

b \cdot \nabla \cdot Π

and

\nabla \cdot (u_{| |} Π_{| |})

are the viscous contributions to the parallel momentum and energy equations, detailed in Sec. II B 3. Q_e and Q_i model the exchange of energy between ions and electrons due to collisions and can be expressed as

Q_{e} = - Q_{i} = 3 \frac{m_{e}}{m_{i}} \frac{n}{τ_{e}} (T_{i} - T_{e}),

(6)

where τ_e is the electron collision time, defined as^10,19

τ_{e} = \frac{6 \sqrt{2} π^{3 / 2} ϵ_{0}^{2} \sqrt{m_{e}} T_{e}^{3 / 2}}{ln Λ e^{4} n},

(7)

\approx \frac{1}{3.64 \times 10^{- 6}} \frac{{(T_{e} [K])}^{3 / 2}}{ln Λ n [m^{- 3}]},

(8)

\approx \frac{1}{2.91 \times 10^{- 12}} \frac{{(T_{e} [eV])}^{3 / 2}}{ln Λ n [m^{- 3}]},

(9)

with

ln Λ

the Coulomb Logarithm. H_P in Eq. is an additional volumetric (charged) particle source term and is an input for the code. It can be used to model a flux of particles entering or leaving the flux tube. Similarly, H_e and H_i in Eqs. and are volumetric energy source terms.

To reduce the number of degrees of freedom, the code can be run under the assumption

T_{i} = \bar{τ} T_{e}

⁠, with

\bar{τ} > 0

an arbitrary constant. Summing Eqs. and , Eqs. can then be rewritten as

\frac{\partial n}{\partial t} + B \frac{\partial}{\partial s} (\frac{n u_{| |}}{B}) = S_{p}^{n} + H_{P},

(10)

\frac{\partial}{\partial t} (m_{i} n u_{| |}) + B \frac{\partial}{\partial s} (\frac{m_{i} n u_{| |}^{2}}{B}) = - (1 + \bar{τ}) \frac{\partial p_{e}}{\partial s} + S_{p, | |}^{u} - b \cdot \nabla \cdot Π,

(11)

\begin{matrix} \frac{\partial}{\partial t} (\frac{3}{2} n (1 + \bar{τ}) T_{e} + \frac{1}{2} m_{i} n u_{| |}^{2}) + B \frac{\partial}{\partial s} (\frac{\frac{5}{2} n (1 + \bar{τ}) T_{e} u_{| |} + \frac{1}{2} m_{i} n u_{| |}^{3} + q_{| |, i}^{cond} + q_{| |, e}^{cond}}{B}) \\ = S_{i}^{E} + S_{e}^{E} + S_{imp}^{E} + H_{e} + H_{i} - \nabla \cdot (u_{| |} Π_{| |}) . \end{matrix}

(12)

SPLEND1D is able to solve either the two fluids model [Eqs. ] or the one-fluid model [Eqs. ], depending on the user inputs.

2. Heat fluxes

The electron and ion parallel conductive heat fluxes are defined as

\frac{1}{q_{| |, e}^{cond}} = \frac{1}{q_{| |, e}^{lim}} + \frac{1}{q_{| |, e}^{S H}},

(13)

and

\frac{1}{q_{| |, i}^{cond}} = \frac{1}{q_{| |, i}^{lim}} + \frac{1}{q_{| |, i}^{S H}},

(14)

where

q_{| |, e}^{S H}

and

q_{| |, i}^{S H}

are defined by

q_{| |, e}^{S H} = - κ_{e}^{S H} \frac{\partial T_{e}}{\partial s}, q_{| |, i}^{S H} = - κ_{i}^{S H} \frac{\partial T_{i}}{\partial s},

(15)

where

κ_{e}^{S H}

and

κ_{i}^{S H}

⁠, the classical Spitzer-Härm electron and ion heat conduction coefficients, are defined as¹⁰^,

κ_{e}^{S H} = 3.16 \frac{n τ_{e} T_{e}}{m_{e}},

(16)

and

κ_{i}^{S H} = 3.9 \frac{n τ_{i} T_{i}}{m_{i}} .

(17)

Here, τ_e is the electron collision time defined in Eq. and τ_i is the ion collision time defined as^10,19

τ_{i} = \frac{12 π^{3 / 2} ϵ_{0}^{2} \sqrt{m_{i}} T_{i}^{3 / 2}}{ln Λ e^{4} n},

(18)

\approx \frac{1}{6.0 \times 10^{- 8}} \sqrt{\frac{m_{i}}{m_{p}}} \frac{{(T_{e} [K])}^{3 / 2}}{ln Λ n [m^{- 3}]},

(19)

\approx \frac{1}{4.8 \times 10^{- 14}} \sqrt{\frac{m_{i}}{m_{p}}} \frac{{(T_{e} [eV])}^{3 / 2}}{ln Λ n [m^{- 3}]},

(20)

where m_p is the proton mass and where we assumed Z = 1. These expressions lead to the well-known Spitzer-Härm conductivity with

κ_{e}^{S H} \propto T_{e}^{5 / 2}

and

κ_{i}^{S H} \propto T_{i}^{5 / 2}

⁠. In low collisionality (long mean free path) regimes, however, the physical heat fluxes may be significantly overestimated by the classical Spitzer-Härm heat fluxes.^1,20–22 In order to avoid this nonphysical divergence, we use so-called “flux limiters”, which limit the maximum value of the heat flux to the free streaming heat flux

q_{| |, {i, e}}^{lim}

⁠, defined as

q_{| |, {i, e}}^{lim} = α_{{i, e}} n T_{{i, e}} v_{{i, e}}^{t h} = α_{{i, e}} n T_{{i, e}} \sqrt{\frac{T_{i, e}}{m_{i, e}}},

(21)

where α_i and α_e are two free parameters, whose typical values are around 0.5,^20,21 although formally one would require kinetic simulations to determine the values of these parameters. By affecting the heat flux, they may influence the predictions of the simulation.^20,21,23 In Sec. V C, we will investigate the sensitivity of the results presented in this paper on these values.

3. Viscosity tensor

SPLEND1D includes the effect of parallel viscosity. The term

b \cdot \nabla \cdot Π

⁠, that appears in Eqs. and , can be written as¹⁰^,

b \cdot \nabla \cdot Π = b \cdot \nabla \cdot [δ p (bb - \frac{1}{3} 1)] = \frac{2}{3} [B^{\frac{3}{2}} \frac{\partial}{\partial s} (δ p B^{- \frac{3}{2}})],

(22)

where

δ p = - η_{| |} (3 \frac{\partial u_{| |}}{\partial s} - B \frac{\partial}{\partial s} (\frac{u_{| |}}{B})),

(23)

and

η_{| |}

is the parallel ion viscosity, defined later in Eq. . This results in

b \cdot \nabla \cdot Π = - \frac{4}{3} [B^{\frac{3}{2}} \frac{\partial}{\partial s} (η_{| |} B^{- 2} \frac{\partial B^{\frac{1}{2}} u_{| |}}{\partial s})] .

(24)

The viscosity also contributes to the ion energy equation. Starting from the contribution of the viscosity to heat generation,

\nabla \cdot (u_{| |} \cdot [δ p (bb - \frac{1}{3} 1)])

⁠, one finds an energy source term of the form

\nabla \cdot (u_{| |} Π_{| |}) = - \frac{4}{3} [B \frac{\partial}{\partial s} (u_{| |} η_{| |} B^{- \frac{3}{2}} \frac{\partial B^{\frac{1}{2}} u_{| |}}{\partial s})],

(25)

where

Π_{| |} = b \cdot Π

⁠. Similar to the classical heat fluxes, the viscous flux can become unphysical at low collisionality and hence requires a flux limiter. This is done by writing

η_{| |} = \frac{η_{| |}^{b r}}{1 + \frac{η_{| |}^{b r} | \frac{\partial u_{| |}}{\partial s} |}{τ_{| |}^{lim}}} + n ν_{num},

(26)

where

τ_{| |}^{lim} = \frac{4}{7} n T_{i}

⁠, and

η_{| |}^{b r}

is the Braginskii ion parallel viscosity, expressed as

η_{| |}^{b r} = 0.96 n T_{i} τ_{i}

⁠. More complete expressions of this flux limiter can be developed,^14,21,24,25 although in this paper we restrict ourselves to this simple form, which is similar to the one implemented in 2D transport codes such as SOLPS-ITER.^14,23 ν_num is a numerical kinematic viscosity that can be employed to facilitate numerical convergence, but that is not enabled by default in the code.

C. Fluid model for the neutrals

This section presents the model retained to describe the dynamics of the neutrals. As the neutrals are not affected by the magnetic field, their dynamics are intrinsically 3D. Furthermore, in typical divertor conditions, the mean free path of the neutrals can be large compared to the system size, leading to a high Knudsen number that in principle requires a kinetic description of the neutral dynamics. This is the approach of EIRENE, one of the main workhorse for neutral dynamics simulations in the fusion community,²⁶ where a Monte Carlo method is used to simulate the behavior of the neutrals. However, such methods are computationally expensive and subject to statistical noise. In recent years, to alleviate these costs, there has been significant work devoted to the development of advanced fluid neutral models, or hybrid neutral models, mixing a kinetic and a fluid description, Refs. 27–29 and references therein. As 1D models aspire for simplicity, in this work, we use fluid neutral models. SPLEND1D implements two different neutral models, the choice of which to use is made by the user: a diffusive neutral model, in which the neutrals diffuse along the θ-direction; and an advective one, where neutrals move in the $θ - φ$ plane (that is, a flux surface), with a velocity that is not necessarily parallel to the magnetic field, similar to the model presented in Refs. 27 and 30. In the following, we further assume that the vectors describing the neutral dynamics (velocities, heat fluxes) do not have a component in the radial direction (⁠ $Ψ$ ⁠). Therefore, while neutrals are not confined to a particular flux tube, they remain confined to a particular flux-surface. This assumption, despite being clearly a simplification, provides the highest fidelity level available for the neutral dynamics description, while maintaining a 1D fluid picture and not requiring further information about the geometry of the problem or the plasma quantities in the divertor. In its current version, SPLEND1D models only a single population of neutral atoms, and does not include contributions from molecules. It is, however, well known from experimental and simulation works that molecules are very strong contributors to the divertor dynamics,^31–33 in particular through elastic collisions with the main ion species. Future works on the SPLEND1D neutral model will therefore encompass the addition of neutral-ion collisions, as well as the addition of molecular effects.

1. Continuity equation for the neutrals

The general form for the continuity equation for the neutral particles is given by

\frac{\partial n_{n}}{\partial t} + \nabla \cdot (n_{n} V_{n}) = S_{n}^{n} + H_{n} - \frac{n_{n}}{τ_{n}},

(27)

where

V_{n}

is the neutral velocity,

S_{n}^{n}

is a neutral particle source (or sink) term due to atomic processes, which will be discussed in Sec. II D, and H_n is an arbitrary neutral source or sink, set as input of the code, for instance to simulate fueling of the plasma by neutrals, or simulate a neutral background, as done in Ref. 16. Since neutrals are not bound to the magnetic field, they may escape the flux tube. This is modeled by an ad hoc sink term,^13,15 characterized by the characteristic neutral residency time τ_N, defined here as a global, constant (in space and time) quantity. In the

(θ, Ψ, φ)

coordinate system,

- \frac{n_{n}}{τ_{n}}

would include the

Ψ

contribution to

\nabla \cdot (n_{n} V_{n})

⁠, such that

- \frac{n_{n}}{τ_{n}} = - \frac{\partial}{\partial Ψ} (n_{n} V_{n}^{Ψ})

⁠.

We now write

V_{n}

in terms of components,

V_{n} (s) = (V_{n}^{θ} (s), V_{n}^{φ} (s))

⁠. After developing the divergence, Eq. can then be rewritten as

\frac{\partial n_{n}}{\partial t} + B \frac{\partial}{\partial s} (\frac{n_{n} V_{n}^{θ}}{B sin α}) = S_{n}^{n} + H_{n} - \frac{n_{n}}{τ_{n}} .

(28)

To solve this equation, a description of

V_{n} (s) = (V_{n}^{θ} (s), V_{n}^{φ} (s))

is required. This is presented in the subsequent sections.

2. Diffusive neutral model

In the diffusive neutral model, we assume that the neutral flux can be written as

n_{n} V_{n} = - D_{n} \nabla n_{n},

(29)

leading to

V_{n}^{θ} = - \frac{D_{n}}{n_{n} sin α} \frac{\partial n_{n}}{\partial s},

(30)

V_{n}^{φ} = 0,

(31)

where D_n is a diffusion coefficient, defined as

D_{n} = \frac{T_{n}}{m_{n} (n {⟨ σ v ⟩}_{c x} + n {⟨ σ v ⟩}_{ion})},

(32)

where T_n is the neutral temperature, assumed, for this diffusive neutral model, to either match the ion temperature, T_n = T_i, or to be a constant (with a value specified by the user).

{⟨ σ v ⟩}_{c x}

is the local charge-exchange reaction rate and

{⟨ σ v ⟩}_{ion}

the local ionization reaction rate, that will be introduced in Sec. II D 1. We remark here that ion-neutral elastic collisions are not considered, as we assume that the neutral and ion populations interact only through atomic processes (ionization, recombination, charge-exchange). We then obtain a diffusion equation for the neutral density n_n,

\frac{\partial n_{n}}{\partial t} = B \frac{\partial}{\partial s} (\frac{D_{n}}{B {sin}^{2} α} \frac{\partial n_{n}}{\partial s}) + S_{n}^{n} + H_{n}^{n} - \frac{n_{n}}{τ_{n}} .

(33)

This neutral model is similar to the one implemented in Refs. 11 and 16, although we retain here a dependence on T_n, whereas Ref. 11 assumes T_n = T_e and B constant along the flux tube.

D. Atomic source terms

Now, we introduce the source terms resulting from atomic interactions between the neutral and plasma populations. We restrict ourselves here to the atomic processes of ionization, charge exchange, excitation, and recombination since we model a single-ion population and consider only atomic neutrals. We remark here that ion-neutral elastic collisions are not considered, as we assume that the neutral and ion populations interact only through these atomic processes (ionization, recombination, charge-exchange). Inclusion of ion-neutral elastic collisions is left for future work.

1. Rates

SPLEND1D's rate coefficients are either obtained from a bi-linear interpolation of the open-ADAS³⁴ database, or evaluated using the fit provided in the AMJUEL manual,³⁵ as decided by the user. Table I presents the sources used for the ionization, excitation, recombination and charge exchange rates as functions of n_e and T_e. The charge-exchange reaction is discussed further below. In the case of AMJUEL data, the excitation rate Q_exc includes both the power radiated by excitation and the potential energy cost E_ion in the case of ionization (⁠ $E_{ion} = 13.6 eV$ in the case of hydrogen).³⁵ In the case of open-ADAS data, the latter is added through an additional sink term proportional to E_ion and the ionization rate.

TABLE I.

Summary of the rate coefficients used in SPLEND1D. The charge-exchange reaction is discussed further in Sec. II D 1.

Reaction	Rate coefficient	AMJUEL³⁵	Open-ADAS³⁴
Ionization	${⟨ σ v ⟩}_{ion}$	Reaction 2.1.5, section 4	SCD
Excitation	Q_exc	Reaction 2.1.5, section 10.2	PLT
Recombination	${⟨ σ v ⟩}_{rec}$	2.1.8, section 4	ACD
Recombination cooling rate	Q_rec	2.1.8, section 10.4	PRB
Charge exchange	${⟨ σ v ⟩}_{C X}$	2.1.9, section 3.1.8	CCD

Reaction	Rate coefficient	AMJUEL³⁵	Open-ADAS³⁴
Ionization	${⟨ σ v ⟩}_{ion}$	Reaction 2.1.5, section 4	SCD
Excitation	Q_exc	Reaction 2.1.5, section 10.2	PLT
Recombination	${⟨ σ v ⟩}_{rec}$	2.1.8, section 4	ACD
Recombination cooling rate	Q_rec	2.1.8, section 10.4	PRB
Charge exchange	${⟨ σ v ⟩}_{C X}$	2.1.9, section 3.1.8	CCD

The charge-exchange rate coefficient

{⟨ σ v ⟩}_{C X}

⁠, when evaluated using the AMJUEL fit, is a function of an effective temperature defined, following the AMJUEL reference manual, as

T_{eff} = \frac{m_{i}}{m_{H}} (T_{i} + β T_{n}),

(34)

where m_i is the ion mass used in SPLEND1D,

m_{H}

the hydrogen (proton) mass, and β is a user flag, set equal to either 0 or 1, determining whether or not to include the finite neutral temperature in the calculation of the effective temperature. Note that in Eq. , we used the assumption m_n = m_i. In the case of the open-ADAS CCD coefficient, no rescaling of T_eff is performed, and the rate is computed from T_e. It is then the responsibility of the user to ensure that the provided rate is indeed adequate.

Since these rates are estimated from either tabulated data (open-ADAS) or from fitted expressions, they are only defined over a certain range of validity. If, during a SPLEND1D simulation, the temperature or density were to go outside of these ranges, the rates are computed using the minimum (or maximum) values for which they are defined.

2. Sources and sinks

The plasma particle source resulting from these atomic reactions is given by

S_{p}^{n} = \underset{ionization}{\underset{︸}{n n_{n} {⟨ σ v ⟩}_{ion}}} \underset{recombination}{\underset{︸}{- n^{2} {⟨ σ v ⟩}_{rec}}} .

(35)

For the neutrals, we have

S_{n}^{n} = - S_{p}^{n}

⁠.

For the energy equation, the source term is separated between the electron and ion contributions. For the electrons, one has^14,27

\begin{matrix} S_{e}^{E} = \underline{\underset{ionization}{\underset{︸}{- n n_{n} {⟨ σ v ⟩}_{ion} E_{ion}}}} \underset{excitation}{\underset{︸}{- n n_{n} Q_{exc}}} \underset{recombination}{\underset{︸}{+ n^{2} (E_{ion} {⟨ σ v ⟩}_{rec} - Q_{rec})}} . \end{matrix}

(36)

The ionization term in Eq. is added only when using open-ADAS coefficients. As mentioned above, it is already included in the Q_exc terms derived from AMJUEL. At this point, it is important to justify the recombination contribution to the electron energy balance, Eq. . The term is split between two contributions,

n^{2} E_{ion} {⟨ σ v ⟩}_{rec}

⁠, that releases the ionization potential energy back to the electrons, and the

- n^{2} Q_{rec}

term, which encompasses the radiative energy losses during recombination, as well as some further Bremsstrahlung losses if Q_rec is taken from the open-ADAS PRB coefficient.³⁴ This is discussed in Ref. 1, Chapter 3, and in Ref. 31, Sec. 4.3. In particular, it was found that in typical TCV conditions, these two contributions approximately balance each other.⁵ This may, however, not be the case when three-body recombination dominates. We remark here that this implementation is consistent with that of EMC3-EIRENE³⁶ and SolEdge-2D (fluid neutrals).³⁷ We further remark that this term can disabled by the user in SPLEND1D. It will however be included in the simulations presented in this paper, although its contribution to the energy balance of the simulations is marginal.

For the ions, we have

\begin{matrix} S_{i}^{E} = \underset{ionization}{\underset{︸}{n n_{n} {⟨ σ v ⟩}_{ion} [\frac{3}{2} T_{n} + \frac{1}{2} m_{n} V_{n}^{2}]}} \\ \underset{recombination}{\underset{︸}{- n^{2} {⟨ σ v ⟩}_{rec} [\frac{3}{2} T_{i} + \frac{1}{2} m_{i} u_{| |}^{2}]}} \\ \underset{charge - exchange}{\underset{︸}{+ n_{n} n {⟨ σ v ⟩}_{C X} [\frac{3}{2} (T_{n} - T_{i}) + \frac{1}{2} (m_{n} V_{n}^{2} - m_{i} u_{| |}^{2})]}} . \end{matrix}

(37)

For the neutral energy equation, we have

S_{n}^{E} = - S_{i}^{E}

⁠, so that the energy that is lost, or gained, by the ions during atomic reactions is transferred to, or from, the neutrals.

Regarding the source terms in the plasma momentum equation, we have

\begin{matrix} S_{p}^{u} = \underset{ionization}{\underset{︸}{m_{i} n n_{n} {⟨ σ v ⟩}_{ion} V_{n}}} \\ \underset{recombination}{\underset{︸}{- m_{i} n^{2} {⟨ σ v ⟩}_{rec} V}} \\ \underset{charge - exchange}{\underset{︸}{+ m_{i} n_{n} n {⟨ σ v ⟩}_{C X} (V_{n} - V)}} . \end{matrix}

(38)

We project Eq. along

b = (sin α, cos α)

and find

\begin{matrix} S_{p, | |}^{u} = \underset{ionization}{\underset{︸}{m_{i} n n_{n} {⟨ σ v ⟩}_{ion} (sin α V_{n}^{θ} + cos α V_{n}^{φ})}} \\ \underset{recombination}{\underset{︸}{- m_{i} n^{2} {⟨ σ v ⟩}_{rec} u_{| |}}} \\ \underset{charge - exchange}{\underset{︸}{+ m_{i} n_{n} n {⟨ σ v ⟩}_{C X} ((sin α V_{n}^{θ} + cos α V_{n}^{φ}) - u_{| |})}} . \end{matrix}

(39)

This implies that, within our model, when a neutral ionizes, the resulting ion inherits only the parallel component of the neutral momentum. The perpendicular component of the neutral momentum is lost, leading to an effective increase in the resulting ion thermal energy, because the parallel momentum of the neutral is transferred to the ion, as well as its total energy. Therefore, the perpendicular energy

\frac{1}{2} m_{n} n_{n} {(V_{n}^{⊥})}^{2}

is redistributed as thermal energy.

For the neutral momentum source terms, we project

S_{n}^{u} = - S_{p}^{u}

onto the (⁠

θ, φ

⁠) basis. This yields

{\begin{matrix} S_{n, θ}^{u} = \underset{ionization}{\underset{︸}{- m_{i} n n_{n} {⟨ σ v ⟩}_{ion} V_{n}^{θ}}} \underset{recombination}{\underset{︸}{+ m_{i} n^{2} {⟨ σ v ⟩}_{rec} u_{| |} sin α}} \\ \underset{charge - exchange}{\underset{︸}{- m_{i} n_{n} n {⟨ σ v ⟩}_{C X} (V_{n}^{θ} - u_{| |} sin α)}}, \\ S_{n, φ}^{u} = \underset{ionization}{\underset{︸}{- m_{i} n n_{n} {⟨ σ v ⟩}_{ion} V_{n}^{φ}}} \underset{recombination}{\underset{︸}{+ m_{i} n^{2} {⟨ σ v ⟩}_{rec} u_{| |} cos α}} \\ \underset{charge - exchange}{\underset{︸}{- m_{i} n_{n} n {⟨ σ v ⟩}_{C X} (V_{n}^{φ} - u_{| |} cos α)}} . \end{matrix}

(40)

The charge exchange reaction introduces a friction term that tends to align the neutral velocity with the magnetic field. In the case of the advective neutral model, it acts as a sink for their perpendicular momentum.

E. Impurity radiation

Impurities are considered to only contribute to volumetric power loss, where all loss mechanisms are grouped into a single term,

S_{imp}^{E}

⁠. We do not consider any dilution of the main plasma ions by the impurity species. For simplicity, we assume that for a given impurity imp, one has

n_{imp} = f_{imp} n

where f_imp is the impurity fraction and n_imp the impurity density. f_imp is assumed to be constant over the flux tube, and therefore we do not self-consistently model the distribution of the impurity density. Furthermore, we assume coronal equilibrium, and the power loss due the impurities is written as

S_{imp}^{E} = - \sum_{imp} f_{imp} n^{2} L_{z, imp} (T_{e}),

(41)

where

L_{z, imp}

is the cooling-rate of the impurity imp, either taken from Ref. 38 or pre-computed using a collisional radiative model (CRM)³⁹ that employs open-ADAS. In this paper, we use the latter. If the electron temperature is outside of the range considered in these two sources, the electron temperature used in the evaluation of

L_{z, imp} (T_{e})

is clamped either to the minimum or maximum value of T_e over which

L_{z, imp} (T_{e})

is defined.

We note here that such impurity model, limited to constant impurity fraction, does not account for the impact of impurity transport and the subsequent occurrence of non-coronal effects. These could be implemented through an ad hoc (free) parameter, τ, a residence time that would allow to form the non-coronal parameter $n_{e} τ$ as done in Ref. 40, together with tabulated values of $L_{z} (T_{e}, n_{e} τ)$ produced with the collisional radiative model (CRM),³⁹ or with the inclusion within SPLEND1D itself of a collisional radiative model for impurities, evolving in time the concentration of impurities' charge states, as done for instance in Ref. 41.

F. Boundary conditions

To solve the equations describing the evolution of the plasma parameters [Eqs. (2)–(5) or (10)–(12)] and the evolution of the fluid neutral parameters [Eqs. (28), (C3), (C4), and (C10)], a set of boundary conditions is required. We distinguish two types of boundaries, depending on whether the considered boundary is a symmetry plane or a “target.” Two kinds of configurations can be simulated with SPLEND1D (Fig. 2):

The flux-tube is assumed to be symmetric around s = 0. Symmetric boundary conditions are applied on the left boundary, while “target” boundary conditions are applied to the right boundary [Fig. 2(a)].
The flux tube connects two targets together, and no symmetry is assumed. In this case, “target” boundary conditions are applied to the left and right boundaries [Fig. 2(b)].

FIG. 2.

View large Download slide

Geometries simulated by SPLEND1D, drawn here for the situation $α = π / 2$ ⁠. In the first case (panel a), the flux-tube is assumed to be symmetric around s = 0. Symmetric boundary conditions are applied at the left boundary, while “target” boundary conditions are applied at the right boundary. In the second case (panel b), the flux tube connects two targets together, and there is not necessarily a symmetry plane. In this case, “target” boundary conditions are applied to the left and right boundaries. In both cases, s is oriented so that it increases toward the right.

The second case can be used, for instance, to study the power sharing between the two sides of the flux tube, which depends on the ratio of the connection lengths.⁴²

In the case of a symmetry plane, at s = 0, we impose the derivative of the scalar fields to be zero

\frac{\partial n}{\partial s} = \frac{\partial n_{n}}{\partial s} = 0, \frac{\partial T_{e}}{\partial s} = \frac{\partial T_{i}}{\partial s} = \frac{\partial T_{n}}{\partial s} = 0,

(42)

while no-flow conditions are imposed for the velocities, such that

u_{| |} = 0, V_{n}^{θ} = V_{n}^{φ} = 0.

(43)

In the case of a “target” boundary condition, we apply the Bohm criterion on the plasma velocity, imposing to either match or exceed the sound speed c_s, as decided by the user. The sound speed is defined as the isothermal sound speed, that is¹^,

c_{s} = \sqrt{\frac{(T_{e} + T_{i})}{m_{i}}},

(44)

with the plus sign if the boundary is at the right side of the domain, and the minus sign if the boundary is at the left side. We remark here that expressing c_s as the isothermal sound speed is an assumption of the model, as various expressions of c_s can be used, depending on the sheath model retained (see Refs. 1 and 43 for a discussion of this topic). The isothermal sound speed has been retained in SPLEND1D as this is the common choice in the community (for instance in the SOLPS-ITER code⁴⁴). In addition to constraining the velocity at the sheath entrance, we also enforce boundary conditions on the heat flux through the sheath, such that

{\begin{matrix} q_{| |, t}^{i} = \frac{5}{2} n T_{i} u_{| |} + \frac{1}{2} m_{i} n u_{| |}^{3} + q_{| |, i}^{cond} + q_{| |, i}^{visc} = γ_{i} n T_{i} u_{| |}, \\ q_{| |, t}^{e} = \frac{5}{2} n T_{e} u_{| |} + q_{| |, e}^{cond} = γ_{e} n T_{e} u_{| |}, \end{matrix}

(45)

where

q_{| |, t}^{i}

and

q_{| |, t}^{e}

are the target ion and electron heat fluxes, and

q_{| |, i}^{visc}

is the viscous contribution to the heat flux. γ_i and γ_e are the sheath transmission coefficients, typically chosen as

γ_{i} = 3.5

and

γ_{e} = 5.5

⁠, see Chapter 2 of Ref. 1. For simplicity, in this paper, we do not retain any dependency of γ_i and γ_e with plasma parameters, although they could in principle depend on local plasma parameters and the Mach number.¹ We also further remark here that the boundary condition is applied on the total heat flux, which includes the viscous contribution

q_{| |, i}^{visc}

⁠. In other codes, this term is often either neglected, or explicitly not included in the heat flux boundary condition. Testing the two formulations with SPLEND1D, we find that, in the conditions of the base case that will be reported in Sec. III, such choice has negligible impact on the simulation outcomes, showing a modest effect only in strongly attached situations.

Finally, for the neutrals, we impose a recycling boundary condition on the neutral flux

Γ_{n} = n_{n} v_{n}

⁠, such that

Γ_{n} \cdot n = - R Γ \cdot n,

(46)

where

Γ

is the ion flux to the target, R is the recycling coefficient, and n the normal vector to the target. In the case of the advective model, the boundary conditions for the neutral velocities and temperature are chosen as

{\begin{matrix} V_{n}^{θ} = \mp \sqrt{\frac{E_{0}}{m_{i}}}, \\ V_{n}^{φ} = 0, \\ T_{n} = \frac{1}{3} E_{0}, \end{matrix}

(47)

where E₀ is the Franck-Condon energy, taken as

E_{0} = 3 eV

⁠.⁴⁵ In the case of the diffusive model,

V_{n}^{θ}

and

V_{n}^{φ}

are set according to Eqs. and , respectively.

G. Numerical implementation

Implemented in Fortran, SPLEND1D employs a finite volume method, solving the plasma and neutrals equations in their conservative forms. The grid is typically chosen to be non-uniform and accumulate toward the boundaries of the computational domain where sheath boundary conditions are enforced, to account for the strong gradients that can develop there. For a symmetric case, the width h_i of a given cell i is typically defined as

h_{i} = q h_{i - 1},

(48)

where

0 < q \leq 1

⁠. A third-order CWENO3 reconstruction,⁴⁶ associated with Rusanov numerical fluxes, is used for the advection terms. A low-order reconstruction, where quantities are assumed constant in a cell, is also implemented. While this further speeds-up the code, this comes at the price of reduced accuracy and will not be discussed in this paper. Ghost cells are used to enforce the boundary conditions. Source terms are typically integrated over a cell using the Simpson's rule.

The code is mostly intended to be used as an initial value problem (IVP) solver. Starting from an arbitrary initial solution, the code evolves the system of equations in time, and converges toward a steady-state, if such a solution exists. Convergence is determined by the user based on the time evolution of various macroscopic quantities, as well as by the norm of the equations' residuals. The code can also directly search for a steady-state by setting the time-derivatives to 0, through a non-linear Newton solver. This is typically run only after a temporal evolution of the equations, to provide the solver a reasonable first guess for the solution.

Several time-stepping schemes are implemented. SPLEND1D can either rely on the time-stepping algorithms implemented in the TS environment of the Portable, Extensible Toolkit for Scientific Computation (PETSc) library,^47,48 using by default the fully implicit Crank-Nicolson method. In schemes that requires the computation of the Jacobian of the system, such as the Crank-Nicolson method, this operation is performed by PETSc using finite-differences with coloring. The Jacobian is typically recomputed every 10–40 iterations, a number that is set at run-time by the user. The time step is free to evolve between user-prescribed minimum and maximum values, with PETSc taking care of the time step adaptivity. Alternatively to the use of PETSc, and not demonstrated in this paper, SPLEND1D is equipped with a “linearized” IMEX (IMplicit-EXplicit) scheme, where the advection and source terms are treated explicitly, and the diffusion and viscosity operators are handled implicitly, using the nonlinear transport coefficients of the previous time step. The implicit part can then be rewritten as a succession of tridiagonal matrix inversions, performed using the Thomas algorithm, which scales linearly with the number of cells. While this scheme comes with a Courant-Friedrichs-Lewy (CFL) condition that restricts the time step, its linear scaling with the resolution and number of equations makes it an interesting solver for cases with a large number of cells (provided the cells are not so small that the CFL condition becomes very strict). Due to the relatively small size of the typical problem solved by SPLEND1D (around ∼500–1000 cells for a maximum of 8 equations, that is, ∼4000–8000 degrees of freedom), the code is currently not parallelized, although extension to an Message Passing Interface (MPI)-OpenMP parallelized code would be relatively straightforward, leveraging the capabilities of PETSc for the MPI parallelization. In Sec. III, after presenting a base case used to demonstrate various features of the code, we will briefly present the performance of SPLEND1D in terms of convergence and speed.

III. PRESENTATION OF BASE CASE

A. Parameters

This section presents the base case investigated in this paper to illustrate in more detail the capabilities of the code, to highlight the role of various modeling parameters and to demonstrate the code's capabilities to unravel the physics behind plasma detachment. The simplified model given by Eqs. (10)–(12) is used, assuming T_i = T_e as opposed to the two-fluid model. The diffusive neutral model is employed [Eqs. (33) and (32)], and we impose T_n = T_i. Symmetric boundary conditions are applied, and the velocity is imposed greater or equal than the sound speed at the right boundary. Table II summarizes the values of the various parameters in these simulations. R, γ_i, γ_e, and τ_n are based on reasonable values, whereas $L_{| |}$ ⁠, c_c, and α are chosen based on typical TCV¹⁷ values. 500 grid cells are simulated, with the grid accumulation chosen so that the width of the cell in front of the target is approximately 0.8 mm in the parallel direction. These choices will be further discussed in Sec. III B.

TABLE II.

Summary of the parameters used for the base case simulations.

Parameter	Definition	Value
$ln Λ$	Coulomb logarithm	Equation (A4)
$\bar{τ}$	$T_{i} = \bar{τ} T_{e}$	1
γ_i	Ion sheath transmission coefficient	3.5
γ_e	Electron sheath transmission coefficient	5.5
α_i	Ion heat flux limiter	0.6
α_e	Electron heat flux limiter	0.6
c_c	Carbon concentration	2%
R	Recycling rate	99%
τ_n	Neutral confinement time	$0.05 ms$
B	Magnetic field	Constant
α	Field-line angle	$4^{°}$
$L_{\| \|}$	Parallel connection length	$25 m$
ν_num	Numerical viscosity	0

Parameter	Definition	Value
$ln Λ$	Coulomb logarithm	Equation (A4)
$\bar{τ}$	$T_{i} = \bar{τ} T_{e}$	1
γ_i	Ion sheath transmission coefficient	3.5
γ_e	Electron sheath transmission coefficient	5.5
α_i	Ion heat flux limiter	0.6
α_e	Electron heat flux limiter	0.6
c_c	Carbon concentration	2%
R	Recycling rate	99%
τ_n	Neutral confinement time	$0.05 ms$
B	Magnetic field	Constant
α	Field-line angle	$4^{°}$
$L_{\| \|}$	Parallel connection length	$25 m$
ν_num	Numerical viscosity	0

Since the system is source-driven, volumetric particle and energy source terms are needed. The energy sources, both for ions and electrons, are defined as Gaussian sources peaked at s = 0 with a standard deviation (characteristic width)

σ = 1.76 m

and amplitudes

{\bar{H}}_{i} = {\bar{H}}_{e} = 2.80 {MWm}^{- 3}

⁠, such that

H_{e} (s) = H_{i} (s) = {\bar{H}}_{i} exp (\frac{- s^{2}}{2 σ^{2}}) .

(49)

{\bar{H}}_{i}

and

{\bar{H}}_{e}

have been chosen so that, in attached conditions, the upstream temperature is approximately 40 eV, which is typical in TCV L-Modes. Considering typical TCV parameters, the obtained value is compatible with that found in moderate plasma current TCV discharges.

The neutral particles' source is defined as a constant source along the flux tube, such that

H_{n} (s) = {\bar{H}}_{n} = 1.58 \times 10^{21} m^{- 3} s^{- 1} .

(50)

The charged particles' source, H_P, is set to 0.

This results in an upstream density of $1.5 \times 10^{19}$ m⁻³ and upstream temperature of $\approx 40$ eV. The full temperature and density profiles are shown in Fig. 3 in green, alongside the velocity and Mach number profiles. A target temperature of 11.6 eV suggests that the base case represents an attached plasma regime. In this scenario, increasing the upstream plasma density via an increase in the particle source H_n leads to a rollover in the target ion flux $Γ_{| |, t}$ ⁠, Fig. 4, indicating the onset of detachment. Example profiles for detached and strongly detached cases are shown in orange and red in Fig. 3.

FIG. 3.

View large Download slide

(a) Plasma density, (b) velocity and Mach number, and (c) temperature profiles along the SOL length, s, for the base case—in attached conditions (green), for a detached case (orange), and for a strongly detached case (red). At the target, the Mach number M is M = 1 for the attached case, $M \approx 1.7$ for the detached and strongly detached cases.

FIG. 4.

View large Download slide

Target ion saturation current as a function of upstream density n_u, with the attached, detached, and strongly detached base cases of Fig. 3 shown by the green, orange, and red dots respectively.

We note that in Fig. 4, each point is the steady-state result of simulations with different values of H_n. They were obtained by simulating $1 s$ of plasma dynamics, and then applying the steady-state solver discussed in Sec. II G. This will be the case in all simulation results presented in this paper, except the time-dependent simulation presented in Sec. VI B.

B. SPLEND1D performance in the base case

In this section, we aim to quantify the performance of SPLEND1D in terms of accuracy in the base-case scenario described previously. In particular, we will use the methodology used in Ref. 16, although one should note that the results are not directly comparable as the considered cases differ. Integrating Eq. from upstream (s = 0, subscript u) to the target (⁠

s = L_{| |}

⁠, subscript t) and neglecting the time derivative, one can define the numerical error of the particle balance, ϵ_part, as

ϵ_{part} = | 1 - \frac{{[\frac{n u_{| |}}{B}]}_{u}^{t}}{\int_{u}^{t} (\frac{S_{p}^{n} + H_{P}}{B}) d s} | .

(51)

Since this paper focuses mainly on steady-state simulations, we will focus on the code performance and accuracy in such conditions. In steady-state conditions, one should ideally find

ϵ_{part} = 0

⁠. Any finite value of ϵ_part comes either from numerical errors or the numerical tolerance, as the steady-state solver has removed the time-derivatives from the system. Similarly, one can define the numerical error of the momentum balance [Eq. ], ϵ_mom, as

ϵ_{mom} = | 1 - \frac{{[\frac{m_{i} n u_{| |}^{2} + p}{B}]}_{u}^{t}}{\int_{u}^{t} \frac{- \frac{p}{B} \frac{\partial B}{\partial s} + S_{p, | |}^{u} - b \cdot \nabla \cdot Π}{B} d s} |,

(52)

and the numerical error of the power balance [Eq. ], ϵ_pow, as

ϵ_{pow} = | 1 - \frac{{[\frac{q_{| |}^{tot}}{B}]}_{u}^{t}}{\int_{u}^{t} \frac{S_{i}^{E} + S_{e}^{E} + S_{imp}^{E} + H_{e} + H_{i}}{B} d s} |,

(53)

where

q_{| |}^{tot} = \frac{5}{2} n (1 + \bar{τ}) T_{e} u_{| |} + \frac{1}{2} m_{i} n u_{| |}^{3} + q_{| |, i}^{cond} + q_{| |, e}^{cond} + q_{| |, i}^{visc} .

(54)

We also define δ_p, as

δ_{p} = | \frac{p - p_{ref}}{p_{ref}} |,

(55)

where

p = 2 n T

is the static pressure and

p_{ref}

the static pressure of a reference simulation. δ_p is defined either upstream (s = 0) or at the target (⁠

s = L_{| |}

⁠).

δ p

quantifies the variation of the static pressure in each simulation to that of a reference simulation. This will be used later in this section to assess the effect of target-cell width and number of cells on the numerical convergence of SPLEND1D. We now evaluate these quantities for the density ramp shown in Fig. 4. Figure 5 shows the numerical error on the particle, momentum and energy balances, ϵ_part, ϵ_mom, ϵ_pow. Two regimes can be identified. For

n_{u} > 1.5 \times 10^{19} m^{- 3}

⁠, SPLEND1D shows excellent convergence properties even in the absence of numerical viscosity ν_num [see Eq. ], with a maximum error on the particle balance of

ϵ_{part}^{max} \approx 0.1 %

⁠, while it is

ϵ_{mom}^{max} \approx 0.0015 %

for the momentum and

ϵ_{pow}^{max} \approx 0.0023 %

for the power balance. For

n_{u} < 1.5 \times 10^{19} m^{- 3}

⁠, the situation is more intricate, and while particle and power balances remain satisfactory, the error on the momentum balance can become important (up to

\approx 30 %

⁠) in the absence of numerical viscosity ν_num. This is due to the very strong velocity gradient that will form just in front of the target to bring the flow from a very low value to the sound speed. This can be alleviated by the addition of a finite ν_num, bringing the momentum error below 1%, Fig. 5(d). This has little impact on the overall outputs of the simulation, Fig. 5(a), where

Γ_{| |, t}

remains virtually unchanged across the various values of ν_num. Similarly, the upstream and target pressure are only affected by up to 3% by the addition of ν_num, Fig. 5(b). Another possibility would be to increase the resolution of the grid near the target, for instance by reducing the width of the cells, as will be shown later in this section. Overall, these results demonstrate that the numerical accuracy of SPLEND1D across the different regimes is satisfactory, from attached to detached regimes.

FIG. 5.

View large Download slide

(a) Target particle flux for different values of ν_num, as a function of the upstream density n_u. (b) Effect of ν_num on the upstream (⁠ $δ_{p, u}$ ⁠) and target (⁠ $δ_{p, t}$ ⁠) δ_p indicator for different values of ν_num, using $ν_{num} = 0$ as a reference, as a function of the upstream density n_u. (c)–(e) Effect of ν_num on the particle balance (ϵ_part), momentum balance (ϵ_mom), and energy balance (ϵ_pow), as a function of the upstream density n_u.

We next explore the performance, in term of computing time, of SPLEND1D, across a density ramp with the same input parameters as that of Fig. 4. For this, we evaluate the simulation time, toward steady-state, for different values of ${\bar{H}}_{n}$ ⁠, applying three different strategies:

Each simulation is run independently, starting from flat, unphysical profiles (⁠ $n = 2.0 \times 10^{19} m^{- 3}, T_{e} = T_{i} = 20 eV, u_{| |} = 0 {ms}^{- 1}$ ⁠).
Each simulation is run independently, starting from “physical” profiles obtained from a steady-state simulation with ${\bar{H}}_{n} = 1.58 \times 10^{21} m^{- 3} s^{- 1}$ ⁠.
The simulations are run sequentially, each starting from the steady-state obtained from the previous simulations.

In all simulations presented in this section, the same convergence parameters for the PETSc solver are used. The time step $Δ t$ is allowed to vary between $Δ t = 5.71 \times 10^{- 11} s$ and $Δ t = 5.71 \times 10^{- 4} s$ ⁠. It is automatically adapted by PETSc based on convergence and error estimates. The simulations are run for $1 s$ of simulated plasma time, after which the steady-state solver is applied to find an exact (within nonlinear solver tolerance) steady-state solution. All simulations were performed on a typical laptop, with an Intel^® Core™i7-8565U CPU. These three simulation strategies provide identical output profiles (within the nonlinear solver tolerance), Fig. 6(a). Figure 6(b) reports the simulation run-times. Starting from flat, unphysical profiles (strategy 1), all simulations converge within $\approx$ 20 s, except for a few outliers that require up to 40 s. Starting from physical profiles of an attached case (strategy 2) yields even lower computation times, approximately tens, except at low densities (lower than the initial simulation), where simulations can take relatively long times or even fail to converge in less than 10 min (after which the simulations were stopped). A similar observation is done for the sequential scan of ${\bar{H}}_{n}$ (strategy 3), restarting from the previous simulation (with a slightly different value of ${\bar{H}}_{n}$ ⁠), where the simulation time can drop to $\approx$ 6 s, except at low density where simulations can struggle to converge in reasonable time. In conclusion, this section demonstrates the numerical performances of SPLEND1D, which can achieve convergence toward a steady-state in less than 30 s across a wide range of regimes, including the strongly detached one.

FIG. 6.

View large Download slide

Results of simulations using the three strategies outlined in Sec. III B. (a) Upstream (solid) and target (dashed) temperature as a function of upstream density n_u, with indistinguishable results across the applied strategies. (b) Simulation run-time required for $1 s$ of simulated plasma time.

We now inspect the dependence of SPLEND1D results and accuracy on the grid resolution, for simplicity in the case of a symmetric domain, as done in the base case. The resolution of the grid is controlled by two parameters, the total number of cells, and the width of the last cell before the target, thereafter denoted h_end, which define entirely the grid, Eq. . In the following, we will assess the effect of both these parameters in the attached base case, varying only one parameter at a time. For simplicity, we group the error on particle balance, momentum balance, and power balance under a single new term, ϵ_tot, defined as

ϵ_{tot} = \sqrt{ϵ_{part}^{2} + ϵ_{mom}^{2} + ϵ_{pow}^{2}} .

(56)

We start by considering the influence of the number of grid cells, N, keeping $h_{end} = 0.8 mm$ ⁠. N is varied from 100 to 1000 by increments of 50, and then from 1000 to 10 000 by increments of 500, Fig. 7. All simulations converged. We first use a simulation with 31 250 cells of 0.8 mm each as a reference for the evaluation of p_ref. As N is increased, we observe that both the upstream and downstream values of δ_p decrease quadratically with N, as one could expect from a second order scheme. ϵ_tot remains largely unaffected. Furthermore, using now a high resolution (125 000 cells, $0.2 mm$ each) simulation as reference, we find that δ_p stagnates when the number of cells is higher than 2000. This is because, as will be shown in the next paragraph, the accuracy of SPLEND1D is largely dictated by the width of the cells close to the target, kept fixed in this scan. Hence, increasing further the number of cells does not lead to an increase in the accuracy. Furthermore, and interestingly, even at fairly low resolution (⁠ $\sim 100$ cells), δ_p for the upstream and target pressure remain lower than 1%.

FIG. 7.

View large Download slide

(Top) δ_p for the upstream pressure (blue) and downstream pressure (red), using either the 31 250 cells case as reference (squares) or the 1 25 000 cells as references (crosses). The scaling relations of δ_p with N are computed using the 31 250 cells case as reference; (bottom) error of the particles, momentum, and energy balances ϵ_tot as a function of the number of grid cells, for a fixed width of the last cell, $h_{end} = 0.8 mm$ ⁠.

We now keep the number of grid cells constant (500 cells) while changing h_end. As h_end is reduced, ϵ_tot and δ_p (both upstream and downstream) strongly decrease, Fig. 8. This indicates that the accuracy of SPLEND1D largely depends on the width of cells near the target. In particular, ϵ_tot scales approximately as $h_{end}^{2}$ ⁠.

FIG. 8.

View large Download slide

(Top) δ_p for the upstream pressure (blue) and downstream pressure (red), using either the 500 cells case as reference (squares) or the 125 000 cells as references (crosses); (bottom) error of the particles, momentum, and energy balances ϵ_tot as a function of the width of the last cell, h_end, for a fixed number of grid cells.

Taken together, the results of Figs. 7 and 8 show that the accuracy of SPLEND1D is largely dictated by the grid resolution near the targets, with the total number of grid cells playing a lesser role. However, as the number of grid cells is increased, or the width of the last cell is decreased, the computational cost increases. From Figs. 6–8, the resolution chosen for the base case (500 grid cells, 0.8 mm target cell width) appears to be a good trade-off between accuracy and numerical cost, and is retained for the rest of the simulations presented in this paper.

C. The two-point model formatting

To evaluate the importance of physical momentum and power loss processes in the base case simulations, we first introduce briefly the two-point model (TPM) formulation. It is a model that does not consider the spatial distribution of plasma parameters along the SOL, but considers quantities at only two points: upstream and target, where upstream can be any point along the flux tube. Furthermore, it considers a steady state situation.

Momentum losses between upstream and target are grouped into a single momentum loss factor, defined as the relative difference in upstream and target total pressures. We remark here that some works¹ define the momentum loss factor as

(1 - f_{mom})

⁠, being simply the ratio of target to upstream pressure

f_{mom} = \frac{p_{tot, u} - p_{tot, t}}{p_{tot, u}},

(57)

= 1 - \frac{n_{t} T_{e, t}}{n_{u} T_{e, u}} \frac{(1 + \frac{T_{i, t}}{T_{e, t}})}{(1 + \frac{T_{i, u}}{T_{e, u}})} \frac{(1 + M_{t}^{2})}{(1 + M_{u}^{2})},

(58)

where M_t and M_u are the target and upstream Mach numbers, respectively:

M = v_{| |} / c_{s}

⁠. Similarly, volumetric power losses along the SOL are described by a single power loss factor,

f_{power}

⁠. However, because the spatial location of the power sources in SPLEND1D can be chosen arbitrarily, care must be taken for the definition of

f_{power}

⁠. We start by integrating, from the target to upstream, the sum of the electron and ion energy equations, Eqs. and , taken at steady-state. Defining the heat flux as the sum of the convective, conductive, and viscous contributions

q_{| |}^{tot} = \frac{5}{2} n (1 + T_{i} / T_{e}) T_{e} u_{| |} + \frac{1}{2} m_{i} n u_{| |}^{3} + q_{| |, i}^{cond} + q_{| |, e}^{cond} + q_{| |, i}^{visc},

(59)

we then have

\begin{matrix} B \frac{\partial}{\partial s} (\frac{q_{| |}^{tot}}{B}) = S_{i}^{E} + S_{e}^{E} + S_{imp}^{E} + H_{e} + H_{i} . \end{matrix}

We now define the total energy loss factor, f_pow, as

f_{pow} = 1 - \frac{q_{| |, t}^{tot} B_{u} / B_{t}}{q_{i n}},

(60)

where

q_{i n} = q_{| |, u}^{tot} + B_{u} \int_{u}^{t} \frac{(H_{e} + H_{i})}{B} d s,

(61)

and

q_{| |, t}^{tot} = q_{| |, t}^{cond, e} + q_{| |, t}^{cond, i} + q_{| |, u}^{conv, t} + q_{| |, u}^{conv, t} + q_{| |, i}^{visc, t},

(62)

= n_{t} u_{| |} (γ_{e} T_{e} + γ_{i} T_{i}) .

(63)

When applied to the SPLEND1D model described in Sec. II B 1, the contribution of each source term to the momentum and power losses can be evaluated such that the important processes can be identified. The breakdown of the momentum and power loss factors is given explicitly in Appendix B. Figure 12 shows the individual contributions to f_pow and f_mom, including atomic sources and viscosity. This will be studied in further detail in Sec. IV B. We remark here that, with this definition of f_mom, geometry effects related to total flux expansion are embedded within volumetric source and sink terms. In order to highlight more clearly the role of total flux expansion on target conditions, it can be preferable to use the f_mom definition proposed in Ref. 49, which is specifically formulated to elucidate the total flux expansion effect. We further remark that, since the viscous heat flux is included in the expression of

q_{| |}^{tot}

and enters the heat flux boundary condition [Eq. ], the contribution of viscosity is not included in f_pow.

Considering particle, momentum, and power balances with the various loss factors, the target temperature and density can be expressed as a function of upstream total pressure,

p_{tot, u}

⁠, and input heat flux q_in [as defined in Eq. ] as follows:

T_{e, t} = m_{i} \frac{{(1 + M_{t}^{2})}^{2}}{M_{t}^{2}} \frac{{(1 - f_{power})}^{2}}{{(1 - f_{mom})}^{2}} {(\frac{B_{t}}{B_{u}})}^{2},

(64)

\times \frac{(1 + \frac{T_{i, t}}{T_{e, t}})}{{(γ_{e} + γ_{i} \frac{T_{i, t}}{T_{e, t}})}^{2}} \frac{q_{i n}^{2}}{p_{tot, u}^{2}},

(65)

n_{t} = \frac{1}{m_{i}} \frac{M_{t}^{2}}{{(1 + M_{t}^{2})}^{3}} \frac{{(1 - f_{mom})}^{3}}{{(1 - f_{power})}^{2}} {(\frac{B_{u}}{B_{t}})}^{2},

(66)

\times \frac{{(γ_{e} + γ_{i} \frac{T_{i, t}}{T_{e, t}})}^{2}}{{(1 + \frac{T_{i, t}}{T_{e, t}})}^{2}} \frac{p_{tot, u}^{3}}{q_{i n}^{2}} .

(67)

This form of the two-point model, labeled as 2-point formatting,^50–52 is a reformulation of the SPLEND1D equations for the target temperature and density, given the power loss and momentum loss factors defined in Eqs. and . Such reformulation of the equations is particularly insightful as it defines global loss parameters (⁠

f_{mom}

and

f_{power}

⁠) that can be used to establish scaling laws, and that can be evaluated from experimental available measurements in a simpler manner than the local sources and sinks.

IV. DETACHMENT ONSET VIA DENSITY RAMP IN THE BASE CASE

A. Observation of a target ion flux roll-over and onset of a total pressure drop

As mentioned earlier, in the base case presented in Sec. III, increasing the upstream plasma density via an increase in neutral particle source H_n leads to a rollover in the target ion flux, Fig. 4, indicating the onset of detachment. Before the roll-over, the target ion flux increases somewhat linearly with increasing n_u, whereas the two-point model would predict $Γ_{| |, t} \propto n_{u}^{3}$ ⁠. Such behavior is reminiscent of the “power starvation” mechanisms, observed in previous simulation works¹⁵ and experiments,⁵ in which the linear evolution of the ion flux and its subsequent roll-over are driven by a reduction of the power that is available for ionization. This is clearly seen in the simulations presented in this paper, where the roll-over coincides with a drop of the ionization source, whereas recombination only plays a role in the reduction of the ion flux at a later stage, Fig. 9. The rollover is accompanied by a reduction of target temperature to less than $1 eV$ [Fig. 10(a)]. In contrast, the upstream temperature T_u is much less sensitive, and only at the highest degree of detachment does it start to degrade significantly, dropping from approximately $\approx 40$ eV to $\approx 28$ eV. The particle flux rollover is accompanied by a total pressure drop, as a signature of detachment, as shown by evaluating f_mom [Eq. (57)], which increases from 0 (no pressure drop) at low density to 0.99 at the highest density achieved in the present simulations, Fig. 10(b). The plasma power loss factor, f_pow [Eq. (60)], exhibits a similar behavior as the momentum loss with increasing density, but reaches saturation (⁠ $f_{pow} \sim 1$ ⁠) at a lower upstream density. This suggests that the onset of detachment is first driven by the increase in power losses, which precedes momentum losses. With the reduction of the total pressure, Fig. 10(b), a reduction of the static pressure at the target is also observed, Fig. 10(c), initially driven by the reduction of the target temperature, and, ultimately, also by a reduction of the target density.

FIG. 9.

View large Download slide

Evolution of the ionization source (blue) and recombination sink (red, absolute value), integrated over the domain, as a function of n_u.

FIG. 10.

View large Download slide

(a) Evolution of target temperature $T_{e, t}$ and upstream temperature $T_{e, u}$ as a function of n_u. (b) Momentum loss factor f_mom and power loss factor f_pow between upstream and target, as defined in Eqs. (57) and (60). (c) Evolution of target static pressure (electrons + ions) and target density n_t as a function of n_u.

We remark in Fig. 10(b) that, at low density, both f_mom and f_pow can be negative. For f_mom, as shown in Figs. 11 and 12, this is due to the action of the viscosity term, which can be relatively large in the initial, attached conditions when the gradient of velocity in front of the target is large. As for f_pow, this is related to the ionization process, Fig. 12, and the particular choice of neutral particle source, H_n, which injects in the system neutrals with a temperature T_n = T_i. Upon ionization, this thermal energy is transferred to the plasma. We note here that these observations are consequences of some of the choices made in the derivation of the model.

FIG. 11.

View large Download slide

Source and sink terms contribution to plasma (top) particle, (middle) momentum, and (bottom) energy equations [Eqs. (10)–(12) respectively], in (left) an attached case, (middle) a detached case, (right) a strongly detached case. In the energy balance, the “recombination” term corresponds to the loss of the ion energy while the “recombination power loss” corresponds to the loss of electron energy.

FIG. 12.

View large Download slide

Source terms contributing to plasma (top) momentum loss and (bottom) power loss factors formulated as in Eqs. (57) and (60), and Appendix B, as a function of upstream electron density n_u. The momentum loss factor includes losses due to ionization (ion), recombination (rec), charge exchange (CX), and viscosity (visc). The power loss factor includes losses due to ionization (ion), recombination (rec), charge exchange (CX), excitation (exc), and impurity radiation (imp). Negative values represent a gain in the plasma momentum or energy.

B. Investigate process at play

The first sign of detachment, Fig. 10, appears to be a target static pressure rollover, along with a target ion flux rollover and a target temperature decrease to below $1$ eV. With a further increase in upstream density, the target electron density rolls over, followed by a saturation in the power losses and then momentum losses. The underlying processes behind these features can be studied with SPLEND1D, which directly outputs each term contributing to the particle, momentum and energy balance equations [Eqs. (10)–(12)], such that the importance of each process can be compared. Figure 11 shows the profile of each of these terms along the flux tube for an attached, a detached, and a strongly detached case. This allows to compute the contribution of each of these terms to the momentum and power losses, Fig. 12. We must note that the choice of neutral model may affect the role of atomic processes in momentum and power losses. In particular, the simulations discussed here employ the diffusive neutral model with T_n = T_i, and so the power loss due to charge exchange may be underestimated. The role of τ_N is discussed in further detail in Sec. V B.

The target electron pressure rolls over as the momentum losses, which are dominated by charge exchange reactions, begin to increase significantly, resulting from a strong increase in electron density in front of the target, as well as from a drop in temperature, which favors charge-exchange reactions over ionization. This leads to the rollover of the target ion flux. The region in front of the target becomes much cooler, and a region of strong temperature gradient moves upstream, with $T_{e, t} < 1$ eV. This front movement is also seen in other atomic processes (ionization, recombination and impurity radiation, Fig. 11), and eventually in the plasma density front as strong detachment is achieved, leading to a rollover in target electron density. The impurity radiation increases as the divertor becomes cooler and denser, becoming the dominant power loss mechanism, until f_pow saturates. As the neutral density continues to increase in the cool divertor, charge exchange momentum losses continue to increase until f_mom also saturates at $f_{mom} \sim 1$ ⁠. Note that although momentum and power losses due to recombination increase with upstream density, they remain negligible compared to other atomic sources throughout the density range explored. However, recombination appears as a significant contributor of the particle balance, Fig. 11.

One of the key assumptions of the standard two-point model (TPM) that is often used as a first model to study SOL physics¹ is that the heat transport is mostly due to (electron) heat conduction. While convection can be enabled in the extended TPM through some ad hoc f_cond parameter, the TPM itself does not provide a self-consistent way to estimate the value of f_cond. In SPLEND1D, we find that convective heat transport in the current simulations is small but non-negligible, Fig. 13. The difference between the SPLEND1D results and the TPM predictions is small at low upstream density (⁠ $n_{u} = 1.6 \times 10^{19} m^{- 3}$ ⁠) but becomes stronger as upstream density is increased, where the convective heat flux becomes significant.

FIG. 13.

Parallel heat flux profiles for the base case with nu=1.5×1019 m−3 (green), for a detached case with nu=3.0×1019 m−3 (orange), and for a strongly detached case with nu=3.5×1019 m−3. The convective (dotted) and conductive (dashed) contributions are plotted alongside the total parallel heat flux (solid).

View large Download slide

Parallel heat flux profiles for the base case with $n_{u} = 1.5 \times 10^{19} m^{- 3}$ (green), for a detached case with $n_{u} = 3.0 \times 10^{19} m^{- 3}$ (orange), and for a strongly detached case with $n_{u} = 3.5 \times 10^{19} m^{- 3}$ ⁠. The convective (dotted) and conductive (dashed) contributions are plotted alongside the total parallel heat flux (solid).

Defining f_cond, the fraction of the heat-flux that is conducted compared to the total heat-flux, averaged over the flux-tube, by

f_{cond} = \frac{1}{L_{| |}} \int (q_{cond} / q_{tot}) d s,

(68)

one observes that the achievable upstream temperature

T_{e, u}

is well correlated with f_cond, Fig. 14. As the convected heat-flux increases and f_cond decreases, the upstream temperature drops, as conduction is essential for sustaining a temperature gradient. Furthermore, the increase in the convected-fraction (decrease in f_cond) is well-correlated with the movement of the ionization peak, as well as the weighted ionization location,

s_{ion}^{weighted}

⁠, defined as

s_{ion}^{weighted} = \frac{\int s (n n_{n} {⟨ σ v ⟩}_{ion}) d s}{\int (n n_{n} {⟨ σ v ⟩}_{ion}) d s} .

(69)

FIG. 14.

View large Download slide

(Blue) Upstream electron temperature $T_{e, u}$ as a function of f_cond. (red, dashed) $s_{ion}^{weighted}$ and (red, dash-dotted) peak ionization location as a function of f_cond.

This is not surprising, since ionization acts as a source term for flows, as shown by Eqs. (2) and (35).

V. EXPLORING THE ROLE OF FREE PARAMETERS

As mentioned in the introduction and in the derivation of the model, 1D codes come with strong assumptions and, as illustrated by Table II, free parameters that may affect the results of the simulations. These become especially important when attempting to use the results of such a model to explain experimental observations (although, in the case of interpretative simulations, some of them can be constrained by experimental data). It can, therefore, be important to understand how sensitive the code results are to these free parameters and assumptions. In this section, we review the impact of some of the main assumptions or parameters on the observations of Sec. IV.

A. Carbon concentration

Let us first examine the role of a parameter strongly affecting the plasma energy sink: the carbon concentration c_c. As expected, c_c has a strong impact on the target parallel particle flux [Fig. 15(a)], with the upstream density required for the $Γ_{| |, t}$ rollover decreasing strongly with c_c. This is due to the strong cooling induced by increased carbon radiation, leading to lower target temperatures for a given upstream density and temperature [Fig. 15(b)], and hence a facilitated access to detachment. We note here that the observed behavior (reduction of the peak particle flux achievable before roll-over) is consistent with the power starvation mechanism described in Sec. IV A. As carbon concentration, and therefore radiation losses, increases, less power is available for ionization, leading to the earlier roll-over of $Γ_{| |, t}$ ⁠.

FIG. 15.

View large Download slide

(Top) Target ion saturation current and (bottom) upstream and target electron temperatures as a function of n_u and for different carbon concentrations.

B. Neutral confinement time

The neutral model implemented in SPLEND1D essentially features three free parameters: the neutral confinement time τ_N [see Eq. (28)], the choice of neutral model and neutral temperature implementation. In this section, we focus on the free parameter τ_N. Figure 16 plots the effect of varying τ_N on the inferred target ion flux and integrated neutral density in the base case density ramp. This reveals that, as expected, τ_N has little effect at low density, where the neutral density is low in all cases. Only at high n_u does the effect of τ_N become significant: the neutral density integrated along the flux tube increases strongly with τ_N, associated with a strong decrease in the target ion flux. The upstream density at which the target ion flux rolls over is not strongly affected by τ_N, except for the most extreme case where τ_N has been decreased by a factor 10.

FIG. 16.

View large Download slide

(Top) Target ion saturation current and (bottom) integrated neutral density along the field line as a function of n_u, for a range of neutral confinement times τ_N.

C. Heat flux limiter

The Spitzer-Harm heat flux can overestimate the physical heat flux at low plasma collisionality, when the electron mean free path, λ_e, is large compared to the electron temperature gradient scale length, $L_{\nabla T} = {(| \nabla T_{e} / T_{e} |)}^{- 1}$ ⁠. Figure 17 shows the ratio of these scale lengths as a function of distance along the SOL, for the base case density and temperature profiles shown in Fig. 3. At the target, the electron collisionality is high enough that the classical Spitzer-Härm heat fluxes are expected to fairly accurately predict the physical heat flux, without the need for a flux limiter. However, the upstream electron mean free path exceeds the temperature gradient scale length and so the effect of enforcing a heat flux limiter should be studied.

FIG. 17.

View large Download slide

Profile of the ratio of the parallel temperature gradient scale length to the electron collisional mean free path¹ along the SOL length, s, for the base case.

The strength of the heat flux limiter can be controlled with the parameter $α_{{i, e}}$ ⁠. Target conditions are compared here for values of $α_{i} = α_{e} = [0.06, 0.3, 0.6, 6, 6 \times 10^{17}]$ ⁠, thus testing values well below and above the typical value of 0.5 used in the literature (Refs. 20 and 21). This includes $α ≫ 1$ ⁠, approaching the situation without heat flux limiters. The results are found to be little affected by this choice. The target ion flux rollover and thus the detachment threshold is weakly affected by the choice of $α_{{i, e}}$ ⁠, Fig. 18, except for the smallest value of $α_{{i, e}} = 0.06$ ⁠. For $α_{{i, e}} \geq 0.3$ ⁠, the magnitude of the target ion flux varies by approximately 10% within the range of $α_{{i, e}}$ studied, Fig. 18, and the pressure drop along the SOL to the target is negligible (not shown). The target density is largely unaffected, and the target temperature is only affected at low upstream density (⁠ $\sim 2 \times 10^{19} m^{- 3}$ ⁠), not shown. These results hold for a range of input power levels.

FIG. 18.

View large Download slide

Target ion flux rollover for the base case with varied heat flux limiter coefficients $α_{i, e}$ ⁠.

Due to the apparent insensitivity of target parameters to the heat flux limiter coefficient for the base case conditions, and for reasonable values of $α_{{i, e}} \geq 0.3$ ⁠, and over a large part of the base case density ramp, $α_{{i, e}}$ was set to 0.6 in all simulations presented in this paper, unless stated otherwise.

D. Bohm boundary condition

SPLEND1D models the plasma along a flux tube up to the entrance of the sheath (in the case of small target angles, this actually corresponds to the entrance of the magnetic pre-sheath), which acts as a perfect sink that absorbs all incoming ions. This results in the so-called Bohm boundary condition, that is, for a purely parallel flow,

u_{| |} \geq c_{s},

(70)

where c_s is the sound speed. In the base case presented so far, we indeed allowed for

u_{| |} \geq c_{s}

(⁠

M \geq 1

⁠) at the sheath entrance. However, other codes sometimes enforce the strict equality

u_{| |} = c_{s}

⁠, thus precluding the presence of supersonic flows at the target. In this section, we briefly discuss how this may influence the various target parameters.

The base case, Sec. III, leads to naturally supersonic flows (⁠ $M \approx 1.7$ ⁠) at the target after roll-over, Fig. 19(a). When enforcing M = 1 at the target in such conditions, the code showed numerical difficulties to converge with satisfactory particle and energy balances, with ϵ_tot [Eq. (56)] reaching up to $\approx 12 %$ ⁠, Fig. 19(d). This is due to the presence of a very sharp velocity gradient required to slow down the flow to M = 1. This is alleviated by the introduction of some numerical viscosity ν_num [Eq. (26)]. Figures 19(b) and 19(c) plot the evolution of various target quantities for three different values of ν_num. The effect of this artificial viscosity on $T_{e, t}$ (not shown), n_t and $Γ_{t, | |}$ is modest, but it greatly improves the code convergence, with ϵ_tot reduced to a maximum of $\approx 1.7 %$ and $\approx 0.23 %$ for $ν_{num} = 1.0 \times 10^{- 5} [a . u .]$ and $ν_{num} = 5.0 \times 10^{- 5} [a . u .]$ ⁠, respectively, Fig. 19(d). We note here that it is not yet clear whether the occurrence of such sharp velocity gradient is a consequence of the choice of parameters for the base case, or a general observation.

FIG. 19.

View large Download slide

(a) Target Mach number M, (b) target parallel ion flux $Γ_{t, | |}$ ⁠, (c) target density n_t, and (d) global accuracy on particle, momentum and energy balances ϵ_tot, plotted as a function of the upstream density n_u, for different forms of the Bohm condition and different values of the numerical viscosity ν_num.

For both target boundary conditions M > 1 and M = 1, $Γ_{t, | |}$ is very similar, Fig. 19(b). Similarly, the target temperature $T_{e, t}$ is unaffected (not shown). However, as expected, a difference arises in the target density n_t and parallel velocity $u_{| |, t}$ ⁠. The case with M = 1 features a lower $u_{| |, t}$ and a higher n_t compared to the M > 1. Since $Γ_{t, | |}$ remains similar across these different cases, we conclude that enforcing the strict equality $u_{| |} = c_{s}$ leads to a redistribution of $Γ_{t, | |}$ between its velocity and density contributions. This could have some implications in simulation codes that enforce $u_{| |} = c_{s}$ ⁠, by promoting higher density at the target. Since many atomic source and sink terms scale with n_e, or even $n_{e}^{2}$ ⁠, this will ultimately influence the particle, momentum, and energy balance of the system.

VI. ADVANCED STUDIES ENABLED BY SPLEND1D

SPLEND1D opens up a large number of possible SOL and detachment studies, such as the investigation of the effect of the total flux expansion, of the parallel connection length, of in-out power sharing, or dynamical behavior. In Secs. VI A and VI B, we present some example studies on the role of different neutral models, ion vs electron heating ratios, and heat pulses effects on the SOL plasma.

A. Independent ion and electron temperatures

SPLEND1D can either solve a single energy equation, Eq. (12), assuming a proportionality relation between the ion and electron temperatures, or two separate energy equations for the electrons and ions, Eqs. (4) and (5). The base case presented in Sec. III has so far employed the T_e = T_i assumption. In this section, we relax this constraint by enabling both temperatures to be independent. It is well known that, at low collisionality, ion and electron temperatures can be different.¹

To set up the simulations, we split the power source between the H_e and H_i terms [Eqs. (4) and (5)], with either H_e = H_i (50/50 split of the total input power between electrons and ions), $H_{e} = 3 H_{i}$ (75/25 split of the input power between electrons and ions), or $3 H_{e} = H_{i}$ (25/75 split of the input power between electrons and ions). In all cases presented in this section, the total power injected in the system is kept constant, as well as all the parameters presented in Table II. Similarly to Sec. IV, we perform density ramps by scanning $\bar{H_{n}}$ ⁠. Allowing $T_{i} \neq T_{e}$ affects the roll-over threshold, Fig. 20(a), and leads to typically higher T_i than T_e, Figs. 20(b) and 20(c). This is easily explained by the lower heat conduction coefficient of the ions, compared to the electrons, Eqs. (16) and (17), while the convective heat-flux and the equipartition term decrease this difference in transported heat flux. If the fraction of input power carried by the electrons is increased, the differences between T_i and T_e decreases, as expected, but remain significant. As density is increased, so does the collisionality and hence the equipartition term [Eq. (6)]. The difference between T_i and T_e becomes negligible near the target (where density is high and temperature low), Fig. 20(b), and reduces at the upstream location, Fig. 20(c). This is also evident when looking at the temperature profile along the flux tubes, Figs. 20(d)–20(f), which show profiles of T_i and T_e at increasing densities. Clearly, while in the attached case [Fig. 20(d)], the two temperatures are strongly different, they get closer to each other as density increases, Figs. 20(e) and 20(f).

FIG. 20.

(a) Γ||,t as a function of upstream density nu, for the base case with Ti = Te (blue), the base case with Ti≠Te and a 75/25 power sharing between ions and electrons (red), the base case with Ti≠Te and a 50/50 power sharing between ions and electrons (yellow) and the base case with Ti≠Te and a 25/75 power sharing between ions and electrons (purple). (b) Target ion (dashed) and electron (solid) temperatures for the same cases as in Panel (a). (c) Upstream ion (dashed) and electron (solid) temperatures for the same cases as in Panel (a). Panels (d)–(f) Ion (dashed) and electron (solid) temperature profiles along the flux tube for nu=1.25×1019 m−3 (d), nu=1.75×1019 m−3 (e), and nu=2.4×1019 m−3 (f).

View large Download slide

(a) $Γ_{| |, t}$ as a function of upstream density n_u, for the base case with T_i = T_e (blue), the base case with $T_{i} \neq T_{e}$ and a 75/25 power sharing between ions and electrons (red), the base case with $T_{i} \neq T_{e}$ and a 50/50 power sharing between ions and electrons (yellow) and the base case with $T_{i} \neq T_{e}$ and a 25/75 power sharing between ions and electrons (purple). (b) Target ion (dashed) and electron (solid) temperatures for the same cases as in Panel (a). (c) Upstream ion (dashed) and electron (solid) temperatures for the same cases as in Panel (a). Panels (d)–(f) Ion (dashed) and electron (solid) temperature profiles along the flux tube for $n_{u} = 1.25 \times 10^{19} m^{- 3}$ (d), $n_{u} = 1.75 \times 10^{19} m^{- 3}$ (e), and $n_{u} = 2.4 \times 10^{19} m^{- 3}$ (f).

B. Time-dependent simulations

As mentioned in Sec. III B, SPLEND1D solves the 1D Braginskii equations as a time-dependent problem. Hence, in addition to finding steady-state solutions, as discussed in earlier sections, it is also possible to use SPLEND1D to explore the dynamics of the system. In this section, we briefly highlight such a possible study enabled by SPLEND1D. We perform a time-dependent simulation, based on a converged, detached simulation of the base-case, with $n_{u} = 3.0 \times 10^{19} m^{- 3}$ ⁠. We then introduce a sequence of short (⁠ $500 μ s$ ⁠) pulses during which the heat-sources, H_e and H_i, are amplified by a factor 10, before being relaxed to their initial values for $10 ms$ ⁠, Fig. 21(a). During each pulse, we observe a strong increase in the target parallel particle flux, $Γ_{| |, t}$ ⁠, Fig. 21(b), together with a strong increase in the target electron temperature $T_{e, t}$ Fig. 21(c), a sign that the plasma is reattaching during these events.

FIG. 21.

(a) Evolution of the target particle flux Γ||,t (left axis, blue) and imposed heat-source integral ∫Hi,e ds as a function of time, for a sequence of 10 heat pulses, with (b) a zoomed version of panel (a). (c) Time-evolution of the target electron temperature ( Te,t) and upstream electron temperature ( Te,u). (d) Time-evolution of the target density (nt), upstream density (nu) and target neutral density ( nn,t), and (e) time-evolution of the integrated plasma ( ∫n ds) and neutral ( ∫nn ds) densities along the flux tube.

View large Download slide

(a) Evolution of the target particle flux $Γ_{| |, t}$ (left axis, blue) and imposed heat-source integral $\int H_{i, e} d s$ as a function of time, for a sequence of 10 heat pulses, with (b) a zoomed version of panel (a). (c) Time-evolution of the target electron temperature (⁠ $T_{e, t}$ ⁠) and upstream electron temperature (⁠ $T_{e, u}$ ⁠). (d) Time-evolution of the target density (n_t), upstream density (n_u) and target neutral density (⁠ $n_{n, t}$ ⁠), and (e) time-evolution of the integrated plasma (⁠ $\int n d s$ ⁠) and neutral (⁠ $\int n_{n} d s$ ⁠) densities along the flux tube.

The target plasma and neutral densities are also strongly affected by the heat pulses. In particular, after a short increase in the target neutral density $n_{n, t}$ (likely due to an increased recycling caused by the increased $Γ_{| |, t}$ ⁠), $n_{n, t}$ drops below its steady-state value, due to the ionization of most of the neutrals present in the system, as evidenced in Fig. 21(e) by the strong decrease in the integrated neutral density in the flux tube, associated with an increase in the integrated plasma density. This leads to an increase in the plasma upstream density n_u and a complex dynamics of the plasma target density n_t, which, after an initial rise and drop, peaks again before relaxing to its steady-state value. These results highlight how SPLEND1D can be used for time-dependent simulations.

VII. CONCLUSION

This paper presented the SPLEND1D code: a 1D plasma fluid model that solves the Braginskii equations projected along a flux tube, together with a simple neutral fluid model. We have highlighted the key features of SPLEND1D, including its flexibility in terms of included SOL physics (magnetic field norm and pitch angle variation along the flux tube, dependent or independent electron and ion temperatures, possibility to enable, disable, or rescale any physical term contributing to the equations, etc.), numerical accuracy, and high computational speed. We believe this makes SPLEND1D a valuable tool for studying the complex dynamics of a divertor plasma interacting with a neutral gas and in contact with a wall. The SPLEND1D code was then used to investigate the physics of detachment in a reference simulation scenario. The code outputs the individual contributions of each term constituting its particle, momentum, and energy equations, allowing for the results to be easily interpreted. The roll-over of the target ion flux, and the onset of detachment, was found to be mainly due to momentum losses owing to charge-exchange reactions, and power losses owing to impurity radiation. We note, however, that the importance of atomic charge-exchange losses is due, in part, to the absence of the molecular charge-exchange loss channel, which is known to be an important player in the simulation of detachment dynamics. To address this, we therefore plan to upgrade the SPLEND1D model to account for molecular species in addition to atomic ones. As with all reduced models, SPLEND1D comes with significant assumptions and free parameters. In this paper, we have shown how they can influence simulation results. However, some important effects are not captured by the current formulation of SPLEND1D. In particular, momentum and energy losses due to radial transport are absent, as well as drifts and SOL currents. These effects, intrinsically relying on perpendicular dynamics, are difficult to capture in a 1D code, likely relying on some ansatz associated with free parameters. Such possibility is left for future work.

SPLEND1D is currently being used to interpret TCV experiments, in particular those related to alternative divertor configurations, where it is being employed to elucidate the role of parallel connection length on the onset of detachment, as well as the impact of total flux expansion on detachment threshold and SOL parallel profiles.

In the future, we plan on continuing the development of the SPLEND1D model. From a numerical perspective, it would be interesting to further improve the capabilities of SPLEND1D, for instance by adaptively refining the mesh depending on local gradients of the solution, which may prove useful in situations where the detachment “front” is moving away from the target. From a physics perspective, we plan to extend SPLEND1D to add the capability to handle multiple plasma and neutral species, in particular molecules, which are known to play a very significant role in the divertor dynamics.^31,32 The impurity model implemented in SPLEND1D, currently limited to constant impurity fraction assuming coronal equilibrium, should also be further enhanced to properly account for the impact of impurity transport and the subsequent occurrence of non-coronal effects. This could be achieved by adding a collisional radiative model for impurities, evolving in time the concentration of impurities' charge states, either with a trace assumption, similarly to Ref. 41, or using a more advanced fluid closure such as the Zdhanov closure.

ACKNOWLEDGMENTS

This work was supported in part by the Swiss National Science Foundation. This work has been carried out within the framework of the EUROfusion Consortium, via the Euratom Research and Training Programme (Grant Agreement No. 101052200 – EUROfusion) and funded by the Swiss State Secretariat for Education, Research and Innovation (SERI). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union, the European Commission, or SERI. Neither the European Union nor the European Commission nor SERI can be held responsible for them.

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

O. Février: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (lead); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). S. Gorno: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal). C. Theiler: Conceptualization (equal); Methodology (equal); Project administration (equal); Supervision (lead); Writing – original draft (equal); Writing – review & editing (equal). M. Carpita: Validation (equal); Writing – review & editing (equal). G. Durr-Legoupil-Nicoud: Validation (equal); Writing – review & editing (equal). M. von Allmen: Software (equal); Validation (equal); Writing – review & editing (equal).

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

APPENDIX A: COULOMB LOGARITHM

The Coulomb logarithm,

ln Λ

⁠, is a slow varying function of density and temperature, and is therefore typically set constant. It is, however, also possible to enforce a local computation of

ln Λ

⁠, based on the local n_e and T_e. The usual definition of

ln Λ

yields,¹⁹ for n_e expressed in

m^{- 3}

and T_e in

eV

⁠,

ln Λ = ln Λ_{low} = 23 - ln ({(n_{e} \times 10^{- 6})}^{0.5} T_{e}^{- 1.5}),

(A1)

for

T_{e} < 10 eV

and

ln Λ = ln Λ_{high} = 24 - ln ({(n_{e} \times 10^{- 6})}^{0.5} T_{e}^{- 1}),

(A2)

for

T_{e} > 10 eV

⁠. Using such a definition of

ln Λ

leads to a discontinuity at

T_{e} = 10 eV

⁠, which causes numerical difficulties to the non-linear solver by introducing a singularity in the Jacobian. To avoid this issue, we introduce a transition parameter, Δ, defined as

Δ = \frac{1}{2} [tanh (\frac{T_{e} [eV] - 10 [eV]}{0.1 [eV]}) + 1],

(A3)

such that

ln Λ = (1 - Δ) ln Λ_{low} + Δ ln Λ_{high} .

(A4)

Equation is a smooth function of n_e and T_e, thus avoiding the numerical difficulties associated with a discontinuous

ln Λ

⁠.

APPENDIX B: CONTRIBUTIONS OF MOMENTUM AND POWER SOURCE TERMS TO LOSS FACTORS

The momentum and power source/sink terms can be calculated from the SPLEND1D model to evaluate the significance of each term, as shown in Fig. 12. Integrating the steady state form of the momentum equation along the flux tube length from the target to any upstream location, we can write the SPLEND1D equations in the form

p_{tot, u} - p_{tot, t} = f_{mom} p_{tot, u}

⁠. We can then separate each of the contributing terms to find

f_{mom} = f_{mom}^{atomic} + f_{mom}^{visc} + f_{mom}^{\nabla_{| |}}

⁠, where

f_{mom}^{atomic} = \frac{1}{p_{tot, u}} \int_{t}^{u} S_{p, | |}^{u} d s,

(B1)

f_{mom}^{visc} = \frac{1}{p_{tot, u}} \int_{t}^{u} \frac{4}{3} [B^{3 / 2} \frac{\partial}{\partial s} (η_{| |} B^{- 2} \frac{\partial B^{1 / 2} u_{| |}}{\partial s})] d s,

(B2)

f_{mom}^{\nabla_{| |}} = \frac{1}{p_{tot, u}} \int_{t}^{u} \frac{\partial B}{\partial s} \frac{m_{i} n u_{| |}^{2}}{B} d s .

(B3)

f_{mom}^{atomic}

can be further decomposed as the sum of each atomic process included (ionization, recombination and charge-exchange, Eq. . Similarly to the momentum equation, each contributing terms in f_pow can be separated, leading to

f_{pow} = f_{power}^{atomic} + f_{power}^{imp}

⁠, where

f_{pow}^{atomic} \frac{q_{i n}}{B_{u}} = \int_{t}^{u} \frac{(S_{e}^{E} + S_{i}^{E})}{B} d s,

(B4)

f_{pow}^{imp} \frac{q_{i n}}{B_{u}} = \int_{t}^{u} \frac{S_{imp}^{E}}{B} d s .

(B5)

f_{pow}^{atomic}

can also be expanded as a sum of each atomic process included [Eqs. and ]: ionization, recombination, charge-exchange, and excitation.

APPENDIX C: ADVECTIVE NEUTRAL MODEL

In this appendix, we present the derivation of an advective model for the neutral dynamics. In this model, the

n_{n} V_{n}

term is described by the addition of a momentum equation. We also add ad hoc loss terms involving τ_n in the momentum and energy equations, that will be developed hereafter. We start from the general form of the momentum equation,

\frac{\partial (m_{n} n_{n} V_{n})}{\partial t} + \nabla \cdot (m_{n} n_{n} V_{n} \otimes V_{n}) = - \nabla p_{n} + S_{n}^{u},

(C1)

- \frac{m_{n} n_{n} V_{n}}{τ_{n}} - \nabla \cdot Π_{n},

(C2)

where

p_{n} = n_{n} T_{n}

is the static neutral pressure,

S_{n}^{u}

is the neutral momentum source term, and

Π_{n}

is the viscous stress tensor. In principle, the viscous stress tensor could be self-consistently evaluated. However, due to the complexity of the resulting expression (see Ref. 30 for instance), we use a simpler form. Projected along the θ and

φ

directions, and using

\frac{\partial}{\partial φ} = 0

⁠, we get

\frac{\partial (m_{n} n_{n} V_{n}^{θ})}{\partial t} + B \frac{\partial}{\partial s} (\frac{n_{n} m_{n} V_{n}^{θ}^{2}}{B sin α}) + \frac{n_{n} m_{n} V_{n}^{φ}^{2}}{B_{φ} sin α} \frac{\partial B_{φ}}{\partial s} = - \frac{1}{sin α} \frac{\partial p_{n}}{\partial s} + S_{n, θ}^{u} - \frac{m_{n} n_{n} V_{n}^{θ}}{τ_{n}} + \frac{\partial}{\partial s} (\frac{η_{n}}{{sin}^{2} α} \frac{\partial V_{n}^{θ}}{\partial s}),

(C3)

\frac{\partial (m_{n} n_{n} V_{n}^{φ})}{\partial t} + B \frac{\partial}{\partial s} (\frac{n_{n} m_{n} V_{n}^{θ} V_{n}^{φ}}{B sin α}) - \frac{n_{n} m_{n} V_{n}^{θ} V_{n}^{φ}}{B_{φ} sin α} \frac{\partial B_{φ}}{\partial s} = S_{n, φ}^{u} - \frac{m_{n} n_{n} V_{n}^{φ}}{τ_{n}} + \frac{\partial}{\partial s} (\frac{η_{n}}{{sin}^{2} α} \frac{\partial V_{n}^{φ}}{\partial s}),

(C4)

where the neutral viscosity η_n is defined as

η_{n} = \frac{n_{n} T_{n}}{(n {⟨ σ v ⟩}_{c x} + n {⟨ σ v ⟩}_{ion})} .

(C5)

This model allows for neutral trajectories that are not aligned to the magnetic field, but constrained within a flux surface. For the neutral energy equation, the procedure is similar. Starting from the general form of the equation, we have

\frac{\partial E_{n}}{\partial t} + \nabla \cdot ([E_{n} + p_{n}] V_{n} + q_{n} + V_{n} \cdot (\nabla \cdot Π_{n})) = S_{n}^{E} - \frac{E_{n}}{τ_{N}},

(C6)

where

S_{n}^{E}

represents the source and sink terms, discussed in Sec. II D. E_n is the total energy, defined as

E_{n} = \frac{3}{2} n_{n} T_{n} + \frac{1}{2} n_{n} m_{n} V_{n}^{2},

(C7)

where

V_{n}^{2} = V_{n}^{θ}^{2} + V_{n}^{φ}^{2}

⁠.

q_{n}

is the neutral conductive heat flux. Denoting κ_n the neutral heat conduction, we write²⁹

κ_{n} = \frac{5}{2} \frac{η_{n}}{m_{n}} = \frac{5}{2} \frac{n_{n} T_{n}}{m_{n} (n {⟨ σ v ⟩}_{c x} + n {⟨ σ v ⟩}_{ion})},

(C8)

following similar notations as for Eq. . The neutral heat flux

q_{n}

is then defined as

q_{n} = - κ_{n} \nabla T_{n} .

(C9)

Developing the divergence, this leads to

\frac{\partial E_{n}}{\partial t} + B \frac{\partial}{\partial s} (\frac{[E + p_{n}] V_{n}^{θ} + q_{n}^{θ} + {[V_{n} \cdot (\nabla \cdot Π_{n})]}_{θ}}{B sin α}) = S_{n}^{E} - \frac{E_{n}}{τ_{n}} .

(C10)

Equations , , , and constitute a set of four equations that are solved to describe the behavior of the neutrals. In Appendix D, we briefly show how Eqs. , , and could be further developed to form a so-called “pressure-diffusion” model, in the spirit of models that have been developed in Ref. 27. The implementation and test of such formulation in SPLEND1D is, however, left for future work.

APPENDIX D: PRESSURE-DIFFUSION EQUATION

In this appendix, we show how Eqs. and can be used to establish a pressure-diffusion equation for

V_{n}

⁠, which can then be incorporated in Eq. . This derivation is similar to the one presented in Ref. 27, and differs only by the presence of

\frac{1}{τ_{N}}

terms. Starting from Eq. , we assume the neutral population to be at steady-state, such that the time-derivative can be removed. Furthermore, we assume the neutral flow to be strongly subsonic, such that the inertia term

n_{n} m_{n} V_{n}^{θ}^{2}

can be neglected when compared to the neutral pressure p_n. We simplify the expression by assuming that the variation of the magnetic field norm along the flux tube can be neglected. Under these assumptions, Eq. can then be simplified as

0 = - \frac{1}{sin α} \frac{\partial p_{n}}{\partial s} - m_{i} n n_{n} {⟨ σ v ⟩}_{ion} V_{n}^{θ} + m_{i} n^{2} {⟨ σ v ⟩}_{rec} u_{| |} sin α - m_{i} n_{n} n {⟨ σ v ⟩}_{C X} (V_{n}^{θ} - u_{| |} sin α) - \frac{m_{n} n_{n} V_{n}^{θ}}{τ_{n}},

(D1)

leading to

n_{n} V_{n}^{θ} = \frac{(n_{n} {⟨ σ v ⟩}_{C X} + n {⟨ σ v ⟩}_{rec}) n u_{| |} sin α - \frac{1}{m_{n} sin α} \frac{\partial p_{n}}{\partial s}}{n {⟨ σ v ⟩}_{ion} + n {⟨ σ v ⟩}_{C X} + \frac{1}{τ_{n}}},

(D2)

which can readily be incorporated in Eq. . Similarly, and under the assumption that

n_{n} m_{n} V_{n}^{θ} V_{n}^{φ}

is small, Eq. can be rewritten as

n_{n} V_{n}^{φ} = \frac{(n {⟨ σ v ⟩}_{rec} + n_{n} {⟨ σ v ⟩}_{C X}) n u_{| |} cos α}{n {⟨ σ v ⟩}_{ion} + n {⟨ σ v ⟩}_{C X} + \frac{1}{τ_{n}}} .

(D3)

The implementation and test of this formulation in SPLEND1D is, however, left for future work.

APPENDIX E: ON THE ROLE OF THE NEUTRAL MODEL

The choice of the neutral model, and the neutral temperature implementation, can affect the plasma dynamics, primarily through the momentum and energy volumetric source terms. The diffusive model [Eqs. (32) and (33)] can be implemented with either T_n = T_i or with a constant, imposed value of T_n (in this section we set $T_{n} = 1$ eV to model cold neutrals), while the advective model [Eqs. (28), (C3), and (C4)] can be implemented with either T_n = T_i or with the evolution of T_n following Eq. (C10).

Figure 22 shows the main differences in the results of the base case simulations for each neutral model outlined above. In the advective model with a self-consistent T_n, and for the diffusive model with cold neutrals (⁠ $T_{n} = 1$ eV), we see an earlier target ion current rollover, along with a lower target plasma pressure, compared to implementing T_n = T_i. The cold neutrals facilitate detachment by increasing both the power loss factor with respect to the cases with T_n = T_i, as a result of increasing the power transferred through charge exchange reactions, by increasing the energy transferred from the ions to neutrals. This, in turn, leads to an increase in the momentum losses by promoting enhanced charge-exchange reaction and reduced ionization. We remark here that the density ramps presented in this paper are performed by ramping up the value of H_n. For different values of H_n, steady-state solutions can have the same upstream density n_u.

FIG. 22.

View large Download slide

(a) Target ion current, (b) momentum loss factor, (c) total plasma pressure, and (d) power loss factor, as a function of n_u for a range of different neutral models/options implemented in SPLEND1D. The base case density ramp is performed using the following neutral models: diffusive model, with T_n = T_i; diffusive model, with $T_{n} = 1$ eV; advective model, with T_n = T_i; advective model, evolving T_n according to Eq. (C10).

The advective model, in comparison with the diffusive model, features an additional neutral transport mechanism, the advective cross-“flux-tube” transport. For a fixed τ_N, as shown in Fig. 22, the advective model displays a much weaker rollover. The total momentum and power loss factors shown in Fig. 22 also differ between the models.

REFERENCES

1.

P.

Stangeby

,

The Plasma Boundary of Magnetic Fusion Devices

, Series in Plasma Physics and Fluid Dynamics (

Taylor & Francis

,

2000

).

Google Scholar

Crossref

2.

S. I.

Krasheninnikov

,

A. S.

Kukushkin

, and

A. A.

Pshenov

,

Phys. Plasmas

23

,

055602

(

2016

).

https://doi.org/10.1063/1.4948273

Google Scholar

Crossref

3.

A. W.

Leonard

,

Plasma Phys. Controlled Fusion

60

,

044001

(

2018

).

https://doi.org/10.1088/1361-6587/aaa7a9

Google Scholar

Crossref

4.

A.

Loarte

,

R.

Monk

,

J.

Martín-Solís

,

D.

Campbell

,

A.

Chankin

,

S.

Clement

,

S.

Davies

,

J.

Ehrenberg

,

S.

Erents

,

H.

Guo

,

P.

Harbour

,

L.

Horton

,

L.

Ingesson

,

H.

Jäckel

,

J.

Lingertat

,

C.

Lowry

,

C.

Maggi

,

G.

Matthews

,

K.

McCormick

,

D.

O'Brien

,

R.

Reichle

,

G.

Saibene

,

R.

Smith

,

M.

Stamp

,

D.

Stork

, and

G.

Vlases

,

Nucl. Fusion

38

,

331

(

1998

).

https://doi.org/10.1088/0029-5515/38/3/303

Google Scholar

Crossref

5.

K.

Verhaegh

,

B.

Lipschultz

,

B.

Duval

,

O.

Février

,

A.

Fil

,

C.

Theiler

,

M.

Wensing

,

C.

Bowman

,

D.

Gahle

,

J.

Harrison

,

B.

Labit

,

C.

Marini

,

R.

Maurizio

,

H.

de Oliveira

,

H.

Reimerdes

,

U.

Sheikh

,

C.

Tsui

,

N.

Vianello

,

W.

Vijvers

,

the TCV Team

, and

the EUROfusion MST1 Team

,

Nucl. Fusion

59

,

100922

(

2019

).

https://doi.org/10.1088/1741-4326/ab4251

Google Scholar

Crossref

6.

P.

Stangeby

and

A.

Leonard

,

Nucl. Fusion

51

,

063001

(

2011

).

https://doi.org/10.1088/0029-5515/51/6/063001

Google Scholar

Crossref

7.

A.

Kallenbach

,

M.

Bernert

,

R.

Dux

,

L.

Casali

,

T.

Eich

,

L.

Giannone

,

A.

Herrmann

,

R.

McDermott

,

A.

Mlynek

,

H. W.

Müller

,

F.

Reimold

,

J.

Schweinzer

,

M.

Sertoli

,

G.

Tardini

,

W.

Treutterer

,

E.

Viezzer

,

R.

Wenninger

,

M.

Wischmeier

, and

the ASDEX Upgrade Team

,

Plasma Phys. Controlled Fusion

55

,

124041

(

2013

).

https://doi.org/10.1088/0741-3335/55/12/124041

Google Scholar

Crossref

8.

B.

Lipschultz

,

B.

LaBombard

,

J. L.

Terry

,

C.

Boswell

, and

I. H.

Hutchinson

,

Fusion Sci. Technol.

51

,

369

(

2007

).

https://doi.org/10.13182/FST07-A1428

Google Scholar

Crossref

9.

S.

Potzel

,

M.

Wischmeier

,

M.

Bernert

,

R.

Dux

,

H.

Müller

, and

A.

Scarabosio

, and

the ASDEX Upgrade Team

,

Nucl. Fusion

54

,

013001

(

2014

).

https://doi.org/10.1088/0029-5515/54/1/013001

Google Scholar

Crossref

10.

S. I.

Braginskii

,

Rev. Plasma Phys.

1

,

205

(

1965

), https://ui.adsabs.harvard.edu/abs/1965RvPP....1..205B/abstract.

Google Scholar

11.

S.

Nakazawa

,

N.

Nakajima

,

M.

Okamoto

, and

N.

Ohyabu

,

Plasma Phys. Controlled Fusion

42

,

401

(

2000

).

https://doi.org/10.1088/0741-3335/42/4/303

Google Scholar

Crossref

12.

R.

Goswami

,

P.

Kaw

,

M.

Warrier

,

R.

Singh

, and

S.

Deshpande

,

Phys. Plasmas

8

,

857

(

2001

).

https://doi.org/10.1063/1.1342028

Google Scholar

Crossref

13.

S.

Togo

,

M.

Nakamura

,

Y.

Ogawa

,

K.

Shimizu

,

T.

Takizuka

, and

K.

Hoshino

,

Plasma Fusion Res.

8

,

2403096

(

2013

).

https://doi.org/10.1585/pfr.8.2403096

Google Scholar

Crossref

14.

E.

Havlíčková

,

W.

Fundamenski

,

F.

Subba

,

D.

Coster

,

M.

Wischmeier

, and

G.

Fishpool

,

Plasma Phys. Controlled Fusion

55

,

065004

(

2013

).

https://doi.org/10.1088/0741-3335/55/6/065004

Google Scholar

Crossref

15.

B. D.

Dudson

,

J.

Allen

,

T.

Body

,

B.

Chapman

,

C.

Lau

,

L.

Townley

,

D.

Moulton

,

J.

Harrison

, and

B.

Lipschultz

,

Plasma Phys. Controlled Fusion

61

,

065008

(

2019

).

https://doi.org/10.1088/1361-6587/ab1321

Google Scholar

Crossref

16.

G. L.

Derks

,

J. P. K. W.

Frankemölle

,

J. T. W.

Koenders

,

M.

van Berkel

,

H.

Reimerdes

,

M.

Wensing

, and

E.

Westerhof

,

Plasma Phys. Controlled Fusion

64

,

125013

(

2022

).

https://doi.org/10.1088/1361-6587/ac9dbd

Google Scholar

Crossref

17.

H.

Reimerdes

,

M.

Agostini

,

E.

Alessi

,

S.

Alberti

,

Y.

Andrebe

,

H.

Arnichand

,

J.

Balbin

,

F.

Bagnato

,

M.

Baquero-Ruiz

,

M.

Bernert

,

W.

Bin

,

P.

Blanchard

,

T.

Blanken

,

J.

Boedo

,

D.

Brida

,

S.

Brunner

,

C.

Bogar

,

O.

Bogar

,

T.

Bolzonella

,

F.

Bombarda

,

F.

Bouquey

,

C.

Bowman

,

D.

Brunetti

,

J.

Buermans

,

H.

Bufferand

,

L.

Calacci

,

Y.

Camenen

,

S.

Carli

,

D.

Carnevale

,

F.

Carpanese

,

F.

Causa

,

J.

Cavalier

,

M.

Cavedon

,

J.

Cazabonne

,

J.

Cerovsky

,

R.

Chandra

,

A. C.

Jayalekshmi

,

O.

Chellaï

,

P.

Chmielewski

,

D.

Choi

,

G.

Ciraolo

,

I.

Classen

,

S.

Coda

,

C.

Colandrea

,

A. D.

Molin

,

P.

David

,

M.

de Baar

,

J.

Decker

,

W.

Dekeyser

,

H.

de Oliveira

,

D.

Douai

,

M.

Dreval

,

M.

Dunne

,

B.

Duval

,

S.

Elmore

,

O.

Embreus

,

F.

Eriksson

,

M.

Faitsch

,

G.

Falchetto

,

M.

Farnik

,

A.

Fasoli

,

N.

Fedorczak

,

F.

Felici

,

O.

Février

,

O.

Ficker

,

A.

Fil

,

M.

Fontana

,

E.

Fransson

,

L.

Frassinetti

,

I.

Furno

,

D.

Gahle

,

D.

Galassi

,

K.

Galazka

,

C.

Galperti

,

S.

Garavaglia

,

M.

Garcia-Munoz

,

B.

Geiger

,

M.

Giacomin

,

G.

Giruzzi

,

M.

Gobbin

,

T.

Golfinopoulos

,

T.

Goodman

,

S.

Gorno

,

G.

Granucci

,

J.

Graves

,

M.

Griener

,

M.

Gruca

,

T.

Gyergyek

,

R.

Haelterman

,

A.

Hakola

,

W.

Han

,

T.

Happel

,

G.

Harrer

,

J.

Harrison

,

S.

Henderson

,

G.

Hogeweij

,

J.-P.

Hogge

,

M.

Hoppe

,

J.

Horacek

,

Z.

Huang

,

A.

Iantchenko

,

P.

Innocente

,

K. I.

Björk

,

C.

Ionita-Schrittweiser

,

H.

Isliker

,

A.

Jardin

,

R.

Jaspers

,

R.

Karimov

,

A.

Karpushov

,

Y.

Kazakov

,

M.

Komm

,

M.

Kong

,

J.

Kovacic

,

O.

Krutkin

,

O.

Kudlacek

,

U.

Kumar

,

R.

Kwiatkowski

,

B.

Labit

,

L.

Laguardia

,

J.

Lammers

,

E.

Laribi

,

E.

Laszynska

,

A.

Lazaros

,

O.

Linder

,

B.

Linehan

,

B.

Lipschultz

,

X.

Llobet

,

J.

Loizu

,

T.

Lunt

,

E.

Macusova

,

Y.

Marandet

,

M.

Maraschek

,

G.

Marceca

,

C.

Marchetto

,

S.

Marchioni

,

E.

Marmar

,

Y.

Martin

,

L.

Martinelli

,

F.

Matos

,

R.

Maurizio

,

M.-L.

Mayoral

,

D.

Mazon

,

V.

Menkovski

,

A.

Merle

,

G.

Merlo

,

H.

Meyer

,

K.

Mikszuta-Michalik

,

P. M.

Cabrera

,

J.

Morales

,

J.-M.

Moret

,

A.

Moro

,

D.

Moulton

,

H.

Muhammed

,

O.

Myatra

,

D.

Mykytchuk

,

F.

Napoli

,

R.

Nem

,

A.

Nielsen

,

M.

Nocente

,

S.

Nowak

,

N.

Offeddu

,

J.

Olsen

,

F.

Orsitto

,

O.

Pan

,

G.

Papp

,

A.

Pau

,

A.

Perek

,

F.

Pesamosca

,

Y.

Peysson

,

L.

Pigatto

,

C.

Piron

,

M.

Poradzinski

,

L.

Porte

,

T.

Pütterich

,

M.

Rabinski

,

H.

Raj

,

J.

Rasmussen

,

G.

Rattá

,

T.

Ravensbergen

,

D.

Ricci

,

P.

Ricci

,

N.

Rispoli

,

F.

Riva

,

J.

Rivero-Rodriguez

,

M.

Salewski

,

O.

Sauter

,

B.

Schmidt

,

R.

Schrittweiser

,

S.

Sharapov

,

U.

Sheikh

,

B.

Sieglin

,

M.

Silva

,

A.

Smolders

,

A.

Snicker

,

C.

Sozzi

,

M.

Spolaore

,

A.

Stagni

,

L.

Stipani

,

G.

Sun

,

T.

Tala

,

P.

Tamain

,

K.

Tanaka

,

A. T.

Biwole

,

D.

Terranova

,

J.

Terry

,

D.

Testa

,

C.

Theiler

,

A.

Thornton

,

A.

Thrysøe

,

H.

Torreblanca

,

C.

Tsui

,

D.

Vaccaro

,

M.

Vallar

,

M.

van Berkel

,

D. V.

Eester

,

R.

van Kampen

,

S. V.

Mulders

,

K.

Verhaegh

,

T.

Verhaeghe

,

N.

Vianello

,

F.

Villone

,

E.

Viezzer

,

B.

Vincent

,

I.

Voitsekhovitch

,

N.

Vu

,

N.

Walkden

,

T.

Wauters

,

H.

Weisen

,

N.

Wendler

,

M.

Wensing

,

F.

Widmer

,

S.

Wiesen

,

M.

Wischmeier

,

T.

Wijkamp

,

D.

Wünderlich

,

C.

Wüthrich

,

V.

Yanovskiy

,

J.

Zebrowski

, and

the EUROfusion MST1 Team

,

Nucl. Fusion

62

,

042018

(

2022

).

https://doi.org/10.1088/1741-4326/ac369b

Google Scholar

Crossref

18.

S.

Gorno

,

O.

Février

,

C.

Theiler

,

T.

Ewalds

,

F.

Felici

,

T.

Lunt

,

A.

Merle

,

F.

Bagnato

,

C.

Colandrea

,

J.

Degrave

,

R.

Ducker

,

G.

Durr-Legoupil-Nicoud

,

B. P.

Duval

,

K.

Lee

,

L.

Martinelli

,

D.

Oliveira

,

A.

Perek

,

H.

Reimerdes

,

L.

Simons

,

G.

Sun

,

B.

Tracey

,

M.

Wischmeier

, and

C.

Wuethrich

, “

X-point radiator and power exhaust control in configurations with multiple X-points in TCV

,”

Phys. Plasmas

31

,

072504

(

2024

).

https://doi.org/10.1063/5.020140

Google Scholar

19.

U.S. Naval Research Laboratory

, see https://www.nrl.navy.mil/News-Media/Publications/NRL-Plasma-Formulary/ for “

Nrl plasma formulary

” (last accessed December 01, 2024).

20.

M.

Day

,

B.

Merriman

,

F.

Najmabadi

, and

R. W.

Conn

,

Contrib. Plasma Phys.

36

,

419

(

1996

).

https://doi.org/10.1002/ctpp.2150360259

Google Scholar

Crossref

21.

W.

Fundamenski

,

Plasma Phys. Controlled Fusion

47

,

R163

(

2005

).

https://doi.org/10.1088/0741-3335/47/11/R01

Google Scholar

Crossref

22.

G.

Ciraolo

,

H.

Bufferand

,

P.

Di Cintio

,

P.

Ghendrih

,

S.

Lepri

,

R.

Livi

,

Y.

Marandet

,

E.

Serre

,

P.

Tamain

, and

M.

Valentinuzzi

,

Contrib. Plasma Phys.

58

,

457

(

2018

).

https://doi.org/10.1002/ctpp.201700222

Google Scholar

Crossref

23.

R.

Schneider

,

X.

Bonnin

,

K.

Borrass

,

D. P.

Coster

,

H.

Kastelewicz

,

D.

Reiter

,

V. A.

Rozhansky

, and

B. J.

Braams

,

Contrib. Plasma Phys.

46

,

3

(

2006

).

https://doi.org/10.1002/ctpp.200610001

Google Scholar

Crossref

24.

E.

Zawaideh

,

F.

Najmabadi

, and

R. W.

Conn

,

Phys. Fluids

29

,

463

(

1986

).

https://doi.org/10.1063/1.865731

Google Scholar

Crossref

25.

E.

Zawaideh

,

N. S.

Kim

, and

F.

Najmabadi

,

Phys. Fluids

31

,

3280

(

1988

).

https://doi.org/10.1063/1.866940

Google Scholar

Crossref

26.

D.

Reiter

,

M.

Baelmans

, and

P.

Börner

,

Fusion Sci. Technol.

47

,

172

(

2005

).

https://doi.org/10.13182/FST47-172

Google Scholar

Crossref

27.

N.

Horsten

,

W.

Dekeyser

,

G.

Samaey

, and

M.

Baelmans

,

Nucl. Mater. Energy

12

,

869

(

2017

).

https://doi.org/10.1016/j.nme.2016.11.007

Google Scholar

Crossref

28.

N.

Horsten

,

M.

Groth

,

W.

Dekeyser

,

W.

Van Uytven

,

S.

Aleiferis

,

S.

Carli

,

J.

Karhunen

,

K.

Lawson

,

B.

Lomanowski

,

A.

Meigs

,

S.

Menmuir

,

A.

Shaw

,

V.

Solokha

, and

B.

Thomas

,

Nucl. Mater. Energy

33

,

101247

(

2022

).

https://doi.org/10.1016/j.nme.2022.101247

Google Scholar

Crossref

29.

W. V.

Uytven

,

W.

Dekeyser

,

M.

Blommaert

,

S.

Carli

, and

M.

Baelmans

,

Nucl. Fusion

62

,

086023

(

2022

).

https://doi.org/10.1088/1741-4326/ac72b4

Google Scholar

Crossref

30.

N.

Horsten

,

W.

Dekeyser

,

G.

Samaey

,

P.

Börner

, and

M.

Baelmans

,

Contrib. Plasma Phys.

56

,

610

(

2016

).

https://doi.org/10.1002/ctpp.201610005

Google Scholar

Crossref

31.

K.

Verhaegh

,

B.

Lipschultz

,

J.

Harrison

,

B.

Duval

,

A.

Fil

,

M.

Wensing

,

C.

Bowman

,

D.

Gahle

,

A.

Kukushkin

,

D.

Moulton

,

A.

Perek

,

A.

Pshenov

,

F.

Federici

,

O.

Février

,

O.

Myatra

,

A.

Smolders

,

C.

Theiler

,

the TCV Team

, and

the EUROfusion MST1 Team

,

Nucl. Fusion

61

,

106014

(

2021

).

https://doi.org/10.1088/1741-4326/ac1dc5

Google Scholar

Crossref

32.

Y.

Zhou

,

B.

Dudson

,

F.

Militello

,

K.

Verhaegh

, and

O.

Myatra

,

Plasma Phys. Controlled Fusion

64

,

065006

(

2022

).

https://doi.org/10.1088/1361-6587/ac6827

Google Scholar

Crossref

33.

O.

Myatra

,

D.

Moulton

,

B.

Dudson

,

B.

Lipschultz

,

S.

Newton

,

K.

Verhaegh

, and

A.

Fil

,

Nucl. Fusion

63

,

076030

(

2023

).

https://doi.org/10.1088/1741-4326/acd9da

Google Scholar

Crossref

34.

H. P.

Summers

,

W. J.

Dickson

,

M. G.

O'Mullane

,

N. R.

Badnell

,

A. D.

Whiteford

,

D. H.

Brooks

,

J.

Lang

,

S. D.

Loch

, and

D. C.

Griffin

,

Plasma Phys. Controlled Fusion

48

,

263

(

2006

).

https://doi.org/10.1088/0741-3335/48/2/007

Google Scholar

Crossref

35.

D.

Reiter

, “

The data file AMJUEL: Additional atomic and molecular data for EIRENE

” (2020), https://eirene.de/Documentation/amjuel.pdf.

36.

H.

Frerichs

,

Y.

Feng

,

X.

Bonnin

,

R. A.

Pitts

,

D.

Reiter

, and

O.

Schmitz

,

Phys. Plasmas

28

,

102503

(

2021

).

https://doi.org/10.1063/5.0062248

Google Scholar

Crossref

37.

M.

Valentinuzzi

, “

Modélisation numérique des flux de puissances sur les composants face au plasma de Tokamak à l'aide de techniques de couplage avancées entre codes fluides et cinétiques

,” Ph.D. thesis (Aix-Marseille Université,

2018

).

Google Scholar

38.

D.

Post

,

R.

Jensen

,

C.

Tarter

,

W.

Grasberger

, and

W.

Lokke

,

At. Data Nucl. Data Tables

20

,

397

(

1977

).

https://doi.org/10.1016/0092-640X(77)90026-2

Google Scholar

Crossref

39.

D.

Wagner

and

J.

Schwartz

, see https://github.com/cfe316/atomic “

Atomic toolbox

” (

2016

).

Google Scholar

40.

A.

Kallenbach

,

M.

Bernert

,

R.

Dux

,

F.

Reimold

, and

M.

Wischmeier

, and

A. U. Team

,

Plasma Phys. Controlled Fusion

58

,

045013

(

2016

).

https://doi.org/10.1088/0741-3335/58/4/045013

Google Scholar

Crossref

41.

Y.

Zhou

,

B.

Dudson

,

T.

Wu

,

Z.

Wang

,

T.

Xia

,

C.

Zhong

,

J.

Gao

,

H.

Du

, and

D.

Fan

,

Plasma Phys. Controlled Fusion

66

,

055005

(

2024

).

https://doi.org/10.1088/1361-6587/ad30fc

Google Scholar

Crossref

42.

R.

Maurizio

,

B.

Duval

,

B.

Labit

,

H.

Reimerdes

,

C.

Theiler

,

C.

Tsui

,

J.

Boedo

,

H. D.

Oliveira

,

O.

Février

,

U.

Sheikh

,

M.

Spolaore

,

K.

Verhaegh

,

N.

Vianello

, and

M.

Wensing

,

Nucl. Mater. Energy

19

,

372

(

2019

).

https://doi.org/10.1016/j.nme.2019.03.020

Google Scholar

Crossref

43.

K. U.

Riemann

,

J. Phys. D: Appl. Phys.

24

,

493

(

1991

).

https://doi.org/10.1088/0022-3727/24/4/001

Google Scholar

Crossref

44.

S.

Wiesen

,

D.

Reiter

,

V.

Kotov

,

M.

Baelmans

,

W.

Dekeyser

,

A.

Kukushkin

,

S.

Lisgo

,

R.

Pitts

,

V.

Rozhansky

,

G.

Saibene

,

I.

Veselova

, and

S.

Voskoboynikov

,

J. Nucl. Mater.

463

,

480

(

2015

).

https://doi.org/10.1016/j.jnucmat.2014.10.012

Google Scholar

Crossref

45.

T.

Rognlien

,

M.

Rensink

, and

D.

Stotler

,

Fusion Eng. Des.

135

,

380

(

2018

).

https://doi.org/10.1016/j.fusengdes.2017.07.024

Google Scholar

Crossref

46.

G.

Puppo

and

M.

Semplice

,

J. Sci. Comput.

66

,

1052

–

1076

(

2016

).

https://doi.org/10.1007/s10915-015-0056-x

Google Scholar

Crossref

47.

S.

Balay

,

S.

Abhyankar

,

M. F.

Adams

,

J.

Brown

,

P.

Brune

,

K.

Buschelman

,

L.

Dalcin

,

V.

Eijkhout

,

W. D.

Gropp

,

D.

Kaushik

,

M. G.

Knepley

,

L. C.

McInnes

,

K.

Rupp

,

B. F.

Smith

,

S.

Zampini

,

H.

Zhang

, and

H.

Zhang

, “

PETSc users manual

,”

Report No. ANL-95/11

(

Argonne National Laboratory

,

2016

).

Google Scholar

48.

S.

Balay

,

W. D.

Gropp

,

L. C.

McInnes

, and

B. F.

Smith

, in

Modern Software Tools in Scientific Computing

, edited by

E.

Arge

,

A. M.

Bruaset

, and

H. P.

Langtangen

(

Birkhäuser Press

,

1997

), pp.

163

–

202

.

Google Scholar

Crossref

49.

M.

Carpita

,

O.

Février

,

H.

Reimerdes

,

C.

Theiler

,

B.

Duval

,

C.

Colandrea

,

G.

Durr-Legoupil-Nicoud

,

D.

Galassi

,

S.

Gorno

,

E.

Huett

,

J.

Loizu

,

L.

Martinelli

,

A.

Perek

,

L.

Simons

,

G.

Sun

,

E.

Tonello

,

C.

Wüthrich

, and

the TCV Team

,

Nucl. Fusion

64

,

046019

(

2024

).

https://doi.org/10.1088/1741-4326/ad2a2a

Google Scholar

Crossref

50.

V.

Kotov

and

D.

Reiter

,

Plasma Phys. Controlled Fusion

51

,

115002

(

2009

).

https://doi.org/10.1088/0741-3335/51/11/115002

Google Scholar

Crossref

51.

D.

Moulton

,

J.

Harrison

,

B.

Lipschultz

, and

D.

Coster

,

Plasma Phys. Controlled Fusion

59

,

065011

(

2017

).

https://doi.org/10.1088/1361-6587/aa6b13

Google Scholar

Crossref

52.

P. C.

Stangeby

,

Plasma Phys. Controlled Fusion

60

,

044022

(

2018

).

https://doi.org/10.1088/1361-6587/aaacf6

Google Scholar

Crossref

SPLEND1D, a reduced one-dimensional model to investigate the physics of plasma detachment

I. INTRODUCTION

II. MODEL

A. Geometry

B. Plasma model

1. Equations

2. Heat fluxes

3. Viscosity tensor

C. Fluid model for the neutrals

1. Continuity equation for the neutrals

2. Diffusive neutral model

D. Atomic source terms

1. Rates

2. Sources and sinks

E. Impurity radiation

F. Boundary conditions

G. Numerical implementation

III. PRESENTATION OF BASE CASE

A. Parameters

B. SPLEND1D performance in the base case

C. The two-point model formatting

IV. DETACHMENT ONSET VIA DENSITY RAMP IN THE BASE CASE

A. Observation of a target ion flux roll-over and onset of a total pressure drop

B. Investigate process at play

V. EXPLORING THE ROLE OF FREE PARAMETERS

A. Carbon concentration

B. Neutral confinement time

C. Heat flux limiter

D. Bohm boundary condition

VI. ADVANCED STUDIES ENABLED BY SPLEND1D

A. Independent ion and electron temperatures

B. Time-dependent simulations

VII. CONCLUSION

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX A: COULOMB LOGARITHM

APPENDIX B: CONTRIBUTIONS OF MOMENTUM AND POWER SOURCE TERMS TO LOSS FACTORS

APPENDIX C: ADVECTIVE NEUTRAL MODEL

APPENDIX D: PRESSURE-DIFFUSION EQUATION

APPENDIX E: ON THE ROLE OF THE NEUTRAL MODEL

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only