Machine learning-assisted early ignition prediction in a complex flow

Combustion and Flame 206 (2019) 451–466

Contents lists available at ScienceDirect

Combustion and Flame journal homepage: www.elsevier.com/locate/combustflame

Machine learning-assisted early ignition prediction in a complex flow Pavel P. Popov a,∗, David A. Buchta a, Michael J. Anderson a, Luca Massa e, Jesse Capecelatro d, Daniel J. Bodony c, Jonathan B. Freund a,b,c a

Coordinated Science Laboratory, University of Illinois at Urbana-Champaign, 1308 W. Main St., Urbana, IL 61801, USA Mechanical Science and Engineering, University of Illinois at Urbana-Champaign, 1206 W. Green St., Urbana, IL 61801, USA c Aerospace Engineering, University of Illinois at Urbana-Champaign, 104 S. Wright St., Urbana, IL 61801, USA d Mechanical Engineering, University of Michigan at Ann-Arbor, 1231 Beal Ave., Ann Arbor, MI 48109, USA e Aerospace and Ocean Engineering, Virginia Polytechnic Institute and State University, Randolph Hall, Blacksburg, VA 24061, USA b

a r t i c l e

i n f o

Article history: Received 26 December 2018 Revised 28 March 2019 Accepted 11 May 2019

Keywords: Ignition Plasma-coupled combustion Ignition prediction Machine learning

a b s t r a c t Machine learning methods are used to improve the efficiency by which turbulence-resolved simulations predict whether a hydrogen jet in air crossflow will successfully ignite. The flush-mounted jet issues perpendicularly from the wall of a low-speed wind tunnel into a turbulent boundary layer wherein a laser-induced optical breakdown (LIB) hotspot is deposited. A detailed hydrogen chemical mechanism is used to model the radicals and any subsequent chemical reactions. A dielectric-barrier discharge actuator generates body forces and hydrogen radicals near the jet orifice. We focus on the success or not of the ignition based on LIB location. A challenge is that definitive determination of this requires long simulations, up to 440 μs after the LIB deposition. This is particularly expensive since multiple simulations are required to find the threshold. To reduce the computational effort, three short-time (91 μs) criteria are proposed, evaluated, and compared: a constructed criterion based on detailed observations of radicals near the stoichiometric surface and two machine learning approaches, each trained on 38 realizations. The constructed criterion provides a low-cost estimate of the ignition boundary that is unambiguous in 45 of the 50 training and test trials, so only 10% of them would need to be simulated longer. The trained neural networks correctly predict outcomes in all cases evaluated, with the more automated procedure— a convolutional neural network (CNN) trained on two-dimensional images—providing the most definitive outcome prediction. From the CNN, a sensitivity analysis is used to determine which kernel features from the two-dimensional input data, as defined by intermediate-layer network weights, are important for identifying ignition. © 2019 The Combustion Institute. Published by Elsevier Inc. All rights reserved.

1. Introduction Plasma-coupled combustion presents numerous challenges for modeling and simulation [1], many of which introduce significant computational cost. We focus on the particular challenge of simulating ignition via laser-induced breakdown (LIB) and mapping the LIB focal point ignition boundaries. In principle, a new simulation is required for each focal point, making it expensive to map the ignition boundary as a function of the laser breakdown location. This is particularly true if there is statistical variation such as in turbulent flow. It is therefore attractive to determine the ultimate success or failure as soon as possible after the LIB deposition, and

∗ Corresponding author at: Coordinated Science Laboratory, University of Illinois at Urbana-Champaign, 1308 W. Main St., Urbana, IL 61801, USA. E-mail address: [email protected] (P.P. Popov).

avoid the significant cost of multiple long simulations to the point of a stabilized flame We consider multiple approaches for estimating successful ignition. A basic metric builds on the well-understood importance of radical precursors to ignition and quantifies their maximum values near the stoichiometric surface. Conceptually, it shares similarities with premixed minimum ignition energy criteria (such as those of Champion et al. [2], Frendi and Sibulkin [3], and Tromans and Furzeland [4]) and techniques such as Computational Singular Perturbation [5] and Chemical Explosive Mode Analysis [6]. The former a priori predict ignition for uniform quiescent mixtures, while the latter characterize the reaction region based on radical compositions. The other two approaches we consider take advantage of the more detailed description of the flow and chemical state available in the simulations and the flexibility of machine-learning procedures. These are trained on a limited number of long simulations

https://doi.org/10.1016/j.combustflame.2019.05.014 0010-2180/© 2019 The Combustion Institute. Published by Elsevier Inc. All rights reserved.

452

P.P. Popov, D.A. Buchta and M.J. Anderson et al. / Combustion and Flame 206 (2019) 451–466

that unambiguously identify whether ignition is sustained. Machine learning has been applied to numerous problems in computational thermo-fluid sciences, such as in the detection of combustion instabilities [7], control of industrial combustors [8], and identifying shocks in turbulent flow fields [9]. In these cases, procedures were crafted primarily for identifying features in data. In our case, we train the machine-learning network on outcomes, and interpret the network weights to understand the early-time features that correspond to these outcomes. To design and demonstrate the approaches in a non-trivial case, we select a configuration that involves multiple mechanisms: the plasma-actuated ignition of a hydrogen jet in a turbulent air crossflow. Plasma actuation in combustion is attractive for improving the likelihood of ignition and sustaining it in challenging environments [1,10,11]. Two different types of plasma applications are used here: a dielectric-barrier discharge (DBD) actuator [12], which generates body forces, hydrogen radicals, and heat in the vicinity of the fuel inlet, and a laser-induced breakdown (LIB) [13], which generates a radical-rich high-temperature ignition seed. The geometry of the present computational study—a jet in crossflow—is relevant to a variety of combustors. Details of the configuration, the simulation model, and its validation are reported elsewhere [12–15] and are therefore only briefly summarized in Section 3 The interaction between the kinetics of plasma generation within the LIB and the neutral chemistry of combustion is of particular importance, though its representation is simplified by significant time-scale sparation between LIB and flow. Nevertheless, the character of the LIB imposes a significant constraint on the choice combustion model: because the LIB region contains significant dissociated oxygen, locally exceeding the concentration of molecular O2 , the chemical mechanism must include O. The 9-species mechanism used is discussed in Section 3. We leverage the fact that the DBD time scales allow a decoupling of its plasma kinetics from the combustion kinetics time scales [12]. The primary goal is to determine the ignition boundary between LIB locations that will and will not lead to sustained ignition. This involves the LIB and DBD models as well as the fuel–oxidizer turbulent mixing, requiring relatively large-scale simulations to represent in detail. It is the expense of these that motivates our efforts to anticipate ultimate outcomes as early as possible in the simulation for any particular LIB placement. By reducing the number of full-length simulations necessary, such an early ignition criterion will speed up future ignition boundary searches. In pursuing this goal, at the outset we emphasize that our purpose is to identify early-time metrics of ignition. Rather, we seek to identify predictive metrics for ultimate ignition. The rest of this paper is organized as follows. In Section 2, we describe the configuration and our specific quantity of interest. A summary of the computational model and procedures, including a detailed assessment of grid resolution, follows in Section 3. Though the present effort focuses on early-time anticipation of the sustained ignition threshold, the behavior of the system leading to a stable flame or blow off is also relevant, and it is described in Section 4. The three ignition criteria are developed: in Section 5 as hand-constructed criterion is introduced for making comparisons, in Section 6 a neural network serves the role of a relatively simple nonlinear fitting function, and in Section 7 a more sophisticated convolution neural network is employed, which requires minimal user input. Though a drawback of this final approach is that it neglects knowledge of combustion, it is simpler to implement and, more interestingly, it can be interrogated to understand what specific features it hones in on that differentiate ultimate outcomes. The predictive capabilities of the models are compared in Section 8.

2. Experimental configuration We define a coordinate system centered at the fuel orifice, with the x-direction parallel to the air crossflow and the y-direction normal to the wall, aligned with the fuel jet (see Fig. 1(a)). Experiments are conducted in the 121 cm × 38 cm × 38 cm test section of a wind tunnel flowing air at either 5 m/s or 15 m/s. The fuel orifice of diameter 4.83 mm is flush with the bottom wall (y = 0) of UL

