5.1 Qualitative comparison of the flow evolution with experiments

A visual comparison of a flow field, obtained through resolved LES and the corresponding experiment of § 2, is shown in figures 17 and 18 for an instant where a triple-point collision occurs for a detonation whose cell size is comparable to the channel height. For the LES, the density evolution is shown for various instants of the collision process for $C_{\unicode[STIX]{x1D705}}=6.7$ . This value of $C_{\unicode[STIX]{x1D705}}$ was found to produce cell sizes and overall qualitative behaviour comparable to experimental observation. Here, the channel height is kept at $H=20$ . For the experiment, schlieren images of each corresponding instant are obtained from high-speed photography and are also shown in figures 17 and 18. For comparison, many features are noted (1–7). In figure 17, feature 1 is a Mach shock with a turbulent reaction zone that follows closely. In both the simulation and the experiment, the reaction zone is slightly decoupled from the shock owing to local unsteadiness of the wave front. This can be seen by the thickening of the shocked unburned gas region as time evolves. In the LES, however, this wave appears to be more irregular, containing what appears to be a smaller and less pronounced cell structure (cells within cells). It should be noted that experimental observations, of § 2, report a very stochastic nature of the propagating wave. In this sense, random bifurcations on wave fronts and changes in cell patterns are expected. In some cases much smaller cells are observed experimentally and numerically for the same quiescent conditions and numerical parameters; see figure 19. Feature 2 in figure 17 is a completely decoupled incident shock reaction zone, as can be seen by the large region of dense unburned gas. As the two triple points approach each other, transverse shock waves compress this gas further. At the focal point of the collision, owing to the increased pressure and temperature, the unburned pocket (2) ignites very rapidly. This ‘explosion’, denoted by feature 3, drives transverse shock waves and a strong Mach shock forward coupled with very rapid burning. Feature 4 is a pocket of unburned gas that is not ignited by compression associated with the original Mach shocks (feature 1). Instead, these pockets mix and burn, on their surfaces, along turbulent shear layers. Owing to shock compression associated with the passage of transverse shock waves, these shear layers experience increased burning rates through turbulent mixing of the burned and unburned gases via Richtmyer–Meshkov instability. These transverse shocks, labelled as feature 5 in figure 18, originate from feature 3. Feature 6 in figure 18 is a new pocket of unburned gas that forms along slip lines associated with the collision process. In fact, the formation of such pockets, in this manner, is typical of irregular detonation propagation (Subbotin Reference Subbotin1975; Oran et al. Reference Oran, Young, Boris and Picone1982; Austin Reference Austin2003; Radulescu et al. Reference Radulescu, Sharpe, Lee, Kiyanda, Higgins and Hanson2005, Reference Radulescu, Sharpe, Law and Lee2007; Kiyanda & Higgins Reference Kiyanda and Higgins2013). Finally, feature 7 in figure 18 is a bifurcation observed on the Mach shock of feature 3, owing to the unstable nature of methane–air detonations. This bifurcation also occurs in the experiment but is less pronounced. Although features 1–7 are captured by the LES, a key difference is observed in the burning rate behind the Mach shock originating from feature 3. In the simulation, the size of this region grows slightly faster compared to experiment, owing to the lower detonation velocity reported in § 2.3.

Figure 20. (a) Numerical density evolution for $C_{\unicode[STIX]{x1D705}}=6.7$ and (b) the corresponding experimental schlieren photographs (Kiyanda & Higgins Reference Kiyanda and Higgins2013) of an irregular detonation propagation in a channel. The time between images is $\unicode[STIX]{x0394}t=0.37$ ( $10.0~\unicode[STIX]{x03BC}\text{s}$ ). The remaining sequence of images is continued in figure 21.

Figure 21. Continuation of figure 20. (a) Numerical density evolution for $C_{\unicode[STIX]{x1D705}}=6.7$ and (b) the corresponding experimental schlieren photographs (Kiyanda & Higgins Reference Kiyanda and Higgins2013) of an irregular detonation propagation in a channel. The time between images is $\unicode[STIX]{x0394}t=0.37$ ( $10.0~\unicode[STIX]{x03BC}\text{s}$ ).

To further compare the performance of the LES with available experimental data (Kiyanda & Higgins Reference Kiyanda and Higgins2013), the simulation was repeated with a channel height of $H=10$ . Both the density evolution from the numerical simulation and the corresponding instants obtained in the experiment through schlieren photography are shown in figures 20 and 21. In this case, all model parameters, including $C_{\unicode[STIX]{x1D705}}$ , are consistent with the previous simulation, except that the reduced channel height permits only half a cell to form. Thus the cell pattern is effectively influenced by the channel walls, or mode-locked into a cell size which corresponds roughly to twice the channel height. In this case, only a single transverse wave is observed. Again, for qualitative comparison, many features (1–7) are noted in figures 20 and 21. Feature 1 in figure 20 is the incident shock wave and decoupled reaction zone, as observed by the region of very dense, shocked and unburned fuel following the shock. From the top wall of the channel, feature 2 is a Mach stem, which initially results from a transverse wave collision with the wall. As observed in figure 17, this Mach stem travels faster than the incident wave, owing to the wave collision process. Feature 3 in figure 20 is a pocket of dense unburned fuel in the wake of the wave, which is consumed by a surface turbulent flame. Furthermore, consumption rates would be enhanced by Richtmyer–Meshkov instability associated with the transverse shock wave, feature 4. In fact, this region burns out only slightly faster than the observed experimental rate. In the experiment, the pocket burns up completely within seven frames ( $70~\unicode[STIX]{x03BC}\text{s}$ ), while the simulation burns up completely within five frames ( $50~\unicode[STIX]{x03BC}\text{s}$ ). Using the method described in § 2.4, the turbulent flame speed of this experiment is found to be ${\hat{S}}_{t}/{\hat{S}}_{L}\approx 6.6$ ( ${\hat{S}}_{t}\approx 110~\text{m}~\text{s}^{-1}$ ). This is comparable to the experimental flame speed found in the current study, § 2, where ${\hat{S}}_{t}/{\hat{S}}_{L}\approx 7.3$ ( ${\hat{S}}_{t}\approx 120~\text{m}~\text{s}^{-1}$ ). Feature 5 in figure 20 is a hotspot that forms behind a bifurcation of the Mach shock and is only clearly visible in the numerical simulation. This bifurcation of the Mach stem near the triple point was also observed in the previous comparison to experiment (figure 17). As the Mach stem, feature 1, continues to evolve, the combustion zone becomes further decoupled from the shock, as observed by an increased thickness of dense unburned fuel behind the shock. Furthermore, a new pocket of dense unburned fuel forms in the wake of the triple-point path, denoted by feature 6 in figure 21. This pocket is initially consumed from its edges via turbulent flame. The passage of the reflected transverse shock wave, feature 7, further contributes to turbulent mixing through Richtmyer–Meshkov instability.

