平行化耦合氣體放電與氣流模擬的發展與其於含雜質常壓氦氣介電質電漿束

全文

(1) . . . . . . . . . . . . . . . .

(2). . . . . . . .

(3). . . . . . . . Development of Parallel Hybrid Simulation of Gas Discharge and Gas Flow and Its Application in the Modeling of Atmospheric-Pressure Helium Dielectric Barrier Discharge Jet Considering Impurities . .

(4). . . . .

(5). . . .

(6)

(7) . . . . . .

(8) . . . . Development of Parallel Hybrid Simulation of Gas Discharge and Gas Flow and Its Application in the Modeling of Atmospheric-Pressure Helium Dielectric Barrier Discharge Jet Considering Impurities. . .

(9). . StudentKun-Mo Lin . AdvisorDr. Jong-Shinn Wu . A Thesis Submitted to Department of Mechanical Engineering National Chiao Tung University in partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy in Mechanical Engineering. . June 2012 Hsinchu, Taiwan. . 2012 !. ".

(10) . . . . . . . .

(11). . . . . . . . . . . . . . . . . . . . . . . . # #. . . # #. #. #. #. #. #. . #

(12). . ##. ## !" . . . . / . . 01. %. &. . 2. 3. I J. . . . -. 4. K. L. M. . . 5. N. et al., 2011]3 a. . . . . . . . 6 +7. . . . &. . '. . 8. 9. . (. :.

(13). . . . . . . ). . .

(14). . . . . . * +,. 0 (domain decomposition) < ;. . . =. . -. >. .. ?. @. (message passing interface, MPI) B 5 6 +C D E : 2 6 F G A. H. #$. . . . . . (O . . . . . . . . . . 5. 6 . \. . 7. P . ]. g. -. Q . S.

(15). . h. R ^. T. U . 5. ' . '!#i. V. WX 'Y W+!#Z 6.

