Proceedings of ITC/ISHW2007

Nonlinear simulation of collapse phenomenon in helical plasma
with large pressure gradient
Naoki MIZUGUCHI1,2) , Yasuhiro SUZUKI1) and Nobuyoshi OHYABU1,2)
1)2)

National Institute for Fusion Science, Toki 509-5292, JapanThe Graduate University for Advanced Studies (SOKENDAI), Toki 509-5292, Japan

Nonlinear magneto-hydrodynamics (MHD) simulations are executed for a helical plasma with large pressure
gradient to investigate the collapsing event that is induced by the MHD instabilities. The simulation results showthat the ballooning-like instabilities with intermediate spatial scale induce both the disordering of the structuresin the barrier region and the drop of the central pressure. It has been revealed that the core pressure fall is related
to the disordering of the magnetic field structure. The simulation results are compared qualitatively with theexperimental observations of the collapse events in the super dense core state of the Large Helical Device (LHD).
Compared with our previous simulation studies for the similar situations in the spherical tokamak (ST) plasma,
the results for the helical cases show milder collapse than those in the ST case.
Keywords: magneto-hydrodynamics, collapse, helical system, simulation, super dense core

1Introduction

1+µ(∇2 v + ∇(∇ · v)),
3In recent experiments in the Large Helical Device (LHD),
high performance confinement with a super-dense core(SDC) inside the internal diffusion barrier (IDB) has beenachieved by using the pellet injection scheme[1]. On there-heating stage in such a discharge, where the plasma beta
gradually increases after ceasing the pelletinjenction, anabrupt flushing of the central density occurs in some cases.
Such a collapse phenomenon, which is named “core density collapse (CDC)”[2, 3], should be avoided since it often
limits the density increase. Although the physical mechanism of CDCs remains unclear up to the present, several
models based on the magneto-hydrodynamic (MHD) instability or the equilibrium limit have been proposed. To
reveal the mechanism of CDC is one of the key issues forfurther development of the confinement toward the highdensity branch.
This article reports a first attempt to carry out a nonlinear MHD simulation in helical plasma with large pressure gradient, which corresponds to the IDB-SDC state ofLHD. The result obtained by the simulation would provide
us with basic understandings for the collapse phenomena,including the CDC, from the point of view of MHD.

2

Simulation model

We solve the time development of the standard set of thecompressive, resistive, nonlinear MHD equations,∂ρ∂t= −∇ · (ρv),∂(ρv) = −∇ · (ρvv) − ∇p + j × B∂tauthor’s e-mail: mizu@nifs.ac.jp

(1)

(2)

∂B
= −∇ × E,(3)∂t∂p= −∇ · (pv) − (γ − 1)(p∇ · v + ηj2 ), (4)∂tin a full-toroidal three-dimensional geometry. The variables, ρ, v, B, and p represent the mass density, the fluid
velocity, the magnetic field, and the pressure, respectively.
The current density j and the electric field E are calculated
from Maxwell’s equations,j = ∇ × B,

(5)

E = −v × B + ηj.

(6)

u1

= r cos (θ − hφ),(7)

= −φ,(9)In (1)-(4), the dissipation terms are included as the resistivity η and the viscosity µ. These terms are assumed to beuniform constants for simplicity. The Ohmic heating term
is evaluated in the evolution of the pressure p, whereas theviscous heating is ignored for simplicity. All the spatial
derivatives are expressed numerically by using the fourthorder central-difference scheme. The time integration issolved by using the fourth-order Runge-Kutta method.
To follow the geometry of the helical devices withcontinuously-wound magnetic coils, we adopt the helicaltoroidal coordinate system which is used in the HINT
code[4]. This coordinate system (u1 , u2 , u3 ) uses a rectangular grid (u1 , u2 ) in the poloidal cross section and atoroidal grid (u3 ) with geometrically regular intervals. By
using the poloidal angle θ, toroidal angle φ, and the minor
radius r, (u1 , u2 , u3 ) is described as[4]u2u3= r sin (θ − hφ),(8)where the number of helical period h is 10/2 to follow theLHD configuration. The boundary condition of the com-

Proceedings of ITC/ISHW2007
case Acase Acase Cenergypressure case Bcase B   major radius[m] Fig. 1 Radial pressure profiles of the initial equilibria for thecase A, B, and C on the horizontally elongated cross section.

