Anomalous transport driven by ion temperature gradient instability in an anisotropic deuterium-tritium plasma

1.   Introduction

The confinement of particle and energy is a key issue in tokamaks. Since the experimental observation of the anomalous ion thermal transport in DIII-D tokamak [1], the research for the heat transport has been developed prosperously. The ion temperature gradient (ITG) instability is considered to play an important role in the anomalies in ion thermal transport [2, 3].

The fusion temperature can be reached through the intense auxiliary heating, such as the neutral beam injection (NBI) heating, ion cyclotron resonance heating (ICRH) as well as electron cyclotron resonance heating. These heating regimes have been widely used in the existent tokamaks such as EAST, JET, DIII-D, and will be adopted in the fusion tokamaks such as ITER and CFETR for the realization of deuterium-tritium (D-T) discharge. Nevertheless, the auxiliary heating regimes could produce strong plasma anisotropy. Generally, tangentially injected NBI drives the anisotropy with T_‖ > T_⊥,and the ICRH typically generates perpendicularly dominant anisotropy with T_⊥ > T_‖,where T_‖ and T_⊥ are the parallel and perpendicular temperatures with respect to the magnetic field. Anisotropic plasma has been experimentally observed. For instance, studies in MAST suggest that the ratio T_⊥ /T_‖ = 1.7 is achieved [4]. In EAST, the degree of the parallel-dominated anisotropy T_‖ /T_⊥ = 1.8 is reached at low electron density with low hybrid wave heating [5].

Influences of the temperature anisotropy on plasma behaviours have been widely researched. The temperature anisotropy plays important roles in the ITG instability [6], the resistive wall mode [7], the low-frequency electromagnetic wave [8], the geodesic acoustic mode [9–11], the zonal flows [12, 13] and so on. In detail, in [6], it is demonstrated that instability of ITG mode is reduced by an ion temperature anisotropy of higher perpendicular temperature or high enough parallel temperature. In [13], it is found that the residual zonal flow driven by the ITG instability shows a strong dependence on the temperature ratio ${\left(T_{\perp }/T_{\Vert }\right)}^{3/2}$ . In [14], it is reported that in the ITG driven turbulence in tokamak plasmas, the different parallel and perpendicular temperatures T_‖ ≠ T_⊥ induced by the ion Bernstein wave heating plays an important role in the generation and dynamics of transport barriers. In [15], the theoretical analysis and numerical evaluation of the impurity flux for the stellarator plasmas suggest that the main ion pressure anisotropy must be taken into account for a correction description of impurity transport in most of the collisionality range.

In this work, we focus on the ITG instability in a temperature anisotropic D-T plasma. The linear dispersion relation is theoretically derived and numerically solved to analyze the effects of temperature anisotropy on the linear behaviour of ITG instability; based on the linear results, the mixing length model approximation is adopted to analyze the anomalous particle and energy transport. The remaining part of this paper is organized as follows. In section 2, dispersion relation of the ITG instability in an anisotropic D-T plasma is derived; in section 3, the linear frequency of ITG instability is numerically evaluated; in section 4, by adopting the mixing length model approximation, the particle flux and energy flux driven by the ITG instability are analyzed; section 5 is the conclusion and discussion.

2.   Dispersion relation

In this section, we derive the dispersion relation in an anisotropic D-T plasma. We consider the case with the electrostatic fluctuation δϕ(r, t), where r denotes the particle position. We introduce the Fourier transformation for the perturbed quantities

where k is the wave vector and ω = ω_r + iγ is the frequency with ω_r and γ the real and imaginary part, respectively. The dispersion relation for ITG instability expressed in the Fourier representation is given by the quasi-neutrality condition

Here, j denotes the ion species, which is adopted as D and T in this work. e is the elementary charge, Z_j is the charge number of j-species. $n_{\mathit{k}}^{\mathrm{e}}$ and $n_{\mathit{k}}^j$ are the perturbed density of electron and j-species, respectively.

Generally, for a ITG-dominated regime, the electron perturbed density is expressed with the adiabatic approximation, namely, with n_e0 and T_e the electron equilibrium density and temperature, respectively. The ion perturbed density is calculated as with $\int \mathrm{d}^3 v$ the velocity space integration and the perturbed distribution function. In the following part, the calculations for the perturbed density of D and T are treated equally, and the label 'j' is omitted for simplicity when no ambiguity is introduced.

By neglecting the parallel effects, the perturbed distribution function is yielded from the gyrokinetic (GK) equations [16, 17]

Here, $\langle A\rangle \equiv \frac{1}{2 \pi} \int_0^{2 \pi} \mathrm{d} \xi A$ denotes the gyro-average. F₀ = F₀(r, μ, E) is the equilibrium distribution function with r the radial coordinate, μ the magnetic moment and E the energy variable. B and b denote the amplitude and the direction of the background magnetic field. denotes the magnetic drift frequency, where k_θ is the poloidal wavenumber, R is the major radius, θ is the poloidal angle, is the magnetic shear, T_⊥ is the perpendicular temperature, and are the parallel and perpendicular velocities normalized by v_{t, ⊥},respectively, where with m the ion mass. is the 0th-order Bessel function, λ_s = k_θρ_s, ρ_s = c_s/Ω with the acoustic speed and Ω = ZeB/m the gyro-frequency. Note that in equation (3), the first term denotes the polarization part, the second term denotes the gyrocenter part. In appendix, we give a brief derivation of equation (3) based on the traditional GK theory [18] for the readers who are interested in.

