Journal Pre-proof

First integrated core-edge fluid simulation of ITER's Limiter-Divertor transition with SolEdge-HDG

M. Scotto d'Abusco, I. Kudashev, G. Giorgiani,
A. Glasser-Medvedeva, F. Schwander, E. Serre, J. Bucalossi,
H. Bufferand, G. Ciraolo, P. Tamain

<figure><img src="image_1.png" /></figure>

PII: S2352-1791(24)00173-X

DOI: https://doi.org/10.1016/j.nme.2024.101750

Reference: NME 101750

To appear in: *Nuclear Materials and Energy*

Received date: 7 June 2024

Revised date: 24 September 2024

Accepted date: 28 September 2024

Please cite this article as: M.S. d'Abusco, I. Kudashev, G. Giorgiani et al., First integrated core-edge fluid simulation of ITER's Limiter-Divertor transition with SolEdge-HDG, *Nuclear Materials and Energy* (2024), doi: https://doi.org/10.1016/j.nme.2024.101750.

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2024 Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
# First integrated core-edge fluid simulation of ITER's Limiter-Divertor transition with SolEdge-HDG

M. Scotto d'Abusco<sup>a</sup>, I. Kudashev<sup>b</sup>, G. Giorgiani<sup>b</sup>, A. Glasser-Medvedeva<sup>b</sup>, F. Schwander<sup>b</sup>, E. Serre<sup>b</sup>, J. Bucalossi<sup>c</sup>, H. Bufferand<sup>c</sup>, G. Ciraolo<sup>c</sup>, P. Tamain<sup>c</sup>

## Abstract

This work explores the Limiter-Divertor transition (L-D) during the current ramp-up of ITER's Q=10 baseline plasma scenario at various central line-integrated density $n_{li}$ values. The analysis, based on transport simulations performed with the latest version of SolEdge-HDG, focuses on the time evolution of heat and ion particle fluxes, revealing regions of elevated temperature on the inner wall and plasma-facing components (PFCs) despite moderate loads. The investigation also delves into the effects of perpendicular convection flux terms on density build-up, comparing different formulations and their interplay with auxiliary heating sources. Furthermore, the paper shows the impact of taking into account the evolution of the parallel neutral momentum on plasma and neutral density at the targets in the context of an ITER steady-state scenario.

**Keywords:** ITER, SolEdge-HDG, Hybridizable Discontinuous Galerkin, ramp-up simulation, integrated core-edge simulation

## 1. Introduction

Achieving robust, hot and dense core plasma while maintaining the target power loads below the operational limits in ITER high fusion performance scenarios will require a well-controlled core-edge synergy Editors et al. (1999). Although several works have focused on fully integrated simulation of ITER scenarios from core to divertor Romanelli et al. (2015); Garzotti et al. (2018), the description of the evolution of the discharge is rudimentary, so that large uncertainties in plasma transport and boundary evolution remain. With the magnetic evolution of the plasma, self-consistent core-edge plasma studies are particularly important to characterize the plasma transport past to the X-point formation prior entering the burning phase, to guarantee the existence of a safe operation scenario. More than the other, the hot and dense ITER edge plasma becomes opaque to hydrogenic neutrals penetration at very low fraction of the Greenwald density limit $n_{Gw}$. Kukushkin et al. (2009), making core fueling with gas puffing almost ineffective. Establishing the correct plasma density build-up throughout the linear ramp-up phase of the plasma current is hence crucial for maintaining the optimal interplay between core fuelling and heating on one hand, with heat loads and radiated power at the scrape of layer (SOL) on the other. One of the tools used for scenario modelling of the whole plasma domain of ITER is the high-fidelity code JINTRAC ROMANELLI et al. (2014), which couples a suite of codes for core and edge plasma. Simulations of the full discharge are performed using 1D core transport codes, assuming averaging over plasma surfaces, coupled with SOL codes, like SOLPS-ITER Wiesen et al. (2015) and EDGE2D Simonini et al. (1994). Although the detailed physics description, this approach

requires the mapping of the solution at the interface to achieve a good match between the two distinct regions. Moreover, code likes SOLPS-ITER and EDGE2D relies on flux aligned mesh which prevent the simulation of non-stationary magnetic configurations. To overcome these limitations, in the present work we apply the new release of SolEdge-HDG Giorgiani et al. (2018, 2020) code to simulate the L-D transition of ITER Q = 10 baseline scenario of a deuterium plasma for different values of the central line integrated density $n_{li}$. SolEdge-HDG is a finite element transport code based on Hybrid Discontinuous Galerkin (HDG) scheme, and employs high-order unstructured meshes to model the whole plasma domain, from the core center down to the PFCs (see Fig. 1), allowing for simulation of non-stationary magnetic equilibrium in realistic tokamak geometry d'Abusco et al. (2022). Therefore the code has been employed to investigate the time evolution of ion particle and heat fluxes and temperatures for different values of $n_{li}$. The model is then enriched with a convective flux terms to simulate anomalous or neoclassical pinch Kudashev et al. (2024a). The influence of this term on the transport of the recycled plasma at the wall and core fueling is shown in three test cases. The plasma recycling at the wall is modelled through neutral transport coupled to the plasma. Neutrals transport plays a crucial role for the correct assessment of volumetric losses of momentum and energy at the divertor target. The neutral dynamic is usually addressed through code coupling with Monte-Carlo codes, which solves the kinetic Boltzmann equation, such as EIRENE D. Reiter and Börner (2005). While providing a rich description of the physics, this method is usually computationally expensive. In order to keep the calculation fast, fluid models are here considered. The transport of the neutral density is assumed to be purely diffusive,