Finally, for the $H=10$ simulation, figure 22 shows how the instantaneous reaction rate obtained numerically compares with the corresponding experimental photographs of Kiyanda & Higgins (Reference Kiyanda and Higgins2013), which recorded the self-luminosity signal. In figure 22, the instantaneous rate of reaction, $\dot{\unicode[STIX]{x1D714}}$ , is superimposed onto two consecutive density evolution plots to show where chemical reactions are occurring, and the intensity. In both the experiment and the simulation, chemical reactions appear to be very intense behind the Mach shock, where short ignition delays are expected due to the overdriven shock. Also, combustion of the unburned gas pockets also becomes intensified with passage of the transverse shock wave. Finally, for the unburned pockets of gas, the chemical reactions occur predominantly on the surfaces, and not uniformly throughout. This suggests that turbulent mixing is the dominant mechanism through which the unburned pockets are consumed. It is therefore of particular interest to quantify the turbulent mixing rates and determine their impact on burnout rates of the pockets. Furthermore, it is of interest to determine the overall impact such burning has on the observed structure and cell pattern of the detonation wave.

Figure 22. (a) Numerical density evolution with superimposed chemical reaction rate ( $\overline{\dot{\unicode[STIX]{x1D714}}}$ ), shown in red, and (b) the corresponding experimental self-luminous images (Kiyanda & Higgins Reference Kiyanda and Higgins2013). The time between images is $\unicode[STIX]{x0394}t=0.74$ ( $20.0~\unicode[STIX]{x03BC}\text{s}$ ).

5.2 Quantitative analysis: turbulent flame speeds

For the detonation simulations presented in § 5.1, where $C_{\unicode[STIX]{x1D705}}=6.7$ for both $H=10$ and $H=20$ , it is of particular interest to estimate and compare the turbulent flame speed at which the pockets of unreacted gas, labelled as feature 4 in figure 17 and as feature 3 in figure 20, burn up. In contrast to the experiments, the turbulent flame speed of the pockets in the simulations is readily available by considering the rates at which reactant mass is consumed within the one-dimensional ‘samples’ of the CLEM subgrid. For the pockets, the local turbulent flame speeds are estimated from

(5.1) $$\begin{eqnarray}S_{t,local}=\frac{{\dot{m}}}{\unicode[STIX]{x1D70C}_{u}A_{c}},\end{eqnarray}$$

where ${\dot{m}}$ is the instantaneous rate of reactant mass consumption, $\unicode[STIX]{x1D70C}_{u}$ is the density of the upstream unburned reactant, located to within $0.5\unicode[STIX]{x1D6E5}_{1/2}$ , and $A_{c}$ is the cross-sectional area through which the flame propagates against the unburned reactant within each ‘sample’. In this case, $A_{c}$ is known and constant for all CLEM domains. In this approach, mass is conserved since the burning rate is balanced by the mass flow rate of unburned reactant. Furthermore, this method was previously applied to calculate turbulent flame speeds in a stand-alone one-dimensional CLEM subgrid formulation (Maxwell et al. Reference Maxwell, Falle, Sharpe and Radulescu2015). Finally, the global turbulent flame speed is obtained by ensemble averaging the flame speeds on the pocket surface in both space and time. Thus

(5.2) $$\begin{eqnarray}S_{t}=\mathop{\sum }_{1}^{n}S_{t,local}/n,\end{eqnarray}$$

where $n$ is the number of samples acquired on the pocket surface.

Using this method, the turbulent flame speeds of the unburned pockets are found to be $S_{t}/S_{L}=3.74$ ( ${\hat{S}}_{t}=61.3~\text{m}~\text{s}^{-1}$ ) and $S_{t}/S_{L}=3.63$ ( ${\hat{S}}_{t}=59.5~\text{m}~\text{s}^{-1}$ ) for the $H=20$ and $H=10$ simulations, respectively. These values differ, by a factor of 2, from those found experimentally, which were $S_{t}/S_{L}\approx 7.3$ ( ${\hat{S}}_{t}=120~\text{m}~\text{s}^{-1}$ ) and $S_{t}/S_{L}\approx 6.6$ ( ${\hat{S}}_{t}\approx 110~\text{m}~\text{s}^{-1}$ ) for $H=20$ and $H=10$ , respectively. The difference in simulation values from the experimental ones are noted in the difficulties of evaluating the true volume and surface area solely from schlieren images in the experiments. The simulation and experiment, however, both predict turbulent flame speeds of roughly half an order to one order of magnitude of the turbulent flame speed.

To gain further insight into the effect of local turbulent mixing rates on turbulent flame speeds, it is possible to obtain average turbulent flame speeds for given ranges of turbulence intensity values, $u^{\prime }$ , since the subgrid kinetic energy is known. In fact, the turbulence intensity is simply

(5.3) $$\begin{eqnarray}u^{\prime }=\sqrt{{\textstyle \frac{2}{3}}k^{sgs}}.\end{eqnarray}$$

Figure 23 shows the average turbulent flame speeds, $S_{t}/S_{L}$ , in the wake of the simulated detonations according to their turbulence intensities, $u^{\prime }/S_{L}$ . For each value of $u^{\prime }/S_{L}$ , the flame speeds are averaged within intervals of $u^{\prime }/S_{L}\pm 0.5$ . Also shown in figure 23 are fan-stirred bomb experiment measurements of Abdel-Gayed, Al-Khishali & Bradley (Reference Abdel-Gayed, Al-Khishali and Bradley1984) for stoichiometric methane–air mixtures at atmospheric pressures. Clearly, the simulation exhibits larger average turbulent flame speeds in regions of higher turbulence intensity. In fact, much like the experiments of Abdel-Gayed et al. (Reference Abdel-Gayed, Al-Khishali and Bradley1984), $S_{t}/S_{L}$ appears to increase proportionally with $u^{\prime }/S_{L}$ . Turbulence intensities were not recorded above the values shown in figure 23, as such high-velocity fluctuations have not had sufficient time to develop on the pocket surfaces. For the turbulence intensities observed, the results are in good agreement with the experiments.

