Two-dimensional particle-in cell/Monte Carlo simulations of a packed-bed dielectric barrier discharge in air at atmospheric pressure

The plasma behavior in a parallel-plate dielectric barrier discharge (DBD) is simulated by a two-dimensional particle-in-cell/Monte Carlo collision model, comparing for the first time an unpacked (empty) DBD with a packed bed DBD, i.e., a DBD filled with dielectric spheres in the gas gap. The calculations are performed in air, at atmospheric pressure. The discharge is powered by a pulse with a voltage amplitude of −20 kV. When comparing the packed and unpacked DBD reactors with the same dielectric barriers, it is clear that the presence of the dielectric packing leads to a transition in discharge behavior from a combination of negative streamers and unlimited surface streamers on the bottom dielectric surface to a combination of predominant positive streamers and limited surface discharges on the dielectric surfaces of the beads and plates. Furthermore, in the packed bed DBD, the electric field is locally enhanced inside the dielectric material, near the contact points between the beads and the plates, and therefore also in the plasma between the packing beads and between a bead and the dielectric wall, leading to values of 4 × 10 8 ?> V m−1, which is much higher than the electric field in the empty DBD reactor, i.e., in the order of 2 × 10 7 ?> V m−1, thus resulting in stronger and faster development of the plasma, and also in a higher electron density. The locally enhanced electric field and the electron density in the case of a packed bed DBD are also examined and discussed for three different dielectric constants, i.e., ϵ r = 22 ?> (ZrO2), ϵ r = 9 ?> (Al2O3) and ϵ r = 4 ?> (SiO2). The enhanced electric field is stronger and the electron density is higher for a larger dielectric constant, because the dielectric material is more effectively polarized. These simulations are very important, because of the increasing interest in packed bed DBDs for environmental applications.


Abstract
The plasma behavior in a parallel-plate dielectric barrier discharge (DBD) is simulated by a twodimensional particle-in-cell/Monte Carlo collision model, comparing for the first time an unpacked (empty) DBD with a packed bed DBD, i.e., a DBD filled with dielectric spheres in the gas gap. The calculations are performed in air, at atmospheric pressure. The discharge is powered by a pulse with a voltage amplitude of −20 kV. When comparing the packed and unpacked DBD reactors with the same dielectric barriers, it is clear that the presence of the dielectric packing leads to a transition in discharge behavior from a combination of negative streamers and unlimited surface streamers on the bottom dielectric surface to a combination of predominant positive streamers and limited surface discharges on the dielectric surfaces of the beads and plates. Furthermore, in the packed bed DBD, the electric field is locally enhanced inside the dielectric material, near the contact points between the beads and the plates, and therefore also in the plasma between the packing beads and between a bead and the dielectric wall, leading to values of 4 10 8 V m −1 , which is much higher than the electric field in the empty DBD reactor, i.e., in the order of 2 10 7 V m −1 , thus resulting in stronger and faster development of the plasma, and also in a higher electron density. The locally enhanced electric field and the electron density in the case of a packed bed DBD are also examined and discussed for three different dielectric constants, i.e., . The enhanced electric field is stronger and the electron density is higher for a larger dielectric constant, because the dielectric material is more effectively polarized. These simulations are very important, because of the increasing interest in packed bed DBDs for environmental applications.