Here, we make some remarks on equation (3). Generally, the trapped ion contribution is important when the frequency of perturbation is close to the ion bounce frequency. Note that the driving force for local ITG instability is due to the diamagnetic drift combined with temperature gradient, the trapped ion contribution can be neglected for simplicity. Besides, only the ion resonance due to the diamagnetic drift is kept and the parallel effect is neglected; the approximation about neglecting the parallel effect is justified when the mode is localized in the outer midplane [19]; correspondingly, we replace the local value of ω_d with the value of θ = 0 in the following part.

To proceed, following the general introduction of anisotropy, we set the equilibrium distribution for each ion species as a bi-Maxwellian form [6, 8, 9, 13]

where n₀ is the equilibrium density. The analyses of the magnetic equilibrium in [9, 20] show that T_‖ = T_‖(ψ ), T_⊥ = BT_‖(ψ )/(B - C₀(ψ )) with ψ(r) the magnetic surface variable, C₀(ψ ) is any function with respect to ψ ; the forms of T_‖ and T_⊥ lead to ∂_θF₀ = 0. Besides, the total plasma current is found to be $I=I_0(\psi)+\mathcal{O}[(\alpha-1) \beta]$, where α = T_⊥ /T_‖ is introduced to represent the temperature anisotropy, β is the ratio of plasma pressure to the magnetic pressure, I₀(ψ ) is the total plasma current in an isotropic plasma. In a low-β plasma, the corrections to the magnetic field are very small and invariably ignored [9, 20]. So the magnetic equilibrium is considered as unaffected by the anisotropy in this work.

After calculating the necessary deviations of F₀ and substituting them into equation (3), we obtain

Here, $\widehat{\omega }=\omega /{\omega }_0$ with ${\omega }_0=-\frac{k_{\theta }{T}_{\perp }}{\mathit{ZeBR}}$ . ν = 1, 2, 3, 4, 5. In order to conveniently analyze the transport coefficients and fluxes in section 4, analogous to the treatment in [21], we have defined a set of driving forces

and a set of weighting functions

Here, κ is the magnetic curvature. In the large-aspect-ratio approximation, $\mathcal{X}^4=\mathcal{X}^5=1 / R$. Note that $\mathcal{X}^4 h_4$ originates from the temperature anisotropy, that is, when α = 1, this term vanishes; $\mathcal{X}^5 h_5$ comes from the diamagnetic drift term.

Next, we calculate the velocity integral of $\left\langle f_k\right\rangle$ to obtain the perturbed density. Before the detailed calculation, we introduce the velocity transformation $\left(\hat{v}_{\|}, \hat{v}_{\perp}\right) \rightarrow\left(\hat{v}_{\|}, \mathcal{E}\right)$ with $\mathcal{E}=\hat{v}_{\|}^2+\hat{v}_{\perp}^2 / 2$ to rigorously treat the magnetic drift [22], then $\int \mathrm{d}^3 \boldsymbol{v}=2 \pi \int_0^{\infty} \mathrm{d} \mathcal{E} \int_{-\sqrt{\mathcal{E}}}^{\sqrt{\mathcal{E}}} \mathrm{d} \hat{v}_{\|}$. After some straightforward calculation, the perturbed density can be formally written as

Here, Γ₀ = I₀(b)e^-b with I₀ the modified 0th-order Bessel function of the first kind, . , and , where κ_n = R/L_n with , with and with . $K_0=\int_0^{\infty} \mathrm{d} \mathcal{E} \frac{V_0}{2 \mathcal{E}-\hat{\omega}}$, $K_1=\int_0^{\infty} \mathrm{d} \mathcal{E} \frac{\mathcal{E} v_0}{2 \mathcal{E}-\hat{\omega}}$ and $K_2=\int_0^{\infty} \mathrm{d} \mathcal{E} \frac{V_2}{2 \mathcal{E}-\hat{\omega}}$, which represent the $\mathcal{E}$-integral in the velocity integral;

which represents the ${\widehat{v}}_{\Vert }$ -integral in the velocity integral.

Note that in equation (8), the factor $\widehat{\omega }$ in $n_{\mathit{k}}^{\mathrm{D}}$ and $n_{\mathit{k}}^{\mathrm{T}}$ is normalized by ω_0D and ω_0T, respectively. In deriving the dispersion relation, we unify the normalization by ω_0D, then the factor $\widehat{\omega }$ in $n_{\mathit{k}}^{\mathrm{T}}$ should be replaced by $C\widehat{\omega }$ with $C=\frac{Z_{\mathrm{T}}{T}_{\perp ,\mathrm{D}}}{Z_{\mathrm{D}}{T}_{\perp ,\mathrm{T}}}$ . By substituting equation (8) into (2), we obtain the dispersion relation for the ITG instability in an anisotropic D-T plasma

Here, we recover the label 'j' to emphasize the role of D and T in the dispersion relation. τ_j = T_⊥,j/T_e. ε_T = n_T0/n_e0 denotes the fraction of T, from which the fraction of D is calculated as ε_D = (1 - Z_Tε_T)/Z_D. By setting ε_T = 0 and α_D = 1, equation (9) reduces to the well-known form with local kinetic limit [19, 23].

3.   Linear version

In this section, equation (9) is numerically evaluated to study the effects of the temperature anisotropy and the fraction of T on the linear frequency. For simplicity, only the temperature anisotropy for D is considered, that is, we set α_T = 1; correspondingly, we denote α_D as α for short in the following. Since the temperature anisotropy observed in the experiments can reach T_⊥ /T_‖ = 1.7 [4] and T_‖ /T_⊥ = 1.8 [5], we vary α from 1/2 to 2 in the following calculation; we mention that the perpendicular temperature is fixed, variation of α means variation of the parallel temperature T_‖ . Besides, we set the perpendicular temperature of D to be the same as that of T, namely, T_⊥,D = T_⊥,T; the ITG anisotropy is neglected for simplicity. κ_{T, D} = κ_{T, T} and κ_{n, D} = κ_{n, T} are assumed. τ_D = τ_T = 1. We consider three typical cases: one is the cyclone base case (CBC) [24], which is widely-used in the GK simulation; the other two are the ITER standard scenario [25] and the CFETR steady-state scenario [26], for the purpose to investigate the ITG instability in the future fusion devices. The key parameters used in the three cases are shown in table 1; in each case, the ITG instability dominates [24–26].

