Dynamics of flames in tubes with no-slip walls and the mechanism of tulip flame formation1

ABSTRACT

A hydrogen/air flame propagation and the development of tulip-shaped flame in 2D tubes of different aspect ratios with both closed ends and in a half-open rectangular channel were studied using high resolution direct numerical simulations of the fully compressible Navier – Stokes equations coupled with a detailed chemistry. Flame propagation in a 3D rectangular channel was studied using large eddy simulations and compared with the results of direct numerical simulations of flame propagation in a 2D rectangular channel with the same aspect ratio. It is shown that the interaction of the rarefaction wave generated by the flame at the deceleration stage with the “positive” flow of unburned gas generated by the flame at the previous accelerating stage leads to a significant decrease of the velocity of the unburned gas flow in the near field zone ahead of the flame front. As a result, the thickness of the boundary layer in the near-field zone ahead of the flame increases significantly, and the profile of the axial velocity of the unburned gas in the near-field zone ahead of the flame front takes the form of a tulip or an inverted tulip, which leads to corresponding changes in the velocities of different parts of the flame front, the flame front inversion, and the formation of a tulip-shaped flame.

KEYWORDS:

• Tulip flame
• flame dynamics
• rarefaction waves
• pressure waves
• boundary layer
• vortices

Introduction

The dynamics of a flame propagating in closed or half-open tubes, is important for understanding combustion processes under confinement, such as explosions and safety issues as well as for industrial and technological applications, e.g. combustion in gas turbines and internal combustion engines. The inversion of the flame front, which propagates from the closed end down the tube, from a convex shape with the tip toward the unburned gas to a concave shape with the tip toward the burned gas, is known as tulip-shaped flame formation and has been observed in many experiments and numerical simulations. It should be noted that the inversion of the flame front can be caused by various processes, such as the interaction of the flame with the shock wave, various hydrodynamic instabilities inherent to the flame, such as Darrieus-Landau (DL) and/or Rayleigh-Taylor (RT) instabilities. This makes the concept of tulip-shaped flame formation not entirely definite and has led to a number of different scenarios in attempts to explain the mechanism of tulip-shaped flame formation since it was first discovered and photographed by Ellis and Robinson (1925), Ellis and Wheeler (1925), Ellis (1928) almost a century ago. A reader can find a long history of experimental investigation and attempts to explain the mechanism of the tulip flame formation in the detailed review by Dunn-Rankin (2009).

In recent years, interest in the flame dynamics in closed and half-open tubes has increased dramatically due to the need for a better understanding of combustion in engines, as well as flame dynamics during the transition from deflagration to detonation, which is required to solve safety problems. Numerous experimental and numerical papers have been published that have studied effects of the aspect ratio, the equivalence ratio of mixtures, the dependence on the ignition location, and the effect of reaction order. However, the physical mechanism leading to the formation of the tulip-shaped flame has only been discussed in a few papers, and despite previous and more recent efforts, there is still no decisive and convincing explanation for the formation of the tulip-shaped flame.

Hariharan and Wichman (2014 and Hariharan and Wichman 2015) explain the formation of the tulip flame as “The stagnation point eventually moves through the flame front and travels at a fixed distance from the tulip cusp at the centerline of the channel as the flame moves toward the far (right) end of the channel. This increases the velocity of the unburned mixture traveling toward the centerline of the flame … ” However, this is just a description of the observed evolution of the velocity stagnation point on the tube axis, which the authors call the saddle point, which is rather a consequence of tulip flame formation but not a physical mechanism explaining the tulip flame formation.

The experimental and numerical simulations were performed by Xiao, Houim, and Oran (2015) to study distorted tulip flames (DTF), which are observed for flames in highly reactive mixtures, and the formation of tulip flames preceding DTF. The conclusion about the mechanism of the tulip-flame formation according Xiao, Houim, and Oran (2015) was “ … our simulations support the mechanism of vortex motion described by Matalon and Metzener. Reverse flow forms and dominates the near-field region due to the vortex motion. This process creates conditions that allow the flame to invert. In addition, … an adverse pressure gradient is imposed on the flame front. The pressure gradient reverses the flow in the unburned region and creates vortices ahead of the flame front near the sidewalls by interacting with the boundary layers.” These are unsubstantiated statements, which are nothing more than a trivial phenomenological description of the figures obtained in simulations. The conclusions “support the mechanism of vortex motion” or “Reverse flame shapes… due to vortex motion” are not supported by any analysis or detailed study of the interaction of vortices with the flame or flow. It is worth noting that Xiao, Houim, and Oran (2015) used the value of hydrogen/air laminar flame velocity $U_f=3.11\mathrm{m}/\mathrm{s}$ and $\gamma =1.17$ for the adiabatic constant. The real values are $U_f=2.36\mathrm{m}/\mathrm{s}$ and $\gamma =1.4$.

In this paper, we use high-resolution 2D numerical simulations to study premixed flames in the stoichiometric hydrogen/air mixture propagating in closed rectangular 2D channels with various aspect ratios and in the half-open rectangular 2D channel, as well as Large Eddy Simulations (LES) to study a flame propagating in closed rectangular 3D channel. The 2D simulations used a high-order numerical method to solve the fully compressible reactive Navier – Stokes equations coupled with a detailed chemical model for H₂/air. The obtained results explain the mechanism of tulip flame formation and show that the tulip flame formation is a purely hydrodynamic phenomenon in agreement with the experimental study by Ponizy, Claverie, and Veyssière (2014). It is shown that the formation of the tulip flame is associated with the deceleration stage of the flame when the decelerating flame creates a rarefaction wave in the unburned gas. This leads to a significant decrease in the unburned gas flow velocity ahead of the flame front, as well as to a significant increase in the thickness of the boundary layer in the proximity region ahead of the flame front. As a result, the axial velocity profile in the near-field region ahead of the flame takes the form of a tulip or inverted tulip, leading to different axial velocities of different parts of the flame front, which in turn leads to in the inversion of the flame front and the formation of a tulip-shaped flame with the location of the tulip petal tips on the inner edges of the boundary layer. It is worth noting that any intrinsic flame front instabilities are not involved in the process in agreement with Ponizy, Claverie, and Veyssière (2014) experimental observations.

It should be emphasized that Guénoche (1964) put forward the fruitful hypothesis that during the deceleration phase the flame generates rarefaction waves, which are the key factor in the tulip flame formation. Unfortunately, 70 years ago, the level of numerical simulations and experimental methods was insufficient to confirm Guénoche’s hypothesis, and later researchers either forgot or ignored his work. Nevertheless, we want to emphasize Guénoche’s priority in describing the mechanism of tulip flame formation. The results of the present work confirm Guénoche’s hypothesis and for the first time not only qualitatively explain the physical mechanism leading to the inversion of the flame front, but also the mechanism of formation of a particular form of tulip flame, such as the location of the tulip petals.

Dynamics of the flame propagating in a channel with no-slip walls

The front of the flame, which is ignited on the center line near the closed end of the channel and propagates to the opposite closed or open end, acquires a convex shape with the tip pointing toward the unburned gas. The dynamics of a premixed flame in a tube depends on the boundary conditions and the tube width. Clanet and Searby (1996) distinguished four stages in the flame evolution: the ignition of a spherical flame, the finger-like flame, the quenching of the lateral flame skirt at the side walls, and the tulip flame formation. The theoretical model developed by Bychkov et al. (2007), using an analytical solution of the Navier – Stokes equations for the two-dimensional case, confirm at least the first two stages of Clanet and Searby classification, but does not explain the deceleration stage and the physical mechanism of the tulip flame formation.

The flow of unburned gas induced by combustion and the boundary layer formed due to no-slip boundary conditions on the side walls are important factors for the tulip-shaped flame formation. After ignition of the flame near the closed end of the tube, the flame front expands rapidly, taking on an almost hemispherical (semi-circular in the 2D case) shape. At the ignition stage the flame propagates with velocity $(\mathrm{\Theta }-1)U_f$ and, according to Clanet and Searby (1996), after a short time $t_{sph}\sim D/2(\mathrm{\Theta }-1)U_f$, begins the next stage of a finger-shaped flame, when the flame accelerates exponentially due to the increase of the flame lateral surface area

