ECE490O: Special Topics in EM-Plasma Simulations JK LEE (Spring, 2006)

1
ECE490O Special Topics in EM-Plasma
SimulationsJK LEE (Spring, 2006)
2
ECE490O Special Topics in EM-Plasma
SimulationsJK LEE (Spring, 2006)

• ODE Solvers
• PIC-MCC
• PDE Solvers (FEM and FDM)
• Linear NL Eq. Solvers

3


??- ??   2005-2
????( ??-??-??)    ?? ??? ????? (3-0-3) ????    ???
????-??(?????)     - 01(EECE695) ????-????    ???? - ??
????/???     ?, ? 1615 1730 LG???-0106   ?, ? 1615 1730 LG???-0106   ?, ? 1615 1730 LG???-0106
???    279-2083  279-2083  279-2083























Computational methods for advanced EM including plasmas and charged particles to be described with hands-on experience projects




-ECE261 EM (II) -ECE490D or others



-Mid-Term Quiz -Homeworks -Projects ECE490O Special Topics on Simulation (Spring 2006)




-(A) CK Birdsall and AB Langdon Plasma Phys. via Computer Simul. (Adam Hilger 1991) -(B) S. Nakamura Applied Numerical Methods in FORTRAN (Prentice 1991) or in C (94) -(C) JD Hoffman Numerical Methods for Engineers and Scientists (McGraw Hill 1992)










Week 1-3 ODE solvers with projects in Fortran Week 4-6 ES1 and PIC-MCC in 1d 2d with projects Week 7-9 FEM Week 10-12 PDE solvers via FDM Week 12-15 Nonlinear eq. Solvers and others




ECE490 (Spring 2006) for entry level simulations recommended

4

Week 1-3 ODE solvers with projects in C/Fortran Week 4-6 ES1 and PIC-MCC in 1d 2d with projects Week 7-9 FEM Week 10-12 PDE solvers via FDM Week 13-15 Nonlinear eq. Solvers and others





5
Computational Eng./Sci.
6
  
"Almost everything will soon be simulated before being made" (John McTague of Ford, et al.) Modeling or computer-intensive simulation for EM problems by using PCs, cluster parallel computers, or supercomputers has become an important tool for difficult problem solving. To be taught in English or in Korean or mixture. For juniors, seniors, and grad. students not offered every year (offered only second time in the 19-year history of POSTECH EE Dept.) to be graded separately for grad. and undergrad. students.
7
"??? ??? ??? ? ?? ???? ? ???? " ??? ???? ???
??? ?????
        '??????? ??'? ?? ???(65) ???? ??? 2? ??? ????? ??? ??. ?? ??? ??? ???? 7? ??? ??? ??? ??? ???? ????. ? ??? "?? ??? ???? ??? ?? ??? ?? ??? ???? ??? ???? ??"? "??? ???? ??? ??? ???? ?? ??? ??? ?? ??? ??? ?"??? ???. ? ??? 1968? '????? ??'?? ????? ???? 20?? ??? ? 88? 8.15??? ?????. ???.???? ?? ???? ??? ?? ??? ??? ????? ??? ??? ? ??? ????. ?? "?? ???? ?? ??? ??? ?? ????? ?? ???? ??? ??? ???? ??? ? ??"? "??? ??? ??? ??? ??? ??? ?? ??, ?? ??? ??"? ??. "???? ??? ??? ? ??? ??? ??? ???? ???? ?? ??? ??? ??"? "???? ?? ??? ???? ???? ?? ???? ???? ?? ??? ??? ??? ??? ??"??? ????. 88? ??? '??????? ??'? ? ??? ???? ? ????? ?? 200? ?? ??? ?? ??? ? ??? ??????. ??? ???? "?????? ???? ?????? ?? ??? ?? ??? ????? ??? ?? ? ??? ???? ??"? ????. ?? ? ??? ?????? "??? ??? ??? ??? ??? ??? ?? ???? ? ?? ?????? ???"? "??? ??? ??? ??? ??? ???? ??? ? ?? '???? ????'? ??? ???"? ???. ??? ?? ltblast_at_joongang.co.krgt         
  2006.03.02 2043 ?? / 2006.03.03 0444 ??
