Stabilizing Effects of Edge Current Density on ... - Semantic Scholar
UW-CPTC 11-14
Stabilizing Effects of Edge Current Density on Peeling-Ballooning Instability P. Zhu, C. C. Hegna, and C. R. Sovinec Department of Engineering Physics and Department of Physics University of Wisconsin-Madison Madison, WI 53706
Abstract Resistive MHD computations using the NIMROD code have found strong dependence of the low-n edge instabilities on edge current density distribution. The low-n edge localized modes can be driven unstable by increasing the edge current density across the peeling-ballooning stability boundary. When the edge peak current density is sufficiently large, the safety factor q-profile develops a region where the magnetic shear becomes zero or negative. In these cases, the lown peeling dominant edge instabilities are partially or fully stabilized. The stabilizing effects of zero or reversed magnetic shear are consistent with analytical theory predictions on the necessary condition for the stability of the peeling mode. Nonlinear simulations indicate that the stabilizing effects of edge current density on the low-n edge instabilities through zero or reverse shear can persist throughout the nonlinear exponential growth phase. Near the end of this nonlinear phase, the filament size in radial direction can well exceed the pedestal width, and disconnected bloblike substructures start to develop within the filaments. Relative pedestal energy loss from these radially extending filaments can reach 20%. Both filament size and pedestal energy loss from the nonlinear low-n edge instabilities can be reduced and regulated by the equilibrium edge current density distribution.
1
I.
INTRODUCTION
Empirical scaling of energy loss due to type-I edge localized modes (ELMs) with collisionality suggests the possible role of edge bootstrap current in pedestal instability [1]. The onset of type-I ELMs is commonly interpreted as the consequence of breaching the marginal stability boundary for the coupled peeling-ballooning modes in pedestal region [2–4]. The peeling (ballooning) component of the edge mode prevails in the low-n (high-n) region of the spectrum, which is primarily driven by the edge parallel current density (pressure gradient). Here n is the toroidal mode number. It is generally speculated that the bootstrap current would have a more dominant presence in the low collisionality regime of edge pedestal region and provide a stronger drive for the peeling component of edge instability. Recently, bloblike structures have been observed to form in the pressure profile during the late nonlinear development stage of peeling-ballooning instability of edge pedestal in MHD simulations [5] using the NIMROD code [6]. These simulation studies suggest that the nonlinear ELM properties may depend upon whether the edge current or the edge pressure gradient is the dominant drive for the edge instability. As an initial step towards understanding the roles of edge current density on ELM dynamics, we need to establish how the edge localized current distribution can affect the characteristics of peeling-ballooning instability. Previous analyses of peeling instability indicate that the linear growth rate of the peeling mode would increase with edge parallel current density [7, 8]. Those analyses mainly consider an edge configuration where the safety factor q profile is monotonically increasing towards the boundary and the local magnetic shear is positive throughout the pedestal region. However, in the edge pedestal region of tokamaks, a bootstrap current can often introduce a local peak in edge current density profile and result in a region with zero and/or reversed magnetic shear. Such a change in the q profile can have a considerable effect on the pedestal stability. Configurations with optimized magnetic shear, either reversed or weak, is known to stabilize curvature driven MHD modes in the core. For high-n ideal MHD ballooning mode, stability boundary calculations as denoted by s − α diagrams show that the gap in critical pressure gradient denoting the first and the second marginal stability regimes becomes smaller as magnetic shear s decreases [9]. For sufficiently small magnetic shear, the second stability regime of ideal MHD ballooning mode can be accessed with a moderate increase
2
in pressure gradient. On the other hand, for a given pressure gradient, reducing magnetic shear to small or even negative (shear reversal) would completely stabilize the high-n ideal MHD-ballooning mode [10], allowing strong peaking of pressure gradient in that region. Global stability analysis carried out for the proposed tokamak TPX found that negative central shear configuration increases the β limit [11]. An n = 1 ideal MHD mode known as the “infernal mode” [12], did develop in an enhanced reversed shear (ERS) discharge on TFTR before terminating in disruption, as predicted by calculation and observed in experiment [13]. In the ideal MHD ballooning stability analyses for the simulation of a proposed second stable core VH-mode (SSC-VH) regime with high-β negative central shear on DIIID [14], as well as for a high performance reversed central shear experiment on DIII-D [15], most of the plasma bulk except the narrow region at edge were found to be in the second stability regime. Stability calculations for both TPX and DIII-D found that low-n ideal MHD kink modes are less stable in negative shear configuration with no conducting wall boundary condition [11, 14] due to more peaked pressure profile and less magnetic shear near the edge than in the configuration with monotonic q profile. In the presence of a conducting vacuum vessel, however, more stabilization is introduced to the reversed central shear case than to the conventional shear case, because of the shifting of the mode center, where q = qmin or s = (r/q)(dq/dr) = 0, closer towards the conducting wall. The computed β limit due to the low-n kink mode for the TPX configuration with a conducting wall boundary condition is higher in the reversed shear case than in the conventional shear case. Calculations for the proposed SSC-VH regime on DIII-D also found that the ideal n = 1 to n = 3 kink modes stable in negative central shear configuration with conducting wall, while retaining high β. Besides high-n ideal MHD ballooning mode and low n ideal MHD kink mode, the toroidal Alfv´en eigenmode (TAE) for the hollow current profile is also reported to be less unstable than that for the monotonic current profile [16]. Earlier analyses of edge stability indicated that the low n kink modes and high n peeling modes can restrict access to the second stability regime of ballooning modes for equilibriums with weak magnetic shear and magnetic wells [3, 17]. Later eigenmode analysis of peelingballooning stability in edge pedestal demontrated that an intermediate range of magnetic shear exists that allows easy access to the second stability regime when the magnetic well of the equilibrium is sufficiently deep [18]. Recently, it was suggested the reversed magnetic 3
shear region in the edge pedestal region due to the bootstrap current can give rise to infernal modes as well, which may explain the observation of edge harmonic oscillations (EHOs) in DIII-D experiments [19]. In this work, we intend to assess the effects of edge current density profile, particularly those of weak and reversed magnetic shear, on peeling-ballooning instability in initial value calculations of edge instability using the NIMROD code. For the following calculations, the TOQ solver [20] is used to construct a range of tokamak equilibra with circular shape and different edge localized current profiles, indicative of bootstrap current profiles present in the tokamak pedestals. It is found from linear NIMROD calculations based on these configurations that the low n components of the linear instability spectrum are sensitive to the edge current density distribution, indicating the dominantly current-driven nature of the instability at low n. Contrary to what would be expected from previous peeling mode theory, the growth rates of the low n components of the linear peeling-ballooning instability actually decreases with increasing edge current density when the peak value of the edge current density becomes sufficiently large. Associated with high edge current density profiles are non-monotonic q profiles with zero and reversed magnetic shear regions located in the pedestal. This leads us to investigate the roles of zero and reversed magnetic shear as a possible explanation for the stabilizing effects of edge current density on the low-n peeling-ballooning instability. The rest of the paper is organized as follows. In Sec. II the preparation and the properties of tokamak equilibriums with different edge-localized current profiles are described. We then present the findings on the edge current effects from linear calculations in Sec. III, where a possible interpretation of the results using a simplified analytical theory is given. The edge current effects on the nonlinear development of the low-n peeling dominant perturbation are briefly investigated in Sec. IV. Finally the main findings are summarized and discussed in Sec. V.
II.
EFFECTS OF EDGE CURRENT DENSITY ON EQUILIBRIUM
In order to focus on the effects of edge profiles, a series of model tokamak equilibria with a simple circular-shaped plasma boundary are constructed with a common pressure profile and different edge parallel current density profiles, using the TOQ solver [20]. A hyperbolic tangent function p(x) = ph tanh [(xp − x)/Lp ] + pp for pressure profile is used to model the 4
steep pedestal region near the edge of an H-mode plasma (Fig. 1). Here, x =
q
Ψp /2π, Ψp is
the normalized poloidal magnetic flux, xp the location of pedestal center, pp the pressure at xp , and ph half the height of the pressure pedestal. For a fixed pressure profile as in Fig. 1 with xp ' 0.72, Lp ' 0.018, pp ' 5094, and ph = 5122, a set of profiles for parallel current density are prescribed with different edge peak values in the pedestal region to reflect various contributions from bootstrap currents (Fig. 2). The total model current density combines contributions from both external current drive and bootstrap current, and is used in an iterative procedure for equilibrium solution as described in Ref. [21]. The flux surface averaged product of the external driven current density JCD and equilibrium magnetic field B is a parameterized function < JCD · B >= cj (1 − xν1 )ν2 . Here the constants cj , ν1 , and ν2 are adjusted to set up the desired current profile. In addition, the calculation of the bootstrap current density is based on a simplified model derived from theory in Ref. [22] that was previously used in Refs. [21, 23], and a constant factor cb is used to scale the overall contribution from the bootstrap current density. For the set of current profiles in Fig. 2, the contributions from external current drive are the same, with cj = 1, ν1 = 1, and ν2 = 1.125; each case has a different contribution from the bootstrap current with a corresponding value of cb , which ranges from 0.1 to 5. Each set of pressure and current density profiles, together with the plasma boundary shape, uniquely determine a tokamak equilibrium configuration. The corresponding safety factor q profiles are mostly similar in the core region (x < 0.7). In the pedestal region, a local minimum in q profile, accompanied with reversed magnetic shear, becomes more prominant as the peak value of current density increases (Fig. 2). Reversed magnetic shear is often observed in experiments at locations of internal transport barriers as a consequence from the combined effects of special current ramping and bootstrap current contribution. Transport simulations of the hybrid and steady state advanced operating modes in ITER have shown that reversed magnetic shear can routinely form in the edge pedestal region due to peaked edge current density profiles in typical parameter regimes and operation scenarios [24]. To study the effects of the edge current density profile on the peeling-ballooning instability, initial value linear calculations are employed to follow the development of perturbation with finite toroidal mode number n in these configurations. The results for these calculations are described in next section.
5
III.
LINEAR CALCULATION RESULTS
The finite element mesh used in NIMROD computation conforms to the equilibrium flux surface (Fig. 3). Since we are primarily interested in the peeling-ballooning instabilities with low to intermediate toroidal mode numbers, our computations use a resistive MHD model without the inclusion of finite Larmor radius effects. A Spitzer resistivity model is used in our study. The temperature of the plasma core (edge) for the equilibrium profiles shown in Figs. (1)-(2) is chosen to be 1673eV (83eV), which corresponds to a Lundquist number of 108 (106 ). For each of the edge current density profiles, the linear growth rates of the MHD modes are obtained from NIMROD calculations for toroidal numbers n = 1 − 12 (Fig. 4). The linear growth rates of low-n (n < 10) peeling-ballooning modes are sensitive to the edge current density. As the peak value of the edge parallel current density increases slightly from Case cb = 0.1 to Case cb = 1, the linear growth rates for the n = 3 − 8 modes show moderate reduction. When the peak edge current density increases from Case cb = 1 to Case cb = 2, the toroidal modes with n = 5 − 12 have a significant decrease in linear growth. Further increases of the peak edge current density from Case cb = 2 to Case cb = 5 only result in a small drop of growth rates in the range of toroidal numbers n ∼ 4 − 7. The linear growth rates of those modes shown in Fig. 4 have reached the desirable numerical convergence. For example, the growth rates in the cb = 1 case have been calculated with increasing spatial and temporal resolutions (Fig. 5). In comparison to the base case with a 5th order polynomial finite-element basis functions and a maximum time step of 5 × 10−8 , which are used to obtain the results in Fig. 4, the growth rates of the modes at higher resolutions stay at almost the same values. Energy principle analysis indicates that the edge current density primarily contributes to the drive for the peeling mode component of the edge localized modes [7]. For a given pressure gradient and toroidal mode number, the edge localized modes can be driven unstable by increasing the edge current density across the peeling-ballooning stability boundary. What’s surprising however is that the growth rates of the low n modes can actually decrease with the peak edge current density, when the magnetic shear becomes zero or negative in the same edge location where the peak edge current density increases, as shown in our calculations. Such a stabilizing effect of edge current may be inferred from the analytic theory on the necessary stability criterion for the localized peeling mode [3, 7].
6
In analogous to the Mercier criterion for localized internal modes, the MHD energy principle for peeling modes was constructed for localized surface modes to determine the necessary condition for stability, which requires [3, 7] q
1 − 4DM
2 I J·B >1+ dl, 2πq 0 R2 Bp3
(1)
where DM is the Mercier index [25, 26], R the major radius, q the safety factor, q 0 = dq/dr, J the total equilibrium current density, Bp the poloidal component of the equilibrium magnetic field B, and dl is the poloidal arc length element on a flux surface. All quantities are evaluated at the plasma surface. To leading order for circular shaped tokamaks in the limit of large aspect ratio, the necessary stability condition for peeling mode in Eq. (1) can be written explicitly as [3] (
1 r 1− 2 α R q
!
Rs + s∆ − ft 2r 0
)
> Rqs
JCD · B B2
(2)
where α = −2µ0 Rq 2 p0 /B 2 , r the minor radius, s = rq 0 /q the magnetic shear, B the magnitude of the equilibrium magnetic field B, p the equilibrium pressure, p0 = dp/dr, and JCD the equilibrium current density from inductive and non-inductive external current drives only, excluding the contribution from the Pfirsch-Schl¨ uter and bootstrap currents. The poloidal plane shaping effects were ignored in this derivation. The Shafranov shift ∆0 in (2) comes from the Pfirsch-Schl¨ uter current contribution, and the term proportional to the trapped particle fraction ft is due to the bootstrap current. Although the condition in (2) was mainly obtained for tokamak configurations with a finite positive magnetic shear s, one may infer from Eq. (2) that decreasing s can reduce the peeling-mode drives from both external and bootstrap currents. Further more, when the magnetic shear is reversed (i.e. s < 0), Eq. (2) suggests that increasing either external or bootstrap current can have a stabilizing effect on the peeling mode. If we interpret the plasma surface considered in the analytical theory on peeling mode as the center of the narrow pedestal region in the model equilibrium used in our numerical calculations, the notion that the zero or reversed shear may be responsible for the stabilizing effect of increasing edge current on edge instability found in our calculations appears to be consistent with the peeling mode theory. Our initial-value calculation indicates that very low n-modes (n ≤ 5) can be completely stabilized by edge current density through zero or reversed magnetic shear. For the same edge curent condition, growth rates of modes with 6 ≤ n ≤ 11 are substantially reduced 7
with increasing edge current. A natural question to ask is whether the stabilizing effects of edge current on these peeling-dominant low-n modes in the range of (6 ≤ n ≤ 11) can persist when they evolve into the nonlinear stage. To address this question, we investigate the nonlinear development of the low-n peeling-ballooning instability in one of the typical equilibria with peaked edge current density in the next section.
IV.
NONLINEAR SIMULATION RESULTS
We use NIMROD simulations to follow the nonlinear evolution of an initial perturbation with toroidal mode number n = 6 in the cb = 1 case of equilibrium configuration, which has an edge region with reversed magnetic shear. The initial perturbation for the nonlinear simulation is dominated by an n = 6 Fourier component, and the simulation includes the dynamics and interactions of 43 Fourier components in the toroidal direction. The growth of the total kinetic energy are mainly driven by that of the n = 6 Fourier component until the late nonlinear stage when the contributions to total kinetic energies from other nonlinearly driven Fourier components become comparable to that from the n = 6 component (Fig. 6). Up to this nonlinear stage, the nonlinear growth rate is similar to and basically determined by the growth rate of the corresponding linear phase. Such a nonlinear phase may be called the nonlinear exponential growth phase. The stabilizing effects of edge current density through zero/reversed shear persist in this nonlinear exponential growth phase of the low-n peeling-dominant modes. The nonlinear exponential growth phase of high-n ballooning instability was previously predicted in theory and confirmed in simulations [27–29]. During this nonlinear phase, a high-n ballooning instability continues to grow at the rate of its linear phase, and the radial √ extent of the ballooning filament reaches the scale of Leq / n, where Leq is the equilibrium scale length and n is the toroidal mode number. Thus in comparison to the higher-n modes, the radial extent of the filament grown from the ballooning mode at a lower n is more prominent in the nonlinear exponential growth phase. Such a property of the nonlinear ballooning instability mostly derives from the toroidal nature of the ballooning mode. For that reason, although the peeling mode is driven by the edge current instead of pressure gradient, the peeling-dominant, low n peeling-ballooning instability may share similar features with the ballooning instability during the early nonlinear phase, such as the exponential growth and 8
the filament size scaling, as suggested by the simulation results reported here. By the end of the nonlinear exponential growth phase, the Rayleigh-Taylor finger or filamentary strucutures of the peeling-ballooning instabilities have extended radially to an amplitude around the outboard midplane about 10% to 20% of the tokamak plasma minor radius, as shown in the total pressure contours at t = 2.30 × 10−4 s (Figs. 7). Radially disconnected blob-like segments start to form within those radially extending filaments that are mostly located in the outboard side of the poloidal plane. The linear growth rate can be a reasonable indicator for the radial extent of the peeling-ballooning filaments during this nonlinear phase. The nonlinear consequences of the mode development can also be quantified by measuring the relative pedestal energy loss. For convenience of definition, we divide the physical domain into three regions in radial direction core
p ≥ ptop ped ,
region:
top
out
pedestal region: pped ≤ p < pped , outer region: p < pout , ped
out where p is the equilibrium pressure, ptop ped the pressure at the top of pedestal, and pped the
pressure of pedestal at the boundary between the pedestal and the outer regions which would normally described by a magnetic separatrix in a divertor tokamak. There is no magnetic separatrix in the model limited tokamak equilibriums considered here, and the value of pout ped top is given by pout ped = pped /2 as an example. The internal energy of plasma in each of the
three regions defined above are measured as a function of time in the nonlinear simulations here for the two equilibrium edge current cases. The relative changes of the internal energy in each region from its initial equilibrium value are plotted in Fig. 8. At the end of the nonlinear exponential growth phase around t = 2.30 × 10−4 s, the pedestal energy loss is about 20%. Due to the persistence of the stabilization effects from edge current density on the nonlinear exponential growth phase, the pedestal energy loss rate, up to the end of the nonlinear exponential growth phase, may be reduced or modulated by the equilibrium edge peak current density through formation of regions with non-monotonic magnetic shear in the pedestal.
9
V.
SUMMARY AND DISCUSSION
In summary, resistive MHD computations using NIMROD code find strong dependence of the low-n peeling dominant edge instabilities on the edge current density distribution. Increasing edge current density does not always increase the linear growth rate of low-n peeling-dominant edge-localized-modes. When the edge peak current density is sufficiently large, the corresponding magnetic shear in the pedestal region becomes zero or negative, and the low-n peeling dominant edge instabilities are partially or fully stabilized. The stabilizing effects of zero or reversed magnetic shear can be inferred from previous analytic theory on the necessary condition for the peeling mode stability [3, 7]. Nonlinear simulations indicate that the stabilizing effects of edge current density on the low-n peeling-dominant modes through zero/reversed shear can persist throughout the nonlinear exponential growth phase. Near the end of this nonlinear phase, the filament size in the radial direction can well exceed the pedestal width, and disconnected blob-like substructures start to develop within the filaments. Relative pedestal energy loss from these radially extending filaments can reach 20%. Both filament size and pedestal energy loss from the nonlinear low-n peeling-dominant instabilities can be reduced and regulated by the equilibrium edge current density distribution. Ideal MHD instabilities known as “infernal modes” can develop in zero and reversed shear region [12, 13]. The relation between the infernal modes and the partially stabilized low-n peeling-dominant edge instabilities in our simulations is not clear. In addition, a circular-shaped, limiter-like tokamak configuration has been used in this work. How the non-circular shape factors such as triangularity and elongation, as well as the separatrix geometry, would affect the stabilization effects of edge current density remains unknown. Finally, a full saturation phase is needed for an evaluation of the pedestal energy loss after each pedestal crash. Additional physics mechanism, numerical schemes, and computational resources are required to meet the challenge of simulating the full pedestal crash process before its recovery. These important issues are beyond the scope of this work and will be addressed in future studies.
10
Acknowledgments
This research is supported by U.S. Department of Energy under Grant No. DE-FG0286ER53218 and No. DE-FC02-08ER54975. This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.
[1] A. Loarte, G. Saibene, R. Sartori, D. Campbell, M. Becoulet, L. Horton, T. Eich, A. Herrmann, G. Matthews, N. Asakura, et al., Plasma Physics and Controlled Fusion 45, 1549 (2003). [2] C. C. Hegna, J. W. Connor, R. J. Hastie, and H. R. Wilson, Phys. Plasmas 3, 584 (1996). [3] J. W. Connor, R. J. Hastie, H. R. Wilson, and R. L. Miller, Phys. Plasmas 5, 2687 (1998). [4] P. B. Snyder, H. R. Wilson, J. R. Ferron, L. L. Lao, A. W. Leonard, T. H. Osborne, A. D. Turnbull, D. Mossessian, M. Murakami, and X. Q. Xu, Phys. Plasmas 9, 2037 (2002). [5] B. Burke, in 52nd Annual Meeting of the Division of Plasma Physics, American Physical Society, November 8-12, 2010, Chicago, IL (American Physical Society, 2010), vol. 55, BAPS.2010.DPP.GI2.3. [6] C. Sovinec, A. Glasser, D. Barnes, T. Gianakon, R. Nebel, S. Kruger, D. Schnack, S. Plimpton, A. Tarditi, M. Chu, et al., J. Comput.Phys. 195, 355 (2004). [7] D. Lortz, Nucl. Fusion 15, 49 (1975). [8] J. Manickam, Phys. Fluids B 4, 1901 (1992). [9] D. Lortz and J. N¨ uhrenberg, Physics Letters 68A, 49 (1978). [10] J. M. Greene and M. S. Chance, Nucl. Fusion 21, 453 (1981). [11] C. Kessel, J. Manickam, G. Rewoldt, and W. M. Tang, Phys. Rev. Lett. 72, 1212 (1994). [12] J. Manickam, N. Pomphrey, and A. M. M. Todd, Nuclear Fusion 27, 1461 (1987). [13] F. Levinton, M. Zarnstorff, S. Batha, M. Bell, R. Bell, R. Budny, C. Bush, Z. Chang, E. Fredrickson, A. Janos, et al., Phys. Rev. Lett. 75, 4417 (1995). [14] A. D. Turnbull, T. S. Taylor, Y. R. Lin-Liu, and H. S. John, Physical Review Letters 74, 718 (1995). [15] E. Strait, L. Lao, M. Mauel, B. W. Rice, T. S. Taylor, K. H. Burrell, M. S. Chu, E. A. Lazarus, T. H. Osborne, S. J. Thompson, et al., Phys. Rev. Lett. 75, 4421 (1995).
11
[16] T. Ozeki, M. Azumi, Y. Ishii, Y. Kishimoto, G. Y. Fu, T. Fujita, G. Rewoldt, M. Kikuchi, Y. Kamada, H. Kimura, et al., Plasma Physics and Controlled Fusion 39, A371 (1997). [17] G. T. A. Huysmans, C. D. Challis, M. Erba, W. Kerner, and V. V. Parail, in Proceedings of the 22nd European Physical Society Conference on Controlled Fusion and Plasma Physics, (Bournemouth, UK, 1995) (European Physical Society, Geneva, 1995), part I, p. 201. [18] H. R. Wilson and R. L. Miller, Phys. Plasmas 6, 873 (1999). [19] L.-J. Zheng, in 2011 International Sherwood Fusion Theory Conference, May 2-4, 2011, Austin, TX (2011). [20] R. L. Miller and J. W. Van Dam, Nucl. Fusion 27, 2101 (1987). [21] R. L. Miller, Y. R. Lin-Liu, A. D. Turnbull, V. S. Chan, L. D. Pearlstein, O. Sauter, and L. Villard, Phys. Plasmas 4, 1062 (1997). [22] S. P. Hirshman, Phys. Fluids 31, 3150 (1988). [23] Y. R. Lin-Liu, A. D. Turnbull, M. S. Chu, J. R. Ferron, R. L. Miller, and T. S. Taylor, Phys. Plasmas 6, 3934 (1999). [24] C. E. Kessel, G. Giruzzi, A. Sips, R. V. Budny, J. F. Artaud, V. Basiuk, F. Imbeaux, E. Joffrin, M. Schneider, M. Murakami, et al., Nucl. Fusion 47, 1274 (2007). [25] C. Mercier, Nucl. Fusion 1, 47 (1960). [26] C. Mercier, Nucl. Fusion Suppl. 2, 801 (1962). [27] P. Zhu and C. C. Hegna, Phys. Plasmas 15, 092306 (2008). [28] P. Zhu, C. C. Hegna, and C. R. Sovinec, Phys. Rev. Lett. 102, 235003 (2009). [29] P. Zhu, C. C. Hegna, C. R. Sovinec, A. Bhattacharjee, and K. Germaschewski, Nucl. Fusion 49, 095009 (2009).
12
p
10000
5000
0 5 4 q
3 2 1
<J||*B>
0 1.5e+06 1e+06 5e+05 0 0
0.2
0.4 0.6 1/2 (Ψp/2π)
0.8
1
Figure 1: Radial profiles for equilibrium pressure (upper panel), safety factor (middle panel), and parallel current density (with a factor from equilibrium magnetic field magnitude) (lower panel) as a function of
q
Ψp /2π, where Ψp is normalized poloidal flux function.
13
<J*B>
3e+06
cb=0.1 cb=1 cb=2 cb=3 cb=5
2e+06
1e+06
0
q
3 2 1 0
0
0.1
0.2
0.3
0.4
0.5 0.6 1/2 (Ψp/2π)
0.7
0.8
0.9
1
Figure 2: Radial profiles from a set of equilibria (as labeled with different values of the factor cb ) for parallel current density (with a factor from equilibrium magnetic field magnitude) (upper panel) and safety factor (lower panel) as a function of
q
Ψp /2π, where Ψp is normalized poloidal
flux function. All equilibria have a common equilibrium pressure profile shown in Fig. 1 (upper panel).
14
-1.5 -1.0 -0.5 0.0 0.5 1.0 1.5
Z
1
2
3
4
5
R Figure 3: The finite element mesh used in the NIMROD computation where the circular lines conform to the equilibrium flux surfaces obtained from the TOQ solver.
15
0.08 cb=0.1 cb=1 cb=2 cb=3 cb=5
γτA
0.06
0.04
0.02
0 0
1
2
3
4
7 8 5 6 toroidal mode number n
9
10
11
12
Figure 4: Linear growth rates (normalized with Alfv´enic time τA = 0.44µs) of edge localized instabilities with toroidal mode number n = 1 − 12 from the set of equilibria (as labeled with different values of factor cb ) for the different edge parallel current density and safety factor profiles shown in Fig. 2.
16
0.08 poly_deg=5, dtm=5e-8 poly_deg=6, dtm=5e-8 poly_deg=7, dtm=5e-8 poly_deg=8, dtm=5e-8 poly_deg=8, dtm=1e-8
γτA
0.06
0.04
0.02
0 0
1
2
3
4
7 8 5 6 toroidal mode number n
9
10
11
12
Figure 5: Linear growth rates of edge localized instabilities with toroidal mode number n = 1 − 12 from a particular equilibrium (Case cb = 1) obtained from NIMROD calculations with different spatial and temporal resolutions. Here poly deg refers to the order of polynomials for finite elements in poloidal plane; dtm refers to the maximum time step in unit of second for time advance.
17
1e+05 1 1e-05 total
kinetic energy
1e-10
n=6
n=42
1e-15
n=18 n=36
n=12
1e-20
n=0
1e-25
n=30
1e-30 n=24
1e-35 1e-40
0
5e-05
0.0001 time (s)
0.00015
n=0 n=6 n=12 n=18 n=24 n=30 n=36 n=42 total 0.0002
Figure 6: Kinetic energies of the toroidal Fourier components and total kinetic energy as functions of time during the nonlinear development of an initial perturbation that is dominated by an n = 6 Fourier component at t = 0 in the cb = 1 equilibrium case. For simplicity, only the energies for the harmonics of the n = 6 component are plotted.
18
Figure 7: Total pressure contours for the cb = 1 equilibrium case at t = 2.30 × 10−4 s (upper panel). Radially disconnected blob-like segments can be seen from the filamentary structure shown in the zoomed-in plot (lower panel).
19
0.3 core pedestal outer
relative change of internal energy
0.2
0.1
0
-0.1
-0.2
-0.3 0.00019 0.000195
0.0002
0.000205 0.00021 0.000215 0.00022 0.000225 0.00023 time (s)
Figure 8: Relative change of internal energies as functions of time from (and in percentage of) equilibrium values in core, pedestal, and outer regions, respectively, for the cb = 1 equilibrium case.
20
Stabilizing Effects of Edge Current Density on Peeling-Ballooning Instability P. Zhu, C. C. Hegna, and C. R. Sovinec Department of Engineering Physics and Department of Physics University of Wisconsin-Madison Madison, WI 53706
Abstract Resistive MHD computations using the NIMROD code have found strong dependence of the low-n edge instabilities on edge current density distribution. The low-n edge localized modes can be driven unstable by increasing the edge current density across the peeling-ballooning stability boundary. When the edge peak current density is sufficiently large, the safety factor q-profile develops a region where the magnetic shear becomes zero or negative. In these cases, the lown peeling dominant edge instabilities are partially or fully stabilized. The stabilizing effects of zero or reversed magnetic shear are consistent with analytical theory predictions on the necessary condition for the stability of the peeling mode. Nonlinear simulations indicate that the stabilizing effects of edge current density on the low-n edge instabilities through zero or reverse shear can persist throughout the nonlinear exponential growth phase. Near the end of this nonlinear phase, the filament size in radial direction can well exceed the pedestal width, and disconnected bloblike substructures start to develop within the filaments. Relative pedestal energy loss from these radially extending filaments can reach 20%. Both filament size and pedestal energy loss from the nonlinear low-n edge instabilities can be reduced and regulated by the equilibrium edge current density distribution.
1
I.
INTRODUCTION
Empirical scaling of energy loss due to type-I edge localized modes (ELMs) with collisionality suggests the possible role of edge bootstrap current in pedestal instability [1]. The onset of type-I ELMs is commonly interpreted as the consequence of breaching the marginal stability boundary for the coupled peeling-ballooning modes in pedestal region [2–4]. The peeling (ballooning) component of the edge mode prevails in the low-n (high-n) region of the spectrum, which is primarily driven by the edge parallel current density (pressure gradient). Here n is the toroidal mode number. It is generally speculated that the bootstrap current would have a more dominant presence in the low collisionality regime of edge pedestal region and provide a stronger drive for the peeling component of edge instability. Recently, bloblike structures have been observed to form in the pressure profile during the late nonlinear development stage of peeling-ballooning instability of edge pedestal in MHD simulations [5] using the NIMROD code [6]. These simulation studies suggest that the nonlinear ELM properties may depend upon whether the edge current or the edge pressure gradient is the dominant drive for the edge instability. As an initial step towards understanding the roles of edge current density on ELM dynamics, we need to establish how the edge localized current distribution can affect the characteristics of peeling-ballooning instability. Previous analyses of peeling instability indicate that the linear growth rate of the peeling mode would increase with edge parallel current density [7, 8]. Those analyses mainly consider an edge configuration where the safety factor q profile is monotonically increasing towards the boundary and the local magnetic shear is positive throughout the pedestal region. However, in the edge pedestal region of tokamaks, a bootstrap current can often introduce a local peak in edge current density profile and result in a region with zero and/or reversed magnetic shear. Such a change in the q profile can have a considerable effect on the pedestal stability. Configurations with optimized magnetic shear, either reversed or weak, is known to stabilize curvature driven MHD modes in the core. For high-n ideal MHD ballooning mode, stability boundary calculations as denoted by s − α diagrams show that the gap in critical pressure gradient denoting the first and the second marginal stability regimes becomes smaller as magnetic shear s decreases [9]. For sufficiently small magnetic shear, the second stability regime of ideal MHD ballooning mode can be accessed with a moderate increase
2
in pressure gradient. On the other hand, for a given pressure gradient, reducing magnetic shear to small or even negative (shear reversal) would completely stabilize the high-n ideal MHD-ballooning mode [10], allowing strong peaking of pressure gradient in that region. Global stability analysis carried out for the proposed tokamak TPX found that negative central shear configuration increases the β limit [11]. An n = 1 ideal MHD mode known as the “infernal mode” [12], did develop in an enhanced reversed shear (ERS) discharge on TFTR before terminating in disruption, as predicted by calculation and observed in experiment [13]. In the ideal MHD ballooning stability analyses for the simulation of a proposed second stable core VH-mode (SSC-VH) regime with high-β negative central shear on DIIID [14], as well as for a high performance reversed central shear experiment on DIII-D [15], most of the plasma bulk except the narrow region at edge were found to be in the second stability regime. Stability calculations for both TPX and DIII-D found that low-n ideal MHD kink modes are less stable in negative shear configuration with no conducting wall boundary condition [11, 14] due to more peaked pressure profile and less magnetic shear near the edge than in the configuration with monotonic q profile. In the presence of a conducting vacuum vessel, however, more stabilization is introduced to the reversed central shear case than to the conventional shear case, because of the shifting of the mode center, where q = qmin or s = (r/q)(dq/dr) = 0, closer towards the conducting wall. The computed β limit due to the low-n kink mode for the TPX configuration with a conducting wall boundary condition is higher in the reversed shear case than in the conventional shear case. Calculations for the proposed SSC-VH regime on DIII-D also found that the ideal n = 1 to n = 3 kink modes stable in negative central shear configuration with conducting wall, while retaining high β. Besides high-n ideal MHD ballooning mode and low n ideal MHD kink mode, the toroidal Alfv´en eigenmode (TAE) for the hollow current profile is also reported to be less unstable than that for the monotonic current profile [16]. Earlier analyses of edge stability indicated that the low n kink modes and high n peeling modes can restrict access to the second stability regime of ballooning modes for equilibriums with weak magnetic shear and magnetic wells [3, 17]. Later eigenmode analysis of peelingballooning stability in edge pedestal demontrated that an intermediate range of magnetic shear exists that allows easy access to the second stability regime when the magnetic well of the equilibrium is sufficiently deep [18]. Recently, it was suggested the reversed magnetic 3
shear region in the edge pedestal region due to the bootstrap current can give rise to infernal modes as well, which may explain the observation of edge harmonic oscillations (EHOs) in DIII-D experiments [19]. In this work, we intend to assess the effects of edge current density profile, particularly those of weak and reversed magnetic shear, on peeling-ballooning instability in initial value calculations of edge instability using the NIMROD code. For the following calculations, the TOQ solver [20] is used to construct a range of tokamak equilibra with circular shape and different edge localized current profiles, indicative of bootstrap current profiles present in the tokamak pedestals. It is found from linear NIMROD calculations based on these configurations that the low n components of the linear instability spectrum are sensitive to the edge current density distribution, indicating the dominantly current-driven nature of the instability at low n. Contrary to what would be expected from previous peeling mode theory, the growth rates of the low n components of the linear peeling-ballooning instability actually decreases with increasing edge current density when the peak value of the edge current density becomes sufficiently large. Associated with high edge current density profiles are non-monotonic q profiles with zero and reversed magnetic shear regions located in the pedestal. This leads us to investigate the roles of zero and reversed magnetic shear as a possible explanation for the stabilizing effects of edge current density on the low-n peeling-ballooning instability. The rest of the paper is organized as follows. In Sec. II the preparation and the properties of tokamak equilibriums with different edge-localized current profiles are described. We then present the findings on the edge current effects from linear calculations in Sec. III, where a possible interpretation of the results using a simplified analytical theory is given. The edge current effects on the nonlinear development of the low-n peeling dominant perturbation are briefly investigated in Sec. IV. Finally the main findings are summarized and discussed in Sec. V.
II.
EFFECTS OF EDGE CURRENT DENSITY ON EQUILIBRIUM
In order to focus on the effects of edge profiles, a series of model tokamak equilibria with a simple circular-shaped plasma boundary are constructed with a common pressure profile and different edge parallel current density profiles, using the TOQ solver [20]. A hyperbolic tangent function p(x) = ph tanh [(xp − x)/Lp ] + pp for pressure profile is used to model the 4
steep pedestal region near the edge of an H-mode plasma (Fig. 1). Here, x =
q
Ψp /2π, Ψp is
the normalized poloidal magnetic flux, xp the location of pedestal center, pp the pressure at xp , and ph half the height of the pressure pedestal. For a fixed pressure profile as in Fig. 1 with xp ' 0.72, Lp ' 0.018, pp ' 5094, and ph = 5122, a set of profiles for parallel current density are prescribed with different edge peak values in the pedestal region to reflect various contributions from bootstrap currents (Fig. 2). The total model current density combines contributions from both external current drive and bootstrap current, and is used in an iterative procedure for equilibrium solution as described in Ref. [21]. The flux surface averaged product of the external driven current density JCD and equilibrium magnetic field B is a parameterized function < JCD · B >= cj (1 − xν1 )ν2 . Here the constants cj , ν1 , and ν2 are adjusted to set up the desired current profile. In addition, the calculation of the bootstrap current density is based on a simplified model derived from theory in Ref. [22] that was previously used in Refs. [21, 23], and a constant factor cb is used to scale the overall contribution from the bootstrap current density. For the set of current profiles in Fig. 2, the contributions from external current drive are the same, with cj = 1, ν1 = 1, and ν2 = 1.125; each case has a different contribution from the bootstrap current with a corresponding value of cb , which ranges from 0.1 to 5. Each set of pressure and current density profiles, together with the plasma boundary shape, uniquely determine a tokamak equilibrium configuration. The corresponding safety factor q profiles are mostly similar in the core region (x < 0.7). In the pedestal region, a local minimum in q profile, accompanied with reversed magnetic shear, becomes more prominant as the peak value of current density increases (Fig. 2). Reversed magnetic shear is often observed in experiments at locations of internal transport barriers as a consequence from the combined effects of special current ramping and bootstrap current contribution. Transport simulations of the hybrid and steady state advanced operating modes in ITER have shown that reversed magnetic shear can routinely form in the edge pedestal region due to peaked edge current density profiles in typical parameter regimes and operation scenarios [24]. To study the effects of the edge current density profile on the peeling-ballooning instability, initial value linear calculations are employed to follow the development of perturbation with finite toroidal mode number n in these configurations. The results for these calculations are described in next section.
5
III.
LINEAR CALCULATION RESULTS
The finite element mesh used in NIMROD computation conforms to the equilibrium flux surface (Fig. 3). Since we are primarily interested in the peeling-ballooning instabilities with low to intermediate toroidal mode numbers, our computations use a resistive MHD model without the inclusion of finite Larmor radius effects. A Spitzer resistivity model is used in our study. The temperature of the plasma core (edge) for the equilibrium profiles shown in Figs. (1)-(2) is chosen to be 1673eV (83eV), which corresponds to a Lundquist number of 108 (106 ). For each of the edge current density profiles, the linear growth rates of the MHD modes are obtained from NIMROD calculations for toroidal numbers n = 1 − 12 (Fig. 4). The linear growth rates of low-n (n < 10) peeling-ballooning modes are sensitive to the edge current density. As the peak value of the edge parallel current density increases slightly from Case cb = 0.1 to Case cb = 1, the linear growth rates for the n = 3 − 8 modes show moderate reduction. When the peak edge current density increases from Case cb = 1 to Case cb = 2, the toroidal modes with n = 5 − 12 have a significant decrease in linear growth. Further increases of the peak edge current density from Case cb = 2 to Case cb = 5 only result in a small drop of growth rates in the range of toroidal numbers n ∼ 4 − 7. The linear growth rates of those modes shown in Fig. 4 have reached the desirable numerical convergence. For example, the growth rates in the cb = 1 case have been calculated with increasing spatial and temporal resolutions (Fig. 5). In comparison to the base case with a 5th order polynomial finite-element basis functions and a maximum time step of 5 × 10−8 , which are used to obtain the results in Fig. 4, the growth rates of the modes at higher resolutions stay at almost the same values. Energy principle analysis indicates that the edge current density primarily contributes to the drive for the peeling mode component of the edge localized modes [7]. For a given pressure gradient and toroidal mode number, the edge localized modes can be driven unstable by increasing the edge current density across the peeling-ballooning stability boundary. What’s surprising however is that the growth rates of the low n modes can actually decrease with the peak edge current density, when the magnetic shear becomes zero or negative in the same edge location where the peak edge current density increases, as shown in our calculations. Such a stabilizing effect of edge current may be inferred from the analytic theory on the necessary stability criterion for the localized peeling mode [3, 7].
6
In analogous to the Mercier criterion for localized internal modes, the MHD energy principle for peeling modes was constructed for localized surface modes to determine the necessary condition for stability, which requires [3, 7] q
1 − 4DM
2 I J·B >1+ dl, 2πq 0 R2 Bp3
(1)
where DM is the Mercier index [25, 26], R the major radius, q the safety factor, q 0 = dq/dr, J the total equilibrium current density, Bp the poloidal component of the equilibrium magnetic field B, and dl is the poloidal arc length element on a flux surface. All quantities are evaluated at the plasma surface. To leading order for circular shaped tokamaks in the limit of large aspect ratio, the necessary stability condition for peeling mode in Eq. (1) can be written explicitly as [3] (
1 r 1− 2 α R q
!
Rs + s∆ − ft 2r 0
)
> Rqs
JCD · B B2
(2)
where α = −2µ0 Rq 2 p0 /B 2 , r the minor radius, s = rq 0 /q the magnetic shear, B the magnitude of the equilibrium magnetic field B, p the equilibrium pressure, p0 = dp/dr, and JCD the equilibrium current density from inductive and non-inductive external current drives only, excluding the contribution from the Pfirsch-Schl¨ uter and bootstrap currents. The poloidal plane shaping effects were ignored in this derivation. The Shafranov shift ∆0 in (2) comes from the Pfirsch-Schl¨ uter current contribution, and the term proportional to the trapped particle fraction ft is due to the bootstrap current. Although the condition in (2) was mainly obtained for tokamak configurations with a finite positive magnetic shear s, one may infer from Eq. (2) that decreasing s can reduce the peeling-mode drives from both external and bootstrap currents. Further more, when the magnetic shear is reversed (i.e. s < 0), Eq. (2) suggests that increasing either external or bootstrap current can have a stabilizing effect on the peeling mode. If we interpret the plasma surface considered in the analytical theory on peeling mode as the center of the narrow pedestal region in the model equilibrium used in our numerical calculations, the notion that the zero or reversed shear may be responsible for the stabilizing effect of increasing edge current on edge instability found in our calculations appears to be consistent with the peeling mode theory. Our initial-value calculation indicates that very low n-modes (n ≤ 5) can be completely stabilized by edge current density through zero or reversed magnetic shear. For the same edge curent condition, growth rates of modes with 6 ≤ n ≤ 11 are substantially reduced 7
with increasing edge current. A natural question to ask is whether the stabilizing effects of edge current on these peeling-dominant low-n modes in the range of (6 ≤ n ≤ 11) can persist when they evolve into the nonlinear stage. To address this question, we investigate the nonlinear development of the low-n peeling-ballooning instability in one of the typical equilibria with peaked edge current density in the next section.
IV.
NONLINEAR SIMULATION RESULTS
We use NIMROD simulations to follow the nonlinear evolution of an initial perturbation with toroidal mode number n = 6 in the cb = 1 case of equilibrium configuration, which has an edge region with reversed magnetic shear. The initial perturbation for the nonlinear simulation is dominated by an n = 6 Fourier component, and the simulation includes the dynamics and interactions of 43 Fourier components in the toroidal direction. The growth of the total kinetic energy are mainly driven by that of the n = 6 Fourier component until the late nonlinear stage when the contributions to total kinetic energies from other nonlinearly driven Fourier components become comparable to that from the n = 6 component (Fig. 6). Up to this nonlinear stage, the nonlinear growth rate is similar to and basically determined by the growth rate of the corresponding linear phase. Such a nonlinear phase may be called the nonlinear exponential growth phase. The stabilizing effects of edge current density through zero/reversed shear persist in this nonlinear exponential growth phase of the low-n peeling-dominant modes. The nonlinear exponential growth phase of high-n ballooning instability was previously predicted in theory and confirmed in simulations [27–29]. During this nonlinear phase, a high-n ballooning instability continues to grow at the rate of its linear phase, and the radial √ extent of the ballooning filament reaches the scale of Leq / n, where Leq is the equilibrium scale length and n is the toroidal mode number. Thus in comparison to the higher-n modes, the radial extent of the filament grown from the ballooning mode at a lower n is more prominent in the nonlinear exponential growth phase. Such a property of the nonlinear ballooning instability mostly derives from the toroidal nature of the ballooning mode. For that reason, although the peeling mode is driven by the edge current instead of pressure gradient, the peeling-dominant, low n peeling-ballooning instability may share similar features with the ballooning instability during the early nonlinear phase, such as the exponential growth and 8
the filament size scaling, as suggested by the simulation results reported here. By the end of the nonlinear exponential growth phase, the Rayleigh-Taylor finger or filamentary strucutures of the peeling-ballooning instabilities have extended radially to an amplitude around the outboard midplane about 10% to 20% of the tokamak plasma minor radius, as shown in the total pressure contours at t = 2.30 × 10−4 s (Figs. 7). Radially disconnected blob-like segments start to form within those radially extending filaments that are mostly located in the outboard side of the poloidal plane. The linear growth rate can be a reasonable indicator for the radial extent of the peeling-ballooning filaments during this nonlinear phase. The nonlinear consequences of the mode development can also be quantified by measuring the relative pedestal energy loss. For convenience of definition, we divide the physical domain into three regions in radial direction core
p ≥ ptop ped ,
region:
top
out
pedestal region: pped ≤ p < pped , outer region: p < pout , ped
out where p is the equilibrium pressure, ptop ped the pressure at the top of pedestal, and pped the
pressure of pedestal at the boundary between the pedestal and the outer regions which would normally described by a magnetic separatrix in a divertor tokamak. There is no magnetic separatrix in the model limited tokamak equilibriums considered here, and the value of pout ped top is given by pout ped = pped /2 as an example. The internal energy of plasma in each of the
three regions defined above are measured as a function of time in the nonlinear simulations here for the two equilibrium edge current cases. The relative changes of the internal energy in each region from its initial equilibrium value are plotted in Fig. 8. At the end of the nonlinear exponential growth phase around t = 2.30 × 10−4 s, the pedestal energy loss is about 20%. Due to the persistence of the stabilization effects from edge current density on the nonlinear exponential growth phase, the pedestal energy loss rate, up to the end of the nonlinear exponential growth phase, may be reduced or modulated by the equilibrium edge peak current density through formation of regions with non-monotonic magnetic shear in the pedestal.
9
V.
SUMMARY AND DISCUSSION
In summary, resistive MHD computations using NIMROD code find strong dependence of the low-n peeling dominant edge instabilities on the edge current density distribution. Increasing edge current density does not always increase the linear growth rate of low-n peeling-dominant edge-localized-modes. When the edge peak current density is sufficiently large, the corresponding magnetic shear in the pedestal region becomes zero or negative, and the low-n peeling dominant edge instabilities are partially or fully stabilized. The stabilizing effects of zero or reversed magnetic shear can be inferred from previous analytic theory on the necessary condition for the peeling mode stability [3, 7]. Nonlinear simulations indicate that the stabilizing effects of edge current density on the low-n peeling-dominant modes through zero/reversed shear can persist throughout the nonlinear exponential growth phase. Near the end of this nonlinear phase, the filament size in the radial direction can well exceed the pedestal width, and disconnected blob-like substructures start to develop within the filaments. Relative pedestal energy loss from these radially extending filaments can reach 20%. Both filament size and pedestal energy loss from the nonlinear low-n peeling-dominant instabilities can be reduced and regulated by the equilibrium edge current density distribution. Ideal MHD instabilities known as “infernal modes” can develop in zero and reversed shear region [12, 13]. The relation between the infernal modes and the partially stabilized low-n peeling-dominant edge instabilities in our simulations is not clear. In addition, a circular-shaped, limiter-like tokamak configuration has been used in this work. How the non-circular shape factors such as triangularity and elongation, as well as the separatrix geometry, would affect the stabilization effects of edge current density remains unknown. Finally, a full saturation phase is needed for an evaluation of the pedestal energy loss after each pedestal crash. Additional physics mechanism, numerical schemes, and computational resources are required to meet the challenge of simulating the full pedestal crash process before its recovery. These important issues are beyond the scope of this work and will be addressed in future studies.
10
Acknowledgments
This research is supported by U.S. Department of Energy under Grant No. DE-FG0286ER53218 and No. DE-FC02-08ER54975. This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.
[1] A. Loarte, G. Saibene, R. Sartori, D. Campbell, M. Becoulet, L. Horton, T. Eich, A. Herrmann, G. Matthews, N. Asakura, et al., Plasma Physics and Controlled Fusion 45, 1549 (2003). [2] C. C. Hegna, J. W. Connor, R. J. Hastie, and H. R. Wilson, Phys. Plasmas 3, 584 (1996). [3] J. W. Connor, R. J. Hastie, H. R. Wilson, and R. L. Miller, Phys. Plasmas 5, 2687 (1998). [4] P. B. Snyder, H. R. Wilson, J. R. Ferron, L. L. Lao, A. W. Leonard, T. H. Osborne, A. D. Turnbull, D. Mossessian, M. Murakami, and X. Q. Xu, Phys. Plasmas 9, 2037 (2002). [5] B. Burke, in 52nd Annual Meeting of the Division of Plasma Physics, American Physical Society, November 8-12, 2010, Chicago, IL (American Physical Society, 2010), vol. 55, BAPS.2010.DPP.GI2.3. [6] C. Sovinec, A. Glasser, D. Barnes, T. Gianakon, R. Nebel, S. Kruger, D. Schnack, S. Plimpton, A. Tarditi, M. Chu, et al., J. Comput.Phys. 195, 355 (2004). [7] D. Lortz, Nucl. Fusion 15, 49 (1975). [8] J. Manickam, Phys. Fluids B 4, 1901 (1992). [9] D. Lortz and J. N¨ uhrenberg, Physics Letters 68A, 49 (1978). [10] J. M. Greene and M. S. Chance, Nucl. Fusion 21, 453 (1981). [11] C. Kessel, J. Manickam, G. Rewoldt, and W. M. Tang, Phys. Rev. Lett. 72, 1212 (1994). [12] J. Manickam, N. Pomphrey, and A. M. M. Todd, Nuclear Fusion 27, 1461 (1987). [13] F. Levinton, M. Zarnstorff, S. Batha, M. Bell, R. Bell, R. Budny, C. Bush, Z. Chang, E. Fredrickson, A. Janos, et al., Phys. Rev. Lett. 75, 4417 (1995). [14] A. D. Turnbull, T. S. Taylor, Y. R. Lin-Liu, and H. S. John, Physical Review Letters 74, 718 (1995). [15] E. Strait, L. Lao, M. Mauel, B. W. Rice, T. S. Taylor, K. H. Burrell, M. S. Chu, E. A. Lazarus, T. H. Osborne, S. J. Thompson, et al., Phys. Rev. Lett. 75, 4421 (1995).
11
[16] T. Ozeki, M. Azumi, Y. Ishii, Y. Kishimoto, G. Y. Fu, T. Fujita, G. Rewoldt, M. Kikuchi, Y. Kamada, H. Kimura, et al., Plasma Physics and Controlled Fusion 39, A371 (1997). [17] G. T. A. Huysmans, C. D. Challis, M. Erba, W. Kerner, and V. V. Parail, in Proceedings of the 22nd European Physical Society Conference on Controlled Fusion and Plasma Physics, (Bournemouth, UK, 1995) (European Physical Society, Geneva, 1995), part I, p. 201. [18] H. R. Wilson and R. L. Miller, Phys. Plasmas 6, 873 (1999). [19] L.-J. Zheng, in 2011 International Sherwood Fusion Theory Conference, May 2-4, 2011, Austin, TX (2011). [20] R. L. Miller and J. W. Van Dam, Nucl. Fusion 27, 2101 (1987). [21] R. L. Miller, Y. R. Lin-Liu, A. D. Turnbull, V. S. Chan, L. D. Pearlstein, O. Sauter, and L. Villard, Phys. Plasmas 4, 1062 (1997). [22] S. P. Hirshman, Phys. Fluids 31, 3150 (1988). [23] Y. R. Lin-Liu, A. D. Turnbull, M. S. Chu, J. R. Ferron, R. L. Miller, and T. S. Taylor, Phys. Plasmas 6, 3934 (1999). [24] C. E. Kessel, G. Giruzzi, A. Sips, R. V. Budny, J. F. Artaud, V. Basiuk, F. Imbeaux, E. Joffrin, M. Schneider, M. Murakami, et al., Nucl. Fusion 47, 1274 (2007). [25] C. Mercier, Nucl. Fusion 1, 47 (1960). [26] C. Mercier, Nucl. Fusion Suppl. 2, 801 (1962). [27] P. Zhu and C. C. Hegna, Phys. Plasmas 15, 092306 (2008). [28] P. Zhu, C. C. Hegna, and C. R. Sovinec, Phys. Rev. Lett. 102, 235003 (2009). [29] P. Zhu, C. C. Hegna, C. R. Sovinec, A. Bhattacharjee, and K. Germaschewski, Nucl. Fusion 49, 095009 (2009).
12
p
10000
5000
0 5 4 q
3 2 1
<J||*B>
0 1.5e+06 1e+06 5e+05 0 0
0.2
0.4 0.6 1/2 (Ψp/2π)
0.8
1
Figure 1: Radial profiles for equilibrium pressure (upper panel), safety factor (middle panel), and parallel current density (with a factor from equilibrium magnetic field magnitude) (lower panel) as a function of
q
Ψp /2π, where Ψp is normalized poloidal flux function.
13
<J*B>
3e+06
cb=0.1 cb=1 cb=2 cb=3 cb=5
2e+06
1e+06
0
q
3 2 1 0
0
0.1
0.2
0.3
0.4
0.5 0.6 1/2 (Ψp/2π)
0.7
0.8
0.9
1
Figure 2: Radial profiles from a set of equilibria (as labeled with different values of the factor cb ) for parallel current density (with a factor from equilibrium magnetic field magnitude) (upper panel) and safety factor (lower panel) as a function of
q
Ψp /2π, where Ψp is normalized poloidal
flux function. All equilibria have a common equilibrium pressure profile shown in Fig. 1 (upper panel).
14
-1.5 -1.0 -0.5 0.0 0.5 1.0 1.5
Z
1
2
3
4
5
R Figure 3: The finite element mesh used in the NIMROD computation where the circular lines conform to the equilibrium flux surfaces obtained from the TOQ solver.
15
0.08 cb=0.1 cb=1 cb=2 cb=3 cb=5
γτA
0.06
0.04
0.02
0 0
1
2
3
4
7 8 5 6 toroidal mode number n
9
10
11
12
Figure 4: Linear growth rates (normalized with Alfv´enic time τA = 0.44µs) of edge localized instabilities with toroidal mode number n = 1 − 12 from the set of equilibria (as labeled with different values of factor cb ) for the different edge parallel current density and safety factor profiles shown in Fig. 2.
16
0.08 poly_deg=5, dtm=5e-8 poly_deg=6, dtm=5e-8 poly_deg=7, dtm=5e-8 poly_deg=8, dtm=5e-8 poly_deg=8, dtm=1e-8
γτA
0.06
0.04
0.02
0 0
1
2
3
4
7 8 5 6 toroidal mode number n
9
10
11
12
Figure 5: Linear growth rates of edge localized instabilities with toroidal mode number n = 1 − 12 from a particular equilibrium (Case cb = 1) obtained from NIMROD calculations with different spatial and temporal resolutions. Here poly deg refers to the order of polynomials for finite elements in poloidal plane; dtm refers to the maximum time step in unit of second for time advance.
17
1e+05 1 1e-05 total
kinetic energy
1e-10
n=6
n=42
1e-15
n=18 n=36
n=12
1e-20
n=0
1e-25
n=30
1e-30 n=24
1e-35 1e-40
0
5e-05
0.0001 time (s)
0.00015
n=0 n=6 n=12 n=18 n=24 n=30 n=36 n=42 total 0.0002
Figure 6: Kinetic energies of the toroidal Fourier components and total kinetic energy as functions of time during the nonlinear development of an initial perturbation that is dominated by an n = 6 Fourier component at t = 0 in the cb = 1 equilibrium case. For simplicity, only the energies for the harmonics of the n = 6 component are plotted.
18
Figure 7: Total pressure contours for the cb = 1 equilibrium case at t = 2.30 × 10−4 s (upper panel). Radially disconnected blob-like segments can be seen from the filamentary structure shown in the zoomed-in plot (lower panel).
19
0.3 core pedestal outer
relative change of internal energy
0.2
0.1
0
-0.1
-0.2
-0.3 0.00019 0.000195
0.0002
0.000205 0.00021 0.000215 0.00022 0.000225 0.00023 time (s)
Figure 8: Relative change of internal energies as functions of time from (and in percentage of) equilibrium values in core, pedestal, and outer regions, respectively, for the cb = 1 equilibrium case.
20