the test section. The Reynolds number is Re = ν jet where U is the air free stream velocity, Ljet = 0.356 m is the distance from the inlet of the test section to the fuel jet, and νair = 1.59 × 10−5 m2 /s is the kinematic viscosity of air. The boundary layer on the bottom wall is tripped turbulent with a strip of 40-grit sandpaper positioned 28 cm upstream of the jet (see Fig. 2). For the 15 m/s case, Re = 3.4 × 105 . The boundary layer turbulence is fully-developed at the fuel jet location [15]. The H2 fuel is injected at T = 298 K, with a volumetric flow rate of 2 L/min. The Froude number, based on  the fuel jet diameter and bulk velocity, is F r = ujet / gDjet = 5.2, so that inertial effects are dominant. The Damköhler number, based ch , and the scalar on the characteristic radical induction time, τind dissipation time, τ χ , both of which are estimated in Appendix B, is τχ Da = ch = 1.25, so that neither reaction or mixing are dominant, τind

as can be expected when we’re interested in marginally igniting flows. A DBD actuator (Fig. 1) is coannular with the orifice. Its exposed electrode is a thin-walled (0.51 mm wall thickness, 4.83 mm outer diameter) tube cylinder that is recessed 4.1 mm below the test section wall. The covered electrode is an annular disk recessed at a depth of 2 mm within the jet’s quartz dielectric material. The two electrodes are driven at 30 kHz, with amplitude of 11 kV. Complete details of the DBD design are reported elsewhere [14]. When the flow is fully-developed, with or without the DBD actuator active, a 39 mJ laser beam is focused at a chosen point, (x0 , y0 , z0 ), in the test section, initiating a laser-induced breakdown (LIB). For the primary simulation cases we discuss and the corresponding experiment, the breakdown is centered at z0 = 0 mm, y0 = 1.7 mm, and variable x0 , and it is slightly angled upstream −1.5◦ from the z-axis due to practical constraints in the design of the experiment. 3. Numerical method The configuration is discretized with three overlapping structured grids and shown in Fig. 2. A Cartesian background grid (grid 1) with 1068 × 257 × 97 points covers a 4.7 cm-wide spanwise section of the test section. The greater spanwise extent of the actual wind tunnel is modeled with periodic boundary conditions. The spanwise extent of the computational domain was determined sufficient based on the fact that the two-point spanwise autocorrelations of fluctuating velocity, between a point with z = 0 and a point with z = ±2.35 cm, are negligible [16]. The other two grids are aligned with the fuel orifice axis. Grid 2, with 266 × 385 × 129 points, is cylindrical, containing the jet walls and the DBD electrodes and grid 3, with 65 × 385 × 65 points, is a Cartesian square prism, overlapping the centerline of grid 2. Grid 3 is included to avoid both the coordinate singularity and the narrow azimuthal spacing of the cylindrical grid. Communication between the grids is performed via cubic interpolation in their overlapping regions. This grid resolution is sufficient to produce accurate mean flow and Reynolds stress profiles [15]. Additional resolution studies pertinent to the current ignition studies are provided in Appendix A. The mass, momentum, energy, and species transport equations are solved in generalized curvilinear coordinates for the conservative dependent variables

Q = [ρ , ρ u1 , ρ u2 , ρ u3 , ρ Etot , ρY1 , . . . , ρYN−1 ],

(1)

P.P. Popov, D.A. Buchta and M.J. Anderson et al. / Combustion and Flame 206 (2019) 451–466

453

Fig. 1. The DBD actuator. (a) a three-dimensional, isometric view, indicating the direction of primary wind-tunnel flow and the coordinate system, and (b) an x − y cross section, showing electrode dimensions, in cm.

Fig. 2. Overset meshes: (a) the full domain, including the trip strip, and (b) the fuel inlet and ignition region.

where ρ is density, ui are the velocity components, the energy Etot includes the kinetic energy, the sensible internal energy and the internal energy of formation, and Yi are species mass fractions. The mass fractions are constructed such that ρY1 , . . . , ρYN−1 are the densities of all species but N2 , with YN2 set by consistency. Finite-difference operators that satisfy the summation by parts (SBP) condition [17] are used for the spatial derivatives: a fourthorder, centered finite difference in the interior, and a second-order, one-sided finite difference on the boundary. The fourth-order Runge–Kutta scheme is used for time integration, with a time step of 3.6 × 10−8 s. Simulataneous-approximation-term (SAT) boundary conditions [18] enforce the no slip, isothermal T = 298 K boundary conditions. Following common practice, conserved variables are mildly low-pass filtered every two time steps [19] to maintain numerical stability without overly degrading resolution [20,21]. The low Mach numbers, with Ma = 0.015 and Ma = 0.044 for the 5 m/s, 15 m/s cases, would significantly restrict time steps. To

avoid this, the Mach number was increased to Ma = 0.3, with reaction rates, transport rates and body forces rescaled to match the Reynolds, Froude and Damköhler numbers. The only nondimensional group that changes is thus the Mach number, and flow with Ma ≤ 0.3 is widely accepted as effectively incompressible. Chemical reactions are modeled with a standard 9-species, 21reaction reduction of the UCSD chemical mechanism, retaining all of the species and reactions that do not involve carbon [22,23]. An ideal gas equation of state closes the equations, with heat capacities from NASA polynomial data fits [24]. Energy and differential species diffusion is implemented through a polyatomic gas extension of Chapman–Enskog theory [25–27]; Soret and Dufour effects are neglected. A modified Debye–Hückel model [28],

  ∂ ∂ξ ξ − ξ0 = , εr ∂ xi ∂ xi λ2N

(2)

454

P.P. Popov, D.A. Buchta and M.J. Anderson et al. / Combustion and Flame 206 (2019) 451–466

Fig. 3. The (a) time-averaged H radical mass fraction generated by the DBD, and (b) the streamwise body force at z = 0.

provides the potential, ξ , for the DBD-induced body force and radical source. In (2), εr = 4.65 is the dielectric constant [29], λN = 1.148 mm is a nominal screening length [15], and the potential offset ξ0 = 7764 V provides a model for the asymmetry of the DBD configuration, with a stronger negative ion sheath close to the dielectric [12], and it is calibrated to match measured body forces [30]. The potential ξ is pre-computed by solving (2); the time-averaged (over the AC cycle) body force and radical production are computed in a separate computation, as reported elsewhere [12]. The resulting steady-state radical concentration and body force fields are shown in Fig. 3. As we can see in Fig. 3(a), the average DBD-generated H radicals are about 104 times below the H concentration present in a typical reacting mixture of H2 and O2 [31], so it is expected at most to affect early stages of ignition. This is discussed in detail elsewhere [12]. Once the flow is statistically stationary, based on velocity fluctuations above the injector at y = 0 m, 0.0017 m and 0.005 m, the LIB hotspot is inserted at the desired target point. The LIB source comes from a separate finite-element solver [13], whose initial condition is a seed concentration of electrons, based on multiphoton ionization (MPI) [32]. These priming electrons absorb laser energy via an inverse Bremsstrahlung mechanism. The resulting plasma is modeled with standard flow equations coupled with a two-temperature thermal non-equilibrium model. The plasma chemistry is modeled with an 11-species plasma kinetics model for air, as reported elsewhere [13]. At 6.3 μs after the 23 ns-long laser pulse starts, the temperature has decayed to T = 4400 K, within the range of applicability of the combustion chemistry model [23], and thermal equilibrium has been achieved. Over the same 6.3 μs, the velocity, density, internal energy, and radical mass fractions of the overset mesh combustion simulation are relaxed to those of the LIB finite-element method solution in the vicinity of the focal point [15], where the temperature is higher than 320 K. This region is fully within a sphere of diameter 1.5 mm, centered on the focal point. The relaxation implicitly assumes that the coupling between the LIB and turbulent combusion solutions is one-way, which is justified due to the timescale separation between the two. Specifically, the largest time scale of the LIB solution—the period from the laser firing to the temperature decaying down to 50 0 0 K, when the San Diego mechanism becomes applicable—is lower than 6.3 μs, whereas, as we shall show subsequently, the smallest time scale of the turbulent combustion solution is the radical induction time, at τind = 30 μs. Following Massa et al. [15], a model for vorticity generation due to shock curvature in the LIB hotspot is also included, which leads to an additional model Mach number parameter, M0 = 1.61, for the strength of the curved shock in the laser direction [15] (Fig. 4)