<footnote><sup>a</sup>Princeton Plasma Physics Laboratory, 100 Stellarator Rd, Princeton, 08540, NJ, US</footnote>

<footnote><sup>b</sup>Aix-Marseille Université, CNRS, Centrale Marseille, M2P2 UMR 7340, Marseille, France,</footnote>

<footnote><sup>c</sup>IRFM, CEA Cadarache, St. Paul-lez-Durance, France,</footnote>
with a non constant diffusion coefficient function of the neutral mean free path as inspired by Horsten et al. (2017). In this work, the model can be enriched by taking into account the evolution of the neutral parallel momentum.

The paper is organized as follows. Section 2 introduces the conservation equations together with the boundary conditions implemented in SolEdge-HDG. The simulations of the ITER L-D transition at different values of $n_{li}$ are analyzed in section 3 with and without convective flux term. The impact of the different neutral fluid models is analyzed in section 4. Finally, concluding remarks and perspectives are summarized in section 5.

## 2. The SolEdge-HDG transport model

### 2.1. The physical model

A set of conservation equations for the mean plasma quantities, i.e. the density, the parallel momentum and the ion and electron energy, known as Braginskii equations, is resolved using an Hybridizable Discontinuous Galerkin method based on high-order finite-elements. The system of equations writes as:

$$ \partial_t n + \nabla \cdot (nu\mathbf{b}) - \nabla \cdot (D\nabla_{\perp}n - v_p n\mathbf{b}_{\perp}) = S_n \quad (1) $$

$$ \partial_t(m_i n u) + \nabla \cdot (m_i n u^2 \mathbf{b}) + \nabla_{\parallel}(k_b n(T_e + T_i)) - \nabla \cdot (\mu \nabla_{\perp}(m_i n u)) = S_{\Gamma} \quad (2) $$

$$ \begin{aligned} & \partial_t \left( \frac{3}{2} k_b n T_i + \frac{1}{2} m_i n u^2 \right) + \nabla \cdot \left( \left( \frac{5}{2} k_b n T_i + \frac{1}{2} m_i n u^2 \right) u \mathbf{b} \right) - \nu e E_{\parallel} - \\ & \nabla \cdot \left( \frac{3}{2} k_b (T_i D \nabla_{\perp} n + n \chi_i \nabla_{\perp} T_i) \right) - \nabla \cdot \left( \frac{1}{2} m_i u^2 D \nabla_{\perp} n + \frac{1}{2} m_i \mu n \nabla_{\perp} u^2 \right) \\ & - \nabla \cdot (k_{\parallel} T_i^5 \nabla_{\parallel} T_i \mathbf{b}) + \frac{3}{2} \frac{k_b n}{\tau_{ie}} (T_e - T_i) = S_{E_i} \end{aligned} \quad (3) $$

$$ \begin{aligned} & \partial_t \left( \frac{3}{2} k_b n T_e \right) + \nabla \cdot \left( \frac{5}{2} k_b n T_e u \mathbf{b} \right) + \nu e E_{\parallel} - \nabla \cdot \left( \frac{3}{2} k_b (T_e D \nabla_{\perp} n + n \chi_e \nabla_{\perp} T_e) \right) - \nabla \cdot (k_{\parallel} T_e^5 \nabla_{\parallel} T_e \mathbf{b}) - \frac{3}{2} \frac{k_b n}{\tau_{ie}} (T_e - T_i) = S_{E_e} \end{aligned} \quad (4) $$

with $\mathbf{b}$ being the unitary vector in the direction parallel to the magnetic field and $\mathbf{b}_{\perp}$ representing the radial direction in the surface perpendicular to $\mathbf{B}$. The governing equations are projected to the magnetic field direction using the differential operator $\nabla_{\parallel} = \mathbf{b} \cdot \nabla, \nabla_{\perp} = -\mathbf{b} \cdot \nabla$. The turbulent transport in the cross field direction is assumed to be diffusive with $D, \mu, \chi_i, \chi_e$ ad-hoc coefficients. $v_p$ is an anomalous convective velocity that will be discussed in a later section. $k_{\parallel i}$ and $k_{\parallel e}$ are the parallel heat conductivity coefficients equal to $k_{\parallel i} = 60[Wm^{-1}V^{-7/2}]$ for the deuterium ion and $k_{\parallel e} = 2000[Wm^{-1}V^{-7/2}]$ for the electrons. $\tau_{ie}$ is the relaxation time for the collisional coupling term between electrons and ions and $S_n, S_{\Gamma}, S_{E_i}, S_{E_e}$ are source terms which model the ionization, recombination and charge-exchange processes caused by neutral-plasma interaction and auxiliary sources (NBI, ECRH, ICRH...) of plasma particle and energy, respectively. An Ohmic heating source term is also calculated, and included in $S_{E_e}$, employing a Spitzer-Harm formulation of plasma resistivity based on density current distribution obtained