(16). . _j. . . k. _. *. l. 3. . b. 6. c. . . d. U. V. I[Hu [ V. 0`. WX. e. . f. 0 (temporal multi-scale. WX. method, TMSM)+, m n o m p q WX k l r s p t +a u v - Q R ] w. x. . . . y. . ". . V. 0. §. . . ®. . . . . . . . . . . . u. `. Z.

(17). . &. . ¨. . . . §. . . . . . . {. . . z. u. . * E. . . ¤

(18). . 0'. '. V. . #. 0 . j. :. . *. . . '. . . Z. . . '#(©. &. .

(19). . «. n ¬. `. . 5 . c. . . *. ¡ . ª. . . ¦. WX '!". . . . . V. + ". u ¥. d. #3. |. C. . }. ~. . '. #. . . . . B . . . . . c. . ". . . b. . <. . . . . . 3. . 6 n. f . c. . ¢. . £. t +!#. &. . . . &. . . . 4. . 6. . ' . . ". &. . 0(gas chromatography)3. u. . (99.99%)¢ +7 P ¯ ¢ 3 i.

(20) . . (t. (quasi-neutral). Â. ¼ . . ` Y. Î. . ². Î. . F. G. . . Ü+Ý. 8. 9 '. ¶. ·. . . ¸. Z . . +¹. ± . E. . . À . µ. º. ". -. . . ¡z. Þ. . É. E. . »¼. ¢. . (»Æ. Ê. ±. .

(21). . c. . . . 1. . ½. X. m. m . Ç. ¾. Á «. ¿. l. ¨. . +E. . . . M. N. ¥. t. . Î. Ä. ¢. . . . m. . . Ö. . Û. . . Ì. '#(Í. . . . ©. m. . Á. (hydroxyl) ×. ÅT. ~. . . . Â. l. r. Á. l. Ï. ppm + y. 'Å©. .

(22). . . áâ (preconditioning) 4 à. (t. '. . . c. He2 *. N2+C Å (»Æ c 1 m 'Ò. . . Y. +Ö. c. 0. Â. Ê. n. Å»Æ. À. . Ã . . L. . 'Ú. . .. ß. m. (t. Ñ . K. Ä. Ñ +H2O+É. Ð. -. Ì. 'Z. c. . . Ë. . Ð. Ù. I J. . . . Ø. &. . »Æ. . . (. f. ®. H. I&. . Ñ . . f. . . Õ. Ð. !# 6. . Ô. . ". . . m. À. . . Ä. . ±. . . Ì. ´'µ. N2+3 É Ê ' E Â (»Æ c (t . . . Ë. ³. (È. ¼ 'Ï. Ó. . ! He2+ ÅÂ. . . ². Ä. . . Â. Z. . . Î. ±. . È. ². -. 'Ã. . s. °. . . * . 5. 6. . <. . . E . : . 2 . Q. 0 (ã ASM; additive Swartz. method) Þ ß ä Ê ; 0 (ã GMRES; generalized minimal residual method) å Â. Y. *. ï. j. . N. . è. æ é. . ý. .

(23). . . R. *. X. . . . a. . . . . . . . ç +è. . . . . . ^. é. . . * +7

(24). . &. . 128 ëáâì. . þ. . . . °. . ¤. &. . . . ù. 6. ñ+E. ¤. . Å. 5. ð. '»ø +ù âx. 4. . P. &. . . . §. z. ¨. 4. . ó. ô. õ. ð. \. . ñú. §. û. ¨. 3. ö. . t. . Ñ. . ÷. 5. . 6. Q. ¦ü. á. . ' #. 0Ñ. 128 ëáâì' í î ê. . . . +!". . 3. ò. . . »j. Ë. ÷ ;. . i d. _ V. . 3

(25). b. . . . . V. 0Z. WX ' *. 3. §. ¨. . Û. . (source terms) Z

(26) ê ¨ Z * ' d V WX 4 0x . .

(27). 0'. m. B. 5. H. 6 n. *. . * . o. Û. m. . & 7. . . p. O. 8 . q. * 9. U. WX. :. . :. 6 l. . d. V +. <. . . . .

(28). . 0(domain decomposition) < ;. 2 k. . r. F s. G p. t. . H. I J. + ii. . ;. n. K o. L m. M. =. >. N 'j . . 4. * ?. k 5. @. l 6. j. k . . WX. (MPI) A. WX . l. H. 0, WX. v.

(29) (time step)W+ · ¸ ° þ ? @ ¼. . . . . . ". v. . 'E. . . . . (&. . j. k. 3 l. . . WX. . ?. 0. @ . . B . E. . . . . u. . . \. 47% d V WX '#(3 i _ V 0ù ¤ 6 . z. WX. . . Á. . . l. û. ·. #»ø. : ¸. #. . *. \ ò. . .

(30) +. !". t +7. . $. . P. #Æ. . n. . .

(31). t c. c . 200 H ¯ ^. :. %. >. ! >. " ª. 3 & +'. §. ¨ (. ). Z ±. '¤Z. . ¯. _. . . _. .'. . . ¡. . ' ¯ +7. . . iii. ª. &. . ,.

(32). . }. ~. 4. -. i.

(33) Development of Parallel Hybrid Simulation of Gas Discharge and Gas Flow and Its Application in the Modeling of Atmospheric-Pressure Helium Dielectric Barrier Discharge Jet Considering Impurities. Student: Lin, Kun-Mo. Advisor: Dr. Wu, Jong-Shinn. Department of Mechanical Engineering National Chiao Tung University Abstract The development of a hybrid algorithm of plasma fluid model (PFM) and gas flow model (GFM) for simulating the atmospheric-pressure dielectric barrier discharge jet (APDBDJ) is reported in this thesis. The gas discharge is modeled by plasma fluid model, discretized through finite-volume method, and parallelized with domain decomposition using message passing interface (MPI) and employed on distributed-memory PC cluster that reduces runtime significantly. The hybrid numerical algorithm is proposed by combining a previously developed parallelized compressible flow equation solver [Hu et al., 2011] to simulate the helium APDBDJ considering impurities in this thesis. A temporal multi-scale method (TMSM), taking advantage of the difference of characteristic timescale between electron and heavy particles, is proposed to further reduce the runtime of PFM dramatically in simulations involving with a large amount of species. The effort of this thesis establishes the foundation in simulating realistic APDBDJs. Chapter 1 of the thesis introduces the research background and motivation, and review of previous studies, and Chapter 2 describes the numerical methods and. iv.

(34) developed algorithm in detail. The results of one-dimensional simulation considering the impurities on the helium APDBDJ are presented in Chapter. 3. Chapter 4 and 5 depict the development of 2D parallel PFM code and the hybrid numerical algorithm of PFM and GFM for simulating APDBDJ respectively. The major findings of the thesis and some recommendations for future work are summarized in Chapter 6. More details of each chapter from Chapter 3 are described as follows in turn. In Chapter 3, the effect of helium impurities (trace amounts of O2, N2 and H2O), measured by the gas chromatography, have been explored and found that the discharge chemistry changes dramatically by considering the impurities. Results have found that the discharges with and without impurities have no quasi-neutral region at the instant of maximum current density. Both discharges with and without impurities have similar levels of electron densities; however, N2+ is found to be the most dominant ion with considering impurities, instead of He2+ in helium discharge without impurities in the breakdown region. In addition, ground-state atomic oxygen is the most dominant neutral species (except the background species) when considering impurities, instead of He2* without considering impurities. The influence of different levels of water vapor is also investigated. The electron densities of helium discharges with various levels of water vapor (1, 5 and 10 ppm) remain at essentially the same level as the amount of water vapor changes. However, the H2O+ replaces the N2+ as the dominant ion as the water vapor increases. Although the ground-state atomic oxygen is still the dominant neutral species, the densities of atomic hydrogen and hydroxyl increase significantly as the water vapor increases. The results show the importance of considering impurities in the helium discharges though the levels of impurities are typical several to tens ppm. In Chapter 4, the thesis reports the development of a two-dimensional plasma fluid modeling code using the cell-centered finite-volume method and its parallel v.

(35) implementation on a distributed-memory PC cluster. Parallel performance of simulating helium APDBDJ resulting from using different degrees of overlapping in the additive Schwarz method (ASM) with preconditioned generalized minimal residual method (GMRES) for different modeling equations is investigated for a small and a large test problem, respectively, employing up to 128 processors. For the large test problem, almost linear speedup can be obtained using 128 processors. Finally, a large-scale realistic two-dimensional APDBDJ problem is employed to demonstrate the capability of the developed fluid modeling code for simulating the low temperature plasma with complex chemical reactions. In Chapter 5, this thesis proposes a hybrid numerical algorithm which couples weakly the PFM and GFM, and two acceleration approaches for simulating the APDBDJ. The weak coupling between gas flow and discharge is introduced by transferring between the results obtained from the steady-state solution of the GFM and cycle-averaged source terms of the PFM respectively. Approaches of reducing the overall runtime include parallel computing of the GFM and the PFM solvers, and employing a TMSM for PFM. Parallel computing of both solvers is realized using the domain decomposition method with message passing interface (MPI) on distributed-memory PC cluster. The TMSM considers only the source and sink terms of chemical reactions by ignoring the transport terms when integrating temporally the continuity equations of heavy species at each time step, and the ignored transport terms are restored only at an interval of several time marching steps. The total reduction of runtime is 47% by applying the TMSM to the example of APDBDJ as presented in this study. Application of the proposed hybrid algorithm is demonstrated by simulating a parallel-plate helium APDBDJ impinging onto a substrate, in which the cycle-averaged properties of the 200th cycle are presented. The distribution patterns of species densities are strongly correlated by the background gas flow vi.

(36) pattern, which shows that consideration of gas flow in APDBDJ simulations is critical. In Chapter 6, major. findings of this thesis are summarized and. recommendations for the future work are outlined.. vii.

(37) /. f. ° &. 0. @. c. ¥. . §. c. 5 . .'. . . . 5. '. E. . . e. . E. . S. c. . ¼. §. ¨. ½. d. . . 5 Ç Ø i. 5 . . E. ç 0. Ù <Â. Ú. }. . . T. ~. . . e. B. . ë. . . . à. ·. . 8. © +].

(38). ². ^. ³. ¢. . . ' y. c. 5 w. \. . . \. ®. . ¼. ´ <µ. . . B. ¶ <·. . . 5 ] . ¸ <µ. \. . . . ¹ <º. b. ¯. . °. @. . S. 3. . WX 'ª. ]. . j. â . . . y. . c. . B. \. û. b. Â. . +f ¦. 3. ¯. . ¥. 6. i. . '. A. . . I. . °. . H <~. P. . t i. . . .

(39). p. . . P. b. G. 7. ¥. 7. M. }. ". . ¤. . . . . 6£. \. . #. J. 9. a. o. @. W. Å`. ". ¦. . _. |. . @. . . H <{. j. °. . <.

(40). °. ^. 5. '. » <¼. +66. C. 2. Í <m. . Á. Â <£. â <. Î <Ï. B. ¤. i. Ð <Ñ. ã 'ä. å. æ. . Ã. Ä. Ò <Ó ç. ¤. i. #. Ô <Õ. <. . ¡:. Å <Æ. ¸ <û. Ö <×. i. , +é. è. T. µ < ,. < ¤. '. í. î. ù d. ¡. . +Â. i. o ' . 3. b. °. 9. <á. ç. ". 3. G C. . V. +. '. @. ß. Wø '½. è. . . . z. I'. B. »Å&. . . U. i. y. A. ^ +R. n. . 3. . °. . . <µ. !#ì. ÷ . À. Ë <Ì. ,. . °. Û Ü Û Ý Þ . :. ç. . . . &. . ¬. . 6. ~. m. ?. . , ]. H <x H. 5. . c. S. «. . . n. D. +. ¿. WX. . j . . ¾. Ê <. »ø f. 6±. 6/. ê. Z. 6. G. T. l. >. ¥. S. . <=. . @. . (. p. . w. 2. I. °. B. ;. °. \. âG. . c. . y i. b. e . . E. i. k. . E. @. 5. ü. . ¡:. H. 3. j. H <v. e. ;. . ¤. j. È <É. . . B. y. . ¨. £. . ªc. ¢. S. ¥ +°. . z. f. ;. 7. . |. 6¢. 3. X. G. H < v. C. 4. u. . Q +R. ;. i. 9 G. P. C. }. 8 F. O. [. . ¥ E. N. 5. . 7 D. 5. . . G. ;. . 6t. . WX. §. 5 6. . »w. D. s. 64 C.

(41). h. 6¡:. . B. ¤Z. g. . +5. . Y. r. . . M. f é. 4. #A. L. . 6q H <. 3. a '¢. e. G. . . X. S. . 6". ;. . d. 2. 5. K ' C. ð. j. 1. . . ú. ï. û. . ?. I+º. -. ü '66Ò f. c. [. a. ½ T »n. viii. d ý. . c. þ. ð. . . ñ. C. c. ò. Â. ¨ '66.

(42). +ö '. f . . . ó l. ô ë. õ y. ö.

(43) Table of Contents. c. ............................................................................................................................i. Abstract ....................................................................................................................iv 6......................................................................................................................viii List of Tables ............................................................................................................xi List of Figures..........................................................................................................xii Nomenclature .......................................................................................................... xv Chapter 1 Introduction ...............................................................................................1 1.1 Classification of Plasmas ..............................................................................1 1.2 Non-Thermal Low Temperature Plasmas (LTP) ............................................2 1.2.1 Direct Current Discharges ..................................................................3 1.2.2 Dielectric-Barrier Discharges (DBD) .................................................3 1.2.3 Radio Frequency (RF) Discharges......................................................4 1.2.4 Microwave Discharges.......................................................................4 1.3 Importance of Helium Discharges.................................................................5 1.4 Numerical Simulation of Gas Discharges......................................................6 1.5 Literature Survey..........................................................................................7 1.5.1 Simulation of Helium Discharges with Impurities ..............................7 1.5.2 Parallel Computing of Fluid Modeling ...............................................8 1.5.3 Parallel Hybrid Numerical Algorithm for Simulating APPJ.............. 10 1.6 Objectives and Organization of This Thesis................................................ 11 Chapter 2 Modeling and Numerical Methods ........................................................... 13 2.1 Governing Equations of Fluid Modeling ..................................................... 13 2.2 Boundary Conditions.................................................................................. 16 2.3 Implementation of Semi-Implicit Schemes ................................................. 17 2.4 Discretization and Numerical Schemes....................................................... 18 2.5 Parallel Implementation of Fluid Modeling................................................. 20 2.6 Parallel Gas Flow Model (GFM) ................................................................ 21 2.7 Hybrid Algorithm for Coupling the GFM and PFM Solvers........................ 23 2.7.1 General Description ......................................................................... 23 2.7.2 Inclusion of Convection Effect in the PFM ...................................... 23 2.7.3 Details of Hybrid Algorithm for Coupling the GFM and PFM Solvers ................................................................................................................. 24 2.8 Temporal Multi-Scale Method for the PFM ................................................ 26 Chapter 3 One-Dimensional Fluid Modeling of Helium Dielectric Barrier Discharge Considering Impurities............................................................................................. 29 3.1 Background and Motivation ....................................................................... 29 ix.

(44) 3.2 Problem Description................................................................................... 30 3.2.1 Plasma Chemistry ............................................................................ 31 3.2.2 Simulation Conditions ..................................................................... 32 3.3 Helium Discharge without Considering Impurities...................................... 34 3.4 Helium Discharge with Impurities .............................................................. 35 3.5 Helium Discharges with Different Levels of Water Vapor ........................... 36 3.6 Summary.................................................................................................... 37 Chapter 4 Development of a Parallel Semi-Implicit 2D Plasma Fluid Modeling Using the Finite-Volume Method ....................................................................................... 39 4.1 Background and Motivation ....................................................................... 39 4.2 Problem Description................................................................................... 40 4.2.1 Plasma Chemistry ............................................................................ 41 4.2.2 Simulation Conditions ..................................................................... 42 4.3 Small Problem Case ................................................................................... 42 4.4 Large Problem Case ................................................................................... 44 4.5 Parallel Performance .................................................................................. 44 4.6 Effect of ASM Overlapping ........................................................................ 45 4.7 Demonstration of the Results of the Large Problem in 5th Cycle ................. 46 Chapter 5 A Parallel Hybrid Numerical Algorithm for Simulating an Atmospheric-Pressure Plasma Jet............................................................................. 49 5.1 Background and Motivation ....................................................................... 49 5.2 Problem Description................................................................................... 50 5.2.1 Plasma Chemistry ............................................................................ 51 5.3 Steady-State Results of Flow Field ............................................................. 52 5.4 Verification of TMSM ................................................................................ 53 5.5 Demonstration of APPJ Simulation............................................................. 55 5.6 Summary.................................................................................................... 56 Chapter 6 Conclusion and Recommendations of Future Work .................................. 58 6.1 Summaries of the Thesis............................................................................. 58 1 One-Dimensional Simulation of Helium Discharge with Impurities.. 58 6.1.2 Development of a Parallel Semi-Implicit 2D Plasma Fluid Modeling Using Finite-Volume Method.................................................................... 59 6.1.3 A Parallel Hybrid Numerical Algorithm for Simulating an Atmospheric-Pressure Plasma Jet ............................................................. 60 6.2 Recommendations for Future Work ............................................................ 62 References ............................................................................................................... 64 Tables ...................................................................................................................... 71 Figures..................................................................................................................... 85 x.

(45) Publication List of Kun-Mo Lin ............................................................................. 122. xi.

(46) List of Tables Table 1-1 Classification of plasmas with their temperatures [3]................. 71 Table 3-1 Measured impurities of helium gas in Taiwan using a gas chromatographer. .............................................................................. 72 Table 3-2 Summary of helium plasma chemistry considering impurities ... 73 Table 4-1 Averaged iteration number per time step of equations with sub-domain of ASM preconditioner solved by ILU and LU methods for small problem case (501 ×310 cells). ................................................ 79 Table 4-2 Averaged iteration number per time step of equations with sub-domain of ASM preconditioner solved by ILU and LU methods for large problem case (1001 ×620 cells). ............................................... 80 Table 4-3 The parallel performance of different level of overlapping ASM preconditioner with sub-domain solved by ILU method on the large problem (1001 ×620 cells) (a) Averaged iteration number per time step of equations (b) The relative runtime per time step. “OL 1”, “OL 2”, and “OL 3” represent one-, two-, and three-level overlapping, respectively....................................................................................... 81 Table 4-4 The parallel performance of different level of overlapping ASM preconditioner with sub-domain solved by LU method on the large problem (1001 ×620 cells) (a) Averaged iteration number per time step of equations (b) The relative runtime per time step. “OL 1”, “OL 2”, and “OL 3” represent one-, two-, and three-level overlapping, respectively....................................................................................... 82. xii.

(47) List of Figures Fig. 1-1 Plasmas classified with the electron density and temperature....... 85 Fig. 1-2 The setup of direct current discharge ........................................... 86 Fig. 1-3 Distributions of physical properties in a direct current discharge . 87 Fig. 1-4 The I-V characteristic of a direct current discharge ...................... 88 Fig. 1-5 Typical arrangements for electrodes and dielectrics in dielectric-barrier discharges. ............................................................. 89 Fig. 1-6. Reactor of (a) capacitively coupled plasma (b) inductively coupled plasma .............................................................................................. 90 Fig. 1-7. Research structure of the thesis................................................... 91 Fig. 2-1 Flowchart of fluid modeling simulation ....................................... 92 Fig. 2-2 The coupling of the GFM and the PFM solvers............................ 93 Fig. 2-3 The flowchart of complete APPJ simulation ................................ 94 Fig. 2-4 Time sequence of the temporal multi-scale method...................... 95 Fig. 3-1 Optical emission spectroscopy measurement of helium APDBDJ.96 Fig. 3-2 Sketch of the planar helium APDBDJ system [39]....................... 97 Fig. 3-3 Schematic diagram of the one-dimensional simulation model of a helium DBD. .................................................................................... 98 Fig. 3-4 Simulated and experimental discharge current densities for helium discharges with and without impurities driven by Vpp = 6.0 kV (peak-to-peak voltage). The frequency of power source is 25 kHz..... 99 Fig. 3-5. Discharge structure at the instant of maximum current density. (a) Case without impurity; (b) Case with impurity (10 ppm O2, 25 ppm N2, and 1 ppm H2O). The legend “Positive” represents the density of all species with positive charge, and the “Negative” represents the density of all species with negative charge. ................................................. 100 Fig. 3-6. The distribution of conduction (of each species), displacement, and total current density at the instant of maximum discharge current density for cases without (a) and with (b) impurities. ...................... 101 Fig. 3-7 Distributions of spatial-averaged number densities of helium discharge without impurity of abundant (a) charged species; (b) neutral species. ........................................................................................... 102 Fig. 3-8 Distributions of spatial-averaged reaction rates of important channels (a) electron; (b) He2 *. The helium discharge contains no impurity. ......................................................................................... 103 xiii.

(48) Fig. 3-9 Distributions of spatial-averaged number densities of helium discharge with impurities (10 ppm O2, 25 ppm N2, and 1 ppm H2O) of abundant (a) charged species; (b) neutral species. ........................... 104 Fig. 3-10 Distributions of spatial-averaged reaction rates of important channels (a) electron; (b) N2+; (c) O. The helium discharge contains 10 ppm O2, 25 ppm N2, and 1 ppm H2O............................................... 106 Fig. 3-11 Distributions of spatial-averaged number densities of helium discharge with impurities (10 ppm O2, 25 ppm N2, and 5 ppm H2O) of abundant (a) charged species; (b) neutral species. ........................... 107 Fig. 3-12 Distributions of spatial-averaged number densities of helium discharge with impurities (10 ppm O2, 25 ppm N2, and 10 ppm H2O) of abundant (a) charged species; (b) neutral species. ........................... 108 Fig. 4-1 Sketch of the simulation domain................................................ 109 Fig. 4-2 The partition of computational domain. ..................................... 110 Fig. 4-3 Averaged time per time step of equations with sub-domain of ASM preconditioner solved by ILU and LU methods for small problem case (501 ×310 cells) .............................................................................. 111 Fig. 4-4 Averaged time per time step of equations with sub-domain of ASM preconditioner solved by ILU and LU methods for large problem case (1001 ×620 cells). ........................................................................... 112 Fig. 4-5 Speedup of cases with sub-domain of ASM preconditioner solved by ILU and LU methods for small problem case (501 ×310 cells). .. 113 Fig. 4-6 Speedup of cases with sub-domain of ASM preconditioner solved by ILU and LU methods for large problem case (1001 ×620 cells).. 114 Fig. 4-7 Proportions of time consumed by different types of equations of the demonstration case (1001 ×620 cells). Note that “others” includes evaluation of transport properties and data communication among processors....................................................................................... 115 Fig. 4-8 Cycle-averaged spatial distributions of (a) potential (b) electron temperature (c) number density of electron (d) number density of N4+ in the 5th cycle................................................................................ 116 Fig. 5-1 Sketch of the computational domain .......................................... 117 Fig. 5-2 The converged results of the GFM solver after three interactions with the PFM. (a) Temperature (b) The Vm = u 2 + v 2 and the streamlines, where u and v are the velocity in the x and y direction xiv.

(49) respectively..................................................................................... 118 Fig. 5-3 The cycle-averaged number densities of species at 10th cycle simulated without (w/o) and with (w) the temporal multi-scale method (a) Electron (w/o), (b) Electron (w), (c) N4+ (w/o), (d) N4+ (w), (e) N2(A3) (w/o), (e) N2(A3) (w).......................................................... 119 Fig. 5-4 Breakdown of runtime for solving equations and other computer operations in one cycle. Note that “Poisson”, “Continuity”, and “Te” represent the Poisson equation, species continuity equations, and the electron energy density equation respectively. The “Others” includes evaluation of transport properties and MPI transmission among processors....................................................................................... 120 Fig. 5-5 Cycle-averaged spatial distributions of (a) electron temperature (Te), (b) number density of electron, (c) number density of N4+, (d) number density of N2(A3) at the 200th cycle.................................... 121. xv.

(50) Nomenclature m. mass. P. pressure. q. charge. kB. Boltzmann constant. n. number density. Γ. particle flux. vth. thermal velocity. γ. factor of secondary electron emission. S. source/sink of continuity equation. Te. electron temperature. ε. energy loss for inelastic electron collision. φ. potential. ε0. vacuum permittivity. εd. dielectric permittivity. E. electric field. f. electron energy distribution function. k. reaction rate coefficient. µ. mobility. D. diffusivity. νm. momentum exchange collision frequency between electron and background neutral particles. ∆x. space interval in x direction. ∆y. space interval in y direction. xvi.

(51) ∆t. time interval. B. Bernoulli function. xvii.

(52) Chapter 1 Introduction 1.1 Classification of Plasmas Plasma is generally considered as the fourth state of matter with the collective behavior of charged and neutral particles. Specifically, plasmas are a collection of quasi-neutral particles containing charged and neutral particles in the form of gaseous or fluid-like mixtures. The plasma state can be produced by heating the substance to the temperature high enough to ionize the neutral atoms, or it can be created by applying the sufficiently high electric field to cause ionizations by energetic free electrons absorbing the electric energy from the applied external electric field. The scales of plasmas observed may be as huge as the intergalactic nebula of cosmic plasmas, or plasmas can be generated in the small devices such as the microplasmas with their scale down to few mm in the laboratory [1]. The plasmas found cover a wide range of electron densities and temperatures as shown in Fig. 1-1 [2]. The plasma density can be as low as 104 cm-3 as that of ionosphere, which stretches from a height of about 50 km to more than 1,000 km of earth’s surface, and the plasma density can be as high as 1018 cm-3 for the thermonuclear reactors. The energy level of electron temperature varies from 0.1 eV of ionosphere to more than 104 eV in the thermonuclear reactors. Plasmas are frequently classified as low- (LTP) and high-temperature plasmas (HTP) as shown in Table 1-1 [3]. The HTP (for example, fusion plasmas) are in their thermal equilibrium, meaning that the temperature of ions is about equal to that of electrons, with the temperature is up to 107 K. The LTP can be subdivided into states in thermal equilibrium and non-thermal equilibrium. The LTP in thermal equilibrium (or thermal LTP) are produced mostly in atmospheric pressure with their temperature 1.

(53) of ions, electrons and neutrals are at the level of 104 K. The thermal LTP, such as arc plasmas and plasma torches, are used in welding, cutting hard materials, and material processing like melting and dissociation of minerals in industries [4]. For LTP with the thermal equilibrium not reached between the electrons and the heavy particles (ions and neutrals) are classified as non-thermal LTP. The temperature of heavy particles of non-thermal LTP is typically close to room temperature, whereas the electron temperature is much greater than that of heavy particles and is as high as 104 ~ 105 K. The non-thermal LTPs are especially important for the industrial applications because the low temperature of heavy particles has no (or tolerable) damage on the materials and the electron temperature is high enough to generate essential reactive species for applications. The non-thermal LTP contribute much to the major processes of microelectronics fabrication such as sputtering; plasma enhanced chemical vapor deposition, plasma etching, ashing, implantation, and surface cleaning. Besides, the non-thermal LTP can also be used for surface modification to improve the surface properties, such as the hardness, resistance of corrosion, dielectric properties of materials, without changing the bulk properties because of the low temperature of heavy particles. Recently, the non-thermal LTP have also been used for many biomedical applications such as sterilization, bio-compatibility of materials, and wound healing as a promising enabling technology [5].. 1.2 Non-Thermal Low Temperature Plasmas (LTP) Non-thermal LTP can be further classified using the frequency of applied power as direct current (DC) discharges, dielectric-barrier discharges (DBD), radio frequency (RF) discharges, and microwave frequency discharges. Each type of discharge is introduced next in turn.. 2.

(54) 1.2.1 Direct Current Discharges The DC discharges are the backbone for understanding the discharge characteristics of different types of non-thermal LTP because the well understood behaviors of DC discharges have been used as analogies for different types of non-thermal LTP. A DC discharge is produced by applying a DC voltage across the parallel-plate system in the low-pressure gas as shown in Fig. 1-2 [2]. The sustaining glow discharge is controlled by the emission of secondary electrons generated mainly by ion impact on the cathode. The structure of DC discharge is characterized with distinct regions exhibiting several bright and dark regions along the discharge tube as shown in Fig. 1-3 [6]. The structures such as negative glow, Faraday dark space, and positive column are usually taken as the analogies to analyze different types of discharges. Details of both bright and dark regions are explained by the variation of electron temperature, electron density, and potential. Fig. 1-4 [2] shows a typical I-V characteristic of a DC discharge. The DC discharge can be operated from the Townsend discharge, characterized by very low current and very high voltage, to the thermal arc discharge results from the thermionic emission.. 1.2.2 Dielectric-Barrier Discharges (DBD) The DBD are generated by using two electrodes of which at least one is covered by a dielectric material. The DBD have been known for more than a century since reported by Siemens (in 1857) concentrating on the generation of ozone [7]. The presence of dielectric serves as the current limiter to avoid formation of sparks and current growth. The typical frequency of applied power of DBD is in the range of 1 ~ 100 KHz. Fig. 1-5 [8] shows various types of arrangement of the DBD. The important advantages of the DBD include the simplicity of their arrangements, stable operation conditions, and scalability from small laboratory reactors to large industrial 3.

(55) installations.. 1.2.3 Radio Frequency (RF) Discharges The RF power can interact with plasmas either inductively or capacitively as shown in Fig. 1-6 [2][9]. The frequency of RF discharges ranges from 1 ~ 100 MHz. Both positive and negative charges are retained between electrodes in the RF discharges rather than quenched at the electrodes as that of DC discharges; therefore, the RF discharges require less breakdown voltage as compared to the DC discharges. RF discharges can be sustained with internal and external electrodes, which is important for some cases with corrosive gases or to reduce contamination of the plasmas with the material of the electrodes, whereas DC discharges require the electrodes to be placed inside the reactor because of the generation of secondary electron emission which is an important factor to sustain the DC discharges. In RF plasmas, the energy of ion bombardment on the substrate can be controlled by another external bias, whereas the DC discharge are exposed to the bombardment of high-energy ions that are accelerated at voltages across the cathode fall leading to the damage on the sensitive substrates. The RF discharges are successfully applied to thin film deposition and etching as well as to the sputtering of insulating materials. A matching box is required to adjust the impedance of the plasma reactor for higher efficiency of power transferred to the discharge.. 1.2.4 Microwave Discharges Microwave discharges are sustained by applied power operating in the range of 300 MHz ~ 10 GHz. The excitation of microwave discharges is similar to the excitation of RF discharges. Microwaves are easily absorbed or reflected by most materials and can not be transmitted by cables without significant losses. The applied 4.

(56) microwave transmits in the waveguide with proper design. The microwave discharges generate usually higher density of plasma than that of RF discharges. However, the wavelength of applied electric power shortens as the frequency increases, which causes to the standing wave problem leading to the serious uniformity issue in practical applications. Many applications of non-thermal LTP, especially in semiconductor fabrication industry, are operated at low pressure for reasons such as to reduce the contaminant of impurities or to eliminate the disturbance of impurities during operation (for example, etching and deposition). Recently, the atmospheric-pressure (AP) non-thermal LTP attract tremendous attention because they do not require the expensive vacuum system for applications, for example, ozone production, pollution control, surface modification, sterilization, and wound healing [5][7] which are less sensitive to the impurities of background gases. Among the sources of non-thermal LTP previously introduced, APDBDJ have been used popularly because of: 1) low cost due to no need of the expensive vacuum pumping system,. 2) its simple implementation and stable. operation conditions, and 3) being a stand-alone module. In this thesis, we are interested in investigating helium APDBDJ numerically and its importance is introduced next.. 1.3 Importance of Helium Discharges Many gases (argon, helium, nitrogen, oxygen, air, etc.) have been used as the working gas to generate discharges. Among these working gases, helium is a commonly selected and studied as the working gas for the APDBDJ because it can be operated in a wide stable operating window at atmospheric pressure [10]. Moreover, helium discharges with addition of other gases for specific purposes are not uncommon for applications. For example, oxygen (up to ~ 2%) are widely added in 5.

