Benchmark studies of BOUT++ code and TPSMBI code on neutral transport during SMBI

JID:PLA

AID:24423 /SCO Doctopic: Plasma and fluid physics

[m5G; v1.211; Prn:28/03/2017; 16:30] P.1 (1-12)

Physics Letters A ••• (••••) •••–•••

1

67

Contents lists available at ScienceDirect

2

68

3

69

Physics Letters A

4

70

5

71

6

72

www.elsevier.com/locate/pla

7

73

8

74

9

75

10

76

11 12 13

77

Benchmark studies of BOUT++ code and TPSMBI code on neutral transport during SMBI

78 79

14 15

80

17

Y.H. Wang ( A.P. Sun (

18

a

16

19 20 21

a,b,c

d

a, c

) , Z.H. Wang ( ) , W. Guo ( ) , Q.L. Ren ( ) , M. Xu ( ) d , A.K. Wang ( ) d , N. Xiang ( ) a, c

a

) ,

d

26 27 28 29 30 31 32 33 34 35 36

83 85 86 87 88

23 25

82 84

Institute of Plasma Physics, Chinese Academy of Sciences, Hefei 230031, China b University of Science and Technology of China, Hefei 230026, China c Center for Magnetic Fusion Theory, Chinese Academy of Sciences, Hefei 230031, China d Southwestern Institute of Physics, Chengdu 610041, China

22 24

81

89

a r t i c l e

i n f o

Article history: Received 15 November 2016 Received in revised form 22 February 2017 Accepted 19 March 2017 Available online xxxx Communicated by F. Porcelli Keywords: Tokamak plasma Transport Fueling SMBI simulation Benchmark

37 38

a b s t r a c t

90 91

SMBI (supersonic molecule beam injection) plays an important role in tokamak plasma fuelling, density control and ELM mitigation in magnetic confinement plasma physics, which has been widely used in many tokamaks. The trans_neut module of BOUT++ code is the only large-scale parallel 3D fluid code used to simulate the SMBI fueling process, while the TPSMBI (transport of supersonic molecule beam injection) code is a recent developed 1D fluid code of SMBI. In order to find a method to increase SMBI fueling efficiency in H-mode plasma, especially for ITER, it is significant to first verify the codes. The benchmark study between the trans_neut module of BOUT++ code and the TPSMBI code on radial transport dynamics of neutral during SMBI has been first successfully achieved in both slab and cylindrical coordinates. The simulation results from the trans_neut module of BOUT++ code and TPSMBI code are consistent very well with each other. Different upwind schemes have been compared to deal with the sharp gradient front region during the inward propagation of SMBI for the code stability. The influence of the WENO3 (weighted essentially non-oscillatory) and the third order upwind schemes on the benchmark results has also been discussed. © 2017 Elsevier B.V. All rights reserved.

92 93 94 95 96 97 98 99 100 101 102 103 104

39

105

40

106

41

107

42

1. Introduction

43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

It is crucial to study effective fuelling of core plasma for the high-performance of the next generation device, such as ITER [1] and DEMO [2]. SMBI is an important plasma fuelling method besides of gas puffing (GP) and pellet injection (PI). SMBI, a pulsed high speed molecular beam formed by a special Laval type nozzle [3], has a higher efficiency than GP and costs much less than PI. SMBI system was firstly applied in HL-1 tokamak device. With unique advantages, it was then installed in HL-1M [4], HL-2A [5, 6] and other fusion devices, such as HT-7 [7], ASDEX Upgrade [8], JT-60U [9], Tore Supra [10], NSTX [11], KSTAR [6], EAST [12], Heliotron J [13] and so on. SMBI could also be further improved and applied as an optional core fueling method for future fusion devices especially for ITER tokamak besides PI in H-mode plasma discharge. SMBI can also be an important tool to control plasma density [9], mitigate ELMs [14], study the transport of impurity [15] and so on.

61 62 63 64 65 66

E-mail addresses: [email protected] (Z.H. Wang), [email protected] (W. Guo). http://dx.doi.org/10.1016/j.physleta.2017.03.036 0375-9601/© 2017 Elsevier B.V. All rights reserved.

Many experiments have been carried out to study edge fluctuation and particle transport [13], ELM mitigation [14], interaction of plasma transport and turbulence [16] by using SMBI. Some physical events observed during SMBI in experiment could also be further used to understand the high fuelling efficiency of SMBI, such as high enhancement of normalized electron temperature gradient [17], a period of degraded confinement following the injection [9]. However, the background physical mechanism of SMBI fueling is rarely studied in both theory and simulation, which is still not well understood. There is only few related theory and simulation [18,19] to study it. Thus, it is urgent and necessary to develop more realistic physical models and codes of SMBI. In simulation, the large-scale 3D parallel simulation code, the trans_neut module of BOUT++ code using fluid simulation method, has been developed. At present, it is the only large-scale 3D code of SMBI in realistic tokamak geometry. It will be the major code used to simulate the SMBI fuelling process and related physics in the following years, especially for self-consistent simulations of ELMs mitigation due to SMBI which is crucial, urgent and challengeable question. It has simulated and found an effective method to increase SMBI penetration depth and fueling efficiency in L-mode plasma dis-

108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:PLA

AID:24423 /SCO Doctopic: Plasma and fluid physics

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39

[m5G; v1.211; Prn:28/03/2017; 16:30] P.2 (1-12)

Y.H. Wang et al. / Physics Letters A ••• (••••) •••–•••

2

charge [20] but not in H-mode yet. Plasma mean profile variation has also been compared between the simulation and the experiment during SMBI in HL-2A [21]. The results are qualitatively consistent with each other except the electron temperature profile in simulation decreases less than that in experiment may due to the radiation energy loss. However, the trans_neut module of BOUT++ has never been benchmarked with other codes (including fluid codes and particle-in-cell codes) on both GP and SMBI transport, which has already limited its credibility and further usage in experiment interpretation and prediction. Thus, it is urgent and important to do benchmark studies of the trans_neut module with some other codes on GP or SMBI transport. Unfortunately, there are few physical models and codes on SMBI except a newly developed TPSMBI code. At present, TPSMBI, a 1D fluid code, has been developed completely in Fortran computation language and also been verified correctly with analytical solutions [22]. TPSMBI is aimed to simulate the transport dynamics during deuterium and impurity SMBI process, such as non-local heat transport, L–H transition, turbulence transport, impurity transport and so on. The advantage of 1D TPSMBI code is as following: a) TPSMBI, aiming at self-consistent multi-scale simulations of SMBI induced ELM mitigation, is easier and more efficient to be revised and improved than the large scale BOUT++ code, especially for ELM mitigation due to SMBI observed in many tokamaks (i.e., HL-2A, EAST, K-STAR.); b) TPSMBI code can be easily coupled with other transport code (i.e., ONETWO, TRANSP) including more transport physics to do integrated simulation during SMBI; c) TPSMBI costs much less computer resource than BOUT++ under the same simulation conditions. Both BOUT++ code and TPSMBI code have their advantages and disadvantages. It will greatly stimulate and improve our understanding on the underlying physics of SMBI with well benchmarked BOUT++ code and TPSMBI code. This paper reports benchmark studies of BOUT++ code and TPSMBI code during fuelling of SMBI in both slab and cylindrical coordinates. The physical model is described in section 2. The numerical methods are shown in section 3. The numerical conditions are illustrated in section 4. The benchmark results are given in section 5. The summary and conclusion are in section 6.

