Acceleration of Premixed Flames in Obstructed Pipes with Both Extremes Open

: Premixed ﬂame propagation in obstructed channels with both extremes open is studied by means of computational simulations of the reacting ﬂow equations with a fully-compressible hydrodynamics, transport properties (heat conduction, di ﬀ usion and viscosity) and an Arrhenius chemical kinetics. The aim of this paper is to distinguish and scrutinize various regimes of ﬂame propagation in this conﬁguration depending on the geometrical and thermal-chemical parameters. The parametric study includes various channel widths, blockage ratios, and thermal expansion ratios. It is found that the interplay of these three critical parameters determines a regime of ﬂame propagation. Speciﬁcally, either a ﬂame propagates quasi-steady, without acceleration, or it experiences three consecutive distinctive phases (quasi-steady propagation, acceleration and saturation). This study is mainly focused on the ﬂame acceleration regime. The accelerating phase is exponential in nature, which correlates well with the theoretical prediction from the literature. The accelerating trend also qualitatively resembles that from semi-open channels, but acceleration is substantially weaker when both extremes are open. Likewise, the identiﬁed regime of quasi-steady propagation ﬁts the regime of ﬂame oscillations, found for the low Reynolds number ﬂames. In addition, the machine learning logistic regression algorithm is employed to characterize and di ﬀ erentiate the parametric domains of accelerating and non-accelerating ﬂames.


Introduction
Semi-open channels and tubes (where one end of the pipe is closed while another end is open; a flame is ignited at the closed end and propagates towards the open extreme), with evenly arranged and closely packed obstacles along the sidewalls, have been reported to experience ultrafast premixed flame acceleration (FA) and eventually an event of deflagration-to-detonation transition (DDT) [1][2][3][4][5]. This acceleration mechanism is conceptually laminar, being devoted to a powerful jet flow along the pipe centerline [1]. The cause of such acceleration is attributed to the cumulative effects of delayed burning in the pockets between the obstacles, with turbulence playing only a supplementary role [2]. On the contrary, industrial and laboratory conduits oftentimes have both extremes open, which attracted the interest for decades, starting with the classical experimental works on chocked flames [6] and quasi-detonations [7] in open obstructed pipes; see also [8,9] and references therein for details. In such a situation of both open extremes of an obstructed pipe, a new gas volume, generated by expansion of the burning matter in the combustion process, flows from the side pockets into two distributed streams moving towards the two open ends, as illustrated in Figure 1. Non-restriction of the flow in axial direction results in a mechanism of flame propagation different from that previously dominate in the flame dynamics [10,11]. Occurrence of these oscillations is due to regularly repeated acceleration-deceleration sequences causing the flame to have varying instantaneous shape, area and velocity [10]. A significant difference in flame propagation observed in these pipe configurations, coupled to potential improvements of practical applications [12] and associated with deeper understanding of flame propagation in such applications, have stirred up a significant interest. In a bid to satisfy the prevailing concerns and gain further insight into the morphology and propagation of premixed flames in obstructed conduits, a number of studies of burning in semi-open obstructed channels have been performed [13][14][15]. The analytical, numerical and experimental investigations of FA in fullyopen channels suggested that hydraulic resistance of the fluid in a pipe is responsible for sudden acceleration and the initial stage of DDT [16,17]. However, a subsequent analytical investigation [18] showed that the hydraulic resistance is not necessary to drive obstacles-based acceleration in an open channel; in fact, it hinders FA at the initial stage of the combustion process. In particular, flame propagation in narrow obstructed channels with open ends has been computationally studied [10], where regular fluctuations of the burning rates around the quasi-steady solutions have been found for the majority of the simulation runs performed. Other geometrical and thermochemical parameters such as the obstacles blockage ratio, the spacing between the obstacles and the thermal expansion ratio have been observed to influence the flame oscillatory behavior in terms of the oscillation periods and amplitudes. However, some pilot studies reported in [10] showed the possibility of FA after the initial oscillations or steady flame propagation.
In this light, exploring the flame behavior and characterizing various propagation regimes that the flame undergoes constitute the basis of the present work. Specifically, we aim to elucidate and quantify the morphology and dynamics of a premixed flame propagating through a comb-shaped array of obstacles in open channels of various blockage ratios, considering fuel mixtures of various thermal expansion ratios. The analysis is conducted by means of computational simulations of the reacting flow equations involving transport properties, a fully-compressible hydrodynamics and an Arrhenius chemical kinetics. A special attention is paid to relatively wide channels, in which the possibility of flame portraying more than one propagation regimes has been indicated. In the present work, it is also of our interest to validate the theoretical formulation [18] by comparing it with the simulation results obtained. Subsequently, we have identified three main stages of flame propagation in a fully-open obstructed pipe, namely: (i) quasi-steady propagation, (ii) exponential FA, and (iii) saturation of the burning velocity. Comparison of the exponential acceleration stage with a respective prediction from [18] shows a good correlation, especially at a later stage of exponential acceleration.

