Numerical Investigation on the Effects of Chemical Reactions on the Discharge Characteristics and Energy Balance of a Nanosecond Repetitive Pulsed Dielectric Barrier Discharge

Numerical investigation on a nanosecond repetitively pulsed dielectric barrier discharge (NS-DBD) in air is a temporal and spatial multi-scale problem involving a large number of species and chemical reactions. To know the effects of the species and chemical reactions on the discharge characteristics and energy balance, a high voltage repetitive plane to plane NS-DBD is numerically studied. Four groups of species and the corresponding chemical reactions are adopted in the investigation. The most complex one has 31 species and 99 chemical reactions that contains all reaction types, in particular, the vibrational-translational relaxation reactions, whereas the simplest one has only 4 species and 4 reactions, which represents the main kinetic processes. The others are in between. The discharge energy reaches to a periodic phase equality state after the second pulse in the repetitive pulses, and the present analysis is focused on the 7th pulse. All the N 2 / O 2 mixture reaction models predict almost the same discharge energies, which are qualitatively similar with that in the simplified 4-species model. The prediction of the discharge energy is determined by the electronic excitation and the energy gain by ions, but the vibrational excitation, negative ions, associative ionization, dissociation of nitrogen and oxygen molecules have very weak effects. The gas heating is determined by the exothermic reactions and the ions. The main processes in the fast and slow gas heating are the energy release of ions and the exothermic reactions, respectively. The negative ions, vibrational excitation, and associative ionization have very weak effects on the gas heating during the high voltage pulse, but they have considerable effects at a larger time scale. The magnitudes of the fast gas heating efficiency ( η G H ) are in the range of 41%∼47% in the N 2 / O 2 mixture reduced kinetic models, but η G H is higher in the 4-reaction model.


Introduction
Dielectric barrier discharge (DBD) actuator in plasma flow control is an active research field in the past few years. The advantages of the DBD actuators are low weight, small size, low energy consumption, without moving parts, fast response, and wide-band control authority. Conventional DBD actuators produce electrodynamic (EHD) flow acceleration (ion wind) of the electric field, which can control the flow separation when the free stream Mach number is less than 0.4. Recently, nanosecond pulsed dielectric barrier discharge (NS-DBD) actuators have been proposed. NS-DBD actuators are driven by nanosecond high voltages repetitively. The discharge is featured by the local ionization at a nanosecond time scale, which induces a fast and slow gas heating at different timescales. The fast gas heating usually occurs at 0.01∼1 µs, which can induce compression waves. However, the slow gas heating occurs at a time scale of about 200 µs. It can introduce wide-band perturbations of temperature, pressure, and density, which can trigger instabilities in the shear flows. NS-DBD plasma actuators can be applied in the high-speed flow controls [1,2].
The air discharge of NS-DBD involves complex reactions, which makes the kinetic numerical modeling of the discharge very complicated. The governing equations of lots of species introduce a very large partial differential equation system. As the reaction rates of different reactions and number densities of different species are much more different, the equation systems are very stiff to solve. Because different chemical reactions of the multi-species contribute much more differently to the energy balance, an appropriate modeling of chemical reactions is very crucial. Kossyi et al. [3] developed a kinetic model of N 2 /O 2 mixtures having 46 species and 395 reactions, which includes electronic excitation, electron impact destruction and ionization, electron attachment and detachment, recombination of electron-ion and ion-ion, associative ionization, ion conversion, and chemical transformations of neutral particles. It is very difficult to adopt all the reactions in a physical modeling of the NS-DBD discharge, due to the very large governing equation system and huge computational cost. Therefore, most of the work adopted different reduced reaction models. Nagaraja et al. [4] adopted a kinetic model that consists 20 species and 176 reactions to do multiscale modeling of a pulsed NS-DBD. Bak and Cappelli [5] developed a 20 species and 49 reactions kinetic model for the low-density plasma and a 19 species and 73 reactions model for the high-density plasma. Bak and Cappelli [6] also adopted an 8 species and 15 reactions kinetic model simulated the plasma formation in a nanosecond pulse plane to plane dielectric barrier discharge at atmospheric pressure and 340 K. Poggie et al. [7] developed a reduced model with 23 species and 50 processes to study the main species and reactions in energy storage and thermalization of an NS-DBD. They [7] also developed a 15-species and 42-process model to simulate a NS-DBD discharge. Zhu et al. adopted two models in plasma simulations, which contain 16-species, 44-process [8], and 13-species, 37-process [9]. Raja et al. [10] simulated an airglow discharge at pressure ∼1 Torr by an 11 species and 21 reactions plasma model. Parent et al. [11] used an 8 species, 28 reactions air chemical model modeled plasmas for computational aerodynamics. The reduced kinetic models neglect some minor species and the corresponding reactions, however, the effects of the neglected species on the preserved major species need to be addressed.
Different models have been proposed to simulate the discharge, which mainly include particle-in-cell (PIC) model [12][13][14][15], Monte Carlo model [16,17], five-moment fluid model [7,18], three-equation drift-diffusion fluid model [4,[6][7][8]10], two-equation drift-diffusion fluid model [19][20][21][22][23][24][25], and some other analytical models according to the fidelity of the physical principles. The PIC and Monte Carlo models based on the principle of particle collisions, typically involve the solution of the Boltzmann equation, which is computation expensive, requiring a long computation time. In a normal NS-DBD plasma discharge under working pressure (>0.5 Torr) [26], the plasma can be assumed as a continuum fluid. It is reasonable to adopt the fluid models, which compromise the fidelity and computational cost.
In the plasma flow control application, the speed of the gas heating and the amount of the discharge energy is transferred to the gas heating are the major concerns. The discharge energy is the total energy of the charged species (electron and ions) gained from the electric field, which is calculated by different methods. In the experiment [27], the discharge power and energy are calculated by the load voltage and the current. In the quasi-one-dimensional analytical model [28], the discharge energy is the sum of energy in the plasma and the energy stored in the capacitor, which is formed by the sheath and the dielectric layers. The gas heating is the energy that transfers from the discharge energy to the thermal energy of the gas, which is also considered by different methods. Nagaraja et al. [4] put all the charged species energy to the fluid. Roy et al. [23] and Che et al. [24] assumed the ions' energy obtained from the electric field fully transfer to the fast gas heating, but a fraction of the electron's energy contributes to the gas heating. Zhu et al. [8] assumed that the electrons contribute to the gas heating through the quenching of the excited neutrals, and the ions' energy is totally transferred into the gas heating instantaneously. Boeuf et al. [20] and Abdollahzadeh et al. [22] predicted the electron energy release to the gas heating through electronic excitation, vibrational excitation, and the elastic and rotational collisions at different timescales. The ions' energy totally contributes to the gas heating instantaneously.
Although the past studies have adopted different chemical reaction models in the NS-DBD investigations, the effects of the chemical reactions on the discharge characteristics, in particular on the discharge energy and gas heating, have not been studied systemically. In addition, there are few numerical studies concerning the vibrational-translational relaxations that are crucial to the slow gas heating. Therefore, the detailed vibrational-translational relaxation reactions are included in the present reaction models. In the present investigation, the objectives are, on the one hand, to figure out the effects of chemical reactions on the discharge characteristics and energy of a plane to plane high-voltage repetitive NS-DBD actuator. On the other hand, foundations are built to simplify chemistry reactions for further researches on the plasma flow control. A three-equation drift-diffusion model is used to reduce the computational cost, the model is adopted in the previous works [29][30][31]. Four sets of the N 2 /O 2 mixture reduced kinetic models (31-species 99-reaction, 28-species 85-reaction, 23-species 50-reaction and 15-species 42-reaction) and a simplified N 2 /O 2 mixture averaged 4-species 4-reaction kinetic model are adopted. The electron-impact reaction-rate coefficients and the electron drift and diffusion coefficients are calculated by the Local-mean-energy approximation (LMEA). Because the initial states, i.e., the species number densities and the electric field, are different in each repetitive pulses, the discharge needs at least two periods to reach a statistical periodic state. The analysis of the discharge characteristics and the effects of chemical reactions is focused on the 7th pulse. The discharge energy is calculated by the current and voltage on the electrodes. The gas heating is predicted by the sum of the heat of the main exothermic reactions and the energy of ions. The conclusions are drawn on the kinetic analysis of the energies in the discharge energy and the gas heating in the different reaction models.