40 41

2. Physical model

44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

A six-field fluid model of SMBI has been obtained according to the original SMBI model [18] and used to do benchmark between TPSMBI and BOUT++ codes on neutral transport in both slab and cylindrical coordinates. Plasma density N e , electron temperature T e , ion temperature T i , atom density N a , molecule density N m , and molecule radial velocity V xm are involved in the model. The particles are molecules, atoms, ions and electrons in this model. It has considered the major particle collision reactions during SMBI including molecule dissociation, atom ionization, ion–atom charge exchange, and electron–ion recombination. Plasma quasi-neutral condition (i.e., N i = N e ) is applied. The neutral transport equations of BOUT++ code in 1D cylindrical coordinate are the same as those of TPSMBI code, as following:

  ∂ 2 Ne ∂ Ne 1 ∂ Ne P = S IP − S rec − D c⊥e + , ∂t r ∂r ∂ r2   ∂2Te ∂ Te 2 1 ∂ Te − χ⊥c e + ∂t 3 r ∂r ∂ r2   2 2 = νrec W rec − v I T e + W I − v diss W diss 3

−

2me T e − T i Mi

τe

c

c

,

(1)

3

(2)

c

67 68 69

(3)



(4)

72 73 74

(5)

75 76

(6)

c



70 71

where D ⊥e , χ⊥e , χ⊥i , and D ⊥a are the perpendicular diffusion coefficients of electron density, electron temperature, ion temperature and atom density, respectively. M i , me M a and M m are ion, electron, atom and molecule mass, respectively. D c⊥a is calculated from atom perpendicular force balance as

77 78 79 80 81 82 83 84

(7)

85

where T a is atom temperature (the atom and ion temperature are assumed equal due to the high charge exchange rate). S IP is the plasma density source due to atom ionization (S P P representing a particle source). S rec is the plasma density sink due P to the electron–ion recombination. S diss is the atom density source due to the molecule dissociation. ν I , νrec , and νdiss are atom ionization rate, plasma recombination rate, and molecule dissociation rate, respectively. W rec = 4.5 eV is the fraction of recombination energy re-absorbed by electron during recombination via some processes (i.e., three body recombination), while W I = 20 eV and W diss = 4.5 eV are electron energy lost per ionization and dissociation, and W bind = 0.5 eV is the binding energy between the two hydrogen atoms in a molecule. Definitions of quantities associated with the sources and rates due to particle collision reactions are as follows:

87

c

D ⊥a = T a /

M a Ca X

ν

,

86 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102

νCa X = N i σC X V th,i ,

(8)

P = N e rec , rec = N i  rec V th,e , S rec P S I = N e I , I = N a  I V th,e , P = N e diss , diss = Nm  diss V th,e , S diss

ν

ν

ν ν ν

42 43

  ∂2Ti ∂ Ti 2 1 ∂ Ti − χ⊥c i + 2 ∂t 3 r ∂r ∂r 2me T e − T i = (νrec − v I ) T i + , Mi τe   ∂ 2 Na ∂ Na 1 ∂ Na P P = − S IP + 2S diss − D c⊥a + + S rec , ∂t r ∂r ∂ r2 ∂ Nm 1 ∂ P + (r V xm Nm ) = − S diss , ∂t r ∂r ∂ Pm ∂ V xm ∂ V xm + V xm = − ∂r , ∂t ∂r Nm Mm

σ

(9)

σ

ν



σ

(10) (11)



where V th,i = kT i / M i and V th,e = kT e / M e are the ion and electron thermal velocity, respectively. Rate coefficients above are calculated by empirical formulas as following:

103 104 105 106 107 108 109 110 111 112

σC X V th,i  = 1.7 × 10−8 + 1.9 × 10−8

113

(1.5T i eV) − (15 eV) cm3 s−1 , (12) (150T i eV)1/3 − (15 eV)1/3   σ I V th,e  = 3 × 10−8 (0.1T e eV)2 / 3 + ( T e eV)2 cm3 s−1 , (13)   σdiss V th,e  = 3 × 10−8 (0.1T e eV)2 / 3 + ( T e eV)2 cm3 s−1 . (14)

114

σrec V th,e  is calculated by using the formulas in reference [23].

120

τe is the electron–ion collision time.

121

In the model, all source and sink terms of particle and heat are included self-consistently. It has also captured the major transport physics during SMBI including: i) neutral molecules and atoms inward propagation due to molecule convection and atom diffusion effect with a constant injection flux of SMBI; ii) formation of both locally peaked plasma density and decreased plasma temperature profiles due to both particle fueling and heat sinking effects of SMBI via neutrals dissociation and ionization, respectively; iii) self-consistent feedbacks of the locally peaked plasma density on following molecule and atom by increasing their losing rates; iv) saturation of SMBI propagating front due to the total

122

×

1/3

1/3

115 116 117 118 119

123 124 125 126 127 128 129 130 131 132

JID:PLA

AID:24423 /SCO Doctopic: Plasma and fluid physics

[m5G; v1.211; Prn:28/03/2017; 16:30] P.3 (1-12)

Y.H. Wang et al. / Physics Letters A ••• (••••) •••–•••

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19

molecule dissociation rate balancing with the molecule injection rate at the SMBI source; v) propagation of plasma density blobs (i.e., source) and ion temperature holes (i.e., sink). Other interesting physics of drift velocities, magnetic shear, turbulent particle and heat transport, particle and heat pinches, etc. will be included later and benchmarked again. The normalization of BOUT++ and TPSMBI is the same, and the characteristic parameters are as following: L 0 = 2.05 m (a reference length of the HL-2A tokamak), T 0 = 10 eV (a reference plasma edge temperature), M 0 = 2M p (the mass of two protons),



m−3 [18]