8
Gonsalves lecture notes (Fall 2005)
9
Gonsalves lecture notes (Fall 2005)
10
ECE490O ODEJK LEE (Spring, 2006)
11
Gonsalves lecture notes (Fall 2005)
12
(No Transcript)
13
Gonsalves lecture notes (Fall 2005)
14
(No Transcript)
15
(No Transcript)
16
Programs of Initial Value Problem Shooting
Method for BVP ODE
Sung Jin Kim and Jae Koo Lee
17
(No Transcript)
18
(No Transcript)
19
Initial Value Problems of Ordinary Differential
Equations
Program 9-1 Second-order Runge-kutta Method
, y(0)1, y(0)0
Changing 2nd order ODE to 1st order ODE,
y(0)1
(1)
z(0)0
(2)
20
Second-Order Runge-Kutta Method
By 2nd order RK Method
21
Program 9-1
k1 hz l1 -h(bmz kmy) k2 h(z
l1) l2 -h(bm(z l1) km(y k1)) y y
(k1 k2)/2 z z (l1 l2)/2 printf(
" 12.6f 12.5e 12.5e \n", time, y, z )
exit(0)
/ CSL/c9-1.c Second Order Runge-Kutta Scheme
(Solving the problem II of Example 9.6) /
include ltstdio.hgt include ltstdlib.hgt include
ltmath.hgt / time t y,z y,y' kount
number of steps between two lines of printing
k, m, b k, M(mass), B(damping coefficient) in
Example 9.6 int main() int kount, n,
kstep0 float bm, k1, k2, km, l1, l2 static
float time, k 100.0, m 0.5, b 10.0, z
0.0 static float y 1.0, h 0.001 printf(
"CSL/C9-1 Second Order Runge-Kutta Scheme \n" )
printf( " t y z\n" ) printf( " 12.6f 12.5e
12.5e \n", time, y, z ) km k/m bm b/m
for( n 1 n lt 20 n ) for( kount 1
kount lt 50 kount ) kstepkstep1 time
hkstep
2nd order RK
result
CSL/C9-1 Second Order Runge-Kutta Scheme
t y z 0.000000
1.00000e00 0.00000e00 0.100000
5.08312e-01 -6.19085e00 0.200000
6.67480e-02 -2.46111e00 0.300000
-4.22529e-02 -1.40603e-01 0.400000
-2.58300e-02 2.77157e-01 0.500000
-4.55050e-03 1.29208e-01 0.600000
1.68646e-03 1.38587e-02 0.700000
1.28624e-03 -1.19758e-02 0.800000
2.83107e-04 -6.63630e-03 0.900000
-6.15151e-05 -1.01755e-03 1.000000
-6.27664e-05 4.93549e-04
22
Various Numerical Methods
h0.1
Exact Solution
h0.01
h0.001
23
Error Estimation
Error estimation
Fourth order Runge-Kutta
24
Program 9-2 Fourth-order Runge-Kutta Scheme
A first order Ordinary differential equation
y(0)1
Fourth-order RK Method
25
Boundary Value Problems of Ordinary Differential
Equations Scharfetter-Gummel method
Sung Soo Yang and Jae Koo Lee
Homepage http//jkl.postech.ac.kr
26
Boundary value problems
Type of Problems Advantages Disadvantages
Shooting method An existing program for initial value problems may be used Trial-and-error basis. Application is limited to a narrow class of problems. Solution may become unstable.
Finite difference method using the tridiagonal solution No instability problem. No trial and error. Applicable to nonlinear problems with iteration. Problem may have to be developed for each particular problem.
27
Program 10-1
Solve difference equation,
With the boundary conditions,
x 0
1
2
9
10
i 0
1
2
9
10
known y(0)1
Especially for i 1,
28
Program 10-1
For i 10,
Summarizing the difference equations obtained, we
write
Tridiagonal matrix
29
Solution Algorithm for Tridiagonal Equations (1)
R2
Based on Gauss elimination
R3
30
Solution Algorithm for Tridiagonal Equations (2)
31
Solution Algorithm for Tridiagonal Equations (3)
void trdg(float a, float b, float c, float
d, int n) / Tridiagonal solution /
int i float r for ( i 2 i lt n
i ) r ai/bi - 1
bi bi - rci - 1 di di -
rdi - 1 dn dn/bn
for ( i n - 1 i gt 1 i-- )
di (di - cidi 1)/bi
return
Recurrently calculate the equations in increasing
order of i until iN is reached
Calculate the solution for the last unknown by
Calculate the following equation in decreasing
order of i
32
Program 10-1
/ CSL/c10-1.c Linear Boundary Value Problem /
include ltstdio.hgt include ltstdlib.hgt
include ltmath.hgt / ai, bi, ci, di
a(i), b(i), c(i), and d(i) n number of grid
points / int main() int i, n
float a20, b20, c20, d20, x void
trdg(float a, float b, float c, float d,
int n) / Tridiagonal solution / printf(
"\n\nCSL/C10-1 Linear Boundary Value Problem\n"
) n 10 / n Number of grid points /
for( i 1 i lt n i ) x
i ai -2 bi 5 ci -2
di exp( -0.2x )
33
Program 10-1
d1 d1 2 dn dn0.5
bn 4.5 trdg( a, b, c, d, n ) d0
1 / Setting d0 for printing purpose /
printf( "\n Grid point no. Solution\n" )
for ( i 0 i lt n i ) printf(
" 3.1d 12.5e \n", i, di )
exit(0)
CSL/C10-1 Linear Boundary Value Problem Grid
point no. Solution 0 1.00000e00
1 8.46449e-01 2 7.06756e-01
3 5.85282e-01 4 4.82043e-01
5 3.95162e-01 6 3.21921e-01
7 2.59044e-01 8 2.02390e-01
9 1.45983e-01 10 7.99187e-02
34
Program 10-3
An eigenvalue problem of ordinary differential
equation
We assume
35
ECE490O PIC FEMJK LEE (Spring, 2006)
36
(No Transcript)
37
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 Maxwells 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.

38
Capacitively Coupled Plasma 1D PIC-MCC

j 1, ? , N


Sheath


Bulk Plasma
Sheath
Substrate

• MCC (Monte-Carlo Collision) Processes
• - Electron-Neutral Collisions
• (Ionization, Scattering, Excitation)
• - Ion-Neutral Collisions
• (Charge-exchange, Scattering)

• 1D Asymmetric Dual-Freq. Voltage-Driven System

39
PIC-MCC Flow Chart

• Particles in continuum space
• Fields at discrete mesh locations in space
• Coupling between particles and fields

I
II
V
IV
III
IV
Fig Flow chart for an explicit PIC-MCC scheme
40
I. Particle Equations of Motion

• Newton-Lorentz equations of motion

• In finite difference form, the leapfrog method

Fig Schematic leapfrog integration
41
III. Electrostatic Field Model

• Possions equation

• Finite difference form in 1D planar geometry

• Boundary condition External circuit

Fig Schematic one-dimensional bounded
plasma with external circuit
42
PDP Structure
AC PDP
Discharge in PDP
43
Striation Profiles in PDP 2D PIC/MCC
Anode
Cathode

• Pressure dependence of striations
• Number of peaks depend on the pressure and
  electrode size.

44
XOOPIC and MAGIC Codes for Electromagnetic Field
S.J. Kim and J.K. Lee
Contents

• Overview of XOOPIC code
• Overview of MAGIC code
• Klystron simulation using XOOPIC code

45
Simulation Domain of Klystron
RF output port
RF input port
9.55 cm
10.05 cm
13.07 cm
E-beam
7.569 cm
6.66 cm
Cylindrical Axis
37.2 cm

• Simulation condition

• Beam emitter I 12 kA, ud 2.48e8 m/s
• Input port Rin2300 ?, R20 ?, f7.69 GHz
• Output port R47.124 ?

46
Example of Klystron Simulation
Phase space
Density
Kinetic energy
uz
47
Simulation Results at 0.5 ns
48
Simulation Results at 10 ns
49
Simulation Results at 20 ns
50
KE as a Function of Beam Current
51
Finite Difference Method (II)
Grid 100x50
52
Full rectangle
100V
-100V
100V
-100V
53
Full circle
100V
-100V
100V
-100V
54
What is FEMLAB?

• FEMLAB a powerful interactive environment for
  modeling and
• solving various kinds of scientific and
  engineering problems based
• on partial differential equations (PDEs).

• Overview

• Finite element method
• GUI based on Java
• Unique environments for modeling
• (CAD, Physics, Mesh, Solver, Postprocessing)
• Modeling based on equations (broad application)
• Predefined equations and User-defined
  equations
• No limitation in Multiphysics
• MATLAB interface (Simulink)

• Mathematical application modes and types of
  analysis

• Mathematical application modes
• 1. Coefficient form suitable for linear or
  nearly linear models.
• 2. General form suitable for nonlinear models
• 3. Weak form suitable for models with PDEs on
  boundaries, edges,
• and points, or for models using terms with
  mixed space and time
• derivatives.
• Various types of analysis
• 1. Eigenfrequency and modal analysis
• 2. Stationary and time-dependent analysis
• 3. Linear and nonlinear analysis

Reference Manual of FEMLAB Software
55
Useful Modules in FEMLAB

• Application areas

• Microwave engineering
• Optics
• Photonics
• Porous media flow
• Quantum mechanics
• Radio-frequency components
• Semiconductor devices
• Structural mechanics
• Transport phenomena
• Wave propagation

• Acoustics
• Bioscience
• Chemical reactions
• Diffusion
• Electromagnetics
• Fluid dynamics
• Fuel cells and electrochemistry
• Geophysics
• Heat transfer
• MEMS

• Additional Modules

1. Application of Chemical engineering Module

• Momentum balances
• - Incompressible Navier-Stokes eqs.
• - Darys law
• - Brinkman eqs.
• - Non-Newtonian flow

• Mass balances
• - Diffusion
• - Convection and Conduction
• - Electrokinetic flow
• - Maxwell-stefan diffusion and convection

• Energy balances
• - Heat equation
• - Heat convection and conduction

2. Application of Electromagnetics Module

• - Electrostatics
• - Conductive media DC
• Magnetostatic
• Low-frequency electromagnetics
• - In-plane wave propagation
• Axisymmetric wave propagation
• Full 3D vector wave propagation
• Full vector mode analysis in 2D and 3D

3. Application of the Structural Mechanics Module

• Plane stress
• Plane strain
• 2D, 3D beams, Euler theory
• Shells

56
FEMLAB Environment
Model Navigator
Pre-defined Equations
57
Monoconical RF Antenna

• Introduction of Conical RF Antenna

• Antennas are used to couple guided
  electromagnetic energy to waves
• in free space.
• Antenna engineering design involves minimizing
  reflection losses and
• obtaining desired directional properties for a
  wide range of frequencies.
• - antenna impedance, radiation pattern,
  frequency analysis.
• Modeling can be extremely useful as it reduces
  the need for building
• prototypes and performing measurements.

• Antenna parameters

Antenna
metal
teflon ?r2.07
Ground plane
Coaxial feed
Coaxial feed rinner1.5mm router4.916mm Z50
58
ECE490O PDE Gonsalves lecture notes (Fall
2005) JK LEE (Spring, 2006)
59
(No Transcript)
60
(No Transcript)
61
(No Transcript)
62
Scharfetter-Gummel method
Tridiagonal matrix

• 2D discretized continuity eqn. integrated by the
  alternative direction implicit (ADI) method

Scharfetter-Gummel method
63
Plasma Display Panel
Plasma Display Panel
Many Pixels
the flat panel display using phosphor
luminescence by UV photons produced in plasma gas
discharge
PDP structure
64
Simulation domain
Sustain 1
Sustain 2
dielectric layer
dielectric and phosphor layer
address
Finite-Difference Method
?
?
?
j
? Electric field, Density
?
?
?
? Potential, Charge
?
?
?
? Flux of x and y
j1
i
i1
? Light, Luminance, Efficiency, Power, Current
and so on
65
Flow chart
fl2p.c
initial.c
pulse.c
charge.c
field.c
flux.c
continuity.c
source.c
history.c
diagnostics.c
Calculate efficiency
time_step.c
current.c, radiation.c, dump.c, gaspar.c,
mu_n_D.c, gummel.c
66
Basic equations

• Continuity Equation with Drift-Diffusion Approx.

and

• Poissons Equation

? surface charge density on the dielectric
surfaces

• Boundary conditions on dielectric surface

for ion
for electron
Mobility-driven drift term
Isotropic thermal flux term
for secondary electron
for excited species
67
Partial Differential Eqs.
General form of linear second-order PDEs with
two independent variables
In case of elliptic PDEs,
?
Jacobi-Iteration method
?
Gauss-Seidel method
?
Successive over-relaxation (SOR) method
In case of parabolic PDEs,
?
Alternating direction implicit (ADI) method
68
Continuity equation (1)
density nsp
Alternating direction implicit (ADI) method
ADI method uses two time steps in two dimension
to update the quantities between t and t?t.
During first ?t/2, the integration sweeps along
one direction (x direction) and the other
direction (y direction) is fixed. The temporary
quantities are updated at t?t/2. With these
updated quantities, ADI method integrates the
continuity equation along y direction with fixed
x direction between t?t/2 and t?t.
1st step
(k means the value at time t) ( means the
temporal value at time t?t/2 )
Discretized flux can be obtained by
Sharfetter-Gummel method.
Spatially discretized forms are converted to
tridiagonal systems of equations which can be
easily solved.
69
Tridiagonal matrix (1)
Based on Gauss elimination
R2
R3
70
Tridiagonal matrix (2)
71
Tridiagonal matrix (3)
/ Tridiagonal solution / void trdg(float a,
float b, float c, float d, int n)
int i float r for ( i 2 i lt n
i ) r ai/bi - 1
bi bi - rci - 1 di di -
rdi - 1 dn dn/bn
for ( i n - 1 i gt 1 i-- )
di (di - cidi 1)/bi
return
Calculate the equations in increasing order of i
until iN is reached.
Ri
Calculate the solution for the last unknown by
Calculate the following equation in decreasing
order of i
72
Continuity equation (2)
2nd step
From the temporally updated density calculated in
the 1st step, we can calculated flux in
x-direction (?) at time t?t/2. Using these
values, we calculate final updated density with
integration of continuity equation in y-direction.
(k1 means the final value at time t?t) (
means the temporal value at time t?t/2 )
(tridiagonal matrix)
From the final updated density calculated in the
2nd step, we can calculated flux in y-direction
(?k1) at time t.
73
Poissons eq. (1)
Poisson equation is solved with a successive over
relation (SOR) method. The electric field is
taken at time t when the continuity equations are
integrated between t and t?t. Time is integrated
by semi-implicit method in our code. The electric
field in the integration of the continuity
equation between t and t?t is not the field at
time t, but rather a prediction of the electric
field at time t?t. The semi-implicit integration
of Poisson equation is followed as
Density correction by electric field change
between t and t?t
The continuity eq. and flux are coupled with
Poissions eq.
This Poissons eq can be discriminated to x and y
directions, and written in matrix form using the
five-point formula in two dimensions.
74
Poissons eq. (2)
j1
ci, j
ai, j
bi, j
j
is the surface charge density
accumulating on intersection between plasma
region and dielectric.
di, j
j-1
i-1
i
i1
Solved using SOR method
75
ECE490O NL Eq. SolversGonsalves lecture notes
(Fall 2005)JK LEE (Spring, 2006)
76
(No Transcript)
77
(No Transcript)
78
(No Transcript)
79
(No Transcript)
80
(No Transcript)
81
Mullers method for solving non-linear equations
FORTRAN features
N. Babaeva
82
Muller_4.c
R0.01
The beam is stronger
Real and imaginary roots of the dispersion
relation
R0.01
83
W.H.Press et al Numerical recipes in C The art
of Scientific Computing
Mullers method
Joe D. Hoffman Numerical Methods for Engineers
and scientists McGraw-Hill, Inc. 1993
In Mullers method, three approximations to the
zero are required. The next approximation is the
zero of the parabola that passes through the
three points. Both real zeros and complex zeros
can be found.Mullers method generalizes the
secant method, but uses quadratic interpolation
among three points instead of linear
interpolation between two. Solving for the zeros
of the quadratic allows the method to find
complex pairs of roots. Given three previous
guesses for the root
And the values of polynomial P(x) at those
points, the next approximation is produces by
the following formulas
Figure, illustrating the Mullers method
Mullers method is also used for finding complex
zeros of analytic functions (not just
polynomials) in the complex plane.