Reduced transport models for a tokamak flight simulator

The ability to predict the evolution of a new plasma discharge can allow one to improve the efficiency of the discharge itself. Moreover, it has become increasingly clear that a fast and reliable predictive ability for the entire evolution of a planned tokamak discharge is required for the operation of reactors like ITER and DEMO, before the plasma discharge is actually performed. This prevents undesired events, like disruptive cases, from occurring in actual experiments, avoiding major damage and the interruption of the reactor operation, with consequent costs of reparations and long delays. Several ways to perform these predictions already exist, and they differ in the level of detail and required computational time. The easiest formulae available are the scaling laws [2, 3], which predict the confinement time of a stationary phase of a discharge. However, they cannot be used to predict the plasma evolution of an entire discharge, including transients. Moreover, they are 0D and miss a physics-based nature, so their applicability and reliability is reduced. In fact, they miss the profile effects, which can be very important for the discharge evolution and to predict certain conditions, like internal transport barriers (ITBs). The predictions can be also made by adopting different suites or codes. For this purpose, more complex models like gyrokinetic or fluid codes can deeply describe the physics, but they are time-consuming. Furthermore, during experimental campaigns some 'last-minute' changes in the discharge program could be applied due to technical reasons, just before running the discharge. Therefore, the need for fast simulations based on first principles has become clear. A predictive tool that meets these prerogatives is the tokamak flight simulator, Fenix [4, 5]. This is a numerical tool that predicts the plasma behavior using the discharge program as input. It is based on the interaction between a control system, which is simulated within Simulink, and a physics model, which mainly includes equilibrium and transport models. It can ensure that either actuator trajectories or plasma parameters satisfy the experimental goals and reduce the probability of plasma disruptions and of exceeding operational limits. While the plasma equilibrium and the control algorithm can already be described analytically, a fully integrated analytical transport model is still missing. This has to be physics-based to be realistic, but also fast enough to be used as an inter-discharge prediction tool. This compromise can be reached by employing analytical models that are derived from first principles theories. An important task in the development of these models is to include as few experimental inputs as possible to strengthen the prediction capability. As a starting point, Fenix was initially tested with models in which the transport was tuned to match some experimental measurements.

An alternative solution based on first principles, which strongly reduces the amount of experimental input, is proposed in this work using a set of four transport physics-based models. Two models are applied in the confined region, depending on the regime of the plasma, while the other two act in the scrape-off layer (SOL). They are linked through the last closed flux surface (LCFS). At this point, the models of the SOL give boundary conditions for the confined region, taking as input engineering parameters, which can be scheduled during the discharge's plan. Moreover, these fast analytical models can be implemented in the flight simulator, so that it can be used inter-discharge to predict the discharge evolution. When coupled with a predictive model for the disruptions, Fenix could also allow the avoidance of disruptive cases. Such a model still does not exist, but its development has been planned for future work. A possible simple approach could be to couple Fenix with a space state description that sets the conditions for which the plasma hits some hard limits (e.g. high density limit, as in [6]).

The models have been developed in the framework of the transport code ASTRA [7, 8]. This framework allows us to couple the different parts of the tokamak simultaneously, giving boundary conditions to each other. Although ASTRA, like any other transport code for confined plasma, does not contain any radial grid beyond the LCFS, additional subroutines are flexibly coupled, which (virtually) simulate the evolution of the SOL concomitantly with the evolution of the confined plasma. The analytical nature of all the models allows for an iteration-free working principle in which all the models proceed concomitantly, interacting with one another through their boundaries.

In section 2 the core model is shown, in section 3 the edge model is presented, while in section 4 models for the SOL are explained. Section 5 shows six first fully-integrated simulations for five H-mode discharges and an L-mode AUG discharge in the Fenix flight simulator. Conclusions are drawn and the outlook described in section 6.

The core of a tokamak plasma is defined here as the region up to $\rho_t = 0.9$, where ρ_t is the toroidal magnetic flux minor radius normalized to the value at the separatrix. This definition of the core comes from the choice of the pedestal width. In fact, this has been fixed equal to 0.1 in ρ_t units. Moreover, almost all the simulations run to build up the core transport model were performed along the radius up to $\rho_t = 0.9$, which is the region where the code used to perform such simulations has already been widely validated in earlier works.

The plasma reaches its highest values of pressure in the core. In this region, when the kinetic pressure is high enough, the fusion reactions can happen. Nevertheless, some dynamics are found to be opposed to the increase of plasma energy, and in some cases they can set a constraint on the maximum reachable pressure and energy. The main phenomena involved in this are the magnetohydrodynamic (MHD) instabilities and the heat and particle transport. The last two are in turn split into neoclassical and anomalous parts. This work focuses on the anomalous (turbulent) transport. Some fast models already exist and they are based on neural networks (NNs), like Qualikiz NN [9]. However, NNs need training on a wide database which can assure that the model will work properly inside the range of variation of the data collected, but it could be unsuccessful outside it, so its use could be compromised. This is particularly important when one extrapolates to larger devices, which is a crucial problem for machines of the next generation. Moreover, even if the NNs can provide a correct prediction of transport coefficients, their structure does not offer physics insights. For all these reasons an alternative model, based on analytical formulae fitted over a smaller database than the one usually needed from a NN, is proposed here. Its theoretical nature, which is based on first principles, offers more transparency and physics insight, as well as an easier extrapolation to larger or different devices. In fact, the coefficients derived in this work are multiplied by a gyroBohm scaling factor which takes into account the dependence on $\rho^{*}$, which is the main parameter to take into account in the extrapolation to bigger machines. In order to define the formulae of the transport coefficients, the main instabilities affecting the core have to be considered. In this region they are the micro-instabilities driven by temperature and density gradients. It has been shown by experiments that this transport has a threshold nature [10, 11]. It is also known that there are some stabilizing effects on these instabilities (e.g. collisionality, magnetic shear, impurity concentration, β, concentration of fast particles and E × B shear) [12–17] and that shape plays a role [18]. Therefore, taking into account all these aspects, some threshold formulae have been adopted for ion temperature gradient (ITG), electron temperature gradient (ETG) and trapped electron mode (TEM) and fitted to a database of TGLF [19–21] simulations, while micro-tearing modes (MTMs) [22–24] have been neglected, because they are usually negligible with respect to the ITG/TEM transport in the core region. This is anyway not always true, especially in spherical tokamaks [25] or in some configurations with high β, which is the main drive of MTM, or in some ITG suppressed advance scenarios. However the database used in the fitting does not include such discharges, but the development of a formula to take into account this instability in an extended database of discharges is planned for the future. TGLF simulations were set in this way: saturation rule 2 (which sets a specific rule for the fitting of the quasilinear fluxes of TGLF on some nonlinear gyrokinetic simulations), three species (electron, D and B/N) and electromagnetic effects. The database consists of some stationary phases of 15 AUG discharges from different scenarios (H-mode, L-mode, I-mode and negative triangularity). Each of these cases has been perturbed in the boundary condition of the kinetic profiles by enhancing/reducing by 10% electron density, electron temperature and ion temperature at $\rho_t = 0.9$. This scan has increased the database and assured that the range of variation of the normalized gradient of the kinetic profiles is covered in a more homogeneous way, with a contiguous spread of the data. This should improve the quality of the fitting. The broad variety of configurations included in the database is largely justified by the general nature of the instabilities treated. In this database, six coordinates in the range $\rho_t = [0-0.9]$ have been considered for the last 20 time steps of the converged TGLF+ASTRA simulations. The simulations included sawteeth, but a strategy to exclude their impact on transport has been used, based on varying the position of the first radial coordinate included in the fitting. Moreover, some coordinates at $\rho_t\gt 0.9$ for one L-mode case are included. Thus, it is possible to use this model in the edge in the L-mode configuration, where turbulent transport is still determining the profiles. The database consists overall of 12 600 values of $\chi_\textrm{e}{}$ and $\chi_\textrm{i}{}$, which are respectively electron and ion heat diffusivity. The database is much smaller than the one usually needed in a NN model, but considering that the maximum number of fitting parameters involved is 15, there are almost 1000 data for each free fitting parameter, and so we have confidence that there are enough data to catch all the dependencies of the formulae and to give a reliable result. The formulae used for the fitting are reported here, in gyroBohm units [26], that is $T^{1.5}_e B^{-2}r^{-1}$, where r is the minor radius:

$$\begin{align} \hat{\chi}_{i,\mathrm{ITG}}=C\cdot H_{\mathrm{ITG}}\left(\frac{R}{L_{Ti}}-\frac{R}{L_{Ti,\mathrm{ITG}}}\right)^{\epsilon_{10}} q^{\gamma_q} e^{-\gamma_{\beta_e}\beta_e} k^{-\gamma_k} e^{-\gamma_\mathrm{imp}(1-c_I)}, \end{align} \tag{ 1 }$$

$$\begin{align} \hat{\chi}_{e,\mathrm{ETG}}=D\cdot H_{\mathrm{ETG}}\left(\frac{R}{L_{Te}}-\frac{R}{L_{Te,\mathrm{ETG}}}\right)^{\epsilon_{20}} q^{\gamma_{q,e}} k^{-\gamma_{k,e}}, \end{align} \tag{ 2 }$$

$$\begin{align} \hat{\chi}_{e,\mathrm{TEM}}=D_2 H_{\mathrm{TEM}}\left(\frac{R}{L_{Te}}-\frac{R}{L_{Te,\mathrm{TEM}}}\right)^{\epsilon_{30}} e^{-\gamma_\nu \nu} e^{-\gamma_s s} e^{\gamma_{\delta,e}\delta}, \end{align} \tag{ 3 }$$

$$\begin{align} \hat{\chi}_{e,\mathrm{ITG}}=\mathrm{max}\left\{1,f_t D_3\frac{L_{Te}}{L_{Ti}}\right\}\hat{\chi}_{i,\mathrm{ITG}}, \end{align} \tag{ 4 }$$

$$\begin{align} \hat{\chi}_e=(1-f_t)\hat{\chi}_{e,\mathrm{ETG}}+f_t\hat{\chi}_{e,\mathrm{TEM}}+\hat{\chi}_{e,\mathrm{ITG}}, \end{align} \tag{ 5 }$$

$$\begin{align} \hat{\chi}_i=\hat{\chi}_{i,\mathrm{ITG}}, \end{align} \tag{ 6 }$$

where R is the major radius, $1/L_T = -\frac{\partial_r T}{T}$, q is the safety factor, β_e is the ratio between electronic kinetic and magnetic pressure, k is the local elongation, δ is the local triangularity, $c_\mathrm{I}$ is the main ion fraction, ν is the collisionality, defined as $R_\mathrm{tor} Z_\mathrm{eff} n_e T^{-2}_e$ with n_e in $10^{19} \mathrm{m}^{-3}$, T_e in keV and $Z_\mathrm{eff}$ equal to effective charge, s is the magnetic shear and f_t , defined as $\sqrt{\frac{r}{R}}$, is analogous of the trapped particle fraction, while C, ₁₀, γ_q , $\gamma_{\beta_e}$, γ_k , $\gamma_\mathrm{imp}$, D, ₂₀, $\gamma_{q,e}$, $\gamma_{k,e}$, D₂, ₃₀, γ_ν , γ_s , $\gamma_{\delta,e}$, D₃ are the fitting parameters. H is a Heaviside function, whose argument is the difference between the normalized temperature gradient and its threshold value. ₁₀, ₂₀ and ₃₀ are the exponents of the differences between the actual normalized gradient and its threshold, representing thus the stiffness, which determines how the kinetic profiles react to the enhanced heating. The electromagnetic effects due to the fluctuations of magnetic fields are included through the β term, as a combination of magnetic flutter and stabilization of ITG modes, resulting in an overall reduction of transport. The shape of the plasma also affects the transport, because it changes the relative length that the particles travel in the high field side (HFS) and in the low field side (LFS), affecting thus the ballooning stability. The shape also changes the local concentration of particles, creating an additive poloidal redistribution. In order to take into account these effects, δ and k are included in the formulae. The collisionality is a crucial parameter to distinguish between different instability regimes. In particular, it reduces the effect of the trapped particle fraction, through collisions with passing electrons. This is modeled in $\hat{\chi}_{e,\mathrm{TEM}}$, through $e^{-\gamma_\nu \nu}$. The magnetic shear reduces the transport, through the tilt of the field lines, which leads to a decorrelation of the neighboring field lines. This effect is taken into account through the terms related to q and s in the formulae. For $\hat{\chi}_{e,\mathrm{ITG}}$ a proportionality with $\hat{\chi}_{i,\mathrm{ITG}}$ has been assumed, through the trapped particle fraction (because the passing electrons are considered adiabatic for this instability) and a ratio between electron and ITG (to simulate the relative contribution of electrons to a ionic instability). The linear threshold formulae are available from the literature [27, 28]. For the ETG and ITG thresholds, some fitting parameters have been added to multiply some dependencies:

$$\begin{align} \frac{R}{L_{Ti,\mathrm{ITG}}}=\mathrm{max}\biggl\{& A_{10}\left(1+B_{10} Z_\mathrm{eff}\frac{T_i}{T_e}\right) \left(1+B_{20}\frac{s}{q}\right)\nonumber\\ & \times \;\;\left(1-1.5{f_t}^2\right)\left[1+0.3\left(k-1\right)\right]; A_{20}\frac{R}{L_{n}}\biggr\}, \end{align} \tag{ 7 }$$