V 0 = kT 0 / M 0 (ion thermal velocity at T 0 ), N 0 = 10 (a reference plasma core density). In simulations, diffusion coc efficients of plasma density and temperature are D c⊥e = χ⊥ e = c 2 χ⊥i = 1 m /s. Atom diffusion coefficient is flux limited to pre13

vent unphysical large values and recalculated as following: D ⊥a =  fl fl fl D c⊥a D ⊥a /( D c⊥a + D ⊥a ), where D ⊥a = V th,a L a , V th,a = kT a / M a , and L a is the atom density gradient length at the steepest gradient region. cl

3

stability to solve the convection equation, the third order upwind scheme is used for the convection term v ∂∂ux . The third order upwind scheme is less diffusive and dispersive than the first and second order upwind schemes especially when a large gradient region appears. It has high enough accuracy to study the transport process of SMBI fuelling. It is able to discretize the time space into uniform intervals of size Dt, denoting t i = i Dt, i = 1, 2, . . . , I . Numerical stability is enhanced by choosing the expansion point so as to realize an implicit differencing scheme. The location of the expansion point is determined by the weights 1 − θ and θ assigned to time points t i and t i +1 [27]. The method is unconditionally stable when θ ≥ 12 . Once θ = 12 , a Taylor series analysis of the scheme (Crank–Nicolson scheme) shows that it is second-order accurate in both time and spatial space. The discrete equations are solved by Gauss elimination with the partial pivoting method. More details of the numerical methods used in TPSMBI can be found in the reference [22].

3.1. Numerical methods of BOUT++

24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53

56 57 58 59 60 61 62 63 64 65 66

70 71 72 73 74 75 76 77 78 79 80 81 82 83 85

3.3. Comparison of the third-order upwind and WENO3 schemes

BOUT++ is a finite-differencing code using Method of Lines (MOL). The numerical methods could be the important part of code. It has provided various numerical methods in source codes for physicists to easily deal with numerical problems with higher code accuracy and stability in both time and spatial solution. In spatial derivatives, there are four kinds of numerical methods such as: the first derivative DDX, the second derivative D2DX2, the upwind derivative VDDX (i.e., for terms like v · ∇ ), and the flux method FDDX (i.e., flux conserving and limiting methods for terms like ∇ · (vf)). The second order and forth order central difference scheme can be used for the first derivative DDX and the second derivative D2DX2. The first, second, third CWENO schemes and upwind schemes in different order are supported in BOUT++ code. In this simulation, the second order central differencing scheme is applied for the first derivative DDX and the second derivative D2DX2. The WENO3 [24] scheme and the third order upwind scheme [25], providing higher accuracy and better stability at steep gradients, are used for the upwind derivative VDDX. In time integration, the general implicit time-integration solver CVODE is used in BOUT++ code [26]. It solves a system of the form, df = g (f, t ), where g (f, t ) is in general a non-linear function not dt involving temporal derivatives of f. To advance the system state in time space, the Newton–Krylov BDF implicit method is used in CVODE. More details of both spatial and temporal methods can be found in developer manual of BOUT++ and in the reference [26]. One of the important advantages of BOUT++ code is that it is an object-oriented, highly integrated and high efficient code. As far as possible, this allows the user to concentrate on the physics, rather than worrying about the numerical methods and problems.

It is crucial to treat the upwind derivative VDDX by using higher order upwind schemes for the code better stability and accuracy, such as the third-order upwind and WENO3 schemes mainly applied in the BOUT++ and TPSMBI codes. It usually leads to some difference of simulation results by using different upwind schemes. To do the benchmark well, it requires explaining previously differences between the third-order upwind scheme and the WENO3 scheme. By using a typical transport equation with convection as following:

(16)

=

u ij



+ −

− +

− Dt v u x + v u x ,

3.2. Numerical methods of TPSMBI

6Dx

u+ x =

2

∂u ∂u ∂ u +v − D 2 = S, ∂t ∂x ∂ x

(15)

where u is an evolving physical quantity, v is a convection velocity, D is a diffusion coefficient and S is a source or sink term. The second order central differencing scheme is applied for the second 2

derivative term D ∂∂ xu2 which is the same as that in BOUT++. By considering that the upwind method has higher accuracy and better

(18)

,

ux =

+ 6u ij +1

− 2u ij −1

u ij +1 − u ij −1 2Dx

− u ij −1

2Dx

93 94 95 97 99 100 101 103 104 105 106 108 109

112

(19)

.

−w

−w

113 114

u ij +2 − 3u ij +1 + 3u ij − u ij −1 2Dx w s +(u ij +2 −2.0u ij +1 +u ij )2 w s +(u ij +1 −2.0u ij +u ij −1 )2

u ij +1

− 3u ij

where w = 1.0/(1.0 + r 2 ), and r =

+ 3u ij −1 2Dx

115 116 117

(20)

,

and w s = 10−8 .

118 119 120 121 122

u+ x in equation (17) is defined differently as

u+ x =

92

111

− 3u ij

6Dx

u ij +1

91

110

For the WENO3 scheme, u − x in equation (17) is defined differently as −

90

107

and u + x is defined as:

where w = 1.0/(1.0 + r 2 ), r =

The equations solved by TPSMBI code generally have the simplified form as following:

(17)

2u ij +1 + 3u ij − 6u ij −1 + u ij −2

−u ij +2

89

102



where v + = max( v , 0) and v − = min( v , 0). For the third-order upwind scheme, u − x in equation (17) is defined as:

u− x =

88

98

it is able to discretize the spatial space into uniform intervals of size Dx, denoting x j = j Dx, j = 1, 2, . . . , J . In general, equation (16) can be written with any type upwind schemes as:

u ij+1

87

96

∂u ∂u +v =0 ∂t ∂x

54 55

69

86

3. Numerical methods

22 23

68

84

20 21

67

123

− u ij −2

124

(21)

,

w s +(u ij −2 −2.0u ij −1 +u ij )2 w s +(u ij +1 −2.0u ij +u ij −1 )2

125 126

.

The difference between the third order upwind scheme and WENO3 scheme is due to the different treatment of the velocities − u+ x and u x . It may cause some difference especially in the largest gradient region.

127 128 129 130 131 132

JID:PLA

AID:24423 /SCO Doctopic: Plasma and fluid physics

[m5G; v1.211; Prn:28/03/2017; 16:30] P.4 (1-12)

Y.H. Wang et al. / Physics Letters A ••• (••••) •••–•••

4

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24 25 26

90

