Proceedings of ITC/ISHW2007

Resistive MHD stability studies for LHD configurations
Luis GARCIA1) , Jacobo VARELA1) and Kiyomasa WATANABE2)
1)

2)

Universidad Carlos III, 28911 Leganés, Madrid, Spain
National Institute for Fusion Science, Toki, Gifu 509-5292, Japan
October 15, 2007

A set of reduced equations is derived without making any assumption on the number of field periods.
The equilibrium used in the equations is the exact three-dimensional equilibrium without any average in thetoroidal angle. Using this set of equations, we study the ideal and resistive MHD properties of different LHD
configurations. The linear results are compared with the ones obtained from the full MHD equations in thepressure-convection limit. The agreement is very good.
Keywords: Resistive MHD, LHD, Reduced equations

1Introduction

Pressure-driven instabilities are a key feature in stellarator stability, since no net toroidal current flows in theplasma, avoiding the appearance of current-driven instabilities. Modes with high poloidal wavenumber m are usually studied using the ballooning formalism, which reduces
the stability problem to finding the eigenvalue of a systemof differential equations to be solved along the magneticlines. The stellarator expansion can be applied to study thestability of modes with low toroidal wavenumber n whosevariation along filed lines is slow compared with the variation of the stellarator terms. A reduced set of equations
expressed in terms of the equilibrium flux coordinates wasderived for those modes [1]. They are formally the sameas the reduced set of MHD equations for a tokamak.
The calculation of global modes using the full threedimensional equilibrium has been based on formulationsof the ideal MHD energy principle in magnetic coordinates [2, 3]. However, these formulations cannot include
the effect of resistivity and are not suitable for nonlinear
calculations. In order to be able of studying the nonlinear evolution, we developed a numerical code (FAR3D),
which solves the full set of resistive MHD equations [4].

2Reduced Equations

For high-aspect ratio configurations with moderate βvalues (of the order of the inverse aspect ratio), we canapply the method employed in Ref. [1] for the derivationof the reduced set of equations without averaging in thetoroidal angle. In this way, we get a reduced set of equations using the exact three-dimensional equilibrium. In thisformulation, we include linear helical couplings between
mode components, which were not included in the formulation developed in Ref. [1].
author’s e-mail: lgarcia@fis.uc3m.es

The main assumptions for the derivation of the set of
reduced equations are high aspect ratio, medium β (of theorder of the inverse aspect ratio ε), small variation of thefields, and small resistivity. With these assumptions, wecan write the velocity and perturbation of the magneticfield as√(1)v = gR0 ∇ζ × ∇Φ, B = R0 ∇ζ × ∇ψ,
where ζ is the toroidal angle, Φ is a stream function proportional to the electrostatic potential, and ψ is the perturbation of the poloidal flux.
The equations, in dimensionless form, are∂ψη= ∇k Φ + Beq
Jζ∂t
S ζ(2)!
√√∂Uβ0 1 ∂ g ∂p ∂ g 1 ∂p= − v · ∇U + 2−∂t∂ρ ρ ∂θ2ε ρ ∂θ ∂ρζ1 ∂  ζ  1 ∂ψ 1 ∂Jeq ∂ψ+ ∇k J ζ −ρJeq+ρ ∂ρρ ∂θ ρ ∂θ ∂ρ(3)d peq 1 ∂Φ
∂p= −v · ∇p +(4)∂tdρ ρ ∂θ √ iζ√ hHere, U = g ∇ × ρm gv , where ρm is the mass density. All lengths are normalized to a generalized minor radius a; the resistivity to η0 (its value at the magnetic axis);
the time to the poloidal Alfvén time τhp = R0 (µ0 ρm )1/2 /B0 ;the magnetic field to B0 (the averaged value at the magnetic axis); and the pressure to its equilibrium value at hemagnetic axis. The Lundquist number S is the ratio of theresistive time τR = a2 µ0 /η0 to the poloidal Alfvén time.
Equilibrium flux coordinates (ρ, θ, ζ) are used. Here, ρis a generalized radial coordinate proportional to the squareroot of the toroidal flux function, and normalized to one at
the edge. The flux coordinates used in the code are those√
described by Boozer [5], and g is the Jacobian of the coordinates transformation. The code uses finite differences

Proceedings of ITC/ISHW2007
1.6

5

1.4

1.2

