Proceedings of ITC/ISHW2007

Growth of eigenmodes with moderate mode numbers in nonlinear
MHD simulation of LHD
Hideaki MIURA1) and Noriyoshi NAKAJIMA1)
1)

National Institute for Fusion Science, Toki 509-5292, Japan(Received 15 October 2007)

Direct numerical simulations of compressible, nonlinear magnetohydrodynamic equations in the fully threedimensional geometry of the large helical device are carried out to study the linear and nonlinear growth of thepressure-driven instability. Two parameter sets are adopted in the simulation study. One parameter set provides
small dissipative coefficients so that the simulation results can be compared to the ideal linear stability analysis.
By the use of this parameter set, the growth rates in the linear stage of our simulation coincides with the ideal
linear growth rates. Another parameter set is provided to study nonlinear evolutions. Excitations of flows parallelto the magnetic field, its role to the nonlinear saturations are studied.
Keywords: direct numerical simulation, MHD instability
DOI: 10.1585/pfr.1.001

1Introduction

In researches of magnetic confinement devices such as
the Large Helical Device (LHD), Magnetohydrodynamics (MHD) instability is one of the key issues to understand plasma behaviors in the devices. Aiming to clarify
the physics of the MHD activities observed in the LHD
experiments[1], many linear analysis and nonlinear simulations have been carried out [2, 3, 4, 5, 6]. Because of
the short wave instability nature of the the pressure driven
instability, small scales must be sufficiently resolved in anumerical computation. In the numerical study of an instability, the linear growth rate is an important index. In linear analysis of an ideal MHD system, high mode behaviors
can be well predicted. However, in numerical simulations,non-ideal (dissipative) MHD equations are studied because
the dissipation is indispensable to avoid numerical instability. Since the resolution available in a numerical simulation
are restricted, the dissipative coefficients in simulations areoften far larger than the values expected to compare numerical results to experimental results. Typically, the adoptionof the large coefficients brings about unfavorable dampingof high modes.
In our recent work on nonlinear MHD simulations ofLHD by the use of the MHD In Non-Orthogonal System(MINOS) code[3], the growth rate of the lowest mode coincides well with the linear ones obtained by the CAS3D
code.[4] However, the growth rates of moderate and highmodes in the simulation are apparently suppressed by theviscosity and the coincidence is lost as the mode number becomes larger. In this paper, the number of the grid
points of our simulations is increased to study the moderate modes behaviors. Numerical simulations are carried out for the equilibrium magnetic field with the vacauthor’s e-mail: miura.hideaki@nifs.ac.jp

