A computational study of a laminar methane–air flame assisted by nanosecond repetitively pulsed discharges

Nanosecond repetitively pulsed (NRP) discharges have been considered a promising technique for enhancing combustion efficiency and control. For successful implementation, it is necessary to understand the complex plasma–combustion interactions involving chemical, thermal, and hydrodynamic pathways. This paper aims to investigate the mechanisms enhancing a laminar methane–air flame assisted by NRP discharges by high fidelity simulations of the jet-wall burner employed in a previous experimental study. A phenomenological plasma model is used to represent the plasma energy deposition in two channels: (1) the ultrafast heating and dissociation of O2 resulting from the relaxation of electronically excited N2 , and (2) slow gas heating stemming from the relaxation of N2 vibrational states. The flame displacement, key radical distribution and flame response under plasma actuation are compared with experimental results in good agreement. The computational model allows a systematic investigation of the dominant physical mechanism by isolating different pathways. It is found that the kinetic effect from atomic O production dominates the flame dynamics, while the thermal effect plays a minor role. Hydrodynamic perturbations arising from weak shock wave propagation appear to be sensitive to burner geometry and is found to be less significant in the case under study.


Introduction
Plasma-assisted combustion (PAC) is a promising technology for enhancing combustion performance across a wide range of * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
The combustion enhancement by NRP discharges primarily occurs through three pathways: thermal, kinetic, and transport.The thermal pathway involves an increase in gas temperature, which accelerates chemical reactions and fuel oxidation [13].The kinetic pathway involves the production of highenergy electrons and ions, which in turn generate active radicals and intermediate species, enabling new reaction routes and accelerating low-temperature fuel oxidation [14].The transport, or hydrodynamic, pathway involves alterations in local flow velocity due to phenomena like ionic wind [15] and shock wave expansion [16].Understanding these pathways is essential for tailoring plasma conditions to specific combustion scenarios and for developing more efficient and cleaner combustion technologies.However, these enhancement pathways often interact, making it challenging to isolate their individual effects to identify the dominant process.For example, Evans et al [17] concluded that an increased chemical reactivity resulting from NRP discharges actuation leads to the upstream displacement of a lean-premixed methane-air stagnation flame.However, this experiment was repeated in [18] where the hydrodynamic effect from plasma was found to be the main reason for the flame displacement.To this end, numerical studies of PAC are well suited in undertaking a parametric study to assess the underlying mechanisms of flame-plasma interaction.Due to the inherent complexities, however, simulations with detailed plasmacombustion coupling are usually limited to 0D [19,20] and 1D [21].Some CFD simulations with simplified plasma models have been carried out to understand the influence of NRP discharges on flame structure.Bak et al [8] applied a reduced combustion-plasma reaction mechanism of methaneair flames and a predefined Gaussian-shaped (temporally) reduced electric field at the discharge region.Tholin et al conducted similar 2D simulations on nanosecond spark discharge [22,23].Recently, Barleon et al [24] developed a fully coupled 2D PAC solver to study the ignition of an atmospheric pressure laminar premixed methane-air mixture by multiple NRP discharges in a pin-pin configuration.Full multi-dimensional simulations with detailed plasma kinetics, however, is still too expensive to predict the flame behavior, especially with NRP discharges involving extremely fast time scales.
Alternatively, a phenomenological model proposed by Castela et al [25,26] circumvents the need for intricate plasma kinetics by incorporating the impact of non-equilibrium plasma on the temperature of gases and the concentration of species within the set of equations governing multi-component reactive flow balance in combustion phenomena, thus lowering the computational burden of direct numerical simulations (DNSs) for high Reynolds number flow and large eddy simulations (LESs) [27].The phenomenological model allows a systematic analysis of individual enhancement pathways by isolating different plasma energy pathways, without facing the physical limitations of experimental design and diagnostics.More recently, Barleon et al [28] proposed an improved phenomenological model for methane-air mixture to describe all important species from plasma reactions in a more comprehensive manner.The model was applied in their LES work of swirled flame stabilization using NRP discharges [29].
The objective of this work is to develop a PAC solver using the phenomenological plasma model in [30] to understand the effects of NRP glow discharges on a laminar methane-air flame, which has been investigated in experiments by Del Cont-Bernard et al [31,32].Unlike other PAC burner setups [17,18] where NRP discharges were produced upstream of the flame, the axisymmetric burner used in [31,32] has a discharge channel generated at the central axis and therefore enables a direct interaction between flame and discharge.Consequently, the chemical effect of plasma on flame front is preserved.It was found in [32] that OH and CH fluorescence intensities in the flame increased by up to 40% and 10%, respectively, during plasma actuation as compared to the base flame.
The paper is organized as follows.Section 2 combines the discussion on the phenomenological plasma model, governing equations, and the axisymmetric computational domain based on the experimental PAC burner, including meshing, boundary conditions, and timestep settings.It also details additional test cases designed to channel all deposited electric energy into a single pathway to assess individual enhancement effects.The results, including analyses of flame displacement, chemical enhancement, hydrodynamic effect, and flame power response, are presented and discussed in section 3. Concluding remarks are provided in section 4.

