Computational Fluid Dynamics Modeling of Low Temperature Ignition Processes From a Nanosecond Pulsed Discharge at Quiescent Conditions

The deflagration waves that consume fuel in spark ignition (SI) engines originate from high temperature, nearly equilibrium, plasma kernels, a.k.a. sparks. Sparks result from dumping stored electrical energy (⁠$∼1$ to 10 mJ) across a minute gaseous gap (⁠$∼0.5$ to 2 mm) that separates two conductors. Once the electric field in this vicinity is high enough, the gas breaks down, becomes conducting, and rapidly transitions to a hot plasma [1]. The hot plasma exceeds the minimum ignition energy (MIE) of the mixture so as to spawn a deflagration wave [2]. A well designed SI system can repeatably produce discharge kernels that transition to flames across a desired range of engine loads.

The need to improve fuel economy has catalyzed research in lean combustion, which currently faces two significant challenges:

1. The operation of three way catalysts is suboptimal at non-stoichiometric conditions. The leaner the mixture, the worse its NO_x remediation [3].

2. Fuel-air mixtures with low adiabatic flame temperatures suffer combustion instabilities: either partial burning or misfire, leading to cycle-to-cycle variability (CCV) [4].

The first point has received widespread consensus: stoichiometry for after-treatment can be achieved by using exhaust gas recirculation (EGR) [5]. The problem of CCV on the other hand is harder and prompts revisiting the electrical and geometrical design of ignition systems because they are directly linked to the occurrence of sporadic flame inception, a cause and symptom of CCV [6]. Ignition technologies different from SI address this problem by increasing the number and volume of ignition sites [7].

One such technology makes use of low temperature plasma discharges (LTP), wherein the plasma is far from equilibrium [8]. Such discharges can be manipulated to occupy much larger volumes than traditional sparks [9]. They can also be repeatedly formed over short periods of time (nanoseconds), unlike traditional sparks (which typically last over microseconds to milliseconds [10]), and can thus obviate cathode fall losses [11] if properly designed.

Low-temperature plasma discharges are better understood today thanks to the availability of software that can numerically simulate them in a variety of configurations [12,13]. Recent computational studies [14,15] have highlighted the role of ionization waves (a.k.a. streamers) in ignition. Research on the occurrence of CCV has also largely benefited from numerical simulations of engine relevant reacting flows [6].

However, the coupling of the chemistry of LTP discharges to that of combustion still remains challenging for the following reasons:

1. The timescales of the underlying processes, particularly electrical breakdown and fast gas heating take place over nanoseconds [16]. This is 2 or 3 orders of magnitude smaller than the typical temporal resolution required to simulate flame ignition and propagation [17].

2. The coupling of the electrical field and electron density, particularly during the transition from a glow phase to an arc phase, is prohibitively expensive to compute due to the extremely small ionization relaxation times encountered [18].

3. The small ionization relaxation times are accompanied by sheath/fall regions much thinner than flames (less than a micrometer) and therefore demand high spatial resolution [1].

4. The electrons in such plasmas are not in local equilibrium and so their energy distribution cannot be assumed to be Maxwellian. The electron's energy distribution function (EEDF) and its many inelastic (excitation, ionization, and attachment) and elastic collision contributions need to be treated differently [1]. This is why a unified infrastructure for Maxwellian and non-Maxwellian kinetics (and transport) is still missing in popular thermochemical software libraries like Cantera [19] and Chemkin [20].

5. The kind of discharge (glow, arc, corona, etc.) that occurs in a gaseous gap is wholly determined by the underlying electrical circuit [1]. Both the voltage across a gaseous gap and the current that flows through it after breakdown are dependent on passive circuit parameters like capacitance, inductance, and resistance. For modeling nanosecond pulses, active circuit components are additionally required since the distance light travels over this duration is comparable to the distance a voltage pulse travels for the same: the connecting cables must be treated as transmission lines [21] that are terminated by known complex valued impedances [22] (for numerical calculations, see Ref. [23]). Simulating this very important electrical aspect is beyond the scope of most CFD software suites.

6. At high pressures, streamers tend to branch chaotically [24]. Discrete particle methods are better suited for their simulation but are more expensive than continuum approaches [25].

The aforementioned difficulties also carry implications in comparing simulated discharges with experimental ones:

1. At high pressures spatial and temporal resolution are very demanding, especially for experiments designed to study both the discharge and ignition processes.