0.8

1.2
1
0.8

2

0.6

10.6
DI

3curvatureiota

DI (iota=0.5)
DI (iota=1)

14

0.4
0.2
00.4

-0.2

000.2

0.4

ρ

0.6

0.8

1Fig. 1 Equilibrium profiles of the rotational transform and curvature for β0 = 2 × 10−3 . The LHD configuration is described in the text.

-0.4
00.01

0.02

0.03
β00.040.05

0.06

Fig. 2 Values of the ideal Mercier criterion vs. β0 for the same
LHD configuration as in Fig. 1 at the radial position ofthe -ι = 1/2 and -ι = 1 rational surfaces.

in the radial direction and Fourier expansions in the twoangular variables. The numerical scheme if semi-implicitin the linear terms. The nonlinear version uses a two semisteps method to ensure (∆t)2 accuracy.

10-1

0.25
0.2-2100.15

The method is especially suitable for the study of LHD
configurations since they verify the assumptions, and thedominant equilibrium modes are (m = 1, n = 0), and(m = 2, n = 10). We start by studying the linear stability
properties of a sequence of zero net-current fixed boundary equilibria with parameters Rax = 3.6 m, Bq = 100%,
γ = 1.25, and a pressure profile p ∼ 1 − ρ2 . The profilesof the rotational transform and curvature for the case with
β0 = 2 × 10−3 are shown in Fig. 1. The rotational transformprofile changes very slightly for the scan in β. The valueof the ideal Mercier stability criterion DI [6] as a functionof β0 at the position of the ι- = 1/2 and -ι = 1 rational surfaces is shown in Fig. 2. From the Mercier criterion we
expect that ideal modes localized at ι- = 1/2 are stable forβ0 < 6 × 10−3 , and those localized at ι- = 1 are stable forβ0 < 0.02. The results of the linear growth rate and thewidth of the dominant Fourier component (m = 2, n = 1)of the n = 1 mode are shown in Fig. 3 for the ideal case andthe case with S = 107 . In these calculations, we include
equilibrium components with n = 0, 10, 20, and, consequently, dynamic components with n = 1, 9, 11, 19, 21.
The mode is ideally stable for β0 < 8 × 10−3 , a value close
to the one given by the Mercier criterion. The resistive andideal modes are very localized for β0 -values below 0.01,
and they are very similar for β0 -values above 0.01.
The effect of the toroidal and helical couplings canbe seen in Fig. 4, where the dominant component of then = 1 family is shown for three different linear calculations: When only one dynamical component is included(Cylinder), when only components with n = 1 are includedγLinear Results

W

30.1

10-30.05
γ (S = 107)γ (Ideal)

W (S = 107)
W (Ideal)0.03
β00.05

10-4
00.01

0.02

0.04

00.06

Fig. 3 Ideal and resistive (S = 107 ) linear growth rate of then = 1 mode vs. β0 for the same LHD configuration as inFig. 1. The width of the dominant (m = 2, n = 1) Fouriercomponent is also represented(Toroidal), and when components with n = 1, 9, 11 are included (Helical). The equilibrium parameters of this configuration are Rax = 3.6 m, Bq = 100%, γ = 1.25, I p < 0,and correspond to an experimental discharge with localized oscillations at the (m = 2, n = 1) rational surface [7].
The Lundquist number S is 8 × 107 , and the growth rateincreases by more than a factor of 2 from the cylindricalto the full calculation. The dominant radial magnetic fieldcomponents are shown in Fig. 5. The importance of thehelical couplings for the magnetic terms is clear from theFigure.

4

Nonlinear Results

We have followed the nonlinear evolution in the cylindrical and toroidal limit. The value of S in these calculations
is reduced in such a way that the linear growth rate andProceedings of ITC/ISHW2007
1.2

10-6

CylinderToroidalHelical1n= 3n= 4
n= 5

10-8

0.8
0.6

Vρ2Vρ2/1

total
n= 1
n= 2

10-7

10-9

10-10

0.4

10-11

0.2

10-12

0
0.45

0.5

0.55
ρ0.6

10-13

0.65

Fig. 4 Comparison of the dominant component of linear eigenfunctions with no couplings, with toroidal couplings, andwith toroidal and helical couplings. The LHD configuration is described in the text.

05000t/τhp

10000

15000