putation for the u1 and the u2 directions is the perfect conductor condition, and that for the u3 direction assumes theperiodicity 2π. The number of the numerical grid points is(N1 , N2 , N3 ) = (62, 146, 500).
The initial condition for the simulation is given by thenumerical solution of the HINT code which models an average experimental configuration of LHD[5]. The simulation geometry includes the region out of the separatrix. Totreat such an external region continuously from the insideof the separatrix, we artificially add low uniform component in pressure initially.
The simulation starts by adding tiny random perturbations in the velocity component of the initial equilibrium.
Thus, the spontaneous time development of the MHD system is solved by using vector-parallel calculation on a supercomputer.

3

time 

Simulation result 

R



Fig. 3 Poloidal mode structure of the perturbations in the pressure on the horizontally elongated cross section.

reference. One can see that the growth rate for case A ismuch larger than that for case B. Also, the resistivity scansshow the independence of the growth rate on the resistivity, which proves the modes to be ideal ones. These modes
have intermediate poloidal mode numbers m ∼ 10. The
poloidal mode structures, shown in Fig.3, clearly exhibit
the ballooning nature in that the fluctuations are localized
in the outer region.
The nonlinear long-term evolution of the energy andthe maximum pressure is shown in Fig.4. One can see thatthe growth and saturation of the energy repeats three timesbefore reaching a relaxed state. The growth shown in Fig.2
corresponds to the first one (t < 100τA ). During this relaxation process, the plasma changes its shape gradually. The
temporal changes in the radial and the poloidal pressure
profile are shown in Fig.5 and Fig.6, respectively. Since theprimarily induced instability is localized in the outer regionof the torus, as described above, the crash of the structures
is first seen in the barrier region (see Fig.6(b)). The plasma
1.0x10-4.22.8.20.18.6.16.4

P max

The MHD equations (1)-(4) are solved in a dimensionlessform by normalizing all the variables with characteristicvalues, such as 1[m] for the spatial scale, 1[T] for the magnetic field. Uniform density ρ = 1 is given for the initial condition for simplicity. The uniform constant for thedissipation coefficient η = 10−4 and µ = 10−3 are usedthroughout this article. The unit of time, i. e., the Alfvèntransit time τA corresponds to the order of ∼ µsec for thestandard parameters of LHD experiment. For the main results described in this article, the initial equilibrium whichroughly models the IDB-SDC state of LHD just before
the CDC events is used. Its pressure profile is shown in
Fig.1 (case A). The maximum and the volume averaged
beta value are β0 = 6.6% and < β >= 1.8%, respectively.
Compared with the standard high-beta discharges of LHD
experiments (case C), the pressure gradient in the barrier
region is markedly large due to the large Shafranov shift.
Figure 2 shows the time evolution of the kinetic energy in solving the linearized system of (1)-(6). One cansee a linear growth of the instability at t > 50 for case A.
In Fig.2, the trace for the case which uses a lower beta
equilibrium with β0 = 2.5% (case B) is also plotted for Fig. 2 Time development of kinetic energy for the linear simulation.

Z 

kinetic energy .14.2.120.10

0

200

time

400

Fig. 4 Time evolution of the total kinetic energy and the maximum pressure.

Proceedings of ITC/ISHW2007

pressure.2

t= 0 τA
200
550.10
2.5

3.0

3.5

4.0

major radius [m]

4.5

5.0

Fig. 5 Temporal change of the radial pressure profile on the horizontally elongated cross section.

surface is deformed, reflecting the linear eigenmode structures in the outer region, and the magnetic field structure in
the inner region, as shown in Fig.3. Part of plasma is lost
due to the disturbance. The resultant pressure gradient in
the barrier region becomes steeper than in the initial state(see Fig.5 at t = 200τA ). It should be noted that the mode
structures of the instability are located only in the barrierregion, whereas the central pressure gradually decreases,
as the lost plasma forms a pedestal pressure in the edge region, as shown in Fig.5 toward t = 550τ A and Fig.6(b)-(d).
If one see the time development of the maximum pressure
at the core as shown in Fig.4, there is an abrupt change inthe trace after t = 300τA . The rapid fall in the core pressure
might be related to the change in the magnetic field structure. Figure 7 shows the time development of the magnetic
field structure which is expressed by the puncture plot of
the field lines. The field line trace is executed by using
a sixth-order Runge-Kutta method together with a fourthorder interpolation of the magnetic field. The magnetic
surface structure is clearly formed entirely from the coreto the edge region initially (see Fig.7(a)). Such a nestedsurface structure is sound during the early stage of the relaxation process, although the edge structure is markedly
deformed, as shown in Fig.7(b). At time t = 300τ A (seeFig.7(c)), the magnetic surface structure abruptly diminishes throughout the whole poloidal plane. Finally, thestructure reappears in the core region as shown in Fig.7(d).
At the moment shown in Fig.7(c), the plasma in the core region at high pressure is linked to the external low-pressure
region with an identical field line. Under this situation, theplasma outward flows due to the pressure imbalances along
the field lines, which might cause the rapid fall of the corepressure, can be induced.
The system reaches a relaxed state within ∼ 600τ A
in this result. The resultant pressure profile becomes a
broader one because of the decrease in the core, and thespread in the edge, as shown in Fig.5 at t = 550τ A .