2. Disentangling incident and reflected voltage traces at the electrodes is non-trivial [21]. The uncertainties in circuit characterization are thus deferred to the computational boundary conditions.

3. The EEDF is difficult (if not impossible) to measure, especially at high pressures. This increases the number of free parameters that can reduce the apparent difference between experimental and numerical observations.

Nevertheless, certain questions regarding the ignitability of fuel-air mixtures that are subject to LTP discharges can be fruitfully asked. This is because:

1. The MIE of a fuel-air mixture is a property of the mixture and to some extent the flow-field [2,26]. The source of the ignition need only meet the MIE criterion in order to ignite. The criterion concerns the temperature and velocity variations across the source as these exert zeroth order control over the local consumption speed (S_c) and stretch rate (Karlovitz number, Ka) of the flame [27].

2. The effects of an LTP discharge manifest chemically and, as a consequence, gas dynamically. An LTP discharge brings about electronic excitation whose relaxation causes fast gas heating [16,28] (an example of such heating can be seen in the rightmost snapshot of Fig. 1). This local and rapid rise in pressure and temperature gives rise to a shock wave whose wake contains an incipient flame kernel [29,30]. The growth/shrinkage of this kernel is effectively decoupled from the electro-chemical character of the plasma owing to the large disparity in timescales: three to six orders of magnitude [14].

Therefore, in order to explain the effects of an LTP discharge as far as ignition is concerned, one need only approximately capture its spatiotemporal gas heating. This can then be modeled as a source term in any suitable form of an energy conservation equation. The source term must comprise parameters that adjust:

1. the width and magnitude of heating per unit distance of streamer propagated.

2. the duration and temporal profile of pulsing (and by extension, gas heating).

3. the pulsing frequency.

The objective of this paper is to introduce such a source term S_LTP, whose parameters are constrained by the properties of LTP discharges. Such models aren't currently implemented in commercial CFD software (as far as the authors are aware), let alone based on the properties of streamers like those shown in Fig. 1. The usefulness of S_LTP is demonstrated here by invoking it in numerical simulations that seek to explain the difference in ignitability between ethylene (C₂H₄)-air and methane (CH₄)-air mixtures.

Experiments were conducted for both stoichiometric C₂H₄-air and CH₄-air mixtures at 343 K and 1.3 bar. A point to point geometry was used, as seen in Figs. 2 and 3. A voltage pulse was generated using a Transient Plasma Systems SSPG-101-HF high voltage pulse generator with $∼12$ ns full-width at half max pulse and an $∼5$ ns pulse rise time. The energy delivered by the system for each pulse was roughly 5 to 8 mJ, but the amounts delivered to the gas could be significantly lower than this owing to transmission line reflections [21]. It was observed that C₂H₄ was capable of igniting with a single pulse (Fig. 2) while for the same voltage profile, CH₄ did not: it was found that CH₄ required two or more pulses in order to ignite (Fig. 3), depending on the dwell time and peak voltage.

S_LTP appears as a nonautonomous source term in the global energy conservation equation. It carries units of energy per unit volume per unit time (⁠$J/m3/s$⁠). Its spatiotemporal character is designed to reflect the fast gas heating that occurs in the wake of streamers. Electrical breakdown of the gas occurs via the propagation of positive (cathode/negative electrode directed) and negative (anode/positive electrode directed) streamers [31]. This comprises the following possibilities for gas heating:

1. The positive streamer is much faster (and therefore more reactive) than the negative streamers so the heating of the wake is significant in a cylindrical region of discharge.

2. The positive streamer is much slower than the negative streamer. This is accompanied by a conical layer of reduced heating.

3. The positive and negative streamers are of similar speeds and collide somewhere between the electrodes. An example of such a case is shown in Fig. 1.

After breakdown, within the duration of the first applied voltage pulse, a surge of current typically follows, raising the electron density of the growingly quasi-neutral plasma [32]. The shape of this current profile versus time is strongly dependent on the underlying circuit: if the real part of the upstream impedance is large, the current is small, and can appear as distorted pulses depending on the magnitude of the imaginary part. Primary and secondary reflections can follow depending on the characteristic impedance and length of the cable connecting the generator to the electrodes. The product of the current and voltage across the electrodes gives the instantaneous power delivered to the gas. This power in a single pulse is split into [1]:

1. Electrical breakdown: enough electrons should acquire kinetic energy in excess of the ionization energy per molecule in order to ionize the gap and make it conducting.