(57) the helium discharges to generate the oxygen reactive species such as O, O*, and O2* which are believed that those species are important for applications like etching [11], surface modification [12][13], surface cleaning [14], and sterilization [15], among others. Therefore, this thesis focuses on the fundamental study of helium APDBDJ.. 1.4 Numerical Simulation of Gas Discharges Despite the gas discharges are the promising technology adopted in many modern applications, the underlying physics and chemistry of many practical discharges have not been understood thoroughly due to their complexity. The design of these discharge devices and process controls of manufacturing mainly depends on time-consuming and expensive trial-and-error approach, which is economically inefficient. It may be limited and difficult, if possible, to comprehend completely these complex physics and chemistry associated with the discharges through the experimental methods. Therefore, numerical simulation provides an alternative approach in revealing the complex physics and chemistry of discharges. Two major approaches, Particle-in-Cell with Monte-Carlo collision (PIC/MCC) and plasma fluid modeling (PFM hereafter), have been widely used for low-temperature discharge simulations. Although the PIC/MCC approach solves the Boltzmann equation directly and statistically, it is very time-consuming for higher pressure condition because of large amount of pseudo particles is needed for obtaining an accurate solution. PFM, which assumes the plasma as a continuum, is often employed to model gas discharges and requires less computational time than the PIC/MCC approach, should the pressure be not too low. Nevertheless, the applications accompanied with large-scale (computational) domain and/or complex chemistry (species and reactions), which is not uncommon in practice, could lead to unacceptable computational time even using PFM. Acceleration of the fluid modeling 6.