4. Characteristics of sustained ignition Before we analyze the ignition threshold, we describe the flame evolution for successful sustained ignition. In these cases, the flame attaches to the upstream lip of the fuel orifice, with its priming reaction zone aligned with the stoichimetric surface above the orifice. This is observed for all cases with sustained ignition, regardless of the LIB placement. In Fig. 5(c), a corresponding thin region of O radicals can be seen, anchored on the front lip, with a H concentration just below it. Several comparisons with corresponding experimental data were reported previously for a similar configuration though with different LIB locations. Massa et al. [15] showed that calculated mean velocity profiles and Reynolds stresses agree with experiments. It was also shown that with the addition of the LIB vorticity generation model, the predicted ignition boundary lies within the uncertainty interval associated with the measured location of laser breakdown. Massa and Freund [12] used more detailed plasma chemistry for cases with DBD actuation, showing an improved prediction. Finally, for the present configuration, Popov et al. [33] make comparisons between computed and experimentally measured H and O radical concentrations. To bracket the ignition threshold position in x for y0 = 1.7 mm, z0 = 0, simulations were run for several values of x0 , until they could be classified as either igniting or non-igniting. This was determined based on the time evolution of the maximum temperature Tmax (t), such as shown in Fig. 6(a). The temperature first decreases due to thermal conduction. In cases that ignite, after 200 μs to 400 μs the heat release of combustion leads to Tmax > 20 0 0 K and a sustained flame. Alternatively, there is extinction, for which Tmax monotonically decreases below 900 K (see Fig. 6(a)). At this temperature, ignition must fail, since the autoignition time of a stoichiometric mixture becomes greater than 10 ms [22]. This considerably exceeds the mean strain time scale, τconv = ( ddyU¯ )−1 , which is τconv = 8.5 × 10−4 s for the 5 m/s mean flow case, and it is still faster for 15 m/s. The scalar dissipation time scale is even faster, per the discussion in Section 5. Using the bracketing procedure in the above paragraph, we obtain the simulated ignition thresholds in agreement with the experiments [33]. The end-to-end uncertainty quantification performed in conjunction with these simulations is not included here but required significant sampling, which was the primary motivation for finding early ignition criteria. In most cases, agreement was within the uncertainty bounds, though these bounds themselves were relatively large. The spacing uncertainties involved details of the LIB dynamics, particularly the LIB-induced hydrodynamics, which are complex [34,35] and beyond the scope to study here. Improving the physical representation of the LIB and accelerating corresponding

P.P. Popov, D.A. Buchta and M.J. Anderson et al. / Combustion and Flame 206 (2019) 451–466

455

Fig. 4. Visualization of the statistically stationary turbulent boundary layer prior to ignition. The model sandpaper trip and fuel jet are shown. An isosurface of the middle eigenvalue of the velocity gradient tensor, colored by the axial velocity, indicates the presence of vortices. The wall boundaries are indicated in dark grey, and the equivalence φ = 0.1 isosurface is shown in bright pink. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 5. Sustained flame visualization for a 15 m/s, no-DBD case: (a) stoichiometric surface, colored by temperature; (b) H radical mass fraction in the x − y at z = 0; and (c) the corresponding O radical mass fraction. In (b) and (c), the white line indicates TALIF measurement locations. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

simulations is another primary motivation for the proposed means of early ignition detection. Examples of the evolution of Tmax (t) for igniting and nonigniting simulations can be seen in Fig. 6(a). There are three alternatives for temperature evolution: asymptotic decrease to a high temperature, decrease to a local minimum followed by an increase to a high temperature, and decrease to below 900 K.

5. Constructed early ignition indicator 5.1. Specification The radical induction time sets the early period of the ignition process [36], during which radical concentrations increase with little heat release. This suggests in indicator based on radicals.

456

P.P. Popov, D.A. Buchta and M.J. Anderson et al. / Combustion and Flame 206 (2019) 451–466

Fig. 6. Evolution of the global maxima of temperature and radical concentrations (normalized against the STP density of air, ρ air,0 ) for the 15 m/s, no-DBD ignition case. Red indicates long-time confirmed successful ignition per the discussion of Section 4. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

For the present chemical mechanism, the temperature and radical species H, O, OH, HO2 , and H2 O2 evolve as shown in Fig. 6. The three simulations closest to the ignition boundary from the bracketing procedure show the global maxima increasing as the LIB approaches the stoichiometric surface. Though the temperature evolves almost identically for the first 100 μs, the radical maxima differ as early as 18 μs after LIB, especially for H. This behavior, however, is not shared by all of the simulation cases (5 m/s and 15 m/s; DBD and no-DBD). For example, it differs significantly for the 15 m/s, DBD-on cases, for which the hydrogen concentration threshold between igniting and non-igniting simulations is 3.7 × 10−5 kmol/m3 (Fig. 6(b)). In order to generalize across simulations, we seek further separation between igniting and non-igniting simulations. Since the ignition delay time is typically lowest near φ = 1 [36], we focus on the metric near the stoichiometric surface: φ ∈ [0.84, 1.16]. Figure 7 confirms that this provides cleaner distinction between igniting and non-igniting cases than the global maxima, and the results are insensitive to the specific neighborhood. Taking φ ∈ [0.75, 1.25] showed no significant differences. These results suggest that an indicator based on the time dependent values of ρ (x,t ) maxφ ∈[0.86,1.14] ( ρX ), for some choice of radicals X, will be efair,0 fective.

To include multiple radical species, we considered both addition and multiplication of their contributions. Multiplication was chosen because it was more straightforward to include the large variation in typical concentration values for different radicals. The indicator function

f (t ) ≡



 X∈R

max

φ ∈[0.86,1.14]

ρX ( x , t ) ρair,0



(3)

accomplishes this, with R being some subset of {H, O, HO2 , H2 O2 , OH}. Different R were tested, with the best found, as quantified subsequently, to be R = {HO2 , H2 O2 }. Finally, due to the large variation in radical concentration levels between igniting and nonigniting cases, we choose to take the natural logarithm of (3):

 c (t ) ≡ log

 max

φ ∈[0.86,1.14]

   ρHO2 (x, t ) ρH 2 O 2 ( x , t ) × max . ρair,0 ρair,0 φ ∈[0.86,1.14] (4)

5.2. Results for the constructed criterion The constructed criterion is designed to anticipate successful ignition for the 38 simulations of Section 4 (19 igniting and 19 nonigniting) closest to the ignition boundary. These include 15 m/s and

P.P. Popov, D.A. Buchta and M.J. Anderson et al. / Combustion and Flame 206 (2019) 451–466

457

Fig. 7. Evolution of the maxima of temperature, and maxima of radical concentrations for φ ∈ [0.84, 1.16], for the 15 m/s, no-DBD ignition case.

5 m/s cases, DBD and no-DBD cases, and y0 = 1.7 mm, 1.1 mm, and 0.562 mm cases. Figure 8 shows the metric (4) for these cases. We define a nominal overlap interval based on how c(t) cases fail to distinguish outcomes:

  (t ) ≡ min c(t ), max c(t ) . ign.

(5)

non-ign.

In making a prediction for a new case, we take the upper boundary of (t) as a threshold for a confident prediction of ignition, and likewise the lower boundary of (t) as a threshold for predicted non-ignition. As such, (t) should be designed to be as narrow as possible. The species selection that led to the specific form (4) is based on minimizing the relative overlap, RO(t), defined as





max c (t ) − min c (t )

RO(t ) =

non-ign.



ign.

max c (t ) − min c (t ) all

.

(6)

all

On average, (t) decreases in time, from  (18 μs) = [−43, −30] to  (383 μs) = [−20.3, −18.3] . As expected, there is a tradeoff between accuracy and utility. Accuracy is low at t = 18 μs and high at t = 383 μs. However, there is less benefit for longer t, since the cost approaches that of running the full simulation. We take t ∗ = 91 μs, which is the highest output time that is still less than a quarter of a complete simulation time interval. This gives us a

Fig. 8. Ignition indicator c(t∗ ) from (4) for the simulations close to the ignition boundary (see text). Igniting cases are colored red and non-igniting cases are colored blue. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