2. Ionization in plasma: when the local electric field is high enough in the plasma, ionization rates can exceed recombination rates, causing a second exponential increase in electron density.

3. Cathode and anode fall heating: where the electric field variations are particularly concentrated.

4. Joule heating: due to elastic collisions of electrons with surrounding molecules.

5. Fast gas heating: due to inelastic collisions of electrons with surrounding molecules.

6. Radical formation: due to inelastic collisions of electrons with surrounding molecules.

Items 1 and 2 have no direct impact on ignition: they only determine the local concentration of electrons. Items 3, 4, 5, and 6 however do, to differing extents.

Item 6 is neglected in S_LTP, i.e., there is no equivalent term that produces radicals in the discharge: these will be produced solely by thermal mechanisms. Although it is known that hot and neutral radicals can decrease ignition delay times τ_ign [33], the extent of this decrease is still impractical to exploit in engine applications as their effect is significant only at low temperatures, where τ_ign is larger than engine reciprocation timescales [34]. In Ref. [34], it was shown that thermal bottlenecks arise during reaction progress: although radicals initially accelerate fuel consumption, the process is swiftly choked due to the participation of thermally sensitive intermediates. τ_ign is a stronger function of temperature and so, for this study, S_LTP is simply a gas heating term. Any fraction of S_LTP that goes to radical formation via thermal decomposition is partially returned when the radicals recombine. This follows rapidly after a real discharge and is the primary reason for large decreases in ignition delay time in the presence of excess radicals [34].

is simply the portion of the energy delivered to the gas that contributes to its increase in temperature.

Of this E_LTP per pulse, only a fraction is available for ignition of the gas as much of it is carried away by a shock wave. For this reason, multiple pulses may sometimes be required. Two situations can arise in a quiescent environment for voltage pulses of similar amplitude:

1. The field applied after breakdown is so low that recombination dominates and thus another breakdown occurs in the subsequent pulse.

2. The field applied after breakdown is so high that ionization dominates and the medium remains conducting when future voltage pulses are applied. Often times, this can lead to a transition to a spark [32].

In nonquiescent conditions, if the rotational frequency of a local eddy is much higher than the frequency of pulsing, it can be assumed that the voltage pulses see fresh nonionized gas. If the opposite holds, a situation akin to quiescent conditions occurs.

z_i is the location along the axis where the positive and negative streamers impinge upon each other. z_a is the axial location of the tip of the anode. H_a is the maximum value of h_a between z_a and z_i.

α and β control the spatial extent of axial heating.

where fl is the floor function and H is the Heavyside function, $τ̂$ is the ratio between the duration of an individual pulse (say of incoming electrical energy) and the duration between pulses τ_p, and N is the number of pulses.

The parameters used in the current study are listed in Table 1 unless stated otherwise.

There is much to be said about the above formulation:

1. Equation (2) covers a large number of possible post-streamer heating profiles. The above parameters give sufficient means to control the cone angles, collision lengths, and amplitude of heating due to streamers.

2. B(t) is a periodic rectangular pulse, so every pulse deposits the same amount of energy. In light of circuit data, this can be appropriately modified. Note that similar E_LTP in time does not imply similar temperature increases in time. For instance, if the first discharge is strong enough, a shock wave lowers the local density and causes the gas to heat by much higher temperatures for the same E_LTP in future pulses. This is in fact what occurs in streamer to arc transitions where the heating occurs due to enhanced reduced electric fields [32].

3. Equation (1) can be handled analytically, but to keep S_LTP general, numerical quadrature was used, namely a MISER Monte–Carlo integration method [35]. For the above parameters, E_LTP evaluates to 0.71 mJ/pulse.

4. The inhomogeneous heating is an important feature of LTP discharges, particularly for large gaps [14]. This influences the shock wave dynamics and ultimately the flame kernel shape and size.

5. Heating close to the electrode tips due to heat transfer from anode and cathode fall regions is incorporated in the function f. An example of such a heating profile can be seen in the rightmost snapshot of Fig. 1: the anode, cathode, and streamer impingement heating zones are clearly visible. Heating can be made arbitrarily concentrated near the electrode tips to mimic the fall regions, but must be appropriately resolved.

6. In the presence of cross-flow at the electrodes, the source can be made to follow the flow as is done in Lagrangian energy deposition models like LESI [36]. The arguments x, y, and z need only be replaced by $x−vxt, y−vyt$⁠, and $z−vzt$ respectively, where v_x, v_y, and v_z are the components of the local gas velocity.