(58) is thus strongly required in simulating realistic gas discharges. Fortunately, this difficulty can be resolved by the rapid advance in computer hardware and the development on the parallel computing. The implementation of parallel computing is capable of reducing the runtime dramatically as discussed in this thesis.. 1.5 Literature Survey 1.5.1 Simulation of Helium Discharges with Impurities APDBDJ with pure helium has been studied for a long time. However, it is very difficult to conduct the experiment with 100 % pure helium because of the limit of industrial production. Different grades of helium gas may be used to conduct the experiments by different research groups. It was shown that discharge current calculated from pure helium data does not quantitatively agree with experimental results [16]. Although the impurity level is typically less than 0.01 %, it was reported that the impurity plays an important role in “ pure” helium discharges [16]-[20]. Effect of impurity was generally modeled with nitrogen because of good agreement with experimental results under reasonable level of impurity (~100 ppm). The simulation indicates that metastable helium atoms generated during a discharge breakdown contribute significantly to the pre-ionization of the gas before the next breakdown through Penning ionization of nitrogen impurity. Although the importance of impurity in the helium APDBDJ has been noticed, there is no report focusing on the study of realistic composition of helium impurities. The existence of impurities change the discharge species chemistry significantly as presented later (Chapter 3) though the impurity level is as low as several to tens ppm for high grade helium gas. Practically, it is important to realize the detailed species of discharge chemistry for deducing the correct conclusion from the experimental observation of applications; therefore, the difference between helium discharge with 7.

