ABSTRACT Title of Document: STABILIZATION OF MHD TURBULENCE BY APPLIED STEADY AND OSCILLATING VELOCITY SHEAR Ching Pui HUNG, Doctor of Philosophy, 2013 Dissertation Directed By: Professor Adil B. Hassam, Department of Physics Some aspects of velocity shear stabilization of magnetized plasma instabilities are considered. In the first part, steady externally forced flow shears are considered. In the second part, resonantly excited oscillating flow shears are considered. The stabilizing effect of steady forced velocity shear on the ideal interchange instability is studied in linear and nonlinear regimes, with a 2D dissipative magnetohydrodynamic (MHD) code. With increasing flow shear V?, the linearly unstable band in wavenumber-space shrinks so that the peak growth results for modes that correspond to intermediate wavenumbers. In the nonlinear turbulent state, the convection cells are roughly circular on the scale of the density gradient. Unstable modes are almost completely stabilized, with the density profile reverting to laminar, when V? is a few times the classic interchange growth rate. The simulations are compared with measurements of magnetic fluctuations from the Maryland Centrifugal Experiment. The spectral data, taken in the plasma edge, are in general agreement with data obtained in higher viscosity simulations. Finally, concomitant Kelvin-Helmholtz instabilities in the system are also examined. Geodesic acoustic modes (GAMs) are axisymmetric electrostatic poloidal oscillations of plasma in tokamaks. It has been proposed to drive GAMs resonantly by external drivers, thus setting up velocity shears to suppress turbulence. Here, we study enhanced damping of GAMs from (1) phase mixing of oscillations and (2) nonlinear detuning of the resonance. It is well-known that phase mixing of Alfven waves propagating in inhomogeneous media results in enhanced damping. The enhancement goes as the 1/3 power of the dissipation. We study this phenomenon for GAMs in tokamaks with temperature profiles. Our analysis is verified by numerical simulation. In addition, the system of nonlinear GAM equations is shown to resemble the Duffing oscillator. Resonant amplification is shown to be suppressed nonlinearly. The results are applied to the proposed GAM excitation experiment on the DIII-D tokamak. STABILIZATION OF MHD TURBULENCE BY APPLIED STEADY AND OSCILLATING VELOCITY SHEAR by Ching Pui HUNG Dissertation submitted to the Faculty of the Graduate School of the University of Maryland, College Park, in partial fulfillment of the requirements for the degree of Doctor of Philosophy 2013 Advisory Committee: Professor Adil B. Hassam, Chair Professor Thomas Antonsen Professor Richard Ellis Professor Raymond Sedwick Dr. Marc Swisdak ? Copyright by Ching Pui HUNG 2013 ii Acknowledgements I greatly appreciate my advisor, Prof. Adil Hassam, for his patience in the past seven years of guiding me through the research projects. I also thank him for correcting numerous typos and grammatical mistakes in the paper and thesis I wrote. I am also grateful to Dr. P. N. Guzdar for beneficial discussions on my scholarly paper and for being a committee member in my candidacy presentation three years ago, so that I can proceed to the final stage of the PhD program at present. This work was supported by the Department of Energy. iii Table of Contents Acknowledgements....................................................................................................... ii Table of Contents......................................................................................................... iii List of Tables ................................................................................................................ v List of Figures .............................................................................................................. vi Chapter 1: Introduction and Overview ......................................................................... 1 1.1 Introduction......................................................................................................... 1 1.2 Previous work ..................................................................................................... 4 1.3 Overview of the contents .................................................................................... 5 Chapter 2: Residual Turbulence from Velocity Shear Stabilized Interchange instabilities .................................................................................................................... 8 2.1 Introduction......................................................................................................... 8 2.2 Equations: ......................................................................................................... 12 2.2.1 2D MHD Equations ................................................................................... 12 2.2.2 Ordering and Reduced Equations .............................................................. 14 2.2.3 Equilibrium ................................................................................................ 14 2.3 Summary of linear analysis............................................................................... 16 2.3.1 Rayleigh-Taylor instability ........................................................................ 16 2.3.2 Kelvin-Helmholtz instability, Rayleigh inflexion criterion ....................... 18 2.4 Rayleigh-Taylor instabilities with flow shear................................................... 18 2.4.1 Equilibrium profiles and examination of KH stability criterion................ 19 2.4.2 RT instability stabilized by velocity shear in the linear regime................. 23 2.4.3 RT instability with velocity shear in the nonlinear regime........................ 25 2.5 Kelvin-Helmholtz instabilities .......................................................................... 41 2.6 Summary of Chapter 2...................................................................................... 43 Chapter 3: Geodesic acoustic modes .......................................................................... 45 3.1 Introduction....................................................................................................... 45 3.2 Equations .......................................................................................................... 49 Chapter 4: Phase mixing of GAM oscillations ........................................................... 53 4.1 Introduction....................................................................................................... 53 4.2 Initial value problem for linearized GAM equations........................................ 54 4.3 Steady state solution of linearized GAM equations with driving term............. 56 Chapter 5: Nonlinear effects on resonantly driven GAMs ........................................ 64 5.1 Introduction....................................................................................................... 64 5.2 Linear expansion of nonlinear GAM equations with small driver amplitude .. 65 iv 5.3 Analysis of the nonlinear GAM equations at resonance................................... 68 Chapter 6: Implications for experiments to drive GAMs in tokamaks...................... 75 6.1 Introduction....................................................................................................... 75 6.2 Effect of phase mixing on damping rates and power requirements.................. 75 6.3 Comparison of the effect of nonlinearities and the effect of damping ............. 79 6.4 Summary of Chapter 3 to Chapter 6 ................................................................. 81 Chapter 7: SUMMARY and CONCLUSION............................................................ 83 Appendices.................................................................................................................. 86 A. Alfven waves propagating in an inhomogeneous, dissipative medium............. 86 A.1 Equations...................................................................................................... 86 A.2 Numerical simulations for undriven and driven standing Alfven waves..... 90 B. Effects of nonlinearities on driving magnetosonic wave in uniform plasma..... 97 B.1 Equations and expansion in small amplitudes.............................................. 97 B.2 Simulation of the system at resonance ....................................................... 100 C. The Duffing Equation ...................................................................................... 104 Bibliography ............................................................................................................. 107 v List of Tables Table 6.1. Estimated damping rates, Q-factors, spatial spreads and power requirements, by classical collision and phase mixing, at the edge and core of tokamak....................................................................................................................... 77 vi List of Figures FIG. 2.1. Schematic description of equilibrium.......................................................... 15 FIG. 2.2a, b. The 1D equilibrium profiles. Note that the density gradient is opposite to gravity, thus unstable to RT instability. FIG. 2.2c. Cross field flow shear profiles for different strengths characterized by coefficient A in (2.26)........................................ 21 FIG. 2.3a. The density profile can be shifted vertically without largely altering 'n by adjusting the initial average on . FIG. 2.3b. For a parabolic shear flow, the generalized Rayleigh criterion is not guaranteed to be satisfied; the criterion can be satisfied by ?upward shifting? the whole density profile. .............................................................. 22 FIG. 2.4. Rayleigh Taylor unstable growth rate in linear regime, as a function of yk , for parabolic flow shear of various strength parameter A........................................... 24 FIG. 2.5a. Laminar density contour plots when the system is unperturbed. FIG. 2.5b. Turbulent steady state of the system without shear; density profile is inverted to a new quasi-stable equilibrium. FIG. 2.5c, a medium strength shear partially restores the initial profile, with residual fluctuations. FIG. 2.5d, RT convection is completely suppressed by strong flow shear with laminar profile recovered. .............................. 27 FIG. 2.6a-d Time evolution of the spectra along the y-direction with shear strength A as parameter in the nonlinear regime. The random perturbation in density is added at 0t = . Note that the amplitudes at steady state decrease as A is increased from 2. For 4A = (FIG. 2.6d), the system is completely stabilized and the perturbation decreases with time. Also, the nonlinear spectrum peaks at 6ym = , which is the same as that of the linear one............................................................................................................... 28 FIG. 2.7. Density fluctuation (at 0.5x = ) for different shear strengths. .................... 29 FIG. 2.8. Density gradient averaged over all space for different shear strengths....... 29 FIG. 2.9. Comparison of steady state density profiles with laminar density at various shear strengths A. ........................................................................................................ 31 FIG. 2.10. Effect of flow shear on particle flux across the 1/ 2x = axis. The transport is towards negative x direction. At strong shear, classical diffusion accounts for nearly all of the particle transport. ......................................................................................... 33 FIG. 2.11. Linear growth rates as a function of ym , with y? , the viscosity in the y- direction, as a parameter. Shear strength is constant at 2A = . Note the suppression of short wavelength modes when y? is increased. .................................................... 34 FIG. 2.12. Nonlinear residual density fluctuations at A = 2 for a) 42 10y? ?= ? , b) 34 10y? ? = ? , c) 36 10y? ?= ? . Notice the increase in wavelengths and the decrease in fluctuation levels, as the viscosity is increased........................................................... 35 FIG. 2.13. Nonlinear residual convection for A = 2 at a) 42 10y? ?= ? , b) 34 10y? ? = ? , c) 36 10y? ?= ? . Notice the decrease in mode numbers (i.e., elongated convection cells) as the viscosity is increased. ........................................................... 36 FIG. 2.14. Spectral analysis of the residual fluctuations of the y-flow in the y direction at A =2 and low viscosity, which is corresponding to FIG. 2.13a. ............................. 39 vii FIG. 2.15. Spectral analysis of the residual fluctuations of the y-flow in the y direction at A =2 and high viscosity, corresponding to FIG. 2.13c............................................ 40 FIG. 2.16. Contour plot of linear growth rate as a function of on and A, for 5ym = . The contours are labeled by the growth rates. ............................................................ 41 FIG. 2.17. Linear growth rates for strong flow shear in the vicinity of the KH criterion boundary with and without gravity. ............................................................................ 43 FIG. 3.1. Schematic illustration of geodesic acoustic mode in tokamak.................... 46 FIG. 4.1. The profile of GAM frequency [ ( ) 1.2o x x? = + ] and the driving frequency ? along x, which will give a resonance around x = 0. ................................................ 59 FIG. 4.2. The steady state of (normalized) poloidal flow, will the GAM frequency and driver frequency profile indicated in FIG 4.1. Here ( )2 40 ( 0) 10xx L? ? ?= = and 0.1od = . Notice the resonance around x = 0 and the phase difference along x at any fixed time t. ................................................................................................................. 59 FIG. 4.3. Power spectrum measured at x = 0 with ? = 10-4, and ( 0)o x? ? ?? = ? = . ..................................................................................................................................... 60 FIG. 4.4. Logarithmic plot of the spectral width a resonance against the viscosity measured at x = 0, the slope is 0.34 0.02? ................................................................. 61 FIG. 4.5. Logarithmic plot of the amplitude at resonance against the viscosity measured at x = 0, which gives a slope of 0.334 0.002? ? . ....................................... 61 FIG. 4.6. Dependence of width of peak of the power spectrum on the slope of GAM frequency. Viscosity is fixed at ? = 10-4. The slope of the graph is 0.65 0.01? . ....... 61 FIG. 4.7. The amplitude of (normalized) poloidal flow as a function of space at ( 0)o x? ?= = and ? = 10-4, 0.1od = .......................................................................... 62 FIG. 4.8. Logarithmic plot of spatial width of resonance [at ( 0)o x? ?= = ] against the viscosity, with a slope of 0.336 0.001? . .............................................................. 63 FIG. 4.9. Dependence of width of peak of the steady state amplitude on the slope of GAM frequency. Viscosity is fixed at ? = 10-4, and ( 0)o x? ?= = . The slope of the graph is 0.334 0.001? ? . ............................................................................................. 63 FIG. 5.1. The amplitude of u as a function of driving frequency (around the natural GAM frequency) for nonlinear GAM equations (? = 0.2). Notice that the width of resonance is finite and the amplification is bounded even no dissipative terms are included in the equations. ........................................................................................... 71 FIG. 5.2. The power spectrum of nonlinear GAM oscillator. Notice the sharp decrease in amplitude when ?? drops below c?? , a property inherited from Duffing oscillator...................................................................................................................... 72 FIG. 5.3. Logarithmic plot of resonant width against flux surface r at ao = 0.01, the slope is 0.65 0.01? ? ................................................................................................... 73 FIG. 5.4. Logarithmic plot of maximum amplification against flux surface r at ao = 0.01, with a slope of 0.73 0.01? ................................................................................. 73 FIG. 5.5. Logarithmic plot of resonant width against driver amplitude at r = 3, the slope is 0.686 0.003? ................................................................................................. 73 FIG. 5.6. Logarithmic plot of maximum amplification against driver amplitude at r = 3, with a slope of 0.69 0.01? ...................................................................................... 74 viii FIG. A.1. The density and magnetic field profile at equilibrium. .............................. 91 FIG. A.2. The Alfven speed profile corresponds to the profiles in FIG. A.1. ............ 91 FIG. A.3. Time evolution of Alfven wave at the center of the system showing dramatic decay when phase mixing occurs. The corresponding oscillation using a flat profile is shown for comparison. ................................................................................ 92 FIG. A.4. A plot of log ( , )E x t versus 3t , which shows the decay goes at ( )3exp 2 ot t? ??? ? before the consistency condition breaks down at ~ 140( / )x At L v and the curve levels off. The parameters used are 510 [ ( 1/ 2)]x x AL x? ?= =v , 1oB = , 2z zk Lpi= . ................................................................................................................. 94 FIG. A.5, Power spectra measured at the center of the system for driven standing Alfven wave, with and without a profile in Alfven speed. Here ( 1/ 2)o z Ak x? = =v is the local resonant frequency at the center, and 510 [ / ( 1/ 2)]x x AL x? ?= =v , 310 [ / ( 1/ 2)]z x AL x? ?= =v and 510oF ?= are used..................................................... 95 FIG. B.1. Schematic view of driving standing magnetosonic wave in a 2D slab....... 98 FIG. B.2. The schematic diagram showing the breaking down of linear expansion as the driver frequency approaches the resonant frequency.......................................... 100 FIG. B.3. The spectra of driven magnetosonic wave at different driver amplitude ? . Note the flattening of the spectra at a neighborhood of the resonant frequency o k? = ( 0? = ), as a result of the nonlinearities. .................................................................. 102 FIG. B.4. A logarithmic plot of spectral width versus driver amplitude ? . The slope of the straight line is 0.496 0.004? . ......................................................................... 103 FIG. B.5. A logarithmic plot of spectral height versus driver amplitude ? . The slope of the straight line is 0.506 0.001? . ......................................................................... 103 FIG. C.1. A plot of amplification against driving frequency of the Duffing oscillator driven by a sinusoidal force, where * 0.2? = . Note that the peak of resonance shift to higher frequency (compared to the natural frequency), which is opposite to the nonlinear GAM system shown in FIG. 5.1. .............................................................. 106 1 Chapter 1: Introduction and Overview 1.1 Introduction The stabilizing effect of velocity shear on MHD instabilities has been well established. Shear flow can suppress linearly growing instabilities, and can distort convective cells nonlinearly and suppress turbulence. It has been shown that shear flow can stabilize both ideal interchange instabilities and drift instabilities. Experimental observations of magnetized plasmas have also shown the stabilization effect of velocity shear. In general, the stabilization effect is significant when 'V ?> , where 'V is the gradient in the shear flow, i.e., shear strength, and ? is the growth rate of the instability in question. This dissertation addresses topics related to the characteristics of turbulence in the presence of velocity shear and to the efficacy of velocity shear stabilization. The topics are motivated by two separate experimental situations, namely, the Maryland Centrifugal Experiment (MCX) at the University of Maryland, and an experiment proposed for the DIII-D tokamak by Hallatschek and McKee. As such, the dissertation is broadly divided into two parts, corresponding to the experimental situation. In the first part, we study the effects of an externally forced shear flow on macroscopic interchange instabilities, motivated by the MCX. The MCX has shown significant stabilization of interchange modes due to flow shear. The experimental 2 setup resembles a magnetic mirror configuration, namely, a long solenoidal magnetic field with axisymmetric mirror end fields. The plasma is rotated azimuthally by a radially applied electric field [1]. Such a system is well-known to be unstable to interchange modes. However, the velocity shear in the azimuthal flow is thought to stabilize the ideal interchange mode in the system. Although the flow shear in MCX could be strong enough to provide significant stabilization, the stabilizing effect is not necessarily complete: there is residual turbulence, accounting for a significant part of particle and heat transport in the quasi-steady state. This motivates us to conduct a broad study of the strength and characteristics of the residual turbulence in the presence of applied shear. In particular, we focus on the effect of increasing shear flow on the wavenumber spectrum of interchange modes, and its effect on the residual turbulence of partially stabilized interchange instabilities. We apply our findings by comparing our corresponding simulation results with the observed magnetic fluctuations in MCX [2,3]. In addition, the presence of velocity shear can itself lead to Kelvin-Helmholtz (KH) instabilities [4]. In our study, we also take into account the destabilizing tendencies of the KH instabilities. In the second part of the dissertation, we study some feasibility aspects of setting up a time-varying velocity shear in tokamaks to achieve turbulence reduction. This work is motivated by a recently proposed experiment on the DIII-D tokamak. The idea is to resonantly excite Geodesic Acoustic Modes (GAMs) in the tokamak [5], so that the oscillating flow shear can stabilize turbulence. The GAM is an ideal MHD normal mode of tokamak plasma, constituting an axisymmetric poloidal oscillation of flux tubes with well-defined frequency. The frequency of GAM oscillations is known 3 to be ~G sc R? , where sc is the sound speed of the plasma and R is the major radius of the tokamak. GAMs have been experimentally observed in tokamaks [6-8]. The proposed experiment will resonantly excite GAMs with externally pulsed magnetic field coils [5]. To have a significant stabilizing effect, the strength of flow shear should be comparable to the growth rate of the instability, i.e., ' ~ mV u x ?? > , where mu and x? are, respectively, the maximum amplitude of the poloidal flow at the resonating flux surfaces, and the spatial spread of the locally driven GAMs in the radial direction. This gives an estimation of the amplitude of poloidal oscillation required to stabilize the plasma instabilities. The power requirement for driving GAMs depends on the damping rate of the mode. In this dissertation, we study two particular aspects of the resonant driving that may impact power requirements; namely, we study the enhanced damping of GAMs from phase mixing between various magnetic surfaces, and we study the detuning of resonance from nonlinear effects. The effect of phase mixing on waves in inhomogeneous media has been studied for Alfven waves by Heyvaerts and Priest [9]. We extend this idea to GAMs, since the radial temperature gradient in a tokamak also leads to a similar phase mixing effect for driven GAMs. We also show that nonlinear detuning of GAMs in tokamaks is highly similar to nonlinear detuning in the driven Duffing oscillator. We then consider the experimental parameters of the DIII-D tokamak as an example to evaluate the power dissipation by collisional and phase mixed damping, and compare our results with Hallatschek and McKee?s estimates [5]. We also compare the strength of nonlinear terms to the damping terms, for the experimental conditions. 4 1.2 Previous work Theoretical analysis of the effects of flow shear on ideal interchange modes in ordinary fluids was first presented by Kuo (1963) [10]. The stabilization effect of flow shear in plasmas was first experimentally demonstrated by Taylor et al (1989) [11], who applied a radial electric field in the tokamak edge, which resulted in a poloidal rotation and a sharp transport barrier formed at the plasma edge. The edge electric field profiles inferred from the poloidal rotation were compared with that in the transport barriers for edge plasmas in JFT-2M tokamak [12]. The highly localized poloidal rotation profile in the edge plasma was also observed in the DIII-D tokamak [13]. The suppression effects of radial electric field on turbulence and on the transition to better confinement were examined theoretically by Shaing et al (1990) [14]. Linear theory also indicated that flow shear would stabilize drift microturbulence and resistive g-modes [15]. Flow shear was shown to suppress turbulence in tokamaks from distortion of convection cells [16], to stabilize the nonlinear Rayleigh-Taylor instability [17], and to stabilize drift microinstabilities in tokamaks [18]. The above results and other work have been summarized by Groebner (1993) [19] and Burell (1997) [20] in review papers. Numerical simulations have shown that flow shear can also stabilize kink and sausage modes in a Z pinch [21]. It was also shown analytically and numerically that flute interchanges in centrifugally confined plasmas can be stabilized with velocity shear [22]. While resonant excitation of GAMs is a recent proposal, the existence of stable geodesic acoustic normal modes in tokamak plasmas has been well-known and shown analytically in both collisional [23,24] and collisionless [25-28] systems. The GAM 5 frequency was shown to be sonic, or ~G T? , where ( )T T ?= is the temperature of plasma at the flux surface ? . GAMs have also been identified experimentally in tokamaks [6-8]. Excitation of GAMs by resonant magnetic field coils and by resonant heating have been demonstrated numerically [29]. The effect of resonantly excited GAMs on resistive ballooning turbulence at the edge of a tokamak discharge was recently simulated by Hallatschek and McKee [5]. 1.3 Overview of the contents An overview of the organization and contents of this dissertation is as follows. In the first part, we study the effects of externally forced flow shear on interchange modes. This is done in Chapter 2, where we present a detailed investigation of the effect of a forced velocity shear on interchange modes, in both linear and nonlinear regimes, using numerical simulations with a dissipative 2D MHD fluid code. A description of the equations used is given in Sec. 2.2. A summary of the results of linear theory is given in Sec. 2.3. We then focus on the effect of increasing flow shear on wavenumber spectra of linearly unstable modes (Sec. 2.4.2), the residual turbulence, and particle transport (Sec. 2.4.3). We compare the results of the nonlinear simulation and the residual turbulence with the spectra of magnetic fluctuations observed in the MCX (Sec. 2.4.3.5). We note that, since we are dealing with macroscopic instabilities, the spatial variation of flow shear should be taken into account. The latter can introduce Kelvin Helmholtz (KH) instabilities. Therefore, we give a derivation of a generalized version of the Rayleigh inflexion criterion for KH 6 instabilities in Sec. 2.5, and follow up with simulation results to show that a complete stabilization of interchange modes is possible only if the Rayleigh criterion is simultaneously satisfied. In the second part, we consider aspects of resonantly driving a time varying velocity shear in magnetized plasma to stabilize turbulence. This work is presented in Chapter 3 to Chapter 6. The two main topics are to study the damping of GAMs from phase mixing, and to study the effects of nonlinearity on the GAM resonance. We begin with an introduction to GAMs in Chapter 3, where we describe and derive the nonlinear equations for the GAMs. In Chapter 4, we focus on the linear GAM system, and study the effect of phase mixing on the power dissipation in the undriven and driven problems. For the undriven (homogeneous) problem, we solve the equation approximately by an expansion in the viscous dissipative coefficient, and we find that the modal decay is faster than exponential (in time), with the effective time constant 1/3 ~? ? , where ? is the viscosity coefficient. For the driven problem, we present an asymptotic analysis of the system of equations, and verify the results with simulation. We show that the width of resonance is broader than that predicted without the phase mixing effect. The resonant peak is correspondingly smaller. Then in Chapter 5, we study the nonlinear GAM equations with a sinusoidal driving term. We expand the equations in order of the small driving amplitude for driving frequencies far and close to resonance. We then verify the results, in particular, the frequency spectra, with a nonlinear simulation. We find that the nonlinear terms can suppress the resonance in the sense that the magnification of the driver amplitude decreases when the driving amplitude increases, making it harder to get a greater output amplitude at steady state. 7 In Chapter 6, we evaluate the power requirements caused by collisional damping from magnetic pumping and dissipation by phase mixing in DIII- tokamak, at the edge and core, respectively. We find that the power requirements are favorable in both cases. We also estimate the effect of the nonlinear terms with the experimental parameters, and show that nonlinear effects are small compared with the dissipative terms. We then conclude the second part of dissertation in the same chapter. Finally, we summarize the work in this dissertation in Chapter 7. We include three supplementary examples in the Appendices related to the second part of the dissertation. In Appendix A, we outline a derivation of the solution of Alfven waves propagating in an inhomogeneous medium, along with the corresponding simulation results. This calculation is provided as a simple example which illustrates the enhanced dissipative effect caused by phase mixing. In Appendix B, we discuss how resonant driving of magnetosonic waves could be suppressed by nonlinearities in the system. An asymptotic analysis and a numerical simulation are described. In Appendix C, we present a brief review of the Duffing equation as a reference and as comparison with the nonlinear GAM system. 8 Chapter 2: Residual Turbulence from Velocity Shear Stabilized Interchange instabilities 2.1 Introduction It is well established that velocity shear can stabilize electrostatic magnetized plasma instabilities. [19,20,30] The shear in the flow distorts growing convective eddies of the unstable mode and, if the shearing rate, V?, exceeds the growth rate, the mode can be stabilized. This stabilization has been shown for both the ideal interchange instability and drift-type instabilities. [14-16,30,31] The former mode has a robust growth rate and spans wavelengths from macroscopic to gyro-radius scale. Drift modes generally have robust growth rates at short, gyro-scale wavelengths. Considerable experimental data from magnetized flowing plasmas is consistent with velocity shear stabilization of turbulence.[19,20]. In this chapter, we examine the effect of velocity shear on the macroscopic, broad bandwidth, ideal interchange instability, with particular attention to the effect of increasing velocity shear on the unstable wave number spectrum and to the effect of increasing shear on residual turbulence and transport. We examine the spectra of both linear instabilities under velocity shear, and the wavenumber spectra for fully developed turbulence in an unstable system in the presence of applied velocity shear. We compare theoretical turbulent spectra with experimental results from the 9 Maryland Centrifugal Experiment [1], an experiment that studied the effect of velocity shear on ideal interchange modes. For the macroscopic interchange mode at long wavelengths, the velocity shear cannot be considered as approximately constant across the mode width (as can be done for the short wavelength highly localized drift modes) since the broadly spread mode also samples varying velocity shear. In other words, the curvature in the flow also needs to be taken into account. [32-34] As is well known, flow curvature precipitates the Kelvin-Helmholtz instability. Thus, in this chapter, we also study how the destabilizing influence of the KH instability on the interchange mode, from flow curvature, competes with the stabilizing influence of the flow shear. In the rest of this section, we provide a brief description of previous work and a summary of results from this present chapter. Ideal interchange modes can be stabilized by transverse magnetic fields or magnetic shear [32]. Velocity shear provides added stabilization and, in fact, can stabilize ideal interchanges even in the absence of transverse magnetic stabilization. The stabilization in the latter case is incremental in that there is residual turbulence for a given amount of velocity shear. [18] Here we focus on the character of this residual turbulence. The effect of velocity shear on ideal interchange modes without magnetic stabilization was studied in the ordinary fluid case by Kuo [10]. Flute instability theory for strongly magnetized plasma in slab geometry has been studied analytically and numerically. [18,17] Velocity shear has been shown to stabilize ideal interchange sausage modes and kink modes in a Z pinch, with mode suppression increasing with the Mach number of the flow [21]. Velocity shear stabilization has also been extended to the case where the 10 effective gravitational force arises from centrifugal forces of plasma rotation; in particular, it has been shown numerically that centrifugally confined plasma in a mirror field system is stable to flutes, at supersonic rotation with shear in angular frequency. [22] The latter study was in agreement with the rotating mirror experiment MCX, which showed suppression of expected flute interchanges for rotation Mach numbers in the range 1.5 to 2.5. [1,35,36]. The present chapter investigates the linear and nonlinear spectra of interchange modes, given varying amounts of flow shear. It is important to note that, for an ordinary Rayleigh-Taylor equilibrium configuration with a cross-field shear flow, the linearized operator of the system is not of Sturm-Liouville form and, thus, the corresponding eigenmodes obtained are not guaranteed to form a complete set [15,18]. Therefore, in this chapter, linear growth wavenumber spectra are obtained by studying growing randomly seeded initial values. Nonlinear wavenumber spectra are studied similarly in time-averaged mode. In our study of linear spectra, we recover the well-known result that the linear growth rate in the absence of flow shear increases monotonically with the transverse wavenumber (in the ideal range of instability) [32]. As the ambient flow shear is increased, we find that this brings down the growth rate overall, but more effectively for modes of either short wavelengths or of long wavelengths, thus leaving the growth rates to peak at midrange wavenumbers which correspond to roughly circular convection cells in our simulations. Spectra are similarly studied in the nonlinear regime, by time-averaging over the turbulence in the nonlinear phase of the simulation. We find that the shape of the mode spectrum remains roughly the same 11 as in the linear regime; in particular, the wavelength of the dominant mode in the nonlinear averaged steady state is approximately the same as that of the mode with the maximum linear growth rate, for modest velocity shear. The foregoing spectral analysis is done for high Reynolds number turbulence. We also study the system for lower Reynolds numbers. We find that increasing viscosity has a stronger stabilization effect on short wavelength modes, compared to the longer ones, which has the effect of shifting the peak of linear growth rates to low mode numbers, ( ym ~ 1 or 2), that is, to cells which are elongated in the y-direction. Experimental observation and analysis from MCX shows the presence of strong nonlinear interaction between the 0m? = and 2m? = interchange modes [2,3], where m? is the azimuthal mode number. For MCX, 2m? = corresponds to highly elongated cells in the azimuthal direction. Thus, this observation is consistent with our high viscosity simulations. The MCX observations were made in the outer periphery of the plasma, where one may expect significant damping from neutrals. Thus, the viscous effect of plasma in MCX is important in edge interchange modes. Because the system is not Sturm-Liouville in nature, velocity shear stabilization does not in general lead to sharp transitions between stable and unstable regions in parameter space. In particular, the possibility of residual turbulence has to be addressed for a given shear. Thus, in this chapter, we examine the behavior of the residual nonlinear fluctuations as a function of increasing V? shear. We measure the size of the residual turbulence and the corresponding degree of profile flattening, in quasi-steady state, for fixed V? shear. The numerical data shows that fluctuations are reduced as flow shear increases and the density profile is restored to laminar state as 12 the Mach number of the flow exceeds about 2, in agreement with previous work. In addition, the degree of turbulence at different flow shears is evaluated quantitatively by comparing the total particle flux with particle transport from classical diffusion [32,37]. We find that the transport is dominated by turbulence at low shear but reduces to classical transport levels at sufficiently large shear. In the next section, we first describe the MHD equations and the ordering used, and describe the equilibrium. In Sec. 2.3, we give a summary of linear Rayleigh-Taylor theory, and a derivation of the stability criterion for the KH instability with non- constant density profile. In Sec. 2.4.1, we describe the equilibrium and parameters for our 2D MHD simulations. The linear spectrum is obtained in Sec. 2.4.2, and nonlinear analysis, including fluctuation measurements, profile flattening, evaluation of turbulent flux, effect of viscosity on nonlinear spectrum, and comparison with the magnetic fluctuation measurements in MCX, is given in Sec. 2.4.3. In Sec. 2.5, we discuss the KH instabilities and analyze the effect of Rayleigh criterion (with density profile) on an RT unstable configuration with flow shear, numerically. We conclude in Sec. 2.6. 2.2 Equations: 2.2.1 2D MHD Equations We consider a strongly magnetized plasma in rectangular geometry with magnetic field ?( , )B B x y z=harpoonrightnosp , plasma density ( , )n x y , and transverse flow ( , )u x y?harpoonrightnosp in the x-y 13 plane. The system is 2D, with axial symmetry / 0z? ? = assumed. Isothermal conditions are assumed. A gravitational field ?g gx= ?harpoonrightnosp is allowed. The axial flow, zu , is assumed to be zero. The governing equations are: ( )n nu S t ? ? ? +? ? = ? harpoonrightnosp , (2.1) 2 2 ( ) 8 z o u B n u u nT ng nu F t ? pi ? ? ? ? ? ? ? ??? ? + ?? = ?? + + + ? +? ?? ??? ? ? ? harpoonrightnosp harpoonrightnospharpoonrightnosp harpoonrightnosp harpoonrightnosp harpoonrightnosp , (2.2) ( ) 2 2 4 z z z B c u B B t ? pi? ? ? ? +? ? = ? ? harpoonrightnosp . (2.3) Here ? , ? are the viscosity and conductivity respectively, particle mass 1M = , S is a particle source, and F harpoonrightnosp is an external force term which is used to introduce a flow shear in the system. We note that B? harpoonrightnosp is not excited in this 2D system, if initially zero, as it only appears as a homogeneous term in Faraday?s Law, ( ) ( )2 24t B B u uB c B? pi? ? ? ? ?? = ?? ??? + ?harpoonrightnosp harpoonrightnosp harpoonrightnosp harpoonrightnospharpoonrightnosp harpoonrightnosp . Thus equations (2.1)-(2.3) form a closed set with variables { n , xu , yu , zB }. We also note that, although the introduction of F harpoonrightnosp in (2.2) is phenomenological, there is correspondence with the experiments. For example, in the MCX experiment (which we will discuss in Sec. 2.4.3), the external electric field is perpendicular to the magnetic field, and the corresponding cross-field ?leakage? current together with the magnetic field gives a j B? harpoonrightnospharpoonrightnosp force which maintains the rotation and shear flow of the plasma. 14 2.2.2 Ordering and Reduced Equations To simplify the analysis in the linear regime in the next section, it will be more convenient to have a reduced closed set from (2.1) - (2.3), with suitable orderings. Consider the ordering, 2~xngL nT B<< and 2 2| ( ) |~nu F B? ?? << ?harpoonrightnosp . Then, to lowest order, (2.2) implies 0zB?? = , and then (2.3) becomes: t z z B u B ? ? ? = ? ? harpoonrightnosp , (2.4) which implies 0u? ?? ? = harpoonrightnosp for strong zB . From the latter, we define the stream function ?u z ?? = ?? harpoonrightnosp . (2.5) As a result, (2.1) - (2.3) are reduced to: n u n S t ? ? ? + ?? = ? harpoonrightnosp , (2.6) ? ? du z n z n g dt ? ? ? ? ? ?? ? = ?? ?? ?? ? harpoonrightnosp harpoonrightnosp . (2.7) Here (2.6) and (2.7) form a reduced closed set for { n , ? }. 2.2.3 Equilibrium In the full set, equations (2.1) - (2.3), suppose S depends on x only. This creates equilibrium profiles with only x dependence to balance dissipation. For given ( )S x , with ?g gx= ?harpoonrightnosp , the equilibrium profiles ( )n x , ( )xu x , ( )zB x are determined. Also, an 15 external force ?( )F F x y=harpoonrightnosp is imposed to create a cross field flow shear in equilibrium. Then, the steady state system is: ( ) ' ( )xnu S x= (2.8) 2( 8 ) ' 0o znT B ngpi+ + = (2.9) ( ) '' ( ) 0ynu F x? + = (2.10) ( ) 2 ' '' 4x z z c u B B? pi = harpoonrightnosp . (2.11) Here, the viscous term ( ) ''xnu? is dropped in (2.9), as, from (2.8), ~x xnu SL is proportional to ? , and so the viscous term is proportional to 2? , which can be neglected. We note that S and xu are ( )O ? and ~ ( )F O ? . Thus, xu , S , and F are small. Furthermore, with hard wall boundary conditions, (2.11) can be integrated once without an integration constant. Since ynu is completely determined by ( )F x in (2.10), the system { n , B , xu } is driven to a 1D equilibrium { ( )n x , ( )B x , ( )xu x } by ( )S x . This system is schematically shown in FIG. 2.1. FIG. 2.1. Schematic description of equilibrium. 16 2.3 Summary of linear analysis In this section, we provide a summary of well known results from linear theory. Introduce a perturbation ( ) exp[ ( )]A a x i ky t?= ?tildenosp tildenosp to the system {(2.6), (2.7)}, where A can be n , ? . The resulting eigenvalue equation is ( )2 2( )[( ') ' ] ( ') ' 'y y yku n nk k nu gk n ku? ? ? ? ? ?? ? + = ?tildenosp tildenosp tildenosp tildenosp . (2.12) This normal mode equation includes both Kelvin Helmholtz and Rayleigh Taylor physics. 2.3.1 Rayleigh-Taylor instability Suppose ( ) 0F x = , thus 0yu = . In this case (2.12) becomes 2 2 2( ') ' 'n nk gk n? ? ? ?? =tildenosp tildenosp tildenosp . (2.13) To model the solutions of this eigenvalue equation at long wavelength, i.e., ' '/n n? ? < >? ?? . (2.14) Then the solution is | |k xce? ?=tildenosp , (2.15) for some normalization constant c, with eigenfrequency [32] 2 1 0 1 0( ) ( )gk n n n n? = ? ? + . (2.16) This is unstable for 1 0n n> , with the growth rate | |? proportional to k . 17 For ' '/n n? ? >>tildenosp tildenosp , the short wavelength limit, consider the Taylor-expanded density profile at the location of maximum '/n n . 2 2 ' (1 ) n n n n x L L= ? , (2.17) with n L the density scale length, such that | ' | n L ? ?>> tildenosp tildenosp . Then " ' 'n n? ?>>tildenosp tildenosp under this ordering and (2.13) can be written as 2 2 2 2 2'' 1 0 n n n g L g L x k L ? ? ? ? ? ? ? + ? =? ?? ? tildenosp tildenosp . (2.18) The eigenfunctions of (2.18) can be written in terms of Hermite polynomials: 2 22 ( ),ox xl l oe H x x? ?=tildenosp (2.19) with the corresponding eigenvalues [38] ( ) ( ) 2 2 2 , 0,1,... 1 2 1 ( 1 2) ( ) 1 2 for , 1 (2 1) n l n n n n g L l l l kL l g L kL l l kL ? ? = =? ?+ + + + ? +? ? ? ? >> + + (2.20) where ( )1/22 2 2(2 1) 4 (2 1) 2o nx l k L l k= + + ? + , which is also the scale length of the eigenfunctions. The condition o nx L<< is required for consistency, which has the implication that only eigenfunctions with n l kL<< are allowed in the expression of ox . From (2.20), for 1nkL >> , i.e., eigenmodes elongated in x direction, we have 2 / n g L? ? ? , (2.21) which is purely growing for density gradient in the opposite direction to gravity. The growth rate is almost independent of k . 18 2.3.2 Kelvin-Helmholtz instability, Rayleigh inflexion criterion We now consider 0g = . Multiply both sides of (2.12) with *?tildenosp , the conjugate of ?tildenosp , and integrate over x, assuming hard wall and either free slip or no slip boundary conditions at x = 0 and xL , the length of the box along x. This results in ( ) 2 *2 2 2 20 0 ( ') ' | | ( )| ' | | | | | x xL L y y y k nu ku n k dx dx ku ? ?? ? ? ? + = ? ? ? tildenosptildenosp tildenosp . (2.22) From (2.22), 2 * 20 ( ') ' | | Im( ) 0 | | xL y y nu k dx ku ? ? ? = ? ? tildenosp . (2.23) Therefore, if ( ') ' 0, 0y xnu x L? ? ? , (2.24) the frequency must be real and the system is stable. If a constant density profile is used, (2.24) becomes '' 0yu ? for 0 xx L? ? , which is the well-known Rayleigh inflexion criterion [4,32]. Thus (2.24) is a generalized version of the latter for system with non-constant density profile. 2.4 Rayleigh-Taylor instabilities with flow shear In what follows, the effect of velocity shear on Rayleigh-Taylor unstable modes will be analyzed in detail by nonlinear numerical simulation. In general, sheared flow can stabilize RT modes, but an introduction of flow shear to the system may be KH unstable if (2.24) is not satisfied. Moreover, the condition is more complicated since 19 (2.24) involves the density profile and there is gravity in the system, and thus the KH and RT physics is coupled. 2.4.1 Equilibrium profiles and examination of KH stability criterion The 2D code used in our simulation was developed by Guzdar et al [39]. The equations used are (2.1) - (2.3). The system is a 2D box with 0 xx L? ? , 0 yy L? ? . The length is normalized to xL , and the speed is normalized to the Alfven speed Av based on some initial reference oB and on . A resolution of 31 and 101 grid points in the x and y directions, respectively, is employed. A time step of 0.002 for each iteration is used. The boundary conditions at {0, }xx L= are hard conducting wall for density and zB (symmetric), and no slip for xu and yu (antisymmetric). For {0, }yy L= , the boundary conditions are periodic for all the variables. In all the simulations, the parameters used are 1xL = , 5yL = (i.e., a box elongated in y direction), 1Av = , 0.3oT = , 32 10x? ?= ? , 42 10y? ?= ? , 32 10x? ?= ? , 42 10y? ?= ? . The gravity parameter 0.1g = and points in negative x direction. The particle source is made up of two Gaussians, one positive and one negative peaked near both x boundaries, as follows: 2 2( 0.1) ( 0.9)( ) ( )a x a xxS x b e e? ? ? ? ?= ? + , (2.25) 20 where the parameters a and b are used to adjust the strength and width of the source (4 and 60, respectively, for our simulations). Notice that the spatial average of ( )S x is zero so that no net particles are introduced or removed from the system, and the source strength scales as the resistivity, that is to say the source strength depends on the resistive diffusion, which is weak. An initial magnetic field 1 0.2 tanh[10( 0.5)]zB x= ? ? and density 2 / (8 )o z on n B Tpi= ? are used, after which the system is allowed to come to resistive equilibrium given the sources. The relaxed profiles are shown in FIG. 2.2a-b. Notice that the constant on can be used to adjust the minimum density in equilibrium, and also to vary the 2( / 8 )nT B? pi= value of the system; in FIG. 2.2a 1.5on = is used as an example. FIG. 2.2a FIG. 2.2b 21 FIG. 2.2c FIG. 2.2a, b. The 1D equilibrium profiles. Note that the density gradient is opposite to gravity, thus unstable to RT instability. FIG. 2.2c. Cross field flow shear profiles for different strengths characterized by coefficient A in (2.26). A parabolic flow shear profile is produced by setting ( ) 2 xF x A?= , so that, from (2.10), (1 )ynu Ax x= ? . (2.26) Here A is the parameter which characterizes the strength of the shear (FIG. 2.2c). From the Rayleigh inflexion theorem, a parabolic flow shear with constant density profile, and having '' 0yu ? everywhere along x is ideally KH stable. (There could exist a viscous KH, with a much smaller growth rate, of order 10-3, but this is not relevant here). However, it does not follow that the same profile could satisfy (2.24) for a system with density variation. Moreover, (2.24) is derived for the system without gravity; thus, for 0g ? , satisfying the condition (2.24) may not guarantee stability. Later, it will be seen that for large flow shear, even if 0g ? , (2.24) is still a good rule of thumb for a prediction of the stability of the system in the case that the 22 RT mode has been suppressed completely by velocity shear. For the parabolic profile (2.26), we calculate that 2 2 2 min ( ') ' ( ) '' ( ) ' '/ ''/ ( '/ ) 12 [ (1 2 ) ' (1 )( ' / '')] 12 (1 2 ) ' (1 )( ' / '') y y y y ynu nu nu n n nu n n nu n n A x n x x n n n n A x n x x n n n n = ? ? + ? ? = ? + ? ? + ? ?? ?? ? ? ? ? ? + ? ? + ? ?? ?? ? . (2.27) If the initial average density on is increased, say from 1n to 1n n+ ? , the minimum density at equilibrium minn will also increase while 'n doesn?t change much if n? is only a fraction of n (see FIG. 2.3a). Thus it is possible to make ( ') 'ynu negative everywhere by increasing on , for all shear parameters A used herein (see FIG. 2.3b). FIG. 2.3a FIG. 2.3b FIG. 2.3a. The density profile can be shifted vertically without largely altering 'n by adjusting the initial average on . FIG. 2.3b. For a parabolic shear flow, the generalized Rayleigh criterion is not guaranteed to be satisfied; the criterion can be satisfied by ?upward shifting? the whole density profile. To stay away from KH instability and investigate the effect of flow shear on the RT only, 1.5on = is used in the rest of this main section of simulation and analysis, in 23 both linear and nonlinear regimes. Lower values of on , for when the KH criterion is marginal ( ~ 0.9on in FIG. 2.3b), will be investigated in Sec. 2.5. 2.4.2 RT instability stabilized by velocity shear in the linear regime We first investigate in detail the effect of velocity shear on the linear mode spectrum. For the 1D equilibrium shown in FIG. 2.2a-c, where 1.5on = and the KH stability criterion is well satisfied, and for a fixed shear parameter A, a perturbation in density 610 cos(2 / 5)yn m ypi?=tildenosp , where ym is an integer, is added into the system at 0t = . After some transients, the system is found to grow exponentially, until the amplitude of the perturbation is increased by ~104 times, at which point it is seen to go nonlinear. The growth rate RT? is measured in the linear regime, and plotted for a range of ym , the harmonics along y, where 2 / 5y yk mpi= . ( )RT ym? is plotted for various shear parameters A and shown in FIG. 2.4. 24 FIG. 2.4. Rayleigh Taylor unstable growth rate in linear regime, as a function of yk , for parabolic flow shear of various strength parameter A. From FIG. 2.4, we first notice that when there is no shear ( 0A = ), the growth rate increases with ym but flattens out eventually. This is consistent with (2.16) and (2.20) in that the growth rate goes like k for small k , but becomes independent of k when 1xkL >> . Secondly, we observed that an increase in flow shear lowers the growth rate for any one particular harmonic. For weak shear ( 0 1A? ? ), higher harmonics are stabilized more effectively than their lower counterparts, i.e., modes of shorter wavelength along y are stabilized before longer ones as A increases. As an explanation, we note that viscosity can suppress the growth rate (with or without flow shear) and the viscosity damping rate ~ k2?, which is proportional to 2ym . Thus the growth rate ( )ym? for A = 0 is expected to drop and the curve rolls off when my is high enough such that 2 ~ 'yk g n n? . This corresponds to , ~ 28y criticalm in our case 25 (not shown in FIG. 2.4). Modes with smaller growth rates are stabilized preferentially by flow shear first; thus, modes with 1nkL << or 2 ~ 'k g n n? [40] will be preferentially stabilized compared with the modes between these extremes. We note that complete stabilization is attained for all k for 3A ? . From the spectrum, for 1A ? the growth rates are seen to peak at about 6ym = or 7 . This corresponds to ~ 15y nk L . Also notice that the system is only linearly unstable in the region where 'n is opposite to gharpoonrightnosp , i.e., from 0.2x = to 0.8 in FIG. 2.2a. The approximate length scale of this unstable region is ~ 1.6effL . 2.4.3 RT instability with velocity shear in the nonlinear regime We now explore the nonlinear regime. For a fixed shear parameter A, a random initial perturbation in density is added and the system is allowed to evolve with time until the nonlinear terms become important. The system reaches a turbulent steady state, as the turbulent transport, mitigated by viscosity and resistivity, enhances the particle flux needed to attain the average steady state, given the sources ( )S x . In general, a flow shear will suppress the RT turbulence. For example, consider a density profile along the x direction against gravity in FIG. 2.5a, which is perturbed in the presence of different shear profiles and the corresponding steady states are shown. (No shear in FIG. 2.5b, medium shear in FIG. 2.5c, strong shear in FIG. 2.5d.) While the profile is ruptured completely by turbulence when no shear is applied, it is 26 maintained partially (with residual fluctuations) in case of medium shear, and complete stabilization is obtained at strong shear. FIG 2.6a-d are other illustrations of the stabilization effect, showing the time evolution of the amplitude of every mode (in the y-direction) as a function of time, with A as parameter, i.e., time evolution of the Fourier spectrum in the y-direction (at mid-x line). In nonlinear steady state, the dominant mode is 6ym = , which is consistent with the linear spectrum shown in the previous section. Moreover, the amplitude decreases as A increases. For 4A = (or above), the amplitudes decay (from the initial perturbation), RT modes are completely stabilized. In a later sub-section, we will also discuss how the nonlinear spectrum is affected by viscosity. FIG. 2.5a FIG. 2.5b 27 FIG. 2.5c FIG. 2.5d FIG. 2.5a. Laminar density contour plots when the system is unperturbed. FIG. 2.5b. Turbulent steady state of the system without shear; density profile is inverted to a new quasi-stable equilibrium. FIG. 2.5c, a medium strength shear partially restores the initial profile, with residual fluctuations. FIG. 2.5d, RT convection is completely suppressed by strong flow shear with laminar profile recovered. FIG. 2.6a FIG. 2.6b 28 FIG. 2.6c FIG. 2.6d FIG. 2.6a-d Time evolution of the spectra along the y-direction with shear strength A as parameter in the nonlinear regime. The random perturbation in density is added at 0t = . Note that the amplitudes at steady state decrease as A is increased from 2. For 4A = (FIG. 2.6d), the system is completely stabilized and the perturbation decreases with time. Also, the nonlinear spectrum peaks at 6ym = , which is the same as that of the linear one. The effect of flow shear on the instabilities, turbulence, and saturated states, can then be measured in terms of the fluctuation in density at saturation, and in the degree of flattening of profiles. Also comparisons can be made between turbulent steady state profiles and the laminar profiles, and between the turbulent particle flux and the corresponding (laminar) classical particle flux (when the system is one-dimensional), at various levels of shear A. We describe these below. 2.4.3.1 Fluctuation and degree of profile flattening When the system reaches the nonlinear regime and becomes steady on average (i.e., saturated), the time evolution of density at the centre of the simulation box is measured and its temporal fluctuation about the average characterized by a standard deviation. In addition, the slope of the density profile is averaged over all simulation space and time, which characterizes the degree of flattening of the profile caused by the RT disruption. The effect of flow shear on the density fluctuation and the level of flattening (quantitatively expressed as standard deviation over time and space, and, 29 time and space average of density gradient, respectively) are shown as functions of shear strength A in FIG. 2.7 and FIG. 2.8. FIG. 2.7. Density fluctuation (at 0.5x = ) for different shear strengths. FIG. 2.8. Density gradient averaged over all space for different shear strengths. From FIG. 2.7 and FIG. 2.8, when there is weak or no shear flow, the average density gradient is negative and the fluctuation is small. This is because the density profile has collapsed due to the strong RT instability and has reached a marginally stable equilibrium, with average slope the same sign as gravity; the small fluctuation is due to the source term ( )S x and the relatively weak residual convection. At strong shear (A ~ 4 to 5), the average density slope reaches a maximum and levels off as the 30 system is completely stabilized. The fluctuation is nearly zero because of the complete stabilization, with the system reverting to laminar conditions consistent with ( )S x and 1-D classical transport (as depicted in FIG. 2.2a). For shear in between, A ~ 1 to 3, the density profile is partially restored, and the fluctuation is large and attains its maximum, since it is in between two stable configurations (i.e., a new average stable equilibrium after RT disruption near ~ 0A and complete flow shear stabilization for 4A > ). 2.4.3.2 Comparison with laminar flow The time averaged density profiles in x in steady state for different flow shears are plotted and compared with the laminar profile and shown in FIG. 2.9. Here the profiles are averaged over 40 Alfven times. This clearly illustrates the restoration of the configuration when A is increased from 0 (complete inversion of the density gradient) to 5 (reverting to laminar density). 31 FIG. 2.9. Comparison of steady state density profiles with laminar density at various shear strengths A. 2.4.3.3 Comparison between turbulence flux and classical particle flux The degree of turbulent transport in steady state can be quantified by measuring the total particle flux across the line / 2xx L= (i.e., the 1D ?plane? in our 2D box), and compared with the corresponding theoretical classical flux of laminar flow as if there were no turbulence. In turbulent steady state, the system is almost static, as if it were laminar. Thus, with the ordering /t x xu L? << , from (2.2) and (2.3), to lowest order, we have 2( 8 ) 0o znT B ngpi?? + ? ?harpoonrightnosp , (2.28) ( ) 2 2 4z z c u B B? pi? ? ? ? ? ? ?harpoonrightnosp , (2.29) 32 which resemble (2.9), (2.10) and (2.11), except (2.28) and (2.29) are describing a quasi-steady state. Integrating (2.29) once, and since in our system u?harpoonrightnosp and zB?? are zero at boundary, we find 2 4z z c u B B? pi? ? ? ?harpoonrightnosp . (2.30) Eliminate zB?? in (2.28) with the expression (2.30), 2 , 2 [ ( ) ]classical o z c u u nT ng B ? ? ? ?= ? ? ? + harpoonrightnosp harpoonrightnosp harpoonrightnosp , (2.31) which is the classical diffusive transport. Here a subscript is used to distinguish the classical flow from the actual measured flow. Thus, the particle flux resulting from classical diffusion across / 2( 1/ 2)xx L= = in steady state is given by ( ) 2 2, , 1/20 01/2 ( ' ) /y y L L classic x x classic o z xx nu dy n c n T ng B dy? == ? ?? = ? ? +? ?? ? . (2.32) On the other hand, the total flux in steady state is evaluated as ( ) , 0 1/2 yL total x x x nu dy = ? = ? . (2.33) Both the classical flux and the total flux are measured in the simulation, using (2.32) and (2.33) respectively in the nonlinear steady state. The time averaged results are shown in FIG. 2.10. The source (2.25) is independent of time, introducing and removing particles at 1/ 2x > and 1/ 2x < respectively; thus the total flux is negative and nearly constant. At zero or weak shear (A < 1), the classical diffusion only accounts for ~ 1/3 of the total flux, while the rest must be due to the turbulence. As velocity shear increases, the classical flux increases and approaches the total flux, showing a relative reduction of turbulent transport. For large shear (A > 4), classical 33 flux is nearly the same as the total flux, thus completely accounting for particle transport as the turbulence is wholly suppressed. FIG. 2.10. Effect of flow shear on particle flux across the 1/ 2x = axis. The transport is towards negative x direction. At strong shear, classical diffusion accounts for nearly all of the particle transport. 2.4.3.4 Effects of cross field viscosity on stability It is also worth noting that at a fixed flow shear, an increase of viscosity in the y- direction imposes a stabilizing effect in both linear and nonlinear regimes. The linear growth rates for different y? at A = 2, are shown in FIG. 2.11. The growth rates decrease as viscosity in the y-direction increases, and the suppression is stronger for shorter wavelength modes. While the long wavelength modes ( 1, 2m = ) which remain have smaller growth rates, the overall result is stabilizing. Correspondingly, in the nonlinear regime, the density fluctuation is smaller as viscosity is increased, and long wavelength convection dominates. 34 FIG. 2.11. Linear growth rates as a function of ym , with y? , the viscosity in the y-direction, as a parameter. Shear strength is constant at 2A = . Note the suppression of short wavelength modes when y? is increased. The level of turbulence, with increasing y? , is shown as contour plots of the density fluctuations. This, along with vector plots of the convection, are shown in FIG. 2.12 and FIG. 2.13 respectively, which compare low y? and high y? situations. For small y? ( 42 10?? ), the convection cells appear to be roughly circular. To see this, we note that from FIG. 2.2a that the density profile has a positive slope from ~ 0.1x to ~ 0.8x ; thus the scale size of the convection cell in the x-direction ~ 0.7, and 6 cells are observed in the nonlinear quasi-steady state, as seen in FIG. 2.12a and FIG. 2.13a. This results in a ratio of the scale in the x direction to that of the y direction ~ 4.2 : 5. A similar nonlinear steady state for zero y? (which is not shown) shows 7 convection cells along the y direction, which gives the corresponding ratio of scale sizes ~4.9 :5. Thus, nonlinearly the preferred turbulent cell shape seems to be 35 circular on the nL scale. Notably, the convection is elongated in y at the lowest allowable mode as y? increases. FIG. 2.12a FIG. 2.12b FIG. 2.12c FIG. 2.12. Nonlinear residual density fluctuations at A = 2 for a) 42 10y? ?= ? , b) 34 10y? ?= ? , c) 36 10y? ? = ? . Notice the increase in wavelengths and the decrease in fluctuation levels, as the viscosity is increased. FIG. 2.13a 36 FIG. 2.13b FIG. 2.13c FIG. 2.13. Nonlinear residual convection for A = 2 at a) 42 10y? ?= ? , b) 34 10y? ?= ? , c) 36 10y? ? = ? . Notice the decrease in mode numbers (i.e., elongated convection cells) as the viscosity is increased. 2.4.3.5 The MCX experiment To conclude this section, it is interesting to compare the above results with magnetic fluctuation data from the Maryland Centrifugal eXperiment (MCX). The MCX basically makes use of centrifugal force of rotating plasma to provide axial confinement. An accompanying flow shear provides stability to the interchange. The setup resembles a mirror machine in that the magnetic configuration is a long solenoid with axisymmetric mirror end fields. There is an axial conducting core which is insulated from the outer wall. A radial electric field is applied by biasing this inner core at high voltage with respect to the outer wall. This results in an E B? drift causing the plasma to rotate azimuthally. Since the magnetic field is curved, there is a nonzero component of the radial centrifugal force which is parallel to the magnetic 37 field and points inward towards the equatorial plane. This balances the parallel pressure gradient so that the plasma is confined to the straight solenoid portion of the magnetic field. For large Mach number rotation, the conventional mirror loss cone is almost completely closed because of the favorable centrifugal potential. At first glance, the MCX setup is unstable to flute-like interchange modes from magnetic curvature and centrifugal force, where both effects are in an interchange- unfavorable direction relative to the density gradient, i.e., 0effg n?? < harpoonrightnosp . As mentioned, velocity shear is expected to stabilize the flutes. Basic evidence of such suppression comes from the fact that the plasma is observed to be relatively quiescent for thousands of instability growth times [35,36]. More evidence is obtained from magnetic fluctuation data as analyzed, for example, by Choi et al [2] and Uzun- Kaymak et al [3], as described below. There are 16 magnetic probes arranged azimuthally with uniform spacing of pi/8 radians on the inner surface of the vacuum vessel, and halfway between the mirror throat and midplane. The probes are oriented to measure fluctuations in the axial magnetic field. The set of 16 probes gives the fluctuation zBtildenosp as a function of azimuthal angle ?, or mode number m in spectral analysis, with a resolution up to 8m = . Before introducing the results of magnetic fluctuation measurements in MCX, we give a brief description of the operation ranges and plasma parameters of the setup first. The experiment has T ~ 30eV, 6 1 , ~ 7.5 10th iv cms?? , ~ 10nL cm , and radial range of the plasma 5 25cm r cm< < (i.e. the size of plasma at radial direction is ~ 20a cm ), which give 2 5 1~ ~ 6.12 10g s nc rL s? ?? . The B-field at midplane is 2kG and 38 14 3 ~ 10in cm ? , so that 2~ 3.02 10? ?? . The corresponding viscous damping rate 2 1 2 2~ ~ ~ 184n n c s L L ? ? ? pi ? . We then obtain the Reynolds number at the core of the plasma 42Re ~ 3.4 10 1 g nL ? ? = ? >> . The plasma in MCX is observed to rotate at a sonic Mach number between 2 and 3. The magnetic fluctuation data suggests that there are 2 ?blobs? that appear to rotate at the local E B? speed in the sheared flow. This corresponds to 2m? = , or in terms of wavenumber, ( ) ~ 1.33n Rk L? . Spectral analysis indicates that modes with shorter wavelengths ( 2m? > ) are of small amplitude. Theoretically, linear, shear free interchange modes with short wavelengths are more unstable than long wavelength ones [refer to Sec. 2.3.1, (2.16), (2.20)]. Therefore, shear free interchange theory cannot explain the spectrum of magnetic fluctuations in MCX. However, the presence of shear flow in MCX would stabilize the shorter wavelengths, as we have shown. Consider our simulation for the RT instability in the linear regime in Sec. 2.4.2. A suppression of unstable modes with short wavelengths is observed in our simulation, as velocity shear increases (FIG. 2.4): while the linear growth rate[ ( )RT ym? ] is relatively flat for 6ym > for shear free conditions, the plots curve downwards as the parameter A increases. Moreover, the unstable region shrinks. As a result, unstable modes are at long wavelengths as flow shear is increased. This has the correct trend compared with the magnetic fluctuation spectrum in MCX. 39 However, in our simulation in the presence of flow shear, the wavelengths observed in residual nonlinear fluctuations are still significantly too short to explain the experimental results in MCX in the region ~ 1Rk L? . In the simulation, the residual fluctuation, shown in FIG. 2.5c, is dominantly of 6ym = , corresponding to ~ 12y effk L . A spectral analysis has been performed of the steady state of the simulation and confirms that the same wavenumber is observed, as shown in FIG. 2.14; here 2A = is used, i.e., a parabolic shear profile with maximum slope 2( / )A xv L is chosen, or a slope of 3.7( / )s xc L in our system (since 0.3? = ). FIG. 2.14. Spectral analysis of the residual fluctuations of the y-flow in the y direction at A =2 and low viscosity, which is corresponding to FIG. 2.13a. A possible explanation for the discrepancy is the effect of viscosity. In the MCX magnetic fluctuation measurements, the magnetic probes are located in the edge region of the plasma where the plasma is more viscous than in the inner region, because of a relatively higher fraction of neutrals. This increase in viscosity corresponds to a larger viscous term than the one used in our simulation for the results shown in FIG. 2.4, FIG. 2.5c, and FIG. 2.14. A larger viscosity would further suppress short wavelength RT modes. Such an effect is illustrated and shown in our simulation in the previous subsection (Sec. 2.4.3.4); there we found that an increase 40 of viscosity in the y direction further stabilizes the short wavelength modes in both linear (FIG. 2.11) and nonlinear (FIG. 2.12, FIG. 2.13) regimes. FIG. 2.15 shows the spectral analysis of the steady state at high viscosity (30 times larger than that in FIG. 2.14) for 2A = , in which the dominant mode is 1ym = , or ~ 2y effk L . The wavelength is greater than that in the weakly viscous case, with better agreement with the results from the MCX (for which the effective ( ) ~ 1.33n Rk L? ). FIG. 2.15. Spectral analysis of the residual fluctuations of the y-flow in the y direction at A =2 and high viscosity, corresponding to FIG. 2.13c. The viscous effect at the edge is likely related to charge exchange friction between ions and neutrals. Based on the higher value of dissipation used in the code (FIG. 2.15), the equivalent damping rate due to charge exchange [38], ,CX CX th iN v? ?= , (2.34) is given by CX? (normalized to /Av a ), is about ( )2 36 10x nL L ?? . Here, N is the density of the neutrals, 15 2~ 4 10CX cm? ?? is the charge exchange cross section, and ~ 2n xL L is the normalized density scale length along x in the simulation during turbulence. Together with the operation ranges of MCX given in our description previously, we can estimate / iN n , the fractional neutral density at the edge needed to 41 give this viscosity. We find this ratio to be 3/ ~ 1.1 10iN n ?? . This is reasonable given the edge conditions on MCX. Also, the corresponding Reynolds number at the edge is Re ~ 12 Reedge core<< . 2.5 Kelvin-Helmholtz instabilities In the following, the minimum density in equilibrium is adjusted (by changing on ) so that (2.24) is either barely satisfied or violated (FIG. 2.3a,b). Thus, we are allowing KH instabilities. A linear analysis of stability is performed by perturbing the density with 6~ 10n ?tildenosp for the harmonic 5ym = . Viscosities 210x? ? = and 42 10y? ? = ? are used. The corresponding growth rates are measured as a function of on and shear strength A, and shown in contours in FIG. 2.16. FIG. 2.16. Contour plot of linear growth rate as a function of on and A, for 5ym = . The contours are labeled by the growth rates. 42 As seen, a flow shear does not necessarily stabilize the system. Complete stabilization with increasing shear is only possible in the region where (2.24) is satisfied, i.e., above the ?KH criterion boundary? in FIG. 2.16 (this includes the regime where the RT effect is analyzed in Sec 2.4). In the region below the boundary, a strong shear in fact destabilizes the system. Consider moving along a horizontal line for 0.8on = where the criterion is violated; the growth rate decreases initially and reaches a minimum at A ~ 1, as the RT effect is partially stabilized by flow shear. However, it rises again for any further increase in shear flow, as the KH effect becomes significant and the system is destabilized. Another remarkable feature is the dramatic drop in growth rate at the criterion boundary for fixed A as on increases. First consider the vertical line A = 0, i.e., no shear. The growth rate decreases slowly when on increases, which is explained by the RT growth rate prediction (2.20) where '/n n is decreasing. Then compare with the vertical line A = 5: a sharp drop in growth rate is observed as the criterion boundary is crossed (indicated by the contour density). The origin of the effect is the KH, since the RT instability is completely stabilized. To illustrate this, the corresponding growth rates in a system without gravity at A = 5 are measured for comparison and plotted in FIG. 2.17, which shows a similar stabilization behavior in the vicinity of the KH criterion boundary. The system is highly unstable when (2.24) is not satisfied as the shear is strong, but is stabilized once the zeros of ( ') ' 0ynu = are removed. (The two curves are not coincident as the equilibrium profiles with and without gravity are slightly different). 43 FIG. 2.17. Linear growth rates for strong flow shear in the vicinity of the KH criterion boundary with and without gravity. 2.6 Summary of Chapter 2 In summary, we have conducted an in depth study of the effect of flow shear on macroscopic ideal interchange modes in magnetized plasmas, using a 2D dissipative MHD code, in both linear and nonlinear regimes. By simulating the system as an externally forced initial value problem, we have shown that in the linear regimes, a cross-field velocity shear V? suppresses growth rates over a broad spectrum of wavenumbers. The suppression is more effective for modes with either very short or very long wavelengths, leaving the unstable region in the k? spectrum to condense and peak at 6ym = , as V? increases. This is in contrast with the shear-free spectrum, which gives k? ?? at long wavelengths and asymptotes to a constant in the short wavelength limit (ignoring viscous damping). We obtain a similar spectrum by simulating the nonlinear turbulent steady state, and the peaking wavenumber corresponds to roughly circular convection cells of size comparable to the density 44 profile scale length. A spatial Fourier transform of residual fluctuations at different shear strengths is taken: again, we show that the spectrum peaks at 6ym = . Further, we measure the average density fluctuation level in steady turbulence, and the flattening of the density profile in turbulent steady state, and find that both decrease as V? increases. Complete stabilization is possible when V? is greater than a few times the classic interchange growth rate, wherein we observe that the original (RT unstable) density profile is recovered and the density and flow profiles revert to laminar. In addition, we measure the turbulent particle flux in turbulent steady state and compare this with the classical diffusive transport. We show that the latter eventually dominates and accounts for all the transport as V? increases, again showing complete stabilization in the strong shear limit. We compare our simulation results with magnetic fluctuation measurements from MCX, which studied interchange modes in the presence of flow shear. By comparing with the experimental spectral data, we show that the experimental results agree with our simulation results if done at higher viscosity. We confirm that this level of viscosity is consistent with the Reynolds numbers at the edge of the plasma in MCX. Both the simulation and the data are consistent with elongated convection cells in the nonlinear regime. Finally, we investigate the effects of Kelvin-Helmholtz on the system, which must be present in any flow shear system, and show from our simulation that we can stabilize an RT unstable configuration with flow shear only when the generalized Rayleigh Inflexion criterion is satisfied. This restriction still allows a broad parameter space in which the RT mode is effectively suppressed. 45 Chapter 3: Geodesic acoustic modes 3.1 Introduction The geodesic acoustic mode (GAM) is well-known as an MHD normal mode in tokamaks. The GAM is an axisymmetric poloidal oscillation of a toroidally confined magnetized plasma (see FIG. 3.1). When a flux tube is displaced poloidally to regions of different magnetic fields, the flux tube gets compressed (decompressed) because of flux freezing. Then, because of pressure changes, the centrifugal force from thermal motions and from the curvature of the toroidal field acts as a restoring force on the displaced flux tube. The result is a poloidal oscillation within a flux surface. The curvature of the magnetic field acts as an effective radially outward gravity on the flux tubes, so the system resembles a nonlinear pendulum. The frequency of GAM oscillations is found to be 2G sc R? = , where sc T M= is the ion sound speed of the plasma, R is the major radius of the tokamak, and T and M are plasma temperature and ion mass respectively. GAMs have been identified experimentally in tokamaks [6-8] as peaks in the broad frequency spectra of plasma turbulence. Theoretical work on GAM oscillations has been well developed, in both collisional fluid [23,24] and collisionless kinetic theory [25-28,41-43]. In this thesis, we will describe the physics of GAMs from the fluid approach. 46 FIG. 3.1. Schematic illustration of geodesic acoustic mode in tokamak. As we have discussed earlier, flow shear can stabilize MHD instabilities. The related theory has been well established, and there is considerable experimental evidence. In general, the stabilization effect is significant when 'V ?> , where 'V is the gradient of the shear flow, i.e., shear strength, and ? is the growth rate of the instability. Insofar as the GAM is an MHD normal mode of toroidal plasma characterized by poloidal oscillations of flux tubes, one may consider the possibility of establishing a time-dependent poloidal flow shear in tokamak by driving GAMs at resonance. This idea was raised by Hallatschek and McKee [5] recently and an experiment has been proposed. The plan is to drive GAMs to resonance from external currents, which would create large amplitude oscillations of poloidal flow localized within the width of the resonating flux surfaces. The goal would be to establish a time varying flow shear with a maximum shear of ' ~ mV u x? , where mu and x? are the maximum amplitude of the poloidal flow of the resonating flux surfaces and the spatial spread of locally driven GAMs in the radial direction, respectively. To have a significant stabilizing effect, the strength of flow shear should be comparable to the growth rate of the instability. For example, for drift instabilities, 47 the criterion is roughly given by ' ~ ~m drift i s nV u x k c L?? ?? > , where i? is the ion Larmor radius and nL is the density scale length in radial direction. This gives an estimation of the strength of flow shear. The power requirements for driving GAMs to the desired flow shear level would depend on the damping rate of the oscillation. The theoretical and computational work in this part of the thesis is motivated by this proposed experiment. In particular, we examine (i) the effect of enhanced damping from phase mixing of GAMs, because of the differential oscillation frequencies in the radial direction; and (ii) the effect of resonant detuning by nonlinearities. For (i), the effect of phase mixing on waves in an inhomogeneous medium has been studied for Alfven waves [9,44]. We will extend this idea to GAMs. We will show that the radial temperature gradient in a tokamak also leads to a similar phase mixing for driven GAM systems. For (ii), we study the effect of nonlinearity from the convective terms in the GAM equations, and we show that nonlinear detuning can occur in GAMs. This behavior is similar to that in the driven Duffing oscillator [45,46]. In Sec. 3.2, we first describe the magnetohydrodynamic (MHD) equations and the ordering used to derive the reduced nonlinear equations for GAMs. In Chapter 4, we consider the linearized GAM equations with viscous damping, and show that phase mixing occurs in both undriven and driven problems. For the undriven problem (Sec. 4.2), we solve the initial value problem in a multiple time scale expansion, with the viscosity ? as the small parameter. We get an oscillatory solution with a slowly decaying amplitude for the poloidal flow. The amplitude decays as 30exp ( / )t t? ??? ? , where 1/3~ot ? ? . For the driven problem (Sec. 4 3), we introduce a sinusoidal 48 momentum source with constant amplitude in the system, and look for a steady state solution. We first analyze the system asymptotically and then verify our analysis by a numerical simulation. We find that in the presence of phase mixing, the frequency spectrum for GAMs is relatively less peaked and the spatial spread of the mode is greater than that predicted by the typical Lorentzian distribution from a linear harmonic oscillator model. In Chapter 5, we keep nonlinear terms in the GAM equations and introduce a sinusoidal momentum source. We then solve for the non-dissipative steady state solution by an expansion in small driving amplitudes, for driving frequencies close to and far from resonant frequency respectively. We verify the results, in particular, the frequency spectra, by a nonlinear simulation. We find that the nonlinear terms can suppress the resonance in the sense that the magnification of the driver amplitude decreases when the driving amplitude increases, making it harder to get a greater output amplitude in steady state. In Chapter 6, we use plasma parameters from the DIII-D tokamak to estimate the power dissipation by classical collisions and by phase mixed damping, separately for the edge and the core of the tokamak. We find that the conditions suggest a cost effective experiment. We then compare our results to a similar estimation done by Hallatschek and McKee [5], who assume that the power is dissipated by plasma turbulence. We also evaluate the nonlinear Duffing-type effects, given the experimental parameters, and find that this effect is small compared with the damping terms. We conclude in Sec. 6.4. 49 3.2 Equations We consider the ideal magnetohydrodynamic equations: ( )|| 0n nu B nu Bt ? ? +? ? + ?? = ? harpoonrightnosp harpoonrightnosp harpoonrightnospharpoonrightnosp , (3.1) u B nM u u T n j t c ?? ? + ?? = ? ? + ?? ??? ? harpoonrightnospharpoonrightnosp harpoonrightnosp harpoonrightnosp harpoonrightnospharpoonrightnosp harpoonrightnosp , (3.2) 0j?? =harpoonrightnosp harpoonrightnosp , (3.3) u B c ?? = ? harpoonrightnospharpoonrightnosp harpoonrightnosp . (3.4) Assume an axisymmetric magnetic field, i.e., ( )B I? ? ? ?= ? ?? + ?harpoonrightnosp harpoonrightnosp harpoonrightnosp harpoonrightnosp , (3.5) where ? is the poloidal flux function, ? is the toroidal angle, and the flux function ( )I RB?? = . In the following, we will assume the magnetic field to be strong, and flux surfaces to be concentric and circular, with radial coordinate r and poloidal angle ? defined with respect to these circular surfaces. We will assume that parallel thermal conduction is large and thus T is approximately constant on each flux surface. The inverse aspect ratio is defined as r R? = , where cosoR R r ?= + , and oR is the major radius of the geometric center of the flux surface. The parallel component of (3.2) gives ( ) 2: ln 0t su B Bu u c B n? ? + ? + ?? =harpoonrightnosp harpoonrightnosp harpoonrightnosp harpoonrightnosp harpoonrightnospharpoonrightnosp harpoonrightnosp harpoonrightnosp , (3.6) where 2sc T M= . A surface integral of the poloidal component of (3.2) over a flux surface gives ( ) ( )2| | 0pds R B nM du dt T n?? ? + ? =? harpoonrightnosp harpoonrightnosp harpoonrightnospharpoonrightnosp , (3.7) 50 where pB ? ?= ? ?? harpoonrightnosp harpoonrightnosp harpoonrightnosp . By expanding the closed set {(3.1), (3.4), (3.6), (3.7)} in orders of the inverse aspect ratio, a set of reduced equations can be obtained. This has been done in Hassam and Drake (1993) [24], with the ordering ( )||~ ~ ~s s s tu c u c r c? ?? . Here, we reproduce their system of equations and our calculations thereafter will be based on this reduced set of equations. Let 0B harpoonrightnosp be the cylindrical magnetic field and assume the safety factor ( )0( )q r rB RB? ?= is of order 1. Within the assumptions of isothermality and low frequency motions (compared with the fast magnetosonic mode), the magnetic field is essentially static and assumed fixed, and the motion comprises E B? harpoonrightnosp harpoonrightnosp flows within each flux surface ( Eu ), flows parallel to the magnetic field ( ||u ), and the associated density fluctuations ( 1n ) from compression of flux tubes from E B? harpoonrightnosp harpoonrightnosp motion and from sound waves along the field. The resulting equations for { 1n , ||u , Eu } follow from (3.1), (3.6), (3.7) are given by ( ) || ||2 sin 0t E Eu r N u r u? ? ?? + ? ? +? =? ?? ? , (3.8) ( ) 2|| || 0t E su r u c N?? + ? + ? =? ?? ? , (3.9) 22 sin 2 s t E c d u N r ? ? ? pi ? = ? ?integralloop , (3.10) where 0 0 ( )( )E d rc u r B dr ? = , 1 || ( )qR ??? = ? , 1 0N n n= and 0 0 ( )n n r= . {(3.8), (3.9), (3.10)} is the closed set obtained in [24]. For our purposes, we further assume || || /Eu u R? << , in which case (3.8) becomes 51 ( ) 2 sin 0t E Eu r N u r? ? ?? + ? ? =? ?? ? , (3.11) and (3.9) is not needed. Next, we multiply (3.11) by sin? and cos? respectively and integrate over the poloidal angle ?, and define new variables sin 2s dN N? ? pi = ?integralloop and cos 2c dN N? ? pi = ?integralloop . From (3.11), we get 0E Et s c u uN N r R ? ? ? = , (3.12) 0Et c s uN N r ? + = , (3.13) 22 s t E s c u N R ? = ? , (3.14) where (3.10) has been rewritten to get (3.14). And for the ansatz to be true, the assumption || || /Eu u R? << is substituted into (3.8) and (3.9), and we find the self consistency condition ~ st G c qR ?? >> . We will see later that 2G sc R? = , so the condition is 1q >> . Our analysis and simulations in the following sections will be based on the closed set of nonlinear equations {(3.12), (3.13), (3.14)}, with variables { Eu , sN , cN }. To see that {(3.12), (3.13), (3.14)} contains the physics of GAMs, consider small perturbation in the variables about a static equilibrium ( )on r . Then, (3.12) becomes 0Et s uN R ? ? = , (3.15) 52 and, (3.14) and (3.15) together give a harmonic oscillation of poloidal flow Eu (and sN ) with frequency 2o sc R? = , the GAM frequency, as a normal mode of the system. In the next chapter we will discuss the effect of phase mixing on GAMs in the linear regime, using the linearized version of the above equations. Then, we will discuss the effects of the nonlinear terms in the section following. 53 Chapter 4: Phase mixing of GAM oscillations 4.1 Introduction In this chapter, we allow a temperature gradient in the ?? harpoonrightnosp direction, i.e., ( )T T ?= , a general feature in tokamaks. (As mentioned in Sec. 3.2, we will assume the system to be isothermal in the sense that the perturbation in temperature on a given flux surface is small compared with its counterparts in density and flow, i.e., 1 0 ~ E sT T n n u c<> ? , i.e., with an expansion parameter 2 ( 0) 1c oL x? ? ?? ?= = <> ? , and find the consistency condition ( )1 'o ot ? ? ?<< . [ 0sN can be obtained by inserting (4.7) into (4.2).] This is self-consistent for large Reynolds numbers. From (4.7), we see that the oscillations decay as ( )3exp ot t? ??? ? wherever 'o? is nonzero, i.e., where there is a profile in sound speed or temperature along r (i.e., x). The time constant 1/3 1/3 2/36 'o ot ? ? ? ? = is defined by (4.7). This decay is significantly faster than a typical exponential decay, since ot , is proportional to 1/3? ? instead of 1? ? . 4.3 Steady state solution of linearized GAM equations with driving term In this section, we study the driven GAM system. We investigate the steady state solution and the corresponding power spectrum. Consider a sinusoidal (in time) momentum source added to (4.2) as an inhomogeneous driving term: 57 ( )( ) 1 sin( ) s t o s r r o s o c r u N r u d t c r r ? ? ?? + = ? ? + , (4.9) where ? is the driving frequency, and do and ? are constants. Correspondingly, (4.4) becomes 2 2 costt o x t ou u u d t? ? ? ?+ = ? + . (4.10) From an examination of (4.10), we expect there will exist a steady state solution such that the amplitude peaks around the flux surfaces with 0 ( )x? close to the driving frequency ? , i.e., at resonance, and that the amplitude will drop as we move to off- resonant surfaces. Analogous to a driven harmonic oscillator, we look for a steady state solution of the form ( ; ) cos ( ; )sinpu A x t B x t? ? ? ?= + , (4.11) so that the amplitude at steady state is 2 2( ; )a x A B? ? = + . Substituting (4.11) into (4.10), we get ( )2 2o xx oA B d? ? ?? ?? ? = , (4.12) ( )2 2 0o xxB A? ? ??? + = . (4.13) First, we examine the off-resonant regions, where 2 2 2o x? ? ??? >> ? in (4.12) and (4.13). In these regions, we find 2 2( ; ) , 0( ) o o dA x B x ? ? ? ? ? ? ? , (4.14) to lowest order. This is consistent with the fact that the oscillation is in phase with the driver for driving frequency well below the GAM frequency, and completely out of 58 phase when driven well above it. To estimate the lowest value of o? ?? such that the above ordering and (4.14) are still valid, we substitute (4.14) into (4.13) to evaluate B to the next order, and then, from (4.12), require that ( )2 2 0 1o xxA B? ? ??? >> for consistency. This implies the condition ( ) 2 2 32 2 8 10 '1 o o o ?? ? ? ? ? >> ? . (4.15) We interpret (4.15) in two ways. (i) If we fix a particular position xo and vary the driving frequency, the power spectrum at xo will peak around ~ ( )o ox? ? with the width of the resonance bounded by 1/3 2/3~ 'o? ? ?? . In contrast, a driven harmonic oscillator damped by a friction term tu? would give a resonance of width 1~? ? ?? ; dissipation due to phase mixing evidently yields a broader spectrum in terms of the dependence on the dissipation parameters. (ii) If we fix a particular driving frequency and consider the amplitude at steady state as a function of x, there is a corresponding spread. Suppose o? ?= at ox x= , i.e., there is a resonance around xo. We expand ( )o x? around xo, i.e., let os x x= ? , '( ) ...o o ox s? ? ?= + + ; then, we find the resonance to occur in a layer of width 1/3 1/3~ 'os ? ? ?? . Having obtained the above asymptotic information, we now solve (4.10) numerically, as an initial value problem, and look at the solution at times 2ct L ?>> , when the solution has come to a steady state. We use Mathematica 6.0 and solve (4.10) in slab coordinates, where 0.5 0.5x? ? ? ( i.e., 1xL = ). A linear profile (in x) of the GAM frequency is used, namely ( ) 1.2o x x? = + . The initial conditions are 59 ( , 0) 0u x t = = , ( , 0) 0tu x t? = = , and we use free slip boundary conditions at both ends. The profiles of the GAM frequency and the driving frequency are illustrated in FIG. 4.1. The oscillation in poloidal flow at steady state, for ( )2 40 ( 0) 10xx L? ? ?= = , 0.1od = , is illustrated in FIG. 4.2 as an example. FIG. 4.1. The profile of GAM frequency [ ( ) 1.2o x x? = + ] and the driving frequency ? along x, which will give a resonance around x = 0. FIG. 4.2. The steady state of (normalized) poloidal flow, will the GAM frequency and driver frequency profile indicated in FIG 4.1. Here ( )2 40 ( 0) 10xx L? ? ?= = and 0.1od = . Notice the resonance around x = 0 and the phase difference along x at any fixed time t. To verify (i), we drive the system with varying frequency ? and measure the corresponding steady state amplitudes at x = 0 (for a fixed ?). We then obtain a power 60 spectrum about the origin with ? as parameter, and measure the full width at half maximum (FWHM). FIG. 4.3 shows the spectrum for ? = 10-4 as an example. We obtain power spectra at different dissipation strength ?, and measure the corresponding FWHMs. The result is shown in FIG. 4.4 as a logarithmic plot, which is a straight line of slope 0.34 0.02? , in agreement with the expected value 1 3 predicted from the asymptotic analysis in (4.15). Moreover, a similar plot for the corresponding height of the resonant peak against ? is shown in FIG. 4.5: the measured slope is 0.334 0.002? ? . Note that a 1D driven harmonic oscillator would show a Lorentzian spectrum with height of the resonant peak 1~? ? . To study the effect of the GAM frequency gradient, we fix ? at 10-4, and use GAM frequency profiles ( ) 1.2o x mx? = + , where m models the inhomogeneity of the GAM frequency or sound speed. As above, we obtain power spectra for different values of m, and plot the corresponding widths (at resonance) against m (or 'o? ). The logarithmic plot is shown in FIG. 4.6, with a slope of 0.65 0.01? , in good agreement with the predicted value 2/3. This indicates a greater inhomogeneity can suppress the resonance by increasing the rate of the phase mixing along x. FIG. 4.3. Power spectrum measured at x = 0 with ? = 10-4, and ( 0)o x? ? ?? = ? = . 61 FIG. 4.4. Logarithmic plot of the spectral width a resonance against the viscosity measured at x = 0, the slope is 0.34 0.02? . FIG. 4.5. Logarithmic plot of the amplitude at resonance against the viscosity measured at x = 0, which gives a slope of 0.334 0.002? ? . FIG. 4.6. Dependence of width of peak of the power spectrum on the slope of GAM frequency. Viscosity is fixed at ? = 10-4. The slope of the graph is 0.65 0.01? . 62 Next, to verify (ii), we fix the driving frequency at ( 0)o x? ?= = , and plot the amplitude as a function of x, at steady state, with ? as parameter. We measure the (spatial) widths of the peak about x = 0 for different ?. The spatial dependence of the amplitude for ( )2 40 ( 0) 10xx L? ? ?= = is shown in FIG. 4.7 as an example, and a logarithmic plot of the width versus ? (FIG 4.8) shows that the spatial width of resonance ~ ??, where the measured 0.336 0.001? = ? , in agreement with the prediction of ? = 1/3. The dependence of the width of resonant flux surfaces on the slope of the profile is similarly investigated by fixing ? = 10-4, and measuring the (spatial) widths of amplitude in steady state for different values of '( 0)o x? = . The slope of the logarithmic plot is found to be 0.334 0.001? ? (FIG. 4.9), in agreement with the prediction in (ii). FIG. 4.7. The amplitude of (normalized) poloidal flow as a function of space at ( 0)o x? ?= = and ? = 10-4, 0.1od = 63 FIG. 4.8. Logarithmic plot of spatial width of resonance [at ( 0)o x? ?= = ] against the viscosity, with a slope of 0.336 0.001? . FIG. 4.9. Dependence of width of peak of the steady state amplitude on the slope of GAM frequency. Viscosity is fixed at ? = 10-4, and ( 0)o x? ?= = . The slope of the graph is 0.334 0.001? ? . The above results show that, in the presence of phase mixing, the dissipative effects in an inhomogeneous driven system are enhanced. The amplification attained in driving a GAM at resonance in an inhomogeneous medium is significantly lower than that in the 1D harmonic oscillator (or in homogeneous medium). It is also more difficult to get a sharp resonance (width of resonant peak ~ ?1/3 instead of ~ ?) , and it is relatively harder to localize the driven modes at a desired flux surface. 64 Chapter 5: Nonlinear effects on resonantly driven GAMs 5.1 Introduction In this chapter, we retain the nonlinear terms in {(3.12), (3.13), (3.14)} and analyze their effect on GAMs. The nonlinear terms arise from the convective terms in the system. We do not include herein the viscous dissipative effect, as we have already dealt with this in the previous section. Our objective is to present a clear illustration of nonlinear effects. Under these conditions, the terms in our set of equations are free of spatial radial derivatives, so we can simply consider GAMs on one particular flux surface at r. We are interested in GAMs that are externally driven by a momentum source. Thus, we add a momentum source in (3.14): 22 ( )sE s cdu N K t dt R = ? + . (5.1) With the normalization 2E su u c= and ( )' 2 st c R t= , and defining ( )* 2( ') ' 22 ss RK t K Rt c c = , the set {(3.12), (3.13), (5.1)} becomes 0 ' s c dN R uN u dt r ? ? = , (5.2) 65 0 ' c s dN R uN dt r + = , (5.3) * ( ') ' s du N K t dt + = . (5.4) In the next section of this chapter we will try to solve the system {(5.2),(5.3),(5.4)} by linear expansion, with the amplitude of driving terms as the small expansion parameter. We show that the expansion scheme breaks down at the vicinity of resonance and obtain the condition with asymptotic analysis. And then we solve the equation nonlinearly and we show that the nonlinearity has a suppression effect on the resonance. Lastly, we verify our results with a numerical simulation. 5.2 Linear expansion of nonlinear GAM equations with small driver amplitude If we consider the homogeneous problem, i.e., * ( ') 0K t = , and examine {(5.2), (5.3), (5.4)} as an initial value problem, the nonlinear effects will be very small for small amplitude perturbations, and the problem can be approximated by the linearized equations. However, for the driven problem, ( ) sinoK t a t?= (ao is a constant), the nonlinear terms become important when the driver frequency approaches the GAM frequency, even at small amplitudes. To illustrate this, we expand u , sN , cN in (5.2), (5.3) and (5.4) in orders of ~ 1oa ? << according to: 1 3 5 1 3 5 2 4 6 ... ... ... s s s s c c c c u u u u N N N N N N N N = + + +?? = + + +? ? = + + +? . (5.5) To the order ? , we have the linear driven GAM equations and the solution is 66 * 1 2 * 1 2 cos ' 1 sin ', 1s a u t aN t ? ? ? ? ? ? =?? ?? ? =? ?? (5.6) where 2 * 2o sa Ra c= . To the order 2? , a nonlinear term drives 2cN , 2 1 1 0 ' c s dN R u N dt r + = , (5.7) which gives ( ) 2 2* 2 22 cos 2 ' 4 1 c o aRN t C r ? ? ? ? ? ? ?= +? ? ?? ? . (5.8) To determine the integration constant 2 oC? , we eliminate u in (5.2) by multiplying (5.2) with sN and substitute (5.3) into the resulting equation, which gives ( )22 2s cN N r R D+ + = , (5.9) where D is a constant. Taking (5.9) to lowest order, we get 2 2 2D r R= . Then considering (5.9) again to order 2? (the equations at order 1? are trivial), and using the expression in (5.6) for 1sN , we get ( ) 2 2 * 224 1 o aC? ? = ? ? . To order 3? , by differentiating (5.4) with respect to 't and using (5.2) to eliminate ' sdN dt in the resulting equation, we see that 3u is driven by the nonlinear term 1 2cu N : 2 3 3 1 22 ' c d u R u u N dt r + = ? , (5.10) and the corresponding solution is 67 ( ) ( ) 2 3 * 3 3 2 22 2 3 2 * 3 3 2 22 cos ' cos3 ' 1 1 98 1 sin ' 3sin 3 ' . 1 1 98 1 s aR t t u r aR t tN r ? ? ? ? ?? ? ? ? ? ?? ? ? ? ? ? = ?? ? ? ? ? ? ?? ? ? ? ???? ? ? ? ?? = ?? ? ? ?? ? ?? ? ? ? ??? (5.11) Since the higher order terms are driven by nonlinear terms, nonlinearities are important as ? approaches unity. As 1? ? , we see that this gives 3 1u u? , indicating that the approximation breaks down. This happens for ( ) 2 3 * * 4 22 18 1 a aR r ? ? ?? ? ?? ? ?? ? ? ? , (5.12) or equivalently, 1/34 2 2 41 ~ o s R a r c ? ? ? ? ? = ? ? ?? ? . (5.13) Recall that for a driven harmonic oscillator without damping, the amplitude goes as 1 ?? and can become infinite at resonance. If there is a small damping term tu? , the amplification at resonance will be finite and goes as 1 ? . Here, we have shown that the usual resonance is altered by nonlinearities beyond a critical smallness for ?? . To understand the behaviour of the nonlinear GAM near resonance, we notice first that for the homogeneous problem ( * ( ') 0K t = ), if we eliminate sN and cN in {(5.2), (5.3), (5.4)}, u satisfies the homogeneous Duffing equation [45,46]: 22 3 2 1 0 ' 2 d u R u u dt r ? ? + + =? ?? ? . (5.14) Therefore, we expect and will show that the solution to the driven problem is parallel to that of the driven Duffing equation. For illustrative purposes, we present the 68 derivation of the solution to the inhomogeneous Duffing equation (with sinusoidal driver, to lowest order in ?? ) in Appendix C for comparison. 5.3 Analysis of the nonlinear GAM equations at resonance We are interested in driving GAMs at resonance, where the nonlinear effects become important. Therefore, in (5.4) we consider a sinusoidal driving term of frequency close to resonance [46]: ( ) ( )2* * *( ') sin 1 ' sin 'K t a t a t? ? ? ?? ?= + ? = + ?? ?? ? ? ? , (5.15) where 1? ?? = ? and 2~ 1? ?? << . Also, we introduce a slow time scale 2 't? ?= , and solve {(5.2), (5.3), (5.4)} with the orderings 1 1~ ~su N ? , 22 ~cN ? , 3* ~a ? , where we replace 'd dt with 2 't ??? + ? in this two time scale treatment. To the order ? , the solution is 1 1 ( )sin ( ', ) ( ) cos ( ', ) ( ) sin ( ', ) ( ) cos ( ', ),s u A t B t N B t A t ? ? ? ? ? ? ? ? ? ? ? ? = +?? = ?? (5.16) where ~ ~A B ? and ( )2( ', ) 't t? ? ? ? ?= + ? . To the order 2? , from (5.3) we get 2cN [ ] [ ] 2 2 2 2 1sin 2 ( ', ) cos 2 ( ', )2 4c R AB B AN t t C r ? ? ? ? ?? ??= + +? ?? ? . (5.17) Again, using (5.9) and the expression for 1sN in (5.16), we get ( )2 2 21 4C A B? = ? + . To the order 3? , the driving term appears; in this order, (5.2) and (5.4) become 2 ' 3 3 1 1 2t s s c RN u N u N r ??? ? = ? ? ? , (5.18) 69 ( )2 2' * sin 't su N u a t?? ? ? ?? ?? + = ? ? + + ?? ? . (5.19) After combining (5.18) and (5.19), we have ( ) ( ) ( ){ ( ) ( )[ ] } 2 2 2 ' 3 3 2 2 2 2 2 2 2 2 2 2 * 2 sin ( ', ) 2 cos ( ', ) cos ( ', ) 2 sin ( ', ) 4 8 sin ( ', ) 4 8 cos ( ', ) sin ( ', ) cos ( ', ) 4 cos ( ', ). t u u B t B t A t A t R AB A B A t r A B B B A t A B A t B t a t ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? + = + ?? ? ? + ? ? ? ? ? ? ? ? ?? ? ? ?? ? ? ?+ ? ?? ? ? + + + (5.20) [In (5.20) we have dropped terms with dependencies of sin 3 ( ', )t? ? and cos3 ( ', )t? ? , stemming from 1 2cu N in (5.18), as they represent driven response at a frequency far off resonance compared with the resonant terms.] If A(?) and B(?) were arbitrary, 3u from (5.20) would be a linear combination of sin ( ', )t? ? and cos ( ', )t? ? , and the ordering would break down as ? approached 1, as in (5.5)-(5.14). Therefore, we seek A(?) and B(?) such that R.H.S. of (5.20) is identically zero, i.e., we remove the secularity. The corresponding conditions are ( ) ( ) ( ) ( ) 22 2 2 22 2 2 * 2 2 8 0 2 2 8 . B A R r A A B A B R r B A B a ? ? ? ? ? ? ? + ? + + =?? ? + ? + + = ??? (5.21) This yields 0A = , 0B? = , and B is found to satisfy ( )2 3 *16 8 0R r B B a?+ ? + = . (5.22) Define 2 * 2 s o cB B a R = and rewrite (5.22) as 70 3 * * 64 32 0B B? ?+ ? + = , (5.23) where the parameter 2 4 2 4 o s a R r c ? = characterizes the strength of nonlinearity. It can then be seen from (5.23) that when ? << 1, * ~ 1B ?? , which recovers the behavior of the undamped driven oscillator. Furthermore, we note that the normalized amplitude for * B effectively represents the amplification of B . FIG. 5.1 shows a plot of * B versus ?? from (5.23) at 0.2? = . The resonant peak shifts to lower frequency (with respect to the natural GAM frequency), and the width of resonance is ( )10/3 1/3 1/33 2 ~ 0.3c? ? ?? = . Here c?? is the critical value of ?? which gives one real root and two repeated real roots for * B in (5.23). Note that 1 c?? increases with the sharpness of resonance, as indicated in FIG. 5.1. Moreover, the maximum amplification is bounded by the height within this layer 7/3 1/3 1/3 *max 2 ~ 5B ? ? ? ?? (outside 1/3~ 0.3? ?? , the ordering in this calculation does not apply, but the linear expansion (5.5)-(5.13) apply and * ~ 1/B ?? ). Also note that the width of the resonance is consistent with the prediction from linear expansion as given in (5.13). 71 FIG. 5.1. The amplitude of u as a function of driving frequency (around the natural GAM frequency) for nonlinear GAM equations (? = 0.2). Notice that the width of resonance is finite and the amplification is bounded even no dissipative terms are included in the equations. Next, to verify the above analysis, we solve {(5.2), (5.3), (5.4)} numerically with Mathematica 6.0, as an initial value problem. In our simulation, we use a sinusoidal driving term ( ) ( )2*( ') 2 sin 1 'o sK t Ra c t?= + ? , and a damping term u?? is added to the RHS of (5.4), so that the system can relax to steady state at ' 1t? >> . We choose 2 sc R? << so that this damping term will have little effect on the steady state. Initially we set all variables to be zero. In the following, we set ( )0.002 2 sc R? = , 1sc = , 10R = , and the amplitudes of the solutions are all measured at 4 ' 7.5 10t = ? . To get a power spectrum similar to the one shown in FIG. 5.1, we fix a particular flux surface r and driver amplitude oa , drive the system with frequency 1 ?+ ? , and then measure the amplitude of u at steady state. We repeat this iteration with different ?? and plot the amplitude (normalized to give the amplification) versus driving frequency. FIG. 5.2 shows the power spectrum at 3r = and 0.01oa = . The spectrum shows a sharp change in amplitude at c? ?? = ? : this is as predicted in FIG. 5.1. Also, the maximum amplification maxAmp is less than 12, in our simulation, 72 while linear theory would predict ~ / 500oAmp ? ? = , suggesting a strong suppression of resonance because of the nonlinearity. FIG. 5.2. The power spectrum of nonlinear GAM oscillator. Notice the sharp decrease in amplitude when ?? drops below c?? , a property inherited from Duffing oscillator. We are also interested to verify the dependencies of the resonant width c?? and maximum amplification maxAmp on the nonlinear parameter ? predicted by (5.23) (i.e., FIG. 5.1). To check this, we i) iterate the system and get power spectra for different flux surfaces at fixed ( 0.01)oa = (we use 1sc = , independent of r), and plot c?? and maxAmp as a function of r; and ii) similarly obtain the spectra for different oa , at a fixed flux surface ( 3)r = , and plot c?? and maxAmp against oa . The results are shown in FIG. 5.3 ? FIG. 5.6. From our calculation [(5.23) and FIG. 5.1], we predict 1/3~c? ?? , 1/3 max ~Amp ? ? , where 2 2~ oa r? ? ; our simulation results shown in FIG. 5.3 ? FIG. 5.6 indicate 0.686 0.003 0.65 0.01~c oa r? ? ? ?? and 0 .69 0.01 0.73 0 .01 m ax ~ oAm p a r ? ? ? , which are in good agreement with theory. 73 FIG. 5.3. Logarithmic plot of resonant width against flux surface r at ao = 0.01, the slope is 0.65 0.01? ? . FIG. 5.4. Logarithmic plot of maximum amplification against flux surface r at ao = 0.01, with a slope of 0.73 0.01? . FIG. 5.5. Logarithmic plot of resonant width against driver amplitude at r = 3, the slope is 0.686 0.003? . 74 FIG. 5.6. Logarithmic plot of maximum amplification against driver amplitude at r = 3, with a slope of 0.69 0.01? From the theory of the 1D linear harmonic oscillator, the amplification is o? ? and the spectral width at resonance ~ ? (Lorentzian). Both are independent of the driving amplitude oa . However, in the nonlinear GAM equations, even in the undamped case, the amplification is limited, and both the amplification and sharpness of resonance depend on oa . The parameter ? increase as 2 oa ; thus increasing oa decreases the amplification and broadens the peak around resonance. This nonlinear detuning effect makes it relatively harder to drive GAMs in a tokamak. 75 Chapter 6: Implications for experiments to drive GAMs in tokamaks 6.1 Introduction In this Chapter, we use tokamak parameters to assess the efficacy of resonantly driving GAMs. We evaluate the power requirement for resonanting GAMs by assuming the dissipation is from classical collisions and evaluating the enhanced damping from phase mixing in Sec. 6.2. Also, we evaluate the strength of the nonlinearities to investigate their significance in tokamaks. Hallatschek and McKee [5] have proposed the driven GAM experiment for DIII-D in Sec 6.3. Thus, we will take the DIII-D tokamak as an example, and we will consider both the outer edge region and the core region. In the former case we take the temperature 200T eV= , and we use 1T keV= for the latter. The number density is typically 19 310n m?= , major radius of the tokamak 1.5R m= , minor radius 0.5a m= , and toroidal field 2B T= . Deuterium plasma is used. Finally we will summarise our study of GAMs, i.e. the second part of this dissertation (Chapter 3 to Chapter 6), in Sec. 6.4 6.2 Effect of phase mixing on damping rates and power requirements GAMs could be excited in tokamaks experimentally by, for example, external time varying magnetic fields created by coils which drive an axisymmetric oscillation of flux tubes in the tokamak. The temperature and, so, the GAM frequency, decreases 76 radially outwards. As a result, the magnetic flux surfaces on which GAMs are to be excited can be chosen by the driver frequency ? of the external field. Consider that we excite GAMs at some flux surfaces at ~ or r , i.e. ( )G or? ?= , with a corresponding radial spread x? . The power required is given by 2 max 2P mu?= , (6.1) where ? is the damping rate, m is the mass of the plasma within the resonantly driven magnetic surfaces, and maxu is the amplitude of the (oscillatory) poloidal flow at resonance. Such an oscillation gives a time-varying velocity shear of amplitude max' ~V u x? . The idea is to excite GAMs to a large amplitude such that the associated shear flow is comparable to the growth rate of the plasma turbulence, so as to obtain a significant stabilization effect, i.e., max' ~ growthV u x ?? > . We assume the growth rate of the instability to be suppressed to be of order the drift frequency, * ~growth i s nk c L?? ? ?= , which is expected to be the maximum of the instability growth rate. Here k? is the poloidal wavenumber of the drift mode, i? is the ion Larmor radius and nL is the radial scale length of the density profile. In the following, we assume ~ 1ik?? , where maximum growth is expected, and ~nL a . This gives an estimate of the required poloidal flow, max ~ su c x a? . The spatial spread x? is related to the width of frequency spectrum as ~ ' ~G Gx a a Q? ? ? ?? = , with the Q-factor GQ ? ?= . Substituting the above scalings into (6.1), and taking ( )2 22 2i im NM nM Ra x api= = ? , we get ( )42 2~ 4 G GP nRa Tpi ? ? ? . (6.2) 77 For collisional damping, the damping rate is known to be ~col ii? ? , the ion-ion collision frequency, which is proportional to 3/2T . This is the rate expected from magnetic pumping. On the other hand, the effective damping rate associated with phase mixing is given by ( ) 1/321/32 2 2 2~ ' ~ i ph i ii G ii G a ?? ? ? ? ? ?? ?? ?? ? [from (4.7)]. We evaluate damping rates and the corresponding Q-factors, spatial spreads and power requirements [from (6.2)], for both the edge ( 200T eV= ) and the core ( 1T keV= ). We summarize our results in Table 6.1. Table 6.1. Estimated damping rates, Q-factors, spatial spreads and power requirements, by classical collision and phase mixing, at the edge and core of tokamak. T 200eV (edge) 1keV (core) Collisional ( )ii? Phase mixing ( )ph? Collisional ( )ii? Phase mixing ( )ph? ( )1s? ? 1800 640 180 860 Q 71 200 1600 330 x? (cm) 0.71 0.24 0.031 0.15 P (W) 25 0.36 0.0010 0.53 From Table 6.1, in the edge region, collisional damping is more important than that from phase mixing. The power requirement for the former is 25W, with a spatial spread of 0.71cm: these are very favorable conditions. For the core, damping by phase mixing is more important, as the power dissipation by collisional damping drops rapidly as temperature increases. We find that the power requirement for phase mixing is less than 1W with a spread less than 1cm, which is very favorable. We note that the scaling of T is important in our estimation, since 5/3~ii ph T? ? ? , and dissipation from phase mixing dominates towards the core. Also, since 78 13/2 ~collisionP T ? and 1/6~phP T , dissipation by collision damping decreases strongly as temperature rises while that of phase mixing is relatively insensitive to temperature change. Therefore, if the power requirement is low in the edge, it may still be low in the core. Similar estimates for power requirements for excitation of GAMs in the edge of the DIII-D tokamak were also done by Hallatschek and McKee. They found a power requirement of 700W. Instead of using ion-ion collision as the basis of damping, they assumed a Q-factor of 20 obtained from the turbulence simulation, for their estimates. Thus their Q-factor is about 1/3 of ours at the edge (based on ii? ). Their smaller Q- factor is possibly due to turbulence broadening. They also used 2x cm? = , apparently estimated directly from the simulation. This gives them a higher estimate of power dissipation. Notice that a wider spatial spread x? implies a thicker transport barrier for the same flow shear, which leads to enhanced confinement. However, as discussed above, a large x? corresponds to a small Q-factor thus the resulting power dissipation [from (6.2)] is higher. In our above calculations for phase mixing, we have used ~nL a , where a = 50cm. This estimate is applicable to the core plasma and also, generally, to edge plasma in L-mode. L-mode in tokamaks refers to shallow gradients at edge, whereas H-mode refers to much steeper gradients wherein nL and TL can be of order 5cm. We have calculated x? and power P also for these steep gradient situations. We find that phase mixing dominates, with a required power of 17W. However, the corresponding ~ 0.1x cm? , which is of order the Larmor radius i? . At these short distances, the effect of finite Larmor radius on velocity shear stabilization should be examined. 79 6.3 Comparison of the effect of nonlinearities and the effect of damping In Chapter 5, we show that the nonlinear convective terms in the nonlinear GAM equations can lower the amplification of resonantly driven GAMs as the driver amplitude oa increases. On the other hand, we know that in the linear model, the amplification is limited by the damping coefficient, with the amplification going as ~ G? ? . To determine whether nonlinearities are more important, we investigate the nonlinear system {(5.2),(5.3),(5.4)} in Sec. 5.3 again, and add a damping term ( )G u? ? on the L.H.S. of (5.4), and then compare the size of the dissipative term and the nonlinear terms in the equations. We assume 2~G? ? ? in the expansion (5.15)-(5.23); then the system can be solved in a similar manner as in Sec. 5.3. (We did not solve the nonlinear system with damping in Sec. 5.3, so as to illustrate the nonlinear effect clearly.) In the presence of the damping term, similar to the linear harmonic oscillator, the solution is not exactly in phase with the driving term, i.e. ( )A ? in the ansatz (5.16) is no longer zero, in fact, ( ) 3~ ~GA ? ? ? ? ? . In this case, it is more convenient to write sinA X ?= and cosB X ?= and solve for the amplitude ( ); ,X ? ? ? and the phase difference (from the driver term) ( ); ,? ? ? ? respectively. Repeating the expansion (5.17)-(5.21) with the introduction of the damping term, we get ( ) ( ) 1/22 22* * *64 32 32 0X X? ? ?? ?+ ? + + =? ?? ? , (6.3) 80 where 2 * 2 s o cX X a R = and * G? ? ?= , and ( ) * 22 2 * * sin 2 32X ?? ? ? ? = ? + + . (6.4) Note that in the limit * 0? ? , the phase ? goes to zero and (6.3) reduces to (5.23). On the other hand, in the limit 0? ? , (6.3) implies ( ) 1/22 2* *4X ? ? ?? ? + , which is the Lorentzian spectrum for a linear harmonic oscillator. We now use (6.3) to show that the nonlinear effect is self-consistently small compared with the dissipative effect. From (6.3), assume 2 * * 32X? ?lessmuch and for 0?? = (i.e., at resonance), * * ~ 1X ? . Then the consistency condition needed to neglect nonlinearities is 3 * 32? ?lessmuch . (6.5) Equation (6.5) can also be interpreted heuristically: we observe from FIG. 5.1 that the nonlinear effect limits the amplification to 1/3~ ? ? , and we know that in the harmonic oscillator the amplification at resonance is G? ? . Therefore, if the latter is more important, which makes the excitation more expensive, we will have 1/3G? ? ? ?<< , which is equivalent to (6.5). We recall 4 2 2 4 o s R a r c ? = , and use 2o P a m ? = for the driver amplitude: then, in the edge region ( 200T eV= ), where classical collisional damping is more important, we find the ratio of nonlinear effects to damping effects to be 3 * : 32 1: 250? ? = . For the estimation by Hallatschek and McKee, which assumed a turbulence scenario with 81 20Q = and 700P W= , the ratio is 1: 390 . In the core, ( ~ 1T keV ), where damping from phase mixing dominates, the ratio becomes 1: 640 . Therefore, we conclude, from the self-consistent check, that nonlinear effects are not important for the experimental parameters considered. 6.4 Summary of Chapter 3 to Chapter 6 In this section, we give a summary of the second part of this dissertation from Chapter 3 to Chapter 6. In summary, we have investigated two factors which may affect the efficacy of driving GAMs resonantly in tokamaks, namely the damping enhancement from phase mixing, and the detuning of resonance by nonlinearities. We began with describing the physics of GAMs from MHD equations and derive a closed set of nonlinear reduced equations as the basis of our analysis. We then linearize the reduced equations, and solve the homogeneous problem, to establish the damping factor ( )3exp ot t? ??? ? , with 1/3~ot ? ? and ? the coefficient of viscosity. This decay is faster than exponential and is a characteristic factor for phase mixing damping, as also shown by Heyvaerts and Priest for Alfven waves. We also study the driven problem, by asymptotic analysis and simulation. From both approaches, we find a spatially localized steady state, with a spatial spread of 1/3 1/3~ 'os ? ? ?? , and the peak of the resonant 1/3max ~u ? . This is broader than that of a typical harmonic oscillator model, implying that phase mixed damping increases the challenge to drive GAMs in a tokamak. 82 We then consider the driven nonlinear non-dissipative equations, by an asymptotic analysis and an expansion in small driving amplitudes in the vicinity of resonance. We find that the system is similar to the Duffing oscillator, with the parameter 4 2 2 4 o s R a r c ? = characterizing the strength of nonlinearity, and we find that the peak and width of the resonance are 1/3~ ? ? and 1/3~ ? respectively. This solution is consistent with our simulations. Thus, the resonance will be suppressed if nonlinearity is important. Finally, we estimate the power requirement for resonating GAMs in a tokamak, using typical experimental parameters for the DIII-D tokamak. We compare our numbers with a similar estimation by Hallatschek and McKee, who assume turbulent scale sizes to estimate the spatial localization. We find that classical collisions dominate at the edge, and that the corresponding spatial spread of GAMs is 3 times smaller than Hallatschek?s estimates. In the core where the temperature is higher, we find phase mixing is more important. In both cases, however, the power requirement is small, thus advocating for the proposed experiment. We also estimate the nonlinear Duffing type terms for the operating parameters of the experiment, and find that nonlinearities will not be important compared with the damping effects. 83 Chapter 7: SUMMARY and CONCLUSION In this dissertation we have considered some aspects of velocity shear stabilization of magnetized plasma instabilities. In the first part, steady externally forced flow shears are considered. In the second part, resonantly excited oscillating flow shears are considered. In the first part of the dissertation, we perform 2D simulations to explore the stabilizing effect of flow shear on ideal interchange modes, in both linear and nonlinear regimes. We show that a cross-field velocity shear V? suppresses growth rates over a broad spectrum of wavenumbers in linear regimes. The effect is more significant for modes with either very short or very long wavelengths, leaving the unstable region in the k? spectrum to condense and peak at 6ym = , as V? increases. This is in contrast with the shear-free spectrum, which gives k? ?? at long wavelengths and asymptotes to a constant in the short wavelength limit (ignoring viscous damping). We also find that the spectrum in the nonlinear turbulent steady state is similar to that of the linear spectra: the peaking wavenumber corresponds to roughly circular convection cells of size comparable to the density profile scale length. A spatial Fourier transform of residual fluctuations at different shear strengths indicates that the spectrum peaks at 6ym = . Moreover, both the average density fluctuation and the fractional flattening of the density profile in turbulent steady state decrease as V? increases. Complete stabilization is shown to be possible when V? is greater than a few times the classic interchange growth rate, wherein we observe that the original (RT unstable) density profile is recovered and the density and flow 84 profiles revert to laminar. In addition, we measure the particle flux in turbulent steady state and compare this with the classical diffusive transport. We show that the latter eventually dominates and accounts for all the transport as V? increases, again showing complete stabilization in the strong shear limit. We compare our simulation results with magnetic fluctuation measurements from MCX. This experiment studied interchange modes in the presence of flow shear. We show that the experimental wavelengths measured agree with simulation results if done at higher viscosity. We confirm that this level of viscosity is consistent with the Reynolds numbers at the edge of the plasma in MCX, where the measurements are taken. Both the simulation and the data are consistent with elongated convection cells in the nonlinear regime. Finally, we investigate the effects of Kelvin-Helmholtz instabilities on the system; the latter instability must be present in any flow shear system. We show that we can stabilize an RT unstable configuration with flow shear only when the generalized Rayleigh Inflexion criterion is simultaneously satisfied. This restriction still allows a broad parameter space in which the RT mode is effectively suppressed. In the second part of the dissertation, we study some factors which may affect the driving of Geodesic Acoustic Modes (GAMs) at resonance to establish significant oscillating flow shear stabilization in tokamaks. These factors are (1) enhanced damping from phase mixing, and (2) the detuning of resonance by nonlinearities. We first investigate phase mixing in Chapter 4, and find that the damping of the mode is faster than exponential, with an effective time constant 1/3~ot ? ? , where ? is the coefficient of viscosity. We then study the driven problem in Sec 4.3, which resembles the situation of a proposed experiment to drive GAMs in tokamaks. We 85 find a spatially localized steady state, with a spatial spread of 1/3 1/3~ 'os ? ? ?? , corresponding to a resonant peak of 1/3max ~u ? , which is broader than the standard harmonic oscillator model. Thus phase mixed damping may broaden the resonance of GAMs in tokamak. We next study nonlinear detuning of GAMs, by nonlinear perturbation theory, and find that the system is similar to the Duffing oscillator, with the parameter 4 2 2 4 o s R a r c ? = characterizing the strength of the nonlinearity. We show that the peak and width of the resonance scale as 1/3~ ? ? and 1/3~ ? , respectively. We show that our analytic calculations are consistent with simulation. Thus, GAM resonance could be suppressed if the nonlinearity is important. Finally, we estimate the power requirements to resonate GAMs in a tokamak for experimental parameters in the DIII-D tokamak. We consider both the collisional and the phase mixing models, and compare with similar estimations done by Hallatschek and McKee which assumed turbulent broadenings. We find that classical collisions dominate at the edge, and the corresponding spatial spread of GAMs is ~1/3 of Hallatschek and McKee?s results. In the core of tokamak, where the temperature is higher, phase mixing is more important. In both cases, the power requirement is small and thus the proposed experiment is cost effective. We also estimate the nonlinear Duffing terms in the operating range of the experiment and find that nonlinearities are not important compared with the damping effects. 86 Appendices A. Alfven waves propagating in an inhomogeneous, dissipative medium In this Appendix, we review the phase mixing of Alfven waves. We consider Alfven waves propagating in an inhomogeneous medium, with the inhomogeneity perpendicular to the magnetic field and the polarization of the Alfven wave. This calculation was done by Heyvaerts and Priest [9]; here we give a slightly different derivation with assumptions that make the problem simpler but still serve our illustrative purpose. We show that the decay from viscous and resistive dissipation is faster than exponential, with the decay factor 3exp ( / )ot t? ??? ? , where the time constant 1/3 ~ot ? ? (or 1/3~? ? ), ? and ? are coefficient of viscosity and resistivity of the plasma respectively. This is another example to illustrate the enhanced dissipative effect due to phase mixing in addition to GAMs in a tokamak, discussed in Sec. 4.2. We then verify the results by simulation in the second part of this appendix, for both the undriven and driven Alfven system. A.1 Equations Consider a magnetized plasma with magnetic field 0 ?( )zB B x z= harpoonrightnosp (rectangular coordinates). Assume zB B? << harpoonrightnosp . The perpendicular components of the momentum equation and Faraday?s law are, respectively, 87 2 ( ) 8o u B n u u nT n u t ? pi ? ? ? ? ? ? ? ??? ? + ?? = ?? + + ?? ?? ?? ??? ? ? ? harpoonrightnosp harpoonrightnosp harpoonrightnosp harpoonrightnosp harpoonrightnospharpoonrightnosp harpoonrightnosp harpoonrightnosp , (A.1) ( ) ( ) 2 24t z z z cB u B u B B u B? pi? ? ? ? ? ? ? ? ? = ? ? ? ? ? ? + ? harpoonrightnosp harpoonrightnosp harpoonrightnosp harpoonrightnosp harpoonrightnospharpoonrightnosp harpoonrightnosp , (A.2) where ? is the coefficient of viscosity and ? is the resistivity of the plasma, and we assume the system is isothermal and ion mass 1M = . Consider a 1D equilibrium given by 2( 8 ) 0x o znT B pi? + = , (A.3) and u? harpoonrightnosp and B? harpoonrightnosp are zero at equilibrium. Now consider perturbations of u? harpoonrightnosp and B? harpoonrightnosp in the y direction, i.e., only consider Alfven wave of y-polarization. (A.1) and (A.2) become ( ) ( )t y x x y o z ynu n u B B?? ? ? ? = ? tildenosptildenosp tildenosp , (A.4) 2 2( / 4 )t x y o z yc B B u? pi? ?? ? ? = ?? ? tildenosp tildenosp , (A.5) where ~ 1y A y zu B B <> ?tildenosp tildenosp and drop the latter term in (A.6). We will check that the solution so obtained will be consistent with this approximation for 2 zt kpi>> Av . We then get a system of coupled wave equations: 88 2 2 [ ] .[ ] t x z y t x y z w B B w ? ? ?? ? ? = ? ??? ? ? = ? ?? tildenosptildenosp tildenosp tildenosp A A v v (A.8) We now solve (A.8) as an initial value problem with the ordering 2 2 ~ 1x t At L? ? ?? ? = <> ?tildenosp tildenosp . From (A.14) again, to lowest order of ? , we have ( )2 2 2 20 0~ ' ''x w k t ik t w? ? +tildenosp tildenospA Av v . If ( ) ( )2 2' ~ ' ~ "n n n nA Av v , the effect of the terms we neglected is small compared with dissipation when 2 zt kpi>> Av . From (A.14), we note that the decay rate of the Alfven wave in the presence of phase mixing is of the form ( )3exp ot t? ??? ? , which is faster than the usual exponential decay, with the time constant 1/3~ot ? ? . This has been shown by Heyvaerts and Priest, except they assumed a solution of constant frequency (in space) and found that the wave decays as ( )3exp z L? ??? ? , where z is the distance from the origin. In a similar way, we can obtain standing Alfven wave solutions when z is confined in an interval from 0 to Lz instead of whole space. Consider the boundary conditions 90 ( 0) ( 0) 0 ( ) ( ) 0 z y z z y z w z B z w z L B z L ?= = ? = = ?? = = ? = = ?? tildenosptildenosp tildenosptildenosp . (A.16) By superimposing Alfven waves travelling in opposite directions in (A.14), we find the standing waves solutions ( ) ( ) 2 2 3 0 0 1 sin 1 " 2 sin( )' exp 3 cos 1 " 2 cos( ) l ll y l l k t t k zw k t B k t t k z ?? ? ? ?+? ?? ? ? ? ? ?? ?= ?? ? ? ?? ? ? +? ?? ?? ? ? ?? ? tildenosp tildenosp A A AA A A A v v vv v v v , (A.17) ( ) ( ) 2 2 3 0 0 2 cos 1 " 2 sin( )' exp 3 sin 1 " 2 cos( ) l ll y l l k t t k zw k t B k t t k z ?? ? ? ?+? ?? ? ? ? ? ?? ?= ?? ? ? ?? ?+? ?? ?? ? ? ?? ? tildenosp tildenosp A A AA A A A v v vv v v v , (A.18) where 2l zk l Lpi= , and the same consistency condition (A.15) applies. A.2 Numerical simulations for undriven and driven standing Alfven waves In this section, we use a 3D dissipative MHD code (developed by Guzdar et al [39]) to illustrate the enhanced dissipation in the presence of phase mixing for standing Alfven waves. Our simulations are based on (A.1) and (A.2) with the equilibrium given by (A.3). The size of the simulation box is x y zL L L? ? , with 1x yL L= = , 5zL = , and the grid size is 401 7 31x y zN N N? ? = ? ? . The viscous coefficients [normalized with ( 2)x A xL x L=v ] are 510x y? ? ?= = , 310z? ?= , and the resistivities [normalized with ( 2)x A xL x L=v ] are 2 2 2 2 54 4 10x yc c? pi ? pi ?= = , 2 2 34 10zc? pi ?= . We use an equilibrium density profile 0.11 cos x x n L pi pi = + , and 10oT = , so that the magnetic field at equilibrium is given by (A.3). We use periodic 91 boundary conditions for all variables at both ends of the y and z direction, and free slip boundary conditions for density and magnetic field, fixed end boundary conditions for the flow, at both ends of the x direction. The equilibrium profiles of density and magnetic field are shown in FIG. A.1, and the corresponding profile in Alfven speed is shown is FIG. A.2. FIG. A.1. The density and magnetic field profile at equilibrium. FIG. A.2. The Alfven speed profile corresponds to the profiles in FIG. A.1. A.2.1 Initial value problem for undriven system We are going to simulate standing Alfven waves along the z and polarized in the y direction, and investigate the decay. The initial conditions are ( 0) 0ynu t = =tildenosp , 92 ( )( 0) ( 1 2) 0.01cos 2y o zB t B x z Lpi= = =tildenosp , i.e., z zL? = . The time evolution of yB at the center of the box is shown in FIG. A.3, and the corresponding plot using uniform density profile (i.e., constant Av ) is shown on the same graph for comparison. The figure indicates a rapid decay when the profile is nonuniform, while for the uniform profile, the time constant for dissipation ( ) 12 21/ ~ 4 633e z z Ac k t? ? pi ? = . FIG. A.3. Time evolution of Alfven wave at the center of the system showing dramatic decay when phase mixing occurs. The corresponding oscillation using a flat profile is shown for comparison. Furthermore, we can measure the decay from the simulation quantitatively. Consider our equation set (A.8) without dissipation, i.e., 0? = . We multiply the first and second one with wtildenosp and yBtildenosp respectively, and sum the two equations; we get ( ) ( )2 212 t y z yw B wB? + = ?tildenosp tildenosptildenosp tildenospAv . (A.19) By averaging (A.19) along the z direction, the R.H.S. of (A.19) vanishes and we define the quantity 2 21( , ) ( , ) 2 y y z E x t nu B x t= + tildenosptildenosp , (A.20) 93 which is time independent for any particular x plane when the medium is non- dissipative. Theoretically, from (A.17), ( , )E x t decays as ( )3exp 2 ot t? ??? ? , where ( )1/32 23 'o lt k?= Av . Therefore, we evaluate log ( , )E x t (at the mid-plane x = 1/2) from our simulation and make a plot of log ( 1/ 2, )E x t= versus 3t is shown in FIG. A.4. The curve is a straight line with a negative slope for 140[ / ( 1/ 2)]x At L x? =v , and levels off after that. The slope of the linear portion is ( ) 5 31.20 0.01 10 [ / ( 1/ 2)]x AL x? ?? ? ? =v , which is in good agreement with the theoretical value 5 31.09 10 [ / ( 1/ 2)]x AL x? ?? =v . The curve levels off as the consistency condition (A.15) breaks down. To illustrate this, we substitute the parameters used in our simulations into (A.15), and we get 3270[ / ( 1/ 2)]x At L x ?<< =v for (A.17) to be valid, thus the fact that our simulations show that (A.17) is no longer valid in describing the solutions to the system for 140[ / ( 1/ 2)]x At L x> =v , is an expected result. 94 FIG. A.4. A plot of log ( , )E x t versus 3t , which shows the decay goes at ( )3exp 2 ot t? ??? ? before the consistency condition breaks down at ~ 140( / )x At L v and the curve levels off. The parameters used are 510 [ ( 1/ 2)]x x AL x? ?= =v , 1oB = , 2z zk Lpi= . A.2.2 Steady state solution for driven system Next, we investigate the effect of phase mixing on Alfven wave driven by a momentum source. We add a momentum source to (A.4), which becomes: ( ) ( ) sin cost y x x x y o z y o znu n u B B F k z t? ?? ? ? ? = ? +tildenosptildenosp tildenosp . (A.21) Here the driving amplitude oF is a constant, and we use 2z zk Lpi= , i.e., we try to drive a standing Alfven wave along the z direction. Other parameters and boundary conditions are the same as our simulation for the undriven problem. We start the simulations with zero initial conditions, and 510oF ? = , and let the system come to a steady state. The amplitudes of the Alfven wave, yBtildenosp , at the center of the simulation box, are measured as a function of driving frequency, and are plotted in FIG. A.5. The spectrum corresponding to flat density profile (and therefore flat Alfven speed profile) is also plotted in the same graph. 95 FIG. A.5, Power spectra measured at the center of the system for driven standing Alfven wave, with and without a profile in Alfven speed. Here ( 1/ 2)o z Ak x? = =v is the local resonant frequency at the center, and 510 [ / ( 1/ 2)]x x AL x? ?= =v , 310 [ / ( 1/ 2)]z x AL x? ?= =v and 510oF ?= are used. When the profile is flat, the height and width of the spectrum [full width at half maximum (FWHM)] are 33.1 10?? and 35 10?? respectively. While a simple harmonic oscillator model predicts the height ( )2 5~ ~ 6 10o z zF k? ?? and width 3 ~ ~ 1.3 10z z Ak? ??v , the simulation shows that in the case of flat profile the system behave as a driven harmonic oscillator. However, for the case where the profile is given by FIG. A.2, the height of the peak is only 42.5 10?? , and the spectrum has a wider spread with FWHM = 0.07. The amplification is smaller and the bandwidth is broadened, compared with the homogeneous case. From the theory of phase mixing, the effective damping rate can be roughly estimated as ( )1/32~ ' ~ 0.03( / )eff x z A x Ak L? ? v v , which gives 4~ 4 10o effF ? ?? and 96 2 ~ 2 10eff o? ? ?? . The deviation from the simulation result is mainly due to the fact that the profiles used are not linear; thus the damping rate will also depend on higher derivatives of the profile to a lesser extent. 97 B. Effects of nonlinearities on driving magnetosonic wave in uniform plasma This section provides an illustration of the effects of nonlinearities in the excitation of a magnetosonic wave by a resonant momentum source. We first show by an expansion in small driving amplitude, that the system behaves similarly to that of linearized equations when the driving frequency is far off the resonant frequency, and that nonlinear terms enter when the system is gets close to resonance, causing the expansion to break down. Next, we investigate spectra in the vicinity of the resonance by simulation, and show that the resonance is suppressed in comparison with the linear system. B.1 Equations and expansion in small amplitudes We start with MHD equations with ?o oB B z= harpoonrightnosp , i.e., in equilibrium, and we assume zB B? << harpoonrightnosp . The density n is also uniform, and the system is isothermal. The governing equations are continuity, momentum equation (in the perpendicular direction; we assume 1M = ) and Faraday?s law (in the direction of the magnetic field): ( ) 0tn nu? +? ? =harpoonrightnosp harpoonrightnosp , (B.1) 2 8o u B n u u nT F t pi ? ? ? ? ? ? ??? ? + ?? = ?? + +? ?? ??? ? ? ? harpoonrightnosp harpoonrightnosp harpoonrightnosp harpoonrightnospharpoonrightnosp harpoonrightnosp , (B.2) ( )t z z zB u B B u? ? ?? = ? ? ?harpoonrightnosp harpoonrightnosp harpoonrightnosp . (B.3) 98 First, to simplify the problem, we consider the low ? limit, i.e., ( )24 1onT Bpi << . In that a case, we only need to consider (B.2) and (B.3). We consider a 2D slab geometry, of size y zL L? . A standing magnetosonic wave is driven in the y direction, with the momentum source ? sin cosoF yF ky t?= harpoonrightnosp , oF = constant. FIG. B.1 shows a schematic view of the system. Since xu and B? harpoonrightnosp are unexcited and remain zero throughout, the y-components of (B.2) and (B.3) form a complete set with variables { yu , zB }. FIG. B.1. Schematic view of driving standing magnetosonic wave in a 2D slab. For clarity, we normalize yu , zB , y and t (and also 1?? ) in (B.2) and (B.3) to Av , oB , yL and y AL v respectively, where 4 o A B npi =v , and we assume 1y AL =v in our normalization. Then we rewrite our equations as ( )2 2 2 sin sint y y y zu u B ky t? ?? + ? + = , (B.4) ( ) 0t z y y zB u B? + ? = , (B.5) 99 where ( )2o y AF L n? = v . We see that the convection of flow and the magnetic energy terms account for the nonlinearities of the system. The y-flow yu and y zB? are zeros at 0y = and yy L= . We solve {(B.4),(B.5)} by a expansion in the driver amplitude ? : 2 1 2 2 1 2 ... 1 ... y z u u u B B B ? ? ? ? ? = + +? = + + +? . (B.6) To order ? , we have 1 1 1 1 sin sin 0; t y t y u B ky t B u ? ?? + ? =?? ? + ? =? (B.7) the corresponding solution is 1 2 2 1 2 2 sin cos cos sin . u ky t k kB ky t k ? ? ? ? ? ? ? ? = ??? ?? ? =? ?? (B.8) To order 2? , we have ( ) ( ) 2 2 2 2 1 1 2 2 1 1 2 ; t y y t y y u B u B B u u B ?? + ? = ?? +?? ? + ? = ???? (B.9) we note that the nonlinear terms in the system become the inhomogeneous term in the second (and higher) order equations. (B.9) gives a second order correction to (B.8): ( ) ( ) ( ) ( ) ( ) 2 2 2 2 32 2 2 2 2 2 2 2 3 2 22 2 3 sin 2 sin 2 8 3 1 cos 2 cos 2 cos 2 . 88 kk u ky t k k k B ky t ky kk ?? ? ? ? ? ? ? ? ?? ? +? = ? ? ??? +? = ? +? ? ??? (B.10) 100 The second order corrections, i.e., the nonlinearities, become important, when 2 1~u u , which implies | 1|k k? ? ?= ? < . (B.11) Thus, the expansion breaks down when the driving frequency is close to the frequency of any one of the harmonics in which condition (B.11) holds. The amplitude of the oscillation of yu (at any particular point other than the end points) is shown schematically in FIG. B.2. In the next subsection we will try to find the amplitude and width of the spectrum by simulation. FIG. B.2. The schematic diagram showing the breaking down of linear expansion as the driver frequency approaches the resonant frequency. B.2 Simulation of the system at resonance Next, we simulate the system {(B.4),(B.5)} in the vicinity of the resonance in order to investigate the behavior of the system in the nonlinear regime, as indicated by our schematic diagram FIG. B.2. Mathematica 6.0 is used. We add a viscous term 2y yu?? 101 and a resistive term 2y zB?? on the R.H.S. of (B.4) and (B.5) respectively, so that the system can relax to the steady state even if the initial conditions are not the steady state solutions. The dissipative coefficients ? and ? are chosen such that 2 2 ~k k? ? ?<< , thus having negligible effects on the steady state. We start from zero initial conditions, and k pi= for the driver term, i.e., we drive a half wavelength mode along the y direction. We use 310? ?= and 310? ?= throughout the simulation. The system is allowed to reach steady state and the amplitude in yu is noted. For a fixed driver amplitude ? , we plot the amplitude of yu against the driver frequency (in ? ), and obtain a frequency spectrum. FIG. B.3 shows the spectra obtained at different driver amplitudes, ? , as parameter. From linear harmonic oscillator theory, the amplification is given by 2 2 ~ 300o k k k? ? ?= , and the peak of spectra would be ~ 300? . However, FIG. B.3 shows that the amplification is significantly lower than the linear theory, and the spectra are not Lorentzian. The resonance is also broadened. 102 FIG. B.3. The spectra of driven magnetosonic wave at different driver amplitude ? . Note the flattening of the spectra at a neighborhood of the resonant frequency o k? = ( 0? = ), as a result of the nonlinearities. FIG. B.4 shows a logarithmic plot of the width (FWHM) versus driver amplitude. The slope of the graph is 0.496 0.004? , which is consistent with the asymptotic analysis (B.11). A similar logarithmic plot of the height of the peak versus driver amplitude (FIG. B.5) shows that the maximum amplitude at resonance is directly proportional to 0.506 0.001? ? . This shows that in the presence of nonlinearities, driving magnetosonic waves is less effective compared with the prediction from the linear oscillator model; increasing the driver amplitude by a factor of f can only increase the output amplitude by f . This suppression effect is independent of dissipation. 103 FIG. B.4. A logarithmic plot of spectral width versus driver amplitude ? . The slope of the straight line is 0.496 0.004? . FIG. B.5. A logarithmic plot of spectral height versus driver amplitude ? . The slope of the straight line is 0.506 0.001? . 104 C. The Duffing Equation In this Appendix, we give a review of the solution of the Duffing equation with a sinusoidal driving term, in the vicinity of resonance [45,46]. [To be more accurate, we are solving a special case of the Duffing equation, which will be defined in (C.1) below.] This work has been done in many references; the purpose of this review is to give a comparison between the driven Duffing equation and the nonlinear GAM equations of Sec. 5.3. We have already shown [see Equation (5.14)] that the homogeneous GAM system is a special form of the homogeneous Duffing equation. The Duffing Equation is given by ( )2 21 ( )ox x x x f t? ? ?+ + + =dotnospdotnosp dotnosp , (C.1) where ( )x x t= and ( )x dx t dt=dotnosp . We can regard this equation as describing a harmonic oscillator driven by the forcing term ( )f t . o? is the natural frequency, with a damping coefficient ? , and the nonlinearity is characterized by the parameter ? . If x is length and t is time, then ? and o? have dimensions of time-1, ? has dimension of length-2, and f has the dimension of acceleration. In this section, we consider only the special case 0? = , i.e., no damping, for simplicity. This will serve our purpose, which is to compare our findings with the results of the nonlinear GAM equations in Sec 5.3. Consider a driving term ( ) sinof t a t?= , such that the driving frequency ? is close to o? , the resonant frequency. We then rewrite ( )f t as [ ] ( )2( ) sin (1 ) sino o o o of t a t a t? ? ? ? ? ? ?? ?= + ? = + ?? ? with 1o? ? ?? = ? , and we introduce a slow time scale 2t? ?= . The orderings ~ 1x ? << , 3~oa ? , and 2~? ?? 105 are assumed. We will write 2tx x x??= +dotnosp in a two time scale expansion. To order ? , (C.1) becomes 2 1 1 0tt ox x?+ = , (C.2) which gives [ ]1 ( )sin ( )ox A t? ? ? ?= + , (C.3) with ~A ? . To order 2? , the system is trivial. To order 3? , we get ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 2 3 2 3 3 1 1 2 2 3 2 2 sin 2 cos sin 4 3sin sin 3 sin . tt o t o o o o o o o o o o o o o x x x x a t A t A t A t t a t ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ?+ = ? ? + + ?? ? = ? + ? +? ?? ? ? + ? +? ?? ? ? ?+ + ?? ? (C.4) This gives 0A? = , ( )2o? ? ? ? ?= ? and ( )2 3 23 4 2 0o o oA A a?? ? ?? ? ? = . (C.5) Rewriting (C.5) with the normalization ( )2* o oA a A?= , we have 2 3 * *4 3 2 1 0 4 o o a A A? ? ? ? ? ? = . (C.6) Defining 2 * 4 3 4 o o a?? ? = as a dimensionless parameter characterizes the nonlinearity. (C.6) is of similar structure to (5.23), and FIG. C.1 shows a plot of * A against ?? from (C.6), for * 0.2? = . The width of the (detuned) resonance 5/3 1/3 * 3 2? ??? = ? , and the amplification is bounded by 2/3 1/3 * 2 ? ? . As ? approaches zero, i.e., the nonlinear effect is very small, ?? is small and the amplification approaches infinity; thus the system reduces to the driven undamped linear harmonic oscillator. Increasing 106 the driving amplitude will increase the width and lower the amplification, which suppresses the resonance. Note the similarity between the driven nonlinear GAM system in Sec 5.3 and the driven Duffing equation. On the other hand, unlike the nonlinear GAM, the peak at resonance for the Duffing oscillator shifts to higher frequency (for positive nonlinear parameter * ? ) compared with the natural frequency, while the corresponding nonlinear parameter of former (? ) is always positive, the peak shifts to lower frequency. This can also be seen by noting the difference in the signs of the ?? terms between (5.23) and (C.6). FIG. C.1. A plot of amplification against driving frequency of the Duffing oscillator driven by a sinusoidal force, where * 0.2? = . Note that the peak of resonance shift to higher frequency (compared to the natural frequency), which is opposite to the nonlinear GAM system shown in FIG. 5.1. 107 Bibliography [1] R. F. Ellis, A. Case, R. Elton, J. Ghosh, H. Griem, A. Hassam, R. Lunsford, S. Messer, and C. Teodorescu, Phys. Plasmas 12, 055704 (2005); doi: 10.1063/1.1896954 [2] S. Choi, P. N. Guzdar, A. Case, R. Ellis, A. B. Hassam, R. Lunsford, C. Teodorescu, and I. Uzun-Kaymak, Phys. Plasmas 15, 042507 (2008); doi: 10.1063/1.2903053 [3] I. U. Uzun-Kaymak, P. N. Guzdar, R. Clary, R. F. Ellis, A. B. Hassam, and C. Teodorescu, Phys. Plasmas 15, 112308 (2008); doi: 10.1063/1.3028312 [4] S. Chandrasekhar, Hydrodynamic and Hydromagnetic Stability (Dover, 1970) [5] K. Hallatschek and G. R. McKee , Phys. Rev. Lett. 109, 245001 (2012), doi: 10.1103/PhysRevLett.109.245001 [6] G. R. McKee, R. J. Fonck, M. Jakubowski, K. H. Burrell, K. Hallatschek et al., Phys. Plasmas 10, 1712 (2003); doi: 10.1063/1.1559974 [7] R. Nazikian et al., Phys. Rev. Lett. 101, 185001 (2008); doi: 10.1103/PhysRevLett.101.185001 [8] J. C. Hillesheim, W. A. Peebles, T. A. Carter, L. Schmitz, and T. L. Rhodes, Phys. Plasmas 19, 022301 (2012); doi: 10.1063/1.3678210 [9] J. Heyvaerts and E. R. Priest, Astron. Astrophys. 117, 220-234 (1983) [10] H. L. Kuo, Phys. Fluids 6, 195 (1963); doi: 10.1063/1.1706719 [11] R. J. Taylor, M. L. Brown, B. D. Fried, H. Grote, J. R. Liberati, G. J. Morales, and P. Pribyl, Phys. Rev. Lett. 63, 2365?2368 (1989), doi: 10.1103/PhysRevLett.63.2365 108 [12] K. Ida and S. Hidekuma, Phys. Rev. Lett. 65, 1364?1367 (1990) , doi: 10.1103/PhysRevLett.65.1364 [13] K. H. Burell et al, Phys. Fluids B 2, 1405 (1990); doi: 10.1063/1.859564 [14] K. C. Shaing, E. C. Crume, and W. A. Houlberg, Phys. Fluids B 2, 1492 (1990); doi: 10.1063/1.859473 [15] A. B. Hassam, Comments Plasma Phys. Controlled Fusion 14, 275 (1991) [16] H. Biglari, P. H. Diamond, and P. W. Terry, Phys. Fluids B 2, 1 (1990); doi: 10.1063/1.859529 [17] A. B. Hassam, Phys. Fluids B 4, 485 (1992); doi: 10.1063/1.860245 [18] A. B. Hassam , T. M. Antonsen, Jr., J. F. Drake, P. N. Guzdar, C. S. Liu, D. R. McCarthy, and F. L. Waelbroeck, Phys. Fluid B 5, 2519 (1993) [19] R. J. Groebner, Phys. Fluids B 5, 2343 (1993); doi: 10.1063/1.860770 [20] K. H. Burrell, Phys. Plasmas 4, 1499 (1997); doi 10.1063/1.872367 [21] S. DeSouza-Machado, A. B. Hassam, and Ramin Sina, Phys. Plasmas 7, 4632 (2000); doi: 10.1063/1.1316086 [22] Yi-Min Huang and A. B. Hassam, Phys. Rev. Lett. 87, 235002 (2001) [23] N. Winsor, J. L. Johnson, and J. M. Dawson, Phys. Fluids 11, 2448 (1968); doi: 10.1063/1.1691835 [24] A. B. Hassam and J. F. Drake, Phys. Fluids B 5, 4022 (1993); doi: 10.1063/1.860622 [25] V. B. Lebedev, P. N. Yushmanov, P. H. Diamond, S. V. Novakovskii, and A. I. Smolyakov, Phys. Plasmas 3, 3023 (1996); doi: 10.1063/1.871638 109 [26] S. V. Novakovskii, C. S. Liu, R. Z. Sagdeev, and M. N. Rosenbluth, Phys. Plasmas 4, 4272 (1997); doi: 10.1063/1.872590 [27] M. N. Rosenbluth and F. L. Hinton, Phys. Rev. Lett. 80, 724?727 (1998); doi: 10.1103/PhysRevLett.80.724 [28] K. Hallatschek and D. Biskamp, Phys. Rev. Lett. 86, 1223?1226 (2001); doi: 10.1103/PhysRevLett.86.1223 [29] R. G. Kleva and A. B. Hassam, Phys. Plasmas 20, 032508 (2013); doi: 10.1063/1.4794837 [30] P. W. Terry, Rev. Mod. Phys. 72, 109 (2000); doi: 10.1103/RevModPhys.72.109 [31] Y. Z. Zhang and S. M. Mahajan, Phys. Fluids B 4, 1385 (1992). doi: 10.1063/1.860095 [32] J. P. Freidberg, Ideal Magnetohydrodynamics (Plenum, New York, 1987) [33] Yi-Min Huang and A. B. Hassam, Phys. Plasmas 11, 2459 (2004); doi: 10.1063/1.1651102 [34] Yi-Min Huang and A. B. Hassam, Phys. Plasmas 11, 3738 (2004); doi: 10.1063/1.1761602 [35] C. Teodorescu, R. Clary, R. F. Ellis, A. B. Hassam, R. Lunsford, I. Uzun- Kaymak, and W. C. Young, Phys. Plasmas 15, 042504 (2008); doi: 10.1063/1.2904903 [36] C. Teodorescu, R. F. Ellis, A. Case, C. Cothran, A. Hassam, R. Lunsford, and S. Messer, Phys. Plasmas 12, 062106 (2005); doi: 10.1063/1.1924391 [37] M. D. Kruskal and R. M. Kulsrud, Phys. Fluids 1, 265 (1958); doi: 10.1063/1.1705884 110 [38] R. J. Goldston and P. H. Rutherford, Introduction to Plasma Physics (Taylor and Francis 2000) [39] P. N. Guzdar, J. F. Drake, D. McCarthy, A. B. Hassam, and C. S. Liu, Phys. Fluids B 5, 3712 (1993) [40] A. B. Hassam, W. Hall, J. D. Huba, and M. J. Keskinen, J. Geophys. Res., 91, 13513 (1986); doi:10.1029/JA091iA12p13513 [41] P. H. Diamond, S.-I. Itoh, K. Itoh and T. S. Hahm, Plasma Phys. Control. Fusion 47 R35 (2005); doi: 10.1088/0741-3335/47/5/R01 [42] H. Sugama and T.-H. Watanabe, J. Plasma Phys. 72, 825 (2006); doi: 10.1017/S0022377806004958 [43] Z. Qiu, L. Chen and F. Zonca, Plasma Phys. Control. Fusion 51, 012001 (2009); doi: 10.1088/0741-3335/51/1/012001 [44] A. Hasegawa and L. Chen, Phys. Rev. Lett. 32, 454?456 (1974); doi: 10.1103/PhysRevLett.32.454 [45] C. M. Bender and S. A. Orszag, Advanced Mathematical Methods for Scientists and Engineers (Springer, 1999) [46] M. H. Holmes, Introduction to Perturbation Methods (Springer, 1998)