Table  1.  Key parameters of three cases. a is the minor radius.

Case	R0 (m)	a (m)	B0 (T)	T⊥ (keV)	κn	κT
CBC	1.67	0.6	1.9	1.97	2.23	6.96
ITER	6.2	2.0	5.3	10.0	3.0	6.0
CFETR	5.7	1.6	5.0	12.5	3.93	8.91

 | Show Table

DownLoad: CSV

Equation (9) is solved as follows: for given (α, ε_T, k_θρ_{s, D}), we scan $({\widehat{\omega }}_{\mathrm{r}},\widehat{\gamma })$ to numerically calculate the left-hand-side expression of equation (9); the increment for ${\widehat{\omega }}_{\mathrm{r}}$ and $\widehat{\gamma }$ are both 5 × 10^-3. When the absolute value of the left-hand-side expression is the smallest, the corresponding $({\widehat{\omega }}_{\mathrm{r}},\widehat{\gamma })$ is considered as the solution of equation (9). For each given (α, ε_T, k_θρ_{s, D}), the absolute values of both the real and imaginary parts of the left-hand-side expression should be smaller than 5 × 10^-3. Figure 1 shows the real frequency (a)–(c) and growth rate (d)–(f) versus k_θρ_{s, D} in the CBC case with three different temperature anisotropic factor α = 1/2, 1 and 2; the T fraction ε_T varied from 0.0 to 0.5. Totally speaking, when α = 1, for a fixed ε_T, tendencies of the dependence of real frequency and growth rate on k_θρ_{s, D} are similar as the simulation result from the gyrokinetic simulation [27], except that the growth rate is overestimated. In [28], it is pointed out that the overestimated growth rate may be induced by the assumption that the ITG instability is localized at the midplane. It is shown that, when α = 1, for a fixed k_θρ_{s, D}, dependence of the real frequency on ε_T is not obvious; dependence of the growth rate on ε_T is not obvious in small k_θρ_{s, D} (0.1–0.2) region, the growth rate decreases as ε_T increases in large k_θρ_{s, D} (0.3–0.5) region. Note that in our modelling, when α = 1, the equilibrium distribution functions of D and T are the same, then the difference between the role of D and T in the dispersion relation shown in equation (9) is reflected by the finite Larmor radius k_θρ_{s, j}. In large k_θρ_{s, D} region, this difference is obvious. Compared to the real frequency, effects of the finite Larmor radius on the growth rate are more obvious.

Figure  1.  Real frequency (a)–(c) and growth rate (d)–(f) versus k_θρ_{s, D} with ε_T varied in the CBC case. α = 1/2, 1 and 2. The normalized frequency c_s/R = 156 kHz; R = R₀ + 0.5a is adopted, where the biggest temperature gradient occurs.

DownLoad: Full-Size Img PowerPoint

When α = 1/2, for a fixed k_θρ_{s, D}, the real frequency decreases while the growth rate increases as ε_T increases; on the contrary, the real frequency increases while the growth rate decreases as ε_T increases when α = 2. In small k_θρ_{s, D} region, dependence of the growth rate on ε_T is more obvious when α = 1/2 than those when α = 1 and 2. It is also indicated that the ITG instability is destabilized (stabilized) when α = 1/2 (α = 2) by increasing ε_T.

Figure 2 shows the real frequency and growth rate versus k_θρ_{s, D} in the CBC case with three different T fractions ε_T = 0.0, 0.1 and 0.5, respectively; the temperature anisotropic factor α is varied from 1/2 to 2. It is shown that, tendencies of the dependence of real frequency on α at different ε_T are similar, that is, the real frequency decreases as α increases; the decrement is obvious in large k_θρ_{s, D} and/or small ε_T region. Tendency of the dependence of growth rate on α is opposite to that of the real frequency. In detail, for a fixed ε_T, the growth rate increases as α increases, especially in large k_θρ_{s, D} and/or small ε_T region. It is interesting to find that, in the region of α < 1, dependence of the growth rate on α is more significant compared to that in the region of α > 1. We take ε_T = 0.5 and k_θρ_{s, D} = 0.4 as an example. When α increases from 1/2 to 1, the growth rate increases from about 0.38 to about 0.48, the increment is about 30%; when α increases from 1 to 2, the growth rate increases from about 0.48 to about 0.52, the increment is about 8%. It is indicated that the ITG instability may be significantly stabilized by controlling α < 1.

Figure  2.  Real frequency (a)–(c) and growth rate (d)–(f) versus k_θρ_{s, D} with α varied in the CBC case. ε_T = 0.0, 0.1 and 0.5.

DownLoad: Full-Size Img PowerPoint

Besides, it is seen that as ε_T increases, the real frequency and growth rate with α ≠ 1 are close to those with α = 1. In another word, effect of the temperature anisotropy of D on the frequency becomes weaker with higher ε_T. Note that this result relates to the assumption that the temperature anisotropy is only from D; as ε_T increases, the fraction of D decreases, as a result, effect of the temperature anisotropy becomes weak.