$$\begin{align} \frac{R}{L_{Te,\mathrm{ETG}}}=\mathrm{max}\biggl\{& F_{10}\left(1+G_{10} Z_\mathrm{eff}\frac{T_e}{T_i}\right) \left(1+G_{20}\frac{s}{q}\right)\nonumber\\ & \times \;\;\left(1-1.5{f_t}^2\right)\left[1+0.3\left(k-1\right)\right]; F_{20}\frac{R}{L_{n}}\biggr\}, \end{align} \tag{ 8 }$$

$$\begin{align} \frac{R}{L_{Te,\mathrm{TEM}}}=0.357\frac{f_t+0.271}{f_t} &\Big[4.9-1.31\frac{R}{L_n}+2.68\mathrm{s}\nonumber\\ & \;\; +\log \left(1+20\nu\right)\Big] \end{align} \tag{ 9 }$$

where $1/L_n = -\frac{\partial_r n_e}{n_e}$, while A₁₀, B₁₀, B₂₀, A₂₀, F₁₀, G₁₀, G₂₀, F₂₀ are other fitting parameters. The ITG threshold has the same shape as the ETG one, except $\frac{T_e}{T_i}$ is inverted, because in the linear limit the two instabilities have the same physical picture [29]. In the threshold formulae the dependence on $\frac{T_e}{T_i}$ appears to take into account the contribution of the ions and the electrons, respectively, to the electronic and ionic instabilities.

Starting from the outputs of the ASTRA simulations, some TGLF standalone simulations have been run for each radial coordinate of each discharge. It has been found that the contribution of transport due to $k_y\rho\gt 1$ instabilities is a very small percentage of the total transport. Considering that this is the only spectral region where ETG appears, the need for another database becomes clear, one that is more focused on the ETG-wavelength instability region. This has been obtained by multiplying the original database of $\chi_\textrm{e}{}$ in output from ASTRA by $\frac{Q_\mathrm{high}}{Q_\mathrm{tot}}$, where $Q_\mathrm{high}$ is the anomalous electron heat coming from $k_y\rho\gt 1$ transport, and $Q_\mathrm{tot}$ is the total transport over all the spectra. This ratio represents the percentage of electron heat transport due to small scale turbulence. Considering that the temperature gradient is constant over the spectra $\chi_e\frac{Q_{e,\mathrm{high}}}{Q_{e,\mathrm{tot}}} = \chi_e\frac{\chi_{e,\mathrm{high}}}{\chi_{e}} = \chi_{e,\mathrm{ETG}}$, because $\chi_{e,\mathrm{high}}$ is assumed to be due to ETG. This should provide a rough distinction between ion scales and mixed-electron scales. This procedure has led to a disentanglement of the microinstabilities, which are fitted separately, leading to a reduction in degrees of freedom, which goes more in the direction of a physics-based fitting rather than an artificial fitting, whose success is a coincidence of numerical optimizations. However, one could in principle contest the rough separation between ETG and TEM+ITG contributions made by splitting spectra in two clear regions. To overcome this problem, a first fitting on the ETG database has been performed by also including $\chi_{e,\mathrm{ITG}}$ formula, and it has been seen to match a branch of high $\chi_{e,\mathrm{TGLF}}$ values. In fact, a clear difference between the $\chi_{e,\mathrm{TGLF}}$ matched by the ITG formula and by the ETG formula has been found in terms of parameter dependencies and order of magnitude. A filter to separate these two branches of coefficients has been put in place by excluding data with both $\frac{R}{L_{T,e}}\lt 20$ and $\chi_\textrm{e}{}\gt 10$ (in gB units) conditions met. Those conditions have been found to belong to ITG-driven cases. Anyway, an unavoidable small residual contribution from TEMs is present in the ITG-filtered database, due to the fact that this instability appears also on electron scales. The ETG threshold has been fitted, because in this way it was found to better match the TGLF output. The results show stronger $\frac{T_e}{T_i}$ and lower s and $\frac{R}{L_n}$ contributions in the linear ETG threshold with respect to the formula from the literature. This effect could be due to a residual presence of $\frac{R}{L_n}$-driven TEM, which is not included in the model. In fact, a threshold formula is present in the literature, but the hypothesis under which it is valid is too constraining (no temperature fluctuations). Another cause of it could be the difference between the linear and nonlinear ETG threshold, where the second is found to produce stronger transport due to the formation of ETG streamers [29–32]. However, TGLF is a quasilinear model; then, the presence of a nonlinear threshold could be just the result of the multiscale spectral calibration over the gyrokinetic database fitted to obtain the saturation rule (i.e. is not a theoretical threshold but the result of a fitting procedure). An additional formulation of the density gradient drive of TEM is planned for the future. The TEM threshold has been chosen not to contain fitting parameters in order to reduce the degrees of freedom during the fitting process and optimize the iteration procedure. One can see that the main stabilizing effects have been taken into account, in the threshold or in the coefficient expressions, while the stiffness is assured by the power law on the difference between normalized gradients and their critical values. Also, cross dependencies, which are effects of a species-related instability on other species, have been included by $\hat{\chi}_{e,\mathrm{ITG}}$. This means that ITG in this model also drives some electronic heat transportin particular, due to the trapped electron contribution in the ITG mode.