Figure 23. Average turbulent flame speeds ( $S_{t}/S_{L}$ ) versus turbulent intensities ( $u^{\prime }/S_{L}$ ), obtained from simulation and compared to experiments of Abdel-Gayed et al. (Reference Abdel-Gayed, Al-Khishali and Bradley1984).

5.3 Quantitative analysis: the reaction zone and hydrodynamic thickness

In order to determine the effect that turbulent mixing intensity has on the reaction zone thickness, the average thickness of the structure ( $\unicode[STIX]{x1D6E5}_{R}$ ) has been measured for various $C_{\unicode[STIX]{x1D705}}$ values for the domain height of $H=20$ . This was achieved by following the procedure to generate the Favre-averaged reactant profiles in figure 16. Figure 24 thus shows the reaction zone thickness ( $\unicode[STIX]{x1D6E5}_{R}$ ), defined here as the distance from the shock wave, $x=x_{s}$ , to where ${\tilde{Y}}(x)<2\,\%$ , as a function of $C_{\unicode[STIX]{x1D705}}$ . Clearly, there is an increasing trend in thickness with $C_{\unicode[STIX]{x1D705}}$ . In fact, the thickness appears to increase linearly with $C_{\unicode[STIX]{x1D705}}$ for the range of values simulated here. A line of best fit (obtained from linear regression of each data point), $\unicode[STIX]{x1D6E5}_{R}=1.19C_{\unicode[STIX]{x1D705}}+4.49$ , is also indicated in figure 24. Also shown are error bars indicating the standard deviation of thickness for each value of  $C_{\unicode[STIX]{x1D705}}$ .

Figure 24. Effect of turbulence intensity (through $C_{\unicode[STIX]{x1D705}}$ ) on the mean reaction zone thickness ( $\unicode[STIX]{x1D6E5}_{R}$ ) for the CLEM-LES (all with resolution $b=1/32$ and $N=16$ elements per LES cell). The error bars indicate the standard deviation in $\unicode[STIX]{x1D6E5}_{R}$ for each simulation. Note that $\unicode[STIX]{x1D6E5}_{R}$ is normalized by $\hat{\unicode[STIX]{x1D6E5}}_{1/2}$ .

Figure 25. Schlieren photograph (this study), which shows that pockets of unburned gas can take up to ${\sim}9\hat{\unicode[STIX]{x1D6E5}}_{1/2}$ behind the leading shock to burn up.

To quantitatively compare the CLEM-LES to experimental observations in terms of the reaction zone thickness, or distance it takes for the pockets of unburned gas to burn, results from § 4.3 are considered. In figure 24, it was observed that the average reaction zone thickness captured by all values of $C_{\unicode[STIX]{x1D705}}$ , for $H=20$ , was in the range $5<\unicode[STIX]{x1D6E5}_{R}<17$ . More specifically, for $C_{\unicode[STIX]{x1D705}}=6.7$ and $C_{\unicode[STIX]{x1D705}}=10.0$ , the average thickness was found to be $\unicode[STIX]{x1D6E5}_{R}=13.0$ and $\unicode[STIX]{x1D6E5}_{R}=16.2$ , respectively. These computed average thicknesses are much higher than previous reports of burnout distances of ${\sim}(4{-}6)\hat{\unicode[STIX]{x1D6E5}}_{1/2}$ (Radulescu et al. Reference Radulescu, Sharpe, Law and Lee2007), which were measured from Euler simulations with significant numerical viscosity. Upon close inspection, however, of the experimental images shown in § 2, burnout distances were estimated to be as high as ${\sim}9\hat{\unicode[STIX]{x1D6E5}}_{1/2}$ , as can be seen in figure 25. Thus, compared to the experiment in figure 25, $C_{\unicode[STIX]{x1D705}}=6.7$ yields an error of ${\sim}44\,\%$ , while $C_{\unicode[STIX]{x1D705}}=10.0$ yields an error of ${\sim}80\,\%$ . Despite these large errors, however, it should be noted that the simulations for $C_{\unicode[STIX]{x1D705}}=6.7$ and $C_{\unicode[STIX]{x1D705}}=10.0$ have standard deviations of $\unicode[STIX]{x1D6E5}_{R}\pm 6.4$ and $\unicode[STIX]{x1D6E5}_{R}\pm 5.7$ , respectively. Thus, a large range of possible expected thicknesses lie anywhere from $6.6<\unicode[STIX]{x1D6E5}_{R}<22$ . Furthermore, the current experiments, which are stochastic in nature, have not quantified the probability and variance of the reaction zone thickness that are possible during propagation.

Finally, the reaction zone thickness has been well correlated with the detonation cell size, where chemical reactions are completed roughly within one cell cycle (Lee Reference Lee1984). From figure 15, the cell size $\unicode[STIX]{x1D706}\approx 10$ when $C_{\unicode[STIX]{x1D705}}=6.7$ . This compares well with the average reaction zone thickness, $\unicode[STIX]{x1D6E5}_{R}=13.0$ . Furthermore, the entire hydrodynamic thickness of the wave ( $\unicode[STIX]{x1D6E5}_{H}$ ), or distance from the leading shock to the trailing sonic plane (or CJ state), is also well correlated with the cell size: $\unicode[STIX]{x1D6E5}_{H}=6.5\unicode[STIX]{x1D706}$ (Lee Reference Lee1984). In figure 26, the ensemble space- and time-averaged pressure profile for $\bar{p}(x)$ is presented for the $C_{\unicode[STIX]{x1D705}}=6.7$ and $H=20$ case and compared to the ZND solution for $p(x)$ . For this case, it takes approximately ${\sim}50\unicode[STIX]{x1D6E5}_{1/2}$ for the pressure to reach the CJ solution state, which corresponds to the sonic plane in the wake of the wave. This is consistent with findings in Lee & Radulescu (Reference Lee and Radulescu2005) for highly unstable detonations involving methane at similar conditions. Furthermore, since $\unicode[STIX]{x1D706}\approx 10$ for $C_{\unicode[STIX]{x1D705}}=6.7$ , then $\unicode[STIX]{x1D6E5}_{H}\approx 5\unicode[STIX]{x1D706}$ . This result is of the same order as that reported by Lee (Reference Lee1984) from previous experimental correlations. To summarize the findings here, increasing the value of $C_{\unicode[STIX]{x1D705}}$ was found to generate a larger reaction zone thickness. This, in turn, is believed to lead to an increased hydrodynamic thickness and cell size, due to their dependence on the reaction zone thickness.