Fig. 1. The evolution of radial profiles of molecule velocity in 1D slab coordinate: TPSMBI results (solid curves) and BOUT++ results (star points). The third order upwind scheme and the WENO3 scheme are used in both TPMSBI and BOUT++ codes, respectively.

27 28

4.1. Flux-driven boundary conditions

31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61

At the core boundary, the evolving quantities N e , T i and T e are given fixed gradient or flux-driven boundary conditions as following:

∂ N e ∂x

= −C N e

core

N0 L0

∂ T e,i T0 = − C T e ,i , ∂ x core L0

(22)

66

94 95 96 97 98 99

5. Benchmark results

100

5.1. Benchmark results in slab coordinate

102 103

(23)

where the coefficients C N e = 10 and C T e,i = 800 are input parameters and can be calculated from C N e = Γe |core L 0 / D c⊥e N 0 and

c C T e,i = Q e,i |core L 0 / 23 χ⊥ N | T [18], where Γe |core and Q e,i |core e ,i e core 0 are inflowing particle and heat fluxes at the core boundary. The inflowing particle flux, Γe |core = − D c⊥e ∂∂Nre |core , is about 5N 0 m s−1 .

∂ T e ,i c The inflowing heat flux, Q e,i |core = − 23 χ⊥ N| | , is about e ,i i core ∂ r core − 2 4.3 kW m . These flux-driven boundary conditions at the core

represent continuously inflowing particle and heat fluxes, or ad hoc effective local sources of particle and heat that are necessary to maintain the plasma profiles at steady state. Dirichlet boundary conditions are applied for electron density, electron and ion temperatures at the edge plasma region as follows:

N e |edge = 0.1N 0 ,

(24)

T i ,e |edge = 0.5T 0 .

(25)

At the core, other evolving quantities except for electron density, electron and ion temperatures are given Neumann radial boundary conditions. At the edge, Neumann boundary condition is applied for atom density. 4.2. SMBI radial boundary conditions

64 65

injected in both toroidal and poloidal directions for the benchmark study on the radial transport dynamics. In simulation, N m0 = N 0 = 1 × 1019 m−3 and V xm0 = −1000 m/s are applied by the relevant calculation according the experiment parameters of SMBI in HL-2A [18].

101

,

62 63

92 93

4. Numerical conditions

29 30

91

In simulations, molecule density and molecule radial velocity are given as constants at the edge. SMBI is assumed symmetrically

The benchmark between TPSMBI code and the trans_neut module of BOUT++ code is first carried out in the slab coordinate, which is convenient to identify differences between them. At the beginning, the simulation results between the trans_neut module of BOUT++ code and TPSMBI code were not nearly as encouraging, such as the evolution of radial profiles of molecule velocity and electron density as shown in Figs. 1 and 2. There are some mismatches especially in the largest region where the SMBI front locates. Finally, it is found this is due to different upwind schemes applied in the molecule density and velocity equations between the two codes. The third order upwind scheme is used in the TPSMBI code while the WENO3 scheme is used in the BOUT++ code. As introduced in the previous subsection 3.3 of Section 3, the two schemes have some different treatment of the convection term. The gradient of molecule velocity at about r = 0.85a in TPSMBI is larger than that in BOUT++ as shown in Fig. 1, since the third order upwind scheme is less diffusive than the WENO3 method in the sharpest gradient region [25,26]. The mismatch of the molecule velocity differences then have some influence on other coupling and evolving physical quantities, such as the evolution of radial profiles of electron density as shown in Fig. 2. The radial inward transport dynamic of SMBI front (i.e., the largest gradient in Fig. 1) and the location of the plasma density peaks (Fig. 2) are consistent with each other though there is some mismatch. It can be solved by applying the same upwind scheme as shown later. The benchmark is then carried out again by keeping all the same except applying the same third order upwind scheme in both codes. Now it finds the TPSMBI code benchmarks very well

104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:PLA

AID:24423 /SCO Doctopic: Plasma and fluid physics

[m5G; v1.211; Prn:28/03/2017; 16:30] P.5 (1-12)

Y.H. Wang et al. / Physics Letters A ••• (••••) •••–•••

5

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24 25 26

90

Fig. 2. The evolution of radial profiles of electron density in 1D slab coordinate: TPSMBI results (solid curves) and BOUT++ results (star points). The third order upwind scheme and the WENO3 scheme are used in both TPMSBI and BOUT++ codes, respectively.

91 92

27

93

28

94

29

95

30

96

31

97

32

98

33

99

34

100

35

101

36

102

37

103

38

104

39

105

40

106

41

107

42

108

43

109

44

110

45

111

46

112

47

113

48

114

49

115

50

116

51

117

52

118

53 54

Fig. 3. The evolution of radial profiles of molecule velocity in 1D cylindrical coordinate: TPSMBI results (solid curves) and BOUT++ results (star points). The third order upwind scheme is applied in both BOUT++ and TPSMBI codes. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

55 56 57 58 59 60 61 62 63 64 65 66

119 120 121

with the trans_neut module of BOUT++ code for all of the evolving physical quantities. For brevity, only the consistent simulation results of molecule velocity and electron density are shown in Figs. 3 and 4. Initially, there is no molecule velocity (the red curve in Fig. 3) in the whole resolution domain except a fixed velocity at the edge. The two codes completely written in two different computer languages, BOUT++ code in C/C++ and TPSMBI in Fortran, have benchmarked well with each other in 1D slab coordinate with the same spatial derivative schemes (the second order central derivative and the third order upwind) but different im-

plicit temporal methods. The molecule radial velocity determined by the pressure gradient driven term is 1000 m/s at edge (i.e., the injection velocity) and then increase rapidly at the SMBI propagating front where the molecule density decreases rapidly due to the dissociation. Once the molecules are all dissociated, the molecule density and its pressure decreases to be zero in the region r < 0.8a. Thus, the molecule radial velocity loses its drive and keeps its previous value right before the molecules are all dissociated. The inflowing molecule flux (i.e., N m V xm ) is zero in the region r < 0.8a where the molecule density is zero, even though

122 123 124 125 126 127 128 129 130 131 132

JID:PLA

6

AID:24423 /SCO Doctopic: Plasma and fluid physics

[m5G; v1.211; Prn:28/03/2017; 16:30] P.6 (1-12)

Y.H. Wang et al. / Physics Letters A ••• (••••) •••–•••

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24

90

25 26 27

91

Fig. 4. The evolution of radial profiles of electron density in 1D cylindrical coordinate: TPSMBI results (solid curves) and BOUT++ results (star points). The third order upwind scheme is applied in both BOUT++ and TPSMBI codes.

92 93

28

94

29

95

30

96

31

97

32

98

33

99

34

100

35

101

36

102