Figures 1 and 2 show the scattering of the fitted coefficients compared to the respective values calculated by TGLF. The red points are the values calculated using the fitting formulae, while the respective TGLF values lay on the black solid line. The coefficients are shown in gyroBohm units. One can see that the general trend is reproduced, and some scattering should be tolerated considering that the coefficients are taken along the entire radius (up to separatrix in few cases) and the simple structure of the analytical formulae. Anyway, some analytical coefficients are completely mismatching the values from TGLF, which is the case of the spot of red points around $\chi_{e,\mathrm{TEM+ITG,fitted}} = 10^{-2}$. These cases were investigated to understand the nature of the mismatch and the safety factor was found to be below 1, so it is argued that the sawteeth model has increased the TGLF coefficients. The means and the standard deviations of the ratios between the fitted and the TGLF transport coefficients for $\chi_{i,\mathrm{ITG}}$, $\chi_{e,\mathrm{ITG+TEM}}$ and $\chi_{e,\mathrm{ETG}}$ are, respectively, $\mu_{i,\mathrm{ITG}} = 1.83$, $\sigma_{i,\mathrm{ITG}} = 2.5$, $\mu_{e,\mathrm{ITG+TEM}} = 1.4$, $\sigma_{e,\mathrm{ITG+TEM}} = 1.226$, $\mu_\mathrm{ETG} = 1.15$, $\sigma_\mathrm{ETG} = 1.27$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 3–5 assure the wide range of variations of most of the parameters and show that most of the parameter dependencies are well captured, proving the successful performance of the fitting routine. However, in figure 3(c) a branch of analytical χ_i could not match the TGLF values. In fact, in this plot, one can clearly see a bifurcation of the dependence of χ_i on $\frac{R}{L_{T_i}}$, where one branch is increasingly weaker and the other is very steep. In order to investigate the mismatching χ_i with strong $\frac{R}{L_{T_i}}$ dependence, the input parameters of two simulations with $9\lt \frac{R}{L_{T_i}}\lt 10$ giving, respectively, $\chi_i\gt 50$ and $\chi_i\lt 40$, were selected, and some scans around the nominal input parameters were performed with TGLF stand alone. After a detailed analysis, the main reasons for the difference between the transport coefficients in the two cases were identified for different values of collisionality and radial position (and trapped particle fraction). In particular, the collisionality of the high χ_i case was three times the one of the low χ_i , while the normalized minor radius of the former was 1.5 times the one of the latter. This is not surprising because high collisionality usually drives transport by shifting upwards the saturation level given by zonal flows [33]. Anyway, the similar values of the other parameters (in particular the high value of magnetic shear) should prevent the triggering of an ITG for both the cases, while the TEM should be reduced by high collisionality. The differences between the two χ_i have been compared to two corresponding linear calculations with the gyrokinetic code GENE, using the same input parameters as the TGLF cases. In the high collisionality case, the growth rate spectrum of GENE is found to be one order of magnitude smaller. Correspondingly, also a quasi-linear estimate of the ion heat flux using the GENE results and applying a quasi-linear rule, which is analogous to that of TGLF ,yields values which are one order of magnitude lower than the TGLF. This comparison with GENE, and the observation that the TGLF predicted ion temperature profiles in these conditions are significantly below the experimental ones, consistent with a TGLF overestimation of the heat conductivities in comparison to the power balance ion heat conductivities, suggested neglecting this branch of TGLF results in the derivation of the analytical model. Anyway, a further investigation with nonlinear gyrokinetic simulation is well beyond the scope of this work. It is worth noticing that the branch of mismatching χ_i in figure 3(c) is reflected also on χ_e in figure 4(c), through the term $\hat{\chi}_{e,\mathrm{ITG}}$ in equation (4).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In figure 5 each of the 12 plots shows a subset of TGLF χ_e which are not matched by the ETG analytical formula. This is probably caused by a minimum residual of ITG-driven transport in the TGLF database used for the ETG fitting. This is supported by the fact that the trend in $\frac{R}{L_{T_i}}$ is captured except for high values of the normalized gradient, which is consistent with the presence of ITG. Moreover, the mismatching transport coefficients are underestimated for low values of β, which reduces ITG. Some other mismatching χ_e could be due to the presence of TEM-driven TGLF transport coefficients, because most of these χ_e mismatch for high values of trapped particle fraction and high values of density gradient, which can drive TEM. The presence of residual TEM transport is also supported by the larger heat conductivities for high $\frac{R}{L_{T_e}}$. Anyway, the ETG contribution to the transport was found to be smaller than 5% for most cases, so its relevance is reduced in this work. A better formulation of the ETG transport is planned for the future.

The means and standard deviations of $\frac{\chi_\mathrm{fitted}}{\chi_\mathrm{TGLF}}$ for different values of the parameters included in the fitting are shown in figures 6–8 respectively for $\chi_{i,\mathrm{ITG}}$, $\chi_{e,\mathrm{TEM+ITG}}$ and $\chi_{e,\mathrm{ETG}}$. Here, one can see the most critical values of the parameters in which the standard deviation is high, then where the fitting procedure loses quality and does not match the TGLF values.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The output fitting parameters are:

$$\begin{align} & C=26.15, \epsilon_{10}=0.207, \gamma_q=0.0546, \gamma_{\beta_e}=1, \gamma_k=0.01,\nonumber\\ & \gamma_{imp}=1, D=0.116, \epsilon_{20}=1.18, \gamma_{q,e}=0.11, \gamma_{k,e}=0.097, \nonumber\\ & D_2=1, \epsilon_{30}=1.74, \gamma_\nu=0.0035, \gamma_s=0.001, \gamma_{\delta,e}=0.031,\nonumber\\ & D_3=0.92\,\, A_{10}=0.227, B_{10}=0.0119, B_{20}=2.645,\nonumber\\ & A_{20}=0.635\, F_{10}=1.8, G_{10}=1.27, G_{20}=0, F_{20}=0.1 .\end{align} \tag{ 10 }$$

The particle diffusivity has been assumed to be equal to $C\cdot\chi_e$, where C is a calibration factor which has been fixed equal to 1 to match the kinetic profiles of a discharge, while particle pinch has been modeled with a heuristic formula which is proportional to diffusivity to assure stationarity and takes into account the effect of $L_{T_e}$, s and ν [34–38]

$$\begin{align} v_p=-D_n \frac{\mathrm{max}\left\{0;0.2 R_\mathrm{tor}\frac{|\partial_r T_e|}{T_e}+0.15\mathrm{s}-\nu_\textrm{ei}{}/15\right\}}{R_\mathrm{tor}} .\end{align} \tag{ 11 }$$

The maximum function has been used here to assume a negative convection, that is a pinch, in order to have stationarity through the balance between particle diffusivity and particle pinch. This is not a really strong assumption for AUG discharges. A physics-based model for particle transport is a task that has not yet been performed, due to the difficulty of modeling the particle convection, whose value can change sign locally. However, an improvement of the model used here is planned for the future. The database shown has been built using TGLF with saturation rule 2, but earlier saturation rule 1 has also been adopted, in order to strengthen the fidelity of the procedure used in the model. Also, in this case , a small scattering of the TGLF coefficients around the fitting values has been found, but the main dependencies are well reproduced by the formulae.