7. Energy deposition models like ISSIM [37] that assume a relationship between current and voltage are generally inapplicable for LTP discharges since the timescales for which such relationships are derived are those of the local flow. Fast gas heating due to electronic de-excitation is independent of the local flow field.

8. Multi-temperature models [38] do not capture all of the relevant streamer physics unless coupled to charge transport. Given the large uncertainties in gas heating rates, it is cheaper to use a simpler source term like S_LTP.

S_LTP was invoked in reacting flow simulations in the converge CFD software [39] suite using the user defined function interface. No turbulence model was used as the environment is quiescent and interest here is only on the initial stages of ignition and propagation when the flow is largely laminar. Adaptive mesh refinement was activated and the minimum cell size was set to 31.25 μm: the flame thickness of the slowest flame in the current study (stoichiometric CH₄-air) is $≈477$μm while it is $≈288$μm for the fastest flame (stoichiometric $C2H4−air$⁠). Coarser grids were observed to lead to extinction, possibly due to excessive numerical dissipation. The chemical kinetic model used was foundational fuel chemistry model (FFCM) [40]. The geometry of the electrodes was the same as that used in experiments and their surfaces were meshed using the Constructive Solid Geometry tool SDF [41]. The domain surrounding the electrodes was a sphere of radius 5 mm, and was assigned an outflow boundary condition in order to avoid pressure increase during flame propagation. converge CFD automatically handles meshing in the gas phase.

where γ is the adiabatic index, L is a suitable length scale, α is the thermal diffusivity, and c is the speed of sound. With $γ=1.4$⁠, $L∼1$ mm, $α=2×10−4$$m2/s$⁠, and c = 850 m/s (calculated for a temperature, T = 2000 K), this value is 10 μs, which is around the time the shock wave is seen to depart from the kernel as indicated by the density of cells allocated by adaptive mesh refinement.

It is possible to estimate the degree of expansion for perfectly cylindrical kernels using a simple Riemann problem solver (as was done in Ref. [14]). Figure 5 shows the results of such a calculation: the final radius of hot gas increases with the initial energy and significantly deviates from isentropic expansion (constant $Tr2(γ−1)$⁠) due to entropy generation from the passage of shock waves. For successful flame propagation, the combination of final temperature and radius needs to be large enough so as to exceed the mixture's MIE.

The nonzero mass fraction of formaldehyde (⁠$CH2O$⁠) in Fig. 4 represents the reaction zone of the flame: successful flame propagation is observed. A stoichiometric CH₄-air mixture pulsed with the same S_LTP is shown in Fig. 6. No flame propagation is observed. The reasons for this can be systematically probed:

1. τ_ign of the CH₄ mixture is larger than C₂H₄ at higher temperatures [34,43]: see Fig. 7. This is notable from Fig. 6 as well: at 10.3 μs, CH₄ is consumed only near the electrode tips.

2. The flame speed of CH₄ (34 cm/s) is lower than C₂H₄ (63 cm/s) (computed using Cantera and FFCM). The flame speed of the C₂H₄ mixture can be made lower by diluting it with N₂. The exact molar proportion of O₂ to N₂ that accomplishes this is $1:5.115$⁠, according to FFCM. Such a dilution has an insignificant effect on τ_ign, as seen in Fig. 7. When the diluted C₂H₄-air mixture is pulsed with the same S_LTP, it indeed fails to propagate. Figure 8 shows that although C₂H₄ is mostly consumed along the axis, the incipient flame is too slow to overcome the heat losses due to entrainment and diffusion. This shows that τ_ign is not the sole limiting factor.

3. Entrainment is the most important fluid mechanical interaction between a discharge kernel and the surrounding gas as it determines the local cooling and growth/stretch rate of the incipient flame [42,45,46]. This can be seen from the streamlines of Fig. 9. The formation of toroidal structures is due to momentum conservation during the transition of a cylindrical shock to a spherical one [47]. The closed streamlines indicate trapped gas that is rapidly mixed along its periphery with the surroundings. Secondary toroids are observed to form when ignition does not take place: this is because small vortex structures are effectively destroyed by flames [48] and is a suggestion that ignition enhancement via turbulence is possible only at scales larger than the discharge kernel. For the CH₄ and diluted C₂H₄ cases, although entrainment increases the size of the kernel, it also cools the mixture, proving detrimental to ignition. The effects of entrainment are best seen in the later stages of Figs. 6 and 8.

