Integrated tokamak modelling with the fast-ion Fokker–Planck solver adapted for transient analyses

Upcoming tokamak devices such as JT-60SA and ITER have to sustain advanced plasmas, which have a high non-inductive current ratio, good confinement characteristics and sufficient magneto-hydrodynamic stabilities. To achieve advanced plasmas, it is essential to control the plasma pressure and the current profile using auxiliary heating and current drive systems. To study the advanced plasmas, integrated modelling codes such as TOPICS [1], CRONOS [2], and ASTRA [3] have been developed, which comprise of a one-dimensional (1D) transport solver and a two-dimensional (2D) equilibrium solver with several heat and current source models. One of the main purposes of integrated modelling is to explore high-performance steady plasma, and establish a robust approach to such plasma. We have extended the capability of TOPICS in the transient phase, such as the current ramp-up or ramp-down phases, through consistent linkage between the transport and free-boundary equilibrium solvers [4].

The current ramp up with reduced central solenoid (CS) flux consumption in JT-60SA was investigated using TOPICS [5]. The plasma current could be ramped up without CS flux consumption and without the ideal magnetohydrodynamic instabilities. The current drive due to the fast ions being produced by a neutral beam (NB) played a key role in the ramp up, but Alfvén instabilities driven by the fast ions might have affected the driven current through fast-ion redistribution. To confirm the accuracy of the obtained results and to model the additional effects of the instabilities, the fast-ion modelling needs to be improved, as explained below. TOPICS offers two main options to evaluate the heating and the current drive produced by fast ions. One solves a 1D Fokker–Planck equation in the energy direction to evaluate heating and applies an analytic model [6] to the current drive, and the other uses the orbit following the Monte Carlo code OFMC [7]. The former was verified in ITPA activity [8] and was used in the above mentioned JT-60SA simulations. Both the former current drive model and the latter OFMC were originally embedded for steady analysis. For the transient analysis, an ad hoc time-dependent treatment based on fast-ion power deposition was used; however, a steady-state solution for the fast-ion distribution function was used in the former model. In addition, the analytic solution has little extendability for further development, e.g. introducing fast-ion redistribution by instabilities. In contrast, the latter is too time-consuming from the perspective of analysing a long-period transient phase, while it naturally considers the effects of magnetic ripples and finite orbit width (FOW).

To enable a self-consistent evolution among the fast ions, the bulk plasma, and the equilibrium magnetic fields, we employs a bounce-averaged Fokker–Planck solver (BAFP), which rapidly solves the fast-ion distribution function in a 2D velocity space, and also integrated the BAFP with TOPICS. Recently, with almost the same motivation, a rapid fast-ion Fokker–Planck solver RISK for CRONOS and the European integrates modelling framework has been reported [9]. The RISK and BAFP codes obtained essentially equivalent solutions because they are based on the same BAFP equation, but there are some differences. The RISK code solves the equation for NB-injected ions by expanding the distribution function in the eigenfunctions of the collisional pitch angle scattering operator. FOW effects are partially considered using an ad hoc orbit broadening algorithm [9]. The BAFP code directly solves the equation for any type of fast ions, NB-injected ions, alpha particles, and other fusion products in two dimensions, namely, particle speed and pitch angle. Although FOW effects are not yet considered, the electric acceleration term is retained in the BAFP equation. The effect of the parallel electric field on ions is usually neglected because it is much smaller than that on electrons based on their mass ratios. Fast ions, however, circulate around the toroidal loop at great speeds and pick up the loop voltage multiple times compared with bulk ions. Therefore, the electric acceleration of fast ions is worthy of evaluation for consistent simulation with bulk plasma and equilibrium. A certain level of electric effect on the fast ions was observed and will be discussed later with the integrated simulation results.

The paper is organised as follows. In section 2, we briefly introduce the BAFP code and its integration with TOPICS. Integrated simulations applied to JT-60SA and ITER plasmas are presented in section 3. Finally, in section 4, we summarise the work and discuss its outlook.

The evolution of the fast-ion distribution function in a state of thermodynamic non-equilibrium in high temperature plasma is governed by the Fokker–Planck equation. The equation generally describes the evolution of the distribution function in 6D phase space. To reduce the computation time for facilitating the transient analysis, three dimensions are reduced via toroidal symmetry, gyro averaging, and bounce averaging; then, we have the bounce-averaged Fokker–Planck equation