Here, we give a phenomenological explanation about the dependence of linear frequency on α and ε_T shown in figures 1 and 2 from the point of diamagnetic drift resonance. k_θρ_{s, D} is fixed in the analysis. The diamagnetic drift resonance occurs when ${\widehat{\omega }}_{\mathrm{r}}=2{\widehat{v}}_{\Vert }^2+{\widehat{v}}_{\perp }^2$ . The equilibrium distribution for D is $F_{0\mathrm{D}}\propto {\mathrm{e}}^{-\alpha {\widehat{v}}_{\Vert }^2}$ , and for T is $F_{0\mathrm{T}}\propto {\mathrm{e}}^{-{\widehat{v}}_{\Vert }^2}$ . When only D exists (ε_T = 0.0), the resonant parallel velocity ${\widehat{v}}_{\Vert }^{\mathrm{res}}$ can be roughly assumed as $1/\sqrt{\alpha };$ when only T exists (ε_T = 1.0), the resonant parallel velocity is around 1. When ε_T increases, the fraction of D decreases and the fraction of T increases. For α < 1 (α > 1), the increment of ε_T means the decrement (increment) of ${\widehat{v}}_{\Vert }^{\mathrm{res}}$ , resulting to the decrement (increment) of the real frequency ${\widehat{\omega }}_{\mathrm{r}}$ , as shown in figures 1(a) and (c). As for the growth rate, as concluded in [6], when the parallel temperature is higher, ITG mode can be more stabilized due to the enhancement of the Landau damping effect. Compared to the case with ε_T = 0.0, the increment of ε_T denotes that more particles' equilibrium parallel temperature varied from T_⊥ /α to T_⊥ . Thus, for α < 1 (α > 1), the increment of ε_T means the decrement (increment) of parallel temperature, which results in the increment (decrement) of the growth rate due to the decrement (increment) of the Landau damping effect, as is shown in figures 1(d) and (f). Analogously, in figure 2, for a fixed ε_T, when α increases, the resonant parallel velocity contributed from D decreases, leading to the decrement of the real frequency; on the other hand, the increment of α means the decrement of equilibrium parallel temperature for D, which leads to the increment of the growth rate due to the reduction of the Landau damping effect.

In figure 3, the real frequency and growth rate versus k_θρ_{s, D} in the CBC, ITER and CFETR cases are shown. It is found that tendencies of the dependence of (normalized) real frequency on α and ε_T in the ITER and CFETR cases are similar with those in the CBC case; for a fixed (k_θρ_{s, D}, α, ε_T), the (normalized) real frequency in three cases are of the same order. Tendencies of the dependence of (normalized) growth rate in the ITER case on α and ε_T are similar with those in the CBC case, except that the magnitude is about 30% smaller. The magnitude of growth rate in the CFETR case is close to that in the CBC case in a large range of α and ε_T, however, the tendency is not quite similar. It is seen that by varying ε_T from 0.0 to 0.5, the growth rate in the CFETR case is nearly unchanged when α = 1/2, and only about 20% increment when α = 1 and 2 in large k_θρ_{s, D} region. Analogously, dependence of the growth rate on α is not quite obvious. When α varies from 1/2 to 2, the increment of the growth rate in the CFETR case is about 50%, which is much smaller than the increments in the CBC and ITER cases. These results show that the ITG instability is the weakest in the ITER case, and is least affected by varying α and ε_T in the CFETR case.

Figure  3.  Real frequency (a)–(c) and growth rate (d)–(f) versus k_θρ_{s, D} in the CBC, ITER and CFETR cases. α = 1/2, 1 and 2, ε_T = 0.0 and 0.5. The normalized frequency c_s/R = 156 kHz (CBC), 96 kHz (ITER), 116 kHz (CFETR); here, R = R₀ + 0.5a for ITER and R = R₀ + 0.6a for CFETR.

DownLoad: Full-Size Img PowerPoint

To end this section, we make some remarks on the results about the linear frequency. In our modelling, only the temperature anisotropy of D is considered; effects of the anisotropy factor α_D and the fraction of T ε_T on the linear frequency depend on this assumption. Obviously, when the temperature anisotropy of T is also considered, the situation becomes complicated. For example, in the case with ε_T increased, though the influence of temperature anisotropy from D decreases, the influence of temperature anisotropy from T increases, and the final results depend on the competition of the temperature anisotropy between D and T. In this work, we do not focus on the analysis of this complicated situation, but emphasize the significant role of the temperature anisotropy on the linear frequency in a relatively easy way.

4.   Quasilinear transport

In this section, based on the linear results in section 3, a quasilinear transport model is adopted to study the effects of temperature anisotropy on the particle and energy fluxes driven by the ITG instability.

We begin with the relation between the driven fluxes and the driving forces [21, 29–31]

Here, η, ν = 1, 2, 3, 4, 5, L_ην is the transport coefficient, which can be explicitly written in the quasilinear theory framework as [32, 33]

It associates with the velocity integral for the weighting factors h_η(ν) and the linear frequency. Obviously, L_ην satisfies the Onsager symmetry [29, 30] automatically.

The particle and energy fluxes are generally defined as

Here, δf is the perturbed particle distribution function with the Fourier component shown in equation (5). δv_r is the radial component of the perturbed E × B drift. $\langle\cdot\rangle_{\varepsilon}$ denotes the ensemble average, which is usually defined as the average over the periodic poloidal and toroidal angle, and average over the time. Note that the phase difference between δv_r and the polarization part in δf is π/2, only the gyrocenter part in δf contributes to the transport fluxes. By substituting equation (5) into equations (12a) and (12b), and by combining the definitions in equations (10) and (11), it is easy to find

from which we can define the heat flux