balance between accuracy and utility—in Fig. 8, we see that all igniting simulations satisfy c (t ∗ ) > −25, and all non-igniting simulations satisfy c (t ∗ ) < −29. Only 4 of 38 cases have c(t∗ ) between these values, which would need to be run longer to conclusively determine the outcome. The time t∗ is comparable to the characteristic mixing time τ χ ≈ 40 μs and the characteristic reaction time

458

P.P. Popov, D.A. Buchta and M.J. Anderson et al. / Combustion and Flame 206 (2019) 451–466

Fig. 9. Two ignition boundary searches, visualized with equivalence ratio, φ : (a) one turbulence realization with full strength LIB, (b) another realization with 65% strength LIB, and (c) a schematic of the intervals of predicted ignition and non-ignition and the actual ignition boundary for both cases. The actual ignition boundaries are within the overlap interval of the c∗ criterion. ch ≈ 30 μs (Appendix B). It should be noted that, while reaction, τind

diffusion and mixing are not explicitly included in the c(t∗ )-based criterion, or the machine learning criteria in the next two sections, each of these critera include data for the increase of radical concetrations over a fixed time period, and thus reaction, diffusion and mixing are accounted for implicitly. 5.3. Generalization We evaluate the portability of the criterion for other parts of the ignition boundary in this same configuration, with different strain rates, mean velocities, mixedness and scalar gradients, and select two different turbulence realizations. The LIB is taken to lie along x0 = 2.4 mm and z0 = 0, and we seek the ignition threshold in y0 . For the first of these realizations, illustrated in Fig. 9(a), a series of 11 simulations up to t ∗ = 91 μs was performed, with y0 varying in even increments of 0.4 mm from y0 = 1.1 mm to y0 = 5.1 mm. All of the cases with y0 ≤ 2.1 mm have c (t ∗ ) > −25, and all of the cases with y0 ≥ 2.5 mm have c (t ∗ ) < −29. Thus, the ignition bracket, based on the c(t∗ ) criterion, is y0 ∈ [2.1, 2.5] mm. To confirm and further refine this bracket, 5 longer simulations were run, at intervals of 0.1 mm between y0 = 2.1 mm and y0 = 2.5 mm. Those with y0 ≤ 2.3 mm all ignited, and the one with y0 = 2.5 mm did not. Therefore, the uncertainty interval for the ignition boundary (defined as the interval between the last igniting simulation and the first non-igniting one) is wholly inside the uncertainty interval of the criterion (−25 < c (t ∗ ) < 29), and hence no false predictions are made for this parameter search. This validates the c(t∗ )-based criterion for this related configuration. A second search used LIB energy reduced by 35%, with the results seen in Fig. 9(b). Despite these changes, the early ignition criterion successfully refines the ignition boundary down to a 0.4 mm-wide interval, which is further refined down to a 0.2 mm interval by three long runs to confirm that the criterion is correct. 6. Simple machine learned ignition criterion Though effective, the constructed c(t∗ ) criterion (4) required detailed examination of the simulation results, and risks also being specific to the particular chemical mechanism and configuration.

For a more complex chemical mechanism or more complex configuration, such an approach may even become intractable or good indicators might be missed. To search for criteria more broadly and avoid the tedious task of seeking them by hand, we use a machine learning procedure based on a neural network [37] trained to take information from the simulation flow fields and provide a scalar ignition indicator. In this section, we consider a simple approach using a two-layer network (the simplest type which can approximate a general non-linear function) with hand-designed inputs. The following section employs a more sophisticated convolutional network to receive more general inputs. As input, we form a 15-dimensional vector of the maxima of all minor species concentrations (H, HO2 , H2 O2 , O, and OH) for φ ∈ [0.84, 1.16] at times t = 0.6t ∗ , t = 0.8t ∗ and t = t ∗ . All are included to afford more flexibility than could be identified via the inspection process that led to the criterion of Section 5. For training, the first attempted target was a binary function: 0 for non-igniting cases and 1 for igniting cases. The best performance achieved had 4 undetermined outcomes of the 50 training and test simulations, which is not an improvement over the constructed criterion. To improve this for the case of ignition anticipation, we regularized the ignition pointer, I p (S ), which depends on the simulation S. It was designed to be: 1. a smooth function of the simulation S, 2. change monotonically (increasing by choice) with its input as S changes from less igniting to more igniting, 3. well-conditioned, on the order of 1 [38], to reduce the number of iterations before the NN starts converging with standard hyperparameters and initialization algorithms, and 4. sensitive for simulations S near the ignition boundary. We take



Ip =

Cnorm

25

Tmax (tk )e

Ce tk

Cpwr − Coffset

(7)

k=1

based on conditions 1 to 4 above. For each simulation in the training set, we have 25 solution fields, evenly spaced at t∗ /5 intervals

P.P. Popov, D.A. Buchta and M.J. Anderson et al. / Combustion and Flame 206 (2019) 451–466

so tk = kt ∗ /5. Since Ip is only evaluated for simulations in the training set (ignition predictions will be made based on the predicted Ip ), we use all 25 stored realizations. Per conditions 1 and 2, a weighted sum distinguishes between simulations that ignite at different times but similar eventual temperature maxima (e.g., the two igniting cases in Fig. 6(a)). The factor eCe tk is chosen to satisfy condition 2: per its first part, the earlier temperatures should affect Ip , but its second part makes the eventual outcome more significant. The value of Cnorm is set so that the quantity in the square brackets has a maximum 1 for the training set, per criterion 3. Finally, the offset Coffset and Cpwr = 1/3 jointly serve to satisfy condition 4: Coffset ensures Ip > 0 for ignition and Ip < 0 for non-ignition, and the most marginal igniting and non-igniting cases are approximately equidistant from I p = 0, so that when we take the 1/3 power the sensitivity of Ip for marginal cases is increased. The net effect of satisfying condition 4 through Coffset and Cpwr is that the trained network has one less undertemined test case, thus their inclusion is justified. Specific values are Ce = 1.1 × 105 s−1 , Cnorm = 2.2 × 10−13 K−1 , Cpwr = 1/3, and Coffset = 0.24. We use a two-layer NN configuration, which is common for simple-function estimation. It includes a hidden layer with 15 tansig neurons [37], which yielded optimal performance—fewer neurons result in a worse fit of the training data (higher mean-squareerror between true and predicted ouput values), whereas more neurons result in more overfitting (higher mean-square error on the validation and test data). Starting from the inputs pi , for i = 1, 2, . . . , 15, the first (hidden) layer’s outputs are,





(1 )

aj

15   ≡ tansig wi,(1j) pi + b(j1) ,

j = 1, 2, . . . , 15,

(8)

i=1

where i and j are respectively the indices of the input vector and hidden layer, and wi,(1j) and b(j1 ) are the corresonding weights and offsets. The output layer approximates Ip in (7) as

INN p

15 ≡ w(j2) a(j1) + b(2) ,

(9)

j=1

where w(j2 ) are the weights of this layer and b(2) is the linear offset. To train the network, to find wi,(1j),(2 ) , b(j1 ) and b(2) , we use the Levenburg–Marquardt algorithm with Bayesian regularization [39,40]. This means that the learning process aims to minimize a loss function F that combines the error and network weights

F = β ED + α EW , where

ED =





(10)

I p − INN p

2

(11)

samples

and



EW =

{

(1 )

w2 . (2 )

w∈ wi, j ,w j

(12)

}

The weights α and β in (10) are positive constants, ED is the sumsquare error, and EW is the sum-square of the weights. Bayesian regularization automatically adjusts α and β during training [40]. The implementation uses MATLAB’s neural network toolbox [41]. The character of the ignition process presents a particular challenge to direct application of standard neural-network tools. Its root is the explosive growth of radicals and the underlying nonlinearities in the kinetic mechanisms that lead to this behavior. This is obvious in Fig. 10(a), which shows a step-like dependence of Ip on one of the radical concentrations. A consequence of this for

459

