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

5. Computational Eng./Sci.

7. ECE490O: ODEJK LEE (Spring, 2006)

8. Plasma Application Modeling Group POSTECH Programs of Initial Value Problem & Shooting Method for BVP ODE Sung Jin Kim and Jae Koo Lee

11. Plasma Application Modeling, POSTECH 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)

12. Plasma Application Modeling, POSTECH Second-Order Runge-Kutta Method By 2nd order RK Method

13. Plasma Application Modeling, POSTECH Program 9-1 k1 = h*z; l1 = -h*(bm*z + km*y); 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 <stdio.h> #include <stdlib.h> #include <math.h> /* 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, kstep=0; 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 <= 20; n++ ){ for( kount = 1; kount <= 50; kount++ ){ kstep=kstep+1; time = h*kstep ; 2nd order RK result CSL/C9-1 Second Order Runge-Kutta Scheme t y z 0.000000 1.00000e+00 0.00000e+00 0.100000 5.08312e-01 -6.19085e+00 0.200000 6.67480e-02 -2.46111e+00 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

14. Plasma Application Modeling, POSTECH Various Numerical Methods h=0.1 Exact Solution: h=0.01 h=0.001

15. Plasma Application Modeling Group POSTECH Error Estimation Error estimation Fourth order Runge-Kutta

16. Plasma Application Modeling, POSTECH Program 9-2 Fourth-order Runge-Kutta Scheme A first order Ordinary differential equation y(0)=1 Fourth-order RK Method

17. Plasma Application Modeling Boundary Value Problems of Ordinary Differential Equations & Scharfetter-Gummel method Sung Soo Yang and Jae Koo Lee Homepage : http://jkl.postech.ac.kr

18. * Three types of boundary conditions 1) : Dirichlet type boundary condition Plasma Application Modeling@ POSTECH 1) : Neumann type boundary condition 1) : Mixed type boundary condition Boundary value problems

19.          Plasma Application Modeling@ POSTECH 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,

20. Plasma Application Modeling@ POSTECH Program 10-1 For i = 10, Summarizing the difference equations obtained, we write Tridiagonal matrix

21. Plasma Application Modeling@ POSTECH Solution Algorithm for Tridiagonal Equations (1) R2 Based on Gauss elimination R3

22. Plasma Application Modeling@ POSTECH Solution Algorithm for Tridiagonal Equations (2)

23. Plasma Application Modeling@ POSTECH 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 <= 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]; for ( i = n - 1; i >= 1; i-- ) { d[i] = (d[i] - c[i]*d[i + 1])/b[i]; } return; } Recurrently calculate the equations in increasing order of i until i=N is reached Calculate the solution for the last unknown by Calculate the following equation in decreasing order of i

24. Plasma Application Modeling@ POSTECH Program 10-1 /* CSL/c10-1.c Linear Boundary Value Problem */ #include <stdio.h> #include <stdlib.h> #include <math.h> /* a[i], b[i], c[i], d[i] : a(i), b(i), c(i), and d(i) n: number of grid points */ int main() { int i, n; float a[20], b[20], c[20], d[20], 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 <= n; i++ ) { x = i; a[i] = -2; b[i] = 5; c[i] = -2; d[i] = exp( -0.2*x ); }

25. Plasma Application Modeling@ POSTECH Program 10-1 d[1] = d[1] + 2; d[n] = d[n]*0.5; b[n] = 4.5; trdg( a, b, c, d, n ); d[0] = 1; /* Setting d[0] for printing purpose */ printf( "\n Grid point no. Solution\n" ); for ( i = 0; i <= n; i++ ) { printf( " %3.1d %12.5e \n", i, d[i] ); } exit(0); } CSL/C10-1 Linear Boundary Value Problem Grid point no. Solution 0 1.00000e+00 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

26. a b    i-1 i i+1 Plasma Application Modeling@ POSTECH Program 10-3 An eigenvalue problem of ordinary differential equation We assume

27. Plasma Application Modeling@ POSTECH Scharfetter-Gummel method Tridiagonal matrix • 2D discretized continuity eqn. integrated by the alternative direction implicit (ADI) method Scharfetter-Gummel method

28. Muller’s method & FORTRAN features Muller’s method for solving non-linear equations & FORTRAN features N. Babaeva

29. Muller’s method & FORTRAN features Muller_4.c R=0.01 The beam is stronger Real and imaginary roots of the dispersion relation R=0.01

30. W.H.Press et al “Numerical recipes in C” The art of Scientific Computing Muller’s method & FORTRAN features Muller’s method Joe D. Hoffman “Numerical Methods for Engineers and scientists” McGraw-Hill, Inc. 1993 In Muller’s 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 Muller’s method Muller’s method is also used for finding complex zeros of analytic functions (not just polynomials) in the complex plane.

31. Electric input power • Discharge • VUV radiation • Phosphor excitation • Visible light in cell • Display light White light emission Discharge Discharge Discharge Front panel bus electrode dielectric Plasma Application Modeling@ POSTECH MgO layer ITO electrode phosphors barrier address electrode Back panel 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

32. y x ny Plasma Application Modeling@ POSTECH nx 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 j+1 i i+1  Light, Luminance, Efficiency, Power, Current and so on

33. 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

34. Plasma Application Modeling@ POSTECH Basic equations • Continuity Equation with Drift-Diffusion Approx. and • Poisson’s 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

35. Plasma Application Modeling@ POSTECH 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

36. Plasma Application Modeling@ POSTECH 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.

37. Plasma Application Modeling@ POSTECH Tridiagonal matrix (1) Based on Gauss elimination R2 R3

38. Plasma Application Modeling@ POSTECH Tridiagonal matrix (2)

39. Plasma Application Modeling@ POSTECH 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 <= 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]; for ( i = n - 1; i >= 1; i-- ) { d[i] = (d[i] - c[i]*d[i + 1])/b[i]; } return; } Calculate the equations in increasing order of i until i=N is reached. Ri Calculate the solution for the last unknown by Calculate the following equation in decreasing order of i

40. Plasma Application Modeling@ POSTECH 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. (k+1 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 (k+1) at time t.

41. Plasma Application Modeling@ POSTECH Poisson’s 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 Poission’s eq. 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.

42. Plasma Application Modeling@ POSTECH Poisson’s eq. (2) j+1 ci, j bi, j ai, j j is the surface charge density accumulating on intersection between plasma region and dielectric. di, j j-1 i-1 i i+1 Solved using SOR method

43. Finite Difference Method (II) Grid 100x50

44. Full rectangle 100V -100V 100V -100V

45. Full circle 100V -100V 100V -100V

46. Plasma Application Modeling, POSTECH 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

47. 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. • - Dary’s 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

48. Plasma Application Modeling, POSTECH FEMLAB Environment Model Navigator Pre-defined Equations

49. Plasma Application Modeling, POSTECH 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 r=2.07 Ground plane Coaxial feed Coaxial feed: rinner=1.5mm router=4.916mm Z=50

50. 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.