$$\begin{eqnarray}&&\frac{\partial \mathcal{F}}{\partial t}=-{{\langle \nabla \cdot \mathbf{\Gamma}\rangle}_{\text{B}}}+{{\langle S\rangle}_{\text{B}}}-{{\langle L\rangle}_{\text{B}}},\end{eqnarray} \tag{ 1 }$$

where $\mathcal{F}$ is the bounce-averaged distribution function; $\mathbf{\Gamma}$ is the particle flux in phase space; S is the particle source that stems from ionisation and charge-exchange for NB species and fusion reactions for α particles; and L is the thermalisation loss of the fast ions which assimilates into the bulk ions. The bounce average operator is defined as

$$\begin{eqnarray}&&{{\langle \ldots \rangle}_{\text{B}}}\equiv \frac{1}{{{\tau}_{\text{B}}}}\mathop{\oint}^{}(\ldots )\frac{\text{d}s}{{{v}_{\parallel}}},\end{eqnarray} \tag{ 2 }$$

where ${{\tau}_{\text{B}}}=\mathop{\oint}^{}\text{d}s/{{v}_{\parallel}}$ is the bounce period; ${{v}_{\parallel}}$ is the speed parallel to the magnetic field; and $\text{d}s$ is the line element along the field line. An integral path is taken along the particle orbit, following the circuit of the poloidal cross-section for the passing particles, and between the two turning points for the trapped particles. Drift motion is neglected in the BAFP. Hence, during a bounce period, staying in same radial position, particles do not feel the radial electric field. The magnetic flux label ρ and the particle speed v are virtually constant, while the pitch angle η of the particle velocity to the magnetic field line varies with the bounce motion: ${{v}_{\parallel}}=v\,\text{cos}\,\eta$ oscillates due to the magnetic well.

In general, the divergence of particle flux in the phase space is written as

$$\begin{eqnarray}&&\nabla \cdot \mathbf{\Gamma}=\frac{1}{\mathcal{J}}\frac{\partial}{\partial {{x}^{i}}}\mathcal{J}\frac{\partial {{x}^{i}}}{\partial {{y}^{j}}}\left({{\mathcal{A}}^{j}}-{{\mathcal{D}}^{jk}}\frac{\partial {{x}^{l}}}{\partial {{y}^{k}}}\frac{\partial}{\partial {{x}^{l}}}\right)\mathcal{F},\end{eqnarray} \tag{ 3 }$$

where ${{\mathcal{A}}^{j}}$ and ${{\mathcal{D}}^{jk}}$ are the advective and diffusive coefficients, respectively; $\mathcal{J}$ is the Jacobian of the transformation from the Cartesian coordinates to the x coordinates; and the summation over repeated indices is understood. To facilitate bounce-averaging, we choose $x\equiv \left({{\rho}_{0}},\theta ,\phi ,{{v}_{0}},{{\eta}_{0}},\zeta \right)$, where θ is the poloidal angle; ϕ is the toroidal angle, ignorable by axi-symmetry; ζ is the gyroangle, irrelevant in the drift-kinetic regime; and the subscript 0 denotes a value at the minimum magnetic field intensity B₀ on a magnetic surface. As mentioned above, $\rho ={{\rho}_{0}}$ and v   =   v₀ during a bounce period. In addition, by adiabatic invariance of the magnetic moment, $\eta =\text{si}{{\text{n}}^{-1}}\left(\sqrt{B/{{B}_{0}}}\text{sin}\,{{\eta}_{0}}\right)$, where B is the local magnetic field intensity. Other coordinates, y, are reserved for describing the local collisional flux and the parallel electric acceleration. A choice of $y=(\rho ,\theta ,\phi ,v,\eta ,\zeta )$ is suitable for describing the Fokker–Planck collisional flux. To couple with a 1D fluid bulk transport solver, we adopt the Fokker–Planck coefficients in case of collision with the Maxwellian background plasma:

$$\begin{eqnarray}&&{{\mathcal{A}}^{v}}=-{{m}_{a}}\underset{b}{\mathop \sum}\,\frac{{{\gamma}_{ab}}U}{{{T}_{b}}} \\ &&{{\mathcal{D}}^{vv}}=\frac{1}{v}\underset{b}{\mathop \sum}\,{{\gamma}_{ab}}U \\ &&{{\mathcal{D}}^{\eta \eta}}=\frac{1}{2{{v}^{3}}}\underset{b}{\mathop \sum}\,{{\gamma}_{ab}}\left(\text{erf}(u)-U\right),\end{eqnarray} \tag{ 4 }$$

with $U=\left(\text{erf}(u)-u\;\text{er}{{\text{f}}^{\prime}}(u)\right)/2{{u}^{2}}$, $u=v\sqrt{{{m}_{b}}/2{{T}_{b}}}$, and ${{\gamma}_{ab}}=\left(e_{a}^{2}e_{b}^{2}/4\pi \epsilon _{0}^{2}m_{a}^{2}\right){{n}_{b}}\ln \Lambda$, where $\text{erf}(u)$ is the error function; $\text{er}{{\text{f}}^{\prime}}(u)$ is its derivative; e is the charge; m is the mass; n is the density; T is the temperature; the subscripts a and b denote the fast ion species and the background plasma species, respectively; and $\ln \Lambda$ is the Coulomb logarithm. By choosing a local coordinate as $y=\left(\rho ,\theta ,\phi ,{{v}_{\parallel}},{{v}_{\bot}},\zeta \right)$, the parallel electric acceleration is expressed in a advective coefficient:

$$\begin{eqnarray}&&{{\mathcal{A}}^{{{v}_{\parallel}}}}=\frac{e}{m}{{E}_{\parallel}},\end{eqnarray} \tag{ 5 }$$

where ${{E}_{\parallel}}$ is the parallel electric field. It should be remarked that the particle flux appears only in a velocity space under our conditions for two reasons: both the collision and the electric field locally changes the particle velocity only, and the projection $\partial {{x}^{i}}/\partial {{y}^{j}}$ contains no cross-spatial component between the real coordinates and the velocity coordinates. If drift motion was reverted, canonical toroidal angular momentum ${{p}_{\varphi}}$ should be chosen rather than the magnetic flux label ${{\rho}_{0}}$ in x coordinates. In that case, because ${{p}_{\varphi}}$ is a real-velocity mixed variable, for example, $\partial {{p}_{\varphi}}/\partial \eta \ne 0$ leads to particle flux in real space as well.

Bounce-averaging the divergence of the particle flux (3) with the coefficients (4) and (5) respectively, we have explicit expressions of the bounce-averaged term in (1):

$$\begin{eqnarray}&&{{\langle \nabla \cdot {{\mathbf{\Gamma}}_{\text{coll}}}\rangle}_{\text{B}}}=\frac{1}{\mathfrak{S}}\left(\frac{\partial}{\partial {{v}_{0}}}\mathfrak{S}\left({{\mathcal{A}}^{v}}-{{\mathcal{D}}^{vv}}\frac{\partial}{\partial {{v}_{0}}}\right)-\frac{\partial}{\partial {{\eta}_{0}}}\mathfrak{T}{{\mathcal{D}}^{\eta \eta}}\frac{\partial}{\partial {{\eta}_{0}}}\right)\mathcal{F},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{{\langle \nabla \cdot {{\mathbf{\Gamma}}_{\text{elec}}}\rangle}_{\text{B}}}=\frac{e}{m}\frac{\mathfrak{U}}{\mathfrak{S}}\left(\frac{\partial}{\partial {{v}_{0}}}\mathcal{V}\,\text{cos}\,{{\eta}_{0}}-\frac{\partial}{\partial {{\eta}_{0}}}\mathcal{V}\frac{\text{sin}{{\eta}_{0}}}{{{v}_{0}}}\right)\mathcal{F},\end{eqnarray} \tag{ 7 }$$