(1) $X_{tip}\propto exp(2\mathrm{\Theta }{U}_{f}t/D),$(1)

where $D$ is the channel width, $\mathrm{\Theta }={\rho }_u/{\rho }_b$ is the ratio of the densities of unburned ${\rho }_u$ and burned ${\rho }_b$ gases, and $U_f$ is the laminar flame speed.

Depending on the size of the initially ignited flame (hot spot) and the width of the tube, the finger-shaped stage lasts a short time and ends when the lateral sides of the flame skirt touch the side walls of the channel and extinguish. From the beginning, the thermal expansion of the high temperature combustion products pushes the unburned gas toward the opposite end of the tube creating a flow of unburned gas ahead of the flame, resulting in the flame acceleration compatible with boundary conditions. The velocity of each point of the flame front in the laboratory reference frame is the sum of the laminar velocity of the flame $U_f$ relative to the unreacted mixture and the local axial flow velocity of the unburned gas in the near-field zone ahead of that point $U_{+}(x,y)$ with which the flame is carried by the flow

(2) $U_{fl}=U_f+U_{+}(x,y).$(2)

The velocity $U_{+}(x,y)$ is almost constant in the bulk cross section of the channel and drops to zero at the side walls inside the boundary layer of thickness ${\delta }_l\approx 5X/\sqrt{Re}$, where $Re=U_{+}X/\nu$ is the Reynolds number, $X$ is the coordinate along the tube, $\nu$ is kinematic viscosity.

The shape of the flame front resembles the velocity profile of the unburned gas, remaining almost flat in the bulk (effectively retaining the convex shape obtained during the ignition stage), with the trailing edges of the flame skirt extending backward in the boundary layer. Therefore, the flame surface area increases, the flame consumes fresh fuel over a larger surface area, resulting in an increase in the combustion wave speed. In the framework of the thin flame model, the increase of the burning rate is proportional to the relative increase of the flame surface area (length in 2D case), which grows linearly in time with accuracy ${\delta }_l/D<<1$ and which in turn gives rise to a larger gradient of the velocity field. This means that the speed of the flame increases exponentially as was for the first time shown by Liberman et al. (2010):

(3) $U_{fl}\propto exp(\alpha \mathrm{\Theta }{U}_{f}t/D)$(3)

where $\alpha$ is a numerical factor of the order of unity. Expression (3) is similar to Equation 1(1) $X_{tip}\propto exp(2\mathrm{\Theta }{U}_{f}t/D),$(1) , so it can be difficult to distinguish between the finger-shaped exponential stage (1) and the exponential growth of the flame speed (3) from one another in the experimental Schlieren images and in numerical studies.

The lateral part of the flame skirt, stretched back along the tube wall in the boundary layer, and the side wall of the tube form a gap (fold) filled with unburned gas. As the flame speed increases, the flow velocity in the unburnt gas also increases and the thickness of the boundary layer decreases. When the thickness of the boundary layer becomes sufficiently small (of order $L_f$), the temperature of the unburned gas inside the fold increases significantly and the unburned fuel inside the fold will be burned out very fast. This leads to the collapse of the lateral part of the flame skirt and a sharp decrease in the flame speed. In the case of isothermal boundary conditions, the scenario is similar: when the boundary layer thickness becomes small enough, the lateral part of the flame skirt comes close to the cold side wall of the tube and extinguishes. In the case of adiabatic boundary conditions, when the thickness of the boundary layer decreases, the part of the unburned gas near the top of the fold is “overheated” first, and the flame skirt at the top of the fold begins to spread much faster toward the side wall, where it collapses. In a sense this resembles a cumulative jet and looks like a small explosion that generates weak pressure waves diverging from the side wall. The effect of these pressure waves is small compared to the rarefaction waves produced by the decelerating flame. The rarefaction wave propagates forward and significantly reduces and can even reverse the velocity of the unburned flow ahead of the flame, resulting in an increase in the thickness of the boundary layer, creating a tulip-shaped of the axial velocity profile in the unburned gas, which leads to the inversion of the flame front. A tulip-shaped axial velocity profile in the unburned gas in the near-field zone ahead of the flame means that according to Equation 2(2) $U_{fl}=U_f+U_{+}(x,y).$(2) the local velocity of the flame front is minimal at the axis ($y=0$), gradually growing toward the walls of the channel, and decreases inside the boundary where $U_{+}(x,y)$ drops to zero on the tube walls at $y=\pm D/2$. This is the mechanism of tulip-shaped flame formation, which also predicts the shape of the tulip flame: the height and location of the tips of the tulip petals.

Dunn-Rankin (2009) noticed that for a tulip flame to form, the width of the tube should be large enough (>25 mm) to obviate the wall-heat transfer and boundary-layer effects, but small enough (<100 mm) to prevent large buoyancy effects. In fact, the wall-heat transfer and boundary-layer effects are important factors for the tulip flame formation. We show that for a tulip flame to form there is another limitation for the smallest width of the channel. The main difference in the dynamics of flame propagation in a “wide” and in a “narrow” tube is due to the different time required for establishing a Poiseuille flow with a parabolic velocity profile in the unburned gas ahead of the flame. According to Liberman et al. (2010) and ZeldovichYa (1947), the characteristic time for establishing a Poiseuille flow ahead of the flame is $t_P\approx D^2/100\nu$. For example, in a hydrogen-air ($\nu =0.22c{m}^2/s$) a parabolic velocity profile in the flow ahead of the flame in the “wide” tube $D=1c\mathrm{m}$ establishes in time $t_P\approx 50\mathrm{m}\mathrm{s}$, which is much longer than the time of the tulip flame formation, which is about (1–5) ms. On the other hand, in a narrow tube $D=1\mathrm{m}\mathrm{m}$, the parabolic velocity profile is established in $t_P\approx 0.5\mathrm{m}\mathrm{s}$. It can be shown, (Liberman 2021), that in this case the flame velocity increases exponentially up to the transition to detonation, bypassing the deceleration stage. It is for this reason that tulip flames have never been observed in narrow tubes. An important feature of flame dynamics in a tube with non-slip walls is that an accelerating flame generates pressure waves and increases the flow velocity ahead of the flame, while a decelerating flame generates rarefaction waves and decreases the flow velocity ahead of the flame. It should be emphasized that, unlike a stationary flame, a flow with an accelerating or decelerating flame is not isobaric. During the acceleration phase, the pressure increases at about the same rate as the flame speed, creating pressure waves and negative pressure gradients that support the flow ahead of the accelerating flame. The rarefaction wave causes a decrease in the flow speed in the unburned gas and creates a reverse (positive) pressure gradient.

Numerical models

Two-dimensional DNS

The two-dimensional computational domains, which were modeled using high resolution DNS, represented rectangular channels with both closed ends of width $D=1\mathrm{c}\mathrm{m}$ and aspect ratios of $L/D=6,12,18$, and a rectangular channel of width 1 cm with an open right end. The numerical simulations solve the 2D time-dependent, reactive compressible Navier–Stokes equations including molecular diffusion, thermal conduction, viscosity and detailed chemical kinetics. The governing equations are

(4) $\frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}+\frac{\mathrm{\partial }(\rho u_i)}{\mathrm{\partial }{x}_i}=0,$(4)

(5) $\frac{\mathrm{\partial }(\rho u_i)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }(P{\delta }_{ij}+\rho {u_i}^2)}{\mathrm{\partial }{x}_j}=\frac{\mathrm{\partial }{\tau }_{ij}}{\mathrm{\partial }{x}_j},$(5)

(6) $\frac{\mathrm{\partial }(\rho E)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }\left[(\rho E+P)u_i\right]}{\mathrm{\partial }{x}_i}=\frac{\mathrm{\partial }({\tau }_{ij}{u}_i)}{\mathrm{\partial }{x}_i}-\frac{\mathrm{\partial }{q}_i}{\mathrm{\partial }{x}_i},$(6)