Phenomenological plasma model
Based on experimental and numerical observations, Castela et al [25] developed a semi-empirical model that captures the essential pathways and their associated time scales through which the electric energy of the discharge is transferred into the reacting gas medium.The discharge pulses predominantly cause electronic and vibrational excitation of nitrogen molecules N 2 .These excited electronic states, N 2 (A,B,C,…), are relaxed through dissociative quenching with O 2 molecules, resulting in a swift surge in gas temperature and the dissociation of O 2 into atomic oxygen.A much slower relaxation of N 2 vibrational states results in a slow gas heating process.The total rate of plasma energy deposit, Ėp , is therefore split into three channels: where Ėp chem and Ėp heat represent the amount of energy transferred into the ultrafast dissociation of O 2 into atomic O and ultrafast gas heating, respectively.Ėp vib is the rate of energy transferred to vibrational states of N 2 .These three channels are modeled as: where Y f O2 is the mass fraction of O 2 in the fresh mixture, e k represents the specific internal energy of species k. α is the fraction of pulse energy contributing to the electronic excitation of N 2 , subsequently leading to the fast gas-heating (α − η) and O 2 dissociation (η).The remaining (1 − α) goes into the slow gas heating of N 2 vibrational states relaxation.Two key assumptions are used in determining constants α and η.First, only air-plasma kinetics is considered in lean premixed flames as air is the abundant component in the gas.Second, it has been reported that for a reduced electric field in the range of 100-400 Td at atmospheric pressure, around 20% of the discharge energy goes into ultrafast heating and 45% goes into slow gas heating [13], even in lean burned gas regime [33].Therefore, the chosen plasma power fractions are approximated as fixed values in this study.Note that in the improved phenomenological plasma model by Barleon et al [28], plasma enhanced processes such as CH 4 → CH 3 + H was also considered, while this study employed the original Castela model.

Governing equations
Following Castela et al [25], in addition to the mass and momentum conservation equations, governing equations for total energy e, vibrational energy e vib and species mass fraction Y k are written as: where Ṙp VT is the relaxation rate of the vibrational energy into gas heating, ωc k and ωp k refer to the reaction rates due to combustion and plasma, respectively.The details of the constitutive relations for these terms can be found in [25].

Numerical configuration
Computational simulations are set up for the experimentally studied jet-wall axisymmetric burner from [31,32], shown in figure 1.A wall-stabilized laminar methane-air flame in V shape is actuated by pin-pin generated NRP glow discharges across the flame front, enabling radicals and excited species from plasma to directly participate in the reaction zone.The applied voltage of 8 kV across the 10 mm gap in the experiment of [32] is estimated to produce a reduced electric field within the validity range (100-400 Td) of the phenomenological plasma energy pathway described above.The phenomenological plasma model was originally implemented in DNS reacting flow solver YWC by Castela at al [25].It was later incorporated into several other reacting flow solvers including OpenFOAM [34] by [35], YALES2 by [27,36] and AVBP by [29,37] for LES.This work also chose As shown in figure 2, the 2D axisymmetric computational domain has three layers of grid size varying from 20 µm to 80 µm.The finest mesh is necessary to cover the narrow cylindrical plasma channel which is applied at the flame center (the axis) and spatial varied F = erfc(r 2 /a) b , where r refers to the radial distance from the discharge, and a and b are fitted geometric parameter taken from [26].The NRP discharges has a frequency of 10 kHz and deposited energy is fixed at 100 µJ per pulse.First-order Euler implicit scheme is used for time integration with time step varying from 0.2 ns to 0.1 µs.The enforced time step adjustment efficiently resolves the multi-timescale plasma-combustion physics observed during an NRP discharge period, which can be segmented into three stages as shown in figure 2(c): (1) nanosecond timescale of plasma energy deposition, (2) weak shock wave expansion from discharge kernel at few microseconds, and (3) laminar flame dynamics responding to previous plasma actuation at milliseconds.The domain is fed with a lean methane-air mixture (ϕ = 0.7) and nitrogen coflow with fixed velocity profile, as marked by the orange line in figure 2(a).Wall and free outlet boundaries are also marked in black and green, respectively.Non-reflective waveTransmissive boundary condition is used for the free outlets to allow the weak shock wave induced by the discharges to propagate outside the domain freely.A 2 mm step for methane/air inlet is included in the computational domain to match the inlet velocity closely and to reduced wave reflections.Details of the configuration are discussed in section 3.4.
Generating a stabilized base flame without plasma is the first step in the PAC simulation.Accurate inlet flow velocity profile is required for numerical setup to match flame morphology and position with experimental observations.Hot wire measurement were available, but was undertaken at  2 mm above the burner outlet due to the physical limitation.Therefore, the measured data cannot be used directly as inlet velocity in the simulation.Alternatively, the numerical inlet velocity was adjusted until the resulting velocity profile at 2 mm above the inlet matches the experimental one, as shown in figure 3. The skeletal mechanism with 30 species and 184 reactions for methane-air [38] is used and the default unity Lewis number assumption in reactingFoam is kept.A steady V shape base flame without plasma, as shown later in figure 4(b), can be generated and set as the initial condition for the simulation.