from experimental data or equilibrium stability calculations. Thus the source of heat is not prescribed value at the core boundary as in many other SOL codes. The basic fluid neutral model uses a simple diffusion equation to evaluate the neutral density Horsten et al. (2017); Kudashev et al. (2024b) which writes as:

$$ \partial_t n_n - \nabla \cdot (D_{n_n} \nabla n_n) = S_{n_n} \quad (5) $$

with $D_{n_n} = \frac{eT_i[eV]}{m_i n(\sigma v_{cx}) + (\sigma v_{c\epsilon})}$ is a non constant diffusion coefficient. $T_n = T_i$ is assumed, and $\langle \sigma v_{cx} \rangle, \langle \sigma v_{c\epsilon} \rangle$ are the ionization and charge-exchange rate coefficient based on ADAS Summers (2004) and AMUJEL Reiter (2000) database, respectively. The extension of this model to include the description of neutral momentum will be presented in Sec. 4. The current physical model implemented in SolEdge-HDG doesn't include drift effects and impurity transport.

### 2.2. Boundary conditions

The plasma wall interactions are modelled through appropriate boundary conditions coupled to the former conservation equations. The plasma density is left free, while an outgoing super sonic velocity $|ub \cdot n| \ge c_s$ is imposed at the wall (Bohm boundary condition) with $c_s = \sqrt{k_b(T_i + T_e)/m_i}$ the plasma sound speed. The total energy fluxes are imposed to the sheath transmission values which leads to:

$$ \begin{aligned} & (nE_i + p_i)u\mathbf{b} - q_{\perp,i} - k_{\parallel i} T_i^5 \nabla_{\parallel} T_i \mathbf{b} = \left( \gamma_i \text{unk}_b T_i + \frac{1}{2} n m_i u^3 \right) \mathbf{b} \\ & (nE_e + p_e)u\mathbf{b} - q_{\perp,e} - k_{\parallel e} T_e^5 \nabla_{\parallel} T_e \mathbf{b} = \gamma_e \text{unk}_b T_e \mathbf{b} \end{aligned} \quad (6) $$

where $E_i = 3/2k_b T_i + 1/2m_i u^2$, $E_e = 3/2k_b T_e$ are the ion and electron energies in $[m^2s^{-1}]$, $q_{\perp,i}, q_{\perp,e}$ are the ion and electron perpendicular heat fluxes defined by the cross field term in Eqs. 3, 4, and $\gamma_i = 2.5, \gamma_e = 4.5$ are the ion and electron sheath transmission factors, respectively. The definition of $\gamma_i$ used in this work differs from the one introduced in Eq (2.92) in $?$. The plasma recycling into neutrals is described by imposing the net neutral flux as:

$$ -D_{n_n} \nabla n_n \cdot \mathbf{n} = [-R(-D\nabla_{\perp} n + nu\mathbf{b}) - \Gamma_{\text{puff}} + \Gamma_{\text{pump}}] \cdot \mathbf{n} \quad (7) $$

where $R$ is the recycling coefficient which states the percentage of the plasma particle flux hitting the wall that is recycled and re-injected as neutrals. $\Gamma_{\text{puff}}, \Gamma_{\text{pump}}$ are two ad hoc boundary fluxes which mimic the gas puff injection and the neutrals extraction rate in the private flux region (PFR).

## 3. Simulation of the ITER L-D transition

The time evolution of the plasma quantities during the L-D transition and primary ramp-up phase is investigated at different fractions of the Greenwald density $n_{Gw}$. Three different scenarios with a central line integrated density $n_{li}$ respectively equal to 10%, 20%, 25% of $n_{Gw}$ are analysed with primary attention on the effects on ion particle flux, heat flux and temperatures in different regions of interest of the wall boundary defined as in Fig. 2a. The possibility to perform L-D transition at larger
<figure><img src="image_2.png" /><figcaption>Figure 1: Snapshot of the ITER plasma ion density and electron temperature at $t = 75$ s of ITER's Q = 10 scenario as computed by SolEdge-HDG.</figcaption></figure>

values of $n_{li}$ is then explored making assumption on the shape of an inward particle pinch term, and of the interplay with auxiliary power sources. Purely diffusive neutrals will be assumed through all the scan.

### 3.1. Numerical setup

The realistic geometry of the ITER cross-section as well as an example of the unstructured mesh are shown in Fig. 2a. The tokamak PFCs are labeled from the pump region, and moving counterclockwise to the outer target (OT), outer wall (OW), top wall (TW), inner wall (IW) and inner target (IT). To carry out the