Research Methodology
The research methodology employed here is based on scrutinizing reacting flows in fully-open, obstructed channels by means of computational simulations of the governing partial differential equations (PDE) for a conservation/balance of mass, momentum, energy and species, including a A significant difference in flame propagation observed in these pipe configurations, coupled to potential improvements of practical applications [12] and associated with deeper understanding of flame propagation in such applications, have stirred up a significant interest. In a bid to satisfy the prevailing concerns and gain further insight into the morphology and propagation of premixed flames in obstructed conduits, a number of studies of burning in semi-open obstructed channels have been performed [13][14][15]. The analytical, numerical and experimental investigations of FA in fully-open channels suggested that hydraulic resistance of the fluid in a pipe is responsible for sudden acceleration and the initial stage of DDT [16,17]. However, a subsequent analytical investigation [18] showed that the hydraulic resistance is not necessary to drive obstacles-based acceleration in an open channel; in fact, it hinders FA at the initial stage of the combustion process. In particular, flame propagation in narrow obstructed channels with open ends has been computationally studied [10], where regular fluctuations of the burning rates around the quasi-steady solutions have been found for the majority of the simulation runs performed. Other geometrical and thermochemical parameters such as the obstacles blockage ratio, the spacing between the obstacles and the thermal expansion ratio have been observed to influence the flame oscillatory behavior in terms of the oscillation periods and amplitudes. However, some pilot studies reported in [10] showed the possibility of FA after the initial oscillations or steady flame propagation.
In this light, exploring the flame behavior and characterizing various propagation regimes that the flame undergoes constitute the basis of the present work. Specifically, we aim to elucidate and quantify the morphology and dynamics of a premixed flame propagating through a comb-shaped array of obstacles in open channels of various blockage ratios, considering fuel mixtures of various thermal expansion ratios. The analysis is conducted by means of computational simulations of the reacting flow equations involving transport properties, a fully-compressible hydrodynamics and an Arrhenius chemical kinetics. A special attention is paid to relatively wide channels, in which the possibility of flame portraying more than one propagation regimes has been indicated. In the present work, it is also of our interest to validate the theoretical formulation [18] by comparing it with the simulation results obtained. Subsequently, we have identified three main stages of flame propagation in a fully-open obstructed pipe, namely: (i) quasi-steady propagation, (ii) exponential FA, and (iii) saturation of the burning velocity. Comparison of the exponential acceleration stage with a respective prediction from [18] shows a good correlation, especially at a later stage of exponential acceleration.

Research Methodology
The research methodology employed here is based on scrutinizing reacting flows in fully-open, obstructed channels by means of computational simulations of the governing partial differential Energies 2020, 13, 4094 3 of 19 equations (PDE) for a conservation/balance of mass, momentum, energy and species, including a fully-compressible hydrodynamics, transport properties (heat conduction, diffusion and viscosity), and an Arrhenius chemical kinetics represented by a one-step irreversible reaction of the first order.

Governing Equations
In a two-dimensional (2D) geometry, considered in the present work, the governing PDEs read: where Y F is a progress variable representing the mass fraction of the fuel mixture; e = QY F + C v T and h = QY F + C p T are the specific internal energy and enthalpy, respectively, with the specific heats at constant volume, C v , and pressure, C P , and the energy release from the reaction Q = (Θ − 1)C P T f . Here T f is the initial temperature of the fuel mixture and Θ ≡ ρ f /ρ b is the thermal expansion ratio. Equation (3) describes a single, first-order Arrhenius reaction, with the activation energy E a and the constant of time dimension τ R . In our code, the cumulative choice (interplay) of the three parameters Q, E a and τ R constitutes an eigenvalue problem, solution to which yields the unstretched laminar burning velocity S L = S L (Q, E a , τ R ), which is also a unity of velocity dimension in the present work. Here we took S L = 3.47 m/s resembling typical hydrogen-air flames as well as oxy-methane or oxy-propane combustion. The speed of sound in this fuel mixture is 100 times larger, c s = 347 m/s such that the hydrodynamics is almost incompressible at the initial stage of burning, with the initial Mach number associated with flame propagation being as small as M 0 ≡ S L /c s = 10 −2 . However, for a corrugated flame front, the real (total) burning velocity U w exceeds S L , and it is calculated as [19]: Both the fuel mixture and burnt matter are assumed to be ideal gases obeying the equation of state: with the same molecular weight of m = 2.9 × 10 −2 kg/mol and the specific heats of C v = 5R p /2m and C p = 7R p /2m, where R u = 8.314 kJ/(kmol·K) is the universal gas constant. The stress tensor γ i,j and the energy diffusion vector q i are given by with the dynamic viscosity ζ = ρν, the Prandtl, Pr, and Schmidt, Sc, numbers. To avoid any impact of the diffusional-thermal instability, here we consider equidiffusive burning, so the Lewis number (defined as the thermal-to-mass diffusivities ratio) is unity, Le = 1. In the simulations, this condition is imposed by taking the Prandtl and Schmidt numbers Pr = Sc = 1 such that Le ≡ Sc/Pr = 1. Then the thermal flame thickness can be defined, conventionally, as: Energies 2020, 13, 4094