Further analyses show that $\mathcal{J}_2=q_{\perp}^r / T_{\perp}$, which denotes the perpendicular heat flux, $\mathcal{J}_3=q_{\|}^r / T_{\|}$, which denotes the parallel heat flux. Note that in the above derivation, we omit the label 'j' to treat D and T equally. In the following part, we use Γ^r/n₀ and Q^r/n₀T_⊥(q^r/n₀T_⊥ ) to describe the transport level for particle and energy (heat).

In the quasilinear transport theory, the mixing length model approximation [32, 34, 35] is widely used to quantitatively analyze the turbulent transport fluxes. It reads as

then the transport coefficients can be rewritten as

Here, $\langle\cdot\rangle_{\hat{v_{\|}}}$ denotes the -integral. and denote the real frequency and growth rate normalized by ω_0D, respectively. χ_GB, defined as , is usually used to normalize the transport coefficients. According to the parameters shown in table 1, χ_GB is evaluated to be 5.83 m² s^-1 (CBC), 2.57 m² s^-1 (ITER), and 5.05 m² s^-1 (CFETR). Generally, only the spectrum which corresponds to the maximum growth rate should be considered in the calculation of transport coefficients [34], which is relevant to the assumption that the nonlinear transport may be dominated by the mode with the maximum linear growth rate. From figure 3, it is seen that the growth rates with k_θρ_{s, D} = 0.3–0.5 are comparable in large range of (α, ε_T); note that the saturation amplitude shown in equation (14) is proportional to , a smaller k_θρ_s corresponds to a larger saturation amplitude for comparable growth rates. Besides, in the CBC case, when only D exists and α = 1, the energy flux is calculated as Q ^r/n_D0T_⊥ = 5.49 m s^-1 with k_θρ_{s, D} = 0.4, from which the effective ion heat diffusivity is estimated as χ = 0.27χ_GB, which agrees well with the GK simulation result [36]. Thus, we take k_θρ_{s, D} = 0.4 in the following calculation.

4.1   Particle flux

Figure 4 shows the particle flux versus ε_T for D $({\mathrm{\Gamma }}_{\mathrm{D}}^r/n_{\mathrm{D}0})$ , T $({\mathrm{\Gamma }}_{\mathrm{T}}^r/n_{\mathrm{T}0})$ and the total flux $({\mathrm{\Gamma }}_{\mathrm{total}}^r/n_{\mathrm{e}0})$ with α varied in the CBC case. The particle fluxes for D and T are calculated according to equation (13a), and the total flux is calculated as

It is seen from figure 4(c) that the total particle flux with different (ε_T, α) is nearly zero, the absolute value is smaller than 5 × 10^-3, the error range we adopted to calculate the left-hand-side of equation (9). Note that in the derivation of equation (9), the adiabatic electron approximation and the quasi-neutrality condition are adopted. Due to the adiabatic approximation, the particle flux for electrons driven by the ITG instability is zero, then the particle flux for ions should also be zero to keep the quasi-neutrality condition. Figure 4(c) is obtained based on the linear frequency and the mixing length model approximation, it proves the accuracy of the linear frequency calculated in section 3 and provides the reliability for the calculation of the quasilinear transport fluxes in section 4.

Figure  4.  Particle flux versus ε_T with α varied in the CBC case. (a) D; (b) T and (c) total. The gray dash line denotes Γ^r = 0.

DownLoad: Full-Size Img PowerPoint

The particle fluxes for D and T are shown in figures 4(a) and (b), respectively. It is seen that direction of the particle flux for T by varying (α, ε_T) is opposite to that for D, which qualitatively agrees with equation (16) by noting that the total particle flux is zero. In detail, for a fixed ε_T, it is shown that ${\mathrm{\Gamma }}_{\mathrm{D}}^r/n_{\mathrm{D}0}<0$ ( ${\mathrm{\Gamma }}_{\mathrm{T}}^r/n_{\mathrm{T}0}>0$ ) in α < 1 region, that is, the particle flux for D(T) is inward (outward), it decreases as α increases; in α ≥ 1 region, ${\mathrm{\Gamma }}_{\mathrm{D}}^r/n_{\mathrm{D}0}>0$ ( ${\mathrm{\Gamma }}_{\mathrm{T}}^r/n_{\mathrm{T}0}<0$ ), that is, the particle flux for D(T) is outward (inward), it increases as α increases. Besides, for a fixed α, as ε_T increases, the particle flux for D(T) increases (decreases), no matter it is inward or outward. When ε_T = 0 (or ε_T = 1.0), only D (or T) exists, the particle flux for D (or T) should be zero for any α. Thus, the lines denoted by different α should intersect at (ε_T, Γ^r) = (0.0, 0.0) and (1.0, 0.0) in figures 4(a) and (b), respectively.

According to the definition in equations (10) and (13a), the particle flux can be divided into three parts: one is driven by the density gradient ($L_{11} \mathcal{X}^1$), the second is driven by the temperature gradient ($L_{12} \mathcal{X}^2+L_{13} \mathcal{X}^3$), and the third is due to the inward pinch ($L_{14} \mathcal{X}^4+L_{15} \mathcal{X}^5$). Totally speaking, as α and ε_T varied, the linear real frequency and growth rate vary, then the transport coefficients calculated from equation (15) vary, leading to the variation of the contribution from these three parts to the particle flux. The analysis about the contribution from these three parts to the particle flux may be helpful to understand the tendency of the particle flux varied with α and ε_T. As an example, we take ε_T = 0.5, α = 1/2, 1 and 2. The values are summarized in table 2. It is shown that with different α, the flux driven by the density gradient is always outward, the flux due to the inward pinch is always inward, the flux driven by the temperature gradient can be inward or outward. For a fixed α, the density gradient parts and the inward pinch parts for D and T are similar, but the temperature gradient parts for D and T are quite different, especially when α = 1/2 and 2. As α increased from 1/2 to 2, the density gradient parts for D and T increase; the inward pinch parts for D and T increase; the temperature gradient part for D(T) changes from inward (outward) to outward (inward), correspondingly, the particle flux for D(T) changes from inward (outward) to outward (inward). The similar results can also be obtained for ε_T = 0.1 to 0.4. It is indicated that the temperature gradient part plays an important role in determining the direction of particle flux for D and T.