37

103

38

104

39

105

40

106

41

107

42

108

43

109

44

110

45

111

46

112

47

113

48

114

49

115

50

116

51

117

52 53

Fig. 5. The evolution of radial profiles of molecule velocity for a higher spatial resolution in 1D slab coordinate: TPSMBI results (solid curves) and BOUT++ results (star points). The third order upwind scheme and the WENO3 scheme are used in both TPMSBI and BOUT++ code, respectively.

54

57 58 59 60 61 62 63 64 65 66

119 120

55 56

118

121

the radial molecule velocity is very high. Since the damping term due to molecule plasma collision has not been considered in the momentum transport equation, the decrease of the molecule velocity due the damping has not been observed here but it will be improved and studied later. In general, the simulation results will also been influenced by the time and spatial resolution of the research domain except the different numerical schemes. Thus, we have also tried to reduce the time interval step Dt in TPSMBI code by increasing its time resolution, while Dt is automatically small until it meets the requirement of absolute and relative solver tolerances defined in

BOUT++ code. It cannot decrease the mismatch. The benchmark is then carried out again by increasing the spatial resolution with a much smaller spatial step Dx in 1D slab coordinate. The third order upwind scheme and WENO3 scheme are used in TPSMBI and BOUT++ codes, respectively. The mismatch of all the variables decreases much once the total number of the spatial discrete points increases from 68 to 260 in both TPSMBI and BOUT++ codes. For brevity, only the simulation results of molecule velocity and electron density are shown in Figs. 5 and 6. The results are consistent well with each other for a higher spatial resolution but different upwind schemes.

122 123 124 125 126 127 128 129 130 131 132

JID:PLA

AID:24423 /SCO Doctopic: Plasma and fluid physics

[m5G; v1.211; Prn:28/03/2017; 16:30] P.7 (1-12)

Y.H. Wang et al. / Physics Letters A ••• (••••) •••–•••

7

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24

90

25 26

Fig. 6. The evolution of radial profiles of electron density for a higher spatial resolution in 1D slab coordinate: TPSMBI results (solid curves) and BOUT++ results (star points). The third order upwind scheme and the WENO3 scheme are used in both TPSMBI and BOUT++ codes, respectively.

91 92

27

93

28

94

29

95

30

96

31

97

32

98

33

99

34

100

35

101

36

102

37

103

38

104

39

105

40

106

41

107

42

108

43

109

44

110

45

111

46

112

47

113

48

114

49

115

50

116

51 52 53

117

Fig. 7. The evolution of radial profiles of atom density in 1D cylindrical coordinate: TPSMBI results (solid curves) and BOUT++ results (star points). The same difference schemes are used in both BOUT++ and TPSMBI codes.

54 55

58 59 60 61 62 63 64 65 66

119 120

5.2. Benchmark results in cylindrical coordinate

56 57

118

The benchmark of the two codes is also carried out in the 1D cylindrical coordinate to increase the reliability. As the benchmark in the slab coordinate, the third order upwind scheme is first applied in TPSMBI code while WENO3 scheme is applied in BOUT++ code in the molecule density and velocity equations. A similar level of discrepancy is found between BOUT++ and TPSMBI results as that in the slab coordinate if the applied upwind schemes are not the same. For brevity, it is not shown again here. The difference of simulation results between the trans_neut module of BOUT++ code and TPSMBI code can be reduced by

using the same spatial derivative schemes. Consistent simulation results of the trans_neut module of BOUT++ code and TPSMBI code have been achieved by applying the same spatial derivative schemes in both codes. The evolution of all the physical quantities is well consistent with each other between the two codes. For brevity, only the simulation results of atom density and ion temperature are shown in Figs. 7 and 8, respectively. The atom density front propagates deeper into r ∼ 0.8a, becomes steeper with the increase of the atom density and then it moves a little bit outwards due to high ionization but having a higher atom density at the edge (Fig. 7). The initial ion temperature is decreased locally on the SMBI fueling path, espe-

121 122 123 124 125 126 127 128 129 130 131 132

JID:PLA

8

AID:24423 /SCO Doctopic: Plasma and fluid physics

[m5G; v1.211; Prn:28/03/2017; 16:30] P.8 (1-12)

Y.H. Wang et al. / Physics Letters A ••• (••••) •••–•••

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24 25 26

90

Fig. 8. The evolution of radial profiles of ion temperature in 1D cylindrical coordinate: TPSMBI results (solid curves) and BOUT++ results (star points). The same difference schemes are used in both BOUT++ and TPSMBI codes.

91 92

27

93

28

94

29

95

30

96

31

97

32

98

33

99

34

100

35

101

36

102

37

103

38

104

39

105

40

106

41

107

42

108

43

109

44

110

45

111

46

112

47

113

48

114

49

115

50

116

51

117

52 53 54

118

Fig. 9. The evolution of radial profiles of molecule density for a lower molecule injection density in 1D cylindrical coordinate: TPSMBI results (solid curves) and BOUT++ results (star points). The same difference schemes are used in both BOUT++ and TPSMBI codes.

55 56 57 58 59 60 61 62 63 64 65 66

119 120 121

cially in the edge region where SMBI penetrates (Fig. 8). Thus, it has proved that the trans_neut module of BOUT++ code and the TPSMBI code are also well benchmarked with each other in the 1D cylindrical coordinate by applying the same spatial derivative schemes but different implicit temporal methods. In order to increase more the credibility, it is still necessary to benchmark the codes in a much broader range of plasma parameters and boundary conditions. It will further do more benchmarks in the 1D cylindrical coordinate as follows in Subsection 5.3, 5.4 and 5.5.

5.3. Benchmark with a lower molecule injection density

122 123 124

The benchmark with an edge boundary condition of a lower molecule density, N m0 = 0.02N 0 , is also made in 1D cylindrical coordinate without changing any other parameters, boundary and initial conditions of the Subsection 5.2 of the Section 5. All the mean profiles obtained in BOUT++ and TPSMBI codes are consistent with each other very well when the same difference schemes are used in both BOUT++ and TPSMBI codes. For brevity, only the simulation results of molecule density and electron density are shown

125 126 127 128 129 130 131 132

JID:PLA

AID:24423 /SCO Doctopic: Plasma and fluid physics

[m5G; v1.211; Prn:28/03/2017; 16:30] P.9 (1-12)

Y.H. Wang et al. / Physics Letters A ••• (••••) •••–•••

9

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24 25 26

90

Fig. 10. The evolution of radial profiles of electron density for a lower molecule injection density in 1D cylindrical coordinate: TPSMBI results (solid curves) and BOUT++ results (star points). The same difference schemes are used in both BOUT++ and TPSMBI codes.

