⇦Chinese Physics Letters, 2019, Vol. 36, No. 4, Article code 045201 Impact of Sheath Boundary Conditions and Magnetic Flutter on Evolution and Distribution of Transient Particle and Heat Fluxes in the Edge-Localized Mode Burst by Experimental Advanced Superconducting Tokamak Simulation * Yan-Bin Wu (吴彦斌)1,2,3, Tian-Yang Xia (夏天阳)3**, Fang-Chuan Zhong (钟方川)1**, Bin Gui (桂彬)3, EAST Team Affiliations 1College of Science, Donghua University, Shanghai 201620 2School of Physics and Electronic Engineering, Anqing Normal University, Anqing 246133 3Institute of Plasma Physics, Chinese Academy of Sciences, Hefei 230031 Received 7 November 2018, online 23 March 2019 *Supported by the National Key R&D Program of China under Grant Nos 2017YFE0301100, 2017YFE0301101, 2017YFE0301104 and 2014GB106001, the National Natural Science Foundation of China under Grant Nos 11275047, 11675217, 11505236 and 11405215, the Youth Innovation Promotion Association Chinese Academy of Sciences under Grant No 2017479, and the K. C. Wong Education Foundation.
**Corresponding author. Email: xiaty@ipp.ac.cn; fczhong@dhu.edu.cn Citation Text: Wu Y B, Xia T Y, Zhong F C, Gui B and EAST Team et al 2019 Chin. Phys. Lett. 36 045201    Abstract To study the evolution and distribution of the transient particle and heat fluxes during the edge-localized modes (ELMs) burst on the experimental advanced superconducting tokamak (EAST), the BOUT$^{++}$ six-field two-fluid model with sheath boundary conditions (SBCs) and magnetic flutter terms in the parallel thermal conduction is used to simulate the evolution of the profiles and growing process of the fluxes at divertor targets. Although SBCs hardly play a role in the linear phase, in the nonlinear phase both SBCs and magnetic flutter can change the dominant toroidal mode. SBCs are able to broaden the frequency distribution of the turbulence. The magnetic flutter increases the ELM size from 2.8% to 8.4%, and it doubles the amplitudes of the radial heat and particle transport coefficients at outer midplane (OMP), at around 1.0 m$^{2}$s$^{-1}$. It is then able to increase the particle and heat flux at the divertor targets and to broaden the radial distribution of the parallel heat flux towards the targets. DOI:10.1088/0256-307X/36/4/045201 PACS:52.55.Fa, 52.55.Dy, 52.55.Rk © 2019 Chinese Physics Society Article Text The H-mode scenario is now considered to be one of the most promising candidates to achieve magnetic confinement fusion in a tokamak. The steep pressure gradient and associated bootstrap current in the H-mode can destabilize the peeling-ballooning (P-B) modes, which may lead to the repetitive unstable events, which is named the edge localized mode (ELM).[1-3] The excessive heat loads on divertor targets during ELM bursts may cause unacceptable erosion of the materials, and then the plasma facing components (PFCs) will be severely damaged.[4] Therefore, to achieve high performance in future tokamaks, it is essential to understand the ELM dynamics and the induced heat flux. BOUT$^{++}$ is a new and highly adaptable, object-oriented $C++$ framework to perform parallel plasma fluid simulations. It is flexible in solving an arbitrary number of equations in 3D curvilinear coordinates using finite-difference methods.[5] Good agreement for an ideal magneto-hydrodynamic (MHD) model is found between BOUT$^{++}$ and other linear MHD eigenvalue codes, such as ELITE and GATO.[3] It has been successfully used to simulate the crash dynamics of ELMs in tokamaks, such as DIII-D, NSTX, C-Mod and the experimental advanced superconducting tokamak (EAST).[3,6,7] For example, as discussed in previous works,[8-10] the electromagnetic turbulence of electrons plays a dominant role in the evolution of divertor heat flux during ELM burst in EAST, C-Mod and DIII-D. Linear analysis shows that the main driving mechanisms are the resistive ballooning mode and drift-Alfven wave (DAWs). The profiles of ion saturation current density ($j_{\rm s}$) are simulated and are comparable to those of experiments by divertor probes, including the peak value and the radial distribution of $j_{\rm s}$.[11] For the simulations of the transient particle and heat fluxes induced by ELMs, the sheath boundary conditions (SBCs) are very important for research of tokamak edge physics and they need to be considered near the divertor targets. SBCs will affect the evolutions of the perturbations, which may have influence on the turbulence in the pedestal region and affect the parallel heat flux during the ELM events.[12,13] Moreover, the magnetic flutter in the thermal conduction is found to be able to increase the total energy loss,[14] which can have a strong impact on the energy deposition on the divertor targets.[15] Therefore, in this study we try to find the details of the impact of SBCs and magnetic flutter on ELMs to understand the physical mechanism of the transient heat flux distributions during an ELM crash. To discover the effects of SBCs and magnetic flutter on particle transport and heat flux evolution and distribution during an ELM event in EAST, three different cases using the BOUT$^{++}$ code are simulated and compared. The six-field model for the plasma edge with a tokamak configuration is constituted of six evolving equations based on Ref.  [14], which are derived in the drift ordering with the flute reduction assumption, $$\begin{alignat}{1} \frac{\partial}{\partial t}\varpi =\,&-\frac{1}{B_{0}}{\boldsymbol b}\times \nabla _{\bot}{\it \Phi} \cdot \nabla \varpi-\frac{1}{B_{0}}{\boldsymbol b}\times \nabla _{\bot}{\it \Phi} \cdot \nabla \varpi_{0}\\ &+B_{0}^{2}\nabla _{\parallel}\Big(\frac{J_{\parallel}}{B_{0}} \Big)+2{\boldsymbol b}\times {\boldsymbol\kappa}\cdot \nabla p\\ &-\frac{1}{2{\it \Omega} _{\rm i}}\Big[\frac{1}{B_{0}}{\boldsymbol b}\times \nabla P_{\rm i}\cdot \nabla (\nabla _{\bot}^{2}{\it \Phi})-Z_{\rm i}eB_{0}{\boldsymbol b}\\ &\times \nabla n_{\rm i}\cdot \nabla \Big(\frac{\nabla _{\bot}{\it \Phi}}{B_{0}} \Big)^2 \Big] \end{alignat} $$ $$\begin{alignat}{1} &+\frac{1}{2{\it \Omega} _{\rm i}}\Big[\frac{1}{B_{0}}{\boldsymbol b}\times \nabla {\it \Phi} \cdot \nabla (\nabla _{\bot}^{2}P_{\rm i})\\ &-\nabla _{\bot}^{2}\Big(\frac{1}{B_{0}}{\boldsymbol b}\times \nabla {\it \Phi} \cdot \nabla P_{\rm i} \Big) \Big]+\mu_{\parallel {\rm i}}\nabla _{\parallel 0}^{2}\varpi,~~ \tag {1} \end{alignat} $$ $$\begin{alignat}{1} \frac{\partial}{\partial t}n_{\rm i}=\,&-\frac{1}{B_{0}}{\boldsymbol b}\times \nabla _{\bot}{\it \Phi} \cdot \nabla n_{\rm i}-\frac{2n_{\rm i}}{B_{0}}{\boldsymbol b}\times {\boldsymbol\kappa}\cdot \nabla {\it \Phi} \\ &-\frac{2}{Z_{\rm i}eB_{0}}{\boldsymbol b}\times {\boldsymbol\kappa}\cdot \nabla P_{\rm i}-n_{\rm i}B_{0}\nabla _{\parallel}\Big(\frac{V_{\parallel {\rm i}}}{B_{0}} \Big),~~ \tag {2} \end{alignat} $$ $$\begin{alignat}{1} \frac{\partial}{\partial t}V_{\parallel {\rm i}}=\,&-\frac{1}{B_{0}}{\boldsymbol b}\times \nabla _{\bot}{\it \Phi} \cdot \nabla V_{\parallel {\rm i}}-\frac{1}{{m_{\rm i}n}_{i0}}{\boldsymbol b}\cdot \nabla P,~~ \tag {3} \end{alignat} $$ $$\begin{alignat}{1} \frac{\partial}{\partial t}T_{\rm i}=\,&-\frac{1}{B_{0}}{\boldsymbol b}\times \nabla _{\bot}{\it \Phi} \cdot \nabla T_{\rm i}-\frac{2}{3}T_{\rm i}\Big[ \Big(\frac{2}{B_{0}}{\boldsymbol b}\times {\boldsymbol\kappa} \Big)\\ &\cdot \Big(\nabla {\it \Phi} +\frac{1}{Z_{\rm i}en_{i0}}\nabla P_{\rm i}+\frac{5}{2}\frac{k_{\rm B}}{Z_{\rm i}e}\nabla T_{\rm i} \Big)\\ &+B_{0}\nabla _{\parallel}\Big(\frac{V_{\parallel {\rm i}}}{B_{0}} \Big) \Big]+\frac{2}{3n_{i0}k_{\rm B}}\nabla _{\parallel}(\kappa_{\parallel {\rm i}}\nabla _{\parallel}T_{\rm i})\\ &+\frac{2m_{\rm e}}{m_{\rm i}}\frac{Z_{\rm i}}{\tau_{\rm e}}(T_{\rm e}-T_{\rm i}),~~ \tag {4} \end{alignat} $$ $$\begin{alignat}{1} \frac{\partial}{\partial t}T_{\rm e}=\,&-\frac{1}{B_{0}}{\boldsymbol b}\times \nabla _{\bot}{\it \Phi} \cdot \nabla T_{\rm e}-\frac{2}{3}T_{\rm e}\Big[ \Big(\frac{2}{B_{0}}{\boldsymbol b}\times {\boldsymbol\kappa} \Big)\\ &\cdot \Big(\nabla {\it \Phi} -\frac{1}{en_{e0}}\nabla P_{\rm e}-\frac{5}{2}\frac{k_{\rm B}}{e}\nabla T_{\rm e} \Big)\\ &+B_{0}\nabla _{\parallel}\Big(\frac{V_{\parallel {\rm e}}}{B_{0}} \Big) \Big]+0.71\frac{2T_{\rm e}}{3en_{e0}}B_{0}\nabla _{\parallel}\Big( \frac{J_{\parallel}}{B_{0}} \Big)\\ &+\frac{2}{3n_{e0}k_{\rm B}}\nabla _{\parallel}(\kappa_{\parallel {\rm e}}\nabla _{\parallel}T_{\rm e})-\frac{2m_{\rm e}}{m_{\rm i}}\frac{1}{\tau_{\rm e}}(T_{\rm e}-T_{\rm i})\\ &+\frac{2}{3n_{e0}k_{\rm B}}\eta_{\parallel}J_{\parallel}^{2},~~ \tag {5} \end{alignat} $$ $$\begin{alignat}{1} \frac{\partial}{\partial t}A_{\parallel}=\,&-\nabla _{\parallel}\phi +\frac{\eta}{\mu_{0}}\nabla _{\bot}^{2}A_{\parallel}+\frac{1}{en_{e0}}\nabla _{\parallel}P_{\rm e}\\ &+0.71\frac{k_{\rm B}}{eB_{0}}\nabla _{\parallel}T_{\rm e}-\frac{\eta_{H}}{\mu_{0}}\nabla _{\bot}^{4}A_{\parallel}.~~ \tag {6} \end{alignat} $$ In this model, all of the variables—such as ion density $n_{\rm i}$, ion and electron temperature $T_{\rm j}$—can be written as $F=F_{0}+F_{1}$, where $F_{0}$ represents for the equilibrium part of arbitrary field quantity, and $F_{1}$ is the perturbed component. Here vorticity $\varpi=n_{i0}\frac{m_{\rm i}}{B_{0}}(\nabla _{\bot}^{2}\phi +\frac{1}{n_{i0}}\nabla _{\bot}\phi \cdot \nabla _{\bot}n_{i0}+\frac{1}{n_{i0}Z_{\rm i}e}\nabla _{\bot}^{2}p_{i1} )$, parallel current $J_{\parallel}=J_{\parallel 0}-\frac{1}{\mu_{0}}\nabla _{\bot}^{2}A_{\parallel}$, parallel electron velocity $V_{\parallel {\rm e}}=V_{\parallel {\rm i}}+\frac{1}{\mu_{0}Z_{\rm i}en_{\rm i}}\nabla _{\bot}^{2}A_{\parallel}$. $A_{\parallel}$ is the perturbed parallel vector potential, and ${\boldsymbol b}$ is the unit vector of the total magnetic field. The curvature vector of the magnetic field is defined as ${\boldsymbol\kappa}={\boldsymbol b}_{0}\cdot{\nabla {\boldsymbol b}}_{0}$. Pressure $P_{\rm j}=n_{\rm j}k_{\rm B}T_{\rm j}$, ${\it \Phi} ={\it \Phi} _{0}+\phi $ is the total electric potential, and ${\it \Omega} _{\rm i}=Z_{\rm i}eB/m_{\rm i}$ is the ion gyro motion frequency. In this study, the equilibrium of a typical ELMy H-mode discharge #62585 with plasma current $I_{\rm p}=600$ kA at simulation time $t=3800$ ms is used. The geometry and profiles used in the simulations are shown in Ref.  [12] and the equilibrium is reconstructed by kinetic EFIT code[16] with the measured electron density and temperature profiles. The pressure profile is calculated from the density and temperature measurements with the assumptions of $T_{\rm i}=T_{\rm e}$ in the edge region, and also the contributions from fast ions. The simulation domain is chosen from the range $0.85\le {\it \Psi}_{\rm N} \le1.05$, which covers both the pedestal region and the scrape-off layer (SOL), where ${\it \Psi}_{\rm N}$ is the normalized magnetic flux. Therefore, all of the density and temperature profiles are measured inside the separatrix and are extrapolated into the SOL. The quasi-neutral condition is applied in this model to keep the self-consistency of the equilibrium because the fast ions are not taken into account in the six-field two-fluid model. Due to the lack of the rotation measurements in this shot, the radial electric field is assumed to be balanced with the ion diamagnetic flow. In this work, only one-fifth of the torus is simulated considering the computational efficiency.
cpl-36-4-045201-fig1.png
Fig. 1. The growth rate of $T_{\rm e}$ with different toroidal mode numbers $n$ in the linear phase for different cases. The black triangle curve represents the case with SBCs and flutter—the dominant toroidal mode number is $n=20$. The green square curve is the result with neither SBCs nor flutter. The red diamond curve, which is overlapped by the green curve, shows the case with SBCs but without flutter. The dominant toroidal mode number is $N=25$ in these two cases.
The SBCs are very important for the research of tokamak edge physics. For example, the dynamics of the effects on the radial electric field by thermal sheath and radial-frequency sheath have been studied in Ref.  [17]. An electrostatic sheath will form at any plasma boundary and acts to filter all but the high energy electrons, while it attracts ions, controls the particle and energy flux leaving the plasma. Therefore in the SOL and private flux region, the divertor plate boundary conditions[13] are $$\begin{align} V_{\rm j}=\,&c_{\rm se}=\sqrt \frac{T_{\rm i}+T_{\rm e}}{M_{\rm j}},~~ \tag {7} \end{align} $$ $$\begin{align} j_{\parallel}^{\rm e}=\,&N_{\rm i}e\Big[c_{\rm se}-\frac{v_{T_{\rm e}}}{2\sqrt \pi}\exp\Big(-\frac{e\phi}{T_{\rm e}} \Big) \Big],~~ \tag {8} \end{align} $$ $$\begin{align} q_{\rm se}=\,&-\kappa_{\parallel {\rm e}}\partial_{\parallel}T_{\rm e}=\gamma_{\rm e}N_{\rm i}T_{\rm e}c_{\rm se},~~ \tag {9} \end{align} $$ $$\begin{align} q_{\rm si}=\,&-\kappa_{\parallel {\rm i}}\partial_{\parallel}T_{\rm i}=\gamma_{\rm i}N_{\rm i}T_{\rm i}c_{\rm se},~~ \tag {10} \end{align} $$ $$\begin{align} \partial_{\parallel}\varpi =\,&0,~~ \tag {11} \end{align} $$ $$\begin{align} \partial_{\parallel}N_{\rm i}=\,&0,~~ \tag {12} \end{align} $$ where $c_{\rm se}$ is the particle velocity, $q_{\rm s}$ is the total heat flux at the sheath edge, $\gamma_{\rm i}\cong 2.5$ and $\gamma_{\rm e}\cong 7$ are sheath energy transmission factors. At the sheath entrance, the sheath potential is applied as the boundary condition. It then quickly diffuses in the parallel direction and also generates the radial electric field. This generates an ${\boldsymbol E}\times {\boldsymbol B}$ drift flow and the velocity difference of the flows between the different flux surface results in vorticity. Diffusion also smooths the local radial electric field. The vorticity and the electrostatic potential are in a steady state when it comes to saturation. In the following, the thermal SBCs are implemented at the divertor targets in the simulations. The results with SBCs are labeled as the w/SBC curves in the relevant figures. For comparison, the results with the Dirichlet boundary condition are labeled by w/o SBC (w, with; w/o, without). The Dirichlet boundary condition means that the evolving variables at the boundaries of divertors are fixed at zero, which is the default in this code. The magnetic flutter has strong impact on the distribution of heat fluxes towards divertor targets. Notice that the thermal conduction in Eqs. (4) and (5) should be written as $$\begin{alignat}{1} &\nabla _{\parallel}(\kappa_{\parallel j}\nabla _{\parallel}T_{\rm j})\\ =\,&\nabla _{\parallel 0}(\kappa_{\parallel j}\nabla _{\parallel 0}T_{j1})\\ &+{\boldsymbol b}_{0}\times \nabla \psi \cdot \nabla (\kappa_{\parallel j}\nabla _{\parallel 0}T_{j1})\\ &+\nabla _{\parallel 0}(\kappa_{\parallel j}{\boldsymbol b}_{0}\times \nabla \psi \cdot \nabla T_{\rm j})\\ &+{\boldsymbol b}_{0}\times \nabla \psi \cdot \nabla (\kappa_{\parallel j}{\boldsymbol b}_{0}\times \nabla \psi \cdot \nabla T_{\rm j} ),~~ \tag {13} \end{alignat} $$ if all of the terms are kept. The second, third and fourth terms are generated due to the direction of perturbed magnetic field, and they are called the magnetic flutter terms in the parallel thermal conduction. The first term on the right-hand side of Eq. (13) is the expression for the case without magnetic flutter. We also discuss the effects induced by these magnetic flutter terms. In the figures, they are labeled as w/flutter. To find the effects of SBCs and magnetic flutter in thermal conduction on the edge instabilities, it is necessary to analyze the growth rate of $T_{\rm e}$ at the peak gradient in the outer midplane (OMP) with different toroidal mode numbers $n$. Resistive ballooning mode and drift-Alfvén wave (DAWs) are two dominant instabilities for the equilibrium and they play important roles in driving ELMs. The growth rate of $T_{\rm e}$ with different toroidal mode numbers $n$ in the linear phase of different cases is shown in Fig. 1. The black triangle curve represents the case with SBCs and flutter, and the peak growth rate is located at the toroidal mode number $n=20$. The green square curve is the result with neither SBCs nor flutter. The red diamond curve shows the case with SBCs but without flutter. It is clear that the red curve is overlapped by the green curve, which means that the linear growth rates are identical. Here $N=25$ is the dominant toroidal mode number in both the cases. Figure 1 shows clearly that SBCs hardly change the linear behavior of the resistive-ballooning modes. This happens because the linear instabilities driven by the resistive-ballooning modes and DAWs only arise in the peak gradient region and do not spread to divertor region. The boundary conditions at targets are unable to affect the solutions in the simulation domain. The amplitudes of $T_{\rm e}$ in these two cases are slightly smaller than the case with the magnetic flutter, which indicates that the magnetic flutter has a contribution to the edge instabilities. The magnetic flutter increases the linear growth rate when $n \leqslant 25$ and decreases it when $n \geqslant30$ for this equilibrium. The time evolutions of $T_{\rm e}$ with different toroidal mode numbers $n$ in the nonlinear phase are plotted in Fig. 2. The time is selected from $t=0.15$ ms to $t=0.45$ ms which covers the whole nonlinear phase for all the three cases. For the case where SBCs and flutter are activated, the simulated toroidal mode with $n=5$ dominates the whole nonlinear phase. For the case without flutter, the simulated toroidal modes with $n=5$ and $n=10$ are almost at the similar amplitude. For the case with neither SBCs nor flutter, the $n=10$ mode is dominant at the early stage, and it is changed to $n=5$ when the fluctuations get saturated.
cpl-36-4-045201-fig2.png
Fig. 2. The time evolution of the toroidal mode numbers $n$ in the nonlinear phases for the three cases: (a) with SBCs and flutter, (b) with SBCs but without flutter, and (c) with neither SBCs nor flutter.
cpl-36-4-045201-fig3.png
Fig. 3. The frequency spectra in the radial direction for the three cases: (a) with SBCs and flutter, (b) with SBCs but without flutter, and (c) with neither SBCs nor flutter.
Figure 3 shows the radial distribution of frequencies for all the three cases. Notice that the signals are all calculated from $T_{\rm e}$ because the simulated divertor heat fluxes mainly come from the electron channel. For the case in Fig. 3(a), the perturbations of $T_{\rm e}$ spread widely in the radial direction with the same frequency and almost go through all the simulation domain inside the separatrix. The fluctuations within this frequency, which exhibit the characteristic of the coherent mode, can reach the separatrix. The heat flux is then able to go into SOL and it then arrives the divertor targets. The peak value of frequency appears at 32 kHz and 95 kHz. For the cases without flutter, similar frequency behavior is found in Figs. 3(b) and 3(c). For these two cases, the amplitude of high frequency is reduced and the peak value of frequency decreases to 10–20 kHz. The radial extensions of the coherent mode are slightly narrower, from the inner boundary to ${\it \Psi}_{\rm N}=1.0$, compared to ${\it \Psi}_{\rm N}= 1.01$ in Fig. 3(a). Although the peak fluctuations are greater inside the pedestal region, they still spread to the separatrix. The frequency distribution is slightly wider in the case with SBCs shown in Fig. 3(b), in the range of 10–20 kHz, compared to 10–15 kHz without SBCs in Fig. 3(c). In conclusion, the magnetic flutter terms in thermal conduction change the peak position of the fluctuations. However, while SBC has less impact on frequency, it does broaden the frequency distribution.
cpl-36-4-045201-fig4.png
Fig. 4. (a) The time evolution of the growth of the pressure perturbation. (b) The time evolution of the ELM size. The black, red, green curves stand for the cases with SBCs and flutter, with SBCs but without flutter, and with neither SBCs nor flutter, respectively. The magnetic flutter leads to larger energy loss but only slightly increases the growth of the perturbation at the linear phase. The dashed-horizontal lines represent the mean ELM size at the saturation of the nonlinear phase.
cpl-36-4-045201-fig5.png
Fig. 5. The radial profiles of transport coefficients of (a) particle flux, (b) electron heat flux and (c) ion heat flux at the OMP. The black, red, green curves stand for the cases with SBCs and flutter, with SBCs but without flutter, and with neither SBCs nor flutter, respectively.
The growing of the pressure fluctuation is plotted in Fig. 4(a). Compared with the cases, both the green and red curves have a similar gradient at the linear phase, which are lower than the black curve. Therefore, the introduction of the magnetic flutter terms in parallel thermal conduction only slightly increases the linear growth rate, which is consistent with the results in Fig. 1. The fluctuation with flutter will reach its saturation value at $t=0.07$ ms and will be much later to $t=0.12$ ms if there is no flutter. The saturated pressure fluctuation of the case with flutter is around 15% larger than that of the case without flutter. During this period, the difference of the energy loss between the two cases is enlarged because the larger fluctuations allows more energy to get out of the pedestal region. Figure 4(b) shows the time evolution of the total energy loss, or ELM size, which is defined in Ref.  [18]. It is clear that at the saturation of the nonlinear phase, the mean ELM size with the flutter is three times larger than the case without flutter, from 2.8% to 8.4%, as shown by the dashed horizontal lines. This result is consistent with Ref.  [12].
cpl-36-4-045201-fig6.png
Fig. 6. Radial profiles of the (a) particle flux and (b) electron heat flux at the upper outboard target in the saturation phase of the simulations. The $x$-axis $R-R_{\rm sep}$ is the distance from the strike point. The black, red, green curves stand for the cases with SBCs and flutter, with SBCs but without flutter, and with neither SBCs nor flutter, respectively.
Within the simulations, the radial heat and particle transport coefficients are easy to calculate based on the definition in Ref.  [19]. The convective components of the fluxes are also dominant and the conductive component is much smaller, similarly to Ref.  [14]. Therefore, the cross-field fluxes discussed here are just induced by the ${\boldsymbol E}\mathbf{\times}{\boldsymbol B}$ drift. Figure 5 shows the radial profiles of the relative transport coefficients at the OMP. The black, red, green curves stand for the cases with SBCs and flutter, with SBCs but without flutter, and with neither SBCs nor flutter, respectively. The coefficients are plotted by the time average during the saturation period. The black curves stand for the period from 0.4004–0.4433 ms, red for 0.5434–0.5863 ms, and green for 0.4147–0.4576 ms. For the cases without flutter, as shown by the red and green curves, the amplitudes of all the three transport coefficients are at the same amplitude, around 0.5 m$^{2}$s$^{-1}$. However, when the flutter is considered, the amplitudes of all the three transport coefficients are doubled, at around 1.0 m$^{2}$s$^{-1}$, as shown by the black curves. This means that the radial transport in the pedestal region is enhanced remarkably by the magnetic flutter terms in the thermal conduction. The transport coefficients have peaks because the gradients of density and temperature have a zero point near ${\it \Psi}_{\rm N} =0.89$.
cpl-36-4-045201-fig7.png
Fig. 7. The evolution of (a) ion particle flux and (b) electron heat flux at the upper outboard target. The black, red, green curves stand for the cases with SBCs and flutter, with SBCs but without flutter, and with neither SBCs nor flutter, respectively.
Figure 6 shows the profiles of the particle and parallel heat fluxes towards the upper outboard (UO) divertor target during the ELM event for the simulated discharge. The $x$-axis $R-R_{\rm sep}$ is the distance from the strike point. To be consistent with the previous legend in Fig. 5, the differently colored curves represent different cases. Notice that the fluxes at the UO divertor target are saturated during the late time of the simulations, which may be inferred from the evolution of ion density flux and electron heat flux shown in Fig. 7. The fluxes are calculated by the time average during the same period shown in Fig. 5. The parallel heat fluxes towards the divertors during the ELM burst event are dominated by electrons, and the peak of the electron divertor heat flux is around 5.5 MW/m$^{2}$, nearly two orders of magnitude larger than that of ions in the simulations. Therefore, the ion heat flux at the divertor is neglected. The heat fluxes at the inner target plate, in spite of ion or electron, are far from saturation because the time scale for ion heat flux towards divertor targets is much larger than electrons, which is far beyond the range of this model. Due to the enhancement of the radial transport, which has been proven in Fig. 5, more particles can be transported from the pedestal region into SOL, and finally arrive at the divertor targets. Therefore, the magnetic flutter terms in thermal conduction are able to increase the particle flux and, consequently, the electron heat flux at the divertor targets, which can be seen from the black curves in Figs. 6(a) and 6(b). Furthermore, the flutter can broaden the radial distributions of parallel heat flux to the targets, and the heat flux width grows from 0.473 cm to 0.493 cm, as the horizontal dashed lines show in Fig. 6(b). Compared with the Dirichlet boundary conditions, the SBCs play important roles in the divertor region. Without the SBCs, the electron heat flux at UO divertor target is much smaller and oscillates near zero, as shown by the green curve in Fig. 6(b). In summary, to find the effects of SBCs and magnetic flutter in parallel thermal conduction on the distribution and the evolution of particle and heat fluxes induced by ELM, simulations in an ELMy H-mode discharge #62585 in EAST are carried out within a six-field two-fluid module in the BOUT$^{++}$ framework. SBCs hardly change the behavior in the linear phase, such as growth rate and dominant mode number. Meanwhile, for the nonlinear phase of the simulations, the SBCs and flutter change the toroidal modes. The magnetic flutter terms in thermal conduction change the peak position of the fluctuations of the edge turbulence. However, SBC has few impacts on frequency but it does broaden the frequency distribution. The pressure fluctuation with flutter reaches its saturation earlier and the saturated value is around 15% larger than the case without flutter. Therefore, the ELM size is three times larger than the case without flutter, from 2.8% to 8.4%. The flutter enhances the radial transport in the pedestal region remarkably. Consequently, the amplitudes of all the radial transport coefficients are doubled, at around 1.0 m$^{2}$s$^{-1}$. The transient heat fluxes on the divertors during the ELM burst event is dominated by electrons, and the peak of the electron divertor heat flux is 5.5 MW/m$^{2}$. The flutter is able to increase the particle and heat flux at the divertor targets and it broadens the radial distribution of parallel heat flux to the targets. Numerical computations were performed on the ShenMa High Performance Computing Cluster in Institute of Plasma Physics, Chinese Academy of Sciences. References Magnetohydrodynamic stability of tokamak edge plasmasEdge localized modes and the pedestal: A model based on coupled peeling–ballooning modesNonlinear ELM simulations based on a nonideal peeling–ballooning model using the BOUT++ codeThe physics of edge localized modes (ELMs) and their role in power and particle exhaustBOUT++: A framework for parallel plasma fluid simulationsNonlinear Simulations of Peeling-Ballooning Modes with Anomalous Electron Viscosity and their Role in Edge Localized Mode CrashesELMy H-mode linear simulation with 3-field model on experimental advanced superconducting tokamak using BOUT ++Edge turbulence and divertor heat flux width simulations of Alcator C-Mod discharges using an electromagnetic two-fluid modelAnalysis of a multi-machine database on divertor heat fluxesPlasma turbulence in the scrape-off layer of tokamak devicesSimulations of particle and heat fluxes in an ELMy H-mode discharge on EAST using BOUT++ codeDivertor heat flux simulations in ELMy H-mode discharges of EASTNonlinear fluid simulation of particle and heat fluxes during burst of ELMs on DIII-D with BOUT++ codeNon-linear MHD simulation of ELM energy depositionMHD Equilibrium Reconstruction in the DIII-D TokamakCalculation of the radial electric field with RF sheath boundary conditions in divertor geometryInfluence of equilibrium shear flow on peeling-ballooning instability and edge localized mode crashSix-field two-fluid simulations of peeling–ballooning modes using BOUT++
[1] Connor J W et al 1998 Phys. Plasmas 5 2687
[2] Snyder P B et al 2002 Phys. Plasmas 9 2037
[3] Xu X Q et al 2011 Nucl. Fusion 51 103040
[4] Zohm H 1996 Plasma Phys. Control. Fusion 38 1213
[5] Dudson B D et al 2009 Comput. Phys. Commun. 180 1467
[6] Xu X Q et al 2010 Phys. Rev. Lett. 105 175005
[7] Liu Z X et al 2012 Phys. Plasmas 19 102502
[8] Chen B et al 2017 Nucl. Fusion 57 116025
[9] Makowski M A, Elder D, Gray T K, LaBombard B, Lasnier C J, Leonard A W, Maingi R, Osborne T H, Stangeby P C, Terry J L and Watkins J 2012 Phys. Plasmas 19 056122
[10] Ricci P and Rogers B 2013 Phys. Plasmas 20 010702
[11] Wu Y B, Xia T Y, Zhong F C, Zheng Z, Liu J B and Team EAST 2018 Plasma Phys. Control. Fusion 60 055007
[12] Xia T Y, Xu X Q, Wu Y B, Huang Y Q, Wang L, Zheng Z, Liu J B, Zang Q, Li Y Y, Zhao D and Team EAST 2017 Nucl. Fusion 57 116016
[13]Xu X Q, Umansky M V, Dudson B D and Snyder P B 2008 Commun. Comput. Phys. 4 949
[14] Xia T Y and Xu X Q 2015 Nucl. Fusion 55 113030
[15] Huijsmans G T A and Loarte A 2013 Nucl. Fusion 53 123023
[16] Lao L L, St H E, Peng Q, Ferron J R, Strait E J, Taylor T S, Meyer W H, Zhang C and You K I 2005 Fusion Sci. Technol. 48 968
[17] Gui B, Xia T Y, Xu X Q, Myra J R and Xiao X T 2018 Nucl. Fusion 58 026027
[18] Xi P W, Xu X Q, Wang X G and Xia T Y 2012 Phys. Plasmas 19 092503
[19] Xia T Y, Xu X Q and Xi P W 2013 Nucl. Fusion 53 073009