2. The Plasma Fluid Modeling 11
2.3 Numerical Methods and Algorithms
2.3.1 Sharfetter-Gummel Scheme for Mass Flux
It is well-known that it required very fine grid sizes for the convection-diffusion type PDE to be numerically stable if standard finite-difference method is employed. The Scharfetter-Gummel (SG) scheme, which originally had been developed for semiconductor device simula-tions, provides an optimum way to discretize the drift-diffusion equation for particle transport [Scharfetter and Gummel, 1969]. The idea of SG scheme is to approximate the flux by a bound-ary value problem along each edge in a given grid, which yields an exponential approximation to the potential distribution. Therefore, it’s also called an exponentially fitted method.
To demonstrate how the SG scheme works, we derive the corresponding formulation using a one-dimensional example. Assume there are two grid points i and i + 1, and the flux between the grids is Γi+1
2. The drift-diffusion equation of positive charged particle can be written as Γp,i+1
where the first and second terms in the RHS represent the diffusion and drift fluxes, respectively.
By the change of variables from position to potential through chain rule, we can obtain
∂φ
in which assume that ∂φ∂x is independent of potential. After rearrangement, we can have the first order ordinary differential equation as follows:
∂
The solution of above ODE is readily obtained as
2 can be found easily as,
C = ni− ni+1
Similarly we can write the electron flux equation as Γe = − 1
The corresponding ODE for electrons can be written as kBTe Again, for one-dimensional case, we can obtain the solution of the above ODE for electrons as
ne(φ) = Ce−
where C is a constant. a and b are defined, respectively, as a = ∇ (kBTe)
kBTe + q kBTe
E,~ (2.78)
b = − Γe
kBTe
meν
. (2.79)
The constant C and the flux Γecan be readily obtained from the boundary conditions. They are written as
. Or equivalently, the electron flux can be expressed as
Γe,i+1
2 =−De,i+1
2
∆x [ni+1B(X)− niB(−X)] (2.82) 2.3.2 Non-dimensionalization
Let n0 denote the background gas density. The dimensionless potential, length, time, and velocity are defined as follows:
where u0is defined as
u0 =
The dimensionless electron energy density is defined as
0 = e
kBT0 (2.88)
, which the is in unit of eV .
x0 = x
λ (2.89)
, where the characteristic length is the mean free path λ = 1
n0σp. (2.90)
The characteristic time scale is thus defined as t0 = u0
λ t. (2.91)
The dimensionless diffusion and mobility coefficient are, respectively, defined as
The ion flux, eq (2.65) , with drift-diffusion approximation can thus be non-dimensionalized as
Similarly, the electron flux, eq (2.74), can be non-dimensionalized as Γ0e,i+1
With the above dimensionless expression for the flux, the general continuity equation can be non-dimensionalized as The dimensionless electron energy density equation can e written as
u0
where the dimensionless electron energy density source term is written as S0 =−0
eSe0 − ΓeE (2.101)
Finally, the dimensionless Poisson’s equation is expressed as 1
or
where the dimensionless electron energy density is expressed as
0 = 0u20mp λ2e2n0
(2.105)
2.3.3 Finite-difference Discretization
Finite differences method (FDM) is widely used in the field of Computational Fluid Dy-namics (CFD), which is based on results of Taylor series expansion. One of the advantages of FDM is that it could be implemented easily, especially in the case of simple geometry. The error caused by the discretization process depends on the number of terms in the Taylor series which are kept. Generally, we can solve a PDE or a set of PDEs using the finite-difference method with dependent variables as functions of spatial coordinates and time in implicit, explicit, or semi-implicit approaches. The explicit approaches are relatively simple to implement; however, it is often suffered by the stability constraint, which is known as the Courant-Friedrichs-Lewy (CFL) condition, resulting in very small computational time step. In contrast, the implicit approaches are stable even for larger values of the time step. However, a system of algebraic equations must be solved at each time step using implicit techniques which is more complicated to implement and possibly more time-consuming. In this thesis, we descritize the fluid modeling equations through the use of finite-difference method. Since continuity equations for different species are similar, we only present the discretizations for the eqs. (2.57), (2.59), (2.61) in the following in turn . By employing the backward Euler scheme for time integration and the Sharfetter-Gummel scheme for representing the particle flux and after tedious algebraic rearrangement, the resulting discretized equation of a species continuity equation on a typical grid point (m, n) can be written as,
where α can be electrons, ion, and uncharged species, xm is the distance from original point to grid point m, ∆xm and ∆ynare the step length at grid (m, n), γ is a factor using to deal with the coordinate system (for cylindrical coordinate is 1 and for Cartesian coordinate is 0) and Γα are the flux defined between grids which can be expressed in Scharfetter-Gummel form as,
Γα,m+1 where X is a non-dimensional variable define as,
Xm+1
and B is the Bernoulli function which is defined as B(x) = x
ex− 1. (2.111)
The electron energy equation is descritized as 3
where Te and i are in unit of eV , νm is the electron collision frequency, and the descretized form energy fluxes Γεare expressed as
Γε,m+1
Finally, the Poisson’s equation is discretized as γεm
2.3.4 Boundary Conditions
Boundary condition for electrons is applied on the electron flux to the wall as Γe,BC = 1
4neve,th− De∇ne− neµeE−∑
i
γiΓi (2.116)
where γi is the secondary emission coefficient when ion species i bombarding on the wall and the electron thermal velocity
ve,th =
√8kBTe
πme . (2.117)
The diffusion and drift terms of electron are considered only when the flux direction is to the wall. The electrons moving to the wall are removed while on the dielectric wall the charges are added to the net charge and are accumulated during runtime. Boundary condition for ions is also applied on the ion flux to the wall
Γi,BC =−Di∇ni− niµiE (2.118) where drift term of ion is considered only when the direction is to the wall. Similarly, the ion moving to the conductor wall are removed while on the dielectric wall the charge are added to the net charge and are accumulated during runtime.
For the non-excited and non-radical neutral species, Neumann boundary condition is em-ployed at all solid walls. However, boundary condition for the excited or radical species is assumed to be
Γp,BC =−Dp∇np (2.119)
where the resulting species on the wall is removed. Finally, the boundary condition of electron energy density equation is expressed as
Γ = 2kBT eΓe. (2.120)
2.3.5 Newton-Krylov-Schwarz (NKS) Algorthm
In this thesis, the parallel fully coupled Newton-Krylov-Swartz(NKS) algorithm [Cai et al., 1998a] was employed to solve the large sparse system of nonlinear discretized equations, which were derived in the previous section. All the functionals for each dependent variables form a
global functional vector. Jacobian matrix is then computed based on this global vector. Nev-ertheless, this treatment results in a fully implicit scheme, which may allow larger time step as compared to semi-implicit or explicit scheme. Another advantage of solving the coupled equa-tions directly, rather than solving the equaequa-tions one by one or some heuristic way of coupling, is that it can have better time accuracy with appreciable time-step size and much better parallel performance since the grain size is much larger. In this method, an inexact Newton method is used to solve the coupled nonlinear discretized equations at each time step. The resulting Jacobian system computed by using a hydrid analytical and numerical (finite difference) for the Newton corrections are solved with a preconditioned Krylov subspace type method, relying di-rectly only on iterative operations. With the use of hybrid analytical and numerical scheme, the time required for evaluating Jacobian matrix is reduced greatly. The Krylov method requires preconditioning for achieving acceptable convergence speed of inner interactions. A good pre-conditioner saves time and memory by allowing small number of iterations in the Krylov loop and smaller storage for the Krylov subspace. In this study, we have utilized a parallel additive Schwarz(AZ) type preconditioner with an inexact or exact solver such as incomplete LU(ILU) or LU factorizations in each subdomain. It was shown that this AS preconditioner can greatly reduce the runtime required for inner iterations in an inexact Newton method [Hwang and Cai, 2005]. This results from that the smaller subdomain blocks maintain better cache residency and shorter convergence time for an approximate solver. In addition, either scheme was used to solve the preconditioned matrix equation at each time step. By combining the AS precon-ditioner with Krylov type subspace method BiCG-STAB(or BCGS) [van der Vorst, 1992] or GMRES [Saad and Schultz, 1986] in the present study within an inexact Newton method leads to a coherent fully parallel solver: Newton-Krylov-Schwarz(NKS) algotithm [Cai et al., 1998b].
In practical implementation, we employed the PETSc package [Balay et al., 2001], which features distributed data structure - index sets, vectors, and matrices - as functional objects. Iter-ative linear and nonlinear solvers are combined modularly, recursively, and extensively through a uniform application programmer interface. In addition, the Jacobian system within the inex-act Newton method can be computed automatically by the finite-difference scheme within the PETSc framework or by the user himself/herself. Portability is achieved through MPI, which is a standard in message passing among processors nowadays, although the details are not requires
in practical coding. All the codes were programmed in C/C++.
Chapter 3
Verifications and Parallel Performance of the Fluid Modeling Code
To validate the developed fluid modeling code, we have conducted several simulations in one-dimensional and two-one-dimensional cases and results were compared with experimental data and previous simulations wherever possible. For the one-dimensional cases, they include: (1) he-lium discharge driven by RF power source, (2) hehe-lium discharge driven by AC power source and (3) nitrogen discharge driven by AC power source. For the two-dimensional case, a GEC chamber-scale simulation was conducted. At the end of Chapter 3, parallel performance of the developed parallel fluid modeling code is presented using the 2D-axisymmetric GEC chamber-scale simulation. These validations and parallel performance were described in the following in turn.
3.1 One-Dimensional Simulation of Helium Discharge Driven by a Radio-Frequency Power Source (13.56 MHz)
3.1.1 Simulation Conditions
A one-dimensional fluid modeling which mimics a parallel-plate atmospheric-pressure driven by a radio-frequency power source (13.56 MHz) was conducted and validated with experimen-tal data obtained in our group. Corresponding test conditions include: two electrodes of area 25 cm2 covered by dielectric layers with thickness of 1 mm each and relative permittivity of 11.63, power electrode applying a sinusoidal voltage waveform having frequency of 13.56 MHz and peak-to-peak voltage of 496 V , discharge gap of 1 mm, and background gas pressure of 760 T orr.
3.1.2 Plasma Chemistry
This simulation employed a “complex” set of Helium reaction channels as summarized in Table 4.1, which includes 7 species (electron, He∗m, He∗∗ex, He∗2, He+, He+2, He) and 27
reaction channels. This set of helium plasma chemistry includes momentum transfer collision (0), electron impact excitations (1)-(7), direct ionization (8)-(9), electron impact de-excitation (10), electron impact dissociation (11), electron-ion recombination (12) (15)-(16), electron-ion dissociative recombination (13)-(14), Hombeck-Molnar associative ionization (17), metastable - metastable associative ionization(18), metastable - metastable ionization (19), ion conversion (20), metastable - induced association (21), metastable - induced dissociative ionization (22), metastable-induced ionization (23), dimer - induced dissociative ionization (24), dimer-induced ionization (25), helium-atom induced dissociation (26).
3.1.3 Validation with Experiment Results
Figure 3.1 shows the comparison of the simulated discharged currents with the measured data that was obtained in our group. It clearly shows that the present fluid modeling code using the complex helium plasma chemistry can predict quantitatively the temporal evolution of discharge current of helium RF discharge very well.
3.2 One-Dimensional Simulation of Helium Dielectric Barrier Discharge Driven by AC Realistic Distorted Sinusoidal Voltages
3.2.1 Simulation Conditions
A one-dimensional fluid modeling which mimics a parallel-plate atmospheric-pressure driven by a distorted sinusoidal AC power source was conducted and validated with experimental data obtained in our group. Corresponding test conditions include: two electrodes of area 25 cm2 covered by dielectric layers with thickness of 1 mm each and relative permittivity of 12, power electrode applying a distorted sinusoidal voltage waveform having frequency of 60 KHz and peak-to-peak voltage of 14000 V, discharge gap of 2.5 mm, and background gas pressure of 760 T orr.
3.2.2 Plasma Chemistry
Again, we have applied the same set of “complex” Helium plasma chemistry as the previous RF helium discharge.
3.2.3 Validation with Experiment Results
Figure 3.2 shows the comparison of the simulated discharged currents with the measured data that was obtained in our group. It clearly shows that the present fluid modeling code using the complex helium plasma chemistry can predict quantitatively the temporal evolution of discharge current of helium DBD very well.
3.3 One-Dimensional Simulation of Nitrogen Dielectric Barrier Discharge Driven by AC Realistic Distorted Sinusoidal Voltages
3.3.1 Simulation Conditions
Similar to those test conditions as described in Section 3.1.1 for the helium RF discharge, one-dimensional fluid modeling is conducted for the same configuration using pure nitrogen gas, which was driven by a distorted sinusoidal power source (60 kHz) with a gap distances of 0.5 mm.
3.3.2 Plasma Chemistry
The nitrogen plasma chemistry employed in the present study includes 9 species (electron, N2, N2+, N4+, N2(X1Σ+g, ν = 1− 8), N2(A3Σ+u), N2(B3Πg), N2(C3Πu), and N2(a01Σ−u)), and 31 reaction channels, which are summarized in Table 3.1. This set of nitrogen plasma chem-istry includes direct ionization (1), excitation into excited, metastable and vibration states (5), de-excitation (6), recombination (3), associative ionization (3), light emission from excited, metastable states (4) and excitation into vibration states (10). Note we have ignored the N+and N3+ in the simulation since they have been found to unimportant in nitrogen plasma simulation [Choi et al., 2006].
3.3.3 Validation with Experiment Data
Figure 3.3 shows that comparison of simulated and measured discharged currents of nitro-gen DBD driven by a quasi-pulsed power (60 kHz) for different gap distances (0.5, 0.7, 1.0 and 1.2 mm) along with the experimental photo images (0.2 sec exposure time) of discharge
on the right. Note the measured relative permittivity and thickness of the ceramic material is 11.63 and 1 mm, respectively. In general, the simulations agree very well quantitatively with the experimental data for the cases of d = 0.5 mm and 0.7 mm, but begin to deviate slightly from the measurements as d = 1.0 mm and exhibit large discrepancy with the measurements as d = 1.2 mm. For the cases of smaller gap (0.5, 0.7 and 1.0 mm), the simulations demon-strate they are typical homogeneous Townsend-like discharges with much fewer electrons than ions (not shown here). For the case of larger gap (d = 1.2 mm), the simulation shows it is a glow-like discharge (quasi-neutral in the bulk) with very high current density during the break-down phase. However, this is obviously against the measurements. This is attributed to the fact that the discharge has transitioned from Townsend-like to filamentary-like (microdischarge), as shown in the photo images in Figure 3.3, which makes the one-dimensional fluid modeling invalid. This shows that one has to be very cautious about the use of one-dimensional fluid modeling for simulating parallel-plate nitrogen DBD. In ref. [Choi et al., 2006], the conclusion of transition from Townsend-like to glow-like discharge was misleading without the validation of experiments. To capture the correct physics with a larger gap, one should employ at least two-dimensional fluid modeling, which is currently in progress and will be reported elsewhere in the near future.
Figure 3.4 shows the spatial-average temporal discharge properties as d = 0.7mm. Results clearly demonstrate that the total number density of ions (N2+ and N4+) is always larger than electron number density throughout a cycle. The simulated electric field across the gap is almost linear without any distortion by the charge density at all times during a cycle (not shown here).
Both the above two phenomena show that this is a typical Townsend-like discharge. In addition, N2+ is found to be most abundant during the breakdown process (same period of those current peaks in Figure 3.3), while N4+ is found to be dominant after the breakdown caused by the associative ionization of excited/metastable nitrogen species.
3.4 Two-dimensional Helium Discharge Simulation in a Gaseous Elec-tronics Conference (GEC) Reference Cell
3.4.1 Simulation Conditions
A GEC discharge simulation is conducted for validating the two-dimensional fluid model-ing code as developed in the thesis. The simulation conditions are summarized as follows. The power and grounded electrodes are both 2 inches in radius, the electrode gap is 1 inch in length, the applied peak-to-peak voltage is 150 V, the radio frequency is 13.56 MHz, and the back-ground gas is pure helium at 500 mTorr. This validation employed a “complex” set of Helium reaction channel listed in Table 4.1, in which the details are discussed in Chapter 4. Gas temper-ature is assumed to be 400 K. A non-uniform grid and the total of 200 time steps per RF cycle are employed. Additive Swartz with a minimal overlapping and LU as a subdomain solver were used. The relative stopping condition of the Newton method and the Krylov subspace method were 10×−5 and 10×−4, respectively.
3.4.2 Validation with Experimental Data
The numerical simulations of the cycle averaged parameters are shown in Figure 3.5, includ-ing (a) electron, (b) He+, (c) He+2, (d) He∗2, (e) He∗, and (f) Hemeta. The results indicate that the maximum value occur near the outer edge between two electrode gaps, where the electric field is the strongest. In addition, the dominant ion species are the molecular ions, rather than the atomic ions. Similar observations were found in the previous study [Martens et al., 2007b].
Figure 3.6 shows a comparison of the simulated peak electron densities with the theoretical prediction and the experimental data [Riley et al., 1994] for various applied voltages. Note that only the 1-D electron Boltzmann equation was used for the theoretical prediction, which was questionable for the present 2-D case. Thus, the data were included for reference only. In gen-eral, the simulation data followed the trend of the measurements reasonably well, except for the lowest (50 V) and the highest (200V) cases. Note that the electron densities were measured by a microwave interferometer probing through the center along the radial direction at the midpoint between the electrodes by integrating the data along the line-of-sight path. The data were only
correct if the densities were uniform throughout the microwave path, which was obviously not true for the current test cases. Nevertheless, the simulated data were in reasonable agreement with the measurements.
3.5 Parallel Performance Test of Fluid Modeling Code
3.5.1 Test Conditions
The simulation conditions used for parallel performance studies are similar to those for the validation of the parallel 2D-axisymmetric fluid modeling code, except that we now consider a 122× 123 uniform grid and 200 V of peak-to-peak value of applied voltage. 400 timesteps per RF cycle were used throughout the simulation. We investigated the parallel performance of the solver by using two types of Krylov subspace method, GMRES and BiCGStab in conjunction with standard AS preconditoners, where the subdomain problems were solved by either the LU decomposition or ILU(0) (incomplete LU decomposition with zero level fill in). All the calculations were done on the V’ger cluster system (Xeon 3GHz, dual core, dual CPU) at the Center for Computational Geophysics, National Central University, Taiwan.
3.5.2 Results and Discussion
Figure 3.7 illustrates the parallel performance including speedup analysis and runtime per time step as a function of the number of processors. In the parallel computing we have applied either GMRES or BiCGStab as the KSP solve and LU of iLU as sub-domain solve. Runtime per timestep was calculated by averaging the total time of 5000 timesteps. The cases of 96 processors using GMRES or BiCGStab with iLU combination performs the best. Generally, resulting speedup shows a typical super-linear behavior when the number of processors is in the range of 32˜144 except the case of BiCGStab linear solve with LU preconditioner.
3.6 Brief Summary of This Chapter
In Chapter 3, several validations of the fluid modeling codes by comparing with experimen-tal data or previous simulations and parallel performance were presented. Results show that the developed fluid modeling code is able to predict one- and two-dimensional
atmospheric-and low-pressure gas discharges very accurately. Superlinear speedup up to 144 processors of the parallel fluid modeling code with proper selection of the simulation parameters was also demonstrated.
Chapter 4
One-Dimensional Simulation of Helium Dielectric Barrier Dis-charge Driven by AC Realistic Distorted Sinusoidal Voltages
This chapter discusses the effects of selecting plasma chemistry in simulations of helium DBD by comparing simulations with experimental results. Simulated temporal discharged currents using the complex plasma chemistry are in excellent agreement with the measurements obtained in the present study, which validates the fluid modeling code. It also indicates that the inclusion of heavy particle related reactions is critical in accurately predicting the helium DBD because of the slow varying electric field in the range of the tens of kilohertz. Detailed temporal variations of spatial-average plasma properties are discussed by partitioning in time.
4.1 Background and Motivation
Atmospheric-pressure plasmas (APP) have attracted tremendous attention in the past two
Atmospheric-pressure plasmas (APP) have attracted tremendous attention in the past two