Numerical study of acoustic coupling in spinning detonation propagating in a circular tube

Combustion and Flame 160 (2013) 2457–2470

Contents lists available at SciVerse ScienceDirect

Combustion and Flame j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m b u s t fl a m e

Numerical study of acoustic coupling in spinning detonation propagating in a circular tube Yuta Sugiyama ⇑, Akiko Matsuo Department of Mechanical Engineering, Keio University, Yokohama, 3-14-1 Hiyoshi, Kohoku-ku, Yokohama, Kanagawa, Japan

a r t i c l e

i n f o

Article history: Received 22 February 2013 Received in revised form 29 April 2013 Accepted 3 June 2013 Available online 25 June 2013 Keywords: Spinning detonation Propagation mechanism Detonation limit CFD

a b s t r a c t Spinning detonations propagating in a circular tube were numerically investigated with a two-step reaction model by Korobeinikov et al. The time evolutions of the simulation results were utilized to reveal the propagation behavior of single-headed spinning detonation. Three distinct propagation modes, steady, unstable, and pulsating modes, are observed in a circular tube. The track angles on a wall were numerically reproduced with various initial pressures and diameters, and the simulated track angles of steady and unstable modes showed good agreement with those of the previous reports. In the case of steady mode, transverse detonation always couples with an acoustic wave at the contact surface of burned and unburned gas and maintains stable rotation without changing the detonation front structure. The detonation velocity maintains almost a CJ value. We analyze the effect of acoustic coupling in the radial direction using the acoustic theory and the extent of Mach leg. Acoustic theory states that in the radial direction transverse wave and Mach leg can rotate in the circumferential direction when Mach number of unburned gas behind the incident shock wave in the transverse detonation attached coordinate is larger than 1.841. Unstable mode shows periodical change in the shock front structure and repeats decoupling and coupling with transverse detonation and acoustic wave. Spinning detonation maintains its propagation with periodic generation of sub-transverse detonation (new reaction front at transverse wave). Corresponding to its cycle, whisker is periodically generated, and complex Mach interaction periodically appears at shock front. Its velocity history shows the fluctuation whose behavior agrees well with that of rapid fluctuation mode by Lee et al. In the case of pulsating mode, as acoustic coupling between transverse detonation and acoustic wave is not satisfied, shock structure of spinning detonation is disturbed, which causes failure of spinning detonation. Ó 2013 The Combustion Institute. Published by Elsevier Inc. All rights reserved.

1. Introduction Detonations are supersonic flow phenomena with leading shock waves that ignite premixed gas. Many researchers have investigated the structure and properties of detonation using experimental, theoretical, and numerical approaches. Detonation has been experimentally shown to propagate in both longitudinal and transverse directions with unstable three-dimensional structure which has not been able to be clarified due to the difficulty of threedimensional visualization by experimental devices. Shock structure of detonation is composed of incident shock wave, Mach stem, and transverse waves that propagate perpendicular to the shock front. A few modes have been observed in detonation within a circular tube, including spinning (single-headed), two-headed, and multi-headed modes, and they are classified according to the number of transverse waves. Spinning detonation in a circular tube, discovered experimentally in 1926 by Campbell ⇑ Corresponding author. Fax: +81 045 566 1495. E-mail address: [email protected] (Y. Sugiyama).

and Woodhead [1] and followed by a succession of papers by Campbell and Woodhead [2] and Campbell and Finch [3], is observed near the detonation limit and is the lowest mode that has only one transverse wave in the circumferential direction, whereas two-headed detonation has two transverse waves along the circumference direction and one transverse wave in the radial direction. Spinning detonation propagates helically on the wall, and transverse detonation rotates around the tube wall. The mean propagating velocity of spinning detonation is reported to be about 0.8–0.9VCJ (VCJ: CJ detonation velocity). A theoretical study by Fay [4] showed that the ratio of spin pitch to tube diameter, 3.13, was derived from an acoustic theory on the wall, and agrees well with that observed from a trajectory of triple point and transverse detonation [5]. This indicates that coupling with transverse detonation and acoustic wave on the wall is one of important mechanisms for the propagation of spinning detonation. The acoustic theory can explain the property of spinning detonation on the wall but cannot explain its structure inside a whole region. Some researches have been conducted in order to understand the propagation mechanism not only on the wall but also inside a

0010-2180/$ - see front matter Ó 2013 The Combustion Institute. Published by Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.combustflame.2013.06.003

2458

Y. Sugiyama, A. Matsuo / Combustion and Flame 160 (2013) 2457–2470 Table 1 Chemical parameters of two-step reaction model in the present study. Q R E1/R E2/R k1 k2

2.33 MJ/kg 397.6 J/kg/K 9850 K 2000 K 3.0  108 m3/kg/s 4.185  105 m4/N2/s

Fig. 1. Smoked track image in the case where initial pressure P1 and channel width L are 1.0 atm and 30 mm, respectively.

Table 2 The relation between initial pressure P1, cell width k by present study, and cell width k from Ref. [23] and empirical critical diameter Dcr estimated by k/p. Cell width by present study k (mm)

Cell width from Ref. [23] k (mm)

Empirical critical diameter Dcr (mm)

1.0 0.6 0.2

1.62 3.33 10.2

10 18 43

0.516 1.06 3.25

whole region. Schott [5] tried to understand the shock structure of spinning detonation, and concluded that the shock front contains a complex Mach interaction. Voitekhovskii et al. [6] measured the Mach configuration by examining smoked disks attached to the end plate of the detonation tube. They experimentally observed that shock front consists of a ‘‘leg’’ and one or two ‘‘whiskers.’’ They also used the term ‘‘leg,’’ as in ‘‘Mach leg.’’ Topchian and Ul’yanitskii [7] investigated the instability of the spinning detonation and found three different pitch modes: stable pitch, periodical unstable pitch, and pitch covered with a cellular pattern. Over two decades, many researchers investigated twodimensional detonation and showed remarkable insight of the propagation behavior and physics. However, detonation is essentially a three-dimensional phenomenon, and three-dimensional calculations should be performed. For detonations in a tube, the initial and boundary conditions play important roles for the propagation of detonation waves and for the limits where detonation waves fail. For small tubes, the boundary layer introduces the heat and momentum losses that cause the detonation velocity deficits. Furthermore, the detonation wave fails to propagate when tube diameter is below a critical value. Diameter is an important parameter to understand the effect of the boundary layer. Theories for detonation limit with these losses have been developed. Zel’dvich included heat and momentum loss terms in the equations for one-dimensional detonation and showed that detonation velocity depends on losses [8] and diameter. Detonation limit should be

