Homepage : jkl.postech.ac.kr

Plasma Application Modeling Boundary Value Problems of Ordinary Differential Equations & Scharfetter-Gummel method (Computational E&M: 490D) Sung Soo Yang, and Jae Koo Lee March 23, 2004 Homepage : http://jkl.postech.ac.kr

* 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

         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,

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

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

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

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

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 ); }

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

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

Plasma Application Modeling@ POSTECH Program 10-3 (Inverse Power method) Step1 : for all i and are set to an arbitrary initial guess. is solved by the tridiagonal solution for Step2 : Step3 : The next estimate for is calculated by the equation: is solved by the tridiagonal solution for Step4 : Step5 : The operations similar to Step 3 and 4 are repeated as the iteration cycle t increases. Step6 : The iteration is stopped when the convergence test is satisfied. Criterion for convergence

Plasma Application Modeling@ POSTECH Program 10-3 while (TRUE) { k = k+1; for (i=1; i<=n; i++) fb[i] = f[i]*ei; for (i=1; i<=n; i++) { a[i] = as[i]; b[i] = bs[i]; c[i] = cs[i]; d[i] = ds[i]*fb[i]; } trdg(a, b, c, d, n); sb = 0; s = 0; for (i=1; i <= n; i++) { f[i] = d[i]; s = s + f[i]*f[i]; sb = sb + f[i]*fb[i]; } eb = ei; ei = sb/s; printf("%3.1d %12.6e \n", k, ei); if (fabs(1.0-ei/eb) <= ep) break; if (k >it) { printf("Iteration limit exceeded.\n"); break; } } main() { • • • • • while (TRUE) { k = 0; n = 10; ei = 1.5; it = 30; ep = 0.0001; for (i=1; i <= n; i++) { as[i] = -1.0; bs[i] = 2.0; cs[i] = -1.0; f[i] = 1.0; ds[i] = 1.0; } printf("\n It. No. Eigenvalue\n"); Step 2,4 Step 1 Step 3 Step 6

Plasma Application Modeling@ POSTECH Program 10-3 It. No. Eigenvalue 1 8.196722e-02 2 8.102574e-02 3 8.101422e-02 4 8.101406e-02 Eigenvalue = 0.0810141 Eigenfunction i f(i) 1 2.84692e-01 2 5.46290e-01 3 7.63595e-01 4 9.19018e-01 5 1.00000e+00 6 1.00000e+00 7 9.19018e-01 8 7.63594e-01 9 5.46290e-01 10 2.84692e-01 ------------------------------------ Type 1 to continue, or0 to stop. z = 0; for (i=1; i <= n; i++) if (fabs(z) <= fabs(f[i])) z = f[i]; for (i=1; i <= n; i++) f[i] = f[i]/z; eigen = ei; printf("Eigenvalue = %g \n", eigen); printf("\n Eigenfunction\n"); printf("i f(i)\n"); for (i=1; i<=n; i++) printf("%3.1d %12.5e \n", i, f[i]); printf("------------------------------------\n"); printf("Type 1 to continue, or0 to stop. \n"); scanf("%d", &kstop); if (kstop != 1) exit (0); } }

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

Plasma Application Modeling@ POSTECH Scharfetter-Gummel method Flux is constant between half grid points and calculated at the grid point   Flux is a linear combination of and Analytic integration between i and i+1 leads to

Plasma Application Modeling@ POSTECH Scharfetter-Gummel method where, The main advantage of SG scheme is that is provides numerically stable estimates of the particle flux under all conditions (Big potential difference) (small potential difference)

Plasma Application Modeling@ POSTECH Scharfetter-Gummel method In case of Drift flux

Plasma Application Modeling@ POSTECH Scharfetter-Gummel method In case of Diffusion flux

Homepage : jkl.postech.ac.kr