Plasma and Fusion Research: Regular Articles

Volume 19, 1403016 (2024)

Eﬀect of Edge Shear Flow on Radial Spreading of Ballooning Mode
Turbulence
Ayumi TAKANO, Seiya NISHIMURA∗ ,1) , Makoto SASAKI2) and Yusuke KOSUGA3)
Department of Electric and Electronic Engineering, Hosei University, Koganei 184-8584, Japan
1)
National Institute for Quantum Science and Technology, Naka 311-0193, Japan
2)
College of Industrial Technology, Nihon University, Narashino 275-8575, Japan
3)
Research Institute for Applied Mechanics, Kyushu University, Fukuoka 816-8580, Japan(Received 3 October 2023 / Accepted 8 February 2024)
Turbulent transport by the ballooning mode in tokamak plasmas with edge pedestals is simulated using a
reduced set of two-fluid equations. In the absence of the equilibrium poloidal flow, global heat transport by
the secondary nonlinear evolution of the resistive ballooning mode turbulence is observed. By examining theeﬀect of the edge shear flow, the global heat transport is suppressed to be almost half, if the edge shear flow
is strong enough. A detailed analysis on the radial profile and the poloidal spectrum of the heat flux is newly
performed. It is revealed that such confinement improvement is caused by the suppression of the formation of
streamer structures that lead to the strong convective heat transport.
c 2024 The Japan Society of Plasma Science and Nuclear Fusion Research⃝

Keywords: tokamak plasma, resistive ballooning mode, nonlinear evolution, turbulence spreading, shear flow
DOI: 10.1585/pfr.19.1403016

1. Introduction
For magnetically confinement fusion plasmas such as
tokamak plasmas, it is an important issue to understand andcontrol magnetohydrodynamic (MHD) turbulence. In particular, the edge localized mode (ELM) due to the ballooning mode has been observed in the edge region of tokamaks, causing partial collapse of the pressure profile andlarge heat load on the divertor [1]. In addition, a recent
experimental observation shows that the radially elongatedflow, i.e. the streamer, generated by the ballooning mode
turbulence plays an important role in the triggering of thetype-III ELM [2].
Theoretical and numerical studies on the ballooning
mode are mainly classified into the following four categories. First, the stability analysis of the ballooning mode
has been conducted taking into account realistic poloidal
geometries and various equilibrium profiles [3–5]. In thisline of research, it has been considered that coupling of thepeeling mode, driven by the edge bootstrap current, andthe ballooning mode is essential for the ELM accompanying the largest collapse, so-called type-I ELM [6]. Second,nonlinear simulations of the ballooning mode have been
performed to understand nonlinear evolution and saturation mechanisms of the ballooning mode turbulence [7–9].
Third, in order to reproduce the repetitive burst transport
due to the ELM, nonlinear simulations of the ballooning
mode turbulence in the presence of the heat source havebeen performed [10, 11]. Fourth, associated with the interaction between the ballooning mode turbulence and mag∗ Corresponding author’s e-mail: nishimura.seiya@qst.go.jp

netic islands, studies have been carried out to examine thephysical mechanism of how bursty heat transport by the
ELM is mitigated or suppressed by the so-called resonant
magnetic perturbation (RMP) [12–15].
In this study, we focus on the nonlinear evolution of
the resistive ballooning mode turbulence in the absenceof the heat source and the peeling mode. In the preceding works, it has been shown that the linear instability ofthe ballooning mode occurs locally, whereas the ballooning mode turbulence spreads globally in the radial direction [7–9]. Originally, the turbulence spreading in the radial direction has been studied in detail in terms of the drift
wave turbulence [16–18]. In the turbulence spreading, thepressure gradient and the turbulence amplitude move radially inward are simultaneously observed, which is also
called the front propagation. However, for the ballooning
mode turbulence, the detailed mechanism of the turbulence
spreading is not fully understood. In particular, it has been
pointed out that the loss of the stored energy due to theballooning mode turbulence is suppressed to a low level,
when the edge poloidal flow exists [8, 9]. In order to clarify such a reduction mechanism, we newly introduce ananalyzing method of the radial profile and the mode number spectrum of the heat flux accompanying the ballooning
mode turbulence.
Remainder of this paper is organized as follows. In
Sec. 2, the reduced two-fluid model of tokamak plasmas
is introduced. We show a law of the conservation of energy and derive a time evolution equation of the stored
energy. In addition, parameters and equilibrium profiles
used in nonlinear simulations are also shown. The nonlin-1403016-1