Discharge Model and Numerical Method
The three-equation drift-diffusion fluid model contains the mass-conservation equations, an energy-conservation equation and a Poisson's equation for the voltage. The electron's mass and energy conservation equations are ∂ ∂t where e and ε are the electron and electron-energy, respectively. n e and n ε are the number density and the energy density of the electrons. R e and R ε are the reaction induced source terms. Γ e and Γ ε are the fluxes with drift-diffusion approximation, which are The drift and diffusion coefficients are calculated by using a Boltzmann equation solver Bolsig+ [32] with the electron energy probability function (EEPF) f 0 [33], where N is the total gas number density and m e is the mass of an electron. σ m is the effective momentum transfer cross section. v c is the external gas velocity. R e and R ε are ∑ (α − η) Γ e + k r Πn r and ∑ (α − η) Γ ε ∆ε + k r ∆rΠn r . α and η are the ionization and the attachment coefficients. k r is the reaction rate and n r is the reactants' number density. ∆ε and ∆r are the energy loss of electronic ionization reactions and other electron impact reactions. The heavy species' mass-conservation equations are where ω k is the mass fraction of k. R k is the source term due to the chemical reactions. Γ k are the fluxes where µ k and D k indicate the drift and diffusion coefficients, respectively. z k and M are the charge number and the average molar mass. The drift coefficients of ions are determined by the electric field [34]. The diffusion coefficients of ions are calculated from the generalized Einstein relation, and diffusion coefficients of the neutrals are calculated by the classical gas kinetic theory [35,36]. ρ is the gas density, which is determined by the equation of state for a perfect gas, R is the universal gas constant. p and T are the pressure and temperature, respectively. The electric potential φ induced by the charge is described by the Poisson's equation − ∇ · (ε∇φ) = e ∑ (−n e + z k n k ). (12) e is the elementary charge and ε is the permittivity. ε is 1.0 in the air and a dielectric constant in the dielectric barrier. n k are the number densities of k, which is n k = ω k N A ρ M k , where N A is the Avogadro constant. The electric field E is calculated from the potential, E = −∇φ. The present governing Equations (1), (2), (9) and (12) are solved by adopting a finite element method [30]. In each element, a second-order Lagrange shape function is used, and the time integration is implemented by using a second-order backward differentiation formula. The boundary conditions and the time step are the same as that in the previous work [30].