Results and discussion
A total of 250 pulses of discharges are applied to the flame in the simulations.The duration of plasma actuation, 25 ms, is sufficient for a new steady state of the flame to be established.One of the most prominent observables as an outcome of NRP discharges is flame displacement [31,32], which is well reproduced in this numerical study.To assess the roles played by individual plasma enhancement pathways, we consider five different cases by allocating discharge energy in specific directions.Table 1 lists the test cases, abbreviated as full-model, all-heat, all-disso, all-vib, and all-heat-20µJ.

Flame displacement
Figure 4 shows flame displacement in the full-model case by comparing planar laser-induced fluorescence (PLIF) images of CH from [32], compared with the computational prediction of the heat release rate (HRR) q.Several observations can be made.Firstly, the morphology and position of the base flame in the simulation satisfactorily match the experimental results.A slight discrepancy in the tailing edge section may be attributed to hot wire measurement errors in the low-velocity region around the nitrogen coflow, leading All Ėp goes into N 2 vibrational relaxation α = 0, η = 0 all-heat-20µJ Same as all-heat but only with 20 µJ per pulse α = 55%, η = 35% to corresponding inaccuracies in the inlet velocity profile within the 5-8 mm range.Secondly, the flame with NRP discharges stabilizes further upstream, exhibiting significant displacement (~2 mm) of the leading edge near the central axis.The magnitude of displacement agrees well with experimental observations.Lastly, the HRR increases across the discharge channel, indicating enhanced reactivity in this zone.
To identify the most important mechanism of plasma actuation affecting the flame displacement, different test cases listed in table 1 are conducted.The flame position, defined as the maximum HRR on the central axis, is compared in figure 5 as a function of time.Set as the benchmark, the black solid curve corresponding to the full-model case considers all plasma energy pathways and reproduces results in figure 4(b).A higher amount of flame displacement is observed in the all-disso case, where all deposited electric energy is spent on O 2 dissociation.In contrast, the flame barely moves if all input energy is used for vibrational excitation of N 2 , withstanding the subsequent slow heating effect (case all-vib).
The all-heat case exhibits a unique response, where the flame initially moves towards the burner but retreats to its original position with a small level of fluctuation.In this case, NRP discharges serve only as pulsed heating sources.The ultrafast energy deposition raises the temperature in the discharge channel and leads to weak shock wave propagation.The flame tip immediately responds to the enhanced reactivity at a timescale of combustion chemistry during the initial 40 pulses.After 5 ms, a potential hydrodynamic effect from gas perturbation tends to play a dominant role for the flame position.This is discussed in detail in section 3.4.

Chemical effect of NRP discharges
The flame position response in the all-disso case stands out among all test cases, indicating that a strong chemical effect from atomic O production is the most relevant factor for the displacement.This is supported by examining the radical production of O and OH in the full-model case.As shown in figure 6, the mass fractions of these two key reactive radicals are significantly enhanced near the flame leading edge, which coincides with the discharge channel.
Experimental results of PLIF images of CH and OH for different excitation lines also show a significant increase in florescence intensity in the flame front compared to the base flame [32].Since the reported OH excitation lines, labeled as OH O12(4) and OH S21(2), only represent a fraction of total OH in the flame, a direct quantitative comparison of OH distribution between simulation and the PLIF image is not possible.Despite this limitation, a qualitative comparison demonstrates a good consistency of OH enhancement and indicates a stronger reactivity with plasma actuation.