The model described in the previous section is, in principle, extendable to the entire confined region. However, this is not always possible due to the different physics which govern the edge of the plasma in some specific configurations. In fact, in the high-confinement regime (H-mode) or the improved confinement regime (I-mode), which are observed when a threshold in heating power is overcome, a pedestal arises respectively in both density and temperature, or only in temperature profiles. This pedestal is formed in a small region of the plasma near the LCFS, or separatrix, but it plays an important role, because it sustains a sizeable fraction of the plasma energy. In this region, steep gradients in the kinetic profiles can be sources of free energy and the presence of the bootstrap current, which modifies the current profile and the safety factor, can affect the turbulence. Moreover, a combination of destabilizing and stabilizing effects, for instance connected with the strong gradients, the high safety factor, the collisionality and the plasma beta, as well as magnetic and rotational shear, imply that the pedestal can host a variety of instabilities and consequent turbulence featuring different thresholds and transport properties. Therefore, a quantitative prediction of the transport in these conditions is not a trivial task. Nevertheless, a simplified approach can be used in the flight simulator, since it will give the boundary condition for the core model even if the pedestal structure is not captured in detail. To distinguish this region from the core another model has been used starting from $\rho_t = 0.9$, which is assumed to be the position of the top of pedestal, outwards. The pedestal width is then assumed. This assumption comes from the absence of a pedestal transport model, that coupled with an MHD stability model could calculate consistently the pedestal width. The choice of $\rho_t = 0.9$ as top of pedestal position has been made to bound smoothly core and edge models. Some tests have been done changing this position up to $\rho_t = 0.96$ and no big impact has been found on the results. Some pedestal transport models to estimate the pedestal width exist (e.g. EPED [39] and IMEP [40]) but it was hard to implement them in a fast flight simulator, because they need to compute MHD stability. A self-consistent calculation of the pedestal width is anyway planned as future work, because its prediction can be crucial for MHD stability in some ELM-free or reactor-relevant scenarios. It is known that in high-confinement mode (H-mode) the dynamic is usually mainly dominated by edge localized modes (ELMs) [41, 42], which are MHD instabilities, while micro-instabilities are often reduced by sheared flows (as ExB shear) [43], even though some residual turbulent transport is still present, as was discussed previously. In this configuration, an ELM average model has been adopted and the diffusivities have been assumed to lay on the marginal stability limit of the MHD peeling–ballooning model [44, 45]. This approach is justified by the fact that it is behind the goal of this project to accurately reproduce the temporal evolution of the instability, as could be done by a high fidelity code or any discrete ELMs model. Furthermore, the temporal scale of these instabilities is much shorter than the duration of the discharge or the reaction time of the control system, so one can assume an average behaviour without leading to big limitations of the model when it is used within a flight simulator. Any discrete ELMs model could be added in the future. By assuming a marginal stability on the MHD peeling–ballooning boundary, a heuristic formula has been derived, which is based on a power law of the ratio between $\beta_{p,\mathrm{top}}$ at the top of pedestal and $\beta_{p,\mathrm{MHD}}$, which is the critical value for the onset of the instability

$$\begin{align} \chi_e=\left(\frac{\beta_{p,\mathrm{top}}}{\beta_{p,\mathrm{MHD}}}\right)^4 .\end{align} \tag{ 12 }$$

The critical value of β_p is

$$\begin{align} \beta_{p,\mathrm{MHD}}=0.686\sqrt{k}(1+\delta)^{1.68}{q}^{1.61}{\beta_{p,\mathrm{top}}}^{0.33}\hat{n}_e^{0.06} w_p^{1.29}, \end{align} \tag{ 13 }$$

which is a scaling that has been obtained from an EPED [39] database, according to [46]. Here, $\hat{n}_e$ is the electron density normalized to the Greenwald limit, while w_p is the pedestal width in the normalized poloidal flux label, which is fixed to 0.1. This is a good assumption for AUG, but it needs to be generalized in order to extrapolate to other devices. This scaling also includes the effect of the shape, which allows us to take into account and also predict the purely peeling or peeling-ballooning limited cases.

$\chi_i = \chi_e+\chi_{i,\mathrm{nc}}$ has been assumed in this case, supported by [40]. The particle diffusivity has been fixed to be equal to $C\cdot \chi_e$, where C has been fixed to 0.03 to match experimental stationary profiles in some discharges. This is still in agreement with [40].

In the low-confinement regime (or L-mode, i.e. when the power threshold is not overcome), the fitting model used for the core has simply been extended to this region, while particle diffusivity is kept equal to $F\cdot \chi_e$, where F is calibrated to match experiments in the database (F = 0.1).

In order to predict the transition between L- and H-mode, criteria based on the ion power crossing the separatrix has been chosen, according to the Schmidtmayr scaling [47]. This scaling has been previously tested successfully on different H-mode discharges. Differences in pedestal widths between density and temperature [48, 49] have been neglected.

The simulation of ELM-free scenarios (e.g. I-mode, EDA H-mode or negative triangularity) can be more challenging, due to the complicated and wide physics involved in the edge. These scenarios are very attractive for future reactors, therefore at least some reduced prediction capability should be included. Then, for these configurations some ad-hoc models are present, which have a less physics-based nature, due to the fact that the dynamic involved is still not fully understood. For the transition from L- to I-mode a scaling derived in [50] has been used. During the transition, a pedestal arises only for the temperatures, while the density remains equal. For the negative triangularity discharges, a correction on the L–H transition has been numerically introduced to modify the power threshold.

A model of the SOL needs to be included to create a fully integrated model. Its role is to furnish boundary conditions at the separatrix for the confined region. An important task is to reduce the amount of experimental input to strengthen the prediction capability. However, a minimum number of experimental fit parameters, which are often machine-related, is necessary due to the more complex physics and the weaker theoretical description of this region. Moreover, the dynamic in the SOL depends on the systems facing the SOL and some parameters can be machine-dependent. Two different models have been implemented to provide electron density and temperature at the separatrix.

4.1. Power exhaust

Several aspects should be taken into account concerning the power exhaust modeling, because this indirectly influences the fusion performances and fixes some technological constraints for large-size future machines. The two-point model has been largely validated across the years [51–56]. Its 0D structure based on a power balance between the outer midplane (OMP) and the divertor target makes it adequate for use in our work. Obviously, its implementation leads to a lack of detail about directional phenomena, radiation and more complicated aspects, but still represents a good candidate to obtain the necessary elements for a fully integrated code, thanks to its physics-based nature.

A simple formula to determine the electron temperature at the separatrix has been derived in [54], which can be used in the attached condition:

$$\begin{align} T_{e,\mathrm{sep}}=\biggl(\frac{7}{2}\frac{q_\parallel l_\parallel^{*}\pi q_\mathrm{cyl}R}{k_z k_0}\biggr) .\end{align} \tag{ 14 }$$

Here, $q_\parallel$ is the parallel electron heat flux, $l^*_\parallel$ is a dimensionless factor used to account for extended field-line lengths in alternative configurations, k₀ is the parallel electronic conductivity for Z = 1 divided by $T_e^{\frac{5}{2}}$, $q_\mathrm{cyl}$ is a safety factor defined as in [54], R is the major radius, and $k_z = \frac{1}{2}(0.672+0.076\cdot Z_\mathrm{eff}^{\,0.5}+0.252\cdot Z_\mathrm{eff})^{-1}$. $T_{i,\mathrm{sep}}$ is assumed equal to $T_{e,\mathrm{sep}}$.