Figure 26. Ensemble space- and time-averaged pressure profile ( $\bar{p}(x)$ ) for the $C_{\unicode[STIX]{x1D705}}=6.7$ and $H=20$ simulation and compared to the ZND solution for $p(x)$ . Note that distances $x$ and $x_{s}$ are normalized by $\hat{\unicode[STIX]{x1D6E5}}_{1/2}$ . The pressure $p$ is normalized by $\unicode[STIX]{x1D6FE}p_{o}$ .

5.4 Qualitative comparison of cell patterns with experiments

To further compare the CLEM-LES with experiment, in terms of cell patterns and irregularity, numerically obtained soot foils for both $H=10$ and $H=20$ at $C_{\unicode[STIX]{x1D705}}=6.7$ were compared to those obtained experimentally for $\text{CH}_{4}+2\text{O}_{2}$ at $\hat{p}_{o}=3.5$  kPa, courtesy of Radulescu et al. (Reference Radulescu, Sharpe, Lee, Kiyanda, Higgins and Hanson2005), and $\text{CH}_{4}+2\text{O}_{2}+0.2\text{air}$ at $\hat{p}_{o}=11$  kPa (Austin Reference Austin2003). For $H=20$ , the numerical soot foil, presented in figure 27, corresponds exactly to the soot foil image of figure 15 for $C_{\unicode[STIX]{x1D705}}=6.7$ . Although the experimental soot foil, also shown in figure 27, was obtained for a smaller channel height (127 mm) and also for a much higher pressure, the two are comparable in terms of scaling through the half reaction length. For the experiment (Austin Reference Austin2003), $\hat{\unicode[STIX]{x1D6E5}}_{1/2}=5.1$  mm. This value was obtained in this study using the Cantera libraries (Goodwin et al. Reference Goodwin, Moffat and Speth2016) and the GRI-3.0 detailed kinetic mechanism (Smith et al. Reference Smith, Golden, Frenklach, Moriarty, Eiteneer, Goldenberg, Bowman, Hanson, Song and Gardiner2016), consistent with the procedure in Maxwell (Reference Maxwell2016). Thus, the experimental soot foil was estimated to be ${\sim}25\hat{\unicode[STIX]{x1D6E5}}_{1/2}$ high. This is comparable to the simulation, whose domain is $20\hat{\unicode[STIX]{x1D6E5}}_{1/2}$ in height. In the portion of the numerical soot foil presented in figure 27, the overall cell size and pattern match well to the experiment. In the experiment, however, a very distinct substructure was observed. This was observed by the presence of much smaller cells within the larger and more predominant cell structure. Although such a fine substructure is not so obvious in the simulation, there are still some streaks on the numerical soot foil, which are indicative of triple-point paths associated with cell bifurcations. The lack of detailed substructure in the simulation could be an artifact of the supergrid resolution, differences in mixture properties, or most likely associated with the limitation of the one-step model in capturing the very short-duration energy release of methane combustion and the associated instabilities that exist on that finer small scale. Furthermore, it is likely that, at the resolution presented here, 32 grids per $\hat{\unicode[STIX]{x1D6E5}}_{1/2}$ , the fine details associated with the expected substructure have effectively been filtered out. It is noted that higher-resolution simulations were not conducted for $C_{\unicode[STIX]{x1D705}}=6.7$ . Despite this lack of fine-scale detail for the substructure, the overall cell behaviour and irregularity are captured well by the simulation.

Figure 27. (a) Numerical soot foil obtained for the CLEM-LES with $C_{\unicode[STIX]{x1D705}}=6.7$ and $H=20$ compared with (b) an experimental soot foil for $\text{CH}_{4}+2\text{O}_{2}+0.2\text{air}$ at $\hat{p}_{o}=11~\text{kPa}$ (Austin Reference Austin2003). Note that the channel height for the experiment is ${\sim}25\hat{\unicode[STIX]{x1D6E5}}_{1/2}$ (127 mm), while the simulation is $20\hat{\unicode[STIX]{x1D6E5}}_{1/2}$ (203 mm).

For $H=10$ , figure 28 shows numerical soot foils obtained for two different portions of the numerical domain and compares them to soot foils of experiments conducted in Radulescu et al. (Reference Radulescu, Sharpe, Lee, Kiyanda, Higgins and Hanson2005). Here, the experimental soot foils were obtained using the same experimental apparatus as described in Kiyanda & Higgins (Reference Kiyanda and Higgins2013). For both experimental soot foils, the numerical simulation captures well the cell size behaviour and irregularity. In general the channel height allows only half a cell to form. The cells, however, occasionally bifurcate, giving rise to the formation of cells of the order of the channel height. This behaviour can be observed in the bottom frames of figure 28. Finally, it is noted that a prominent substructure was not observed experimentally or numerically for $H=10$ , as was observed experimentally in figure 27. This could be attributed to differences in the mixture or conditions (i.e. the pressure) at which tests were performed.

Figure 28. Numerical soot foils obtained for the CLEM-LES with $C_{\unicode[STIX]{x1D705}}=6.7$ and $H=10$ compared with two experimental soot foils for $\text{CH}_{4}+2\text{O}_{2}$ at $\hat{p}_{o}=3.5~\text{kPa}$ , courtesy of M. Radulescu. Note that the channel height for the experiment is also $10\hat{\unicode[STIX]{x1D6E5}}_{1/2}$ .

5.5 Quantitative analysis: velocity of the wave