Introduction
Dielectric barrier discharges (DBDs) at atmospheric pressure [1][2][3][4] usually occur in the filamentary regime in most gases [5]. These filaments are developed by electron avalanches at high voltage [6]. A typical filament in a DBD may be a few hundred micrometers in diameter and can exist for several nanoseconds [6]. Filamentary discharges are characterized by a high electron concentration in narrow filaments and they produce high concentrations of chemically active species and radicals through electron-impact processes driven by energetic electrons [7]. This type of discharge is of special interest in the field of gas pollutants and hydrocarbon conversion [8][9][10].
Several numerical studies have been published in literature, dealing with filaments in DBDs. Steinle et al studied a filamentary DBD in air at atmospheric pressure through a two-dimensional (2D) fluid model with general finite-element software [11]. A 2D fluid model was also used to analyse the filamentary mode of a DBD in nitrogen at atmospheric pressure with a discharge gap of 1 mm [12]. Kushner and coworkers applied a 2D fluid model to investigate the influence of dielectric spheres, with 10s of micrometer in diameter, on the dynamics of a positive streamer in humid atmospheric pressure air [13]. A multi-fluid hydrodynamics simulation of the ion energy and angular distributions on polymer surfaces produced by atmospheric pressure DBD streamers in air was reported by the same group [14]. A particle-in cell/Monte Carlo collision (PIC/MCC) code was developed to study the streamer formation process in air from the initial stage of a single free electron in a background electric field [15], which was the first study of streamer dynamics with a particle model in electric discharges without dielectric. Akishev et al reported the development of surface streamers for negative pulsed surface DBDs in argon at atmospheric pressure in pin-to-plane and strip-on-plane geometries using a 1D fluid model [16]. Recently, Kushner and coworkers have investigated plasma filaments in high voltage (10-20 kV) DBDs in humid air at atmospheric pressure, and displayed self-organized patterns of the plasma density during a single filament and a single discharge pulse, through a 2D fluid model with their nonPDPSIM computer code [17].
Most of the above studies focused on the filaments in an unpacked DBD. Only in [13] a kind of packed bed DBD was studied. A packed bed DBD reactor, in which dielectric beads or pellets are inserted in the gap between the electrodes, is gaining increasing attention for environmental applications, such as air pollution control and greenhouse gas conversion, because a much better energy efficiency can be reached [18][19][20][21][22][23]. Moreover, when inserting a catalytically active packing material (or when coating the dielectric beads with a catalytic material), i.e., yielding so-called plasma catalysis, the selective production of specific compounds can be targeted [24][25][26]. In a packed bed DBD, the plasma properties will be modified, for example, the electric field in the vicinity of the beads will be enhanced [27]. However, Jo et al have reported that the presence of a Pt catalyst can result in a lower electric field than in a bare Al 2 O 3 packing, due to the enhanced surface conductivity [28]. Furthermore, pollutants will be adsorbed on the dielectric surface [29], and surface discharges will be developed along the surface of the packing beads, which leads to higher electron impact ionization and dissociation rates, and thus a higher concentration of reactive species around the insulating surface compared to volumetric streamers [30]. It was demonstrated that surface discharges can occur for a higher volume fraction of the packing material [31]. Furthermore, it was reported that a typical filamentary discharge could transit to a combination of spatially limited microdischarges and surface discharges [21].
In spite of the growing interest in packed bed DBDs for environmental applications, a detailed understanding of the plasma behavior and the fundamental mechanisms of the filaments in a packed bed DBD from a physical perspective is almost non-existent. Only a few, very approximate, numerical studies have been performed for packed bed DBD reactors [27,32,33]. Chang and Takaki [32] developed a simplified time averaged 1D numerical model based on Poissonʼs equation and transport equations by assuming a spherical void between the pellets to study a ferroelectric packed bed non-thermal N 2 plasma, which shows that all the plasma parameters increase with increasing applied voltage and pellet dielectric constant. Kang et al [27] used a 2D fluid model to study a ferroelectric packed bed DBD. Only the surface discharge with high-electric field was observed in a nanosecond timescale for stacked double pellets. Finally, Russ et al [33] implemented a 2D fluid model to simulate transient microdischarges in a packed bed reactor for dry exhaust gas (80% N 2 , 20% O 2 and 500 ppm NO). However, this work was limited to a short one-directional discharge with a constant driving voltage. Recently, within our group PLASMANT, a more comprehensive 2D fluid model for a packed bed DBD reactor in helium is being developed , but it does not incorporate the description of the filaments. In fact, for an accurate description of the filamentary behavior in a packed bed DBD, a PIC/MCC model is much more appropriate.
Therefore, in the present paper, we investigate, for the first time, the characteristics of the filaments in a packed bed DBD with a 2D PIC/MCC simulation, and we compare the discharge properties between an unpacked and a packed bed DBD. The outline of the paper is as follows. In section 2, the simulation model is explained. In section 3.1, our numerical result for the formation of an individual filament obtained by the PIC/ MCC method is compared to the result obtained with a fluid model by Kushner and coworkers [17], for an unpacked DBD, to benchmark our model. Section 3.2 deals with the description of the filamentary behavior in a packed bed DBD and the comparison with an unpacked DBD. The effect of the dielectric constant of the dielectric barriers and beads on the filamentary discharge behavior is investigated in section 3.3. Finally, our concluding remarks are made in section 4.

