1 / 55

Plasma Simulations using Particle Codes

J.A. Heikkinen Valtion Teknillinen Tutkimuskeskus (VTT), Finland. Plasma Simulations using Particle Codes. Contributors: C.-S. Chang 1 , S. Janhunen 2 , T.P. Kiviniemi 2 , S. Leerink 2 , T. Makkonen 2 , M. Nora 2 , F. Ogando 9 , S. Sipilä 2

lysandra
Download Presentation

Plasma Simulations using Particle Codes

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. J.A. Heikkinen Valtion Teknillinen Tutkimuskeskus (VTT), Finland Plasma Simulations using Particle Codes Contributors: C.-S. Chang1, S. Janhunen2, T.P. Kiviniemi2, S. Leerink2, T. Makkonen2, M. Nora2, F. Ogando9, S. Sipilä2 1Courant Institute, New York University, USA, 2Teknillinen korkeakoulu, Espoo, Finland, 3UNED, Madrid, Spain CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008

  2. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 Outline • Introduction to Particle-in-Cell (PIC) Codes - first principles fluid dynamics, plasma and beam physics, and astrophysics and cosmology - solid state device design, accelerator physics, galaxy formation, fusion reactors, fluid flows, and magnetohydrodynamics (MHD) - PIC applications are now at a threshold where they can be used for precision prediction rather than just as quantitative indicators of system behaviour. • Gyrokinetic Particle Simulation of Magnetized Plasmas - required for microturbulence simulation - neoclassical effects, macroscopic plasma evolution?

  3. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 The challenge for modelling: span the large separation in length and time scales I LASER FUSION PLASMA IMPLOSION COLLISIONS ION MOTION PLASMONS, LASER Particle simulations Hybrid methods Hydrodynamic νDT << DT fusion ναe << Alpha slowing-down ωpi << Ion oscillations ωpe Electron oscillations 101 102 103 104 105 106 107 108 Rate of variation (nanoseconds-1) Separation of time scales requires long, large-scale PIC simulations → Cells, EF-scale machines

  4. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 Fast ignition integrated simulation code from K. Mima, Osaka University

  5. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 The challenge for modelling: span the large separation in length and time scales II TOKAMAK PLASMA MICROTURBULENCE MHD EVOLUTION TRANSPORT Fluid transport Extended MHD Gyrokinetic particle simulations τA << Alfvén transit time τs << Sound transit time τevol << MHD evolution time τtransτres Transport time Resistive time 10-6 10-5 10-4 10-3 10-2 10-1 100 101 Problem time (sec.) Separation of time scales requires long, large-scale PIC simulations → Cells, EF-scale machines

  6. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 Multiscale code Integration in CPES in Kepler framework C.S. Chang, New York Univ. XGC1 Short time gyrokinetic edge turbulence/neocl Collab. GA ELITE, etc Linear ELM Stabilitycheck XGC0 Experimental time evolution of kinetic edge equilibrium M3D & NIMROD Nonlinear ELM crash Joint, CEMM Core Gyrokinetic code Collab. GPS + RF-edge interaction, energetic particles, impurity, radiation, 3D magnetic perturbation physics, and core MHD effects are to be built into XGC0 DEGAS2 Neutrals Atomic physics Plasma-Surface Interaction code

  7. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 Macroscopic flows and turbulence

  8. XGC Gyrokinetic particle motions in the plasma edge (poloidal plane) ions electrons V|| >0 V|| <0 CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008

  9. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 Basis of PIC plasma simulation Acceleration and increment of velocity Calculation of forces from fields and velocity Computation of electric field. Magnetic is given. Initial step with =0 Displacements and new positions. Boundary conditions. Calculation of density. Current profile fixed. Resolution of Poisson equation for the electrostatic potential.

  10. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 1D Electrostatic case with Poisson Eq GRID j’th particle (q,m) Leap-frog algorithm for particle co-ordinates cell d i-2 i-1 i i+1 CIC Poisson solver NGP

  11. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 Superparticles • It is not practical to simulate all particles in laboratory (1020-1033 m-3) or in space (1018 km-3) plasmas. • Therefore, each simulation particle (a superparticle) represents many particles of a real plasma (grid spacing ~ λD) Collisional effects due to Coulomb field are not modelled correctly with superparticles Remedy: Use finite size particles to smear out collisions: With the size ~ several λD, this method retains collective interaction

  12. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 Finite number of simulation particles creates noise 3D solution, potential and density fluctuation amplitudes at a given position • Physical solution in any PIC simulation is always accompanied by a discrete particle noise • In the picture, regression shows almost perfect N-1/2 scaling, indicating that results in this case are dominated by noise. • Discreteness may, however, be beneficial in finding proper asymptotics, e.g., in nonlinear Landau damping

  13. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 From noise to turbulence spectrum •Wave number spectrum from different time instants demonstrates how one moves from white noise to physical spectrum in modes.

  14. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 Finite size filtering • CIC provides filtering which creates a transfer function to noise, too.

  15. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 Stability and accuracy • Differencing gives time constraint ωpeΔt<2, while retaining accuracy for the plasma waves requires ωpeΔt<0.2. • Numerical instabilities may arise due to (aliasing) - discrete spatial grid - finite-sized time steps • High frequency modes appear not correctly on the mesh and modes k and k+2πn/Δx cannot be distinguished. • If k+2πn/Δx resonates with any of the natural modes of the plasma, system can become numerically unstable. • Remedy: Make either aliased modes Landau damped (λD>Δx/π) or use finite size particles (size ≥ 0.7Δx).

  16. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 Landau damping • Truely dissipationless Vlasov-Poisson simulation cannot be realized due to the effects of numerical gridsize • With Δx<λD, numerical collision time however longer than ωpe-1nλD • When the size of phase-space filaments reaches the interparticle separation, envelope oscillations disappair, and wave goes on propagating at constant amplitude as in O’Neil’s analysis (1965) Discretness in PIC simulations helps to describe the true asymptotics of nonlinear Landau damping. The continuum Eulerian simulations in principle fail in this task Carbone, Valentini, & Veltri, 2007

  17. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 A B PIC-MCC: Binary collisions I (Takizuka, 1977) • Scattering between two test particles randomly selected from a grid (left). • After that, the relative velocity vector (top-right) is rotated • The new velocity vectors for both particles are recalculated so that momentum is conserved

  18. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 Species A Species A Species B Species B Binary collisions II • Pairing of particles for binary collisions is done within a grid cell • For A-B collision, loop all the particles of species A, and for every particle, find a random partner from species B • This works if the particles have a same weight, i.e., each simulation particle corresponds to the same amount of true particles. A-B collisions vs. B-A collisions

  19. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 A A B 1. 1. B 2. 1. 2. 3. 1. Binary Collisions 2. B 2. A 3. Binary collisions III • The problem of different weight factors can be solved by applying the binary collision operator between groups of test particles that represent an equal amount of actual particles KANA= KBNB NA, NB are the weights All particles from the first group collide with an imaginary particle representing the state of the second group, and all particles from the second group collide with and imaginary particle representing the state of the first group

  20. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 Collision effects easy to test Maxwellian energy distribution; all the test particles started with an energy of 10 keV Both particles started with a very low kinetic energy in an environment of average energy of 10 keV (electron thermalizes faster)

  21. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 Laboratory plasmas are collisional • ES, PD (Plasma Device), and OOPIC (relativistic, electromagnetic, 2d3v, collisions) code family (C.K. Birdsall, PTSG, US Berkeley) • Plasma-material boundaries and SOL an important application in magnetized fusion BIT1 code result (Tskhakaya, 2005) for the Mach number and charge density along the SOL magnetic field up to the plasma-wall interface • MP and DS regions cannot be modelled in fluid or in gyrokinetic approximation. • Binary collisions for plasma particles. • Null-collision model for charged-neutral particle collisions. Self-consistent neutrals (Matyash, 2003) • Either gyro-averaged or Boltzmann electrons • Either 1d3v or 2d3v • Sometimes PSI included

  22. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 Simplifications for enhancing code performance • Take electromagnetic field as given (standard particle code) • Solve only for perturbations δf • Advance only orbit invariants, not track particle orbits • Either freeze ions or use adiabatic approximation for electrons or ions • Solve for quasineutrality • Advance particle co-ordinates and fields implicitly • Take drift-kinetic limit (guiding-center code) • Apply gyrokinetic approximation

  23. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 GYROKINETIC PARTICLE SIMULATION

  24. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 Gyrokinetic (GK) method Standard nonlinear GK ordering: Distribution function F can evolve arbitrary far from its initial value.

  25. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 Gyrokinetic model • With the proper averaging process, FLR effects are still kept. • Coordinates change to gyrocenter. Angle of gyration goes by averaging. • This introduces changes in equations that use real space coordinates (i.e. Poisson)

  26. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 Larmor ring discretization • Four-point averaging procedure (W.W. Lee, 1987) often used for estimates of gyroaverages <..> • Sampling of density ñi on the grid from each particle uses the same four-point averaging ∂f/∂μ is difficult to evaluate with particle methods; usually f is taken as a Maxwellian ∫<Φ> ∂f/∂μ dv can then be evaluated with either one, two, or three Larmor rings with the corresponding weights from the Maxwellian

  27. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 Direct implicit solver for ion polarization • Ion polarization density can also be calculated directly from simulation particles (requires different guiding-center transformation, JCP 2008) • Isolate implicitly ion polarization drift contribution to density. Larmor circle i’th subparticle cloud y k i py+ px0 px- xp,yp px+ py0 ds k py- x p’th point on the Larmor orbit of the k’th ion gyro-center of the k’th ion

  28. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 Implicit solver for electron parallel nonlinearity • E field in ∆xa is calculated at advanced time but at the position after free streaming • We demand that |∆xfs|>>|∆xa|. It is a constraint in Δt A.B Langdon et al (1983)

  29. Gyrocentre motion • Gyrokinetic model deals with particle gyrocenters, which present smoother transversal trajectories. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008

  30. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 Magnetic background

  31. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 Initialization • Loading particles for local Maxwellian velocity distribution and for wished density and temperature profiles usually satisfactory • However, in tokamaks neoclassical equilibrium conditions difficult to initialize • Finite orbit effects can be accounted for by initialising particles according to orbit invariants (Hively, 1981) • Electrons to be loaded on ion gyro-orbits • In full f calculation, use equal weights for particles among a given species

  32. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 Simplified boundary conditions

  33. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 Plasma-wall boundary problems

  34. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 Nonlinear GK (kinetic or gyro-fluid) studies • Zonal flows, geodesic acoustic modes, neoclassical equilibrium • Most drift-wave-type fluctuations like ITG, TEM, TIM, ETG, universal and dissipative drift instabilities. • Interchange turbulence • Tearing and internal kink instabilities • Microtearing and drift-tearing mode • Drift-Alfvén and kinetic shear-Alfvén waves • Compressional Alfvén waves

  35. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 Gyrokinetic particle codes for simulation of magnetized toroidal plasmas

  36. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 ∂f and full f particle codes for magnetized fusion • A number of ∂f PIC codes (GEM,GT3D, GYGLES, GTC, EUTERPE, ORB5, UCAN, and many others) have paved the way over almost 20 years for precision prediction for turbulent transport in tokamaks and stellarators • Full f PIC codes (ELMFIRE, TPC, XGC1, and some others) and full f semi-Lagrangian code (GYSELA) approach precision prediction for turbulent transport and macroscopic plasma evolution in tokamaks • A number of full f particle codes (ASCOT, SELFO, XGC0, and many others) capable for consistent calculations of fields and particles in tokamaks and stellarators

  37. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 3d2v particle-in-cell full f code ELMFIREJ.A. Heikkinen, S. Janhunen, T. Kiviniemi, S. Leerink, M. Nora, and F. Ogando (TKK,VTT,UNED) • Global full-F GK code • Drift-kinetic electrons and gyrokinetic ions and impurities. • Heat (RF, NBI, Ohmic) and particle sources + recycling • Full binary collision operator

  38. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 R=1.1 m, a=0.08 m B=2.1 T, I=22 kA, n(0)=41019 m-3 Ti,e(0)=180 eV Neoclassical radial electric field in turbulent FT-2 tokamak Simulation can agree with the neoclassical radial electric field in L-mode turbulent simulations (Vφfrom simulation)

  39. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 ELMFIRE benchmarking

  40. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 Influence of noise on results • Particle simulation noise not only masks physical density fluctuations, but it produces undesirable fluxes that perturbate the neoclassical equilibrium. • Associated diffusivity can be estimated from the radial particle shift during decorrelation time. • Radial ion conductivity can be calculated from mixing-lenght estimate of the physical level of fluctuations, being also proportional to T3/2.

  41. Effects on calculated conductivity • Image shows influence of collisions and potential averaging on ion radial conductivity. • Collisionless cases show residual noise conductivity. • Noise is filtered out by averaging potential over flux surface. • So far noise is reduced by “brute force” ... higher N! Scaled Cyclone Base Case with kinetic electrons χiNC CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008

  42. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 Saturation with adiabatic electrons I • Evolution of i is studied with nonlinear runs of Cyclone Base (Dimits et al., PoP ’00). • Measured at r=a/2 (q=1.4). R/LT=8.28. No collisions; T(a/2)=2000 eV, n=5.5*1019 m-3. • Adiabatic electrons: ne=n0(r)*(1+e(Φ-<Φ>)/kT(r)) χi (m2/s) Why turbulent transport saturates at such a low value? <χ>i << χgB

  43. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 Saturation with adiabatic electrons II

  44. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 Saturation with adiabatic electrons III Due to turbulence growth, ion temperature profile gets relaxed. In regions where Ti decreases, ion orbits narrow, with an opposite behaviour where Ti increases. As a consequence, ion polarization has to compensate the resulting perturbation in ion density → Electrostatic potential gets a bipolar perturbed structure This leads to suppressed turbulence by sheared zonal flows. Cooled Heated

  45. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 FT-2 simulations in L-mode Turbulent damping of zonal flows

  46. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 Steep pressure profile causes spin-up Strong heating (100 kW) causes radial electric field to deviate from neoclassical Sheared poloidal flow suppresses transport R=1.1 m, a=0.08 m B=2.1 T, I=22 kA, n(0)=41019 m-3 Ti,e(0)=250 eV

  47. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 3d2v Particle-in-cell full f code XGC1 C.S. Chang (New York University, Center for Plasma Edge Simulation (CPES)) Kinetic XGC1 - Global full-f GK code - short time simulation of turbulence transport - consistently with neoclassical physics Kinetic XGC0 - long time evolution of neoclassical ion-electron- neutral equilibrium (full-f, PIC) • Includes 3D magnetic perturbation effects • Turbulence transport to be imported from XGC1 coupling

  48. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 First 3D electrostatic solution across entire edge Early time 3D  (Neoclassical + Fluctuations) XGC1 Equilibrium  Negative potential well is forming in H-layer (dark blue) Positive potential hill in scrape-off (yellow/brown)

  49. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 Strongly sheared ExB flow profiles in the entire H-mode edge from XGC-1 Wall (eV) N GAM

  50. CEA-EDF-INRIA school in Numerical Analysis; Numerical Models for controlled Fusion – Nice, France, 8-12 September, 2008 Bootstrap current calculation across separatrix with XGC Good agreement with analytical models Agreement not so good with steeper pedestal n = 1.5 cm n = 0.75 cm Simulation with 150000 ions and electrons in XGC0

More Related