Fig. 2. Time variations of density distribution in two-dimensional periodic channel. (a) Initial condition, (b–d) propagation of one-sided detonation, and (e) cellular detonation with two triple points shifted from one-sided detonation.

3.5

D* = Dcr*/2 = λ/2π Steady mode Unstable mode Pulsating mode

3.0

Diameter D [mm]

Initial pressure P1 (atm)

2.5 2.0 1.5 1.0 0.5 0.0 0.0

0.2

0.4

0.6

0.8

1.0

1.2

Initial pressure P1 [atm] Fig. 3. Calculation conditions and propagation modes in a plane of initial pressure P1 and diameter D; steady (s), unstable (4), and pulsating mode (), for spinning detonation in a circular tube. Red circles ( ) show D Dcr =2 ¼ k=2p (k: simulated cell width).

2459

Y. Sugiyama, A. Matsuo / Combustion and Flame 160 (2013) 2457–2470

Contact surface

Mach leg

Acoustic wave Transverse detonation

(a)

(b)

(c)

Fig. 4. (a) Shock front, (b) density distribution on the wall, and (c) exothermic reaction variable b on the wall in the case of steady mode.

Table 3 Calculation conditions and propagation modes; steady (s), unstable (4), and pulsating (), for detonations in a circular tube. Propagation mode

0.20 0.26 0.30 0.40 0.60

  s s s

0.6 0.6 0.6 0.6

0.51 0.60 0.75 0.90

  s s

0.2 0.2 0.2 0.2

1.50 1.95 2.25 3.00

 4 4 s

considered with the boundary conditions and instabilities of detonation. However, in numerical simulations, conducting three-dimensional calculations with full Navier–Stokes equations is difficult, because of high grid resolution for the boundary layer near the wall. Since calculation with Euler equations can artificially eliminate the loss of heat and momentum, the effect of initial conditions, tube dimension, and geometry could be revealed in the simulations separately. Three-dimensional calculation requires a large number of grid points. Therefore, very powerful computer hardware is required, and some researchers have conducted such calculation. Williams et al. [9] studied the two-headed mode in a square tube, and Deledicque and Papalexandris [10] studied multi-headed mode in a rectangular tube using a one-step chemical reaction model of Arrhenius form. Eto et al. [11] and Tsuboi et al. [12] investigated the detailed shock structures in a rectangular tube of single- and two-headed modes using a detailed reaction model. Tsuboi et al. [13,14] and Virot et al. [15] investigated spinning detonation in a circular tube, and their simulated results agreed well with experimental data. Kasimov and Stewart [16] numerically investigated hydrodynamic instability of spinning detonations using a three-dimensional linear perturbation to find the heat release and activation energy dependence on the spinning mode. The present authors [17] studied the influence of activation energy by one-step reaction model and showed that an increase in activation energy renders the simulated spinning detonation unstable. A few study [14] reported shock structure and propagation mechanism of spinning detonation not only on the wall but also inside a whole region. This paper clarifies the propagation behavior and stability of spinning detonation in a circular tube by three-dimensional numerical simulations, and the effect of initial pressure and diameter with Euler equations. The previous papers introduced here indicate that a coupling with transverse detonation and acoustic wave is one of important mechanisms

Fig. 5. Maximum pressure history on the wall in the case of steady mode.

10.0

Maximum pressure of transverse detonation P/PvN

Diameter D (mm)

1.0 1.0 1.0 1.0 1.0

8.0 6.0 4.0 2.0 0.0 15

16

17

18

19

20

Time [μs]

(a) 2.0

Normalized detonation velocity V/VCJ

Initial pressure P1 (atm)

1.5

1.0

0.5

0.0

12

14

16 18 Time [μs]

20

(b) Fig. 6. (a) Maximum pressure history of transverse detonation on the wall, and (b) time history of detonation velocity.

for the propagation of spinning detonation. Therefore, detailed discussion is carried out to explain the unsteady propagation behavior and the shock structure of spinning detonation with time evolutions of the simulated results not only on the wall but also inside

2460

Y. Sugiyama, A. Matsuo / Combustion and Flame 160 (2013) 2457–2470

  da 1 E1 ¼  ¼ k1 q exp  tig dt RT

a whole region from the viewpoint of acoustic coupling between a transverse detonation and an acoustic wave.

xa 

2. Physical model and numerical method

db xb  ¼ dt

α

β1

M1

unburned gas region

(track angle)

contact surface between burned and unburned gas

incident shock wave

θ1

M2

M2 0

acoustic wave

β3

h  E  i 2 k2 P2 b2 exp  RT  ð1  bÞ2 exp  E2RTþQ ða 6 0Þ

