Streamer Discharge Modeling for Plasma-Assisted Combustion

1. Introduction

In the last two decades many positive results have been reported in which non-thermal plasma-assisted combustion has shown a reduction of ignition time delay, extension of flammability limits, and flame stabilization at various fuel-air mixture flows including supersonic and fuel composition including hydrocarbons and hydrogen [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19]. Plasma assisted combustion is an extreme non-equilibrium excitation of the gas induced by electrical discharges in contrast to common combustion which is a thermal process requiring high temperature for ignition. In the plasma, electrons gain energy from an external electric field and through collisions, transfer this energy into processes that create excited radicals. Due to the slow exchange of energy to heat during the electrical discharge phase, the temperature of the heavy particles remains low compared to electrons leading to non-equilibrium conditions. The essence of such plasma is the nonequilibrium population of excited states which leads to an increase in the reactivity and facilitates ignition and flame propagation. From this point of view, the most important question for plasma-enhanced combustion chemistry is the partitioning of the electrical discharge energy into heat and vibrational, rotational, and electronic excitation. The rate of molecule excitation by electron impact in discharge plasma depends on the electron energy. Thus, the possibility to control the electron energy means the possibility to control the direction of energy deposition and selective excitation of different states of molecules.

Improvement in fuel efficiency and reduction of emission from hydrocarbon combustion used for automotive and gas turbine platforms can be achieved by improved control of lean flame blowout which will lead to reduction in the emission of greenhouse gases [20]. The low temperature combustion was originally investigated to reduce emissions, mainly NOx and particulate matter. More recently improving efficiency and lowering fuel consumption has also been a motivating factor [20,21,22]. An electrical discharge can be the source of energy which creates a plasma and channels energy into processes that can be used to control combustion. It has been demonstrated that nonthermal plasmas have the potential to control ignition/combustion.

The existing means of creating the plasma in the combustion zone include injection from an external source, spark, dielectric barrier discharges, near surface electric discharge, microwave discharges, repetitive pulsed discharges, and DC discharges [3,4,5,6,23]. The repetitive pulsed discharges have been shown to be the most effective in experimental demonstration in controlling ignition and combustion characteristics. In such a discharge, the energy coupling can be controlled by the applied voltage, pulsed width and repetition rate [24]. Plasma assisted combustion is a promising way to stabilize ignition at high altitudes and at low dynamic pressures and temperatures [7,8,9].

There is a need for fundamental understanding of the discharge physics that leads to the onset of ignition of the air fuel mixture. There are numerous reports of modeling and simulation of plasma assisted combustion in hydrogen, methane and ammonia as fuels. The majority of the reported results use zero-dimensional models to estimate radical species production. There are several reports on modeling of the discharge phase as streamers to estimate radical production [25,26,27]. Bouwman et al. discuss the streamer characteristics such as branching and the suppression of photoionization, and reactive species that are generated due to direct electron impactof 3D particle simulations in air-methane streamer [26]. Nagaraja et al. investigated the ignition characteristics of hydrogen-air mixtures using pulsed nanosecond dielectric barrier discharge [27]. The model was one-dimensional at low pressures which limits its validity at or near atmospheric pressure where at least two dimenstional models are required to correctly estimate fields due to space charge. Barleon et al. investigated the impact of nanosecond repetitive pulses frequency on the ignition characteristics of a lean methane-air mixture using a fully coupled plasma-combustion model in cylindrical geometry. They conclude the minimal energy to ignite decreases with increasing frequency because the plasma channel does not have the time to diffuse, and the necessary radical concentrations remain confined [25].

In this paper we investigate the modeling and simulation of the transient phase of the electrical discharge followed by adiabatic combustion to quantify a few macroscopic properties which are critical for the chemical conversion phase. Proper modeling of the discharge phase would lead to better predictable outcomes as this phase determines the critical radical densities for chemical conversion reactions leading to ignition. As noted earlier, there are numerous reports of streamer simulations for combustion in various gas mixtures. However, this paper reports on some properties of the steamers which can be generalized for combustion studies. This will lead to predictive models for plasma enhancement with energy addition, such as, the effect of specific plasma-produced species, efficient energy coupling into a high-speed and highly turbulent reactive flow, and in situ production of elevated temperatures. [12]. The results are presented for hydrogen/air fuel mixture. However, the proposed approach is valid for other fuel/air systems. For the successful application of this methodology to a specific application, several fundamental data sets are required including the electron impact cross section of the constituent gas molecules/atoms and chemical reaction rates of the neutral particles with the radicals produced by the discharge.