The main limit of this model is that it does not include a consistent radiation model and it is derived for the attached condition, i.e. without pressure losses along the magnetic field lines. However, detachment modeling has to be included in the future because it is a necessary condition that has to be reached by machines of the next generation in order to prevent exceeding maximum heating loads and damage to the plasma facing components (PFCs). The modularity of the integrated model allows the future updating of the formula under detachment conditions. The formula shown needs $Z_\mathrm{eff}$ at the separatrix, which is not provided by the model for the confined region. This 'external' parameter is provided by the particle balance model in the SOL, which gives the value of the density at the separatrix for all the ion species. This model is discussed in the next section.

4.2. Particle balance

The proper treatment of particle diffusion in the SOL is a complicated task that can be pursued by the use of time-consuming codes, which include neutral treatment, such as SOLPS-ITER [57, 58]. This is beyond the scope of this work. Nevertheless, a simplified model can be used to translate engineering parameters into boundary conditions for the separatrix. Some formulas have already been used in the past for AUG, for example, in [59]. However, these are often related to the technological aspects of a specific experiment, so they miss general validity. Moreover, they are not time-dependent, which makes them inadequate for use in a flight simulator. This is the reason why a model that takes into account the minimum possible amount of experimental observations and more detailed geometrical aspects, has to be developed. This can be done by using a multispecies 0D particle balance that takes into account only the main features related to different regions included in the SOL, but ignores any directional detailed dynamic. For this purpose, in the model developed, the SOL has been split into six different zones that communicate with each other. These are the two divertor legs, the PF, the two upper sides of the tokamak cross section, and an additional zone close to the wall on the LFS. The splitting of the SOL region into the different zones is shown in figure 9. Such splitting can possibly be modified according to different divertor configurations by manipulating the regions (e.g. including another upper private flux region (PFR) in the case of a double null configuration, i.e. DN). Every region has different aspects to be taken into account, even though the main physics is the same. This allows a modular treatment of every zone, without losing the generality of the background structure. The model consists of a particle balance between all the confining regions, which includes diffusive terms, ionization sources, recombination sinks and sources of plasma coming from LCFS. The recombination sinks are assumed to be equal to zero, which is consistent with the assumption of the attached divertor. Gas puff and vacuum pump are also included as source and sink in the regions where they work. These are examples of the possibility of treating local effects in the respective zones. One other example is the recycling factor, which is only modeled in the wall zone (i.e. low recycling regime), even though a more precise description of it, including both the divertor legs, is planned for the future. This parameter is equal to 0.99 for D and H. Moreover, a diffusivity between this last region and outside the tokamak has been assumed (leakage), while for all the other zones, perfect reflectivity of the PFCs has been considered. Some coefficients have been added to the equations to keep concentration gradients between confining regions unbalanced, in order to mimic convection. These are called enrichment factors and they multiply the density of species:

$$\begin{align} \frac{\partial N_k}{\partial t}=S_k+P_k+\sum\limits_{j=1,j\neq k }^6 D_{jk}\left(\epsilon_{jk}\frac{N_j}{V_j}-\frac{N_k}{V_k}\right) .\end{align} \tag{ 15 }$$

Zoom In Zoom Out Reset image size

In the previous equation S_k represents the source term, which also contains the plasma coming from the confined region through the separatrix, P_k is the sink, D_jk is a 0D diffusion coefficient defined as an exchange rate between the regions j and k (in m s⁻³), N is the number of particles, V is the volume and _jk is the enrichment factor. From these balance equations, the temporal evolution of particle density in every region can be found. The average between the densities of the two upper parts of the SOL (SOLL and SOLR) is assumed to be equal to density at the separatrix, according to the 1.5D structure of ASTRA. In the model, a multispecies treatment is included up to eight species, but this number is arbitrary because it does not greatly impact the computational time, due to the analytical nature of the equations. One can possibly build neutral modeling in a similar way, by considering recombinations as sources and ionizations as sinks. This model is necessary to simulate detachment conditions. Considering six equations like (15) for six regions, the following system can be derived for each species:

$$\begin{align} \overline{N}^{\,j}=\overline{N}^{\,j}_0+\overline{\overline{D}}\cdot\overline{N}^{\,j}, \end{align} \tag{ 16 }$$

where $\overline{N}^{\,j}$ is the vector of the particle number for every region, j index indicates a specific species, $\overline{N}^{\,j}_0$ is the vector that includes all the local effects (ionizations, recombinations, recycling, vacuum pump or gas puff) and $\overline{\overline{D}}$ is the matrix of the diffusivities normalized by some factors, which also include the enrichment factors. The system is calculated implicitly.

This model allows us to calculate the density at the separatrix with gas puff and vacuum pump velocity as inputs, reducing so the necessary experimental data to the engineering parameters, which are controllable and needed anyway for the discharge plan. However, a crucial problem is the determination of the diffusivities between confining zones. These have been derived heuristically, considering the main contributions to the parallel and perpendicular particle transport:

$$\begin{align} & D_\parallel=0.1\frac{M}{L_\mathrm{par,sep}} \quad \left(\frac{m^4}{s}\right), \quad D_\perp=\frac{0.05\nu_\mathrm{sep}}{(R_{w}-R_\mathrm{omp})^2} \quad \left(\frac{m^5}{s}\right), \nonumber\\ & \qquad D_{2,3}=D_{3,4}=D_{2,6}=1000 D_\parallel, \nonumber\\ & \quad D_{1,0}=0.03 D_\perp,\qquad D_{1,2}=0.5 D_\perp \nonumber\\ & \qquad D_{4,5}=D_{5,6}=0.1 \quad \left(\frac{m^3}{s}\right),\quad D_{1,6}=0 \quad \left(\frac{m^3}{s}\right) .\end{align} \tag{ 17 }$$