Effect of slow gas heating
HRR from relaxation of vibrational states of nitrogen N 2 (v) is calculated by where e eq vib (T) is the equilibrium vibrational energy at T and τ VT is the time scale of vibrational-transitional (VT) relaxation.The VT-relaxation process is dominated by oxygen atoms O( 3 P) [39]: τ VT is therefore determined by where n O is number density of atomic O.
Figure 7 shows the comparison of the temperature profile along the central y axis for different test cases.All cases involving energy deposition exhibit an increase in temperature compared to the base flame without plasma.Unlike all other PAC test cases, which present a hotter upstream fresh gas, the all-vib case does not increase the upstream flow temperature at all.This is expected since Ṙp VT is positively correlated with gas temperature and atomic O concentration, as can be inferred by inserting equations ( 9) and ( 10) into (8).Visualization of Ṙp VT in figure 8(case full-model) confirms that the slow gas heating effect from NRP discharges mainly enhances the temperature in the reactive leading edge  of the flame, where the highest concentration of atomic O from combustion chemistry resides in the discharge channel.Without the additional plasma-induced O 2 dissociation in the all-vib case, the slow heating process would be even more challenging and sluggish.

Hydrodynamic effect
The hydrodynamic effect of NRP primarily results in ionic wind or gas expansion leading to shock wave propagation [15,16,22].Since the mathematical model used in this study does not solve for the transport equations of electrons and ions as in [40], possible electrohydrodynamic phenomena cannot be assessed and is left for future work.
Fast gas heating on a nanosecond timescale increases temperature and pressure within the discharge channel, forming a weak shock wave that propagates at the speed of sound [22].In an ideal case without any reflection, the wave expansion quickly dissipates and the pressure field restores to its original state, as long as the timescale of the sound speed is smaller than the period of NRP discharges (1/f).This is illustrated in a test case, shown in figure 9(a), where 10 kHz NRP heating sources (100 µJ per pulse) are deposited at the center of an open square domain composed of stagnant nitrogen at 300 K and 1 bar.The waveTransmissive boundary condition set for all sides is able to filter out the circular wave induced by the high-pressure kernel.The pressure field around the heating source restores to a uniform value as early as 9 µs, and the entire computational domain completely returns to the initial 1 bar condition at 100 µs.As a result, the velocity field perturbation in the test case is negligible, even after multiple pulses.
However, wave reflection is unavoidable in a jet-wall PAC burner, not only physically from the top wall and inlet tube, but also numerically from the inlet boundaries, as marked in figure 1.This poses a challenge in the setup of the inlet boundary condition.A fixed velocity profile close to the flame ensures a well-positioned flame morphology, but it enhances wave reflection.Alternatively, extending the inlet tube in the computational domain can dampen wave propagation and therefore reduce its reflective effect from the inlet, with the downside of resulting in a less satisfactory flame shape.As a compromise, the methane/air inlet tube is kept 2 mm long in the computational domain.As seen in figure 9(b), reflective boundaries restrain the damping of the pressure propagation, resulting in a significant pressure imbalance at the end of a pulse period (100 µs).This computational artifact translates into velocity field modification, as presented in figure 10, where vertical velocities U y along the flame centerline are compared for different cases.A significantly higher U y in the upstream flame is observed, thus pushing the flame downstream even when reactivity can be increased from extra heat deposition.This explains why flame displacement is much reduced in the all-heat case.To further validate this explanation, the all-heat-20µJ case is set to inject only 20 µJ per heating pulse, an amount corresponding to the 20% fractional power for fast gas heating in the phenomenological model.Although being less heated, a higher flame displacement (see figure 5) and marginal upstream velocity modification can be achieved, simply because the hydrodynamic effect diminishes with much weaker wave propagation and reflection.Without reflection from physical inlet boundaries, it is expected that the hydrodynamic effect from shock wave expansion is even smaller in this PAC configuration.It should also be noted that gas expansion consumes energy and therefore reduces the effective thermal effect from heating source deposition.
The above results contradicts the findings of a previous study [41] which concluded that hydrodynamic effect of NRP discharges dominates over thermal and kinetic pathways in [41].In this experimental study using a similar PAC configuration, the authors argued that the NRP discharges act similarly to a flow blockage, producing a low-velocity region and inducing flame displacement.However, the mechanism of such a hydrodynamic effect was simply attributed to the ionic wind effect without clear validation.Moreover, the individual contributions from thermal and kinetic pathways were not examined in [41], making the claimed hydrodynamic dominance inconclusive.In contrast, the velocity decrease at the discharge vicinity resulting from fast gas perturbation was not observed in this numerical study.In summary, the potential significance of the ionic wind effect requires further investigation.