of 19
It is however noted that the quantity (7) is only a thermal-chemical parameter of length dimension, while the real size of the internal flame structure can be an order of magnitude larger [20]. Indeed, according to [20], a real width of the burning zone, determined by the characteristic length scale of the temperature profile gradient, is about ( 4 ∼ 5)L f , while it is the thickness of the active reaction zone, which is actually of the order of L f [20]. Nevertheless, L f is one of the key parameters in our studies, for the following reasons. First, it controls the choice of the computational mesh size making sure that the grid is sufficiently fine to resolve the flame front, having several grid point inside it; see a discussion below. Second , L f is a characteristic flame length scale. Along with a characteristic flame velocity scale S L , they determine the characteristic flame time scale t f ≡ L f /S L . Moreover, in addition to L f , another important length scale in this work is the radius (half-width) of the channel R, but it will also be measured in terms of L f , because their ratio constitutes the Reynolds number associated with flame propagation, Re = RS L /ν = R/PrL f = R/L f . The ratio of R and S L allows introducing a convenient dimensionless time τ ≡ tS L /R. Consequently, in our investigations, we are monitoring the instantaneous scaled burning velocity, U w (τ)/S L , as well as the instantaneous scaled flame tip positions and velocities, Z tip (τ)/R and U tip (τ)/S L , versus the scaled time τ, for the given values of Finally, the quantity of L f may be needed to impose appropriate ignition conditions. We usually imitate the initial flame structure by the classical Zeldovich-Frank-Kamenetskii (ZFK) solution for the profiles of the temperature, fuel mass fraction, pressure and velocities in the domain before and after a planar flame front, which reads [21]: such that Y F = 0, and T = T b = ΘT f in the burnt matter (z < 0), while Y F → 1 , and T → T f in the fresh fuel mixture (z L f > 0). The original pressure P throughout the domain equals the initial pressure of the fuel mixture, P f , and both velocity components are set to zero at the very beginning.

Computational Platform
We use an in-house computational fluid dynamics (CFD) fully-compressible, finite-volume Navier-Stokes code, which solves, in particular, the reacting flow Equations (1)- (6). The embryo of this code was originally developed in the middle of 1990s by Eriksson at Volvo Aero. In the following decades, this computational platform has been comprehensively updated and tested by numerous research groups worldwide, including (but not limited to) researchers from Chalmers University of Technology, Royal Institute of Technology, Uppsala and Umeå Universities in Sweden, Tsinghua University in China, Princeton and West Virginia Universities in the U.S. The current version of the solver is adapted for parallel computations and is available in 2D versions (Cartesian and cylindrical axisymmetric ones) as well as in a fully three-dimensional (3D) Cartesian version (though solely a 2D Cartesian set of Equations (1)-(6) is modelled in the present work). In order to save the computational time while resolve the regions of particular interest (leading pressure and shock waves, detonation or flame fronts), a self-adaptive structured computational grid is employed, as will be detailed below. This adaptive grid makes the computational package perfect, in particular, for fundamental studies of reacting and non-reacting flows in conduits of high aspect ratios such as used in the present work.
Finally, it should be mentioned in passing that modelling of DDT often faces difficulties. In particular, no real mechanism of detention triggering is available for a simplified Arrhenius chemical kinetics employed here, while detailed chemistry would affect the results, in particular, the timing of the detonation onset. However, in the present work, we do not actually focus on DDT and detonation but consider mostly FA at the early (subsonic) stage of the combustion process. Moreover, we analyze the scaled values U tip /S L , Z tip /R, τ, Re, Θ etc.-rather than the dimensional quantities of velocities, lengths, time and energies. This makes our consideration rather qualitative than quantitative. In other words, it is conventionally adopted that the chemistry is "immersed" in the chosen parameters such as E a , τ R , Q and, thereby, S L , while the evolution of the scaled variables should not change much. In fact, in this work we scrutinize FA in fully-open obstructed channels as compared with the situations considered in our previous works [1][2][3][4][5]10]. Since all these works employed the same computational platform, with the same one-step Arrhenius kinetics, there is no background to use anything else in the present work, which justifies the usage of this computational platform for this particular problem.