Table  2.  Contribution from the density gradient ($\left(L_{11} \mathcal{X}^1\right)$), the temperature gradient ($L_{12} \mathcal{X}^2+L_{13} \mathcal{X}^3$), and the inward pinch ($L_{14} \mathcal{X}^4+L_{15} \mathcal{X}^5$) to the particle flux (total). The unit for the flux is m s^-1. ε_T = 0.5, and α = 1/2, 1, 2. ' ← ' and ' → ' denote the inward and outward particle flux, respectively.

	α	$L_{11} \mathcal{X}^1$	$L_{12} \mathcal{X}^2+L_{13} \mathcal{X}^3$	$L_{14} \mathcal{X}^4+L_{15} \mathcal{X}^5$	Total	Direction
	1/2	1.453	-0.720	-1.108	-0.375	←
D	1	2.543	-0.117	-2.261	0.166	→
	2	2.987	0.795	-2.900	0.882	→
	1/2	1.365	0.324	-1.315	0.374	→
T	1	2.339	-0.463	-2.039	-0.163	←
	2	2.804	-1.431	-2.256	-0.883	←

 | Show Table

DownLoad: CSV

Figure 5 shows the particle flux versus ε_T in the CBC, ITER and CFETR cases. α = 1/2, 1 and 2 are adopted. Totally speaking, tendencies of the particle flux for D and T versus (ε_T, α) in the ITER and CFETR cases are similar with those in the CBC case. In large range of chosen (ε_T, α), the particle fluxes for D and T are about 0.5 m s^-1, 0.03 m s^-1 and 0.2 m s^-1 in the CBC, ITER and CFETR cases, respectively, that is, the particle flux is the largest in the CBC case and smallest in the ITER case. These results can be qualitatively explained by combining equation (15) and figure 3. The transport coefficient shown in equation (15) is proportional to γ³ and χ_GB, both γ and χ_GB in the ITER case are about half of those in the CBC and CFETR cases, resulting that the particle flux in the ITER case is about an order of magnitude smaller. It is seen that the particle flux in the CFETR case is smaller than that in the CBC case, which is closely related to the contribution from the temperature gradient part and the inward pinch part. We take ε_T = 0.5 as an example. When α = 2, the particle fluxes for D in the CBC and CFETR case are 0.882 m s^-1 and 0.321 m s^-1, respectively; however, the temperature gradient parts in the CBC and CFETR case are 0.795 m s^-1 and -0.712 m s^-1, respectively. The bigger outward particle flux for D in the CBC case is mainly due to the positive contribution from the temperature gradient part. When α = 1/2, the particle fluxes for D in the CBC and CFETR case are -0.375 m s^-1 and -0.162 m s^-1, respectively; however, the inward pinch parts in the CBC and CFETR case are -1.108 m s^-1 and -0.460 m s^-1, respectively. The bigger inward particle flux for D in the CBC case is mainly due to the positive contribution from the inward pinch part.

Figure  5.  Particle flux versus ε_T in the (a) CBC, (b) ITER and (c) CFETR cases. α = 1/2, 1 and 2. Solid: D, dash: T, dot: total.

DownLoad: Full-Size Img PowerPoint

Generally, the inward or smaller outward particle flux means the better particle confinement. In figure 5, it is shown that for a fixed ε_T, the particle fluxes for both D and T are closer to zero when α approaches to 1 than that when α is far away from 1, indicating that the confinement is better for the isotropic case than that for the anisotropic case. On the other hand, for the anisotropic case, it is found that for a certain α < 1, larger ε_T (for example, α = 1/2 and ε_T = 0.5) is beneficial for the particle confinement; for a certain α > 1, smaller ε_T (for example, α = 2 and ε_T = 0.1) is beneficial for the particle confinement.

4.2   Energy flux

Figure 6 shows the heat and energy fluxes for D versus ε_T in the CBC case. The heat and energy fluxes are calculated according to equations (13c) and (13b), respectively. It is seen that for a fixed ε_T, both the heat and energy fluxes increase as α increases, dependences of the heat and energy flux on α are more sensitive in α < 1 region than that in α > 1 region. For ε_T = 0.5, when α increases from 1/2 to 1, the heat flux increases from about 2.3 to about 4.5 m s^-1, the increment is about 96%; when α increases from 1 to 2, the heat flux increases from about 4.5 to about 5.4 m s^-1, the increment is about 20%. As shown in equation (13c), the heat flux contains two parts, which are denoted as the perpendicular heat flux ($\mathcal{J}_2$) and the parallel heat flux ($\mathcal{J}_3$). A careful analysis shows that the perpendicular part dominates the heat flux (in our cases, the ratios of the parallel part to the total heat flux are 1%, 14% and 20% for α = 1/2, 1 and 2 with ε_T = 0.5, respectively). Furthermore, the perpendicular heat flux is dominated by the temperature gradient part ($L_{22} \mathcal{X}^2+L_{23} \mathcal{X}^3$). It is found that dependence of the heat flux on α has the similar tendency to that of the temperature gradient part.