(7) $\frac{\mathrm{\partial }\rho Y_k}{\mathrm{\partial }t}+\frac{\mathrm{\partial }\rho u_i{Y}_k}{\mathrm{\partial }{x}_i}=\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_i}\left(\rho V_{ik}{Y}_k\right)+{\dot{\omega }}_k,$(7)

(8) $P=\rho R_{B}T\left({\sum }_{i=1}^{N_s}\frac{Y_i}{W_i}\right),$(8)

Here $\rho$, $u_i$, $T$, $P$, $E$, ${\tau }_{ij}$, $q_i$ are density, velocity components, temperature, pressure, specific total energy, viscosity stress, heat flux mass, $R_B$ is the universal gas constant. $Y_i,W_i,V_{i,j}$ are mass fraction, molar mass and diffusion velocity of species $i$. The viscosity and thermal conductivity of pure-species and mixture are calculated by the Chapman-Enskog expression Byron Byron Bird and Stewart Lightfoot (2002) and a semi-empirical formula Wilke (1950).

The reaction rate of species $i$ is determined as

(9) ${\dot{\omega }}_k=W_i\sum _{j=1}^{N_r}({v{ }^{\mathrm{\prime }\mathrm{\prime }}}_{jk}-{v{ }^{\mathrm{\prime }}}_{jk})\cdot \left(k_{f,j}\prod _{k=1}^{N_s}{\left(\frac{\rho Y_k}{W_k}\right)}^{{v{ }^{\mathrm{\prime }}}_{jk}}-k_{b,j}\prod _{k=1}^{N_s}{\left(\frac{\rho Y_k}{W_k}\right)}^{{v{ }^{\mathrm{\prime }\mathrm{\prime }}}_{jk}}\right)$(9)

where $v{{ }^{\mathrm{\prime }}}_{jk}$ and $v{{ }^{\mathrm{\prime }\mathrm{\prime }}}_{jk}$ are stoichiometric coefficients of species $k$ of the reactant and product sides of reaction $j$. A detailed chemical kinetic model for hydrogen/air developed by Kéromnès et al. (2013) was implemented in the simulations. This chemical model has been extensively tested against experimental data and found to accurately predict ignition delay and laminar flame velocity over a wide range of pressures (1–70 bar), temperatures (900–2500 K), and equivalence ratio (0.1–4.0). The initial conditions are $P_0=1atm$, $T_0=298K$, and during the formation of a tulip flame, the pressure increases from an initial value of $P_0=1atm$ up to $\approx 1.5\mathrm{a}\mathrm{t}\mathrm{m}$ for high aspect ratio channels and up to $\sim 2.0\mathrm{a}\mathrm{t}\mathrm{m}$ in a short channel $L/D=6$.

Three-dimensional LES

The three-dimensional DNS resolving properly all characteristic scales is currently very expensive. An efficient approach is to use Large Eddy Simulation (LES). A comparison of the flame dynamics obtained with DNS for a flame propagating in a two-dimensional channel with the flame dynamics obtained with LES is shown in Appendix A. It is seen that qualitatively the results of DNS and LES are reasonably close.

In LES, all quantities are filtered in the spectral space by applying filter to the governing equations. Since the flow is compressible, the mass-weighted Favre is introduced. The filtered governing equations are

(10) $\frac{\mathrm{\partial }\bar{\rho }}{\mathrm{\partial }t}+\frac{\mathrm{\partial }\bar{\rho }{\tilde{u}}_i}{\mathrm{\partial }{x}_i}=0,$(10)

(11) $\frac{\mathrm{\partial }\bar{\rho }{\tilde{u}}_i}{\mathrm{\partial }t}+\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}\left(\bar{\rho }{\tilde{u}}_i{\tilde{u}}_j+\bar{P}{\delta }_{ij}\right)=\frac{\mathrm{\partial }{\bar{\tau }}_{ij}}{\mathrm{\partial }{x}_j}-\frac{\mathrm{\partial }{\tau }_{ij}^{sgs}}{\mathrm{\partial }{x}_j},$(11)

(12) $\frac{\mathrm{\partial }\bar{\rho }\tilde{E}}{\mathrm{\partial }t}+\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_i}(\bar{\rho }\tilde{E}+\bar{P}){\tilde{u}}_i=\frac{\mathrm{\partial }({\tilde{u}}_i{\bar{\tau }}_{ij})}{\mathrm{\partial }{x}_i}-\frac{\mathrm{\partial }{\bar{q}}_i}{\mathrm{\partial }{x}_i}-\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_i}{H}_i^{sgs}+\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_i}{\sigma }_{ij}^{sgs},$(12)

(13) $\frac{\mathrm{\partial }\bar{\rho }{\tilde{Y}}_k}{\mathrm{\partial }t}+\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_i}(\bar{\rho }{\tilde{Y}}_k{\tilde{u}}_i)=\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_i}\left(\bar{\rho }{\tilde{D}}_i\frac{\mathrm{\partial }{\tilde{Y}}_k}{\mathrm{\partial }{x}_i}\right)-\frac{\mathrm{\partial }{\phi }_i^{sgs}}{\mathrm{\partial }{x}_i}+{\bar{\omega }}_k$(13)

The unresolved fluxes ${\sigma }_{ij}^{sgs}$ are negligible, and the unresolved Reynolds stresses ${\tau }_{ij}^{sgs}$ are modeled according to the Boussineq assumption

(14) ${\tau }_{ij}^{sgs}-\frac{1}{3}{\delta }_{ij}{\tau }_{kk}^{sgs}=-\rho v_t\left(\frac{\mathrm{\partial }{\tilde{u}}_i}{\mathrm{\partial }{x}_j}+\frac{\mathrm{\partial }{\tilde{u}}_j}{\mathrm{\partial }{x}_i}-\frac{2}{3}{\delta }_{ij}\frac{\mathrm{\partial }{\tilde{u}}_k}{\mathrm{\partial }{x}_k}\right).$(14)

The sub-grid viscosity $v_t$ is calculated by wall-adapting local eddy viscosity model and the sub-grid heat flux is defined by Garnier, Adams, and Sagaut (2009)

(15) $H_i^{sgs}=-\bar{\rho }\frac{v_t{c}_p}{{Pr}_t}\frac{\mathrm{\partial }\tilde{T}}{\mathrm{\partial }{x}_i},$(15)

where $\underset{t}{Pr}=0.75$ sub-grid Prandtl number.

In simulations of turbulent reactive flow, the most important is the combustion model, which is used to approximate the filtered reaction rates. In this paper, we employ the so-called thickened flame model Colin et al. (2000). The flame front is thickened artificially without changing the laminar flame velocity. For a premixed laminar flame, the laminar flame velocity $U_f$ and the flame thickness $L_f$ can be expressed as

(16) $U_f\propto \sqrt{D_{th}A},L_f\propto \frac{D_{th}}{U_f}=\sqrt{D_{th}/A},$(16)

where $D_{th}$ and $A$ are thermal diffusivity and pre-exponential factor. Supposing that $D_{th}$ and $A$ are replaced by $F{D}_{th}$ and $A/F$, it is easy to check that $L_f$ increases by a factor $F$ while $U_f$ remains unchanged. The wrinkling factor function $\epsilon$ is introduced to simulate a flame front wrinkled by vortices at the subgrid scale.

The power-law wrinkling factor function derived by Charlette, Meneveau, and Veynante (2002) is used in simulations. Appling the thickened flame model, the filtered Equation 13(13) $\frac{\mathrm{\partial }\bar{\rho }{\tilde{Y}}_k}{\mathrm{\partial }t}+\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_i}(\bar{\rho }{\tilde{Y}}_k{\tilde{u}}_i)=\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_i}\left(\bar{\rho }{\tilde{D}}_i\frac{\mathrm{\partial }{\tilde{Y}}_k}{\mathrm{\partial }{x}_i}\right)-\frac{\mathrm{\partial }{\phi }_i^{sgs}}{\mathrm{\partial }{x}_i}+{\bar{\omega }}_k$(13) reads