To gain insight into why increasing  $C_{\unicode[STIX]{x1D705}}$  has the effect of generating larger hydrodynamic structures, the subgrid kinetic energy associated with random velocity fluctuations ( $k^{sgs}$ ) is Favre- and ensemble-averaged for each $C_{\unicode[STIX]{x1D705}}$ case using the same procedure as was done for the reactant profiles, i.e. equations (4.2)–(4.4). The averaged profiles obtained for $k^{sgs}$ at various values of $C_{\unicode[STIX]{x1D705}}$ are thus shown in figure 29. The principal observation made here is that, as $C_{\unicode[STIX]{x1D705}}$ is increased, more subgrid kinetic energy is generated. This is not surprising since $k^{sgs}$ is a function of $C_{\unicode[STIX]{x1D705}}$ , as demonstrated in Maxwell (Reference Maxwell2016). This increase in $k^{sgs}$ with $C_{\unicode[STIX]{x1D705}}$ thus generates more frequent stirring events on the CLEM subgrid, and therefore faster burning rates on pocket surfaces. Despite this faster mixing and burning, larger pockets of unburned gas are able to form in the wake. This suggests that an overall velocity deficit develops on the wave front, allowing for increased ignition delays.

Figure 29. Effect of $C_{\unicode[STIX]{x1D705}}$ on the averaged subgrid kinetic energy ( $k^{sgs}$ ) profile for the CLEM-LES. Note that distances $x$ and $x_{s}$ are normalized by $\hat{\unicode[STIX]{x1D6E5}}_{1/2}$ .

To quantify the statistical distribution of velocities experienced by the wave front, and how such a distribution was affected by $C_{\unicode[STIX]{x1D705}}$ , p.d.f.s were constructed from the numerical simulations and are shown in figures 30 and 31 for $H=20$ and $H=10$ , respectively. Also shown in these figures are the experimental p.d.f.s previously shown in figure 6. To construct the p.d.f.s from the numerical simulations, velocity measurements were taken on the top and bottom walls for a total of 500 data points in each simulation. The p.d.f.s were then evaluated at $D\pm 0.5$ intervals and normalized by $D_{avg}$ accordingly. For the CLEM-LES, conducted at various $C_{\unicode[STIX]{x1D705}}$ values, and also an Euler simulation at the same supergrid resolution ( $b=1/32$ ), the average propagation speed is always within 1 % error of the CJ value, where $D_{CJ}=6.30$ ( $2240~\text{m}~\text{s}^{-1}$ ).

Figure 30. The p.d.f. of a detonation wave, in a channel with $H=20$ , having a certain velocity ( $D/D_{avg}$ ) at any given moment and location. Also shown is a p.d.f. compiled from experiments (this study). Note that $D_{avg}=5.19$ ( $1850~\text{m}~\text{s}^{-1}$ ) for the experiment and $D_{avg}=6.30$ ( $2240~\text{m}~\text{s}^{-1}$ ) for the simulations.

Figure 31. The p.d.f. of a detonation wave, in a channel with $H=10$ , having a certain velocity ( $D/D_{avg}$ ) at any given moment and location. Also shown is a p.d.f. compiled from Kiyanda & Higgins (Reference Kiyanda and Higgins2013). Note that $D_{avg}=5.53$ ( $1970~\text{m}~\text{s}^{-1}$ ) for the experiment and $D_{avg}=6.30$ ( $2240~\text{m}~\text{s}^{-1}$ ) for the simulations.

Clearly, from figures 30 and 31, as $C_{\unicode[STIX]{x1D705}}$ increases, the likelihood or probability of the detonation wave to travel at speeds below the average CJ value also increases. This is especially true for $C_{\unicode[STIX]{x1D705}}\geqslant 6.7$ . At these high values of $C_{\unicode[STIX]{x1D705}}$ , the detonation tends to favour a greater chance of having wave speeds below the average propagation speed value, i.e. when $(D/D_{avg})<1$ . Conversely, at lower values of $C_{\unicode[STIX]{x1D705}}$ , and the Euler simulation, velocities above the average propagation speed are favoured. This would lead to much shorter ignition delays, thus explaining the relatively quick burning times associated with the observed shortened reaction zone thickness, previously shown in figure 24. For higher $C_{\unicode[STIX]{x1D705}}$ values, since the wave spends more time at below CJ velocities, most of the unburned gas that is shocked by the wave front has a much longer ignition delay compared to the ZND model. For this reason, unburned pockets of gas are able to form in the wake, which eventually burn up through turbulent mixing. This favouring of velocities below the average propagation speed thus lengthens the overall hydrodynamic structure of the wave. Consequently, the cells also increase in size and become more irregular in appearance, as observed in the cell patterns of figure 15. Since higher $C_{\unicode[STIX]{x1D705}}$ values inherently generate more random fluctuations, the cell irregularity would also be expected to increase.

In general, the p.d.f.s collected for $C_{\unicode[STIX]{x1D705}}\geqslant 6.7$ compared well to experiments, in terms of the most probable velocity expected on the wave front. The experiments and the simulations, for $C_{\unicode[STIX]{x1D705}}\geqslant 6.7$ , both show that the most probable velocity on the wave front, at any given instant or location, is actually below the average propagation value, $(D/D_{avg})<1$ . The simulations for $C_{\unicode[STIX]{x1D705}}\geqslant 6.7$ also reproduce well the decaying behaviour of the wave velocity, where the probability of the wave speed exhibits a power-law dependence on wave speeds above the favoured value. For $H=20$ , although the numerical simulations for $C_{\unicode[STIX]{x1D705}}\geqslant 6.7$ do not collapse onto this correlation, the same decaying trend is observed. For $H=20$ , the most probable wave speed favoured, experimentally, is $(D/D_{avg})=0.93\pm 0.04$ with a peak p.d.f. value of 3.3. For the simulation at $C_{\unicode[STIX]{x1D705}}=6.7$ , the most probable velocity was $(D/D_{avg})=0.96\pm 0.04$ with a peak p.d.f. value of 1.74. This favoured velocity value has 8.6 % error compared to experiment. For $C_{\unicode[STIX]{x1D705}}=10.0$ , the most probable velocity was $(D/D_{avg})=0.85\pm 0.04$ with a peak p.d.f. value of 2.25. This agrees well for the favoured velocity value, with only 3.2 % error. Despite this, however, the range of expected velocities on the wave front is much larger in the simulations compared to experiment. These differences in peak p.d.f. values and expected ranges of velocities are believed to be influenced by several factors. First, the numerical simulations do not account for energy losses that are present in the experiments. Therefore, the detonation, in the simulations, has higher velocities. Furthermore, cell patterns for a channel height of $H=20$ are very stochastic and irregular. The experimental p.d.f. in figure 30 was compiled for only four different experiments. It is likely that many more experiments are required in order to capture the events in the tails of the p.d.f., beyond the velocity limits currently obtained. Finally, velocity measurements are captured experimentally every $11.53~\unicode[STIX]{x03BC}\text{s}$ , which may filter out any high-speed and short-lived velocity fluctuations, which may occur beyond the current p.d.f. limits. It is likely that, in reality, the experimental p.d.f. should have a larger range with a smaller peak p.d.f. value. Despite the differences in p.d.f. values between experiment and simulation, it should be noted that $C_{\unicode[STIX]{x1D705}}\geqslant 6.7$ does much better at capturing the peak value location $(D/D_{avg})$ compared to Euler or low-value $C_{\unicode[STIX]{x1D705}}$ simulations. This location, when $(D/D_{avg})<1$ , is the key feature that influences pocket formation, and thus increased hydrodynamic thickness, owing to increased ignition delays at most locations on the wave front.

