The prediction of power fluxes and plasma-wall interactions impacted by MHD processes during ITER operation [disruption, Edge Localized Modes (ELMs), 3D magnetic fields applied for ELM control, etc.] requires models that include an accurate description of the MHD processes themselves, as well as of the edge plasma and plasma-wall interaction processes. In this paper, we report progress on improving the edge plasma physics models in the nonlinear extended MHD code JOREK, which has capabilities to simulate the MHD response of the plasma to the applied external 3D fields, disruptions and ELMs. The extended MHD model includes E × B drifts, diamagnetic drifts, and neoclassical flows. These drifts can have large influences, on e.g., divertor asymmetries. Realistic divertor conditions are important for impurity sputtering, transport, and their effect on the plasma. In this work, we implemented kinetic and fluid neutral physics modules, investigated the influence of poloidal flows under divertor conditions in the ITER PFPO-1 (1.8T/5MA) H-mode plasma scenario, and compared the divertor plasma conditions and heat flux to the wall for both the fluid and kinetic neutral model (in JOREK) to the well-established 2D boundary plasma simulation code suite SOLPS-ITER. As an application of the newly developed model, we investigated time-dependent divertor solutions and the transition from attached to partially detached plasmas. We present the formation of a high-field-side high-density-region and how it is driven by poloidal E × B drifts.

Simulation of transient phenomena driven by MHD in tokamak plasmas and their control have so far mostly focused on the MHD processes while the associated plasma-edge and plasma–wall interaction aspects of these processes have been addressed with specialized edge plasma codes in open-loop.2,3 This consists in taking the magnetic geometry as given by the MHD codes as fixed in time and using the edge-plasma codes to model the plasma in these conditions. This has two drawbacks: one is that transient modeling of edge-plasma and plasma wall interactions driven by MHD processes cannot be accurately done and the second is that the varying behavior of the edge plasma and plasma–wall interactions does not feedback on the MHD. The work presented here aims to address these two drawbacks by improving the edge plasma physics description in the non-linear extended MHD code JOREK, which has demonstrated capabilities to model the MHD and core plasma behavior during Edge Localized Mode (ELM) suppression by external 3D fields,4,5 ELMs6,7 and disruptions8 to cite few important examples.

JOREK is a fully implicit 3D non-linear extended MHD code for realistic tokamak X-point plasmas.1,9,10 The JOREK code has typically been developed and used for studying MHD instabilities such as ELMs and disruptions.6,11,12 The present version of JOREK includes some divertor physics in terms of sheath boundary conditions (Mach-one velocity and sheath transmission factors) and a simple neutral fluid model (reflection of outgoing ions as neutrals, ionization, recombination, and radiation). This model has been extensively applied for simulations of ELMs12 and ELM control by RMPs5 and by pellet triggering.13 

The aim of this paper is to describe, benchmark, and apply an improved physics model for the neutrals in the JOREK code. Previously, a kinetic particle framework14 has been added to the JOREK fluid MHD models. This particle model has been used to study the transport of Tungsten in the 3D ergodic fields due to ELMs.15 Kinetic particles have also been used to simulate the fast particle destabilization of TAE modes.16 This particle framework will be used to replace the neutral fluid model with a fully kinetic neutral model (not including molecular physics at this time).

In this paper, we first summarize the essentials of JOREK in Sec. II. Then, in Sec. III, the atomic neutral-plasma interaction physics of interest are derived and presented in Euler equations form (and the form as used for JOREK in Sec. III A) in a three-fluid model; both forms are particle, momentum, and energy conserving. Then, it is reduced with JOREK's assumptions to a two-fluid model. The focus here lies on volumetric ionization and recombination. Charge-exchange (CX) will only be used in the kinetic model. After the derivation, Sec. IV describes the implementation of the kinetic neutral model in JOREK. For the fluid neutral model, we refer to Fil et al.17 In Sec. V, we further present a benchmark with SOLPS-ITER simulations of a ITER PFPO-1 5MA H-mode scenario18,19 and the difference between utilizing fluid or kinetic neutral models within JOREK. We report on the energy conservation properties in Sec. V F. As a first application, the effect of poloidal E × B drifts on the heat and particle fluxes in this ITER scenario is investigated. This is investigated in the conduction-limited regime for both steady-state solutions in Sec. VI and in dynamic simulations in Sec. VII. In Sec. VIII, we present the dynamic transition toward plasma detachment and the formation of a High-Field-Side High-Density (HFSHD)-region driven by poloidal E × B drifts. Section IX is the conclusion.

The importance of E × B drifts in the SOL and divertor has been identified in Refs. 20–22. It is typically very difficult to measure the electric fields and the E × B drifts velocity directly in experiments. Hence, there is a need of understanding the effect of electric fields by means of modeling and simulations. Since the implementation of E × B drifts in codes such as SOLPS-ITER23 and SOLEDGE2D,24,25 the E × B drifts have attracted increasingly more attention in recent years; also in the MHD plasma modeling community.25–32 These edge plasma studies highlight the impact of drifts on divertor asymmetries and SOL flows in multiple machines. Changes to heat fluxes and to detachment behavior by drifts might impact the effective plasma operation space with acceptable power loads on the wall in future fusion devices. The interplay between E × B drifts and radiating impurities has been identified as a necessary condition for the onset of the experimental “detachment cliff” in DIII-D with a nonlinear dependence.33, E × B drifts further have been found to play a fundamental role in the development of the high-field-side high-density (HFSHD)-region in ASDEX-Upgrade.34,35 The HFSHD-region with radiating impurities can be a precursor for a X-point radiator plasma, which if not controlled can lead to disruptions. However, if controlled well can lead to ELM-suppressed regimes.36 

This section shortly summarizes the base MHD model in JOREK and the reduced MHD form to provide a background for the derivations of the plasma-neutral interactions derived in Sec. III. We start with the base MHD model, then state the reduced MHD assumptions and how to obtain the weak form of the reduced MHD equations. A more complete explanation of the model can be found in the study by Hoelzl et al.1 

The MHD model is formulated as a set of normalized equations for the time evolution of the magnetic potential ( A) Eq. (1), the fluid velocity ( v) Eq. (2), the mass density (ρ) Eq. (3), and the total pressure (p) Eq. (4). The normalization of the equations is performed such that parameters as the vacuum permeability μ0 and the Boltzmann constant kb do not appear in the equations,
A t = E Φ ,
(1)
ρ v t = ρ v · v p + J × B + · τ ¯ + S ̃ v ,
(2)
ρ t = · ( ρ v ) + · ( D ¯ ρ ) + S ̃ ρ ,
(3)
p t = v · p γ p · v + · ( κ ¯ T ) + ( γ 1 ) τ ¯ : v + ( γ 1 ) ( ( E + v × B ) · J ) + S ̃ p .
(4)
The magnetic field vector, B, and the current, J, are defined as
B = F ϕ + × A and J = × B .
(5)
Here, F ( ψ ) R B ϕ is the toroidal flux function. The electric field E is defined by Ohm's law including resistivity, η, and drift-ordered (diamagnetic) terms,
E = v × B + η ( J J * ) + F 0 δ * ρ ( p i p e ) .
(6)
Φ is the electric potential, τ ¯ the viscous stress tensor, D ¯ the diffusivity, and κ ¯ the heat conductivity. S ̃ ρ , S ̃ v, and S ̃ p in Eqs. (2)–(4) denote the source/sink terms. Due to the non-Eulerian form of the fluid equations, some extra terms appear in the source terms to keep the equations conservative. Source terms due to plasma–neutral interactions are derived and presented in Sec. III.
In order to reduce computational requirements, JOREK can also employ a reduced MHD model, which eliminates fast magnetosonic waves while retaining the relevant physics. Reduced MHD has fewer unknowns compared to full MHD, which decreases the computational costs and memory requirements for simulations. The reduced MHD model used in JOREK is derived following the ansatz-based approach.1 The ansatz considered here assumes that the time dependent part of the magnetic potential is dominated by the toroidal component. The ansatz for the magnetic field is deduced by approximating B ϕ as the vacuum F 0 / R toroidal field. After some calculations, we can obtain an expression of the poloidal component of the velocity ( v pol = ( e ϕ × v ) × e ϕ), which captures E × B effects and the ion diamagnetic drift velocity; and the parallel component of the velocity field v = v · B,
B = F 0 R e ϕ + 1 R ψ × e ϕ A · e ϕ = ψ ϕ , v = v pol + v = v E × B + v dia , i + v = R u × e ϕ ( δ * R / ρ ) p i × e ϕ + v B .
(7)

The velocity stream function u is defined as Φ / F 0. In this paper, the diamagnetic drift velocity is not used.

In JOREK, the equations are solved in the weak form using a finite element method. To obtain the weak form, we multiply Eqs. (1)–(4) with test functions as shown in Eqs. (8)–(10) and then perform integration by parts,
A t · B A * d V 1 R 2 ψ t A * d V .
(8)
With the reduced MHD assumptions, the weak form of the induction equation Eq. (8) is reduced from a vector potential to the scalar poloidal magnetic flux. Equation (2) is split into a parallel and poloidal momentum equation in the following way:
v t · B v * d V and v t · v E * d V .
(9)
For the pressure and mass continuity equations, we obtain
p t p * d V and ρ t ρ * d V .
(10)
As the current and vorticity are higher order derivative terms, these are introduced as separate parameters for numerical reasons. These two parameters and the previously presented reductions together give us seven equations with the seven variables: poloidal magnetic flux ψ, velocity stream function u, toroidal current density j, toroidal vorticity w, mass density ρ, temperature T = T i + T e, and parallel velocity v . There is also the possibility to run with a two-temperature model in JOREK, but is not used is this work. Any other time-evolved parameter (e.g., the neutral density for the fluid neutral model17) requires another equation with another unknown to solve.

The physics variables are Discretized using finite elements based on bicubic Bézier surfaces10 in the poloidal plane and by a Fourier series in the toroidal direction. These higher order elements provide G 1 continuity, i.e., continuity of the variables and their gradients in the R,Z plane. Even higher order continuity is available.37 The iso-parametric finite element grid is flux-aligned after running a Grad–Shafranov solver on the initial equilibrium. Furthermore, the grid can be fine-tuned to add resolution in regions of interest (e.g., the pedestal region and close to the strike-points).

For the time-evolution of the MHD fluid, the fully implicit time-evolution schemes, Crank–Nicholson or Gears38 (i.e., a BDF2 scheme), are utilized. These second-order accuracy schemes are unconditionally stable; therefore, the limitations of the time steps are determined by timescales of physical processes present in the simulation.

We describe two different neutral models used based on the same physics. First the fluid neutral model, and a new approach for the kinetic neutrals that are two-way coupled to the plasma fluid. For now, we are only looking at atomic species as a first step. The plasma-neutral reactions we take into account are electron-impact ionization, effective radiative recombination reactions, and exclusively for the kinetic neutrals, charge exchange (CX). For the neutral fluid model, neutral transport is described as a diffusive process with a spatially constant diffusivity.17 

We first derive a more general three-fluid model and thereafter reduce it with approximations to a two-fluid model and write it in the fluid equations form that JOREK uses. This is then the form implemented for the fluid neutral model. For the kinetic neutral model, this gives us a basis for what coupling terms are needed.

Let us first look at ionization and recombination reactions and their energy exchange in a stationary frame. Here, we look at the initial and final energy of the reaction participants. As a shorthand in this article, T ( i , e , n ) is written in terms of energy (eV). γ is the ratio of specific heats, ϕ ion is the ionization potential, and X is an arbitrary initial energy higher than the ionization potential. The reaction equation plus the energy terms are shown in the following equation:
e ( X ) + n 0 ( T n γ 1 ) ionisation i + ( m i m n T n γ 1 ) + 2 e ( m e m n T n γ 1 + X ϕ ion ) , e ( T e γ 1 ) + i + ( T i γ 1 ) recombination n 0 ( m i m n T i γ 1 + m e m n T e γ 1 ) + photon ( m e m n T i γ 1 + m i m n T e γ 1 + ϕ ion ) .
(11)

With ionization, the ion gains m i / m n times the neutral energy, the electron gains the rest of the energy but loses the ionization potential. In an electron fluid description, it does not matter how the energy is distributed between the two electrons. If the recombination results in the neutral being in an excited electronic state e + i + rec n * + photon, called dielectronic recombination, the excited state will cascade down to the ground (or metastable) state. This will result in most of the ionization potential energy being converted into photons. Generally, we may assume that all the radiatively recombined electrons cascade to the ground state.39 For implementation in JOREK, we should take the plasma fluid macroscopic velocity and kinetic energy into account as well. This results in an extra 1 2 m v 2 in the energy terms.

For mass sources and losses due to ionization and recombination, we use effective rate coefficients from the OpenADAS database.40 The reactions rates are then
Γ n rec = n e n i S rec , Γ ( i , e ) ion = n n n e S ion ,
(12)
where S ( rec , ion ) ( m 3 / s ) are the recombination and ionization rate coefficients (i.e., σ v ) and n n / i / e ( m 3 ) are the number densities of neutrals, ions and electrons. Γ ( m 3 / s ) is the number density source. Note that the rate coefficient depends on the electron density and temperature S = S ( n e , T e ). The ionization and recombination reactions involve one neutral, one electron, and one ion and, thus, it holds that Γ i ion = Γ e ion = Γ n ion and Γ n rec = Γ i rec = Γ e rec. Let us write just the partial time derivative of the continuity equation and the source terms of interest. The other terms do not influence the result. Throughout the paper, we use the subscripts ion and rec or ir to indicate whether the partial derivative is associated with ionization, recombination or both,
( n i t ) i r = Γ i ion Γ n rec , ( n e t ) i r = Γ e ion Γ n rec , ( n n t ) i r = Γ n ion + Γ n rec .
(13)
For the momentum, the sources and sinks for the different species arise formally from taking the first moments of ionization and recombination collision operators ( m i v C i ion d v Γ i ion m i v n) (Ref. 41). The momentum sources in the Eulerian form of the momentum equation are then
( m i n i v i t ) i r = Γ i ion m i v n Γ n rec m i v i , ( m e n e v i t ) i r = Γ i ion m e v n Γ n rec m e v e , ( m n n n v n t ) i r = Γ i ion m n v n + Γ n rec ( m i v i + m e v e )
(14)
with subscripts i, e, n denoting the property belonging to species ions, electrons, and neutrals.
For the energy equation, the sink and source terms correspond to the second moment of the ionization and recombination collision operators. As we are dealing with a “drifting” plasma, we need to include the macroscopic kinetic energy in the energy equation and the sink and source moments. Taking the second moment of these collision operators to obtain the energy exchange takes the form of 1 / 2 m α v 2 C α ion d v m α m n ( Γ α ion 1 / 2 m n v n 2 + Γ α ion T α γ 1 ) (Refs. 41 and 42). The energy sources in the Eularian form of the energy equation are then
( E α t ) = ( ( 1 2 m α n α v α 2 + n α T α γ 1 ) t ) ,
(15)
which gives us terms from the continuity and momentum equation in our energy/pressure equation. It is important to note that these terms are also multiplied by ( γ 1 ) as we are now solving for a pressure equation. In Secs. III C 1 and III C 2, the ionization and recombination terms are presented.