(17) $\frac{\mathrm{\partial }\bar{\rho }\tilde{Y}}{\mathrm{\partial }t}+\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_i}(\bar{\rho }\tilde{Y}{\tilde{u}}_i)=\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_i}\left(\bar{\rho }\epsilon {\tilde{D}}_i\frac{\mathrm{\partial }\tilde{Y}}{\mathrm{\partial }{x}_i}\right)+\frac{\epsilon \bar{\omega }}{F}$(17)

In both 2D and 3D simulations we used adiabatic no-slip reflecting boundary conditions at the walls of the tube: $\vec{u}=0,\mathrm{\partial }T/\mathrm{\partial }\vec{n}=\mathrm{\partial }{Y}_k/\mathrm{\partial }\vec{n}=0$, where $\vec{n}$ is the normal to the wall.

Numerical schemes and modeling parameters

The two-dimensional direct numerical simulations were performed using the DNS solver, which uses the fifth order weighted essentially non-oscillatory (WENO) finite difference schemes (Jiang and Shu 1996) to resolve the convection terms of the governing equations. The advantage of the WENO finite difference method is the capability to achieve arbitrarily high order accuracy in smooth regions while capturing sharp discontinuity. To ensure the conservation of the numerical solutions, the fourth order conservative central difference scheme (Rango and Zing 1999) is used to discretize the non-linear diffusion terms. The time integration is third order strong stability preserving Runge – Kutta method by Gottlieb and Shu (1998). The resolution and convergence (a grid independence) tests similar to that in our previous publications Liberman et al. (2019) and Wang et al. (2018) were thoroughly performed to ensure that the resolution is adequate to capture details of the problem in question and to avoid computational artifacts.

The reliable simulations of reactive flows require a proper resolution of the inner structure of a flame. According to Liberman et al. (2019), a resolution of 8 grids per flame front is sufficient to predict the correct laminar flame structure and the flame velocity. Since the thickness of the laminar flame decreases with increasing pressure, higher resolution is required to determine the flame structure at maximum pressure. The maximum pressure during the formation of the tulip flame is about 2 bar, which leads to decrease in the thickness of the flame front from $L_f=350\mu \mathrm{m}$ to $L_f\simeq 130\mu \mathrm{m}$. In two-dimensional simulations we used uniform mesh with resolution $\mathrm{\Delta }x\simeq 12.5\mu \mathrm{m}$, which corresponds to 28 grid points over the flame width at the beginning of the process, and over 14 grid points at the maximum pressure. In the three-dimensional large eddy simulations, we used a coarse mesh of width $\mathrm{\Delta }x\simeq 125\mu \mathrm{m}$ defined as a function of pressure which ensures that the flame front is well resolved at all pressures during the tulip flame formation.

Two-dimensional and three-dimensional numerical simulations were performed to clarify what physical processes lead to flattering and inversion of the flame front during the formation of a tulip flame. Two-dimensional calculations were performed for rectangular channels of the width $D=1\mathrm{c}\mathrm{m}$ with both closed ends and various aspect ratios $L/D=6,12,18$ and for the half-open channel with the right opened end. The flame was ignited at the center of the left closed end of the tube by a small hot spot of the burnt gas behind a half-circular flame of the radius 3 mm in the 2D case and by a hemispherical flame of the radius 3 mm in the 3D case with a hot, burnt gas behind the flame front. The no-slip adiabatic boundary conditions were used at the tube walls.

Tulip flame formation in a two-dimensional rectangular channel L/D = 6

Figure 1 shows the time evolution of the flame skirt surface area $F_f$, local velocities of the flame front along the center line $U_{fl}(y=0)$ and near the side wall $U_{fl}(y=0.4\mathrm{c}\mathrm{m})$ and the combustion wave velocity (burning rate) $S_f$. It can be seen in Figure 1 that the flame propagation speed in the laboratory reference frame is closely related to the variation of the flame surface area. From the beginning and up to 0.5 ms the flame surface area $F_f$ as well as combustion wave velocity $S_f$ and the local velocities of the flame front $U_{fl}(y=0)$ and $U_{fl}(y=0.4\mathrm{c}\mathrm{m})$ increase. When the lateral parts the flame skirt begin to collapse at the side wall, the flame surface area $F_f$ begins to decrease at 0.5 ms, at the same time the combustion wave velocity $S_f$ begins to decrease, but the velocity of the flame tip (the point on the flame front at the center line) begins to decrease after 0.55 ms. However, the speed of the flame front near the wall, at $y=0.4\mathrm{c}\mathrm{m}$, continues to increase and after 0.9 ms it exceeds $U_{fl}(y=0)$. The change in the flame acceleration rate at (0.3–0.4) ms is due to the collision of the flame with the pressure wave reflected from the right end of the channel. This collision slows down the flame and, together with the increased pressure in the unburnt gas mixture, explains the higher rate of flame deceleration in a short channel compared to channels with a higher aspect ratio. Figure 2 shows a sequence of numerical images of Schlieren patterns and the stream lines for selected times, which show the evolution and dynamics of the flame and flows in the channel $\mathrm{D}=1\mathrm{c}\mathrm{m}\mathrm{a}\mathrm{n}\mathrm{d}$$\mathrm{L}/\mathrm{D}=6$ during the tulip-shaped flame development.

Figure 1. Time evolution of the flame surface area $F_f$ (dashed-dotted), local velocities of the flame front at $y=0$ and at $y=0.4cm$ and combustion wave velocity $S_f$; tube $L/D=6$.

Figure 2. (a, b). Time sequence of computed Schlieren images and streamlines for the flame propagating in the tube $L/D=6$ with both ends closed.

It is seen in Schlieren images (Figure 2) that oblique pressure waves generated by a convex flame front at the stage of flame acceleration, propagate through the tube, repeatedly colliding and reflecting from the sidewalls, reflect from the right end of the tube and run back and forth along the tube. This results in slight fluctuations in pressure and flame speed. The decelerating flame acts like a piston, which moves out of the pipe creating a rarefaction wave in the unburned gas ahead of the flame front. It can be seen from Figure 2 that the reverse flow in the burnt gas is formed in a short time after the beginning of flame deceleration, which is due to the rarefaction wave created by the decelerating flame in the unburned gas.

An important consequence is that the flow velocity in the near-field zone ahead of the flame front and, consequently, the local velocity in the flame front decrease more at the center line y=0 and lesser toward the channel wall up to the inner edge of the boundary layer, at y=0.4cm.. It should be emphasized that the rarefaction wave reduces the velocity and changes the structure of the flow in the near-field zone ahead of the flame, which leads to a significant increase in the thickness of the boundary layer and changes the axial velocity profile, which takes the shape of a tulip or an inverted tulip (see Figure 5). Due to the reverse (negative) flow created behind the rarefaction wave in the near-field region of the unburnt gas, the difference between velocities of the unburnt flow near the wall and at the axis increases.

Figures 3(a,b) show the time evolution of the local velocities of the flame front on the center line $U_{fl}(y=0)$ and near the boundary layer $U_{fl}(y=0.4\mathrm{c}\mathrm{m})$, and the flow velocity in the unburned gas at 1 mm and 3 mm ahead of the flame front at $y=0$ and at $y=0.4\mathrm{c}\mathrm{m}$. It can be seen that the local velocities of the flame front change synchronously with the change in the axial flow velocities in the near-field zone ahead of the flame. As a result, the central part of the flame front begins to lag behind compared to the lateral parts of the flame.

Figure 3. A) Flame front velocity at $y=0$ (dashed line) and flow velocity at 1 mm (▲) and 3 mm (●) ahead of the flame front. b) Flame front velocity at $y=0.4\mathrm{c}\mathrm{m}$ (dashed line) and flow velocity 1 mm (▲) and 3 mm (●) ahead of the flame front.