with $\mathfrak{S}=\mathcal{V}\mathop{\oint}^{}\frac{\text{cos}{{\eta}_{0}}}{\text{cos}\eta}\text{d}s$, $\mathfrak{T}=\mathcal{V}\mathop{\oint}^{}\frac{\text{cos}{{\eta}_{0}}}{\text{cos}\eta}\frac{\text{ta}{{\text{n}}^{2}}{{\eta}_{0}}}{\text{ta}{{\text{n}}^{2}}\eta}\text{d}s$, $\mathfrak{U}=\mathop{\oint}^{}{{E}_{\parallel}}\text{d}s$, and $\mathcal{V}=v_{0}^{2}\text{sin}{{\eta}_{0}}$. To correspond with arbitrary magnetic configuration, $\mathfrak{S}$ and $\mathfrak{T}$ are numerically integrated under a given equilibrium. The electric acceleration of trapped particles cancels out during a back-and-forth motion: $\mathfrak{U}=0$ for trapped particles. Passing particles pick up a poloidal loop voltage during a transit period: $\mathfrak{U}=\sigma q{{V}_{\text{loop}}}$, where σ is the sign representing the go-around direction; q is the safety factor; ${{V}_{\text{loop}}}$ is the toroidal loop voltage. When these field related coefficients $\mathfrak{S}$, $\mathfrak{T}$, and $\mathfrak{U}$, simultaneously with the plasma related coefficients ${{\mathcal{A}}^{v}}$, ${{\mathcal{D}}^{vv}}$, and ${{\mathcal{D}}^{\eta \eta}}$, are all given, the evolution of the distribution function $\mathcal{F}$ can be computed by discretising (6) and (7). Implementing the BAFP code, we consulted [10], wherein the bounce-averaged equation is derived and the numerical methods, the trapped and passing boundary handling, and other necessary matters to implement the code are described.

Figure 1 shows the flowchart of the TOPICS calculation combined with the BAFP solver. The time-dependent calculation starts with zero fast ions and the given initial profile of bulk plasma. There are three time steps corresponding to each objective to solve: $\Delta {{t}_{\text{FI}}}$ for the fast ions, $\Delta {{t}_{\text{TR}}}$ for the bulk transport, and $\Delta {{t}_{\text{EQ}}}$ for the equilibrium. $\Delta {{t}_{\text{FI}}}$ and $\Delta {{t}_{\text{TR}}}$ are automatically adjusted during the time-dependent calculation so that the time-derivative term $\partial /\partial t$ remains sufficiently small. The equilibrium solver, as the name suggests, involves no time-derivative term; $\Delta {{t}_{\text{EQ}}}$ merely represents the update interval, chosen small enough so as not to significantly vary the metric. Using the metric and bulk plasma profiles, the BAFP solver proceeds with the innermost iteration $\Delta {{t}_{\text{TR}}}/\Delta {{t}_{\text{FI}}}$ times. Subsequently, using the BAFP solution as sources, the transport solver proceeds with the middle iteration $\Delta {{t}_{\text{EQ}}}/\Delta {{t}_{\text{TR}}}$ times. As indicated in figure 1, the sources are updated at every time-step of the middle iteration. In this regard, $\Delta {{t}_{\text{TR}}}=\Delta {{t}_{\text{FI}}}$ in most time-steps at a quasi-steady state. Afterward, the equilibrium is solved consistent with the total pressure and the current profiles in the outermost iteration.

Zoom In Zoom Out Reset image size

The collisional flux ${{\mathbf{\Gamma}}_{\text{coll}}}$ in (6) is the sum of the fluxes caused by Coulomb collisions with the thermal ions and electrons, and it contains the slowing down, energy diffusion, and pitch-angle scattering processes, whose coefficients are given as ${{\mathcal{A}}^{v}}$, ${{\mathcal{D}}^{vv}}$, and ${{\mathcal{D}}^{\eta \eta}}$ in (4), respectively. They depend on the electron density ${{n}_{\text{e}}}$, the ion densities n_i, the electron temperature ${{T}_{\text{e}}}$, and the ion temperatures T_i. The transport solver of TOPICS obtains the profiles of these plasma parameters. The flux ${{\mathbf{\Gamma}}_{\text{coll}}}$ of fast ions along the energy direction in the phase space is inextricably associated with the heating power supplied to the electrons ${{P}_{\text{e}}}$ and the bulk ions P_i. The pitch angle component of the flux determines the ratio between the numbers of trapped and passing particles, and the fast-ion current ${{j}_{\text{FI}}}$.