<figure><img src="image_3.png" /><figcaption>Figure 2: (a) Positions of the separatrix at different instants together with the mesh in the ITER cross-section. The main regions of interest plus pump and puff location are also labeled. (b) Temporal evolution of the perpendicular transport coefficients (red curves) and plasma current $I_p$ (blue curve).</figcaption></figure>

simulation of the ITER L-D transition, the plasma equilibrium coming from an equilibrium stability calculation (IMAS shot #135011) performed at different time steps, are imported and interpolated on the mesh. At $t = 0$ s the code is run in steady mode using the first magnetic equilibrium and current density distribution. As initial condition we seek for a steady state solution where $n_{li} = 100\% n_{Gw}$. At each following time step, the new magnetic equilibrium as well as the current distribution to determine the Ohmic power are uploaded, and the code is run over a time $dt$ that is equal to the time interval in between

two different magnetic configurations, using the transport solution obtained at the previous step as initial condition. Here $dt = 0.05$ s is chosen. The Ohmic power deposition varies with plasma resistivity modification as $Z_{eff}/T^{3/2}$ where $Z_{eff} = 1.6$ is assumed. In absence of experimental evidence, the perpendicular transport coefficients are assumed to vary proportionally to $\propto 1/I_p$ to simulate the improvement of confinement time as predicted by global scaling law Editors et al. (1999). At $t = 0$ s the transport coefficients are set all equal to $2.68 m^2/s$, and then decrease to $D = \mu = 0.3 m^2/s$, $\chi_i = \chi_e = 1 m^2/s$ when the steady state phase of the discharge is reached. The time evolution of the plasma perpendicular transport coefficients is shown in Fig. 2b. The recycling coefficient at the wall boundary R is set equal to 0.99 and It's kept fixed in time and space for all the density scenario analyzed. The pump out of neutrals in the PFR is described by setting $\Gamma_{pump} = (1-R)\Gamma_{nn} \cdot n$, with $\Gamma_{nn} = -D\nabla n_n$ denoting the neutral flux vector.

### 3.2. Plasma density scan

The time evolution of plasma quantities at the plasma facing components (PFCs) is here analyzed for three scenarios where $n_{li}$ is linearly raised to fraction of [10%, 20% 25%] of $n_{Gw}$ with the same evolution of the magnetic configuration. We implemented a feedback filter on the gas puff rate to control the plasma density rises at the prescribed fraction of $n_{Gw}$. For computational reasons, the analysis is limited to the first 20 s of the discharge. Simulations of the entire discharge, including the steady state and ramp-down phases, will be performed in a future work. The time evolution of $n_{li}$, the gas puff rate *Puff*,

<figure><img src="image_4.png" /><figcaption>Figure 3: From top to bottom, time evolutions for the scan in density ramp-up of the line integrated density, puff rate, injected total power and on axis ion and electron temperatures. Time interval where L-D transition occurs is highlighted in red</figcaption></figure>

the injected power *P* and the on axis values of $T_i$, $T_e$ are shown in Fig. 3. In absence of auxiliary sources, here the total injected power is equal to the Ohmic heating source. Typical 2D spatiotemporal maps of ion particle and heat fluxes and temperatures as a function of the distance along the wall $S_W$ are shown in Fig. 4 a,b. The 0 coordinate is located in the middle of the pump below the dome, and moves poloidally counterclockwise
<figure><img src="image_5.png" /><figcaption>Figure 4: Spatio-temporal maps for the ion particle flux (a) and the heat flux (b) projected in the direction normal to the wall, as a function of the wall coordinates $S_W$ and time. The case with $n_{li} = 20\% n_{Gw}$ is here shown.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: Spatio-temporal map of the ion and electron temperatures along the OT target location for $n_{li} = 10\% n_{Gw}$ (first row) and $n_{li} = 25\% n_{Gw}$ (second row).</figcaption></figure>

along the tokamak PFCs, as labeled in Fig. 2a, towards the high field side. Results show the inner wall, and the outer and inner targets are the PFCs mostly interested by plasma fluxes. One can also note there are no flux at the divertor targets before the L-D transition ($t < 8$ s), while the fluxes peak at the inner and outer target when the X-point is established. Even if the magnitude of the fluxes present a quite moderate intensity in this phase of the discharge below the engineering limit of 10 MW/m⁻²), temperatures at the target are still critical for tungsten sputtering.