The rarefaction wave, which is generated by the decelerating flame in the unburned mixture, consists of the head which propagates toward the right end of the channel at the speed of sound, with the flow behind, which moves in the opposite (negative) x-direction. Contrary to the accelerating flame which creates negative pressure gradients supporting the unburned flow ahead of the flame, the refraction wave creates a positive pressure gradient.¹ During the development of the tulip flame, the velocity and accelerations in the different parts of unburned flow change from positive to negative values. This leads to the formation of a complex unsteady flow ahead of the flame, when the flow velocity $U_{+}$ in the near-field zone ahead of the flame front can become negative with the maximum value at the axis of the tube.

Figure 4(a, b) shows the difference between the axial flow velocities at 1 mm ahead of the flame front at $y=0;0.25\mathrm{c}\mathrm{m};0.4\mathrm{c}\mathrm{m}$. After 0.6 ms, $\mathrm{\Delta }{U}_{+}=U_{+}(y=0.4)-U_{+}(y=0)$ becomes positive, and since then, the parts of the flame front that are closer to the center line propagate in the laboratory system at lower velocities compared to the local velocities of the parts of the flame front that are closer to the side walls.

Figure 4. Velocities of the flow at 1 mm ahead of the flame front at $y=0;0.4;0.25\mathrm{c}\mathrm{m}$, and the difference $\mathrm{\Delta }{U}_{+}$ between (a): $U_{+}(y=0.4)$ and $U_{+}(y=0)$, and between (b): $U_{+}(y=0.25)$ and $U_{+}(y=0)$.

It should be emphasized two important factors involved in the initiation of reverse flow behind the rarefaction wave. The reverse flow behind the head of the rarefaction wave interacts with the previously created during the acceleration stage “positive” flow, which leads to the significant decrease in the flow velocity and to the increase of the boundary layer thickness, which is maximal in the near-field zone ahead of the flame front. Figure 5 shows the temporal evolution of the axial velocity profiles, shown by solid lines (the axial velocity distribution in y-direction) in the unburned gas mixture ahead of the flame, in the cross section of the tube shown by the dashed lines. After 0.70 ms the velocity profile acquires the shape of a tulip, thus creating conditions for the flame front inversion and the tulip-shaped flame formation.

Figure 5. Axial velocity profiles (solid lines) in the flow ahead of the flame in different cross sections of the tube (L/D=6) shown by the dashed lines.

Since the local velocity of the flame front is the sum of the laminar flame velocity relative to the unburned mixture and the flow velocity in the near-field zone ahead of the flame front, the change in the local velocity of various parts of the flame front in the laboratory reference system is as greater as closer the part of the flame front to the tube axis. In fact, the tulip flame has already formed at 1.20ms, and at 1.35 ms, the surface area of the flame begins to increase, and this is so-called distorted tulip flame.

It can be seen in Figure 2 that the stage of flame front flattening, which begins at 0.75 ms, is accompanied by the formation of two vortices that appeared in the burnt gases near the sidewalls of the tube. According to the theoretical model developed by Matalon and Metzener (1997, 2001) vortices play an important role in the formation of a tulip flame. The role of a pair of vortices, created by the highly curved parts of the flame near the sidewalls, apparently due to the baroclinic effect, have been discussed by many authors e.g. by Dunn-Rankin, Barr, and Sawyer (1986) and Dunn-Rankin and Sawyer (1998). Matalon and Metzener (1997) and (2001Metzener and Matalon 2001) put forward a seemingly plausible scenario that: “The highly curved segments of the flame near the walls generate a pair of intense vortices in the burned gas that advects the flame into the tulip shape.” We cannot confirm this scenario and believe that the formation of vortices is only a consequence of the flow in the unburned and in the burned gases which arises due to the rarefaction wave during the development of the tulip-shaped flame, but not a cause of a tulip flame formation. Moreover, the height, $h={\int }_{0.8ms}^{1.35ms}\{U_{+}(y=4mm,t)-U_{+}(y=0,t)\}dt=1.56mm$(the difference in x coordinates of the flame front at $y=0.4mm$ at the tulip petal tip location and at the central line, $y=0$) is completely determined by the axial velocity profiles of the unburned flow in the close proximity ahead of the flame front. The height of the tulip petal mesured in computed Schlieren images in Fig. 2 is $h=1.58\pm 0.01mm$.

It is known (Cohen and Kundu 1989) that the baroclinic effect leads to the generation of vorticity when a curved interface between two matters of different densities is exposed to a pressure gradient. The baroclinic displacement is determined by the inviscid two-dimensional vorticity equation: $d\omega /dt=\left[\mathrm{\nabla }\rho \times \mathrm{\nabla }P\right]/{\rho }^2$.

When the flame front begins to flatten and the curvature of the flame front near the side walls becomes maximum, the conditions become favorable for baroclinic vortex formation. We emphasize that the boundary conditions require the occurrence of recirculation in the reverse flow of combustion products near the left closed end of the channel. These large-scale vortices subsequently drift in the combustion products toward the flame front, forming a focus of the flow velocity streamlines near the axis of the tube, that apparently extend in both directions from the flame front, as can be seen in Figure 2 at 0.85 ms and later. Obviously, this is a purely hydrodynamic process, which does not involve the instabilities of the flame front. Thus, we conclude that the physical process that explains and fully describes the formation of a tulip-shaped flame is a rarefaction wave generated by the flame during the deceleration stage. Flow features, such as vortices or a stagnation point (saddle point) along the center line, are formed when the flame front is inverted due to the flow created by rarefaction waves, and are a consequence, not a cause, of the formation of a tulip flame.

Figures 6(a, b) show pressure $P(x,y=0)$ and flow velocity $U_{+}(x,y=0)$ profiles along the center line and near the boundary layer $U_{+}(x,y=0.4cm)$, for selected times during the tulip flame formation. The arrows at the Figures 6 show the direction of waves propagation. To trace the temporal evolution of the pressure gradient, additional pressure profiles are shown shortly before and shortly after the time shown in the upper right corner of each slide. It is seen that from the beginning ($t<0.6ms$), the flow velocity at the center line, in the near-field zone of unburned gas ahead of the flame is larger than the velocity near the sidewall. Therefore, the local speed of the flame front at the central line is larger than the flame front velocities near the wall.

Figure 6. (a) Sequences of pressure and velocity distribution along the center line, and the flow velocity profiles $U_{+}(x,y=0)$ and $U_{+}(x,y=0.4)$ before the flame front flattering; (b) Sequences of pressure and velocity distribution along the center line, and the flow velocity profiles $U_{+}(x,y=0)$ and $U_{+}(x,y=0.4)$ during the tulip flame formation.

After 0.6 ms, the rarefaction wave creates a reverse flow to the left end of the channel in the unreacted gas, resulting in an increase in the thickness of the boundary layer in the near-field flow ahead of the flame as can be seen Figure 5. This, in turn, causes a reverse flow of the burned gas, which is maintained by creating a reverse (positive) pressure gradient at 0.7 ms in Figure 6 (а). The formation of a reverse (positive) pressure gradient during the inversion of the flame front is a natural consequence of the rarefaction wave structure (Landau and Lifshitz 1989). Figure 6(b) shows the time evolution of the pressure $P(x,y=0)$ and the flow velocity $U_{+}(x,y=0)$ profiles along the center line and near the boundary layer $U_{+}(x,y=0.4cm)$ during the late times of the tulip flame formation. It can be seen that after 0.85 ms a positive pressure gradient is almost “permanently” established, although sometimes it is smoothed out by reflected pressure waves. At the same time, the flow velocity near the boundary layer $U_{+}(x,y=0.4cm)$ in the near-field zone of the unburned flow exceeds the flow velocity at the center line $U_{+}(x,y=0)$ all the way up to the formation of a tulip flame.

Effect of aspect ratio L/D = 18

The velocity of the flame tip at the center line $y=0$ computed for the flames propagating in tubes of various aspect ratios $L/D=6,12,18$ with both closed ends and for the flame in the half-open tube (open right end) are shown in Figure 7.

Figure 7. The velocity of the flame tip along the center line $y=0$, computed for tubes $L/D=6,12,18$ and the half-open tube.