c 2024 The Japan Society of Plasma⃝
Science and Nuclear Fusion Research

Plasma and Fusion Research: Regular Articles

Volume 19, 1403016 (2024)

ear simulation results and analyses on them are shown inSec. 3. The reduction of the stored energy loss due to theedge poloidal shear flow is observed. The heat flux due tothe ballooning mode turbulence is analyzed, and the reduction of the convective heat flux by the shear flow is foundto be essential. Then, time evolution of the radial profileand the mode number spectrum of the convective heat flux
are analyzed. It is concluded that the edge shear flow affects the nonlinear characteristics of the ballooning mode
turbulence, which avoids both the radial turbulence spreading and the global heat transport. Section 4 is devoted fora summary.

2. Simulation Model
2.1 Reduced two-fluid model
Our simulation model is based on the conventional
Hazeltine’s reduced two-fluid model [19], where tokamak
plasmas with large aspect ratios are modeled by cylindricalplasmas with eﬀective toroidal curvatures. In the model,
a quasi-neutral condition between ion and electron densities is considered, and the cold ion approximation is used.
For simplicity, we neglect the ion diamagnetic drift, theion parallel velocity, and the electron inertia. The vorticity
equation, the generalized Ohm’s law, and the pressure evolution equation (the electron continuity equation for constant temperature) are given by[ ]dU= ∇∥ j∥ − x, p + µ∇2⊥ U,
dt∂A∥= −∇∥ (ϕ − δp) − η∥ ( j∥ − j∥0 ),∂tdp= 2δβ∇∥ j∥ + 2β[x, ϕ − δp] + χ∇2⊥ p,
dt

(1)
(2)

2.2

Energy conservation and heat flux

In the cylindrical coordinates, we assume that an arbitrary variable f (r, θ, z, t) can be expanded into complex
Fourier series asf (r, θ, z, t) = f0 (r) + f˜(r, θ, z, t)∑= f0 (r) +f˜m,n (r, t) exp [i(mθ − nz)],
m,n

(10)
(3)

respectively. The definition of U and the Ampère’s law are
given byU = ∇2⊥ ϕ,(4)

j∥ = −∇2⊥ A∥ ,(5)

respectively. The operators are defined byd∂=+ [ϕ, ],dt ∂t∂
− [A∥ , ],∇∥ =∂z[ f, g] = ẑ · ∇⊥ f × ∇⊥ g,∂
1 ∂∇⊥ = r̂ + θ̂
,
∂r
r ∂θ

the electron pressure, respectively, and j∥0 is the equilibrium parallel current. The transport coeﬃcient {µ, η∥ , χ}denote the perpendicular momentum diﬀusivity, the parallel resistivity, and the perpendicular thermal diﬀusivity,
respectively. The major radial position measured from theplasma center is defined by x = r cos θ, and the eﬀective
curvature of the toroidal magnetic field field is given byκ0 = −∇⊥ x.
The time and length are normalized as t/τA → t,
r/a → r, and z/R0 → z, where τA = R0 /vA , vA is the
Alfvén velocity, a is the minor radius, and R0 is the major radius. The normalized parameters are defined as follows: β is the ratio between the kinetic pressure and themagnetic pressure measured at the plasma center, and δ
is the ion inertia length normalized by 2a. The other variables and operators are normalized as follows: a∇⊥ → ∇⊥ ,R0 ∇∥ → ∇∥ , R0 κ0 → κ0 , (cτA /a2 B0 )ϕ → ϕ, A∥ /ϵaB0 →A∥ , (4πR0 /cB0 ) j∥ → j∥ , (8π/ϵB20 )p → p, (τA /a2 )µ → µ,(c2 τA /4πa2 )η∥ → η∥ , and (τA /a2 )χ → χ, where c is the velocity of light, ϵ = a/R0 is the inverse aspect ratio of torus,
and B0 is the representative magnitude of the toroidal magnetic field.

(6)
(7)
(8)
(9)

