Effect of toroidal rotation on impurity transport in tokamak improved confinement

The TASK code is one of the integrated codes developed in Japan, and has a modular structure combining many codes such as equilibrium, transport, and wave analyses. The multiple modules exchange data with each other through the data interface BPSD [30], which enables self-consistent simulations of time evolutions of plasma quantities.

In the transport calculations of plasma main component using the TASK/TR code (transport module for main plasma species), the following one-dimensional diffusive equations for particle and heat transport and magnetic diffusion are solved.

$$\begin{align} & \frac{{\text{1}}}{{V^{\prime}}}\frac{\partial }{{\partial t}}\left( {{n_s}V^{\prime}} \right)\nonumber\\ & \quad = - \frac{{\text{1}}}{{V^{\prime}}}\frac{\partial }{{\partial \rho }}\left( {V^{\prime}\left\langle {\left| {\nabla \rho } \right|} \right\rangle {V_s}{n_s} - V^{\prime}\left\langle {{{\left| {\nabla \rho } \right|}^2}} \right\rangle {D_s}\frac{{\partial {n_s}}}{{\partial \rho }}} \right) + {S_s},\end{align} \tag{ 1 }$$

$$\begin{align} & \frac{{\text{1}}}{{{{V^{\prime}}^{{\text{5/3}}}}}}\frac{\partial }{{\partial t}}\left( {\frac{{\text{3}}}{{\text{2}}}{n_s}{T_s}{{V^{\prime}}^{{\text{5/3}}}}} \right) \nonumber\\ & \quad = - \frac{{\text{1}}}{{V^{\prime}}}\frac{\partial }{{\partial {{\rho }}}}\left( V^{\prime}\left\langle {\left| {\nabla \rho } \right|} \right\rangle \frac{{\text{3}}}{{\text{2}}}{n_s}{T_s}{V_{{E_s}}} - V^{\prime}\left\langle {{{\left| {\nabla \rho } \right|}^{\text{2}}}} \right\rangle {\chi _s}\frac{{\partial {n_s}{T_s}}}{{\partial \rho }}\right.\nonumber\\ & \left.\qquad -\ \left\langle {{{\left| {\nabla \rho } \right|}^{\text{2}}}} \right\rangle \left( {\frac{{\text{3}}}{{\text{2}}}{D_S} - {\chi _s}} \right){T_s}\frac{{\partial {n_s}}}{{\partial \rho }} \right) + {P_s},\end{align} \tag{ 2 }$$

$$\begin{align} & \frac{\partial }{{\partial t}}\left( {\frac{{\partial \psi }}{{\partial \rho }}} \right)\nonumber\\ & \quad = - \frac{\partial }{{\partial \rho }}\left[ \frac{{{\eta _\parallel }}}{{{\mu _0}}}\frac{I}{{V^{\prime}\langle {R^{ - 2}}\rangle }}\frac{\partial }{{\partial \rho }}\left( {\frac{{V^{\prime}}}{I}\left\langle \frac{{{{\left| {\nabla \rho } \right|}^2}}}{{{R^2}}}\right\rangle \frac{{\partial \psi }}{{\partial \rho }}} \right) \right.\nonumber\\ & \left.\qquad -\ \frac{{{\eta _\parallel }}}{{I\langle {R^{ - 2}}\rangle }}\langle \left( {{J_{CD}} + {J_{BS}}} \right){B_t}\rangle \right].\end{align} \tag{ 3 }$$

Here, n_s, T_s, S_s, P_s, D_s, V_s, χ_s, V_Es, B_t, η_//, J_CD, J_BS, I, ρ = r/a, are the plasma density, the plasma temperature, the particle source, the heat source, the particle diffusion coefficient, the particle convective velocity, the heat diffusion coefficient, the heat convective velocity, the toroidal magnetic field, the plasma resistivity along the magnetic field lines, the induced current density, the bootstrap current, B_tR and the normalized minor radius, respectively, and V' = dV/dρ is the radial derivative of plasma volume and < > denotes a magnetic flux surface average. In the particle transport equation, the quasi-neutrality condition is assumed. Transport coefficients include neoclassical and turbulent contributions. The heating source of neutral beam injection (NBI) is given by

$$\begin{align}{P_{{\text{NB}}}} = {P_{{\text{NB}},0}}\exp \left[ { - {{\left( {\frac{{\rho - \rho _{\text{0}}^{{\text{NB}}}}}{{\rho _{\text{w}}^{{\text{NB}}}}}} \right)}^2}} \right].\end{align} \tag{ 4 }$$

Here, P_NB,0 is the intensity of NBI power, $\rho _{\text{0}}^{{\text{NB}}}$ is the radial position of NBI, $\rho _{\text{w}}^{{\text{NB}}}$ is the radial width of NBI. The lower hybrid (LH) heating profile is also given with the same form of the equation. A H-mode pedestal profile is given by a modified tanh function expressed as [31]

