530 likes | 669 Views
J.A. Heikkinen 2 , S.J. Janhunen 1 , T.P. Kiviniemi 1 , S.Leerink 1 , M. Nora 1 , and F. Ogando 1,3 , 1 Teknillinen Korkeakoulu (TKK), Finland 2 Valtion Teknillinen Tutkimuskeskus (VTT), Finland 3 Universidad Nacional de Educación a Distancia (UNED), Spain.
E N D
J.A. Heikkinen2, S.J. Janhunen1, T.P. Kiviniemi1, S.Leerink1, M. Nora1, and F. Ogando1,3, 1 Teknillinen Korkeakoulu (TKK), Finland 2 Valtion Teknillinen Tutkimuskeskus (VTT), Finland 3 Universidad Nacional de Educación a Distancia (UNED), Spain Kinetic Simulation of Sheared Flows in Tokamaks Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 Outline • Background • Gyrokinetic full f approach for toroidal fusion plasmas. • Examples of Vlasov codes. • Description of a PIC full f code. • Testing: • Comparison to neoclassical theory. • Linear and nonlinear benchmarks of unstable modes. • Influence of noise on results. • Transport simulations in toroidal plasmas. • Future challenges of full f
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 Section IBackground
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 Turbulence and zonal flows • Self-organization and zonal flow generation with turbulence is important not only in fusion plasmas but also in other fluids like planetary atmoshpheres and solar plasma • Transport transients and related transport barriers are important, e.g., for ITER fusion reactor design • Empirical evidence from a number of tokamak facilities indicates a complex scaling law for transport transients • The scaling laws for transport barriers in tokamaks are today to a large extent not understood by physics principles
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 Sheared flows and turbulence appear in toroidal plasma simulation
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 Kinetic simulation required • All trials up to date around the world to see the spontaneous confinement transition in an experimental-like first principles plasma simulation have failed so far. • Reason to this may be that not all relevant physics like kinetic details of particle motion have been included in the simulation codes, lack of self-consistency with background evolution, or model simplifications like using fluid codes have been too crude.
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 Section IIGyrokinetic full f approach for simulation of toroidal plasmas
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 Full f simulation of tokamak plasmas • In full f simulation, one solves for the whole particle distribution function in phase space and in time either non-perturbatively, with gyrokinetics or drift-kinetically • Pioneering non-perturbative 5D particle-in-cell (PIC) simulation of magnetically confined toroidal plasmas (unrealistic me/mi) : Cheng & Okuda, ’78; Sydora, ’86; LeBrun, ’93, Kishimoto, ’94 • Eulerian and semi-Lagrangian Vlasov 3D, 4D, and 5D simulation: Cheng & Knorr, ’76; Manfredi, ’96; Sonnendruecker, ’99; Xu, ’06; Scott, ’05; Grandgirard, ’06 • Gyrokinetic 5D particle simulation: Furnish, ’99; Heikkinen, ’00, 03’; Chang, ’03.
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 Calculation flux in particle codes 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.
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 Delta f vs. full f • Delta f calculates perturbations from an assumed background distribution f0. • Powerful for small f-f0 • Linear mode analysis • “Snapshot” transport analysis • Path-breaking global transport studies for large toroidal installations • Full f calculates the whole particle distribution. • Fitting processes that perturbate strongly the particle distribution • Strong transient or long time scale transport in core or edge plasmas • Strong particle/energy sources • Full neoclassical equilibrium • Edge plasma (sheaths, wall losses, recycling, separatrix, flows) • Large MHD (sawteeth, ELMs)
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 Delta f vs. full f • Automatic extension of Delta f to full f calculation may be possible • Non-conservative effects; collisions, sources • Loading optimization • Automatic increase of the number of markers when required • Few particles (~10-100) per cell are needed for good results with small f-f0. • Full f can be realized both in PIC and Vlasov (Eulerian, semi-Lagrangian) • Vlasov simulation is free of noise and has controlled numerical resolution in phase space. • PIC is numerically flexible and straightforward to implement • Needs many particles per cell (~1000) or fine grid discretization of velocity space for an acceptable noise level or accuracy.
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 Section IIIExamples of full f Vlasov codes
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 dtF=0 t t- t x,v Z/i Interpolation (R-R0)/i Full f 5D gyrokinetic code GYSELA Grandgirard, ’06 (CEA, LPMIA-Univ, IRMA-Univ, LSIIT-Illkirch) • Semi-Lagrangian GYSELA code: - fixed grid, follows trajectories backwards, - global code - 5D GAM and ITG Cyclone tests successful
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 Full f 5D gyrokinetic code TEMPEST Xu ’06 (LLNL, Calif-Univ, GA, LBNL, PPPL) • - Based on modified gyrokinetics valid for large long-wavelength and small short-wavelength variations (Qin, 06) - Fixed grid, equations solved via a Method-of-Lines approach and an implicit backward-differencing scheme using iteration advances the system in time - Developed for circular core or divertor edge geometry - 4D neoclassical tests successful. 5D drift wave and ITG benchmarking ongoing
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 Section IVThe ELMFIRE code; an example of a full f PIC code
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 The ELMFIRE group Founded in 2000 at Euratom-Tekes Contributors from Finland Spain Holland Main affiliations VTT TKK ... but also ... CSC Åbo Akademi UNED (Spain)
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 ELMFIRE code • Full f nonlinear gyrokinetic particle-in-cell approach for global plasma simulation (present version electrostatic). • Magnetic coordinates (ψ,,) Boozer ’81. • Guiding-center Hamiltonian White & Chance, ’84. • Gyrokinetics is based on Krylov-Boholiubov averaging method in description of FLR effects (P. Sosenko, ‘01). • Adiabatic or kinetic electrons with impurities. • Parallelized using MPI with very good scalability. • Based on free software: PETSc and GSL for math calc.
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 ELMFIRE full f features • An initial canonical distribution function avoids the onset of unphysical large scale ExB flows (Heikkinen, ‘01) • Direct implicit ion polarization (DIP) and electron acceleration (DEP) sampling of coefficients in the gyrokinetic equation • Quasi-ballooning coordinates to solve the gyrokinetic Poisson equation • Versatile heat (RF, NBI, Ohmic) sources and particle sources/ recycling • Full binary collision operator
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 Poisson equation • W.W. Lee proposed “standard” model with polarization drift included in equation operator. • Ion density evaluated from ion motion without polarization drift • P. Sosenko proposes including polarization in the ion density. • Ion density evaluated from ion motion with polarization drift. Circular gyro-orbits.
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 Implementation into ELMFIRE • Solve by isolating ion polarization drift contribution to density. • That contribution is calculated implicitely every timestep using also previous values of . • The gyroaveraged electric field is interpolated from grid potential values for the ion polarization drift. 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
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 Section VBenchmarking of ELMFIRE
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 Testing ELMFIRE • Comparison to neoclassical theory in the presence of turbulence. • Frequency of GAM and Rosenbluth residual. • Neoclassical radial electric field. • Comparison to other codes has been done in the Cyclone Base cases. • Linear growth of unstable modes and their phase. • Nonlinear saturation of transport in both adiabatic and kinetic-electron case. • Comparison to experimental results. • Collaboration with IOFFE Institute and St. Petersburg Polytechnic working with the FT-2 tokamak.
Available resources • Gyrokinetic full f computation is very demanding computationally. Parallel computation is a need. • CSC (The Finnish IT Center for Science) provides shared use of high-end parallel computers. • IBM eServer cluster 1600. 512 processors with 2.2TFlops, 384GB RAM and High Performance Switch communication. • Cluster of 768 AMD OpteronTM processors up to 3.2TFlops, 1600GB RAM, Infiniband network. • Cray XT4 (Hood): 70TFlop, 70TB RAM; HP ProLiant Supercluster, 10.6 Tflop, 100 TB Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007
Geodesic Acoustic Modes • Neoclassical theory predicts GAM frequency and Rosenbluth residual. • Results show good wide agreement with theory. • Simulations done on a plasma annulus. • R=0.3-0.9 m, a=0.08 m, B=0.6-2.45 T, q=1.28-2.91, Ti=90-360 eV, ni=5.1×1019m-3 (r/a=0.75) Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 R=1.1 m, a=0.08 m, B=2.1 T, I=22 kA, parabolic ion heated n,Ti,e profiles (r=0.04m). Neoclassical radial electric field • Neoclassical radial electric field is well followed in conventional (L-mode) turbulent simulations both in radius and in time
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 Neoclassical benchmarking Parameters and boundary conditions • FT-2 like parameters, R0=0.55 m and a=0.08 m, • Itot = 55 kA, T,n,j ~ (1-(r/a)2) , • n0=5*1019 1/m3 • T0=300 eV in high T case, T0=120 eV in low T • no heating, no loop voltage • relaxing profiles, cooling by recycling on outer radius
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 Neoclassical benchmarking Number of particles must be sufficient as Er may depend on noise Er radial dependence fairly well predicted by standard neoclassical theory. However, Reynolds stress and poloidal Mach number can be important.
Linear growth of unstable modes • Test based on adiabatic Cyclone Base Case (Dimits PoP '00) • Red points from ELMFIRE, blue line: fit from article. • Figures show growth rates and typical time evolution for a mode with ki=0.3 Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007
Linear growth of unstable modes • Test based on adiabatic Cyclone Base Case (Dimits PoP '00) • Red points from ELMFIRE, blue line: fit from article. • Figures show growth rates and typical time evolution for a mode with k i=0.3 Region of linear growth Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007
Linear growth of unstable modes • Test based on kinetic Cyclone Base Case (Chen NF '03) • Filled circles and squares from ELMFIRE • Dashed line: fit for the growth rate from Chen NF ‘03. Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007
Evolution of thermal conductivity • Evolution of i is studied with nonlinear runs of Cyclone Base. • Measured at r=a/2 (q=1.4). Using kinetic electrons. R/LT=10. Weak collisionality; T(a/2)=2000 eV, n=5*1017 m-3. • Convergence requires a large number of particles per cell. Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 Section VIInfluence of noise on results
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 Influence of noise on results • Particle simulation not only covers physical density fluctuations, but it produces undesirable noise with fluxes that perturbate the solution. • Associated diffusivity can be estimated from the radial particle shift during decorrelation time. • Physical radial ion heat conductivity can be estimated from mixing-length estimate of the physical level of fluctuations, being also proportional to T3/2.
Effects on calculated conductivity • Image shows influence of strong collisionality and potential averaging on ion radial heat 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; T=100 eV, n=4.5*1019m-3 Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 Contribution of noise to the heat flux The convective noise flux prediction gives a fairly good estimate of the unphysical ion heat conductivity in the simulations
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 From noise to turbulence spectrum •Wave number spectrum from different time instants demonstrates how one moves from white noise to physical spectrum in modes.
Problematic levels of noise Fluctuations @ r=a/2 • Figure show exceptionally bad case regarding noise effects. • It is a kinetic cyclone base case with scaled parameters and low T=100eV and n=4.5*1017 m-3. • Density fluctuations remain almost constant in time. • Regression shows almost perfect N-1/2 scaling, indicating that results are dominated by noise. Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007
but not always problematic... • In FT-2, fluctuation levels are much higher (10-40%) than in scaled Cyclone Base Case (~1%). • So high perturbation level warrants the use of a full f scheme. • Image shows density fluctuations relative to flux surface average. • Relative importance of noise values can be seen in videos of both cases. Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007
Probability distribution functions • Fluctuation levels can be represented by the PDF graphs. • PDFs show fluctuation distributions over a middle flux surface. • Density values are averaged over time. • Turbulence in FT-2 takes fluctuations up to 40% in start (up) and 15% after 50µs (down). Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 Section VIITransport simulations
Case 1 under study: LH heated FT-2 • Parameters from 100 kW LH heated 22 kA FT-2 tokamak. • Case shows a rapid growth of potential and electric field and reduction of thermal conductivity interpreted as ITB formation. • Top figure diagram (r,t) of flux surface average of potential. • Bottom figure ion diffusivity at mid radius. Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 Evolution of profiles • Strong ExB flow shear is created at r=0.05 m, where a knee-point in Ti profile is found
Calculated spectra • Parameters of the FT-2 case, devised to cause the formation of Internal Transport Barrier. • B=2.2T, I=22kA, n=3.5×1019m-3, Ti=250eV, Te=300eV S(k) at r=5.1cm vs. time S(k) at r=7.5cm Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007
Case 2 under study: LH heated FT-2 • Heating phase for 100 kW LH heated 22 kA FT-2 tokamak (O8+ impurities included). Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 Evolution of large scale fluctuations • Density fluctuations plots show the formation and further destruction of macroscopic structures
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 Evolution of diffusivity • Both particle diffusivity and heat conductivity drop drastically when poloidal flow shear destroys the turbulent structures • The figures show values from the middle radius
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 Case 3 under study – Ohmic FT-2 Reproduce the FT-2 reflectometer signal I(t) with ELMFIRE I(t)=∫ w(r,θ) δn(r,θ,t) r dr dθ • W(r, θ ) = Weighting function calculated by beam tracing code using exp. data • δn(r,θ,t) = Density fluctuations simulated by ELMFIRE code Results • Poloidal velocity calculated by the shift of the power spectrum of I(t) is in reasonable agreement with the poloidal velocity measured by the Doppler reflectometer for an Ohmic discharge • Width of the power spectrum is too narrow compared to the experimental results.
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 Frequency broadening is too narrow in the simulation
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 Section VIIIFuture challenges
Numerical Flow Models for Controlled Fusion – Porquerolles, France, 16-20 April, 2007 Unresolved issues in full f • How to ensure a sufficient number of particles in all grid cells when density varies strongly in global PIC simulation • Adaptation of good grid resolution for strongly varying f in Vlasov simulations • Full collision operator in Vlasov simulations • SOL plasma simulation with sheath boundaries • Gyrokinetics with strong fluctuations