The 2D spatio-temporal maps of ion and electron temperature for $n_{li} = 10\% n_{Gw}$ in Fig. 5, show peaks up to 300 eV for $T_i$ and 100 eV for $T_e$ in the time interval right prior and after the onset of the X-point ($7 s < t < 9 s$). During the L-D transition, both density and temperature profiles follow the evolution of the magnetic equilibrium. However, the large parallel conduction rapidly flatten the temperature profile allowing for very large temperature values for fast transitory events. The peak values reduce by a factor three in the case $n_{li} = 25\% n_{Gw}$ highlighting a more favorable regime to perform the L-D transition. The tungsten sputtering rate increases significantly with the energy of impinging ions once it surpasses the sputtering threshold?. For a deuterium ion, this threshold is reached for values of $T_i$, $T_e \ge 45$ eV. For instance, SOLPS-ITER reactor simulations indicate substantial tungsten sputtering in regions where $T_i$ is approximately 100 eV, despite the moderate levels of $q \cdot n$ and $\Gamma \cdot n$?. Thus, further investigations into tungsten erosion in areas with high $T_i$ and $T_e$ are essential. Fig. 6 shows the comparison of the time evolutions of the maximum value achieved for the ion particle flux, heat flux and temperatures at the OT (first row), IW (second row) and IT (third row) location for the three density scenarios analyzed. The beneficial effect of the L-D transition at higher density is confirmed. For lager values of $n_{li}$ the ion particle flux remains almost unchanged at the IT while decreases at the OT, the heat flux decreases, and the temperatures reduce more than a factor two. However, unless for temperature, the ion particle flux and heat flux continue to increase at the IW. The simulations of the L-D transition at larger values of $n_{li}$ lead to numerical instabilities and to the computations overflow due to very low temperatures in the divertor area.

### 3.3. Anomalous pinch

Diffusion driven perpendicular transport inevitably leads to flat density profiles in the core region, which disagree with experimental evidences Kudashev et al. (2024a). In the simulations presented here we investigate the possibility of simulating anomalous or neoclassical pinch effects adding a convective flux term $\nabla \cdot (-v_p n \mathbf{b}_\perp)$ to Eq. 1. We compare the development of a particle pinch defined as a function of the normalized magnetic flux $\psi_N$, with a collisionality dependent one as used in Asp et al. (2022). The 1D profile of $v_p$ at midplane is shown in Fig. 7 and defined as:

$$v_p(\psi_N) = v_0 \times \left( \psi_N^2 + \psi_N^2 \times \tanh \left( \frac{0.92 - \psi_N}{0.02} \right) \right) \quad (8)$$

$$v_p(v^*) = f(v^*)0.5D \sqrt{\psi_N} H(\psi_N = 0.9) \quad (9)$$

$v_0$ is the constant value used to tune the magnitude of $v_p$ in order to have comparable values for the two formulations. D is the particle perpendicular diffusion coefficient, $v^* = 6.921e-18 \times \frac{R_{0qzn_{eff}}\Lambda_e}{\epsilon^2 T_e^2}$ Sauter et al. (1999) is the normalized electron collisionality, $\Lambda_e = 31.3 - \ln(\frac{\sqrt{\bar{n}}}{\bar{v}^*})$ is the Coulomb logarithm, q the safety factor profile and $f(v^*) = \min(1, \exp(1 - v^*/v_{th}))$. To take into account the larger collisionality at low values of the injected power, the constant $v_{th} = v_{th}(I_p)$ has here a tanh decay from the plasma current and assumes value in the range [0.04, 1].
<figure><img src="image_7.png" /><figcaption>Figure 6: Time evolutions of the ion particle flux (a), heat flux (b) and ion and electron temperatures (c) achieved for the scan in density ramp-up. Each row corresponds to the maximum value achieved at different location along the wall as defined in Fig. 2a.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: Plasma density and anomalous pinch midplane profiles for different formulations of the latter (a), and the respective ion and electron temperature profiles (b),</figcaption></figure>

as computed by SolEdge-HDG at $t = 20$ s.

<figure><img src="image_9.png" /><figcaption>Figure 8: Time evolutions of the line integrated density, the puff rate, the total power injected and the on axis ion and electron temperatures in the case a pinch term is included. Density ramp-up is now fixed at ≈ 35% $n_{Gw}$.</figcaption></figure>

We propose here the comparison between three different cases: one with $v_p = v_p(\psi_N)$, $v_0 = 0.3$, and two with $v_p = v_p(v^*)$, where in one of this we add an external source of power to

