A self-consistent multi-component model of plasma turbulence and kinetic neutral dynamics for the simulation of the tokamak boundary
A Coroado 1, P Ricci 11 École Polytechnique Fédérale de Lausanne (EPFL), Swiss Plasma Center (SPC), CH-1015
Lausanne, Switzerland.
[email protected]
Abstract
A self-consistent model is presented for the simulation of a multi-component plasma in the tokamak boundary. A deuterium plasma is considered, with the plasma species that include electrons, deuterium atomic ions and deuterium molecular ions, while the deuterium atoms and molecules constitute the neutral species. The plasma and neutral models are coupled via a number of collisional interactions, which include dissociation, ionization, charge-exchange and recombination processes. The derivation of the three-fluid drift-reduced Braginskii equations used to describe the turbulent plasma dynamics is presented, including its boundary conditions. The kinetic advection equations for the neutral species are also derived, and their numerical implementation discussed. The first results of multi-component plasma simulations carried out by using the GBS code are then presented and analyzed, being compared with results obtained with the single-component plasma model.
June 2021
1 Introduction
The boundary of a tokamak plays a crucial role in determining the overall performance of the device, as it sets the confinement of particles and heat, determines the heat exhaust to the vessel walls and controls the impurity level in the core [1]. The boundary is also the region where the plasma is fueled and helium ashes generated by fusion reactions are removed.
The tokamak boundary is characterized by the presence of several ion and neutral species that interact through a complex set of collisional processes [2, 1]. In particular, neutral atoms and molecules are relevant in the boundary as they result from processes such as the plasma recycling at the vessel walls and gas puffs. Recycling occurs because ions and electrons, transported along the magnetic-field lines by the parallel flow or across them by turbulent motion, eventually end at the vessel, where they recombine and re-enter the plasma as neutral particles, either being reflected, in which case they keep the energy of the original ion or, following an absorption process, being reemitted at the wall temperature. In case of absorption, a significant fraction of the atoms may associate to form molecules before being reemitted back to the plasma [3]. The exact probability of reflection or reemission, as well as the probability that atoms associate into molecules, depends on the physical properties of the limiter or divertor plate material [4]. At the same time, external injection of neutral molecules can be used to fuel the plasma, reduce the heat load on the vessel wall (e.g., by reducing the temperature of the plasma and hence inducing volumetric recombination processes), or diagnose the plasma.
The neutral atoms and molecules interact with the plasma in the tokamak boundary. Indeed, neutral atoms and molecules can be ionized, thus generating atomic and molecular ions, leading to a multi-component plasma. Molecular species also undergo dissociative processes that break them into mono-atomic species. Recombination, charge-exchange, elastic and inelastic collisions are also at play. These collisional interactions convert neutral particles into ions and electrons and vice versa, affect the temperature of the plasma species because of the energy required to trigger ionization and dissociation processes and modify the plasma velocity. As a result, since the dynamics of the plasma in the boundary are strongly influenced by its interaction with neutral species, it is important that simulations of the plasma dynamics in the tokamak boundary take into account its multi-component nature and the interactions between the different species to provide reliable quantitative predictions.
The description of a multi-component plasma is usually addressed by means of a fluid-diffusive model, which typically considers a version of the Braginskii fluid equations for the plasma species simplified by modelling cross-field transport through empirical anomalous transport coefficients. This approach is used by codes such as B2 [5, 6], EDGE2D [7], EMC3 [8], SOLEDGE-2D [9] and TECXY [10]. Sometimes, neutral particle species are also modelled using a diffusive fluid approach [11], for example in the UEDGE code [12]. However, diffusive models are no longer valid when the neutral mean free path is large, i.e. of the order of the plasma gradient scale length, which is often the case in the tokamak boundary. For this reason, neutrals are more commonly modelled by using a kinetic description valid for all ranges of mean free path. These models, typically based on Monte Carlo methods for the numerical solution, are implemented in the DEGAS2 [4], EIRENE [13], GTNEUT [14] and NEUT2D [15] codes. As a matter of fact, heat exhaust studies strongly rely on integrated neutral-plasma simulations of the tokamak boundary, which are most often based upon the coupling of the aforementioned fluid-diffusive models for the multi-component plasma and Monte Carlo-based models for the several neutral species, such as B2-EIRENE [13], EDGE2D-EIRENE [16], EMC3-EIRENE [17] and SOLPS [18]).
In order to shed light on perpendicular transport processes, simulations of plasma turbulence in the tokamak boundary have been carried out since a decade by using fluid and gyrofluid models, implemented in codes such as BOUT++ [19] (and its module Hermes [20]), FELTOR [21], GBS [22, 23], GDB [24], GRILLIX [25], HESEL [26] and TOKAM3X [27]. Kinetic models, implemented in other codes such as Gkeyll [28] and XGC1 [29, 30], have also been used. These simulations have allowed remarkable progress in the understanding of the mechanisms underlying turbulence and cross-field transport in a single-ion species boundary plasma. On the other hand, multi-component plasma simulations that include turbulent transport processes are still in their early days. Recent progress was made thanks to the synergy between the SOLEDGE2D [31] and TOKAM3X [27] codes. Based on the EIRENE Monte Carlo code for the kinetic simulation of the neutral species, multi-component plasma simulations are now enabled by the SOLEDGE3X code. The investigations carried out with SOLEDGE3X focused on the study of the dynamics of carbon impurities in the tokamak boundary [32]. Progress was also made by coupling the two-dimensional fluid code HESEL [26] with a one-dimensional fluid-diffusive model for the neutral particles, which accounts for both atomic and molecular species. The resulting nHESEL [33, 26, 34] code thus allows for the simulation of a single-ion plasma including the interactions with three neutral species: cold hydrogen molecules puffed into the system, warm atoms resulting from the dissociation of the hydrogen molecules and hot hydrogen atoms generated by charge-exchange processes. Such a model was used to study the plasma fueling in the presence of gas puffs and the formation of a density shoulder in the tokamak boundary at a high gas puffing rate.
In the present work, we describe the development and numerical implementation in the GBS code of a multi-component model that addresses the turbulent multi-ion species plasma dynamics through a set of fluid drift-reduced Braginskii equations, while each multiple neutral species are simulated by solving a kinetic equation. This work generalizes the implementation of the neutral-plasma interaction in GBS described in [35] for single-component plasmas. Single-component GBS simulations were used to study the electron temperature drop along the magnetic field [36] and to determine the influence of neutrals on gas puff imaging diagnostics [37]. The model has been improved recently through the implementation of mass-conservation by taking toroidal geometry consistently into account and by making use of particle-conserving boundary conditions to properly describe the recycling processes [38].
While the methodology presented in this work has the potential to include an arbitrary number of particle species and the corresponding complex scenarios, we consider a deuterium plasma, composed of five different particle species: three charged particle species, namely electrons (), monoatomic deuterium ions () and diatomic deuterium ions (), and two neutral species, namely deuterium atoms () and molecules ().
The model constitutes the first implementation of a kinetic multi-species model that avoids the statistical noise from the Monte Carlo method. In fact, the neutral kinetic equation, valid for any neutral mean free path, is solved by discretizing the kinetic equation integrated along the neutral path. The model has the potential to provide the fundamental elements necessary for the description and understanding of the key mechanisms taking place in the boundary, such as the fueling or gas puff imaging, where molecular species play an important role.
The results of the first simulation carried out with the multi-component model are also described in the present work, shedding some light on the processes underlying plasma fueling. In particular, for the limited configuration and sheath limited regime considered, we show that molecular dissociation processes have an impact on the location of the ionization source and plasma profiles, with respect to single-component simulations.
The outline for the present paper is as follows. After the Introduction, the collisional processes at play within the multi-component deuterium plasma model now implemented in GBS is presented in Sec. 2. In Sec. 3, we guide the reader through the derivation of the set of drift-reduced Braginskii equations used to describe a multi-component plasma, extending the approach previously followed by the single-component version of GBS. In Sec. 4, we present the boundary conditions we apply at the tokamak wall. The kinetic model for the neutral species is discussed in Sec. 5, where the numerical approach is also described, based on discretizing the kinetic equation integrated along the neutral path, in a generalization of the approach developed in [35] for a single neutral species model. Finally, in Sec. 6 we present and discuss the preliminary results of the first multi-component plasma GBS simulations, analyzing the impact of the molecules on the plasma dynamics, also presenting comparison to results from previous single-ion plasma simulations of GBS. The summary follows. In App. A, we derive the average energy of the reaction products and the average electron energy loss for the dissociative processes considered in the model. App. B presents the derivation of the friction and thermal force terms in the velocity and temperature equations used for the multi-component plasma description, following the Zhdanov closure [39] and considering the approach described in [32]. App. C features the list of kernel functions used to express the system of equations solved for the neutral species, while App. D presents the neutral system of equations in the matrix form implemented in GBS.
2 Collisional processes in multi-component deuterium plasmas
In the present paper we aim at describing an experimentally relevant multi-component deuterium plasma. Similarly to the works in [11], [40] and [41], the plasma we consider is composed of the , and species and we consider the D and neutral species. The molecules are present as the result of the association of atoms at the vessel walls and external injection. The molecules can be ionized, thus giving rise to ions, while dissociative processes are responsible for generating mono-atomic ions, , and neutrals, D, the later being possibly further ionized into . In general, the presence of ions is negligible in deuterium plasmas. Additionally, in the typical conditions of the tokamak boundary considered here, also the concentration of species that might be present in a deuterium plasma, such as the and ions, is negligible [11, 42, 43]. Considering five different species contrasts with the three species model used in the previous GBS simulations of a single-ion species plasma [35], where only mono-atomic deuterium ions and neutrals are evolved. We highlight that, by introducing the tools necessary to deal with the fundamental processes at play in multi-component plasmas, the present model can be further extended to describe more complex scenarios that include a larger number of plasma and neutral species.
The charged particle and neutral species are coupled by means of collisional processes, which include ionization, recombination, charge-exchange and dissociation processes, as well as electron-neutral collisions. These processes appear both in the neutral and plasma species model as particle and heat sources or sinks, as well as friction terms.
We henceforth list the collisional processes considered in our multi-component model, as well as their respective reaction rates, in Table 1. We remark that we neglect the distinction between fundamental and excited states for atoms, molecules and ions. In particular, we use the total cross section for each process considering the sum over the accessible electronic states of the reactants and products, following [42] and [43]. Based on momentum and energy considerations, we also compute the values of the velocity and energy of the collision products. Since these also depend on the electronic states of the reactants and products, we perform an average over the states relevant to a given reaction, taking into account the cross section of each state.
We denote with , and the modulus of the electron, and velocities, while , and represent their densities. The cross sections and refer to the collisions leading to the ionization of D and respectively, and are the cross sections for recombination of and with electrons, and stand for the cross sections of elastic collisions between electrons and D and respectively, and represent the dissociation cross sections of and , and are the cross sections for dissociative ionization of and , is the cross section of dissociative recombination of ions and, finally, , , and represent the cross sections for , , and charge-exchange interactions.
By considering Krook collision operators, the collision frequencies for ionization, recombination, elastic collisions and dissociative processes are computed as the average over the electron velocity distribution function, neglecting therefore the velocity of the colliding massive particle (D, , or ) when computing the relative velocity between the electron and the other particle. In fact, electrons have significantly larger thermal velocity than ions or neutrals. As for the charge-exchange interactions between ions and the neutral species D and , since the dependence of the cross section upon the ion-neutral relative velocity is weak [1], we neglect the velocity of the neutral particles (D or ) when evaluating the relative velocity of the colliding particles (we note that the velocity of a neutral particle is typically smaller than the velocity of the ions). Thus, we compute the reaction rates and by averaging over the distribution function of the species, which we assume to be a Maxwellian with temperature . Following the same approach when computing the cross section of charge-exchange interactions between ions and the and D neutrals, we average the cross sections and over the distribution function, assumed to be a Maxwellian of temperature .
Table 1: Collisional processes considered and their respective reaction rates.
Collisional process
Equation
Reaction Frequency
Ionization of D
Recombination of and
elastic collisions
Ionization of
Recombination of and
elastic collisions
Dissociation of
Dissociative ionization of
Dissociation of
Dissociative ionization of
Dissociative recombination of
Charge-exchange of
Charge-exchange of
Charge-exchange of
Charge-exchange of
We highlight that the values of the product for most of the reactions considered in Table 1 are obtained from the AMJUEL [42] and HYDEL [43] databases (precise references for each cross section are listed in Table 1 of [41]). While these databases list the cross sections for ordinary hydrogen plasmas, we assume here that they apply also to deuterium. More precisely, the cross section for the elastic collisions is obtained from [44] (page 40, Table 2), while for the elastic collision we use [45] (page 917, Table 13). The cross section for the charge-exchange reaction is taken from the HYDEL database (H.4, reaction 4.3.1), while for the charge-exchange we use the cross section values in the ALADDIN database [46], which are obtained from [47, 48]. For all the other reactions, we use the cross sections from the AMJUEL database [42]. The product for the collisional processes considered in this work is plotted as a function of the temperature in Fig. 1.
Figure 1: product for the collisional processes considered in this work. Ionization processes, elastic collisions and charge-exchange processes are displayed on the top panel, dissociative reactions on the bottom panel. The product is plotted as a function of the temperature of the colliding particle.
Having listed the collisional processes, we now focus on the velocity and energy of their products. For charge-exchange interactions of the kind , we assume that, while and exchange an electron, their velocities are not affected and energy is conserved. As a consequence, the ion is released from the charge-exchange collision with the velocity of , and is released with the velocity of . For the elastic collisions, given the large electron to deuterium mass ratio, we consider that the D velocity is not affected by the collision, while the electron is emitted isotropically in the reference frame of the massive particle according to a Maxwellian distribution function, , centered at the velocity of the incoming D particle, . The temperature is established by energy conservation considerations. Precisely, we observe that the average energy of the incoming electrons consists of the sum of the kinetic energy associated with the electron fluid velocity, , and the thermal contribution, given by . On the other hand, the energy of the outcoming electrons has a contribution given by the collective re-emission velocity, , and a thermal contribution, . It follows that satisfies the following balance, . The elastic collisions between electronds and can be described similarly. The re-emitted electrons have a distribution , with obtained from an analogous conservation law, .
We now consider the electrons generated by ionization of D. We assume that they are described by the Maxwellian distribution function centered at the fluid velocity of the D atom with that takes into account the ionization energy loss, , whose value is presented in Table 2. More precisely, satisfies the energy conservation law, , as the reaction gives rise to two electrons with the same properties. The same approach is followed for the ionization of , with the two released electrons being described by a Maxwellian centered at the velocity of the molecules, , and with temperature obtained from , with the average energy loss due to ionization of (see Table 2). We highlight that we neglect multi-step ionization processes when computing the cross section for ionization of D and , and we do the same for all other electron impact-induced reactions, such as the dissociative processes considered here.
We apply the procedure used for ionization processes to describe the properties of the electrons resulting from dissociative processes, with the electron generated by dissociation of being described by the Maxwellian centered around and with temperature obtained from . Regarding dissociation of , the resulting electron is similarly modelled by a Maxwellian centered at the velocity of the ion and with temperature given by energy conservation, . On the other hand, dissociative ionization of generates two electrons, whose Maxwellian distribution function, , is centered around the velocity, , and characterized by a temperature , obtained from . Similarly, the electrons generated by dissociative ionization of are assumed to follow a Maxwellian centered at and with temperature obtained from the corresponding energy conservation law, .
The evaluation of the temperature of the D atoms and ions released from dissociative reactions are based on modelling these reactions as Franck-Condon dissociation processes. These temperatures are summarized in Table 2 and rely on data from [43]. The detailed calculations are presented in the App. A. We also remark that these particles are emitted isotropically in the frame of the centre of mass of the incoming or particle. Thus, the D atoms generated in dissociation of molecules, for instance, are assumed to have a Maxwellian distribution . Analogously, we describe the neutral D atoms and ions generated by dissociative-ionization of molecules by the Maxwellian distributions and respectively, with the temperature listed in Table 2 and evaluated in App. A. Similarly, and are the Maxwellian distributions of D atoms and ions generated by dissociation of ions, where is the fluid velocity of the ion population that includes the leading order components (see Sec. 3). Finally, dissociative-ionization of generates ions that are described by a Maxwellian distribution . To conclude, we note that the D atoms and generated by dissociative-recombination of are described by the Maxwellian distributions and respectively, with the average thermal energy of the reaction products.
Table 2: Average electron energy loss and average energy of reaction products for the ionization and dissociative processes included in the model.
Collisional process
Electron energy loss
Reaction product temperature
Ionization of D
————————————
Ionization of
————————————
Dissociation of
Dissociative-ionization of ()
Dissociative-ionization of ()
Dissociation of
Dissociative-ionization of
Dissociative-recombination of
————————————
3 The three-fluid drift-reduced Braginskii equations
The kinetic equations for , and , which include the terms associated with the neutral-plasma interactions, are the starting point for the derivation of the Braginskii set of equations, used here to model the plasma dynamics. These equations generalise the ones considered in the single-ion species model described in [23, 35], by adding the new collisional neutral-plasma interaction terms listed in Table 1, as well as an equation for the description of molecular ions, . The kinetic equations are
(1)
(2)
and
(3)
In Eqs. (1-3), is the particle velocity, is the particle acceleration due to the Lorentz Force, is the gradient in real space and in the velocity space. The , and terms represent Coulomb collisions between charged particles affecting the , and distribution functions, respectively.
The Braginskii equations for the 3-species plasma (, and ) are then obtained by taking the first three moments of the kinetic equations for each species in the limit , with the cyclotron frequency ( denotes the ion mass and is the elementary charge) and the characteristic Coulomb collision time for ions. The Braginskii equations, including the neutral-plasma interaction terms, can be derived by following the steps presented in [49], and take the following form
(4)
(5)
(6)
(7)
(8)
(9)
(10)
(11)
(12)
where is the component of the stress tensor along the and directions, is the friction force acting on the electrons, is the electron heat flux density, is the electron heat generated by Coulomb collisions and is the electron advective derivative. The equivalent notation is used for the and species.
The drift-limit of the Braginskii equations is finally derived by applying the ordering, valid in typical conditions of the tokamak boundary. Only leading order components in are retained in the electron perpendicular velocity, i.e. , with the drift and the electron diamagnetic drift, thus neglecting electron inertia. Similarly, the perpendicular velocity is decomposed as , where the leading order perpendicular velocity,
(13)
is the sum of the drift and the diamagnetic drift, . The polarization drift,
(14)
is of higher order than in the expansion, as shown in [49]. Similarly, the drift arising from friction between ions and other species,
(15)
is also of higher order in . This term includes contributions from collisions of with D, and particles. Assuming , and , and noticing that , one obtains . For this reason, and are approximated with their leading order components, i.e. and , in Eq. (15). In Eqs. (14) and (15) we introduce the giroviscous term for ions , the viscosity , the magnetic field curvature vector , the gradient along the magnetic field , the gradient perpendicular to the magnetic field , and the magnetic field unit vector .
For the derivation of the drift-limit of the velocity, we follow a similar approach, as the ordering is also valid in typical tokamak boundary conditions. The ions perpendicular velocity is thus given by , with
(16)
the leading order component, with . The velocity denotes the polarization drift and stands for the drift velocity arising from friction between ions and other species. Their expressions are given by
(17)
and
(18)
with the giroviscous term and the related viscosity. The approximation is used in Eq. (18).
To obtain an expression for the parallel friction forces and parallel heat fluxes and close the Braginskii equations, we use the collisional closure proposed by Zhdanov in [39], leveraging the formulation presented in [32, 50], more suitable for numerical implementation. The application of this procedure to the particular case of the multispecies plasma considered here is described in App. B, where we take advantage of the fact that the density is considerably smaller than the density, i.e. , for typical tokamak boundary conditions, which leads to because of quasi-neutrality. On the other hand, the contributions from the perpendicular heat flux arising from and in the and equations, respectively, can be evaluated as in the single-ion species model [22, 23], in particular following the derivation presented in [49]. This approach can be generalised to evaluate the term arising from the perpendicular component of in the equation. Thus, the drift-reduced Braginskii system of equations is composed of the continuity equation for the electron species, the continuity equation for the species, the vorticity equations that ensures quasi-neutrality, , and the equations for the parallel velocities and temperature of all species. They take the form
(19)
(20)
(21)
(22)
(23)
(24)
(25)
(26)
and
(27)
In Eqs. (19-27) we introduce , and the plasma vorticity , with the contribution given by and an analogous contribution, . The system is thus closed by the generalized Poisson equation, which is obtained by inverting the definition of the plasma vorticity, , yielding
(28)
We remark that the electron giroviscous term in Eq. (22) is defined similarly to the ion giroviscous terms, . Eq. (21) is written avoiding the Boussinesq approximation and taking into account all components of the velocity of the ion species and , including the higher order polarization and friction contributions. On the other hand, in order to express the advective derivative for the ion species, and , we only consider the leading order components of the perpendicular velocity, and , therefore neglecting and . Similarly, we neglect the friction and polarization drifts in the continuity equation for . We also remark that the terms of higher order in in the perpendicular velocity of ions are neglected when writing in the temperature equations, Eqs. (26) and (27), which is a necessary assumption in order to avoid explicit time derivatives arising from the polarization drift velocity, . Nevertheless, all terms are considered in the divergence of the perpendicular velocity of ions in Eq. (26), as we make use of to write in terms of and . Finally, when taking the divergence of these terms, we make use of to neglect the contribution of the velocity of D atoms, which is valid since (with the sound Larmor radius of ions, the ions sound speed and the mean free path of a D atoms). This relation can be generalized to the other neutral and ion species, namely molecules and ions. Thus we neglect the contribution of the divergence of neutral particle velocities when compared to the divergence of ion velocities.
We note that dimensionless units are used in Eqs. (19-27) and in the rest of the paper. The densities, , and , are normalized to the reference value , while temperatures, , and , are normalized to the respective reference values, , and , which are related by the dimensionless quantity . Conversely, lengths parallel to the magnetic field are normalized to the tokamak major radius, , lengths perpendicular to the magnetic field are normalized to the ion sound Larmor radius, , where is the normalized ion sound speed and is the ion cyclotron frequency at the magnetic axis, and time is normalized to . All other normalizations follow, namely the parallel velocities, , and , are normalized to , the plasma vorticity is normalized to , perpendicular diffusion coefficients and conductivities are normalized to , while the parallel diffusion coefficients and conductivities are normalized to . Normalized quantities are used in the rest of the paper, except when explicitly mentioned. The parameter is the ratio between the ion sound Larmor radius and the tokamak major radius . We also note that is the dimensionless resistivity given by , with the parallel conductivity defined in terms of the electron characteristic time as .
We conclude with a few final remarks on Eqs. (19-27). We first note that the parallel conductivity appearing in the temperature equations for electrons is expressed in the form , where we retain the Spitzer temperature dependence while we neglect the weaker space and time variation of the factor, similarly to the approach followed in the single-component plasma model previously implemented in GBS ([22, 23]). A similar approach is followed for and . This is not expected to impact the simulation results in the sheath-limited regime, where conductivity-related contributions are small. Finally, we point that, since the density may drop to a very low value, numerical issues may arise in the equations for and due to terms featuring a dependence. For a robust numerical approach, we evolve the parallel flux and pressure of the ion species, and , instead of and . The equations for the time evolution of and are
(29)
and
(30)
with , and given, respectively, by Eqs. (20), (24) and (27). We focus on the parallel flux, , and pressure, , when presenting the simulation results.
4 Boundary conditions
The boundary conditions implemented in the previous GBS models for single-ion species plasma are extended in the present work to include the molecular ion species . In the case considered here, of a plasma with a toroidal limiter, the domain boundary includes the limiter plates, the outer wall and the interface with the core, where the low plasma collisionality questions the application of a fluid model.
We first consider the boundary conditions at the limiter plates, where most of the plasma ends by flowing along the magnetic field lines. Those are the most important boundary conditions to impact the simulation dynamics. The boundary conditions are imposed at the interface between the collisional pre-sheath (CP) and the magnetic pre-sheath (MP), derived from the Bohm-Chodura boundary conditions, following the approach described in Ref. [51] in the cold ion limit and generalized in Ref. [52] to account for finite ion temperature. Here, we further extend this procedure to the case of a multi-ion species plasma. For this purpose, we use the coordinates, with the direction of the magnetic field, the direction perpendicular to the magnetic field and parallel to the limiter surface, and the direction perpendicular to both and (here all spatial coordinates are normalized to , while the other quantities are normalized as in Eqs. (19-27)). We also introduce , the coordinate perpendicular to the limiter plate, with the angle between the magnetic field line and the plane of the limiter.
As a first step in the derivation of the boundary conditions, we note that the steady-state dynamics of the multispecies plasma in the CP is described by means of the continuity equation for the and species (quasi-neutrality provides the electron density) and the parallel momentum equations for , and . In steady state, these can be written as
(31)
(32)
(33)
(34)
and
(35)
with , and the particle sources for and , and , and the momentum sources for , and .
From Eqs. (31-35) and following the approach described in [52], a system of five equations for , , and , is hence obtained for the interface between the CP and the MP border, considering the limit and isothermal ions and electrons. For this purpose, at the MP entrance, gradients along the direction are assumed weaker than gradients along by a factor , with , and respectively the scale lengths of , and along the direction. In addition, finite Larmor radius (FLR) effects are neglected and, to express the and components of the velocity of each ion species, and , we consider only the leading order terms in (see Eqs. (13) and (16)). This yields
(36)
(37)
(38)
and
(39)
where and are respectively the and components of the drift velocity, and are the and components of the diamagnetic velocity and and are the and components of the diamagnetic velocity. Finally, we define the velocity of the ions along the direction as . We also introduce the velocity of the ions along the direction that excludes the diamagnetic contribution, that is and for the ions. The system in Eqs. (31-35) yields
(40)
(41)
(42)
(43)
and
(44)
where , and . We then make use of the quasi-neutrality condition, , to obtain a system of five linear equations that we express in matrix form as , with
(45)
(46)
and
(47)
Following [51, 52], we observe that, while the source terms are important in the CP, they are small at the MP entrance with respect to the gradient terms. This allows one to assume at the MP entrance. Thus, the linear system reduces to at the MP entrance. We solve with respect to to obtain the non-trivial solution at the MP entrance. For this purpose, following [53], the parallel velocity of the ion species, , is related to ,
(48)
In addition, we assume (and therefore ) and keep only zero order terms in , neglecting therefore all derivatives along the direction. The condition then yields
(49)
where the signs refer to the magnetic field lines entering/leaving the vessel and we have defined . We now note that , since we neglect , to obtain the boundary condition for at the limiter,
(50)
The expressions of the boundary conditions for the other plasma quantities then follow. In fact, Eq. (42) can be inverted to express in terms of , which yields
(51)
We then use Eq. (44) to express in terms of , that is
(52)
and, making use of , we also obtain
(53)
Regarding the density of the ions, we use Eq. (41), deriving the following boundary condition
(54)
In order to derive the boundary conditions for , and , we notice that temperature gradients along the direction perpendicular to the wall are small compared to the gradients of the other physical quantities. In fact, [51, 52] show that . In the present work, we follow this prescription and assume (we note that our tests show that imposing does not affect the simulation results noticeably).
To obtain the boundary condition for at the MP entrance, we start from its definition, . We write the second order derivatives in the directions perpendicular to the magnetic field retaining only derivatives along the direction, since . Since at the limiter, the factor can be considered constant when the derivatives defining are evaluated. We then write the derivatives along the direction in terms of derivatives along and consider (for simplicity). This yields
(55)
We now take advantage of Eqs. (52) and (54) to express and in terms of and use Eq. (51) to obtain the final expression of the boundary condition for , that is
(56)
Finally, the boundary condition for the electron parallel velocity is obtained from the analysis of the electron kinetic distribution function at the MP entrance. As discussed in [51], this gives
(57)
where .
At the vessel outer wall and the core interface, boundary conditions are considered, similarly to the approach used in previous models of GBS [51, 52, 23]. In fact, a set of first-principles boundary conditions is yet to be derived for such boundaries. The impact of these boundary conditions upon the simulation results is controlled by extending radially the simulation domain towards the wall and the core. The conditions we impose include homogeneous Neumann boundary conditions to , , , , , , and . Since the density of ions is expected to be very low at the core-edge interface (no particles outflowing from the core), we use Dirichlet boundary conditions at the core interface for , setting it to a residual value, while homogenous Neumann boundary conditions are considered at the vessel outer wall. We also use Dirichlet boundary conditions for the vorticity, setting at both the wall and the core interface. Regarding the boundary conditions, we follow the approach presented in Ref. [54], where is considered at the vessel wall. Finally, is considered at the core interface, where is a constant value chosen to prevent large gradients of .
5 The kinetic model for the neutral species and its formal solution
In order to compute the neutral distribution functions of D and , and , we consider a set of two coupled kinetic equations, that is
(58)
and
(59)
The formal solution of Eqs. (58) and (59) can be obtained by using the method of characteristics, assuming that the plasma quantities are known. This yields
(60)
and
(61)
The solutions presented in Eq. (60-61) show that the distribution functions of D and at position , velocity and time result from the neutrals generated at a location , in the plasma volume or at the boundary, and at time , with the unit vector aligned with the neutral velocity and the distance measured from to (the subscript "b" denotes the intersection point between the domain boundary and the characteristic starting at with direction ). Since the neutrals are solved on the coordinate system, with the distance from the torus axis and the vertical coordinate measured from the equatorial midplane, the integral includes the Jacobian corresponding to the coordinate system . The volumetric source due to the collisional processes in Eq. (60) is
(62)
since D ions can be generated in the plasma volume by and charge-exchange interactions, recombination of ions with electrons, dissociation of molecules into two D atoms, dissociative ionization of into D and , dissociation of ions into D and , and dissociative recombination of into two D atoms.
Similarly, molecules can be generated in the plasma by and charge-exchange interactions or recombination of ions with electrons. Therefore, the volumetric source term in Eq. (61) is
(63)
We remark that is a Maxwellian distribution function describing the ion population, centered at the ion velocity , which includes only the leading order components, i.e. , and based on the temperature, . Similarly, is a Maxwellian distribution that describes the ions and is defined analogously. We remark that, when evaluating the average velocity of the Maxwellian distributions describing neutrals generated from and , we assume that and can be neglected, i.e. and . Regarding the dissociative processes, we recall that the temperature is the average thermal energy of D atoms generated by dissociation of , presented in Table 2 and calculated in App. A. The energy of the neutral D atoms generated by other dissociative processes is evaluated using a similar approach.
The effective frequencies for depletion of neutral particles are given by
(64)
and
(65)
since the volumetric sinks of D atoms are due to ionization or charge-exchange with or , while are depleted by ionization, charge-exchange with or , dissociation or dissociative ionization.
A contribution to the neutral distribution functions in Eqs. (60) and (61) is related to the neutral recycling at the boundary walls. Therefore, we now focus on the neutral processes that take place there. A fraction, , of the ions that reach the boundary walls, after recombination with electrons and formation of neutrals is reflected back into the plasma. The remaining fraction, , is absorbed and reemitted at wall temperature as , also following a recombination process. Analogous considerations hold when describing neutrals that reach the boundary. The molecules can, in fact, be reflected or reemitted with the same probability as .
Regarding the atomic species, since the wall temperature is low, a fraction, , of the and D particles absorbed at the walls associate and are reemitted back into the plasma as molecules. The ions and D neutrals reaching the boundaries that do not associate undergo reflection and reemission processes similar to the ones described for ions and particles, the probability of reflection, , being the same. As a consequence, the distribution functions at the vessel, and , for (with the unit vector normal to the boundary) yield
(66)
and
(67)
We first analyse the contributions of reflected particles in Eqs. (66-67). The reflected D and are described by the distribution functions and , since is the velocity of the reflected neutrals as they flow towards the wall, with the velocity along the direction normal to the wall surface. On the other hand, the contribution from the reflected and is modelled by considering the projection of the flux of outflowing and along the direction normal to the boundary surface, given respectively by and . These fluxes include the contributions of the plasma parallel flow and the leading order perpendicular drifts, i.e. the and diamagnetic drifts, yielding
(68)
and
(69)
We assume that the velocity distribution of the D neutrals generated by reflection of ions is described by a Maxwellian centered at the velocity, , with , and with temperature of the incoming ions, , given by . Analogously, the neutrals generated by reflection of ions are described by a Maxwellian distribution of velocities, , being , with and the temperature of the incoming ions.
We now focus on the contributions in Eqs. (66-67) accounting for reemission of neutrals from the boundary. These are written in terms of
(70)
and
(71)
In addition to the projections of the boundary ion fluxes, and , Eqs. (70) and (71) take into account the projections along the direction normal to the boundary of the fluxes of D and outflowing to the limiter and walls, and . These are defined based on the neutral fluxes directed towards the boundary (i.e. for ) as
(72)
and
(73)
We assume that the velocity distribution of reemitted particles follows the Knudsen cosine law for a given wall temperature, . This yields, for the D neutrals,
(74)
while the expression for molecules is analogously given by
(75)
We now follow the same approach described in [35] to obtain a set of time-independent two-dimensional integral equations for the D and densities, making the numerical implementation of the formal solution in Eqs. (60) and (61) feasible. More precisely, we first make use of the fact that the neutral time of flight is typically shorter than the characteristic timescales of turbulence, , a condition that we denote as the neutral adiabatic regime. This allows us to approximate in Eqs. (60-67) or, equivalently, and in Eqs. (58-59). Second, we note that the neutral mean free path is typically smaller than the characteristic elongation of turbulence structures along the magnetic field, . Therefore, our description of neutral motion is reduced to the analysis of a set of independent two-dimensional planes perpendicular to the magnetic field, approximately coincident with the poloidal planes. Then, integrating Eqs. (60-61) over the velocity space, a system of two coupled equations for the densities of D and is obtained,
(76)
and
(77)
where the same geometrical arguments presented in [35] is used when considering the integral along the neutral path and the integral along the perpendicular velocity angle, that is
(78)
where is the area element in the 2D poloidal plane and is a generic function. In addition, we use
(79)
with being a line element along the boundary of D, denoted as , and .
We now express the volumetric source terms appearing in Eqs. (62) and (63), and , in terms of and , and the distribution functions of the neutral species at the boundary appearing in Eqs. (66) and (63), and , in terms of , , and . For , this yields
(80)
while for one has
(81)
Replacing in Eqs. (73) and (72), the normal projections of the fluxes of and D can be written respectively as and . By replacing and by their expressions as given in Eqs. (66) and (67), these fluxes can be rewritten in terms of , , , , and as
(82)
and
(83)
We note that that the neutral particle densities and fluxes in Eqs. (80-83) are multiplied by a factor in order to account only for the contribution of particles that are reemitted at the boundary, hence excluding reflection. Neutral reflection is included, in the definition of the kernel functions that appear in Eqs. (76-77).
We now turn to the definition of the kernel functions appearing in Eqs. (80-83). These are defined as integrals over velocity space. For instance, quantifies the amount of neutrals found at a location in the plasma volume () being generated from collisions involving neutralization of ions at a location inside the plasma volume (). Its expression is given by
(84)
which separates the contributions to arising from the direct path of length connecting and , , and the path corresponding to the trajectory of neutrals that are reflected at the boundary, . Both and have the same expression,
(85)
where and is the distance between and measured along the path (for the direct trajectory is given by the distance between the two points along a straight line, while for the reflected trajectory is the sum of the distances between and the boundary and the distance from the boundary to ). We remark that is the integral along the parallel velocity of the Maxwellian distribution function, . We also remark that in Eq. (85) is valid in case the points are optically connected, i.e. if the straight line connecting the two points does not cross the core region nor the limiter plates. Otherwise, if the points are not connected, one has . As for , in the present work we assume no reflection at the outer walls, while reflection of ions and neutrals may take place at the limiter plates. The other kernels appearing in Eqs. (76-77) have the same structure as , and they take into account possible direct and reflected paths connecting the two points. These kernels are presented in detail in App. C.
We now turn to the evaluation of the non-homogeneous terms appearing in Eqs. (80-83), i.e. the terms that are not proportional to nor . For instance, these terms include the contribution of the ions recycled at the wall. In fact, the reflection and reemission of ions that outflow to the boundary and recombine with electrons contribute to the density of neutral D atoms, through the term
(86)
where is defined in Eq. (68). Similarly, the recombination of ions with electrons at the walls that are then either reflected or reemitted as , and the recombination of ions with electrons at the walls and the following association into molecules contribute to the density of the species. These contributions can be expressed as
(87)
and
(88)
We also define the non-homogeneous terms appearing in Eqs. (82-83), that provide the contributions to the flux of neutrals at the boundary, and , given by the ions outflowing to the wall. Following a similar approach to the one presented above, this can be expressed as
(89)
(90)
and
(91)
We now turn to the evaluation of the contributions to the neutral particles appearing in Eqs. (80-83) caused by volumetric processes that involve the ion species and . The contribution to the density as a result of recombination processes is given by
(92)
and the contribution to the flux of to the boundary, also associated to recombination events, is expressed as
(93)
Similar contributions from volumetric recombination processes are considered for the D neutral species. The contribution to the D density as a result of recombination yields
(94)
while an analogous definition is used for the flux of D,
(95)
Finally, the contribution of the dissociation of ions to appearing in Eq. (81) is evaluated as
(96)
Similarly, a dissociation process of ions results in a contribution to in Eq. (83) given by
(97)
For their numerical solution, the system of kinetic equations for the neutral species is discretized on a regular cartesian grid in the plasma and then written in matrix form. The details of the numerical implementation of the neutral model are discussed in App. D.
6 First simulation of a multi-component plasma with the GBS code
We present the first results from simulations of turbulence in the tokamak boundary carried out by using the multi-component plasma model described in Secs. II-IV and implemented in the GBS code. Similarly to [22, 35], we consider a tokamak with an infinitesimally thin toroidal limiter at the HFS equatorial midplane, with major radius , and we simulate a three-dimensional domain with an annular cross section that includes the edge and the open-field line region of the device. The radial size of the domain is and a the poloidal size is at the core interface. Since the limiter has a radial width of ), both the open and closed field-line regions have a radial extension of , corresponding to half the size along the radial direction. The parameters chosen for the present simulation are , , , , , , , , , , , , , , , , , , , and , . Regarding the reflection probability at the limiter, we remark that it depends strongly on the particle energy and the wall material (see [1]). In this simulation, reflection of ions and neutrals takes place at the limiter plates with a given probability , constant along the limiter surface. The fraction of reflection at the boundary is therefore defined as
(98)
We choose to consider metallic boundaries and hence we assume , a value similar to the one adopted in [35]. We also assume , which is consistent with the usual assumption that most D atoms associate into molecules at the boundary (see e.g. [3, 55]).
Regarding the numerical parameters, we note that the plasma grid resolution is while neutral grid resolution is . The time step is and the neutral quantities are evaluated every . Although we have not carried out convergence studies with the multispecies model presented in this paper, convergence on plasma and neutral grid refinement has been studied within the single component framework. The conclusions presented in [56], which we expect to remain valid in the multispecies model presented here, show that our results are converged with respect to the frequency of neutral calculation.
For the description of the simulation results we focus on the quasi-steady state regime, established after a transient, when the plasma and neutral profiles fluctuate around constant values. We take toroidal and time averages of the plasma quantities evolved by Eqs. (19-27) and (28) over a time interval of . These quantities are shown in Fig. 2 on a poloidal cross section. In Fig. 3, we present the density of the neutral species, and , and the neutral-plasma collisional interaction terms taken into account in our model. The results of the multispecies simulations are compared with the one of a single-component plasma, with corresponding parameters. The time and toroidal averages of the plasma and neutral main quantities for the single species simulation are shown in Fig. 4.
Figure 2: Cross section plots of the electron density (), density (), density (), electron parallel velocity (), parallel velocity (), parallel velocity (), electron temperature (), temperature (), temperature () and electrostatic potential (), toroidal and time-averaged over an interval of from the quasi-steady state of the multi-component plasma simulation described in Sec. 6.Figure 3: Cross section plots of the neutral species densities and source terms resulting from the neutral-plasma interaction, toroidal and time-averaged over an interval of from the quasi-steady state of the multi-component plasma simulation described in Sec. 6.
We first focus on some general considerations on the plasma and neutral densities. The plots in Figs. 2 reveal that the density of the molecular ion species is three to four orders of magnitude smaller than the density of the main ion species , a result in agreement with the assumption used in Eqs. (22-27) for the derivation of the parallel friction and heat flux terms and in Eqs. (40-44) to obtain the boundary conditions at the limiter. We highlight that the density of peaks just inside the LCFS next to the limiter, since most of the molecules cross the open-field line region without interacting and are then dissociated and/or ionized by the denser and warmer plasma inside the LCFS. As a matter of fact, exhibits a similar behavior to the profile of the molecular ionization source presented in Fig. 3, which also peaks in the edge near the limiter. On the other hand, Fig. 2 shows that and are comparable to near the limiter plates, while they are about one order of magnitude smaller than in the rest of the SOL and up to two orders of magnitude smaller inside the LCFS. Furthermore, regarding the relative importance of D and , Fig. 3 shows that is larger than by a factor between two and three in the open-field line region around the limiter, while is larger than inside the LCFS at the HFS, as a consequence of the higher plasma densities and temperatures that lead to the dissociation of molecules in that region.
As a second set of observations, we focus on the asymmetry of the plasma density and flow. An up-down asymmetry in the edge region is shown by the profiles of and , which are noticeably larger below the equatorial midplane than above it. The underlying reason of this asymmetry can be inferred from the , and profiles. In fact, the and parallel flows are directed in the counterclockwise direction in the edge region. Therefore, the ionization of neutrals inside the LCFS, which occurs mostly in the proximity of the limiter at the HFS, leads to plasma particles subject to a downward flow. Albeit being small, this flow leads to a slightly larger density of and below the equatorial midplane of the device. The parallel flux of ions is also directed counterclockwise in the edge at the HFS, which further enhances this mechanism, even though densities are small compared to the other species. We highlight that the and profiles are slightly different when the single-component model of GBS is considered, as illustrated in Fig. 4. In this case, although it is also observed an up-down asymmetry in the profile, this is related to the fact that the ionization source, , is larger in the edge region below the limiter than above it, due to larger recycling rates at the lower limiter plate. In fact, contrary to the multispecies case, the is characterized by a counterclockwise parallel flow of ions in the edge below the midplane, while above it the parallel flow is directed clockwise.
In the multi-component plasma simulation we also observe a larger parallel flow of plasma in the open-field line region towards the upper side of the limiter when compared to the lower side. This can be observed in Fig. 2, that shows larger and above the limiter plates than below it, ultimately leading to higher recycling rates and hence larger and densities in the region above the limiter, as shown in Fig. 3. The reason behind this behaviour is again related to the profile. In fact, while the ions flow counterclockwise along the magnetic field lines in the edge, they undergo cross-field transport towards the SOL. As a result of the counterclockwise parallel flow and related asymmetry of in the edge region, most ions cross the LCFS above the equatorial midplane, while flowing along the magnetic field lines towards the upper side of the limiter.
Figure 4: Cross section plots of plasma density , ion parallel velocity , ion temperature and ionization source term , toroidal and time-averaged over an interval of from a quasi-steady state single-component plasma simulation. The grid sizes and simulation parameters are the same as the ones considered in the multi-component simulations, except for the wall re-emission temperature, which is set to , to mimick Franck-Condon dissociation processes, and .
It is also observed that is slightly larger in the HFS compared to the LFS, which is due to the existence of sources in the HFS around the midplane. Similarly, also in the single-species simulation, is larger in the HFS as a consequence of the ionization source, .
Focusing on the temperature of the plasma components, we observe that the profile presents a similar behaviour to the one observed in single-component plasma simulations. A clear asymmetry between the HFS and the LFS is observed for , which is qualitatively similar to the results for a single-component simulation in Fig. 4. As a matter of fact, the temperature is considerably lower on the HFS compared to the LFS, which is related to the generation of cold ions inside the LCFS due to ionization of D atoms, dissociative processes and charge-exchange interactions. This effect is particularly important above the limiter, where the recycling rates are larger. On the other hand, the profile of exhibits a maximum inside the LCFS at the HFS, where the majority of the ions are generated by ionization of molecules coming from the limiter. The up-down asymmetry of the pressure around the limiter is also due to the asymmetry of the recycling rates. As an aside note, we remark that, since it is strongly related to the profile [57], the electrostatic potential profile revealed by the multi-fluid simulations is similar to the one observed in the single-component plasma model.
Analyzing the neutral-plasma interaction terms presented in Fig. 3, we first notice that ionization processes tend to be more important in the edge region at the HFS, with atomic and molecular ionization rates exhibiting similar profiles. However, peaks in the vicinity of the LCFS, while peaks further inside the LCFS and has a larger radial spread. In fact, molecules generated in the open-field line region and are dissociated and/or ionized in the proximity of the LCFS, where the plasma is warmer and denser. In contrast, although most D atoms are generated in the open-field line region, they are also created by dissociation of molecules in the edge. This shifts the maximum of radially inwards and makes the ionization source spread across a wider area. As a result of being larger than in the edge, the maximum of is also almost two times larger than .
Focusing on the electron-neutral collisions, we note that the reactions involving occur more often in the open-field line region, mainly in the area surrounding the limiter plates where the majority of the neutral molecules are generated. Reactions with become less important in the edge, since most molecules are dissociated and/or ionized due to the higher densities and temperatures. On the other hand, electron-atom collision reactions involving the D species peak inside the LCFS, because the cross sections of these reactions are larger in the edge region due to the higher plasma density and temperature and because of the presence of D atom resulting from dissociative processes. We also highlight that elastic collisions and charge-exchange reactions are more frequent on the upper side of the limiter, in agreement with the strong up-down asymmetry discussed above. Regarding charge-exchange reactions, we observe that they are spatially localized similarly to the electron-neutral collisions. The reactions between the two molecular species ( collisions) occur less often than the charge-exchange between mono-atomic species ( collisions) by three to four orders of magnitude, which is a result of the to ratio. In addition, the terms arising from charge-exchange interactions between molecules and ions ( collisions) are found to be two orders of magnitude smaller than the ones between the atomic species ( collisions) in the region of the domain where these interactions are important. In turn, charge-exchange between ions and D atoms is three orders of magnitude smaller than charge-exchange, which is due to the fact that is proportional to .
Finally, we analyse the dissociative processes, which represent a sink of molecular species and and sources of D atoms and ions. Simple dissociation of and , described by the terms and respectively, which do not involve ionization nor recombination processes, are found to be dominant dissociation processes, and occur with a frequency similar to that of the ionization of D and . We remark that dissociation of molecules peaks just above the limiter plate (where most molecules are generated) and in the edge region, in the vicinity of the LCFS, and then it is significantly smaller in the core, since drops rapidly across the edge. In contrast, dissociation of ions is very small in the open-field line region, where the density of is negligible (at the typical electron temperature of the SOL, the ionization cross section is small), and is important only inside the LCFS, where ions are generated. The profile therefore closely follows the profile, with a larger radial spread when compared with the dissociation of . As for dissociative ionization of and , and respectively, we observe that the rates are smaller by one to two orders of magnitude with respect to the simple dissociation of and and peak in the edge region a bit further inside. This is due to the fact that the energy required to trigger dissociative ionization processes is considerably larger than the one needed to dissociate the particles without triggering an ionization process, as shown in Table II. Hence, these processes are only relevant in the edge region, where densities and temperatures are sufficiently high to make these cross sections significant. This is particularly the case of , since this term is also proportional to the density of ions, which is relevant only inside the LCFS. Nevertheless, we highlight that these reactions become considerably less important towards the core, as very few and cross the edge region without being dissociated. As for dissociative-recombination of particles, , its amplitude is also smaller than that of simple dissociation by one to two orders of magnitude and follows very closely the profile, since there is no energy threshold to trigger the reaction, unlike dissociative ionization processes.
These results allow us to draw a global picture of the main processes determining the dynamics of neutrals in the boundary. Although some molecules are dissociated in the SOL region, most of them cross the LCFS and are dissociated into D atoms within a short distance as they get in contact with the warmer and denser plasma of the edge. The remaining molecules penetrate further towards the core and are ionized by the increasingly warmer and denser plasma, giving rise to ions, which in turn are quickly dissociated into ions and D atoms.
We remark that, in the multi-component as well as in the single-component simulations, given the low plasma density of the SOL, a significant amount of D atoms generated in the open-field line region (emitted at the limiter or created by dissociation of molecules) penetrate in the edge, where ionization takes place due to the higher plasma density and temperature. However, the presence of the D sources inside the LCFS in the multi-component simulations shifts the ionization processes, , towards the core with respect to the results with respect to single-component simulations, as shown in Fig. 4.
Figure 5: Radial profiles of the ions and neutrals species densities, averaged over the toroidal and poloidal directions, evaluated over an interval of from a quasi-steady state simulation described in Sec. 6.
To conclude, we present radial plots of the particle densities (Fig. 5) and radial fluxes (Fig. 6), obtained by evaluating the time, toroidal and poloidal average of these quantities. In Fig. 6, we discriminate the contributions of the , diamagnetic and polarization drifts to the flux of the plasma ion species, and . The results from the single-component simulations are shown in Fig. 7. The profile in Fig. 5 is similar to the one observed within the single-component plasma simulation in Fig. 7, with a large density gradient region near the LCFS and a density shoulder appearing in the far SOL. In turn, the density of is small in the whole domain and peaks in the edge, across the LCFS, where most molecules are ionized and decreases rapidly towards the core, due to the small penetration of molecules in the warmer and denser plasma of the region. On the other hand, the ions observed in the open-field line region result from charge-exchange interactions between and (see Fig. 3) and the ionization of molecules reemitted from the limiter and vessel wall.
Figure 6: Radial profiles of the radial flux for ions (top), ions (middle) and neutral species D and (bottom), averaged over the toroidal and poloidal directions, evaluated over an interval of from the quasi-steady state multi-component plasma simulation described in Sec. 6. The components of the and radial flux are discriminated.
Focusing on the neutral species, we note that peaks in the open-field line region, in contrast to the single-component plasma simulation. This is the result of the molecules dissociated into D atoms in the edge and near SOL. On the other hand, we observe that decreases monotonically from the outer wall to the core interface, since molecules are generated in the open-field line region as the result of recycling processes are lost due to dissociation and ionization processes which take place mostly in the edge and near SOL.
The dissociation of molecules also impacts , the radial flux of D, presented in Fig. 6. In contrast with the single-component plasma simulation presented in Fig. 7, points radially inwards in the edge, but reverses sign in the SOL region, a consequence of the release of D atoms because of the dissociation of molecules, particularly important close to the LCFS. In addition, the D atoms reaching the outer wall associate and are reemitted as molecules, thus contributing to the outward flux of D. The multi-component simulation shows that peaks in the edge region, while for a single-component model is maximum at the LCFS. This is due to the D atoms that are generated in the edge region close to the LCFS in a multi-component model, compensating their ionization. At the same time, we note that , the radial flux of molecules, points radially inwards in the whole domain (see Fig. 6). More precisely, is approximately constant in the SOL, because the loss of molecules due to dissociation is compensated by the molecules recycled at the limiter. Then, decreases in the edge as a consequence of the molecules being dissociated and/or ionized because of the larger temperatures and densities in this region and becomes negligible towards the core.
Figure 7: Radial profiles of density (top) and radial flux (bottom) for the and D species, averaged over the toroidal and poloidal directions, evaluated over an interval of from a quasi-steady state single-component plasma situation. The components behind the radial ion flux are discriminated. Plasma and neutral grid resolution, as well as simulation parameters, are the same considered in Fig. 4.
Turning to the dynamics of the ion species, we note that the radial flux of ions points radially outwards across the whole domain and is mostly determined by the dominant flux except near the core, where the diamagnetic flux, dominates over the flux. The polarization drift contribution is negligible in the whole domain. We also remark that the flux increases across the edge region from the core to the separatrix, having a maximum in the near SOL, and then decreases gradually across the open-field line region. This contrasts with the behavior of the ion flux in the single-component plasma simulation (see Fig. 7), where the flux peaks at the LCFS. This difference is related to the location of the ionization source . Indeed, while the source has a smooth profile and peaks at the LCFS in the single-component model, the ionization source peaks further inside the edge in the multi-component model, accounting for a sharp increase of the flux in the edge close to the LCFS.
Fig. 6 shows that the radial flux of ions points radially outwards in the SOL and radially inwards in the edge. This is a consequence of the fact that most are generated in the vicinity of the LCFS, where the molecules are ionized by the warmer and denser plasma. The radial flux is determined by the balance between the inward pointing and outward pointing diamagnetic drift components in the SOL, by the flux in the edge close to the LCFS, and by the diamagnetic component towards the core.
We also note that the inward pointing is sharply peaked in the edge, close to the LCFS. This is because most ions are generated by ionization of molecules in that region and are then dissociated after traveling a short distance. Indeed, the location of the peak of corresponds to the one of the profile in Fig. 5. The flux of associated with the polarization drift is not represented in Fig. 6 because it is neglected in our model. We note that is three to four orders of magnitude smaller than , which is a consequence of the ratio . Since the polarization drift component is expected to be small compared to the total molecular ion flux, , we conclude that neglecting the polarization drift terms in Eqs. (19-27) has indeed a negligible impact on the simulation results.
7 Conclusions
In this work we present a multi-component model for the self-consistent description of the neutral and plasma dynamics in the tokamak boundary. This model is implemented in the GBS code, allowing for the simulation of a deuterium plasma in the edge and SOL regions of a tokamak, including electrons, and ions, D atoms and molecules. The neutral and the plasma models are coupled through a number of collisional processes, which give rise to neutral-plasma interaction terms in the plasma and neutral equations. The reactions considered include ionization, electron-neutral elastic collisions, charge-exchange and dissociative processes. The multi-component plasma model relies on the Braginskii fluid equations derived in the drift limit, being an extension of the single ion species model to account for ions and closed by following Zhdanov approach. As for the neutral species, we extend the approach considered in the single neutral species model of GBS [35] to include the molecular species, . The neutrals are computed by solving two coupled kinetic equations for the D and species, which is carried out by using the method of characteristics. The resulting system of linear integral equations are then discretized and solved for the and densities.
The results from the first simulation carried out using the multi-component model are described in the sheath-limited regime in a toroidally limited plasma. The results exhibit some noticeable differences with respect to the single-ion component implemented in GBS. We observe an up-down asymmetry in the and density, which are larger below the equatorial midplane. This is related to the counterclockwise parallel flow of the plasma in the edge, observed in the profiles of , and . This feature also leads to larger recycling rates and a higher density of neutral particles in the upper side of the limiter, compared to the lower side. Moreover, the simulation shows that the density of the neutral species, , is about one order of magnitude smaller than in the open-field line region and two orders of magnitude smaller in the edge, while is about three to four orders of magnitude smaller than , even in the edge close to the LCFS, where peaks.
By taking into account the molecular dynamics, the first simulations based upon the multi-component model also shed some light on the role played by molecules on the plasma fuelling. As a matter of fact, particles are generated close to the LCFS. A large fraction of molecules reach the closed field line region, where they are most often dissociated into atomic D by the warmer and denser plasma. The resulting D atoms and the remaining molecules are then ionized inside the edge, with the ions being quickly dissociated as a consequence of the high electron densities and temperatures. The simulation results therefore show that the peak of the ionization of D atoms is shifted radially inwards with respect to the results from the single-species simulations.
The radial profiles of the densities and radial fluxes are also impacted by the presence of molecular species. We observe that the radial flux of increases sharply in the edge close to the LCFS as a result of the peak of the ionization source observed in that region. The flux of then remains high in the vicinity of the LCFS, and decreases sharply again in the near SOL, where the sources of are outweighed by the sinks at the limiter. This is a major difference with respect to the flux observed in the single-ion species simulation, which is maximum at the LCFS. On the other hand, the D density peaks in the SOL due to the ions dissociated there. This also explains why the D radial flux reverses sign, pointing radially outwards in the far SOL. On the other hand, the inward flux of D atoms in the edge increases radially inwards in the vicinity of the LCFS, since D atoms are also generated in that region as a result of dissociation of molecules.
Ultimately, our results show that the multi-component model for the self-consistent description of the neutral-plasma interaction can provide a description of a deuterium plasma that captures the main features of the molecular dynamics and its overall impact. While describing the turbulent phenomena that lead to cross-field transport, it is possible to address a multi-component plasma and more than one neutral species at a kinetic level. The procedure described here can be extended to include additional plasma and neutral species, as well as additional collisional processes.
Appendix A: Evaluation of average electron energy loss and reaction product energies in collisional processes
The Franck-Condon principle [58, 59] states that electronic excitation occurs over a timescale considerably shorter than the characteristic timescale associated with vibration or dissociation of the diatomic species. In turn, the vibration or dissociation timescales are much shorter than the electron deexcitation timescale. As a result, when an electron impacts a molecule or a ion, an electronic excitation is observed with no significant change in the inter-atomic distance (vertical transition). If the excited state is not stable, the molecule dissociates before deexcitation takes place. In this case, the difference between the excitation energy and the dissociation energy is converted into kinetic energy of the products (ionization and dissociative energies are discussed in [60]). We note that the exact energies of the products of dissociation reactions depend on the vibrational level of the molecule or ion. Considering the excitation of a molecule in a given initial state, the set of vibrational levels accessible for the molecule in the final state are the ones lying within the region of the potential energy surface accessed by that particular vertical transition, known as the Franck-Condon region. The mean energy of the reaction products is thus the average over the Franck-Condon region, taking into account all accessible vibrational states.
In the present work, we model the products of dissociative reactions by considering that they are reemitted isotropically in the reference frame of the incoming massive particle ( or ), thus approximating their velocity distribution as a Maxwellian centered at the velocity of the incoming or . The temperature of the Maxwellian, together with the average electron energy loss for each process, are obtained from the values presented in [43]. Since these energies depend on the intermediate excited state of the or particle, different values are found for different channels within the same dissociative process. This requires that an average is performed over all possible excited states, taking into account the respective cross section of each process. We present these calculations in detail for each process, following [43].
The energy loss and the energy of the reaction products may depend on the electronic levels () and sub-levels () of the reaction products, on the molecular orbital (MO) of the intermediate state, if bonding or antibonding, and on the energy of the incident electron. The energy values are experimentally determined for all relevant dissociation channels. These quantities are then averaged over all vibrational states of the molecules or ion and over the Franck-Condon region, from [43].
We start by considering the dissociation of molecules, i.e.
(99)
For this reaction, the values of the electron energy loss, , and reaction product energies, , depend significantly on the electronic state of the products. Hence, considering that there are electronic states of the reaction products and, associated, different sub-processes contributing to the dissociation of , the average electron energy loss is obtained by performing a weighed average of , the energy loss for the sub-process , based on the reaction rate, yielding
(100)
For simplicity, we evaluate all quantities at the reference temperature, . Similarly, the average value for the energy of the reaction products is obtained as
(101)
with the average energy of the products for the sub-process .
The values of , , are presented in Table 3 for all sub-processes. The additional information between brackets refers to the minimum and maximum of the range of energies accessible to and , following the values listed in [43]. We highlight that denotes a D atom in the fundamental state (electron at the lowest orbital ), while and denote an atom in the excited state with the electron in an orbital of type or , respectively, and represents an atom in the excited state . Following [43], we assume that the energy is equally distributed over the reaction products, regardless of the fact that their electronic states are the same. Based on the values in Table 3, from Eqs. (100) and (101), we obtain and , respectively, at . These are the values mentioned in Table 2.
Reaction
Table 3: product, average electron energy loss and average energy of reaction products for each sub-process of dissociation.
Focusing now on the dissociative-ionization of ,
(102)
we consider three cases. If the incoming electron has an energy , with , no dissociation takes place. If , with , the electron can ionize the molecule, resulting in an unstable ion, which then dissociates into a D atom and a ion. The short-lived has the electron in a bonding molecular orbital (MO) with -symmetry, thus exhibiting (g) symmetry (German for even) state, denoted as . If , the intermediate ion has the electron in a higher-energy antibonding MO with -symmetry, which exhibits (u) symmetry (German for odd), thus denoted as . As a result of the different energy levels of the intermediate ion, the energy of the final products will also be different, as well as the average electron energy loss. According to the results presented in [43], these energies still depend on the energy of the incoming electron within each sub-process. To simplify the evaluation of the and the energy of the products, we consider the energy to be evenly distributed by the reaction products (D and ) and we consider the two cases separately. For , all dissociative-ionization events originate an intermediate state , while for all events generate an intermediate state . The values for the electron energy loss and reaction product energies being considered for each case are evaluated for [43] and listed in Table 4. We note that this is just an approximation, as even with there are electrons with energies superior to the threshold that will generate a ion in a state, and vice-versa. Nevertheless, this approximation avoids us to evaluate and at every single value of .
Table 4: Average electron energy loss and average energy of reaction products for the two cases of dissociative-ionization of .
Reaction
For the dissociation of , i.e.
(103)
different sub-processes are taken into account, following an approach similar to the one adopted to treat the dissociation of . We perform a weighed average of the electron energy loss and the reaction products energy by using Eqs. (100) and (101), respectively. The values of , and for each sub-process are presented in Table 5. The weighed averaged values for the electron energy loss and reaction products energy at the reference temperature, , yield and , as listed in Table 5.
Reaction
Table 5: product, average electron energy loss and average energy of reaction products for each sub-process of dissociation.
Regarding the dissociative-ionization of , i.e.
(104)
we follow [43], where the average energy of the resulting ions is obtained from an average performed over all vibrational states () of the ion and over the Franck-Condon region. This yields , while the average electron energy loss is .
We finally focus on the dissociative-recombination of , which generates a D atom in the fundamental state (electron in orbital ) and a D atom in an excited state (electron with principal quantum number ), i.e.
(105)
We assume that the energy of the products is evenly distributed among the two D atoms and is given by
(106)
with the Rydberg unit of energy (corresponding to the electron binding energy in a hydrogen atom in the fundamental state). Since this expression depends on the energy of the incoming electron, , and the electronic level of the excited atom, , we assume an energy of the incident electron of , the typical value in the region around the LCFS at the HFS, and consider that these atoms are most likely in the accessible state of lowest energy (considering a higher excited state would not change the value of the energy of the products by a significant amount). Under these assumptions, we get .
Appendix B: Zhdanov collisional closure
We focus on the derivation of the parallel friction forces and the parallel heat fluxes, denoted respectively by and for a given species , with and , where we introduce the thermal component of the velocity, , with the fluid velocity of the species, and the collision operator , with describing collisions of species with species . We consider the collisional closure derived by Zhdanov in [39], relying on the approach proposed in [32] and discussed in [50] for its numerical implementation.
Following [39], the parallel component of the friction forces and heat fluxes of the species is related to the parallel gradients of the temperature and parallel velocity of all species through
(107)
where denotes the temperature of plasma species and is the parallel component of the fluid velocity of species with respect to the center of mass of the plasma, , with . The matrix relates the parallel heat fluxes and friction forces with the parallel gradients of temperature and parallel velocity. We remark that Eq. (107) simplifies the general result obtained by Zhdanov [39] to the case of singly-ionized states, neglecting possible multiplicity of charge states for the chemical species present in the plasma.
In order to compute the matrix , we consider the -moment approximation of the distribution function [39], thus including the moments up to the fifth order moment. We first express and in terms of these moments of the distribution function, namely the first order moment, , the third order moment, , and the fifth order moment, , where we introduce the velocity with respect to the center of mass of the plasma, , and the parameter , with . Since only the expressions for the parallel component of the friction forces and heat fluxes are needed, we consider only the parallel component of these equations. The heat flux, , simply corresponds to the third order moment, , while the friction forces, , are obtained in terms of , and [61], yielding
(108)
(109)
where and are respectively the mass and density of species , is the reduced mass, and are polynomial functions of the local plasma density and temperature, their exact expressions being presented in [39] (chapter 8.1, pp. 163-164). Eqs. (108) and (109) can then be written in matrix form as
(110)
where the matrices and are defined to satisfy Eqs. (108) and (109). We now aim at expressing the moments and in terms of and . This can be achieved by solving a system of moment equations similar to the one presented in [39] (chapter 8.1, pp. 162-163), including the time evolution of the moments (, and ) and the time evolution of basic thermodynamic variables (, and ). We neglect time derivatives and nonlinear terms. For simplicity, we also assume that, for two massive particle species and , the condition is fulfilled, which allows us to write . Moreover, as long as is verified, can also be replaced by , following [39] (the simulation results shown in Fig. 2 meet these conditions). We therefore impose , while no assumption is made on the temperature and pressure gradients, i.e. temperature gradients can be different from species to species [39].
The parallel projection of the system of moment equations can then be written as (see [61])
(111)
(112)
where is the pressure of species . Rewriting Eqs. (111-112) in matrix form, one obtains
(113)
which can be inverted to express the parallel third and fourth order fluid moments in terms of the parallel gradient of temperature and relative parallel velocity as
(114)
Finally, making use of Eq. (114) to express and in Eq. (110) in terms of the parallel temperature gradients and relative velocities, one obtains the expressions for the parallel heat flux and friction forces in the matrix form presented in Eq. (107), that is
(115)
Since the matrices , , and are fully determined by Eqs. (108), (109), (111) and (112), the expressions of the parallel heat flux and friction forces can be found. Following Zhdanov [39], these matrices can be expressed in terms of the local values of plasma quantities, namely densities , and and temperatures and (we again assume , mass ratios and characteristic time scales and , with defined as the inverse of the collision frequency for momentum transfer between electrons and ions, and the ion timescale defined as the inverse of the collision frequency for momentum transfer between ions. We retain only terms of leading order in , while terms proportional to the fast electron timescale are neglected when compared to terms proportional to , which considerably simplifies the final expressions. We also highlight that, besides imposing the quasi-neutrality relation , we take into account the fact that the density of the molecular ion species is much smaller than the density of the main ion species for typical tokamak boundary conditions, i.e. , keeping therefore only leading order terms in . As a result, the friction forces between molecular ions and other species are neglected, as well as molecular ion temperature gradient terms, while friction and thermal force contributions involving and species are kept in the expressions of the parallel components of the heat fluxes and friction forces. The expressions obtained for the friction forces and heat fluxes finally yield
(116)
The expressions in Eqs. (116) can be simplified by applying the relation between the electron and ion characteristic times,
(117)
having again assumed . This enables one to write appearing in Eq. (116) in terms of . Following Braginskii’s approach [62] and considering that the the electron characteristic time is , we then write Eqs. (116) in terms of the resistivity, defined as [22, 23]
(118)
The parallel friction forces and heat fluxes, as they appear in Eqs. (22-24) and Eqs. (25-27), respectively, are therefore written in normalized units as
(119)
We note that, similarly to the single-ion species model implemented in GBS [22], the ohmic heating terms are neglected.
Appendix C: List of kernel functions
The kernels used in Eqs. (80-83) for , , and are defined as
(120)
(121)
(122)
(123)
(124)
(125)
(126)
(127)
(128)
(129)
(130)
(131)
(132)
(133)
(134)
(135)
(136)
(137)
(138)
(139)
(140)
(141)
where the kernel functions for a given are defined as
(142)
(143)
(144)
(145)
(146)
(147)
(148)
(149)
(150)
(151)
(152)
(153)
(154)
(155)
(156)
(157)
(158)
(159)
(160)
(161)
(162)
(163)
We remark that all velocity distributions given by a Maxwellian or a Knudsen cosine law are integrated along the parallel velocity, that is
(164)
(165)
(166)
(167)
(168)
(169)
(170)
(171)
(172)
(173)
Appendix D: Numerical solution of the neutral equations
The coupled neutral equations for and D, Eqs. (80-83), may be discretized as a linear matrix system, , with the unknown representing the density and boundary flux of the and D species. Indicating with the number of points that discretize the poloidal plane and the number of points discretizing the boundary, is a vector of size , is a matrix and is a vector that includes all contributions not proportional to the neutral density or flux, namely the effect of recombination of and with electrons, the effect of dissociative processes to which ions are subject and the contributions from the flux of and ions to the boundary.
The matrix , and the vectors and can then be written as
(174)
where is a matrix of size ,
(175)
that discretizes the kernel defined in Eq. (84) at the spatial points where is evaluated. The matrix
(176)
has size and discretizes the kernel defined in Eq. (134) at the points where is evaluated. The other matrices appearing in the definition of are defined similarly,
(177)
(178)
(179)
(180)
(181)
(182)
(183)
(184)
(185)
(186)
(187)
(188)
(189)
(190)
The vector is defined through the vectors and of size ,
(191)
(192)
and the vector and of size ,
(193)
(194)
It is remarked that the vector can also be written as , where refers to the densities and boundary fluxes of the and ion species,
(195)
and the matrix can be expressed as
(196)
with entries
(197)
(198)
(199)
(200)
(201)
(202)
(203)
(204)
(205)
(206)
(207)
(208)
(209)
(210)
(211)
(212)
We remark that a convergence study to estimate the error introduced by the discretization of the neutral equation was carried out for a single neutral species model and it is reported in [56].
ACKNOWLEDGEMENTS
We would like to thank C. Theiler, C. Wersal, D. Mancini, K. Verhaegh and M. Wensing for useful discussions. The simulations presented herein were carried out in part at CSCS (Swiss National Supercomputing Center) under the project ID s882 and in part on the CINECA Marconi supercomputer under the GBSedge project. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Fond National Suisse de la Recherche Scientifique and from the Euratom research and training programme 2014-2018 and 2019-2020 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding authors upon reasonable request.
REFERENCES
References
[1]
Peter C Stangeby.
The Plasma Boundary of Magnetic Fusion Devices.
IOP, Bristol, 2000.
[2]
Charles Spencer Pitcher.
Tokamak Plasma Interaction with Limiters.
PhD thesis, University of Toronto, 1988.
[3]
Chen Zhang, Chaofeng Sang, Liang Wang, Mingyu Chang, Daoyuan Liu, and Dezhen
Wang.
Effect of carbon and tungsten plasma-facing materials on the
divertor and pedestal plasma in EAST.
Plasma Physics and Controlled Fusion, 61:115013, 2019.
[4]
Daren Stotler and Charles Karney.
Neutral Gas Transport Modeling with DEGAS 2.
Contributions to Plasma Physics, 34(2-3):392–397, 1994.
[5]
B. J. Braams.
Radiative Divertor Modelling for ITER and TPX.
Contributions to Plasma Physics, 36:276, 1996.
[6]
M. Baelmans, D. Reiter, and R. R. Weynants.
New Developments in Plasma Edge Modeling with Particular Emphasis on
Drift Flows and Electric Fields.
Contributions to Plasma Physics, 36:117, 1996.
[7]
R. Simonini, G. Corrigan, G. Radford, J. Spence, and A. Taroni.
Models and Numerics in the Multi-Fluid 2-D Edge Plasma Code
EDGE2D/U.
Contributions to Plasma Physics, 34:368, 1994.
[8]
Y. Feng, F. Sardei, and J. Kisslinger.
3D fluid modelling of the edge plasma by means of a Monte Carlo
technique.
Journal of Nuclear Materials, 266:812–818, 1999.
[9]
H. Bufferand, G. Ciraolo, L. Isoardi, G. Chiavassa, F. Schwander, E. Serre,
N. Fedorczak, Ph Ghendrih, and P. Tamain.
Applications of SOLEDGE-2D code to complex SOL configurations and
analysis of Mach probe measurements.
Journal of Nuclear Materials, 415:S589–S592, 2011.
[10]
R. Zagórski, H. Gerhauser, and H. A. Claaßen.
Numerical simulation of the TEXTOR edge plasma including drifts and
impurities.
Contributions to Plasma Physics, 38:61, 1998.
[11]
Keiji Sawada and Takashi Fujimoto.
Effective ionization and dissociation rate coefficients of molecular
hydrogen in plasma.
Journal of Applied Physics, 78:2913, 1995.
[12]
A. Yu Pigarov, S. I. Krasheninnikov, W. P. West, T. D. Rognlien, J. A. Boedo,
D. G. Whyte, C. J. Lasnier, T. W. Petrie, M. J. Schaffer, and J. G. Watkins.
DIII-D edge plasma simulations with UEDGE code including
non-diffusive anomalous cross-field transport.
Journal of Nuclear Materials, 313-316(SUPPL.):1076–1080, 2003.
[13]
D. Reiter, M. Baelmans, and P. Börner.
The EIRENE and B2-EIRENE codes.
Fusion Science and Technology, 47(2):172–186, 2005.
[14]
John Mandrekas.
GTNEUT: A code for the calculation of neutral particle transport in
plasmas based on the transmission and escape probability method.
Computer Physics Communications, 161(1-2):36–64, 2004.
[15]
K. Shimizu, T. Takizuka, S. Sakurai, H. Tamai, H. Takenaga, H. Kubo, and
Y. Miura.
Simulation of divertor detachment characteristics in JT-60 with
superconducting coils.
Journal of Nuclear Materials, 313-316(SUPPL.):1277–1281, 2003.
[16]
K. D. Lawson, M. Groth, D. Harting, S. Menmuir, D. Reiter, S. Brezinsek,
G. Corrigan, P. Drewelow, C. F. Maggi, A. G. Meigs, J. Simpson, M. F. Stamp,
S. Wiesen, and Contributors.
A study of the atomic and molecular power loss terms in
EDGE2D-EIRENE simulations of JET ITER-like wall L-mode discharges.
Nuclear Materials and Energy, 12:924–930, 2017.
[17]
J. D. Lore, J. M. Canik, Y. Feng, J. W. Ahn, R. Maingi, and V. Soukhanovskii.
Implementation of the 3D edge plasma code EMC3-EIRENE on NSTX.
Nuclear Fusion, 52(5):054012, 2012.
[18]
R. Schneider, X. Bonnin, K. Borrass, D. P. Coster, H. Kastelewicz, D. Reiter,
V. A. Rozhansky, and B. J. Braams.
Plasma edge physics with B2-Eirene.
Contributions to Plasma Physics, 46(1-2):3–191, 2006.
[19]
B. D. Dudson, M. V. Umansky, X. Q. Xu, P. B. Snyder, and H. R. Wilson.
BOUT++: A framework for parallel plasma fluid simulations.
Computer Physics Communications, 180(9):1467–1480, 2009.
[20]
B. D. Dudson and J. Leddy.
Hermes: Global plasma edge fluid turbulence simulations.
Plasma Physics and Controlled Fusion, 59(5), 2017.
[21]
Matthias Wiesenberger, Lukas Einkemmer, Markus Held, Albert Gutierrez-Milla,
Xavier Sáez, and Roman Iakymchuk.
Reproducibility, accuracy and performance of the FELTOR code and
library on parallel computer architectures.
Computer Physics Communications, 238:145–156, 2019.
[22]
P. Ricci, F. D. Halpern, S. Jolliet, J. Loizu, A. Mosetto, A. Fasoli, I. Furno,
and C. Theiler.
Simulation of plasma turbulence in scrape-off layer conditions: The
GBS code, simulation results and code validation.
Plasma Physics and Controlled Fusion, 54(12):124047, 2012.
[23]
F. D. Halpern, P. Ricci, S. Jolliet, J. Loizu, J. Morales, A. Mosetto,
F. Musil, F. Riva, T. M. Tran, and C. Wersal.
The GBS code for tokamak scrape-off layer simulations.
Journal of Computational Physics, 315:388–408, 2016.
[24]
M. Francisquez, B. Zhu, and B. N. Rogers.
Global 3D Braginskii simulations of the tokamak edge region of IWL
discharges.
Nuclear Fusion, 57(11), 2017.
[25]
Andreas Stegmeir, David Coster, Alexander Ross, Omar Maj, Karl Lackner, and
Emanuele Poli.
GRILLIX: A 3D turbulence code based on the flux-coordinate
independent approach.
Plasma Physics and Controlled Fusion, 60(3), 2018.
[26]
A. S. Thrysøe, M. Løiten, J. Madsen, V. Naulin, A. H. Nielsen, and
J. Juul Rasmussen.
Plasma particle sources due to interactions with neutrals in a
turbulent scrape-off layer of a toroidally confined plasma.
Physics of Plasmas, 25(3), 2018.
[27]
Camille Baudoin, Patrick Tamain, Hugo Bufferand, Guido Ciraolo, Nicolas
Fedorczak, Davide Galassi, Philippe Ghendrih, and Nicolas Nace.
Turbulent heat transport in TOKAM3X edge plasma simulations.
Contributions to Plasma Physics, 58:484–489, 2018.
[28]
E. L. Shi, G. W. Hammett, T. Stoltzfus-Dueck, and A. Hakim.
Gyrokinetic continuum simulation of turbulence in a straight
open-field-line plasma.
Journal of Plasma Physics, 83(3), 2017.
[29]
C. S. Chang and S. Ku.
Spontaneous rotation sources in a quiescent tokamak edge plasma.
Physics of Plasmas, 15(6), 2008.
[30]
S. Ku, C. S. Chang, and P. H. Diamond.
Full-f gyrokinetic particle simulation of centrally heated global
ITG turbulence from magnetic axis to edge pedestal top in a realistic tokamak
geometry.
Nuclear Fusion, 49(11), 2009.
[31]
H. Bufferand, C. Baudoin, J. Bucalossi, G. Ciraolo, J. Denis, N. Fedorczak,
D. Galassi, Ph Ghendrih, R. Leybros, Y. Marandet, N. Mellet, J. Morales,
N. Nace, E. Serre, P. Tamain, and M. Valentinuzzi.
Implementation of drift velocities and currents in
SOLEDGE2D–EIRENE.
Nuclear Materials and Energy, 12:852–857, 2017.
[32]
H. Bufferand, P. Tamain, S. Baschetti, J. Bucalossi, G. Ciraolo, N. Fedorczak,
Ph Ghendrih, F. Nespoli, F. Schwander, E. Serre, and Y. Marandet.
Three-dimensional modelling of edge multi-component plasma taking
into account realistic wall geometry.
Nuclear Materials and Energy, 18:82–86, 2019.
[33]
Alexander S. Thrysøe, Laust E.H. Tophøj, Volker Naulin, Jens Juul
Rasmussen, Jens Madsen, and Anders H. Nielsen.
The influence of blobs on neutral particles in the scrape-off
layer.
Plasma Physics and Controlled Fusion, 58(4), 2016.
[34]
A. S. Thrysøe, J. Madsen, V. Naulin, and J. Juul Rasmussen.
Influence of molecular dissociation on blob-induced atom density
perturbations.
Nuclear Fusion, 58(9), 2018.
[35]
C. Wersal and P. Ricci.
A first-principles self-consistent model of plasma turbulence and
kinetic neutral dynamics in the tokamak scrape-off layer.
Nuclear Fusion, 55(12):123014, 2015.
[36]
C. Wersal, P. Ricci, and J. Loizu.
A comparison between a refined two-point model for the limited
tokamak SOL and self-consistent plasma turbulence simulations.
Plasma Physics and Controlled Fusion, 59(4):044011, 2017.
[37]
C. Wersal and P. Ricci.
Impact of neutral density fluctuations on gas puff imaging
diagnostics.
Nuclear Fusion, 57(11):116018, 2017.
[38]
A. Coroado and P. Ricci.
Moving toward mass-conserving simulations of plasma turbulence and
kinetic neutrals in the tokamak boundary with the GBS code.
Physics of Plasmas, 28(2):022310, 2021.
[39]
V M Zhdanov.
Transport Processes in Multicomponent Plasma.
Plasma Physics and Controlled Fusion, 2002.
[40]
D. P. Stotler, C. H. Skinner, R. V. Budny, A. T. Ramsey, D. N. Ruzic, and R. B.
Turkot.
Modeling of neutral hydrogen velocities in the tokamak fusion test
reactor.
Physics of Plasmas, 3:4084, 1996.
[41]
M. Wensing, B. P. Duval, O. Février, A. Fil, D. Galassi, E. Havlickova,
A. Perek, H. Reimerdes, C. Theiler, K. Verhaegh, and M. Wischmeier.
SOLPS-ITER simulations of the TCV divertor upgrade.
Plasma Physics and Controlled Fusion, 61:085029, 2019.
[42]
D Reiter.
The data file AMJUEL: Additional Atomic and Molecular Data for
EIRENE.
Technical report, FZ, Forschungszentrum Julich GmbH FRG, 2011.
[43]
R.K. Janev and Langer W.D.
The Janev-Langer Hydrogen-Helium database, 1987 Contents.
Springer, 1987.
[44]
R. K. Janev.
Atomic and Molecular Processes in Fusion Edge Plasmas.
Plenum Press, New York, 1995.
[45]
Jung Sik Yoon, Mi Young Song, Jeong Min Han, Sung Ha Hwang, Won Seok Chang,
Bongju Lee, and Yukikazu Itikawa.
Cross sections for electron collisions with hydrogen molecules.
Journal of Physical and Chemical Reference Data, 37:913, 2008.
[46]
H. K. Chung.
Data for Atomic Processes of Neutral Beams in Fusion Plasma, Summary
Report of the First Research Coordination Meeting.
Technical report, IAEA Nuclear Data Section, 2017.
[47]
Predrag S. Krstić.
Inelastic processes from vibrationally excited states in slow
H^+ + H_2 and H+ H_2^+ collisions: Excitations and charge
transfer.
Physical Review A - Atomic, Molecular, and Optical Physics,
66(4):042717, 2002.
[48]
Predrag S. Krstić.
Vibrationally resolved collisions in cold hydrogen plasma.
In Nuclear Instruments and Methods in Physics Research, Section
B: Beam Interactions with Materials and Atoms, pages 58–62, 2005.
[49]
A. Zeiler, J. F. Drake, and B. Rogers.
Nonlinear reduced Braginskii equations with ion thermal dynamics in
toroidal plasma.
Physics of Plasmas, 4(6):2134–2138, 1997.
[50]
Madhusudan Raghunathan, Yannick Marandet, Hugo Bufferand, Guido Ciraolo,
Philippe Ghendrih, Patrick Tamain, and Eric Serre.
Generalized Collisional Fluid Theory for Multi-Component,
Multi-Temperature Plasma Using The Linearized Boltzmann Collision Operator
for Scrape-Off Layer/Edge Applications.
Plasma Physics and Controlled Fusion, 29:1–29, 2021.
[51]
J. Loizu, P. Ricci, F. D. Halpern, and S. Jolliet.
Boundary conditions for plasma fluid models at the magnetic
presheath entrance.
Physics of Plasmas, 19(12):122307, 2012.
[52]
Annamaria Mosetto, Federico D. Halpern, Sébastien Jolliet, Joaquim Loizu,
and Paolo Ricci.
Finite ion temperature effects on scrape-off layer turbulence.
Physics of Plasmas, 22(1), 2015.
[53]
D. Tskhakaya, S. Kuhn, and Y. Tomita.
Formulation of boundary conditions for the unmagnetized
multi-ion-component plasma sheath.
Contributions to Plasma Physics, 46(7-9):649–654, 2006.
[54]
Paola Paruta, P. Ricci, F. Riva, C. Wersal, C. Beadle, and B. Frei.
Simulation of plasma turbulence in the periphery of diverted tokamak
by using the GBS code.
Physics of Plasmas, 25(11):112301, 2018.
[55]
D. Galassi, H. Reimerdes, C. Theiler, M. Wensing, H. Bufferand, G. Ciraolo,
P. Innocente, Y. Marandet, and P. Tamain.
Numerical investigation of optimal divertor gas baffle closure on
TCV.
Plasma Physics and Controlled Fusion, 2020.
[56]
M Giacomin, P Ricci, A Coroado, G Fourestey, D Galassi, E Lanti, and D Mancini.
The GBS code for the self-consistent simulation of plasma turbulence
and kinetic neutral dynamics in the tokamak boundary.
submitted, 2021.
[57]
J. Loizu, P. Ricci, F. D. Halpern, S. Jolliet, and A. Mosetto.
On the electrostatic potential in the scrape-off layer of magnetic
confinement devices.
Plasma Physics and Controlled Fusion, 55(12):124019, 2013.
[58]
J. Franck.
Elementary processes of photochemical reactions.
Transactions of the Faraday Society, 21:536–542, 1926.
[59]
Edward U. Condon.
Nuclear motions associated with electron transitions in diatomic
molecules.
Physical Review, 32(6):858–872, 1928.
[60]
Jinjun Liu, Edcel J. Salumbides, Urs Hollenstein, Jeroen C.J. Koelemeij,
Kjeld S.E. Eikema, Wim Ubachs, and Fŕdŕic Merkt.
Determination of the ionization and dissociation energies of the
hydrogen molecule.
Journal of Chemical Physics, 130:174306, 2009.
[61]
V. M. Zhdanov and P. N. Yushmanov.
Diffusion and heat transfer in a multicomponent completely ionized
plasma.
Journal of Applied Mechanics and Technical Physics, 1980.
[62]
S. I. Braginskii.
Transport processes in a plasma.
In Reviews of Plasma Physics, volume 1. Consultants Bureau, New
York, 1965.