(59) and without the realistic impurities is discussed in this thesis.. 1.5.2 Parallel Computing of Fluid Modeling One-dimensional plasma fluid model is a good tool to investigate the temporal and spatial variations of plasma physics and chemistry for a homogeneous discharge. However, one-dimensional simulations are not capable of revealing the species distributions associated with the flow dynamics in the discharge and afterglow region which determines the effect of discharge treatment for applications. This difficulty can be resolved by using two-dimensional simulations for most of the applications if the discharge behavior can be simplified as two-dimensional planar or axisymmetric simulations. Fluid modeling generally requires less computational time as compared to PIC/MCC approach. Nevertheless, it is still an issue for large-scale two-dimensional problems with many species and complex chemistry, which could lead to unacceptable computational time. Fortunately, the development of parallel computing has been proved to be able to reduce effectively the computational runtime dramatically for simulations. Recently, a representative plasma simulation package Plasimo using the finite-volume method, developed by van Dijk et al. [21], has been demonstrated as a parallel version using symmetric multi-processing (SMP) with OpenMP protocol. However, OpenMP allows only data to be shared within a single node with multi-processor, which limits the problem size or chemistry complexity of simulations. To speed up the computation of large-scale plasma fluid modeling, parallel computing using very popular distributed memory machines with message passing interface (MPI) is required. In the past, there are very few studies focusing on parallel implementation of low-temperature plasma fluid modeling on distributed memory machines, albeit the 8.