4Discussion

The simulation result described in Sec.3 can be compared
qualitatively with the experimental observations on theFig. 6 Temporal change of the pressure profile. The iso-contour
maps of the pressure on the horizontally elongated crosssection are drawn. The time equals to (a)0, (b)130,(c)300, and (d)550 τA .

CDC event in LHD. The time scale of the whole simulation process roughly agrees with that of the crash phaseof CDCs. There remains some several questions for themechanism of occurrence of CDCs. The simulation results might provide us with interpretation for them. Firstly,
the reason why the core density decreases is a basic question for CDCs. The fact that only the density decreases,
keeping the core temperature unchanged, implies that thecollapse is governed not by the conductive processes, butby the convective ones. In our simulation result, the coreplasma is extracted convectively through the transiently
disordered or reconnected field lines as shown in Fig.7.
This convective loss mechanism is valid even if there isno significant unstable mode in the core region. The possibility of the coexistence of the instability in the barrier
region and the significant drop in the core pressure is also
comparable with another question why the edge fluctuation
is observed just before the density collapse in the experiment. The simulation result gives a reasonable scenario forit. More detailed comparison would need systematic scans
for the edge stability problem both in the simulation andexperiment in future.
Now we compare the simulation result with the previous simulation which studies the nonlinear dynamics ofST plasma on an ELM crash[6]. Although both the simulations are initiated by the ballooning mode instability, andProceedings of ITC/ISHW2007Fig. 8 Three-dimensional structure of plasma pressure at t =365τA . The iso-contour of the pressure together with an
iso-pressure surface are drawn.

5Fig. 7 Temporal change of the magnetic field. The puncture
plots of the magnetic field lines on the horizontally elongated cross section are plotted. The time equals to (a)0,(b)130, (c)300, and (d)550 τA .

have several points in common, there is also qualitative difference between them. To compare both simulations would
help us understand the key properties of the process. InRef.[6], the crash process of the pressure profile is divided
in two stages. The edge region is perturbed in the formerstage, whereas the core pressure falls in the latter one, as isalso observed in our simulation. However, for the ST case,the core collapse is caused by a secondary induced internal low-n instability, where n is the toroidal mode number,in contrast to the direct parallel loss in the helical case described in Sec.3. No prominent secondary low-n structure
can be observed in the helical case simulation, as shown in
Fig.8. The resultant temporal change of the core pressure
profile shows milder activity for the helical case than for
the ST case. The large internal plasma current of the ST
would make difference in the core instability from the helical system, because the induced internal mode is a current
driven one.

Summary

We have executed MHD simulations and revealed the nonlinear dynamics of the ballooning mode in helical systems
with large pressure gradient. The simulation result has
shown a physical mechanism that the pressure collapse in
the barrier region is followed by the convective loss fromthe core region through the disturbance of the magnetic
field. The simulation is qualitatively comparable with theexperimental observations of the CDC events in LHD. Although the simulation result shows a closed scenario initiated by the ballooning instability, it is necessary to carryout the systematic stability analyses under the experimental situation before discussing the cause of the CDC in future.
The authors would like to thank Drs. S. Ohdachi,
J. Miyazawa, R. Sakamoto, Y. Narushima, K. Itoh, K.
Ichiguchi, Y. Todo, H. Miura, and N. Nakajima for fruitful discussion. This work is performed with the supportand under the auspices of the NIFS Collaborative Research
Program (NIFS07KNXN097).
[1] N. Ohyabu, T. Morisaki, S. Masuzaki, et al., Phys. Rev. Lett.
97 55002 (2006).
[2] R. Sakamoto, et al., this conference.
[3] J. Miyazawa, et al., this conference.
[4] K. Hatafuji, T. Hayashi, and T. Sato, J. Comp. Phys. 81 169(1989).
[5] Y. Suzuki, N. Nakajima, K. Watanabe, et al., Nucl. Fusion
46 L19 (2006).
[6] N. Mizuguchi, R. Khan, N. Nakajima, and T. Hayashi, Nucl.
Fusion 47 579 (2007).