Computational Approach, Parameters and Domain
The numerical approach is based on a cell-centered, finite-volume numerical scheme, which is of the second-order accuracy in time and the fourth-order in space for the convective terms, and of the second-order in space for the diffusive terms. The code employs the balance Equations (1)-(3) in a unified form [25] ∂G ∂t where G stands for any of the variables ρ, ρu x , ρu z , ρY F and e, while E G and F G designate the related axial and radial fluxes, respectively, and H G is the source term. The spatial discretization is obtained by integrating any of the conservation laws (1)-(3) in the form (11) over a given grid cell. More details about the numerical method are available, for instance, in [2,5,10,25]. In the present work, a premixed flame propagates in a long 2D channel of width 2R, with both extremes open, where the fraction of width 2Rα, near the walls, is blocked by the obstacles with the spacing (the distance between two neighboring obstacles) ∆Z. In this light, the major dimensionless parameters of this study are the thermal expansion ratio Θ, the blockage ratio α, the scaled spacing ∆Z/R, and the Reynolds number associated with flame propagation Re = R/L f . To be specific, the following thermal-chemical and configurational parameters have been taken: various Θ = 5, 6, 8, 10; α = 1/3, 1/2, 2/3; R/L f = 24, 36, 48; and fixed ∆Z/R = 1/4. The initial temperature, pressure and density of the fuel mixture are T f = 300 K, P f = 100 kPa, and ρ f = 1.16 kg/m 3 , respectively.
In order to save the computational cost, similar to [5,10], we are looking for a symmetric solution such that only half of the channel is actually modeled by employing the symmetry condition at the bottom as illustrated in Figure 1. To keep the flame segments as planar as possible, we employ the walls of the channels and the obstacles to be adiabatic, n·∇T = 0, and free-slip, n·u = 0, where n is a Energies 2020, 13, 4094 6 of 19 normal vector at the surface of the wall/obstacle. The absorbing (non-reflecting) characteristic boundary conditions are set at both open extremes to prevent the reflection of the sound waves and weak shocks from these outlets. Such non-reflecting boundary conditions have been widely tested in our previous works [19,20,25]. Initially, the parameters at the right (unburnt) end are ρ = ρ f , T = T f and u z = −S L , while those at the left (burnt) extreme are ρ = ρ f /Θ, T = ΘT f , and u z = −ΘS L . The initial flame structure was imitated by the ZFK solution for a planar flame front, Equations (8)- (10). Such a planar ZFK flame was initiated at a distance 50 L f from the left (burnt) end of the channel.
The computational domain is being reconstructed depending on an instantaneous location of a flame front and the leading pressure wave, so the length of the domain grows with time. Such a reconstruction of the grid is achieved by using the third-order splines to interpolate the flow variables during the grid reconstruction and preserves the second-order accuracy of the numerical scheme. The dynamic grid consists of a structured rectangular grid, with the walls of the grid parallel to each direction. This adaptive mesh has the minimal cell size of 0.2 L f , which has been justifed in the previous works employing this solver; see, for instance, [10,20] (and references therein). Indeed, it is recalled that L f correlates with the width of the active reaction zone, while the thermal (temperature-based) flame thickness is about ( 4 ∼ 5)L f [20]. Consequently, taking the grid size of 0.2 L f , we obtain about five grid points inside the active reaction zone and 20-25 grid points inside the flame front. Therefore, the grid used in the present work has a sufficiently fine resolution as the thickness of the reaction zone appears distributed over several computational cells.