(

ða > 0Þ ð2Þ

0

As discretization methods, Yee’s Non-MUSCL Type 2nd-Order Upwind Scheme [20] is used for the spatial integration, and the Point-Implicit Method that treats only the source term implicitly is used for the time integration. Grid resolution is defined as the number of grid points in induction length Lind calculated by onedimensional steady solution. We conducted a grid convergence study in two-dimensional channel calculation and obtained converged cell widths described in Section 3.1 when 17 or more grid points in induction length are set at each initial pressure (P1 = 0.2, 0.6 and 1.0 atm). Therefore, in all two- and three-dimensional calculations, 17 grid points in induction reaction length Lind are set in all directions. The computational grid for channel flow is an orthogonal system. Channel widths L depend on the initial pressure and are approximately L = 200Lind (L = 30 mm at P1 = 1 atm, L = 45 mm at P1 = 0.6 atm, and L = 112.5 mm at P1 = 0.2 atm) in Section 3.1 and L = 5.9 Lind (0.94 mm at P1 = 1 atm) in Section 3.2. For two-dimensional calculations of channel flow, adiabatic and slip wall in the case of 200Lind and periodic boundaries in the case of 5.9Lind are adapted for the lower and upper boundary conditions, respectively.

100

Shock angle β1 [deg.]

The governing equations are the compressible and reactive twoand three-dimensional Euler equations. The fluid is modeled as an ideal gas with constant specific heat ratio, and all diffusive effects are neglected. The Zel’dvich–von Neumann–Doring (ZND) wave structure model is useful for understanding a CJ detonation. The ZND model consists of a shock discontinuity, followed by a zone of chemical reaction after a certain ignition delay. A two-step reaction model by Korobeinikov et al. [18] is used for chemical reaction, instead of all of the elementary chemical reactions occurring behind a leading shock wave. This simplified model represents the reaction mechanism with two phases, induction and exothermic periods, whose reaction rates are shown in Eqs. (1) and (2), based on the ZND model. In Eqs. (1) and (2), k1 and k2 are reaction rate constants, E1 and E2 are activation energies, and Q and R are heat release and gas constant, respectively. The parameters of this model in the present work are the same as those of Inaba [19] and listed in Table 1. Premixed gas is modeled as stoichiometric hydrogen–air. P, q, and T are pressure, density, and temperature, respectively. In this model, induction progress variable a changes from 1 to 0 (exothermic progress variable b is constant at 1) in the induction period, and exothermic progress variable b changes from 1 to beq (a is constant at 0) in the exothermic period. Initial pressure is chosen as a parameter (P1 = 0.2, 0.6, and 1.0 atm), and initial temperature is fixed at 293 K.

ð1Þ

80

60 33º

40

20

Θ

Incident shock wave

transverse detonation

0 -180

Mach stem

180

Fig. 8. Relation between shock angle b1 and detonation coordinate H from triple point in the circumferential direction on the wall.

60

7

55

6.5

50

6

45 5.5 40 5

35 30

Fig. 7. (a) Schematic diagram of shock structure of spinning detonation with acoustic coupling and definitions of shock angle b, deflection angle h, and detonation coordinate H from triple point in the circumferential direction, and (b) density distribution near triple point on the wall with streamline. M1 and M2 denote Mach numbers of incoming flow before and after incident shock wave in transverse detonation attached coordinate, respectively.

90

0

0.2

0.4

0.6

0.8

1

4.5

Mach number of premixed gas into incident shock wave M1

(b)

0

Angle from triple point Θ [deg.]

Incident shock angle β1 [deg.]

(a)

-90

Mach stem

Distance from tube axis r/R Fig. 9. Relation between distance from tube axis r/R, incident shock angle b1, and Mach number into incident shock wave M1 in the transverse detonation attached coordinate.

Y. Sugiyama, A. Matsuo / Combustion and Flame 160 (2013) 2457–2470

Surface A

Surface A

(a) r/R = 0.2 Unburned gas pocket

(b) r/R = 0.4 Unburned gas pocket

(c) r/R = 0.6 Unburned gas pocket

(d) r/R = 0.8

(e) r/R = 1 Fig. 10. Circumferential distributions of density and reaction progress variables at various radii.

2461

2462

Y. Sugiyama, A. Matsuo / Combustion and Flame 160 (2013) 2457–2470

transverse detonation transverse wave

(a)

reaction front by transverse detonation

(b)

Fig. 11. Distributions of (a) density and (b) mass fraction of reactant at surface ‘‘A’’ described in Fig. 10. Radii r/R of dotted and dashed circles are 0.40 and 0.56, respectively.

The computational grids for three-dimensional detonation are a cylindrical system, whose diameters D are varied between 0.2 and 3.0 mm. The present computational grid of a circular tube has a singular point at a center, where physical values are an average around it. The wall boundary in a circular tube has adiabatic and slip conditions. In all two- and three-dimensional calculations, the premixed gas velocity of incoming flow is 2000 m/s, which is slightly (about 3%) overdriven with respect to CJ velocity in the present chemical parameters. The axial length in the computational grid is more than 100Lind, in order to avoid disturbance from the outflow boundary proposed by Gamezo et al. [21]. The results of one-dimensional steady simulation are used as an initial condition. Sheets of two- and three-dimensional unburned gas mixture behind the detonation front are artificially added in order to create initial disturbances.

3. Results and discussions

Symbols with and without asterisk () indicate the data from the present study, and the previous experimental data [23]. Figure 1 shows the smoked track images in the case of P1 = 1 atm and L = 30 mm. Since pressure near a triple point shows a higher value in the calculation region, the black tracks give trajectories of the triple points that form cellular structures. Table 2 shows the relation between initial pressure P1, cell width k by present study, and cell width k from Ref. [23] and empirical critical diameter Dcr estimated by k/p. Cell widths k are obtained by grid convergence study. Previous paper [24] showed that a cell width by numerical simulation tends to be smaller than that by an experiment even if a detailed chemical reaction model is adopted. Cell widths listed in Table 2 by the present study are also a factor of 4–6 times smaller than those k from Ref. [23], but they are proper values in our numerical method. It is because that, in our paper, spinning detonation observes at conditions of pD  k (k; cell width, D; diameter) as well as those by a previous experiment. Therefore, the present numerical simulations using a simplified chemical reaction model are properly conducted for spinning detonation in a circular tube.

3.1. Cell width k A previous paper [22] showed the existence of critical diameter, which is empirically obtained as pDcr = k (k: cell width). In order to discuss the effect of diameter, cell width k and empirical critical diameter Dcr are estimated by two-dimensional calculations.

Mach number of premixed gas behind incident shock wave M2

3.5

3

2.5

1.841

2

1.5

1

0

0.2

0.4

0.6

0.8

1

Distance from tube axis r/R Fig. 12. Relation between Mach number of premixed gas behind incident shock wave M2 in the transverse detonation attached coordinate and distance front tube axis r/R.

3.2. Acoustic coupling for the propagation of two-dimensional detonation Before we investigate the propagation behavior of spinning detonation, we show that acoustic coupling between transverse wave and acoustic wave behind it plays an important role for the propagation of detonation. In this section, two-dimensional calculation is conducted. First, we calculate the spinning detonation in the case where diameter D and initial pressure P1 are 0.3 mm and 1.0 atm. In this condition, spinning detonation maintains its stable propagation. Using the physical values on the wall, two-dimensional calculation with periodic condition at upper and lower boundaries is conducted in order to understand the importance of acoustic coupling. This simulation shows the propagation of one-sided detonation with one transverse wave in a two-dimensional periodic channel whose width is 0.94 mm. Figure 2 shows the time variations of density distribution of two-dimensional calculation. Figure 2a shows the initial condition from the result on the wall of spinning detonation in a circular tube at initial pressure P1 = 1.0 atm and diameter D = 0.3 mm. Time intervals between Fig. 2b–d are one rotation of the acoustic wave in the periodic channel, and one-sided detonation propagates in a periodic channel. Red and black arrows in Fig. 2b to d denote the positions of the transverse wave and the acoustic wave behind it, respectively. After that, cellular detonation with two triple points appears and maintains stable propagation as shown in Fig. 2e. Indeed, the

2463

Y. Sugiyama, A. Matsuo / Combustion and Flame 160 (2013) 2457–2470

Rear whisker

Front whisker

Reaction front attachment with rear whisker

Precedence of reaction front near wall

(D) 128.31 μs (A) 126.96 μs New whisker

Reaction front attachment with new whisker

Leading shock front pushed by reflected wave

(B) 127.44 μs

(E) 128.50 μs

New whisker

Reaction front attachment with new whisker

(C) 127.83 μs

(a)

(F) 128.69 μs

(b)

Fig. 13. Time evolutions A–F of (a) the shock front from the front side, (b) the surface of induction valuable a = 0 from the front side, (c) density on the wall, (d) absolute gradient of density on the wall, and (e) exothermic reaction variable on the wall.

‘‘acoustic wave’’ far behind triple point is just a weak shock wave. We define ‘‘acoustic wave’’ propagating at burned gas region as a wave with nearly the same of sonic velocity estimated by acoustic theory as shown below. As a state at the burned gas region is assumed as a CJ state, sonic velocity is calculated as 2 aCJ c cMCJ þ 1 c ¼  cþ1 V CJ c þ 1 cM 2CJ

ð3Þ

where MCJ and c are CJ Mach number and specific heat ratio, respectively. After restarting the calculation of two-dimensional periodic channel from the solution on the wall of the stable spinning detonation in a circular tube, a long acoustic wave appears in Fig. 2b. The positions of the transverse wave and the acoustic wave are in-phase in Fig. 2b but shift to out-of-phase in Fig. 2c and d because the velocity of the transverse wave is faster than that of the acoustic wave. Acoustic theory states that the acoustic wave behind the transverse wave propagates at sonic velocity in CJ state in a twodimensional channel. In the simulated case, the averaged speed of the acoustic wave between Fig. 2b and d is 1197 m/s and agrees well with sonic velocity in CJ state (1142 m/s from Eq. (3)). Averaged velocity of the transverse wave between Fig. 2b and d is 1589 m/s, which is much larger than that of the acoustic wave. Therefore, the transverse wave is located ahead of the acoustic

(a)

(b) Fig. 13 (continued)

wave in Fig. 2c and d. Acoustic coupling considers acoustic vibrations in detonation products driven by disturbances on the detonation front. This assumes that there is a finite-amplitude transverse wave attached to an acoustic wave propagating in the detonation products at the CJ state. In the present case, due to the different velocities of transverse and acoustic waves, acoustic coupling between them is not satisfied for the propagation of one-sided detonation. Shock front structure is disturbed and, one-sided detonation cannot steadily maintain its propagation. As a result, a new disturbance propagating in the opposite direction is generated and evolves into a new triple point as shown in Fig. 2e. After that, an averaged velocity of transverse wave (1119 m/s) in the cellular detonation agrees with that of acoustic waves (1142 m/s), and the cellular detonation maintain its propagation. This result indicates that the acoustic coupling between a transverse wave and an acoustic wave is important mechanism to maintain propagation of a detonation. The above-mentioned simulations imply that the propagation of spinning detonation should be discussed while focusing on the acoustic coupling between the transverse wave and acoustic wave. 3.3. Propagation behaviors of spinning detonations in a circular tube Table 3 and Fig. 3 show the calculation conditions and propagation modes: steady (s), unstable (4), and pulsating mode (), for

2464

Y. Sugiyama, A. Matsuo / Combustion and Flame 160 (2013) 2457–2470

Front whisker on the wall

1 2 (A) 126.96 μs

detach 1’ Unburned gas pocket

sub- transverse detonation (B) 127.44 μs

detach 1’ sub-transverse detonation

Unburned gas pocket (C) 127.83 μs

(c)

(d)

(e)

Fig. 13 (continued)

spinning detonation in a circular tube. Red circles ( ) show D ¼ Dcr =2 ¼ k =2p. As diameter decreases, spinning detonation tends to be unstable from steady mode to pulsating mode, as shown in Fig. 3. Symbols with asterisk () indicate the data from the present study. In our simulations, the simulated critical diameter D of spinning detonation is around half the empirical critical diameter Dcr =2 ¼ k =p. Our calculations show three propagation modes; steady mode, unstable mode, and pulsating mode, and their propagation mechanisms are discussed from the viewpoint of acoustic coupling between transverse and acoustic waves in Sections 3.3.1, 3.3.2, and 3.3.3, respectively. 3.3.1. Steady mode This section shows the propagation behavior of steady mode. Since propagation behavior is essentially the same as in other observed steady modes, one condition (P1 = 1.0 atm, D = 0.40 mm) is used to discuss the shock structure of steady mode. Figure 4 shows (a) shock front, (b) density distribution on the wall, and (c) exothermic reaction variable b on the wall of steady mode. Fay’s acoustic theory [4] states that sonic velocity on the wall of circular tube is estimated as 2124 m/s by

V acou ¼ 1:841a

ð4Þ

Here, 1.841 is the first zero of the first derivative of the Bessel function J1. As shown in Fig. 4b, the acoustic wave is rotating in the burned gas region. Therefore, we apply sonic velocity Vacou of 2104 m/s in Eq. (4) as that at CJ state aCJ of 1142 m/s calculated from Eq. (3). We define the wave propagating with the equivalent velocity of Eq. (4) in the burned gas region as ‘‘acoustic wave’’ as shown in Fig. 4b. Our simulations in steady mode show that velocities of acoustic wave are calculated between 0.93Vacou and Vacou at various initial pressures and diameters, and agree with that predicted by Eq. (4) of Fay’s acoustic theory. A white line in Fig. 4c indicates the

shock front on the wall. In the case of steady mode, shock front structure and physical distributions are independent of time in a transverse detonation attached coordinate, and Fig. 4 shows representative images of steady mode. A Mach leg, which stands perpendicular to the wall, exists on the shock front, and no whisker from the edge of Mach leg appears in Fig. 4a. Density distribution shows transverse detonation on the wall, which completes the reaction of the induction zone behind the shock wave in Fig. 4c. In the case of steady mode, the acoustic wave, connecting to the transverse detonation at the contact surface between burned and unburned gases, always maintains rotation with transverse detonation in phase, which indicates that spinning detonation maintains stable propagation when acoustic coupling between transverse detonation and acoustic wave is always satisfied. Figure 5 shows the maximum pressure history on the wall imitating the soot track images in an experiment. The black belt represents the trajectory of transverse detonation on maximum pressure history, and is of constant width, which also indicates that detonation structure maintains the same shape. Track angle a is defined as arctangent of the pitch of the spin divided by circumferential length. The track angle in the experimental observation is about 45° under various experimental conditions, and this value is also derived by an acoustic theory as shown in the following equation

tan a ¼ 1:841

aCJ V CJ

ð5Þ

where aCJ and VCJ are sonic velocity at CJ state and CJ velocity, respectively. In our calculations, the simulated track angle in Fig. 5 is 46.9°, which agrees well with the theoretical value of 47.0° calculated by Eq. (5). Figure 6 shows (a) the maximum pressure history of transverse detonation on the wall and (b) time history of detonation velocity. Pressures and velocities are normalized with respect to von Neumann spike of CJ condition

Y. Sugiyama, A. Matsuo / Combustion and Flame 160 (2013) 2457–2470

2465

Penetration and reflection at contact surface

1’ Penetrated wave

2’ Reflected wave (D) 128.31 μs

Generation of a new whisker Penetrated wave

Reflected wave (E) 128.50 μs

Penetrated wave

Reflected wave (F) 128.69 μs

(c)

(d)

(e)

Fig. 13 (continued)

of one-dimensional steady simulation and CJ velocity, respectively. When we get some physical quantities such as the detonation velocity, we use discrete data at grid points such as displacement of shock front in a small time period. In this case, error may be observed, but is limited in the order of distance between several grid points. The black line in Fig. 6b indicates CJ velocity. Figure 6a shows that the maximum pressure at transverse detonation is constant around 4.8PvN. Longitudinal velocity of steady mode is always constant for a given CJ detonation velocity. We concluded that the spinning detonation in steady mode shows very stable propagation. Figure 7 shows (a) a schematic picture of shock structure of spinning detonation with acoustic coupling and definition of shock angle b, deflection angle h, and detonation coordinate H from triple point in a circumferential direction on the wall, and (b) density distribution near triple point on the wall with streamline. M1 and M2 denote Mach numbers of incoming flow before and after incident shock wave in transverse detonation attached coordinate, respectively. Huang et al. [25] conducted some experiments with two mixtures (H2– O2–Ar and C2H2–O2–Ar) and tried to analyze the shock structure of spinning detonation on the wall. Their experiments showed track angle a and distribution of shock angle b1. At Mach stem, shock angle b1 becomes 90°, gradually decreases far away from the triple point, and approaches around 30° at incident shock wave. From experimental data [25], they investigated the shock structure of spinning detonation, which agrees with the analyzed data using the Rankine–Hugoniot relation, shock angle b1 and deflection angle h1. We also obtain the relation between shock

angle b1 and angle from triple point H as shown in Fig. 8. The present simulated data also show that shock angle b1 at Mach stem (H = +0° in Fig. 7) is 90° and approaches 33° at incident shock wave (H = 0° in Fig. 7) in Fig. 8. This analysis will be utilized in Section 3.3.3 from the viewpoint of the reason why transverse detonation fails in the case of pulsating mode. We analyze the effect of acoustic coupling in the radial direction using the acoustic theory and the extent of Mach leg. Ul’yaniskii [26] carried out experimental measurements and proposed an extent of the Mach leg of 0.5R. Voitsekhovskii et al. [6] reported that the extent of the Mach leg was 0.6R and that his experiments show a uniform spinning detonation. Tsuboi and Hayashi [13] showed that the extent of the Mach leg changes between 0.2R and 0.5R. Here, we apply the transverse detonation attached coordinate as shown in Fig. 7a. In order to estimate Mach number M2 behind an incident shock wave, we get Mach number M1 incoming flow before the incident shock wave in transverse detonation attached coordinate and incident shock angle b1 not only on the wall but also inside a whole region. We assume that acoustic coupling between transverse detonation and acoustic wave on the wall is satisfied and that their velocities on the wall are calculated by 1.841aCJ of Fay’s acoustic theory. In the shock front described in Fig. 4a, the velocity Vtran of the transverse wave and Mach leg in the circumferential direction at radius r is calculated by Eq. (6), because it varies linearly with radius r.

V tran ¼ 1:841aCJ

r R

ð6Þ

Y. Sugiyama, A. Matsuo / Combustion and Flame 160 (2013) 2457–2470

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  r 2 V ¼ V 2CJ þ 1:841aCJ R R r 

ð7Þ

Incident shock angle b1 is described as discrete data at H = 0° in 20 grid points of radii 0.05 6 r/R 6 1 since they are obtained by the simulated data with the Rankine–Hugoniot relation of simulated pressure jump by the incident shock wave and Mach number of incoming premixed gas M1. Figure 9 shows the relation between distance from tube axis r/R, incident shock angle b1 at H = 0°, and Mach number into incident shock wave M1 = V(r/R)/a1 in Fig. 7. Here, a1 is sonic velocity of incoming flow. Mach number into incident shock wave M1 is calculated from Eq. (7), and therefore data is described as a curve. Since the velocity of premixed gas increases with increment of r/R as described in Eq. (7), Mach number of incoming premixed gas M1 becomes large at higher r/R. Incident shock angle b1 decreases with increasing r/R. Figure 10 shows circumferential distributions of density and reaction progress variables at various radii r/R. As r/R becomes small, transverse detonation decays and an unburned gas pocket behind the transverse wave is observed at r/R = 0.2, 0.4 and 0.6. Chemical reaction of the unburned gas is induced by auto-ignition behind the leading shock wave. Figure 11 shows distributions of (a) density and (b) reaction progress variables at surface ‘‘A’’ described in Fig. 10. White dotted and dashed lines denote r/R = 0.40 and 0.56, respectively. In Fig. 11, transverse detonation rotates in the circumferential direction at 0.56 6 r/R 6 1 because the reaction front is attached to it. In contrast, at 0.40 6 r/R 6 0.56, a transverse wave, which does not induces instantaneous chemical reaction, is observed. When r/ R is smaller than approximately 0.40, no transverse wave is observed in Fig. 11. Mach leg of shock front in Fig. 4 is a line of triple point that is an intersection of Mach stem, incident shock wave, and transverse wave. Therefore, the extent of the Mach leg corresponds to that of the transverse wave at r/R = 0.4 as shown in Fig. 11 in the present study. Generally, shock wave can be generated only when flow is supersonic (M P 1), because the effective region is limited, and disturbance can be accumulated. Therefore, sonic condition of M = 1 is an important parameter to discuss the existence of a shock wave. In the case of disturbance propagating in the circumferential direction in a circular tube, acoustic theory shows that it rotates at the velocity of 1.841a (a; local sonic velocity). This indicates that a shock wave rotating in the circumferential direction of a circular tube can be generated only when the Mach number is larger than 1.841. Figure 12 shows Mach number of premixed gas behind incident shock wave M2 in the transverse detonation attached coordinate in Fig. 7a. A red line indicates M2 = 1.841 and intersects with the simulated Mach number M2 at r/R = 0.4. As shown in Fig. 12, a transverse wave appears only when r/R > 0.40 and M2 > 1.841. This indicates that the critical Mach number M2 ¼ 1:841 behind the incident shock wave in the transverse detonation attached coordinate is an important parameter to discuss the extent of Mach leg and shock front structure of spinning detonation in a circular tube. The extent of Mach leg in the present study (r/R = 0.4) is comparable to the previous paper by Ul’yaniskii [26] (r/R = 0.5), Voitsekhovskii et al. [6] (r/ R = 0.6) and Tsuboi and Hayashi [13] (0.2 6 r/R 6 0.5). Our analysis shows that acoustic theory determines the shock structure of spinning detonation inside a whole region. Further investigations, for example, by changing CJ Mach number M1 and specific heat ratio c, can assist our results.

3.3.2. Unstable mode Unstable mode is obtained when the diameter is smaller than that of steady mode, as shown (4) in Fig. 3 and Table 3. The unstable mode in the present study is equivalent to the previous numerical results by Tsuboi et al. [13,14] because of its periodicity. Here, we discuss the propagation mechanism from the viewpoint of acoustic coupling between transverse and acoustic waves. This section shows the propagation behavior of unstable mode, using the shock structure under one set of conditions (P1 = 0.2 atm and D = 2.25 mm). As unstable mode has periodicity during propagation of transverse detonation, one cycle of unstable mode is used to explain the propagation behavior. Figure 13 is time evolutions A–F of (a) the shock front from the front side, (b) the surface of induction valuable a = 0 from the front side, (c) density on the wall, (d) absolute gradient of density on the wall, and (e) exothermic reaction variable on the wall. A sequence in Fig. 13 schematically and quantitatively shows the generation and decay of the complex Mach interaction and whiskers. The images of ‘A’ are almost the

Fig. 14. Maximum pressure history on the wall in the case of unstable mode.

Maximum pressure of transverse detonation P/PvN

where aCJ is the sonic velocity at CJ state, and R is radius of a circular tube. Therefore, the velocity of premixed gas entering into the incident shock wave in the transverse detonation attached coordinate is obtained by velocity composition of VCJ and Vtran, and is a function of r/R, as

10

8

6

4

2

0

40

45

50

55

Time [μs]

(a) 2

Normalized detonation velocity V/VCJ

2466

1.5

1

0.5

0

25

30

35

40

45

50

55

Time [μs]

(b) Fig. 15. (a) Maximum pressure history at transverse detonation on the wall, and (b) time history of detonation velocity.

2467

Y. Sugiyama, A. Matsuo / Combustion and Flame 160 (2013) 2457–2470

same as the images of ‘F’, and therefore the duration between ‘A’ and ‘F’ is one cycle of unstable mode. A series of events is explained as follows with denotation in Fig. 13. (A) One Mach leg and two (front and rear in (a)) whiskers and complex Mach interaction are observed on the shock front. The front whisker is observed in (c) and (d) on the wall. Transverse detonation completes the reaction of premixed gas behind the shock wave in exothermic reaction variable of (e). Two acoustic waves 1 and 2 exist behind the transverse detonation, and they connect at the contact surface between burned and unburned gas in density distribution of (c). (B) Front whisker disappears because of the collision with the wall in (a). Triple point collides with the front whisker on the wall. This induces generation of a sub-transverse detonation (new reaction front at transverse wave) that will propagate to the unburned gas pocket in (e). Acoustic waves 1 and 2 are connected as shown in 1’ and are detached from transverse detonation at the contact surface between burned and unburned gas in ‘detach’ of (c). Between (A) to (C), velocity of acoustic wave 10 is 2104 m/s, which agrees well with that of 2124 m/s predicted by Fay’s acoustic theory, whereas velocity of transverse detonation is 1951 m/s. Therefore, an acoustic wave exists ahead of transverse detonation on the wall. This indicates that acoustic coupling between them is not satisfied. (C) The shock front structure is similar to that of ‘B’ from the front side (a) and (b). The sub-transverse detonation is propagating to the unburned gas pocket described in (c) and (e). The transverse detonation and Mach leg on the wall are instantaneously accelerated from 1951 m/s to 2103 m/s in the circumferential direction. (D) After completion of exothermicity of

the unburned gas pocket, precedence of reaction front is observed near the wall in (b). Sub-transverse detonation reflects on and penetrates the contact surface described in (c). Reflected wave propagates to the leading shock front, whereas penetrated one becomes an acoustic wave 20 . It propagates with the velocity of 2104 m/s, which agrees well with that predicted by Fay’s theory and couples with the transverse detonation. Acoustic coupling is satisfied between the transverse detonation and the penetrated wave 20 at contact surface of burned and unburned gas in (c). (E) The reflected wave pushes the leading shock front. This makes wavy contact surface in (e). The intersection of the leading shock front and the reflected wave becomes a triple line called as new whisker, and it propagates to the center of a circular tube in (a) and (d). Since a chemical reaction behind the new whisker is prompted, the reaction front attaches to a new whisker near the wall as shown in (b). (F) Shock front structure is the same as that of moment (A). It is observed approximately in a half cycle of one round, and its period is about 1.7 ls under these conditions (P1 = 0.2 atm, D = 2.25 mm) as shown in Fig. 13. Compared with the steady mode, spinning detonation becomes slightly unstable due to the repetition of coupling and decoupling between transverse detonation and acoustic wave. If sub-transverse detonation, generated by successive collisions of a triple point and whisker on the wall, does not appear, unburned gas region enlarges, and chemical reaction is inhibited. As a result, spinning detonation fails, and a detonation becomes pulsating mode. Figure 14 shows maximum pressure history on the wall resembling the soot track images in the experiment. The black belt is the

Entrance of acoustic wave into unburned gas region

(A) 31.66 μs

Entrance of acoustic wave into unburned gas region

(B) 32.87 μs

(C) 33.51 μs

(a)

(b)

(c)

Fig. 16. Time evolutions A–C of (a) density, (b) absolute gradient of density, and (c) exothermic reaction variable on the wall around the moment of the failure of spinning detonation. White lines in (c) denote the shock front.

2468

Y. Sugiyama, A. Matsuo / Combustion and Flame 160 (2013) 2457–2470

trajectory of transverse detonation on maximum pressure history. The black belt repeatedly becomes wide and narrow in a half cycle of one round, which is the same behavior as that of the generation and decay of complex Mach interaction as shown in Fig. 13. Averaged track angle 47.6° agrees well with the theoretical value of 47.0° yielded by Eq. (5). Figure 15 shows (a) the maximum pressure history of transverse detonation on the wall and (b) time history of detonation velocity. Pressures and velocities are normalized with respect to von Neumann spike of CJ condition of one-dimensional steady simulation and CJ velocity, respectively. The black line in Fig. 15b indicates a CJ velocity. Since the collision of a triple point and a whisker and the generation of sub-transverse detonation occur intermittently, the strength of transverse detonation shows the periodical change in Fig. 15a. Its averaged period is 1.7 ls, which agrees well with that of generation and decay of complex Mach interaction as discussed in Fig. 13. As calculation conditions move toward lower initial pressure and smaller diameter, spinning detonations become unstable and show the periodical propagation behavior. Longitudinal velocity of unstable mode in Fig. 15b shows fluctuation around a CJ velocity. Lee et al. [27] reported velocity fluctuation near extinction limit and five modes of propagations: stable detonation mode, rapid fluctuations mode, stuttering mode, galloping mode, and low-velocity stable mode (fast flame). Longitudinal velocity is constant near VCJ in stable detonation mode, fluctuating around VCJ in rapid fluctuation mode, periodically changes from 0.6VCJ to VCJ with repetition of failure and re-ignition in stuttering mode, and periodically changes from 0.4VCJ to 1.5VCJ with periodical local explosion in galloping mode and is constant near 0.5VCJ in low-velocity stable mode. In the case of unstable mode, since the velocity fluctuation appears around CJ detonation velocity, it is classified as rapid fluctuation mode. 3.3.3. Pulsating mode Pulsating mode is obtained as shown () in Table 3 and Fig. 3. See our previous work [28] for detailed discussion on the propagation behavior of pulsating mode. Pulsating mode shows cyclic behavior in the longitudinal direction, including three regions of multi-headed, single-headed spinning, and no detonation. In this paper, we focus on the failure mechanism of single-headed spinning detonation from the viewpoint of acoustic coupling between transverse detonation and acoustic wave behind it, using one set of simulated conditions (P1 = 1.0 atm, D = 0.20 mm). Figure 16 shows the time evolutions A–C of (a) density, (b) absolute gradient of density and, (c) exothermic reaction variable on the wall around the moment of the failure of spinning detonation. White lines in (c) denote the shock front. The velocity of acoustic wave is 2104 m/s, which agrees well with that estimated by Fay’s acoustic theory, whereas the velocity of transverse detonation is 1770 m/s and much smaller than that of acoustic wave. This indicates that acoustic coupling between transverse detonation and acoustic wave is not satisfied. At (A) and (B), acoustic wave is ahead of transverse detonation and enters the unburned gas region behind the incident shock wave. At (C), unlike unstable mode, a sub-transverse detonation is not generated. Chemical reaction behind transverse detonation is inhibited, transverse detonation suddenly decays, and an unburned gas pocket is observed behind the transverse wave. In the case of spinning detonation, transverse detonation has an important role to maintain propagation because it completes the potential exothermicity. However, the transverse detonation cannot maintain stable rotation in the pulsating mode and, eventually, spinning detonation fails. Figure 17 shows (a) schematic picture of shock structure of spinning detonation without acoustic coupling and definitions of shock angle b, deflection angle h, and detonation coordinate H from triple point in the circumferential direction on the wall and (b) density distribution near triple point on the wall with

streamline at 32.87 ls (close-up view of Fig. 16a-B). M1 and M2 denote Mach numbers of incoming flow before and after incident shock wave, respectively, and M3 denotes Mach number after the secondary shock wave in transverse detonation attached coordinate. In the enlarged picture in Fig. 17a, there are defined zone2: a flow after incident shock wave, zone-3: a flow after the secondary shock wave from zone-2, zone-4: flow after transverse detonation from zone-2, zone-5: after transverse detonation from zone-3. Just before spinning detonation fails, the flow is almost steady in the circumferential direction, and streamline is described in Fig. 17b. Premixed gas is deflected by the secondary shock wave. In the case of pulsating mode, since the acoustic wave is faster than transverse detonation, the acoustic wave is located in front of transverse detonation, and enters the unburned gas region behind incident shock wave as shown in Figs. 16 and 17. Since the velocity of the acoustic wave (2104 m/s) is larger than a sonic velocity of unburned gas (a2  800 m/s), the acoustic wave becomes a secondary shock wave that affects discontinuous changes of physical values such as pressure and temperature. Here, we analyze the effect of the secondary shock wave. We make some assumptions below. These are: that track angle a is 47° and Mach number into the leading shock wave is 1.47MCJ (=7.05), that shock angle at incident wave b1 near triple point on the wall (H  0°) is 33° as well as that in Fig. 8a, that transverse detonation angle b4 in Fig. 17a is 90°, that velocity vectors in zones 4 and 5 are headed in the same direction, and that the secondary shock angle b2 is a function in

Incident shock wave M1 β 1 α (track angle) M2

M2

θ2

acoustic wave 0

M2

β4

contact surface between burned and unburned gas

θ1

β2

2

unburned gas region M2

β3

Θ

M3

transverse detonation

3 Mach stem

4

5

contact surface between flows with and without the secondary shock wave

(a) β 2 ~ 37º

(b) Fig. 17. (a) Schematic diagram of shock structure of spinning detonation without acoustic coupling and definitions of shock angle b, deflection angle h, and detonation coordinate H from triple point in a circumferential direction on the wall and (b) density distribution near triple point on the wall with streamline at 32.87 ls (close-up view of Fig. 16a-B). M1 and M2 denote Mach numbers of incoming flow before and after incident shock wave, respectively, and M3 denotes Mach number after the secondary shock wave in transverse detonation attached coordinate.

2469

1.4

1.6

1.2

1.4

Density ρtran/ρ*

Pressure Ptran/P*

Y. Sugiyama, A. Matsuo / Combustion and Flame 160 (2013) 2457–2470

1 0.8 0.6 0.4 0.2 0 15

1 0.8 0.6 0.4

20

25

30

35

0.2 15

40

20

25

30

35

40

Secondary shock wave angle β2 [deg.]

Secondary shock wave angle β2 [deg.]

(a)

(b)

0.9

100

Induction time tind/tind*

Temperature Ttran/T*

1.2

0.8 0.7 0.6 0.5 0.4 0.3 15

20

25

30

35

40

10

1 15

20

25

30

35

40

Secondary shock wave angle β 2 [deg.]

Secondary shock wave angle β2 [deg.]

(c)

(d)

Fig. 18. Relation between the secondary shock wave angle b2 and (a) pressure, (b) density, (c) temperature, and (d) induction time in zone 5. Plots indicate the condition that pressures at zone-4 of P and zone-5 of Ptran are continuous at contact surface in the enlarged picture of Fig. 17.

this analysis. From the above assumptions, the secondary shock wave angle b3, and physical values of pressure Ptran, density qtran, temperature Ttran and induction time tind after the transverse wave and before chemical reaction are estimated in Fig. 18. Here, induction time tind is calculated by zero-dimensional analysis of Eq. (8) with qtran and Ttran after transverse wave.

t ind ¼

1   1 k1 qtran exp  RTEtran

ð8Þ

In Fig. 18, physical values are normalized with respect to those without the secondary shock wave in zone-4 (P = 184 atm, q = 15.4 kg/m3, T = 3070 K and tind ¼ 0:0053 ls), and plots indicate the condition that pressures in zone-4 and zone-5 are continuous at the contact surface described in dashed line of the enlarged picture of Fig. 17. As shown in the streamline of Fig. 17b, the secondary shock wave angle b2  37°, which agrees well with plots in Fig. 18. Therefore, the estimated data in Fig. 18 are applicable to estimate physical values in zone-5 after the secondary shock wave. At all secondary shock wave angles b2, the induction time becomes larger than that in zone-4, and the secondary shock wave has an effect on deceleration of the chemical reaction behind the transverse wave. In the plot in Fig. 18d, induction time becomes 1.90 times larger than that in zone-4. The failure mechanism of spinning detonation is explained as follows. As acoustic coupling is not satisfied, the acoustic wave is located ahead of transverse wave and enters into unburned gas region. It becomes the secondary shock wave at unburned gas region. The shock structure of spinning detonation is disturbed and chemical reaction is decelerated by the secondary shock wave. As a result, transverse detonation fails. This is to say

that changes of physical values and flow angle deflection by the secondary shock wave are very sensitive to the propagation of transverse detonation, and that acoustic coupling should be satisfied in order to maintain the propagation of spinning detonation. 4. Conclusions Spinning detonations in a circular tube at various initial pressures were numerically investigated using three-dimensional Euler equations with a two-step chemical reaction model proposed by Korobeinikov et al. Three distinct propagation modes were obtained: namely, steady mode, unstable mode, and pulsating mode. Steady mode showed stable propagation without change in the shock structure. The maximum pressure history of transverse detonation on the wall and velocity history of detonation remained nearly constant, and a Mach leg always existed on the shock front and rotated at a constant speed. Coupling with transverse detonation and acoustic wave was always satisfied on the wall. We analyzed the effect of acoustic coupling in the radial direction using acoustic theory and the extent of Mach leg. Acoustic theory stated that in the radial direction transverse wave and Mach leg could rotate in the circumferential direction when Mach number behind the incident shock wave in transverse detonation attached coordinate was larger than 1.841. The extent of Mach leg in the present study was comparable to the previous papers. Unstable mode showed periodical change in the shock structure. Coupling and decoupling with transverse detonation and acoustic wave were repeated, inducing the periodical flow field. Sub-transverse detonation was repeatedly generated by the collision of triple point and whisker on the wall.

2470

Y. Sugiyama, A. Matsuo / Combustion and Flame 160 (2013) 2457–2470

Reflected and penetrated waves were generated by collision of sub-transverse detonation and contact surface between burned and unburned gas. The penetrated wave decelerated to the acoustic wave, which coupled with the transverse detonation, whereas the reflected wave intersected with the leading shock front, and as a result, intersection line became a triple line called as a new whisker on a shock front. The velocity history showed fluctuation around VCJ, which qualitatively agreed with that of rapid fluctuation mode by Lee et al. In the case of pulsating mode, spinning detonation could not maintain its propagation due to decoupling of acoustic wave and transverse detonation. The acoustic wave was located ahead of transverse wave and enters into unburned gas region. It becomes the secondary shock wave at unburned gas region. The shock structure of spinning detonation is disturbed and chemical reaction is decelerated by the secondary shock wave. As a result, transverse detonation fails. Acknowledgments This work has been supported in part by Grant in Aid for the Global Center of Excellence Program for ‘‘Center for Education and Research of Symbiotic, Safe and Secure System Design’’ from Ministry of Education, Culture, Sport, and Technology in Japan. References [1] C. Campbell, D.W. Woodhead, J. Chem. Soc. (1926) 3010–3021.

[2] [3] [4] [5] [6]

[7] [8] [9] [10] [11] [12] [13] [14] [15]

[16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28]

C. Campbell, D.W. Woodhead, J. Chem. Soc. (1927) 1572–1578. C. Campbell, A.C. Finch, J. Chem. Soc. (1928) 2094–2106. J.A. Fay, J. Chem. Phys. 20 (1952) 942–950. G.L. Schott, Phys. Fluids 8 (1965) 850–865. B.V. Voitsekhovskii, V.V. Mitrofanov, M.E. Topchian, Izd-vo Sibirsk, Odetl. Adak. Nauk SSR, Novosibirsk (1963) (Translation 1966, Wright-Patterson Air Force Base, Report, FTD-MT-64-527 (AD-633-821)). M.E. Topchian, V.Y. Ul’yanitskii, Acta Astronaut. 3 (1976) 771–779. Y.B. Zel’dvich, NASA Technical Memorandum 1261 (1950). D.N. Williams, L. Bauwens, E.S. Oran, Proc. Combust. Inst. 26 (1996) 2991– 2998. V. Deledicque, M.V. Papalexandris, Combust. Flame 144 (2006) 821–837. K. Eto, N. Tsuboi, A.K. Hayashi, Proc. Combust. Inst. 30 (2005) 1907–1913. N. Tsuboi, M. Asahara, K. Eto, A.K. Hayashi, Shock Waves 18 (2008) 329–344. N. Tsuboi, A.K. Hayashi, Proc. Combust. Inst. 31 (2007) 2389–2396. N. Tsuboi, K. Eto, A.K. Hayashi, Combust. Flame 149 (2007) 144–161. F. Virot, B. Khasainov, D. Desbordes, H. Presles, IN: 21st International Colloquium on the Dynamics of Explosions and Reactive Systems, 2007 on CD-ROM. A.R. Kasimov, D.S. Stewart, J. Fluid Mech. 466 (2002) 179–203. Y. Sugiyama, A. Matsuo, Proc. Combust. Inst. 32 (2009) 2331–2337. V.P. Korobeinikov, V.A. Levin, V.V. Markov, G.G. Chernyi, Astronaut. Acta 17 (1972) 529–537. K. Inaba, PhD Thesis, Keio University, 2004. H.C. Yee, NASA Technical, Memorandum 89464, 1987. V.N. Gamezo, D. Desbordes, E.S. Oran, Shock Waves 9 (1999) 11–17. A.A. Vasil’’ev, Shock Waves 18 (2008) 245–253. A.A. Vasil’ev, J. Propul. Power 22 (6) (2006) 1245–1260. B.D. Taylor, D.A. Kessler, V.N. Gamezo, E.S. Oran, Proc. Combust. Inst. 34 (2) (2013) 2009–2016. Z.W. Huang, M.H. Lefebvre, P.J. Van Tiggelen, Shock Waves 10 (2000) 119–125. V.Y. Ul’yanitskii, Fiz. Goreniya Vzryva 16 (1980) 105–111. J.J. Lee, G. Dupre, R. Knystautas, J.H. Lee, Shock Waves 5 (1995) 175–181. Y. Sugiyama, A. Matsuo, Combust. Flame 159 (2012) 3646–3651.