For $H=10$ , much better agreement with experiment is observed when $C_{\unicode[STIX]{x1D705}}=6.7$ . The most probable wave speed favoured, experimentally (Kiyanda & Higgins Reference Kiyanda and Higgins2013), is $(D/D_{avg})=0.814\pm 0.005$ with a peak p.d.f. value of 6.15. In this case, the simulation at $C_{\unicode[STIX]{x1D705}}=6.7$ has a favoured velocity of $(D/D_{avg})=0.84\pm 0.04$ with a peak p.d.f. value of 2.55. This also agrees well for the favoured velocity value, with only 1.0 % error. This better agreement can probably be attributed to the fact that a domain height of $H=10$ effectively mode-locks the detonation in such a way that much less variation in cell size was observed. This was seen in the soot foils presented in figure 28. Clearly, the analysis conducted here supports the need to calibrate $C_{\unicode[STIX]{x1D705}}$ in order to ensure the correct velocity probability trend and qualitative features are obtained numerically.

Figure 32. Numerical soot foils obtained for the 3D CLEM-LES for $C_{\unicode[STIX]{x1D705}}=6.7$ and $H=10$ .

5.6 Three-dimensional effects and the validity of the two-dimensional LES approach

Finally, in order to further validate the two-dimensional (2D) numerical results obtained throughout this work, a thin-channel three-dimensional (3D) simulation has been conducted for the $H=10$ case, with $C_{\unicode[STIX]{x1D705}}=6.7$ and a domain width of $2\hat{\unicode[STIX]{x1D6E5}}_{1/2}$ . This corresponds roughly to the experimental set-up of Kiyanda & Higgins (Reference Kiyanda and Higgins2013). The grid resolution used for the 3D simulation was $b=1/16$ with $N=16$ subgrid elements within each cell. Figure 32 shows a portion of soot foils which were obtained on the side and top walls of the domain. More specifically, these are obtained at $z=0$ in the $x{-}y$ plane and $y=10$ in the $x{-}z$ plane, respectively. Here, $x$ is the direction along the channel length, $y$ is along the channel height and $z$ is along the channel width. Upon comparing the sidewall soot foil with that of the 2D LES simulation in figure 28, the cell size obtained is comparable. In both the 2D and 3D simulations, the resulting cell size is one-half cell per channel height throughout most of the domain. Furthermore, the soot foil obtained in the 3D simulation does not exhibit a visible substructure, as was the case for the 2D simulation in figure 28. Finally, the vorticity streaks on the top wall soot foil suggest the resulting cell pattern is predominantly 2D for the thin channel.

Figure 33. Profile sections of density extracted in the $x{-}y$ , $z{-}y$ and $x{-}z$ planes for the 3D CLEM-LES.

Figure 34. Average reactant profiles for both 2D and 3D CLEM-LES. Note that distances $x$ and $x_{s}$ are normalized by $\hat{\unicode[STIX]{x1D6E5}}_{1/2}$ .

To further investigate the resulting flow pattern of the 3D simulation, various sections of density profiles are extracted in the $x{-}y$ , $z{-}y$ and $x{-}z$ planes for a single instant in time and presented in figure 33. This was done for an instant in time just before a triple point reflects on the top wall, where a pocket of unburned gas has nearly formed in the wake of the wave. In three sections taken in the $x{-}y$ plane, at $z=0.25$ , $z=1.0$ and $z=1.75$ , it is clear that the dense pocket of unreacted gas does not have a uniform shape along the width of the channel. Upon investigating the sections along various $x$ locations in the $z{-}y$ plane and $y$ locations in the $x{-}z$ plane, it is clear that the burning of the pocket of unreacted gas indeed occurs in all three dimensions. The presence of unburned gas, however, varies much more significantly in the $x$ and $y$ directions compared to the channel width in the $z$ direction. Significant burn-up along the channel width is only observed in frames where the remaining pocket is nearly consumed, sufficiently far from the leading shock. This can be seen in frames in the $z{-}y$ plane for $x\leqslant 734.5$ and the $x{-}z$ plane for $y\leqslant 5.5$ . Finally, despite the non-uniformities in the pocket shape observed across the channel width, and the three-dimensional structure of the pocket burnout, the positions of the incident and Mach shocks remain uniform throughout the channel width.

