Transcript Document

ECE586: Advanced E&M Simulation (2004)
On PDX1 Program
HyunChul Kim and J.K. Lee
Plasma Application Modeling, POSTECH
2004. 9. 16
 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)
A Series of XPDX1*
 XPDx1: X window (using xgrafix library),
Plasma Device, 1 Dimensional (1d3v), Bounded
(with external circuit drive), Electrostatic Code
• XPDP1 (x=P) : Planar Configuration
• XPDC1 (x=C) : Cylindrical Configuration
• XPDS1 (x=S) : Spherical Configuration
r
~
LRC
Computation Space
* Developed by PTS group, UC Berkeley
All are available at http://ptsg.eecs.berkeley.edu
PIC Overview
 PIC Codes Overview
• Plasma behavior of a large number of charges
particles are simulated by using a few
representative “super particles”.
• PIC codes solve fundamental equations, the
Newton-Lorentz equation of motion to move
particles in conjunction with Maxwell’s
equations (or a subset) with few approximations.
• PIC codes are quite successful in simulating
kinetic and nonlinear plasma phenomenon like
ECR, stochastic heating, etc.
Computer Simulation of Plasma
Kinetic
Description
Vlasov,
Fokker-Planck
Codes
Particle
Codes
Fluid
Description
Hybrid
Codes
 Particle codes
• The particle-particle model
• The particle-mesh model
• The particle-particleparticle-mesh model
Fluid
Codes
XPDx1 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
Subcycling : ti slow  kt fast , k  1
IV
III
IV
Fig: Flow chart for an explicit PIC-MCC scheme
V
I. Particle Equations of Motion
 Newton-Lorentz equations of motion
d
mv  F  q (E  v  B)
dt
d
xv
dt
 In finite difference form, the leapfrog method
vt  t / 2  vt t / 2 q t vt  t / 2  vt t / 2
 (E 
 Bt )
t
m
2
x t  t  x t
 vt  t / 2
t
• Second order accurate
• Stable for wp t  2
I. Particle Equations of Motion
• Boris algorithm

v v
t  t / 2
qtE t

2m
t
q

t
E
v   vt  t / 2 
2m
v  v
q


( v   v  )  Bt
t
2m
v
v  v

v  v
v

t
q
B
t
t  bˆ tan( ) 
2
2m
t
I. Particle Equations of Motion
Finally,
v'  v   v   t t
2t t
v  v  v' 
1 tt  tt


2t t
v' 
1 tt  tt
v

v'
v

v  t t
XPDx1 Flow Chart
I
II
V
IV
III
IV
Fig: Flow chart for an explicit PIC-MCC scheme
II. Particle Boundary
 Absorption
• Conductor : absorb charge, add to the global σ
 Secondary electron emission
+ io n
– 
electron
–
se
• Ion impact secondary emission
• Electron impact
secondary emission
XPDx1 Flow Chart
I
II
V
IV
III
IV
Fig: Flow chart for an explicit PIC-MCC scheme
III. Electrostatic Field Model
 Maxwell’s equation in vacuum
 E  
B
,   D   , D  0 E
t
D
 H  J 
,   B  0 , B  0 H
t
• In electrostatics,
  E  0  E  
  2  

0
(Poisson’s equation)
Or Gauss’ law
 D  d s 
S
V
 dV Qenclosed
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
J
EJ  

0
E0 

From Gauss’s law,
 t0  1t  0t   0t x 2
E1/ 2 

x

t
Q t  Q t  t
t
t  t
 0   0   J plasma dt 
t  t
A
• Short circuit
0 (t ) is specified, J (t )  0
• Open circuit
 
t
0
t  t
0

t
t  t
J plasma dt
III. Electrostatic Field Model
• 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 )
― One second order difference equation
where
XPDx1 Flow Chart
I
II
V
IV
III
IV
Fig: Flow chart for an explicit PIC-MCC scheme
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
x  S ( x  xi )
1.0
Quadratic spline
NGP
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
 ( X j )   qi S ( X j  xi )
i
Fi  qi x E j S ( X j  xi )
j
IV. Coupling Fields to Particles
 Weighting in 1D
• For linear weighting in cylindrical coordinates,
(0<j<N)
XPDx1 Flow Chart
I
II
V
IV
III
IV
Fig: Flow chart for an explicit PIC-MCC scheme
Collisions
 Electron-neutral collisions
• Elastic scattering (e + A → e + A)
• Excitation (e + A → e + A*)
• Ionization (e + A → e + A+ + e)
 Ion-neutral collisions
• Elastic scattering (A+ + A → A+ + A)
• Charge exchange (A+ + A → A + A+)
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 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 T ( i )i
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.
Since XPDx1 deals with only one collision in the
timestep, the total number of missed collisions