91 92

27

93

28

94

29

95

30

96

31

97

32

98

33

99

34

100

35

101

36

102

37

103

38

104

39

105

40

106

41

107

42

108

43

109

44

110

45

111

46

112

47

113

48

114

49

115

50

116

51 52 53

117

Fig. 11. The evolution of radial profiles of electron density for the edge Neumann boundary condition of electron density in 1D cylindrical coordinate: TPSMBI results (solid curves) and BOUT++ results (star points). The same difference schemes are used in both BOUT++ and TPSMBI codes.

54 55 56 57

60 61

in Figs. 7 and 8. It indicates that both of the codes are also reliable for the long enough time simulation and even a lower molecule density condition. See also Figs. 9 and 10. 5.4. Benchmark with the edge Neumann boundary condition of electron density

62 63 64 65 66

119 120

58 59

118

The benchmark with the edge Neumann boundary condition for electron density is also made in 1D cylindrical coordinate without changing any other parameters, boundary and initial conditions of N the Subsection 5.2 of the Section 5. ∂∂Nre |edge = −20.0 L 0 is used 0

in simulation. All the mean profiles of evolving physical quantities obtained in BOUT++ and TPSMBI codes agree very well with each other when the same difference scheme are used in both codes. For brevity, only the radial profiles of electron density and temperature are shown in Figs. 11–12. It indicates that the two codes are reliable in both of Neumann and Dirichlet boundary conditions.

121 122 123 124 125 126 127

5.5. Benchmark with a much lower molecule injection density and a lower initial electron temperature profile

128 129 130

The benchmark is done again with a much lower molecule injection density and a lower initial electron temperature profile in

131 132

JID:PLA

10

AID:24423 /SCO Doctopic: Plasma and fluid physics

[m5G; v1.211; Prn:28/03/2017; 16:30] P.10 (1-12)

Y.H. Wang et al. / Physics Letters A ••• (••••) •••–•••

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24 25 26

90

Fig. 12. The evolution of radial profiles of electron temperature for the edge Neumann boundary condition of electron density in 1D cylindrical coordinate: TPSMBI results (solid curves) and BOUT++ results (star points). The same difference schemes are used in both BOUT++ and TPSMBI codes.

91 92

27

93

28

94

29

95

30

96

31

97

32

98

33

99

34

100

35

101

36

102

37

103

38

104

39

105

40

106

41

107

42

108

43

109

44

110

45

111

46

112

47

113

48

114

49

115

50

116

51

117

52

118

53 54

Fig. 13. The evolution of radial profiles of molecule density for a much lower molecule injection density and a lower initial electron temperature profile in 1D cylindrical coordinate: TPSMBI results (solid curves) and BOUT++ results (star points). The same difference schemes are used in both BOUT++ and TPSMBI codes.

55

58 59 60 61 62 63 64 65 66

120 121

56 57

119

122

1D cylindrical coordinate without changing any other parameters, boundary and initial conditions of the Subsection 5.2 of the SecT tion 5. N m0 = 0.005N 0 and the core boundary ∂∂Tre = −600 L 0 are 0 used in the simulation. All the radial profiles of evolving physical quantities obtained in BOUT++ and TPSMBI codes agree very well with each other. The profiles of molecule density and electron temperature are shown in Figs. 13–14. It demonstrates that the numerical methods of TPSMBI and BOUT++ code are still reliable when a much lower molecule injection density and a different initial plasma temperature profile is used.

A lot of consistent benchmarks have been obtained and only a part of them are shown above for simplicity. In conclusion, the benchmarks between the trans_neut module of BOUT++ code and the TPSMBI code have been done successfully in both 1D slab and cylindrical geometries in L-mode discharge under various conditions. It indicates that BOUT++ and TPSMBI codes have the ability to give the same simulation results of the radial SMBI transport process in L-mode discharge without the SOL region. The H-mode dynamics (i.e., pedestal formation, turbulence suppression, ELM induced transport) is much more complicated and has not

123 124 125 126 127 128 129 130 131 132

JID:PLA

AID:24423 /SCO Doctopic: Plasma and fluid physics

[m5G; v1.211; Prn:28/03/2017; 16:30] P.11 (1-12)

Y.H. Wang et al. / Physics Letters A ••• (••••) •••–•••

11

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24 25 26

90

Fig. 14. The evolution of radial profiles of electron temperature for a much lower molecule injection density and a lower initial electron temperature profile in 1D cylindrical coordinate: TPSMBI results (solid curves) and BOUT++ results (star points). The same difference schemes are used in both BOUT++ and TPSMBI codes.

27 28 29 30 31 32 33

the key physics of ELM cycles [16], turbulent transport [28], profile bifurcation [29], ExB velocity shear [30] etc. It is an important verification of the trans_neut module of BOUT++ and TPSMBI code. It will greatly accelerate the related research of SMBI, especially for SMBI induced ELM mitigation and its fueling efficiency increment in H-mode.

6. Summary and conclusions

Acknowledgements

The benchmark of the trans_neut module of BOUT++ code and TPSMBI code on radial transport dynamic during SMBI has been first successfully achieved in both 1D slab and cylindrical coordinates. It has covered a much broader range of simulation conditions relevant and crucial for SMBI modeling. The simulation results, obtained from the trans_neut module of BOUT++ code and the TPSMBI code, are very well consistent with each other by applying the same spatial derivative schemes (the second order central derivative and the third order upwind) but different implicit temporal methods. The physical results revealing by the two codes are also consistent well with each other during SMBI such as the SMBI front propagation, the mean profile evolution and so on. Different upwind schemes have been compared and used in molecule density and molecule velocity equations to deal with the large gradient front region during the inward propagation of SMBI. The third order upwind scheme is found to be less diffusive than the WENO3 scheme in the largest gradient region. The influence of WENO3 and the third order upwind schemes on the benchmark results has also been discussed. More benchmarks with other fluid or particle codes (i.e., UEDGE, DEGAS2, SOLPS) on neutral transport dynamics during GP will be done later. It demonstrates that TPSMBI code is appropriate for the radial transport dynamic studies of the SMBI front penetration and mean profile evolution inside of the separatrix in plasma. It is also appropriate to interpret the radial mean profile variation during SMBI in the experiment by making some reasonable assumptions. TPSMBI is much more efficient in computation than the large scale simulation code. It is still challenging for TPSMBI to benchmark with BOUT++ or validate with the experiment in the H-mode plasma, since it requires the modification of the physical model to include