To further compare the 2D and 3D simulations, in terms of the rate at which the unreacted gas is consumed, the Favre-averaged reactant profiles are presented in figure 34. It is clear that the reactant profiles of the 3D simulation match closely the profiles obtained from the 2D simulation for $C_{\unicode[STIX]{x1D705}}=6.7$ . In fact, both the 2D and 3D simulations have an average hydrodynamic thickness of $\unicode[STIX]{x1D6E5}_{R}=8.2$ with standard deviations of $\unicode[STIX]{x1D6E5}_{R}\pm 2.7$ and $\unicode[STIX]{x1D6E5}_{R}\pm 2.0$ , respectively. To further compare 2D and 3D simulation results, the p.d.f.s of wave propagation velocity for both simulations are presented in figure 35, and compared to the experiment (Kiyanda & Higgins Reference Kiyanda and Higgins2013). It is clear that the 2D and 3D simulations both exhibit a similar p.d.f. trend. The 3D simulation however favours a slightly higher velocity $(D/D_{avg})=0.87\pm 0.04$ with a peak p.d.f. value of 3.30 compared with the 2D simulation which favoured a velocity of $(D/D_{avg})=0.84\pm 0.04$ and a peak p.d.f. value of 2.55. The two results, however, are within measurement error of each other. Since these peak p.d.f. locations are below $(D/D_{avg})<1$ , as was observed in the experiment, pockets of unburned fuel are able to form in the wake, owing to increased ignition delays at most locations on the wave front, thus lengthening the hydrodynamic structure as observed in figure 34. The results presented in figures 34 and 35 suggest that, despite the 3D nature of local burning on unburned pocket surfaces, the use of 2D LES is justified for flows with high aspect ratios, where cell sizes are sufficiently large in height compared to their channel width. This can be attributed to the large-scale flow patterns, which are predominantly 2D. Furthermore, the use of 2D LES, in this study, is further justified by the fact that the LEM subgrid accounts for the small-scale 3D shear instabilities, on the unburned gas pocket surfaces, through the stirring terms, $\dot{F_{T}}$ and $\dot{F_{Y}}$ , in (3.10) and (3.11) (Kerstein Reference Kerstein1991a ). To this effect, the LEM subgrid simulates 3D isotropic turbulent mixing according to the large-scale flow conditions. The principal caveat here is that, owing to the low-Mach-number approximation on the subgrid, shock-driven instabilities, such as Richtmyer–Meshkov, are not captured on the subgrid itself. This is currently an area of active research in the framework of one-dimensional turbulence, an alternative grid-within-a-grid approach for LES (Jozefik, Kerstein & Schmidt Reference Jozefik, Kerstein and Schmidt2016). For the high-aspect-ratio flows considered here, the large-scale 2D flow instabilities are sufficient to generate the bulk of the subgrid kinetic energy, which influences stirring events on the subgrid scales accordingly. The second caveat is that the LES scale is assumed to be sufficiently resolved to the inertial range of the Kolmogorov cascade (Frisch Reference Frisch2000). If this were the case, one might expect isotropic shear-driven mixing to dominate on the subgrid scales. However, it is noted that, while the integral scales of turbulence have been resolved here, it is not clear to what extent the turbulence characteristics of the shear layers have been captured, with respect to the Kolmogorov cascade. Furthermore, vortex breakdown, according to the Kolmogorov cascade, depends on the ability of vortices to stretch in all three dimensions. Flows that are predominantly 2D lack this ability and tend to experience significant backscatter as a result (Kraichnan Reference Kraichnan1967; Leith Reference Leith1968; Batchelor Reference Batchelor1969). It is also noted that Kolmogorov’s theory of turbulence was developed for incompressible flow. It is not yet understood how these assumptions, regarding isotropic shear-driven mixing, apply to highly compressible flows.

Figure 35. The p.d.f.s of a detonation wave having a certain velocity ( $D/D_{avg}$ ) at any given moment and location. These p.d.f.s have been compiled from both 2D and 3D simulations where $C_{\unicode[STIX]{x1D705}}=6.7$ and $H=10$ . Also shown is a p.d.f. compiled from Kiyanda & Higgins (Reference Kiyanda and Higgins2013). Note that $D_{avg}=5.53$ ( $1969.1~\text{m}~\text{s}^{-1}$ ) for the experiment and $D_{avg}=6.30\pm 0.05$ ( $2243.3\pm 17.8~\text{m}~\text{s}^{-1}$ ) for the simulations.

Figure 36. ZND temperature profiles computed for $\text{CH}_{4}+2\text{O}_{2}$ using the one-step reaction mechanism (3.13) compared to the SD Toolbox ZND solver (Kao & Shepherd Reference Kao and Shepherd2008) using the GRI-3.0 mechanism (Smith et al. Reference Smith, Golden, Frenklach, Moriarty, Eiteneer, Goldenberg, Bowman, Hanson, Song and Gardiner2016).

5.7 Validity and limitations of the global one-step reaction mechanism

In the current study, the global one-step ignition reaction mechanism, given by (3.13), although simplistic, was calibrated to reproduce the correct ZND half reaction length for a detonation travelling at the CJ velocity. Figure 36 shows the theoretical ZND temperature profile, computed using both the one-step model and the GRI-3.0 reaction mechanism (Smith et al. Reference Smith, Golden, Frenklach, Moriarty, Eiteneer, Goldenberg, Bowman, Hanson, Song and Gardiner2016), for a detonation travelling at the CJ velocity of $D_{CJ}=6.30$ ( $2240~\text{m}~\text{s}^{-1}$ ). Clearly, the correct reaction length is recovered by the one-step model, where the chemical reactions terminate around $1\unicode[STIX]{x1D6E5}_{1/2}$ (10 mm) downstream from the shock location, shown at $x=0$ . However, notable differences in temperature are observed, between the two reaction mechanisms, in the post-shock and fully burned states. Differences in the state variables of these two locations arise from the perfect-gas assumptions adopted in the governing equations, where $\unicode[STIX]{x1D6FE}$ was assumed constant. In reality, however, $\unicode[STIX]{x1D6FE}$ varies across the shock and reaction zone. Despite this, the amount of heat release, $Q$ , is consistent between the two models, yielding the correct propagation velocity. Another notable difference between the two reaction mechanisms is the temperature profile within the reaction zone itself. While the GRI-3.0 mechanism has a very distinct induction zone, with very little heat release until $x<-0.5$ ( $-$ 5 mm), the one-step model begins releasing heat as soon as the gas is shocked. This is observed by a gradual increase in temperature from $x<0$ onwards. As a result, the exothermicity of the GRI-3.0 model is confined to a much smaller layer towards the end of the reaction phase compared to the one-step model. Thus, there may be implications on the turbulence–chemistry interactions, and therefore the required value of $C_{\unicode[STIX]{x1D705}}$ , since the one-step model does not accurately capture the very short-duration energy release of methane combustion. Furthermore, this may explain the lack of substructure observed in the numerical soot foil records shown in § 5.4. Despite this deficiency, the correct post-shock ignition delay is captured by the one-step model for the range of shock velocities expected during detonation propagation. Figure 37 shows how the constant-volume ignition delays, computed from the one-step and GRI-3.0 reaction mechanisms, compare for shock speeds in the range $D=4{-}9$ (1400– $3200~\text{m}~\text{s}^{-1}$ ). Although there exists a slight underprediction of ignition times at low shock speeds, and overprediction at high speeds, the one-step model performs well at capturing the correct order of magnitude for the ignition delay times across the range of values. In fact, the one-step model works well for predicting ignition delays of methane chemistry, since there are no chain-branching cross-over effects (Browne, Liang & Shepherd Reference Browne, Liang and Shepherd2005). The sensitivity of the mixture to changes in temperature and pressure remain relatively unchanged, for the range of post-shock conditions expected.