Air Chemistry
N 2 /O 2 mixture is assumed as the initial neutral air, and the molar fraction of N 2 and O 2 is 4:1. There are 46 species and 395 chemical reactions [3] and the chemical reactions are very complex even for this simplified air. Kinetic analysis conducted by Popov et al. [37] reveals that the electrons impact dissociation of O 2 molecules, the quenching of electronic excited N 2 molecules by O 2 molecules and the quenching of excited O atoms by N 2 molecules mainly contribute to the fast gas heating when the reduced electric field E/N ≤ 200 Td. The electron impact dissociation of N 2 molecules and the processes of the charged species are the major processes when E/N > 400 Td. A reduced chemistry model is considered by containing the main neutral and charged species involved processes. The reduced kinetic model consists of 31 species and 99 reactions, which is denoted as R99 hereafter. The species and the reactions are shown in Table A1 in Appendix A. The abbreviations of the types of reactions are shown in Table A2. The details of the reaction are shown in Table A3.
As the fractions of the contribution of the dissociation reactions are smaller than the corresponding quench and charged species processes on the gas heating in low and high electric field zones [37]. A further reduced kinetic model is considered based on R99 with all the dissociation processes and the associated species excluded. The reduced kinetic reaction model consists of 28 species and 85 reactions and is denoted as R85 hereafter.
Based on the sensitivity analysis conducted by Uddi et al. [18], the electronic excited state molecules The model is the same as the one used by Poggie et al. [7], which consists 23 species and 50 reactions. It is denoted as R50. The details of R50 can be referred from [7].
Vibrational excited nitrogen contributes to the gas heating through the vibration-translation (VT) relaxations at a large timescale. However, it contributes no more than 10% of the total energy released at the relatively high reduced field (E/N > 200 Td) [38]. Therefore, a kinetic model is further reduced based on R50 by excluding all the vibrational excited nitrogen species N 2 (υ = [1][2][3][4][5][6][7][8]. It is denoted as R42. The details of the reactions are also be referred from [7] but excluding all the 8 vibrational excitation reactions. Another simple kinetic model is considered by averaging the N 2 and O 2 species as an "Air" molecule, which incorporates 4-species and 4-reaction mechanism [29][30][31], as shown in Table A1. It is denoted as R4 hereafter. Details of the reactions are shown in Table A4. The model is the same as the one used in the previous works [20,39,40].

Energy Transfer Models
The energy transfer process of a DBD plasma actuator is analyzed from the electric power source to the gas heating, as shown in Figure 1. The discharge energy supplied by the power source transfer to the electrons and ions through the electric field. The energy gained by the electrons dissipates to the electron-energy change and electron impact endothermic reactions. The electron-impact endothermic reactions produce excited neutral species and ions. This part of electron energy transfers to the chemical energy of the resultants. The chemical energy stored in the excited neutrals and ions transfer to the gas heating through the exothermic reactions. As the molar masses of ions are comparable to the one of the neutral species, it is assumed that the energy of ions is totally transferred into the gas heating instantaneously.
The discharge energy is mainly determined by the current and the voltage, which is calculated by multiplying the current and the voltage on the powered electrode. The model is denoted as A_ui and the discharge energy is where T p is the integrating time.
The gas heating power is denoted as P GH and the corresponding efficiency is denoted as η GH . The contribution of ions to P GH is P ion = J ion · E. The contribution of reactions to the gas heating is P R . The predictions of P R and Q R are based on the kinetic analysis of the kinetic models. Q R is calculated by the sum of the energies of the dominant reactions that contribute to gas heating. The power P GH and energy Q GH released for the gas heating are shown in Table 1, where ∆E and r are the reaction energy loss and the reaction rate, respectively.

Name
Power Energy The gas heating efficiency is calculated as η GH = Q GH Q A_ui .

Simulation Setup
A low-pressure plane to plane DBD is adopted from the experiments conducted by Takashima et al. [27] to investigate the effects of kinetic reactions on the discharge characteristics and energy. The referred experiment contains abundant experimental data. The actuator consists of two rectangular copper plate electrodes with each one covered by a 1.75 mm thickness quartz layer as a dielectric barrier, where ε = 4.8. The electrode is 14 mm wide and 65 mm long. The length of the discharge gap between the two dielectric barriers is 1 cm. In the experiment, the discharge was sustained in the dry air. The pressure and temperature of the dry air are 60 Torr and 300 K, respectively. The discharge is powered by a Chemical Physics Technologies (CPT) generator, which generates a bipolar negative-positive wave shape with peaks of −22.9 kV and +18.5 kV, respectively. The pulse duration is 85 ns with 45 ns and 40 ns full widths in the negative and positive pulses, respectively. The pulse repetition rate is 10 kHz. The experimental estimated and the fitted voltage waveforms are shown in Figure 2. Exp denotes the experimental waveform. The Gaussian functions fitting to Exp is denoted as Gaussian. The decrease and oscillations of the voltage after 75 ns in Exp is neglected in Gaussian. Time durations of 10∼33 ns, 33∼45 ns and 45∼55 ns correspond to the fall and rise phases in the negative pulse and ascending phase in the positive pulse. Time duration from 55 ns to 85 ns corresponds to the decaying phase of the positive pulse in Gaussian. The decaying phase after 55 ns in Exp voltage lasts more than 125 ns.
According to the transition conditions [41,42], the mode of the present discharge is considered as a quasi-uniform mode, which is validated by the ICCD camera images in the experiment [27]. Because of the quasi-uniform characteristics, a two-dimensional computation model is adopted. The computational domain is shown in Figure 3   The structured mesh is adopted in the simulation. The mesh in the y-axis is uniform, with a size of 100 µm. The maximum mesh sizes in the x-axis in the discharge gap and the dielectric barriers are 5 µm and 10 µm, respectively. The minimum mesh is located at x = 0 mm and x = 10 mm, and the size is 2 µm. The mesh across the entire domain is generated from the minimum mesh to the maximum mesh by a linear growth function, and the growth rate is 1.3. The total number of mesh elements is 175,464.
Because the discharge is quasi-uniform, quasi-one-dimensional simulation models are used to investigate the effects of chemical reactions to reduce the computational costs, which is the same as the models used by Poggie et al. [7] in a similar simulation. The one dimensional simulation is established in the x-axis, and the grid distributions are the same as the one in the two-dimensional model.
In R99 and R85, the initial number densities of the electrons and the negative ions are 5 × 10 11 m −3 , and the sum of the number densities of positive ions is 1 × 10 12 m −3 . In R50 and R42, the initial electron number density and the sum of ions' number density are all set to be 1 × 10 12 m −3 to keep an initial zero net charge. The initial O 2 mole fraction is taken to be 0.2. The mole fraction of other neutrals is 1 × 10 −10 . The mole fraction of N 2 is calculated by the law of mass conservation. In R4, the initial number densities of the electrons and the negative ion are 1 × 10 12 m −3 . To keep the initial balance of the charge, the positive ion number density is 2 × 10 12 m −3 . The low level of pre-ionization has a negligible influence on the discharge energy [7].

Model Validations
To validate the present numerical methods, the current on the electrode and the discharge energy in the referred experiment [27] are referred for comparisons. The one-dimensional model is validated by comparing the currents in the experiment and the two-dimensional model. Then the results in the one-dimensional models are validated further.
To show the reasonability of the quasi-one-dimensional model, the comparison of the currents on the electrode powered by Exp waveform in R4 in the one-dimensional, the two-dimensional models and that in the experiment is shown in Figure 4, the voltage waveform is also shown in the figure.
Because the Exp waveform is formed by piecewise linear interpolations and the gradient due to time is discontinues, the current oscillates during the pulse and in the post-discharge.  The currents in the one-dimensional and two-dimensional models both agree that the experimental current in the fall phase of the negative pulse and the currents in the two models are nearly the same. However, the current in the one-dimensional model is larger than that in the two-dimensional model, and the current in the one-dimensional model has a closer agreement with the experimental current. The currents in both models agree well with that in the experiment in the ascending and decay phase of the positive pulse, and they are nearly the same. It can be concluded that the one-dimensional model and two dimensional model are both feasible for the prediction of the current in the present plane to plane NS-DBD discharges. Therefore, one-dimensional models are used to study the effects of reactions on the discharge characteristics hereafter.
The past numerical simulations [4,8,20,22] of a repetitive NS-DBD focused on the first pulse. The experimental results [27] show that the energy in every pulse is nearly the same in the whole discharge. However, the energy in the first pulse differs from that in the other pulses [27]. In order to know when the discharge energy can become steady, the discharges of the first seven pulses are performed in all simulations. Q A_ui in R99 in the 1st, 2nd, 3rd, 5th, and 7th pulses of the Gaussian voltage waveform are shown in Figure 5. Q A_ui after the 1st pulse are all translated to 0 ns for clarity. It can be seen that Q A_ui are almost the same as each other in each phase of the pulse after the 2nd pulse. It indicates that Q A_ui reaches to periodic state following the discharge pulse, which is similar to the experimental results [27]. Hereafter, the results of the 7th pulse are adopted in all simulations.
To study the effects of voltage waveforms on the current, the currents on the electrode powered by Exp and Gaussian waveforms in R99 are compared, as shown in Figure 6. The start time of the experiment curve at 180 ns [27] is translated to 0 ns for clarity. The magnitude of the current using the Gaussian waveform is smaller than the experimental data near the first negative peak. In the time duration after 75 ns the result of the Gaussian waveform rapidly reduces to zero without oscillates due to the neglect of voltage oscillations after 75 ns in the Gaussian waveform. The predicted currents using Exp and Gaussian waveforms both agree well with the experimental data during the high voltage pulses. Therefore, the present numerical model predicts good results for the discharge current.
The predicted discharge power P A_ui and energy Q A_ui are also compared with the ones in the referred experiment [27], which are shown in Figure 7a and Figure 7b, respectively. The predicted power under the Exp and Gaussian waveforms also both agree well with experimental data, except that the predicted power under the Exp waveform oscillates during the fall phase in the negative pulse (10∼33 ns). The discharge energy predicted by both waveforms is nearly equivalent during the first negative pulse (10∼45 ns). The peaks of Q A_ui in the Exp and Gaussian waveforms at 33 ns are 3.23 mJ and 3.09 mJ, which are 3% and 7% lower than the experimental value 3.33 mJ. During the positive pulse, the two calculated energies are also nearly the same as each other and both are higher than the experimental results. The peaks of Q A_ui in the Exp and Gaussian waveforms at 55 ns are 3.30 mJ and 3.42 mJ, which are 5% and 9% larger than the experimental value 3.13 mJ. The discharge energy under the Exp waveform oscillates after 75 ns, but the result of the Gaussian waveform rapidly reaches to a constant without any oscillations. Although Q A_ui in the simulations are very different from the experimental value, the agreement of the discharge energy during the pulse is well. It is concluded that the numerical model in this work calculates good discharge energy. Because the discharge power and energy predicted by the Exp and Gaussian waveforms are almost identical to each other during the high voltage pulse and the power after the pulse is very small, the Gaussian waveform is used subsequently for simplification.
To investigate the effects of grid resolution, a finer mesh with half the size in the whole computational domain is used for comparison. Q A_ui predicted by R99 in both resolutions are shown in Figure 8. The predicted results are nearly the same, indicating that the grid convergences. Therefore, the coarse mesh is used subsequently.

Discharge Characteristics and Energy
In this section, the discharge characteristics are discussed in R99, and then the investigation on the energy transfer process is provided.

Discharge Characteristics
In this section, the discharge characteristics of the NS-DBD actuator in the high voltage pulse and post-discharge periods are investigated. The energy is investigated by analyzing the magnitudes, the time scales, the species number densities, the discharge energy, and gas heating.
The spatial-temporal distributions of the number densities of the electrons n e and N + 2 ion n N + 2 are shown in Figure 9a,b, respectively. Before the discharge breakdown at a time of about 23 ns, the electric field in the whole domain is low and nearly distributes uniformly, as shown in Figure 10. As a result, the ionization process is weak in this phase, and the low densities of the electrons and ions are shown. After the breakdown, the number densities of the electrons and the ions increase as a result of the electron avalanche. E/N, n e , and N + 2 at time 25 ns after the breakdown are shown in Figure 11. A sheath forms because of the charge accumulation on the dielectric surface of the negatively biased driven electrode. The sheath is featured by a strong electric field and the violation of quasi-neutrality. The sheath thickness is characterized by the length of the high-electric field area according to the Bohm criterion [43]. The sheath edge is considered at the location where the electric field magnitude is decreased by 99%. The thickness of the sheath is 0.179 mm and the maximum E/N in the sheath is about 7290 Td at time 25 ns. In the sheath, an electric field builds up until an equilibrium of electron and ion fluxes towards the surface is established. Because the drift and diffusion coefficients of electrons are much higher than that of the ions, n e is lower than n N + 2 in the sheath. n e and n N +  Figure 11. The sheath has a length of about 0.1 mm, where E/N is larger than 2000 Td. The sheath is nearly free of electrons. n e and n N + 2 display high values in the edge of the sheath, which reach magnitudes of 6.5 × 10 19 m −3 and 4 × 10 19 m −3 , respectively. According to the research of Poggie et al. [7], the error of the two-term approximation of the Boltzmann equation becomes substantial above E/N = 3000 Td. However, n e is very low in the sheath. Therefore, the present model for ionization at the edge of the sheath is accurate enough to study the discharge characteristics.  To study the characteristics of n e locally and globally, n e at the point x = 0.09 mm (Point "A") at the edge of the sheath and the total n e in the whole domain n etot are shown in Figures 12 and 13, respectively. Because the computational model is one dimensional, the number densities in the other two orthogonal dimensional directions are homogeneous. The total number density n etot per unit area is calculated by integrating the number density n e in the computational domain, the unit of n etot is consequently m −2 . The gradient of n e , n etot , and the voltage are also shown in the figures. n e and n etot both increase rapidly in the fall phase in the negative pulse and the ascending phase in the positive pulse, corresponding to two discharges. In the rise phase of the negative pulse, n e and n etot increase slowly, while n e and n etot start to decrease from the decay phase of the pulse. The time of the maximum dn e /dt and dn etot /dt are about 24 ns and 23 ns, when the breakdown happens. Therefore, the breakdown voltages are predicted as 11.9 kV and 10.3 kV. The breakdown voltage predicted by n e locally is in a closer agreement with the experimental estimation 13 kV [27].

Energy Transfer Process
The characteristics of the energy are analyzed by separating each parts of the energy following the energy transfer path, as shown in Figure 1. Q A_ui is the sum of the energies of the electrons Q electrons and the ions Q ions gain from the electric field, and Q electrons is dissipated to the electron-energy change ∆En e and the energies of the electron impact endothermic reactions Q Ee , Q Ei , and Q Ve .
First, because Q Ee , Q Ei , Q Ve , and Q ions are determined by the number densities of the corresponding species, the total number densities of N 2 (A 3 Σ), N 2 (B 3 Π), N 2 (C 3 Π), N 2 (a 1 Σ), the ions N + 2 , O + 2 , and the vibrational excited species N 2 (υ = 1-8) are shown in Figures 14-16, respectively. As the number densities are adopted from the 7th pulse, the initial number densities of the species are not same to that in the initial values. The number densities of the nitrogen excited molecules increase mainly in the two discharges during the pulse, and the increments in the two discharges are nearly the same. The total number density of N 2 (A 3 Σ) increases up to 4.9 × 10 16 m −2 , which is about 4 times larger than the maximum density of N 2 (B 3 Π) and 6 times larger than the maximum densities of N 2 (C 3 Π) and N 2 (a 1 Σ). After the pulse, N 2 (C 3 Π) and N 2 (a 1 Σ) decay rapidly with nearly the same rate, however, N 2 (A 3 Σ) and N 2 (B 3 Π) decay more slowly. The number density of N 2 (A 3 Σ) decreases in about 10 4 ns, which correlates with the time scale of the energy release of the quench of the nitrogen excited molecules, as shown in Figure 21a. The number densities of N + 2 and O + 2 also increase in the two discharges, and the increments in the two discharges are almost the same. The maximum density of O + 2 is about 2 times larger than that of N + 2 . The time scales of the decay of N + 2 and O + 2 are about 500 ns and 8.7 × 10 5 ns, which also correlates to the time scale of the change of the chemical energy stored in the ions, as shown in Figure 20. The number densities of N 2 (υ = 1-8) also increase mainly in the two discharges during the pulse, which is the same as that of the excited species and the ions. The number densities of the vibrational excited nitrogen decrease with the increase of the vibrational quantum number, and the number density of N 2 (υ = 1) is about two orders larger than that of N 2 (υ = 8) after the pulse, which are 1.1 × 10 17 m −2 and 1.3 × 10 15 m −2 , respectively. Because the VT relaxation rates are low, the decay of N 2 (υ = 1-8) starts at a time scale of around 1 ms, which is close to the time scale 0.5 ms in the experiments conducted by Montello and the co-workers [44]. The time scale correlates with the time scale of the energy release of N 2 (υ = 1-8), as shown in Figure 21a, which corresponds to the slow gas heating. According to the analysis of the number densities, the comparison of Q Ee , Q Ei , and Q Ve is shown in Figure 17, and ∆En e , Q ions are shown in Figure 18.  The magnitudes of Q Ee , Q Ei , and Q Ve all increase in the two discharges during the pulse, which is the same as that in the change of the number densities, and they keep nearly constant after the discharges. Since the number densities of nitrogen excited species are higher than that of the ions, the magnitude of Q Ee is higher than that of Q Ei , and they reach the magnitudes up to 0.85 mJ and 0.09 mJ. As the number densities of N 2 (υ = 1-8) are lower than the nitrogen excited species, the maximum magnitude of Q Ve is about 0.1 mJ. ∆En e mainly increases to 0.015 mJ during the pulse in the discharges and decays rapidly after the discharge. Q ions increases to 0.55 mJ in the two discharges and keeps almost constant after the discharges. Therefore, Q A_ui is mainly transferred to Ee and the energy of the ions. However, the energy transfer to Ei, Ve, and ∆En e is neglectable.
Second, according to the energy transfer path, the electron impact endothermic reactions transfer a part of energy from the electrons to the chemical energy stored in the heavy species, which includes N 2 (υ = 1-8), the electronic excited state species (Ee species), the ground excited state species (Ges species) and the ions. Since the chemical energy in the heavy species correlates to the corresponding number densities, the total number densities of N 2 (υ = 1-8), the Ee species and the ions are shown in Figures 14-16, respectively. The total number densities of the Ges species are shown in Figure 19. It can be seen from the figure that O 3 and O are the dominate species. O is primarily formed by the electron impact dissociation (8) and the quenching of electronically excited nitrogen (33), (39), (44) in Table A3, and O 3 is primarily formed by the Ges reactions (31)  The comparison of the chemical energy stored in N 2 (υ = 1-8), the Ee species, the Ges species, and the ions is shown in Figure 20. The chemical energy stored in the heavy species mainly increases during the pulse and the energies stored in Ee and Ges are larger than the ones in the other species. The energy stored in Ee and Ges decreases due to the quench reactions and Ges reactions, respectively, which are also closely related to the change of the corresponding species number densities. The energy stored in Ee decreases in about 3000 ns after the pulse, while the energy stored in Ges still increases slowly after the pulse and decays from 10 5 ns. Therefore, the time scale of the energy release of the quench of electronically excited reactions is much lower than the ones of the other reactions.
Third, to investigate the energies of the exothermic reactions, the comparisons of the energies of Q_N 2 , Q_N, Q_O 2 , Q_O, VT relaxations, Dis_N 2 , Dis_O 2 , Re_ei, Re_ii, Pic, Ges, and Ai are shown in Figure 21a-c, respectively. During the high voltage pulse, Q_O, Q_N 2 , and Dis_O 2 are the main exothermic reactions, which is consistent with the conclusions of Popov et al. [37], when E/N ≤ 200 Td. The energy of other reactions is at least one order of magnitude smaller and can be neglected. The increase of Q_N 2 mainly happens at two time scales 85 ns and 10 4 ns, which correspond to the quench of N 2 (C 3 Π), N 2 (a 1 Σ) and N 2 (A 3 Σ), N 2 (B 3 Π), respectively, as can be seen in Figure 14. VT increases from the timescale 1 ms, which is mainly the result of the quench of N 2 (υ = 1 − 2), as shown in Figure 16a. It is mainly responsible for the slow gas heating. Q_O 2 and Q_N release slowly after the pulse, which starts at the time scale 1 µs and 1 ms, respectively. At a long time after the pulse, Q_O, Q_N 2 , and Dis_O 2 are still the dominant exothermic reactions. Q_O 2 and Q_N are comparable to that of the dominant reactions, which also contribute to the slow gas heating. However, Popov et al. [37] found that the gas heating of the dissociations of the nitrogen molecules and the processes involving the charged particles are the main precesses when E/N > 400 Td.  To study the energies of different reactions at different electric field regions in detail, the distributions of the number densities of the typical species that involved in the entomic reactions at time 35 ns during the pulse are shown in Figure 22a,b, respectively. In the low E/N (<200 Td) region, O, N 2 A 3 Σ and O 3 P are the dominant species, which are associated with the Q_O, Q_N 2 and Dis_O 2 reactions, as shown in Figure 23a. In most parts of the sheath, where E/N is higher than >400 Td, Dis_N 2 and Pic are the dominant exothermic reactions, as shown in Figure 23a. However, the dominant species are O, N, and O 3 P there. The disparity is the result of the different reaction heats of the corresponding reactions. At extremely high E/N (>12,300 Td) in the sheath, all the power densities are very small, and the power density of Ges reactions become the largest one. However, Popov et al. [37] revealed that the Dis_N 2 and the processes involving charged particles are still dominant. The difference may be caused by the limitations of the two-term approximation of the Boltzmann equation and the drift-diffusion approximation at a very high E/N region [33]. It can be concluded that the energies of the dominant reactions are determined by the reactions power densities at low E/N region, as the length of the high E/N sheath is small and the power densities are quite low there.
Finally, according to the energies of the dominant reactions, the total gas heating P GH , Q GH in the whole domain and the efficiency η GH in R99 are calculated and shown in Figure 24. The gas heating releases mainly in two time scales 85 ns and 7000 ns, which are corresponding to the fast and slow gas heating. The power P GH in the fast and slow gas heating processes are higher than 10 5 W and lower than 10 3 W, respectively. The comparison of the contributions of the reaction heat and the ions are shown in Figure 25. According to the magnitudes of P R and P ion , the fast gas heating is at timescale 85 ns, which is dominated by P ion . The slow gas heating is at timescale 7000 ns, which is dominated by P R . It can be seen that the fraction of the contribution of ions is larger than the one of the reactions in the fast gas heating. However, the contributions of the ions are lower than that of the reactions for the slow gas heating at a long time scale. The slow gas heating is the result of the exothermic reactions, which are mainly the quench of the vibrational excited nitrogen molecules.    E n e r g y ( m J ) Figure 25. Comparison of the powers and energies of reaction heat P R , Q R and the ions P ion , Q ion in the gas heating. The corresponding fast and slow energy release times are also labeled.

Effects of Chemical Reactions
In this section, the comparison of the discharge energy and the gas heating predicted by the kinetic models are conducted, the energies and the densities of the preserved dominant processes are also compared to study the effects of the reductions of the kinetic models.
To investigate the effects of chemical reactions on the discharge characteristics and energy, the energy transfer of the electrons and the ions are analyzed. The comparison of P A_ui and Q A_ui predicted by using all the kinetic models is shown in Figure 26a Figure 26. Comparison of the P A_ui (a) and Q A_ui (b) predicted by all the kinetic models and the ones in the experiment [27].
The N 2 /O 2 mixtures reaction models R99, R85, R50, and R42 predict almost the same P A_ui , which are qualitatively similar to that of R4. However, the positive peaks in R4 are lower than that in the other models, which result in a lower Q A_ui in R4 than that in the experiment and the other models. Therefore, the reductions from R99 to R42 have very weak effects on the prediction of the discharge energy. During the pulse, R4 under-predicts the discharge energy, but R99, R85, R50, and R42 over-predict the discharge energy, which agrees better with the experiment. Because P A_ui is calculated by the current on the electrode, the currents in the N 2 /O 2 mixtures reaction models are also nearly the same, which are also different from that in R4.
According to the energy transfer path, the energies of the electron impact exothermic reactions, the electron-energy change and the energy gained by ions in all the kinetic models are compared, which are shown in Table 2. As the energies are qualitatively similar, the values at time 35 ns in the pulse and 85 ns after the pulse are compared quantitatively. The energies of Ei increase from R99 to R42 and the maximum increments at 35 ns and 85 ns are about 29% and 18%. However, because Ee and Ve are not considered in R4, the energies are mainly consumed by Ei. The energy of Ei is much higher than that in the N 2 /O 2 mixtures reaction models. The energies of Ee from R99 to R42 also increase, the maximum increments at 35 ns and 85 ns are 12% and 4%. Similarly, the increments of Ve are 15% and 13%. Compared with the difference of the energies in the N 2 /O 2 mixtures reaction models and R4, the energies of the electron impact exothermic reactions in R99, R85, R50, and R42 can be considered nearly the same. Similar conclusions can be made as to the energies of ∆En e and Q ion . Table 2. Comparison of the energies of Ee, Ei, Ve, the electron-energy change ∆En e , and the energies of the ions Q ions in all the kinetic models.

Ee (mJ)
Ei ( The main exothermic reactions are Q_O, Q_N 2 , Dis_N 2 , and Dis_O 2 . Because Q_O and Q_N 2 are all included in the N 2 /O 2 mixtures reaction models, to investigate the effects of the neglected reactions on the preserved reactions, the comparison of the energies of Q_O, Q_N 2 in the N 2 /O 2 mixtures reaction models is shown in Figure 27a,b, respectively. In all the models, the energies of Q_N 2 are nearly the same. Because Dis_N 2 and Dis_O 2 are reduced in R85, R50, and R42, O( 1 S) and O( 3 P) are also reduced there, the energy of Q_O in R99 is higher than that in the other models. However, the energy of Q_N 2 is not affected a lot by the dissociation reactions. Although the negative ion, the associative ionization, the electronic excited state nitrogen atoms and oxygen molecules are removed from R85 to R50, the energies of Q_O and Q_N 2 in R85 and R50 are nearly the same. Therefore, the reduction from R85 to R50 has no effects on the dominant exothermic reactions. The corresponding total number densities of O( 1 D) and N 2 (A 3 Σ) are shown in Figure 28a,b, respectively. The total number densities of N 2 (A 3 Σ) are almost the same. Therefore, all the reductions from R99 to R42 have no effects on the reactions related to the excited nitrogen molecules. The total number densities of O( 1 D) in R99 is higher than the ones in all the other models, but they are nearly the same in R85, R50, and R42.   To study the effects of different kinetic models on the gas heating, the comparison of the total gas heating energies and the corresponding efficiencies in all the kinetic models are shown in Figure 29a,b, where the gas-heating energy in all the kinetic models are Q GH _R99, Q GH _R85, Q GH _R50, Q GH _R42, and Q GH _R4, and the corresponding efficiencies are η GH _R99, η GH _R85, η GH _R50, η GH _R42, and η GH _R4. The sequence of the gas-heating energy is Q GH _R99 > Q GH _R85 > Q GH _R42 > Q GH _R50. The η GH also follows the same relation. Q GH and η GH in R4 are much more different than that in the other models. The gas heating Q GH , η GH , Q R , and Q ion at the end of pulse (85 ns) and the long time after the pulse (up to 10 3 ns and 10 5 ns) are shown in Table 3.
At 85 ns, from R99 to R42, Q GH ranges from 0.63 mJ to 0.74 mJ, and η GH ranges from 41% to 47%. Q GH predicted by R4 is lower (0.42 mJ) and the corresponding η GH is 52%. It has been concluded that the fast gas heating efficiency is in the range of 30%∼50% [37,38,45,46]. Therefore, the predictions of gas heating are reasonable. Q GH and η GH from R99 to R42 at 10 3 ns and 10 5 ns are all larger than the ones at 85 ns mainly due to the increase of Q R . However, due to the different Q R predicted in the models, Q GH and η GH are very different. The difference of η GH is up to 20% at 10 5 ns. It can be concluded that the reduction of the kinetic models have very weak effects on the prediction of fast gas heating. However, because the reductions have effects on the prediction of Q R at a long time scale, the reduction has considerable effects on the prediction of the slow gas heating. R4 predicts much larger Q GH and η GH at 10 3 ns and 10 5 ns, which is due to the neglect of the main exothermic reactions there. Table 3. Comparison of Q GH , η GH , Q R , and Q ion at 85 ns, 10 3 ns and 10 5 ns in all the kinetic models.

Conclusions
The effects of the chemical reactions on the discharge characteristics and energy of a plane to plane NS-DBD actuator have been studied by adopting a three-equation drift-diffusion model. Four sets of the N 2 /O 2 mixture reduced kinetic models (31-species 99-reaction, 28-species 85-reaction, 23-species 50-reaction and 15-species 42-reaction) and a simplified N 2 /O 2 mixture averaged 4-species 4-reaction kinetic model are adopted. The electron-impact reaction-rate coefficients and the electron-transport coefficients are obtained by using LMEA. The validity of the numerical model on the prediction of the discharge energy and gas heating also has been proved.
The discharge energy in the repetitive NS-DBD reaches a periodic phase equality state after the second pulse following the discharge pulse, which is mainly balanced by the gas heating and energy stored in the heavy species. All the N 2 /O 2 mixture kinetic models over-predict almost the same discharge energy during the pulse. The simplified 4-species 4-reaction kinetic model under-predicts the discharge energy, but it is qualitatively similar to that in the other models. The discharge energy is mainly affected by the electronic excitation and the energy gain by ions from the electric field. The vibrational excitation, negative ions, associative ionization, dissociation of nitrogen and oxygen molecules have weak effects on the prediction of the discharge energy.
The gas heating is mainly contributed by the ions and exothermic reactions. The contribution of ions is larger than that of the exothermic reactions. The dissociation of oxygen molecules, the quench of electronically excited oxygen atoms and nitrogen molecules are the dominant exothermic reactions in the low electric field zone (E/N ≤ 200 Td), and the dissociation of nitrogen molecules and positive ion conversation reactions are the main exothermic reactions in the high electric field region (E/N > 400 Td). The reduction of the negative ions, the vibrational excitation, and associative ionization reactions have very weak effects on the prediction of fast gas heating, but it has considerable effects on the slow gas heating. The gas heating mainly releases during the high voltage pulse and the fast gas heating efficiencies are in the range of 41%∼47% predicted by all the reaction models, which is consistent with the results in the past. However, the simple 4-reaction model predicts unreasonable high gas heating efficiencies after the pulse due to the unreasonable exothermic reaction in the model.
Based on the contributions of the different types of chemical reactions on discharge physics, specifically on discharge energy and gas heating, the present results can be used to reduce the chemical reactions for NS-DBD when different gas heating stages are specified. They can also be used to simplify chemical reactions in the flow-control applications, when fast gas heating is dominated, compromising the fidelity and the computational cost.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. The Details of the Reaction Mechanisms of all the Kinetic Models
The species and the reactions involved in the kinetic models are shown in Table A1, where the abbreviations of the species are listed in Table A2. The details of the reactions of R99 are shown in Table A3. In the tables, the reaction rates of the heavy species are calculated by the Arrhenius equation k = AT exp(E/T). The reaction rates of certain electron impact reactions are also obtained by the Arrhenius equation k = AT e exp(E/T e ), where A and E are given in the Tables A3 and A4. T and T e are the ambient temperature and the electron temperature, respectively. The reaction rates in R99 and R85 are referred from the literature [7,37,45,[47][48][49]. The reaction rates in R50 and R42 are the same as that in the references [7,8], respectively. The reaction rates of the other electron impact reactions are k (ε) = γ ∞ 0 εσ f 0 dε, where γ = (2e/m e ) 1/2 , ε and m e represent the electron energy and the electron mass, respectively. The mean electron energy probability function (EEPF) f 0 is calculated by solving a two-term approximation of the Boltzmann equation [33], and the collision cross section σ is obtained from the LXCAT project data [50]. α and η in Table A4 are the functions of ε [30].
Re_ei. N 2 (υ = 1-8      Table A3. Details of reactions of R99. Units are consistent with number densities in cm −3 , temperature in K, one-body rates in s −1 , two-body rates in cm −3 s −1 , and three-body rates in cm 6 s −1 . The medium M can be neutral molecules N 2 and O 2 .