simulate ECRH heating. The heating source is described by an ad-hoc Gaussian profile which intensity linearly ramp-up to 2 MW in the time interval $10 \text{ s} < t < 15 \text{ s}$ (see third row in Fig. 8). The L-D transition is here performed at $n_{li} = 35\% n_{Gw}$ and the perpendicular transport coefficients assume the same values as the ones used in section 3.2. The inward pinch velocity speeds up the transport of the recycled plasma source toward the core center, allowing for larger values of plasma core density compatible with the SOL and target ones. Fig. 8 shows the time evolution of $n_{li}$, Puff, P, $T_{i,e}$ and of the applied ECRH source, while Fig. 9 shows the comparison of the time evolutions of the maximum ion particle flux, the heat flux and the temperatures at the OT, IW, IT location. The time evolution of $n_{li}$ (see first row Fig. 8) does not show significant changes proving the code is able to match the required target density for all the three cases. The gas puff rate slightly decreases when the ECRH heating is applied due to the lower collisionality, and hence a larger pinch velocity. The time signals at the PFCs mostly overlap up to $t = 11 \text{ s}$, then the applied ECRH source leads to larger fluxes and temperatures almost everywhere, except at the IW location, where a steepening of the density and energy gradients is reported at the separatrix, which reduces the heat and particle loads at the IW. $T_e$ increases by a factor 1.5, while $T_i$ by a factor larger than 2. Since $v_p(\psi_n)$ reaches its maximum nearby the separatrix, this causes larger density values and more energy dilution in the plasma edge, leading to a generally smaller fluxes and temperatures peaks at the PFCs, especially during the L-D transition. The mid-plane plasma density profiles, and the respective ion and electron temperature profiles, at $t = 20 \text{ s}$ for the three cases presented are shown in Fig. 7. The collisional pinch leads to more peaked density profile providing a self-consistent description of the density peaking as a function of the local collisionality. This one is influenced by the ECRH heating source which enlarges the steepness width and intensity. A broadening of the temperature profiles is also reported when Eq. 8 instead of Eq. 9 is used to compute $v_p$.
<figure><img src="image_10.png" /><figcaption>Figure 9: Temporal evolution of the ion particle flux (a), the heat flux (b) and the ion and electron temperatures (c) for different anomalous pinch velocities. Each row corresponds to the maximum value achieved at different locations along the ITER wall as defined in Fig. 2a.</figcaption></figure>

## 4. Impact of neutral modelling

In the context of developing more reliable reduced and fast fluid models, the purely diffusive transport model for the density is enriched here by an equation for the neutral parallel momentum. The velocity at which neutrals depart from the wall determines the mean free path and the location of the ionization front. Moreover, as the system approaches attached or partially

<figure><img src="image_11.png" /><figcaption>Figure 10: Plasma and neutral density profiles comparison at the IT and OT employing the $N_n$ and the $N_n \Gamma_{nn}$ neutral fluid models. The neutral pressure distribution in the divertor region is also shown.</figcaption></figure>

detached regimes, the large plasma density source in front of the target, in combination with low temperatures, leads to non physical values of the neutral diffusion coefficient $D_{nn}$ as used in Eq. 5. The proposed model is a simplified version of the pressure diffusion and parallel momentum equation derived in Horsten et al. (2017). Dependency of the neutral particle flux from the neutral temperature and the off diagonal losses of parallel neutral momentum are neglected. With this assumption the system of equation writes as:

$$ \begin{aligned} \partial_t n_n + \nabla \cdot (n_n u_n \mathbf{b} - D_{nn} \nabla n_n) &= S_{nn} \\ \partial_t (n_n u_n) + \nabla \cdot (n_n u_n^2 \mathbf{b} - \eta \nabla (u \mathbf{b})) + \nabla_{\|} (n_n k_b T_n) &= S_{\Gamma_{nn}} \end{aligned} \quad (10) $$

where $D_{nn}$ is the same as in Eq. 5, and $\eta = \frac{n_n k_b T_n}{m_n n(\sigma v)_{cx}}$ is the neutral parallel viscosity, which is a function of the neutral pressure and charge-exchange reaction rate. $T_n = T_i$ is assumed. The recycling condition is then imposed by replacing the left hand side of Eq. 7 with the new expression of the neutral flux $\Gamma_{nn} \cdot \mathbf{n} = [n_n u_n \mathbf{b} - D_{nn} \nabla n_n] \cdot \mathbf{n}$ and constraining the neutral parallel speed to be equal and opposite to the plasma one $u_n = -u \mathbf{b} \cdot \mathbf{n}$. One can note that in the case of small values of $D_{nn}$ the neutral particle flux becomes mostly convective.

We set-up a reference test-case where we compare the solutions obtained using Eq. 5 ($N_n$ model) or Eqs. 10 ($N_n \Gamma_{nn}$ model) as neutral model. The domain of simulation corresponds to the whole ITER poloidal cross section and we focus on the steady state using the magnetic equilibrium at $t = 75$ s shown in Fig. 2. The total injected power is equal to the Ohmic heating source as described in section 3. The perpendicular transport coefficients are constant in space and time, and set all equal to $1 \text{ m}^2/\text{s}$, $R = 0.98$ everywhere at the boundary, $\Gamma_{puff} = 1.25 \times 10^{21} \text{ p/s}$ and we set an outgoing neutral flux $\Gamma_{pump} = (1 - R) \Gamma_{nn} \cdot \mathbf{n}$ in the private flux boundary to mimic the presence of a pump. Fig. 1 shows the 2D map of the plasma density and the electron temperature as computed by using the $N_n \Gamma_{nn}$ model. The comparison of the plasma and neutral density profiles at the IT and OT in Fig. 10 shows comparable density values for the two models. However we report sharper profiles with faster exponential decay out the separatrix when neutral parallel momentum equation is solved. One can also note a larger asymmetry and a remarkable shift of the peak location, which is now closer to the strike point, when the $N_n \Gamma_{nn}$ model is employed. The 2D distribution of the neutral pressure in the lower divertor region confirms a reduced transport of neutrals in the PFR and towards the extremity corner of the divertor legs. The neutral model is based on the charge-exchange/ionization processes dominance, and thus does not describe correctly the neutral dynamics in the regions without plasma, such as the PFR. A full kinetic model should be used there. Future studies are forthcoming to investigate the reason of these different behaviours.
## 5. Summary and conclusions