1. Ionization

For the ionization, the ions and electrons gain the energy of the neutral atom, based on their mass ratios. However the cost of ionization is only paid by the electrons. This gives us Eq. (16),
( E i t ) ion = Γ i ion m i m n T n γ 1 + 1 2 Γ i ion m i v n 2 , ( E e t ) ion = Γ i ion m e m n T n γ 1 + 1 2 Γ i ion m e v n 2 Γ i ion ϕ ion , ( E n t ) ion = Γ i ion T n γ 1 1 2 Γ i ion m n v n 2 .
(16)

2. Recombination

For recombination, both the ion and electron fluid lose their total energy per particle. However, the neutral only gains the energy based on the mass fractions of the ions and electrons as shown in Eq. (17). The rest of the energy (i.e., almost the entire electron energy) is emitted as photons and is lost from the system
( E i t ) rec = Γ n rec T i γ 1 Γ n rec 1 2 m i v i 2 , ( E e t ) rec = Γ n rec T e γ 1 Γ n rec 1 2 m e v e 2 , ( E n t ) rec = Γ n rec ( m i m n T i γ 1 + m e m n T e γ 1 ) + Γ n rec ( 1 2 m e v e 2 + 1 2 m i v i 2 ) .
(17)

To use the model in JOREK, we can reduce the equations and remove some terms by implementing assumptions made in JOREK. We reduce the ion and electron continuity equation by regarding only one plasma fluid. Let the center-of-mass density be ρ m i n i + m e n e, and ne = ni in the absence of impurities. We then assume the limit of m e 0. Then, ρ = m i n i and m i / m n = 1. This further results in the electron momentum equations being zero and all terms with electron mass disappearing.

Furthermore, note that in Eq. (A6), there is only a loss term for Te, and nothing for the energy to go to. This is the energy lost due to radiation. In the ADAS database the Continuum and line power driven by recombination is combined with Bremsstrahlung of dominant ions (PRB). We can remove Γ rec T e from electron energy equation if and only if we explicitly take radiation into account. We have to be careful here, as the PRB recombination is the sum of radiation from free-electron recombination, the dielectronic-cascade + Bremsstrahlung power coefficients. This includes the dielectronic-cascade power loss from an excited energy state to the ground state. This energy is not lost from the electron fluid energy equation. When we explicitly write the electron energy loss in terms of radiation power density, we should write P rad = ( P PRB Γ n rec ϕ ion ) to correct for the ionization potential energy. In terms of radiated power coefficients,
P rad = ( P PRB Γ n rec ϕ ion ) = L PRB n e n i S rec n e n i ϕ ion = ( L PRB S rec ϕ ion ) n e n i .
(18)
This all together results in Eqs. (19)–(21), our neutral physics interaction used in the JOREK models,
( m i n i t ) i r = m i Γ n ion m i Γ n rec , ( m n n n t ) i r = m n Γ n ion + m n Γ n rec ,
(19)
m i n i ( v i t ) i r = Γ i ion m i ( v n v i ) , m n n n ( v n t ) i r = m i Γ n rec ( v i v n ) ,
(20)
( ( n i T i ) t ) i r = Γ i ion T n + ( γ 1 ) Γ i ion 1 2 m i ( v n v i ) 2 Γ n rec T i , ( ( n e T e ) t ) i r = ( γ 1 ) ( Γ i ion ϕ ion ( L PRB S rec ϕ ion ) n e n i ) , ( ( n n T n ) t ) i r = Γ n ion T n + Γ n rec T i + ( γ 1 ) Γ n rec 1 2 m i ( v i v n ) 2 .
(21)
Note that there is no absorption term for radiation. This model assumes the plasma is transparent to radiation. In JOREK, there is the choice to use a one- or two-temperature model. As n e = n i = n, and partial pressure is additive, p p e + p i = n T e + n T i = n ( T e + T i ). We can add the electron and ion energy equations and write it as function of T T e + T i. In the benchmark presented in Sec. V, the single temperature model is used.

The kinetic neutral model is part of the kinetic particle framework developed in JOREK.14 The particle framework is very flexible for the implementation of many other physical models. In this section, we first present the general way of coupling particles to the fluid description of the plasma, then focus on fully kinetic neutrals.

The neutral particles utilize the same pusher as fully kinetic ions where the full orbits of the particles are followed using the well-known Boris method43 in toroidal geometry in the time varying electric and magnetic fields of the MHD fluid. Particles are followed in both ( R , Z , ϕ) space and in the local ( s , t , ϕ) space of the cubic Bezier finite elements. The ionization of the neutrals gives rise to a source of ions, momentum, and energy in the fluid. To add these sources to the fluid equations, the discrete particle distribution needs to be converted to the continuous, in space, description of the finite element basis. For example, the ion source S ρ ( x ) = m f ( x , v ) d v, where f, the distribution of ionized neutrals is represented in the finite element H i j ( s , t ) and Fourier basis: S ̃ ρ ( x ) = ijk p ρ , ijk H i j ( s , t ) H ϕ , k ( ϕ ), with pijk being the expansion coefficients. In the weak form, this yields a set of equations,
v * ( x ) S ̃ ρ ( x ) d x = v * S ρ ( x ) d x ; v * ( s , t , ϕ ) ijk p ρ , ijk H i j ( s , t ) H ϕ , k ( ϕ ) d x = n v * m w n δ ( x x n ) .
(22)
Since it depends only on the finite element geometry, the LHS of the system of equations needs to be factorized only once. The system needs to be solved at every fluid time step. The time step of the kinetic particles is determined by the underlying physics time-scales in order to get an accurate spatial distribution of the source. The particles are considered superparticles with a certain weight wn. At every ionization event this weight is reduced until it reaches a critical lower value at which it is removed. For the conservation properties of the combined fluid–particle system, it has been found essential to collect the contributions to the right-hand side particle sums at every particle time step. To reduce the noise in the finite element representation (or as regularization), the following equation is solved instead of Eq. (22),
v * ( x ) ( 1 λ 2 ζ 4 ) S ̃ ρ ( x ) d x = v * S ρ ( x ) d x ,
(23)
where λ and ζ are the regularization coefficients and v * the test function of the finite element method (weak form). The regularization coefficients can be chosen with different values for the n = 0 harmonic (as opposed to higher harmonics) and in the perpendicular and parallel directions. This method of projection allows us to project an arbitrary number of any moments of the particles while only calculating one matrix factorization.
Plasma–wall interaction is essential for a high-recycling plasma, modeling detachment and sputtering of the wall material. A general sputtering module has already been implemented14 based on the Eckstein sputtering coefficient.44 Both the fluid and kinetic particles can reflect and sputter on the wall. For recycling (i.e., plasma hitting the boundary and being re-introduced in the domain as a neutral), the TRIM-code database45 is used to determine the probability of fast reflection or thermal desorption based on the electron temperature and wall material. For the plasma fluid, the boundary values are known and we can integrate the total reflected or sputtered flux on the boundary,
N α refl / sputter = Y ( E 0 ( T ) , θ ) Γ plasma d A .
(24)
A fixed number of superparticles with combined weight N α refl / sputter is sampled where the location is weighted by the yield times plasma flux. The energy of a reflected neutral is either a fraction of the incoming sampled particle from the plasma, or a fixed value for thermal desorption based on an estimated fixed wall temperature. Kinetic particles impacting the boundary are evaluated and recycled individually.

As the incident plasma fluence is only defined based on the MHD fluid time stepping scheme, sputtering and reflection events are performed only once per fluid time step. This means, after a fluid time step, superparticles are introduced in a sputtering event, then in the kinetic particle loop, they are evolved per kinetic time step just like the other superparticles that already existed. Although for the kinetic particles this a pulsed particle source, as it is coupled back to the plasma fluid only every fluid time step, it does not result in oscillating or pulsed solutions.

Depending on the type of kinetic particles, different moments and coupling schemes are necessary to couple the right information to the JOREK fluid plasma. Evaluating the plasma–neutral interaction physics as written in Sec. III, three moments are needed to describe the coupling, mass density, momentum and energy. Here, all the effects of the plasma-neutral interactions are summed to obtain one source term per equation; in the conservative form the source are written as S ρ , S v, and SE. This results in the following source terms in the JOREK form,
( ρ t ) = S ̃ ρ = S ρ , ρ ( v t ) = S ̃ v = S v v S ρ , ( ( ρ T ) t ) = S ̃ E = ( γ 1 ) ( S E v · ( S v ) + 1 2 v 2 ( S ρ ) ) ,
(25)
which is derived in the same manner as in Eq. (15).

Although the moments of kinetic particles are coupled as explicit terms, they are multiplied by implicit MHD terms [e.g., v in Eq. (15)] to maintain the highest order accuracy. In these coupling schemes, when kinetic particles are the dominant species, the accuracy of the time-evolution scheme might be of a lower order. Since we have knowledge on the dependency of rate coefficients on the plasma density and temperature (ne, Te), it is possible to improve the order of accuracy by making the coupling terms (semi-) implicit in future work.

As the kinetic particle framework is a Monte Carlo model, we need to sample from probabilities to determine which interaction(s) will happen. Ionization, CX, and line radiation are dependent on the kinetic neutrals population and it should be calculated based on the kinetic time step. Phenomena that create neutrals from the fluid plasma are only well-defined over an entire fluid time step ( Δ t fluid) and are not an interaction that can be based on the existing superparticles in the plasma domain. We, thus, use separate routines to initialize the superparticles once per fluid time step for plasma–boundary interaction, puffing, and recombination. However, in principle, it is possible to initialize these phenomena more than once per tfluid. Section IV D 1 describes the implementation of the plasma–neutral interactions.

1. Ionization

The ionization probability P ion is calculated from the density and temperature of the background plasma at the location of the superparticle. If a uniformly distributed pseudo-random number U rng [ 0 , 1 ] is lower than P ion, then the weight of the superparticle is updated according to the probability,
P ion = ( 1 e S ion ( n e , T e ) n e Δ t kinetic ) if U rng [ 0 , 1 ] P ion , then w n t + Δ t kinetic = w n t P ion .
(26)
The ionized weight w n ( 1 P ion ) is then the ionization source term for the plasma fluid with its according momentum and energy exchange,
S ρ ion = m i w n ( 1 P ion ) , S v ion = S ρ ion v or S v ion = S ρ ion B · v , S E ion = S ρ ion ( 1 2 v n 2 ϕ ion / m i ) ,
(27)
for the one temperature model.

2. Charge-exchange

Since CX between ionized hydrogen and neutral atomic hydrogen is a symmetric collision, there are no radiation losses. Based on the CX reaction rate probability, P C X, we can sample an ion from the background plasma shifted Maxwellian with temperature Te and macroscopic drift velocity v drift and exchange its momentum with the superparticle,
P C X = ( 1 e S C X ( n e , T e ) n e Δ t kinetic ) if U rng [ 0 , 1 ] P C X , then v n = N rng T e m i + v drift ,
(28)
where N rng ( 3 ) is a normally distributed random number, and v drift ( 3 ) the plasma macroscopic drift velocity. The plasma sources due to CX then become
S ρ C X = 0 , S v C X = w n ( v n old ( N rng T e m i + v drift ) ) , S E C X = 1 2 w n ( ( v n old ) 2 ( v n new ) 2 ) .
(29)

3. Line radiation

Line radiation due to excitation of atoms is an underlying mechanism, which the code does not resolve. Therefore, only the effective radiated power is calculated at a corresponding temperature,
P rad , line = n e w n S PLT ( n e , T e ) Δ t kinetic ,
(30)
where S E = P rad. There are no further mass density or momentum losses involved.

4. Recombination

For recombination of the background plasma into kinetic neutrals, a different strategy is devised. As the neutrals from recombination come from everywhere in the plasma, we need to account for the recombined mass and all lost energy. Therefore, all the losses of the plasma due to recombination are calculated in the fluid, and then sampled from, per finite-element, to initialize superparticles. The losses from the plasma are then identical as the fluid neutral model and as in Sec. III. From each finite-element, the recombined weight, lost momentum and energy are integrated over time Δ t fluid. The velocity given to the superparticle is the lost plasma momentum divided by the recombined weight as shown in Eq. (31),
w n = 1 m i elm ρ rec d V = elm S rec n e n i Δ t fluid d V , v = rec ρ rec v drift d V m i w n .
(31)
This strategy ensures mass and energy conservation at the cost of initializing at least one superparticle per finite-element. Specifically, at hot regions of the plasma this would be a negligible superparticle. This, however, does not cost the simulation much computational effort as these superparticle will be ionized in one kinetic time step.

The location of the recombined superparticle is sampled in the poloidal plane from the local element space (s,t) with two uniform distributed random numbers U rng [ 0 , 1 ], and in the toroidal direction, ϕ = ϕ plane + 1 2 U rng [ 0 , 1 ] Δ ϕ to place it in the sector associated with the poloidal plane with coordinate ϕ plane and distance between planes Δ ϕ.

