Princeton Plasma Physics Laboratory PPPL-4075 PPPL-4075 Numerical Study of Field-reversed Confi gurations: The Formation and Ion Spin-up E.V. Belova, R.C. Davidson, H. Ji, M. Yamada, C.D. Cothran, M.R. Brown, and M.J. Schaffer June 2005 PPPL PRINCETON PLASMA PHYSICS LABORATORY Prepared for the U.S. Department of Energy under Contract DE-AC02-76CH03073. PPPL Report Disclaimers Full Legal Disclaimer This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, nor any of their contractors, subcontractors or their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or any third party’s use or the results of such use of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise, does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof or its contractors or subcontractors. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof. Trademark Disclaimer Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise, does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof or its contractors or subcontractors. PPPL Report Availability This report is posted on the U.S. Department of Energy’s Princeton Plasma Physics Laboratory Publications and Reports web site in Fiscal Year 2005. The home page for PPPL Reports and Publications is: http://www.pppl.gov/pub_report/ Office of Scientific and Technical Information (OSTI): Available electronically at: http://www.osti.gov/bridge. Available for a processing fee to U.S. Department of Energy and its contractors, in paper from: U.S. Department of Energy Office of Scientific and Technical Information P.O. Box 62 Oak Ridge, TN 37831-0062 Telephone: (865) 576-8401 Fax: (865) 576-5728 E-mail: [email protected] National Technical Information Service (NTIS): This report is available for sale to the general public from: U.S. Department of Commerce National Technical Information Service 5285 Port Royal Road Springfield, VA 22161 Telephone: (800) 553-6847 Fax: (703) 605-6900 Email: [email protected] Online ordering: http://www.ntis.gov/ordering.htm Numerical Study of Field-Reversed Configurations: The Formation and Ion Spin-up E. V.Belova 1),R. C. Davidson1),H.Ji 1),M. Yamada1),C. D. Cothran2), M.R. Brown2), M. J. Schaffer3) 1)PrincetonPlasmaPhysics Laboratory, PrincetonNJ,USA,2)SwarthmoreCollege,Swarth- more PA, USA, 3)General Atomics, SanDiegoCA, USA E-mail: [email protected] Shorttitle: NumericalStudyofFRCs PACS numbers: 52.55.Lf, 52.65.Kj, 52.65.Rr, 52.35.Vd 1 Abstract Resultsofthree-dimensionalnumericalsimulationsoffield-reversedconfigurations(FRCs)are presented. Emphasis of this work is on the nonlinear evolution of magnetohydrodynamic(MHD) instabilities in kinetic FRCs, and the new FRC formation method by counter-helicity spheromak merging. Kinetic simulations show nonlinear saturation of the n = 1 tilt mode, where n is the toroidal mode number. The n = 2 and n = 3 rotational modes are observed to grow during the nonlinear phase of the tilt instability due to the ion spin-up in the toroidal direction. The ion toroidalspin-up is shown to be related to the resistive decay of the internal flux, and the resulting loss of particleconfinement. Three-dimensional MHD simulations ofcounter-helicityspheromak merging and FRC formation show good qualitative agreement with results from the SSX-FRC experiment. The simulationsshowformationofanFRC inabout20−30 Alfve´ntimesfortypical experimental parameters. The growth rate of the n = 1 tilt mode is shown to be significantly reduced compared to the MHD growth rate due to the large plasma viscosity and field-line-tying effects. 2 1. Introduction The field-reversed configuration (FRC) is a compact toroid with little or no toroidal field. It offers a unique fusion reactor potential because of its compact and simple geometry, translation properties, and high plasma beta [1]. At present, the most importantissues are FRC stabilitywith respecttolow-n(toroidalmodenumber)MHDmodes,andthedevelopmentofnewFRCformation andcurrent-drivemethods. The traditionaltheta-pinchformationmethodusually produceshighlykineticFRCs with rel- ativelylowflux, smallS(cid:3) (theFRCkinetic parameterS(cid:3) isthe ratiooftheseparatrixradius tothe ionskin depth)and largeelongationE. Atheoreticalunderstandingoftheobserved FRCstability properties has proven to be elusive, due to the complicated interplay of several non-ideal MHD effects[2–4],includingfinite-Larmor-radius(FLR)effects,theHallterm,andplasmafloweffects. Advancednumericalsimulationsare requiredtodescribe theself-consistent stabilitypropertiesof kineticFRCs [3,5,6]. Results ofaset ofsuch simulationsarepresented inthispaper. A slow FRC formationtechnique, based on counter-helicityspheromak merginghas demon- strated the advantage of this approach compared with traditional theta-pinch formation methods [7]. The counter-helicity spheromak merging method allows formation of the configuration with large S(cid:3), thus permitting experimental studies of large-S(cid:3) FRC stability properties. The SSX- FRC experiment[8]is designed tostudy FRC formationbythe counter-helicityspheromakmerg- ing method, and to examine the general issue of FRC stability properties at large S(cid:3). Three- dimensional MHD simulations have been performed in support of the SSX-FRC experiment, and showgoodagreementwiththe experimentalresults. 2. StudyofFRC nonlinearstabilityproperties Numerical studies of the nonlinear evolution of magnetohydrodynamic (MHD) instabilities 3 in kinetic, prolate FRCs (theta-pinch-formedFRCs) have been performedusing the 3D nonlinear hybridandMHDsimulationcodeHYM [5]. The stabilitypropertiesofMHDmodeswithtoroidal mode numbers n (cid:21) 1 have been investigated, including finite-ion-Larmorradius (FLR) and rota- tional effects. It has been demonstrated that due to the strong FLR stabilization of the higher-n modes, the n = 1 tilt mode is the most unstable mode for nearly all experimentally-relevantnon- rotatingFRC equilibria[3]. AnempiricalFLR scaling ofthe tiltmodelineargrowthratehas been obtained: γ = CV =R Eexp(−3E(cid:26) =R ),whereγ = CV =R E istheMHDgrowthrate,R A s i s mhd A s s is theseparatrixradius,E is theseparatrixelongation,(cid:26) is ionthermalLarmorradius,and C (cid:25) 2 i is aconstant. NonlinearkineticsimulationsperformedforasetofFRCequilibriawithE = 4−6andS(cid:3) = 10−80 show that the n = 1 tiltmode saturates nonlinearlywithoutdestroying the configuration, provided the FRC kinetic parameter is sufficiently small, S(cid:3) (cid:24) < 20. In addition to the saturation of the tilt mode, nonlinear hybrid simulations show that the ions spin-up toroidally in the ion diamagneticdirection, and the n = 2 modegrows inthe nonlinearphase ofthe simulation. Initial conditionsforthe simulations have been set at t = 0 so that allof the equilibriumtoroidalcurrent is carriedby the electrons, andthe ions have a non-rotatingMaxwellian distribution. These initial conditions are consistent with experimental observations just after FRC formation. However, as the simulation proceeds, the ions graduallybegin to rotate, and near the end of the simulationrun theiontoroidalflowvelocitybecomes comparableto theiondiamagneticvelocity. Therefore, the saturation of the tilt instability occurs in the presence of a significant ion toroidal rotation, and it is accompanied by the growth of the n = 2 rotational mode, which is often seen in experiments [1]. The saturation of the n = 1 tilt mode and the growth of the n = 2 rotationalmodecanbeseen inFig.4ofRef.[3],wheretheresultsofnonlinearhybridsimulations 4 are shown for a configuration with E = 6:25 and S(cid:3) (cid:25) 20 (for comparison, the parameters for typicaltheta-pinchFRC experiments [1]areE = 5−8 and S(cid:3) (cid:24) < 20). Inthe nonlinearphase, the ionrotationrateiscomparabletothelineargrowthrateofthetiltmode. Therefore,theiontoroidal spin-up, in addition to driving the rotational instability, is likely to contribute significantly to the saturationofthetiltinstability. 2.1. Ionspin-up A set of 2D (axisymmetric) nonlinear hybrid simulations has been performed in order to study theresistive evolutionofthe kineticFRC, andinvestigate the mechanismof theion toroidal spin-up. The initial ion distribution function is taken to be f = f(") = Aexp(−"=T ), where 0 " = m v2=2+e(cid:30) is the ion energy, T is the (uniform)ion temperature, and (cid:30) is the electrostatic i 0 potential. Since the ion pressure is scalar, a kinetic equilibrium for the hybrid runs has been constructedsimilartotheidealMHDequilibrium,i.e., bysolvingtheGrad-Shafranovequationfor a chosen pressureprofile. The cold-fluiddescriptionis used forthe electrons, andquasi-neutrality is assumed. The numerical simulations have been performed for an FRC with E = 4, S(cid:3) = 20, and an elliptical separatrix shape. The value of the normalized resistivity at the field null is (cid:17) = o 1=S = 10−4, where S = V R =(cid:17) is the Lindquist number, and a resistivity profile has been used A c with(cid:17) inversely proportionalto theplasmadensity. Thesimulationsshowthatthereisasignificantparticlelossassociatedwiththeresistivedecay of the poloidal flux. Figure1 shows the time evolutionof the trapped poloidalflux (t), and the 0 number of ions inside the separatrix region < 0. It can be seen that the decay of the poloidal flux results in the loss of a significant fraction of the particles at t (cid:24) > 50t . Here, the Alfve´n A time is defined as t = R =V , where R is the radius of the flux conserving shell, and V is the A c A c A characteristic Alfve´n velocity. Analysis of the particle phase-space shows that the particles with 5 initial values of p such that 0 < p (cid:24) < (cid:1) are lost from the closed-field-line region when the (cid:30) (cid:30) absolute value of the trapped poloidal flux is reduced by (cid:1) . Here p = m =eRv − is the 0 (cid:30) i (cid:30) canonical toroidal angular momentum of the ion, and the signs are chosen such that p > 0 is (cid:30) the approximate confinement condition. Since most of the particles with small p have negative (cid:30) toroidal velocity, the particle loss results in a net flux of the negative momentum away from the separatrix region, and therefore there is a net positive ion rotation inside the separatrix (where the positive direction corresponds to the current direction). Therefore, the ion toroidal spin-up in the simulations is related to the resistive decay of the internal flux, and the resulting loss of weakly-confinedparticles. Figures 2a shows the time evolution of the toroidal angular momentum of all ions, L = R R R v f d3vd3x, and the angularmomentumofthe ions inside the separatrixregion obtainedin (cid:30) the simulations shown in Fig. 1. It can be seen that the net angular momentum is approximately conservedduetotheimposedperiodicboundaryconditionsinz-direction. Incontrast,theangular momentum of the part of the plasma confined inside the separatrix is positive, and increases in time as the configuration decays. (Due to the imposed periodic boundary condition, particles which leave the closed-field-line region remain on the open field lines, so that the plasma outside the separatrix spins-up in the negative direction.) The maximum value of the ion flow velocity is plotted in Fig. 2b. The peak value of V is about 0.2-0.3 V , which is comparable to the ion (cid:30) A diamagnetic velocity for S(cid:3) (cid:25) 20. In the final state, the ions carry a significant fraction of the total current. Similarvalues ofthe ion toroidalflowvelocity have also been obtained in nonlinear three-dimensional(3D)simulations[3]. Both 2D and 3D simulations with zero initial ion rotation demonstrate the formation of an approximaterigid-rotorprofileinsidetheseparatrixinabout40-60Alfve´ntimes,dependingonthe 6 plasma resistivity. Poloidal contour plots and radial profiles of the ion toroidal flow velocity are showninFig.3. Figure3shows thattheionflowvelocitychanges signoutside theseparatrix,and that there is a significant velocity gradient close to the separatrix at the plasma edge. The n = 2 and n = 3 rotational instabilities are found to be localized in the vicinity of maximum velocity shear near the edge, and have a similar structure to the external modes. The growth rates of the rotational modes are found to be larger in smaller-S(cid:3), more kinetic configurations. The details of theiontoroidalspin-updeterminethe nonlinearevolutionoftheseinstabilities. Toroidal ion spin-up has always been observed during the quasi-steady-state decay phase of FRC experiments [1]. The measured ion rotation frequency is comparable to the ion diamagnetic frequency (cid:11) = Ω =Ω (cid:24) 1 by the time when about one-half of the plasma has been lost due to i di decay, and the n = 2 rotational instability begins to grow. In addition, recent observations in the TS-3andTS-4experiments[9]suggestthattheiontoroidalspin-upplaysanimportantroleduring the observed nonlinear stabilization of the n = 1 tilt mode. These observations are in a good agreementwiththeresults of2Dand3D hybridsimulationsdescribed aboveand inRef. [3]. Asignificant numberoftheoreticalandexperimentalstudies have beenperformedtoinvesti- gate the plasma spin-upin FRCs [1]. Two majorphysical mechanisms which have been proposed to explain the observed ion rotation are: particle loss [10], and the end-shortening of the radial electric field [1,11]. The simulation results presented here provide good agreement with the ob- served ion rotation,even though theend-shorteningmechanism is not includedin oursimulations due to the imposed periodic boundary conditions. The experimental estimates of the particle loss time and the spin-up time are also consistent with the assumption of spin-up resulting from the particle loss, at least for smaller FRC devices [12]. The conclusion of Reference [11] that the particle loss will result in spin-up in the wrong direction when electric-field efects are taken into 7 account is incorrect, because (1) it neglects the magnetization current of the lost ions; and (2) it considersrotationofthelostionsrelativetothelaboratoryframe,ratherthantheirrotationrelative tothe bulkionpopulation. The linear velocity profile (V (cid:24) R) obtained in 2D simulations (Fig. 3) suggests that the (cid:30) distributionfunctionf ofthe ionsinside theseparatrixevolves towardsanexponentialrigid-rotor i distributionfunction, whichcorresponds toa shiftedlocal Maxwelliandistribution. The evolution off towardsaMaxwelliandistributioninacollisionless modelandintheabsence of3D instabil- i ities mayindicatethat thestochasticity oftheionorbitsplays a significantroleinFRC relaxation. Earlierone-dimensional(neglectingaxialvariation)hybridsimulationsofthe iontoroidalspin-up have modeled the particle loss by removing the simulation particles whose gyrocenter was close to the separatrix [13]. The resulting plasma rotationwas localized near the separatrix, in contrast to the nearly-rigid-rotation profiles found in the 2D simulations described above (Fig. 3). This differencein the rotation profiles can be explained by the complexity of the ion orbits in the two- dimensionalpoloidalfield,and (perhaps)themodelchosen fortheparticleloss inRef. [13]. 3. Counter-helicityspheromak mergingsimulations FRC formationby the counter-helicitysperomak merging method has been developed in the TS-3 experiments in Japan [7,9]. These experiments have shown that an FRC is formed after the opposing toroidal magnetic fields of two merging spheromaks are annihilated, and the plasma is heated to form an FRC-like pressure profile. This method allows the slow formation of an FRC withlargepoloidalflux, as wellas studies ofmagneticreconnectionand MHD relaxationin high- betaplasmas [14,15]. The SSX-FRC experiment is designed to study the new FRC formation method by counter- helicity spheromak merging, and the FRC stability properties for large values of S(cid:3). In addition, 8 theeffectsoftheresidual(axially-antisymmetric)toroidalfieldonmacroscopicstabilityproperties are being studied. Experimental results demonstrate the formation of an FRC-like configuration with S(cid:3) > 35, and indicate the presence of the global n = 1 instability, consistent with the tilt- mode instability [8]. However, the observed growth rate of this instability is smaller by a factor of 6-8 than that of the ideal MHD growth rate. In addition, the experimental formation studies consistently show the presence of an axially-antisymmetric toroidal field, which does not com- pletelyannihilateduringthe reconnection. The underlyingreasons forthisare notyetunderstood. NumericalsimulationsusingtheHYM codehavebeen performedtoinvestigatethisissues, andto study3D spheromakmergingforexperimentally-relevantparameters. Two- and three-dimensional MHD simulations of counter-helicity spheromak merging have been performed for the SSX-FRC geometry and parameter range. The simulations have been performedusingthe3DnonlinearresistiveMHDversionoftheHYMcode[5]withhighresolution up to 127 (cid:2) 513 grid points in the poloidal (R;Z) plane. The boundary conditions are taken to correspond to a cylindrical flux conserver which is perfectly conducting on the perturbation (fast)time-scale,butallowingforequilibriumpoloidalmagneticfieldpenetration,thustakinginto account the field-line-tying effects present in the experiment. The initial spheromak formation by plasma guns have not been simulated. Instead, the initial conditions for the numerical studies have been chosen tocorrespondto the experimentalconditions at the beginningofthe spheromak mergingprocess. Thustheinitialconditionsinthesimulationscorrespondtotwo,low-beta,nearly force-freespheromakswithoppositetoroidalmagneticfields, placedclosetothedevicemidplane. The SSX-FRC experimental parameters are as follows: flux conserver radius and length are R = 20 cm and L = 60 cm, respectively, the edge magnetic field is B = 1 kG, the plasma c 0 densityis1015cm−3,theplasmatemperature(aftertheFRCformation)isT (cid:25) 30eV,andtheFRC 9 kinetic parameter is S(cid:3) = R =(cid:21) = 28. The Alfve´n time is defined as t = R =V , which can c i A c A be estimated as t = 2:8(cid:22)s. It has been observed that the poloidal flux at the midplane rises to A its maximum value in t (cid:24) 10 − 15t (after the spheromaks are ejected at t = 20(cid:22)s), and that a A high-betaFRC-likeconfigurationisformedatthemidplaneatt (cid:24) > 20−25t . Growthofthen = 1 A instability is seen at t (cid:24) > 18t , when the perturbation amplitude becomes larger than the level of A the n = 1 turbulence, i.e., at about 10% of poloidal magnetic energy density. The characteristic growth time of the n = 1 mode is estimated to be 10 − 14t , and the n = 1 component of the A magneticfieldenergybecomesthedominantcomponentatt (cid:24) > 30t . Theexperimentalparameters A andresults aredescribed ingreaterdetail elsewhere[8]. 3.1. Axisymmetric simulations A set of axisymmetric simulations of counter-helicity spheromak merging have been per- formedinordertostudythedependenceofthereconnectionrateandthetoroidalfieldannihilation onvaluesoftheplasmaresistivityandviscosity. Forsimplicity,theresistivityandviscosityprofiles have been assumed to be uniform in the simulations. The Lindquist number S = V R =(cid:17) based A c on classical resistivity can be estimated as S (cid:25) 103 for T (cid:24) T (cid:25) 12eV. For these temperatures, e i the ions are collisional with ! (cid:28) (cid:24) 1, and there are several estimates for the normalized plasma ci i viscosity coefficient (cid:23) = 1=Re, where Re is the Reynolds number. The estimates of Braginskii’s unmagnetized and weakly-magnetized ion viscosity are (cid:23) (cid:25) 8 (cid:2) 10−3 and (cid:23) (cid:25) 2 (cid:2) 10−3, re- 0 1 spectively, and the gyroviscosity is (cid:23) (cid:25) 5 (cid:2) 10−3. Due to the strong dependence on the ion gyro temperature, there is a large uncertainty in the above values. Nevertheless, it can be seen that the SSX-FRCplasmais intheviscosity-dominatedregimewith(cid:23) > (cid:17). Numericalsimulations for(cid:17) = 10−3 and(cid:22) = 10−3 −4(cid:2)10−3 showformationofanFRC in about20-30Alfve´ntimes,where(cid:22) = n^(cid:23) isthedynamicviscosity,andn^ isthenormalizedplasma 10 density. Evidently, large toroidal and poloidal flows with flow velocity up to (cid:24) 0:5−1V (based A on the edge field) are generated during the reconnection phase (Fig. 4b). The plasma pressure is significantlyincreased byOhmicandviscous heating,andtheFRC-likepressureprofileisformed due to convective transport by poloidal flows (Fig. 4a). The flow velocity reduces by an order- of-magnitude after the FRC formation is complete. These results agree with previous numerical studies ofaxisymmetricspheromakmerging[16,17]. Figure 5 shows the results of five MHD simulation runs performed for values of (cid:17) = 2:5 (cid:2) 10−4 − 4 (cid:2) 10−3, and (cid:22) = 10−3. It can be seen that the time for a complete reconnec- tion depends weakly on the resistivity (this time can be approximately determined by the strong reduction in the kinetic energy in Fig. 5a). Thus the reconnection is complete and the FRC forms inabout20t forlargervaluesof(cid:17),andinabout25-30t forsmaller(cid:17). Thisisconsistentwithpre- A A vious theoretical studies of driven magnetic reconnection [16,18]. On the other hand, the details of the time evolution of the plasma kinetic energy (Fig. 5a) and magnetic energy (Fig. 5b) vary strongly with (cid:17). For larger values of (cid:17), the kinetic energy peaks at t (cid:24) 5t , and the flow velocity A (mostlytoroidal)increaseswith(cid:17). Thisinitialincreaseinthekineticenergyiscaused bytheforce imbalance present in the initial conditions. For values (cid:17) (cid:21) 0:004, the configuration decays faster than the FRC forms. In the smaller resistivity cases ((cid:17) (cid:20) 0:002), there is a maximum in the flow energy which occurs at t (cid:24) 20 −25t . This peak in the kinetic energy correlates with the large A reconnection rate and the fast reductionin the magnetic energy seen in Fig. 5b. The toroidalflow velocityinthese cases iscomparabletothe characteristicAlfve´nvelocity,V (cid:24) V . (cid:30) A Figure 6 shows the results of a set of MHD simulation runs performed to investigate effects of plasma viscosity for (cid:17) = 5(cid:2) 10−4 . As expected, larger viscosity results in a strong reduction of the plasma flow (Fig. 6a). The reduced flow velocity, on the other hand, is related to a slower 11 reconnectionofthetoroidalfield(Fig.6b). Figure7shows contourplots ofthe toroidalfieldfrom the simulations with (cid:17) = 0:001, and different values of viscosity (cid:22) = 0:001, and (cid:22) = 0:004. It canbeseenthatlargervaluesofplasmaviscosityresultinincompletereconnectionandsignificant residualtoroidalfields,whicharepresentatt (cid:24) 20−30t . Figure7alsoshowsthereversalofthe A initialtoroidalfieldontheouterfluxsurfacesat t (cid:24) > 20t associated withtheso-called“sling-shot A effect” [14,15]. For the case with (cid:17) = (cid:22) = 0:001 (Fig. 7a), the magnitude of the reversed field is about 10% of the maximum toroidal field at t = 25t . This effect is weak in the large viscosity A case(Fig.7b)duetothereductioninthetoroidalflowvelocity. Thewidthofthereconnectionlayer estimatedfromtheaxialprofileofJ is(cid:14) (cid:25)0.8cmforthelow-viscositycasewith(cid:22) = 0:001,butit R increasesto(cid:14) (cid:25)1.6cmfor(cid:22) = 0:004. Forcomparison,theexperimentallymeasuredreconnection layeris about2cm to3cmwide [19]. Therefore, the numerical results indicate that the relatively large values of plasma viscos- ity in the experiments ((cid:22) (cid:24) > 0:002) may be responsible for the incomplete reconnection of the toroidalmagnetic field which is observed in SSX-FRC. The simulations predict a reduced “sling- shot” effect in the SSX-FRC experiment, and a reconnection layer width of about 2 cm, which is comparablewiththeexperimentallymeasuredone. 3.2. Three-dimensionalsimulations Three-dimensional MHD simulations have been performed to study the effects of plasma viscosity, self-generatedflows and magnetic field-line-tyingeffectson the unstable globalmodes. A random initial perturbationis applied at t = 0. The n = 1 tilt mode is found to be a dominant mode in all cases, and the configuration is tilted at the end of the simulations. Figure 8a shows the time evolution of the n = 1 mode obtained in four simulation runs with (cid:17) = 10−3 and (cid:22) = 5(cid:2)10−4 −4(cid:2)10−3. In the linear phase, t (cid:24) < 20t , the growth rate is largest forsmall viscosity, A 12 and the growth rate reduces as (cid:22) increases. The simulations demonstrate that large values of viscosity have a stabilizing effect on the n = 1 tilt mode, as well as the higher-n MHD modes. Thenonlinearslow-downofthetiltinstability,seen att (cid:25) 20−25t inlow-viscosityruns,occurs A when the amplitude of the n = 1 mode becomes comparable to that of the n = 0 mode. The slowing-downofthe tiltmotionin othercases ((cid:22) = 0:002 and (cid:22) = 0:004) observed att > 20t is A due to the evolutionof the backgroundquasi-equilibrium, and the resistive decay of the magnetic field,which reducesthe instabilitydrive. Figure 8b shows the time evolution of the total kinetic energy obtained for the same set of simulations as shown in Fig. 8a. The growth of the n = 1 tilt mode is not seen until t > 20t , A when the amplitude of this mode becomes comparable to the n = 0 component of the kinetic energy. Note thatbythattimethegrowthrate ofthe n = 1 modeis reducedcomparedtoitslinear value. The calculated growth rates obtained from the simulation shown in Fig. 8a for (cid:22) = 0:004 are: γ (cid:25) 0:33 (cid:1) V =R for t < 20t , and γ (cid:24) < 0:15 (cid:1) V =R for t > 20t . For comparison, an A s A A s A estimatefortheidealMHD growthrateisγ (cid:24) (0:7−1:3)V =R . Severalfactorsmaycontribute 0 A s to the reduction of the instability growth at later times (t (cid:24) > 20t ), including nonlinear mode A interactionsandan increaseofthe separatrixelongation. Additionaltime-evolutionplotsareshowninFig.9wherethepeakvaluesofthepoloidaland toroidalfields, the radial currentdensity, the axial flowvelocity, and the toroidalflowvelocity are shown for 3D simulations with (cid:17) = 0:001 and (a) (cid:22) = 0:001 and (b) (cid:22) = 0:004. The growth of the tilt instability can be seen in these plots at t > 25t for (cid:22) = 0:001, and at t > 35t for A A the (cid:22) = 0:004 simulation. The toroidal magnetic field reconnection rate is proportional to the radial current J , which is significantly smaller in case (b) for t (cid:24) < 22t . Figure9 shows that the R A maximum toroidalflow velocity is reduced approximatelyby a factorof twowhen (cid:22) is increased 13 from0.001 to0.004. The possible stabilizing effect of the toroidal flows, generated during the reconnection pro- cess, appearsto beless significantthanthe stabilizingeffectofviscosity, because thegrowthrates for the larger V cases (i.e., smaller (cid:22) ) are larger than the growth rates for smaller V (larger (cid:22)). (cid:30) (cid:30) The finite residual toroidalfield, on the other hand, may contributeto the reduction of the growth rateofthen = 1 modeinthe high-viscositycases ((cid:22) > 10−3). ThesimulationsshowninFigs.8and9havebeenperformedforrealisticboundaryconditions, includingtheeffectsofmagneticfieldline-tying. Anotherset ofsimulationshave beenperformed with different boundary conditions, i.e., by neglecting these effects. A comparison between the runswithdifferentboundaryconditionsshowthatline-tyingeffectsincreasethereconnectiontime by5−10t ,and reducethen = 1modegrowthrate. Ithas alsobeen foundthatthepeak toroidal A flowvelocityis reducedbya factoroftwodue toline-tyingeffects. Insummary,thesimulationsshowthatboththelargeplasmaviscosityandthefield-line-tying boundary conditions reduce the growth rate of the tilt mode. In addition, there is a nonlinear reduction of the instability drive at the time when the mode amplitude becomes experimentally observable. These results may provide an explanation for the experimentally-measured growth ratesthat areabout6−8 times smallerthantheideal MHDgrowthrate. 4. Conclusions The hybrid simulations presented here show that, while ion FLR effects determinethe linear stability properties of non-rotatingFRCs, the inclusion of nonlinear and ion-toroidal-floweffects is necessary for a satisfactory description of plasma behavior in low-S(cid:3) FRC experiments. In particular,ithasbeenshownthattheiontoroidalspin-upplaysanimportantroleinFRCnonlinear evolution,includingthatofthen = 1 tiltmode. 14 The 3D hybrid simulations have been able to reproduce all major experimentally observed stabilitypropertiesofkinetic(theta-pinch-formed)FRCs. Namely,thescalingofthelineargrowth rateofthen = 1tiltinstabilitywithS(cid:3)=E parameterhasbeendemonstratedforaclassofelongated ellipticFRCs[3];andiontoroidalspin-up,thenonlinearsaturationofthetiltmode,andthegrowth ofthen = 2rotationalmodehavebeen demonstrated. Ithasbeen shownthattheloss ofionswith a preferentialsign of toroidal velocity due to the resistive decay of the poloidal flux results in the ion toroidal spin-up, which reproduces very well the experimentally-observed ion rotation. The timescale oftheionspin-upis determinedbythefluxdecay time. The MHD version of the HYM code has also been used to study FRC formation by the counter-helicityspheromakmerginginsupportoftheSSX-FRCexperiment[8],andcontributedto interpretationofseveralpuzzlingexperimentalobservations. Aparameterscan,whichisnoteasily accessible experimentally, has been carried out using numerical simulations, and it has provided a better qualitative understanding of the several of the experimental results. In particular, the persistence of the residual toroidal field, and the slower-than-MHD growth of the tilt instability have beenshowntoberelatedtothe largeplasmaviscosityand line-tyingeffectsin theSSX-FRC experiments. 15 Acknowledgments This research was supported by DOE contract DE-AC02-76CH03073. Calculations were performedatthe U.S. NationalEnergyResearch SupercomputingCenter. 16 References [1] TuszewskiM.1988Nucl.Fusion282033 [2] IshidaA.,MomotaH.andSteinhauerL.C.1988Phys.Fluids313024 [3] BelovaE.V.,DavidsonR.C.,JiH.andYamadaM.2004Phys.Plasmas112523 [4] BarnesD.C. 2002Phys.Plasmas9560 [5] BelovaE.V.etal2000Phys.Plasmas74996 [6] OhtaniH.,HoriuchiR.andSatoT.2003Phys.Plasmas10145 [7] OnoY.etal1997Phys.Plasmas41953 [8] CothranC.D.etal2003Phys.Plasmas101748 [9] OnoY.etal2003Nucl.Fusion43649 [10] EberhagenA.andGrossmannW.1971Z.Phys.248139 [11] SteinhauerL.C.2002Phys.Plasmas93851 [12] TuszewskiM.etal1982Phys.Fluids251696;1988Phys.Fluids31946 [13] HarnedD.S.andHewettD.W.1984Nucl.Fusion24201 [14] Yamada M.etal1990Phys.Rev.Lett.65721 [15] OnoY.etal1996Phys.Rev.Lett.763328 [16] Watanabe T.-H., Sato T. and Hayashi T. 1997 Phys. Plasmas 4 1297; Sato T., Oda Y. and Otsuka S. 17 1983Phys.Fluids263602 [17] LukinV.S.etal2001Phys.Plasmas81600 [18] BiskampD.andWelterH.1980Phys.Rev.Lett.441069 [19] KornackT.W.,SollinsP.K.andBrownM.R.1998Phys.Rev.E58R36 18 Figure captions Fig.1. Time evolution of (a) the normalized value of trapped poloidal flux, and (b) the nor- malized number of ions inside the separatrix obtained from 2D hybrid simulations with S(cid:3) = 20 andE = 4. Fig.2. (a)Timeevolutionofthenormalizedangularmomentumofallions,andtheionsinside the separatrix; and (b) the maximum value of the ion toroidal flow velocity obtained from same simulationsas inFig. 1. Fig.3. (a) Contour plots of the ion toroidal velocity in the r − z plane at t = 40t and A t = 80t obtainedfrom2DhybridsimulationswithS(cid:3) = 20andE = 4;and(b)Radialprofilesof A the ion toroidal flow velocity at the FRC midplane at t = 20, 40, and 80t . The separatrix radius A is R =R (cid:25) 0:6. s c Fig.4. Contourplotsof(a)the pressure, and (b)the toroidalvelocity at t=t =10, 14, 20, 25 A obtained from 2D MHD simulations of counter-helicity spheromak merging for (cid:22) = (cid:17) = 0:001. The velocity is normalized to the characteristic Alfve´n velocity; the contourvalues correspond to V = 0:4 andV = −0:2fort = 10−20, and V = 0:2and V = −0:1fort = 25. max min max min Fig.5. Time evolution of (a) the normalized kinetic energy and (b) the magnetic energy ob- tained from 2D MHD simulations of counter-helicityspheromak merging for (cid:22) = 0:001 and sev- eralvalues ofresistivity. Fig.6. Time evolution of (a) the normalized kinetic energy and (b) the magnetic energy (to- tal and toroidal) obtained from 2D MHD simulations for (cid:17) = 5 (cid:2) 10−4 and different values of viscosity: (cid:22) = 0:0005 (dottedline),(cid:22) = 0:001 (solid),and(cid:22) = 0:002 (dashed). Fig.7. Contour plots of the toroidal magnetic field at the poloidal plane obtained from 2D MHD simulations with (cid:17) = 0:001 and (a) (cid:22) = 0:001 (t=t = 10;20;30); and (b) (cid:22) = 0:004 A 19 (t=t = 20;30;40). The maximum toroidal magnetic field value (normalized to the initial edge A field)isshown foreach plot. Fig.8. Timeevolutionof(a)then = 1 modeenergy, and(b)thetotalkineticenergyobtained from four 3D MHD simulation runs of counter-helicity spheromak merging for (cid:22) = 5 (cid:2) 10−4 − 4(cid:2)10−3. Fig.9. Time evolution of the maximum values of the poloidal and toroidal magnetic fields, the radial current density, the axial flow velocity, and the toroidal flow velocity obtained in 3D simulationswith(cid:17) = 0:001 and (a)(cid:22) = 0:001, and (b)(cid:22) = 0:004. 20 40 35 (a) 30 y 0 25 20 15 1.0 (b) N (y <0) 0.75 0.5 0 10 20 30 40 50 60 70 80 t / t A Figure1. 21 4 (a) 3 L in 2 L 1 0 L total -1 0 10 20 30 40 50 60 70 80 0.3 (b) 0.2 V / V j A 0.1 0 0 10 20 30 40 50 60 70 80 t / t A Figure2. 22 (a) t= 40 t= 80 0.3 (b) 0.2 t= 80 t= 40 V / V 0.1 j A t= 20 0.0 - 0.1 0.0 1.0 R / Rc Figure3. 23 (a) (b) t= 10 t= 14 t= 20 t= 25 Figure4. 24 0.004 (a) 0.012 0.01 0.008 0.0005 E 0.002 K 0.006 0.00025 0.004 0.001 0.002 0 0 10 20 30 40 (b) 0.6 0.4 E M 0.00025 0.2 0.0005 0.001 0.004 0.002 0 0 10 20 30 40 t / t A Figure5. 25 0.0005 (a) 0.01 0.008 E 0.006 K 0.001 0.004 0.002 0.002 0 0 10 20 30 40 0.6 total (b) 0.4 E M toroidal 0.2 0 0 10 20 30 40 t / t A Figure6. 26 (a) (b) t=10, B m a x = 0.76 t=20, B m a x = 0.54 t=20, B m a x = 0.50 t=30, B m a x = 0.32 t=30, B = 0.02 t=40, B = 0.15 max max Figure7. 27 0.01 (a) m =0.0005 m =0.001 0.0001 m =0.002 m =0.004 |V |2 1e-06 1 1e-08 1e-10 0 10 20 30 40 50 0.02 (b) m =0.0005 0.015 E K 0.01 m =0.001 m =0.002 0.005 m =0.004 0 0 10 20 30 40 50 t / t A Figure8. 28 (a) (b) 1.4 1.4 1.2 1.2 1 B 1.0 pol B B 0.8 B 0.8 pol 0.6 0.6 B tor B 0.4 0.4 tor 0.2 0.2 0 0.0 0.6 0.6 0.5 0.5 0.4 0.4 J J R 0.3 R 0.3 0.2 0.2 0.1 0.1 0 0.0 0.5 0.4 Vj 0.2 Vj 0.3 V V 0.2 0.1 V V 0.1 Z Z 0 0.0 0 10 20 30 40 0 10 20 30 40 50 t / t t / t A A Figure9. 29 External Distribution Plasma Research Laboratory, Australian National University, Australia Professor I.R. Jones, Flinders University, Australia Professor João Canalle, Instituto de Fisica DEQ/IF - UERJ, Brazil Mr. Gerson O. Ludwig, Instituto Nacional de Pesquisas, Brazil Dr. P.H. Sakanaka, Instituto Fisica, Brazil The Librarian, Culham Science Center, England Mrs. S.A. Hutchinson, JET Library, England Professor M.N. Bussac, Ecole Polytechnique, France Librarian, Max-Planck-Institut für Plasmaphysik, Germany Jolan Moldvai, Reports Library, Hungarian Academy of Sciences, Central Research Institute for Physics, Hungary Dr. P. Kaw, Institute for Plasma Research, India Ms. P.J. Pathak, Librarian, Institute for Plasma Research, India Dr. Pandji Triadyaksa, Fakultas MIPA Universitas Diponegoro, Indonesia Professor Sami Cuperman, Plasma Physics Group, Tel Aviv University, Israel Ms. Clelia De Palo, Associazione EURATOM-ENEA, Italy Dr. G. Grosso, Instituto di Fisica del Plasma, Italy Librarian, Naka Fusion Research Establishment, JAERI, Japan Library, Laboratory for Complex Energy Processes, Institute for Advanced Study, Kyoto University, Japan Research Information Center, National Institute for Fusion Science, Japan Professor Toshitaka Idehara, Director, Research Center for Development of Far-Infrared Region, Fukui University, Japan Dr. O. Mitarai, Kyushu Tokai University, Japan Mr. Adefila Olumide, Ilorin, Kwara State, Nigeria Dr. Jiangang Li, Institute of Plasma Physics, Chinese Academy of Sciences, People’s Republic of China Professor Yuping Huo, School of Physical Science and Technology, People’s Republic of China Library, Academia Sinica, Institute of Plasma Physics, People’s Republic of China Librarian, Institute of Physics, Chinese Academy of Sciences, People’s Republic of China Dr. S. Mirnov, TRINITI, Troitsk, Russian Federation, Russia Dr. V.S. Strelkov, Kurchatov Institute, Russian Federation, Russia Kazi Firoz, UPJS, Kosice, Slovakia Professor Peter Lukac, Katedra Fyziky Plazmy MFF UK, Mlynska dolina F-2, Komenskeho Univerzita, SK-842 15 Bratislava, Slovakia Dr. G.S. Lee, Korea Basic Science Institute, South Korea Dr. Rasulkhozha S. Sharafiddinov, Theoretical Physics Division, Insitute of Nuclear Physics, Uzbekistan Institute for Plasma Research, University of Maryland, USA Librarian, Fusion Energy Division, Oak Ridge National Laboratory, USA Librarian, Institute of Fusion Studies, University of Texas, USA Librarian, Magnetic Fusion Program, Lawrence Livermore National Laboratory, USA Library, General Atomics, USA Plasma Physics Group, Fusion Energy Research Program, University of California at San Diego, USA Plasma Physics Library, Columbia University, USA Alkesh Punjabi, Center for Fusion Research and Training, Hampton University, USA Dr. W.M. Stacey, Fusion Research Center, Georgia Institute of Technology, USA Director, Research Division, OFES, Washington, D.C. 20585-1290 05/16/05 The Princeton Plasma Physics Laboratory is operated by Princeton University under contract with the U.S. Department of Energy. Information Services Princeton Plasma Physics Laboratory P.O. Box 451 Princeton, NJ 08543 Phone: 609-243-2750 Fax: 609-243-2751 e-mail: [email protected] Internet Address: http://www.pppl.gov