Here, the 'sep' subscript stands for separatrix, M is the sound velocity, calculated as $\sqrt{\frac{T_{i,\mathrm{sep}}}{m_\mathrm{D}}}$, where $m_\mathrm{D}$ is the deuterium density, $L_\mathrm{par,sep}$ is an estimate of the distance traveled by the particles along the field lines and it is calculated as $\pi q_\mathrm{sep}R$, R_w and $R_\mathrm{omp}$ are the major radius at the wall and OMP positions. $D_{2,3}$, $D_{3,4}$ and $D_{2,6}$ were calculated assuming convection as the dominating parallel transport. $D_{4,5}$ and $D_{5,6}$ were fixed to 0.1 to mimic a lower transport between the divertor legs across the PFR with respect to the parallel transport along the field lines, which is a reasonable assumption for ions [60] in L-mode. The high gradients between divertor legs and PFR ensures a higher value of diffusion across the PFR in H-mode [60]. $D_{1,6}$ has been fixed to 0 to neglect the transport from the wall directly to the divertor (i.e. perpendicular transport is the only dynamic assumed in the thin layer of plasma near the wall). $D_{1,0}$ and $D_{1,2}$, which represent the perpendicular transport in the LFS, have been assumed proportional to the collisionality, according to [61], neglecting poloidal asymmetries. In particular, $D_{1,2}$ has been multiplied by 0.5 which is the standard value used in SOLPS calculations [62, 63], while $D_{1,0}$ has been multiplied by 0.03 to take into account a certain level of D wall retention, which counteracts the outwards diffusion. The previous equations contain some numerical factors which represent the calibration of the model to match an experimental case.

These coefficients do not take into account the drifts, which have been modeled by employing some enrichment factors, that keep the densities of confining regions unbalanced and ensure a certain background level of diffusion. In this way, the compression factors can be taken into account. The enrichment factors are

$$\begin{align} \begin{aligned} & \epsilon_{1,}=\epsilon_{23}=\epsilon_{45}=\epsilon_{16}=\epsilon_{56}=1, \qquad \epsilon_{34}=10, \\ & \qquad \epsilon_{26}=5\cdot \mathrm{max}\left\{1,\mathrm{min}\left\{20,0.2{n_\mathrm{D}}^{0.67}\right\}\right\} \end{aligned} \end{align} \tag{ 18 }$$

in favourable configuration (i.e. $B_tI_p\gt 0$), while in unfavorable configuration $\epsilon_{23} = \frac{10}{R^2(1+2r R^{-1}))}$. This last coefficient has been used to model a $\nabla B\times B$ drift correction, which is assumed to be inversely proportional to the size of the machine and includes an aspect ratio dependence. ₃₄ has been fixed to 10 to modulate a certain background level of HFS high density front [62, 64, 65]. ₂₆ is a derivation of the steady-state formula of [59], which has also been successfully implemented in IMEP [40], for a time-dependent situation. $n_\mathrm{D}$ is the density of D or H at the divertor. Also, here the numerical factors are calibrated to match the separatrix electron density of an H-mode AUG discharge in the flight simulator. A better formulation of the enrichment factors including detachment, X-point radiator [66], $\nabla B\times B$ and E × B drifts is planned for the future. All the particle diffusivities and enrichment factors shown here are used for D or H, but similar formulae are also applied for impurities, taking into account some technological aspects, like W sputtering, B coming from the boronization of the machine and N seeding. In particular, the W sputtering is roughly modeled by adding a simple constant input of this species in the SOL. A similar approach is used for the residual B coming from the boronization.

The temporal evolution of D ions in a Fenix simulation of the H-mode #40446 discharge is shown in figure 10, together with the D gas puff. In this figure, one can see that the densities above the X-point are one order of magnitude lower than the regions below. This is consistent with many experimental observations and theoretical predictions, especially if the source coming from the divertor region (e.g. by ionizations) is bigger than the one from the OMP. One can also notice that the PFR density is lower than those of the divertors for $t\lt 2\mathrm{s}$, while it lays between the inner and outer divertor values for $t\gt 2\mathrm{s}$. These two temporal ranges are before and after the onset of H-mode, and it is consistent with observations from [60]. It is worth mentioning that the only parameters that depend on the temperature are the ionization coefficients of H, D and He. It is still complicated to run simulations with both power exhaust and particle balance models because they depend on each other's output.

Zoom In Zoom Out Reset image size

The main limit of the model is the presence of a specific divertor configuration, which is lower single null (LSN). As already mentioned, other configurations (e.g. upper single null, i.e. USN, and DN) can still be included by splitting the SOL into more zones, while maintaining a similar approach in describing the particle content in each zone. Anyway, the model has been validated only against LSN configurations, thus further validation with other configurations is needed. Another limitation is the modeling of the X-point radiator and detachment, which are still not included, but they could be implemented by adding the treatment of neutrals in the particle balance and modifying the enrichment factors, without changing the background structure. Finally, there is no description of transport and heat and particle loads deriving from ELMs. This is really important to prevent any damage to the divertor, but should not have much affect on the evolution of the global parameters, which are of primary importance in a flight simulator. Moreover, the W influx should be independent of the conditions of the divertor, even in the presence of ELMs, because the region which rules such quantity in AUG is the main chamber [67]. Then the same model with the same tuning coefficients is used for both L- and H-mode. We are aware that such limitations reduce strongly the physical details and range of applicability of the model, but they still represent a good starting point for future developments. Moreover, the amount of predicted operation-relevant aspects (e.g. divertor heat loads) will be increased.