(60) importance of reduced computational time cannot be overemphasized for practical applications. Among the very few, the parallel fully implicit Newton-Krylov-Schwarz (NKS) algorithm was employed to solve the coupled large sparse, algebraic nonlinear system of the discrete governing equations of fluid modeling derived from the fully implicit scheme [22]-[23]. Although the speedup of parallel computing is scalable up to hundreds of processors, the overall computational time is too large for realistic large-scale multidimensional problems. This obstacle could be overcome by the so-called semi-implicit method, which solves the fluid modeling equations independently with proper linearization of the source terms of the Poisson equation [24] and the electron energy density equation [25]. In this approach, the coupled nonlinear system of plasma fluid modeling equations become linear and decoupled so that they can be solved sequentially. It was shown that much larger time step could be employed to greatly shorten the computational time in sequential implementation. Thus, one of the major objectives of this thesis is to develop a parallel two-dimensional plasma fluid modeling code using the cell-centered finite-volume method with the semi-implicit approach. The resulting linear systems of discretized equations are solved by the parallel generalized minimal residual method (GMRES) [26] in conjunction with the parallel additive Schwarz method (ASM) [27] as the preconditioner to accelerate its convergence. The Schwarz type methods have been proved to be theoretically optimal for many types of problems, and practically powerful for solving large problems on computers with thousands of processors. The preconditioner is decomposed into several sub-domains by domain decomposition for parallel computing. The computational time could be dramatically reduced with the combination of preconditioning and linear matrix solvers for various modeling equations as presented in Chapter 4.. 9.

(61) 1.5.3 Parallel Hybrid Numerical Algorithm for Simulating APPJ Gas discharge simulation considering diffusion without convection for heavy species (i.e., ignoring fluid dynamics effect) is considered to be valid at the low-pressure condition. However, fluid dynamics is expected to have a strong impact on the gas discharge at high pressure condition such as the atmospheric-pressure plasma jet (APPJ). Thus, it is necessary to properly model and integrate the gas flow and gas discharge simultaneously for a better understanding the APPJ. However, it is found that only very few studies in the literature have focused on this subject [21][28][29]. This is mainly because that a complete simulation of the APPJ coupling the fluid dynamics and gas discharge often takes from weeks up to months of runtime. Specifically speaking, simulation of gas discharge often takes generally about 90% of the overall runtime for the simulation of the APPJ because of the very small time step limited by the very light electron. In other words, the bottleneck for speeding up the APPJ simulation is to shorten the runtime consumed by the modeling of gas discharge. For the PFM employed in this study, the evolution of gas discharge is modeled by the self-consistent solution of the Poisson equation, the charged and neutral species continuity equations, and the electron energy density equation. It is known that electrons respond extremely fast to the temporal variation of the electric field, leading to very large transport properties (i.e., mobility and diffusivity), whereas ions respond relatively slow to the temporal variation of the electric field. The neutrals transported by diffusion are even slower if compared with the drift of charged species induced by the electric field. The time step size used in solving the PFM is generally constrained by the electron motion, in the order of 10-10 seconds, and must be small enough to resolve the electron dynamics for a faithful simulation of gas discharge. This leads to the possibility of neglecting the transport of heavy species as compared to chemical 10.

(62) reaction at each time step without losing numerical accuracy of the simulation. Realization of this concept can further shorten the runtime for solving the PFM by using the temporal multi-scale method (TMSM) as presented in Chapter 5.. 1.6 Objectives and Organization of This Thesis The major objective of this thesis is to develop an efficient hybrid numerical algorithm to couple the PFM and GFM utilizing parallel computing for simulating the helium APDBDJ with impurities. To achieve this objective, several specific goals are summarized as follows: 1. To study the impact on chemistry of the helium APDBDJ with and without impurities by one-dimensional fluid modeling; 2. To develop a parallel semi-implicit two-dimensional plasma fluid modeling code using finite-volume method; 3. To develop a hybrid numerical algorithm to simulating gas flow and gas discharge of the APPJ. The thesis is organized as follows: Chapter 2 describes the modeling and numerical method developed in this thesis including the governing equations of PFM and GFM, discretization of equations, the hybrid algorithm for coupling the PFM and GFM solvers, and the development of TMSM for PFM. Chapter 3 studies the helium APDBDJ with and without considering impurities simulated by one-dimensional PFM. The distribution of species of helium APDBDJ with and without considering impurities are compared, and the reaction paths for generating dominant species are presented and discussed. The effect of different levels of water vapor as impurities is also presented and discussed at the end of this chapter. Chapter 4 investigates the parallel performance of the developed PFM code 11.

(63) through a small and a large problem with the configuration of numerical domain and size equivalent to the experimental setup performed by our group. The detailed performance of different types of equations is discussed, and a demonstration of helium APDBDJ without including the flow dynamics is given at the end of this chapter. Chapter 5 addresses on the development of an efficient hybrid numerical algorithm to couple the developed PFM and GFM solvers for simulating a realistic parallel-plate APDBDJ. It is shown that the runtime of PFM can be reduced tremendously for the problem tested using the proposed TMSM. The cycle-averaged species distributions at 200th cycle are presented which shows that it is important to consider fluid dynamics for the APPJ simulation. Chapter 6 concludes by summarizing the major findings found in this thesis and outlining the recommendations for the future work. Fig. 1-7 shows the research structure of the thesis.. 12.

(64) Chapter 2 Modeling and Numerical Methods 2.1 Governing Equations of Fluid Modeling Fluid modeling is suitable for low-temperature plasmas in a wide range of pressures (from low pressure to atmospheric pressure). Generally, there are two types of approximations used in the fluid modeling: 1) local field approximation (LFA) and 2) local-mean-energy approximation (LMEA). The former assumes the locally absorbed electric power is fully balanced by the power dissipated through ionization, while the latter solves the electron energy density equation to obtain the electron temperature which is related to the evaluation of reaction rate constants and other transport properties associated with electrons. LMEA has been shown to be more accurate than LFA in fluid modeling of low-pressure gas discharges [30]. For wider future applications, the LMEA is adopted in the current study to consider non-local effect of electron energy distribution that LFA generally lacks. The governing equations of fluid modeling with LMEA include the continuity equations of each species, the Poisson equation for calculating the electric field, and the electron energy density equation for evaluating the electron temperature subject to the non-local effect of electron kinetic energy. The general continuity equation for ion species can be written as ∂n p ∂t. r r + ∇ ⋅Γ p=. ∑S rp. i =1. pi. p=1,…,K. (2-1). where n p is the number density of ion species p, K is the number of ion species, rp is the number of reaction channels that involve the creation and destruction of ion r species p and Γ p is the particle flux that is expressed, based on the drift-diffusion. approximation, as. 13.

(65) r r r Γ p = sign(q p ) µ p n p E − D p∇n p r r E = −∇φ. (2-2) (2-3). r where q p , E , φ , µ p , and D p are the ion charge, the electric field, the electric potential, the ion mobility, and the ion diffusivity respectively. Note that the form of the source term S pi can be modified according to the modeled reactions describing how the ion species p is generated or destroyed in reaction channel i. The continuity equation for electron species can be written as. ∂ne r r + ∇ ⋅Γ e= ∂t. ∑S re. i =1. (2-4). ei. where ne is the number density of electrons, re is the number of reaction channels r that involve the creation and destruction of electrons and Γ e is the corresponding. particle flux that is expressed, based on the drift-diffusion approximation, as r r r Γ e = − µe ne E − De ∇ne. (2-5). where µe and De are the electron mobility and electron diffusivity, respectively. These two transport coefficients can be readily obtained as a function of the electron temperature from the solution of a publicly available computer code for the Boltzmann equation, named BOLSIG+ [31]. Similar to S Pi , the form of S ei can also be modified according to the modeled reactions that generate or destroy the electron in reaction channel i. The continuity equation for neutral species can be written as. ∂nuc r r + ∇ ⋅ Γ uc = ∂t. ∑S ruc. i =1. uci. uc=1,…,L. (2-6). where nuc is the number density of uncharged neutral species uc, L is the number of neutral species, ruc is the number of reaction channels that involve the generation. 14.