where f and g are arbitrary variables.
In the cylindrical coordinates (r, θ, z), the coordinate
variables {r, θ, z} are the minor radial position, the poloidal
phase angle, and the toroidal position, respectively, and{r̂, θ̂, ẑ} are unit vectors. The variables {U, ϕ, A∥ , j∥ , p} denote the vorticity, the electrostatic potential, the parallel
component of the vector potential, the parallel current, andwhere {m, n} are the poloidal and toroidal mode numbers,
respectively. Hereafter, variables with the subscript 0 denote equilibrium quantities, and those with the tilde denote
perturbation quantities.
In our preliminary study, it has been observed thatnumerical instability occurs around r = 0 when a finite
value of δ is considered. In order to avoid such instability, we limit the domain of the r as r0 ≤ r ≤ 1,
where r0 = 0.1 is used in this study. The boundary
conditions for the perturbation amplitudes are given suchthat f˜m,n (r0 , t) = f˜m,n (1, t) = 0 for (m, n) , (0, 0) and∂r f˜0,0 |r=r0 = f˜0,0 (1, t) = 0.
Then, we write down Eqs. (1)-(3) for the (m, n) component of the complex Fourier amplitude as[]∂Ũm,n = − ikθ ϕ′0 Ũm,n + ikθ U0′ ϕ̃m,n − ϕ̃, Ũ
m,n
∂t
[
]
′
+ ik∥ j̃∥m,n + ikθ j∥0 Ã∥m,n − Ã∥ , j̃∥
m,n
(
)
[ ]
2
− x, p̃ m,n + µ ∇⊥ Ũ
,
(11)
m,n
[
]
∂
Ã∥m,n = − ik∥ ϕ̃m,n − ikθ ϕ′0 Ã∥m,n + Ã∥ , ϕ̃
m,n
∂t
[
]
′
+ iδk∥ p̃m,n + iδkθ p0 Ã∥m,n − δ Ã∥ , p̃1403016-2m,n− η∥ j̃∥m,n ,(12)

Plasma and Fusion Research: Regular Articles

Volume 19, 1403016 (2024)[]∂p̃m,n = − ikθ ϕ′0 p̃m,n + ikθ p′0 ϕ̃m,n − ϕ̃, p̃
m,n
∂t
(
[
] )
+ 2δβ ik∥ j̃∥m,n + ikθ j′∥0 Ã∥m,n − Ã∥ , j̃∥
m,n
[
]
(
)
2
+ 2β x, ϕ̃ − δ p̃
+ χ ∇⊥ p̃
,
(13)
m,n

m,n

∫r
Operating r dr r to Eq. (22) gives0

1 ∂− 2 ∆W = hconv + hEM + hcurv + hdiﬀ ,
4π ∂t
with

where kθ = m/r, k∥ = m/q − n, q is the safety factor (1/q =−(1/r)(dA∥0 /dr) in our normalization), the prime indicate
the radial derivative, and (∇2⊥ f˜)m,n = (1/r)∂r (r∂r f˜m,n ) −kθ2 f˜m,n .
The energy conservation law of Eqs. (1)-(3) in termsof the complex Fourier amplitude is given by
dE
=L,
dt

∫ r[
]
ϕ̃, p̃ rdr,
0,0
r0
∫ r[
]
hEM = −2δβÃ∥ , j̃∥ rdr,
0,0r
∫ r [0]x, ϕ̃ − δ p̃ rdr,
hcurv = 2β0,0r
∫ r (0)hdiﬀ = χ∇2⊥ p̃ rdr,
hconv = −

(14)

0,0

r0

with
E = E ϕ + E A∥ + E p ,(15)

L = S ϕ′0 + S j′∥0 + S p′0 + S µ + S η∥ + S χ ,(16)where the components of each perturbation energy are defined by∑∑1∫ 1 (
) 2
Eϕ =Eϕm,n =
∇⊥ ϕ̃
rdr,
(17)
m,n
2 r0
m,n
m,n
∑
∑1∫ 1 (
) 2∇⊥ Ã∥
rdr, (18)E A∥ =E Am,n
=
∥
m,n
2 r0
m,n
m,n
∑
∑ 1 ∫ 1
2
Ep =
E m,n
=
p̃m,n rdr,
(19)
p
4β
r0
m,n
m,n

(23)

(24)
(25)
(26)
(27)

where {hconv , hEM , hcurv , hdiﬀ } are the heat fluxes associatedwith the convection, the electromagnetic perturbation, thecurvature, and the thermal diﬀusion, respectively.
For later use, we also derive an expression for extracting[ the] modes that drive the convective heat transport.
First, ϕ̃, p̃
in Eq. (24) can be rewritten as0,0[]ϕ̃, p̃

0,0

=−

∑

ikθ