Figure  6.  Heat (a) and energy (b) flux for D versus ε_T with α varied in the CBC case.

DownLoad: Full-Size Img PowerPoint

As shown in equation (13b), the energy flux contains the particle flux and the heat flux. According to figures 4(a) and 6(a), the particle flux is much smaller than the heat flux, thus the heat flux dominates the energy flux. It is seen that dependence of the energy flux on α is similar to that of the heat flux. In addition, for a fixed α, dependence of the heat flux on ε_T is just opposite to that of the particle flux, which leads to that the energy flux moderately decreased as ε_T increases. When ε_T varies from 0.1 to 0.5, the biggest change of energy flux occurs at α = 1, which is dropped from about 5.3 to about 4.7 m s^-1. It indicates that the energy transport can be significantly reduced in the situation where the parallel temperature dominates (α < 1), and effectively aggravated in the situation where the perpendicular temperature dominates (α > 1).

Figure 7 shows the heat and energy fluxes for T versus ε_T in the CBC case. It is seen that dependence of the heat flux for T on (α, ε_T) and its value are similar to those for D. The heat flux for T increases as α increases, and the increment is more significant in α < 1 region. Analogously, the perpendicular heat flux dominates the heat flux, and the temperature gradient part dominates the perpendicular heat flux. Unlike the situation for the heat flux, dependence of the energy flux on (α, ε_T) shows different characters. Since the particle flux for T is comparable to the heat flux for T, due to the modulation of particle flux, the energy flux is not sensitive to ε_T, and its value obviously increases (decreases) in α < 1 (α > 1) region, where the particle flux is outward (inward). It is seen that the biggest energy flux appears at α = 1. When α is far away from α = 1, the energy flux for T can be moderately suppressed by about 20%–40% in the chosen parameter ranges.

Figure  7.  Heat (a) and energy (b) flux for T versus ε_T with α varied in the CBC case.

DownLoad: Full-Size Img PowerPoint

Figure 8 shows the energy flux versus ε_T in the CBC, ITER and CFETR cases. α = 1/2, 1 and 2 are adopted. Tendencies of the energy flux for D and T versus (ε_T, α) in the ITER and CFETR cases are similar with those in the CBC case. In each case, qualitatively, dependence of the energy flux for D on α is more significant than that for T, dependences of the energy flux for both D and T on ε_T are not obvious, which are due to the modulation of the particle flux, as illustrated in the following figures 6 and 7. Quantitatively, analogous to the particle flux shown in figure 5, the energy flux in the ITER case is an order of magnitude smaller than that in the CBC case, which is due to the smaller growth rate and χ_GB, as is discussed following figure 5. The energy flux in the CFETR case is about half of magnitude smaller than that in the CBC case. Note that the energy flux can also be divided into three parts, which are driven by the density gradient, the temperature gradient and the inward pinch. As an example, in table 3, we show the contribution from these three parts to the energy flux for D with ε_T = 0.5 and α = 1/2, 1 and 2. It is shown that in the CBC case, the density gradient part and the inward pinch are mostly cancelled, the value of the total energy flux is close to that of the temperature gradient part. In the CFETR case, values of these three parts are smaller than those in the CBC case. Though the transport coefficients described by equation (15) are close in CBC and CFETR cases, the fluxes are not close, since the the fluxes are calculated by the product of transport coefficients and driving forces ($\mathcal{X}^\nu$). For example, in our modelling, $\mathcal{X}^2=\mathcal{X}^3=\kappa_{\mathrm{T}} / R$, the bigger R results in the smaller energy flux in the CFETR case.

Figure  8.  Energy flux for D (solid) and T (dash) versus ε_T. (a) CBC; (b) ITER; (c) CFETR.

DownLoad: Full-Size Img PowerPoint

Table  3.  Contribution from the density gradient part, the temperature gradient part and the inward pinch part to the energy flux for D in the CBC and CFETR cases. The unit for the flux is m s^-1. ε_T = 0.5, and α = 1/2, 1, 2.

	α	Density gradient	Temperature gradient	Inward pinch	Total
	1/2	2.473	1.418	-2.331	1.560
CBC	1	3.777	5.208	-4.287	4.699
	2	3.867	7.441	-4.825	6.483
	1/2	1.807	0.320	-0.908	1.220
CFETR	1	2.205	1.372	-1.294	2.283
	2	2.109	1.911	-1.339	2.681

 | Show Table

DownLoad: CSV

Besides, it is seen that the optimal confinement for the energy of D and T appears at α = 1/2, regardless of the value of ε_T.

5.   Conclusion and discussion

In conclusion, the particle and energy transports driven by the ITG instability are investigated in an anisotropic D-T plasma.

At first, the linear dispersion relation is derived in the local kinetic limit. Three sets of typical parameters, denoted as CBC, ITER and CFETR case, respectively, are chosen to numerically solve the dispersion relation. Dependences of the linear mode frequency on parameters, such as the temperature anisotropy factor of D α, the fraction of T particle ε_T and the poloidal spectrum k_θρ_{s, D}, are studied. Based on the linear results, the quasilinear transport coefficients are derived by adopting the mixing length model approximation, and the particle and energy fluxes for D and T are evaluated.

The main results are concluded as follows:

(1) Dependences of the linear frequency and the quasilinear transports on the parameters in the three cases show the similar tendency. Values of the real frequency in the three cases are nearly the same, the growth rate is the smallest in the ITER case, and largest in the CBC case. The quasilinear particle and energy transports are the largest in the CBC case and the smallest in the ITER case, which indicates the better particle and energy confinement in large devices.