Results and Discussion
The morphology and propagation of premixed flames in obstructed channels with both ends open demonstrate unique and distinct characteristics as compared to semi-open channels. Initially, an oscillating trend has been observed for R/L f = 12 [10], with the oscillation characteristics being dependent on the thermal expansion ratio Θ, the blockage ratio α, and the obstacle spacing ∆Z. The series of snapshots in Figure 2 shows the flame development and the flow distributions for burning with Θ = 8 in a channel having R/L f = 48 and α = 1/2. An initially planar flame front undergoes folding as it propagates due to the influence of thermal expansion across the flame front, while the entire flow is directed towards the left extreme as seen in Figure 2a. Eventually, after propagating through circa two channel widths, Figure 2b, the flow is distributed between the two extremes and the flame has gradually attained a tulip-like shape. The momentum released from the left extreme, along with the accompanied heat losses, moderate flame propagation significantly. As seen in Figure 2, the pockets of the unburnt fuel mixture are generated as the flame moves further, thereby creating the effect of delayed burning (thermal expansion) from each of these pockets. As the unburnt gas in these pockets starts burning, the expanding gases flow into the core of the channel and are distributed between the two flows moving towards two extremes. With further propagation away from the left extreme, the unburned pockets are further away from the left outlet, Figure 2c,d, and the new gas volume, generated by thermal expansion in the burning process, flows in the direction causing least resistance, i.e., in the direction of the flame. In this regime, the flame starts undergoing significant acceleration, which is sustained as the flame front travels further. This is accompanied by continuous formation and isolation of the new pockets from the left extreme, thereby propelling the flame further through formation of a jet flow in the core path of the channel. This regime of burning is observed in Figure 2e-g. As burning progresses further, the flame velocity grows, and the flame shape develops into a highly corrugated structure as seen in Figure 2g. With increasing Mach number, the influence of gas compressibility becomes significant, subsequently causing moderation of FA and saturation of the burning velocity to a near-constant value. Interestingly, at the saturated velocity state, the flame front appears more stable such that the small-scale wrinkles disappear as seen in Figure 2h. Energies 2020, 13, x FOR PEER REVIEW 7 of 18 temperature snapshots as well as those for vorticity and shock waves are shown. It is seen that the leading vortices, inherent during the accelerating phases, are absent at the saturated phase, with strong vortices confined behind the flame front. Likewise, the leading shock waves in the accelerating phase are not apparent when the flame velocity saturates to a steady one.  These observations provide the insights into the physics of premixed flames in obstructed channels with both ends open. Namely, a flame undergoes initial quasi-steady propagation, which dominates in the region of the simultaneous (near equal) flows of the expanding gases towards both extremes of the channels. Subsequently, the onset of FA occurs, which is then supplanted by the state of saturated propagation velocity. Figure 3a shows the evolution of the flame tip position scaled by the channel half-width R. To visualize these regimes clearly, it is important to analyze the derivative of the flame tip position (i.e., the flame tip velocity). This quantity, scaled by the laminar flame speed S L , is presented in Figure 3b, indicating the regimes of quasi-steady propagation, exponential acceleration and the saturated flame velocity. The flame tip variations with time are qualitatively similar to the experimental results of Middha and Hansen [36] as well as that of Yanez et al. [37]. The exponential acceleration regime mimics the trend seen in the semi-open channel with closely packed obstacles [37], namely, exponential FA without initial quasi-steady propagation. However, such initial quasi-steady propagation is quite similar to the oscillatory regime identified in narrower channels [10]. Apparently, the oscillations can correlate with prolonged quasi-steady propagation as the flame tries to balance the complex flow dynamics and gas expansion in such a highly confined space.  Figure 4) to that at the saturated phase ( Figure 5). To be specific, the characteristic instantaneous temperature snapshots as well as those for vorticity and shock waves are shown. It is seen that the leading vortices, inherent during the accelerating phases, are absent at the saturated phase, with strong vortices confined behind the flame front. Likewise, the leading shock waves in the accelerating phase are not apparent when the flame velocity saturates to a steady one.     propagation trend is strongly influenced by the blockage ratios. Namely, as grows, the duration of quasi-steady propagation also increases. For the lowest blockage ratio considered, = 1/3, a flame propagates quasi-steadily without acceleration. For low α , the core (unobstructed) flow area is increased as such the volume of the gas flowing out increases accordingly, while the new gas volume, generated in the pockets, is reduced. Consequently, the contribution of delayed burning to propelling the flame front decreases, thereby moderating the flame velocity tremendously. We next study the impacts of the blockage ratio α, the channel half-width R and the thermal expansion ratio Θ on the flame dynamics. Specifically, Figures 6-8 show the scaled flame tip position Z tip /R and the scaled burning velocity U w /S L versus the scaled time τ = tS L /R. The flame propagation trend is strongly influenced by the blockage ratios. Namely, as α grows, the duration of quasi-steady propagation also increases. For the lowest blockage ratio considered, α = 1/3, a flame propagates quasi-steadily without acceleration. For low α, the core (unobstructed) flow area is increased as such the volume of the gas flowing out increases accordingly, while the new gas volume, generated in the pockets, is reduced. Consequently, the contribution of delayed burning to propelling the flame front decreases, thereby moderating the flame velocity tremendously.
We next study the impacts of the blockage ratio , the channel half-width and the thermal expansion ratio Θ on the flame dynamics. Specifically, Figures 6-8 show the scaled flame tip position / and the scaled burning velocity / versus the scaled time = / . The flame propagation trend is strongly influenced by the blockage ratios. Namely, as grows, the duration of quasi-steady propagation also increases. For the lowest blockage ratio considered, = 1/3, a flame propagates quasi-steadily without acceleration. For low α , the core (unobstructed) flow area is increased as such the volume of the gas flowing out increases accordingly, while the new gas volume, generated in the pockets, is reduced. Consequently, the contribution of delayed burning to propelling the flame front decreases, thereby moderating the flame velocity tremendously.   Moreover, increasing the channel half-width reduces the duration of quasi-steady propagation for = 1/2 and 2/3. This effect is observed more apparently in Figure 9a. Here, as the channel halfwidth grows from 24 to 48 , the scaled time , taken before acceleration starts, decreases by a factor of 10, from 6 to 0.6. Consequently, this initial time delay (prior to acceleration) diminishes as  Moreover, increasing the channel half-width reduces the duration of quasi-steady propagation for = 1/2 and 2/3. This effect is observed more apparently in Figure 9a. Here, as the channel halfwidth grows from 24 to 48 , the scaled time , taken before acceleration starts, decreases by a factor of 10, from 6 to 0.6. Consequently, this initial time delay (prior to acceleration) diminishes as the channel widens. Therefore, a growth of the channel width provides an impact similar to that of growing α. For a wider channel, = 48 , the large blockage ratios, = 1/2 and 2/3, provide the  Moreover, increasing the channel half-width reduces the duration of quasi-steady propagation for α = 1/2 and 2/3. This effect is observed more apparently in Figure 9a. Here, as the channel half-width grows from 24 L f to 48 L f , the scaled time τ, taken before acceleration starts, decreases by a factor of 10, from 6 to 0.6. Consequently, this initial time delay (prior to acceleration) diminishes as the channel widens. Therefore, a growth of the channel width provides an impact similar to that of growing α. For a wider channel, R = 48 L f , the large blockage ratios, α = 1/2 and 2/3, provide the similar effects on flame propagation and the difference in the durations of the stages of quasi-steady propagation is minor, see Figure 8, while no acceleration is seen for α = 1/3 in this case. Moreover, increasing the channel half-width reduces the duration of quasi-steady propagation for = 1/2 and 2/3. This effect is observed more apparently in Figure 9a. Here, as the channel halfwidth grows from 24 to 48 , the scaled time , taken before acceleration starts, decreases by a factor of 10, from 6 to 0.6. Consequently, this initial time delay (prior to acceleration) diminishes as the channel widens. Therefore, a growth of the channel width provides an impact similar to that of growing α. For a wider channel, = 48 , the large blockage ratios, = 1/2 and 2/3, provide the similar effects on flame propagation and the difference in the durations of the stages of quasi-steady propagation is minor, see Figure 8, while no acceleration is seen for = 1/3 in this case.  Similar to the effect of the blockage ratio, an increase in the thermal expansion ratio Θ shortens the stage of quasi-steady propagation, and at the subsequent stage of acceleration, the acceleration rate is increased as shown in Figure 9b. For the fuel mixtures of large Θ, the expansion produced across the flame is obviously enhanced, and the flame surface area will grow accordingly, leading to faster transition to the accelerating phase and faster FA. In contrast, the flames of lower thermal expansion propagates quasi-steadily without undergoing any significant acceleration.
In order to scrutinize the exponential regime of FA, which follows the quasi-steady phase observed in this scenario, we next recall a theory of FA in fully-open obstructed channels [18]. While [18] considered both inviscid and viscous formulations, here we chose the inviscid theory only in order to avoid considering the length on the conduit. Then, within the inviscid formulation, the flame tip velocity U tip obeys the evolution equation [18]: with the solution: where the scaled exponential acceleration rate for the fully open channel, σ o , is given by [18]: The indexes o and s in Equation (14)  and to understand the exponential acceleration regimes clearly, we plot Equation (13) and compare it with the simulation results. For this purpose, a quasi-steady stage of flame propagation is extracted from the simulation by shifting the flame tip position Z tip , accordingly, such that the instant t = 0 indicates that exponential FA starts. This comparison is shown in Figures 10-12 for the accelerating flame regime only. Figure 10 presents the trend in the R = 24 L f channel; in Figure 10a, the theory predicts the later stage of exponential acceleration quite well, while in Figure 10b, α = 2/3, the theory diverges from the simulation later in time but shows better prediction at an early propagation stage.
The indexes o and s in Equation (14) are devoted to the fully-open and semi-open channels, respectively. It is seen form Equation (14) that the exponential acceleration rate in fully-open channels is much less than that in semi-open channels, for the same fuel. Obviously, this is because the flow is split in two parts when both ends of a pipe are open. To validate this theory by our simulation results and to understand the exponential acceleration regimes clearly, we plot Equation (13) and compare it with the simulation results. For this purpose, a quasi-steady stage of flame propagation is extracted from the simulation by shifting the flame tip position , accordingly, such that the instant = 0 indicates that exponential FA starts. This comparison is shown in Figures 10-12 for the accelerating flame regime only. Figure 10 presents the trend in the = 24 channel; in Figure 10a, the theory predicts the later stage of exponential acceleration quite well, while in Figure 10b, = 2/3, the theory diverges from the simulation later in time but shows better prediction at an early propagation stage.   In Figure 11a, representing the channel of half-width 36 and = 1/2, the theory captures the later trend better for Θ = 10. However, for the = 2/3 case, Figure 11b, the later acceleration trend is better captured when Θ = 8. The trends for = 1/2, 2/3 and Θ = 8, 10 in the channel of half-  In Figure 11a, representing the channel of half-width 36 and = 1/2, the theory captures the later trend better for Θ = 10. However, for the = 2/3 case, Figure 11b, the later acceleration trend is better captured when Θ = 8. The trends for = 1/2, 2/3 and Θ = 8, 10 in the channel of halfwidth 48 are shown in Figure 12. Here, the later trend of exponential acceleration is better captured for = 2/3 and Θ = 8, 10, see Figure 12b. Although the theory yields a good approximation of  In Figure 11a, representing the channel of half-width 36 L f and α = 1/2, the theory captures the later trend better for Θ = 10. However, for the α = 2/3 case, Figure 11b, the later acceleration trend is better captured when Θ = 8. The trends for α = 1/2, 2/3 and Θ = 8, 10 in the channel of half-width 48 L f are shown in Figure 12. Here, the later trend of exponential acceleration is better captured for α = 2/3 and Θ = 8, 10, see Figure 12b. Although the theory yields a good approximation of exponential acceleration, it narrowly captures the trend in the earlier acceleration phase. In fact, theoretically, the exponential acceleration rates have definite value even for the low blockage and thermal expansion ratios. For instance, for α = 1/3 with Θ = 6, we have σ 0 = 2.17. However, no acceleration has actually been observed with such a set of parameters. The simulation results for the combination of such a low α and Θ do not show FA. In fact, acceleration is only noticeable when both α and Θ are large. Specifically, at α = 1/2 with Θ = 8, FA becomes noticeable when σ 0 = 3.66. Therefore, initiation of actual exponential acceleration will happen when 2.17 ≤ σ 0 ≤ 3.66, which corresponds to 1/3 ≤ α ≤ 1/2 and 6 ≤ Θ ≤ 8.
The resulting saturated velocity occurring due to gas compression is expected to be lower than the Chapman-Jouget (CJ) deflagration velocity, which is defined as [38]: Figure 13 presents the scaled flame tip velocity versus the scaled flame tip position (such a phase diagram), for Θ = 8 and Θ = 10 in Figure 13a,b, respectively. The solid horizontal lines show the scaled CJ deflagration velocity, Equation (15), which is thereby compared to the saturated flame tip velocity. It is seen that Equation (15) exceeds the computational saturated flame velocity for all the cases considered in Figure 13. In the case of Θ = 8, Figure 13a, weak R-dependence is observed: the saturated flame velocity is slightly higher for R = 24 L f than that for R = 36 L f . However, in the case of Θ = 10, Figure 13b, the saturated velocity curves collapse to the same value for all widths for each blockage ratios α = 1/2 and α = 2/3. We can therefore conclude from Figure 13 that the saturated flame tip velocity for α = 1/2 exceeds that for α = 2/3, and this value is almost independent of the channel width, particularly for Θ = 10.
Energies 2020, 13, x FOR PEER REVIEW 13 of 18 cases considered in Figure 13. In the case of Θ = 8, Figure 13a, weak -dependence is observed: the saturated flame velocity is slightly higher for = 24 than that for = 36 . However, in the case of Θ = 10, Figure 13b, the saturated velocity curves collapse to the same value for all widths for each blockage ratios = 1/2 and α = 2/3. We can therefore conclude from Figure 13 that the saturated flame tip velocity for = 1/2 exceeds that for = 2/3, and this value is almost independent of the channel width, particularly for Θ = 10. While the present simulations are limited to a 2D Cartesian geometry, it is of interest to discuss how the computational results might change in a 3D geometry, either cylindrical axisymmetric or 3D Cartesian, as both situations are closer to the practical reality than the 2D planar channel studied here. We start the discussion with the analogies observed in unobstructed conduits with both ends open, where computational simulations demonstrated regular flame oscillations in the 2D channels [20] but flame acceleration in the cylindrical pipes [32]. We expect the same qualitative trend in the obstructed conduits because a 3D geometry is generally associated with higher increment in the flame surface area and, thereby, faster flame propagation. The exponential acceleration rate grows more than twice.
Quantitatively, the cylindrical counterpart of the evolution Equation (12) will take the form: While the present simulations are limited to a 2D Cartesian geometry, it is of interest to discuss how the computational results might change in a 3D geometry, either cylindrical axisymmetric or 3D Cartesian, as both situations are closer to the practical reality than the 2D planar channel studied here. We start the discussion with the analogies observed in unobstructed conduits with both ends open, where computational simulations demonstrated regular flame oscillations in the 2D channels [20] but flame acceleration in the cylindrical pipes [32]. We expect the same qualitative trend in the obstructed conduits because a 3D geometry is generally associated with higher increment in the flame surface area and, thereby, faster flame propagation. The exponential acceleration rate grows more than twice.
Quantitatively, the cylindrical counterpart of the evolution Equation (12) will take the form: with the solution: (see [2,18] to understand the origin of these equations). According to Equation (17), the exponential acceleration rate in the cylindrical geometry exceeds its 2D counterpart more than twice indeed. The same trend has been seen in semi-open conduits, obstructed as well as unobstructed, when switching from the 2D planar to the cylindrical-axisymmetric geometries. In addition, turbulence is an essentially 3D phenomenon, so it will provide an extra effect in the 3D geometry. In all these respects, we expect much stronger flame acceleration in the 3D case, with a little (if any) chance of non-accelerative trends of flame propagation.
We next develop a predictive model for the propagation trends observed in obstructed channels with both ends open by employing a machine learning (ML)-based technology. The goodness of fit of the model has been calculated using the confusion matrix method for logistic regression method [39]. It is presented as the prediction accuracy for the models and it represents how well the model predicts the observed accelerating/non-accelerating flame trends. Specifically, a logistic regression algorithm is adopted to train these data sets in order to classify the accelerating versus non-accelerating flame trends considering Θ, α and Re = R/L f as the parameters. Logistic regression has been chosen to be employed because this is one of the most widely used ML algorithms for the classification problems, which nevertheless is relatively simple to implement. To justify such a choice, it is noted that logistic regression has already been employed in a similar complex problem of flame propagation, which involved fast flame detection to address accidental fire scenarios [40]. Figure 14 is a schematic of the ML model. The procedure involves feeding the learning algorithm with the training sets: Θ, α and/or Re and the "classification" variable. Here, the classification variable indicates accelerating or not accelerating event. So, the algorithm uses this data set to learn a hypothesis function, h, that maps input x to the predicted output y.
Energies 2020, 13, x FOR PEER REVIEW 14 of 18 employed because this is one of the most widely used ML algorithms for the classification problems, which nevertheless is relatively simple to implement. To justify such a choice, it is noted that logistic regression has already been employed in a similar complex problem of flame propagation, which involved fast flame detection to address accidental fire scenarios [40]. Figure 14 is a schematic of the ML model. The procedure involves feeding the learning algorithm with the training sets: Θ , and/or and the "classification" variable. Here, the classification variable indicates accelerating or not accelerating event. So, the algorithm uses this data set to learn a hypothesis function, ℎ, that maps input to the predicted output . Specifically, we consider a hypothesis function in logistic regression in the form: where is the sigmoid function, which is depicted in Figure 15 and is defined as: Specifically, we consider a hypothesis function in logistic regression in the form: where g is the sigmoid function, which is depicted in Figure 15 and is defined as: Specifically, we consider a hypothesis function in logistic regression in the form: where is the sigmoid function, which is depicted in Figure 15 and is defined as: Figure 15. The sigmoid function g(z) of Equation (19).
As a result, the function ℎ ( ) describes the probability of = 1 for a given value. As such, if ℎ ( ) ≥ 0.5 then = 1, and = 0 for ℎ ( ) < 0.5. From the simulation results, the cases undergoing acceleration are classified as 1, and the cases with only quasi-steady propagation are 0. The training sets consist of the input variable (Θ, α, ) and . The results obtained for the hypothesis function have been subsequently used to classify and develop the regime diagrams depicted in Figures 16-17. As a result, the function h θ (x) describes the probability of y = 1 for a given x value. As such, if h θ (x) ≥ 0.5 then y = 1, and y = 0 for h θ (x) < 0.5. From the simulation results, the cases undergoing acceleration are classified as 1, and the cases with only quasi-steady propagation are 0. The training sets consist of the input variable x(Θ, α, Re) and y. The results obtained for the hypothesis function have been subsequently used to classify and develop the regime diagrams depicted in Figures 16 and 17.     Figure 17 shows this result, for Θ = 8 and Θ = 10 in Figure 17a  Specifically, Figure 16a shows the classification results obtained for Re = 24, with the solid line indicating the decision boundary represented by the variable z. Essentially, the equation satisfying the accelerating phase is −131.69 + 9.36Θ + 112.29α ≥ 0, with a training accuracy of 91.67 %. For Re = 36 and 48 in Figure 16b,c, the predicted hypothesis is the same, and the decision boundary is represented by −111.73 + 10.19Θ + 69.85α ≥ 0, with the training accuracy of 100 %.
Finally, the data are classified based on the flame propagation Reynolds number Re, and the developed model predicts if a flame accelerates or not based on Re. Figure 17 shows this result, for Θ = 8 and Θ = 10 in Figure 17a,b, respectively. The model prediction is therefore −17.01 + 0.19Θ + 22.62α ≥ 0, with 83.3 % prediction accuracy for Θ = 8, and −14.4 + 1.14Θ + 39.9α ≥ 0, with a 100 % prediction accuracy for Θ = 10.
It should be noted that the prediction accuracies reported here are solely based on the training data set. However, in the future works, our research plan is to generate more data to improve and test the present model more extensively.

