Numerical Simulation of Pin-Plane Discharge for Plasma-Assisted Combustion of Hydrogen

1. Introduction

The electrical breakdown of gases caused by electron impact ionization results in plasma formation. Depending on various electrical discharge conditions, the resulting plasmas can be thermal or non-thermal. In thermal plasmas the gas and electron temperature are equal and high. Whereas in most non-thermal plasmas the electron temperatures are high and the gas temperature is close to the room temperature. Non-thermal plasma is an enabling technology which forms the basis for many mature and emerging technologies including semiconductor processing, manufacturing and plasma medicine. [1,2,3,4]. The radicals and chemical species in non-thermal plasmas provide a highly reactive environment under ambient conditions.

At low pressures the electric discharge will produce a glow like bulk discharge, However, when the pressure is near atmospheric or higher, electrical breakdown leads to arc formation as they reach steady state [5]. If the electrical discharge can be terminated by some means during the transient phase, low temperature non-thermal plasmas can be generated. In a dielectric-barrier discharge (DBD), due to charging of the dielectric during the transient phase, the electric field decreases in the gap leading to the extinction of the micro-discharges [2]. The pulse width of nano second repetitive pulsed discharge can be lowered so that the discharge does not transition to the arc mode. Atmospheric pressure discharges start off as an avalanche, and as the space charge build-up starts to dominate the electric field, it transitions to a streamer discharge [6,7,8,9]. Due to the enhancement of the electric field in certain regions of the discharge, the transient phase of the electrical breakdown at atmospheric pressure is very effective in generating radicles which forms the basis of plasma chemical conversion in most technical applications [10]. The methods to create non-thermal plasma near atmospheric pressure have evolved over the last several decades with emergence of new applications. The choice of power source to create non-thermal plasmas is primarily dictated by the gas composition and gas pressure with radio frequency and dc discharges used for low pressure and near atmospheric pressure the sources include high voltage repetitive pulses and AC sources for DBD. [1,11,12,13]].

With a DC or pulsed voltage, the commonly used electrode geometries include plane-plane, pin-plane, and pin-pin. The type of electrode geometry used in the discharge determines the spatial profile of the electric field due to the applied voltage and space charge, and that affects the temporal evolution of the electrical discharge. Both experimental and modeling studies of the formation of streamers in plane-plane geometry have been reported extensively [14,15,16,17]. In our previous work we reported on fluid models appropriate for electrical discharges near atmospheric pressure, and simulation results of streamer propagation and radical generation characteristics for plane-plane, pin-plane and pin-pin electrodes [17,20]. At or near atmospheric pressure the pin-plane electrode geometry is frequently used for non-thermal plasma generation in particular where a liquid surface acts as the plane electrode [18,19]. Several modeling efforts can be found in the literature for the pin-plane geometry [21,22,23,24]. The effect of electrode geometry on the characteristics of streamers in pin-plane electrode geometry has been reported by Bourdon et al.[22]. They found that the radius of curvature of the tip of the pin electrode or its shape had little impact on the discharge properties of the ionization front propagation. They also found that the propagation speed of the ionization front was high close to the tip of the electrode pin and decreased as the front propagated further away from the tip. In a pin-plane electrode geometry, the effect of voltage polarity on the streamer characteristics in atmospheric pressure argon was reported by Sato et al. [23]. They concluded that stable high electron density plasma can be obtained at low-voltage operations with positive polarity applied to the pin electrode. The scaling law for pin-plane electrode setup was studied by Fu et al. [24]. They conclude that electrode geometry which produces local field enhancements is critical in determining the breakdown characteristics. Using finite element method, Yamazawa and Yamashita have done detailed calculations of the Laplacian electric field for pin-plane electrodes. [22].

