Transcript Computational E&M: 490D

ECE490O: Special Topics
in EM-Plasma Simulations
JK LEE (Spring, 2006)
ECE490O: Special Topics
in EM-Plasma Simulations
JK LEE (Spring, 2006)
•
•
•
•
ODE Solvers
PIC-MCC
PDE Solvers (FEM and FDM)
Linear & NL Eq. Solvers
ECE490O: PDE
Gonsalves’ lecture notes (Fall 2005)
JK LEE (Spring, 2006)
Plasma Display Panel
Plasma Display Panel
Many Pixels
the flat panel display using phosphor luminescence
by UV photons produced in plasma gas discharge
Discharge
PDP structure
bus electrode
dielectric
MgO layer
barrier
address
electrode
White light emission
Discharge
Discharge
Front panel
ITO electrode
phosphors
Back panel
(1)
Electric input power
(2)
Discharge
(3)
VUV radiation
(4)
Phosphor excitation
(5)
Visible light in cell
(6)
DisplayPlasma
light Application
Modeling @ POSTECH
Simulation domain
y
Sustain 1
Sustain 2
dielectric layer
ny
dielectric and phosphor layer
x
address
nx



i






j
Finite-Difference Method
 Electric field, Density
 Potential, Charge
j+1
 Flux of x and y
i+1
 Light, Luminance, Efficiency, Power, Current and so on
Plasma Application
Modeling @ POSTECH
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
Plasma Application
Modeling @ POSTECH
Basic equations
• Continuity Equation with Drift-Diffusion Approx.
nsp
t
Γ p   D p n p   p n p E
   Γ sp  S sp
Γ e   Dene  e ne E
and
Γ ex   Dex nex
• Poisson’s Equation
  ( V )    ql nl  
 : surface charge density on the dielectric surfaces
• Boundary conditions on dielectric surface
Γi  n  i ni E  n  0.25 ni vi
Γe  n  e ne E  n  0.25 ne ve  Γ se  n
Mobility-driven drift
term
Γ se  - se,iΓ i
Isotropic thermal flux
term
for ion
for electron
 th 
8k B T
m
for secondary electron
i
Γex n  0.25 nex vex
for excited species
Plasma Application
Modeling @ POSTECH
Partial Differential Eqs.
General form of linear second-order PDEs
with two independent variables
auxx  buxy  cuyy  dux  euy  fu  g  0
b 2  4ac  0, Hyperbolic(ex. Waveeq.)
 2