To test the implementation of the neutrals code development, JOREK simulations are compared to SOLPS-ITER simulations of an ITER (pre-fusion power operation) PFPO-1 H-mode scenario without impurities (first published in Park et al.18) In this section, we compare the plasma power and peak heat flux to the divertor walls as a function of the upstream density (i.e., the outer mid-plane separatrix density). Here, the simulations are performed for a range of different neutral fueling rates in the divertor. The details of the scenario are B T = 1.85 T , I p = 5 M A, Pheating = 20 MW, D = 0.3 m2s−1, χ i , = χ e , = 1.0 m2s−1. Parallel heat transport is based on the Braginskii closure.46,47 The fluid neutral diffusivity was chosen to match the SOLPS-ITER neutral density decay from the strike point following the field lines at low puffing rates. This resulted in D n 2 × 10 2 m2s−1. For the kinetic neutrals, such an assumption is not needed. To reach a quasi-steady-state, JOREK simulation were run between 11 and 55 ms depending on the neutral gas puffing rate.

Although SOLPS-ITER solves plasma transport differently and has hydrogen molecules included, we still expect to see similar behaviors with JOREK. It is known that with temperatures under 1 eV, molecules are important for the correct detachment behavior and the exact behavior of detachment front.18 However, to benchmark as a first step let us look at the power incident on the inner and outer divertor targets and the peak heat fluxes.

Although more advanced options are available in JOREK, to implement and study neutral physics we do not use everything at once. Here, the visco-resistive reduced MHD is used as model where the time-evolved parameters are ψ , j , u , w , ρ , T = T i + T e , v , and ρn, the neutral density for the addition of fluid neutrals.

Figure 1 shows the JOREK FE (finite element) grid extending up to the first wall/divertor and the SOLPS-ITER plasma and kinetic grids. SOLPS-ITER does not simulate the full core, and the plasma does not reach the wall except for the divertor. The red lines are possible puffing and pumping surfaces. The JOREK grid does include the core. SOLPS-ITER uses two separate grids, one for the plasma fluid, one for the kinetic particles. The JOREK kinetic particles exist in real space (R, Z, ϕ) and local finite element coordinates.

FIG. 1.

(Left) Representation of the JOREK finite element grid cells. The grid is extended up to the first wall. The dome structure and below are not included in the grid. (Right) SOLPS-ITER plasma grid in blue, Eirene kinetic particle grid in green, extended up to the first wall, dome structure included. (SOLPS-ITER grid figure altered from Ref. 19).

FIG. 1.

(Left) Representation of the JOREK finite element grid cells. The grid is extended up to the first wall. The dome structure and below are not included in the grid. (Right) SOLPS-ITER plasma grid in blue, Eirene kinetic particle grid in green, extended up to the first wall, dome structure included. (SOLPS-ITER grid figure altered from Ref. 19).

Close modal

As the fueling rate scan in the study by Park et al.18 was performed with SOLPS-ITER without drifts, we have reduced the JOREK models here to obtain driftless behavior where also the magnetic field is frozen. For this, the time-evolution of the parameters ψ , j , u, and w is removed. This fixes the magnetic field and removes the E × B flows.

For this benchmark, the neutral gas puffing location was chosen to be in the divertor. In steady-state and without impurities, the puff location has no significant impact in the near-SOL and divertor region.18 However, this could have an effect on time-dependent behavior and the effectiveness of the fueling the plasma.

The pumping surface in the JOREK grid is in the private region boundary. Here, diffusively transported neutrals are pumped out.

Figure 2 shows the resolution of the FE grid on top of the plasma density. To resolve the sharp (density) gradients close to the strike points, a resolution on the order of 1 mm is need in the direction along the field. Insufficient resolution will alter the divertor solution of the simulation. The figure shows the gradients to be well resolved. Monte Carlo noise of the projections of the kinetic particles depends on the number of superparticles per finite element. A finer grid, thus, requires more superparticles for similar noise levels.

FIG. 2.

Plasma density (a.u.) shown on the finite elements. The spatial resolution along the field 0.74 mm, cross field 8.1 mm. Strong gradients along the field over 5.7 mm can be smoothly resolved.

FIG. 2.

Plasma density (a.u.) shown on the finite elements. The spatial resolution along the field 0.74 mm, cross field 8.1 mm. Strong gradients along the field over 5.7 mm can be smoothly resolved.

Close modal

In this subsection, we present the results of the fueling scan done in JOREK with fluid neutrals compared to the SOLPS-ITER results. Only the results without drifts are directly comparable between JOREK and SOLPS-ITER. Here, only a few metrics are compared as a global overview.

First, the heat load to the wall as function of the neutral gas puffing rate was investigated in Ref. 48. This showed the dependence of the total plasma heat load on the fueling rate are qualitatively in good agreement.

Figure 3 shows the plasma heat load to the divertor as function of the outer mid-plane separatrix density nupstream. Although JOREK obtains the same power ranges; in a stationary state, nupstream at a given fueling rate is lower in JOREK simulations compared to SOLPS-ITER. This initial drop in nupstream is related to the difficulty of keeping the shape of the initial density profile, which is assumed ad hoc rather than consistent with gas fueling.

FIG. 3.

The inner, outer, and total divertor plasma heat load as function of the outer mid-plane separatrix electron density for SOLPS-ITER and JOREK with fluid neutrals (fluid). All results in this figure are without drifts.

FIG. 3.

The inner, outer, and total divertor plasma heat load as function of the outer mid-plane separatrix electron density for SOLPS-ITER and JOREK with fluid neutrals (fluid). All results in this figure are without drifts.

Close modal

The use of diffusive fluid neutrals with a spatially constant Dn is not sufficient to provide the correct plasma mechanics. However the addition of fluid neutrals allows JOREK simulation to transition from a typically convective (sheath-limited) regime to the conduction limited regime. This relatively simple fluid neutral model captures enhanced energy losses with increased neutral gas puffing, which leads to decreased divertor heat loads. Given the simplicity of the model, we can consider this as reasonable agreement. However, further study for this model is needed to identify the reason for the discrepancies with the SOLPS-ITER heat loads at the low end of upstream densities ( n upstream 1.5 × 10 19 m 3).

The fueling scan using the JOREK model with kinetic neutrals compared to the SOLPS-ITER results is shown in Fig. 4. Only the results without drifts are directly comparable between JOREK and SOLPS-ITER. As opposed to fluid neutrals, kinetic neutrals do not need to assume a diffusivity for the neutral transport. In the kinetic neutral model, the neutrals do have a velocity and can exchange momentum and energy with the plasma via a combination of CX, ionization and recombination.

FIG. 4.

The inner, outer, and total divertor plasma heat load as a function of the outer mid-plane separatrix electron density for SOLPS-ITER and JOREK with kinetic neutrals (kinetic). All results in this figure are without drifts.

FIG. 4.

The inner, outer, and total divertor plasma heat load as a function of the outer mid-plane separatrix electron density for SOLPS-ITER and JOREK with kinetic neutrals (kinetic). All results in this figure are without drifts.

Close modal

Figure 4 presents the same scenarios as in Sec. V B. At the lower- to mid-range of the upstream density ( n upstream 2 × 10 19 m 3), there is good agreement in total power of JOREK with SOLPS-ITER. The inner–outer divertor power balance is slightly more symmetric than SOLPS-ITER. Beyond n upstream 2 × 10 19 m 3, the SOLPS-ITER plasma transitions from high-recycling to a fully detached divertor plasma. While this happens nupstream starts to saturate. In the no-drift JOREK simulation, the transition to detachment occurs at about 10% higher nupstream than SOLPS-ITER. The transition as function of nupstream is less sharp in JOREK. The slope ( Power / n upstream) is decreasing for higher nupstream This does indicate a similar but less strong form of density saturation from n upstream 2.4 × 10 19 m 3.

Since our kinetic neutral model now only consists of atoms, and neutral–neutral collisions are neglected, one can expect the comparison with SOLPS-ITER to diverge when detaching the plasma from the divertor. We suspect the absence of molecules and neutral–neutral collision together with the exact geometry of the simulation, which is not reproduced in JOREK, can explain the difference in our results from the ones of SOLPS-ITER.

Neutral–neutral collisions combined with the geometry of our simulation grid act to redistribute neutrals in the divertor region. The lack of neutral–neutral collisions allow hot neutrals to be ballistic and reach the wall with a lot of energy. Neutral–neutral collisions would disperse this energy through a lot of neutrals before reaching the wall. The pressure created in the neutrals then drive transport under the divertor dome. This redistribution of neutral particles can effectively change the power distributed between the inner and outer divertor.49 In our simulations, there is no grid under the divertor dome.

For plasma detachment, molecules are key around a few eV. They play an important role in getting the detachment behavior right.18,50 The molecules create strong (momentum) losses for the plasmas at these low temperature; this creates a stronger ion flux rollover and decreases heat load at the wall. Stronger momentum losses along the field lines in the divertor will alter the transport behavior. This might result in stronger upstream density saturation as fueling from the divertor becomes less efficient.18,50 The implementation of molecules in JOREK is left for future work.

A key divertor metric that is important to determine the compatibility of the plasma power fluxes with the plasma facing components at the divertor is the peak heat flux ( q p k , ). Figure 5 shows q p k , on the inner and outer divertor target. Note that for both SOLPS-ITER and JOREK a toroidally averaged wall is considered (i.e., the tilting and beveling of the actual divertor tiles are not taking into account). To include the shape of the monoblocks, a corrected angle can be taken into account.51 Generally, we observe JOREK having a higher q p k , . For the same plasma heat load, this indicates that JOREK profiles are narrower than SOLPS-ITER. The slope ( q p k , / n upstream) is similar between the two codes. However, as with the total power, the reduction of the peak heat flux with an increased fueling rate in JOREK occurs at higher upstream densities compared to SOLPS-ITER. In both simulations at high nupstream, the inner and outer target peak heat fluxes become the same.

FIG. 5.

The inner and outer, peak plasma heat fluxes as function of the outer mid-plane separatrix electron density for SOLPS-ITER and JOREK with kinetic neutrals (kinetic). All results in this figure are without drifts.

FIG. 5.

The inner and outer, peak plasma heat fluxes as function of the outer mid-plane separatrix electron density for SOLPS-ITER and JOREK with kinetic neutrals (kinetic). All results in this figure are without drifts.

Close modal

Generally, the inner target q p k , is lower than the outer target. At the low end of nupstream, JOREK has a higher inner q p k , . Why there is a strong difference between the behavior between JOREK and SOLPS-ITER in the inner target peak heat flux at low upstream density is not entirely clear. For JOREK, the ion flux is higher on the inner target, and the density is not yet high enough to really decrease the temperature. So the high ion flux with higher temperature combines in a higher heat flux.

Differences in the peak heat fluxes could have several causes. For the same power, a high peak heat flux indicates JOREK has a lower effective cross field heat diffisivity κ , effective. There are two notable differences in the simulation that could account for this. For this benchmark, JOREK employs the one-temperature model, while SOLPS-ITER used a two-temperature model, which is important since the same gradients along the field line are implicitly assumed for ion and electrons, which is not the case when two temperatures are considered. Furthermore, molecule–plasma interaction can provide efficient momentum-exchange channels that effectively enhances the cross field particle and energy transport.

We have discussed the physics results of the used fluid neutral and kinetic neutral model. The kinetic neutral model obtains quantitatively better results. The fluid neutral model with a spatially constant diffusion coefficient does not suffice to obtain the desired divertor solutions. Although the fluid model is much simpler to implement, i.e., add one extra equation and variable to the system of equations, it is typically not computational more advantageous. On typically four nodes, with 36 cpus per node, the kinetic neutral model using 10 6 simulation particles using a time step of 2 × 10 8 s, is as fast using the fluid neutral model. This means that in axisymmetric simulations, there is no computational penalty for using the kinetic neutrals over the fluid neutrals. For bigger computational problems (e.g., 3D), the kinetic neutral model scales almost ideally with the number of compute nodes and scales more favorably than the fluid neutral model. To reach steady-state, it takes 9000 10 000 cpuh, which equates to 60 70 h wall-time on said four nodes.

The initial conditions of some of the variables are defined by (re-)constructed magnetic equilibrium. This then gives us an initial profiles for ψ, ρ, and j. From experimental profiles (or from other codes in our case), the pressure p is separated into the mass density ρ and temperature T. If experimentally the parallel velocity is known in the core and edge, we can use it as initial condition for v . These initial conditions are all functions of the flux, asymmetries will arise naturally. The potential and vorticity u and w are initialized to zero but arise naturally within a few time steps.

For the boundary conditions, in this case, we use a fixed magnetic boundary, ψ thus has a Dirichlet boundary condition arising from the (re-)constructed equilibrium. The potential and vorticity u and w are Dirichlet conditions set to the value of 0. ρ has a Neuman boundary condition. The boundary conditions for the v and the parallel heat flux (i.e., a function T) are described by the usual magnetized sheath boundary conditions with a sheath heat transmission factor.52 

The kinetic particles Boris pusher used for advancing the particles in time preserves the particle energy conservation down to machine precision.43 However, for kinetic-fluid interaction not all properties are perfectly numerically conserved. The conservation further depends on the chosen time-stepping scheme. The coupling of moments of the kinetic particle distribution result in explicit source terms in the MHD equations. The time stepping of the JOREK MHD fluid utilizes the fully implicit Gears scheme38 (i.e., a BDF2 scheme). The Gears scheme does not conserve energy exactly but is second order accurate in the time step.

The recombination process is calculated from the finite elements of the plasma fluid. This results in an average relative error in momentum 10 4. However, the relative error in energy is 10 3. The ionization, CX, and radiation errors will depend on the error in the projection to the finite elements. The error should scale with the amount of superparticles per finite element. As the neutral–plasma interaction is localized on the centimeter scale, it is non-trivial to assess this a priori. However, most neutrals “live” close to the strike points where plasma–neutral interaction is high. To ensure adequate statistical coverage of the plasma–neutral or plasma-wall processes, simulations were performed with 10 6 10 7 superparticles with an average weight of 5 × 10 13 5 × 10 14 atoms. The weight of each superparticle is self-governed by the physical processes.