The authors wish to thank X.Q. Xu, B.D. Dudson and M.V. Umansky for their contributions to BOUT++ framework. This work was supported by the National Natural Science Foundation of China under Grant Nos. 11575055, 11475219, 11375053 and the National Magnetic Confinement Fusion Science Program of China under Grant Nos. 2013GB111005, 2014GB108004, 2015GB110001. The authors would also like to acknowledge the ShenMa High Performance Computing Cluster at the Institute of Plasma Physics, Chinese Academy of Sciences.

38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

94 95 96 97 98 99 100

36 37

92 93

been included in the physical model. Benchmark of BOUT++ and TPSMBI codes in H-mode discharge will be done later with an improved physical mode of SMBI including the H-mode dynamics. The benchmark of neutral transport dynamics during GP in 2D and 3D geometries will be done between BOUT++ code and other codes (such as UEDGE, DEGAS2, SOLPS) later.

34 35

91

101 102 103 104 105 106 107 108 109 110 111 112

References

113 114

[1] L. Garzotti, P. Belo, G. Corrigan, F. Köchl, J. Lönnroth, V. Parail, G. Pereverzev, S. Saarelma, G. Tardini, M. Valoviˇc, I. Voitsekhovitch, S. Wiesen, Simulations of density profiles, pellet fuelling and density control in ITER, Nucl. Fusion 52 (2012) 013002. [2] P. Vincenzi, F. Koechl, L. Garzotti, D.B. King, E. Tindale, T. Bolzonella, P.T. Lang, B. Pegourié, M. Romanelli, R. Wenninger, Fuelling and density control for DEMO, Nucl. Fusion 55 (2015) 113028. [3] L.H. Yao, N.Y. Tang, Z.Y. Cui, D.M. Xu, Z.C. Deng, X.T. Ding, J.l. Luo, J.F. Dong, G.C. Guo, S.K. Yang, C.H. Cui, Z.G. Xiao, D.Q. Liu, X.P. Chen, L.W. Yan, D.H. Yan, E.Y. Wang, X.W. Deng, Plasma behaviour with molecular beam injection in the HL-1M tokamak, Nucl. Fusion 38 (1998) 631. [4] J.F. Dong, N.Y. Tang, W. Li, J.l. Luo, Y. Liu, Z.G. Xiao, F.B. Yao Lianghua, Li Bo, Qin Yunwen, Research on the penetration characteristics of three refuelling methods on HL-1M, Plasma Phys. Control. Fusion 44 (2002) 371–379. [5] L.H. Yao, B.B. Feng, C.Y. Chen, Z.B. Shi, B.S. Yuan, Y. Zhou, X.R. Duan, H.J. Sun, J. Lu, Y.M. Jiao, G.Q. Ni, H.Y. Lu, W.W. Xiao, W. Li, Y.D. Pan, W.Y. Hong, H. Ran, X.T. Ding, Y. Liu, Plasma behaviour with hydrogen supersonic molecular beam and cluster jet injection in the HL-2A tokamak, Nucl. Fusion 47 (2007) 1399–1410. [6] W.W. Xiao, P.H. Diamond, W.C. Kim, L.H. Yao, S.W. Yoon, X.T. Ding, S.H. Hahn, J. Kim, M. Xu, C.Y. Chen, B.B. Feng, J. Cheng, W.L. Zhong, Z.B. Shi, M. Jiang, X.Y. Han, Y.U. Nam, W.H. Ko, S.G. Lee, J.G. Bak, J.W. Ahn, H.K. Kim, H.T. Kim, K.P. Kim, X.L. Zou, S.D. Song, J.I. Song, Y.W. Yu, T. Rhee, J.M. Kwon, X.L. Huang,

115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:PLA

12

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38

AID:24423 /SCO Doctopic: Plasma and fluid physics

[m5G; v1.211; Prn:28/03/2017; 16:30] P.12 (1-12)

Y.H. Wang et al. / Physics Letters A ••• (••••) •••–•••

D.L. Yu, K.D. Lee, S.I. Park, M. Jung, S. Zoletnik, M. Lampert, G.R. Tynan, Y.S. Bae, J.G. Kwak, L.W. Yan, X.R. Duan, Y.K. Oh, J.Q. Dong, ELM mitigation by supersonic molecular beam injection: KSTAR and HL-2A experiments and theory, Nucl. Fusion 54 (2014) 023003. [7] X. Gao, Y.X. Jie, C.Y. Xia, M.S. Wei, Y. Yang, S.Y. Zhang, J.Y. Zhao, L.Q. Hu, Y.B. Zhu, J.R. Luo, Y.P. Zhao, N. Qiu, J.G. Li, B.N. Wan, G.L. Kuang, X.D. Zhang, X.N. Liu, X.Z. Gong, Y. Bao, B.L. Lin, Z.W. Wu, Y.D. Li, Y.J. Shi, M. Song, P. Fu, X.M. Zhang, M. Zeng, A.G. Xie, N.Z. Cui, H.L. Ruan, L. Wang, B. Sheng, S. Liu, W.W. Ye, K. Yang, J.K. Liu, Y.F. Cheng, H.Y. Fan, S.X. Liu, X.D. Tong, J.S. Mao, X.M. Gu, J.K. Xie, Y.X. Wan, High density operation on the HT-7 superconducting tokamak, Nucl. Fusion 40 (2000) 1875–1883. [8] P.T. Lang, J. Neuhauser, J. Bucalossi, A. Chankin, D.P. Coster, R. Drube, R. Dux, G. Haas, L.D. Horton, S. Kalvin, G. Kocsis, M. Maraschek, V. Mertens, V. Rohde, V. Rozhansky, R. Schneider, I. Senichenkov, I. Veselova, E. Wolfrum, A.U. Team, Impact of a pulsed supersonic deuterium gas jet on the ELM behaviour in ASDEX upgrade, Plasma Phys. Control. Fusion 47 (2005) 1495–1516. [9] H. Takenaga, Particle control study towards burning plasma control in JT-60U, J. Nucl. Mater. 390–391 (2009) 869–875. [10] B. Pégourié, E. Tsitrone, R. Dejarnac, J. Bucalossi, G. Martin, J. Gunn, D. Frigione, D. Reiter, P. Ghendrih, C. Clément, Supersonic gas injection on Tore Supra, J. Nucl. Mater. 313–316 (2003) 539–542. [11] V.A. Soukhanovskii, R. Maingi, C.E. Bush, R. Raman, R.E. Bell, R. Kaita, H.W. Kugel, C.J. Lasnier, B.P. LeBlanc, J.E. Menard, S.F. Paul, A.L. Roquemore, Divertor heat flux reduction and detachment experiments in NSTX, J. Nucl. Mater. 363–365 (2007) 432–436. [12] X. Zheng, J. Li, J. Hu, J. Li, R. Ding, B. Cao, J. Wu, Comparison between gas puffing and supersonic molecular beam injection in plasma density feedback experiments in EAST, Plasma Phys. Control. Fusion 55 (2013) 115010. [13] L. Zang, T. Mizuuchi, N. Nishino, S. Yamamoto, K. Kasajima, K. Hashimoto, M. Sha, M. Takeuchi, K. Mukai, H.Y. Lee, N. Kenmochi, K. Nagasaki, S. Kado, H. Okada, T. Minami, S. Kobayashi, N. Shi, S. Konoshima, Y. Nakamura, F. Sano, Effect of supersonic molecular-beam injection on edge fluctuation and particle transport in Heliotron J, Phys. Plasmas 21 (2014) 042308. [14] W.W. Xiao, P.H. Diamond, X.L. Zou, J.Q. Dong, X.T. Ding, L.H. Yao, B.B. Feng, C.Y. Chen, W.L. Zhong, M. Xu, B.S. Yuan, T. Rhee, J.M. Kwon, Z.B. Shi, J. Rao, G.J. Lei, J.Y. Cao, J. Zhou, M. Huang, D.L. Yu, Y. Huang, K.J. Zhao, Z.Y. Cui, X.M. Song, Y.D. Gao, Y.P. Zhang, J. Cheng, X.Y. Han, Y. Zhou, Y.B. Dong, X.Q. Ji, Q.W. Yang, Y. Liu, L.W. Yan, X.R. Duan, Y. Liu, ELM mitigation by supersonic molecular beam injection into the H-mode pedestal in the HL-2A tokamak, Nucl. Fusion 52 (2012) 114027. [15] X.W. Cui, Z.Y. Cui, B.B. Feng, Y.D. Pan, H.Y. Zhou, P. Sun, B.Z. Fu, P. Lu, Y.B. Dong, J.M. Gao, S.D. Song, Q.W. Yang, Investigation of impurity transport using supersonic molecular beam injected neon in HL-2A ECRH plasma, Chin. Phys. B 22 (2013) 125201.

