1 Introduction
It is well established that while detonation waves are inherently unstable and possess a three-dimensional structure, globally they maintain a steady average velocity that is remarkably close to the Chapman–Jouguet (CJ) velocity. Without the presence of cells, detonation waves can be classified as being stable, weakly unstable and highly unstable. Corresponding globally cellular detonations have also been categorized as being regular, mildly irregular and highly irregular by the regularity of the fish-scale cells that are formed via the trajectories of the triple points, as observed in Soloukhin (Reference Soloukhin1966). Eckett, Quirk & Shepherd (Reference Eckett, Quirk and Shepherd2000) studied theoretically and numerically the direct initiation of cylindrical detonations with smooth front for stable or weakly pulsating mixtures, and found that unsteadiness in the induction zone is the primary reason for the initiation failure. Kasimov & Stewart (Reference Kasimov and Stewart2005) studied detonation initiation and failure by using asymptotic analysis of weakly curved detonations, while Sow, Chinnayya & Hadjadj (Reference Sow, Chinnayya and Hadjadj2014) investigated the dynamics of detonations with smooth fronts, and found that in the presence of friction the detonation-velocity deficit increases as the mean reaction zone becomes thicker than that of the ideal Zel’dovich–von Neumann–Doring (ZND) model, and that pulsating instability enhanced by friction renders the detonation even more prone to failure. Faria & Kasimov (Reference Faria and Kasimov2015) discussed the instability of detonation in the presence of significant curvature and found that curvature plays a destabilizing role.
Furthermore, Radulescu, Sharpe & Law (Reference Radulescu, Sharpe and Law2007a ) numerically studied the direct initiation of globally planar and cellular detonations, and showed that cellular instability lengthens the global reaction zone structure and hence renders the initiation more difficult. The cell bifurcation mechanisms have also been studied experimentally and numerically (Vasil’ev & Trotsyuk Reference Vasil’ev and Trotsyuk2003; Vasil’ev, Vasiliev & Trotsyuk Reference Vasil’ev, Vasiliev and Trotsyuk2010), showing that the formation of new transverse waves on the expanding detonation front is preceded by cellular instability. Vasil’ev et al. (Reference Vasil’ev, Vasiliev and Trotsyuk2010) observed the bifurcation properties of the detonation cells, and demonstrated that the presence of small-scale cells and additional hot spots at the points of collision of transverse waves do not lead to any noticeable enhancement in the initiation. On the other hand, Ng et al. (Reference Ng, Kiyanda, Morgan and Nikiforakis2015) observed the opposite trend on the initiation of cellular instability by considering different wavelengths in the initial perturbation. As such, the role of cellular instability on detonation initiation remains unclear.
We next note that most previous theoretical (Lee & Stewart Reference Lee and Stewart1990; Bourlioux & Majda Reference Bourlioux and Majda1992; Short & Stewart Reference Short and Stewart1998; Short Reference Short2001; Clavin & Denet Reference Clavin and Denet2002; Clavin & Williams Reference Clavin and Williams2012), experimental (Lee Reference Lee2008) and computational (Sharpe & Falle Reference Sharpe and Falle2000; Ng et al. Reference Ng, Higgins, Kiyanda, Radulescu, Lee, Bates and Nikiforakis2005; Deledicque & Papalexandris Reference Deledicque and Papalexandris2006; Henrick, Aslam & Powers Reference Henrick, Aslam and Powers2006; Radulescu et al. Reference Radulescu, Sharpe, Law and Lee2007b ; Sharpe & Radulescu Reference Sharpe and Radulescu2011) studies were mainly concerned with globally planar detonations in channels and square tubes, while in many situations the detonation fronts are globally curved in open cylindrical and spherical geometries. It is therefore of interest to assess the role of the front curvature on the pulsating and cellular instabilities of the front. Specifically, Watt & Sharpe (Reference Watt and Sharpe2003) studied the effect of pulsating instability of weakly curved detonations through linear instability analysis, and subsequently (Watt & Sharpe Reference Watt and Sharpe2005) showed computationally that the presence of curvature associated with a radially expanding detonation can induce pulsating instability on the front which is otherwise stable in the planar propagation. Experimental investigations (Lee et al. Reference Lee, Knystautas, Guirao, Bekesy and Sabbagh1972) also showed that cylindrical and spherical detonation fronts have cellular structures, which sustain self-propagating diverging detonation by matching the increasing front area with the correspondingly increasing triple-shock units as the detonation expands. Recently, numerical investigations on two-dimensional cylindrical detonations (Jiang et al. Reference Jiang, Han, Wang and Zhang2009; Asahara, Tsuboi & Hayashi Reference Asahara, Tsuboi and Hayashi2010) demonstrated that cell divergence is caused by the local explosions of unreacted fuel pockets between the shock front and the reaction layer. There is also considerable interest in the role of diffusion in the cellular cylindrical detonation structure. Specifically, while the role of diffusion in the detonation wave with regular cell structure is usually weak because almost all the gas passing through a globally planar detonation is burnt, unreacted pockets (Gamezo, Desbordes & Oran Reference Gamezo, Desbordes and Oran1999) are produced in the highly irregular detonation propagating in channels, and as such their burning rate can be enhanced by the turbulent mixing of burnt and fresh gases at their boundaries via vortices produced by the Richtmyer–Meshkov (R–M) and Kelvin–Helmholtz (K–H) instabilities. Radulescu et al. (Reference Radulescu, Sharpe, Law and Lee2007b ) demonstrated numerically and experimentally that in highly irregular detonations, turbulent mixing promotes consumption of the unreacted gases, as suggested in Gamezo et al. (Reference Gamezo, Desbordes and Oran1999) and Mahmoudi & Mazaheri (Reference Mahmoudi and Mazaheri2012).
Following these previous worthwhile contributions, the present work aims to further identify the separate and coupled roles of cellularity and global curvature on the structure as well as the global propagation velocities of planar and cylindrical detonations, by simulating and comparing these properties with the corresponding situations in which the cells are suppressed. We shall also address the role of diffusion in the weakly unstable cylindrical detonation. In the following we shall first state the governing equations and the numerical method, and then present and discuss the results.
2 Governing equations
In two-dimensional (2-D) coordinates, the equations of mass, momentum and energy conservation as well as the heat transfer equation are as follows:
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\boldsymbol{U}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}\boldsymbol{F}(\boldsymbol{U})}{\unicode[STIX]{x2202}x}+\frac{\unicode[STIX]{x2202}\boldsymbol{G}(\boldsymbol{U})}{\unicode[STIX]{x2202}y}=S,\end{eqnarray}$$
where the conserved variable vector
$\boldsymbol{U}$
, the flux vectors
$\boldsymbol{F}$
and
$\boldsymbol{G}$
and the source term
$S$
are given, respectively, by
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{U}=(\unicode[STIX]{x1D70C},\unicode[STIX]{x1D70C}u,\unicode[STIX]{x1D70C}v,E,\unicode[STIX]{x1D70C}Y)^{T}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{F}(\boldsymbol{U})=(\unicode[STIX]{x1D70C}u,\unicode[STIX]{x1D70C}u^{2}+p,\unicode[STIX]{x1D70C}uv,u(E+p),\unicode[STIX]{x1D70C}uY)^{T}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{G}(\boldsymbol{U})=(\unicode[STIX]{x1D70C}v,\unicode[STIX]{x1D70C}vu,\unicode[STIX]{x1D70C}v^{2}+p,v(E+p),\unicode[STIX]{x1D70C}vY)^{T}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle S(\boldsymbol{U})=(0,0,0,0,\unicode[STIX]{x1D70C}\unicode[STIX]{x1D714})^{T}. & \displaystyle\end{eqnarray}$$
Here
$u$
and
$v$
are the Cartesian components of the fluid velocity in the
$x$
and
$y$
directions respectively and
$\unicode[STIX]{x1D70C}$
is the density,
$p$
the pressure,
$T$
the temperature given by the ideal gas law,
$p=\unicode[STIX]{x1D70C}R_{u}T$
and
$E$
the total energy per unit volume, given by
$$\begin{eqnarray}\displaystyle E=\frac{p}{(\unicode[STIX]{x1D6FE}-1)}+\frac{1}{2}\unicode[STIX]{x1D70C}(u^{2}+v^{2})+\unicode[STIX]{x1D70C}YQ, & & \displaystyle\end{eqnarray}$$
with
$Y$
being the reactant mass fraction,
$Q$
the heat of reaction per unit volume and
$\unicode[STIX]{x1D6FE}$
the specific heat ratio, which are respectively taken as 50 and 1.2. The reaction rate per unit mass is given by the Arrhenius law,
$\unicode[STIX]{x1D714}=-KY\exp (-E_{a}/R_{u}T)$
, where
$E_{a}$
is the activation energy,
$R_{u}$
the gas constant,
$M$
the molecular weight and
$K$
the constant pre-exponential factor. The equations are made dimensionless based on the state of the unburned gas, as:
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \overline{\unicode[STIX]{x1D70C}}=\unicode[STIX]{x1D70C}/\unicode[STIX]{x1D70C}_{0},\quad \overline{p}=p/p_{0},\quad \overline{T}=T/T_{0},\quad \overline{u}=u/u_{0},\\ \overline{p}=v/u_{0},\quad \overline{x}=x/L_{1/2},\quad \overline{t}=t/t_{0},\\ \displaystyle \overline{E}=E/p_{0},\quad \overline{K}=Kx_{0}/u_{0},\quad \overline{Q}=Q/u_{0}^{2},\quad \overline{E}=E_{a}/R_{u}T_{0},\end{array}\right\} & & \displaystyle\end{eqnarray}$$
where
$u_{0}=\sqrt{R_{u}T_{0}}$
and
$t_{0}=L_{1/2}/u_{0}$
. The reference length
$L_{1/2}$
is half of the reaction zone length. Because of self-similarity of the Euler equations, the dimensionless and the original dimensional forms are identical except the reaction rate becomes
$\unicode[STIX]{x1D714}=-KY\exp (-E_{a}/T)$
. Consequently, for simplicity in notation the overbars on the variables are dropped in the following.
For steady planar detonation, the CJ detonation velocity is given by
$$\begin{eqnarray}\displaystyle D_{CJ}=\sqrt{(\unicode[STIX]{x1D6FE}^{2}-1)Q/2}+\sqrt{(\unicode[STIX]{x1D6FE}^{2}-1)Q/2+\unicode[STIX]{x1D6FE}p_{0}/\unicode[STIX]{x1D70C}_{0}}. & & \displaystyle\end{eqnarray}$$
For a curved front of radius,
$R$
, the detonation speed is given by He & Clavin (Reference He and Clavin1994), He (Reference He1996), Yao & Stewart (Reference Yao and Stewart1996) and Law (Reference Law2006).
In the following simulation, an activation energy
$E_{a}=24$
is selected, with neutral stability boundary on the ZND model (Lee Reference Lee2008), and the corresponding detonation is weakly unstable. We note herein that Radulescu & Lee (Reference Radulescu and Lee2002) experimentally observed two major classes of detonation behaviour: weakly and highly unstable detonations. Bourlioux & Majda (Reference Bourlioux and Majda1992) have shown experimentally and numerically that the degree of irregularity of the detonation cell depends on the activation energy. The importance of the physical or numerical diffusion in the detonation structure is closely related to the regularity associated with the activation energy (Oran et al.
Reference Oran, Weber, Stefaniw, Lefebvre and Anderson1998; Gamezo et al.
Reference Gamezo, Desbordes and Oran1999; Arienti & Shepherd Reference Arienti and Shepherd2005), while Radulescu et al. (Reference Radulescu, Sharpe, Law and Lee2007b
) and Mazaheri, Mahmoudi & Radulescu (Reference Mazaheri, Mahmoudi and Radulescu2012) have shown that physical diffusion is not very important in stable and weakly unstable detonations. Therefore, inviscid reactive flow adopted in the present study should be sufficient to describe the dynamics of weakly unstable cylindrical detonations.
3 Numerical method
To numerically solve the governing equations, we apply the fifth-order local characteristics based on the weighted essentially non-oscillatory (WENO) conservative finite-difference scheme (Jiang & Shu Reference Jiang and Shu1996) with third-order total variation diminishing (TVD) Runge–Kutta time discretization (Shu & Osher Reference Shu and Osher1988; Balsara & Shu Reference Balsara and Shu2000). Relative grid resolution and preserving-positivity scheme are verified for our numerical scheme so that it can reach the fifth order (Wang et al. Reference Wang, Zhang, Shu and Ning2012, Reference Wang, Shu, Han and Ning2013; Zhang & Shu Reference Zhang and Shu2012).
For the simulation of the expanding cylindrical detonation without cellular structure, a grid resolution of 20 points per half-reaction length is used, which is sufficient for convergence (Watt & Sharpe Reference Watt and Sharpe2003) and captures the ZND structure and peak pressure (Wang et al.
Reference Wang, Shu, Han and Ning2013). Reflective boundary conditions are set at
$R=0$
; free outflow condition at
$R=\infty$
. The initial conditions are given by
$\unicode[STIX]{x1D70C}=1$
,
$u=0$
and
$Y=1$
everywhere, and a high-energy source
$E_{s}$
is set at the origin with
$p=p_{source}$
for
$R=R_{source}$
and
$p=1$
for
$R>R_{source}$
, where
$R_{source}=5L_{1/2}$
is the half-width of the reaction zone,
$L_{1/2}$
.
For the 2-D simulation allowing for cellular structure, Cartesian coordinates are used. We take the quarter area of the numerical domain (
$1500\times 1500$
;
$x\geqslant 0,y\geqslant 0$
) and set a symmetric boundary condition at the lines:
$x=0$
and
$y=0$
, while the other sides have free outflow conditions, and radius
$R^{2}=x^{2}+y^{2}$
. Initial conditions are similar to those of the 1-D simulation, with the high-energy source
$E_{s}$
in the region
$R<R_{source}$
. The total number of cells with grid resolution of 20 points is approximately 900 million. Due to the very large size of the computational domain and the large number of grid points for the 2-D simulations, it is essential to implement the high-resolution and efficient parallel code on a massive parallel platform, as is done in this work by 900 processors.
4 Results and discussion
To facilitate the following discussions, we shall designate the planar and cylindrical detonations with smooth fronts as 1-D planar and 1-D cylindrical detonation, respectively and the corresponding propagation with cellular fronts as 2-D planar and 2-D cylindrical detonations.
Pulsating evolution of the shock pressure (a) and detonation velocity (b) with radius.
4.1 Role of curvature on pulsating instability
To identify the effects of curvature on pulsating instability, we first compare results of the propagation of the 1-D planar and 1-D cylindrical detonations. With the 1-D planar detonation being stable for
$E_{a}=24$
, results shown in figure 1 confirm the earlier study (Watt & Sharpe Reference Watt and Sharpe2005) that (positive) curvature promotes pulsating instability in that the 1-D cylindrical detonation exhibits pulsation for a relatively large range of
$E_{s}$
near the critical source energy. Furthermore, the 1-D cylindrical detonation becomes more stable with increasing radius and hence decreasing curvature, with the shock pressure and velocity approaching the CJ state. The mechanism of longitudinal pulsation has been explained for the 1-D planar detonation (McVey & Toong Reference McVey and Toong1971; Alpert & Toong Reference Alpert and Toong1972; Law Reference Law2006). Figure 2 shows evolution of the pressure profile for the 1-D planar and 1-D cylindrical detonations, in the upper and lower panels respectively. Obviously, pressure jumps (see solid line) and rapidly decays through the thin reaction layer denoted by the mass fraction of the reactant in both planar and cylindrical detonations. Furthermore, for the 1-D cylindrical detonation, the pressure profile transforms into a square wave at
$R\sim 200$
–290 as the detonation decays from the initial overdriven state, demonstrating that the detonation structure consists of a thick induction zone followed by a short reaction layer. Because of the heat release during this period, the leading edge of the pressure profile steepens and the induction thickness (referring to the thickness of the Neumann zone) decreases, indicating an increase in the reaction rate. Eventually, the detonation pressure peaks and then decays as the detonation is weakened by the downstream rarefaction wave. The decay, however, increases the downstream flow velocity, which would eventually reach the sonic state and thereby cut off the downstream influence. This then leads to re-intensification of the reaction rate, and the process repeats itself, hence constituting the observed pulsation (Alpert & Toong Reference Alpert and Toong1972; Law Reference Law2006). This therefore confirms the results of Watt & Sharpe (Reference Watt and Sharpe2005) and Faria & Kasimov (Reference Faria and Kasimov2015) obtained with similar simulations, and in addition provides the mechanistic explanation for the pulsating instability. As a corollary of this understanding, planar detonations do not exhibit pulsating instability because their induction thickness is constant.
Pulsating instability can also suppress successful detonation initiation, as observed in the 1-D simulation of Watt & Sharpe (Reference Watt and Sharpe2005) and explained in the studies of He & Clavin (Reference He and Clavin1994) and Ng & Lee (Reference Ng and Lee2003). Specifically, He & Clavin (Reference He and Clavin1994) theoretically showed that there exists a critical initiation radius below which generalized Chapman–Jouguet solutions do not exist, although subsequent study (Eckett et al. Reference Eckett, Quirk and Shepherd2000) showed a non-uniqueness of the critical initiation energy in some cases and as such suggested that the quasi-steady assumption could have limited direct initiation. Furthermore, Eckett et al. (Reference Eckett, Quirk and Shepherd2000) found that, for direct initiation with near-critical initiation energy, the detonation initially forms as the blast wave approaches the quasi-steady solution, and then is destroyed quickly by the pulsating instability. Consequently, while initiation of the 1-D cylindrical detonation without the cellular structure is affected by the pulsating instability, the role of cellular instability in the initiation and propagation of 2-D cylindrical detonations need to be identified, to be discussed later.
Evolution of pressure profile for the planar detonation (a) and cylindrical (b) 1-D detonation. The same source energy (
$E_{s}=12\,000$
) is used for both cases.
Figure 3 shows the effects of curvature and source energy on the shock pressure and propagating velocity. It is seen that, for the source energy of
$E_{s}=30\,000$
, decay of the average velocity of the planar detonation is slower than that of the cylindrical detonation case, while they both remain strongly overdriven, without exhibiting any pulsation, and decay to the CJ state. However, as the source energy is reduced to
$E_{s}=12\,000$
, while pulsation remains inhibited for the planar detonation shown in figure 2(a), it appears in the cylindrical detonation, with progressive damping of the oscillation amplitude, as indicated by the black solid lines in figure 3(a,b). The detonation velocity initially rapidly decays below the CJ state and then gradually approaches it as the curvature decreases. It is also noted that some kinks are present during early decay for both the planar and cylindrical detonations in the 1-D simulation, shown in figure 3(a,b). These are caused by disturbances from the rear boundary and their formation was explained in Short, Kapila & Quirk (Reference Short, Kapila and Quirk1999) and Han, Gao & Law (Reference Han, Gao and Law2015). Specifically, formations of these kinks, or more exactly small pulsations in the rapid shock decay period, have similar mechanisms with ordinary 1-D detonation pulsations. However, they subsequently disappear due to the presence of the transverse flow in 2-D simulations.
Shock pressure (b) and average detonation velocity (a) versus propagating distance for 1-D detonation.
Evolution of cellular structure and cylindrical front: (a) maximum pressure; (b) and (c) the pressure and temperature contours at
$t=55.9$
, 108.3, 159.3 and 209.7; (d) zoom-in of the dashed box, for the pressure, temperature and mass fraction contours.
4.2 Role of curvature on cellular instability
Since pulsation does not occur for the 1-D cylindrical propagation with
$E_{s}=30\,000$
, we shall use this value to study the role of curvature on the cellular instability. Figure 4(a) shows that evolution of the cellular structure for the 2-D cylindrical detonation consists of three stages: no cell, growing cell size and diverging cells. Specifically, initially no cell forms due to the strong overdrive controlled by the high source energy. As the detonation decays, the effect of curvature amplifies the transverse disturbance, and leads to the formation of cells. For
$R\sim 120$
–500, the cells grow larger as the detonation front expands. For
$R>500$
, as the overdrive and the global curvature further decrease, the cells diverge in order to match the area of the expanding front (Jiang et al.
Reference Jiang, Han, Wang and Zhang2009), as indicated by the increasing number of triple points. During the entire detonation evolution, the global curvature mainly affects the cell growth and the initial cell-diverging stages. Its effect becomes weak in the late cell-diverging stage as the detonation front expands and the transverse waves become dominant, with unreacted pockets appearing behind the Mach shock of the detonation structure, as respectively shown in figure 4(b,d). This is particularly prominent for relatively small radius, and hence large curvature, due to the reduced residence time behind the short Mach shock in the larger cells, as for
$R\sim 400$
–600 shown in figure 4(c). For larger radius and hence smaller curvature, such as
$R\sim 800$
–1100, the unreacted pockets mainly appear behind the highly curved front regions. At
$R\sim 1450$
, the number of unreacted pockets becomes even less due to weakening of the curvature. Figure 4(d) shows the pressure, temperature and mass fraction contours by enlarging the boxed regime in figure 4(b). It is seen that globally the cellular structure consists of the incident shock (IS), the Mach shock (MS) and the transverse wave (TW) extending to the downstream and the shear layers behind the MS, which resembles the structure of unstable cellular detonation in Radulescu et al. (Reference Radulescu, Sharpe, Law and Lee2007b
). It is further seen that the induction time of the gas shocked by the IS is much longer than that of the gas shocked by the MS, leading to reaction delay; see the temperature contour in figure 4(d). Consequently, the frontal cellular structure is sustained by the rapid heat release behind the MS and the delayed heat release behind the IS. Furthermore, finer structures are also observed for the MS, as shown in figure 4(d). This finer structure originates mainly from the downstream instabilities at earlier times, e.g. hydrodynamic instabilities led by an entropy-vortex roll created after the collision of two triple points or the explosion in the downstream unreacted pocket. As the two triple points connecting the IS approach each other, the excess unreacted gas enters into the shear layer. After their collision, the IS transforms to a strong Mach shock, leaving the unreacted pocket behind it due to the short residence time, as substantiated by the unreacted pockets observed behind the collisions 1 and 2, shown in figure 4(d).
Effect of the curvature on the average velocity of the planar and cylindrical detonations. Red line and diamond are the average velocity of 1-D and 2-D cylindrical shock wave, respectively.
Figure 5 compares the 1-D and 2-D average velocities of the planar and cylindrical detonations, with the same source energy. The average velocity is determined by calculating the ratio of the distance difference between two average shock radii, which is obtained by integrating across the entire cylindrical front to the corresponding time difference. It is seen that while the 1-D and 2-D propagation velocities largely agree with each other for the planar detonation that still is strongly overdriven, the 1-D propagation velocities are uniformly higher than those of the 2-D propagation for the cylindrical detonation. Furthermore, for the cylindrical detonation, it is seen that while the 1-D and 2-D propagation velocities are largely similar during the initial stage when they are mainly controlled by the initial overdrive, and during the latter stage as they decrease and approach the CJ value, the 2-D velocities are smaller than the corresponding 1-D velocities in the intermediate period. However, the numerical simulation of the inert shock expansion in cylindrical geometries shows that the 1-D and 2-D velocities globally follow a uniform decay with minor disagreement that is caused by weak pulse in 1-D case and weakened significantly due to absence of the reaction, see red diamond and line in figure 5. This demonstrates that for the cellular cylindrical detonation the cellular instability causes more deficit of the average velocity in the regime with large curvature and hence more difficulty to attain direct initiation than in the 1-D cylindrical detonation.
The above discussions show that unreacted pockets caused by the cellular instability produce more deficit of the 2-D cylindrical detonation. It is noted that Mahmoudi, Mazaheri & Parvar (Reference Mahmoudi, Mazaheri and Parvar2013) and Mahmoudi & Mazaheri (Reference Mahmoudi and Mazaheri2015) also numerically observed unreacted pockets behind the double-Mach shock near the wall for the global planar detonation with highly irregular cellular structure. In the present simulation, curvature enhances cellular instability by modifying the globally cylindrical detonation structure and leads to unreacted pockets. These take on a tongue-like pattern surrounded by the sheared-layer boundary, which resembles that of Mazaheri et al. (Reference Mazaheri, Mahmoudi and Radulescu2012), while a small-scale vortex caused by the R–M instability is also observed along the shear-layer boundaries. Nevertheless, the simulations of Mahmoudi et al. (Reference Mahmoudi, Mazaheri and Parvar2013), Mahmoudi & Mazaheri (Reference Mahmoudi and Mazaheri2015) also showed that burning unreacted gas receding from the shock compression is related to physical diffusion that modifies the shear-layer structure. Hence, we have also assessed the role of the diffusion in the detonation structure of the 2-D cylindrical detonation in the present simulation, presented in appendix B.
Average detonation velocities of planar detonations with and without cellular front.
For the 2-D planar detonation in figure 5, with the domain width being the same as that of Radulescu et al. (Reference Radulescu, Sharpe, Law and Lee2007b
), it is seen that the upper and lower periodic boundary condition weakens the wall confinement effect. A lower source energy (
$E_{s}=3000$
) is used to further identify the effect of the cellular instability with the small overdrive. Figure 6 shows that the average velocity of the 2-D simulation still is close to that in the 1-D case, although the instant velocity oscillates around the CJ state as the overdriven detonation decays to the CJ state. It is however noted that Radulescu et al. (Reference Radulescu, Sharpe and Law2007a
) suggested that cellular instabilities can affect the initiation due to unreacted pockets. The difference is that their case involves a highly irregular detonation with
$E_{a}=27$
such that some large unreacted pockets are generated, which is also enhanced by the confinement effect of the wall boundary condition, while the present investigation is on the weakly irregular detonation, with
$E_{a}=24$
. Consequently, more deficit of the 2-D velocity is not observed for the planar detonation with a weak confinement as the detonation decays to the CJ state (figure 6), and the cellular instability has no effect on direct detonation initiation. On the other hand, in the 2-D cylindrical detonation, since the effective activation energy increases due to the global curvature effect, cellular instability is enhanced consequently and has significant effect on the initiation, as further discussed next.
Average 1-D and 2-D velocities of the cylindrical detonations at different source energies.
4.3 Effect of cellular instability on initiation of 2-D cylindrical detonation
Eckett et al. (Reference Eckett, Quirk and Shepherd2000), Watt & Sharpe (Reference Watt and Sharpe2005) and Sow et al. (Reference Sow, Chinnayya and Hadjadj2014) have shown that successful initiation of the 1-D cylindrical detonation is associated with pulsating instability. However, for the 2-D cylindrical detonation, its initiation is more difficulty than predicted by the 1-D case due to the cellular instability; see figure 5. We conducted 2-D simulations by lowering the source energy to examine further the negative effect of the cellular instability on the 2-D initiation. Figure 7 shows the average velocity of the 2-D cylindrical detonation at different source energies. It is seen that for
$E_{s}=27\,500$
, the 2-D cylindrical detonation has faster decay of the average velocity than that of the 1-D case, see figure 7. The cellular detonation forms and remains sustaining, shown in figure 8(a), while its velocity is slightly lower than the CJ value, as indicated by the blue symbols and line in figure 7. For
$E_{s}=25\,000$
, while the 1-D detonation still propagates at the overdriven speed (red line in figure 7), the 2-D case (red symbols in figure 7) decays rapidly and subsequently quenches at
$x\sim 400$
, which also is substantiated by the absence of cells in the far field of figure 8(b). As
$E_{s}$
is reduced to 20 000, the 2-D detonation (black symbols in figure 7) fails rapidly at
$x\sim 100$
while the 1-D detonation (black line in figure 7) can still propagate, being at the CJ state at
$x\sim 650$
. Consequently, cellular instability tends to weaken the initiation of the 2-D cylindrical detonation, as shown in figure 8(b), while initiation also fails easier for reduced source energy.
It is also of interest to assess the dependence of the initiation and diverging propagation on the initial perturbation. Figure 9 shows the decay of the detonation velocities with and without perturbation. In the simulation, a sinusoidal pressure profile is imposed initially along the circumference of the front of the source energy in order to promote the development of the cellular instability. It is seen that in the strongly overdriven stage (
$x<100$
), the average velocities of the two simulations all have consistent decay rate for
$E_{s}=27\,500$
(figure 9). After the cell forms at
$x\sim 100$
, the two velocities continue to have global consistence with slight differences, as shown in figure 8(a). For
$E_{s}=25\,000$
, the initially imposed perturbation still does not have any obvious effect on the decay of the average velocity, even though globally the detonation undergoes transition from overdriven detonation to quenching detonation. Consequently, for the present simulation, initiation or quenching of the detonation is almost independent of the initial condition.
Cell evolution of cylindrical detonations at different source energies: (a)
$E_{s}=27\,500$
; (b)
$E_{s}=25\,000$
.
Detonation velocities for cases with (diamonds) and without (squares) initial perturbations for the 2-D cylindrical detonation.
4.4 Propagation and transverse instability caused by global curvature in 2-D cylindrical detonation
The detonation structures in the cell-growing and cell-diverging stages shown in figure 4 are shown in figure 10, which explain the difference between the 1-D and 2-D average velocities of the cylindrical detonations, as well as further reveal effects of curvature on the structure of the 2-D cylindrical detonation. It is seen that the cylindrical front is initially smooth and small perturbations constantly propagate upstream into the detonation front. The triple points appear regularly at the expanding front for
$R\sim 322$
(
$x,y\sim 228$
), followed by the uniform, unreacted gas caused by the large curvature, as shown in figure 10. The unreacted pockets are rapidly consumed by the large reaction rate of the source energy, which renders the front structure relatively thin. As the degree of overdrive decreases, effects of the global curvature increase through thickening of the detonation structure so as to accommodate the fresh gas entering the expanding front. In the thickened detonation structure, unreacted pockets behind the Mach stem over the expanding front emerges at
$t=108.3$
–125.5, which in turn leads to more velocity deficit in the 2-D than the 1-D cases; see figures 5 and 10. At
$t=209.7$
, the volume of the unreacted pockets is small as the effect of the curvature weakens and the globally overdriven cylindrical detonation approaches the CJ state (figure 10). The average velocity of the 2-D cylindrical detonation is also close to that of the 1-D cases at
$R\sim 1450$
(figure 5). Consequently, since the 1-D and 2-D average velocities of the cylindrical detonations would eventually be identical as the radius approaches infinity, the cause of the increased deficit of its averaged velocity, as compared to the 1-D case, is the unreacted pocket behind the 2-D detonation front due to the global curvature.
Unreacted pockets behind globally cylindrical expanding front in the cell-growing and cell-diverging stages.
Recognizing that curvature increases longitudinal instability in the range of radius with relatively large curvature (Ng et al.
Reference Ng, Radulescu, Higgins, Nikiforakis and Lee2005), we next study the physical variables, spatially integrated along the cylindrical circumference (
$0<R<900$
), to identify whether the curvature still has an effect on the averaged detonation structure. Figure 11 shows evolution of the average pressure at
$R\sim 100$
–900. Globally, figures 10 and 11 show that the detonation structure (
$350<R<700$
) becomes thicker. The average pressure, temperature and Mach number for L1–L6 shown in figure 11 are calculated to closely identify the subtle changes in the average longitudinal thickness (Radulescu et al.
Reference Radulescu, Sharpe, Law and Lee2007b
).
Evolution of pressure averaged along the cylindrical circumference: L1–L6 denote the locations of the detonation front.
Figure 12 shows the average Mach number, temperature and pressure for locations L1–L6 in figure 11. The average Mach number relative to the detonation front is calculated by using the average sound speed. Defining the average distance between the detonation front and the average sonic plane where the Mach number is unity as the hydrodynamic thickness (Radulescu et al.
Reference Radulescu, Sharpe, Law and Lee2007b
), figure 12(a) then shows that compared to the case of Radulescu et al. (Reference Radulescu, Sharpe, Law and Lee2007b
), the location of the sonic plane is delayed further downstream, being affected by the global curvature. For example, for the case of L1, the location is delayed from
${\sim}20$
in the planar case to
${\sim}45$
in the cylindrical case. As such, the average structure consists of a main reaction layer followed by a long tail in which reaction gradually tends to the equilibrium state at the sonic plane. It is further found that the hydrodynamic thickness is approximately
${\sim}50$
for L1, 53 for L2, then decreases to 40 for L3 due to reduction in the degree of the overdrive, and increases to
${\sim}50$
for L4 and L5 because of unreacted pockets, and finally shortens to
${\sim}40$
at L6 with divergence of the transverse wave and reduction of the curvature; see the small plot in figure 12(a). A closer examination of the thickness of the main reaction layer, namely the detonation front structure, shows that its thickness also changes, from approximately 2 for L1, then increases to
${\sim}4$
for L2,
${\sim}8$
for L3 and 14 for L4 and L5, and eventually shortens to
${\sim}6$
for L6; see the small plots in figure 12(a). This change of the thickness of the main reaction layer is reasonable in that, from L2 to L5, the fresh mixture entering the overdriven detonation front with the large expanding curvature cannot be consumed rapidly by the finite reaction rate, hence leading to a larger thickness at
$R\sim 350$
–700; the thickness subsequently decreases with the reduction of curvature and the diverging transverse wave, as shown for L6. The study of Gamezo et al. (Reference Gamezo, Desbordes and Oran1999) showed that for the globally planar detonation with irregular cellular structures, the average reaction zone is larger and the maximum reaction rate is lower than those in the 1-D case, which also is substantiated by comparison of the 1-D and 2-D results in figure 12(a,b). This effect also has been shown to increase with the activation energy due to the larger unreacted gas pockets. In the present results, the curvature can increase the effective activation energy and thus enhance the cellular instability. Consequently, elongation of the average reaction zone also is much more obvious than the ZND profile of the 1-D simulation shown in figure 12(a), resembling those of the highly unstable cases of Gamezo et al. (Reference Gamezo, Desbordes and Oran1999) and Radulescu et al. (Reference Radulescu, Sharpe, Law and Lee2007b
).
The above discussion demonstrates that on average the curvature extends the induction length which increases the extent of the instability. Together with the transverse detonation structure, these two factors result in unreacted pockets. Furthermore, within the extended detonation structure, temperature and pressure perturbations due to the downstream local explosion of unreacted gas can propagate upstream and eventually lead to the appearance of new transverse waves. This is further considered next.
Profiles of (a) Mach number, (b) temperature and pressure averaged along the cylindrical circumference: (a) dash line denotes
$M_{a}=1$
; (b) solid lines in the upper half of panel are temperature corresponding that in the lower half of panel. Dot-dot-dash line is the ZND profile of 1-D simulation, and reaction region (R) and induction region (I) are labelled for the 1-D simulation.
Figure 13 shows the formation of unreacted pockets and new transverse waves for the cellular cylindrical detonations in the region (
$x\sim 600$
–750,
$y\sim 0$
–300). From the pressure contour in figure 13(a), the triple-point structure, consisting of the IS, MS and TW, appears at the curved expanding front. Referring to
$t=107.6$
, it is seen that the detonation globally has a thicker front structure due to the curvature, as observed from the main reaction layer in figure 12(a). With the transverse wave movement, the MS with the curved expanding front now covers more fresh gas. At
$t=114$
, the strong local explosion induced by the collision between TP1 and TP2 causes the front to abruptly expand, and thus is unable to rapidly consume the excess fresh gas in the thickened structure. The bulk of the unreacted pocket appears downstream and its temperature (
${\sim}4$
) is very low relative to that in the main reaction zone, as observed from the temperature and mass fraction contours in figure 13(b,c), which leads to the average velocity deficit of the 2-D cylindrical detonation.
By inspecting the curved shock front, we find that cusps appear at the front and as such induce new transverse waves; see figure 13(c). Observed from the mass fraction and temperature contours, it is seen that a small unreacted pocket fills the cusp and its local explosion generates new transverse waves, as shown in the track of figure 13(a). Furthermore, the expanding front of the MS reveals the cusp formation. For location 1 in figure 13(c), the MS front consists of small perturbations, see figure 13(b,c), which are mainly caused by the downstream unreacted pockets because the flow in the thick structure is subsonic, with average
$M_{a}<1$
from the front to the sonic surface, as shown in figure 12(a). With expansion of the curved front, perturbation develops and a small kink appears at the locations 2–4. At location 5, a small cusp with unreacted pocket forms and explodes locally, such that eventually new transverse waves appear at location 6. The generation of new transverse waves at the cusp does not add new unreacted pockets in the downstream flow, while globally it reduces the thickness of the detonation front structure. This also substantiates the result that the thickness of the main reaction layer and the average detonation structure decrease as the cells diverge steadily, as shown in figures 4 and 12.
In summary, the lengthening of the average detonation structure caused by the global curvature, coupled with the transverse structure, results in the appearance of unreacted pockets as the fresh mixture entering the detonation front cannot be completely burned in this region, as shown in figures 10 and 13.
Divergence of transverse wave and formation of un-reacted pocket: dashed lines denote tracks of new transverse waves. Incident Wave (IS) and Mach Shock (MS), and TP1 and TP2 are triple points; 1–6 denote locations of MS front.
5 Conclusions
In this work, effects of the global curvature on the cylindrically propagating detonation are computationally investigated. Results from the 1-D simulation show that pulsation is mainly caused by the increase in the instability factor from the induction thickness due to the global curvature. For the 2-D simulation, cell evolution approximately follows the three stages of no cell, growing cells and diverging cells. By comparing the differences in the average 1-D and 2-D velocities, it is noted that the latter has more velocity deficit in the cell-growing stage, and that this deficit is mainly due to the incomplete reaction behind the detonation front. From the average detonation structure, it is shown that the global curvature increases the average instability factor by extending the induction zone, creating unreacted pockets and hence large deficits of the averaged velocity. In the stage of diverging cells, local explosion in the cusp of the expanding shock front creates new transverse waves and consequently reduces the detonation thickness. With continuous propagation the downstream unreacted pocket is progressively consumed, leading to the eventual convergence of the average velocity determined from the 1-D and 2-D simulations.
The present study further demonstrates that the cellular instability renders the initiation more difficult in the cylindrical detonation, based on the simulation of the weakly unsteady detonation with the activation energy of
$E_{a}=24$
, for which the diffusion effect in the burning of unreacted pockets is very weak, as shown in figure 17. This leads to unreacted pockets and hence the observed difficulty in initiation. The tendency, however, may be reversed if the mixture is sufficiently unstable (cf.
$E_{a}>30$
as in many real hydrocarbon mixtures) such that the cellular instability is easier to be amplified because significant flow instability with high levels of compressible turbulence (Kessler, Gamezo & Oran Reference Kessler, Gamezo and Oran2011) or the embedded hydrodynamic fluctuation can promote cell formation (Qi & Chen Reference Qi and Chen2016). Its simulation requires much higher numerical resolution to resolve small-scale turbulence and capture shock–shock and shock–vortex interactions, which can break down the unburnt pockets (Kessler et al.
Reference Kessler, Gamezo and Oran2011) and diffusion effect should also play a prominent role. This problem merits further study.
Acknowledgements
This work was supported by the National Science Foundation of China under grant nos 91541206, 91441131, 51206088 and 11325209 and the National Key Basic Research Program of China (no. 2014CB239603).
Appendix A. Verification of grid resolution convergence
In the overdriven regime, the reaction length scale is much thinner than that of the CJ wave and hence a higher grid resolution is required. In the stage of initial decay of strongly overdriven detonation, the cell development depends significantly on the grid size. Watt & Sharpe (Reference Watt and Sharpe2005) simulated the 2-D cylindrical detonation with
$E_{a}=24$
by using three levels of grid resolutions of 4, 8 and
$16~\text{points}/L_{1/2}$
(
$\text{pts}/L_{1/2}$
). Their simulations used a third-order monotonic upstream-centred scheme (MUSCL)-type TVD scheme to discretize the convective fluxes, and identified that for
$16~\text{pts}/L_{1/2}$
there is still grid dependency. We simulated the decay process of the strongly overdriven cylindrical detonation with three levels of grid resolution of 10, 20 and
$40~\text{pts}/L_{1/2}$
to verify grid convergence. Results show that solutions obtained with
$20~\text{pts}/L_{1/2}$
converge to those obtained with
$40~\text{pts}/L_{1/2}$
, especially the present simulations used a fifth-order WENO scheme. As such, accuracy of the present WENO scheme is verified, and fifth-order resolution can indeed be achieved for a 3-D detonation (Wang et al.
Reference Wang, Shu, Han and Ning2013). Figure 14 shows the cell evolution at the three levels of the grid resolutions. For
$10~\text{pts}/L_{1/2}$
, the cellular instabilities start to appear at
$x\sim 20$
and is not uniform along the expanding front. For
$20~\text{pts}/L_{1/2}$
and
$40~\text{pts}/L_{1/2}$
, cellular instability all appears at
$x\sim 60$
and the cells are more uniform, demonstrating that the dependence of the cellular instability on the perturbation from the use of Cartesian grid for curved wave is very weak. Overall, the decay rate and key physical feature in the cell shown in figure 14 are not affected significantly by the grid resolution. From figure 15, it is seen that the average detonation velocities of simulations with all three resolutions follow the same global trend of evolution with minor differences. Therefore, grid resolution of
$20~\text{pts}/L_{1/2}$
selected as a compromise to describe the detonation frontal structure is reasonable.
Cellular structures at grid resolutions: (a)
$10~\text{pts}/L_{1/2}$
, (b)
$20~\text{pts}/L_{1/2}$
, (c)
$40~\text{pts}/L_{1/2}$
.
Decay rate of average detonation velocity: (a)
$10~\text{pts}/L_{1/2}$
(red line), (b)
$20~\text{pts}/L_{1/2}$
(blue line), (c)
$40~\text{pts}/L_{1/2}$
(black line).
Detonation velocities during initiation stage for the 2-D cylindrical detonation with and without physical diffusion: Euler (black line), N–S (red line).
Detonation structure with unreacted pocket in the 2-D simulation. (a) Euler (b) NS.
Appendix B. Effects of diffusion on initiation of 2-D cylindrical detonation
Previous studies (Mahmoudi & Mazaheri Reference Mahmoudi and Mazaheri2011, Reference Mahmoudi and Mazaheri2012) have shown that burning of the unreacted pockets is associated with diffusion (Mazaheri et al.
Reference Mazaheri, Mahmoudi and Radulescu2012). It is therefore necessary to assess the role of physical diffusion in the detonation structure, and to further identify whether initiation is significantly affected by diffusion in our study. Consequently, simulations were carried out by solving the reactive N–S equations so that the results can be used to compare with those obtained with the Euler equations. The N–S equations were made dimensionless by the quantities in (10) such that the effects of viscous, thermal and mass diffusion in the dimensionless equations can be reflected respectively by the non-dimensional Reynolds (Re), Prandtl (Pr) and Schmidt (Sc) numbers. Furthermore, we have set
$Pr=1$
,
$Sc=1$
and Lewis numbers
$Le=Pr/Sc=1$
, with
$Re=\sqrt{\unicode[STIX]{x1D6FE}}\unicode[STIX]{x1D70C}_{0}u_{0}L_{1/2}/\unicode[STIX]{x1D707}_{0}=104$
, where
$\unicode[STIX]{x1D707}_{0}$
is assumed to be a constant. Re defined here is only applicable to the overall structure of the detonation, while the local laminar flame consuming the pockets is related to the local Reynold number. From figure 7, it is seen that for
$E_{s}=25\,000$
the 2-D detonation experiences transition from strongly overdriven to quenching, with the formation and then absence of the cell. Therefore, this case is selected to investigate the effect of the physical diffusion on the burning of unreacted pocket and the decay of the average detonation velocity. Figure 16 shows the decay of the velocities obtained by the Euler and N–S equations at
$x\sim 0{-}270$
. Globally, the two velocities have a consistent trend of the decay. Furthermore, it is seen from the enlarged plot in figure 16 that the velocity obtained by the N–S equation is slightly higher than that obtained by the Euler equation. It is noted that the two decay rates are the same at
$x\sim 0{-}50$
, while some minor differences are observed at
$x\sim 50{-}270$
since the cell appears. Indeed, while the front structures are different for the two cases, the global characteristics are qualitatively similar at the same time, as shown in figure 17; for the two cases the detonation structures are globally stretched and the reaction layer almost decouples from the leading shock. Nevertheless, at this time, the detonation obtained by the N–S simulation has globally higher shocked temperature than that by the Euler simulation, which is also observed more obviously at the triple point shown in figure 17.
Figure 17(a) shows the temperature and mass fraction for the Euler simulation. It is seen that the curved leading shock detaches globally from the main reaction zone with corrugated surface and unreacted pockets surrounded by the shear-layer boundaries; see locations 1–4. Small-scale vortices produced by hydrodynamic instability appear at the shear-layer boundaries. At the triple-point location, large amount of the unreacted gas slips into the layer while recedes from the shocked zone, and is still not burnt out and remains in the downstream far from the front. It is noted that the vortices are smeared due to diffusion for the N–S solution, shown in figure 17(b). The unreacted gases entrained into this vortex are also more rapidly reacted, and thereby damped; as a result, the volume of the unreacted pockets becomes globally small. Nevertheless, the modified burning rate of the unreacted gas does not change significantly the process of initiation. Hence, for the 2-D cylindrical detonation that are weakly unstable, diffusion does not change the failure of the initiation or the successful development of the diverging detonation, while slightly affects the detonation structure and the unreacted pockets. A series of N–S computations (Bourlioux & Majda Reference Bourlioux and Majda1992; Oran et al. Reference Oran, Weber, Stefaniw, Lefebvre and Anderson1998; Gamezo et al. Reference Gamezo, Desbordes and Oran1999; Radulescu & Lee Reference Radulescu and Lee2002; Arienti & Shepherd Reference Arienti and Shepherd2005) have been performed to identify the role of diffusion in 2-D planar detonations and found that its importance depends on the regularity associated with the activation energy. Oran et al. (Reference Oran, Weber, Stefaniw, Lefebvre and Anderson1998) suggested that the small-scale structure eliminated in the Euler computation affects neither the overall the detonation structure nor the reactive dynamics, as indicated by Mazaheri et al. (Reference Mazaheri, Mahmoudi and Radulescu2012). Radulescu et al. (Reference Radulescu, Sharpe, Law and Lee2007b ) made the quantitative comparisons between experiments and computation to identify the limitations of inviscid reactive flow in capturing the structure of globally planar detonations with highly irregular cell (Mazaheri et al. Reference Mazaheri, Mahmoudi and Radulescu2012). The role of the diffusion in burning unreacted pocket is found to be relatively weak for the weakly unstable detonation (Mahmoudi & Mazaheri Reference Mahmoudi and Mazaheri2011; Mazaheri et al. Reference Mazaheri, Mahmoudi and Radulescu2012). Consequently, the present simulation results obtained using the inviscid system in the context of weakly unstable detonations are viable.