b  4ac  0, Parabolic (ex.Continuityeq.)
b 2  4ac  0, Elliptic (ex.Poisson's eq.)

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
Plasma Application
Modeling @ POSTECH
Continuity equation (1)
nsp
t
   Γ sp  S sp
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
ni*, j  nik, j
t / 2

x*,i 0.5, j  x*,i 0.5, j
x

Sik, j

yk,i , j 0.5  yk,i , j 0.5
y
(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.
x,i0.5,j   x,i,j ni , j   x,i , j ni 1, j
Spatially discretized forms are converted to tridiagonal systems of equations which can
be easily solved.
*
*
*
Ai, j ni 1, j  Bi, j ni, j  Ci, j ni 1, j  Di, j
Plasma Application
Modeling @ POSTECH
Tridiagonal matrix (1)
Based on Gauss elimination
 B1
A
 2
 0
R2
 A2
 B B1
 1
 A2
 0


 A2
A
 2
 0
C1
B2
A3
0  1   D1 
C2  2    D2 
B3  3   D3 
 A2
0

 0
 A2
0

 0
B2
A3
0  1   R2 D1 
C2  2    D2 
B3  3   D3 
R2C1
B2
B2  R2C1
A3
0  1   R2 D1 
C2  2    D2  R2 D1 

B3  3  
D3
A3
0
 A2
0

 0
A2

 A2 
C1 0  
D
 1   B 1
B1

 1 
B2
C2  2    D2 
A3
B3  3   D3 






R2C1
R2C1
 A2

0

 0
R3
R2C1
A3
B2
B2
A3
0  1   R2 D1 
R3C2  2    R3 D2 
B3  3   D3 
 1 
 
A3 
B3 R3C2  2  
0
B3  R3C2  3 
 R2 D1 
 R D 
 D3 3 2 
 D3  R3 D2 
R2C1
0
0     R2 D1 
1
A3     A3 
C2  2   D2 
B2     B2 
 
B3  3   D3 
D2
Plasma Application
Modeling @ POSTECH
Tridiagonal matrix (2)
 A2
0

 0
R2C1
A3
0
0  1   R2 D1 
R3C2  2    R3 D2 
B3  3   D3 
A32  R3C23  R3 D2
 2 
 2 
D3
B3
R3
A3
( D2  C23 ) , R3 
A3
B2
1
( D2  C23 )
B2
A21  R2C12  R2 D1  1 
 1 
 3 
R2
A
( D1  C12 ) , R2  2
A2
B1
1
1
( D1  C12 ) 
( D1  C12 )
B1
B1
 i  ( Di  Cii 1 ) Bi
Plasma Application
Modeling @ POSTECH
Tridiagonal matrix (3)
Ai, j ni*1, j  Bi, j ni*, j  Ci, j ni*1, j  Di, j
/* Tridiagonal solution */
void trdg(float a[], float b[], float c[], float d[], int n)
{
int i;
float r;
for ( i = 2; i <= n; i++ )
{
r = a[i]/b[i - 1];
b[i] = b[i] - r*c[i - 1];
d[i] = d[i] - r*d[i - 1];
}
d[n] = d[n]/b[n];
Ri
Bi
Di
Calculate the equations in increasing
order of i until i=N is reached.
Calculate the solution for the last
unknown by
 BN

N  DN
Calculate the following equation in
decreasing order of i
 i  ( Di  Cii 1 ) Bi
for ( i = n - 1; i >= 1; i-- )
{
d[i] = (d[i] - c[i]*d[i + 1])/b[i];
Plasma Application
Modeling @ POSTECH
Continuity equation (2)
Ai, j ni*1, j  Bi, j ni*, j  Ci, j ni*1, j  Di, j
z
t
Ai , j  
2x 2
t
Ci , j  
2x 2
z
1
i , j
2
z 1
e
1
i , j
2
z 1
e
i
2
i
2
Di 1, j
z
e
, j
i
1
, j
2
Bi ,
1
j
t
1 
Di ,
2x 2
Di 1, j
, j
Di , j  nik, j 
1
j
 z
1
z
i
, j
1

i 
,
2
2
e
 z
1
 e i 2 , j 1

z
j

z
e
i
i 



1

1
, j
2
1
, j
2

t k
t  k
Si , j 
 1  k 1 
i, j  
2
y  i , j  2
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.
1
nik, 
 ni*, j
j
t / 2

1
k 1
yk,
i , j  0.5  y ,i , j 0.5
y

Sik, j

x*,i 0.5, j  x*,i 0.5, j
x
(k+1 means the final value at time t+t)
( * means the temporal value at time t+t/2 )
1
'
k 1
1
'
Ai', j nik, 
 Ci', j nik, 
j 1  Bi , j ni , j
j 1  Di , j (tridiagonal matrix)
From the final updated density calculated in the 2nd step, we can calculated flux in ydirection (k+1) at time t.
Plasma Application
Modeling @ POSTECH
Poisson’s eq. (1)
  ( 0 V )   e n p  ne 
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
  ( V
k 1
)  
 n p n p 
e  k
k
n

n


t

p
e
0 
t




Density correction by electric field change
between t and t+t
The continuity eq. and flux are coupled with Poission’s eq.


e
e
k k 
k 1


    
t  nl  l  V   
0
0


l

 (sign) nlk
l
This Poisson’s eq can be discriminated to x and y directions, and written in matrix form
using the five-point formula in two dimensions.
ai , jVi 1,
j
 bi, jVi1, j  ci , jVi , j1  d i , jVi , j1  ei , jVi , j  f i , j
Plasma Application
Modeling @ POSTECH
Poisson’s eq. (2)
ai , jVi 1,
ai ,

j
bi , j 
ci ,
di , j 
ei ,
j
 bi, jVi1, j  ci , jVi , j1  d i , jVi , j1  ei , jVi , j  f i , j

1
 i,
2x 2 

j
 i,
j 1


e
0

t  ni , j  i , j  ni ,
j
i ,


j 1

l

1
 i,
2y 2 

j
  i 1,
j


e
0

t  ni , j  i ,
j
 ni 1,
j
 i  1,
j


l



1
e

 i , j 1   i 1, j1  
t  ni , j 1i , j 1  ni 1, j 1 i 1, j 1 
2 
2y 
0
l

  (ai ,
 bi, j  ci ,
j
 di, j )
j
j+1
f i,
j 1


1
e

 i  1, j  i  1, j 1 
t  ni 1, j i 1, j  ni 1, j 1 i 1, j 1 
2 
2x 
0
l


j
j
 
N
q
l
l 1
 (
d iek .
)
0
j
 diek
is the surface charge density
accumulating on intersection between
plasma region and dielectric.
ci, j
bi, j
ai, j
di, j
j-1
i-1
i
i+1
Solved using SOR method
Plasma Application
Modeling @ POSTECH
Scharfetter-Gummel method
n
  Γ  S
t
Γ   Dn  sgn(q)nE
Ai, j ni*1, j  Bi, j ni*, j  Ci, j ni*1, j  Di, j
Ai', j nik, j 11  Bi', j nik, j 1  Ci', j nik, j 11

Di', j
Tridiagonal matrix
• 2D discretized continuity eqn. integrated by the alternative direction implicit
(ADI) method
ni*, j  nik, j
t / 2
nik, j 1  ni*, j
t / 2


x*,i 0.5, j  x*,i 0.5, j
x
yk,i ,1j 0.5  yk,i,1j 0.5
y
 Sik, j 
 Sik, j 
yk,i, j 0.5  yk,i, j 0.5
y
x*,i 0.5, j  x*,i 0.5, j
x
x,i0.5,j   x,i,j ni, j   x,i, j ni1, j Scharfetter-Gummel method
Plasma Application
Modeling @ POSTECH
Gonsalves’ lecture notes (Fall 2005)
Gonsalves’ lecture notes (Fall 2005)
Gonsalves’ lecture notes (Fall 2005)
Gonsalves’ lecture notes (Fall 2005)
Gonsalves’ lecture notes (Fall 2005)
ECE490O: NL Eq. Solvers
Gonsalves’ lecture notes (Fall 2005)
JK LEE (Spring, 2006)