Pi 2
r   Pi 
.
1  Pi
k
k
Hence, XPDx1 is constrained by vmaxt  1
for accuracy.
where vmax  ng max ( T ( ) )

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. Nc  N  PT
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.
Numerical Parameters
 Choose Δx and Δt to resolve the smallest important
physical feature
 Require Δx < Debye length, sheath length, wave
length, Larmor radius, boundary feature, etc.
 Require t  x / max for all species (“particle
Courant”) for accurate sampling of fields.
 Require t  0.2 / wp for accuracy of explicit leap
frog mover or for accuracy when space charge forces are
important.
 Require t  1 / when collisions are important.
 Require # of superparticles per cell > 10. It should
be larger in simulations where particles remain trapped
for long times.
Example of XPDP1 Input File
RF DISCHARGE(IN MKS UNITS)
Voltage-driven with electron-neutral collisions (Argon atom)
-nsp---nc---nc2p---dt[s]---length[m]--area[m^2]--epsilonr---B[Tesla]---PSI[D]-2 400 1.8e6 8e-12 0.03
0.01
1.0
0 .0
0.0
-rhoback[C/m^3]---backj[Amp/m^2]---dde--extR[Ohm]--extL[H]---extC[F]---q0[C]0.0
0.0
0.0 0.0
0.0
1.0
0.0
-dcramped--source--dc[V|Amp]--ramp[(V|Amp)/s]---ac[V|Amp]---f0[Hz]--theta0[D]0
v
0.0
0.0
100
13.56e6
0
--secondary--e_collisional---i_collisional---reflux---nfft--n_ave--nsmoothing--ntimestep-1
1
2
0
0
276549 6
0
--seec(electrons)---seec(ions)---ion_species----Gpressure[Torr]---GTemp[eV]---imp-0.0
0.2
2
100e-3
0.026
0
---GAS----psource--nstrt-1
0
0
SPECIES 1
----q[C]-------m[Kg]---j0L[Amp/m^2]---j0R[Amp/m^2]----initn[m^-3]----k--1.602e-19 9.11e-31
0.0
0.0
5e14
1
--vx0L[m/s]---vxtL[m/s]--vxcL[m/s]---vxLloader(0=RNDM,1=QS)-0.0
4.19e5
0.0
1
--vx0R[m/s]---vxtR[m/s]--vxcR[m/s]---vxRloader
0.0
4.19e5
0.0
1
--v0y[m/s]---vty[m/s]---vyloader---v0z[m/s]---vtz[m/s]--vzloader-0.0
4.19e5
1
0.0
4.19e5
1
--nbin----Emin[eV]----Emax[ev]---maxnp—
200
0.0
20.0
300000
-For-Mid-Diagnostic---nbin----Emin[eV]---Emax[eV]----XStart--XFinish—
300
0.0
20.0
0.0
0.03
SPECIES 2
----q[C] ------m[Kg]---j0L[Amp/m^2]---j0R[Amp/m^2]----initn[m^-3]----k1.602e-19 6.68e-26
0.0
0.0
5e14
1
-vx0L[m/s]---vxtL[m/s]--vxcL[m/s]---vxLloader(0=RNDM,1=QS)-0.0
97.8
0.0
1
--vx0R[m/s]---vxtR[m/s]--vxcR[m/s]---vxRloader
0.0
97.8
0.0
1
--v0y[m/s]---vty[m/s]---vyloader---v0z[m/s]---vtz[m/s]--vzloader-0.0
97.8
0
0.0
97.8
1
--nbin----Emin[eV]----Emax[ev]---maxnp-100
0.0
100.0
300000
-For-Mid-Diagnostic---nbin----Emin[eV]---Emax[eV]----XStart--XFinish-200
0.0
0.3
0.0
0.03
Some Input Parameters
 nsp : Number of species.
 nc: The number of spatial cells. Δx=length/nc
 nc2p: Superparticle to actual particle weight. The initial
number of superparticles is N=initn·area·length/nc2p.
 dt: Timestep for simulation in seconds.
 length: The length of the system (distance between
electrodes) in meters.
 B: Applied homogeneous magnetic field in Tesla
 PSI: Angle of the B-field in degrees
 extC: The external circuit capacitance in Farads. Used in
conjuction with extL, extR and the driving source.
 source: Either a voltage (v) or current (i) source
 f0: AC frequency of the source.
 GAS: The type of gas, important when collisions are
turned on. Helium = 1, Argon = 2, Neon = 3, Oxygen = 4.
 Gpressure : Background gas pressure in Torr.
 q: Charge of the particle in Coulombs.
 m: Mass of the particle in Kgs.
 initn: Initial particle number density
 For details, refer the source code itself or the manual
inside the package of source file.
Example of Result (driven by RF)
Vx vs. x for electrons
Density vs. x
Ion flux vs. Ion Energy
Vx vs. x for ions
Potential vs. x
Electron Temperature vs. x