The arrows in Figure 7 show the change of the flame tip velocity caused by the collision of the flame with the pressure wave reflected from the right closed end of the tube, which occurs at$t\approx 0.29\mathrm{m}\mathrm{s},0.62\mathrm{m}\mathrm{s},0.84\mathrm{m}\mathrm{s}$, for $L/D=6,12,18$, respectively.

As already mentioned, the most important is the deceleration stage, which results in a rarefaction wave, that reduces the velocity of the unburned flow, which leads to an increase in the width of the boundary layer in the near-field zone ahead of the flame and the tulip-shaped form of the axial velocity profile, thereby initiating the formation of a tulip-shaped flame. Interesting is to compare the intensity of the rarefaction wave produced by a decelerating flame in channels with different aspect ratios, in a semi-open channel and in a three-dimensional channel. For the classical problem (Landau and Lifshitz 1989), when the rarefaction wave is created by a piston, which starts to move out of the tube at a constant acceleration $U=-at$, the intensity of the rarefaction wave is characterized by the value of $a=\mathrm{c}\mathrm{o}nst.$. In our case the “piston” acceleration is not constant, and the intensity of the rarefaction wave can be characterized by the magnitude of the “negative” flame acceleration at the deceleration stage, the absolute values of which at the beginning are given in Table 1.

Table 1. Absolute values of the flame acceleration.

$L/D=6$	$a\left(0.7ms\right)\hspace{0.17em}=\hspace{0.17em}9\times {10}^4\mathrm{m}/{\mathrm{s}}^2$	$a\left(0.9ms\right)\hspace{0.17em}=\hspace{0.17em}1.65\times {10}^5\mathrm{m}/{\mathrm{s}}^2$
$L/D=12$	$a\left(0.7ms\right)\hspace{0.17em}=\hspace{0.17em}8\times {10}^4\mathrm{m}/{\mathrm{s}}^2$	$a\left(0.9ms\right)\hspace{0.17em}=1.3\times {10}^5\mathrm{m}/{\mathrm{s}}^2$
$L/D=18$	$a\left(0.7ms\right)\hspace{0.17em}=5.5\times {10}^4\mathrm{m}/{\mathrm{s}}^2$	$a\left(0.9ms\right)\hspace{0.17em}=2.1\times {10}^5\mathrm{m}/{\mathrm{s}}^2$
Half-open channel	$a\left(0.7ms\right)\hspace{0.17em}=5.5\times {10}^4\mathrm{m}/{\mathrm{s}}^2$	$a\left(1.2ms\right)\hspace{0.17em}=2.1\times {10}^5\mathrm{m}/{\mathrm{s}}^2$
$L/D=6$, 3D channel	$a\left(0.5ms\right)\hspace{0.17em}=1.5\times {10}^6\mathrm{m}/{\mathrm{s}}^2$

Table 1 shows that the magnitude of the flame deceleration and the corresponding intensity of the rarefaction wave are greater when the aspect ratio of the channel is smaller, since for a channel with a smaller aspect ratio, the pressure in the unburned gas grows faster. As will be shown in Section 4.4, the acceleration and deceleration stages in the case of a three-dimensional flame occur much faster than in the case of a two-dimensional flame, and, accordingly, the intensity of the rarefaction wave in the three-dimensional problem is much greater.

Figure 8 shows the time evolution of the flame front surface area $F_f$, the combustion wave velocity $S_f$, and local velocities of the flame front at the center line, $U_{fl}(y=0)$ and near the side wall, $U_{fl}(y=0.4\mathrm{c}\mathrm{m})$. It is seen that the flame surface begins to decrease at 0.55 ms when the trailing edge of the flame skirt touched the side wall. The combustion wave velocity $S_f$ and the velocity of the flame tip $U_{fl}(y=0)$ begin to decrease at 0.57 ms. As explained earlier, the rate of the flame deceleration in this case is lower than in the shorter channel. Therefore, the rarefaction wave is weaker and the flame front begins to flatten later. After 0.95 ms, the deceleration rate increases as the flame surface area decreases. The local speed of the flame front at the center line $U_{fl}(y=0)$ decrease, while the local speed of the flame front near the boundary layer $U_{fl}(y=0.4\mathrm{c}\mathrm{m})$ increases. After 1.0 ms, the local velocity of the flame front at y = 0.4 cm exceeds the velocity of the flame front at y = 0, and from this moment the inversion of the flame front begins. The entire scenario of the flame front inversion and the tulip flame formation is qualitatively similar to that considered in previous section for the short channel $L/D=6$. But since there is no collisions of the pressure wave reflected from the right closed end with the flame front, the key role of the rarefaction wave is better seen for channels with larger aspect ratio.

Figure 8. Evolution of combustion wave velocity $S_f$, the flame surface area $F_f$, and local velocities of the flame front at $y=0$ and at $y=0.4\mathrm{c}\mathrm{m}$ for the tube $L=18\mathrm{c}\mathrm{m}$.

Figure 9 shows sequences of calculated schlieren images and streamlines at selected times during the development of a tulip flame for the tube with aspect ratio $L/D=18$.

Figure 9. (a,b) Schlieren images and stream lines for the premixed hydrogen/air flame propagating in the tube with both closed ends, aspect ratio $L/D=18$.

Figure 10 (a, b) shows the evolution of axial velocities in the unburned gas at the distance of 1 mm ahead of the flame front and the difference $\mathrm{\Delta }{U}_{+}$ between the velocity near the side walls and the velocity at the center line. The values of $\mathrm{\Delta }{U}_{+}$ become positive after 0.95 ms, and reach maximum, when the rarefaction wave reverse the flow of the unburned gas ahead of the flame.

Figure 10. (a, b) Velocities of the flow at 1 mm ahead of the flame front at $y=0,0.4,0.25\mathrm{c}\mathrm{m}$ and the difference $\mathrm{\Delta }{U}_{+}$ between $U_{+}(y=0.4),U_{+}(y=0.25)$ and $U_{+}(y=0)$; $L/D=18$.

As mentioned earlier, in this case the streamlines ahead and behind the flame front are “focused” to the center line. The convergence of streamlines is a feature of the flow, which develops due to the tulip-like shape of the axial velocity profile in the unburned gas and large-scale vortex motion in the flow of burnt products. Figure 11 shows the time evolution of the axial velocity profiles in different cross section of the channel in the unburned flow ahead of the flame. It can be seen that due to the rarefaction waves generated by the decelerating flame the axial velocity profile in close vicinity ahead of the flame front acquires a tulip-like form creating conditions for the flame front inversion and the formation of the tulip flame.

Figure 11. Axial velocity profiles (solid lines) in the flow ahead of the flame in different cross sections of the tube (dashed lines).

Tulip flame formation in a half-open tube

In the case of a half-open tube, when the flame is ignited at the left closed end and propagates to the right open end, there are no reflected pressure waves, and the pressure ahead of the flame remains almost constant. As can be seen in Figure 7, the acceleration stage in a half-open tube is almost the same as for flames in closed tubes, but the rate of deceleration in the case of a half-open tube is noticeably lower than in tubes with both closed ends. In tubes, with both ends closed, the rate of the flame deceleration increases due to the increasing pressure. On the contrary, in a half-open tube only the rarefaction waves affect the dynamics of the flow ahead of the flame and, consequently, the formation of the tulip flame.

Figure 12 shows the time evolution of the flame surface area $F_f$, the speed of the combustion wave $S_f$, and the local velocities of the flame front at the central line and near the sidewall, at$y=0.4\mathrm{c}\mathrm{m}$. It can be seen that the flame deceleration stage in a half-open tube lasts almost twice longer than in a tube with both ends closed. The inversion of the flame front and the formation of a tulip-shaped flame also last about two times longer compared to tubes with both ends closed, and the decrease in the flame surface area is monotonous and smoother.

Figure 12. Calculated time evolution of combustion wave velocity $S_f$, the flame surface area $F_f$ (dashed-dotted) and local velocities of the flame front at $y=0$ and near the sidewall at $y=0.4\mathrm{c}\mathrm{m}$ for the half-open tube.