Fig. 6 Nonlinear evolution of the integral of the mean squared
value of V ρ for the case of Fig. 4.
1 10-4
ToroidalCylinder
Linear-58 10

3
2/1
1/1
3/1
4/11
3/11
5/11
1/9

Bρ

1Vρ26 10-54 10-5
2 10-5000-1-2
00.2

0.4

ρ

0.6

0.8

1

0.2

0.4

ρ

0.6

0.8

1Fig. 7 Comparison of the time-averaged root mean squared
value of V ρ during the stationary phase for the nonlinear
cylindrical and toroidal cases. The linear eigenfunctionis also plotted.

Fig. 5 Dominant radial magnetic field components for the caseof Fig. 4 with all the couplings included.

5width of the dominant component are similar to those of
the linear calculation including helical couplings. The calculation for the cylindrical limit is single-helicity, that is,
only components with m/n = 2 are included. As can beseen in Fig. 6, the nonlinear evolution for the cylindricalcase leads to saturation with bursting activity. The pressure profile flattens around the ι- = 1/2 rational surface. For
the toroidal limit, the components with the same n are linearly coupled. The saturation level increases with respect
to the single-helicity case, and the profile of the root meansquared value of the radial velocity widens with respect tothe linear eigenfunction, as can be seen in Fig. 7. In both
calculations, the (m = 2, n = 1) component dominates thespectrum. The calculation of the nonlinear evolution including linear helical couplings is under way.

Pressure-convection Limit

The natural generalization of the set of reduced equations
would be to write the velocity and magnetic field asv=√ 
g R0 ∇ζ × ∇Φ + ∇θ × ∇ (ρΛ) ,B = R0 ∇ζ × ∇ψ + ∇θ × ∇ (ρχ) ,(5)(6)

where Φ and Λ are velocity stream-functions, and ψ and χ,
the perturbations of the poloidal and toroidal flux, respectively. In this formulation, the incompressibility conditionis approximate(higher order in the reduced equations),√ ∇ · v/ g = 0. It corresponds to the pressure-convection
limit of Ref. [4]. To ascertain the validity of the reduced
set of equations for LHD, we have compared the resultsof linear calculations with those obtained using the fullMHD equations in the pressure convection limit. The results practically do not change, and are consistent with theProceedings of ITC/ISHW2007
0.80.02Φ (Press. conv.)
Φ (Reduced eqs.)

0.60.015

Λ (Press. conv.)

0.01

0.2

0.005

Λ1/1

Φ1/1

0.4

00-0.2-0.005
00.2

0.4

ρ

0.6

0.8

1Fig. 8 Comparison of the dominant Fourier component of thestream-function Φ obtained from the reduced and pressure convection limit equations. The stream-function Λis also plotted using a different scale. The corresponding
LHD configuration is described in the text.

approximations made. This is illustrated in Fig. 8 where
we have plotted the dominant component (m = 1, n = 1)
obtained from both calculations for a configuration withaspect ratio 8.3 and γ = 1.13.

6Summary

We have derived a set of reduced equations without any assumption on the number of field periods. For the calculation of the linear growth rates and nonlinear evolution, weuse Boozer coordinates. The linear results agree very wellwith the ones obtained from the full MHD equations in
the pressure-convection limit. We have studied the nonlinear evolution of the fluctuations in the cylindrical (singlehelicity) and toroidal (only n = 0 equilibrium modes) limits. The saturated level widens with respect to the linear
eigenfunction. The calculation of the nonlinear evolution
including linear helical couplings is under way.
One of the authors (L.G.) thanks the staff at NIFS for
their support during his stay when part of this work wasdone.
[1] L. Garcia, et al., Phys. Fluids B 2, 2162 (1990).
[2] D.V. Anderson, et al., Proceedings of Joint VarennaLausanne International Workshop on Theory of Fusion Plasmas (Compositori, Bologna, 1988) p. 93.
[3] C. Schwab, Phys. Fluids B 5, 3195 (1993).
[4] L. Garcia, Proceedings of the 25th EPS International Conference, Prague, 1998, VOL. 22A, Part II, p. 1757.
[5] A.H. Boozer, Phys. Fluids 25, 520 (1982).
[6] C. Mercier, Nucl. Fusion 1, 47 (1960).
[7] A. Isayama, et al., Plasma Phys. Control. Fusion 48, L45(2006).