Princeton Plasma Physics Laboratory PPPL-4024 Numerical Study of the Formation, Ion Spin-up and Nonlinear Stability Properties of Field-reversed Configurations E.V. Belova, R.C. Davidson, H. Ji, M. Yamada, C.D. Cothran, M.R. Brown, and M.J. Schaffer November 2004 PRINCETON PLASMA PHYSICS LABORATORY PPPL Prepared for the U.S. Department of Energy under Contract DE-AC02-76CH03073. PPPL Report Disclaimers Full Legal Disclaimer This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, nor any of their contractors, subcontractors or their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or any third party's use or the results of such use of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise, does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof or its contractors or subcontractors. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof. Trademark Disclaimer Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise, does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof or its contractors or subcontractors. PPPL Report Availability This report is posted on the U.S. Department of Energy's Princeton Plasma Physics Laboratory Publications and Reports web site in Fiscal Year 2005. The home page for PPPL Reports and Publications is: http://www.pppl.gov/pub_report/ Office of Scientific and Technical Information (OSTI): Available electronically at: http://www.osti.gov/bridge. Available for a processing fee to U.S. Department of Energy and its contractors, in paper from: U.S. Department of Energy Office of Scientific and Technical Information P.O. Box 62 Oak Ridge, TN 37831-0062 Telephone: (865) 576-8401 Fax: (865) 576-5728 E-mail: [email protected] National Technical Information Service (NTIS): This report is available for sale to the general public from: U.S. Department of Commerce National Technical Information Service 5285 Port Royal Road Springfield, VA 22161 Telephone: (800) 553-6847 Fax: (703) 605-6900 Email: [email protected] Online ordering: http://www.ntis.gov/ordering.htm IC/P6-34 Numerical Study of the Formation, Ion Spin-up and Nonlinear Stability Properties of Field-Reversed Configurations E. V. Belova 1), R. C. Davidson 1), H. Ji 1), M. Yamada 1), C. D. Cothran 2), M. R. Brown 2), M. J. Schaffer 3) 1) Princeton Plasma Physics Laboratory, Princeton NJ, USA 2) Swarthmore College, Swarthmore PA, USA 3) General Atomics, San Diego CA, USA E-mail: [email protected] Abstract. Results of three-dimensional numerical simulations of field-reversed configurations (FRCs) are presented. Emphasis of this work is on the nonlinear evolution of magnetohydrodynamic (MHD) instabilities in kinetic FRCs, and the new FRC formation method by the counter-helicity spheromak merging. Kinetic simulations show nonlinear saturation of the n = 1 tilt mode, where n is the toroidal mode number. The n = 2 and n = 3 rotational modes are observed to grow during the nonlinear phase of the tilt instability due to the ion spin-up in the toroidal direction. The ion toroidal spin-up is shown to be related to the resistive decay of the internal flux, and the resulting loss of particle confinement. Three-dimensional MHD simulations of counter-helicity spheromak merging and FRC formation show good agreement with results from the SSX-FRC experiment. Simulations show formation of an FRC in about 30 Alfvén times for typical experimental parameters. The growth rate of the n = 1 tilt mode is shown to be significantly reduced compared to the MHD growth rate due to the large plasma viscosity and field-line-tying effects. 1. Introduction The field-reversed configuration (FRC) is a compact toroid with little or no toroidal field. It offers a unique fusion reactor potential because of its compact and simple geometry, translation properties, and high plasma beta. At present, the most important issues are FRC stability with respect to low-n (toroidal mode number) MHD modes, and the development of new FRC formation and current drive methods. The traditional theta-pinch formation method usually produces highly kinetic FRCs with relatively low flux, small S* (the FRC kinetic parameter, S*, is the ratio of the separatrix radius to the ion skin depth) and large elongation, E. A theoretical understanding of the observed FRC stability properties has proven to be elusive, due to the complicated interplay of several non-ideal MHD effects [1,2], including finite Larmor radius (FLR) effects, the Hall term, and plasma flow effects. Advanced numerical simulations are required for the self-consistent investigation of the stability properties of kinetic FRCs. Results of such simulations are presented in this paper. A slow FRC formation technique, based on counter-helicity spheromak merging demonstrated the advantage of this approach compared with traditional theta-pinch formation methods [3]. The counter-helicity spheromak merging method allows forming the configuration with large S*, thus permitting experimental studies of large-S* FRC stability properties [3,4]. The SSX-FRC experiment [4] is designed to study FRC formation by the counter-helicity spheromak merging method, and to examine the general issue of FRC stability at large S*. Three-dimensional MHD simulations have been performed in support of the SSX-FRC experiment, and show very good agreement with the experimental results. 2. Study of FRC nonlinear stability properties Numerical studies of the nonlinear evolution of magnetohydrodynamic (MHD) instabilities in kinetic FRCs has been performed using the 3D nonlinear hybrid and MHD simulation code HYM [1]. The stability properties of MHD modes with toroidal mode numbers n ≥ 1 are investigated, including finite ion Larmor radius (FLR) and rotation effects, and finite electron pressure. It has been demonstrated that due to strong FLR stabilization of the higher-n modes, the n = 1 tilt mode is most unstable mode for nearly all experimentally-relevant non-rotating FRC equilibria [1]. An empirical FLR scaling of the tilt mode growth rate has been obtained: γ = C(VA/RsE) exp(−3Eρi/Rs), where γmhd = C(VA/RsE) is the MHD growth rate, E is the separatrix elongation, ρi is ion thermal Larmor radius, and C ≈ 2 is a constant. [Figure 1: Time evolution of the n = 0 - 4 Fourier harmonics of the ion kinetic energy obtained from hybrid simulations with E = 6.25 and S* = 20.] [Figure 2: Contour plots of the ion density perturbation at the FRC cross section of the maximum perturbation amplitude at t = 60tA (|δn/n| ≤ 0.2) and t = 76tA (|δn/n| ≤ 0.4).] Nonlinear kinetic simulations performed for a set of FRC equilibria with E = 4 - 6 and S* = 10 - 80 show that the n = 1 tilt mode saturates nonlinearly without destroying the configuration, provided the FRC kinetic parameter is sufficiently small, S* ≲ 20. In the initial configurations, all of the equilibrium toroidal current is assumed to be carried by the electrons, and the ions have a non-rotating Maxwellian distribution, which is consistent with experimental conditions just after FRC formation. However, as the simulation proceeds, the ions gradually begin to rotate toroidally, and near the end of the simulation run the ion toroidal flow velocity becomes comparable to the ion diamagnetic velocity. Therefore the saturation of the tilt instability occurs in the presence of ion toroidal spin-up, and it is accompanied by the growth of the n = 2 rotational mode, which is often seen in experiments [2]. The saturation of the n = 1 tilt mode and the growth of the n = 2 rotational mode can be seen in Figs. 1 and 2, where the results of nonlinear hybrid simulations are shown for a configuration with E = 6.25 and S* ≈ 20 (for comparison, the parameters for typical theta-pinch FRC experiment are E = 5 - 8 and S* ≲ 20 [2]). A set of 2D (axisymmetric) nonlinear hybrid simulations has been performed in order to study the resistive evolution of the kinetic FRC. Analysis of the simulation results show that the ion toroidal spin-up is related to the resistive decay of the internal flux, and the resulting loss of weakly confined particles. An approximate confinement condition can be written as: p_φ > 0, where p_φ = R v_φ - e/m ψ is the ion canonical toroidal angular momentum, and ψ is the poloidal flux (ψ < 0 inside the separatrix). Weakly confined particles, which have small initial toroidal angular momentum p_φ ≲ |Δψ| will be lost from the closed-field-line region, when the magnitude of the trapped flux ψ0 decreases by Δψ. Since most of the particles with small p_φ have negative toroidal velocity, the particle loss results in the net flux of the negative momentum away from the separatrix region, and therefore there is a net positive ion rotation inside the separatrix (with the ion flow velocity in the current direction). [Figure 3: (a) Time evolution of the normalized angular momentum of all ions, and the ions inside the separatrix; (b) maximum value of the ion toroidal flow velocity obtained from 2D hybrid simulations with S* = 20 and E = 4.] Figures 3a shows the time evolution of the toroidal angular momentum of all ions, and the ions inside the separatrix region. It can be seen that the net angular momentum is approximately conserved due to imposed periodic boundary conditions in z-direction. In contrast, the angular momentum of the part of the plasma confined inside the separatrix is positive and increases in time as the configuration decays. Here the Alfvén time is defined as tA = Rc/VA, where Rc is the radius of the flux conserving shell, and VA is the characteristic Alfvén velocity. The maximum value of the ion flow velocity is plotted in Fig. 3b. The peak value of Vφ is about 0.2-0.3 VA, which is comparable to the ion diamagnetic velocity for S* ≈ 20. In the final state, the ions carry a significant fraction of the total current. Both 2D and 3D simulations with zero initial ion rotation demonstrate the formation of an approximately rigid-rotor rotation profile inside the separatrix in about 40-60 Alfvén times, depending on the plasma resistivity. Poloidal contour plots and radial profiles of ion toroidal flow velocity are shown in Fig. 4. The linear velocity profile (Vφ ~ R) suggests that the distribution function fi of the ions inside the separatrix evolves towards an exponential rigid-rotor distribution function, which is a shifted local Maxwellian distribution. The evolution of fi towards a Maxwellian distribution in the collisionless model in the absence of 3D instabilities indicates that the stochasticity of the ion orbits plays a significant role in the FRC relaxation. [Figure 4a: Contour plots of the ion toroidal velocity in the r - z plane at t = 40 and t = 80tA from 2D hybrid simulations with S* = 20 and E = 4.] [Figure 4b: Radial profiles of the ion toroidal flow velocity at the FRC midplane at t = 20, 40, and 80tA. The separatrix radius is Rs/Rc ≈ 0.6.] Figure 4 shows that the ion flow velocity changes its sign outside the separatrix, and there is significant velocity gradient close to the separatrix at the plasma edge. The n = 2 and n = 3 rotational instabilities are found to be localized in the vicinity of maximum velocity shear near the edge, and have a similar structure to the external modes. Growth rates of the rotational modes are found to be larger in the smaller-S*, more kinetic configurations. The details of the ion toroidal spin-up determine the nonlinear evolution of these instabilities. 3. Counter-helicity spheromak merging simulations The SSX-FRC experiment is designed to study the new FRC formation method by counter-helicity spheromak merging, and the FRC stability properties for large values of S*. In addition, the effects of the residual (axially-antisymmetric) toroidal field on macroscopic stability properties are being studied. Experimental results demonstrate the formation of the FRC-like configuration with S* > 35, and indicate the presence of the global n = 1 instability, consistent with the tilt-mode instability [4]. However, the observed growth rate of this instability is smaller by a factor of 6-8 than that of the ideal MHD growth rate. In addition, the experimental formation studies consistently show the presence of an axially-antisymmetric toroidal field, which does not completely annihilate during the reconnection. The underlying reasons for this are not understood. Numerical simulations using the HYM code have been performed to investigate this issues, and study 3D spheromak merging for experimentally-relevant parameters. 3.1. Axisymmetric simulations Two- and three-dimensional MHD simulations show good agreement with experimental results. Simulations show formation of an FRC in about 30 Alfvén times for typical experimental parameters. Large toroidal and poloidal flows with flow velocity up to ~ 0.5 - 1VA (based on the edge field) are generated during the reconnection phase. The plasma pressure is significantly increased by Ohmic and viscous heating, and the FRC-like pressure profile is formed due to convective transport by poloidal flows [6]. The flow velocity reduces by an order-of-magnitude after the FRC formation is completed. Axisymmetric simulations of counter-helicity spheromak merging have been performed in order to study the dependence of the reconnection rate and the toroidal field annihilation on values of plasma resistivity and viscosity. The SSX-FRC experimental parameters for Ti = 12eV, n = 10^15 cm^-3 and VA = 7 · 10^4 m/s, are as follows: S ≈ 900, so that the resistivity is η ≈ 10^-3, and an estimate for the viscosity coefficient is μ = 2 · 10^-3 - 10^-2. (Here the MHD normalization is used, i.e. η = 1/S and μ = 1/Re, where S is the Lundquist number, and Re is the Reynolds number.) [Figure 5: Time evolution of (a) normalized kinetic energy and (b) magnetic energy obtained from 2D MHD simulations of counter-helicity spheromak merging for μ = 0.001 and several values of resistivity.] [Figure 6: Time evolution of (a) normalized kinetic energy and (b) magnetic energy (total and toroidal) obtained from 2D MHD simulations for η = 0.001 and several values of viscosity.] Figure 5 shows the results of five MHD simulation runs performed for values of η = 2.5 · 10^-4 - 4 · 10^-3 and μ = 10^-3. It can be seen that the time for a complete reconnection depends weakly on the resistivity (this time can be approximately determined by the strong reduction in the kinetic energy). Thus the reconnection is completed and the FRC forms in about 20tA for larger values of η, and in 25-30tA for smaller η. This is consistent with previous theoretical studies of driven magnetic reconnection [5,6]. On the other hand, the details of the time evolution of the plasma kinetic energy (Fig. 5a) and magnetic energy (Fig. 5b) vary strongly with η. For larger values of η, the kinetic energy peaks at t ~ 5tA, and the flow velocity (mostly toroidal) increases with η. This initial increase in the kinetic energy is caused by the force imbalance present in the initial conditions. In smaller resistivity cases (η ≤ 0.002), there is a maximum in the flow energy which occurs at t ~ 20 - 25tA. This peak in the kinetic energy correlates with the large reconnection rate and the fast reduction in the magnetic energy seen in Fig. 5b. The toroidal flow velocity in these cases is comparable to the characteristic Alfvén velocity, Vφ ~ VA. Figure 6 shows the results of a set of MHD simulation runs performed to investigate effects of plasma viscosity for η = 5 · 10^-4. As expected, larger viscosity results in a strong reduction of plasma flows. Reduced flow velocity, on the other hand, is related to a slower reconnection of the toroidal field (Fig. 6b). This effect is more pronounced for smaller values of resistivity. Figure 7 shows contour plots of the toroidal field from the simulation runs with μ = 0.001 and μ = 0.004. It can be seen that large values of plasma viscosity result in incomplete reconnection, and significant residual toroidal fields are present at t ~ 20 - 30tA. Numerical results thus indicate that the large values of plasma viscosity observed in the experiments (μ ≳ 0.002) may be responsible for the incomplete reconnection observed in the SSX-FRC. [Figure 7: Contour plots of toroidal magnetic field at the poloidal plane obtained from 2D MHD simulations with η = 0.001: (a) μ = 0.001, and (b) μ = 0.004. Time is normalized to the Alfvén time.] 3.2. Three-dimensional simulations Three-dimensional MHD simulations have been performed to study the effects of plasma viscosity, self-generated flows, and magnetic field-line-tying effects on the unstable global modes. A random initial perturbation was applied at t = 0. The n = 1 tilt mode was found to be a dominant mode in all cases, and the configuration was tilted at the end of the simulations. Figure 8a shows the time evolution of the n = 1 mode from four simulation runs with η = 10^-3 and μ = 5 · 10^-4 - 4 · 10^-3. The linear growth rate is largest for small viscosity, and it reduces as μ increases. Therefore, large values of viscosity have a strong stabilizing effect on the n = 1 tilt mode, as well as the higher-n modes in the MHD regime. Figure 8b shows the time evolution of the total kinetic energy obtained for the same set of simulations as shown in Fig. 8a. The growth of the n = 1 tilt mode is not seen until t > 20tA, when the amplitude of the mode becomes comparable to that of the n = 0 component of the kinetic energy. Note that by that time the growth rate of the n = 1 mode is reduced compared to its linear value (Fig. 8a). The calculated growth rates obtained from the simulation shown in Fig. 8a (for μ = 0.004) are: γ ≈ 0.33 · VA/Rs for t < 20tA, and γ ≲ 0.15 · VA/Rs for t > 20tA. Several factors may contribute to the reduction of the instability growth at later times (t ≳ 20tA), including the nonlinear mode interactions and an increase of the separatrix elongation. For comparison, an estimate for an ideal MHD growth rate is γ0 ~ (0.7 - 1.3)VA/Rs. Possible stabilizing effect of the toroidal flows, generated during the reconnection process, appears to be less significant than that of viscosity, because the growth rates for smaller-μ cases (and, hence, larger Vφ) are larger, than those for larger μ (smaller Vφ). The finite residual toroidal field, on the other hand, may contribute to the reduction of the growth rate of the n = 1 mode in high-viscosity cases (μ > 10^-3). [Figure 8: Time evolution of (a) n = 1 mode energy, and (b) total kinetic energy from 3D MHD simulations of counter-helicity spheromak merging for μ = 5 · 10^-4 - 4 · 10^-3.] The simulations shown in Fig. 8 have been performed for realistic boundary conditions, including the effects of the line-tying. Another set of simulations have been performed with different boundary conditions, neglecting these effects. A comparison between the runs with different boundary conditions have demonstrated that line tying effects increase the complete reconnection time by 5 - 10tA, and reduce the n = 1 mode growth rate. It has also been found that the peak toroidal flow velocity is reduced by a factor of 2 due to line-tying effects. Simulations show that both the large plasma viscosity and the field-line-tying boundary conditions have strong stabilizing effect on the tilt mode. In addition, there is a nonlinear reduction of the instability drive at the time when the mode amplitude becomes observable. These results provide an explanation for the experimentally measured growth rates that are about 6 - 8 times smaller than the ideal MHD growth rate. 4. Conclusions The hybrid simulations presented here show that, while ion FLR effects determine the linear stability properties of non-rotating FRCs, the inclusion of nonlinear and ion toroidal flow effects is necessary for a satisfactory description of plasma behavior in low-S* FRC experiments. All major experimentally observed stability properties of kinetic (theta-pinch-formed) FRCs have been reproduced and explained with the help of numerical simulations. Namely, the scaling of the linear growth rate of the n = 1 tilt instability with S*/E parameter has been demonstrated for a class of elongated elliptic FRCs [1]; the ion toroidal spin-up, the nonlinear saturation of the tilt mode, and the growth of the n = 2 rotational mode have been demonstrated. It has been shown that the loss of ions with a preferential sign of toroidal velocity due to the resistive decay of the poloidal flux results in the ion toroidal spin-up, which reproduces very well the experimentally observed ion rotation. The time scale of the ion spin-up is determined by the flux decay time. The HYM code has also been used to study FRC formation by the counter-helicity spheromak merging in support of the SSX-FRC experiment [4], and contributed to interpretation of several puzzling experimental observations. In particular, the persistence of the residual toroidal fields and slower-than-MHD growth of the tilt instability have been shown to be related to the large plasma viscosity and the line-tying effects in the SSX-FRC experiments. Acknowledgments This work was supported by DOE contract DE-AC02-76CH03073. Calculations were performed at the U.S. National Energy Research Supercomputing Center. References [1] BELOVA, E.V., et al., Phys. Plasmas 7, 4996 (2000); 11, 2523 (2004). [2] TUSZEWSKI, M., Nucl. Fusion 28, 2033 (1988). [3] ONO, Y., M. INOMOTO, T. OKAZAKI, Y. UEDA, Phys. Plasmas 4, 1953 (1997). [4] COTHRAN, C. D., et al., Phys. Plasmas 10, 1748 (2003). [5] BISKAMP, D., H. WELTER, Phys. Rev. Lett. 44, 1069 (1980). [6] WATANABE, T.-H., T. SATO, T. HAYASHI, Phys. Plasmas 4, 1297 (1997); T. SATO, Y. ODA, S. OTSUKA, Phys. Fluids 26, 3602 (1983). 07/07/03 External Distribution Plasma Research Laboratory, Australian National University, Australia Professor I.R. Jones, Flinders University, Australia Professor João Canalle, Instituto de Fisica DEQ/IF - UERJ, Brazil Mr. Gerson O. Ludwig, Instituto Nacional de Pesquisas, Brazil Dr. P.H. Sakanaka, Instituto Fisica, Brazil The Librarian, Culham Laboratory, England Mrs. S.A. Hutchinson, JET Library, England Professor M.N. Bussac, Ecole Polytechnique, France Librarian, Max-Planck-Institut für Plasmaphysik, Germany Jolan Moldvai, Reports Library, Hungarian Academy of Sciences, Central Research Institute for Physics, Hungary Dr. P. Kaw, Institute for Plasma Research, India Ms. P.J. Pathak, Librarian, Institute for Plasma Research, India Ms. Clelia De Palo, Associazione EURATOM-ENEA, Italy Dr. G. Grosso, Instituto di Fisica del Plasma, Italy Librarian, Naka Fusion Research Establishment, JAERI, Japan Library, Laboratory for Complex Energy Processes, Institute for Advanced Study, Kyoto University, Japan Research Information Center, National Institute for Fusion Science, Japan Dr. O. Mitarai, Kyushu Tokai University, Japan Dr. Jiangang Li, Institute of Plasma Physics, Chinese Academy of Sciences, People's Republic of China Professor Yuping Huo, School of Physical Science and Technology, People's Republic of China Library, Academia Sinica, Institute of Plasma Physics, People's Republic of China Librarian, Institute of Physics, Chinese Academy of Sciences, People's Republic of China Dr. S. Mirnov, TRINITI, Troitsk, Russian Federation, Russia Dr. V.S. Strelkov, Kurchatov Institute, Russian Federation, Russia Professor Peter Lukac, Katedra Fyziky Plazmy MFF UK, Mlynska dolina F-2, Komenskeho Univerzita, SK-842 15 Bratislava, Slovakia Dr. G.S. Lee, Korea Basic Science Institute, South Korea Institute for Plasma Research, University of Maryland, USA Librarian, Fusion Energy Division, Oak Ridge National Laboratory, USA Librarian, Institute of Fusion Studies, University of Texas, USA Librarian, Magnetic Fusion Program, Lawrence Livermore National Laboratory, USA Library, General Atomics, USA Plasma Physics Group, Fusion Energy Research Program, University of California at San Diego, USA Plasma Physics Library, Columbia University, USA Alkesh Punjabi, Center for Fusion Research and Training, Hampton University, USA Dr. W.M. Stacey, Fusion Research Center, Georgia Institute of Technology, USA Dr. John Willis, U.S. Department of Energy, Office of Fusion Energy Sciences, USA Mr. Paul H. Wright, Indianapolis, Indiana, USA The Princeton Plasma Physics Laboratory is operated by Princeton University under contract with the U.S. Department of Energy. Information Services Princeton Plasma Physics Laboratory P.O. Box 451 Princeton, NJ 08543 Phone: 609-243-2750 Fax: 609-243-2751 e-mail: [email protected] Internet Address: http://www.pppl.gov