J. Plasma Fusion Res. SERIES, Vol. 5 (2002) 495499

Nonlinear MHD Simulation of Pressure Deformation
Dynamics in the Large Helical Device
MIURA Hideaki, HAYASHI Takaya and SATO Tetsuya
National Institute for Fusion Science, Toki 509-5292, Japan(Received: 11 December 2001 / Accepted: 17 May 2OO2)

Abstract
Nonlinear behaviors of a magnetohydrodynamic (MHD) plasma in the Large Helical Device (LHD)
are investigated by means of a numerical simulation. Computations for a configuration with the radius ofthe vacuum magnetic axis R* = 3.6 m are conducted under the stellarator symmetry. The MHD plasma isobserved to be dominated by a ballooning instability. Simulations with several cases of resistivities showthat tho saturation levels of the plasma kinetic energy are nearly proportional to the resistivity. The
numerical results are compared with our previous simulation results on Ru= 3.7 m. The saturation levels
are fotrnd to be almost comparable between the two configurations, although the R"* = 3.7 m
configuration is more magnetically hilly and Mercier-unstable. Pressure budget analysis shows thatviscous heating plays a key role in making pressure deformations saturated.

Keywords:
nonlinear MHD simulation, pressure deformation, viscous heating

In this paper, we aim to investigate full 3D

1. lntroduction
Understanding of complex nonlinear behaviors of

nonlinear behaviors of an MHD plasma in the LHD withthe magnetic axis at R* = 3.6 m and compare with theresults of the previous simulations with R"* = 3.7 m.
While it is known that at LHD configuration becomes
more magnetically hilly and Mercier-unstable as thean MHD plasma is important for achievement of
successful confinement of plasmas in fusion devices.
Numerical simulation is one of the most powerful andsuccessful approaches to investigate the subject. Thoughmagnetic axis is shifted inwardly [3], a series ofexperiments with Ru* = 3.6 m reports a better
confinement than a Ru* = 3.7 m configuration [4]. We
expect that a comparison of Ru^ = 3.6 m and 3.7 m
simulations allows us to investigate why and how ahuge memory and long computational time arenecessary for a nonlinear simulation on a helical device

with fully three-dimensional (3D) geometry, recent
development of high-perforrnance supercomputers has
enabled us to conduct such a large computation.
Recently, we have developed a numerical code tosimulate full 3D, compressible and dissipative MHD
equations in a helical device [1]. Our previous

good confinement is achieved with Ro = 3.6 m magnetic

axis from a viewpoint of nonlinear MHD.

In Sec. 2, the outline of our simulations isdescribed. In Sec. 3, numerical results by our nonlinear

simulations on a LHD configuration with the magnetic
axis at Ru* = 3.7 m [2], under the stellarator symmetry,simulations are shown. We investigate growth of thekinetic energy as an indicator of plasma fluctuation. We
pay also our special attention to pressure deformation

has shown excitation of unstable modes, which aredominated by the resistive ballooning instability, andsaturation of the excited modes,@2O02 by The Japan Society of Plasma

Corresponding author's e-mail: miura@toki. theory. nifs. ac. jp495

Science and Nuclear Fusion ResearchMiura H. et al., Nonlinear MHD Simulation of Pressure Deformation Dynamics in the Large Helical Device

toroidal flux. Figures l(a) and l(b) are the Poincar6plots of magnetic lines on horizontally- and vertically-because the pressure gradient should be the main source

to drive an instability in the system studied in this paper.
Summary is shown in Sec. 4.

elongated poloidal section, respectively. The outer sideof the torus is in the right-hand side. The profile of therotational transform, l, is shown in Fig. 1(c) as afunction of the mean plasma radius. From the2. Outlines of Simulations
We solve the fully-nonlinear, compressible andresistive MHD equations in a full 3D geometry. The

calculation of the specific volume, the core part of theconfiguration has a property of magnetic well while theedge region is magnetic hill. In Figs. 1(d) and 1(e),bird's eyes views of the pressure and the toroidal current
on a poloidal section Q = nl2O (@ is the toroidal angle)
are shown, respectively. The central beta, Bs, is about
3.6 7o. Note that an initial equilibrium which is used inref. 2 is roughly described by Ru* = 3.7 m, the initialstellarator symmetry is imposed in the toroidal direction.
We set the position of the vacuum magnetic axis to beRu* = 3.6 m. In our numerical code, we adopt 4th-orderaccuracy discretization schemes both in time and space.
Number of grid points are 97 x 97 on a poloidal section
and 32 in the toroidal direction for a half-pitch period.
Refer to ref. 2 on details of our simulation code. (Seealso ref. 5 on normalization of the MHD equations.)
In Fig. 1, profiles of an initial ideal equilibrium areshown. This equilibrium is obtained by the use of theHINT code [5]. An initial pressure profile is given byp(Vr) = p(0) (l - y)2 where y is the initial normalized