To quantify the error in energy conservation of the combined fluid-kinetic neutrals model, the power balance of one of the simulations is shown in Fig. 6. It presents the entire thermal and kinetic power balance. In this 2D equilibrium, there is no change in the magnetic energy and no Ohmic heating losses are used to feed magnetic energy into thermal energy. The x axis shows the simulation time steps of the JOREK plasma fluid. Each time step equals 2.246 × 10 6 s. In this figure, volumetric and surface processes are combined. Essentially, the equation that is compared is
ϵ + V P heating d U d t d V = V P rad + P plasma neutral d V + A P plasma wall + P neutral wall d A .
(32)
FIG. 6.

Power balance between energy sources and sinks for the steady-state phase of the simulation. The results are from a high-recycling but still attached simulation.

FIG. 6.

Power balance between energy sources and sinks for the steady-state phase of the simulation. The results are from a high-recycling but still attached simulation.

Close modal

The volumetric losses of the plasma and neutrals are the plasma–neutral (pn) interaction and corresponding radiation losses. The surface losses at the boundary are the plasma–wall interaction and neutral–wall interaction. The effective plasma heat load on the wall is the outgoing plasma heat flux minus the energy retained after wall-assisted recombination. The effective neutral heat load is the outgoing kinetic neutral heat flux minus the energy the neutrals retain after colliding with the wall. P heating d U d t is how much energy is supplied to the plasma minus the change in internal energy. In steady-state, the power balance is conserved with a relative error ε 10 3 10 2. The error depends on the time step size and the amount of simulation particles. For these applications this is a sufficient accuracy.

In the benchmark with SOLPS-ITER, V, we focused solely on JOREK results without drifts. By design, the JOREK MHD models include E × B flows as an essential ingredient. In this section, we compare standard JOREK reduced MHD model (including the evolution of magnetic flux, and vorticity) with JOREK simulations without drifts.

Figure 7 presents the results of the same scenarios as in Secs. V B and V C. At the lower end of nupstream, the drifts do not influence the solutions in the divertor. With increasing nupstream, the E × B drifts start redistributing the power balance such that more power goes to the outer divertor target. Simultaneously the inner target heat load decreases faster than in the simulations without drifts. Drifts cause a generally stronger decrease in power with increasing nupstream. Specifically, the E × B flow from the outer target to the inner target through the private region causes the power to drop to the inner target. Around n upstream 2 × 10 19 m 3, the plasma is already detached from the inner target. While detachment behavior (roughly T e , target < 5 eV) without drifts starts at much higher nupstream.

FIG. 7.

The inner, outer, and total divertor plasma heat load as a function of the outer mid-plane separatrix electron density for JOREK without drifts and JOREK with drifts. All results in this figure are with kinetic neutrals.

FIG. 7.

The inner, outer, and total divertor plasma heat load as a function of the outer mid-plane separatrix electron density for JOREK without drifts and JOREK with drifts. All results in this figure are with kinetic neutrals.

Close modal

Around n upstream 2 2.2 × 10 19 m 3, the upstream density rises with increased fueling while the power to the divertor almost remains constant. This is a phase-like transition of the plasma and coincides with the formation of a so-called high-field-side high-density region (HFSHD-region). The dynamic formation of the HFSHD-region is presented and explained in Sec. VIII.

Figure 8 presents the steady-state peak heat fluxes (qpk) as function of nupstream. Although at low nupstream the total power is the same with and without drifts, the heat flux profiles have changed. At the low nupstream, the inner target peak heat flux q p k , I T is increased while q p k , O T is decreased. q p k , I T decreases more quickly due to the E × B drifts. Counter-intuitively with increasing nupstream q p k , O T increases up until n upstream 1.5 × 10 19 m 3. At higher nupstream, q p k , O T decreases similar to its no-drift counterpart until the phase-like transition of the plasma.

FIG. 8.

The inner and outer, peak plasma heat fluxes as a function of the outer mid-plane separatrix electron density for JOREK with drifts (w/drift) and JOREK without drifts (w/o drifts).

FIG. 8.

The inner and outer, peak plasma heat fluxes as a function of the outer mid-plane separatrix electron density for JOREK with drifts (w/drift) and JOREK without drifts (w/o drifts).

Close modal

The ion flux through the boundary in JOREK is Γ ρ v · n ̂, where the boundary parallel velocity is set to the sound speed because of the Bohm sheath criterion, v = c s. This then gives us Γ ρ T 1 / 2. Figure 9 shows the (quasi-) steady-state ion flux profiles over the inner and outer divertor target, with and without drifts, for a range of outer mid-plane plasma densities, nupstream. The upper row shows JOREK results (with drifts) and the lower row in the figure shows ion flux results of JOREK without drifts. The results from Fig. 9 correspond to the points in Figs. 7 and 8. As the plasma responds differently to the same amount of gas puffing when the drifts are removed, the comparison cannot easily be done for exactly the same nupstream.

FIG. 9.

Quasi-steady-state ion flux profile for the runs corresponding to Fig. 7 on the inner and outer target vs plasma density as function of the distance to the separatrix. The coloring indicates the corresponding upstream outer mid-plane plasma density. (Top left) Inner divertor target with drifts. (Top right) Outer divertor target with drifts. (Bottom left) Inner divertor target without drifts. (Bottom right) Outer divertor target without drifts.

FIG. 9.

Quasi-steady-state ion flux profile for the runs corresponding to Fig. 7 on the inner and outer target vs plasma density as function of the distance to the separatrix. The coloring indicates the corresponding upstream outer mid-plane plasma density. (Top left) Inner divertor target with drifts. (Top right) Outer divertor target with drifts. (Bottom left) Inner divertor target without drifts. (Bottom right) Outer divertor target without drifts.

Close modal

The total plasma heat load to the wall (as shown in Fig. 7) shows that at the low end of nupstream, the drifts do not have a big influence. Yet/However with increasing nupstream, the power redistributes differently between the IT and OT. At higher n upstream 1 × 10 19 m 3, the location of and height of the peak ion flux is different. With increasing nupstream, we see the qualitative behavior more diverging, i.e., the location and height of the peak, single- vs double-peaked, the rollover and detachment.

In all results (with and without drifts, and at the IT and OT), the ion flux moves upwards away from the strike point (to higher s s sep) with increasing nupstream as partial detachment sets in. The E × B drift compresses the ion flux on the OT, where the profile remains closer to the separatrix than its “no-drift” counterpart. On the inner target, the E × B drift pushes the profiles upwards, further away from the separatrix.

The no-drift simulations show a clear double-peaked profile. Observing the plasma profiles 1 cm from the wall, the left peak is due to a combination of peaked ne and low v , and the right peak is at lower ne, but v sharply increases. This accelerates the plasma differently to the expectations from the simple sheath and we obtain a double-peak ion flux profile. In the steady-state profile with drift, we do not see a strong double-peaked profiles, only more flattened profiles, which is due to the same phenomena, but less strong. Furthermore, the “drift” IT receives significantly lower ion fluxes than the no drift result. The asymmetry between drift OT and IT is due to the poloidal E × B transport from the OT to the IT and from the IT upwards. The upward E × B flows the IT upwards are an important phenomenon while detaching the plasma, which will be discussed in more depth in Sec. VIII.

The transition from the high-recycling regime toward detachment is typically marked with the “rollover” of the ion flux (i.e., decreasing ion flux with increasing nupstream). The “no drift” rollover is fairly symmetrical between the OT and IT. However, as the E × B drift drives asymmetry, it creates a more asymmetrical rollover. On the drift IT rollover begins around n upstream 1.8 × 10 19 m 3 while the OT start seeing rollover around n upstream 2.15 × 10 19 m 3.

Remarkably, at around n upstream 2.3 × 10 19 m 3 in the drifts simulation, there are two different profiles at almost equal nupstream ( 2.28 × 10 19 and 2.3 × 10 19 m 3). This coincides with the formation and stabilization of the HFSHD-region. This non-linear transition has a strong effect on the power distributed to the IT and OT and on the fueling of the plasma from the divertor. This is discussed in more depth in Sec. VIII. During this transition, the ion flux profile first moves downwards (lower s s sep) then moves 1 cm upwards again. This is due to the strong E-field created by the HFSHD-region. This transition only occurs with drifts. Simulation without drifts Fig. 9 show a much smoother behavior on the IT and OT heat fluxes with increasing density.

The plasma heat flux through the boundary in JOREK, q γ sheath ρ T v . With v = c s, this gives q γ sheath ρ T 3 / 2. Figure 10 shows the (quasi-) steady-state heat flux profiles over the inner and outer divertor targets, with and without drifts, for a range of outer mid-plane plasma densities, nupstream. Again, this is only the plasma heat flux. The upper row shows JOREK results (with drifts) and the lower row shows heat flux results of JOREK without E × B drifts. The results from Fig. 10 correspond to the points in Figs. 7–9.

FIG. 10.

Quasi-steady-state heat flux profile for the runs corresponding to Fig. 7 on the inner and outer targets vs plasma density as function of the distance to the separatrix. The coloring indicates the corresponding upstream outer mid-plane plasma density. (Top left) Inner divertor target with drifts. (Top right) Outer divertor target with drifts. (Bottom left) Inner divertor target without drifts. (Bottom right) Outer divertor target without drifts.

FIG. 10.

Quasi-steady-state heat flux profile for the runs corresponding to Fig. 7 on the inner and outer targets vs plasma density as function of the distance to the separatrix. The coloring indicates the corresponding upstream outer mid-plane plasma density. (Top left) Inner divertor target with drifts. (Top right) Outer divertor target with drifts. (Bottom left) Inner divertor target without drifts. (Bottom right) Outer divertor target without drifts.

Close modal

Figures 7 and 8 shows that although the total heat load to the wall is not strongly influenced by drifts at the low end of nupstream, the peak heat fluxes are different. As the heat flux q ρ T 3 / 2 and the ion flux Γ ρ T 1 / 2, it is expected to see stronger changes with heat flux profiles with increasing nupstream.

The no drifts results in Fig. 10 show a similar response on the IT and OT. q n o drift , I T decreases faster with increasing nupstream than q n o drift , O T; however, qualitatively they behave similarly, and qpeak moves progressively upwards (higher s s sep). As with with the ion flux in Fig. 9, the E × B drifts compresses and pushes the OT heat flux profile downwards (lower s s sep) and pushes the IT heat flux upwards. Due to the T 3 / 2 dependency and higher temperature in the far-SOL, the heat flux peak is always at higher s s sep than the ion flux.

At n upstream 2 × 10 19 m 3 , q I T , n o drift is still significant, however, q I T , drift has already completely collapsed and is in a detached state (taking detachment as T e , target 5 eV). Due to the E × B drifts, the heat flux profiles on the IT and OT behave strongly asymmetrically. The earlier q I T , drift decay seems to be a result of two factors, both caused by E × B drifts. E × B drifts cause a different distribution of the transport from plasma core through the edge to the wall. Second, E × B drifts cause transport from the OT to IT via the private flux region. The higher density at the IT lowers the temperature leading to the low heat fluxes on the IT. Double-peaked structures are seen, similarly to those in the ion flux, due to a combination of regions with high ρ and low T and low ρ and high T.

The no drift results show a steady progressive decay of the heat flux. There is no sign of phase or regime change into detachment. With drifts, the results are asymmetric and around n upstream 2.3 × 10 19 m 3, the heat flux profile moves first downwards (to lower s s sep) and upwards. The same as for the ion flux results (Fig. 9).

So far, 0D and 1D profiles have been presented. To understand how these changes occur due to E × B flows, let us look at the 2D poloidal flow profiles. 2D flows and divertor geometries have been extensively studied in experiments and simulations to design next generation divertors (such as the ITER divertor).53 Even without drifts, neutral recycling associated with the vertical divertor geometry has a deep impact on plasma flows. Herein strong recycling leads to flow reversal away from the target in the separatrix region while it causes strong flows in the outer SOL. Here, we will focus more on the effect of drifts on top of the already existing flows due to the divertor geometry. Figure 11 (left) shows the electron density profile ne, with on top of that, the 2D poloidal flows visualized with colored arrows. The color represents the relative importance of the E × B flows compared to the rest of the poloidal flows (the parallel component in the poloidal plane), i.e., E × B factor | v E × B | | v pol v E × B |. The length of the arrow represents magnitude of v pol.

FIG. 11.

2D electron density profiles for an high-recycling ( n upstream 1.95 × 10 19 m 3) case for ITER PFPO-1. (Left) 2D poloidal flows visualized on top of the 2D electron density profile. The color of the quivers shows the importance of the E × B term. The importance is here defined as E × B factor | v E × B | | v pol v E × B |; the magnitude of the E × B velocity divided by the magnitude of poloidal velocities excluding the E × B velocity. The thin white line is the separatrix. (Right) Electric potential (white) contour on top of the electron density profile. E × B drifts are along the (white) potential contour lines. The red contours are the electron temperature Te for 1, 10, and 50 eV.

FIG. 11.

2D electron density profiles for an high-recycling ( n upstream 1.95 × 10 19 m 3) case for ITER PFPO-1. (Left) 2D poloidal flows visualized on top of the 2D electron density profile. The color of the quivers shows the importance of the E × B term. The importance is here defined as E × B factor | v E × B | | v pol v E × B |; the magnitude of the E × B velocity divided by the magnitude of poloidal velocities excluding the E × B velocity. The thin white line is the separatrix. (Right) Electric potential (white) contour on top of the electron density profile. E × B drifts are along the (white) potential contour lines. The red contours are the electron temperature Te for 1, 10, and 50 eV.

Close modal

The quasi-steady-state simulation shown in Fig. 11 is in a high recycling regime ( n upstream 1.95 × 10 19 m 3). It is still attached, but clearly there is already some density transport upwards helped by E × B flows. Figure 11 (right) shows a contour plot of the electric potential; and the E × B flow direction is along these contours. It is assumed that the wall structures are grounded and we can treat the electric potential with a Dirichlet boundary condition. The plasma tends to form three E × B vortex regions (i.e., with a closed contour in a small region), all relatively close to the X-point. One on the LFS (low-field-side), one on the HFS (high-field-side), and one in the private region just below the X-point toward the HFS. When an E × B vortex becomes dominant over non-E × B flows, it can form a high-density region, in the attached plasma case, these vortices are not dominant.