The first fully-integrated modeled discharge has been tested inside the Fenix flight simulator, including core, edge and SOL models. For the neoclassical transport, some coefficients from the literature [68–70] have been used. From $\rho_t = 0$ to $\rho_t = 0.9$, the fitting model has been used for $\chi_{i,\mathrm{an}}$ and $\chi_{e,\mathrm{an}}$, while the particle diffusivity has been assumed to be equal to $\chi_{e,\mathrm{an}}$. The particle pinch has been modeled with the heuristic formula mentioned in section 2, which takes into account the main parameters involved. From $\rho_t = 0.9$ to $\rho_t = 1$, the model used depends on the specific mode of the plasma: in L-mode the fitting core model is extended to this region because the nature of the transport is the same; in H-mode the ELMs averaged model is used; in I-mode an H-mode-like model is used for heat diffusivities and an L-mode-like one is used for particle diffusivity. For this regime, a power threshold criterion similar to the Schmidtmayr scaling has been adopted. The discharge analyzed (#40446) is a standard H-mode, so the $\chi_{e,\mathrm{an}}$ and $\chi_{i,\mathrm{an}}$ are calculated by a temporal combination of the various models, which depends on the temporal evolution of the regime of the plasma, which is predicted by the ion power threshold criteria (Schmidtmayr scaling). The particle diffusivity in the edge has been fixed to $C\cdot\chi_{e,\mathrm{an}}$, where C varies between 0.1 in L-mode and 0.03 in H-mode during the simulation. This suffers from a lack of predictive power of the particle transport in the edge in L-mode, but it was necessary to validate the heat diffusivities model. An improvement inthe formula used for this phase is planned for the future.

The trajectories of some engineering parameters and actuators of the discharge #40446 are shown in figures 11–16. The feedback mechanisms in Fenix for this discharge act on the average density, the position and the shape. The feedback on β_p has been switched off and feedforward is used in order to validate the transport model with the experimental NBI time trace. In figure 11, gyr1, gyr2, gyr3 and gyr4 panels show the ECRH gyrotrons with not null time traces, while ICRH pairs 1 and 2 show the couples of ICRH sources. These heating systems do not actuate control on any parameter. In figure 12, one can see instead the time trace of the NBI sources, which is generally used to control β_p . In this simulation, the β feedback was switched off to validate the transport model. The shape of the plasma and its position are controlled by the magnetic coils. In figure 14 the time traces of such coil currents are shown, while the position of the plasma (i.e. some geometrical parameters) is shown in figure 15. In figure 13 the geometry of the magnetic coils with the associated nomenclature is shown, to help the reader identify them in figure 14. More details about the connection between the current coils, the position and the shape are given in the legends of the figure. The feedback on the average density is shown in figure 16, where the average density is shown in the upper plot and the divertor gas puff is shown in the lower plot. One can see in the upper plot a sensitive difference between the actual value of the density and its target value during the feedback phase. This is due to the fact that the inward particle transport is not strong enough to increase the inner core density on the right timescale, despite the strong gas puff shown in the lower plot. This is due to the inaccurate particle transport model, whose limitations are more evident during the ramp-up (i.e. in the feedback phase for this discharge). Information about the control system implemented in Fenix and its feedback algorithms can be found in [71, 72]. Figure 17 shows the time trace of some global parameters, while in figure 18 some time traces related to the integrated model are shown. These are the electron density and temperature at the separatrix (respectively figures 18(a) and (b)), the ion power crossing the separatrix (figure 18(c)) and the comparison between β_p at the top of the pedestal and its critical value calculated with the scaling (figure 18(d)), as in section 3. One can see that when the power at the separatrix exceeds the critical value (around 2 s) $\beta_{p,\mathrm{top}}$ approaches the critical value. In figure 17(c) one can see that β_p follows the experimental trajectory for most of the duration of the discharge, so a first validation of the model is given, considering also that the profiles have been matched. By looking at β_p , one can see that the early phase has not been matched. This is due to the fact that during the ramp-up a TEM-dominated regime at high collisionality is found [73], and our model is not able to reproduce it. Moreover, the core particle transport model underestimates the density in the ramp-up phase, because the predicted inwards diffusion and its timescale are too small.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In figure 19 the kinetic profiles are shown for different phases of the discharge. Here, one can see that the profiles from the simulation are not within the experimental error bars, but for ρ < 0.9 the maximum error is $20\%$ and for ρ > 0.9 it is $7\%$. These errors show that the transport model can be improved, but we consider it a good starting point, considering that the task of the flight simulator is not a detailed analysis of the physics, but the prediction of the evolution of some global parameters. Considering that the discharge #40446 is an H-mode where the transition takes place at the beginning of the flat-top, one could argue that the mismatch in β_p found in the first 1.5 s is due to the fact that during this time range the plasma is in L-mode. In order to disprove this and disentangle the ramp-up from the L-mode, an additive Ohmic heated L-mode discharge (#38998) has been run in Fenix. The time traces of this discharge are shown in figure 20. Figure 21 shows the kinetic profiles between 2.4 and 2.7 s. Here on the right the density profile from Fenix (solid line) shows a non-negligible mismatch with the experimental one (dots). Such a disagreement deserves justification. The particle transport model used in the core is poor. An attempt to develop a more detailed and physics-based formula for the particle pinch matching TGLF in the simulations has been done, but it was not successful, therefore, the peaking factor cannot always be guaranteed. Moreover, this discharge is an L-mode, which is harder to predict, due to the absence of a pedestal model in the edge. It is worth mentioning that there is no inconsistency between the line averaged density time trace and the underestimation of the density profile, because there is no density feedback in this discharge.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A simulation including the match of the ramp-up phase is planned for the future. Finally, regarding the discharge #40446, in figure 22 the different contributions to χ_e via the ETG, ITG and TEM channels are shown along the entire discharge. One can see that, for this case, ITG is the dominant mode. In order to strengthen the validation of the integrated transport model, another four H-mode discharges with different particle sources and heating power time traces have been run. The trajectories of the plasma current and β_p for the various discharges are shown in figure 23. Moreover, the time traces of the heating power and the D particle source for these discharges are compared respectively in figures 24(a) and (b), while the time trace of the average density of the core (for ρ < 0.3) is shown in figure 25. In figure 24(b) there is no gas puff signal for discharge #39977, because there was no diagnostic data available. One can notice that in figures 25(a), (c) and (d) there is an underestimation of the average density in the inner core (for ρ < 0.3). This is attributable to the particle transport model, which is not able to reproduce the peaking factor of the density. Moreover, the slow timescale of the simulated density reaction to the actuators does not allow us to reproduce a strong increase in the density in short times, as in figures 25(c) and (d). The feedback on the line averaged density for these discharges is switched off at around 1.5 s. All these results suggest that the model can predict the temperature profiles and the global parameter evolution, while some improvement is necessary to model the particle transport, especially in the presence of a strong concentrated increase in density and its peaking factor. Nevertheless, the main physics during the flat-top + ramp-down phase has been predicted, despite the differences between the various transients [73, 74]. This suggests that the transport model is appropriate for use in the usual AUG conditions in the Fenix suite.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A set of analytical formulae that describe the transport in different plasma regions has been derived from theory-based models. For the core, some formulae from the literature have been fitted to TGLF+ASTRA simulations, to model temperature gradient-driven micro-instabilities. In the edge, an averaged ELM model has been used for the H-mode, coupled with Schmidtmayr scaling for the L-H transition, while in L-mode the core model has been extended to this region. For the SOL, two different models have been built to compute species densities and temperatures at the separatrix, respectively, based on a particle balance and the two-point model. Six discharges, five H-mode (#40446, #40009, #40254, #39977, #39967) and an L-mode (#38898), have been successfully simulated with the fully-integrated transport model in the Fenix flight simulator framework, and a match of the experimental trajectories during the flat-top and the ramp-down has been found, proving the ability to predict plasma energy and kinetic profile evolution with few experimental inputs. However, the ability to predict the peaking factor of the density profile is limited, thus some additional improvements to the particle transport must be investigated in the future. The analytical physics-based nature of the models allows us to use them in a realistic and fast way as a predictive tool. Their modularity opens up the possibility of constant updating by including new physical ingredients (e.g. detachment). The next step is the validation of the model against multiple AUG discharges in a wide range of configurations, and the extrapolation of it to machines of different sizes (e.g. ITER and DEMO). Long-term development also includes a physics-based modeling of the particle transport and the inclusion of a detachment model.

All data that support the findings of this study are included within the article (and any supplementary files).

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.