Figure 13(a, b) shows the sequences of calculated Schlieren images and streamlines during the development of a tulip flame for the flame propagating in the half-open tube.

Figure 13. (a, b) Sequences of computed schlieren images and the stream lines for the premixed hydrogen–air flame propagating in a half-open tube.

Figure 14 shows the flow velocities along the center line and at $y=0.25cm$ and at $y=0.4cm$ in the unreacted gas at 1 mm ahead of the flame front and the difference in velocities $\mathrm{\Delta }{U}_{+}=U_{+}(y=0.25cm)-U_{+}(0)$ and $\mathrm{\Delta }{U}_{+}=U_{+}(y=0.4cm)-U_{+}(0)$.

Figure 14. (a, b): Velocities of the flow at 1 mm ahead of the flame front at $y=\hspace{0.17em}0.4\hspace{0.17em}\mathrm{c}\mathrm{m}$ (a) and $y=\hspace{0.17em}0.25\hspace{0.17em}\mathrm{c}\mathrm{m}$ (b) and the corresponding difference between $U_{+}(y=0.4cm),U_{+}(y=0.25cm)$ and $U_{+}(y=0)$; half open tube.

It can be seen that after 1 ms, the flow velocities in the unburned gas closer to the side walls exceed the flow velocity at the center line, and the further away from the center line the greater the difference $\mathrm{\Delta }{U}_{+}$. This trend continues all the way up to the “inner” edge of the boundary layer, where the flow velocity closer to the wall becomes maximum, decreases inside the boundary layer and vanishes at the side wall of the channel. Accordingly, the local velocity of the flame front $U_{fl}=U_f+U_{+}$ becomes minimum at the central line, gradually increases toward the boundary layer where it becomes maximum and decreases inside the boundary layer.

Figure 15 shows the axial velocity profiles in the unburned flow in various cross sections of the tube ahead of the flame front. It can be seen that the boundary layer thickness gradually increases, and when the flame front inverts to form a tulip flame, the boundary layer thickness increases up to approximately 1 mm, so that the location of the tips of the tulip petals is determined by the thickness of the boundary layer and by the tulip-shaped axial velocity profile in the near-field flow ahead of the flame front.

Figure 15. Axial velocity profiles (solid lines) in the flow ahead of the flame in different cross sections of the half-open tube (dashed lines).

Figure 16 shows the dynamics of pressure profiles and velocities at the center line and closer to the boundary layer in the flow ahead of the flame during the formation of the tulip flame. It can be seen that a negative pressure gradient exists only during the acceleration stage, up to 0.65 ms. After 0.65 ms, when the flame surface begins to decrease and the decelerating flame generates rarefaction waves, a reverse (positive) pressure gradient is formed, which persists until the tulip flame is formed. Since the flame deceleration is weaker in the case of a half-open tube, the rarefaction waves created by the decelerating flame are also weaker than in the case of tubes with both ends closed. Therefore, rarefaction waves can reverse the pressure gradient and create a positive pressure gradient in a half-open channel, but their intensity is not sufficient for creating a reverse flow in the unburnt mixture.

Figure 16. Pressure and the flow velocity profiles ahead of the flame during the development of the tulip flame for half open tube.

3D simulations of tulip flame formation in a channel L/D = 6

The flame dynamics and the formation of a tulip flame in a three-dimensional rectangular channel is qualitatively similar to those in a two-dimensional channel. The most significant difference is the higher rates of flame acceleration and deceleration, so that the maximum flame speed achieved in the acceleration stage is almost two times higher in 3D case compare to the 2D case. Accordingly, the duration of the acceleration and deceleration stages is almost 1.5 times shorter in 3D case compare to the 2D case. The greater acceleration of a three-dimensional flame in a channel with no-slip walls compared to a two-dimensional flame is explained by the fact that the acceleration (deceleration) of the flame is determined by the increase (decrease) in the surface area of the flame, which is a line in the two-dimensional case and a surface in the three-dimensional case. The difference in the dynamics of the 2D and 3D flames was first demonstrated and explained by Ivanov et al. (2013) and by Liberman (2021) in connection with modeling the transition from deflagration to detonation in 2D and 3D rectangular channels.

Figure 17 shows the time evolution of the flame surface area $F_f$, the combustion wave velocity $S_f$, and the velocities of the flame front $U_{fl}(y=0)$ and $U_{fl}(y=0.3\mathrm{c}\mathrm{m})$ for the flame propagating in the rectangular channel of length $L=6cm$ and cross section $D\times D=1\mathrm{c}{\mathrm{m}}^2$. It can be seen that, qualitatively the flame dynamics shown in Figure 17 is similar to the flame dynamics shown in Figure 1 for the two-dimensional case. In the two-dimensional case, the surface area of the flame begins to decrease at 0.5 ms, and the deceleration stage begins at 0.55 ms, while in the three-dimensional case for a channel with the same aspect ratio these times are 0.33 ms and 0.38 ms. In particular, there are no oscillations caused by the reflected pressure wave collision with the flame front. The intensity of the rarefaction wave is greater for flames in three-dimensional channels compared to the two-dimensional case. As shown in the Appendix, the LES model somewhat underestimates the effect of pressure waves, which is noticeable for channels with a small aspect ratio, but weaker or even absent for channels with a large aspect ratio and for a semi-open channel.

Figure 17. Time evolution of the flame surface area $F_f$ (dashed-dotted), local velocities of the flame front at $y=0$ and at $y=0.3cm$, and combustion wave velocity $S_f$,for the flame propagating in tube with both ends closed, $L/D=6$, $D\times D=1\mathrm{c}{\mathrm{m}}^2$.

Figure 18(a,b) shows the flow velocities at the cross section $(x,y,z=0)$ at 1 mm ahead of the flame front at the central line $y=z=0$, and at $y=0.15\mathrm{c}\mathrm{m}$ and $y=0.3\mathrm{c}\mathrm{m}$, and the differences $\mathrm{\Delta }{U}_{+}=U_{+}(y=0.15cm)-U_{+}(y=0)$ and $\mathrm{\Delta }{U}_{+}=U_{+}(y=0.3cm)-U_{+}(y=0)$. It is seen that the local velocities of the flame front away from the central line begin to exceed the flame front velocity along the axis ($y=z=0$) after 0.42 ms. After 0.5 ms, the unburned flow velocity closer to the walls exceeds the flow velocity along the tube axis, the convex flame front begins to flatten, and then a tulip-shaped flame develops very quickly. The difference between the flow velocities closer to the wall and the velocity at the axis increases in the direction from the channel axis toward to side walls and becomes maximum at the distance 0.3 cm from the axis, resulting in a tulip-shaped profile of the axial velocity in the unburned mixture close ahead of the flame front. Eventually, this leads to the formation of a tulip-shaped flame, with the tips of the petals forming at a distance of 0.3 cm from the channel axis.

Figure 18. (A, b) the flow velocities at 1mm ahead of the flame front at the cross section $(x,y,z=0)$, at $y=0$, $y=0.15\mathrm{c}\mathrm{m}$, and $y=0.3\mathrm{c}\mathrm{m}$; a) $\mathrm{\Delta }{U}_{+}=U_{+}(y=0.3cm)-U_{+}(y=0)$, b) $\mathrm{\Delta }{U}_{+}=U_{+}(y=0.15cm)-U_{+}(y=0)$.

Figure 19(a,b) shows numerical Schlieren images with the stream lines calculated for selected times for the flame propagating in three-dimensional channel $L=6\mathrm{c}\mathrm{m}$, $D\times D=1\mathrm{c}{\mathrm{m}}^2$

Figure 19. (a, b) Schlieren images and stream lines; premixed hydrogen/air flame propagating in the 3D tube with both closed ends $L=6\mathrm{c}\mathrm{m}$, cross section $1\times 1\mathrm{c}{\mathrm{m}}^2$, shown is cross section $(x,y,z=0)$.