The E × B flows are dominant around the X-point, both up- and downstream. The E × B drift causes the plasma to favor distribution of the upstream plasma toward the OT. The plasma transport from the OT to the IT via the private region is strongly enhanced by E × B flows. At the IT strike-point, the plasma is partially transported upwards. The low temperature allows more neutrals to travel further upwards before being ionized; allowing for more density above the separatrix strike point, leading to the broadening of the ion flux. The E × B flows then transport the plasma further along the target upwards. With increased puffing the plasma transports more and more to the HFS E × B vortex, leading to a transition in the plasma (see Sec. VIII).

From Fig. 11 we can also observe how the core plasma is being fueled from the divertor plasma. In Fig. 11, along the separatrix the plasma is flowing downwards. In the SOL, both on the IT and OT, the plasma moves upwards. However, the flow is not exactly aligned with the flux surfaces; they are slightly pointing inward. The flow inversion along the separatrix occurs around the OMP for this scenario. The flow inversion point varies with different amounts of gas puffing. It typically moves upwards with more puffing.

In Secs. V and VI, steady state divertor solutions have been presented. The JOREK steady state solutions are in fact the result of time dependent simulations, with an accurate description of the dynamic evolution of the system. Depending on the timescale of the change in parameters (such as the gas puffing rate), the time dependent solution will not be the same as a series of steady state solutions. In the following, we will discuss some aspects of the time dependent behavior, including the dynamics of detachment and of the formation of the HFSHD region.

JOREK simulations (with kinetic neutrals) are typically started from a initial scenario without neutrals and the divertor legs being fully developed with a attached divertor. As we proceed with puffing neutrals, a dynamic plasma response is obtained from an attached to conduction-limited, to high-recycling and eventually detachment, depending on the puffing rate.

As an example, Fig. 12 shows the time evolution of the plasma heat load to the wall, plotted as function on n upstream ( t ) at the OMP. As soon as the neutral gas is introduced (at t = 0 ms) in the plasma, we observe a quick decay of 8.3 MW of power in the first 5 ms. Simultaneously, the upstream density increases from 1 to 1.5 × 10 19 m−3. In the following 36.6 ms, the power further decreases by only another 1.3 MW, whereas the upstream density slowly increases, reaching saturation after 30–40 ms. Thus, for this ITER scenario, at least 40 ms (in real time) has to be simulated to obtain a steady-state solution. The figure also includes the power to the divertor from the steady state solutions (red dots). It clearly illustrates the difference between the dynamic and steady state solutions: in the dynamic solution, the power to the divertor decreases about twice as fast (as a function of upstream density) compared to the steady state solution.

FIG. 12.

Time-dependent total plasma heat load to the divertor as function of the upstream density at the outer mid-plane. The vertical dashed lines indicate the tens of ms have passed after the gas puff reached the plasma in the divertor. The points shown are each 100 JOREK fluid time steps ( Δ t = 0.2246 ms) apart. The red dots indicate the steady-state results as shown in Fig. 7.

FIG. 12.

Time-dependent total plasma heat load to the divertor as function of the upstream density at the outer mid-plane. The vertical dashed lines indicate the tens of ms have passed after the gas puff reached the plasma in the divertor. The points shown are each 100 JOREK fluid time steps ( Δ t = 0.2246 ms) apart. The red dots indicate the steady-state results as shown in Fig. 7.

Close modal

The initially rather quick response of the heat load is expected as two simultaneous effects work on decreasing the heat load. First, the plasma-neutral interaction leads to energy losses in the divertor by ionization and radiation. Second, due to the fueling, the plasma becomes more high-recycling, leading to an strong increase in the plasma density at the target. The plasma “tries” to equilibrate the pressure over a magnetic field line, and thus, as the density increases, the plasma temperature decreases.

The time and manner in which this initial power decay occurs is very similar for simulations with different rates of puffing neutrals in the divertor. To which degree the power decreases strongly depends on this rate. From our simulations, typically 50 % of the change in nupstream occurs in the 5 ms after the neutrals reach the plasma. This behavior might change depending on the waveform used for puffing neutrals and is associated with the propagation of the increased density source upstream toward the mid-plane separatrix whose timescale is dictated by the ion transport timescale along the field line from the divertor to the mid-plane.

If we remove the evolution of the electric potential and vorticity, (i.e., no E × B drifts), this general pattern remains, although the path Power ( n upstream ) is less concave. Thus, the power decreases slower with an increase in nupstream than with E × B drifts included.

As the dynamic progression of the power as a function of nupstream and time does not follow the same path as the steady-state path, one might expect the heat and ion fluxes to differ as well. Figure 13 shows the time evolution of the plasma heat flux profile within one simulation along the inner and outer divertor target, with (upper row) and without E × B drifts (lower row), for a range of outer mid-plane plasma densities.

FIG. 13.

Heat flux profile evolution on the inner and outer target over time as function of the distance to the separatrix. The coloring indicates the corresponding upstream outer mid-plane plasma density. Every time slice is a 4000 fluid time steps (0.8984 ms) apart. (Top left) Inner divertor target with drifts. (Top right) Outer divertor target with drifts. (Bottom left) Inner divertor target with drifts. (Bottom right) Outer divertor target with drifts.

FIG. 13.

Heat flux profile evolution on the inner and outer target over time as function of the distance to the separatrix. The coloring indicates the corresponding upstream outer mid-plane plasma density. Every time slice is a 4000 fluid time steps (0.8984 ms) apart. (Top left) Inner divertor target with drifts. (Top right) Outer divertor target with drifts. (Bottom left) Inner divertor target with drifts. (Bottom right) Outer divertor target with drifts.

Close modal

The results with drifts from Fig. 13 are from the same simulation as Fig. 12. Without drifts, the dynamic peak heat flux decreases faster as function of nupstream than seen in the steady-state results (i.e., for the same nupstream, the dynamic peak heat flux is lower). The dynamic heat flux profiles also do follow the same shapes as seen in the steady-state results. The height, width, and location of the heat flux peaks are slightly altered.

With drifts, the changes are slightly bigger. The same general behavior where the q p k ( t , n upstream ) decreases faster as function of nupstream. In addition the shape of the profile—e.g., singly or doubly peaked profile—is qualitatively different. In the dynamic profiles, much stronger secondary peaks are observed.

In the simulations with and without drifts, the heat flux is decreasing, but without drifts, we do not really see the plasma heat flux or temperature “detaching” from the IT. With drifts, it quickly almost reduces to zero. One could call this a partially detached state, because although the plasma temperature on the IT is below 5 eV, the ion flux just started to decrease.

Figure 14 shows the ion flux profile evolution over the inner and outer divertor target, with and without drifts, for a range of outer mid-plane plasma densities, nupstream. Again, this is only the plasma ion flux. The upper row shows JOREK results (with drifts) and the lower row in Fig. 13 show heat flux results of JOREK without E × B drifts.

FIG. 14.

Ion flux profile evolution on the inner and outer target over time as function of the distance to the separatrix. The coloring indicates the corresponding upstream outer mid-plane plasma density. Every time slice is a 4000 fluid time steps (0.8984 ms) apart. This corresponds with the time slices of Fig. 13 (Top left) Inner divertor target with drifts. (Top right) Outer divertor target with drifts. (Bottom left) Inner divertor target with drifts. (Bottom right) Outer divertor target with drifts.

FIG. 14.

Ion flux profile evolution on the inner and outer target over time as function of the distance to the separatrix. The coloring indicates the corresponding upstream outer mid-plane plasma density. Every time slice is a 4000 fluid time steps (0.8984 ms) apart. This corresponds with the time slices of Fig. 13 (Top left) Inner divertor target with drifts. (Top right) Outer divertor target with drifts. (Bottom left) Inner divertor target with drifts. (Bottom right) Outer divertor target with drifts.

Close modal

The results with drifts from Fig. 14 are from the same simulation as Figs. 12 and 13. At low nupstream, with and without drifts, the profiles are similar for both the IT (inner target) and OT (outer target). When increasing the neutral fueling; and thus, also nupstream, a strong increase occurs into the ion flux as the plasma starts developing to a high-recycling regime.

In the range below n upstream 2 × 10 19, for equal nupstream, the steady-state peak ion flux is lower than the dynamic peak ion flux. The dynamic ion flux reaches earlier higher recycling stage than would be expected from steady-state solutions due to the dynamic SOL transport, with parallel transport being limited by the parallel sound speed, discussed above. Furthermore, the dynamic ion flux profiles are qualitatively similar to the steady-state solutions.

With enough gas puffing and a high enough upstream density, the losses in the divertor region become large enough such that the plasma will detach from the wall. In literature, many definitions have been used for describing and declaring the state of plasma detachment (e.g., Ref. 54) In this work, we shall not take a very strict definition of the detachment regime: the plasma is considered detached when the target electron temperature T e , target 5 eV. In this section, we explore mostly more qualitative effects of detaching the plasma.

In Sec. VI, we have already presented some results of plasmas in detached states (Figs. 7–10). Including drifts, Fig. 7 shows a transition in the power as function of nupstream. Figure 9 shows the shift of the ion flux around nupstream indicating a movement of the divertor plasma density in the divertor and at higher nupstream a decay of the ion flux. Figure 10 shows the almost complete decay of the heat flux with high nupstream. In all these figures, we observe a strong dependence on E × B drifts, which eases the detachment on the IT. This section will focus in depth on the JOREK simulations of detached plasmas; and how E × B drifts cause the high-field-side high-density region and a phase-like transition in the plasma when detaching.

During a simulation, the heat load to the wall decreases as expected with increasing nupstream. This continues until we reach a critical nupstream, as shown in Fig. 15. The power toward the wall continues decreasing, however nupstream suddenly starts decreasing. nupstream drops slightly and then increases again with almost a constant plasma heat load to the wall. In this particular simulation, the transition lasts for about 6 ms. The stagnation and drop of nupstream indicate a change in the fueling behavior from the divertor. The plasma flows from the divertor going upstream through the sol must have been changed. This then settles, and again start allowing for plasma fueling. By this point, the power on the IT has already been diminished, so this transition is mostly an effect on the OT power, as shown in Fig. 7.

FIG. 15.

Time-dependent total plasma heat load to the divertor walls as function of the upstream density at the outer mid-plane. Progression toward detachment, at around n upstream 2.4 × 10 19 m 3 the plasma undergoes a transition. The power keeps slowly decreasing while nupstream suddenly decreases, and then starts increasing again. The vertical dashed lines indicate the amount of ms that have passed after the gas puff reached the plasma in the divertor. The points shown are each 100 JOREK fluid time steps ( Δ t = 0.2246 ms) apart.

FIG. 15.

Time-dependent total plasma heat load to the divertor walls as function of the upstream density at the outer mid-plane. Progression toward detachment, at around n upstream 2.4 × 10 19 m 3 the plasma undergoes a transition. The power keeps slowly decreasing while nupstream suddenly decreases, and then starts increasing again. The vertical dashed lines indicate the amount of ms that have passed after the gas puff reached the plasma in the divertor. The points shown are each 100 JOREK fluid time steps ( Δ t = 0.2246 ms) apart.

Close modal

Note, that the stagnation and drop of nupstream in this transition is not seen in the steady-state results from Fig. 7. Again, as also shown in Fig. 12, the evolution of the power to the targets as function of nupstream in these dynamic simulations does not follow the same path as the steady-state results. This is due to the short timescale in which the gas puffing is scaled up, and equilibration time between up- and downstream. The neutral gas fueling rate increases more quickly in Fig. 15 than in Fig. 12, this is the reason the power to the divertor is lower for the same nupstream.

This transition (i.e., the stagnation, drop, and re-increase in nupstream) occurs over many JOREK time steps and is well resolved in both space and time. If we study 2D results of the simulations, we observe the nupstream back-and-forth transition coincides with the formation and movement of a high-field-side high-density region.

Figure 16 shows a comparison of the electron density in the poloidal plane with drifts (left) and without drifts (right) for comparable nupstream. Both simulations are just detached (i.e., T e , target 5 eV). Including the E × B drifts lead to qualitatively very different results. E × B drifts are essential for transporting plasma across the field and redistributing the plasma flow from the core and divertor. In Fig. 16 (right), without E × B drifts, the density is continuously building up; however, there is no mechanism of transporting it away from the divertor legs. This behavior is quite symmetrical and agrees qualitatively well with SOLPS-ITER without drifts.18 As explained in Sec. VI, the E × B drifts push the OT plasma leg downwards and transport plasma via the private region from the OT to the IT. On the IT, the E × B drifts push the plasma upwards enhancing the cross field transport. With more fueling from the divertor, the HFS E × B vortex becomes dominant and starts trapping increasingly more density forming the so-called HFSHD-region.

FIG. 16.

2D poloidal plane of n e ( 10 20 m 3 ) for JOREK simulation with (left) and without (right) drifts. The central electron density is 1.2 × 10 19 m 3 Both runs have identical fueling amounts but stabilize at different upstream electron densities.

FIG. 16.

2D poloidal plane of n e ( 10 20 m 3 ) for JOREK simulation with (left) and without (right) drifts. The central electron density is 1.2 × 10 19 m 3 Both runs have identical fueling amounts but stabilize at different upstream electron densities.

Close modal

The HFSHD-region is not a passive structure, it influences the flows in the divertor region while transitioning to a detached plasma state. This can be observed more clearly by examining the dynamic ion flux profiles toward the wall. We note that SOLPS simulations of ITER divertor plasmas with drifts at high SOL power and with impurities show a similar behavior with the appearance of an HFSHD-region because of similar flow patterns.29,31 In our simulations, however, we find that this phenomenology is driven by drifts and main plasma flow dynamics not specifically requiring high power operation and the injection of impurities, as described below.

