ECENT years, non-equilibrium plasma discharges have drawn great attention due to its potential to enhance and stabilize combustion in internal combustion engines, gas turbines, and scramjet engines [1-4]. Experimental studies have demonstrated that plasmas can shorten the ignition delays [5-7], extend extinction limits [8], improve flame stabilization [8], increase flame speed [9] and suppress soot formation [10, 11]. However, the understanding of the mechanism of plasma/combustion interaction is still lacking due to the complicated thermal, kinetic and transport coupling between plasma and combustion. In particular, it is still unclear which reaction pathways are dominant, and what reaction pathways are still missing. In order to answer this question, a combination of both computational and experimental efforts with detailed kinetics is required.

In contrast to experimental efforts, there have been only a few detailed numerical studies of plasma assisted combustion (PAC). Numerical simulation can be immensely advantageous in complementing the experimental efforts and in providing significant insights into the mechanism of plasma enhancement of combustion. The multi-scale nature of PAC, however, creates significant difficulties of stiffness for comprehensive modeling studies. As a result, many studies have resorted to simplified 0D models to gain understanding of the plasma kinetic and thermal effects in fuelair mixtures [12-15]. In such models, plasma discharge was assumed to be uniform in the entire volume during each discharge pulse. The electric field and electron density were pre-tuned such that the coupled energy can match that in experiments. These models did not consider sheath formation, and charge accumulation on the dielectric layers. In addition, species and thermal diffusion effects were ignored. In the past few years, Nagaraja et al.[6, 7, 16, 17] developed a self-consistent, 1D parallel model to simulate pulsed, nanosecond discharges in fuel-air mixtures. The framework is capable of resolving electric field transients from pulse to pulse, as well as the cumulative effects of multiple pulses on fuel oxidation and combustion. However, the comprehensiveness of this model results in significant computational difficulties even on our supercomputer capable of peta-scale and beyond. For example, the simulation of 150 discharge pulses for C2H4/O2/Ar mixtures took more than 3 weeks. One of the major computational challenges comes from the large number of species and the dramatic stiffness of the detailed plasma-combustion kinetics which need to be calculated on every grid point and every time step. In a plasma-combustion system, the characteristic timescales of different phases can range from millisecond to picosecond and even beyond. In most past works, the variable-coefficient ODE solver (VODE)[18] is usually applied to solve the stiff ODE system. However, the computation time of VODE solver increases as cubic of the number of species due to the iterative “Jacobian matrix decomposition”. For this reason, it’s almost impossible to involve large detailed mechanisms or moderate mechanisms in simulations with large number of grid points, especially in the high dimensional computational domain. In order to enable the detailed plasma-combustion kinetics in comprehensive modeling, people can either reduce the stiffness of the ODE system, such as the computational singular perturbation (CSP) method [19, 20], the intrinsic low-dimensional manifold (ILDM) method [21], the dynamic stiffness removal method [22] and the hybrid multitimescale (HMTS) method [23], or reduce the size of the kinetic mechanism, such as the visualization …