training a network is that the gradients in the two regimes (igniting and non-ignition) are not comparable. Standard training procedures depend on this [37]. They fail in the regularization procedure of the weights. In the low-gradient portion of the range of behaviors, the weights’ contribution to (10) is negligible, which leads to overfitting in that region. The weights end up being slavishly linked to local gradient fluctuations, as demonstrated in Appendix C, and are therefore ineffective in extrapolation. A logarithmic regularization (Fig. 10(b)) reduces sensitivity to the exponential growth of radicals during ignition. This improves performance—the spread seen in Fig. 10(c) (which shows the correlation of the actual and predicted Ip ) is removed in Fig. 10(d). To estimate the uncertainty in the machine learning predictions, and to reduce its variance, we use the “bagging” ensemble method [42], which is based on independently training the network weights several times (here, 10 times), with different random initial weights, sample order, and training/validation split. The final prediction is the average of the predictions from the ensemble:





N N :k INN , p = mean I p

k = 1, 2, . . . , 10,

(13)

N :k denotes the prediction given by the k-th set of network where IN p weights in the ensemble. The variation within this ensemble provides an estimate of the 95% confidence interval for the ensembleaveraged prediction. NN is From the mean estimate, INN p , ignition is predicted if I p above 0.17 with 95% confidence, and non-ignition is predicted if INN is similarly below −0.195. If neither are satisfied we continue p the simulation. Results will be compared with other indicators and discussed in Section 8. For now, we compute the sensitivity of INN to perp turbations in each of the pi , to assess which of the network’s inputs affect the prediction most. The sensitivities were computed ∂I with a simple second-order finite difference approximation for ∂ pp . i

Figure 11 shows these sensitivities for 4 of our training cases. The sensitivities for the marginal cases are greatest for t = 0.8t ∗ and t = t ∗ . We also see that all radical species contribute to the NNtrained indicator. We also note that HO2 has the greatest overall sensitivity at t∗ indicating that it was a good choice for the constructed criterion. However, the other, H2 O2 , has the greatest sensitivity for the marginal simulations at the earlier time, 0.6t∗ . Additionally, we can see that OH is the least sensitive. For marginally igniting cases, OH concentrations increase later, as can be seen in Fig. 7(e). The peculiar negative sensitivity to O is a consequence of the LIB. The LIB dissociates most of the oxygen. This O reassociates but is also consumed faster as the reactions progress. Thus, less O at later times is consistent with more chemical reactions leading to ignition. Though successful, this standard ML approach for fitting a nonlinear function has some significant deficiencies. Most notably, the dependence on a user selected ignition pointer Ip is a limitation, especially since it requires the appropriate value for Ce , Coffset and Cnorm based on an analysis of the training set, and also to run every simulation in the training set for at least 5t∗ . To avoid these limitations, in the next section we deploy a convolutional neural network (CNN), whose performance with a binary target (ignition vs. non-ignition) is comparable without the ignition pointer Ip . It can be viewed as having sufficient flexibility to determine an effective equivalent of Ip . This reduces both the computational and human time required for the development of a predictive criterion, and it should be more portable to related applications. In essence, the current simple NN can be considered to have a hybrid character— a computer-trained NN with an essential hand-designed component (the ignition pointer, Ip ). This places it midway between the machine-learning reliant CNN and the previous hand-designed

460

P.P. Popov, D.A. Buchta and M.J. Anderson et al. / Combustion and Flame 206 (2019) 451–466

Fig. 10. Input regularization: (a) raw concentration inputs, and (b) their logarithms, and corresponding correlations in (c) and (d) respectively.