Figure 17 shows the ion flux profiles as function of s s sep progressing over time for comparable simulations with and without drifts. The color indicates nupstream. This dynamic response shows how the plasma transition toward detachment. Figure 17 (bottom row) shows the ion flux profiles without drifts. The ion flux first increases with increasing nupstream. The peak at the IT occurs at a slightly lower nupstream than the OT peak. The progression of the ion flux over time and nupstream is qualitatively the same for the IT and OT. Both sides develop double peaked profiles. One could argue the plasma is not fully detached as the ion flux has not fully collapsed; however, for the goal of comparing the physics behavior, this is sufficient.

FIG. 17.

Ion flux profile evolution on the inner and outer target over time as a function of the distance to the separatrix. The ion flux profile evolution shows a dynamic ion flux rollover into detachment for JOREK simulation with and without drifts. The colors indicate the upstream density. Each successive line is separated by 100 JOREK fluid time steps (i.e., 0.2246 ms real time).

FIG. 17.

Ion flux profile evolution on the inner and outer target over time as a function of the distance to the separatrix. The ion flux profile evolution shows a dynamic ion flux rollover into detachment for JOREK simulation with and without drifts. The colors indicate the upstream density. Each successive line is separated by 100 JOREK fluid time steps (i.e., 0.2246 ms real time).

Close modal

Figure 17 (top row) shows the ion flux profiles with E × B drifts. The dynamic response and qualitative results differ strongly from the case without drifts. E × B drifts favor distribution of ion flux toward the OT creating higher peak ion fluxes. The ion flux rollover point on the OT is around the same nupstream for both simulations with and without drifts. The ion flux on the IT with drifts [Fig. 17 (top left)] has lower peak ion fluxes than its no-drift counterpart. With increasing nupstream, the IT ion flux diminishes more quickly and more upwards (to higher s s sep). At around n upstream 2.4 × 10 19 m−3, the IT ion flux sharply increases around 9 cm above the separatrix strike point. At even higher nupstream, the IT ion flux decays rather quickly again. During this transition in the plasma, the HFSHD-region starts forming and moving and this has a strong impact on the resulting ion flux at the IT.

The IT ion flux disappears completely at the high end of nupstream, thus reaching full detachment, while this is not the case for the ion flux without drifts. In this 20 MW heating scenario, full detachment at the IT always goes paired with the formation of the HFSHD-region.

For the ITER PFPO-1 H-mode 5MA 20 MW scenario, the JOREK simulations with E × B drifts always develop a HFSHD-region when going from the conduction-limited regime toward detachment. To fully understand the detachment behavior, we now look more closely on the HFSHD-region, its formation and its interaction with divertor and core plasma.

1. Importance of the HFS E × B vortex

Figure 11 has shown that in an attached high-recycling scenario three E × B vortices form; one on the HFS, one the LFS and one in the private region, all in the vicinity of the X-point. In the lower range of nupstream, the E × B drifts due to the vortices are weaker than the other plasma flows associated with a vertical divertor configuration and the associated recycling. With increasing nupstream, especially the HFS E × B vortex builds up in combination with the increased density coming from the IT strike point and this strongly impacts the divertor plasma flow.

Figure 18 shows the detached scenario with a fully formed HFSHD-region. It shows the electron density profile with the poloidal flow vectors colored with the E × Bfactor (the ratio of the E × B flow amplitude divided by the poloidal component of the parallel flow | v E × B | | v , pol |) on Fig. 18 (left) and the electric potential contours on Fig. 18 (right). Specifically, the HFS E × B vortex has grown to become dominant as the E × B factor 1. The HFS E × B vortex draws in plasma density from the IT, and this becomes trapped in the vortex. Once the HFSHD-region is stabilized, the electric potential forms a plateau where the density is highest. We observe the potential contour patterns from Fig. 11 become more exaggerated toward detachment. The LFS E × B vortex is lower and does not create a vortex flow pattern. The private region E × B vortex contributes to enhanced transport from the X-point toward the OT and decreases transport from the X-point to the IT. The private region E × B vortex is the dominant actor in transporting the plasma from the OT to the IT via the private region.

FIG. 18.

2D electron density profiles for a detached case for ITER PFPO-1. (Left) 2D poloidal flows visualized on top of the 2D electron density profile. The color of the quivers shows the importance of the E × B term. The importance is here defined as E × B factor | v E × B | | v pol v E × B |; the magnitude of the E × B velocity divided by the magnitude of poloidal velocities excluding the E × B velocity. The quivers are visualized every 50th gauss point of the grid to not saturate the image. The thin white line is the separatrix. (Right) Electric potential (white) contour on top of the electron density profile. E × B drifts are along the potential contour lines. The red contours are the electron temperature Te for 1, 10, and 50 eV.

FIG. 18.

2D electron density profiles for a detached case for ITER PFPO-1. (Left) 2D poloidal flows visualized on top of the 2D electron density profile. The color of the quivers shows the importance of the E × B term. The importance is here defined as E × B factor | v E × B | | v pol v E × B |; the magnitude of the E × B velocity divided by the magnitude of poloidal velocities excluding the E × B velocity. The quivers are visualized every 50th gauss point of the grid to not saturate the image. The thin white line is the separatrix. (Right) Electric potential (white) contour on top of the electron density profile. E × B drifts are along the potential contour lines. The red contours are the electron temperature Te for 1, 10, and 50 eV.

Close modal

The transition toward a dominant HFS E × B vortex is a rather dynamic phenomenon. This transition is shown as a series of five 2D ne profiles in Fig. 19. The arrows here indicate a general movement of the density. There is an initial slow buildup phase of density close to the wall. This phase is followed by a density front building up, moving from the wall toward the X-point. The E × B vortex is significant enough to move this high-density front. Once the front is close to the X-point, it quickly builds up more and more density (till n e 10 × 10 20 m−3). It then reaches a critical point when the high density region close the X-point suddenly moves away from the X-point and installs itself at the stable location shown in the last profile in Fig. 19 and also in Fig. 18. Specifically, during this fast movement away from the X-point the flow patterns everywhere are altered. This is most notable on the HFS of the divertor but can also be seen on the OT and in the entire SOL.

FIG. 19.

The formation of the high-field-side high-density region presented as a series of the 2D electron density profiles in the poloidal plane. The arrows indicate general movement of the high density. Note that due to the E × B vortex, the precursor and HFSHD-region rotate as well. Thus, it does not solely moves in- and outwards.

FIG. 19.

The formation of the high-field-side high-density region presented as a series of the 2D electron density profiles in the poloidal plane. The arrows indicate general movement of the high density. Note that due to the E × B vortex, the precursor and HFSHD-region rotate as well. Thus, it does not solely moves in- and outwards.

Close modal

The buildup and movement of the HFSHD-region naturally influences its surroundings. This is why we see this dynamic phenomenon also in the ion flux profiles (Fig. 17) and even in the upstream density (Fig. 15). Due to the strong potential difference the HFSHD-region pushes and pull the plasma flow in the divertor legs.

The exact equilibrium position on which the HFSHD-region settles is somewhat dependent on the gas fueling rate in the divertor. Altering the gas fueling rate results in horizontal motion. Exactly how this relation works is not yet fully understood in our simulations.

These simulations only consider as source atomic neutrals, and no molecular species and/or impurities. Since molecules can change the momentum loss along the SOL and, together with impurities, increase energy losses, the formation and position of the HFSHD-region could change. However, the E × B vortex regions always form with plasma density build-up in the divertor and are, thus, almost unavoidable. Second, this simulation does not contain the diamagnetic drift (i.e., p drift). The transport in the SOL will be altered with the addition of diamagnetic drifts. However, as the HFSHD-region forms during detachment in regions of low pressure, the diamagnetic drift is not likely to alter the transport significantly. The HFSHD-region itself, although containing very high densities, is a very low pressure region.

2. Fueling from the divertor

The poloidal velocity profile changes strongly in the presence of the HFSHD-region. There is a clear qualitative difference in the v pol profiles between Figs. 11 and 18. This means, that with enough divertor gas puff, the flows in the divertor, and from the divertor to the core plasma will become significantly different. The ion mass flux in the poloidal plane Γ pol = ρ v pol, then tells us in this scenario how the core is fueled from the divertor. This subsection highlights the changes in Γ pol after the formation of the HFSHD-region.

Figure 20 shows two schematics of the ion flux directions. Figure 20 (left) represents the attached plasma case, and Fig. 20 (right) the detached plasma after the formation of the HFSHD-region. Although the ion flux profiles are more complicated, this representation illustrates the general behavior. The color of the arrows indicates from where the flows originate, the size of the arrows do not indicate the magnitude; thus, they only trace the flow direction. The flux surface chosen are the separatrix and an approximate flux surface to indicate separation between the SOL and far-SOL regions.

FIG. 20.

Schematic illustration of the ion flux direction in the poloidal plane. Orange indicates flows from the core plasma. Cyan indicates flows from the divertor plasma close to the wall. Purple indicates the HFSHD-region or HFS E × B vortex. Green indicates some flows that arise due to multiple flow inversions. (Left) The poloidal mass flows in the attached case corresponding to Fig. 11. (Right) The poloidal mass flows in the detached case corresponding to Fig. 18.

FIG. 20.

Schematic illustration of the ion flux direction in the poloidal plane. Orange indicates flows from the core plasma. Cyan indicates flows from the divertor plasma close to the wall. Purple indicates the HFSHD-region or HFS E × B vortex. Green indicates some flows that arise due to multiple flow inversions. (Left) The poloidal mass flows in the attached case corresponding to Fig. 11. (Right) The poloidal mass flows in the detached case corresponding to Fig. 18.

Close modal

Along the separatrix the ion flow from the core flows through (and slightly around) the X-point toward the strike point. Figure 20 (left) shows that due to E × B drift the ion flux is mostly distributed to the outer divertor leg. In the private flux region, a small part is redirected to the inner divertor leg. Another part of the private flux region transport comes from recycling on the outer target. The blue arrows indicate flows from the divertor plasma. These arise from plasma recycling and ionization of the gas puff. Figure 20 (left) shows the flux in the SOL comes mostly from the divertor plasma. Higher up along the SOL the flows are almost parallel to the separatrix, but it has a small component crossing the flux surfaces ( Γ pol · ψ 0). Depending on the exact gas puffing, in the attached case, the flow inversion point on the separatrix is somewhere between the OMP and the top of the plasma.

After the HFSHD-region has formed, the Γ pol profile has qualitative strongly changed. Figure 20 (right) shows in purple the HFSHD-region and with it the HFS E × B vortex. The recycled flux (cyan) from the IT is interrupted traveling upwards and circles back. The flow from the core (orange) travels from the X-point to halfway along the divertor plasma leg at the OT. Along the leg it gets deflected into the private region toward the IT. Due to multiple flow inversions the result is a flux along the inner divertor plasma leg toward the X-point that then goes into the private region. Generally, the flow patterns become more complicated.

These changes due to the HFSHD-region are not only confined to the divertor region. Due to the interrupted flow from the IT, there is no plasma traveling upstream in the HFS SOL. The flow reversal point along the separatrix is now approximately where the blue arrow crosses the separatrix in Fig. 20 (right). That entails that the flows in the SOL come for the OT and travel all along the SOL to finally reach the X-point from the HFS and the HFSHD-region. It should be noted that strong changes of plasma flows with the onset of detachment have been long seen in divertor tokamaks.55,56

To observe whether the plasma retains this behavior on a longer timescale, some simulations have been extended to 150 ms. From the typical timescales covered by JOREK, this simulation is completely steady-state. However, as in the simulation the upstream density was ramping up due to increased fueling, it takes of the order of a confinement time before this is fully equilibrated with the core of the plasma at the magnetic axis. The flows and effective fueling will still change on the timescale of a confinement time when considering equilibration of the core plasma density or the applied fueling.

3. Ionization patterns in the presence of HFSHD-region

As a final dive investigating detachment and the HFSHD-region with our kinetic neutral model, we investigate the ionization patterns and the neutral density distribution. In the method, how the kinetic particles are coupled to JOREK, moments of the particle distribution can be projected to the same finite element grid. The ionization rate moment in Fig. 21 (left) is normalized and, for visibility reasons, both plots only show four orders of magnitude. The neutral density is presented as the number density (atoms/m3).

FIG. 21.

Results taken from the same simulation and time stamp as Fig. 18. (Left) The normalized ionization rate in the divertor on a log-scale. The thin white line is the separatrix. (Right) The neutral atomic density m−3 on a log-scale. Both log-scales span four orders of magnitude more clearly show the patterns.

FIG. 21.

Results taken from the same simulation and time stamp as Fig. 18. (Left) The normalized ionization rate in the divertor on a log-scale. The thin white line is the separatrix. (Right) The neutral atomic density m−3 on a log-scale. Both log-scales span four orders of magnitude more clearly show the patterns.

Close modal

As this simulation only fuels kinetic neutrals from the divertor region, there are only two ways by which neutrals can be found outside the divertor region. One is by recombination of the core plasma, and more importantly by transport across the divertor and SOL plasma. In the attached plasma case, nearly all neutrals are ionized and do not cross the divertor plasma leg. E × B drifts transport the plasma higher up the inner target where they undergo wall-assisted recombination. This provides the main source of neutrals outside of the private region. Figure 21 (right) show that in the detached case there is a significant neutral density close to the inner target at the height of the X-point. The ionizing region is still close to the wall on IT. Te at the IT strike point is approximately 1 eV; this is low enough such that neutrals can transport past this region. In this particular case the outer target strike point temperature is approximately 10 eV, thus, not detached. On the OT, the ionization region is strongest close the wall and, due to the even higher temperature in the OT, SOL transport do not transport neutrals significantly upwards on the OT.

It is expected that when we introduce molecular hydrogen to the kinetic neutral model, we will obtain a more realistic ionization front. Since molecular hydrogen–plasma interaction introduces many more energy loss channels, it is expected that the plasma will detach at lower nupstream. Specifically close to the strike points, molecules should lower Te for the 20 MW scenario since this is a rather low power level for these ITER plasmas. For the ITER Q = 10 scenarios with 100 MW, solely adding neutrals shall not be sufficient to detach the plasma. For studying detachment in these conditions, the puffing of impurity species is needed such as Neon.

