- Email: [email protected]
A finite element approach to the computation of the MHD spectrum of straight noncircular plasma equilibria
Computer Physics Communications 11(1976) 313—323 © North-Holland Publishing Company
A FINITE
ELEMENT APPROACH TO THE COMPUTATION OF THE MHD SPECTRUM OF STRAIGHT NONCIRCULAR PLASMA EQUILIBRIA D. BERGER, R. GRUBER and F. TROYON Centre de Recherches en Physique des Plasmas, Ecole Polytechnique Fédérale de Lausanne, Switzerland Received 27 July 1976
A numerical scheme to compute the MHD spectrum of a straight plasma column of arbitrary cross section is described. It is based on a finite element expansion in terms of a carefully chosen set of basis functions. Numerical tests show that kink modes are well reproduced for plasmas of arbitrary low /3 but that the weakly growing internal modes are difficult to obtain in the low-Ø regime.
1. Introduction The stability of a straight cylindrical plasma column of circular cross section in which the equilibrium quantities depend only on one variable (the radius) is fairly well understood. Codes have been written which compute, for a given equilibrium, growth-rates of instabilities and corresponding eigenfunctions. These codes use either a shooting method, solving the Ham— LUst equation [1] directly, or a finite element method [2—4]. The finite element method has the advantage that it gives quickly an overall view of the whole spectrum and thus permits a better qualitative understanding of the dependence of the eigenfrequencies and the eigenmodes on the parameters of the configuration. The main interest is now shifting to axisymmetric toroidal configurations of small aspect ratio with non-circular cross-section. The equilibrium quantities for these plasmas are functions of two variables (radius and angle) so that the stability problem becomes two-dimensional. Stability criteria have been derived, such as the Mercier criterion [5]. The influence of toricity on the internal kink mode [6] and the effect of non-circularity [7] have been studied by expanding the potential energy in powers of the inverse aspect ratio or in powers of the ellipticity. A major disadvantage in all these analytic calculations
is that growth-rates and eigenmodes are obtained only qualitatively. The only possibility to obtain the actual growth-rates and the eigenfunctions is to solve the MHD stability equations numerically. We distinguish between two main approaches: the evolutionary method and the variational method. The first successful codes to solve the ideal MHD stability equations were the time evolutionary codes [8,9]. In these calculations the linearized equations of motion are solved numerically by a finite difference scheme. An initial equilibrium state is perturbed randomly and the asymptotic time behaviour of this initial noise is determined. This technique works well in the high-j3 regime and gives reasonable accuracy for the most unstable kink mode. The time evolution. ary codes are inadequate for treating Tokamak.like plasmas with a low 13. For such plasmas the proximity of the nearly marginally stable modes of the continuum makes it possible to extract the unstable mode in a reasonable time. One also finds that localized modes cannot be obtained with the time evolutionary codes. This is probably due to the choice of the coordinate system. A cartesian system does not fit the magnetic flux surfaces. As a consequence, a diffusion of the fluid elements through the flux surfaces takes place which results in a stabilizing mode coupling. The main advantage of the time evolutionary method is that the inclusion of dissipative effects is straightfor.
314
D. Berget et aL/Computation of the MHD spectrum of a straight plasma column
ward. It also constitutes the first step towards obtaining a fully nonlinear code [10,11]. In the variational method the displacement vector is expanded in terms of a set of basis functions [12] and substituted into the lagrangian of the system. Different choices of basis functions can lead to quite different methods. Kerner et al. [13,14] use a global Fourier expansion in the poloidal angle and a sophisticated modified Bessel-expansion in the radial direction. The Princeton group also tried global functions in the two dimensions, and in the last version of their code [15,16] a combination of a global Fourier expansion in the poloidal angle and a finite element expansion in the “radial” direction is taken. The global expansion in the poloidal angle seems to be well suited to their nonorthogonal, so-called “natural”, coordinate system [17] and to the straight elliptical problem that they solve. This work describes an alternative approach which uses a finite element expansion in both directions. The basis functions are chosen within a class of functions which can satisfy V = 0 over an extended region, without imposing additional non-physical constraints. This condition has been shown to be a necessary condition in the one-dimensional case [4] to represent correctly the unstable part of the spectrum at low /3. It is also a sufficient condition to remove the generally occurring “spectrum pollution” [18]. In the two-dimensional case there exists no analytic work on the spectrum, but from the structure of the 8W, it is clear that this incompressibility condition becomes an exact constraint in the limit /3 = 0, and the elements should then be chosen such that this limit can be reproduced. This choice does not mean that V ~ will vanish, but only that it can vanish. A code has been written along these lines and applied to the Gajewski equilibrium which is a straight elliptical equilibrium. The Princeton group [15] has computed analytically the spectrum for this equilibrium in the long wavelengths low-j3 limit (kink limit), assuming the plasma to be surrounded by a vacuum region without any conducting wall. We find good agreement between our numerical results and the exact results. We then investigate the dependence on ellipticity of the unstable modes in the range 0
elements. We also find that ellipticity destabilizes the long wavelength modes. However, the internal modes which are clustered around nq = 1 have growth-rates two orders of magnitude smaller than the kink. Therefore, even for the circular case, we need a large number of elements to represent these modes. As the ellipticity increases the situation becomes rapidly worse and with an ellipticity of 1 .5 no clear convergence law has emerged with the largest number of elements we could use. The difficulty is traces to a poor representation of the B V operator, inherent to the finite element method. This problem is enhanced as the ellipticity increases. This is due to the stronger angular variation of the eigenmodes and to the singular behaviour of the orthogonal coordinate system. In conclusion, the method described is perfectly adequate to represent the fast growing kinks, making it a good tool for studying the global stability of noncircular configurations. There remains a problem of representation of the slower growing internal modes. •
2. The physical problem We consider a straight plasma column with a noncircular cross section confined by a magnetic field B(x, y) and surrounded by a vacuum region. Choosing the magnetic axis as the z axis we denote the axial cosi-iponent of B as Bz and the projection on the (x,y) plane as B~.p, p and / are the density, the plasma pressure and the axial current density respectively. The flux function i~(x,y) is introduced in the usual way through =
Vqi X Ez.
In equilibrium Bz(~’)and p(j’) are functions of only, and i/i is a solution of the equation 2 ~2 d
i/i
(1) 3x 3y Instead ofx andy we use an orthogonal coordinate system P, x- The surfaces x = const. are orthogonal to the flux surfaces j’ = const. The volume element in this system is ~‘
dr~dxdydzJdb dxdz
-
(2)
The labeling of the x surfaces is fixed by imposing the condition
D. Berger et al/Computation ofthe MHD spectrum of a straight plasma column
(3)
aJ(1i~,x)Iax = 0
The vector ~ is related to the physical displacement by /X\
where the index s refers to the plasma surface. As a result the x coordinate coincides, on the surface, with the nonorthogonal natural poloidal coordinate used by the Princeton group [171,but everywhere else the coordinate systems are different. The equations of evolution of a perturbation of the form
315
~as(
/JB~~~\
Y)=( \~z/
J.
J~
(10)
\~~/BpJ
The variation in eq. (5) has to be performed with the two constraints (ii
X A) 5
~(x) = Re[~/i,X)eit1*~]
=
(11) (4)
A(x) = Re[A(x,y)e~*~]
where ~(x) is the displacement in the plasma and A(x) the potential vector in the vacuum region,are equivalent to the variational principle (5) 6(Wp+Wv_W2WK)=0,
(nXA)~=0,
where s refers to the plasma vacuum interface and c to the conducting shell which surrounds the entire system. Eqs. (5) lead toand an eigenvalue problemkin 2 isand the(11) eigenvalue ~the eigenvector; is which w the axial wave number of the perturbation.
where 3. The regular finite elements method
w~=~jr ffd~Pdx(~-!~-IF[x]I2
In a previous paper [19] we have shown in the onedimensional case (straight circular cross-section plasma column) that a bad choice of elements leads to spectral ++IF[Y] —B 2 pollution [18] and to abad representation of the un5V~I stable part of the spectrum. These difficulties have their origin in the form of the plasma potential energy 2 W~,and they have been solved altogether by an adequate +JB~IF[ZJ_V. jX~ choice of elements [4]. The criterion is that the basis functions for the three components of the displace+~p~ V~I2 2~IXI2), / ment ~ must be such that V = 0 can be satisfied identically over an extended region without imposing an additional unphysical constraint on the motion WV=~7T ffdxdyIVXAI2, [4,20]. We apply the same receipe to the two-dimenvacuum sional case! x 12 1y12 Performing the transformation (10), V~~ becomes WK=~1T ffd~dxpJ(l~—I +IZBpI2)~ V~~=(ax/a~,)+iky÷(az/ax). (12) plasma (8) plasma
JB~
4~
—
‘-
~I7I
and F[fl V~ =~+ikY+~,
In the variational form (5) the component X appears as X, axia~and ax/ax. The other two components Yand Z only appear as Y, ay/ax and Z, az/ax with no ~I’-derivatives. This implies that in eq. (5) the variation of say 8~,belongs to a space U(&~)defined by -
~,
ax K=
—i
4
(9)
U(~7)~H~(&2)C {log(JB~)}
where
2(fl)C2(&2),
(13)
D. Berger et aL/Computation of the MHD spectrum of a straight plasma column
316
af af H~2)” f,~,~—EL2(~2),f(0,x)0
e~ ,
(14)
~,
when
or x~x0
or ~
~JI~JI~
-‘
or x~x1_1
C2(~2)= [f~~iiEL2(~2)}.
(15)
~i~~L::’
when ~
~
and
x1_1~x~x1
‘I’i—v)i_I X/—X/_l
L2(~7)is the space of square integrable functions on the domain ~l = [0 ~ ~ ~ 0 ~ 2ir] such that f(,(i, 2ir) =f(iji, 0). The notation in eq. (13) means that X belongs to J4(~2)and Y, Z to C2(~l).To approximate eq. (5) we choose a finite dimensional subspace V C Udefined by ~
~p—v1i_1X—XJ+1
,
when
lP—~~Lit+j
X~Xj_i
~i_~i+1
~f~~j~l’
‘p—’pi+1 x—x1+1 tI)i—’I,i+i
where V1 is a finite dimensional subspace of H~(~2) containing the “finite elements”; V2 and V3 are two finite dimensional subspaces of C2(&~)which have the following properties: for all X E V1 there exists YE V2 and Z E V3 such that V as expressed in eq. (12), vanishes everywhere. Among all possible subspaces V satisfying this condition, we choose the most simple one: for X we choose linear basis functions in ~Dand x. We then require that each term in eq. (12) be piecewise constant in Li and linear in x. As a result Yand Z are both piecewise constant in ~1i.In the direction Yvaries linearly and Z quadratically. Wex use this subspace to expand both the displacement ~ and its variation ö~. in order to effectively construct the elements of V, we cover the domain ~2with a rectangular mesh (~,, x) of N~,intervals in the i~i direction and N~ intervals in the x direction. An element of V 1, X, can then be expressed in terms of a set of basis functions e~by
and x1~x~x1~1
=
(16)
V=V1XV2XV3,
.iJi~~
‘~Ii___~1i—1 Xj_~Xj+i
xi_xi+i’
0,
when
and
x,_1~x~x1
when
and
x1~x~x1~1
when
or
or X>XNX or x~’x1~1.
(18)
~,
The dimension of V1 is N~ N~. An element of V2 can be analogously expressed in as 2~ terms of basis functions e~~” N~,—1 ,-
—
i=0
i+1/20
j+1/20~Y
Nxl
1
+
Y
2~ 1~112~ e~”
/=1
X~
where
Nx_l
(Xio(e~+e~X)+~x~14}~
(19)
.
~~‘ith 0,
when
or
x~x0or x~x1_1
x—x1_1 when
N~
i+1/2N~
+ey
e~2~
and
x,1~x~x1 (20)
~
(17)
‘~, ~
0,
when
andxx~~1
or X~XNxor
x~xç11
The dimension of V2 is also Nq, ‘Nx. The elements of V3 are higher order in x. An element Z can be expanded in terms ofastwo sets of basic func214~I’2 tions e1~1,‘2j, and e~hI’
D. Berger et al/Computation of the MHD spectrum ofa straight plasma column
317
N~—1~
~jZj+i12o(e~~+e~i’2Nx)
Z
e5
1=0 ~ Nx_l
/
~Z~÷
+
11~1e~1I21) /=1
N~,—1N,<~~l +
2~hh’2
i=0
~i;z~+1121+112e~h/ /=0
(21)
where
/
eYL~~\\I
e~/21=
0, when
or x
2(x—x 1_1)(x—x1_1 /2) 2 (x/—x/_l)
=
when
p1
and
x1_1 ~x
e 51
2(x—xj~
e~2
1)(x—x1+112)
(x1—x,~1)2
when
and
x~x~x,~1
when ~p>~p~÷1 or X>XNx Ot
0,
—‘p
(22)
and 1~1~’2 = e~~2
0, when
‘~p<~p1 orx
Fig. 1. The regular basis functions. The maximum heights of
ex, ey, ez1 and ez2 when
<~<~p~ and
x,~x
are
1.
x
0,
when
Xij.i
and x>x11.1
Yi.l/2j.1 2j.1
Z
(23)
•lj
•l/
• Z~v2j.l/2
The dimension of V3 is 2N~, N,~.The elements for I, Y and Z are shown in fig. 1. Note that e~ 2’ ahas nd a support extending over four mesh cells, e~~ e~11’21are non-zero in two mesh cells, while e~~1’2/+h/2 extends only over one mesh cell. The points, where X = X~ 1,Y = ~ Z = Z.÷112,, Z = Z~1121~.112, are the nodal points. Their locations in a mesh cell are shown in fig. 2. It also explains the indexing of the elements. All that has been said until now concerns exclusively the plasma contribution to 6W. In the vacuum region the variation of A, in eq. (5), sayöA, must be
Zi*I/2 1 Xj
Yi.l/2j
Xj.ij
—
Fig. 2. Position of the nodal points in mesh cell.
such that V X ÔÁ belongs to C2(&2) and that the two boundary conditions (11), with &~replacing~, be satisfied. There is no natural coordinate system in vacuum and each configuration has to be considered independently. But since A does not appear in the kinetic energy term, for a given configuration, the variation on A can be done once and for all, considering
D. Berger et a!./Computation of the MHD spectrum of a straight plasma column
318
and ö~as given functions. The 6 Wv in eq. (5) can then be rewritten as
The surfaces of constant magnetic flux ~~ti(x,y)are similar ellipses
~irffdx dx’ Gk(x, x?)(FE6Xs])*Fl[X~] (24) where Gk(x, x’) is a real symmetric function. It can
~i(x,y)
~~{x2 +~(y2/e2)~,
with ~ti~=
e2//2(e2 +
=
,
also be considered as the variation of a potential energy W,,,
=
dx’ Gk(x, ~‘)~F[X
~-irffdx
)*FI 5] ~X~] ,
(25)
(27)
1). The poloidal field B~becomes (28)
B~~x,y) = 2~S\/~2 + ~2/~4)~
The location of the conducting shell is given by the 2 + (y/C)3)2 = 1 (29) ~x/g~) where ,
which replaces expression (7) and the boundary conditions The new potential nofunction longer contains the(11). potential vector A. Theenergy Green’s Gk(x, x’) contains all the needed information on the vacuum region and on the location of the conducting shell. It does not depend on the equilibrium magnetic field or on any other equilibrium parameter except the geometry of the vacuum region. The discretization of 6W is made by expanding X~ and 6X 5 in terms of hat functions in x. The function Gk(x, x’) can only be calculated analytically in simple cases. In general it will have to be computed numerically. The discretized variational problem reduces to a matrix eigenvalue problem of the form .
AE~—w2B~0
26~ ‘-
-‘
in which the matrices A and B have a particular block structure. A program which solves this problem, and which utilizes the block structure of the matrices, has been developed [21] It is an essential part of the code. Without it, we could never have run the cases presented below on our computer (CDC 6500).
.s~ =
—
~2 +
[(1+ e2)2 +4e2(r~
—
1)]h/2}1~’2
,-
w’
in which r~,is the ratio of the cross sectional area of the shell to the cross-sectional area of the plasma. The shape of the shell is chosen only for its conviencence in the calculation of the vacuum contribution. The orthogonal coordinate x is defined by .
(y/ea
sin x)
=
x/(a cos x)
(30)
With this choice, the plasma surface corresponds to
x = cos x~y = ea sin x showing that x coincides there with the angular coordinate chosen by Dewar et al. [15]. For all other values of /i * i~i~ they naturally differ. The jacobian J is then given in semi-implicit form by e2
.
.
4. Application to an elliptical equilibrium
{1
=
1 2xy e 2 cos 2 x + sin 2 x / sin X ~ 2 ~4 2 .
+
(31)
On the plasma surfacef5 = (~2+ l)/(e/). Green’s function G(~,x’) which characterizes the vacuum contribution, eq. (25), is calculated numerically by expanding it in terms of finite elements e~(x):
4.1. The long wavelength, low-current regime G(~,x’) = ~G11e1(~)e1(~’) To test the code we choose an equilibrium proposed by Gajewski [221 which has an elliptical cross-section and a constant current density /. The axial magnetic field B~and the mass density p also are taken to be homogeneous. The elliptical plasma surface is defined by the minor radius of the plasma a, and the excentricity of the ellipse, e = b/a where b is the major plasma radius. We assume a vacuum region around the plasma and a confocal conducting wall enclosing the elliptical plasma and the vacuum region.
.
(32)
“.‘
The element e1(x) is chosen as a hat function centred at x = The calculation of G.1 is described in the Appendix. For e = 1 the equilibrium reduces to the one-dimensional code [23]. For e ~ 1 the spectrum has been determined analytically by Dewar et al. [15] in the long wavelength, low /3-limit, characterized by B~ /Bz ka ~ 1 (so-called kink ordering) and assuming r~= oo• The modes can be divided into classes which ~.
D. Berger et al./Computation of the MHD spectrum of a straight plasma column
are the analogue of the separation in azimuthal numbers m in the circular case. Kink modes are the only unstable modes, the internal modes remaining stable at this order. 2, Thewhich square of the growth-rate of the becomes the m = 1 kink in the kink mode, r’ circular limit is given by F2(q)c~.,~/A —A2 1 (‘~JA1+A2—~..JA1—A2),(33) where
319
The crosses and the dots are the numerical results obtained with rw = 20. The small discrepancy in the circular case is due to the presence of the conducting shell. For rw there is always stability for very long wavelengths (q very small). This stabilizing effect seems to become stronger for the same value of r~as the ellipticity increases. We have verified by taking successively r~= 10, 20 and 40 keeping q = 0.1 fixed, that the discrepancy decreases as r~.Extrapolating to rw = gives the point shown as a square n in fig. 3, which does indeed lie on the analytical curve. The mode is correctly represented in the circular case and has the right behaviour as described in Dewar et a!. [15] for the elliptical case although we did not make a detailed comparison in this last case. We found that as the ellipticity increases, the number of mesh points required increases also. We believe this is due to the pathological behaviour of the orthogonal coordinate system near the axis. In conclusion, we find that the code reproduces correctly the kink modes in the low beta limit. It should be pointed out that a choice of the finite elements which would violate what we have called the divergence condition, could never do this. As the current/ is lowered, the number of mesh points has to increase as an inverse power off to reproduce the kink. 00
=
__________
p(l
+
A
e2)2a2
2,
2
A1
q=
2~i +q2) (1 + 2 k(l +e2~~
(34)
ef
l—q
Note, that the growth-rate scales asj2, for a given safety factor q. The mode is only unstable in the range 0
4.2. The
strong
current regime
As a second test we choose the same equilibrium, but we take a larger current / = 0.4 and a conducting shell close to the plasma, r~ 2. We look for the unstable part of the spectrum in the range 0
q Fig. 3. Comparison between the numerical and analytical values of the growth-rates of the “m = 1” external kink for = 1 and e = 1.5. The solid lines are the analytic results for r~= =. The X and 0 represent the numerical results obtained with rw = 20 for e = 1 and a = 1.5, respectively. The 0 represents the result of an extrapolation to rw results at r~= 10, 20 and 40.
=
from the
320
D. Berger et al./Computation of the MHD spectrum of a straight plasma column
r2a -r
0 .25
.5
75
9
q
Fig. 4. Unstable part of the spectrum in the range 0
~-
1,
in the strong current regime, with the conducting shell at rw = 2.0. The solid lines represent the six most unstable modes for the circular case. The a represent the numerical results for an ellipticity of a = 1.5.
axis. Note that the kink itself becomes also internal near q = 1, the change in behaviour of the mode coinciding with the kink in the curve visible in the figure. There are in fact an infinite number of unstable modes accumulating at q = 0, F2 = 0. All these features are not visible in fig. 3 because of the scale used and because of the smallness of the current. The dotted circles represent the numerical results obtained with the 2-D code for an elongation ~ = 1.5, keeping N~= 12 fixed. The points are obtained by an extrapolation to N~= from N~= 12, 16 and 20 intervals, assuming a quadratic convergence. We see that most of the dependence of the growth-rates on the ellipticity is eliminated by normalizing F2 to w~.At long wavelengths the ellipticity destabilizes, while the conducting wall stabilizes. The logarithmic scale tends to attenuate the difference between the elliptical and circular case. Near q = 1 we only find two unstable modes and the results suggest a stabilizing effect of the ellipticity. This contradicts the conclusions of Lavel [7] which shows, at least in a
particular case, that the internal modes are strongly destablized by the ellipticity. To understand this apparent discrepancy we do a full convergence study of the growth-rate of the second mode, which shows already a strong stabilization for a fixed value of q = 1 and for two different elongations = i and e = 1.5. The results are shown in fig. 5 in which the square of the growth-rate is plotted versus N~2for a fixed N 11 = 12. In the circular case the convergence is indeed quadratic above N~ 20, but in the elliptic case no clear convergence law has emerged yet with the largest number of cells we can use, namely N~= 48. From the behaviour of the curve we see that a quadratic extrapolation underestimates the final value which will most probably be higher than the circular value. This shows that we should not trust the results for the slow growing modes. Even in the circular case, which looks good in fig. 5, there are difficulties in representing internal modes. For example, with the same current! = 0.4, no internal modes are found near q = 2. Also for smaller values of the current even the modes near q = 1 become impossible to find. The ellipticity only magnified the difficulty. The explanation of this difficulty is found by looking at the expression for Wi,, eq. (6). Internal modes are characterized by FX FY FZ 0. With polynomial finite elements F cannot vanish identically, resulting in a stabilizing effect analogous to the effect of a bad
-r Wa
00
Fig. 5. Plot of the square of 2,kthe eeping growth-rate N of an internal mode as a function of N~ 11 fixed. The • represents the value calculated with a 1-D code with the same number of elements.
D. Berger at al./Computation of the MHD spectrum of a straight plasma column
representation of the divergence, but smaller by a factorB~/B~. This explains that at the order of the kink we only need to satisfy the divergence condition while at the order of the internal modes we need to satisfy
321
both.
greatfully acknowledged. F. Debonneville, R. Schreiber and the staff of the Computing Center have been a great help in preparing and running the code. This work is supported in part by the Swiss National Science Foundation.
5. Conclusions
Appendix
From the examples presented above we can draw some interesting conclusions. With its judicious choice of regular basis functions, the code is capable of well reproducing kink modes for arbitrary /3, with a reasonable number of mesh points. This result has been obtained with a shearless equilibrium, but from our experience with the 1—D case, shear only affects the distribution of mesh points in the ~,D direction. In all our calculations we never observed spectral pollution [18]. For elongations larger than two the orthogonal coordinate system becomes very singular near the axis and it is difficult to represent accurately the kink modes. This could be taken care of by a change of coordinate system. The toroidal effects will not bring any new problems. The V ~ condition [seeeq. (12)] must be replaced by a similar condition on thefrom quantity 2)[24] where r denotes the distance the V (~/raxis. toroidal However, the weakly unstable internal modes are strongly stabilized. This comes from a bad representation of the operator F in eq. (6) which has to vanish identically on a rational surface. This condition can never be fulfilled with the described method, since the two terms in the F-operator have different functional dependences. This gives rise to a similar stabilizing effect as recognized in the one-dimensional code when we did not filfull the divergence condition [19]. Whereas the divergence condition affects the kink mode, a bad representation of the operator F only stabilizes the internal modes and does not influence the kink. This problem is specific to the two-dimensional case, since in the one-dimensional case Fis a number which can vanish exactly.
We limit ourselves to the case k ~ 0. The Green’s function Gk(x, x’) can be computed by reformulating the vacuum problem in terms of a scalar potential. We write =
V~(x,y,z),
(Al) ~uRe[~(x,y)e”~] where t~i(x,y)satisfies the Helmholz equation
a2
a2 ax2 ay2
(A2)
with the boundary conditions dtP/dflIcO, (A3) d~/dnI 5= F[X5]
/(.J5B~),
in which F is the operator defined in eq. (9), which reduces on the surface to
J5ax + ikB
F = -~-
-~-
(A4)
The vacuum contribution to the potential energy W~, eq. (7), can be written as an integral over the plasma— vacuum interface 1i 2ir Wv= 7Tf J5B~dx(~)~.
(A5)
We replace the x,y coordinate system by the elliptical system ~i, x defined by x
=
a
sh p cos x~
Acknowledgement
y
=
a ~.,/~TTi
ch p sin
Discussions with Professor J.L. Descloux and Dr. J. Rappaz on the mathematical implications of our choice of elements and on spectral pollution are
The plasma surface and the shell are surfaces of constant p, called ji~ and ~ respectively. On the surface, x as defined by eq. (A6) coincides with the x introduced inside the plasma.
(A6)
x.
D. Berger et al/Computation of the MHD spectrum of a straight plasma column
322
Eqs. (A2) and (A3) become
W~=
ap2 ax2 2
(A7)
acb/apI
aØ/apIc=0, a2(2
5=F[x5]
—
where A A 5 =J5B~. These equations are equivalent to the weak variational principle: find a function ~(p,x) periodic in x such that
fdxfd~
(p ~p
+~ ~
+ k2A2~6ø}
6~(p 5,x)O,
—fdxF[X5]
‘A8’ ‘
‘
1 (vacuum). for Eq. any(A8) 6cb EisHsolved by a finite element expansion of ~ and 60: =
~0Oi, x)
(A9)
~60 1e1(p,
x)
The same basis is used for the two variables. Substituting in eq. (A8) we obtain the inhomogeneous system of equations
(AlO) where S. is a vector given by S1
=
f
dxF[X51
e1(p5,
x).
(Al l~
[8]
0
0 is obtained in terms of S by solving eqs. (A10). We are interested only in the value of ~ x) which can be written
0(~5,x)= ~
G11S1e1(p5,x),
(A12)
1E5,JES
where the sum is carried only over elements which contribute on the surface. G~,is a real symmetric matrix. Replacing in eq. (AS)
G~
f
d~A5(F[X5]
)* e1(p5, x)
2ir
X
1) (sh2 p + cos2x). Note that
~
f
d~’A~F’[X~] ej(p5,
x’).
(Al 3)
This expression is identical to eq. (24) provided we identify Gk(x, x’) with the expression (32). To carry the calculation we cover the domain 0 ~ ~ 2ir, ~ ~ ~ ~ with a rectangular mesh (p1, x) and choose as elements pyramids centred on the mesh points with an hexagonal base (see ref. [12] for details). The elements e~(p5,x) are then indeed hat functions in x~ which achieves the identification with eq. (32). In practice we found that the mesh should be concentrated near the plasma surface in order to obtain the best results. The numerical approximation of Gk(x, x’), which is a singular integrable function, by regular functions, leads to a value of W~which is too small for large m, thus underestimating the stabilizing vacuum contribution.
References [1] K. Ham and R. Lust, Z. Naturforschung 13a (1958) 936. [2] T. Takeda, Y. Shimomura, M. Ohta and M. Yoshikawa, Phys. Fluids 15 (1972) 2193. [31 T.J.M. Boyd, G.A. Gardner and L.R.T. Gardner, Nucl. Fusion 13 (1973) 229. [41 K. Appert, D. Berger, R. Gruber and J. Rappaz, J. Comput. Phys. 18 (1975) 284. [51 C. Mercier, Nuel. Fusion 1 (1960) 47. [6] M.N. Bussac, D. Edery, G. Laval, R. Pellat and J.L. Soulé, 7th Europ. Conf. on Contr. Fusion and Plasma Phys., Vol.1 (1975) p. 100. [7] G. Laval, Phys. Rev. Letters 34 (1974) 1316. A. Sykes and J.A. Wesson, Nucl. Fusion 14 (1974) 645. [9] G. Bateman, W. Schneider and W. Grossmann, Nucl. Fusion 14(1974)669. [10] J. Wooten, H.R. Hicks, G. Bateman and R.A. Dory, ORNL TM 4784 (1974). [Ill A. Sykes and J.A. Wesson, 7th Europ. Conf. on Contr. Fusion and Plasma Phys., Vol.1(1975), p. 111. [12] O.C. Zienkiewicz, The Finite Element Method in Structural and Continuum Mechanics (McGraw-Hill, New York, 1971). [13] W. Kerner and H. Tasso, IAEA-CN-33/A13-l (Tokyo, 1975), p. 475. [14] W. Kerner, MAlT-i 229 (1976).
D. Berger et aL/Computation of the MHD spectrum of a straight plasma column [151 R.L. Dewar, R.C. Grimm, J.L. Johnson, E.A. Frieman, J.M. Greene and P.H. Rutherford, Phys. Fluids 17
(1974) 930. [16] R.L. Dewar, J.M. Greene, R.C. Grimm and J.L. Johnson, J. Comput. Phys. 17 (1975). [17] R.C. Grimm, J.M. Greene and J.L. Johnson, in: Methods of Computational Physics, Vol. 16, eds. B. Alder, S. Fernbach, J. Killeen and M. Rotenberg (Academic Press, New York, 1976, ch. 4. [18] J. Rappaz, PhD Thesis, EPF-Lausanne (1976).
323
[19] K. Appert, D. Berger, R. Gruber, F. Troyon and J. Rappaz, ZAMP 25 (1974) 229. [20J R. Gruber, PhD Thesis, EPF-Lausanne (1976). [21] R. Gruber, Contributed at the 2nd Europ. Conf. on Comp. Phys., Garching, 1976, Paper F2. [22] R. Gajewski, Phys. Fluids 15 (1972) 70. [23] K. Appert, D. Berger, R. Gruber, F. Troyon and K.V. Roberts, Comput. Phys. Commun. 10(1975) 11. [24] J.L. Soulé, private communication.
A FINITE
ELEMENT APPROACH TO THE COMPUTATION OF THE MHD SPECTRUM OF STRAIGHT NONCIRCULAR PLASMA EQUILIBRIA D. BERGER, R. GRUBER and F. TROYON Centre de Recherches en Physique des Plasmas, Ecole Polytechnique Fédérale de Lausanne, Switzerland Received 27 July 1976
A numerical scheme to compute the MHD spectrum of a straight plasma column of arbitrary cross section is described. It is based on a finite element expansion in terms of a carefully chosen set of basis functions. Numerical tests show that kink modes are well reproduced for plasmas of arbitrary low /3 but that the weakly growing internal modes are difficult to obtain in the low-Ø regime.
1. Introduction The stability of a straight cylindrical plasma column of circular cross section in which the equilibrium quantities depend only on one variable (the radius) is fairly well understood. Codes have been written which compute, for a given equilibrium, growth-rates of instabilities and corresponding eigenfunctions. These codes use either a shooting method, solving the Ham— LUst equation [1] directly, or a finite element method [2—4]. The finite element method has the advantage that it gives quickly an overall view of the whole spectrum and thus permits a better qualitative understanding of the dependence of the eigenfrequencies and the eigenmodes on the parameters of the configuration. The main interest is now shifting to axisymmetric toroidal configurations of small aspect ratio with non-circular cross-section. The equilibrium quantities for these plasmas are functions of two variables (radius and angle) so that the stability problem becomes two-dimensional. Stability criteria have been derived, such as the Mercier criterion [5]. The influence of toricity on the internal kink mode [6] and the effect of non-circularity [7] have been studied by expanding the potential energy in powers of the inverse aspect ratio or in powers of the ellipticity. A major disadvantage in all these analytic calculations
is that growth-rates and eigenmodes are obtained only qualitatively. The only possibility to obtain the actual growth-rates and the eigenfunctions is to solve the MHD stability equations numerically. We distinguish between two main approaches: the evolutionary method and the variational method. The first successful codes to solve the ideal MHD stability equations were the time evolutionary codes [8,9]. In these calculations the linearized equations of motion are solved numerically by a finite difference scheme. An initial equilibrium state is perturbed randomly and the asymptotic time behaviour of this initial noise is determined. This technique works well in the high-j3 regime and gives reasonable accuracy for the most unstable kink mode. The time evolution. ary codes are inadequate for treating Tokamak.like plasmas with a low 13. For such plasmas the proximity of the nearly marginally stable modes of the continuum makes it possible to extract the unstable mode in a reasonable time. One also finds that localized modes cannot be obtained with the time evolutionary codes. This is probably due to the choice of the coordinate system. A cartesian system does not fit the magnetic flux surfaces. As a consequence, a diffusion of the fluid elements through the flux surfaces takes place which results in a stabilizing mode coupling. The main advantage of the time evolutionary method is that the inclusion of dissipative effects is straightfor.
314
D. Berget et aL/Computation of the MHD spectrum of a straight plasma column
ward. It also constitutes the first step towards obtaining a fully nonlinear code [10,11]. In the variational method the displacement vector is expanded in terms of a set of basis functions [12] and substituted into the lagrangian of the system. Different choices of basis functions can lead to quite different methods. Kerner et al. [13,14] use a global Fourier expansion in the poloidal angle and a sophisticated modified Bessel-expansion in the radial direction. The Princeton group also tried global functions in the two dimensions, and in the last version of their code [15,16] a combination of a global Fourier expansion in the poloidal angle and a finite element expansion in the “radial” direction is taken. The global expansion in the poloidal angle seems to be well suited to their nonorthogonal, so-called “natural”, coordinate system [17] and to the straight elliptical problem that they solve. This work describes an alternative approach which uses a finite element expansion in both directions. The basis functions are chosen within a class of functions which can satisfy V = 0 over an extended region, without imposing additional non-physical constraints. This condition has been shown to be a necessary condition in the one-dimensional case [4] to represent correctly the unstable part of the spectrum at low /3. It is also a sufficient condition to remove the generally occurring “spectrum pollution” [18]. In the two-dimensional case there exists no analytic work on the spectrum, but from the structure of the 8W, it is clear that this incompressibility condition becomes an exact constraint in the limit /3 = 0, and the elements should then be chosen such that this limit can be reproduced. This choice does not mean that V ~ will vanish, but only that it can vanish. A code has been written along these lines and applied to the Gajewski equilibrium which is a straight elliptical equilibrium. The Princeton group [15] has computed analytically the spectrum for this equilibrium in the long wavelengths low-j3 limit (kink limit), assuming the plasma to be surrounded by a vacuum region without any conducting wall. We find good agreement between our numerical results and the exact results. We then investigate the dependence on ellipticity of the unstable modes in the range 0
elements. We also find that ellipticity destabilizes the long wavelength modes. However, the internal modes which are clustered around nq = 1 have growth-rates two orders of magnitude smaller than the kink. Therefore, even for the circular case, we need a large number of elements to represent these modes. As the ellipticity increases the situation becomes rapidly worse and with an ellipticity of 1 .5 no clear convergence law has emerged with the largest number of elements we could use. The difficulty is traces to a poor representation of the B V operator, inherent to the finite element method. This problem is enhanced as the ellipticity increases. This is due to the stronger angular variation of the eigenmodes and to the singular behaviour of the orthogonal coordinate system. In conclusion, the method described is perfectly adequate to represent the fast growing kinks, making it a good tool for studying the global stability of noncircular configurations. There remains a problem of representation of the slower growing internal modes. •
2. The physical problem We consider a straight plasma column with a noncircular cross section confined by a magnetic field B(x, y) and surrounded by a vacuum region. Choosing the magnetic axis as the z axis we denote the axial cosi-iponent of B as Bz and the projection on the (x,y) plane as B~.p, p and / are the density, the plasma pressure and the axial current density respectively. The flux function i~(x,y) is introduced in the usual way through =
Vqi X Ez.
In equilibrium Bz(~’)and p(j’) are functions of only, and i/i is a solution of the equation 2 ~2 d
i/i
(1) 3x 3y Instead ofx andy we use an orthogonal coordinate system P, x- The surfaces x = const. are orthogonal to the flux surfaces j’ = const. The volume element in this system is ~‘
dr~dxdydzJdb dxdz
-
(2)
The labeling of the x surfaces is fixed by imposing the condition
D. Berger et al/Computation ofthe MHD spectrum of a straight plasma column
(3)
aJ(1i~,x)Iax = 0
The vector ~ is related to the physical displacement by /X\
where the index s refers to the plasma surface. As a result the x coordinate coincides, on the surface, with the nonorthogonal natural poloidal coordinate used by the Princeton group [171,but everywhere else the coordinate systems are different. The equations of evolution of a perturbation of the form
315
~as(
/JB~~~\
Y)=( \~z/
J.
J~
(10)
\~~/BpJ
The variation in eq. (5) has to be performed with the two constraints (ii
X A) 5
~(x) = Re[~/i,X)eit1*~]
=
(11) (4)
A(x) = Re[A(x,y)e~*~]
where ~(x) is the displacement in the plasma and A(x) the potential vector in the vacuum region,are equivalent to the variational principle (5) 6(Wp+Wv_W2WK)=0,
(nXA)~=0,
where s refers to the plasma vacuum interface and c to the conducting shell which surrounds the entire system. Eqs. (5) lead toand an eigenvalue problemkin 2 isand the(11) eigenvalue ~the eigenvector; is which w the axial wave number of the perturbation.
where 3. The regular finite elements method
w~=~jr ffd~Pdx(~-!~-IF[x]I2
In a previous paper [19] we have shown in the onedimensional case (straight circular cross-section plasma column) that a bad choice of elements leads to spectral ++IF[Y] —B 2 pollution [18] and to abad representation of the un5V~I stable part of the spectrum. These difficulties have their origin in the form of the plasma potential energy 2 W~,and they have been solved altogether by an adequate +JB~IF[ZJ_V. jX~ choice of elements [4]. The criterion is that the basis functions for the three components of the displace+~p~ V~I2 2~IXI2), / ment ~ must be such that V = 0 can be satisfied identically over an extended region without imposing an additional unphysical constraint on the motion WV=~7T ffdxdyIVXAI2, [4,20]. We apply the same receipe to the two-dimenvacuum sional case! x 12 1y12 Performing the transformation (10), V~~ becomes WK=~1T ffd~dxpJ(l~—I +IZBpI2)~ V~~=(ax/a~,)+iky÷(az/ax). (12) plasma (8) plasma
JB~
4~
—
‘-
~I7I
and F[fl V~ =~+ikY+~,
In the variational form (5) the component X appears as X, axia~and ax/ax. The other two components Yand Z only appear as Y, ay/ax and Z, az/ax with no ~I’-derivatives. This implies that in eq. (5) the variation of say 8~,belongs to a space U(&~)defined by -
~,
ax K=
—i
4
(9)
U(~7)~H~(&2)C {log(JB~)}
where
2(fl)C2(&2),
(13)
D. Berger et aL/Computation of the MHD spectrum of a straight plasma column
316
af af H~2)” f,~,~—EL2(~2),f(0,x)0
e~ ,
(14)
~,
when
or x~x0
or ~
~JI~JI~
-‘
or x~x1_1
C2(~2)= [f~~iiEL2(~2)}.
(15)
~i~~L::’
when ~
~
and
x1_1~x~x1
‘I’i—v)i_I X/—X/_l
L2(~7)is the space of square integrable functions on the domain ~l = [0 ~ ~ ~ 0 ~ 2ir] such that f(,(i, 2ir) =f(iji, 0). The notation in eq. (13) means that X belongs to J4(~2)and Y, Z to C2(~l).To approximate eq. (5) we choose a finite dimensional subspace V C Udefined by ~
~p—v1i_1X—XJ+1
,
when
lP—~~Lit+j
X~Xj_i
~i_~i+1
~f~~j~l’
‘p—’pi+1 x—x1+1 tI)i—’I,i+i
where V1 is a finite dimensional subspace of H~(~2) containing the “finite elements”; V2 and V3 are two finite dimensional subspaces of C2(&~)which have the following properties: for all X E V1 there exists YE V2 and Z E V3 such that V as expressed in eq. (12), vanishes everywhere. Among all possible subspaces V satisfying this condition, we choose the most simple one: for X we choose linear basis functions in ~Dand x. We then require that each term in eq. (12) be piecewise constant in Li and linear in x. As a result Yand Z are both piecewise constant in ~1i.In the direction Yvaries linearly and Z quadratically. Wex use this subspace to expand both the displacement ~ and its variation ö~. in order to effectively construct the elements of V, we cover the domain ~2with a rectangular mesh (~,, x) of N~,intervals in the i~i direction and N~ intervals in the x direction. An element of V 1, X, can then be expressed in terms of a set of basis functions e~by
and x1~x~x1~1
=
(16)
V=V1XV2XV3,
.iJi~~
‘~Ii___~1i—1 Xj_~Xj+i
xi_xi+i’
0,
when
and
x,_1~x~x1
when
and
x1~x~x1~1
when
or
or X>XNX or x~’x1~1.
(18)
~,
The dimension of V1 is N~ N~. An element of V2 can be analogously expressed in as 2~ terms of basis functions e~~” N~,—1 ,-
—
i=0
i+1/20
j+1/20~Y
Nxl
1
+
Y
2~ 1~112~ e~”
/=1
X~
where
Nx_l
(Xio(e~+e~X)+~x~14}~
(19)
.
~~‘ith 0,
when
or
x~x0or x~x1_1
x—x1_1 when
N~
i+1/2N~
+ey
e~2~
and
x,1~x~x1 (20)
~
(17)
‘~, ~
0,
when
andxx~~1
or X~XNxor
x~xç11
The dimension of V2 is also Nq, ‘Nx. The elements of V3 are higher order in x. An element Z can be expanded in terms ofastwo sets of basic func214~I’2 tions e1~1,‘2j, and e~hI’
D. Berger et al/Computation of the MHD spectrum ofa straight plasma column
317
N~—1~
~jZj+i12o(e~~+e~i’2Nx)
Z
e5
1=0 ~ Nx_l
/
~Z~÷
+
11~1e~1I21) /=1
N~,—1N,<~~l +
2~hh’2
i=0
~i;z~+1121+112e~h/ /=0
(21)
where
/
eYL~~\\I
e~/21=
0, when
or x
2(x—x 1_1)(x—x1_1 /2) 2 (x/—x/_l)
=
when
p1
and
x1_1 ~x
e 51
2(x—xj~
e~2
1)(x—x1+112)
(x1—x,~1)2
when
and
x~x~x,~1
when ~p>~p~÷1 or X>XNx Ot
0,
—‘p
(22)
and 1~1~’2 = e~~2
0, when
‘~p<~p1 orx
Fig. 1. The regular basis functions. The maximum heights of
ex, ey, ez1 and ez2 when
<~<~p~ and
x,~x
are
1.
x
0,
when
Xij.i
and x>x11.1
Yi.l/2j.1 2j.1
Z
(23)
•lj
•l/
• Z~v2j.l/2
The dimension of V3 is 2N~, N,~.The elements for I, Y and Z are shown in fig. 1. Note that e~ 2’ ahas nd a support extending over four mesh cells, e~~ e~11’21are non-zero in two mesh cells, while e~~1’2/+h/2 extends only over one mesh cell. The points, where X = X~ 1,Y = ~ Z = Z.÷112,, Z = Z~1121~.112, are the nodal points. Their locations in a mesh cell are shown in fig. 2. It also explains the indexing of the elements. All that has been said until now concerns exclusively the plasma contribution to 6W. In the vacuum region the variation of A, in eq. (5), sayöA, must be
Zi*I/2 1 Xj
Yi.l/2j
Xj.ij
—
Fig. 2. Position of the nodal points in mesh cell.
such that V X ÔÁ belongs to C2(&2) and that the two boundary conditions (11), with &~replacing~, be satisfied. There is no natural coordinate system in vacuum and each configuration has to be considered independently. But since A does not appear in the kinetic energy term, for a given configuration, the variation on A can be done once and for all, considering
D. Berger et a!./Computation of the MHD spectrum of a straight plasma column
318
and ö~as given functions. The 6 Wv in eq. (5) can then be rewritten as
The surfaces of constant magnetic flux ~~ti(x,y)are similar ellipses
~irffdx dx’ Gk(x, x?)(FE6Xs])*Fl[X~] (24) where Gk(x, x’) is a real symmetric function. It can
~i(x,y)
~~{x2 +~(y2/e2)~,
with ~ti~=
e2//2(e2 +
=
,
also be considered as the variation of a potential energy W,,,
=
dx’ Gk(x, ~‘)~F[X
~-irffdx
)*FI 5] ~X~] ,
(25)
(27)
1). The poloidal field B~becomes (28)
B~~x,y) = 2~S\/~2 + ~2/~4)~
The location of the conducting shell is given by the 2 + (y/C)3)2 = 1 (29) ~x/g~) where ,
which replaces expression (7) and the boundary conditions The new potential nofunction longer contains the(11). potential vector A. Theenergy Green’s Gk(x, x’) contains all the needed information on the vacuum region and on the location of the conducting shell. It does not depend on the equilibrium magnetic field or on any other equilibrium parameter except the geometry of the vacuum region. The discretization of 6W is made by expanding X~ and 6X 5 in terms of hat functions in x. The function Gk(x, x’) can only be calculated analytically in simple cases. In general it will have to be computed numerically. The discretized variational problem reduces to a matrix eigenvalue problem of the form .
AE~—w2B~0
26~ ‘-
-‘
in which the matrices A and B have a particular block structure. A program which solves this problem, and which utilizes the block structure of the matrices, has been developed [21] It is an essential part of the code. Without it, we could never have run the cases presented below on our computer (CDC 6500).
.s~ =
—
~2 +
[(1+ e2)2 +4e2(r~
—
1)]h/2}1~’2
,-
w’
in which r~,is the ratio of the cross sectional area of the shell to the cross-sectional area of the plasma. The shape of the shell is chosen only for its conviencence in the calculation of the vacuum contribution. The orthogonal coordinate x is defined by .
(y/ea
sin x)
=
x/(a cos x)
(30)
With this choice, the plasma surface corresponds to
x = cos x~y = ea sin x showing that x coincides there with the angular coordinate chosen by Dewar et al. [15]. For all other values of /i * i~i~ they naturally differ. The jacobian J is then given in semi-implicit form by e2
.
.
4. Application to an elliptical equilibrium
{1
=
1 2xy e 2 cos 2 x + sin 2 x / sin X ~ 2 ~4 2 .
+
(31)
On the plasma surfacef5 = (~2+ l)/(e/). Green’s function G(~,x’) which characterizes the vacuum contribution, eq. (25), is calculated numerically by expanding it in terms of finite elements e~(x):
4.1. The long wavelength, low-current regime G(~,x’) = ~G11e1(~)e1(~’) To test the code we choose an equilibrium proposed by Gajewski [221 which has an elliptical cross-section and a constant current density /. The axial magnetic field B~and the mass density p also are taken to be homogeneous. The elliptical plasma surface is defined by the minor radius of the plasma a, and the excentricity of the ellipse, e = b/a where b is the major plasma radius. We assume a vacuum region around the plasma and a confocal conducting wall enclosing the elliptical plasma and the vacuum region.
.
(32)
“.‘
The element e1(x) is chosen as a hat function centred at x = The calculation of G.1 is described in the Appendix. For e = 1 the equilibrium reduces to the one-dimensional code [23]. For e ~ 1 the spectrum has been determined analytically by Dewar et al. [15] in the long wavelength, low /3-limit, characterized by B~ /Bz ka ~ 1 (so-called kink ordering) and assuming r~= oo• The modes can be divided into classes which ~.
D. Berger et al./Computation of the MHD spectrum of a straight plasma column
are the analogue of the separation in azimuthal numbers m in the circular case. Kink modes are the only unstable modes, the internal modes remaining stable at this order. 2, Thewhich square of the growth-rate of the becomes the m = 1 kink in the kink mode, r’ circular limit is given by F2(q)c~.,~/A —A2 1 (‘~JA1+A2—~..JA1—A2),(33) where
319
The crosses and the dots are the numerical results obtained with rw = 20. The small discrepancy in the circular case is due to the presence of the conducting shell. For rw there is always stability for very long wavelengths (q very small). This stabilizing effect seems to become stronger for the same value of r~as the ellipticity increases. We have verified by taking successively r~= 10, 20 and 40 keeping q = 0.1 fixed, that the discrepancy decreases as r~.Extrapolating to rw = gives the point shown as a square n in fig. 3, which does indeed lie on the analytical curve. The mode is correctly represented in the circular case and has the right behaviour as described in Dewar et a!. [15] for the elliptical case although we did not make a detailed comparison in this last case. We found that as the ellipticity increases, the number of mesh points required increases also. We believe this is due to the pathological behaviour of the orthogonal coordinate system near the axis. In conclusion, we find that the code reproduces correctly the kink modes in the low beta limit. It should be pointed out that a choice of the finite elements which would violate what we have called the divergence condition, could never do this. As the current/ is lowered, the number of mesh points has to increase as an inverse power off to reproduce the kink. 00
=
__________
p(l
+
A
e2)2a2
2,
2
A1
q=
2~i +q2) (1 + 2 k(l +e2~~
(34)
ef
l—q
Note, that the growth-rate scales asj2, for a given safety factor q. The mode is only unstable in the range 0
4.2. The
strong
current regime
As a second test we choose the same equilibrium, but we take a larger current / = 0.4 and a conducting shell close to the plasma, r~ 2. We look for the unstable part of the spectrum in the range 0
q Fig. 3. Comparison between the numerical and analytical values of the growth-rates of the “m = 1” external kink for = 1 and e = 1.5. The solid lines are the analytic results for r~= =. The X and 0 represent the numerical results obtained with rw = 20 for e = 1 and a = 1.5, respectively. The 0 represents the result of an extrapolation to rw results at r~= 10, 20 and 40.
=
from the
320
D. Berger et al./Computation of the MHD spectrum of a straight plasma column
r2a -r
0 .25
.5
75
9
q
Fig. 4. Unstable part of the spectrum in the range 0
~-
1,
in the strong current regime, with the conducting shell at rw = 2.0. The solid lines represent the six most unstable modes for the circular case. The a represent the numerical results for an ellipticity of a = 1.5.
axis. Note that the kink itself becomes also internal near q = 1, the change in behaviour of the mode coinciding with the kink in the curve visible in the figure. There are in fact an infinite number of unstable modes accumulating at q = 0, F2 = 0. All these features are not visible in fig. 3 because of the scale used and because of the smallness of the current. The dotted circles represent the numerical results obtained with the 2-D code for an elongation ~ = 1.5, keeping N~= 12 fixed. The points are obtained by an extrapolation to N~= from N~= 12, 16 and 20 intervals, assuming a quadratic convergence. We see that most of the dependence of the growth-rates on the ellipticity is eliminated by normalizing F2 to w~.At long wavelengths the ellipticity destabilizes, while the conducting wall stabilizes. The logarithmic scale tends to attenuate the difference between the elliptical and circular case. Near q = 1 we only find two unstable modes and the results suggest a stabilizing effect of the ellipticity. This contradicts the conclusions of Lavel [7] which shows, at least in a
particular case, that the internal modes are strongly destablized by the ellipticity. To understand this apparent discrepancy we do a full convergence study of the growth-rate of the second mode, which shows already a strong stabilization for a fixed value of q = 1 and for two different elongations = i and e = 1.5. The results are shown in fig. 5 in which the square of the growth-rate is plotted versus N~2for a fixed N 11 = 12. In the circular case the convergence is indeed quadratic above N~ 20, but in the elliptic case no clear convergence law has emerged yet with the largest number of cells we can use, namely N~= 48. From the behaviour of the curve we see that a quadratic extrapolation underestimates the final value which will most probably be higher than the circular value. This shows that we should not trust the results for the slow growing modes. Even in the circular case, which looks good in fig. 5, there are difficulties in representing internal modes. For example, with the same current! = 0.4, no internal modes are found near q = 2. Also for smaller values of the current even the modes near q = 1 become impossible to find. The ellipticity only magnified the difficulty. The explanation of this difficulty is found by looking at the expression for Wi,, eq. (6). Internal modes are characterized by FX FY FZ 0. With polynomial finite elements F cannot vanish identically, resulting in a stabilizing effect analogous to the effect of a bad
-r Wa
00
Fig. 5. Plot of the square of 2,kthe eeping growth-rate N of an internal mode as a function of N~ 11 fixed. The • represents the value calculated with a 1-D code with the same number of elements.
D. Berger at al./Computation of the MHD spectrum of a straight plasma column
representation of the divergence, but smaller by a factorB~/B~. This explains that at the order of the kink we only need to satisfy the divergence condition while at the order of the internal modes we need to satisfy
321
both.
greatfully acknowledged. F. Debonneville, R. Schreiber and the staff of the Computing Center have been a great help in preparing and running the code. This work is supported in part by the Swiss National Science Foundation.
5. Conclusions
Appendix
From the examples presented above we can draw some interesting conclusions. With its judicious choice of regular basis functions, the code is capable of well reproducing kink modes for arbitrary /3, with a reasonable number of mesh points. This result has been obtained with a shearless equilibrium, but from our experience with the 1—D case, shear only affects the distribution of mesh points in the ~,D direction. In all our calculations we never observed spectral pollution [18]. For elongations larger than two the orthogonal coordinate system becomes very singular near the axis and it is difficult to represent accurately the kink modes. This could be taken care of by a change of coordinate system. The toroidal effects will not bring any new problems. The V ~ condition [seeeq. (12)] must be replaced by a similar condition on thefrom quantity 2)[24] where r denotes the distance the V (~/raxis. toroidal However, the weakly unstable internal modes are strongly stabilized. This comes from a bad representation of the operator F in eq. (6) which has to vanish identically on a rational surface. This condition can never be fulfilled with the described method, since the two terms in the F-operator have different functional dependences. This gives rise to a similar stabilizing effect as recognized in the one-dimensional code when we did not filfull the divergence condition [19]. Whereas the divergence condition affects the kink mode, a bad representation of the operator F only stabilizes the internal modes and does not influence the kink. This problem is specific to the two-dimensional case, since in the one-dimensional case Fis a number which can vanish exactly.
We limit ourselves to the case k ~ 0. The Green’s function Gk(x, x’) can be computed by reformulating the vacuum problem in terms of a scalar potential. We write =
V~(x,y,z),
(Al) ~uRe[~(x,y)e”~] where t~i(x,y)satisfies the Helmholz equation
a2
a2 ax2 ay2
(A2)
with the boundary conditions dtP/dflIcO, (A3) d~/dnI 5= F[X5]
/(.J5B~),
in which F is the operator defined in eq. (9), which reduces on the surface to
J5ax + ikB
F = -~-
-~-
(A4)
The vacuum contribution to the potential energy W~, eq. (7), can be written as an integral over the plasma— vacuum interface 1i 2ir Wv= 7Tf J5B~dx(~)~.
(A5)
We replace the x,y coordinate system by the elliptical system ~i, x defined by x
=
a
sh p cos x~
Acknowledgement
y
=
a ~.,/~TTi
ch p sin
Discussions with Professor J.L. Descloux and Dr. J. Rappaz on the mathematical implications of our choice of elements and on spectral pollution are
The plasma surface and the shell are surfaces of constant p, called ji~ and ~ respectively. On the surface, x as defined by eq. (A6) coincides with the x introduced inside the plasma.
(A6)
x.
D. Berger et al/Computation of the MHD spectrum of a straight plasma column
322
Eqs. (A2) and (A3) become
W~=
ap2 ax2 2
(A7)
acb/apI
aØ/apIc=0, a2(2
5=F[x5]
—
where A A 5 =J5B~. These equations are equivalent to the weak variational principle: find a function ~(p,x) periodic in x such that
fdxfd~
(p ~p
+~ ~
+ k2A2~6ø}
6~(p 5,x)O,
—fdxF[X5]
‘A8’ ‘
‘
1 (vacuum). for Eq. any(A8) 6cb EisHsolved by a finite element expansion of ~ and 60: =
~0Oi, x)
(A9)
~60 1e1(p,
x)
The same basis is used for the two variables. Substituting in eq. (A8) we obtain the inhomogeneous system of equations
(AlO) where S. is a vector given by S1
=
f
dxF[X51
e1(p5,
x).
(Al l~
[8]
0
0 is obtained in terms of S by solving eqs. (A10). We are interested only in the value of ~ x) which can be written
0(~5,x)= ~
G11S1e1(p5,x),
(A12)
1E5,JES
where the sum is carried only over elements which contribute on the surface. G~,is a real symmetric matrix. Replacing in eq. (AS)
G~
f
d~A5(F[X5]
)* e1(p5, x)
2ir
X
1) (sh2 p + cos2x). Note that
~
f
d~’A~F’[X~] ej(p5,
x’).
(Al 3)
This expression is identical to eq. (24) provided we identify Gk(x, x’) with the expression (32). To carry the calculation we cover the domain 0 ~ ~ 2ir, ~ ~ ~ ~ with a rectangular mesh (p1, x) and choose as elements pyramids centred on the mesh points with an hexagonal base (see ref. [12] for details). The elements e~(p5,x) are then indeed hat functions in x~ which achieves the identification with eq. (32). In practice we found that the mesh should be concentrated near the plasma surface in order to obtain the best results. The numerical approximation of Gk(x, x’), which is a singular integrable function, by regular functions, leads to a value of W~which is too small for large m, thus underestimating the stabilizing vacuum contribution.
References [1] K. Ham and R. Lust, Z. Naturforschung 13a (1958) 936. [2] T. Takeda, Y. Shimomura, M. Ohta and M. Yoshikawa, Phys. Fluids 15 (1972) 2193. [31 T.J.M. Boyd, G.A. Gardner and L.R.T. Gardner, Nucl. Fusion 13 (1973) 229. [41 K. Appert, D. Berger, R. Gruber and J. Rappaz, J. Comput. Phys. 18 (1975) 284. [51 C. Mercier, Nuel. Fusion 1 (1960) 47. [6] M.N. Bussac, D. Edery, G. Laval, R. Pellat and J.L. Soulé, 7th Europ. Conf. on Contr. Fusion and Plasma Phys., Vol.1 (1975) p. 100. [7] G. Laval, Phys. Rev. Letters 34 (1974) 1316. A. Sykes and J.A. Wesson, Nucl. Fusion 14 (1974) 645. [9] G. Bateman, W. Schneider and W. Grossmann, Nucl. Fusion 14(1974)669. [10] J. Wooten, H.R. Hicks, G. Bateman and R.A. Dory, ORNL TM 4784 (1974). [Ill A. Sykes and J.A. Wesson, 7th Europ. Conf. on Contr. Fusion and Plasma Phys., Vol.1(1975), p. 111. [12] O.C. Zienkiewicz, The Finite Element Method in Structural and Continuum Mechanics (McGraw-Hill, New York, 1971). [13] W. Kerner and H. Tasso, IAEA-CN-33/A13-l (Tokyo, 1975), p. 475. [14] W. Kerner, MAlT-i 229 (1976).
D. Berger et aL/Computation of the MHD spectrum of a straight plasma column [151 R.L. Dewar, R.C. Grimm, J.L. Johnson, E.A. Frieman, J.M. Greene and P.H. Rutherford, Phys. Fluids 17
(1974) 930. [16] R.L. Dewar, J.M. Greene, R.C. Grimm and J.L. Johnson, J. Comput. Phys. 17 (1975). [17] R.C. Grimm, J.M. Greene and J.L. Johnson, in: Methods of Computational Physics, Vol. 16, eds. B. Alder, S. Fernbach, J. Killeen and M. Rotenberg (Academic Press, New York, 1976, ch. 4. [18] J. Rappaz, PhD Thesis, EPF-Lausanne (1976).
323
[19] K. Appert, D. Berger, R. Gruber, F. Troyon and J. Rappaz, ZAMP 25 (1974) 229. [20J R. Gruber, PhD Thesis, EPF-Lausanne (1976). [21] R. Gruber, Contributed at the 2nd Europ. Conf. on Comp. Phys., Garching, 1976, Paper F2. [22] R. Gajewski, Phys. Fluids 15 (1972) 70. [23] K. Appert, D. Berger, R. Gruber, F. Troyon and K.V. Roberts, Comput. Phys. Commun. 10(1975) 11. [24] J.L. Soulé, private communication.