$$\begin{align} & X\left( \rho \right) = \frac{{{X_{{\text{ped}}}} - {X_{{\text{sol}}}}}}{2}\left\{ {{\text{mtanh}}\left( {\frac{{{\rho _{{\text{mid}}}} - \rho }}{{2{\rho _{\text{w}}}}},{\alpha _{{\text{ped}}}}} \right) + 1} \right\} + {X_{{\text{sol}}}},\end{align} \tag{ 5 }$$

$$\begin{align} & {\text{mtanh}}\left( {\rho ,{\alpha _{{\text{ped}}}}} \right) = \frac{{\left( {1 + {\alpha _{{\text{ped}}}}\rho } \right){e^\rho } - {e^{ - \rho }}}}{{{e^\rho } + {e^{ - \rho }}}}.\end{align} \tag{ 6 }$$

Here, X represents each plasma quantity, and X_ped, X_sol, ρ_mid, ρ_w and α_ped are the height of pedestal, the height of scrape-off layer, the middle position of pedestal, the radial width of pedestal, the radial profile fitting parameter, respectively.

The diffusion coefficients are given by linear summation of neoclassical and turbulent diffusion coefficients. For the neoclassical transport model, NCLASS [32] is applied. For the turbulent transport model of main plasma, CDBM (current diffusive ballooning mode) model is applied. The diffusion coefficients by the CDBM model are given as follows [33];

$$\begin{align}{\chi _{{\text{CDBM}}}} = {\chi _e} = {\chi _i} = CF\left( {s,\alpha } \right){F_\kappa }{F_{E \times B}}{\alpha ^{3/2}}\frac{{{c^2}}}{{\omega _{pe}^2}}\frac{{{v_A}}}{{qR}},\end{align} \tag{ 7 }$$

where, c is the speed of light, ω_pe is the electron plasma frequency, v_A is the Alfvén velocity, q is the safety factor, s ≡ (r/q)(dq/dr) is the magnetic shear, α ≡ −(2μ₀ q² R/B²)(dp/dr) is the normalized pressure gradient. The fitting parameter C = 12.0, which is correction between the CDBM confinement time and the $\tau _{{\text{89}}}^{{\text{ITER}}}$ confinement scaling [34]. The each fitting factor F(s, α), F_κ [35] and F_{E × B} [36] is defined as