uum magnetic axis position 3.6m and the pressure profile p(ψ) = (1 − ψ2 ) where the ψ is the toroidal magnetic
flux. Growth of the moderate mode numbers, generationof toroidal flows and nonlinear saturations of the instability are discussed. In this articles, two parameter sets areadopted: (A) the isotropic heat conductivity κ = 1 × 10−6(that is, the parallel heat conductivity κ// and the perpendicular heat conductivity κ⊥ is identical), the resistivity η = 1 × 10−6 , the viscosity µ = 1 × 10−6 , and (B)κ// = 1 × 10−2 , κ⊥ = 1 × 10−6 , η = 1 × 10−6 , µ = 1 × 10−4 .
These parameters are already normalized by some typicalquantities (see Ref.[5]) and hence the reciprocal of the dissipative number can be considered as some similarity parameters such as the Reynolds number Re = 1/µ and theLundquist number S = 1/η. The number of grid points
are 193 × 193 on a poloidal cross-section and 640 in thetoroidal direction. The number of grid points is doubled in
the two directions of a poloidal cross-section compared toour earlier works.[5, 6] The parameter set (A) provides ahigh Reynolds number and the behaviors of the solutions
may be easily compared to an ideal linear analysis. Though
the simulation with this parameter set cannot be completed
because the numerical resolution (number of grid points)
is not sufficient for such a high Reynolds number flow, it
serves to study linear stages. The parameter set (B) provides the moderate Reynolds number and nonlinear saturation within the numerical resolution, although the linear
stage is influenced by the dissipative coefficients µ and κ// .
Thus we first study the former case, parameter set (A) sothat the linear growth are clearly seen. Then the nonlinear
growth is studied for the parameter set (B).

Proceedings of ITC/ISHW2007

2Comparison to linear ideal analysis

Here we study the growth of Fourier amplitudes of the velocity vector. The velocity vector field V is decomposed
into the normal (V∇ψ ), the parallel (Vb ) and the binormalcomponents (V∇ψ×b ) Then each of the three vector components are decomposed into the Fourier modes m and n on
the Boozer coordinate. By the use of the Fourier
coeffiP Rψ 2cient Vα,mn (ψ), the power spectrum Aα,n = mVi,mn dψ
is defined.[5] In Fig.1(a), the time evolution of the volumeintegrated amplitudes of the three vector components ofthe velocity vector field are shown for the parameter set(A). Hereafter, the abscissa t is normalized by the toroidal
Alfven time τA whenever a time series is plotted. InFig.1(b), the time evolutions of Aα,n are shown. Here werecall that a single Fourier mode on the Boozer coordinate
is not the linear eigen-function of the MHD equations byitself. A linear eigen-function consists of multiple Fouriermodes, as are shown by Nakajima[8]. However, as isshown in the reference, the toroidal coupling in the LHD isso weak that the poloidal Fourier wavenumber n can be approximately recognized as a good quantum number. Basedon this understanding, we see only the n-wavenumber dependence of the Fourier amplitude growth by summing mover all wavenumbers. Fig.1(c) is the magnification of thelinear stage of this simulation. Observe in Figs.1(b) and(c) that some of the Fourier amplitudes achieve exponential growth while some show further acceleration. Becauseof the high Reynolds number 1/µ = 106 , the computation
is terminated in the midst of the growth. However, observations of the growth of the energy faster than an exponential growth in Figs.1(a) and (b) reveal that the final time 55
is already in the nonlinear stage of the evolution. We alsofind that the growth of the parallel energy coincides withthose of the other two components only in short period at40 < t < 45, or rather faster than them. As is commented inthe above, the final period of this simulation is influenced
by nonlinear coupling in the MHD equations. Deviationsof the parallel energy growth from the other two components should be attributed to the nonlinear growth. Basedon studies of the growth of separate m/n Fourier modes( figures of which are omitted here), the linear growth isconsidered in between t ' 30 and t ' 40.
In Fig.2, the growth rates obtained in this simulation
are compared to the ideal linear analysis. The CAS3D3[?]
computation provides estimates of the growth rates. The
growth rates for n ≤ 4 in the simulation shows reasonablecoincidence to the CAS3D computation. Compared to ourprevious work[5], the number of grid points are doubled
into the two directions of a poloidal cross-section. Although the parameters are the same as the ones in the previous work, the artificial dissipations due to the truncation is
drastically decreased thanks to the high resolution properties of the compact scheme.[10] Nonetheless, we still haveto recognize that the moderate Fourier modes n > 7 is under the influence of the dissipative coefficients. We also

note that µ = 1 × 10−6 is the minimum value which makes
sense as a physical viscosity for the current numerical resolution. Although we could reduce the viscosity below
1 × 10−6 , the growth of unstable modes becomes rather
insensitive to the change of µ, suggesting that the numerical viscosity associated with the truncation errors become
gradually dominant.

3

Nonlinear evolution

Next, the nonlinear evolutions of the instability is studied
by the simulation with the parameter set (B). In Fig.3, thetime evolution of the volume-integrated amplitudes of thethree vector components of the velocity vector field areshown. The growth rate of the parallel flow component is
the same as the other two component, making a difference
to the observation in Fig.Fig:Run553total. It shows that thethree components obey to the same kind of the instability.
Among the three components of the velocity vector, theparallel component becomes dominant finally even though
it is initially much smaller than the other two components.
Fourier amplitudes for some toroidal Fourier modes nof (a) the normal, (b) the binormal and (c) the parallel components of the velocity vector are shown in Figs.4(a)-(c),respectively. In Figs.4(a) and (b), Fourier amplitudes associated with the odd toroidal numbers n = 1, 3, 5, 7, 9, and
11 grows the fastest. The growth of the Fourier modes aremostly contributed by the Fourier mode set m/n associatedwith the ι/2π = 0.5 rational surface. Growth rates of these
Fourier modes are almost identical. Although the separation between odd and even wavenumbers are not clarified
yet, the exponential growth of the even toroidal numbers
are considered as a consequence of nonlinear couplings ofthe odd wavenumbers rather than their own linear instability, since the growth rates of these even numbers are aslarge as the twice of the growth rates of the odd numbers.
Among the four dissipative coefficients, the parallel heat
conductivity κ is considered as the most influencing one.
For example, once the unstable Fourier modes which areresonant to the ι/2π = 0.5 rational surface grow, the pressure distribution is modified because of the large κ// . Since
eigen-functions of low n modes cover wide are across theflux surfaces, the κ// effect can also modify the pressure
profile in wide area. In Fig.4(c) the time evolutions of theFourier amplitudes of the parallel velocity component areshown. In comparison to Figs.4(a) and (b) in which thetoroidal modes n = 1, 3 and 5 have almost the same amplitudes, it is clearly seen that the n = 1 Fourier modeachieves the largest amplitude. Furthermore, there is noclear distinction between even and odd mode numbers.
In Fig.5, the pressure isosurfaces and contours on apoloidal cross-sections are shown at t = 110τ A (upper)and 150τA (lower). At t ' 110τA , many of the Fouriermodes are saturated or going to be saturated. The contourlines on the poloidal cross-section in Fig.5(a) exhibits clear

Proceedings of ITC/ISHW2007
mushroom-like structures. The mushroom-like structures
have been repeatedly reported in our earlier works.[4, 5, 6]
While the structures in the earliest work consist mainlyof m/n = 2/1 Fourier modes, the structure in this article consist of many more Fourier modes which are resonant to the ι/2π = 0.5 surface. The overlapping of theresonant modes provides quite large pressure deformation.
Consequently, the contour plots of the pressure (and therefore the Poincare plots of the magnetic field lines, which
are not shown here, too) become quite chaotic after thenonlinear saturations. However, at t = 150τ A , the pressure contour associate with the outermost closed surface is
not deformed very much through the linear and nonlinear
evolution, suggesting that the confinement is not criticallydamaged in spite of the strong instability. Furthermore,
contour lines tend to form the concentric profiles again,
showing a possibility of organizing relatively well confinedstate again.

-4(a) 10-5

1010-610-7"Run551.spec_total.dat" u 1:2"Run551.spec_total.dat" u 1:3"Run551.spec_total.dat" u 1:410-8010203040(b)-410-510-6

10

10-7

n=0
n=1
n=2
n=3
n=4
n=5
n=6
n=7
n=8n=9n = 10n = 11

10-8-910-10

104Concluding Remarks

10-1110-12Direct numerical simulation study is conducted to study
linear and nonlinear evolution of the instability in LHD.
A simulation with relatively small dissipative coefficientsshows good coincidence of the linear growth rates with theideal analysis. In nonlinear simulations, the parallel flow
becomes dominant after the nonlinear saturation. It is noteworthy that the parallel flow tends to grow much larger
than the other two components of the velocity vector in thesimulation with smaller dissipative coefficients than in theone with larger coefficients. It suggests that the parallel
flow dominance is much stronger in simulations with verysmall dissipative coefficients. The numerical simulationsin this work was carried out on NEC SX-7 “Plasma Simulator” in National Institute for Fusion Science under theNIFS program NIFS06KTA033.
[1] K. Watanabe et al., Nucl. Fusion 45, 1247 (2005).
[2] H. Miura et al., Phys. Plasmas 8, 4870 (2001)[3] K. Ichiguchi et al., Nucl. Fusion 43 1101 (2003).
[4] H. Miura et al., IAEA-CSP-25/CD/TH/2-3 (2004).
[5] H. Miura et al., AIP Conference Proceedings 871 (2006).
[6] H. Miura et al., Fusion Science and Technology 51, 8 (2007).
[7] C.Nührenberg, Phys. Plasmas 6 137 (1999).
[8] N. Nakajima, Phys. Plasmas 3 4556 (1996).
[9] C. Nührenberg, Phys. Plasmas, 6 137 (1999).
[10] S.K. Lele, J. Compu. Phys. 103 16 (1992).

50t

10-1301020304050t(c)
10-4
10-5
10-6-710n=0n=1n=2n =3n=4n=5n=5n=7n=8n=9n = 10n = 10

10-8
10-9
10-10
10-11
10-12
10-1330

35404550tFig. 1 Time evolutions of (1)the volume-integrated energy of
the three components of the velocity vector, (b)theFourier amplitudes of the normal velocity componentA∇φ,n , and (c) an magnification of (b), in the parameterset (A) simulation.
0.3

0.25

0.2

0.15

0.1

0.05
CAS3DMINOS param set (A)00246810t

Fig. 2 Comparison of the growth rates of the Fourier amplitudes
in the parameter set (A) simulation to those of the CAS3Dcomputations.

Proceedings of ITC/ISHW2007

Run553 : total power-1

1010-2-3

1010-410-510-610-7"Run553.spec_total.dat" u 1:2"Run553.spec_total.dat" u 1:3"Run553.spec_total.dat" u 1:410-8-9100204060

80

100120140t

Fig. 3 Time evolutions of the Fourier amplitudes of the normal
velocity component in the parameter set (B) simulation.

(a)
10-2
10-4
10-6

n=0n=1n=2n=3n=4n=5n=6n=7n=8n=9n = 10n = 11

10-8-10

10-12

10-14

10

020

40

60

80

100

120

140

t(b)
10-2
10-4
10-6

n=0
n=1
n=2
n=3
n=4
n=5n=6n=7n=8n=9n = 10n = 11

10-8-10

10

10-12
10-14
0

20

40

60

80

100

120

140

t(c)
10-2
10-4
10-6

n=0
n=1
n=2
n=3
n=4
n=5n=6n=7n=8n=9n = 10n = 11

10-810-1010-1210-14020406080100120140t

Fig. 4 Time evolutions of the Fourier amplitudes of (a)the normal, (b) the binormal and (c) the parallel velocity components in the parameter set (B) simulation.

Fig. 5 Isosurfaces and contours of the pressure on a verticallyelongated poloidal cross-sections at t = 110τA (upper)and 150τA (lower). At the midst of the two times, thecontour lines of the pressure are much more chaotic. After the destruction of the closed contour lines, concentric
profiles of the pressure tend to be recovered.