Transcript Document

EECE695: Computer Simulation (2005)
Particle-in-Cell Techniques
HyunChul Kim and J.K. Lee
Plasma Application Modeling Group,
POSTECH
 References:
• Minicourse by Dr. J. P. Verboncoeur (PTS Group of
UC Berkeley) in IEEE International Conference on
Plasma Science (2002)
• “Plasma Physics via Computer Simulation” by C.K.
Birdsall and A.B. Langdon (Adam Hilger, 1991)
PIC Overview
 Applications of PIC model
• Basic plasma physics: waves and instabilities
• Magnetic fusion
• Gaseous discharges
• Electron and ion optics
• Microwave-beam devices
• Plasma-filled microwave-beam devices
PIC Overview
 PIC Codes Overview
• PIC codes simulate plasma behavior of a large
number of charges particles using a few
representative “super particles”.
• These type of codes solve the Newton-Lorentz
equation of motion to move particles in
conjunction with Maxwell’s equations (or a
subset).
• Boundary conditions are applied to the
particles and the fields to solve the set of
equations.
• PIC codes are quite successful in simulating
kinetic and nonlinear plasma phenomenon like
ECR, stochastic heating, etc.
PIC-MCC Flow Chart
• Particles in continuum space
( x,v)i
• Fields at discrete mesh locations in space (E,B) j
• Coupling between particles and fields
I
II
V
IV
III
IV
Fig: Flow chart for an explicit PIC-MCC scheme
I. Particle Equations of Motion
 Newton-Lorentz equations of motion
d
mu  F  q ( E  v  B )
dt
u  v
d
xv
dt
 In finite difference form, the leapfrog method
Fig: Schematic leapfrog integration
I. Particle Equations of Motion
ut  t / 2  ut t / 2 q t ut  t / 2  ut t / 2
t
 (E 

B
)
t
t
m
2
xt  t  xt ut  t / 2
 t  t / 2
t

• Second order accurate
• Requires minimal storage
• Requires few operations
• Stable for wp t  2
I. Particle Equations of Motion
• Boris algorithm

u u
t  t / 2
qtEt

2m
t
q

t
E
u   u t  t / 2 
2m
u  u
q

 t (u   u  )  Bt
t
2 m
u
u  u

u  u
u

t
q
B
t
t  bˆ tan( ) 
2
2 t m
t
I. Particle Equations of Motion
Finally,
u'  u   u   t t
2t t
u  u  u' 
1 tt  tt


2t t
u' 
1 tt  tt
u

u'
u

u  tt
II. Particle Boundary
 Absorption
• Conductor : absorb charge, add to the global σ
• Dielectric : deposit charge, weight q locally to mesh
 Reflection
• Physical reflection
reverse xbc  x  x bc  x
 vx  vx
• Specular reflection
 1st order error
v
t  t / 2*
 v
t  t / 2
qEt 2 | x t  xbc |

(
 t )
t  t / 2
m
|v
|
 Thermionic Emission
 Fowler-Nordheim Field Emission
 Child’s Law Field Emission
 Gauss’s Law Field Emission
II. Particle Boundary
 Secondary electron emission
+ io n
photon ,  ex
–  electron
– se
• Photoemission
• Ion impact secondary emission
• Electron impact secondary emission
 Important in processes related to high-power
microwave sources
III. Electrostatic Field Model
 Possion’s equation
  (x, t )   (x, t ),
• Finite difference form in 1D planar geometry
 j 1  2 j   j 1
x 2
j

,

 Boundary condition : External circuit
Fig: Schematic one-dimensional bounded plasma
with external circuit
III. Electrostatic Field Model
E0 
0

EJ  
J

From Gauss’s law,
E1/ 2
 t0  1t  0t   0t x 2


x

• Voltage driven series RLC circuit
From Kirchhoff’s voltage law,
d 2Q(t )
dQ(t ) Q(t )
L

R

2
dt
dt
C
 V (t )   J (t )   0 (t )
 0t   0t  t
Q t  Q t  t
  J plasma dt 
t  t
A
t
• Short circuit
0 (t ) is specified, J (t )  0
• Open circuit
 
t
0
t  t
0

t
t  t
J plasma dt
IV. Coupling Fields to Particles
 Particle and force weighting
: connection between grid and particle quantities
• Weighting of charge to grid
• Weighting of fields to particles
grid point
a point charge
IV. Coupling Fields to Particles
• Nearest grid point (NGP) weighting
 fast, simple bc, noisy
• Linear weighting
: particle-in-cell (PIC) or cloud-in-cell (CIC)
 relatively fast, simple bc, less noisy
• Higher order weighting schemes
 slow, complicated bc, low noisy
Quadratic spline
NGP
1.0
Linear spline
Cubic spline
0.5
0.0
xi  2x
Position (x)
xi  x
xi
xi  x xi  2x
Fig: Density distribution function of a particle at xi
for various weightings in 1D
IV. Coupling Fields to Particles
Areas are assigned to grid
points; i.e., area a to grid
point A, b to B, etc
Fig: Charge assignment for linear weighting in 2D
V. Monte-Carlo Collision Model
• The MCC model statistically describes the
collision processes, using cross sections for each
reaction of interest.
• Probability of a collision event
Pi  1  exp[ng (x) T ( i )i t ]
where  T ( i )   j  j ( i )
• For a pure Monte Carlo method, the timestep is
chosen as the time interval between collisions.
ti  
ln(1  R)
ng (x) T ( i )i
where 0< R< 1is a uniformly distributed random
number.
However, this method can only be applied when
space charge and self-field effects can be neglected.
V. Monte-Carlo Collision Model
• There is a finite probability that the i-th particle will
undergo more than one collision in the timestep.
Thus, the total number of missed collisions
(error in single-event codes)

2
P
r   Pi k  i
1  Pi
k
Hence, traditional PIC-MCC codes are constrained by
vmaxt  1 for accuracy.
where vmax  max (ng ( x )) max ( T ( ) )
x

V. Monte-Carlo Collision Model
• Computing the collision probability for each particle
each timestep is computationally expensive.
→ Null collision method
1. The fraction of particles undergoing a collision each
time step is given by
PT  1  exp[ maxt ].
2. The particles undergoing collisions are chosen at
random from the particle list.
3. The type of collisions for each particle is determined
by choosing a random number, 0  R   max .
Null collision
Collision
type 3
Collision type 2
Collision type 1
Fig: Summed collision frequencies for the null collision
method.