support whose parameters are trainable. Pooling layers collapse neighboring values into one. As an example, here we use two-dimensional fields for input, with 2 × 2 convolutional and max pooling layers (average pooling layers were also tested, yielding slightly higher validation error). In ) general, the input of the z-th layer, ai,(zj,k , is an nc -valued field over

an nx × ny domain and has the shape nx × ny × nc . Then, for a 2 × 2 convolutional layer with nf filters (output fields), the output is



(z+1 )

ai, j,k

= ReLU bk +

l,m∈[0,1],n∈[1,2,...,nc ]



(z )

ql,m,n,k ai+l, j+m,n

(14)

for i = 1, 2, . . . , nx − 1, j = 1, 2, . . . , ny − 1, k = 1, 2, . . . , n f , where ReLU(x ) = max(0, x ) is the rectified linear unit function. For a 2 × 2 max pooling layer the output is Fig. 11. Sensitivities of INN to the network inputs. p

criterion. In terms of how much of the prediction task it can take on, there is no expectation that it will perform so well as the following CNN. 7. Convolutional neural network indicator Convolutional neural networks (CNN) are widely used for image recognition and have been used for post-processing [7] and identifying features [9] of combustion applications. A CNN has a deep structure (i.e., multiple layers), which can represent nonlinear functions more accurately. It can operate at similar computational cost on larger data sets and is well-suited to recognize features (spatial correlations). Its output is also translationally invariant [43], consistent with the translational invariance of physics. Convolution layers convolve the input with a function of small

) ai,(z+1 = j,k

max a2(zi+) m,2 j+n,k

m,n∈[−1,0]

(15)

for i = 1, 2, . . . , nx /2 , j = 1, 2, . . . , ny /2 , k = 1, 2, . . . , nc . Successive convolutional and max pooling layers identify spatial features; at the same time, the intermediate output becomes translationally independent of the input [44]. This loss of information is embodied in the max pooling layers, which do not distinguish the origin of the max entry within the 2 × 2 stencil [43]. Such networks are relatively computationally tractable, because they reduce nx and ny for successive layers (note that the range of the output indices i, j in (14) and (15) is lower than that of the input indices), and they have a low number of trainable parameters— zero in the case of max pooling, and (4nc + 1 )n f for the 2 × 2 convolutional layer—which do not scale with nx and ny . In contrast, dense fully-connected layers (such as those used in the last section) with nx × ny × nc input and output fields would have approximately n2x n2y n2c trainable parameters.

P.P. Popov, D.A. Buchta and M.J. Anderson et al. / Combustion and Flame 206 (2019) 451–466

461

Fig. 12. The convolutional neural network used.

Our CNN takes as inputs 15 two-dimensional fields at the center plane z = 0, covering x ∈ [−0.75 cm, 0.75 cm], y ∈ [0, 0.5 cm]. While most widely used successfully for image analysis, CNNs can also be applied to three-dimensional input data [45,46]. We have opted for 2D input fields for the sake of minimizing the computational cost of training and evaluating the CNN criterion, and it will be shown sufficient for this purpose. Scaling up to three-dimensional inputs in the future is not expected to be a problem, since many networks have been trained to operate on inputs of significantly higher dimension than the present 150 × 50 × 15 [47]. The CNN output will be a sigmoid classification ICNN ∈ [0, 1]. The CNN is implemented in TensorFlow, using the Keras API [48]. Details of these algorithms are well-documented [48,49]. The network is trained using binary crossentropy [50] as the loss function F

F = −[Itrue × ln (ICNN ) + (1 − Itrue ) × ln (1 − ICNN )],

(16)

where Itrue is the true outcome. The specific training algorithm is RMSProp [51] with a learning rate (similar to the step in a fractional Newton’s method, see [37]) of 0.001 and a decay rate (which controls learning rate decay, see [51]) of 0.02. Denoting a full loop over all of the training and validation data as an epoch, each training sequence consists of 150 epochs and takes 120 s on a Core i7 modern CPU. After training is complete, we accept the weights of the epoch with the lowest validation loss, provided its value is below 0.1, for use as our predictor. Since every predicion with binary crossentropy loss has the same validation error function, the cutoff value of 0.1 used here can be used in general. Figure 12 shows the layer-by-layer structure of the CNN used, following Keras terminology [48]. The input layer accepts the 150 × 50 × 15 input array, which is the discretized two-dimensional field for 15 variables. The input passes through three convolution and pooling layer pairs, each of which consists of a 2 × 2 convolutional layer with a ReLU activation function (and 16,32 and 64 filters, respectively), and a 2 × 2 max pooling layer. Next, we flatten the resulting 17 × 5 × 64 tensor (into a 5440dimensional vector) and add 2 fully-connected dense layers with 512 units each and ReLU activation. Each of these dense layers is followed by a dropout layer with a dropout coefficient of 0.85. Dropout, which entails randomly disabling a certain fraction (in our case, 85%) of the neurons in a layer, has been shown to counter overfitting [52]. In the present case, the dropout keeps the validation error decreasing together with the training error

for a larger number of epochs, up to about 150, versus about 20 without dropout. This can be seen on Fig. 13—with dropout, the validation loss decreases with the training loss until around the 120-th epoch, and attains a minimum value below 0.1, whereas without dropout the network starts to overfit (validation loss increases while training loss increases) at around 20 epochs, and its lowest attained value is ≈ 0.2. The final layer is a single-unit, fully-connnected layer, with a sigmoid activation function, which results in a single-valued output ICNN ∈ [0, 1]. Choosing effective input variables is vital. Based on the sensitivity study in Fig. 11, we choose to only use data from t = 0.8t ∗ and t = t ∗ . Thus, 10 of the input variables have the form





Inpl = log10 Yα fφ (t ) + 



,

l = 1, 2, . . . , 10.

(17)

For improved performance, we add fields tracking the increments between 0.8t∗ and t∗ , based on the principle that not only the radical concentrations, but also their rate of increase, may be a strong predictor of ignition:









Inp11–15 = log10 Yα fφ (t ∗ ) +  − log10 Yα fφ (0.8t ∗ ) + 



. (18)

Adding the fields in (18) is called feature engineering [38]—they are not independent of the first 10 fields, but their addition improves the validation accuracy of the trained network. As we shall see in a sensitivity analysis performed at the end of this section, the trained network is just as sensitive to (18) as to (17). The inclusion of the horizontal and vertical velocities in addition to the above 15 input fields was also tested, and determined to not provide significant improvement of the predictive accuracy. Thus, we proceed with the radicals only. Sensitivity of the network to the input data is tested by training and evaluating it on alternative data sets in which all pointwise radical concentrations are multiplied by a normally distributed random factor F ∼ N (1, σN2 ). The network maintained the same performance for values of σ N up to 0.075. Given time-developing character of the ignition kernel, it would be possible to employ a neural network architectures that is designed for time series, such as Long-short Term Memory (LSTM) [53] or Recurrent Neural Networks (RNN) [54]. The reasons behind this decision are threefold. Firstly, both LSTM and RNN are usually applied to longer time series [55,56] than the 2 (in this section) to 3 (in the previous section) time points considered here. Secondly, LSTM and RNN are usually applied to series with a large degree of translational symmetry, whereas it is unclear that the

462

P.P. Popov, D.A. Buchta and M.J. Anderson et al. / Combustion and Flame 206 (2019) 451–466

Fig. 13. Comparison of the training and validation losses during training of the CNN: (a) using dropout 0.85, and (b) without dropout.

Fig. 14. An example of the convolutional network’s 15 input fields, for a marginal igniting simulation. The region plotted is x ∈ [−0.75 cm, 0.75 cm], y ∈ [0, 0.5 cm], z = 0. The magenta line shows the LIB axis, and the magenta dot shows the initial LIB location. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

present ignition case is fundamentally time transient, given the known time for the laser pulse. Finally, the present simpler networks are fully successful for this configuration; other configuration or other prediction objectives might warrant more specialized networks. In (17) and (18) above, α ∈ {H, HO2 , H2 O2 , O, OH }, t ∈

) {0.8t ∗ , t ∗ },  = 10−9 and fφ = exp(− (φ2−1 ), with df being a d2 2

f

parameter that controls weighting of the fields near the stoichiometric sufrace, similarly to the procedure in Section 6. The parameter  is used to avoid the singularity of log10 at 0. The value  = 10−9 was chosen to be at least 10−3 times the lowest radical mass fraction maxima for marginal simulations at t = t ∗ , and the value d f = 0.53 was determined to produce the lowest error on the validation subset of the training set. For a broad range of these values, the results are insensitive to the specific choice of  and df , giving similar conclusions as in Section 8 for  = 10−8 and df ∈ [0.3, 0.8]. Example input fields for a marginal igniting case are shown in Fig. 14. Outputs of the three convolutional layers are given in Fig. 15. We can see that, with successive convolutional layers, the outputs isolate either the interior or the edges of the field: after the first convolutional layer, there remain contributions both of the interior and edges of the kernel, whereas by the third all are focused on one or the other. Note also that all three layers contain trivial (zero) filters, as a result of the regularization, which reduces the weights of uninformative neurons.

This split between interior and edge components of the flame kernel suggests that the network output depends on multiple qualitatively different metrics. To further explore this split, we compute the sensitivities of ICNN with respect to each of the 64 output filters of the third convolutional layer. Each sensitivity with respect to a filter is a single-valued scalar computed essentially by a finite difference operation: the entire filter field is proportionally increased a small increment and the change in outcome is quantified. The results are shown in Fig. 15(c); various alternative marginal cases have been interrogated, with similar results. Figure 15(c) shows all 64 filters ordered in order of increasing sensitivity, from negative to positive. We can see that almost all filters to which ICNN has a positive sensitivity identify edge-type features, while the majority of negative-sensitivity filters have interior-type features. This can be seen more clearly in Fig. 15(d), which shows the sparial gradients of the squares of the filters, for those filters to which ICNN is most positively and most ngatively sensitive. The negative filters have their largest gradients in the interior of the kernel, and the positive filters have their largest gradients on the kernel’s edges. This can be interpreted in terms of the different shapes of igniting and non-igniting kernels. Since the LIB is fired away from the fuel jet, an important feature of an igniting kernel is radical production at the region of the kernel which is closest to the stoichiometric surface: this region will be on the kernel’s edge. On the other hand, the evolution of a non-igniting kernel more closely resembles a decaying Gaussian profile, similar to the filters to which ICNN is negatively sensitive to. The sensitivity ordering in Fig. 15(c)

P.P. Popov, D.A. Buchta and M.J. Anderson et al. / Combustion and Flame 206 (2019) 451–466

463

Fig. 15. Representative CNN output fields after the (a) first and (b) second convolutional layer of Fig. 12, and (c) the full third convolutional layer, ordered by increasing sensitivity of ICNN to each filter (the sensitivity value is shown above each plot), and (d) amplitudes of the gradients of the squares of the filters, for the 8 filters to which ICNN is most negatively sensitive, and the 8 filters to which ICNN is most positively sensitive. As these are intermediate layers of the CNN, the colormap ranges are not significant, and are omitted for better figure visibility. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

indicates that the CNN uses this distinction as part of the learned ignition prediction criterion. Finally, before we quantify the predictive performance in the next section, we compute the sensitivity to perturbations in its inputs to quantify which input fields the network are most important. Figure 16 shows the sensitivity of ICNN with respect a rescaling of each of the 15 input fields—pointwise sensitivities were also computed, but due to the max pooling layers these sensitivities are noisy and uninformative. We calculate sensitivities for the same 4 simulations as those for the simple NN in Fig. 11. The marginal simulations exhibit the greatest sensitivity. Also of note is that the network is most sensitive to both species in the constructed criterion, HO2 and H2 O2 . Additionally, we can see again that OH has the lowest sensitivity, similarly to last section’s findings. We can also see that almost all sensitivities are positive, whereas some of those in the last section are negative. This can be explained by the

fact that the CNN accounts not only for the maximum of a species’ concentration, but its extent, which is generally larger for more strongly igniting cases. Additionally, we have similar sensitivities for t = 0.8t ∗ and t = t ∗ , whereas for the simple NN the latter sensitivities were considerably higher. This can be explained by the additional criteria (such as the center/edge distinction from above) which the CNN prediction is based on—the radical profile shapes at t = 0.8t ∗ may be informative even if the radical maximal values are not. 8. Ignition criteria as a priori predictors All three ignition prediction criteria were trained and evaluated on the 38 simulations introduced in Section 5.2 and Fig. 9. For the test set, we use the 12 simulations introduced in Section 5.3. The results are in Figs. 17–19, respectively. In Fig. 17, we can see that

464

P.P. Popov, D.A. Buchta and M.J. Anderson et al. / Combustion and Flame 206 (2019) 451–466

Fig. 16. Sensitivities of ICNN to the network inputs.

Fig. 17. Constructed early ignition criterion from Section 5. The horizontal lines indicate the ignition and non-ignition thresholds.

Fig. 18. The simple NN early ignition criterion of Section 6.

Fig. 19. The CNN predictive criterion of Section 7. Ignition is predicted for ICNN > 0.5 and non-ignition for ICNN < 0.5.

for the constructed criterion, 1 test set simulation is classified as ambiguous, and the remainder are correctly classified as either igniting or non-igniting. Though the criterion gives no false predictions, it is inconclusive for 4 of the training set simulations. In contrast, both of the machine-learning criteria (Figs. 18 and 19) correctly classify all of the training and test set cases. In this sense, the machine learning criteria perform better than the constructed criterion, and similarly to each other, though the CNN criterion has the advantage over the two-layer NN criterion in that it does not require a training-set-specific indicator function Ip or its regularization. It should be clear, however, that the machinelearning methods still leveraged the researchers’ physical insight. This is evident both in the choice of inputs—the minor radicals at 0.6t∗ , 0.8t∗ , t∗ —and from the input manipulation (logarithmic preconditioning and preferential weighting close to the stoichiometric surface), which was required for their success. Thus the machine learning methods are probably best viewed as a tool to facilitate the analysis of multi-dimensional data, not as a general classifier into which to feed data indiscriminately. These results indicate a route to generating an early ignition criteria for more complex chemical mechanisms. The simulations over which the present criteria were designed include two different free-stream flow velocities, resulting in either laminar or turbulent boundary layers. Additionally, LIB hotspots are at different locations, with significantly different local flow conditions: the local velocity strain rate varies between 500 s−1 and 12, 0 0 0 s−1 . Furthermore, simulations both with and without DBD, and simulations with two different strengths of the LIB hotspot are included in the training validation and test sets for all three criteria. Therefore, we expect the criteria to generalize well provided that the flow conditions (jet in crossflow, in similar flow regime) and chemistry are similar. There is an expectation that the samples on which we trained do not span the entire parameter space for ignition by a hotspot. Much larger or much smaller LIB hotspots, or different flow fields will likely require some re-training of the networks. The possibility that a pre-training on the current samples might accelerate learning in a novel configuration is an exciting research direction. We also anticipate multiple additional opportunities for improved training and generality. Due to the local character of LIB ignition, the present criteria only uses local information. It should be possible to obtain relevant yet computationally cheaper training samples by simulating ignition on an appropriate ersatz configuration, such as a cube-shaped domain of comparable size to the ignition region. Additionally, it may be possible to automate smallscale experiments so that they can span a large parameter range,

P.P. Popov, D.A. Buchta and M.J. Anderson et al. / Combustion and Flame 206 (2019) 451–466

465

thus providing experimental data on which to train ignition criteria. The prediction goals of the present network are specific to the needs of cost control is large-scale simulation, which does not have a direct analog in experimental since the sustained flame would be easily achieved. However, within a combustion system, there are other needs and opportunities. For examples, CNNs could be trained for blowoff prediction and combustion control. Though the implementation of such a control mechanism might pose significant practical challenges, advances in training such as we present and increased availability of computational hardware specialized to machine learning applications will make evaluation sufficiently quick for many applications. For the present configuration, we can make specific time-scale estimates for a blowoff prediction. The mean strain time scale is 8.5 × 10−4 s for the 5 m/s case and 2.8 × 10−4 s for the 15 m/s. Without any specific effort to accelerate it, the current CNN evaluation time is 1.3 ms (using a CPU/GPU combination), so as implemented the present network design is too slow for control. However, pruning of the trained CNN network can speed up the inference time significantly [57]. Additionally, faster chip architectures specialized for machine learning, such as Google’s Tensor Processing Units (TPU) [58] could provide extra speedup. The specific two-dimensional field loss functions we introduce would require online optical diagnostics, but other less invasive choices might be sufficient depending on the goals of the study.

radical profile, and filters which are focused on its edge. The former are well-suited for measuring radical maxima, while the latter are appropriate for measuring radical spread and compositions at the flame kernel’s boundary. It was found that the prediction outcome is positively sensitive to the edge-filters, and negatively sensitive to the interior-filters, indicating that the network assicoates ignition with features on the outer edge of the kernel, and extinction with features on its interior.

9. Summary

References

In this study, we examined the ignition of a hydrogen jet in a crossflow, with LIB as the ignition source and DBD as a control mechanism. The multi-physics simulation procedure is based on a compressible reactive flow solver with an integrated electric field solver for the DBD, which provides the basis of the DBD body force and radical generation. LIB was included from an offline simulation, taking advantage of the separation of scales A set of three criteria for early ignition prediction were developed and tested as a priori predictors. The constructed criterion performs well, with only five inconclusive and no incorrect classifications, though it required significant hand tuning. Performance was improved for two machine-learning approaches. A simple two layer NN model used Bayesian regularization-enhanced training and a logarithmic regularization of inputs. This method provided, in essence, a more complete and seemingly reliable ignition indication. However, it still relied on a designed ignition pointer, whose calculation for the training set samples requires the simulation of each case for a set long interval, even if ignition (or not) is already conclusively observed. A more involved convolutional neural network (CNN) approach obviated the need to define an ignition pointer, yielding the same accuracy when trained to predict a simple binary target (ignition/no ignition). It should be noted that the somewhat superior performance of the CNN over the NN criterion should in no way be taken as evidence that CNNs are inherently superior: in the present case, we included the NN to demonstrate that a simple architecture can also yield good performance, with a bit more manual work on the part of the researcher. Overall, we conclude that machine learning can be a very useful tool for ignition prediction, allowing the researcher to more effectively analyze the available high-dimensional data. At the same time, the researcher’s physical insight is essential for formulating the problem in such a way that the networks can converge to the data without overfitting. The behavior of the CNN’s intermediate layers suggests that the final evaluation is made based on multiple fundamentally different internal metrics. Specifically, the CNN’s layers progressively split the input into filters which are focused on the interior of the

Acknowledgments We thank Prof. Greg Elliott and Drs. Ryan Fontaine and Munetake Nishihara for providing and discussing experimental results. We also thank Profs. Carlos Pantano and Marco Panesi, and Drs. Alessandro Munafo and Pilbum Kim for discussions and work on the integration of LIB and DBD in the overall simulation procedure. We thank Mr. Houyi Du for discussions on machine learning training strategies. This material is based in part upon work supported by the Department of Energy, National Nuclear Security Administration, under Award Number DE-NA0 0 02374. Supplementary material Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.combustflame.2019.05. 014.

[1] A. Starikovskiy, N. Aleksandrov, Plasma-assisted ignition and combustion, Prog. Energy Combust. Sci. 39 (2013) 61–110. [2] M. Champion, B. Deshaies, G. Joulin, Relative influences of convective and diffusive transports during spherical flame initiation, Combust. Flame 474 (1988) 161–170. [3] A. Frendi, M. Sibulkin, Dependence of minimum ignition energy on ignition parameters, Combust. Sci. Technol. 73 (1990) 395–413. [4] P.S. Tromans, R.M. Furzeland, An analysis of lewis number and flow effects on the ignition of premixed gases, Proc. Combust. Inst. (1986) 1891–1897. [5] S.H. Lam, Using CSP to understand complex chemical kinetics, Combust. Sci. Technol. 89 (1993) 375–404. [6] Z. Luo, C.S. Yoo, E.S. Richardson, J.H. Chen, C.K. Law, T. Lu, Chemical explosive mode analysis for a turbulent lifted ethylene jet flame in highly-heated coflow, Combust. Flame 159 (2012) 265–274. [7] A. Akintayo, K.G. Lore, S. Sarkar, S. Sarkar, Prognostics of combustion instabilities from hi-speed flame video using a deep convolutional selective autoencoder, Int. J. Progn. Health Manag. [8] E. Schaffernicht, V. Stephan, K. Debes, H.M. Gross, Machine learning techniques for self-organizing combustion control, 32nd Annual Conference on Artificial Intelligence, Springer, Paderborn, Germany (2009), pp. 395–402. [9] M. Monfort, T. Luciani, J. Komperda, B. Ziebart, F. Mashayek, G.E. Marai, A deep learning approach to identifying shock locations in turbulent combustion tensor fields, Modeling, Analysis, and Visualization of Anisotropy, Springer (2017), pp. 375–392. [10] S.M. Starikovskaya, Plasma assisted ignition and combustion, J. Phys. D: Appl. Phys. 39 [11] H. Do, M.A. Cappeli, M.G. Mungal, Plasma assisted cavity flame ignition in supersonic flows, Combust. Flame 157 (2010) 1783–1794. [12] L. Massa, J.B. Freund, Plasma-combustion coupling in a dielectric-barrier discharge actuated fuel jet, Combust. Flame 184 (2017) 208–232. [13] A. Alberti, A. Munafo, A. Sahai, C. Pantano, M. Panesi, FEM simulation of laser-induced plasma breakdown experiments for combustion applications, 55th AIAA Aerospace Sciences Meeting, Grapevine, TX (2017). AIAA 2016–1111 [14] J.E. Retter, R.A. Fontaine, J.B. Freund, N.G. Glumac, G.S. Elliott, Coaxial DBD actuator design for control of a hydrogen diffusion flame, 54th AIAA Aerospace Sciences Meeting, San Diego, CA (2016). AIAA 2016–0199 [15] L. Massa, J. Capecelatro, D.J. Bodony, J.B. Freund, An integrated predictive simulation model for the plasma-assisted ignition of a fuel jet in a turbulent crossflow, 54th AIAA Aerospace Sciences Meeting, San Diego, CA (2016). AIAA 2016–2154 [16] J. Capecelatro, L. Massa, J.B. Freund, XPACC: simulation-1, The Center for Exascale Simulation of Plasma-Coupled Combustion (XPACC), Urbana, IL, 2015. Technical report [17] B. Strand, Summation by parts for the finite difference approximations for d/dx, J. Comput. Phys. 110 (1994) 47–67. [18] D.J. Bodony, Accuracy of the simultaneous-approximation-term boundary condition for time-dependent problems, J. Sci. Comput. 43 (2010) 118–133. [19] S.K. Lele, Compact finite-difference schemes with spectral-like resolution, J. Comput. Phys. 103 (1992) 16–42.

466

P.P. Popov, D.A. Buchta and M.J. Anderson et al. / Combustion and Flame 206 (2019) 451–466

[20] M.R. Visbal, D.V. Gaitonde, High-order-accurate methods for complex unsteady subsonic flows, AIAA J. 37 (10) (1999) 1231–1239. [21] D.V. Gaitonde, M.R. Visbal, High-order schemes for navier-stokes equations: algorithm and implementation into FDL3DI, 45433, Air Force Research Laboratory, Wright-Patterson Air Force Base, Ohio, 1998. Technical Report AFRL– VA-WP-TR-1998-3060, Air Vehicles Directorate [22] A.L. Sanchez, F.A. Williams, Recent advances in understanding of flammability characteristics of hydrogen, Prog. Energy Combust. Sci. 41 (2013) 1–55. [23] Chemical-kinetic mechanisms for combustion applications, http://web.eng. ucsd.edu/mae/groups/combustion/mechanism.html, accessed: 2018-04-30 [24] B.J. McBride, S. Gordon, M.A. Reno, Coefficients for calculating thermodynamic and transport properties of individual species, 4513, NASA, 1993. NASA Technical Memorandum. [25] R.J. Kee, M.E. Coltrin, P. Glarborg, Chemically reacting flow: theory & practice, John Wiley & Sons, 2003. [26] G. Dixon-Lewis, I. Transport, Flame structure and flame reaction kinetics i phenomena in multicomponent systems, Proc. R. Soc. A 307 (1968) 111–135. [27] J.O. Hirschfelder, C.F. Curtiss, R.B. Bird, Molecular theory of gases and liquids, Wiley, New York, 1954. [28] D.M. Orlov, Modelling and simulation of single dielectric barrier discharge plasma actuators, Ph.d. thesis. University of Notre Dame, 2006. [29] M.R. Stuart, Dielectric constant of quartz as a function of frequency and temperature, J. Appl. Phys. 26 (1955) 1399–1404. [30] N. Benard, A. Debien, E. Moreau, Time-dependent volume force produced by a non-thermal plasma actuator from experimental velocity field, J. Phys. D: Appl. Phys. 46. [31] C.J. Butler, A.N. Hayhurst, Measurements of the concentrations of free hydrogen atoms from observations of ions: Correlation of burning velocities with concentrations of free hydrogen atoms, Combust. Flame 115 (1998) 241–252. [32] C.D. Michelis, Laser induced gas breakdown: a bibliographical review, IEEE J. Quantum Electron. 5 (1969) 188–202. [33] P.P. Popov, M. Nishihara, D.A. Buchta, M.J. Anderson, K. Tang, L. Massa, J. Capecelatro, D.J. Bodony, G. Elliott, J.B. Freund, Computational prediction and experimental measurements of H and O radicals in a turbulent jet in crossflow, The Center for Exascale Simulation of Plasma-Coupled Combustion (XPACC), Urbana, IL, 2018. Technical report [34] T.X. Phuoc, Laser-induced spark ignition fundamental and applications, Opt. Lasers Eng. 44 (5) (2006) 351–397. [35] H. Torikai, Y. Soga, A. Ito, Schlieren visualization of blast extinguishment with laser-induced breakdown, Proc. Combust. Inst. 36 (2) (2017) 3297–3304. [36] G.D. Alamo, F.A. Williams, A.L. Sanchez, Hydrogen-oxygen induction times above crossover temperatures, Combust. Sci. Technol. 176 (2004) 1599–1626. [37] M.T. Hagan, H.B. Demuth, M.H. Beale, Neural network design, PWS Publishing, Boston, MA, 1996. [38] G. Montavon, G.B. Orr, K.R. Müller, Neural networks: tricks of the trade, Springer, 1998.

[39] D.J.C. MacKay, Bayesian interpolation, Neural Comput. 4 (1992) 415–447. [40] F.D. Foresee, M.T. Hagan, Gauss-newton approximation to Bayesian learning, International Joint Conference on Neural Networks (1997). [41] Neural network toolbox, https://www.mathworks.com/help/nnet/, accessed: 2018-08-29. [42] Z.H. Zhou, Ensemble methods: foundations and algorithms, CRC Press, Boca Raton, FL, 2012. [43] T. Wiatowski, H. Bölcskei, A mathematical theory of deep convolutional neural networks for feature extraction (2017). [44] C.M. Bishop, Pattern recognition and machine learning, Springer, 2009. [45] D. Nie, X. Cao, Y. Gao, L. Wang, D. Shen, Estimating CT image from MRI data using 3d fully convolutional networks, Deep Learning and Data Labeling for Medical Applications, Springer, Cham, Switzerland (2016). [46] Q. Dou, H. Chen, Y. Jin, L. Yu, J. Qin, P.A. Heng, 3d deeply supervised network for automatic liver segmentation from ct volumes, Medical Image Computing and Computer-Assisted Intervention, MICCAI, Springer, Cham, Switzerland (2016). [47] A. Karpathy, G. Toderici, S. Shetty, T. Leung, R. Sukthankar, L. Fei-Fei, Large-scale video classification with convolutional neural networks, IEEE Conference on Computer Vision and Pattern Recognition, Columbus, OH (2014). [48] Keras documentation, https://keras.io/, accessed: 2018-08-29. [49] ML practicum: image classification, https://developers.google.com/machinelearning/practica/image-classification/, accessed: 2018-08-29. [50] C.M. Bishop, Neural networks for pattern recognition, Oxford University Press, 1995. [51] G. Hinton, N. Srivastava, K. Swersky, Neural networks for machine learning: lecture 6a, overview of mini-batch gradient descent, https://www.cs.toronto. edu/tijmen/csc321/slides/lecture_slides_lec6.pdf. accessed: 2018-09-25. [52] N. Srivastava, G. Hinton, A. Krizhevsky, I. Sutskever, R. Salakhutdinov, Dropout: a simple way to prevent neural networks from overfitting, J. Mach. Learn. Res. 15 (2014) 1929–1958. [53] S. Hochreiter, J. Schmidhuber, Long short-term memory, Neural Comput. 9 (1997) 1735–1780. [54] J.J. Hopfield, Neural networks and physical systems with emergent collective computational abilities, Proc. Natl. Acad. Sci. 8 (1982) 2554–2558. [55] M. Sundermeyer, R. Schlüter, H. Ney, Lstm neural networks for language modeling, Thirteenth Annual Conference of the International Speech Communication Association, Portland, OR (2012). [56] O. Vinyals, A. Toshev, S. Bengio, D. Erhan, Show and tell: a neural image caption generator (2015). [57] S. Han, H. Mao, W.J. Dally, Deep compression: compressing deep neural networks with pruning, trained quantization and Huffman coding (2015). [58] MLPerf v0.5 results, https://mlperf.org/results/, accessed: 2019-04-18.