(66) r and destruction of uncharged species uc and Γuc is the corresponding particle flux. which can be expressed as r r Γuc = − Duc∇nuc. (2-7). where Duc is the diffusivity of neutral species. It is noted that the convective effect is neglected in this study. Similarly, the form of Suci can also be modified according to the modeled reactions that generate or destroy the neutral species in reaction channel i. The electron energy density equation can be expressed as. r r ∂nε r r + ∇ ⋅ Γ nε = − eΓ e ⋅ E − ∂t. where.  3  nε  = ne k BTe   2 . ∑ε k n n sc. i =1. i i i e. − 3. me ne k B vm ( Te − Tg ) M. is the electron energy density, Te. (2-8). is the electron. temperature, ε i and ki are the energy loss and rate constant for the ith inelastic electron collision respectively, ni is the number density of species related to the ith inelastic electron collision, Sc is the number of reaction channels of inelastic electron collision, kB is the Boltzmann constant, ν m is the momentum exchange collision frequency between the electron (mass me ) and the background neutral (mass M), Tg r is the background gas temperature. Γ nε is the corresponding electron energy density. flux and can be expressed as r r r 5 5 Γ nε = k BTe Γ e − De ne∇ ( k BTe ) 2 2. (2-9). The second term on the right-hand side of Eq. (2-8) represents the sum of the energy losses of the electrons due to inelastic collision with other species. The last term on the right-hand side of Eq. (2-8) can be ignored for low-pressure gas discharges, while it is important for medium-to-atmospheric pressure discharges. The Poisson equation for electrostatic potential can be expressed as 15.

(67) K r r ∇ ⋅ (ε∇φ ) = −∑ ( qn )i. (2-10). i =1. where φ is the potential, K is the total number of charged species and the permittivity ε , is a function of position, is written as. ε = ε rε 0. (2-11). where ε 0 is the vacuum permittivity, and ε r is the relative permittivity of each region. Several different regions including discharge, dielectric materials (such as alumina, substrate, and Teflon), and conductors are considered simultaneously in this study. All of these regions are meshed and solved to obtain the electrostatic potential distribution by using the Poisson equation.. 2.2 Boundary Conditions The flux-type boundary conditions for the ions, electrons, and neutral species are employed on the solid surfaces (dielectric or electrode) as. r r r Γ p = a ⋅ sign(q p ) µ p n p E − Dp ∇n p r r r 1 Γ e = − a ⋅ µene E − De∇ne + ne vth 4 r r Γuc = − Duc∇nuc. (2-12) (2-13) (2-14). r where a = 1 if drift velocity ( sign(q p ) µ p E ) points toward the dielectric surface, and a = 0 otherwise. We assume that the ions and electrons accumulate and the neutral species quench at the dielectric surface in the present study. The thermal velocity of electron is vth =. 8k BTe π me. (2-15). where me is the electron mass. Note that the effect of secondary electron emission is neglected. For all species, the fluxes at the boundaries of computational domain. 16.

(68) (except the dielectric surfaces) are assumed to be zero. The boundary conditions of electron energy density flux at the dielectric surfaces are r r Γ nε = 2k BTe Γe. (2-16). For the Poisson equation, the potentials of powered and grounded electrode are assigned with applied voltage and zero potential respectively. Neumann boundary conditions with zero gradients are applied to the other boundaries of the computational domain for the Poisson equation.. 2.3 Implementation of Semi-Implicit Schemes It was reported that explicit evaluation of the source term of the Poisson equation leads to a very small time step due to the restriction of dielectric relaxation time [24]. The so-called semi-implicit treatment is thus applied on the source term of the Poisson equation to expand the time step by a Taylor’s series expansion (TSE) in time. With some derivations based on a TSE in time and approximations, the Poisson equation, Eq. (2-10), can be rewritten as K r   r  ∇ ⋅   ε + ∆ t ∑ ( q µ n )i  ∇ φ  = − i =1    . ∑q n K. i =1. i i. (2-17). Note the number densities and mobilities of semi-implicit term in equation (2-17) are approximated from the values of previous time level. Similar constraint on time step size can be found on the source term of the electron energy density equation, Eq. (2-8), and the energy source term is linearized by a TSE in electron energy with some approximations for increasing the time step size of the simulation [25]. Thus, the electron energy density equation can be rewritten as. 17.

(69) sc r r ∂nε r r m + ∇ ⋅ Γ nε = − eΓ e ⋅ E − ne ∑εi ki ni − 3 e nekBvm (Te − Tg ) ∂t M i =1 (2-18) r r s  e r ∂Γ e ∂µe e r ∂Γ e ∂De c ∂ki  me ∂vm E⋅ ni + 3 kB −  E⋅ + + ∑εi (Te − Tg ) nε − neε M ∂ε  ne ∂µe ∂ε ne ∂De ∂ε i=1 ∂ε  r r ∂Γe ∂µe ∂Γ e ∂De 3 has been derived where ε = 2 kBTe . The discretization form of + ∂µe ∂ε ∂De ∂ε. (. ). by Hagelaar et al. [25], and the finite difference method is applied to evaluate and. ∂ki ∂ε. ∂vm . Details of the implementation can be found in references [24][25], and are ∂ε. not described here for brevity.. 2.4 Discretization and Numerical Schemes In the present study, the above equations are discretized using the collocated cell-centered finite-volume method [32] as. ∂η Fi +1/ 2, j − Fi −1/ 2, j Gi , j +1/ 2 − Gi , j −1/ 2 + + = Si , j ∂t ∆ xi , j ∆ yi , j. with. Γp    np   Γ    n   e   e η =  nuc  , F =  Γ uc  , G =        Γε    nε   Γφ    0    x  i, j. (2-19). Γp  Sp    S  Γe   e Γuc  , S =  Suc     Γε   Sε   Sφ  Γφ  y   i, j. where the subscripts i and j represent the indices of cell in x- and y-direction respectively. For simplicity of presentation, the rectangular computational domain is assumed and a set of regular grids is considered. ∆x and ∆y are the cell width in xand y-direction respectively. The fluxes in the continuity equations and the electron energy density equation are calculated with the Scharfetter–Gummel (SG) scheme [33]. After the backward Euler method is employed as a time-integrator, the discretized form of continuity equation can be written as nik, +j 1 − nik, j. ∆ t. +. Γik++11. 2. ,j. − Γik− +11. ∆ xi , j. 2. ,j. 18. +. Γ ik,+j1+ 1 − Γik,+j1− 1 2. ∆ yi , j. 2. = Sik, j. (2-20).