Although the HFSHD-region is a high density region, the pressure in this region is low. Ionization occurs most significantly on the edge of the HFSHD-region, where the temperature sharply increases again. This allows the HFSHD-region as well to gather neutral density.

Note that if a significant part of the gas fueling would be provided by the top gas valve (as shown in Fig. 1), this would cooldown the (far-) SOL temperature. This could then allow for stronger neutral transport across the OT strike point, although this is not found to have a major impact on divertor behavior according to Ref. 18.

Realistic divertor conditions are important to predict impurity sputtering and transport and to study the consequences due transient MHD activity on the divertor. We have derived a plasma-atomic neutral interaction model based on the effective rate coefficients (From OpenADAS); this model has been implemented in JOREK as both a fluid neutral and a kinetic neutral model. The fluid neutral model evolves only the neutral mass density equation. The kinetic neutral model evolves individual superparticles in 3D geometry and includes ionization, recombination, and charge-exchange. The kinetic model further includes plasma and neutral boundary interactions based on the SDTRIM coefficients.

To benchmark our new models, we have performed a neutral gas puffing scan for an ITER early operations H-mode scenario (ITER PFPO-1, 5MA, 20 MW, H-mode) and compared it to SOLPS-ITER simulations (without poloidal drifts).18 To compare relatively similar physics models, we reduced the JOREK physics model by removing the evolution the magnetic fields and by removing E × B drifts. For all simulations, a divertor gas puff was used.

Using the fluid neutral model, we observe a power decay with increased gas puffing, but quantitatively, it does not match well. A one-variable fluid neutral model with spatially constant diffusivity is not sufficient for realistic divertor plasma modeling.

Using the kinetic neutral model, at low-to-medium upstream densities ( n upstream = 1 × 10 19 2 × 10 19 m 3), the total plasma heat load of the no-drift JOREK models agrees well (difference < 10 %) with SOLPS-ITER simulation.18 At higher upstream densities ( n upstream > 2 × 10 19 m 3), the transition to detachment occurs at a 10 % higher upstream density and is less sharp than the detachment transition in SOLPS-ITER. The JOREK kinetic neutral model only includes atomic species. The extra energy loss channels provided by molecular species are not included, and thus, JOREK obtains detachment condition ( T target < 5 eV) at higher upstream densities than SOLPS-ITER. The JOREK peak heat fluxes in the conduction-limited regime are approximately twice as high as SOLPS-ITER simulation. The exact underlying reason for this is not fully identified, but it is expected to be due to the use of a single temperature model that does not distinguish electrons and ions for parallel transport in the SOL.

The kinetic neutral model using 10 6 simulation particles is as fast as the neutral fluid model and scales more favorably with increasing computational nodes.

As a first application of the developed kinetic neutral model, we have looked at the effects of E × B drifts on the divertor solutions and the plasma heat load. The electric field forms vortices (on either side of the X-point and in the private flux region), leading to predominantly poloidal E × B flows. In the steady-state, E × B drifts result in a faster, compared to no-drifts, decrease in power to the wall as a function of the outer mid-plane separatrix density. This comes paired with a more asymmetric distribution of power between the inner and outer targets and allows the inner target to detach more easily with a rapidly diminishing heat flux. At higher upstream density, n upstream 2.2 × 10 19 m 3, a phase-like transition in the plasma develops. This coincides with the formation of a high-field-side high density (HFSHD)-region. This phase-like transition only develops with E × B drifts.

Understanding this transition requires time-dependent simulations since the formation of the HFSHD-region is a transient phenomenon at a timescale of 10 20 ms. While the plasma density in the divertor is ramping up, an E × B vortex close to the X-point on the HFS is becoming dominant. Simultaneously, E × B flows enhance cross field transport upwards at the IT. Plasma density is collected in the center of the HFS E × B vortex and the newly formed HFSHD-region moves horizontally toward the X-point and back, reaching a stable position. The HFS E × B vortex alters the poloidal flows not only in the divertor region, but in the SOL as a whole. This may have consequences for plasma fueling.

Next steps to further improve the divertor solutions in JOREK include the addition of hydrogenic molecules. JOREK models already exist including two temperatures, diamagnetic flows, and radiative kinetic impurities, but have not been used here. The effects of these will be left for further studies. The molecules increase the energy and momentum losses in the divertor changing the precise upstream density at which detachment occurs. Since the HFSHD-region is a low pressure region, the diamagnetic drift is not expected to have a strong influence there. On the contrary, for the H-mode edge plasmas in which large pressure gradients are expected, diamagnetic flows should play an important role.

With the implementation of a kinetic neutral model and, thus, more realistic divertor solutions, we have improved the capability to study the consequences of MHD activity (ELMs and disruptions) on the divertor plasmas and associated power fluxes including the application of 3D fields for ELM control which render the divertor plasma 3D as well.

This work was part of the research programme (No. 19KE01) of the Foundation for Nederlandse Wetenschappelijk Onderzoek Instituten (NWO-I), which was part of the NWO and has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No. 101052200 – EUROfusion) and by the ITER Organization under Implementing Agreement No. 43-2174 with the Eindhoven University of Technology (TU/e) and with TU/e work order No. 107640. ITER is the Nuclear Facility, INB No. 174. Computations were performed using computer clusters Marconi-Fusion and the ITER Scientific Data and Computing Center (SDCC). The views and opinions expressed herein do not necessarily reflect those of the European Commission nor the ITER Organization nor NWO-I. Discussions with members of the ITER Organization (R. A. Pitts, S. D. Pinches, and X. Bonnin) are acknowledged.

The authors have no conflicts to disclose.

Sven Quirijn Korving: Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (lead). Guido T. A. Huijsmans: Funding acquisition (equal); Methodology (equal); Resources (equal); Software (equal); Supervision (equal); Writing – review & editing (equal). Jae-Sun Park: Data curation (equal); Methodology (equal); Resources (equal). Alberto Loarte: Funding acquisition (equal); Methodology (equal); Resources (equal); Supervision (equal); Writing – review & editing (equal).

Raw data were generated at the ITER Scientific Data and Computing Center (SDCC) and Marconi-Fusion large scale facility. Derived data supporting the findings of this study are available from the corresponding author upon reasonable request.

This appendix contains the JOREK form of the conservation equations used in Sec. III. JOREK uses a different form of the conservation equations than the Eularian form. As we rewrite the terms the JOREK form of the equations are still conservative. The assumptions used in JOREK are not included in this sections. The projections of the momentum equations to obtain the reduced MHD model are not included either.