Electrode geometry determines the spatial contours of the electric field which impacts the electrical breakdown of gases and partitioning of energy into processes leading excited species generation. In this paper we have extended our previous study on electrode geometries and report the results of effect of diameter of the pin electrode on streamer formation and propagation [17]. As with the previous study, this was done for hydrogen air mixture at atmospheric pressure which is of interest for plasma assisted control of combustion [23]. The combustion kinetics of hydrogen/air mixture have been studied extensively and are well documented [25,26,27,28,29,30,31,32,33]. In order to control the combustion characteristics, efficient production of relevant radicals is important. For this study a first order drift-diffusion fluid model for electrons was used along with Poisson’s equation to find the electric field for self-consistency [25,34]. The purpose of this study is to understand the effect of pin diameter in a pin-plane electrode geometry on radical generation.

2. Fluid Model for Streamer Discharge Simulation

Since the electron mass is small compared to ions, they usually have orders of magnitude higher drift velocities and diffusion coefficients [35]. In the ns timescale the heavy ions do not gain significant energy from the electric field and can be considered stationary. So, the convective part in the ion conservation equation can be neglected. Due to the significant mass difference between the electrons and neutrals, the energy transfer due to elastic collision is inefficient. Although a significant amount of energy is stored in the vibrational states, the translational relaxations (V-T) processes are slow, and in the ns time frame it does not lead to gas heating [35]. In streamer discharges there is rapid growth of charges due to ionization caused by space charge field enhancement which leads to sharp density gradients both in the electric field and electron density [36,37]. The fluid models suitable for simulation of the ionization fronts are the conservation equations for electron and ion particle density, electron momentum and electron energy. These are derived from the first three moments of the Boltzmann’s kinetic transport equations [38,39]. Since each moment equation depends on the higher-order moment tensor, some assumptions such as isothermal plasma and ideal gas law are made to close the system of equations. In the first order model all transport and rate coefficients are parametrized using the local reduced electric field (E/N). In this model the electron energy is assumed to be in equilibrium with local electric field. In the second order fluid model, the electron energy conservation equation is explicitly solved to obtain electron energy, and the transport parameters and electron impact rate coefficients are functions of the local mean electron energy [39]. The choice of order depends on several factors such as gas pressure and composition. We found that the first order model to be suitable for the study of electrical breakdown in atmospheric pressure of gases [20,34]. In attaching gasses negative ions are formed, therefore the first order drift-diffusion fluid model includes three particle conservation equations: electrons, positive ions and negative ions [6,9,20,36].

$${\displaystyle \frac{\partial n_e}{\partial t}}=-\nabla .{\mathsf{\Gamma }}_{\mathit{e}}\left(\mathrm{E}/\mathrm{N}\right)+n_e{\nu }_I\left(\mathrm{E}/\mathrm{N}\right)-n_e{\nu }_a\left(\mathrm{E}/\mathrm{N}\right)+S$$

(1)

$${\displaystyle \frac{\partial n_i}{\partial t}}=n_e{\nu }_I\left(\mathrm{E}/\mathrm{N}\right)+S$$

(2)

$${\displaystyle \frac{\partial n_n}{\partial t}}=n_e{\nu }_a\left(\mathrm{E}/\mathrm{N}\right).$$

(3)

where n_e, n_i and n_n are the electron, positive ion, and negative ion densities respectively, ν_I is the ionization frequency, and ν_a is the attachment frequency. Since the gas composition of interest in this study has oxygen, it is necessary to include equation 3. The quantity “S” in equations 1 and 2 represents various ion/electron source or sink mechanisms such as electron ion recombination, formation of electron ion pair due to photoionization, or remnant space charge in repetitive discharges. Assuming isothermal plasma and ideal gas law the momentum conservation equation can be simplified to obtain electron flux which is given by [9,36,37]

$${\mathsf{\Gamma }}_{\mathit{e}}\left(\mathrm{E}/\mathrm{N}\right)=-n_e{\mu }_e(\mathrm{E}/\mathrm{N})\mathit{E}-D_e(\mathrm{E}/\mathrm{N})\nabla n_e,$$

(4)