Figure 37. Ignition delay times for various shock speeds, computed using the one-step reaction mechanism (3.13) compared to the GRI-3.0 mechanism (Smith et al. Reference Smith, Golden, Frenklach, Moriarty, Eiteneer, Goldenberg, Bowman, Hanson, Song and Gardiner2016).

5.8 Kolmogorov constant $C_{\unicode[STIX]{x1D705}}$ as a tuning parameter

From the numerical soot foils presented figure 15, it is clear that the $C_{\unicode[STIX]{x1D705}}$ parameter requires tuning in order to match the desired cell size and behaviour observed in the experiments of this study and Kiyanda & Higgins (Reference Kiyanda and Higgins2013). From a quantitative perspective, figures 30 and 31 have shown that numerically obtained p.d.f.s of detonation velocity distributions matched well the trends obtained through experiments when $C_{\unicode[STIX]{x1D705}}\geqslant 6.7$ . In general, higher values of $C_{\unicode[STIX]{x1D705}}$ allowed below-average velocities to be favoured along the wave front, as was the case experimentally. As a result, it is clear that $C_{\unicode[STIX]{x1D705}}$ requires tuning in order to ensure the correct velocity distribution, and consequently the correct expected reaction zone thickness. Furthermore, as observed in figures 17–28, a value of $C_{\unicode[STIX]{x1D705}}=6.7$ was found to be adequate to form qualitative cells and burning patterns comparable to experiments. For smaller values of $C_{\unicode[STIX]{x1D705}}$ , only much smaller cells appear. For larger values of $C_{\unicode[STIX]{x1D705}}$ , the cells increase in size until they are mode-locked by the channel height. It should be noted, however, that specific ‘matching’ of $C_{\unicode[STIX]{x1D705}}$ values can prove difficult due to the stochastic nature of both the experiments and the numerical simulations. For this, better-converged statistics would be required. It should be further noted that, even though the value of $C_{\unicode[STIX]{x1D705}}=6.7$ was found to give good agreement with experimental observation, $C_{\unicode[STIX]{x1D705}}=1.4{-}2.0$ is a more universally accepted value (Chasnov Reference Chasnov1991; Bradley Reference Bradley1992). This difference may be due to several reasons.

First of all, low values of $C_{\unicode[STIX]{x1D705}}$ correspond to Kolmogorov’s theory of the turbulence cascade in the incompressible limit (Frisch Reference Frisch2000). The detonation phenomenon under investigation here, on the other hand, is highly compressible. It is possible that highly compressible regimes generate much more turbulent fluctuations than incompressible theory predicts. This is especially true for Richtmyer–Meshkov instability, where turbulent motion is generated from shock waves interacting with material interfaces. Secondly, the large-scale motions of the experimental flow fields are essentially 2D due to the thin channel cross-section widths. In fact, the aspect ratios of the channel heights to widths ( $H/W$ ) in the experiments are ${\sim}10$ (this study) and ${\sim}4$ (Kiyanda & Higgins Reference Kiyanda and Higgins2013). In both cases, the channel widths are much smaller than the observed cell sizes. Thus, straining and decay of turbulent motions along the cross-sectional width may not be as prominent as would be expected in a full-scale 3D experiment. Furthermore, turbulence in two dimensions has been shown, theoretically, to exhibit backscatter (Kraichnan Reference Kraichnan1967; Leith Reference Leith1968; Batchelor Reference Batchelor1969), and hence generate more turbulent motions. This has also been assumed to be the case for atmospheric flows, which are also considered essentially 2D at planetary scales (Lilly Reference Lilly1966). In this sense, small-scale turbulence is able to feed into, and amplify, larger-scale turbulent motions. Although 3D burning does occur in the wake of the detonation wave at smaller scales, a 3D simulation conducted in the previous section suggested that the rate of fuel burn-up is largely governed by the large-scale fluid motions along the channel length and height. Since the experiments of this study and Kiyanda & Higgins (Reference Kiyanda and Higgins2013) are quasi-2D, it is possible that there are more turbulent motions behind the detonation front than would be expected in a naturally occurring 3D detonation wave propagation. In this sense, it is possible that the corresponding 2D simulation would also require a higher $C_{\unicode[STIX]{x1D705}}$ value to artificially generate more turbulence. Furthermore, Kraichnan (Reference Kraichnan1971) predicts a theoretical value of $C_{\unicode[STIX]{x1D705}}=6.69$ for 2D turbulence, compared to 1.4 for 3D turbulence. Another review of 2D turbulence (Danilov & Gurarie Reference Danilov and Gurarie2000) indicates that a $C_{\unicode[STIX]{x1D705}}$ value of 5.8–7 should be adopted due to the increased turbulence associated with backscatter. It should be noted, however, that in the current work this agreement is simply treated as a coincidence and has not been investigated further.

Other sources of discrepancy may arise from discrepancies associated with experimental losses, and errors associated with limitations of the chemical model used. Turbulence–flame interactions of the artificially stretched reaction zone structure, obtained from the one-step model chemistry, may have an effect on the distribution of hotspots during subgrid stirring. This should be investigated further in future work. Despite this potential deficiency, the one-step model was found to produce the correct post-shock ignition delays for the range of shock speeds expected on the wave front. Therefore, while the exact value of $C_{\unicode[STIX]{x1D705}}$ may not be captured correctly by the one-step model, the methodology is able to capture the underlying physics of how turbulence–reaction layer interactions affect changes in the detonation structure for the unsteady wave front, with varying shock speeds.

Finally, in terms of calibration, it should be noted that it is possible to implement a dynamic procedure to obtain the $C_{\unicode[STIX]{x1D705}}$ value, i.e. the localized dynamic kinetic energy model (LDKM) method (Schumann Reference Schumann1975; Menon & Calhoon Reference Menon and Calhoon1996; Menon & Kerstein Reference Menon and Kerstein2011). Although this is possible, validation and verification would still be required. Furthermore, the method may be more appropriate for 3D simulations and would therefore increase computational expense significantly. Also, the dynamic procedure may not adequately account for increased turbulent fluctuations due to compressibility effects. For now, the static method of prescribing $C_{\unicode[STIX]{x1D705}}$ proves advantageous in two dimensions for examining the effect of turbulence intensity on detonation propagation patterns and behaviours.