The electric flux ${{\mathbf{\Gamma}}_{\text{elec}}}$ in (7) depends on q and ${{V}_{\text{loop}}}$ through $\mathfrak{U}$. Both q and ${{V}_{\text{loop}}}$ are obtained using the transport solver. The radial profile of ${{V}_{\text{loop}}}$ is provided in addition to the q profile rather than the scalar value at the surface. The equilibrium solver of TOPICS obtains the poloidal flux function ψ and the toroidal flux function F for calculating the field related coefficients $\mathfrak{S}$ and $\mathfrak{T}$, using the total pressure ${{p}_{\text{tot}}}$ and the safety factor q. The total pressure ${{p}_{\text{tot}}}$ is composed of bulk pressures and the fast-ion pressure ${{p}_{\text{FI}}}$. The safety factor q is obtained including the contribution of the external current ${{j}_{\text{ext}}}$. One of the external currents is the NB-driven current, which is the product of the fast-ion current ${{j}_{\text{FI}}}$ and the shielding factor [11]. The fast-ion current ${{j}_{\text{FI}}}$ and its pressure ${{p}_{\text{FI}}}$ are, respectively, obtained as the first and the second velocity moments of the BAFP solution $\mathcal{F}$. The zeroth moment, i.e. the density of the fast ions ${{n}_{\text{FI}}}$, is also obtained to maintain quasi-neutrality; however, ${{n}_{\text{FI}}}$ usually makes a negligible contribution to the total plasma density, even if ${{p}_{\text{FI}}}$ accounts for a considerable portion of ${{p}_{\text{tot}}}$. Details of the transport and equilibrium solvers of TOPICS are described in [4].

To verify the BAFP code and our implementation, we compared it with a calculation of the OFMC code [7]. The calculation was applied to an ITER-class plasma (plasma current ${{I}_{\text{p}}}=15$ MA, toroidal magnetic field density ${{B}_{\text{tor}}}=5.2$ T, ${{n}_{\text{e}}}=1.0\times {{10}^{20}}$ m⁻³) with the tangential NB injection having an accelerating energy of 1 MeV and a power of 33 MW. The effects of the ripple and the electric field were switched off in both calculations to ensure identical conditions. Figure 2 shows the steady-state radial profiles of the heating power and the fast-ion current caused by the NB. The BAFP calculation was in good agreement with the OFMC calculation especially in the case without the FOW effects. We suppressed the FOW effects setting drift velocity to zero; the fast ions deposited all their energy on the magnetic surface they ionised. Because both the plasma current and the toroidal magnetic field faced clockwise in this calculation, the actual orbits of the co-injected ions traced inside the magnetic surface of birth. This orbital effect leads to the slight inward shift of the power deposition observed in figure 2.

Zoom In Zoom Out Reset image size

The computational costs of these two methods are very different. These calculations were performed on a parallel computer comprising of Intel Xeon CPUs (X5680 or above). The OFMC calculation required 1 h with 80 cores using 20 000 test particles, while the BAFP calculation required 5 s with 10 cores. The BAFP code solved equation (1) for each magnetic flux surface in parallel with the multiple processing cores. In the steady-state calculation mode of the BAFP, the time-derivative term (left hand side of (1)) serves as a numerical convergence factor instead of its original meaning, and the convergence iteration is typically conducted five to ten times. In terms of calculation time, a convergence iteration in a steady-state calculation is equal to a temporal evolution in a transient calculation, which is completed within a second.

3.1. JT-60SA simulation

First, we applied the integrated code TOPICS with the BAFP to JT-60SA [12] plasma in the ramp-up phase. Operational conditions imposed on the simulation were as follows. Volume expansion and transition to the divertor configuration were achieved in 2.5 s. The plasma current was ramped up at the rate of ${{\dot{I}}_{\text{p}}}=0.4$ MA s⁻¹ until 13.0 s so that a flattop current of ${{I}_{\text{p}}}=5.5$ MA was attained. The average electron density was ramped up to $2.5\times {{10}^{19}}$ m⁻³ along the current ramp-up as shown in figure 3.

Zoom In Zoom Out Reset image size