[16] P. Tamain, G. Bonhomme, F. Brochard, F. Clairet, C. Gil, J. Gunn, P. Hennequin, G. Hornung, J.L. Segui, L. Vermare, P. Ghendrih, Interaction of plasma transport and turbulence on particle fuelling, J. Nucl. Mater. 438 (2013) S148–S154. [17] D.L. Yu, C.Y. Chen, L.H. Yao, B.B. Feng, X.Y. Han, L.M. Yang, W.L. Zhong, Y. Zhou, K.J. Zhao, Y. Huang, Y. Liu, L.W. Yan, Q.W. Yang, J.Q. Dong, X.R. Duan, Penetration characteristics of supersonic molecular beam injection on HL-2A tokamak, Nucl. Fusion 50 (2010) 035009. [18] Z.H. Wang, X.Q. Xu, T.Y. Xia, T.D. Rognlien, 2D simulations of transport dynamics during tokamak fuelling by supersonic molecular beam injection, Nucl. Fusion 54 (2014) 043019. [19] Y.L. Zhou, Z.H. Wang, X.Q. Xu, H.D. Li, H. Feng, W.G. Sun, Comparisons between tokamak fueling of gas puffing and supersonic molecular beam injection in 2D simulations, Phys. Plasmas 22 (2015) 012503. [20] Y.L. Zhou, Z.H. Wang, M. Xu, Q. Wang, L. Nie, H. Feng, W.G. Sun, Investigation of molecular penetration depth variation with SMBI fluxes, Chin. Phys. B 25 (2016) 095201. [21] Z.H. Wang, X.Q. Xu, T.Y. Xia, D.L. Yu, G.Y. Zheng, S.Y. Chen, W.L. Zhong, Z.B. Shi, A.P. Sun, J.Q. Dong, M. Xu, T.T. Sun, L.H. Yao, Simulations and validations of transport during fueling by SMBI in HL-2A tokamak, IAEA, DOI (2014) TH/P7-30. [22] Y.H. Wang, W.F. Guo, Z.H. Wang, Q.L. Ren, A.P. Sun, M. Xu, A.K. Wang, N. Xiang, Radial transport dynamics studies of SMBI with a newly developed TPSMBI code, Chin. Phys. B 25 (2016) 106601. [23] M.J. Seaton, Radiative recombination of hydrogenic ions, Mon. Not. R. Astron. Soc. 119 (1959) 81–89. [24] G.S. Jiang, C.W. Shu, Efficient implementation of weighted ENO schemes, J. Comput. Phys. 126 (1996) 202–228. [25] M. Dehghan, On the numerical solution of the one-dimensional convection– diffusion equation, Math. Probl. Eng. 2005 (2005) 61–74. [26] B.D. Dudson, M.V. Umansky, X.Q. Xu, P.B. Snyder, H.R. Wilson, BOUT++: a framework for parallel plasma fluid simulations, Comput. Phys. Commun. 180 (2009) 1467–1480. [27] E.S. Gross, V. Casulli, L. Bonaventura, J.R. Koseff, A semi-implicit method for vertical transport in multidimensional models, Int. J. Numer. Methods Fluids 28 (1998) 157–186.

67

[28] Z.H. Wang, P.H. Diamond, Ö.D. Gürcan, X. Garbet, X.G. Wang, Turbulence propagation in heat flux-driven plasmas: implications for temperature profile structure, Nucl. Fusion 51 (2011) 073009. [29] K. Miki, P.H. Diamond, O.D. Gürcan, G.R. Tynan, T. Estrada, L. Schmitz, G.S. Xu, Spatio-temporal evolution of the L → I → H transition, Phys. Plasmas 19 (2012) 092306. [30] K.H. Burrell, E.J. Doyle, P. Gohil, R.J. Groebner, J. Kim, R.J. La Haye, L.L. Lao, R.A. Moyer, T.H. Osborne, W.A. Peebles, C.L. Rettig, T.H. Rhodes, D.M. Thomas, Role of the radial electric field in the transition from L (low) mode to H (high) mode to VH (very high) mode in the DIII-D tokamak Phys. Plasmas 1 (1994) 1536.

96

68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 97 98 99 100 101 102 103 104

39

105

40

106

41

107

42

108

43

109

44

110

45

111

46

112

47

113

48

114

49

115

50

116

51

117

52

118

53

119

54

120

55

121

56

122

57

123

58

124

59

125

60

126

61

127

62

128

63

129

64

130

65

131

66

132