This work reports the first core-edge non steady calculations of different ITER's Q = 10 scenarios. This is made possible using the transport code SolEdge-HDG which relies on an high-order finite-elements method that guarantees an accurate description of the tokamak geometry including the variation of the magnetic equilibrium. The ion particle flux, the heat flux and temperatures at the PFCs have been investigated simulating a linear ramp-up of the central line integrated density at different fraction of the Greenwald density limit. In the case of purely diffusive perpendicular transport, the density ramp-up at fraction larger than 25% of $n_{Gw}$ leads to numerical instabilities. However, convective flux terms prove to be significant for core refuelling, extending this interval up to value of $n_{li} \approx 35\% n_{Gw}$. Despite the moderate particle and heat loads at the PFCs, the temperatures are still critical, with large spots reached at the inner wall and at the inner and outer target respectively prior and after the X-point formation. The implementation of a conservation equation for the parallel neutral momentum shows the neutral speed is a critical factor for generating non-symmetric plasma density sources between the inner and outer targets.

This work demonstrates the promising capabilities of SolEdge-HDG to investigate transport during the ramp-up phase in ITER using integrated core-edge fluid simulations. The present results underscore the importance of considering non-steady transitory events to accurately assess the impact of tungsten erosion on plasma performances.

## Acknowledgements

This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No 101052200 EUROfusion). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them. This work has been carried out thanks to the support of the AMIDEX project (ANR-11IDEX-0001 02, SAPHIR project) funded by the "Investissements d'Avenir" French Government program, managed by the French National Research Agency (ANR). This work also received support by the US DoE under DE-AC02-09CH11466.

## References

Asp, E.M., Corrigan, G., da Silva Aresta Belo, P., Garzotti, L., Harting, D., Köchl, F., Parail, V., Cavinato, M., Loarte, A., Romanelli, M., Sartori, R., 2022. Jintrac integrated simulations of iter scenarios including fuelling and divertor power flux control for h, he and dt plasmas. Nuclear Fusion 62, 126033. URL: https://dx.doi.org/10.1088/1741-4326/ac90d4, doi:10.1088/1741-4326/ac90d4.

D. Reiter, M.B., Börner, P., 2005. The eirene and b2-eirene codes. Fusion Science and Technology 47, 172-186. URL: https://doi.org/10.13182/FST47-172, doi:10.13182/FST47-172, arXiv:https://doi.org/10.13182/FST47-172.

d'Abusco, M.S., Giorgiani, G., Artaud, J., Bufferand, H., Ciraolo, G., Ghendrih, P., Serre, E., Tamain, P., 2022. Core-edge 2d fluid modeling of full tokamak discharge with varying magnetic equilibrium: from west start-up to ramp-down. Nuclear Fusion 62, 086002. URL: https://dx.doi.org/10.1088/1741-4326/ac47ad, doi:10.1088/1741-4326/ac47ad.

Editors, I.P.B., Chairs, I.P.E.G., Co-Chairs, Team, I.J.C., Unit, P.I., 1999. Chapter 1: Overview and summary. Nuclear Fusion 39, 2137. URL: https://dx.doi.org/10.1088/0029-5515/39/12/301, doi:10.1088/0029-5515/39/12/301.

Garzotti, L., Belo, P., Corrigan, G., Harting, D., Köchl, F., Loarte, A., Asp, E.M., Parail, V., Ambrosino, R., Cavinato, M., Mattei, M., Romanelli, M., Sartori, R., Valovič, M., 2018. Integrated core-sol modelling of fuelling, density control and divertor heat loads for the flat-top phase of the iter h-mode d-t plasma scenarios. Nuclear Fusion 59, 026006. URL: https://dx.doi.org/10.1088/1741-4326/aaf2f3, doi:10.1088/1741-4326/aaf2f3.

Giorgiani, G., Bufferand, H., Ciraolo, G., Ghendrih, P., Schwander, F., Serre, E., Tamain, P., 2018. A hybrid discontinuous galerkin method for tokamak edge plasma simulations in global realistic geometry. Journal of Computational Physics 374, 515-532. URL: https://www.sciencedirect.com/science/article/pii/S0021999118304893, doi:https://doi.org/10.1016/j.jcp.2018.07.028.

Giorgiani, G., Bufferand, H., Schwander, F., Serre, E., Tamain, P., 2020. A high-order non field-aligned approach for the discretization of strongly anisotropic diffusion operators in magnetic fusion. Computer Physics Communications 254, 107375. URL: https://www.sciencedirect.com/science/article/pii/S0010465520301612, doi:https://doi.org/10.1016/j.cpc.2020.107375.