Simulation Model and Conditions
PIC simulations [34] take into account the detailed kinetic behavior of charged particles, which is not achievable in fluid simulations. The simulation model used in this investigation is a 2D PIC/MCC method developed within the VSim simulation package [35], using a fully explicit algorithm.
The geometry used in the model is shown in figures 1(a) and (b). Note that, in order to clearly show the dielectric spheres in space dimension, only the central part with packed beds is indicated in figure 1(b), while the simulation region also extends from 0 to 10 mm in the x direction, i.e., the same as in the unpacked DBD, but with only beads in the central region. The dielectric spheres in figure 1(b) are not in strict spherical packing here, because there must be some spacing for the streamers to propagate. Indeed, in reality, the geometry is threedimensional, and there will still be some spacing between the beads even in perfect spherical packing. The computational domain is 10 mm × 1.65 mm. The gas gap is 1 mm and is bounded by two dielectric plates of 0.3 mm thickness with relative permittivity, r  . As is well known, in a strict spherical packing (two rows with two beads in the top row and three beads in the bottom row) in a discharge gap d = 1 mm, the dielectric beads have a radius r d 2 3 ( ) = + mm and the distance between adjacent sphere centers is r 2 in the x direction and r 3 in the y direction [36]. In this case, however, there is no gap between the beads, and thus no possibility for streamer propagation. In order to simulate streamer propagation in this 2D geometry, we define a gap between the beads of r 0.2 in both x and y directions. Thus, in this case the dielectric beads have a radius of r d 2. = mm in the y direction, which is not a strict spherical packing, as shown in figure 1(b).
In the Poisson solver, we included the surface charge accumulation on all dielectrics by self-consistently accounting for the deposited surface charges, as well as the dielectric coefficients. Here we consider perfect dielectrics, which means that no conductivity and no charge leaking are included. The top dielectric plate is bounded by a planar metal electrode with thickness of 0.05 mm. The bottom dielectric boundary is connected to the grounded electrode. In the simulation, a voltage of −20 kV was assumed, with a rise time of 0.1 ns and then kept constant for the duration of the calculation. This is the same as used in the fluid model of Babaeva et al [17], in order to easily compare our results with the results from literature for the unpacked DBD case. The maximum calculation time is 0.35 ns for the unpacked case (section 3.1), while it is 0.85 ns for the packed bed DBD case, when assuming a dielectric constant of 30 at the top plate, 3 at the bottom plate and 9 in the beads (i.e., section 3.2), and 0.75 ns when assuming a dielectric constant of 22, 9 and 4 for the top and bottom plate and the beads (section 3.3). These maximum calculation times are based on the different streamer formation times for the different dielectric conditions. The number of mesh points is 1000 500 in the whole simulation region  mm). Subsequently, the plasma sustains itself through volumetric/surface ionization and secondary electron emission. The secondary electron emission coefficient is assumed to be 0.1 for N 2 + and O 2 + ions on the dielectric surface [37]. The choice of this value is not so critical, since the secondary electron emission has little effect on this filamentary discharge for such a short time (up to 0.85 ns). Although photoionization might be important in the mechanism of streamer formation [17], it is not taken into account in this model, because in a packed bed DBD, the photons cannot travel freely like in an unpacked DBD, as the gas gap spacing is small and the geometry is complex, so most photons will more likely be lost at the surface of the dielectric beads instead of contributing to ionization. The MCC procedure is employed to account for elastic, excitation, ionization and attachment collisions of electrons with N 2 and O 2 gas molecules, as shown in table 1. We include only three types of ions, which are the most important ones, as they have the lowest ionization threshold and thus they have a larger density compared to other ions. Furthermore, as our simulations only run for less than 1 ns, to describe the filament development, this is so short that electron-ion and negative ion-positive ion recombination will not affect the results, because the recombination reactions have a larger relaxation time (i.e., in the order of 100 ns, considering the small cross sections/reaction rate constants). Finally, some electron impact dissociation and vibrational excitation processes are also omitted here to speed up the simulations, and because these reactions will only lead to electron energy losses in the current model, and their cross sections are relatively small, thus omitting these reactions will almost not affect the electron kinetics. All cross sections are adopted from the lxcat database, and were taken from [38][39][40]. As the number of electrons and ions rapidly increases due to the ionization avalanches, a particle merging algorithm is used when the number of super-electrons or super-ions exceeds 10 in each grid. In this case, four particles are combined into two particles with both conservation of momentum and energy. The simulation time-step is fixed at 10 13 s and the simulations will run for several tens of thousands of time-steps to reach a total simulation time up to 1 ns.  Figure 2 shows the formation and propagation of a single filament in an unpacked DBD, by means of the time evolution of the electron density. Here 30 r  = in the top dielectric plate and 3 r  = (a typical value for polymers) in the bottom dielectric plate. We assumed these values in order to validate our model, by comparing with results with literature [17]. For the conditions under study, i.e., an applied voltage of −20 kV with rise time of 0.1 ns and then kept constant for the duration of the calculation, the electrons are multiplied from a few initial seed electrons at the middle of the top dielectric surface and they form an anode-directed avalanche, due to ionization collisions. The electron density inside the filament is of the order of 10 10 21 23 m −3 , but a transition to an arc does not occur, due to the dielectric barriers, stabilizing the discharge. Locally, the electron density can reach values of 10 22 -10 23 m −3 , which is a factor 100 higher than the maximum values reported in literature [17] at 0.35 ns for the same conditions. The reason is that a fluid model [17] treats the electrons as a continuum and calculates the average electron density, so it can probably not capture the very localized peak densities as shown in figure 2. The majority of the electrons, however, have a somewhat lower density, as is clear from figure 2, similar to the fluid results in literature [17]. The negative filament originating from the electron avalanche propagates to the bottom dielectric in 0.15 ns, as is clear from figure 2. Thus, we can roughly estimate that the average propagation speed of the streamers is about 10 6 m s −1 , which is similar to the values reported in the literature [17]. However, it is nearly an order of magnitude higher than the gas-phase streamer velocity reported in [41], due to the higher electric field ( 10 8 V m −1 ) compared to the value of 10 7 V m −1 in [41] for a gas-phase streamer. Until 0.1 ns, the filament is weak and the speed of the streamer is low, due to the limited avalanches from the few seed electrons. However, after 0.1 ns, the avalanche becomes stronger because of the high applied voltage with amplitude of −20 kV. When the negative streamer reaches the bottom dielectric surface, the volume streamer changes to a surface discharge along the bottom dielectric surface, due to the lateral component of the electric field, induced from the electron charges on the surface. Thus, the electron avalanche continues in the lateral direction, over the surface of the bottom dielectric. In general, the characteristics of the streamer formation and propagation are similar as in literature [17], although photoionization is not included in our work. Figure 3 shows the electron density distribution in the packed bed DBD, for the same moments in time   figure 2(a). In the time frame of 0.35 ns, the streamer propagates over a distance of only 0.3 mm, which is much shorter than in the unpacked DBD case (see figure 2). Thus, the electrons are only located between the first and second dielectric beads. The negative streamer continues to propagate towards the bottom dielectric surface when time progresses, as is illustrated in figure 4. The maximum electron densities are 2.2 10 23 , 2.1 10 23 and 5.8 10 23 m −3 at 0.5, 0.6 and 0.85 ns, respectively. These values are a factor 2 higher than the maximum values of the unpacked DBD case at 0.25 and 0.35 ns, i.e., when the negative streamer has also crossed the discharge gap. The reason why the maximum electron density is slightly lower at 0.6 ns than at 0.5 ns is because of the formation of a positive streamer at 0.55 ns. Indeed, the streamer in the gap between beads 1 and 2 propagates back towards the top dielectric surface, while the streamers between beads 3 and 4 and between beads 4 and 5 also move upwards, starting from the bottom dielectric, as shown in figure 4(b).