µ_e is the electron mobility, D_e is electron diffusion coefficient, and E is the electric field vector. For a self-consistent solution of the streamer propagation, the electric field must be calculated at any instant in time from the charge particle densities at that time. The spatial distribution of the electric field E can be obtained from the solution of the Poisson equation [3,37].

$$\left\{\begin{array}{c}{\nabla }^2\Phi =-e(n_i-n_e-n_n)/{ϵ}_o\\ \mathit{E}=-\nabla \Phi \end{array}\right.,$$

(5)

where $\mathsf{\Phi }$ is the electric potential, e is the unsigned electron charge and ${ϵ}_o$ is the free space permittivity. The energy (in eV) gained by the electrons from the local electric field during time interval $∆t$ can be obtained from the equation shown below.

$$∆ℇ\left(t\right)=∆t\iiint |\mathit{E}{|}^2{n}_e{\mu }_e{d}^{3}r.$$

(6)

Once the spatial electron density distribution is found, at each time step the rate equation shown in equation 7 can be used to find the generation of radicals of interest [17].

$$d{N}_k\left(t\right)=∆t\iiint n_e\left(r,z,t\right){\nu }_k\left(\mathrm{E}/\mathrm{N}\right)d^{3}r$$

(7)

where the excitation frequency, ${\nu }_k\left(\mathrm{E}/\mathrm{N}\right)$, is a function of the local reduced electric field, E/N.

A two-dimensional finite difference scheme assuming azimuthal symmetry was used to self-consistently solve the set of equations 1-5. Numerical diffusion is a concern due to very steep gradients in electron and ion densities, and electric field at the streamer tip. In our work, the method for handling the convective part of the electron conservation is based on Flux-Corrected Transport (FCT) [40,41]. This was proposed by Boris and Book for modeling shock fronts and uses a high order and low order scheme to correct the particle flux. In our implementation the high order flux was obtained from an eight-order interpolation, and the low order flux was derived from the donor cell method. The FCT method has been used by several authors and has been found to be very effective in dealing with ionization fronts [6,7,8]. The electric potential is obtained from the solution of Poisson’s equation at every timestep.

For the solution of the Poisson’s equation, we used the successive over relaxation (SOR) method [3]. Once the streamer forms the charge density perturbation is small compared to the previous step and the SOR method converges rapidly. The electrode configuration investigated in this study consists of pin-plane geometry of varying pin diameter and electrode separation. For setting the boundary condition in the radial direction for an applied voltage V, and an electrode separation of d. we chose a large radial distance R and the potential at that radius $\mathsf{\Phi }\left(\mathrm{R},\mathrm{z}\right)={\displaystyle \frac{\mathrm{V}}{\mathrm{d}}}\mathrm{*}\mathrm{z}$ and the on axis radial component of the electric field ${\left.{\displaystyle \frac{\partial \mathsf{\Phi }}{\partial r}}\right|}_{r=0}=0$. The pin electrode potential for a surface contour s₁ is set to $\mathsf{\Phi }\left(\mathrm{r},{\mathrm{s}}_1\right)=\mathrm{V}$, and at the plane electrode the potential is zero. The pin electrode consists of a cylinder with a semi hemispherical top with the same diameter as the cylinder. The pin diameters investigated were from 0.5 mm to 1.5 mm in steps of 0.25 mm. Beyond 1.5 mm we did not see significant difference in the results. In all cases the pin protruded from a plane electrode. The protrusion of the pin electrode from the plane surface is 0.75 mm for all pin diameters investigated. Similar electrode configurations have been reported by other investigators [22].

3. Results and Discussions

To study the effect of the pin geometry we varied the pin diameter and the gap distance between the pin tip and the plane electrode. The two gap distances (plane to plane) investigated were 5 mm and 10 mm with an applied voltage of -20 kV and – 40 kV respectively to set the applied field in the gap equal. The mesh used to discretize the space in both the z and r directions is non uniform with a finer mesh near the tip of the pin-electrode in the z direction and fine radial mesh near the axis. The smallest grid used is 5 μm and the outer radial grids are relatively large. As discussed earlier, the Pin-Plane electrode geometry with similar parameters have been reported by other investigators [21,22,23]. Bourdon et al. found that the radius of curvature of the pin electrode tip (range of 100-1000 μm) and the shape of the pin electrode, whether rod or hyperbola, had negligible impact on ionization front characteristics. We have done our simulation over a limited set of parameters primarily to get a qualitative understanding of the pin-electrode diameter effect on streamer formation and propagation and excited specie generation efficiency.

According to the Rather and Meek criteria, the streamer develops from the avalanche phase of the discharge when the electric field due to the space charge is comparable to the Laplacian electric field [42]. The simulation of the avalanche phase takes a significant computational time, and we bypass the avalanche growth by placing an electrically neutral plasma on the pin electrode. It has been shown that once the streamer forms after the adjustment of the initial charge, the streamer characteristics are not affected by the shape or magnitude of the neutral plasma [7,36]. A Gaussian-shaped neutral plasma is placed on the electrode. The shape of the initial neutral plasma (electrons and positive ions) is given as n_e(z, r) = 10¹⁸exp(−((z-z_o)²/0.001)exp(−r²/0.001)/m³ where the tip of the pin electrode is at z= z₀. We stop the simulations when the streamer bridges the gap.

Electron ion recombination in the ns time frame is not important. Therefore, the source term “S” in equations 1 and 2 is due to photoionization or remnant charge from the previous discharge. In an air discharge the photons emitted by excited nitrogen molecules in the wavelength of 98-102.5 nm have sufficient energy to ionize oxygen molecules to generate the electron-ion pairs. [21]. Photoionization is very sensitive to impurities and depends on the gas composition, and calculating the rate would require the cross-section data which are not readily available. The commonly used cross-section data is a report by Penny and Hummert where they experimentally determined photoionization rates in air, oxygen and nitrogen [43]. Most reports using this photoionization data are concentrated on understanding details such as streamer branching [44,45,46]. Wei et al. have concluded that the photoionization had negligible influence on the propagation of streamer in humid air since the rate is about three orders of magnitude lower than electron impact ionization [45]. In this study the photoionization term “S” in equations (1) and (2) is simulated by including a fixed background charge of 10¹¹/m³ distributed uniformly in space. Most reported simulations of streamers avoid the details of photoionization by adding a background charge and it has been shown that this assumption has very little impact on streamer characteristics and radical generation over a wide range of conditions and gas composition [13,14,16].

Our main interest is to understand radical production for an application. In plasma assisted combustion of hydrogen/air mixture, the radicals of interest are the electron impact dissociated products. The results presented here are for hydrogen-air mixture for complete stoichiometric combustion which has an equivalence ratio of ϕ=1 (N₂ 58%, O₂ 14%, H₂ 28%). The impact of varying gas composition on radical production for plane-plane geometry has been reported earlier [18]. For this given gas composition, the transport parameters and the rate coefficients as a function of reduced electric field (E/N) were obtained from an open-source Boltzmann equation solver, BOLSIG [47]. The electron impact cross-sections needed for the Boltzmann equation solver are available at lxcat.net which is the most comprehensive and up-to-date database [48].

Streamer development and propagation are largely influenced by the Laplacian field for the given electrode geometry and applied voltage. In figure 1 the axial Laplacian field for two different gap distances are shown. The influence of the pin diameter on the field is significant for pin diameters in the range of 0.25 to 1.5 mm. For diameters greater than this the order of field enhancement at the tip is very similar. This is true for both the gap distances studied. For distances of about 1mm away from the tip the electric field is about the same for all pin diameters. The radial extent of the axial field is shown in figure 2 at two different distances from the cathode. Close to the tip of the pin (figure 1 (a)), the low pin diameters have a narrower region of field enhancement along the radial direction. However farther away from the tip (figure 2 (b)) the radial profile of the axial field for all pin diameters is very similar. Our goal is to get a qualitative understanding of the pin geometry effect on streamer properties.

Figure 3 shows the images of the magnitude of the electric field and electron density when the streamer head was approximately at 7mm from the cathode. The time instances are slightly different since the simulation picks a time interval according to the Courant-Friedrichs–Lewy (CFL) condition [50]. The initial charge was placed on the pin electrodes which is the cathode. This results in the formation of an anode directed streamer with electron drift in the same direction as the propagation of the ionization front. The differences in density and electric field for different pin diameters are pronounced near the tip, but further away from the tip the profiles look very similar. This is to be expected since the Laplacian field shows more enhancement near the tip for smaller pin diameter, and at larger distances from the tip the field is the same as shown in figure 1.

Figure 4 shows the streak image of the streamer electron density and axial electric field for two different pin diameters. The y-axis is the distance from the cathode, and the change in the quantity of interest is illustrated by the color bar. As the streamer head propagates along the gap, the sharp changes in the electron density and the electric field at the streamer front are observed. The parabolic curvature of the front indicates the streamer is gaining speed as it propagates towards the opposite electrode. The axial electric field at the streamer front shows a sharp peak which can also be seen in figure 3 (b). The electric field in the bulk of the plasma is shielded, and near the tip the density increases with time, and the extent of the shielding also increases. There is very little noticeable difference between the two diameters shown in the figure. This is to be expected as the Laplacian field past a fraction of mm from the tip is about the same for all tip diameters we tested.

Shown in figure 5 is the cross-sectional view of the axial electric field and the electron density along the radial direction at a fixed distance of 2 mm from the cathode for pin diameters of 0.5 mm and 1.5 mm. Once the streamer head reaches the axial position of 2 mm, the radial profile for both the diameters shown have very similar characteristics both in terms of the magnitude and radial expansion. The main difference between the two pin diameters shown is that the streamer head with a lower pin diameter arrives at this axial position before the streamer head with a higher diameter. Since the Laplacian field enchantment is higher at the tip of the pin electrode with smaller diameter, the streamer formation is faster than pin electrodes with higher diameter. At an axial distance of 2 mm from the cathode the Laplacian fields are very similar for all pin diameters studied. Therefore, the temporal growth of radial profiles of electron density and electric field are similar for all pin diameters.

At every timestep, the energy gained from the electric field is calculated using equation 6. Cumulative energy over time is obtained by summing the energy gained for each timestep. The plots in figure 6 show the total energy deposited in the discharge as a function of time for different pin diameters. The energy deposition trend is very similar for different pin diameters. As alluded to earlier, the field enhancement at the tip of the pin electrode is inversely related to the pin diameter (figure 1). The higher the field enchantment, the quicker the streamer formation is. Therefore, the energy deposited at any given time is highest for the smallest pin diameter as shown in Fig 6.

In most applications of non-thermal plasmas there are certain radicals of interest which promote the relevant plasma chemical reactions. For plasma-assisted combustion of hydrogen the production of oxygen and hydrogen atoms are important in controlling the combustion chemistry. These products are created by dissociating the molecules by electron impact excitation to a repulsive electronic state which leads to highly reactive dissociated products. The important reactions leading to the generation of the atoms is shown in equations 8-11 [50,51].

$$e+O_2\underset{6.0eV}{\to }{O}_2{(A}^3{\mathsf{\Sigma }}_u^{+})+e\to O{(}^{3}P)+O{(}^{3}P)+e$$

(8)

$$e+O_2\underset{8.4eV}{\to }{O}_2{(B}^3{\mathsf{\Sigma }}_u^{-})+e\to O{(}^{3}P)+O{(}^{1}D)+e$$

(9)

$$e+H_2\underset{8.9eV}{\to }{H}_2\left(b^3{\mathsf{\Sigma }}_u^{+}\right)+e\to H{(}^{1}S)+H{(}^{1}S)+e$$

(10)

$$e+H_2\underset{11.8eV}{\to }{H}_2{(B’}^1{\mathsf{\Sigma }}_u^{+})+e\to H{(}^{1}S)+H{(}^{1}S)+e$$

(11)

The gain of electron kinetic energy is from the electric field, and energy losses are through elastic and inelastic collision. Due to the large mass difference, the electron energy loss due to elastic collisions is negligible compared to the inelastic losses. The rate coefficients of inelastic processes are dependent on the mean electron energy, and for the first order fluid model it is a function of E/N. For comparison purposes, the G-factor is often used in plasma chemical applications which is defined as the number of radicals produced by electron impact for every 100 eV of electrical energy from the power source [25]. In the simulations, equation 7 is used to find the number density of excited species from the spatial distribution of electron density and electric field of the streamer. Figure 7 shows the G-factors for oxygen and hydrogen atom production as a function of time for pin diameter of 1 mm. Also shown is the G-factor of vibrational excitation of the nitrogen molecule which is high since the threshold for this process is low compared to the electronic excitation. A significant amount of the energy goes into vibrational excitation which will eventually lead to gas heating. The G-factors depend on the electron impact rates which are higher for low threshold processes. Since hydrogen dissociation requires higher energy compared to oxygen, the G-factor for hydrogen is lower. The electric field in the bulk of the streamer is low compared to the front due to shielding. As the streamer propagates the electrons in the bulk have low mean energy due to lower electric field. This results in energy partitioning which favors low threshold processes such as rotational and vibrational modes of molecules. As a result, we observe a slight decrease in the G-factors for dissociated products as time progresses. The streamer is fully formed from the initial charge at about 0.2 ns. The G-factor of H and O atoms for different pin diameters are shown in figure 8 (a) and 8 (b) respectively for a gap distance of 10 mm. The pin diameter for this gap distance has very little effect on the G-factors and the efficiency of radical production. However, when the gap distance is 5 mm, there are some differences in the G-factor, and as the pin diameter increases, the G-factor improves. The Laplacian field for 5 mm gap shown in figures 1 and 2 show that the field enhancement is high near the pin tip for the shorter distance. However, as shown in figure 2 (a), for smaller diameter pin the extent of the radial field drops faster compared to the higher diameter pins. The larger extent of the field for higher pin diameters helps with the production of electron impact species leading to higher G-factors.

Figure 9. The time evolution of G-factors for hydrogen (a) and oxygen (b) in an air/hydrogen gas mixture at atmospheric pressure for an applied voltage of -40 kV. The plots are for varying pin diameter from 0.25 mm to 1.5 mm and 5 mm electrode gap distance.

Figure 9. The time evolution of G-factors for hydrogen (a) and oxygen (b) in an air/hydrogen gas mixture at atmospheric pressure for an applied voltage of -40 kV. The plots are for varying pin diameter from 0.25 mm to 1.5 mm and 5 mm electrode gap distance.

4. Conclusions

In a pulsed discharge, the early phase of the breakdown can be modelled as a streamer which is a space charge dominated ionization front. The dependence of the streamer characteristic on the pin diameter in a pin-plane electrode geometry was studied. We presented the electron density and electric field of a propagating streamer for different pin diameters. The Laplacian field is enhanced near the tip of the electrode, but the radial extent of this enhancement is lower for smaller gap distances. This impacts the efficiency of radical production. Since the Laplacian field does not differ significantly due to pin size for the larger gap (10 mm) studied, the G-factors showed little dependence on the pin diameter. However, for the smaller gap (5 mm) distances there were some differences with higher diameters pins with higher G-factors. The pin diameter does influence the time it takes to form the streamer which is quicker for smaller pin diameters due to larger field enhancement at the tip. But once formed the streamer properties and G-factors are very similar for all pin-diameters studied particularly for larger gap distances.