Conclusions
While apparent flame oscillations are found to dominate in the relatively narrow fully-open obstructed channels with Re = R/L f = 12, for the wider channels such as having Re = 24, 36, 48, the propagation trends are observed to be quite distinctive. To be specific, when the parameters Re, α and Θ are beyond their respective critical values, we have identified the three distinctive phases/ regimes of burning such as (i) initial quasi-steady flame propagation, followed by (ii) exponential FA and (iii) the final saturation of the flame velocity. Typically, all three regimes occur when Θ = 8, 10, α = 1/2, 2/3 and Re = 24, 36, 48. It is noticed that all three parameters (Re, α and Θ) significantly impact the duration of quasi-steady propagation. Specifically, while in the Re = 24 channel, quasi-steady propagation lasts as long as ∼ 6 R/S L , this time reduces to 0.6 R/S L when Re = 48. The regime of quasi-steady propagation is also influenced by the blockage ratio α and the thermal expansion ratio Θ. Namely, as α decreases, the volume of the pockets of unburned gas, left behind the leading flame segment, also reduces, causing drastic moderation of the effect of delayed burning, which is the major factor for the Bychkov mechanism of ultrafast FA [1]. Likewise, thermal expansion plays a similar role: a lower gas expansion, Θ, means reduced propulsion forces and, as a result, moderated propagation of the flame front.
The nature of exponential acceleration has been explored by comparing the simulation results with the theory [18] developed for FA in obstructed channels with closely packed obstacles and both extremes open. According to the plots of this comparison, the theory [18] predicts the later stage of exponential acceleration better than it predicts the earlier stage. The simulation results helped to understand the applicability domain of this theory to realistic flames. This domain is identified to be somewhere between {α = 1/3; Θ = 6} and {α = 1/2; Θ = 8}. Furthermore, the regime of the saturated burning velocity has been scrutinized in the present instigation. Specifically, the growing flame velocity saturates when in starts approaching the speed of sound. It has been identified that the main parameter influencing the final saturated flame tip velocity is the blockage ratio α. In particular, with α = 1/2, the flame velocity saturates to a higher value than that in the case of α = 2/3. It is also noted that all the saturated flame velocities obtained in the present simulations appeared lower than the respective CJ deflagration velocity.
Finally, a logistic regression ML algorithm was utilized to learn about the hypothesis functions that predict when a flame undergoes acceleration and when it continuously propagates in a quasi-steady state without FA. As a result, the models, identifying when the flames would experience all three regimes (quasi-steady propagation, acceleration, and saturation) and when it only propagates with near-constant velocity, are developed. Such an identification is based on the interplay between the Reynolds number associated with flame propagation Re, the blockage ratio α, and the thermal expansion ratio Θ, as demonstrated by the diagrams of Figures 16 and 17.