(2) As α increases, the real frequency decreases, while the growth rate increases; in α > 1 (α < 1) region, the real frequency increases (decreases) as ε_T increases, while the growth rate decreases (increases). The ITG instability may be significantly stabilized by controlling α < 1.

(3) Due to the adiabatic electron approximation and the quasi-neutrality condition, the total particle flux for ions is found to be zero; this result proves the accuracy for the calculation of the linear frequency and provides the reliability for the calculation of the quasilinear transport fluxes. The direction of the particle flux for D is opposite to that of T; the temperature gradient driven particle flux plays an important role in determining the direction of the particle flux for D and T. It is shown that the particle confinement is better for the isotropic case than that for the anisotropic case when ε_T is fixed. On the other hand, for the anisotropic case, it is found that for a certain α < 1 (α > 1), larger ε_T (smaller ε_T) is beneficial for the particle confinement.

(4) As α increases, the heat fluxes for D and T increase, the increment is more sensitive in α < 1 region than that in α > 1 region. Due to the modification of the particle flux, dependences of the energy flux for D and T on α show different tendencies. As α increases, the energy flux for D increases, while the energy flux for T increases in α < 1 region and then decreases in α > 1 region. The largest energy flux for T appears at α = 1; as the deviation of α is larger from 1, the energy flux becomes smaller. It is concluded that smaller α is beneficial for energy confinement.

(5) Totally speaking, in the ITG-dominated regime, the better particle and energy confinement occur by choosing smaller α and larger ε_T.

To the end, it should be pointed out that a significant source of anisotropy in high performance machines is NBI, which may also produce a sheared flow in the plasma. Consideration of the consistency between temperature anisotropy and sheared flows may be helpful for a better estimation of turbulent transport level in devices such as ITER and CFETR. This issue will be left as a future work.

Appendix. Derivation of the perturbed distribution function from the traditional GK theory.

In this appendix, we give a brief derivation of the perturbed distribution function shown in equation (3) based on the traditional GK theory [18].

In terms of the gyrocenter (GY, following the terminology of [18], we denote the guiding center (GC) in perturbed fields as gyrocenters) coordinates $\overline{\mathit{Z}}=(\overline{\mathit{X}},{\overline{v}}_{\Vert },\overline{\mu },\overline{\xi })$ , with $\overline{\mathit{X}}$ the GY position, ${\overline{v}}_{\Vert }$ the parallel velocity, $\overline{\mu }$ the magnetic moment and $\overline{\xi }$ the gyro-angle, the GK equation is written as

where $\overline{F}(\overline{\mathit{Z}};t)$ is the GY distribution function. Linearizing equation (A.1) yields

where

Here, the subscript '0' and '1' denote the unperturbed and perturbed quantities, respectively. $\langle\delta \phi\rangle$ denotes the gyro-average of δϕ. with .

By introducing the coordinate transformation $(\overline{\mathit{X}},{\overline{v}}_{\Vert },\overline{\mu })\to (\overline{\mathit{X}},\overline{\mu },E)$ with $E=\frac{1}{2}m{\overline{v}}_{\Vert }^2+\overline{\mu }B(\overline{\mathit{X}})$ , equation (A.2) can be rewritten in terms of $(\overline{\mathit{X}},\overline{\mu },E)$ as

By using the identities and $\dot{\bar{X}}_1 \cdot \bar{\mu} \nabla B+m \bar{v}_{\|} \dot{\bar{v}}_{\|, 1}=-\dot{\boldsymbol{\bar{X}}}_0 \cdot e \nabla\langle\delta \phi\rangle$, equation (1.5) can be reduced as

We further introduce the adiabatic transformation $\bar{F}_1=\bar{H}_1+e\langle\delta \phi\rangle \partial_E \bar{F}_0$, where denotes the non-adiabatic part of the perturbed distribution function, then (1.6) can be transformed to be

We then solve equation (1.7) in the Fourier space. Note that equation (1) shows the Fourier transformation for a function with respect to the particle position r. We introduce the GC coordinates Z = ( X, v_‖,μ, ξ ), where X, v_‖,μ and ξ denote the GC position, parallel velocity, magnetic moment and gyro-angle, respectively. By using r = X + ρ(μ, ξ ), where ρ is the gyro-radius, equation (1) can be rewritten in the GC coordinates as

where J_n is the nth order Bessel function with argument λ = k_⊥ ρ. In the GK theory, the perturbed quantities shown in equation (1.7) are directly evaluated at $\overline{\mathit{Z}}$ , namely, $\delta A{|}_{\mathit{Z}=\overline{\mathit{Z}}}$ . Easily, equation (1.7) is solved in the Fourier space as

where J₀ originates from the gyro-average of δϕ. Note that the parallel dynamics has been neglected for simplicity.

In the GK theory, the GC distribution function F( Z) can be obtained from the GY distribution function $\overline{F}(\overline{\mathit{Z}})$ by using the pull-back transformation [37, 38], which can be written in the linear version as

Here, G₁'s denote the generating vector fields, which are explicitly expressed as [18, 33]

Here, S₁ is the 1st order gauge function, which is calculated as

where $\widetilde{\delta \phi}=(\delta \phi-\langle\delta \phi\rangle)$. Then the generating vector fields are calculated as

The particle distribution function f(r) can be obtained from the GC distribution function F(X) by noting the transformation X = r - ρ, namely, f(r) = F(r - ρ) = e^-iρ·∇F(r). By using the identities

the gyro-averaged particle distribution function can be obtained according to equations (1.9)–(1.13d) as

which agrees with equation (3). We mention that the key points in the derivation are (i) the derivation of the linear GK equation and (ii) the transformation among the particle, GC and GY distribution functions.