2. Modeling of the Transient Electrical Discharge and Combustion

In streamer discharges the ionization front propagates at speeds several times higher than the local drift velocity and experience steep spatial gradients due to rapid growth of ionization [28,29]. The energetic electrons in a streamer discharge at or near atmospheric pressure produce excited species and charged particles which are very reactive. The electrons see a spatially varying electric field due to the space charge created by the ionization front. This process typically is in the few nano-second time scale. Following the electron impact process, several plasma chemical reactions including recombination and neutral excited species reactions in the time frame of 100s of nano-second produces the precursors for slower chemical conversion reactions. The final stage consists of the chemical conversion for the desired applications. The time frame for a repetitive pulsed nanosecond discharge for plasma-assisted combustion is shown Figure 1, where the final stage consists of the exothermic reactions initiated by the fuel and oxygen precursors.

The first step in the simulation process is to determine the radicals produced by the streamer discharge under the given condition of the gas composition and the applied discharge voltage. Since the transport and rate coefficients are dependent on the gas composition, these quantities are determined from the solution of the Boltzmann equation with the appropriate electron impact cross sections.

During the transient phase the heavy particles do not gain energy in the short period and the neutral gas and ion are at or near room temperature. Also, in the time scale of interest (few ns) the ions can be considered to be stationary compared to the lighter electrons. We have previously reported that at atmospheric pressure for molecular gases, the first order model gives satisfactory results. The first order fluid model for an attaching gas includes the following particle conservation equations. [30,31,32,33]

