AN INTRODUCTION TO THE NUMERICAL MODELING OF FUSION PLASMAS Herv´ e Guillard INRIA Sophia-Antipolis & LJAD, UMR 7351, 06100 Nice, France CEMRACS 2013 Herv´ e Guillard (INRIA & LJAD) 22/08/2013 1 / 36
Outline Fusion Plasmas 1 Kinetic Models 2 Fluid models 3 The MHD limit 4 Herv´ e Guillard (INRIA & LJAD) 22/08/2013 2 / 36
Outline Fusion Plasmas 1 Kinetic Models 2 Fluid models 3 The MHD limit 4 Herv´ e Guillard (INRIA & LJAD) 22/08/2013 3 / 36
Deuterium-Tritium Reaction Herv´ e Guillard (INRIA & LJAD) 22/08/2013 4 / 36
Nuclear reactions Herv´ e Guillard (INRIA & LJAD) 22/08/2013 5 / 36
Fusion on earth ??? Essentially two studied technologies : Inertial confinement Magnetic confinement Herv´ e Guillard (INRIA & LJAD) 22/08/2013 6 / 36
Inertial confinement technology relies on fast and violent heating of fusion targets. Herv´ e Guillard (INRIA & LJAD) 22/08/2013 7 / 36
Target Ablator heating Ablator is illuminated by powerful - X rays and reaches a plasma state Herv´ e Guillard (INRIA & LJAD) 22/08/2013 8 / 36
Ablator expansion & DT Compression Ablator is accelerated outwards By rocket effect DT is accelerated inwards (v ∼ 300 km/s) Herv´ e Guillard (INRIA & LJAD) 22/08/2013 9 / 36
Deuterium - Tritium Compression DT ice is compressed up to 1000 times its initial density ( 20 time lead density) Herv´ e Guillard (INRIA & LJAD) 22/08/2013 10 / 36
Hot spot Ignition Central gas is heated up to several 10 6 of degrees When ρ R ∼ 0 . 3 g / cm 2 the gas is self heated by alpha particles Fusion reactions start and propa- gate into DT ice Note : 1 mg D-T → 340 MJ : Fu- sion is equivalent to combustion of 10 kg coal Herv´ e Guillard (INRIA & LJAD) 22/08/2013 11 / 36
The other way to fusion : magnetic confinement Objective : confine hot D-T plasma by strong magnetic fields Poloidal coils Toroidal coils Herv´ e Guillard (INRIA & LJAD) 22/08/2013 12 / 36
Tokamaks : magnetic field lines vacuum chamber Tokamaks design in the 1950s by Igor Tamm and Andrei Sakharov Herv´ e Guillard (INRIA & LJAD) 22/08/2013 13 / 36
ITER Tokamak ITER (International Thermonuclear Experimental Reactor) in construction in Cadarache, construction begins in 2010, first plasma in 2020, first fusion plasma in 2027 840 m 3 150M O K 13 Tesla = 200 000 x earth magnetic field Herv´ e Guillard (INRIA & LJAD) 22/08/2013 14 / 36
Numerical simulations Necessary for Equilibrium computation : steady state and control of the machines prediction of the possible occurrence of instabilities magnetic instabilities hydrodynamic instabilities determination of the value of key parameters e.g transport coefficients due to turbulence understand and explain physical phenomena Very large number of different numerical models Herv´ e Guillard (INRIA & LJAD) 22/08/2013 15 / 36
Outline Fusion Plasmas 1 Kinetic Models 2 Fluid models 3 The MHD limit 4 Herv´ e Guillard (INRIA & LJAD) 22/08/2013 16 / 36
Kinetic models Vlasov or Boltzmann eq for each particle (ions, electrons, neutral) ∂ f s � ∂ t + div x ( v f s ) + div v ( F f s ) = C s = C ss ′ Force field F = e s m s [ E + v × B ] (note div v F = 0) E , B given by Maxwell equations: ∂ B Υ = � e s f s ( x , v , t ) dv 3 ∂ t + curl E = 0 � J = � e s − 1 ∂ E f s v ( x , v , t ) dv 3 ∂ t + curl B = µ 0 J � c 2 ε 0 div E = Υ div B = 0 Herv´ e Guillard (INRIA & LJAD) 22/08/2013 17 / 36
Kinetic models Self-consistent and closed model extremely heavy from computational point of view 6 D model : 3 space dimensions, 3 velocity dimensions covers huge range of time and space scales some simplifications possible : Ampere law, Quasi-Neutral assumption, Electrostatic assumption, gyrokinetic theory Used for some very specific tasks Herv´ e Guillard (INRIA & LJAD) 22/08/2013 18 / 36
Kinetic models : one example of applications Anomalous transport in tokamaks Radial diffusion of mass and tem- perature exceed by order of magni- tude “laminar” values due to micro-turbulence and micro- instabilities : Estimate the value of turbulent transport coefficients by direct sim- ulation Development of gyrokinetic codes (e.g GYSELA code of CEA) Herv´ e Guillard (INRIA & LJAD) 22/08/2013 19 / 36
Outline Fusion Plasmas 1 Kinetic Models 2 Fluid models 3 The MHD limit 4 Herv´ e Guillard (INRIA & LJAD) 22/08/2013 20 / 36
From kinetic to fluid Idea of fluid models : instead of computing the whole distribution function f s compute only a small number of its moments, i.e : Fluid variables are moments of the distribution function f s • fluid density � f s ( x , v , t ) dv 3 n s ( x , t ) = • fluid velocity � f s ( x , v , t ) v dv 3 n s ( x , t ) u s ( x , t ) = Herv´ e Guillard (INRIA & LJAD) 22/08/2013 21 / 36
From kinetic to fluid • Pressure tensor � f s ( x , v , t )( v − u s ) ⊗ ( v − u s ) dv 3 P s ( x , t ) = m s • scalar pressure = 1/3 trace of pressure tensor p s ( x , t ) = m s � f s ( x , v , t ) | ( v − u s ) | 2 dv 3 3 • temperature T s ( x , t ) = p s ( x , t ) n s ( x , t ) • Energy flux Q s ( x , t ) = m s � f s ( x , v , t ) | v | 2 v dv 3 2 More complex fluid models, e.g 14 Moment closure by D. Levermore Herv´ e Guillard (INRIA & LJAD) 22/08/2013 22 / 36
From kinetic to fluid Mass conservation First velocity moment of the Boltzmann (Vlasov) equation [ ∂ f s � ∂ t + div x ( v f s ) + div v ( F f s ) = C s ] gives f s dv 3 def v f s dv 3 def � = n s � = n s u s C s dv 3 = 0 � div v ( F f s ) = 0 � ∂ n s ∂ t + div x ( n s u s ) = 0 Herv´ e Guillard (INRIA & LJAD) 22/08/2013 23 / 36
Fluid models proceeding in the same way for the other moments : m s v × Boltzmann and integrating, gives momentum (velocity) equation 1 2 m s | v | 2 × Boltzmann and integrating, gives energy (temperature) equation Herv´ e Guillard (INRIA & LJAD) 22/08/2013 24 / 36
Fluid models s ∈ ion , electron , neutral particle ∂ n s ∂ t + div( n s u s ) = 0 m s ( ∂ n s u s � + div( n u s ⊗ u s )) + ∇ p s + div x Π s − n s e s ( E + u s × B ) = R ss ′ ∂ t s ′ � = s ∂ n s 2 | u s | 2 + 3 ∂ t (m s 2 p s ) + div Q − e s u s . E = Q s Herv´ e Guillard (INRIA & LJAD) 22/08/2013 25 / 36
The 2-fluid model Much simpler than kinetic model (3D) Clear mathematical structure (1 compressible Navier-Stokes system per species coupled by force terms) But needs closure assumptions Questionable in tokamaks very large disparity in length and time scales still costly e.g for a ion-electron model 10 dof per mesh point not really used but good starting point for subsequent approximations Herv´ e Guillard (INRIA & LJAD) 22/08/2013 26 / 36
From 2-fluid to one Fluid models Normalization units time τ obs velocity u ∗ lenght L ∗ = u ∗ τ obs ion gyro-period τ c i = m i / Z i eB electron gyro-period τ c e = m e / eB plasma period τ pe = ( ε 0 m e / n e e 2 ) 1 / 2 τ c e /τ pe ∼ 1 collision τ coll ∝ (Λ / lnΛ) τ pe Non dimensional parameters Λ = 4 π 3 n λ 3 ε i = τ c i /τ obs M i = u ∗ /ν T i , D , µ = m e / m i Relevant asymptotic regimes ε i → 0 Λ → ∞ µ = m e / m i → 0 etc Give a huge set of one fluid models Herv´ e Guillard (INRIA & LJAD) 22/08/2013 27 / 36
Outline Fusion Plasmas 1 Kinetic Models 2 Fluid models 3 The MHD limit 4 Herv´ e Guillard (INRIA & LJAD) 22/08/2013 28 / 36
From 2-fluid to one Fluid models The MHD scalings MHD : violent instabilities affecting the whole plasma region : x ∗ = e . ga minor radius fast event : speed of MHD waves : Alfen velocity B v A = √ µ 0 ρ u ∗ = O ( ν T i ) ∼ O ( v A ) MHD asymptotic limit ε i = τ c i τ → 0 , M i = O (1) Herv´ e Guillard (INRIA & LJAD) 22/08/2013 29 / 36
THE MHD MODELS Ideal MHD Zero order model in ε i ( ε i = 0) ∂ ∂ t ρ + div ( ρ u ) = 0 ∂ ∂ t ρ u + div ( ρ u ⊗ u ) + ∇ p − ( J × B ) = 0 E + u × B = 0 need the current J → Ampere’s law J = curl B need the magnetic field B → Faraday’s law ∂ B ∂ t + curl E = 0 ⇒ autonomous system for a One Fluid model Herv´ e Guillard (INRIA & LJAD) 22/08/2013 30 / 36
THE MHD MODELS Ideal MHD - summary ∂ ∂ t ρ + div ( ρ u ) = 0 ∂ t ρ u + div ( ρ u ⊗ u ) + ∇ p + 1 ∂ ( B × curl B ) = 0 µ 0 ∂ B ∂ t − curl ( u × B ) = 0 + energy equation Herv´ e Guillard (INRIA & LJAD) 22/08/2013 31 / 36
One example of hydrodynamic instabilities in magnetized plasmas 2 scalar variables : ion density n ; ion parallel velocity u || u = u || b + u ⊥ ∂ ∂ t ρ + div ( ρ u ) = 0 ∂ t ρ u || + div ρ u || u + ∇ || p = ρ u . D b ∂ Dt with (this is the ideal Ohms law solved for the ⊥ velocity) u ⊥ = E × B B 2 B given, E computed by the adiabatic assumption Herv´ e Guillard (INRIA & LJAD) 22/08/2013 32 / 36
Recommend
More recommend