Flame power response to plasma
Figure 11 shows the flame power response in different cases after applying NRP discharges.The base flame without plasma input has a thermal power about 223 W, which well matches the reported flame power estimation of 230 W in [32].In all cases, there is power overshooting at around 8 ms (80 pulses) ranging from 10% (case all-vib) to 22% (case all-disso).After 250 pulses, the flame power restores the initial intensity of 223 W. Such a return to the quasi-steady state is expected since the plasma input energy is less than 1 W, corresponding to less than 0.5% of the thermal power of the flame.Furthermore, as the lean combustion is complete, the flame power cannot sustain the overshooting value based on energy conservation.This phenomenon was observed experimentally in Lacoste et al [42], where HRR fluctuations under a positive step of plasma actuation were reported.Note that the flame power is calculated by integrating the volumetric HRR q (W m −3 ) across the entire computational domain.Although the discharge zone is always enhanced, either thermally and/or chemically, the flame surface morphology changes and therefore affects the resulting integration.In other words, a power drop below the base value of 223 W does not imply that the plasma assistance vanished.
The highest amount of flame power surge in the all-disso case echoes the dominant role of the kinetic effect of NRP discharges in reactivity enhancement and flame displacement.The noticeable discrepancy between a decent power increment (14%) and a negligible steady displacement for the all-heat case further indicates that there must be an additional factor that suppresses the flame movement.This factor is found to be the numerically enhanced hydrodynamic effect from pulsed waves, as discussed in section 3.4.However, the flame power response suggests that even without the modeling limitation of wave reflection, the fast gas heating effect cannot overtake the chemical effect to be the primary plasma enhancement pathway in the PAC configuration under study.

Conclusions
The enhancement pathways of PAC for a laminar methane/air flame under NRP discharges actuation were investigated by computational simulations using a phenomenological plasma model.The study provided valuable insights into the complex interplay of thermal, kinetic, and hydrodynamic effects in PAC.
The results demonstrated that the kinetic pathway, particularly the production of atomic oxygen, plays a dominant role in flame displacement and reactivity enhancement.This finding disputes the conclusion of hydrodynamic effect dominance reported in [41], where the chemical effect of NRP discharges was not incorporated in detail.In contrast, our study showed that the sole kinetic pathway from atomic oxygen production is capable of significantly enhancing reactivity, modifying flame position, and boosting flame power fluctuation.
The thermal pathway, while contributing to the overall enhancement, was found to be less influential in this specific PAC configuration.The hydrodynamic effects, primarily arising from shock wave propagation, were shown to be sensitive to the specific burner geometry and the numerical setup of the simulation, particularly the boundary conditions.
The study also highlighted the challenges in isolating the individual effects of these pathways due to their strong interdependence.However, by carefully controlling the plasma energy pathways in the numerical model, it was possible to isolate and individually assess the plasma enhancement pathways.
This work contributes to a deeper understanding of the complex plasma-flame interactions in PAC and provides a guidance for further studies aimed at optimizing PAC systems for improved combustion performance.Future work should focus on further exploring the interplay of the enhancement pathways in different PAC configurations and operating conditions.

Figure 1 .
Figure 1.Schematic illustration of the plasma-assisted combustion (PAC) burner examined in this study.Reproduced from [31].© The Author(s).Published by IOP Publishing Ltd.CC BY 4.0.

Figure 2 .
Figure 2.An overview of numerical configuration in this study: (a) 2D axisymmetric computational domain with meshing and boundary conditions; (b) plasma region represented by spacial function F and covered by finest mesh; (c) adaptive integration timestep applied for a pulse period x 0.2 ns, y 1-50 ns and z 100 ns.

Figure 5 .
Figure 5. Flame position evolution as a function of time for different test cases.

Figure 6 .
Figure 6.Production of O and OH radicals in mass fraction for the full-model case.

Figure 7 .
Figure 7. Temperature along central axis for different cases at t = 25 ms.

Figure 8 .
Figure 8. Relaxation rate of the vibrational energy into gas heating in case full-model.

Figure 9 .
Figure 9. Pressure field following an NRP discharge.(a) Test case, heating source (100 µJ per pulse) deposited in quiescent flow filed without reflection.(b) Case all-heat with several reflective boundaries.Note that the annotation data range is not same.

Figure 10 .
Figure 10.Vertical velocity profile along central axis from difference cases.

Table 1 .
Abbreviations of test cases and the parameters of discharge energy pathway based on equations (1)-(4).