Horsten, N., Samaey, G., Baelmans, M., 2017. Development and assessment of 2d fluid neutral models that include atomic databases and a microscopic reflection model. Nuclear Fusion 57, 116043. URL: https://dx.doi.org/10.1088/1741-4326/aa8009, doi:10.1088/1741-4326/aa8009.

Kudashev, I., Glasser, A.M., d'Abusco, M.S., Serre, E., 2024a. Impact of variable perpendicular transport coefficients in west simulations using soledge-hdg. IEEE Transactions on Plasma Science, 1-6doi:10.1109/TPS.2024.3384031.

Kudashev, I., M. Scotto d'Abusco, A. M. Glasser, E.S., F. Schwander, H.B., G. Ciraolo, P.G., Tamain, P., 2024b. Global particle buildup simulations with gas puff scan. application to west discharge. Frontier Physics, Fusion Plasma Physics 12. doi:10.3389/fphy.2024.1407534.

Kukushkin, A., Pacher, H., Loarte, A., Komarov, V., Kotov, V., Merola, M., Pacher, G., Reiter, D., 2009. Analysis of performance of the optimized divertor in iter. Nuclear Fusion 49, 075008. URL: https://dx.doi.org/10.1088/0029-5515/49/7/075008, doi:10.1088/00295515/49/7/075008.

Reiter, D., 2000. The data file amjuel: Additional atomic and molecular data for eirene.

ROMANELLI, M., CORRIGAN, G., PARAIL, V., WIESEN, S., AM-BROSINO, R., BELO, P.D.S.A., GARZOTTI, L., HARTING, D., KÖCHL, F., KOSKELA, T., LAURO-TARONI, L., MARCHETTO, C., MATTEI, M., MILITELLO-ASP, E., NAVE, M.F.F., PAMELA, S., SALMI, A., STRAND, P., SZEPESI, G., Contributors, E.J., 2014. Jintrac: A system of codes for integrated simulation of tokamak scenarios. Plasma and Fusion Research 9, 3403023-3403023. doi:10.1585/pfr.9.3403023.

Romanelli, M., Parail, V., da Silva Aresta Belo, P., Corrigan, G., Garzotti, L., Harting, D., Koechl, F., Militello-Asp, E., Ambrosino, R., Cavinato, M., Kukushkin, A., Loarte, A., Mattei, M., Sartori, R., 2015. Modelling of plasma performance and transient density behaviour in the h-mode access for iter gas fuelled scenarios. Nuclear Fusion 55, 093008. URL: https://dx.doi.org/10.1088/0029-5515/55/9/093008, doi:10.1088/0029-5515/55/9/093008.

Sauter, O., Angioni, C., Lin-Liu, Y.R., 1999. Neoclassical conductivity and bootstrap current formulas for general axisymmetric equilibria and arbitrary collisionality regime. Physics of Plasmas 6, 2834-2839. URL: https://doi.org/10.1063/1.873240, doi:10.1063/1.873240, arXiv:https://pubs.aip.org/aip/pop/article-pdf/6/7/2834/19054899/2

Simonini, R., Corrigan, G., Radford, G., Spence, J., Taroni, A., 1994. Models and numerics in the multi-fluid 2-d edge plasma code edge2d/u. Contributions to Plasma Physics 34, 368-373. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/ctpp.2150340242, doi:https://doi.org/10.1002/ctpp.2150340242, arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/ctpp.2150340
Summers, H., 2004. The adas user manual, version 2.6. http://www.adas.ac.uk/.

Wiesen, S., Reiter, D., Kotov, V., Baelmans, M., Dekeyser, W., Kukushkin, A., Lisgo, S., Pitts, R., Rozhansky, V., Saibene, G., Veselova, I., Voskoboynikov, S., 2015. The new solps-iter code package. Journal of Nuclear Materials 463, 480–484. URL: https://www.sciencedirect.com/science/article/pii/S0022311514006965. doi:https://doi.org/10.1016/j.jnucmat.2014.10.012. pLASMA-SURFACE INTERACTIONS 21.
## Highlights

* Core-edge non-steady calculations of ITER's Q = 10 L-D transition have been performed with SolEdge-HDG including the variation of the magnetic equilibrium

* The scan in plasma density shows, despite moderate particle and heat fluxes at the PFCs, areas where large temperature are achieved prior and after the X-point formation

* Convective flux terms are important for core refuelling allowing for fast transport of the recycled plasma sources at the wall

* Preliminary comparison of different fluid neutral models implemented in SolEdge-HDG shows the critical role neutrals modelling plays for a correct assessment of plasma quantities at the divertor target

* Future studies will aim to investigate the W erosion and impurity contamination for fast transitory event like L-D transition
Declaration of interests

☐ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

☒ The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

Manuel Scotto d'Abusco reports financial support, equipment, drugs, or supplies, and writing assistance were provided by Aix-Marseille University. Manuel Scotto d'Abusco reports a relationship with Aix-Marseille University that includes: employment. No If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.