m,n

)
∂ (
ϕ̃m,n p̃∗m,n ,
∂r

(28)

where p̃−m,−n = p̃∗m,n is used, and f ∗ is the complex conjugate of arbitrary f . Substituting Eq. (28) into Eq. (24)gives∑hconv =
hm,n(29)
conv ,
m,nwith |(∇⊥ f˜)m,n |2 = |∂ f˜m,n /∂r|2 + kθ2 | f˜m,n |2 for arbitrary f˜.
In addition, {S ϕ′0 , S j′∥0 , S p′0 } are the energy sources associated with radial gradients of the the electrostatic potential,
the parallel current, and the pressure, respectively, involvedin the linear terms in Eqs. (11)-(13), and {S µ , S η∥ , S χ } are
the energy sinks associated with each transport coeﬃcient.
Note that contributions of various nonlinear terms andquasilinear terms associated with the toroidal curvature arecancelled out each other in the final form of the energy conservation law. For later use, we define the pressure energy
perturbation summed over n as
∑
E m,n
(20)
E mp =
p .
n

with
∗
hm,n
conv =imϕ̃m,n p̃m,n .

(30)

Summing over n and considering the complex conjugation
property, the convective heat transport is also separated as∑mhconv =
Hconv
,
(31)
m>0

with
m
Hconv
=

∑(
−m,−n )
hm,nconv + hconv
n

=

∑

(
)
2Re imϕ̃m,n p̃∗m,n .

(32)

n

In this study, the stored energy is defined by∫ r()
p0 + p̃0,0 rdr = W0 − ∆W,
W =4π2r0

2.3 Simulation condition(21)∫rwhere W0 = 4π2 r p0 rdr is the initial stored energy and0∫r∆W = −4π2 r p̃0,0 rdr is the loss of the stored energy. In
0order to derive the time evolution of the loss of the stored
energy, we write down the (m, n) = (0, 0) component of
Eq. (13) as
[
]
[
]
∂
p̃0,0 = − ϕ̃, p̃ − 2δβ Ã∥ , j̃∥
0,0
0,0
∂t
[
]
(
)
+ 2β x, ϕ̃ − δ p̃ + χ ∇2⊥ p̃
0,0

0,0

.

In this study, the equilibrium profiles are given by
q = 2 + 2r2 ,( r − r )]
)[
β (
s
p0 =
1 − r2 1 − tanh,2ϵdϕ0 = αδp0 ,

(33)
(34)
(35)

where {rs , d, α} are arbitrary. The profile of j∥0 is determined by the cylindrical MHD equilibrium in our normalization:ϵjθ Bz − ϵ 2 j∥0 Bθ = p′0 ,
(36)
(22)
2
1403016-3

Plasma and Fusion Research: Regular Articles

Volume 19, 1403016 (2024)

Fig. 1 Radial profile of the equilibrium pressure p0 and safety
factor q.

Fig. 3 Radial profile of the equilibrium poloidal flow velocity
vE0 for α = 0, 0.2, 0.4, 0.6, 0.8, and 1.

poloidal flow velocity due to nonlinear eﬀects. In order tostudy the eﬀect of the edge shear flow on the evolution of
the ballooning mode turbulence, we conduct the parameter
survey for diﬀerent values of α in Eq. (35). Figure 3 shows
the radial profile of the equilibrium poloidal E × B flow
velocity. It is observed that sharp radial gradients of theflow velocity are formed in the edge region. In addition,
the absolute value of the velocity for the case with α = 1
agrees with that of the electron diamagnetic drift velocity
defined by v∗ = −δp′0 .
Fig. 2 Radial profile of the equilibrium parallel current j∥0 .

3. Simulation Results and Analyses

where jθ = −(dBz /dr), j∥0 = (1/r)d(rBθ )/dr, and Bz =qBθ /r. Then, the determination of j∥0 is reduced to solve
the following diﬀerential equation for B2θ :
)
( ′
ϵ p0
dB2θ
2
2
=− 2 2+ f Bθ ,dr(q /r ) + ϵ 2 2
(
)
d q2
ϵ2
f =
+ .
2
dr 2r
r

(37)
(38)

In the following simulations, the parameters are chosen such that β = 0.01, δ = 0.04, ϵ = 0.3, µ = 10−5 ,η∥ = 10−5 , χ = 10−5 , rs = 0.85, r0 = 0.1, d = 0.05.
The time step is ∆t = 0.01, and the radial resolution is∆r = (1 − r0 )/256, and Fourier modes with −40 ≤ m ≤ 40
and −20 ≤ n ≤ 20 are solved.
Figure 1 shows the radial profile of the equilibriumpressure and the safety factor. There exists a sharp radial
gradient of the pressure around r = 0.85, where the resistive ballooning mode is initially destabilized. Similarly,
Fig. 2 shows the radial profile of the equilibrium parallel
current, which is the solution of the cylindrical MHD equilibrium. The poloidal E × B flow velocity is given byvE = vE0 + ṽE0,0 ,(39)

where vE0 = ϕ′0 is the equilibrium poloidal E × B flow
velocity and ṽE0,0 = ϕ̃′0,0 represents the modification of the

In the following, simulation results of the nonlinear
evolution of the resistive ballooning mode turbulence areshown, and eﬀects of the edge shear flow on the radial turbulence spreading are analyzed in detail.
Figure 4 shows the color contours of the pressure perturbation in the poloidal cross section in the absence of theedge shear flow (α = 0) at diﬀerent times, where Z is thevertical position normalized by the minor radius. The ballooning mode with m ∼ 30 is in the linear growth phase at
t = 100 (Fig. 4 (i)), and the mode is excited near the radial
position of the peak of the equilibrium pressure gradientr = 0.83. It is observed that the mode is in the nonlinear
evolution phase at t = 150 (Fig. 4 (ii)), where the ballooning structure collapses due to the energy cascade of the unstable modes and the transition to a turbulent state with acomplex structure. In more nonlinearly developed phases
at t = 400 and t = 600 (Figs. 4 (iii) and 4 (iv), respectively),
the secondary turbulence spreading occurs in the radial direction, where intermediate-scale structures with m ∼ 10
modes are formed.
Here, the intermediate-scale structures with m ∼ 10
observed in Figs. 4 (iii) and 4 (iv) are considered to be socalled streamer structures. The group of the streamers has
the ballooning structure, where multiple poloidal/toroidal
modes are combined. It is confirmed that m ∼ 10 modes in
the initial growing phase as Fig. 4 (i) are marginally stable.
Therefore, the development of the streamers is a nonlinear

1403016-4

Plasma and Fusion Research: Regular Articles

Volume 19, 1403016 (2024)

Fig. 6 Radial profiles of the pressure p0 + p̃0,0 for α = 0 and
α = 1.

Fig. 4 Color contours of the pressure perturbation in thepoloidal cross section in the absence of the edge shear
flow (α = 0 in Fig. 3) at (i) t = 100, (ii) t = 150, (iii)t = 400, and (iv) t = 600.

Fig. 7 Radial profiles of the poloidal flow velocity vE = vE0 +ṽE0,0 for α = 0 and α = 1, and vE0 for α = 1.

Fig. 5 Color contours of the pressure perturbation in thepoloidal cross section in the presence of the strong edge
shear flow (α = 1 in Fig. 3) at (i) t = 100, (ii) t = 150,(iii) t = 400, and (iv) t = 600.

phenomenon and occurs simultaneously with the deformation of the pressure profile. Since the streamer structureis less susceptible to the diﬀusion damping compared withthe micromodes, its radially extending structure is suitablefor driving global heat transport. As shown below, it isconsidered that the system spontaneously generates such
structures for eﬃcient heat transport.
Figure 5 shows the color contours of the pressure perturbation in the poloidal cross section in the presence ofthe strong edge shear flow (α = 1) at diﬀerent times. Evenin the presence of the strong edge shear flow, similar to
Fig. 4, the ballooning structure at t = 100 (Fig. 5 (i)) andthe transition to the turbulent state at t = 150 (Fig. 5 (ii))
are also observed. On the other hand, in more nonlinearly
developed phases at t = 400 and t = 600 (Figs. 5 (iii) and5 (iv), respectively), the secondary turbulence spreading as
observed in Fig. 4 becomes very weak, and the fluctuations
are mainly localized in the edge region. Therefore, in ourparameter regime, it is found that the edge shear flow does
not strongly aﬀect the transition process of the ballooning
mode to the turbulent state, but tends to suppress the radial
spreading of the turbulence.
Figure 6 shows the radial profile of the pressure p0 +p̃0,0 at t = 600. Compared with the initial state p0 , it is observed that the pressure decreases in the inner region andincreases in the outer region, caused by the global turbulent transport. Comparing the cases of α = 0 and α = 1,
such global transport tends to be reduced by the edge shear
flow, which is caused by the suppression of the radial turbulence spreading. Corresponding to Fig. 6, Fig. 7 shows
the radial profile of the poloidal flow velocity vE0 + ṽE0,0
at t = 600. For the case with α = 0, the poloidal flow isformed even in the core region, which is caused by the ra-

1403016-5

Plasma and Fusion Research: Regular Articles

Volume 19, 1403016 (2024)

dial turbulence spreading. On the other hand, for the casewith α = 1, deformation of the poloidal flow profile is suppressed. Therefore, the strong edge shear flow is found tohave reduction eﬀects of the radially global heat and momentum transport.
Figure 8 shows the time evolution of the energy loss
rate for diﬀerent values of α. Here, from Eq. (21), the energy loss rate is defined as follows:ηloss =

∆W,
W0 r=0.83(40)

where r = 0.83 is the radial position of the peak of theequilibrium pressure gradient. Since the linear growth rateof the ballooning mode is high compared with the shearing rate of the edge shear flow in our parameter regime,
the energy loss rate is not sensitive to the value of α in thelinear regime 0 ≤ t ≲ 150. In the first nonlinear saturation phase, 150 ≲ t ≲ 250, the energy loss rate weakly
decreases as the value of α increases. In the subsequent
nonlinear evolution phase, 250 ≲ t, the energy loss ratestrongly depends on the value of α. In particular, the energy loss rate for α = 1 is roughly half of that for α = 0. In
addition, from comparison of the small α cases (δ = 0.04)
with the δ = 0 case (v∗ = 0 and vE0 = 0), it is observedthat the electron diamagnetic eﬀect strongly enhances theenergy loss in the nonlinear phase. This implies that theelectron diamagnetic eﬀect plays an important role in theformation of the streamers.
Here, we consider the suppression mechanism of thestreamer structures by the edge shear flow. The time required for the streamer to extend globally in the radial direction (for example, at least, radial length 0.5) is definedas the characteristic time of the streamer τstream . From
Figs. 4 (iii) and 4 (iv), we can estimate τstream ∼ 200. If theshear flow suﬃciently twists the streamer in the poloidal
direction during the time that the streamer extends in theradial direction, it is considered that the development of thestreamer can be suppressed. The length of the twist in thepoloidal direction is r∆θ, where ∆θ is the typical twist an-

Fig. 8 Time evolution of the energy loss rate for α = 0, 0.2,
0.4, 0.6, 0.8, and 1. The case with δ = 0 is also shown
as a reference.

gle. Considering that ∆θ is the angle that suﬃciently twists
the m ∼ 10 mode, that is, the angle that mixes the phases of
the neighboring streamers, gives ∆θ = 2π/m ∼ π/5. Therefore, the characteristic time of twisting due to the shear
flow is calculated as follows:
τE0 =

r∆θ
100
∼
.
|vE0 | r=0.83
α(41)

On the time scales longer than τE0 , it becomes diﬃcult to
produce structures cross the edge shear flow. This is consistent the simulation results, Figs. 4, 5 and 8, that α ∼ 0.5
is suﬃcient to prevent the formation of the streamer structures. Therefore, it is found that the strong edge shear flow
suppresses the development of the streamer structures inthe ballooning mode turbulence.
In order to understand the detailed physical mechanism of the reduction of the energy loss rate due to the edge
shear flow in Fig. 8, we separate the heat flux into contributions of the convection, the electromagnetic perturbation,
the curvature, and the thermal diﬀusion as shown in theright-hand side of Eq. (23). Figure 9 shows the radial profiles of the decomposed heat fluxes. In the absence of theedge shear flow (Fig. 9 (a)), it is observed that the contribution of the convection is dominant. This result indicatesthat the radial turbulence spreading is mainly caused by theconvective energy transport. On the other hand, the contribution of the convection becomes drastically small in thepresence of the strong edge shear flow (Fig. 9 (b)), which

Fig. 9 Radial profiles of the heat fluxes associated with the convection hconv , the electromagnetic perturbation hEM , thecurvature hcurv , and the thermal diﬀusion hdiﬀ for (a)α = 0 and (b) α = 1 at t = 400.

1403016-6

Plasma and Fusion Research: Regular Articles

Volume 19, 1403016 (2024)

Fig. 10 Color contours of the convective heat flux spectrummHconv(r) in r-m space for α = 0 at (i) t = 150 and (ii)t = 400.

results in the reduction of the energy loss in Fig. 8. Comparing Figs. 9 (a) and 9 (b), it is observed that the contribution of the electromagnetic perturbation is negligible, andthe contribution of the thermal diﬀusion tends to partiallycancel out that of the convection.
Next, in order to decompose the convective heat fluxminto poloidal mode components, we use Hconv(r) definedmby Eq. (32). Figure 10 shows the color contours of Hconv(r)
in r-m space in the absence of the edge shear flow. In the
beginning of the nonlinear phase at t = 150 (Fig. 10 (i)), it
is observed that various modes contribute to the convective
heat flux. On the other hand, in the phase of the radial
turbulence spreading at t = 400 (Fig. 10 (ii)), contributionsof the low-m modes are important, which is associated withthe streamer structures in Fig. 4. Similarly, Fig. 11 showsmthe color contours of Hconv(r) in r-m space in the presenceof the strong edge shear flow. Comparing Figs. 10 and 11,
it is found that the convective heat flux spectrum is almost
the same in both cases at t = 150, but the contributionsof the low-m modes are less important in the case with thestrong edge shear flow.
The results in Figs. 10 and 11 imply that the strongedge shear flow tends to avoid the inverse energy cascadein the mode number space. In order to verify this mechanism, we use E mp defined by Eq. (20). Figure 12 shows thepoloidal mode number spectrum of the pressure perturbation energy in the phase of the radial turbulence spreading at t = 400. In the absence of the edge shear flow
(Fig. 12 (a)), the peak in the range of 10 ≲ |m| ≲ 20 isobserved, which corresponds to the streamer structures inFig. 4. On the other hand, such prominent peak in theFig. 11 Color contours of the convective heat flux spectrummHconv(r) in r-m space for α = 1 at (i) t = 150 and (ii)t = 400.

Fig. 12 Poloidal mode number spectrum of the pressure perturbation energy E mp for (a) α = 0 and (b) α = 1 at t = 400.

poloidal mode number space is absent in the case with thestrong edge shear flow (Fig. 12 (b)).

4. Summary

1403016-7

In this study, we perform nonlinear simulations ofPlasma and Fusion Research: Regular Articles

Volume 19, 1403016 (2024)

the ballooning mode turbulence using a reduced two-fluid
model, and investigate how the edge shear flow aﬀects thenonlinear evolution of the ballooning mode turbulence. In
particular, from the perspective of the radial spreading ofballooning mode turbulence, the loss of the stored energy
due to the ballooning mode turbulence is analyzed in detail. Based on the evolution equation of the change of thestored energy, a method to analyze the heat flux in the radial and mode number space is newly developed. Then, it
becomes possible to quantitatively investigate the physicalmechanism behind the relationship between the turbulence
spreading and the edge shear flow.
In our simulation parameters, for the case without theequilibrium edge shear flow, the ballooning mode turbulence leads to the global heat transport by forming the secondary streamer structures. On the other hand, it is observed that, when the edge shear flow becomes stronger,
such structures are not prominent and the global heat transport is suppressed. More specifically, it is found that theloss of stored energy roughly becomes half when the edge
shear flow is as strong as the diamagnetic drift. Decomposing the heat flux into its components reveals that theconvective heat transport is responsible for the global heattransport. Furthermore, by extracting the poloidal mode
components that drive the convective heat transport, it isconfirmed that global heat transport is mainly caused by
the streamer structures. By comparing the time scale forthe development of the streamer structures and the characteristic time scale of the edge shear flow, it is newly revealed that there is a physical mechanism that the edge
shear flow prevents the formation of the streamer structures.
In the following, we compare our simulation resultswith the preceding works. It is shown in Ref. [8] that theradial spreading of the ballooning mode turbulence is suppressed by the edge shear flow. In this study, by analyzing the heat flux in the radial direction and mode number
space, it is newly revealed that the edge shear flow suppresses not only the local convective heat transport but the
global convective heat transport by aﬀecting the secondary
nonlinear evolution of the ballooning mode turbulence. In
addition, the interaction between the spontaneously generated zonal flow and the ballooning mode turbulence in theabsence of the equilibrium edge shear flow is investigatedin Ref. [9]. Our simulation results imply that the formationof the zonal flow and its influence on the ballooning mode
turbulence are less important in the presence of the strong
equilibrium edge shear flow.
In this study, we neglect the ion diamagnetic drift, theparallel velocity, and the parallel thermal transport for simplicity, but those may aﬀect the physics of the turbulence
spreading. It is also important to consider a more realistic poloidal cross-section shape for better understanding of
the ballooning mode turbulence in the advanced tokamaks.
In our simulations, we use the relatively high resistivity.
Therefore, it is also important to check how the choice ofthe resistivity aﬀects the physical mechanism of the eﬀectof the edge shear flow on the radial spreading of the ballooning mode turbulence. Those are left as future works.

Acknowledgments
The simulation code has been mainly developed using the facilities of the Center for Cooperative Work on
Computational Science, University of Hyogo. We would
like to acknowledge a Collaborative Research Program of
Research Institute for Applied Mechanics, Kyushu University.
[1] G.T.A. Huijsmans, C.S. Chang, N. Ferraro, L. Sugiyama,
F. Waelbroeck, X.Q. Xu, A. Loarte and S. Futatani, Phys.
Plasmas 22, 021805 (2015).
[2] J. Cheng, J.Q. Dong, K. Itoh, S.-I. Itoh, L.W. Yan, J.Q. Xu,
M. Jiang, Z.H. Huang, K.J. Zhao, Z.B. Shi, W.L. Zhong, S.
Inagaki, T. Kobayashi, K. Ida, Y. Kosuga, M. Sasaki, M.K.
Han, Z.X. Wang, M. Xu, Y. Xu and HL-2A team, Nucl.
Fusion 60, 046021 (2020).
[3] M. Furukawa and S. Tokuda, Phys. Rev. Lett. 94, 175001(2005).
[4] N. Aiba, C. Giroud, M. Honda, E. Delabie, S. Saarelma,
L. Frassinetti, I. Lupelli, F.J. Casson, S. Pamela, H. Urano,
C.F. Maggi and JET Contributors, Nucl. Fusion 57, 126001(2017).
[5] C.K. Sun, X.Q. Xu, C.H. Ma and B. Li, Phys. Plasmas 25,
082106 (2018).
[6] J.W. Connor, R.J. Hastie, H.R. Wilson and R.L. Miller,
Phys. Plasmas 5, 2687 (1998).
[7] X.Q. Xu, B.D. Dudson, P.B. Snyder, M.V. Umansky, H.R.
Wilson and T. Casper, Nucl. Fusion 51, 103040 (2011).
[8] P.W. Xi, X.Q. Xu, X.G. Wang and T.Y. Xia, Phys. Plasmas
19, 092503 (2012).
[9] H. Seto, X.Q. Xu, B.D. Dudson and M. Yagi, Phys. Plasmas
26, 052507 (2019).
[10] P. Beyer, S. Benkadda, G.F. Chaudier, X. Garbet, Ph.
Ghendrih and Y. Sarazin, Plasma Phys. Control. Fusion 49,
507 (2007).
[11] A. Thyagaraja, M. Valovič and P.J. Knight, Phys. Plasmas
17, 042507 (2010).
[12] M. Leconte, P. Beyer, X. Garbet and S. Benkadda, Phys.
Rev. Lett. 102, 045006 (2009).
[13] S. Nishimura and M. Yagi, Plasma Fusion Res. 6, 2403119(2011).
[14] M. Bécoulet, F. Orain, G.T.A. Huijsmans, S. Pamela,
P. Cahyna, M. Hoelzl, X. Garbet, E. Franck, E.
Sonnendrücker, G. Dif-Pradalier, C. Passeron, G. Latu, J.
Morales, E. Nardon, A. Fil, B. Nkonga, A. Ratnani and V.
Grandgirard, Phys. Rev. Lett. 113, 115001 (2014).
[15] D. Chandra, A. Thyagaraja, A. Sen and P. Kaw, Nucl. Fusion 57, 076001 (2017).
[16] X. Garbet, L. Laurent, A. Samain and J. Chinardet, Nucl.
Fusion 34, 963 (1994).
[17] F. Hariri, V. Naulin, J.J. Rasmussen, G.S. Xu and N. Yan,
Phys. Plasmas 23, 052512 (2016).
[18] Y. Kosuga, S.-I. Itoh, P.H. Diamond and K. Itoh, Phys. Rev.
E 95, 031203(R) (2017).
[19] R.D. Hazeltine, M. Kotschenreuther and P.J. Morrison,
Phys. Fluids 28, 2466 (1985).

1403016-8