- Email: [email protected]
Numerical simulation of Reversed-Field Pinch dynamics
Computer Physics Communications 43 (1986) 17—28 North-Holland, Amsterdam
17
NUMERICAL SIMULATION OF REVERSED-FIELD PINCH DYNAMICS D.D. SCHNACK, D.C. BARNES, Z. MIKIC Science Applications International Corporation, P.O. Box 2351, La Jolla, CA 92038, USA
D.S. HARNED Courant Institute of Mathematical Sciences, New York University, 251 Mercer Street, New York, NY 10012, USA
E.J. CARAMANA and R.A. NEBEL Los Alamos National Laboratoiy, Los Alamos, NM 87545, USA
The Reversed-Field Pinch (RFP) is presented as a driven nonlinear dynamical system that evolves on time scales long compared with the fastest normal modes of the system. Numerical methods for the simulation of these processes, and examples of their application to experimentally observable phenomena, are presented.
1. Introduction The reversed-field pinch (RFP) is a toroidal Z-pinch that offers a promising alternative to the tokamak as a controlled fusion reactor. In the tokamak, confinement (radial force balance) is achieved by the balance between radial pressure gradients and stresses produced by the poloidal
The RFP derives its name from the observation that such plasmas exhibit robustly long lifetimes when the axial magnetic field at the outer boundary reverses direction relative to the field on axis (r 0). Experimentally, such states arise naturally whenever the applied axial electric field supports sufficient plasma current to raise the pinch parameter 9 = B9waii/Bz~v~~ 1.3. Another useful parame-
magnetic field; a large toroidal (axial, in the cylindrical approximation) magnetic field is required to stabilize the m = 1 ideal MHD mode. The size of this field must be large enough to assure that the safety factor q = rB2/RB9 (where r is the minor radius and R is the major radius of the torus) is greater than unity everywhere within the device, 2 Thus BZ/BO R/r>>confinement 1 and $ ~p/B~ (r/R) a ~ 1. In the RFP, results from balance between the stresses in the poloidal and toroidal (axial) magnetic fields. Thus BZ/BG 1, and q r/R << 1. Finite plasma pressure arises from any small imbalance in these stresses. Linear stability against ideal kink modes is achieved by a combination of high shear (q’/q) and the presence of an outer boundary that inhibits plasma displacements. Thus, to a good approximation, the RFP is characterized by the presence of nearly force-free magnetic fields (J X B = 0).
ter in the description of the RFP is the field reversal parameter F = BZ/BZ. These two parameters are dependent on the plasma current, the axial flux and the value of the axial magnetic field at the outer boundary. It is a remarkable fact that all RFP experiments are observed to seek values ofof F and 0approximately (F —0.1, the 0 same 1.5)operating independent initial toroidal flux or applied voltage (provided, of course, that sufficient voltage is supplied to achieve field reversal). In a landmark paper [1], Taylor showed that this behavior could be explained in a global sense by assuming that the plasma naturally seeks its lowest energy state (W) subject to the constraint that the axial flux (~)and the magnetic helicity (K), an integral measuring the mutual linkage of the poloidal and toroidal fluxes, remain constant. This simple variational calculation yields force
0010-4655/86/$03.50 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)
18
D.D. Schnack et at
/
Simulation of reversed-field pinch dynamics
free magnetic fields that exhibit axial field reversal whenever 0> 1.2. (States at their minimum value of W/K are said to be relaxed.) This is a simple and elegant theory that describes the global properties of the RFP. However, being of a variational nature, it does not address the dynamical mechanisms responsible for the relaxation process. Nor does it address what is perhaps the most interesting and important feature of RFP plasmas: their anomalously long lifetime in the presence of finite plasma resistivity, Once field reversal is achieved (and it occurs naturally with sufficient axial voltage) the amount of positive axial flux in the plasma is limited to that contained between the field reversal surface and the magnetic axis; only negative axial flux can be supplied to (or withdrawn from) the plasma by the external circuit. Thus the lifetime of the RFP should be limited by the resistive diffusion time of this positive flux. However, as noted above, it is an experimental fact that such plasmas can be maintained almost indefinitely at fixed values of F and 0 in the presence of an applied axial voltage, which supplies poloidal flux to the plasma; the lifetime of the discharge is limited primarily by the volt-seconds available in the external circuit. Poloidal flux is apparently dynamically converted into axial flux in sufficient amounts to maintain F and 0 at their proper values. This is referred to as the RFP dynamo. The operation of the RFP dynamo is best exhibited in a mode of operation called “ramp shots” [2]: after initial field reversal, the applied axial voltage is gradually increased (“ ramped”) to several times it original value. Both the current and the positive axial flux are observed to increase to many times their original values such as to maintain constant values of F and 0. By dynamical means the plasma creates positive and negative flux in equal amounts, and expels the excess negative flux to the external circuit. The process is accompanied by strong m 0 and 1 magnetic field oscillations, indicating that long wavelength MHD activity may be responsible for the dynamo process. Recently theoretical progress has been made in understanding the dynamics of the relaxation and dynamo processes. First, it was recognized [3] that =
transport phenomena (e.g., resistive diffusion) serve to increase the value of W/K relative to a completely relaxed state, (i.e.. K tends to dissipate faster than W). These studies further indicated that such one-dimensional processes could not produce field reversal or profile sustainment [4]. Concurrently, linear stability considerations [3] led to the conclusion that, because of the monotonically decreasing q(r) profile characteristic of RFP plasmas, many m I resistive kink (tearing) modes with different axial mode numbers were likely to be simultaneously active. (This is in contrast to the tokamak, where a single m I mode exists only if q(0) < 1.) These modes were thought to be important because of their known robust nonlinear behavior. Two-dimensional (single helicity) simulations [5] showed that the nonlinear evolution of one of these modes was sufficient to produce a helically deformed state with average field reversal. However, such helical deformations are not experimentally observed. Studies of the simultaneous nonlinear evolution of many such modes for experimentally relevant time scales required the development of efficient three-dimensional numerical algorithms [6,7]. The application of these techniques to the RFP showed that the nonlinear interaction of many m I resistive kink modes leads to a reduction of W/K (relaxation), to an increase in positive axial flux, and to the attainment of a (time averaged) symmetric reversed state [7,8]. (Resistive modes were invoked by Taylor to eliminate all but one of the infinite set of helicity invariants that exist in ideal MHD, although the details of the dynamics were not considered.) Such calculations have subsequently been coupled with one-dimensional transport to provide an explanation of experimentally observed sawtooth oscillations [8,9]. They have also directly simulated “ramp shot” operation. We are thus led to the following picture of the sustained RFP as a driven nonlinear dynamical system: the externally applied axial voltage provides poloidal flux (axial current) that is a continual source of free energy. Transport due to resistive diffusion causes the axial flux to decay and the plasma to evolve away form its desired minimum value of W/K. This serves to destabilize many long wavelength m 1 resistive kink modes =
=
=
=
D.D. Schnack et al.
/
Simulation of reversed-field pinch dynamics
whose subsequent nonlinear evolution and interaction rearrange the magnetic fields in such a way as to increase the positive axial flux, decrease W/K, and sustain field reversal. (Many sideband modes, principally m 0, are excited by this process and can be experimentally observed.) Transport processes then increase W/K, causing the cycle to repeat. The exact behavior of the pinch is thus determined by a balance between transport driving the system away from its relaxed state and nonlinear dynamics restoring it. At low 0( ~ 1.6) the timescales for these processes are in balance and a “quiescent” discharge (no strongly correlated signals) is observed. At higher 0 their timescales become more disparate, and coherent oscillations, such as the “sawtooth”, are observed. We emphasize that, except for the earliest stages of startup, the m 1 modes are always nonlinear. The role of linear theory is merely to identify the important modes. Crucial to the attainment the theoretical picture of RFP dynamics has been the development of numerical algorithms for the long-time-scale simulation of three-dimensional MHD systems. The resistive modes thought to be responsible for the RFP dynamo evolve on time scales long compared to the Alfvén transit time; the relaxation process may require a significant fraction of a resistive diffusion time. Of particular importance has been the emergence of codes based on the semi-implicit algorithm [7], in which the equation of motion is modified by the addition of a term, proportional to the time step, consisting of a linear spatial operator acting on the time derivative of the velocity. With the proper choice of this operator, a consistent implicit algorithm results that is unconditionally stable with respect to normal modes, and yet requires essentially the same computation time per timestep as an explicit method. Using this technique, computations have been performed in only a few hours, instead of hundreds of hours required that would have been required with previous methods. Thus fully three-dimensional, long time parameter studies can be routinely performed. In this paper we present several examples of the computation of nonlinear processes in the RFP. =
=
Section 2 contains a brief discussion of the
19
mathematical model, and of the spatial and temporal approximation required to obtain accurate, long time-scale nonlinear computations, including both the pseudospectral and semi-implicit methods. Specific examples are presented in section 3. These include the attainment and sustainment of field reversal for times exceeding the resistive diffusion time, the simulation of ramp shots, and a discussion of a possible mechanism for sawtooth oscillations. Section 4 contains the summary and conclusions.
2. Numerical methods 2.1. Mathematical model It is now well established that multi-dimensional nonlinear resistive magnetohydrodynamics is an excellent model for the description of the macroscopic dynamics of present magnetic fusion experiments. When the magnetic energy density in the confining magnetic field is large compared with the internal energy density, as is the case in the RFP, the low-frequency, long wavelength motions described by this model are governed primarily by the interaction between the total Lorentz body force and the changing magnetic field induced by fluid motion. It is then appropriate to ignore the pressure force and to remove the energy equation from the model. Then the remaining equations (which remain fully compressible) can be written (in a convenient set of nondimensional variables) as
SV x B
-~—- =
~ ~
=
—
—
(la)
i~J,
s~v ç7V + SJ x B +
v
~2
(ib)
~,
where A is the vector potential, B is the magnetic field, J is the current density, V is the velocity, p is the mass density, and we have chosen a gauge in which the gradient of the electrostatic potential vanishes. In these equations time is measured in units of the resistive diffusion time TR 4iia2/ c2 ~ and velocity is measured in units of the Alfvén velocity VA B 0/~ These equations are further simplified by assuming p to be uni=
=
20
D.D. Schnack et al.
/ Simulation
of reversed-fieldpinch dynamics
form in space and time, and setting p 1. Then S T~/T~ is the Lundquist number (the ratio of the resistive diffusion time to the Alfvén transit time TA a/VA), and r ‘rR/T.,iSC is the magnetic Prandtl number, Eqs. (1) are fully compressible and describe both shear and compressional Alfvén waves, as well as resistive instabilities and resistive and viscous damping. Their nonlinear evolution is particularly appropriate to the description of forcefree relaxation in the reversed-field pinch (RFP), and to other phenomena whose origin is primarily electromagnetic. Because of the widely separated time scales supported by these equations when S is large, their solution presents an extremely difficult test of any numerical method. The mathematical model described here supports a wide variety of solutions ranging from global kink and tearing modes to high frequency, short wavelength turbulence. We emphasize that our goal is to simulate the long time scale, nonlinear, three-dimensional evolution of long wavelength, low frequency motions for which geometnc effects such as magnetic shear, field line curvature and finite boundaries are important. We therefore employ temporal approximations in cylindrical geometry that permit time steps far in excess of those allowed by explicit methods. We make no attempt to simultaneously model turbulence, for which high resolution, short timestep explicit methods are appropriate [II].
the RFP as well, there is reason to believe that many modes may be equally important [3]; thus a large number of mode interactions are possible and no a priori mode selection process can be introduced. Direct evaluation of convolution sums in this case would be impractical, requiring e2(N2) operations (where N is the number of Fourier modes; typically N> 100). However, use of the Fast Fourier Transform (FFT) allows these terms to be evaluated in ~(N ln N) operations, thus becoming competitive in speed with finite differences, and with algorithms that retain only a few modes. These pseudospectral methods are briefly described in this section. We choose cylindrical coordinates (r, 0, z), 0 ~ r ~ ~ 0 0 ~ 2~,0 z L, and define the scale length a such that rmax I. We further assume that the solutions of (1) are periodic in with L 2’rrR, where R is the aspect ratio. Then introducing the discrete mesh 0~. 2i’r(j I)/M, j I, 2 M; Zh 2i’r(k I)R/N, k 1, 2 N: any quantity can be represented as a finite Fourier series
2.2. Spatial approximation: method
~
=
=
=
=
the pseudo-spectral
=
=
=
=
=
=
M/2
f~r,
0~,Z~)
— —
~-‘
fv,,(r)
~
rn=
—
‘~i/2+ I n= —M/2-1- I
>~ (2) where the complex finite Fourier coefficients f,~, are defined as .
•~ ~ r)
k1
—
f( r, 8,.
~)
/
Xe
The periodic nature of the poloidal (azimuthal) and toroidal (axial) directions in many fusion devices allows solutions to be represented by Fourier series in these coordinates. This offers several advantages over finite difference representations, including reduced phase error and cxponential (or high order) convergence properties [12]. Experience has shown that for the simulation of tokamak plasmas only a relatively small number of the possible Fourier modes are important in the dynamics [13]; thus a mode selection procedure may be used that allows for the direct evaluation of convolutions introduced by quadratic nonlineanities. While this may eventually bear out in
—
—
(3
i(,nO-f ti~/R) ‘
-
and the reality of f(r, 0, z) requires that [,,_,, f,,~’,,.The representation (2) is chosen because of the periodic boundary conditions in 0 and z, and the rapid convergence properties for smooth solutions [12]. Introducing (2) into (1) yields =
=S(vxB) 8t
—(i~
~7X ~7XA)
tn,,
““n,
(4a) V,n,n
at
=
—
S(v• ~zv),,,
+
S(J X B)
,
+v(~2V),,,,,,
(4b)
D.D. Schnack et al,
/
Simulation of reversed-field pinch dynamics
where A m,n and Vm,,, are the Fourier coefficients of the vector potential and velocity, and the subscript ( )m,n on a nonlinear term or operator represents the finite Fourier transform of that term. Eqs. (4) are a set of MN nonlinear partial differential equations in the radial coordinate, r, and the time, t, that describe the evolution of the complex Fourier coefficients ~ and ~ They are the equations that we solve numerically, The nonlinear terms appearing on the right hand side of (4) are evaluated with the aid of FFTs. For example: to evaluate (VX B)mn, transform the complex arrays V,,~and Bmn to configuration space to obtain the real arrays Vj,k and BJk, obtain the real product VJk X BJ,k in configuration space by simple multiplication, and then transform this product back to Fourier space. This process introduces aliasing errors, but they can be simply removed [6]. The resulting algorithm retains all the convergence and energy conservation properties of fully spectral methods; indeed, properly dealiased pseudospectral methods are mathematically equivalent to spectral methods [12]. Since the solutions of eqs. (4) are not periodic in the radial coordinate r, the techniques discussed above are not directly applicable. Similar problems are encountered in viscous hydrodynamic problems in nonperiodic domains, such as flow in pipes or rectangular channels. In these cases approximation in the nonperiodic coordinate is accomplished by expansion in Chebyshev polynomials [12], which have the desirable property that their zeros (collocation points) are densely spaced near the outer boundary giving a naturally concentrated mesh there. This is important since hydrodynamic boundary layers form in this region. In resistive MHD, the important boundary layers are associated with filamentation of current, occur internally, and may form spontaneously at nonpredetermined locations. Thus a uniform mesh, or one that is internally dense, is desirable. Expansion in Bessel functions, the natural functions for analytic investigation, is flawed by nonuniform convergence at the outer boundary and lack of a fast transform algorithm [12]. Thus we choose finite differences as the most efficient representation of the radial coordinate,
21
2.3. Temporal approximation: the semi-implicit method As discussed in section 1, the relevant dynamics that determine the evolution of the RFP evolve on time scales that are long compared to those associated with the fastest normal modes of the system. Stability restrictions placed on explicit temporal approximations can result in uneconomically small time steps, while the fully implicit treatment of nonlinear terms either requires iteration (as in finite difference methods [14]), or results in unacceptably large matrices (as in spectral methods). Recently, a new class of methods for solving the time-dependent MHD equations based on an algorithm developed for long range weather simulation [15,16] has been introduced [17,18]. These semi-implicit algorithms allow very large time steps, yet avoid the complexity and large memory requirements associated with implicit methods. In the semi-implicit method for MHD new terms are added to the time-discretized equations that do not affect the consistency of the solution, yet provide a simple and efficient means of enhancing stability. This technique was first applied to the MHD equations in order to eliminate the fast mode time step constraint [17]. Subsequently, this procedure was extended to eliminate the shear Alfvén time step constraint as well [18]. This method is unconditionally stable with respect to all Alfvén modes and consequently permits such large time steps that accuracy becomes the most important consideration in the choice of step size. The semi-implicit method may be implemented using a variety of different time advance schemes and forms for the semi-implicit terms. The choice of the particular semi-implicit method, both in the order of accuracy and the form of the semi-implicit terms, can affect the size of the allowable time step. The proper choices can lead to improved efficiency and accuracy for long-time scale three-dimensional nonlinear MHD computations. The essence of the semi-implicit method can be illustrated by considering the simple first order equation [16] a~a F 5 —.
/
t
—
i ~a
where iw is a frequency or, more generally, an
22
D, D. Schnack et at.
/
Simulation
operator that determines the normal modes of the system. A simple time difference is introduced on the left hand side of (5), i.e., aF,/at (F” ~ F”)/~t. Evaluating the right hand side of (5) at the old time step (F”) results in an explicit method; evaluation at the advanced time step (F”~) yields an implicit method. In the semi-implicit method the frequency (operator) iw is split as i~ i(a + /3), with part evaluated at the new time level and part at the old, i.e., =
=
F”~—F” __________
i(aF”~’+ $F”).
=
(6)
The algorithm introduces no damping, and is unconditionally stable provided a2> ,32• In practice, this may be accomplished by setting /3 w a, with the resulting algorithm =
F”~ F” —
_________
—
,
=
i~aF”+ iaF”~
—
ictF”
(7)
unconditionally stable provided a> w/2. Thus a term iaF is added and subtracted from the right hand side of (5), with one treated at the new time level and the other at the old; the original frequency (operator) is always treated at the old time level, The MHD equations (1) consist of two coupled equations of the form (5). For heuristic purposes, we consider the simple linear system [7] au,/at=Bv, (8a) av/at Au, (8h) =
where A and B are linear operators independent of time, and u and v are vectors. Eqs. (8a, b) can be combined to yield a wave equation for v
a2v/at2
=
Lv,
(9)
where L AB. An identical equation holds for u if A and B commute. Eqs. (8a, b) are analogous to (la, b), and eq. (9) is analogous to the linearized MHD wave equation =
a2 v/at2
=
of
reversed-fieldpinch dynamics
where G is an operator yet to be defined. Combining (ha, b) we find, analogous to (9), Lv, (12) (1 a z~tG)a2v/at which is the wave equation for the modified systern. It can be shown [7] that if the semi-implicit operator G is chosen such that, for the long wavelength modes of interest, its eigenvalues and eigenvectors closely approximate those of L. then the solutions of (12) are also a close approximation to those of (9). Furthermore, the coefficient a appearing in (lib) and (12) can he chosen to yield an unconditionally stable algorithm [7,17.18]. The power of the semi-implicit method lies in choosing G such that it is easily inverted, and in choosing a to yield unconditional stability. In =
—
MHD a particularly useful form for G is obtained from the linearized operator LV=c’xv~’x(VxB 0)xB0 (13) by replacing B0 with a constant vector C0, and then retaining only terms proportional to where a r, 0, z [18]. This yields a simple block tn-diagonal system that is easily inverted. If a modified leap frog algorithm is used for the cxplicit terms [6], unconditional stability then results provided 2/4 1 (14) S~B~k,~,,~t is simultaneously satisfied for a r, 0, z. This serves to determine the coefficients C~.When the right hand side of (14) is negative the algorithm is explicitly stable and C~is set to zero. The semi-implicit method is a simple, fast and =
—
=
accurate algorithm for three-dimensional MHD simulations. It is simple because it is extremely easy to implement, even in existing codes. It is fast because the computation time per time step is dominated by the evaluation of nonlinear convolutions. Since the semi-implicit method requires
ç~x ç~’x (V>< B
0) x B0.
(10)
only the solution of linear tridiagonal systems, it
In applying the semi-implicit method, we introduce the modified system
requires essentially the same CPU time per time step as explicit methods, while allowing time steps that may be orders of magnitude larger. It is accurate because the semi-implicit operator can be chosen to mimic the MHD operator for modes of
au/at By, av,/at Au + a ~1tGav/at, =
=
(ha) (lIb)
D.D. Schnack et al.
/
Simulation of reversed-field pinch dynamics
23
io~
interest. The abillity to simulate driven systems for times that are significant fractions of the resistive diffusion time is essential to understanding RFP dynamics.
io-~ 1 o-~ ‘-I
C.,
3. Examples of RFP dynamics
3.1. Sustained field reversal, t ~
~ T~
[7] c.z~
The essence of the RFP is the sustainment of field reversal and toroidal flux in the presence of an applied toroidal voltage for times exceeding the resistive diffusion time of the system. We simulate this effect by posing an initial value problem in which the initial conditions consist of a one-dimensional equilibrium characterized by the safety factor profile (is) q(r) 0.4(1 1.8748r2 + 0.83232r4), =
—
which we perturb with both the (m 1, n —1) and (m 1, n —2) modes. We sustain the current by applying an axial electric field at the wall that is constant in space and time. The average poloidal electric field at the wall is set to zero so as to conserve toroidal flux. The resistivity is constant in space and time, with Lundquist numben S iO~.With unit aspect ratio, the spatial mesh is 65 radial, 8 poloidal and 16 axial points, This case has been the subject of a previous study [19], where its nonlinear evolution was computed for up to 200 Alfvén times (0.2 resistive diffusion times) with several nonlinear 3-D MHD codes. There it was found that the nonlinear growth and saturation of the unstable nonresonant (1, —2) mode led to the generation and sustainment of an average reversed axial field (Br) at the wall, as is observed in reversed-field pinch (RFP) experiments. These simulations indicated the attainment of a steady field-reversed state dominated by an m 1, n —2 helical perturbation that was maintained against resistive diffusion by finite flow. Using the semi-implicit method we have cornputed the nonlinear evolution of this case for times in excess of one resistive diffusion time (> 1000 Alfvén times). In fig. 1 we plot the energy in various modes as a function of time (measured =
=
LAB
M
N
LAB
M
N
1013
01 32
0 1 02 1 -1 1-2
5 6 87
1 -4 2-1 2 .2 2-3
1015
4
1
9
2
.2
.4
-3
.6
.8
-4
1.0
t/TR
Fig. 1. Energy in selected (m, n) helical modes as a function of time for a driven RFP simulation. Time is measured in resistive diffusion times (TR). Note the nonlinear transition and modelocking at (‘r O.5Tpj.
=
=
=
=
10h1
=
in resistive times). The initially unstable nonresonant (m 1, n —2) mode grows and saturates in about 0.02 resistive times, producing a helical state in agreement with previous results. However, this state is in turn unstable to the resonant (1, 3) mode, which saturates in about 0.1 resistive times at an amplitude approximately equal to that of the =
=
—
.66 .64
-______________________
__________
.62 .60 .58 56
~ 52 50 .48 .46 .42 .40 1.0 tlTR Fig. 2. Safety factor at r = 0 as a function of time for the driven RFP simulation shown in fig. 1. ‘
.2.4.6.8
24
D. D. Schnack ci at.
/
Simulation of reversed-field pinch dynamics
(1, —2). The (1, —2) and (1, —3) then proceed to periodically exchange energy for 0.4 resistive times, resulting in a time-dependent, field-reversed state. At t O.STR a further nonlinear transition occurs, with the (1. —3) and (1, —4) modes inter___________________________________ 200~ ‘i’ l7OkA
changing roles. For t O.7T~ a fully three-dimensional resistive steady state results. The nonlinear periodic energy exchange between the (1, —2) and (1, —3) modes occurring for O.ISTR ~ t ~ O.ST~ can be understood by considering q(0), the value of the safety factor at r 0. In fig. 2 this quantity is plotted as a func=
a
150
RAMP
100
60 kA 50
,I,,,’fl
0 0.15
__
‘‘‘~‘‘~‘‘‘l’
2.0
b
I
I
170 kA
0.10
I~IIb
LA~
-J U-
~kA~
0.05
H
0 I
_____________________________
0.08
I
I
I
-
C
0.06
RAMP
0.04 60 kA~
0.02
.15 25 .05
0 -0.02 _______________________________________
0
5
10 TIME (ms)
15
-
20
Fig. 3. RFP ramp shot experimental data taken from ref. [2]: (a) toroidal current versus time; (b) average toroidal flux versus time; (c) toroidal magnetic field at the wall (vacuum liner) versus time,
0.0
.2
.4 TIME (TR)
Fig. 4. Numerical simulation of RFP ramp shot: (a) toroidal current versus time; (Is) toroidal flux versus time; (c) field reversal parameter F versus time,
D.D. Schnack et a!.
/ Simulation
of reversed-fieldpinch dynamics
tion of time. Initially, q(0) <0.5, so that the unstable (1, —2) mode is nonresonant from above, As is well known [3], the growth of such a mode raises q(0) (“second reconnection”), increasing the shear, and stabilizing the mode. This behavior is seen at early times (t < O.1TR) in fig. 2. However, this profile modification is sufficient to destabilize the resonant (1, —3) mode, whose evolution is such as to lower q(0), i.e., to remove the unstable
25
resonance (“first reconnection”) [3]. This is seen in fig. 2 for O.1TR ~ t ~ O.2T~.The achievement of q(0) <0.5 again destabilizes the (1, —2), causing q(0) to increase, which in turn destabilizes the (1, 3), causing q(0) to decrease. This periodic oscillation about q(0) 0.5 is apparent in fig. 2 for O.2Tp, ~ I ~ O.STR. Finally, at t O.STR the ifl terchange of the (1, —3) and (1, —4) causes q(0) to increase abruptly, followed by the attainment —
=
t=10
t=20 ~“
“•:‘~~‘~
‘~~“
:~
,~
•
20
f
20
‘~
16 12
I
,,~
8
•
0.0
.20
.~ ••~
z
,
12
,J”
.60
.80
1,00
‘5’”
~“
0.0
.20
.40
t=30
‘~
1
20~
1.00
24 20
16~
I
c
16 ‘~‘“•‘‘~ ‘.
8
0 0.0
.80
28
24
12
.60
t=40
28
I
•
—r——.
—r—-—+
~
~
,,
,i.,.
8
.40
f
‘~
16
z
I
. B,,
,‘
‘•
..,‘
‘:
,
“
,‘
, ,.
12
‘~
8
~
~
.20
.40
.60
—r—-—
.80
1.00
0 0.0
:,
‘
,,‘
‘
‘~
‘:
“‘
‘‘:..
‘
‘
,
—‘
.20
.40
.60
.80
1.00
—r———~
Fig. 5. Field line plots (surfaces of section) in the (r, z) plane at various times for resistive mode evolution in a low 9 RFP discharge. Note that good flux surfaces remain near the wall (r = 1).
D.D.
26
/ Simulation
Schnack ci a!.
of a resistive steady state. This second nonlinear transition and mode-locking is presently being studied. 3.2. Simulation of ramp-shots [10]
As mentioned in the introduction, one of the most interesting reversed-field pinch (RFP) cxperimental results occurs in discharges where the toroidal current is slowly raised after toroidal field reversal is established. These discharges exhibit the RFP dynamo effect in its most unambiguous
of reversed-field pinch dynamics
form showing the strong nonlinear coupling of the external poloidal and toroidal magnetic field circuits via the plasma—magnetic-field configuration. In this section we present the results of the application of the numerical methods described in section 2 to this problem. In fig. 3 ZT-40M data from ref. [2] are shown. Parts a—c give the toroidal current, average toroidal flux and the toroidal magnetic field at the wall (vacuum liner), respectively, as a function of time for three different discharges; two have relatively flat toroidal current at 60 and 170 kA the third is
t=10
t=20
28
~‘
I’,
•,i,’
,‘
‘ ‘..~
‘‘‘‘
28
,
—
24
I,
‘C’
, ‘
‘‘
‘
, ‘
•:~
‘
‘~
.,,,~,s;i
24 “
“•
‘‘4’’,,
,
20
.
‘‘
.
‘.
‘..
,
20
__I
_
0.0
.20
.40
.60
t=30
.80
1,00
0.0
.20
.40
‘
.
.“
,
~.
•‘:
.
~‘,‘
,,
28
24
24
20
20
16
16
‘
‘
,
.1.
.
‘. ,
,
.~
~‘
1’ ‘
‘‘.‘.
~
~
‘,
,.
~.
*5
I]
..
.‘,,
S
12
I
8
. ,
1.00
.
•
12
.80
t=45_________________
_______
‘.5 5,
28
.60
,
8
,,
,•,,,,
,
,
•‘.~‘.
,‘~:,
, •.‘
.~
~,
~
.:‘.
4 1 .80 1.00 0,0 .20 .40 .60 .80 1.00 Fig. 6. Field line plots (surfaces of section) in the (r, z) plane at various times for resistive mode evolution in a high 8 RFP discharge. 0.0
.20
.40
.60
Note that outer flux surfaces are destroyed, and then reform,
D.D. Schnack et a!. / Simulation of reversed-field pinch dynamics
a “ramp” discharge where the current goes between these two values in about 15 ms. Notice the deeper reversal of the ramp case in part c; enhanced MHD fluctuations are clearly visible in both parts b and c of fig. 3 for the rising currçnt case. In fig. 4 we display the results of the numerical simulation of such a ramped current discharge [10]. The initial conditions consisted of a one-dimensional equilibrium with magnetic perturbations added for modes m = 1 and n = —2 to —5 inclusive. The Lundquist number S is 2300. The applied axial voltage at the wall is linearly increased from its initial value with a doubling time of 0.5 TR. The poloidal voltage is adjusted to keep the pinch parameter 0 approximately constant at its initial value of 1.58. Fig. 4a shows the toroidal current I4~~ fig. 4b the toroidal flux TF, and fig. 4c the field reversal parameter F, as functions of time (measured in units of TR). Notice that F remains negative and approximately constant throughout the run, even though the toroidal flux has increased significantly, indicating that the plasma has produced positive and negative flux in equal amounts and expelled the excess negative flux to the external circuit. The modes responsible for this enhanced dynamo effect have been identified as long-wavelength m = 1 helical resistive modes. Thus significant RFP physics lies at long wavelength. Further details and discussion of these results can be found in ref. [10]. 3.3. Sawtooth oscillations Experimental observations have demonstrated that long-wavelength, low-frequency signals are the dominant structures in ZT-40M [20,21]. For 0 1.55 the discharge exhibits large-amplitude sawtooth oscillations in the soft X-ray flux, The sawtooth crash is associated with large m = 0 and 1 magnetic perturbations measured at the wall [20]. Prior to the crash, m = 1 precursor oscillations are observed in the soft X-ray flux and sometimes in the m = 1 edge magnetic probe. For 0 1.55, the sawtooth oscillations are smaller, the m = 0 magnetic perturbation becomes imperceptible, and a continuous m = 1 oscillations appears.
27
Previous numerical simulations [8] have demonstrated that the m = 1 modal activity such as that described in section 3.1, leads to destruction of flux surfaces and the presence of stochastic field lines in the core of the plasma (this is shown in fig. 5). These m = 1 modes couple nonlinearly to m = 0 modes that are resonant at the field reversal surface. At large 0, these driven m = 0 modes grow large enough to overlap with the m 1 modes active in the plasma core. (The m = 0 islands can become large because, at high 9, the field reversal surface is relatively farther from the outer boundary.) The discharge then rapidly becomes completely stochastic (see fig. 6), with an implied subsequent rapid outward transport of thermal energy. At low 0 (fig. 5) the stochastic region resulting from m = 1 interactions is limited to the core, with good flux surfaces, and hence confinement, remaining at the edge. This picture of “edge confinement” is consistent with the cxperimentally observed sensitivity of the RFP to field errors. Toroidal flux generation (dynamo) is experimentally observed to be associated with the sawtooth crash [19]; similar flux generation has been observed in the numerical simulations at high 0 [8]. At low 0, confinement is maintained and a quasicontinuous level of m = 1 turbulence remains. Thus, the following picture of high 0 RFP sawtooth oscillations emerges: transport processes cause temperature and current to gradually peak on-axis, thus lowering q(0). When the q profile becomes sufficiently flattened, several m = 1 resistive modes are destabilized. Their nonlinear interaction causes a stochastic core to grow. Simultaneously, m = 0 magnetic islands are driven at the reversal surface. When these overlap with the m = 1 modes the entire discharge becomes stochastic causing rapid heat transport to the wall, a lowering of core temperature, flattening of current, raising of q(0), and a subsequent stabilization of the m = 1 activity. The process then repeats.
4. Summary and conclusions A picture of the RFP as a driven, nonlinear, dynamical, three-dimensional system has been
28
D.D. Schnack er at. /
Simulation of reversed-field pinch dynamics
given. Numerical methods appropriate for the accurate simulation of these long time scale processes have been presented. We have presented examples of the application of these methods to essential experimentally observed processes: long time driven sustainment, ramped current operation and sawtooth oscillations. These calculations have served to demonstrate that long wavelength, lowfrequency helical instabilities are suffictent to cxplain these observed characteristics of the RFP. Acknowledgements The authors wish to thank Drs. B. Carraras, J, Holmes, A. Mirin, W. Matthaeus, D. Montgomery, C. Rowdyshrub, L. Turner and J. Wesson for advice, collaboration, discussion and criticism throughout all phases of the work reviewed in this paper. This work was performed under contract DE-ACO3-83ER53150 with the US Department of Enerov References [1] J.B. Taylor, Phys. Rev. Lett. 33 (1974) 139. [2] J.A. Phillips et al., Los Alamos National Laboratory Report LA-10060-MS (May 1984). [3] E.J. Caramana, R.A. Nebel and D.D. Schnack. Phys. Fluids 26 (1983) 1305. [4] E.J. Caramana and D.A. Baker, Nucl. Fusion 24 (1984) 423.
[5] A.Y. Aydemir and D.B. Barnes, Phys. Rev. Lett. 52 (1984) 930. [6] D.D. Schnack, D.C. Baxter and E.J. Caramana, J. Cornput. Phys. 55(1984) 485. [7] D.D. Schnack, D.C. Barnes, Z. Mikic, D.S. Harned and E.J. Caramana, SAIC Report SAIC-86/3022-APPAT-83; J. Comput. Phys. (1986) submitted. [8] D.D. Schnack, E.J. Caramana and R.A. Nebel, Phys. Fluids 28 (1985) 321. [9] K.A. Werley, R.A. Nebel and GA. Wurden, Phys. Fluids 28 (1985) 1450 [101E.J. Caramana and D.D. Schnack, Los Alamos National Laboratory Report LA-UR-86-965, Phys. Fluids (1986) submitted. [11] D.G. Fox and S.A. Orszag, J. Comput. Phys. 11 (1973) 612. [12] D. Gottlieb and S.A. Orszag, Numerical Analysis of Spectral Methods: Theory and Application (SIAM. Philadelphia, 1977). [13] HR. Hicks, BA. Carraras, J.A. Holmes, D.K. Lee and B V. Waddell, J. Comput. Phys. 44 (1981) 46 [14] C.H. Finan III and 1. Killeen. Comput. Phys. Commun. 24 (1981) 441. [15] A.J. Robert, in: Proc. WMO/IUGG Symp. on Numencal Prediction, Tokyo Methods (1969). Used in Atmospheric [16] Weather A.J. Robert, in: Numerical Models. WMO GARP Publication Series, no. 17. vol. 2 (1979). [17] D.S. Harned and W. Kerner, J. Comput. Phys. 60 (1985) 62. [181 D.S. Harned and D.D. Schnack. J. Comput. Phys. 65 (1986) 57. [19] A.Y. Aydemir, D.C. Barnes, El. Caramana, A.A. Mirin, R.A. Nebel, D.D. Schnack and AG. Sgro, Phys. Fluids 28 (1985) 898. [20] R.G. Watt and R.A. Nebel, Phys. Fluids 26 (1983) 1168. [21] GA. Wurden, Phys. Fluids 27 (1984) 551.
17
NUMERICAL SIMULATION OF REVERSED-FIELD PINCH DYNAMICS D.D. SCHNACK, D.C. BARNES, Z. MIKIC Science Applications International Corporation, P.O. Box 2351, La Jolla, CA 92038, USA
D.S. HARNED Courant Institute of Mathematical Sciences, New York University, 251 Mercer Street, New York, NY 10012, USA
E.J. CARAMANA and R.A. NEBEL Los Alamos National Laboratoiy, Los Alamos, NM 87545, USA
The Reversed-Field Pinch (RFP) is presented as a driven nonlinear dynamical system that evolves on time scales long compared with the fastest normal modes of the system. Numerical methods for the simulation of these processes, and examples of their application to experimentally observable phenomena, are presented.
1. Introduction The reversed-field pinch (RFP) is a toroidal Z-pinch that offers a promising alternative to the tokamak as a controlled fusion reactor. In the tokamak, confinement (radial force balance) is achieved by the balance between radial pressure gradients and stresses produced by the poloidal
The RFP derives its name from the observation that such plasmas exhibit robustly long lifetimes when the axial magnetic field at the outer boundary reverses direction relative to the field on axis (r 0). Experimentally, such states arise naturally whenever the applied axial electric field supports sufficient plasma current to raise the pinch parameter 9 = B9waii/Bz~v~~ 1.3. Another useful parame-
magnetic field; a large toroidal (axial, in the cylindrical approximation) magnetic field is required to stabilize the m = 1 ideal MHD mode. The size of this field must be large enough to assure that the safety factor q = rB2/RB9 (where r is the minor radius and R is the major radius of the torus) is greater than unity everywhere within the device, 2 Thus BZ/BO R/r>>confinement 1 and $ ~p/B~ (r/R) a ~ 1. In the RFP, results from balance between the stresses in the poloidal and toroidal (axial) magnetic fields. Thus BZ/BG 1, and q r/R << 1. Finite plasma pressure arises from any small imbalance in these stresses. Linear stability against ideal kink modes is achieved by a combination of high shear (q’/q) and the presence of an outer boundary that inhibits plasma displacements. Thus, to a good approximation, the RFP is characterized by the presence of nearly force-free magnetic fields (J X B = 0).
ter in the description of the RFP is the field reversal parameter F = BZ/BZ. These two parameters are dependent on the plasma current, the axial flux and the value of the axial magnetic field at the outer boundary. It is a remarkable fact that all RFP experiments are observed to seek values ofof F and 0approximately (F —0.1, the 0 same 1.5)operating independent initial toroidal flux or applied voltage (provided, of course, that sufficient voltage is supplied to achieve field reversal). In a landmark paper [1], Taylor showed that this behavior could be explained in a global sense by assuming that the plasma naturally seeks its lowest energy state (W) subject to the constraint that the axial flux (~)and the magnetic helicity (K), an integral measuring the mutual linkage of the poloidal and toroidal fluxes, remain constant. This simple variational calculation yields force
0010-4655/86/$03.50 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)
18
D.D. Schnack et at
/
Simulation of reversed-field pinch dynamics
free magnetic fields that exhibit axial field reversal whenever 0> 1.2. (States at their minimum value of W/K are said to be relaxed.) This is a simple and elegant theory that describes the global properties of the RFP. However, being of a variational nature, it does not address the dynamical mechanisms responsible for the relaxation process. Nor does it address what is perhaps the most interesting and important feature of RFP plasmas: their anomalously long lifetime in the presence of finite plasma resistivity, Once field reversal is achieved (and it occurs naturally with sufficient axial voltage) the amount of positive axial flux in the plasma is limited to that contained between the field reversal surface and the magnetic axis; only negative axial flux can be supplied to (or withdrawn from) the plasma by the external circuit. Thus the lifetime of the RFP should be limited by the resistive diffusion time of this positive flux. However, as noted above, it is an experimental fact that such plasmas can be maintained almost indefinitely at fixed values of F and 0 in the presence of an applied axial voltage, which supplies poloidal flux to the plasma; the lifetime of the discharge is limited primarily by the volt-seconds available in the external circuit. Poloidal flux is apparently dynamically converted into axial flux in sufficient amounts to maintain F and 0 at their proper values. This is referred to as the RFP dynamo. The operation of the RFP dynamo is best exhibited in a mode of operation called “ramp shots” [2]: after initial field reversal, the applied axial voltage is gradually increased (“ ramped”) to several times it original value. Both the current and the positive axial flux are observed to increase to many times their original values such as to maintain constant values of F and 0. By dynamical means the plasma creates positive and negative flux in equal amounts, and expels the excess negative flux to the external circuit. The process is accompanied by strong m 0 and 1 magnetic field oscillations, indicating that long wavelength MHD activity may be responsible for the dynamo process. Recently theoretical progress has been made in understanding the dynamics of the relaxation and dynamo processes. First, it was recognized [3] that =
transport phenomena (e.g., resistive diffusion) serve to increase the value of W/K relative to a completely relaxed state, (i.e.. K tends to dissipate faster than W). These studies further indicated that such one-dimensional processes could not produce field reversal or profile sustainment [4]. Concurrently, linear stability considerations [3] led to the conclusion that, because of the monotonically decreasing q(r) profile characteristic of RFP plasmas, many m I resistive kink (tearing) modes with different axial mode numbers were likely to be simultaneously active. (This is in contrast to the tokamak, where a single m I mode exists only if q(0) < 1.) These modes were thought to be important because of their known robust nonlinear behavior. Two-dimensional (single helicity) simulations [5] showed that the nonlinear evolution of one of these modes was sufficient to produce a helically deformed state with average field reversal. However, such helical deformations are not experimentally observed. Studies of the simultaneous nonlinear evolution of many such modes for experimentally relevant time scales required the development of efficient three-dimensional numerical algorithms [6,7]. The application of these techniques to the RFP showed that the nonlinear interaction of many m I resistive kink modes leads to a reduction of W/K (relaxation), to an increase in positive axial flux, and to the attainment of a (time averaged) symmetric reversed state [7,8]. (Resistive modes were invoked by Taylor to eliminate all but one of the infinite set of helicity invariants that exist in ideal MHD, although the details of the dynamics were not considered.) Such calculations have subsequently been coupled with one-dimensional transport to provide an explanation of experimentally observed sawtooth oscillations [8,9]. They have also directly simulated “ramp shot” operation. We are thus led to the following picture of the sustained RFP as a driven nonlinear dynamical system: the externally applied axial voltage provides poloidal flux (axial current) that is a continual source of free energy. Transport due to resistive diffusion causes the axial flux to decay and the plasma to evolve away form its desired minimum value of W/K. This serves to destabilize many long wavelength m 1 resistive kink modes =
=
=
=
D.D. Schnack et al.
/
Simulation of reversed-field pinch dynamics
whose subsequent nonlinear evolution and interaction rearrange the magnetic fields in such a way as to increase the positive axial flux, decrease W/K, and sustain field reversal. (Many sideband modes, principally m 0, are excited by this process and can be experimentally observed.) Transport processes then increase W/K, causing the cycle to repeat. The exact behavior of the pinch is thus determined by a balance between transport driving the system away from its relaxed state and nonlinear dynamics restoring it. At low 0( ~ 1.6) the timescales for these processes are in balance and a “quiescent” discharge (no strongly correlated signals) is observed. At higher 0 their timescales become more disparate, and coherent oscillations, such as the “sawtooth”, are observed. We emphasize that, except for the earliest stages of startup, the m 1 modes are always nonlinear. The role of linear theory is merely to identify the important modes. Crucial to the attainment the theoretical picture of RFP dynamics has been the development of numerical algorithms for the long-time-scale simulation of three-dimensional MHD systems. The resistive modes thought to be responsible for the RFP dynamo evolve on time scales long compared to the Alfvén transit time; the relaxation process may require a significant fraction of a resistive diffusion time. Of particular importance has been the emergence of codes based on the semi-implicit algorithm [7], in which the equation of motion is modified by the addition of a term, proportional to the time step, consisting of a linear spatial operator acting on the time derivative of the velocity. With the proper choice of this operator, a consistent implicit algorithm results that is unconditionally stable with respect to normal modes, and yet requires essentially the same computation time per timestep as an explicit method. Using this technique, computations have been performed in only a few hours, instead of hundreds of hours required that would have been required with previous methods. Thus fully three-dimensional, long time parameter studies can be routinely performed. In this paper we present several examples of the computation of nonlinear processes in the RFP. =
=
Section 2 contains a brief discussion of the
19
mathematical model, and of the spatial and temporal approximation required to obtain accurate, long time-scale nonlinear computations, including both the pseudospectral and semi-implicit methods. Specific examples are presented in section 3. These include the attainment and sustainment of field reversal for times exceeding the resistive diffusion time, the simulation of ramp shots, and a discussion of a possible mechanism for sawtooth oscillations. Section 4 contains the summary and conclusions.
2. Numerical methods 2.1. Mathematical model It is now well established that multi-dimensional nonlinear resistive magnetohydrodynamics is an excellent model for the description of the macroscopic dynamics of present magnetic fusion experiments. When the magnetic energy density in the confining magnetic field is large compared with the internal energy density, as is the case in the RFP, the low-frequency, long wavelength motions described by this model are governed primarily by the interaction between the total Lorentz body force and the changing magnetic field induced by fluid motion. It is then appropriate to ignore the pressure force and to remove the energy equation from the model. Then the remaining equations (which remain fully compressible) can be written (in a convenient set of nondimensional variables) as
SV x B
-~—- =
~ ~
=
—
—
(la)
i~J,
s~v ç7V + SJ x B +
v
~2
(ib)
~,
where A is the vector potential, B is the magnetic field, J is the current density, V is the velocity, p is the mass density, and we have chosen a gauge in which the gradient of the electrostatic potential vanishes. In these equations time is measured in units of the resistive diffusion time TR 4iia2/ c2 ~ and velocity is measured in units of the Alfvén velocity VA B 0/~ These equations are further simplified by assuming p to be uni=
=
20
D.D. Schnack et al.
/ Simulation
of reversed-fieldpinch dynamics
form in space and time, and setting p 1. Then S T~/T~ is the Lundquist number (the ratio of the resistive diffusion time to the Alfvén transit time TA a/VA), and r ‘rR/T.,iSC is the magnetic Prandtl number, Eqs. (1) are fully compressible and describe both shear and compressional Alfvén waves, as well as resistive instabilities and resistive and viscous damping. Their nonlinear evolution is particularly appropriate to the description of forcefree relaxation in the reversed-field pinch (RFP), and to other phenomena whose origin is primarily electromagnetic. Because of the widely separated time scales supported by these equations when S is large, their solution presents an extremely difficult test of any numerical method. The mathematical model described here supports a wide variety of solutions ranging from global kink and tearing modes to high frequency, short wavelength turbulence. We emphasize that our goal is to simulate the long time scale, nonlinear, three-dimensional evolution of long wavelength, low frequency motions for which geometnc effects such as magnetic shear, field line curvature and finite boundaries are important. We therefore employ temporal approximations in cylindrical geometry that permit time steps far in excess of those allowed by explicit methods. We make no attempt to simultaneously model turbulence, for which high resolution, short timestep explicit methods are appropriate [II].
the RFP as well, there is reason to believe that many modes may be equally important [3]; thus a large number of mode interactions are possible and no a priori mode selection process can be introduced. Direct evaluation of convolution sums in this case would be impractical, requiring e2(N2) operations (where N is the number of Fourier modes; typically N> 100). However, use of the Fast Fourier Transform (FFT) allows these terms to be evaluated in ~(N ln N) operations, thus becoming competitive in speed with finite differences, and with algorithms that retain only a few modes. These pseudospectral methods are briefly described in this section. We choose cylindrical coordinates (r, 0, z), 0 ~ r ~ ~ 0 0 ~ 2~,0 z L, and define the scale length a such that rmax I. We further assume that the solutions of (1) are periodic in with L 2’rrR, where R is the aspect ratio. Then introducing the discrete mesh 0~. 2i’r(j I)/M, j I, 2 M; Zh 2i’r(k I)R/N, k 1, 2 N: any quantity can be represented as a finite Fourier series
2.2. Spatial approximation: method
~
=
=
=
=
the pseudo-spectral
=
=
=
=
=
=
M/2
f~r,
0~,Z~)
— —
~-‘
fv,,(r)
~
rn=
—
‘~i/2+ I n= —M/2-1- I
>~ (2) where the complex finite Fourier coefficients f,~, are defined as .
•~ ~ r)
k1
—
f( r, 8,.
~)
/
Xe
The periodic nature of the poloidal (azimuthal) and toroidal (axial) directions in many fusion devices allows solutions to be represented by Fourier series in these coordinates. This offers several advantages over finite difference representations, including reduced phase error and cxponential (or high order) convergence properties [12]. Experience has shown that for the simulation of tokamak plasmas only a relatively small number of the possible Fourier modes are important in the dynamics [13]; thus a mode selection procedure may be used that allows for the direct evaluation of convolutions introduced by quadratic nonlineanities. While this may eventually bear out in
—
—
(3
i(,nO-f ti~/R) ‘
-
and the reality of f(r, 0, z) requires that [,,_,, f,,~’,,.The representation (2) is chosen because of the periodic boundary conditions in 0 and z, and the rapid convergence properties for smooth solutions [12]. Introducing (2) into (1) yields =
=S(vxB) 8t
—(i~
~7X ~7XA)
tn,,
““n,
(4a) V,n,n
at
=
—
S(v• ~zv),,,
+
S(J X B)
,
+v(~2V),,,,,,
(4b)
D.D. Schnack et al,
/
Simulation of reversed-field pinch dynamics
where A m,n and Vm,,, are the Fourier coefficients of the vector potential and velocity, and the subscript ( )m,n on a nonlinear term or operator represents the finite Fourier transform of that term. Eqs. (4) are a set of MN nonlinear partial differential equations in the radial coordinate, r, and the time, t, that describe the evolution of the complex Fourier coefficients ~ and ~ They are the equations that we solve numerically, The nonlinear terms appearing on the right hand side of (4) are evaluated with the aid of FFTs. For example: to evaluate (VX B)mn, transform the complex arrays V,,~and Bmn to configuration space to obtain the real arrays Vj,k and BJk, obtain the real product VJk X BJ,k in configuration space by simple multiplication, and then transform this product back to Fourier space. This process introduces aliasing errors, but they can be simply removed [6]. The resulting algorithm retains all the convergence and energy conservation properties of fully spectral methods; indeed, properly dealiased pseudospectral methods are mathematically equivalent to spectral methods [12]. Since the solutions of eqs. (4) are not periodic in the radial coordinate r, the techniques discussed above are not directly applicable. Similar problems are encountered in viscous hydrodynamic problems in nonperiodic domains, such as flow in pipes or rectangular channels. In these cases approximation in the nonperiodic coordinate is accomplished by expansion in Chebyshev polynomials [12], which have the desirable property that their zeros (collocation points) are densely spaced near the outer boundary giving a naturally concentrated mesh there. This is important since hydrodynamic boundary layers form in this region. In resistive MHD, the important boundary layers are associated with filamentation of current, occur internally, and may form spontaneously at nonpredetermined locations. Thus a uniform mesh, or one that is internally dense, is desirable. Expansion in Bessel functions, the natural functions for analytic investigation, is flawed by nonuniform convergence at the outer boundary and lack of a fast transform algorithm [12]. Thus we choose finite differences as the most efficient representation of the radial coordinate,
21
2.3. Temporal approximation: the semi-implicit method As discussed in section 1, the relevant dynamics that determine the evolution of the RFP evolve on time scales that are long compared to those associated with the fastest normal modes of the system. Stability restrictions placed on explicit temporal approximations can result in uneconomically small time steps, while the fully implicit treatment of nonlinear terms either requires iteration (as in finite difference methods [14]), or results in unacceptably large matrices (as in spectral methods). Recently, a new class of methods for solving the time-dependent MHD equations based on an algorithm developed for long range weather simulation [15,16] has been introduced [17,18]. These semi-implicit algorithms allow very large time steps, yet avoid the complexity and large memory requirements associated with implicit methods. In the semi-implicit method for MHD new terms are added to the time-discretized equations that do not affect the consistency of the solution, yet provide a simple and efficient means of enhancing stability. This technique was first applied to the MHD equations in order to eliminate the fast mode time step constraint [17]. Subsequently, this procedure was extended to eliminate the shear Alfvén time step constraint as well [18]. This method is unconditionally stable with respect to all Alfvén modes and consequently permits such large time steps that accuracy becomes the most important consideration in the choice of step size. The semi-implicit method may be implemented using a variety of different time advance schemes and forms for the semi-implicit terms. The choice of the particular semi-implicit method, both in the order of accuracy and the form of the semi-implicit terms, can affect the size of the allowable time step. The proper choices can lead to improved efficiency and accuracy for long-time scale three-dimensional nonlinear MHD computations. The essence of the semi-implicit method can be illustrated by considering the simple first order equation [16] a~a F 5 —.
/
t
—
i ~a
where iw is a frequency or, more generally, an
22
D, D. Schnack et at.
/
Simulation
operator that determines the normal modes of the system. A simple time difference is introduced on the left hand side of (5), i.e., aF,/at (F” ~ F”)/~t. Evaluating the right hand side of (5) at the old time step (F”) results in an explicit method; evaluation at the advanced time step (F”~) yields an implicit method. In the semi-implicit method the frequency (operator) iw is split as i~ i(a + /3), with part evaluated at the new time level and part at the old, i.e., =
=
F”~—F” __________
i(aF”~’+ $F”).
=
(6)
The algorithm introduces no damping, and is unconditionally stable provided a2> ,32• In practice, this may be accomplished by setting /3 w a, with the resulting algorithm =
F”~ F” —
_________
—
,
=
i~aF”+ iaF”~
—
ictF”
(7)
unconditionally stable provided a> w/2. Thus a term iaF is added and subtracted from the right hand side of (5), with one treated at the new time level and the other at the old; the original frequency (operator) is always treated at the old time level, The MHD equations (1) consist of two coupled equations of the form (5). For heuristic purposes, we consider the simple linear system [7] au,/at=Bv, (8a) av/at Au, (8h) =
where A and B are linear operators independent of time, and u and v are vectors. Eqs. (8a, b) can be combined to yield a wave equation for v
a2v/at2
=
Lv,
(9)
where L AB. An identical equation holds for u if A and B commute. Eqs. (8a, b) are analogous to (la, b), and eq. (9) is analogous to the linearized MHD wave equation =
a2 v/at2
=
of
reversed-fieldpinch dynamics
where G is an operator yet to be defined. Combining (ha, b) we find, analogous to (9), Lv, (12) (1 a z~tG)a2v/at which is the wave equation for the modified systern. It can be shown [7] that if the semi-implicit operator G is chosen such that, for the long wavelength modes of interest, its eigenvalues and eigenvectors closely approximate those of L. then the solutions of (12) are also a close approximation to those of (9). Furthermore, the coefficient a appearing in (lib) and (12) can he chosen to yield an unconditionally stable algorithm [7,17.18]. The power of the semi-implicit method lies in choosing G such that it is easily inverted, and in choosing a to yield unconditional stability. In =
—
MHD a particularly useful form for G is obtained from the linearized operator LV=c’xv~’x(VxB 0)xB0 (13) by replacing B0 with a constant vector C0, and then retaining only terms proportional to where a r, 0, z [18]. This yields a simple block tn-diagonal system that is easily inverted. If a modified leap frog algorithm is used for the cxplicit terms [6], unconditional stability then results provided 2/4 1 (14) S~B~k,~,,~t is simultaneously satisfied for a r, 0, z. This serves to determine the coefficients C~.When the right hand side of (14) is negative the algorithm is explicitly stable and C~is set to zero. The semi-implicit method is a simple, fast and =
—
=
accurate algorithm for three-dimensional MHD simulations. It is simple because it is extremely easy to implement, even in existing codes. It is fast because the computation time per time step is dominated by the evaluation of nonlinear convolutions. Since the semi-implicit method requires
ç~x ç~’x (V>< B
0) x B0.
(10)
only the solution of linear tridiagonal systems, it
In applying the semi-implicit method, we introduce the modified system
requires essentially the same CPU time per time step as explicit methods, while allowing time steps that may be orders of magnitude larger. It is accurate because the semi-implicit operator can be chosen to mimic the MHD operator for modes of
au/at By, av,/at Au + a ~1tGav/at, =
=
(ha) (lIb)
D.D. Schnack et al.
/
Simulation of reversed-field pinch dynamics
23
io~
interest. The abillity to simulate driven systems for times that are significant fractions of the resistive diffusion time is essential to understanding RFP dynamics.
io-~ 1 o-~ ‘-I
C.,
3. Examples of RFP dynamics
3.1. Sustained field reversal, t ~
~ T~
[7] c.z~
The essence of the RFP is the sustainment of field reversal and toroidal flux in the presence of an applied toroidal voltage for times exceeding the resistive diffusion time of the system. We simulate this effect by posing an initial value problem in which the initial conditions consist of a one-dimensional equilibrium characterized by the safety factor profile (is) q(r) 0.4(1 1.8748r2 + 0.83232r4), =
—
which we perturb with both the (m 1, n —1) and (m 1, n —2) modes. We sustain the current by applying an axial electric field at the wall that is constant in space and time. The average poloidal electric field at the wall is set to zero so as to conserve toroidal flux. The resistivity is constant in space and time, with Lundquist numben S iO~.With unit aspect ratio, the spatial mesh is 65 radial, 8 poloidal and 16 axial points, This case has been the subject of a previous study [19], where its nonlinear evolution was computed for up to 200 Alfvén times (0.2 resistive diffusion times) with several nonlinear 3-D MHD codes. There it was found that the nonlinear growth and saturation of the unstable nonresonant (1, —2) mode led to the generation and sustainment of an average reversed axial field (Br) at the wall, as is observed in reversed-field pinch (RFP) experiments. These simulations indicated the attainment of a steady field-reversed state dominated by an m 1, n —2 helical perturbation that was maintained against resistive diffusion by finite flow. Using the semi-implicit method we have cornputed the nonlinear evolution of this case for times in excess of one resistive diffusion time (> 1000 Alfvén times). In fig. 1 we plot the energy in various modes as a function of time (measured =
=
LAB
M
N
LAB
M
N
1013
01 32
0 1 02 1 -1 1-2
5 6 87
1 -4 2-1 2 .2 2-3
1015
4
1
9
2
.2
.4
-3
.6
.8
-4
1.0
t/TR
Fig. 1. Energy in selected (m, n) helical modes as a function of time for a driven RFP simulation. Time is measured in resistive diffusion times (TR). Note the nonlinear transition and modelocking at (‘r O.5Tpj.
=
=
=
=
10h1
=
in resistive times). The initially unstable nonresonant (m 1, n —2) mode grows and saturates in about 0.02 resistive times, producing a helical state in agreement with previous results. However, this state is in turn unstable to the resonant (1, 3) mode, which saturates in about 0.1 resistive times at an amplitude approximately equal to that of the =
=
—
.66 .64
-______________________
__________
.62 .60 .58 56
~ 52 50 .48 .46 .42 .40 1.0 tlTR Fig. 2. Safety factor at r = 0 as a function of time for the driven RFP simulation shown in fig. 1. ‘
.2.4.6.8
24
D. D. Schnack ci at.
/
Simulation of reversed-field pinch dynamics
(1, —2). The (1, —2) and (1, —3) then proceed to periodically exchange energy for 0.4 resistive times, resulting in a time-dependent, field-reversed state. At t O.STR a further nonlinear transition occurs, with the (1. —3) and (1, —4) modes inter___________________________________ 200~ ‘i’ l7OkA
changing roles. For t O.7T~ a fully three-dimensional resistive steady state results. The nonlinear periodic energy exchange between the (1, —2) and (1, —3) modes occurring for O.ISTR ~ t ~ O.ST~ can be understood by considering q(0), the value of the safety factor at r 0. In fig. 2 this quantity is plotted as a func=
a
150
RAMP
100
60 kA 50
,I,,,’fl
0 0.15
__
‘‘‘~‘‘~‘‘‘l’
2.0
b
I
I
170 kA
0.10
I~IIb
LA~
-J U-
~kA~
0.05
H
0 I
_____________________________
0.08
I
I
I
-
C
0.06
RAMP
0.04 60 kA~
0.02
.15 25 .05
0 -0.02 _______________________________________
0
5
10 TIME (ms)
15
-
20
Fig. 3. RFP ramp shot experimental data taken from ref. [2]: (a) toroidal current versus time; (b) average toroidal flux versus time; (c) toroidal magnetic field at the wall (vacuum liner) versus time,
0.0
.2
.4 TIME (TR)
Fig. 4. Numerical simulation of RFP ramp shot: (a) toroidal current versus time; (Is) toroidal flux versus time; (c) field reversal parameter F versus time,
D.D. Schnack et a!.
/ Simulation
of reversed-fieldpinch dynamics
tion of time. Initially, q(0) <0.5, so that the unstable (1, —2) mode is nonresonant from above, As is well known [3], the growth of such a mode raises q(0) (“second reconnection”), increasing the shear, and stabilizing the mode. This behavior is seen at early times (t < O.1TR) in fig. 2. However, this profile modification is sufficient to destabilize the resonant (1, —3) mode, whose evolution is such as to lower q(0), i.e., to remove the unstable
25
resonance (“first reconnection”) [3]. This is seen in fig. 2 for O.1TR ~ t ~ O.2T~.The achievement of q(0) <0.5 again destabilizes the (1, —2), causing q(0) to increase, which in turn destabilizes the (1, 3), causing q(0) to decrease. This periodic oscillation about q(0) 0.5 is apparent in fig. 2 for O.2Tp, ~ I ~ O.STR. Finally, at t O.STR the ifl terchange of the (1, —3) and (1, —4) causes q(0) to increase abruptly, followed by the attainment —
=
t=10
t=20 ~“
“•:‘~~‘~
‘~~“
:~
,~
•
20
f
20
‘~
16 12
I
,,~
8
•
0.0
.20
.~ ••~
z
,
12
,J”
.60
.80
1,00
‘5’”
~“
0.0
.20
.40
t=30
‘~
1
20~
1.00
24 20
16~
I
c
16 ‘~‘“•‘‘~ ‘.
8
0 0.0
.80
28
24
12
.60
t=40
28
I
•
—r——.
—r—-—+
~
~
,,
,i.,.
8
.40
f
‘~
16
z
I
. B,,
,‘
‘•
..,‘
‘:
,
“
,‘
, ,.
12
‘~
8
~
~
.20
.40
.60
—r—-—
.80
1.00
0 0.0
:,
‘
,,‘
‘
‘~
‘:
“‘
‘‘:..
‘
‘
,
—‘
.20
.40
.60
.80
1.00
—r———~
Fig. 5. Field line plots (surfaces of section) in the (r, z) plane at various times for resistive mode evolution in a low 9 RFP discharge. Note that good flux surfaces remain near the wall (r = 1).
D.D.
26
/ Simulation
Schnack ci a!.
of a resistive steady state. This second nonlinear transition and mode-locking is presently being studied. 3.2. Simulation of ramp-shots [10]
As mentioned in the introduction, one of the most interesting reversed-field pinch (RFP) cxperimental results occurs in discharges where the toroidal current is slowly raised after toroidal field reversal is established. These discharges exhibit the RFP dynamo effect in its most unambiguous
of reversed-field pinch dynamics
form showing the strong nonlinear coupling of the external poloidal and toroidal magnetic field circuits via the plasma—magnetic-field configuration. In this section we present the results of the application of the numerical methods described in section 2 to this problem. In fig. 3 ZT-40M data from ref. [2] are shown. Parts a—c give the toroidal current, average toroidal flux and the toroidal magnetic field at the wall (vacuum liner), respectively, as a function of time for three different discharges; two have relatively flat toroidal current at 60 and 170 kA the third is
t=10
t=20
28
~‘
I’,
•,i,’
,‘
‘ ‘..~
‘‘‘‘
28
,
—
24
I,
‘C’
, ‘
‘‘
‘
, ‘
•:~
‘
‘~
.,,,~,s;i
24 “
“•
‘‘4’’,,
,
20
.
‘‘
.
‘.
‘..
,
20
__I
_
0.0
.20
.40
.60
t=30
.80
1,00
0.0
.20
.40
‘
.
.“
,
~.
•‘:
.
~‘,‘
,,
28
24
24
20
20
16
16
‘
‘
,
.1.
.
‘. ,
,
.~
~‘
1’ ‘
‘‘.‘.
~
~
‘,
,.
~.
*5
I]
..
.‘,,
S
12
I
8
. ,
1.00
.
•
12
.80
t=45_________________
_______
‘.5 5,
28
.60
,
8
,,
,•,,,,
,
,
•‘.~‘.
,‘~:,
, •.‘
.~
~,
~
.:‘.
4 1 .80 1.00 0,0 .20 .40 .60 .80 1.00 Fig. 6. Field line plots (surfaces of section) in the (r, z) plane at various times for resistive mode evolution in a high 8 RFP discharge. 0.0
.20
.40
.60
Note that outer flux surfaces are destroyed, and then reform,
D.D. Schnack et a!. / Simulation of reversed-field pinch dynamics
a “ramp” discharge where the current goes between these two values in about 15 ms. Notice the deeper reversal of the ramp case in part c; enhanced MHD fluctuations are clearly visible in both parts b and c of fig. 3 for the rising currçnt case. In fig. 4 we display the results of the numerical simulation of such a ramped current discharge [10]. The initial conditions consisted of a one-dimensional equilibrium with magnetic perturbations added for modes m = 1 and n = —2 to —5 inclusive. The Lundquist number S is 2300. The applied axial voltage at the wall is linearly increased from its initial value with a doubling time of 0.5 TR. The poloidal voltage is adjusted to keep the pinch parameter 0 approximately constant at its initial value of 1.58. Fig. 4a shows the toroidal current I4~~ fig. 4b the toroidal flux TF, and fig. 4c the field reversal parameter F, as functions of time (measured in units of TR). Notice that F remains negative and approximately constant throughout the run, even though the toroidal flux has increased significantly, indicating that the plasma has produced positive and negative flux in equal amounts and expelled the excess negative flux to the external circuit. The modes responsible for this enhanced dynamo effect have been identified as long-wavelength m = 1 helical resistive modes. Thus significant RFP physics lies at long wavelength. Further details and discussion of these results can be found in ref. [10]. 3.3. Sawtooth oscillations Experimental observations have demonstrated that long-wavelength, low-frequency signals are the dominant structures in ZT-40M [20,21]. For 0 1.55 the discharge exhibits large-amplitude sawtooth oscillations in the soft X-ray flux, The sawtooth crash is associated with large m = 0 and 1 magnetic perturbations measured at the wall [20]. Prior to the crash, m = 1 precursor oscillations are observed in the soft X-ray flux and sometimes in the m = 1 edge magnetic probe. For 0 1.55, the sawtooth oscillations are smaller, the m = 0 magnetic perturbation becomes imperceptible, and a continuous m = 1 oscillations appears.
27
Previous numerical simulations [8] have demonstrated that the m = 1 modal activity such as that described in section 3.1, leads to destruction of flux surfaces and the presence of stochastic field lines in the core of the plasma (this is shown in fig. 5). These m = 1 modes couple nonlinearly to m = 0 modes that are resonant at the field reversal surface. At large 0, these driven m = 0 modes grow large enough to overlap with the m 1 modes active in the plasma core. (The m = 0 islands can become large because, at high 9, the field reversal surface is relatively farther from the outer boundary.) The discharge then rapidly becomes completely stochastic (see fig. 6), with an implied subsequent rapid outward transport of thermal energy. At low 0 (fig. 5) the stochastic region resulting from m = 1 interactions is limited to the core, with good flux surfaces, and hence confinement, remaining at the edge. This picture of “edge confinement” is consistent with the cxperimentally observed sensitivity of the RFP to field errors. Toroidal flux generation (dynamo) is experimentally observed to be associated with the sawtooth crash [19]; similar flux generation has been observed in the numerical simulations at high 0 [8]. At low 0, confinement is maintained and a quasicontinuous level of m = 1 turbulence remains. Thus, the following picture of high 0 RFP sawtooth oscillations emerges: transport processes cause temperature and current to gradually peak on-axis, thus lowering q(0). When the q profile becomes sufficiently flattened, several m = 1 resistive modes are destabilized. Their nonlinear interaction causes a stochastic core to grow. Simultaneously, m = 0 magnetic islands are driven at the reversal surface. When these overlap with the m = 1 modes the entire discharge becomes stochastic causing rapid heat transport to the wall, a lowering of core temperature, flattening of current, raising of q(0), and a subsequent stabilization of the m = 1 activity. The process then repeats.
4. Summary and conclusions A picture of the RFP as a driven, nonlinear, dynamical, three-dimensional system has been
28
D.D. Schnack er at. /
Simulation of reversed-field pinch dynamics
given. Numerical methods appropriate for the accurate simulation of these long time scale processes have been presented. We have presented examples of the application of these methods to essential experimentally observed processes: long time driven sustainment, ramped current operation and sawtooth oscillations. These calculations have served to demonstrate that long wavelength, lowfrequency helical instabilities are suffictent to cxplain these observed characteristics of the RFP. Acknowledgements The authors wish to thank Drs. B. Carraras, J, Holmes, A. Mirin, W. Matthaeus, D. Montgomery, C. Rowdyshrub, L. Turner and J. Wesson for advice, collaboration, discussion and criticism throughout all phases of the work reviewed in this paper. This work was performed under contract DE-ACO3-83ER53150 with the US Department of Enerov References [1] J.B. Taylor, Phys. Rev. Lett. 33 (1974) 139. [2] J.A. Phillips et al., Los Alamos National Laboratory Report LA-10060-MS (May 1984). [3] E.J. Caramana, R.A. Nebel and D.D. Schnack. Phys. Fluids 26 (1983) 1305. [4] E.J. Caramana and D.A. Baker, Nucl. Fusion 24 (1984) 423.
[5] A.Y. Aydemir and D.B. Barnes, Phys. Rev. Lett. 52 (1984) 930. [6] D.D. Schnack, D.C. Baxter and E.J. Caramana, J. Cornput. Phys. 55(1984) 485. [7] D.D. Schnack, D.C. Barnes, Z. Mikic, D.S. Harned and E.J. Caramana, SAIC Report SAIC-86/3022-APPAT-83; J. Comput. Phys. (1986) submitted. [8] D.D. Schnack, E.J. Caramana and R.A. Nebel, Phys. Fluids 28 (1985) 321. [9] K.A. Werley, R.A. Nebel and GA. Wurden, Phys. Fluids 28 (1985) 1450 [101E.J. Caramana and D.D. Schnack, Los Alamos National Laboratory Report LA-UR-86-965, Phys. Fluids (1986) submitted. [11] D.G. Fox and S.A. Orszag, J. Comput. Phys. 11 (1973) 612. [12] D. Gottlieb and S.A. Orszag, Numerical Analysis of Spectral Methods: Theory and Application (SIAM. Philadelphia, 1977). [13] HR. Hicks, BA. Carraras, J.A. Holmes, D.K. Lee and B V. Waddell, J. Comput. Phys. 44 (1981) 46 [14] C.H. Finan III and 1. Killeen. Comput. Phys. Commun. 24 (1981) 441. [15] A.J. Robert, in: Proc. WMO/IUGG Symp. on Numencal Prediction, Tokyo Methods (1969). Used in Atmospheric [16] Weather A.J. Robert, in: Numerical Models. WMO GARP Publication Series, no. 17. vol. 2 (1979). [17] D.S. Harned and W. Kerner, J. Comput. Phys. 60 (1985) 62. [181 D.S. Harned and D.D. Schnack. J. Comput. Phys. 65 (1986) 57. [19] A.Y. Aydemir, D.C. Barnes, El. Caramana, A.A. Mirin, R.A. Nebel, D.D. Schnack and AG. Sgro, Phys. Fluids 28 (1985) 898. [20] R.G. Watt and R.A. Nebel, Phys. Fluids 26 (1983) 1168. [21] GA. Wurden, Phys. Fluids 27 (1984) 551.