A negative-ion-based NB (N–NB) with an accelerating energy of 500 keV and a power of 5 MW was injected at 5.0 s. After the current flattop started at 13.0 s, the plasma was still in the transient phase owing to continuous N–NB injection and current diffusion. Although the JT-60SA has various external heating and current drive systems with a total power of 41 MW, only the single N–NB was injected in this calculation to clearly observe the effects of the heat and current sources provided by the BAFP. The NB-absorbed power ${{P}_{\text{NB}}}$ was dominant in the total heating power ${{P}_{\text{tot}}}$. According to the scaling law in [13], ${{P}_{\text{tot}}}$ exceeds the L–H transition power soon after the NB injection; for simplicity we did not impose any model of pedestal formation in this calculation. The ohmic heating power ${{P}_{\text{OH}}}$ decreased after the NB injection owing to a decrease in resistivity with an increase in ${{T}_{\text{e}}}$. The density profiles were prescribed as a function of ${{n}_{\text{e}}}(\rho )={{n}_{\text{e},\text{ax}}}\sqrt{1-{{\rho}^{2}}}$ with the time evolution of electron density at the axis ${{n}_{\text{e},\text{ax}}}$. Meanwhile, the ${{T}_{\text{e}}}$, T_i and q profiles generated by the transport solver were consistent with the heating power and the driven current determined using the BAFP solution. In the transport solver, the thermal transport and magnetic diffusion equations were solved using CDBM [14, 15], a theory-based turbulent transport model, and the matrix-inversion method [16], for calculating the neoclassical properties of plasma, as in [4]. The CDBM model has been compared with other transport models and several discharges in JT-60U and JET. The CDBM has given results close to the experimental data of normal and week shear discharges [17]. Although not yet published, it has been confirmed that the CDBM results are also in good agreement with the experimental data of reversed shear discharges with an internal transport barrier.

The N–NB energy is higher than the critical energy, where the slowing-down friction by the electrons becomes equal to that caused by the bulk ions. The NB heating power supplied to the electrons ${{P}_{\text{NB},\text{e}}}$ was greater than that supplied to the ions ${{P}_{\text{NB},i}}$. The electron temperature ${{T}_{\text{e}}}$ exceeded the ion temperature T_i. The NB driven current ${{I}_{\text{NB}}}$ accounted for 13% of the total plasma current ${{I}_{\text{p}}}$ at t   =   13.0 s. The bootstrap current sustained 8% of ${{I}_{\text{p}}}$. The ohmic current ${{I}_{\text{OH}}}$ decreased, thus, the total plasma current ${{I}_{\text{p}}}$ was maintained at the prescribed value.

Figure 4 shows the time evolution of the NB heating power ${{P}_{\text{NB}}}$, the NB driven current ${{I}_{\text{NB}}}$, and the loop voltage ${{V}_{\text{loop}}}$, in addition to the contour of the fast-ion distribution function $\mathcal{F}$ via a comparison between the cases with and without electric acceleration. The tangential N–NB-induced ions circulated around the toroidal loop at greater speeds than the bulk ions; hence, some level of electric acceleration was observed. In the fast-ion distribution function, considerably more co-circulating particles existed than counter-circulating particles because they came from the co-tangential N–NB. The electric acceleration effect was more noticeable on the current ${{I}_{\text{NB}}}$ than on the power ${{P}_{\text{NB}}}$. The parallel electric field shifted the distribution function towards the positive parallel direction in the phase space. The co-circulating particles were accelerated and the counter-circulating particles were decelerated. The larger amount of accelerated particles compared to the decelerated particles caused a slight increase in the heating power. In contrast, both the acceleration of the co-circulating particles and the deceleration of the counter-circulating particles contributed to the increase in driven current. The shift of the distribution function led to the difference in ${{P}_{\text{NB}}}$ and ${{I}_{\text{NB}}}$ as shown in figure 4. The reduction in ${{V}_{\text{loop}}}$ owing to the NB current drive was observed after t   =   5.0 s, and it was most apparent at $\rho =0.4$, where the NB deposition peak was located. The effect of the electric field on ${{I}_{\text{NB}}}$ remained at 10% or more after current ramp-up because the loop voltage ${{V}_{\text{loop}}}$ sustained by the ohmic current ${{I}_{\text{OH}}}$. The driven current ${{I}_{\text{NB}}}$ reached its peak at $t\approx 8$ s. Fast ions did not exist until the NB injection started at t   =   5.0 s, and then the amount of fast ions increased owing to the particle source S provided by the NB. The born fast ions with high energy slowed down and their energy approached the bulk thermal energy. Thermalisation and assimilation into the bulk ions took place in the form of the loss of fast ions. The thermalisation loss $L={{\nu}_{\text{loss}}}\mathcal{F}$ was proportional to the amount of the fast ions. The amount of the fast ions initially increased and then the loss L exceeded the source S. This led to a decrease in ${{I}_{\text{NB}}}$ after $t\approx 8$ s.