4. Both stoichiometric CH₄-air and C₂H₄-air flames possess positive Markstein lengths [49], i.e. their consumption speed decreases as their stretch rate increases. This effect can be calculated using a low Mach one-dimensional reacting flow code like LTORC [50] and is shown in Fig. 10. At low values of stretch, the drop in consumption speeds is roughly the same for both CH₄ and C₂H₄. At higher stretch rates, the effect on consumption speed is highly non-linear: high enough Karlovitz numbers (a non-dimensional measure of stretch [43]) can cause local extinction. The stretch rate is directly proportional to the curvature of the flame, thus implying a minimum kernel radius r_c for successful flame propagation. It is this radius that is proportional to the MIE [2] and is often easier to compute than measure. Using LTORC, and taking the limit of $ELTP→0$⁠, the values of r_c for stoichiometric C₂H₄-air, CH₄-air, and C₂H₄-air-diluted are found to be 0.36 mm, 0.82 mm, and 0.61 mm respectively. The first of these mixtures is the one that requires the smallest MIE.

5. Since r_c is larger for CH₄, its MIE is higher. E_LTP was thus increased from 0.71 mJ to 1.41 mJ by increasing σ_f from $3σg$ to $12σg$⁠. Here too, much like in the experiments, flame propagation fails for a single pulse. So S_LTP was extended to produce a second pulse of 1.41 mJ 300 μs after the first. This too failed, as can be seen in Fig. 11. Maly and Vogel [10] once pointed out that the success of nanosecond pulsed discharges lies in their ability to form kernels away from electrodes, thereby reducing their heat loss and promoting ignition. But for the LTP discharge in question, the resulting entrainment is observed to extinguish the flame.

6. For E_LTP = 1.41 mJ/pulse, decreasing the dwell time to 100 μs however does lead to ignition (Fig. 12). This implies that multi-pulse ignition is very sensitive to pulse timing: if enough cooling has occurred between pulses, more energy is required per pulse to achieve the required MIE.

A simple model that captures the fast gas heating effects of LTP discharges was implemented and demonstrated in three-dimensional ignition simulations with two different fuels. The source term S_LTP and its integral E_LTP numerically quantify the most relevant aspects of LTP discharges pertaining to ignition. The causes for success (and failure) of ignition were traced to fundamental properties of the mixtures studied (ignition delay time and flame speed) and the energy deposited (shock wave losses and entrainment). Specifically:

1. Shock wave expansion was found to alter the conditions of the kernel. For stoichiometric C₂H₄-air mixtures, this alteration satisfied the MIE criterion while it did not for stoichiometric CH₄-air mixtures.

2. The ignition delay times are a necessary indicator of ignition success.

3. Entrainment cools the kernel while also increasing its size. These two phenomena oppose each other and control the MIE of the mixture.

4. A higher energy per pulse alone was insufficient to transition a CH₄-air kernel to a flame: a second pulse was required.

5. The dwell time between pulses plays a decisive role: longer dwell times tend to require more energy per pulse to guarantee flame propagation.

S_LTP can be improved in light of more exhaustive measurements of energy delivered to the gas and the variation of this energy from pulse to pulse. In turbulent environments, where implementing S_LTP may prove too expensive, a flame initialization based on the gas-dynamic outcome of S_LTP (like that shown in Fig. 5, for example,), will need to be used. Furthermore, the flame initialization will need to mirror the stochasticity of the turbulence, so the interaction of terms like S_LTP with the surrounding flow-field will need to be studied more. This will be the subject of a future paper.

Argonne National Laboratory's work was supported by the U.S. Department of Energy, Office of Energy Efficiency and Renewable Energy, Office of Vehicle Technology under contract DE-AC02-06CH11357. The authors wish to thank the DOE Technology Managers Michael Weismiller and Kevin Stork, and the DOE Program Manager Gurpreet Singh, for funding this research. The authors would also like to thank the LCRC community at ANL for providing the HPC capabilities to run these simulations.

• Vehicle Technologies Program (Award No. DE-AC02-06CH11357; Funder ID: 10.13039/100006177).

E_LTP =

integral of S_LTP over space and time

FFCM =

foundational fuel chemistry model

LTP =

low temperature plasma

MIE =

minimum ignition energy

SI =

spark ignition

S_LTP =

source term for modeling gas heating due to LTP

τ_ign =

ignition delay time