$$\begin{align}F\left( {s,\alpha } \right) & = \left\{ \begin{array}{*{20}{c}} {\frac{1}{{\sqrt {2\left( {1 - s^{\prime}} \right)\left( {1 - 2s + 3{s^{^{\prime} 2}}} \right)} }}{\text{ }}\left( {s^{\prime} = s - \alpha < 0} \right)} \\ {\frac{{\left( {1 + 9\sqrt 2 {{s^{\prime}}^{3/2}}} \right)}}{{\sqrt 2 \left( {1 - 2s^{\prime} + 3{s^{^{\prime} 2}} + 2{{s^{\prime}}^3}} \right)}}{\text{ }}\left( {s^{\prime} = s - \alpha > 0} \right)} \end{array}\right.,\end{align} \tag{ 8 }$$

$$\begin{align}{F_\kappa } & = {\left( {\frac{{2\sqrt \kappa }}{{1 + {\kappa ^2}}}} \right)^{3/2}},\end{align} \tag{ 9 }$$

$$\begin{align}{F_{E \times B}} & = \frac{1}{{1 + \omega _{E \times B}^2/\gamma _{{\text{CDBM}}}^2}},\end{align} \tag{ 10 }$$

where ω_{E × B} = (RB_p /B_t )(dE_r /RB_p ) is the E × B flow shearing rate and γ_CDBM = |α|^0.5(v_A /qR)F(s, α) is the growth rate of CDBM [37]. The radial electric field from the radial force balance is expressed as

$$\begin{align}{E_r} = \frac{1}{{{Z_i}e{n_i}}}\frac{{{\text{d}}{p_i}}}{{{\text{d}}r}} + {V_{{\text{tor}}}}{B_p} - {V_{{\text{pol}}}}{B_t},\end{align} \tag{ 11 }$$

and the poloidal velocity V_pol is evaluated with NCLASS. An additional thermal diffusion coefficient χ_e = χ_e + 0.1exp[−(ρ/0.075a)²] is given for numerical stabilization, since the electron temperature can easily diverge due to the excessive radiation losses of impurities [38]. The turbulent diffusion coefficient is set as D_CDBM = χ_CDBM/2 and turbulent convective velocity is set as ${V_{{\text{turb}}}}\,{{ = }}\, - {\text{0}}{\text{.5}}{D_{{\text{turb}}}}r{\text{/}}a$ [39].

In the impurity transport calculations using the TASK/TI code (transport module for impurities), the following particle transport equations including atomic processes are solved for k-valent ionized impurity ions.

$$\begin{align}\frac{{\text{1}}}{{V^{\prime}}}\frac{\partial }{{\partial t}}\left( {V^\prime{n_k}} \right) & = - \frac{{\text{1}}}{{V^{\prime}}}\frac{\partial }{{\partial {{\rho }}}}\left( {V^{\prime}\langle \left| {\nabla \rho } \right|\rangle {V_k}{n_k} - V^{\prime}\langle {{\left| {\nabla \rho } \right|}^2}\rangle {D_k}\frac{{\partial {n_k}}}{{\partial \rho }}} \right)\nonumber\\ & \quad + {\gamma _{k - 1}}{n_e}{n_{k - 1}} - {\gamma _k}{n_e}{n_k} - {\alpha _k}{n_e}{n_k}\nonumber\\ & \quad + {\alpha _{k + 1}}{n_e}{n_{k + 1}}{\text{ }}\left( {k = 2 - 74} \right),\end{align} \tag{ 12 }$$

$$\begin{align}\frac{{\text{1}}}{{V{^{^{\prime}}}}}\frac{\partial }{{\partial t}}\left( {V{^{^{\prime}}}{n_{\text{1}}}} \right) & = - \frac{{\text{1}}}{{V{^{^{\prime}}}}}\frac{\partial }{{\partial {{\rho }}}}\left( {V{^{^{\prime}}}\langle \left| {\nabla \rho } \right|\rangle {V_1}{n_1} - V{^{^{\prime}}}\langle {{\left| {\nabla \rho } \right|}^2}\rangle {D_1}\frac{{\partial {n_1}}}{{\partial \rho }}} \right)\nonumber\\ & \quad - {\gamma _{\text{1}}}{n_e}{n_{\text{1}}} + {\alpha _{\text{2}}}{n_e}{n_{\text{2}}} - {\alpha _{\text{1}}}{n_e}{n_{\text{1}}} + {S_0}{ }\left( {k = 1} \right),\end{align} \tag{ 13 }$$

where, γ_k is the ionization rate, α_k is the recombination rate and S₀ is the impurity source. The impurity temperatures are assumed to be equivalent to the bulk ion temperature. The γ_k, α_k and radiation loss power of impurities are evaluated by using OPEN-ADAS database [40]. Impurity densities at the boundary are given to be fixed value with the charge distribution determined by the coronal equilibrium. The time evolution of bulk plasmas and impurities is solved by linking TASK/TR and TASK/TI modules. The data of plasma density and temperature profiles are transferred from TR to TI, while the data of radiation loss power and effective charge number are transferred from TI to TR.

Neoclassical transport coefficients for impurities are calculated by introducing geometrical factors (P_A and P_B ) into NCLASS in order to consider the poloidal asymmetry due to the toroidal rotation. The impurity flux is expressed in the following form [11, 41–45],

$$\begin{align}\langle {{\Gamma }}_Z^{{\text{Neo}}} \cdot \nabla r\rangle & = \langle {n_Z}\rangle {q^2}{D_c}Z\left[ {P_A}\left( {\frac{1}{Z}\frac{1}{{{L_{{n_Z}}}}} - \frac{1}{{{L_{{n_i}}}}} + \left( {C_0^Z - 1} \right)\frac{1}{{{L_{{T_i}}}}}} \right)\right.\nonumber\\ & \quad \left. -\ {P_B}\left( {C_0^Z + {k_i}} \right)\frac{1}{{{L_{{T_i}}}}} \right],\end{align} \tag{ 14 }$$

with

$$\begin{align}{P_A} & = \frac{1}{{2{ \epsilon ^2}}}\frac{{\left\langle {{B^2}} \right\rangle }}{{\left\langle {{n_Z}} \right\rangle }}\left[ {\left\langle {\frac{{{n_Z}}}{{{B^2}}}} \right\rangle - {{\left\langle {\frac{{{B^2}}}{{{n_Z}}}} \right\rangle }^{ - 1}}} \right],\end{align} \tag{ 15 }$$

$$\begin{align}{P_B} & = \frac{1}{{2{ \epsilon ^2}}}\frac{{\left\langle {{B^2}} \right\rangle }}{{\left\langle {{n_Z}} \right\rangle }}\left[ {\frac{{\left\langle {{n_Z}} \right\rangle }}{{\left\langle {{B^2}} \right\rangle }} - {{\left\langle {\frac{{{B^2}}}{{{n_Z}}}} \right\rangle }^{ - 1}}} \right],\end{align} \tag{ 16 }$$

where the gradient length of each plasma quantity 1/L_X = −∇X/X, and the aspect ratio = r/R. Absence of the poloidal asymmetry makes P_A = 1 and P_B = 0. The coefficients k_i and $C_{\text{0}}^Z$ are related to the collisionality, the neoclassical ion flow coefficient and the friction coefficient, respectively. For the case where the impurity is in the Pfirsch–Schlüter (PS) regime and the main ion is in the banana regime, $- \left( {C_{\text{0}}^Z - {\text{1}}} \right)$ becomes close to −0.5, and a numerical formula is often used [46, 47]. More recently, extended formulas have been proposed with ion-electron collisional heat exchange for impurities, and $C_{\text{0}}^Z$ is given as follows [48, 49]:

$$\begin{align}C_0^Z = \frac{{\frac{{1.5}}{{1 + {f_1}\left( {\frac{{{A_i}}}{{{A_A}}}} \right)}} - \frac{{0.29 + 0.68\alpha }}{{0.59 + \alpha + \frac{{1.34 + {f_2}}}{{{g^2}}}}}}}{{1 + {f_3}{\mu _{ie}}{g^2}}}.\end{align} \tag{ 17 }$$

On the other hand, $C_{\text{0}}^Z$ + k_i is 0.33 takes the limit of ≪ 1, and the expression for k_i is given by [50–52], which give the following expression including the effect of rotation [49],

$$\begin{align}{k_i} & = \left( {\frac{{{k_{i,0}} + {l_1}{{\left( {{f_t}\nu _i^*} \right)}^{0.5}} + {l_2}\nu _i^{*0.25}}}{{1 + {l_3}\nu _i^{*0.5}}} - {l_4}{l_5}\nu _i^{*2}f_t^{\,6} + {l_6}\nu _i^{*0.25}} \right)\nonumber\\ & \quad \times \left( {\frac{1}{{1 + {l_4}\nu _i^{*2}f_t^{\,6}}}} \right).\end{align} \tag{ 18 }$$

The variables f_i, l_i, g, α and μ_ie in equations (17) and (18) are detailed in [48] and [49]. Here, the trapped particle fraction f_t is evaluated with Lin–Liu and Miller model [53];

$$\begin{align}{f_t} = 0.75\left( {1.5\sqrt \epsilon } \right) + 0.25\left( {\frac{{3\sqrt 2 }}{\pi }\sqrt \epsilon } \right).\end{align} \tag{ 19 }$$

Equation (14) assumes a convection coefficient K = 1 multiplied to the impurity density gradient and main ion density gradient terms. This assumption is known to be valid for the case where the impurity is in the PS regime and the main ion is in the banana regime. Therefore, rewriting equation (14) to equation (7) in [43]. Using the temperature screening coefficient H, the neoclassical flux of the impurity is given as

$$\begin{align}\frac{{\langle {{\Gamma }}_Z^{{\text{Neo}}} \cdot \nabla r\rangle }}{{\langle {n_Z}\rangle }} & = {q^2}{D_c}Z\left[ {P_A}\left( {K\frac{1}{Z}\frac{1}{{{L_{{n_Z}}}}} - K\frac{1}{{{L_{{n_i}}}}} - H\frac{1}{{{L_{{T_i}}}}}} \right)\right.\nonumber\\ & \quad \left. -\ {P_B}\left( {C_0^Z + {k_i}} \right)\frac{1}{{{L_{{T_i}}}}} \right].\end{align} \tag{ 20 }$$

As shown in [43], the above analytical flux formula combined with the NCLASS code describes impurity transport including the poloidal asymmetry. Using the diffusion coefficient $D_{{\text{asym}}}^{{\text{ NCLASS}}}$ and convection velocity $V_{{\text{asym}}}^{{\text{ NCLASS}}}$ with the poloidal asymmetry (labeled with 'asym'), equation (20) can be simply expressed as

$$\begin{align}\frac{{\langle {{\Gamma }}_Z^{{\text{Neo}}} \cdot \nabla r\rangle }}{{\langle {n_Z}\rangle }} = D_{{\text{asym}}}^{\,{\text{NCLASS}}}\frac{1}{{{L_{{n_Z}}}}} + V_{{\text{asym}}}^{\,{\text{NCLASS}}},\end{align} \tag{ 21 }$$

where $V_{{\text{asym}}}^{{\text{ NCLASS}}}$ is defined as

$$\begin{align}V_{{\text{asym}}}^{\,{\text{NCLASS}}} = V_{N,{\text{asym}}}^{\,{\text{NCLASS}}} + V_{T,{\text{asym}}}^{\,{\text{NCLASS}}},\end{align} \tag{ 22 }$$

with the poloidally asymmetric convection velocity $V_{N,{\text{asym}}}^{{\text{ NCLASS}}}$ and $V_{T,{\text{asym}}}^{{\text{ NCLASS}}}$. Here the poloidally asymmetric transport coefficients can be written as follows using the symmetric (labeled with 'sym') transport coefficients $D_{{\text{sym}}}^{\,{\text{NCLASS}}}$, $V_{N,{\text{sym}}}^{{\text{ NCLASS}}}$ and $V_{T,{\text{sym}}}^{{\text{ NCLASS}}}.$

$$\begin{align}D_{{\text{asym}}}^{\,{\text{NCLASS}}} & = {q^2}{D_c}K{P_A} = D_{{\text{sym}}}^{\,{\text{NCLASS}}}{P_A},\end{align} \tag{ 23 }$$

$$\begin{align}V_{N,{\text{asym}}}^{\,{\text{NCLASS}}} & = ZD_{{\text{sym}}}^{\,{\text{NCLASS}}}\frac{1}{{{L_{{n_i}}}}}{P_A} = V_{N,{\text{sym}}}^{\,{\text{NCLASS}}}{P_A},\end{align} \tag{ 24 }$$

$$\begin{align}V_{T,{\text{asym}}}^{\,{\text{NCLASS}}} & = - ZD_{{\text{sym}}}^{\,{\text{NCLASS}}}\frac{1}{{{L_{{T_i}}}}}\left( {\frac{H}{K} + \frac{{C_0^Z + {k_i}}}{K}\frac{{{P_B}}}{{{P_A}}}} \right)\nonumber\\ & = V_{T,{\text{sym}}}^{\,{\text{NCLASS}}}\left( {{P_A} + \frac{{C_0^Z + {k_i}}}{H}{P_B}} \right).\end{align} \tag{ 25 }$$

Therefore, the geometric factor that quantifies the poloidal asymmetry is ${P_A}{\text{ = }}D_{{\text{asym}}}^{{\text{ NCLASS}}}/D_{{\text{sym}}}^{{\text{ NCLASS}}}{\text{ = }}V_{N,{\text{asym}}}^{{\text{ NCLASS}}}/V_{N,{\text{sym}}}^{{\text{ NCLASS}}}$, which increases the inward pinch driven by the main ion density gradient in the presence of strong asymmetry. Since H is $- \left( {C_{\text{0}}^Z - {\text{1}}} \right)$ from equation (14) and H is negative value except when the main ions are highly collisional regime, the screening effect driven by the main ion temperature gradient is reduced by the P_B term in equation (25) [42]. To evaluate the geometric factor P_A and P_B , a two-dimensional impurity distribution considering the poloidal asymmetry using flux surface averaged (FSA) density is given as follows [54]:

$$\begin{align}{n_Z}\left( {r,\theta } \right) = \langle {n_Z}\left( {r,\theta } \right)\rangle \frac{{\exp \left[ {M_Z^2\frac{{{R^2}\left( {r,\theta } \right)}}{{R_0^2}}} \right]}}{{\left\langle \exp \left[ {M_Z^2\frac{{{R^2}\left( {r,\theta } \right)}}{{R_0^2}}} \right]\right\rangle }}.\end{align} \tag{ 26 }$$

Here, only the centrifugal effects due to toroidal rotation with NBI is considered, and the poloidally varying electrostatic potential caused by toroidal rotation-induced electrostatic potential, ion cyclotron resonance heating (ICRH) and electron cyclotron resonance heating (ECH) is not considered. Assuming that the temperatures of the impurity species and bulk ions are comparable, the Mach number of impurities M_Z is determined by ${M_Z}{\text{ = }}\sqrt {{m_Z}/{m_i}} {M_i}$, and Mach number of the main ion M_i is determined by ${M_i}{\text{ = }}{V_{{\text{tor}}}}\sqrt {{m_i}/2{T_i}}$. The toroidal rotation velocity V_tor is a flux surface averaged quantity. The FSA value of the physical quantity X is calculated with

$$\begin{align}\langle X\rangle = \frac{{{{\mathop \int \nolimits}}\frac{X}{{{B_p}}}{\text{d}}{l_p}}}{{{{\mathop \int \nolimits}}\frac{1}{{{B_p}}}{\text{d}}{l_p}}},\end{align} \tag{ 27 }$$

where dl_p is the line element along the poloidal magnetic field.

In this simulation, the following turbulent diffusion coefficient, a combination of constant and parabolic ones, is used

$$\begin{align}{D_{{\text{W}},{\text{turb}}}} = \left\{ \begin{array}{*{20}{c}} {C_1^{{\text{W,turb}}}} \\ {D_1^{{\text{turb}}}{\rho ^2} + D_2^{{\text{turb}}}} \\ {C_3^{{\text{W,turb}}}} \end{array}{ }\begin{array}{*{20}{c}} {\left( {0.0 \unicode{x2A7D} { }\rho \unicode{x2A7D} \rho _1^{{\text{W,turb}}}} \right)} \\ {\left( {\rho _1^{{\text{W,turb}}} \unicode{x2A7D} { }\rho \unicode{x2A7D} \rho _2^{{\text{W,turb}}}} \right)} \\ {\left(\rho _2^{{\text{W,turb}}} < { }\rho \unicode{x2A7D} 1.0\right)} \end{array}\right.,\end{align} \tag{ 28 }$$

where,

$$\begin{align}D_1^{{\text{W,turb}}} & = \frac{{C_2^{{\text{W,turb}}} - C_1^{{\text{W,turb}}}}}{{{{\left( {\rho _2^{{\text{W,turb}}}} \right)}^2} - {{\left( {\rho _1^{{\text{W,turb}}}} \right)}^2}}}{ },\end{align} \tag{ 29 }$$

$$\begin{align}D_2^{{\text{W,turb}}} & = \frac{{C_1^{{\text{W,turb}}}{{\left( {\rho _2^{{\text{W,turb}}}} \right)}^2} - C_2^{{\text{W,turb}}}{{\left( {\rho _1^{{\text{W,turb}}}} \right)}^2}}}{{{{\left( {\rho _2^{{\text{W,turb}}}} \right)}^2} - {{\left( {\rho _1^{{\text{W,turb}}}} \right)}^2}}}{ },\end{align} \tag{ 30 }$$

and the turbulent convection component of impurity is ignored.

We have introduced the impurity transport model described in section 2, and for the first step, the validity of the impurity transport simulation model is checked for the H-mode case in [42]. For that purpose, calculations with fixed background plasma profiles are carried out. As shown in figure 3, the plasma density, temperature and toroidal rotation velocity profiles are used for two time slices at t= 9.32 s, when the W is not peaking in the center of the plasma (case A), and at t= 12.37 s, when the W is peaking in the plasma center region (case B). Here, it is assumed that the electron and ion densities are the same, and the electron and ion temperatures are the same. In case A, the plasma density is hollow in the plasma center, and has positive gradient. In case B, the plasma density is peaked in the plasma center, and has negative gradient. The toroidal rotation velocity profiles are flat in the plasma center with the Mach numbers of the main ion at the plasma center M_i = 0.19 in case A and M_i = 0.30 in case B. The W transport analyses using NEO and GKW show that the neoclassical transport is dominant in the plasma center, while the turbulent contribution is also significant in the plasma periphery [42]. The turbulent diffusion component is nearly equal or larger than the neoclassical component for ρ> 0.45, and the turbulent convection component is larger and directing inward for ρ> 0.8. The turbulent diffusion coefficient of W is estimated from D_W,GKW/ χ_{i, NEO} in [42] to be set as in figures 4(a) and (b) by using equation (28), neglecting convective velocity. The neoclassical thermal diffusivity profile of the main ion obtained by NCLASS is used for the estimation assuming χ_{i, NEO}= χ_{i, NCLASS}. The density of W at the boundary is set as n_w,edge = 2.5 × 10¹⁴ m⁻³ in case A and n_w,edge = 1.0 × 10¹⁴ m⁻³ in case B.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Impurity profiles in the steady states in cases A and B are compared by using TASK/TI. The diffusion coefficient profiles are shown in figures 4(a) and (b). The neoclassical diffusion coefficient D^NC is larger for ρ < 0.6 in case B because the higher density of main ion and W results in higher collisionality. In particular, the increase of D^NC is significant for ρ < 0.15, where the W density profile is strongly peaking. As shown in figure 5, W does not accumulate in case A at the central region, whereas a strong accumulation occurs in case B. In case B, n_WZ at the plasma center is about 1.5% of the electron density, so the quasi-neutrality condition is not significantly affected in fixed background plasma density. The densities at the low field side (LFS) in figures 6(c) and (d) are larger than the magnetic flux surface-averaged (FSA) W density, indicating the W density distribution shifts towards the LFS side due to centrifugal effects as in figures 6(a) and (b). Figure 7 shows that the normalized density gradient profiles of W are strongly depending on the magnitude and direction of the neoclassical convection velocities shown in figures 4(c) and (d). In both A and B cases, the neoclassical convective components $V_{T,{\text{asym}}}^{{\text{ NCLASS}}}$ driven by the main ion temperature gradient are positive in all the radial region. Furthermore, in case A, $V_{N,{\text{asym}}}^{{\text{ NCLASS}}}$ is positive for ρ < 0.6 due to the positive bulk density gradient, resulting in the absence of W density peaking at the plasma center. On the other hand, $V_{T,{\text{asym}}}^{{\text{ NCLASS}}}$ remains positive for ρ > 0.6, whereas $V_{N,{\text{asym}}}^{{\text{ NCLASS}}}$ becomes negative, and the difference of these terms gives negative $V_{{\text{asym}}}^{{\text{ NCLASS}}}$ for ρ > 0.7. As a result, inward convection is driven only in the outer region of the plasma, where a relatively small peak of the W density is formed. In case B, $V_{N,{\text{asym}}}^{{\text{ NCLASS}}}$ is negative in all the radial region, and $V_{{\text{asym}}}^{{\text{ NCLASS}}}$ is negative. In the outer region of the plasma, since the Mach numbers have similar value in both cases, the geometrical coefficients P_A are similar, but at the plasma center, P_A is increased to 40 in case B, while 10 in case A. Therefore, in case B, a sharp peak of the W density is formed at the plasma center. As shown in figures 6(a) and (b), the W density peak is shifted to the LFS. The poloidal asymmetry appears clearer in case A, when the central W density is lower. Note that in case B, the structure is also asymmetric. From figure 12 in [42], the reconstructed experimental values from SXR in case A are n_W,LFS ∼ 5.6–9.2 × 10¹⁵ m⁻³ at ρ ∼ 0.7 and n_W,LFS ∼ 0.5–1.4 × 10¹⁵ m⁻³ at ρ ∼ 0.5, and are in good agreement with the results in figures 6(c) and (d). Similarly, the experimental values in case B are n_W,LFS ∼ 0.3–1.1 × 10¹⁶ m⁻³ at ρ ∼ 0.7 and n_W,LFS ∼ 0.6–1.5 × 10¹⁶ m⁻³ at ρ ∼ 0.15, and are also in good agreement. As shown in figure 7, the normalized W density gradients are similar to the results obtained with NEO and GKW in [42]; the difference around ρ = 0.3 in case A is due to the larger turbulent diffusion coefficient evaluated by GKW.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Next, the impurity transport calculation is coupled with that of main plasma species. The cases with improve confinement plasmas in section 3 are analyzed to include the effects of impurity. In the H-mode case, the same W turbulent diffusion coefficient as in case B in section 4.1 are used. In the hybrid scenario case, the turbulent diffusion coefficient is estimated from [55] to have $\rho _{\text{1}}^{{\text{W,turb}}}$ = 0.4, $\rho _{\text{2}}^{{\text{W,turb}}}$ = 0.8, $C_{\text{1}}^{{\text{W,turb}}}$ = 0.15, $C_{\text{2}}^{{\text{W,turb}}}$ = 1.9, and $C_{\text{3}}^{{\text{W,turb}}}$ = 1.5. In the advanced scenario plasma, $\rho _{\text{1}}^{{\text{W,turb}}}$ = 0.6 is set from the previous conditions. The W density at the boundary n_w,edge = 1.0 × 10¹⁴ m⁻³, which is same with the case in section 4.1, is used for all the cases. The toroidal rotation velocity profiles shown in figure 2(d) give central Mach numbers M_Z ∼ 2.5 (H-mode), 3.5 (Hybrid), and 2.9 (Advanced).

In the H-mode case, peaked W density is obtained as show in figure 8. The normalized W density gradient in the center region of the H-mode case is similar to R/L_nW ∼ 50 comparable to the case with the fixed background plasma profile, which reproduces the result of the central accumulation. With the inclusion of W, the temperature gradient of the plasma is significantly reduced at ρ ∼ 0.1, due to stronger radiation losses with the impurity accumulation. Compared to figure 2(b), the electron temperature at ρ = 0.1 is reduced by about 4% and the normalized electron temperature gradient is reduced by about 40% due to the W accumulation. In this state, the dominant convection component is inward one driven by the bulk density gradient, as shown in figure 9(b). Thus, in this case of the W accumulation is large enough to cause a strong change in the electron temperature profile, and interactions occur between the main plasma and the impurity as follows: (1) increase radiation losses, (2) reduce main plasma temperature gradient and increase main plasma density gradient, (3) increase neoclassical inward flow, and (4) further W accumulation. A loop from (1) to (4) is followed until a steady state, resulting in enhanced and sustained W density peaking. Furthermore, the temperature screening effect is reduced by the term proportional to P_B in equation (25) in the presence of a strong asymmetry in the W density.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the hybrid scenario case, the W density also has a peak at the center of plasma. For ρ < 0.4, as shown in figure 10, the value of both $V_{T,{\text{asym}}}^{{\text{ NCLASS}}}$ and $V_{N,{\text{asym}}}^{{\text{ NCLASS}}}$ increase toward the center of plasma due to an increase in the main ion temperature gradient by the reduced turbulent diffusion and an increase in the main ion density gradient by the NBI particle source. The convective component driven by the main ion density gradient is stronger, resulting in an inward particle pinch for ρ < 0.35. The temperature screening effect is reduced in the center of the plasma, where the W density peaking is large, as in the H-mode case. Compared to the H-mode case, the W density is lower, whereas the W normalized density gradient is of similar value in the plasma center region. This difference is due to the strong inward neoclassical particle pinch driven in the H-mode case in the plasma periphery 0.85 < ρ < 1.0, which causes the inflow into the plasma confinement region with a significant amount of tungsten.

Zoom In Zoom Out Reset image size

In the advanced scenario case, as shown in figure 11, the particle pinch driven by the main ion density gradient is stronger due to the effects of reduced turbulent diffusion and the NBI particle source, as in the hybrid scenario case. Note that the neoclassical convection components are larger than in the hybrid scenario case, even though the Mach number is smaller, because the plasma transport is improved, and the plasma gradient is stronger than in the hybrid scenario case. The normalized W density gradient changes sharply at ρ ∼ 0.4, where the ITB is formed, and the W density peaks becomes wider in the radial direction than the hybrid scenario case. This result indicates the acceleration of W accumulation in the ITB region, as confirmed by previous experiments [14–19].

Zoom In Zoom Out Reset image size

The dependence on Mach number of the main ion is investigated in the advanced scenario case to clarify the effect of toroidal rotation on the impurity profile formation. To avoid violation of the impurity trace limit condition α_trace = n_ZZ²/n_i ≪ 1 and breakdown of calculation due to excessive radiation losses, the boundary condition of the impurity density is set to be smaller value n_w,edge = 1.0 × 10¹³ m⁻³. Mach number M_i0 at the plasma center is used here as a parameter, and the central toroidal velocity is set by ${V_{{\text{tor}},{\text{0}}}}{{ = }}\sqrt {{\text{2}}{T_i}/{m_i}} {M_{i{\text{0}}}}$. The profile of the toroidal velocity is given by using the function ${V_{{\text{tor}}}}{{ = }}{V_{{\text{tor}},{\text{0}}}}{\left( {{\text{1}} - {\rho ^{\text{3}}}} \right)^{\text{2}}}$ to be similar to that of the hybrid case in figure 2(d).

In the small M_i0 case as M_i0 = 0.1, the W density profile flattens due to the turbulent diffusion without the enhacement by the toroidal rotation. The case with M_i0 = 0.1 gives P_A = 3 at the plasma center, while the case with M_i0 = 0.5 gives P_A = 150, indicating formation of a strong asymmetric structure at the LFS. Figure 12 shows the dependence of the W density profile. The W density profiles in the high Mach number cases have the same kind of the structure as in the previous section. The W density profile flattens due to the turbulent diffusion, decreasing the neoclassical transport in the M_i0 = 0.1 case. In addition, in the region near ρ ∼ 0.6, where outward W convection is driven, the enhancement by the toroidal rotation gives impurity ejection. The poloidal asymmetry of the W density distribution as in equation (26) scales to $M_Z^{{\text{ 2}}}$, and the W density is localized in the central region, so the value of P_A is considerably larger at high Mach numbers. Figure 13 shows the dependence of the W normalized density gradient on the Mach number of W, plotting at ρ= 0.2, where the W density peaking is strong, ρ= 0.4, where it is close to the ITB-foot, and ρ= 0.6, where the outward particle pinch is driven. The value of the normalized W density gradient tends to decrease in accordance with decreased toroidal rotation velocity, especially at ρ = 0.2. Except for the higher Mach numbers, the values of $V_{N,{\text{asym}}}^{{\text{ NCLASS}}}$ and $V_{T,{\text{asym}}}^{{\text{ NCLASS}}}$ decrease with decreasing Mach number. At ρ = 0.6, the value of the W normalized gradient decreases from M_Z ∼ 3–1 due to the reduction of enhanced temperature screening effect, but at ρ = 0.4, the value of the W normalized density gradient does not decrease due to a delicate balance of $V_{N,{\text{asym}}}^{{\text{ NCLASS}}}$ and $V_{T,{\text{asym}}}^{{\text{ NCLASS}}}$. On the other hand, the normalized W density gradient is unchanged or even decreased at higher Mach numbers. This is due to the suppression of the temperature screening coefficient H/K via the M_i dependence of k_i as shown in [49]. In [58], it is also theoretically shown that the W transport is suppressed in the high Mach number regime.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The dependence on ITB position is investigated in the advanced scenario case. The plasma parameters in section 4.2 are used, and the LH heating profile is changed with heating position $\rho _{\text{0}}^{{\text{LH}}}$ = 0.4–0.7. The obtained plasma profiles show formation of the ITB, and the radial position of the ITB-foot is consistent with the LH heating position $\rho _{\text{0}}^{{\text{LH}}}$ (figure 14). In the central region of the plasma, the ion temperature profile is peaking, but the electron temperature is flat due to the effect of radiation losses caused by W accumulation. The global plasma performance is improved to have H-factor H_98y2= 1.39 and normalized-β value β_N = 2.59 for $\rho _{\text{0}}^{{\text{LH}}}$ = 0.6, compared to H_98y2= 1.25 and β_N = 2.39 for $\rho _{\text{0}}^{{\text{LH}}}$ = 0.4. As shown in figure 15, the electron density at the plasma center is changed little with that at the ITB position, while the volume-averaged electron density increases from <n_e > = 3.79 × 10¹⁹ m⁻³ to <n_e > = 3.85 × 10¹⁹ m⁻³ for $\rho _{\text{0}}^{{\text{LH}}}$ = 0.4–0.7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The W accumulation inside the ITB is suppressed in the cases of the ITB formed in the outer region. In the cases where the LH heating position $\rho _{\text{0}}^{{\text{LH}}}$ are at the outer region, the particle source due to NBI is relatively small, so the outward $V_{T,{\text{asym}}}^{{\text{ NCLASS}}}$ driven by the main ion temperature gradient suppresses the inward particle pinch. On the other hand, in the cases where the $\rho _{\text{0}}^{{\text{LH}}}$ is close to the center of plasma, the inward $V_{N,{\text{asym}}}^{{\text{ NCLASS}}}$ causes a sharp increase in the W density gradient at the ITB region (figure 8). The concentration of W, c_W= <n_W>/<n_e >, also decreases from 1.37 × 10⁻⁵ (for $\rho _{\text{0}}^{{\text{LH}}}$ = 0.4) to 0.89 × 10⁻⁵ (for $\rho _{\text{0}}^{{\text{LH}}}$ = 0.7). The W accumulation over the entire plasma can be suppressed by increasing the temperature screening effects.