Packed bed DBD and comparison with unpacked DBD
We can identify four different characteristics for the filaments in the packed bed DBD compared to the unpacked DBD. First, the negative streamer needs 0.55 ns to cross the discharge gap (1 mm), which is much longer than in the unpacked DBD case. The reason is that the dielectric spheres in the gap hinder the streamer propagation by enhancing the local electric field. Second, when the negative streamer reaches the bottom dielectric surface, the electrons propagate along the surface, but they are restricted by the dielectric beads. Third, the electrons also propagate along the surface of the dielectric beads. Finally, a cathode-directed positive streamer is developed at 0.55 ns in the gap between beads 1 and 2, as well as between beads 3 and 4 and between beads 4 and 5, due to the large negative potential created at the beads as a result of surface charging. In short, the introduction of the dielectric beads leads to a transition in discharge behavior from a typical filamentary discharge (or negative streamer) to a combination of a predominant positive streamer and limited surface discharges on the dielectric surfaces.

Effect of dielectric constant of the packed bed DBD
The above results were obtained for a top dielectric plate with relative permittivity of 30, to allow comparison with the results of [17]. However, this value does not represent typical dielectric materials used in a DBD, which are mostly ZrO 2 ( It is clear that the presence of the dielectric beads in the discharge greatly enhances the electric field inside and near the dielectrics, both between the packing beads, and between the beads and the dielectric surfaces. This is especially true inside and at the surface of the bottom dielectric plate and around the contact points between the bottom dielectric surface and the beads in the bottom row, due to the increased charge deposition on the bottom dielectric surface and between the packing beads in the bottom row, as shown in figure 5(d). Indeed, the electrons are accumulated and trapped in the gap between the packing beads in the bottom row and the bottom dielectric plate, which is clearly seen in figures 6(c)-(d). A negative streamer is developed from the initial seed electrons at the middle surface of the top dielectric plate and dominates the discharge until 0.5 ns. After the negative streamer reaches the bottom dielectric surface at 0.5 ns, the discharge will change to a combination of surface discharges on the dielectric surfaces (including both the beads and the plate surfaces) and a positive streamer, as shown in figures 6(c)-(d). The surface discharge will lead to more reactive species around the dielectric surfaces compared to a volumetric streamer. We can deduce from figure 7 that the electric field enhancement is most pronounced for the largest dielectric constant considered here ( 22 r  = ). Indeed, the maximum electric field amplitude increases with increasing dielectric constant, at the different times investigated, as shown in table 2. The reason is that, both the dielectric beads and the dielectric plates are more effectively polarized for a larger dielectric constant. Thus, a stronger polarization leads to a stronger locally enhanced electric field. Again, it is well known that packing materials with larger dielectric constants increase the discharge intensity due to the electric field enhancement, which has also been reported in some early papers [42][43][44]. As a result, the electron density is also higher for larger dielectric constants, as demonstrated in figure 8 at 0.5 ns, and in table 3 at different times. Finally, it is clear from figure 8 that the streamer propagation through the packed bed DBD proceeds faster when the beads have a higher dielectric constant. This can be explained because the avalanche is faster for a stronger electric field as a results of the larger dielectric constant.

Conclusion
We have presented a 2D PIC/MCC model, obtained with the VSim simulation code, to study the properties of both negative and positive streamers in an atmospheric pressure DBD operating in air, comparing an unpacked and a packed bed DBD reactor in the filamentary regime, with an applied voltage of −20 kV and a gap size of 1 mm. The calculated electron density reaches maximum values in the order of 10 20 -10 24 m −3 inside the filaments, and there is no transition to an arc discharge due the dielectric barriers, stabilizing the discharge.
The presence of the dielectric packing leads to a transition in discharge behavior from a combination of negative streamer and unlimited surface streamers on the bottom dielectric surface to a combination of a predominant positive streamer and limited surface discharges on the dielectric surfaces of the beads and plates.
We have also compared the plasma behavior in the packed bed DBD for typical dielectric constants used for the packing materials, i.e.,      density are illustrated for these three cases. The electric field is locally enhanced inside the dielectric material, near the contact points between the beads and the plates, as well as inside the plasma, in the gap between the packing beads, and between the bead and the dielectric walls. Furthermore, the electric field enhancement is stronger for a larger dielectric constant, due to the more effective polarization. As a result, the electron density is also higher for a larger dielectric constant, which will also lead to more reactive species in the vicinity of the dielectrics. In our future work, we plan to consider radical formation and its transport to the surface, as well as the surface conditions (such as conductivity, porosity, and the effect of adsorbed gas molecules), by coupling the PIC model with a fluid model. Indeed, the PIC method alone is too time consuming to consider these processes.