pressure profile being p(r4) = p(0) (l- {)2 and B6= 4

Vo. Thls differences in the results of two series ofsimulations compared in the next section are mainly
brought about by the difference in Ru" for the initialequilibria.

(c)

(e)

Fig. 1 An initial equilibrium computed by the HINT code. The Poincar6 plot of the magnetic field lines on (a) the horizontally-elongated poloidal section, (b) the vertically-elongated poloidal section. The outer side of the torus is in theright-hand side. (c) The rotational transform plot. The abscissa is the mean plasma radius. Bird's eyes views of (d)the pressure and (e) toroidal current on A = n/20 poloidal section.
496

Mir-rra H. er al., Nonlinear MHD Simulation of Pressure Deformation Dynamics in the Large Helical Device

3. Nonlinear Simulations and Pressure

1e-06

Deformations
In this section, we see numerical results of ournonlinear full 3D simulations. In the MHD equations,
there are three intrinsic parameters: the conductivity 4
resistivity 4 and viscosity 1t.ln this paper, we fix r= 1x 10-6 and F = 2 x l0-3 throughout this article and study
the influence of the resistivity by varying 4 from I x
10-6tolxlo-:.
In Fig. 2, we see time evolution of the averaged
kinetic energy per unit density Er= l/2 ('u; o;) (or is theLe-ut

1e-08

1e-09

1e-10

1000

l-th component of the velocity vector and (.) represents
the volume average) for five values of the resistivity 4-1200 1400 1600 1800
t2000

Fig. 2 Time evolution of the averaged kinetic energy.
Thick dashed, thick dotted, thick solid, thin dotted,
thin dashed lines represent simulations with 4 = 1

I x 10-6, loua x 10-6, 10r/2 x 10-6, 103/a x 10-6 and I x
l0-5. It is observed that Ey grows almost exponentiallyx 10-6, 101/4 x 10-6, 101/2 x 10$ and 1034 x 10-6, and 1
x 10's, respectively.

first, and begins to decay gradually for all of the fiveresistivities. We shall see later that the exponential
growth in Fig. 2 is dominated by the ballooning

1e-06

instability.

(a)

In Fig. 3(a), saturation levels o of the kinetic
energy growth are plotted as a function of the resistivity

4. Here we define the saturation level as the maximum
value of the kinetic energy. A least-square fitting showso n no'e. Thus o is almost proportional to 4. In Fig.
3(b), time evolutions of E11 with R* = 3.6 m (solid line)and3.7 m (dashed line) are shown. The resistivity is 4 =101/2 x 10-6 for both of the two simulations. We find that
E* of the Ru^ :: l.f m simulation grows faster than thatof the /?u^ = 3.J m simulation. It is a direct result of theformer configuration being more Mercier-unstable thanthe latter. Note that the two kinetic energies saturate atalmost the same level despite the difference of thegrowth rates. We checked the levels of o for threevalues of resisr.ivity 4 - I x 10-6, l0ua x 10-6, and 10r/2 x
10-6 of R* = 3.7 m simulations and verified that o of R*= 3.7 m simulations are comparable with those of R* =3.6 m simulatrons. Furthermore, as we shall see in thenext paragraph, an instability observed in this article is
not destructive and a plasma confinement is recoveredin the latter slage of the simulations. These two facts
suggest that the saturation level of an instability is notsensitive to a detailed value of Ru*, because of anonlinear property of an MHD plasma.
In Figs. 4(a) and (b), contour plots of the pressure

1e-07

1e-08
1e-06

1200 1400 1600 1800

2000

t

on horizontally- and vertically-elongated poloidal
sections at t .= l40O tA (tA is the Alfv6n time unit),
obtained by a simulation with \ = 101t2 x 10-6, areFig.3 (a) Saturation level oas a function of the resistivity 4. Solid line represents o- 4ot. (b) A comparison of the kinetic energy growth between 4" = 3.6m and 8", = 3.7 m simulations. Solid and dashedshown, respectively. The outer side of the torus is in theright-hand side. Contours are deformed in the course oflines represent 8", = 3.6 m and 8", = 3.7 m simulations, respectively.

the time evoltrtion at around t = 2/3 resion where the497

Miura H. et a/., Nonlinear MHD Simulation of Pressure Deformation Dynamics in the Large Helical Device

4e-09(a)3e-09I2e-09

IIIIt\\\\,,II

1e-09
0-5e-101000 1200 1400 1600 1800 2000

tFig. 4 Contour plots of the pressure on horizontally- andvertically-elongated poloidal sections. Outer sideof the torus is in the right-hand side to the paper.
The time stamp is t= 1200 ro in (a) and (b).

8e-09
6e-09
4e-09

2e-09
0-2e-09-4e-09-6e-09-8e-09pressure gradient is the steepest, especially in the outerside of the torus. Poloidal and toroidal Fourier mode
numbers are estimated to be 15 and 10, respectively.
Furthermore, pressure deformations are clearer on thehorizontally-elongated poloidal section than thevertically-elongated poloidal section. An analysis which
is similar to that seen in Miura et al. [2] reveals that thepressure deformations are dominated by the resistive
ballooning instability. It is also observed that plasma
confinement, for which the pressure profile becomes1000 1200 1400 1600 1800 2000

t

Fig.5 Pressure budget evaluation tor a q = 1012 x 10-6
simulation. The pressure flux, compression, thermal conduction, ohmic heating and viscous heating terms are represented by solid line, dottedbroader, is recovered after a sufficiently long timedevelopment.

Next, in order to clarify a physical mechanism

line, dashed line, cross and triangles, respectively.
(a) Volume integration of the five terms in theright-hand side of the pressure budget equation
over the system. (b) Volume integration of the fivewhich drives or suppresses pressure deformations
observed in Fig. 4, we calculate volume integrals of therhs. terms of the pressure-budget equation as has beenterms for their positive and negative value re-done in ref. 2. In Fig. 5(a), the pressure flux,gions.

compression, conduction, ohmic heating and viscousheating terms are shown for the period 1000 zo ( I (2000 11. It is clear that the viscous heating term (adashed line) dominates the pressure budget. Next, thefive terms are integrated for positive and negative values
separately and plotted in Fig. 5(b). It is clear that thepressure flux term has the largest amplitude. It isconsistent with a result shown in Ref. 2, in which it isshown that the pressure flux is a profile in phase withthe pressure fluctuation and drives pressure

4. Concluding Remarks
We have conducted a series of nonlinear MHD
simulations on the LHD system with the position of themagnetic axis R* = 3.6 m. Our simulations have shownthat the saturation level of the kinetic energy is nearly
proportional to 4: o* 4o'e. It implies that the saturationlevel becomes sufficiently small when the resistivity
tends to a small value available in the LHD experiments,suggesting that an instability observed in theexperiments can be quite mild. Although paying too
much attention to the precise value of the scalingexponent 0.9 does not make sense because the numberof samples are small and and range of 4 is too shortdeformations. Note that the viscous heating term grows

rapidly just before t = l4OO fa, when the kinetic energy

becomes maximum. A rapid growth just before thenonlinear saturation suggests that the viscous heatingterm plays a crucial role in the saturation.

(one decade), it still seems that the scaling exponent 0.9498

Miura H. er a/., Nonlinear MHD Simulation of Pressure Deformation Dvnamics in the Larse Helical Device

is much greater than a value expected from a growthrate ofthe resistive ballooning instability 1/3. It suggests
that, in order to explain the behaviors of o, thenonlinearity of the MHD equations (the viscous heatingeffect, for example) should be taken into account. It isuse of the supercomputer NEC SX-4/64M2 in the theory

and computer simulationcenter of NIFS.

Referencest 1l H. Miura et al., to be published in J. Plasma Fusion

worth emphasizing here again that the saturation levelsin the R* = 3.6 m simulations are almost comparablewith Ru, =3.7 m cases. So far as a system under thestellarator syrnmetry is concerned, these two facts
suggest that a plasma can be confined even with a small(shifted to inward) R* because of the nonlinear propertyof aplasma.
This simulation research was conducted bv makins

Res. SERIES 4.

I2l H. Miura et a/., Phys. Plasmas 8, 4870 (2001).
[3] Ichiguchi et al., Nucl. Fusion 33,481 (1993).
t4l N. Ohyabu et al.,Phys. Rev. Lett. 84, 103 (1999).
t5l T. Hayashi, A. Takei and T. Sato, Phys. Fluids B 4,1539 (1992); Harafuji et al.,J. Comput. Phys. El,
169 (1989).

499