Zoom In Zoom Out Reset image size

Figure 5 shows the time evolution of the radial profiles for the prescribed electron density ${{n}_{\text{e}}}$, the safety factor q, the electron temperature ${{T}_{\text{e}}}$, the plasma current density ${{j}_{\text{p}}}$, the NB heating power density to the electrons ${{P}_{\text{NB},\text{e}}}$ and the NB-driven current density ${{j}_{\text{NB}}}$. The NB heating created a temperature gradient, as can be clearly observed at $\rho \approx 0.4$ and $t\approx 6$ s in figure 5(c). The NB current drive raised the plasma current at the same location $\rho \approx 0.4$ as shown in figure 5(d), cooperatively with the bootstrap current due to the temperature gradient. The local current ramp contributes to the formation of the reversed magnetic shear profile as shown in figure 5(b). The negative shear at $\rho <0.5$ lowered the turbulent transport level in CDBM, leading to a further temperature increase inside the transport barrier. At $\rho <0.4$, the neoclassical transport level increased owing to the low inverse aspect ratio and the relatively high q value, playing a key role in the strong flattening as seen in figure 5(c) instead of the lowered turbulent transport. The q profile was quite transient after the total current reached the flattop. It gradually approached the steady-state at the time scale of the current diffusion ${{\tau}_{\text{curr}}}$. The evolution of q was closely related to that of ${{j}_{\text{p}}}$ and, thus, the driven current ${{j}_{\text{NB}}}$. The driven current ${{j}_{\text{NB}}}$ increased ${{j}_{\text{p}}}$ at $\rho \approx 0.4$, and induced the counter electric field, which tended to cancel out the rapid variation in current. The counter field diffused and caused the counter current shown at $\rho \approx 0.2$ in figure 5(d). This evolution of ${{j}_{\text{p}}}$ led to variations in q, as shown in figure 5(b).

Zoom In Zoom Out Reset image size

3.2. ITER simulation

Next, the application of the integrated simulation to the ITER scenario is shown. Calculation conditions were set to be the same as in the previously reported ITER task [18] for the standard scenario, except for the NB injection conditions. Early NB injection before the current flattop has some merit in terms of current profile control and flux consumption reduction. Conventional integrated codes, however, have not had a transient fast-ion distribution function; and hence, we could not study such scenarios with sufficient applicability. Moreover, it is important to study the interaction between the fast ions and the Alfvén waves. The degradation of heating and the current drive efficiency by fast ions due to Alfvén instabilities are important subjects in ITER research. An introduction of the fast-ion redistribution model to the integrated code will be performed in future studies; nevertheless, obtaining the fast-ion distribution function within the integrated modelling framework is an indispensable step to that end.

In the present simulation, N–NB (1 MeV, 33 MW) was injected from t   =   60 s in the middle of the current ramp-up phase. The time evolution of the plasma parameters is shown in figure 6. After t   =   60 s, the plasma heating power ${{P}_{\text{tot}}}$ increasingly exceeded the external input power ${{P}_{\text{NB}}}$ because the fusion reaction was initiated. Additional resonant wave heating (7.0 MW of ion cyclotron range) was injected from t   =   90 s aiming at L–H transition. Further increase in ${{T}_{\text{i}}}$ led to an increase in alpha particle production. The production and subsequent collisional processes of the alpha particles were also calculated using BAFP in addition to the calculation of NB-induced fast ions. The fusion reaction between thermal deuteron and thermal triton was dominant at the flattop. At the beginning of the NB injection, the reaction between beam-induced deuteron and thermal triton occurred more often than the former. The beam-thermal reaction was calculated using the fast-ion distribution function. The resultant alpha heating power ${{P}_{\alpha}}$ is shown in figure 6. The unified handling of fast ions, that is, the NB-injected ions and the alpha particles, is adequate for integrated simulation in terms of unified implementation of transport processes and fusion reactions.