(70) with Γ Γ. k +1 i+ 1 , j 2. Γ. k +1 i−1 , j 2. Γ. k +1 i, j + 1 2. k +1 i, j −1 2. = − = −. ). (. (. ).  B −X nik++1,1j − B Xi+ 1 , j nik, +j1  i+ 1 , j 2 2  xi+1, j − xi, j  Dik+ 1 , j 2. k i− 1 , j 2. D. xi, j − xi−1, j. ). (. (. ). Xi + 1 , j =. ). (. (. Xi − 1 , j =. ).  B −X nik, +j+11 − B Xi, j+ 1 nik, +j1  i, j + 1  2 2 yi, j+1 − yi, j . = −.  B −X nik, +j1 − B Xi, j− 1 nik, +j−11  i, j − 1  2 2 yi, j − yi, j−1  . ). (. k i, j−1 2. D. (. ). D. sign(q)µik− 1 , j 2. k i− 1 , j 2. D. 2. = −. D. 2. k i+ 1 , j 2. 2. B −X nik, +j 1 − B Xi− 1 , j nik−+1,1j  i− 1 , j 2 2  . k i, j + 1 2. sign(q)µik+ 1 , j. Xi, j + 1 =. sign(q)µik, j + 1. Xi, j − 1 =. (φ. 2. (φ. D. sign(q)µik, j − 1 k i, j− 1 2. D. 2. (φ. 2. k i, j+ 1 2. 2. (φ. i +1, j. i, j. −φ. i, j +1. i, j. −φ. i, j. i−1, j. −φ. −φ. i, j. i , j −1. ) ). k +1. k +1. ) ). k +1. k +1. where the superscripts k and k+1 represent properties of the previous and current time X . e −1. levels respectively and the Bernoulli function B ( X ) =. X. Similarly, the discretized form of electron energy density equation (2-8) can be written as nεk(+i1, j ) − nεk(i , j ) ∆ t. Γεk +i1+ 1. (. +. with. 2. ,j. ). − Γεk +i1− 1. ∆ xi , j. (. 2. ,j. ) +. ). (. Γεk +i1, j + 1 − Γεk +i1, j −. (. (. 2. ). ∆ yi , j. (. 1. B − X nεk(+i1+1, j ) − B X i + 1 , j nεk(+i1, j )  i+ 1 , j   2 2. k. B − X nεk(+i1, j ) − B X i − 1 , j nεk(+i1−1, j )  i− 1 , j   2 2. k. B − X nεk(+i1, j +1) − B X i , j + 1 nεk(+i1, j )  i, j + 1   2 2. 5 Di − 12, j Γ ε i− 1 , j = − ( 2 ) 3 xi , j − xi −1, j k +1. 5 Di , j + 12 Γ ε i, j + 1 = − ( 2) 3 yi , j +1 − yi , j k +1. (. ). (. ). (. ). (. (. (. ). ). ). = Sεk(i , j ). (2-21). ). k. 5 Di + 12, j Γ ε i+ 1 , j = − ( 2 ) 3 xi +1, j − xi , j k +1. 2. ). 5 Di , j − 12  Γ ε i, j − 1 = − B − X i , j − 1 nεk(+i1, j ) − B X i , j − 1 nεk(+i1, j −1)   ( 2) 2 2 3 yi , j − yi , j −1   k. k +1. S. k ε (i , j ). (. v v Se m = −eΓe ⋅ E − ∑ ε m kmk (i , j ) nmk ( i , j ) nek(i , j ) + 3 neν m Tek(i , j ) − Tg (i , j ) M m =1. (. = −e Γ. k e , x (i , j ). ⋅E. k x(i , j ). +Γ. k e, y(i , j ). ⋅E. k y(i , j ). ) − ∑ε Se. m =1. ). k k k m m ( i , j ) m ( i , j ) e( i , j ). k. n. n. +3. The Poisson equation (2-10) is discretized in a similar method as 19. (. m k ( neν m ) Tek(i, j ) − Tg (i, j ) M. ).

(71) 1 ∆ xi, j 1 ∆ yi, j = − . . . . εi', jεi'+1, j. . ' '  εi, j ∆ hx(i+1, j) + εi+1, j ∆ hx(i, j). εi', jεi', j +1. . ' '  εi, j ∆ hy( i, j +1) + εi, j +1∆ hy(i, j). ∑( qn)l  m. l =1. . . (φ. (φ. k +1 i +1, j. k +1 i , j +1. −φ. −φ. k +1 i, j. k +1 i, j. )−ε. )−ε. εi'−1, jεi', j. ' i −1, j. ∆ hx(i, j) + εi', j ∆ hx( i−1, j ) εi', j−1εi', j. ∆ hy( i, j ) + εi', j ∆ hy( i, j−1). ' i , j −1. (φ. (φ. k +1 i, j. k +1 i, j. −φ. −φ. k +1 i −1, j. k +1 i , j −1. . ) . . ) . +. (2-22). k. where ∆hx ( i , j ) and ∆hy (i , j ) represents the half cell width of cell (i,j) in the x- and ydirection respectively. Note the effective local permittivity is defined as. ε. ' i, j. . = ε i, j + ∆ t  . ∑( q µ m. l =1. n. ). . i, j i, j l. . k. (2-23). where the semi-implicit treatment is included.. 2.5 Parallel Implementation of Fluid Modeling At each time step, the resulting algebraic linear systems are solved equation by equation using parallel preconditioned Krylov subspace method provided by PETSc library [34] through domain decomposition technique on top of the MPI protocol. Fig. 2-1 shows the proposed flowchart of simulation. After the evaluation of transport properties and rate constants of reaction channels, the discretized governing equations are solved sequentially with acceptable time step size benefiting from the use of semi-implicit scheme. The computational domain is decomposed with vertex-based partition [35] into several horizontal (or vertical) sub-domains along the y- (or x-) direction. In our implementation, each sub-domain is assigned to a single processor. Such partition does not distinguish different types of physical regions such as electrodes, dielectric materials, and discharge region. Hence, the sub-domain of each processor may or may not contain a region with multi-physics. The Poisson equation is an elliptic partial differential equation (PDE) while the 20.

(72) continuity. equations. and. the. electron. energy. density. equation. are. convection-diffusion-reaction equations either parabolic or hyperbolic types of PDE depending on the Peclet number (ratio of drift to diffusion fluxes). The continuity equations of species can be further classified into charged (such as electron and ions) and neutral species. The continuity equation of charged species consists of the mobility, the diffusivity, and the local distribution of electric field, which are varied both temporally and spatially. The coefficients of the matrices for these continuity equations need to be updated at each time step. It follows that the corresponding preconditioners of the continuity equation of charged species need to be reconstructed at each time step. On the other hand, the continuity equations of neutral species are diffusive equations and their diffusivities are treated as constant for most neutral species. Thus, the coefficients of these matrices are unchanged at each time step, leading to a constant preconditioning matrix for neutral species is sufficient.. 2.6 Parallel Gas Flow Model (GFM) The GFM employed in the present study is a two-dimensional planar and axisymmetric flow solver developed in our group [36], which simulates the background gas flow as a continuum by solving a set of governing equations including the continuity, Navier-Stokes (N-S), energy, species transport equations, and the equation of state for ideal gases. The general form of two-dimensional planar governing equations is written in the Cartesian tensor as. ∂ ( ρϕ ) ∂ ∂  + ( ρViϕ ) =  µ ∂t ∂X i ∂X i . ϕ. ∂ϕ   + Sϕ ∂X i . (2-24). where t is the time, X i = ( x , y ) is the position vector, Vi = (u, v ) is the velocity vector, µϕ is an effective diffusion coefficient, Sϕ is the source term, ρ is the. 21.

(73) fluid density, and ϕ = (1, u, v , ht , Yi ) represents for the variables for the mass, momentum,.  ht  ht = . energy,. ∑ Yi hi + i. and. mass. 1 2 2 ∑ V j and hi = 2 j =1. ∫C. fraction. . p ,i. dT . . of. ith. species,. respectively.. is the total enthalpy, where C p,i is the. specific heat capacity of ith species at constant pressure and T is the mixture temperature of the background gas flow. The temperature distribution of solids (e.g., electrode and dielectric material) is obtained by solving the steady state heat conduction equation since we are simulating the experiment when the solids are in their steady states. The steady-state heat conduction equation is written as ∂. . ∂ Xi .  k. ∂ T   + Sth, solid = 0 ∂ Xi . (2-25). where k is the thermal conductivity, Sth,solid is the heat source/sink term of the solid material. Conjugate heat transfer is considered by applying the heat flux continuity at the gas-solid interfaces. The governing equations of the GFM are solved using a cell-centered finite-volume method, with an extended SIMPLE (Semi-Implicit Method for Pressure Linked Equations) scheme. A second-order upwind scheme with linear reconstruction is used to evaluate the inviscid flux across the cell interface. A flux limiter is used to prevent the occurrence of local extrema from being introduced by the data reconstruction. Pressure smoothing is employed to avoid the pressure oscillations on a collocation grid. The use of above numerical approaches allows the GFM solver to simulate both compressible and incompressible flows with a wide range of speeds. The computation performed by the flow solver is also parallelized using domain decomposition approach. Detailed numerical implementation and validations of the. 22.