It can be seen in Figure 19 that the convex flame front begins to flatten at 0.5 ms, which is accompanied by the formation of a pair of vortices behind the flame front near the side walls, and the reverse flow of the burned gas along the side walls. The interaction of the reverse flow near the side walls with the remaining “positive” flow near the tube axis creates a large-scale recirculation – a pair of vortices near the left closed end of the tube.

Figure 20 shows the axial velocity profiles in the unburned mixture ahead of the flame in various cross sections of the tube. It is seen that the rarefaction wave causes the increase in the thickness of the boundary layer in the unburned flow close ahead of the flame.

Figure 20. 3D LES simulations: axial velocity profiles in the flow ahead of the flame in different cross sections of the tube.

After 0.5 ms, the axial velocity profile in the near-field zone ahead of the flame takes the shape of a tulip. As a result, the local velocity of the flame front $U_{fl}(y)=U_f+U_{+}(y)$, becomes minimal at the tube axis, $y=0$, gradually increases toward the side walls until it reaches a maximum near the boundary layer at $y=0.3\mathrm{c}\mathrm{m}$, and then decreases in the boundary layer and vanishes at the side wall at $y=0.5\mathrm{c}\mathrm{m}$, and the tulip flame is formed.

Figure 21 shows the time evolution of the pressure gradient and the velocities in the unburned flow at the tube axis ($x,y=0,z=0$) and near the side walls $(x,y=0.3\mathrm{c}\mathrm{m},0)$. It can be seen that the initially negative pressure gradient formed during the acceleration stage at $t<0.4\mathrm{m}\mathrm{s}$ becomes positive at 0.42 ms when the decelerating flame generates a rarefaction wave, and remains positive until 0.58 ms when the flame surface area begins to increase and the previously formed tulip-shaped flame begins to accelerate. From the beginning, at $t<0.42\mathrm{m}\mathrm{s}$, the flow velocity at the central line $U_{+}(x,y=0,0)$ in the near-field zone ahead of the flame front exceeds the flow velocity close $U_{+}(x,y=0.3\mathrm{c}\mathrm{m},0)$near the side wall, but after 0.42 ms the velocity $U_{+}(x,y=0.3\mathrm{c}\mathrm{m},0)$ exceeds the flow velocity at the central line.

Figure 21. Sequences of pressure profiles along the center line and the flow velocities ahead of the flame at $y=0$ and $y=0.3\mathrm{c}\mathrm{m}$; $L=6\mathrm{c}\mathrm{m}$, cross section $1\times 1\mathrm{c}{\mathrm{m}}^2$.

Summary and conclusions

In this paper, we used high resolution direct numerical simulations and Large eddy simulations coupled with a detailed chemical model to study the propagation of a stoichiometric hydrogen-air flame and the tulip flame formation in two-dimensional rectangular channels with different aspect ratios and both closed ends, in a half-open channel, and in a three-dimensional rectangular channel with closed ends. The simulation showed good qualitative agreement between the numerical and experimental results available in the literature, which indicates that the simulation correctly reproduces the main physical processes involved in the formation of the tulip flame. It was found that the formation of a tulip flame is a pure hydrodynamic process, which does not involve the intrinsic instabilities of the flame front in agreement with the conclusion reached by Ponizy, Claverie, and Veyssière (2014) based on an experimental study of the tulip flame formation in a stoichiometric propane-air flame propagating in a cylindrical tube with closed ends.

It was shown that the rarefaction waves generated by the flame during the deceleration stage, when the flame surface area begins to decrease due to the collapse of the lateral parts of the flame skirt on the side walls of the channel, play a key role in the inversion of the flame front and the formation of a tulip flame. The flame front flattening is accompanied by the creation of the first pair of vortices near the side walls and the reversal of the combusted gas flow. Later-formed large scale circulations in the reversed flow of the burned gas are required by the continuity and boundary conditions. However, neither the vortical motion nor the intrinsic flame front instabilities are the physical processes that cause the tulip flame formation.

The interaction of the (“negative”) flow generated by the rarefaction wave with the positive unburned gas flow generated during the acceleration stage leads to a significant reduction in the velocity of the unburned gas and, consequently, to a significant increase in the thickness of the boundary layer in the near-field zone ahead of the flame. This leads to the formation of a tulip-shaped form of the axial velocity profile ahead of the flame, which in turn leads to the inversion of the flame front and the formation of a tulip flame, because the local velocity of any point on the flame front represents the sum of the laminar flame speed relative to the fresh mixture and the flow speed ahead of that point. It can be seen from the mechanism of tulip-shaped flame formation described above that the tips of the tulip-shaped flame petals are located at a distance from the side walls equal to the thickness of the boundary layer in the near region ahead of the flame, and the height of the tulip-shaped flame petal is determined by the profile of the axial velocity of unburned gas in the immediate vicinity ahead of the flame front.

The formation of the inverse pressure gradient is a natural effect of the rarefaction wave since the front of a rarefaction wave propagates forward with the speed of sound while the gas behind the head of the rarefaction wave propagates in the opposite (negative) direction, therefore the pressure gradient in the rarefaction wave is positive, contrary to the negative pressure gradient, which was formed in the flow ahead of the flame during the accelerating stage.

The 3D LES modeling of tulip flame formation is qualitatively similar to the 2D DNS results. The main difference between the two-dimensional and three-dimensional cases is the significantly shorter characteristic times of the acceleration and deceleration stages in the three-dimensional case compared to the two-dimensional case. This means that two-dimensional modeling of the tulip flame formation can provide only qualitative, but not quantitative agreement with the experimental study. Three-dimensional modeling of the tulip flame formation in highly reactive (hydrogen/air) and slow-reactive (methane/air) mixtures and comparison with the corresponding experimental studies will be published in a separate paper.

Acknowledgements

This work was supported by National Natural Science Foundation of China under grant 11732003 (C.Q and C.W.). M.L. would like to thank Axel Brandenburg the kind assistance to attend ICDERS. The helpful discussions with Grisha Sivashinsky and Paul Clavin are greatly acknowledged. No specific funds were received for this particular work; research of Nordic Institute for Theoretical Physics (NORDITA) is partially supported by Nordforsk. (M.L.)

Disclosure statement

No potential conflict of interest was reported by the authors.

Notes

¹ Formally, this is true for a piston that starts moving with negative constant acceleration in a gas at rest (a piston that begins to move out from a tube filled with gas at rest, Landau and Lifshitz 1989). However, this is correct at least qualitatively in our case also.

Appendix A.

Comparison of DNS and LES

Since three-dimensional simulation of realistic size problems is very expensive and practically impossible even with powerful parallel computers, it has become common to use Large Eddy Simulation (LES) to model such kind problems. Although it is generally believed that DNS, used with a sufficiently high resolution and detailed chemical model, correctly describes combustion processes, the question arises as to whether the different numerical schemes provide an accurate model of the problem under consideration. To answer this question, we compare results of DNS obtained for the flame propagating in a two-dimensional channel $L/D=6$ in Sec.4.1, with results obtained using LES for the same two-dimensional problem. Figure A1 shows the time evolution of the flame skirt surface area $F_f$, local velocities of the flame front along the center line $U_{fl}(y=0)$ and near the side wall $U_{fl}(y=0.4\mathrm{c}\mathrm{m})$ and the combustion wave velocity $S_f$ obtained using large eddy simulations with the mesh size 50μm for the same conditions as in Figure 1.

Figure A1. Time evolution of the flame surface area $F_f$ (dashed-dotted), local velocities of the flame front at $y=0$ and at $y=0.4cm$ and combustion wave velocity $S_f$; tube $L/D=6$ , obtained using large eddy simulations.

In general, the time evolution of the flame surface area and the velocities of the flame front shown in Figure A1 is qualitatively and quantitively very similar to the results of 2D DNS in Figure 1. However, since in the LES method used an artificially thickened flame front, the flame front is significantly lesser sensitive to the pressure waves so that the profiles of flame front velocity are smoother and the intensity of the rarefaction wave, the magnitude of the negative acceleration is somewhat smaller during the tulip flame formation.