Zoom In Zoom Out Reset image size

Figure 7 shows the time evolution of the distribution function $\mathcal{F}$ of the NB-induced fast ions. The distribution function immediately after the start of the NB injection nearly represents the beam ionisation profile in the velocity space. The newborn fast ions were monoenergetic but had quite a broad pitch angle with respect to the magnetic field lines because of the large cross-section of the beam. In the high-energy region, the fast ions predominantly slowed down with slight pitch angle scattering. The fast ions reached the trapped/passing boundary jump to the counterpart of the boundary. After a thermalisation time of about 1 s, the fast ions reached the energy level of the bulk temperature. Energy and pitch angle diffusion processes in the high-energy region were slower than the thermalisation. The distribution function transformed to a quasi-steady state about 10 s after the start of the NB injection. The background plasma profiles and magnetic fields varied within the time required to complete those processes. The fast-ion distribution function is impacted by the history of the background variations. This feature is omitted in steady-state modelling and is thus an important advantage of the present integration of fast-ion modelling.

Zoom In Zoom Out Reset image size

Fast ions play a key role in terms of heating and current drive for the profile formation of advanced plasmas. A transient analysis via integrated modelling is required to explore a robust approach to such advanced plasmas. To simulate the self-consistent evolution of fast ions, bulk plasma and equilibrium magnetic field, we combined a BAFP code for the fast ions with the integrated code TOPICS. The BAFP modelling and our implementation were verified by comparison with the results of OFMC calculations.

The integrated simulation was applied to JT-60SA and ITER plasmas in the ramp-up phase to confirm the capability and effectiveness of the integrated code for transient analyses. It showed the coupled evolution of the profiles of NB heating, the current drive, the plasma pressure, and the safety factor. The alpha particles were handled by the BAFP simultaneously with the NB-injected ions, and the resultant alpha heating was considered in the ITER case. The reasonable results yielded by the integrated simulation gives us the confidence to proceed with predictive simulations, experiment analyses, and further model development.

We have considered the electric acceleration of fast ions in the integrated modelling framework. In recent-day highly non-inductive scenarios, the loop voltage at the flattop tends to be small; however, the electric field always appears in the transient phase, that is, not only in the ramp up/down phases but also in the transient events caused by instabilities and controls during the flattop phase. Moreover, as discussed in [5], a countering electric field would be imposed to suppress flux consumption during the ramp-up phase. The developed integrated modelling code can be used to study those advanced scenarios under a toroidal electric field.

In the present integration, the bounce averaged-version of the Fokker–Planck code is employed as a transient fast-ion solver. Its computational lightness, which allows it to solve a time-step evolution of the 2D distribution function within a second, facilitates consistent linkage with the bulk transport solver. Evolution of the fast-ion distribution function is shown in the integrated simulation. Its characteristic time is several seconds, the interval within which plasma profiles and magnetic fields vary. For further development towards the introduction of the fast-ion interaction with Alfvén waves, it is necessary to solve coupled evolution among the fast-ion distribution function, the plasma profiles, and the equilibrium magnetic field. Therefore, the present consistent linkage is an important improvement in integrated modelling for transient analyses.

In a comparison with the BAFP under identical simulation conditions, the conventional NB-driven current model in TOPICS, mentioned in the introduction, is found to be reasonable at least from a qualitative perspective. Conventional evaluation of the NB-driven current is in good agreement with the BAFP result at the flattop, but the conventional evaluation underestimates the NB-driven current soon after the NB injection, because a slowed down and pitch-angle-scattered steady-state solution is assumed, while the distribution function does not yet take this shape in the transient phase. In terms of computation time, in case of the use of ten processing cores, integrated calculation with the BAFP code requires almost the same time as that with the conventional model. The favourable computation speed encourages us to use the BAFP version in future integrated simulations.

The authors are grateful to Drs Y Kamada, K Shinohara and J Shiraishi for their valuable suggestions and discussions. They also thank Messrs. I Kamata and M Suzuki for their numerical assistance. This work was partially supported by JSPS KAKENHI Grant Number 15K18311.