1. Mass conservation
Instead of the continuity equation for number density, JOREK uses mass density. We obtain this form simply by multiplying Eq. (13) by the mass of the particles. We then obtain
( m i n i t ) i , r = m i ( Γ i ion Γ n rec ) , ( m e n e t ) e , r = m e ( Γ e ion Γ n rec ) , ( m n n n t ) i , r = m n ( Γ n ion + Γ n rec ) .
(A1)
2. Momentum equation
To write the momentum conservation to the JOREK form for any species, we can expand the derivative with the product rule and write the equation as
m α n α ( v α t ) i , r = ( m α n α v α t ) m α v α ( n α t ) .
(A2)
The second term in the right-hand side of Eq. (A2) makes sure that the equation is still in a conservative form; and that the code will still have good conservation properties. Writing Eq. (14) in the form of Eq. (A2) gives us the following result:
m i n i ( v i t ) i , r = Γ n rec m i v i m i v i ( Γ n rec ) + Γ i ion m i v n m i v i ( Γ i ion ) = Γ i ion m i ( v n v i ) , m e n e ( v e t ) i , r = Γ n rec m e v e m e v e ( Γ n rec ) + Γ i ion m e v n m e v e ( Γ i ion ) = Γ i ion m e ( v n v e ) , m n n n ( v n t ) i , r = Γ n rec ( m i v i + m e v e ) m n v n ( Γ n rec ) Γ i ion m n v n m n v n ( Γ n ion ) = Γ n rec ( m i v i + m e v e m n v n ) .
(A3)
Note that the recombination terms for the ions and electrons disappear, and the ionization term disappears for the neutrals.
3. Energy equation
To obtain the JOREK form of the energy conservation, we can expand the kinetic energy part of the derivative and rewrite the equation as follows:
( E α t ) = ( ( 1 2 m α n α v α 2 + n α T α γ 1 ) t ) , ( ( n α T α ) t ) = ( γ 1 ) [ ( E α t ) v α ( m α n α v α t ) 1 2 m α v α 2 ( n α t ) ] ,
(A4)
which gives us terms from the continuity and momentum equation in our energy/pressure equation. It is important to note that these terms are also multiplied by ( γ 1 ) as we are now solving for a pressure equation. In the next sections, the ionization and recombination terms are presented.
a. Ionization
Writing Eq. (16) to the JOREK form and rewriting some vector terms gives us
( ( n i T i ) t ) ion = Γ i ion m i m n T n + ( γ 1 ) Γ i ion 1 2 m i ( v n v i ) 2 , ( ( n e T e ) t ) ion = Γ i ion m e m n T n + ( γ 1 ) Γ i ion ( 1 2 m e ( v n v e ) 2 ϕ ion ) , ( ( n n T n ) t ) ion = Γ n ion T n ,
(A5)
where ( v n v i ) 2 comes from v n 2 2 v i · v n + v u 2. There is further no physical meaning to this term as it purely arises due to the JOREK form of the equation.
b. Recombination
Writing Eq. (17) to the JOREK form gives us the following Eq. (A6), which simplifies the ion and electron terms. In the neutral equation, the quadratic velocity terms in the neutral pressure equation could not be reduced,
( ( n i T i ) t ) rec = Γ n rec T i , ( ( n e T e ) t ) rec = Γ n rec T e , ( ( n n T n ) t ) rec = Γ n rec ( m i m n T i + m e m n T e ) + ( γ 1 ) Γ n rec ( 1 2 m n v n 2 + 1 2 m i v i 2 + 1 2 m e v e 2 m i v n · v i m e v n · v e ) .
(A6)
1.
M.
Hoelzl
,
G. T. A.
Huijsmans
,
S. J. P.
Pamela
,
M.
Bécoulet
,
E.
Nardon
,
F. J.
Artola
,
B.
Nkonga
,
C. V.
Atanasiu
,
V.
Bandaru
,
A.
Bhole
et al, “
The JOREK non-linear extended MHD code and applications to large-scale instabilities and their control in magnetically confined fusion plasma
,”
Nucl. Fusion
61
,
065001
(
2021
).
2.
I.
Vasileska
,
X.
Bonnin
, and
L.
Kos
, “
Kinetic-fluid coupling simulations of ITER Type I ELM
,”
Fusion Eng. Des.
168
,
112407
(
2021
).
3.
H.
Du
,
C.
Sang
,
L.
Wang
, and
X.
Bonnin
, “
Numerical simulation of the energy deposition evolution on divertor target during type-III ELMy H-mode in EAST using SOLPS
,”
Fusion Eng. Des.
61
,
2461
2466
(
2014
).
4.
M.
Bécoulet
,
F.
Orain
,
G. T. A.
Huijsmans
,
S.
Pamela
,
P.
Cahyna
,
M.
Hoelzl
,
X.
Garbet
,
E.
Franck
,
E.
Sonnendrücker
et al, “
Mechanism of edge localized mode mitigation by resonant magnetic perturbations
,”
Phys. Rev. Lett.
113
,
115001
(
2014
).
5.
M.
Becoulet
,
G. T. A.
Huijsmans
,
C.
Passeron
,
Y. Q.
Liu
,
T. E.
Evans
,
L. L.
Lao
,
L.
Li
,
A.
Loarte
,
S. D.
Pinches
,
A.
Polevoi
,
M.
Hosokawa
,
S. K.
Kim
,
S. J. P.
Pamela
,
S.
Futatani
, and
JOREK Team
, “
Non-linear MHD modelling of edge localized modes suppression by resonant magnetic perturbations in ITER
,”
Nucl. Fusion
62
,
066022
(
2022
).
6.
G. T. A.
Huijsmans
,
S.
Pamela
,
E.
van der Plas
, and
P.
Ramet
, “
Non-linear MHD simulations of edge localized modes (ELMs)
,”
Plasma Phys. Controlled Fusion
51
,
124012
(
2009
).
7.
A.
Loarte
,
G.
Huijsmans
,
S.
Futatani
,
L. R.
Baylor
,
T. E.
Evans
,
D. M.
Orlov
,
O.
Schmitz
,
M.
Becoulet
,
P.
Cahyna
,
Y.
Gribov
et al, “
Progress on the application of ELM control schemes to ITER scenarios from the non-active phase to DT operation
,”
Nucl. Fusion
54
,
033007
(
2014
).
8.
F. J.
Artola
,
A.
Loarte
,
M.
Hoelzl
,
M.
Lehnen
,
N.
Schwarz
, and
JOREK Team
, “
Non-axisymmetric MHD simulations of the current quench phase of ITER mitigated disruptions
,”
Nucl. Fusion
62
,
056023
(
2022
).
9.
G. T. A.
Huijsmans
and
O.
Czarny
, “
MHD stability in X-point geometry: Simulation of ELMs
,”
Nucl. Fusion
47
,
659
(
2007
).
10.
O.
Czarny
and
G. T. A.
Huijsmans
, “
Bézier surfaces and finite elements for MHD simulations
,”
J. Comput. Phys.
227
,
7423
7445
(
2008
).
11.
M.
Hoelzl
,
A.
Cathey
,
S.
Futatani
,
G. T. A.
Huijsmans
,
F.
Orain
,
M.
Dunne
,
S.
Pamela
,
M.
Becoulet
,
F. J.
Artola
,
D. C.
van Vugt
,
S.
Smith
,
N.
Schwarz
,
F.
Liu
,
S. Q.
Korving
et al, “
Simulations of edge localized mode (ELM) cycles and ELM control
,” in
28th IAEA Fusion Energy Conference (FEC)
(
2020
).
12.
A.
Cathey
,
M.
Hoelzl
,
K.
Lackner
,
G. T. A.
Huijsmans
,
M. G.
Dunne
,
E.
Wolfrum
,
S. J. P.
Pamela
,
F.
Orain
,
S.
Günter
,
JOREK Team, ASDEX Upgrade Team
, and
EUROfusion MST1 Team
, “
Non-linear extended MHD simulations of type-I edge localised mode cycles in ASDEX upgrade and their underlying triggering mechanism
,”
Nucl. Fusion
60
,
124007
(
2020
).
13.
S.
Futatani
,
A.
Cathey
,
M.
Hoelzl
,
P. T.
Lang
,
G. T. A.
Huijsmans
,
M.
Dunne
,
JOREK Team, ASDEX Upgrade Team
, and
EUROfusion MST1 Team
, “
Transition from no-ELM response to pellet ELM triggering during pedestal build-up-insights from extended MHD simulations
,”
Nucl. Fusion
61
,
046043
(
2021
).
14.
D. C.
van Vugt
, “Nonlinear coupled MHD-kinetic particle simulations of heavy impurities in tokamak plasmas,” Ph.D. thesis (
Eindhoven University of Technology
,
2019
).
15.
D. C.
van Vugt
,
G. T. A.
Huijsmans
,
M.
Hoelzl
,
A.
Loarte
,
ASDEX Upgrade Team
, and
Eurofusion MST1 Team
, “
Kinetic modeling of ELM-induced tungsten transport in a tokamak plasma
,”
Phys. Plasmas
26
,
042508
(
2019
).
16.
A.
Dvornova
, “
Hybrid fluid-kinetic MHD simulations of the excitation of toroidal Alfvén eigenmodes by fast particles and external antenna
,” Ph.D. thesis (
Eindhoven University of Technology
,
2021
).
17.
A.
Fil
,
E.
Nardon
,
M.
Hoelzl
,
G. T. A.
Huijsmans
,
F.
Orain
,
M.
Becoulet
,
P.
Beyer
,
G.
Dif-Pradalier
,
R.
Guirlet
,
H. R.
Koslowski
,
M.
Lehnen
,
J.
Morales
et al, “
Three-dimensional non-linear magnetohydrodynamic modeling of massive gas injection triggered disruptions in JET
,”
Phys. Plasmas
22
,
062509
(
2015
).
18.
J.-S.
Park
,
X.
Bonnin
, and
R. A.
Pitts
, “
Assessment of ITER divertor performance during early operation phases
,”
Nucl. Fusion
61
,
016021
(
2021
).
19.
J.-S.
Park
,
X.
Bonnin
, and
R.
Pitts
, “
Assessment of ITER W divertor performance during early operation phases
,” in
APS Division of Plasma Physics Meeting
(
2019
).
20.
I. H.
Hutchinson
,
B.
LaBombard
,
J. A.
Goetz
,
B.
Lipschultz
,
G. M.
McCracken
,
J. A.
Snipes
, and
J. L.
Terry
, “
The effects of field reversal on the Alcator C-Mod divertor
,”
Plasma Phys. Controlled Fusion
37
,
1389
(
1995
).
21.
P. C.
Stangeby
and
A. V.
Chankin
, “
Simple models for the radial and poloidal E × B drifts in the scrape-off layer of a divertor tokamak
,”
Nucl. Fusion
36
,
839
(
1996
).
22.
A. V.
Chankin
, “
Classical drifts in the tokamak SOL and divertor: Models and experiment
,”
J. Nucl. Mater.
241–243
,
199
213
(
1997
).
23.
S.
Wiesen
,
D.
Reiter
,
V.
Kotov
,
M.
Baelmans
,
W.
Dekeyser
,
A. S.
Kukushkin
,
S. W.
Lisgo
,
R. A.
Pitts
,
V.
Rozhansky
,
G.
Saibene
,
I.
Veselova
, and
S.
Voskoboynikov
, “
The new SOLPS-ITER code package
,”
J. Nucl. Mater.
463
,
480
484
(
2015
).
24.
H.
Bufferand
,
G.
Ciraolo
,
Y.
Marandet
,
J.
Bucalossi
,
P.
Ghendrih
,
J.
Gunn
,
N.
Mellet
,
P.
Tamain
,
R.
Leybros
,
N.
Fedorczak
,
F.
Schwander
, and
E.
Serre
, “
Numerical modelling for divertor design of the WEST device with a focus on plasma-wall interactions
,”
Nucl. Fusion
55
,
053025
(
2015
).
25.
H.
Bufferand
,
C.
Baudoin
,
J.
Bucalossi
,
G.
Ciraolo
,
J.
Denis
,
N.
Fedorczak
,
D.
Galassi
,
P.
Ghendrih
,
R.
Leybros
,
Y.
Marandet
,
N.
Mellet
,
J.
Morales
et al, “
Implementation of drift velocities and currents in SOLEDGE2D-EIRENE
,”
Nucl. Mat. Energy
12
,
852
857
(
2017
).
26.
A. V.
Chankin
,
G.
Corrigan
,
M.
Groth
,
P. C.
Stangeby
, and
JET Contributors
, “
Influence of the E × B drift in high recycling divertors on target asymmetries
,”
Plasma Phys. Controlled Fusion
57
,
095002
(
2015
).
27.
E. T.
Meier
,
R. J.
Goldston
,
E. G.
Kaveeva
,
M. A.
Makowski
,
S.
Mordijck
,
V. A.
Rozhansky
,
I.
Senichenkov
, and
S. P.
Voskoboynikov
, “
Comparing N versus Ne as divertor radiators in ASDEX-upgrade and ITER
,”
Plasma Phys. Controlled Fusion
58
,
125012
(
2016
).
28.
E.
Sytova
,
R. A.
Pitts
,
E.
Kaveeva
,
X.
Bonnin
,
D.
Coster
,
V.
Rozhansky
,
I.
Senichenkov
,
I.
Veselova
,
S.
Voskoboynikov
, and
F.
Reimold
, “
Comparing N versus Ne as divertor radiators in ASDEX-upgrade and ITER
,”
Nucl. Mater. Energy
19
,
72
78
(
2019
).
29.
E.
Kaveeva
,
V.
Rozhansky
,
I.
Senichenkov
,
E.
Sytova
,
I.
Veselova
,
S.
Voskoboynikov
,
X.
Bonnin
,
R. A.
Pitts
,
A. S.
Kukushkin
,
S.
Wiesen
, and
D.
Coster
, “
SOLPS-ITER modelling of ITER edge plasma with drifts and currents
,”
Nucl. Fusion
60
,
046019
(
2020
).
30.
I.
Veselova
,
E.
Kaveeva
,
V.
Rozhansky
,
I.
Senichenkov
,
A.
Poletaeva
,
R.
Pitts
, and
X.
Bonnin
, “
SOLPS-ITER drift modeling of ITER burning plasmas with narrow near-SOL heat flux channels
,”
Nucl. Mat. Energy
26
,
100870
(
2021
).
31.
V.
Rozhansky
,
E.
Kaveeva
,
I.
Senichenkov
,
I.
Veselova
,
S.
Voskoboynikov
,
R. A.
Pitts
,
D.
Coster
,
C.
Giroud
, and
S.
Wiesen
, “
Multi-machine SOLPS-ITER comparison of impurity seeded H-mode radiative divertor regimes with metal walls
,”
Nucl. Fusion
61
,
126073
(
2021
).
32.
C. K.
Tsui
,
J. A.
Boedo
,
O.
Février
,
H.
Reimerdes
,
C.
Colandrea
,
S.
Gorno
, and
TCV Team
, “
Relevance of E × B drifts for particle and heat transport in divertors
,”
Plasma Phys. Controlled Fusion
64
,
065008
(
2022
).
33.
H.
Du
,
G.
Zheng
,
H.
Guo
,
A. E.
Jaervinen
,
X.
Duan
,
X.
Bonnin
,
D.
Eldon
, and
D.
Wang
, “
SOLPS analysis of the necessary conditions for detachment cliff
,”
Nucl. Fusion
60
,
046028
(
2020
).
34.
F.
Reimold
,
M.
Wischmeier
,
S.
Potzel
,
L.
Guimarais
,
D.
Reiter
,
M.
Bernert
,
M.
Dunne
, and
T.
Lunt
, “
The high field side high density region in SOLPS-modeling of nitrogen-seeded H-modes in ASDEX Upgrade
,”
Nucl. Mater. Energy
12
,
193
199
(
2017
).
35.
P.
Manz
,
S.
Potzel
,
F.
Reimold
,
M.
Wischmeier
, and
ASDEX Upgrade Team
, “
Stability and propagation of the high field side high density front in the fluctuating state of detachment in ASDEX Upgrade
,”
Nucl. Mater. Energy
12
,
1152
1156
(
2017
).
36.
M.
Bernert
,
F.
Janky
,
B.
Sieglin
,
A.
Kallenbach
,
B.
Lipschultz
,
F.
Reimold
,
M.
Wischmeier
,
M.
Cavedon
,
P.
David
,
M. G.
Dunne
,
M.
Griener
,
O.
Kudlacek
et al, “
X-point radiation, its control and an ELM suppressed radiating regime at the ASDEX Upgrade tokamak
,”
Nucl. Fusion
61
,
024001
(
2021
).
37.
S. J. P.
Pamela
,
G. T. A.
Huijsmans
,
M.
Hoelzl
, and
JOREK Team
, “
A generalised formulation of G-continuous Bezier elements applied to non-linear MHD simulations
,”
J. Comput. Phys.
464
,
111101
(
2022
).
38.
C. W.
Gear
, “
The numerical integration of ordinary differential equations
,”
Math. Comp.
21
,
146
156
(
1967
).
39.
H. P.
Summers
and
R. W. P.
McWhirter
, “
Radiative power loss from laboratory and astrophysical plasmas. I. Power loss from plasmas in steady-state ionisation balance
,”
J. Phys. B
12
,
2387
(
1979
).
40.
H.
Summers
, see http://www.adas.ac.uk for “
The ADAS user manual
” (
2004
).
41.
E. T.
Meier
and
U.
Shumlak
, “
A general nonlinear fluid model for reacting plasma-neutral mixtures
,”
Phys. Plasmas
19
,
072508
(
2012
).
42.
C.
Dunlea
,
C.
Xiao
, and
A.
Hirose
, “
A model for plasma-neutral fluid interaction and its application to a study of CT formation in a magnetized Marshall gun
,”
Phys. Plasmas
27
,
062103
(
2020
).
43.
H.
Qin
,
S.
Zhang
,
J.
Xiao
,
J.
Liu
,
Y.
Sun
, and
W. M.
Tang
, “
Why is Boris algorithm so good
,”
Phys. Plasmas
20
,
084503
(
2013
).
44.
W.
Eckstein
, “
Sputtered energy coefficient and sputtering yield
,”
IPP-Report No. 17/29
,
2011
.
45.
See http://www.eirene.de/ for “
EIRENE—A Monte Carlo transport solver
” (
2019
).
46.
S. I.
Braginskii
, “
Transport processes in a plasma
,” in
Reviews of Plasma Physics
, edited by
M. A.
Leontovich
(
Consultants Bureau
,
New York
,
1965
), p.
1
.
47.
P.
Hunana
,
T.
Passot
,
E.
Khomenko
,
D.
Martínex-Gómez
,
M.
Collados
,
A.
Tenerani
,
G. P.
Zank
,
Y.
Maneva
,
M. L.
Goldstein
, and
G. M.
Webb
, “
Generalized fluid model of the Braginksii type
,”
Astrophys. J. Suppl. Ser.
260
,
26
(
2022
).
48.
S. Q.
Korving
,
G. T. A.
Huijsmans
, and
J. S.
Park
, “
The influence of drifts in the ITER divertor in nonlinear MHD JOREK simulations with fluid and kinetic neutrals
,” in
47th EPS Conference on Plasma Physics
(
2021
).
49.
V.
Kotov
, “
Numerical study of the ITER divertor plasma with the B2-EIRENE code package
,” Ph.D. thesis (
Institut Fur Energieforschung 4 - Plasmaphysik, Forschungszentrum
,
2007
).
50.
A. S.
Kukushkin
,
H. D.
Pacher
,
V.
Kotov
,
D.
Reiter
,
D.
Coster
, and
G. W.
Pacher
, “
Effect of neutral transport on ITER divertor performance
,”
Nucl. Fusion
45
,
608
(
2005
).
51.
S.
Carli
,
R. A.
Pitts
,
X.
Bonnin
,
F.
Subba
, and
R.
Zanino
, “
Effect of strike point displacements on the ITER tungsten divertor heat loads
,”
Nucl. Fusion
58
,
126022
(
2018
).
52.
P. C.
Stangeby
,
The Plasma Boundary of Magnetic Fusion Devices
(
Taylor & Francis Group
,
2000
), Secs. 2.8 and 2.10, pp.
92
95
, 98–105.
53.
A.
Loarte
, “
Effects of divertor geometry on tokamak plasmas
,”
Plasma Phys. Controlled Fusion
43
,
R183
(
2001
).
54.
A.
Loarte
,
R. D.
Monk
,
J. R.
Martín-Solís
,
D. J.
Campbell
,
A. V.
Chankin
,
S.
Clement
,
S. J.
Davies
,
J.
Ehrenberg
,
S. K.
Erents
,
H. Y.
Guo
,
P. J.
Harbour
,
L. D.
Horton
et al, “
Plasma detachment in JET Mark I divertor experiments
,”
Nucl. Fusion
38
,
331
(
1998
).
55.
S. K.
Erents
,
R. A.
Pitts
,
W.
Fundamenski
,
J. P.
Gunn
, and
G. F.
Matthews
, “
A comparison of experimental measurements and code results to determine flows in the JET SOL
,”
Plasma Phys. Controlled Fusion
46
,
1757
(
2004
).
56.
N.
Asakura
,
H.
Takenaga
,
S.
Sakurai
,
G. D.
Porter
,
T. D.
Rognlien
,
M. E.
Rensink
,
K.
Shimizu
,
S.
Higashijima
, and
H.
Kubo
, “
Driving mechanism of SOL plasma flow and effects on the divertor performance in JT-60U
,”
Nucl. Fusion
44
,
503
(
2004
).