$${\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. The quantity S represents various ion/electron source or sink mechanisms such as photoionization, recombination, or remnant space charge in repetitive discharges. In the absence of magnetic field and assuming the velocity of the electrons is large compared to the slow species, and the plasma is isothermal, the particle flux can be obtained from the momentum conservation equation and is given by [34,35]

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

(4)

where µ_e is the electron mobility and D_e is electron diffusion coefficient, In the first order fluid model, the transport and rate coefficients are determined by the local reduced electric field, E/N. In slowly varying electric field where the magnetic field can be neglected the electric field E is obtained from the solution of the Poisson equation [30,31].

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

(5)

where q_e is the unsigned electron charge and ${ϵ}_o$ is the free space permittivity.

The discharge current can be calculated from the electron density and the drift velocity of electron in the gap from the following equation [30].

$$I_g={\displaystyle \frac{q_e}{d}}∰-v_{ez}{n}_e{d}^{3}r$$

(6)

where d is the gap distance and $v_{ez}$ is the electron drift velocity in the axial direction. The cumulative energy deposited in the gap which is used to determine the energy partitioning is obtained from the gap voltage, $V_g$, and the discharge current.

$$ℇ\left(t\right)={\int }_0^t{V}_g{(\tau )I}_g(\tau )d\tau$$

(7)

A generic electron impact excitation process is shown in equation 8 where the rate coefficient $k_x$ is spatially and temporally dependent on the local electric field, E/N, obtained from the streamer simulation. At each time step, equation 8 is used to calculate the generation of the excited species ”x” from the spatial electron density distribution and the rate constants which is a function of local electric field [36].

$${\displaystyle \frac{d{N}_x\left(r,z,t\right)}{dt}}=n_e(r,z,t)k_x(\mathrm{E}/\mathrm{N})M_x$$

(8)

where N_x is the excited species concentration, and $M_x$ and $k_x$ are the concentration of the neutral species and rate constant for the electron impact process respectively.

In the results presented, we bypass the avalanche phase of the streamer development by placing a neutral plasma at either electrode. The determination of spatial photoionization requires the cross-section data which is not readily available for gas mixtures. To overcome this problem, the photoionization term “S” in equations 1-2 is simulated by a very low density spatially uniform neutral plasma. These assumptions have been widely studied and its impact on the generation and propagation of streamers are not critical for plasma generation for combustion [30,31,32,33].

The set of equation 1-5 was solved numerically using finite difference method. The model is two dimensional with azimuthal symmetry. The Flux-Corrected Transport (FCT) method proposed by Boris and Book was used for the convective term of the electron density equation [37,38]. This method is particularly suitable for handling the steep density and field gradients encountered in streamer propagation. The details of the method as applied to streamers have been expensively reported [30,31,32,33]. The Poisson’s equation is solved for the electric potential by Successive Over Relaxation (SOR) method [39]. This iterative method converges rapidly as there is only small perturbations in the space charge density for the time steps used in the simulation. A detailed description of the numerical method and boundary conditions can be found in reference 39.

For combustion processes, the important precursors produced in the discharge are the dissociated products resulting from electron collision excitation to repulsive electron states. The important products include oxygen and hydrogen atoms and OH radicals. A subset of the reactions leading directly to dissociation is listed below [40,41,42,43].

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

(9)

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

(10)

$$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$$

(11)

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

(12)

$$e+H_{2}O\underset{9.13eV}{\to }H+OH$$

(13)

$$e+H_{2}O\underset{11.5eV}{\to }H+H+O$$

(14)

There are other pathways for the formation of H, O and OH through the reaction of discharge generated reactive species with neutrals which is in the second time duration of afterglow plasma chemistry in the 100 ns time range shown in Figure 1. A significant fraction of the energy during the discharge phase goes into exciting the nitrogen molecular vibrational and electronic excited states, and some of these excited species can dissociate molecules. A partial list of H and O formation by the excited nitrogen molecules produced in the discharge is shown below. This group of excited radicals and neutral consist of a total of 28 reactions [44,45,46].

$$N_2\left(B^3{\mathsf{\Pi }}_g\right)+O_2\underset{3.{·10}^{-10}c{m}^3/s}{\to }{N}_2\left(\nu \right)+2O{(}^{3}P)$$

(16)

$$N_2\left(\mathrm{a}\mathrm{’}1\mathsf{\Sigma }\mathrm{u}\right)+H_2\underset{2.6c{m}^3/s}{\to }{N}_2\left(\nu \right)+2H{(}^{1}S)$$

(17)

The dissociated products such as O, H and OH undergo reactions with O₂ and H₂ during the combustion phase_. This leads to a set of chain reactions shown below [44,45,46].

$$H+O_2\to O+OH$$

(18)

$$O+H_2\to H+OH$$

(19)

$$OH+H_2\to H+H_{2}O$$

(20)

$$O+H_{2}O\to OH+OH$$

(21)

The exhaustive list of chemical conversion reactions which consists of 22 forward and 22 reverse reactions along with the rate constants can be found in reference 43. The chemical conversion code consists of coupled first order nonlinear differential equations which are solved simultaneously. The reactions rates are dependent on the temperature and generally take the form of the Arrhenius equation. Since these reactions change the gas temperature, small time steps are required to calculate the change in temperature from the enthalpy of the reactants and products. These reactions are temperature dependent and are either endothermic or exothermic. The final combustion product is H₂O and the reactions leading to H₂O formation are exothermic resulting in an overall increase in temperature and eventual ignition of the hydrogen/air mixture. These are relatively slow reaction compared to the discharge phase which eventually form water as a byproduct.

3. Results and Discussion

The discharge configuration for this study consists of parallel plate electrodes separated by 5mm. A Gaussian shaped plasma (equal number density of electrons and positive ions) is placed on one of the electrodes. The shape of the initial distribution is given as n_e(z, r) = 10¹⁸exp(-z²/0.001)* exp(-r²/0.001) m^-3. A constant background charge of 10¹¹ m^-3 is added to simulate the photoionization. This configuration produces a single streamer discharge whose properties are characterized for this study. The simulation is stopped when the streamer bridges the gap.

Most of the results presented here is for an applied step voltage of -25 kV and 25 kV (186 Td) for cathode and anode directed streamers respectively. The gas composition is for hydrogen gas in air for an equivalence ratio of ϕ=1, which is the ratio of hydrogen to the oxygen content in the gas that is required for complete stoichiometric combustion. When the initial charge is placed on the anode, the streamer that is formed propagates towards the cathode. Figure 2 shows the cross-sectional image of the electron density, and the electric field of a cathode directed streamer for several instances in time as the streamer propagates. The density and the electric field show sharp gradients at the tip of the streamer. When the initial charge is placed on the cathode, an anode directed streamer develops. The plots of the electron density and axial electric field along the axis of an anode directed streamer is shown in Figure 3. These plots are typical of streamer propagation where the streamer tip shows high gradients for electron density and electric field. In the streamer bulk, away from the tip, the electric field are nearly constant. The extent to which the field enhancement extends beyond the tip of the streamer depends on the radii of the tip. These plots illustrate the ionization fronts which is a characteristic of streamers.

The rate constants of electron impact processes (partial list: equations 9-14) depend on the reduced electric field (E/N) and gas composition. These rates are calculated from the electron impact cross sections using the BOSIG+ Boltzmann equation solver [47,48]. A useful measure of the efficiency of radical production is the G-factor which is defined as the number of radicals produced for 100 eV of electrical power energy input to the discharge. Numerically the number of radicals produced in the discharge is calculated from equation 8 and the energy input is obtained from equation 7. This method of obtaining the G-factors gives an average estimate over the volume of the discharge. The electron impact G-factors of dissociation products as a function of time is shown Figure 4. The plot is shown starting about 0.4 ns when the streamer is fully formed from the initial charge distribution. The G-factor of radicals is almost constant during the streamer propagation. This has also been observed in other fuel air mixtures. This may sound counter intuitive since the electric field enhancement increases with time as shown in Figure 3(b). However, at the high reduced electric field at the tip the rates for most electron impact processes start to flatten. Also, due to shielding, bulk of the electrons are under a low electric field which favors low threshold energy process such as vibrational excitation. This is an important result which suggests that for a given gas composition and applied voltage, the energy partitioning is fixed, and the species density primarily depends on the energy deposited in the volume of gas. Therefore, the G-factor can be defined as a macroscopic property of the streamer which once evaluated can be used to determine the radical density from the electrical energy density of the discharge for chemical conversion simulations.

In most plasma-assisted combustion modeling in the literature, a zero-dimensional calculation is used to estimate the excited species density. This is similar to a Townsend type discharge at low pressures where spatial applied field does not change with time. At atmospheric pressure all discharges quicky turns into a streamer discharge from the Townsend like avalanche phase. For comparison, the G-factors for several combustion relevant species are shown in Figure 5 for a streamer and Townsend-type (the space charge field is not included to the applied field) discharge for different species as a function of the applied field. For comparison, the G-factors for a Townsend type discharge is also shown in the figure. A Townsend discharge has no space-charge field and is the same as zero-dimensional estimation of radical production. The space charge modifies the applied field significantly as shown in Figure 2 and the radical production is spatially dependent on the local electric field. As shown in Figure 5, the zero-dimensional results can deviate significantly and overestimate the high threshold energy radical production in the discharge phase. As the applied electric field is increased, the energy partitioning moves towards processes with higher threshold energy. As shown in the plot, the nitrogen vibrational and oxygen dissociation which have lower threshold decrease with increasing E/N, whereas the hydrogen and nitrogen dissociation which have higher energy threshold increase with increasing E/N.

The next set of simulations following the discharge phase was solving the rate equation of the plasma chemistry reactions (partial list: equations 16 and 17) in the time scale of 100 ns. This is done for zero-dimension time dependent solution till steady-state is reached. This involves the solution a set of simultaneous first order nonlinear differential equation for stiff problems. MATLAB has several functions for this, and we used the ode15s. This gives the final concentration of H, O and OH radicals for the next phase of simulation.

To understand the role of electrical energy deposition which is responsible for creating the precursors necessary for combustion, we did combustion simulations for three different energy density coupled to the discharge. The simulations were done for zero dimension at constant pressure and the specific heats at constant pressure were used to determine the temperature. The combustion simulations were done by solving the rate equations discussed earlier for the combustion phase which consists of 22 forward and 22 reverse reactions (partial list: equations 18-21). In practice, the electrical energy density coupled to the gas can be controlled with the applied voltage and the pulse duration of a nano-second repetitive pulse excitation. The duration of the pulse is limited to 10 ns or less to prevent arcing. These simulations were done by initializing the densities of H and O obtained from the G-factors for different energy densities. There have been some reports on the temperature rise due to fast heating from relaxation of excited electronic states [49]. Although a slower process, but faster than chemical conversion reaction, in streamers we find that any significant temperature rise from electrical energy input comes from V-T relaxation. As shown in Figure 4, at an applied field of 186 Td approximately 20% of the energy goes into nitrogen vibrational excitation. If all the vibrational (N₂, O₂ and H₂) excitations are added up it accounts for approximately 30% of the electrical energy input. Assuming a specific heat under constant pressure for the gas mixture to be 29.2 J mole^-1K^-1, for an energy density of 0.1 J cm^-3, and if 30% of the energy goes int gas heating, we determine the temperature rise to be about 23^oC. This rise in temperature is takin into account but is not significant in determining the combustion characteristics.

The results of the temperature rise as a function of time is shown in Figure 6. This was done for a starting temperature of 682 K which is close to the gas temperature during the compression phase of internal combustion engine [20]. The plasma significantly reduces the ignition delay as shown in Figure 6. The higher the input energy density, the faster is the ignition. The 0.2 J/cm³ and 0.3 J/cm³ energy input resulted in ignition in a few μs, but the 0.1 J/cm³ energy input did not proceed to ignition in the time range shown. The time dependance for an energy density of 0.3 J/cm³ of the reactants H₂ and O_2, and the product H₂O is shown in Figure 7. As expected, the time dependence of density follows the temperature rise shown in Figure 5.

Repetitive nano second pulses and dielectric barrier discharges have been studied for plasma-assisted combustion, and it has proven to be one of the most effective methods of coupling electrical energy into the combustion gas. In repetitive pulsed discharges the gas composition changes with time which will change the electron energy distribution function. Therefore, the electron transport and rate coefficient will change, and the radical generation will be impacted. Similarly, the micro discharges in a dielectric barrier discharge can be characterized by a streamer discharge prior to the voltage collapse due to dielectric charging [39]. To correctly capture the changes as the combustion proceeds, the G-factors were calculated for partially combusted (burnt) gas compositions. Figure 8 shows the G-factors as a percentage of the completed combustion. This is determined by the fraction of H₂ remaining compared to the original concentration. Also, the burnt hydrogen is used to determine the water concentration. As the hydrogen and oxygen are depleted during combustion, the dissociation fraction also decreases. However, as the combustion proceeds, the water concentration increases and therefor the dissociation products associated with it also increases. Since nitrogen concentration remains the same its dissociation remains constant,

We studied the effect of pulse repetition rate on the ignition characteristics. Figure 9 shows the onset of ignition for three different frequencies with each pulse delivering an energy density of 0.01 J/cm³ at an applied reduced field of 186 Td. In these discharges the gas composition changes as the combustion proceeds. For these simulations, during each pulse, the G-factors used for calculating the radical densities were estimated for the composition at that instant in time from the data shown in Figure 8.

From the temperature plots (Figure 6 and Figure 9) as a function of time, as the temperature reaches about 1000 K, the combustion proceeds rapidly towards the final temperature. To get to a temperature of about 1000 K, the energy deposited is roughly about 0.3 J/cm³ for frequencies shown in Figure 9. The slower rep rates take proportionally longer time to get to this threshold energy density. This suggests that in the zero-dimensional case the ignition takes please when the energy deposited reaches a certain threshold. If flow is included this would give a different result as the repetition rate and residence time would determine the energy density coupled to a particular volume of gas.

4. Conclusions

The transient phase of pulsed nano second repetitive discharge at atmospheric pressure is best characterized by streamers. The ionization front of streamers has field enhancement, but the bulk of the streamer is shielded and has low electric field. The gross properties are different from zero-dimensional assumption of a uniform Townsend type discharge. For predictive models the discharges phase should be modelled as streamer to get better results. The energy density coupled to the gas is the most significant macro property of the discharge which determines the radicals produced by the discharge. The ignition characteristic is sensitive to the energy density. For a given applied field, the density of radicals produced per unit energy density input is fairly constant during the discharge phase. The energy density which in turn depends on the pulse length can be used to control the combustion characteristics. Therefore, an accurate estimation of the radicals produced in the discharge phase is important for predictive models. For repetitive discharges, the ignition time shows dependence on the total energy coupled to the gas during the ignition delay period.