Detailed structure of spinning detonation in a circular tube

Combustion and Flame 149 (2007) 144–161 www.elsevier.com/locate/combustflame

Detailed structure of spinning detonation in a circular tube N. Tsuboi a,∗ , K. Eto b , A.K. Hayashi b a Space Transportation Engineering Department, Institute of Space and Astronautical Science, Japan Aerospace Exploration Agency, Yoshinodai 3-1-1, Sagamihara, Kanagawa 229-8510, Japan b Department of Mechanical Engineering, Aoyama Gakuin University, Fuchinobe 5-10-1, Sagamihara, Kanagawa 229-8558, Japan

Received 1 December 2005; received in revised form 2 November 2006; accepted 11 December 2006 Available online 15 February 2007

Abstract A single spinning detonation wave propagating in a circular tube, discovered experimentally in 1926, is simulated three-dimensionally with a detailed chemical reaction mechanism. The detonation front obtained numerically rotates periodically with a Mach leg, whiskers, and a transverse detonation. A long pressure trail, which is distributed from the transverse detonation to downstream, was reproduced, clearly showing that the pressure trail also spins synchronously with the transverse detonation. The formation of an unburned gas pocket behind the detonation front was not observed in the present simulations because the rotating transverse detonation completely consumed the unburned gas. The calculated profiles of instantaneous OH mass fraction have a keystone shape behind the detonation front. The numerical results for pitch, track angle, Mach stem angle, and incident shock angle on the tube wall agree well with the experimental results. © 2006 The Combustion Institute. Published by Elsevier Inc. All rights reserved. Keywords: Detonation; Single spin; 3D structure; Numerical simulation; Detailed chemical reaction model

1. Introduction Detonation is a shock-induced combustion wave propagating through a reactive mixture or pure exothermic compound and has been studied from the safety engineering point of view, such as for coal mine explosions, and from the scientific point of view of astrophysics, as in star explosions. The detailed structures and properties of the detonation have been studied using the experimental and numerical methods as well as through a theoretical approach. Experimental studies have been carried out by many re* Corresponding author.

E-mail address: [email protected] (N. Tsuboi).

searchers, who provide useful fundamental and practical data. However, the detailed three-dimensional structure of detonations has remained unclear because of its fast propagation speed and of the difficulty of three-dimensional visualization by experimental measurement devices. Detonations were measured to observe a few modes, for example, spinning, two-headed, and multiheaded modes. Detonations comprise transverse waves, which move along the circumference and a radius. The single spinning detonation has only a transverse wave along the circumference, whereas the twoheaded detonation has two transverse waves along the circumference and one transverse wave along a radius, respectively. The spinning detonation is the lowest mode in the limit mixtures in a circular tube to propagate with a helical track on the wall and to rotate

0010-2180/$ – see front matter © 2006 The Combustion Institute. Published by Elsevier Inc. All rights reserved. doi:10.1016/j.combustflame.2006.12.004

N. Tsuboi et al. / Combustion and Flame 149 (2007) 144–161

around the tube axis. Campbell and Woodhead first observed the reproducible striations produced on the high-speed photographic records by a spinning detonation in a stoichiometric mixture of carbon monoxide and oxygen in 1926 [1–3]. Bone and co-workers systematically investigated detonation phenomena to conclude that spin was connected with the mode of coupling between the leading shock front and the reaction zone [4–6]. A theoretical study was also started in 1946 [7–9] to show that spin detonation occurs via an acoustic wave moving transversely across the main front, which propagates axially along the tube. Fay [8] showed that the ratio of the pitch of the spin to the diameter of the tube was derived by an acoustic theory and that this ratio is 3.13, which is close to the measured value of 3 [3,10]. However, the acoustic theory cannot explain the shock structure of the spinning mode in detail. Duff [11], Schott [10], Voitsekhovskii et al. [12–14], Denisov and Troshin [15], and Huang et al. [16] tried to find the shock structure of the spinning mode and they concluded that the wave front of the spinning mode contains a complex Mach interaction. Topchian and Ul’yanitskii [17] investigated the instability of the spinning detonation to divide into three different types: the stable pitch, the periodical unstable pitch, and the pitch covered with a cellular pattern. Bykovskii et al. [18] experimentally studied a continuous spin detonation in an annular tube with an expanding cross section to find the existence region of the spin mode and the structure of transverse detonation. The effects of the composition of the mixture on spin were described by Gordon et al. [19] and Barthel [20]. Gordon showed the effect of the initial pressure and composition of highly diluted mixtures of hydrogen and oxygen with a monatomic gas to conclude that the phenomenon of single spin for mixtures close to fuel-lean and fuel-rich mixtures was not influenced by the initial pressure of the mixture. Kasimov and Stewart [21] 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. Achasov and Penyazkov [22] also investigated the evolution of the cellular structure of a gaseous detonation, including a spinning detonation in a circular tube, as a function of the initial pressure. Spinning detonation is also observed in a tube of noncircular cross section as well as in a circular tube. Bone et al. carried out experiments in triangular, square, and oblong tubes to measure the spin pitch [6] and Oppenheim and co-workers presented a smoke-film record of spinning detonation in a square cross-section tube [23]. However, experimental observations for spinning detonation are mostly presented for circular tubes.

145

Recent computational fluid dynamics (CFD) has yielded a remarkable insight in these problems. One of the earliest numerical detonation studies was presented in 1978 by Taki and Fujiwara [24]. They simulated two-dimensional gaseous detonation using a two-step reaction model. Oran, Kailasanath, Gamezo, and co-workers [25,26] have studied the detonation structure over two decades. Three-dimensional simulations have been performed by Williams et al., who studied the rectangular mode in a rectangular tube [27,28], and by the present authors [29–31]. Computation with a detailed chemical reaction model requires high resolution, on the order of 1 µm, near the detonation front, and the physical size of the computational domain is on the order of 1 mm for a hydrogen/air gas mixture at a pressure of 0.1 MPa. Simulations for detonation with cells require a very high resolution and a large number of computational grid points in all computational domains. Therefore, simulation of three-dimensional detonation requires very powerful computer hardware and few researchers have conducted it. As a result of this restriction, two-dimensional simulations of the detonation have mainly been carried out to understand its flow structure and properties in detail. Our group has been working on simulation of three-dimensional hydrogen/air detonation with a detailed chemical reaction model in a rectangular tube for about four years and this work has revealed its two modes, namely a rectangular mode in phase and a diagonal mode [29–32]. These modes have been experimentally observed in the past by many researchers. There are two types of rectangular modes: in phase and partially out of phase. These rectangular modes consist of two two-dimensional waves. The threedimensional cell length for these modes is approximately the same as the two-dimensional cell length. The diagonal mode shows three-dimensional diagonal motion of the triple-point lines and cannot be simulated by two-dimensional simulations. On the other hand, good three-dimensional numerical simulation in a circular tube has not been achieved so far, because of the existence of a singular point in a tube center. Washizu and Fujiwara [33] simulated a spinning detonation in a coaxial tube using a twostep chemical reaction model; however, this result does not correspond to a pure spinning detonation in a circular tube. Some authors computed a threedimensional heterogeneous detonation in a cylindrical tube, which revealed a structure made of periodic two-headed detonation [34]; however, a spinning detonation was not observed. The authors succeeded in simulate a gaseous single spinning detonation in a circular tube [35]; however, more computational resolutions are necessary to simulate the spinning detonation.

146

N. Tsuboi et al. / Combustion and Flame 149 (2007) 144–161

The aim of this work is to reproduce the physical phenomenon and to investigate and characterize the detailed structure of the spinning detonation numerically. The three-dimensional numerical approach enables us to provide a clear unsteady flow structure including transverse detonation, complicated shock structure, and a pressure trail.

is given as Cpi /Ri = a1i + a2i T + a3i T 2 + a4i T 3 + a5i T 4 . (6) The coefficients in Eq. (6) are obtained from the data in the JANAF tables [36]. The production rate of each chemical species, ω˙ i , is given by combining the elementary chemical reactions in the kinetic model so that

2. Physical model and numerical method The governing equations are the three-dimensional Euler equations with 9 species (H2 , O2 , H, O, OH, HO2 , H2 O2 , H2 O, and N2 ) and 18 elementary reactions in rectangular coordinates, ∂Q ∂E ∂F ∂G + + + = S, ∂t ∂x ∂y ∂z where ⎛ ρu ⎞ ⎛ ρ ⎞ ⎜ ρu2 + p ⎟ ⎜ ρu ⎟ ⎜ ⎜ ⎟ ⎟ ⎜ ρuv ⎟ ⎜ ρv ⎟ Q=⎜ ⎟, ⎟, E = ⎜ ⎜ ρuw ⎟ ⎜ ρw ⎟ ⎝ ⎝ ⎠ ⎠ e (e + p)u ρi ρi u ⎛ ρv ⎞ ⎜ ρvu ⎟ ⎜ 2 ⎟ ⎜ ρv + p ⎟ F=⎜ ⎟, ⎜ ρvw ⎟ ⎝ ⎠ (e + p)v ρi v ⎛ ρw ⎞ ⎛ 0⎞ ρwu ⎜ ⎜ 0⎟ ⎟ ⎜ ρwv ⎟ ⎜ ⎟ ⎜ ⎜ 0⎟ ⎟ G=⎜ ⎟, S = ⎜ ⎟, 2 ⎜ ρw + p ⎟ ⎜ 0⎟ ⎝ ⎝ ⎠ ⎠ 0 (e + p)w ω˙ i ρi w

Ns  i=1

1 ρi hi − p + ρ u2 + v 2 + w2 , 2

(1)

(2)

(3)

(4)

where Ns is total number of species. The thermal equation of state for a perfect gas is given by p =  Ns i=1 ρi Ri T , where Ri denotes the gas constant for each species. The specific heat at constant pressure for each species, Cpi , which is used to evaluate the total energy through the specific enthalpy T hi = T0

Cpi dT + h0i ,

Nr 

 − ν  νik ik

(5)



k=1

× kf,k

Ns  i=1

where ρ is density, u, v, and w are velocities in the x, y, and z directions, e is total energy, ρi is density of the ith species, p is pressure, and ω˙ i is production rate of the ith species by chemical reaction, respectively. The total energy is defined as e=

ω˙ i = Wi



cχi νik − kb,k

Ns 





cχi νik ,

(7)

i=1

where ν  are the stoichiometric coefficients of the reactants, ν  are the stoichiometric coefficients of the products, Wi is molecular weight, Nr is the total number of elementary reactions, kf is forward specific reaction-rate constant, kb is backward specific reaction-rate constant, and cχ is concentration. These equations are explicitly integrated by the Strang-type fractional-step method. The chemical reaction source terms are treated in a linearly pointimplicit manner in order to avoid a stiff problem. A second-order Harten–Yee non-MUSCL-type TVD scheme is used for the numerical flux in the convective terms [37]. The averaged state on a computational cell boundary is given by the generalized Roe average [38] to evaluate the numerical flux in the convective terms. In the present simulation, the Petersen and Hanson model is used for chemical kinetics to solve detonation problems; this model was proposed by Petersen and Hanson [39] as a new detailed chemical reaction model. This model contains 18 reactions and 9 species, and it is based on the H2 /O2 submechanism of the RAMEC/Gas Research Institute GRIMech 1.2 methane-oxidation mechanism. A feature of this model is the pressure dependence on a forward reaction coefficient included in the collision reaction with a third body; in other words, HO2 and H2 O2 chemistries near the second and third explosion limits are considered because they are necessary for ignition at extremely high pressure but are lacking in certain finite-rate chemical models currently in use. The effects of pressure dependence in the reaction model on the detonation have to be studied as well as the comparison with simplified reaction models. Our previous simulation [40] showed that the present model including the pressure dependence increases, depending on the conditions, detonation cell size more than the other detailed reaction models. The present simulation neglects viscous effects such as turbulence and boundary-layer interactions. Local vortex or local scales may be affected by the

N. Tsuboi et al. / Combustion and Flame 149 (2007) 144–161

Fig. 1. Computational grid system for GS1. A grid point for every four points in the axial and circular directions is plotted. The plotted axial length is one-half.

viscous effects; however, the largest scales generated by the shock wave might not be affected [41]. The computational mesh is a cylindrical system with 1201 (axial direction) × 41 (radial direction) × 213 (circumferential direction) grid points (GS1) and 601 × 41 × 319 (GS2). The computational grid system for GS1 is shown in Fig. 1. The grid sizes are 5 µm in the propagating direction, 10 µm near the wall and 20 µm close to the center in the radial direction, and 15 µm for GS1 or 10 µm for GS2 in the circular direction along the tube wall. The grid size of 5 µm gives a resolution of 32 grid points in the theoretical half-reaction length [42], which equals 160 µm for H2 at atmospheric pressure [43]. Therefore, the computational domains are 6 or 3 mm in length and 1 mm in diameter. The difference between GS1 and GS2 is the circular resolution and the length of the computational domain. The other mesh systems with 601 × 41 × 213 (GS3) and 601×41×151 (GS4) are also used in order to estimate the grid convergence for the spinning detonation. These meshes and boundary conditions at the downstream (leftside boundary) are summarized in Table 1. The present computational domain is small in order to maintain high resolution to reveal the detailed three-dimensional structure of the spinning detonation. As mentioned previously, numerical simulation with a detailed reaction model requires high grid resolution in space. In our study, circular grid resolution as well as axial grid resolution is found to be a very im-

147

portant feature for simulating a spinning detonation. The spinning detonation has a transverse detonation wave rotating along the circular wall. The transverse detonation must also be resolved in the finer circular grid. Our previous two-dimensional study confirmed that the dependence of the computational domain on the detonation structure is negligible when the ratio of the channel width to the length of the propagating direction is about 3–4 times or more. The boundary conditions are as follows: the upstream conditions are a pressure of 0.1 MPa and a temperature of 300 K, and the inflow gas is a stoichiometric H2 /air gas mixture; the wall boundary conditions are adiabatic, slip, and noncatalytic; and the downstream condition is that a nonreflection boundary, proposed by Gamezo et al. [26], is imposed for GS1 or that pressure is fixed at the boundary and other variables are extrapolated for GS2–GS4. The pressure at the downstream boundary equals 1.78 MPa for GS2–GS4, which is the Chapman–Jouguet (CJ) theoretical value. The detonation velocity in the rectangular tube becomes the CJ velocity easily, however, that in the circular tube requires many iterations. It is thought that a long pressure trail interacts with the downstream boundary; therefore the simulation for GS1 with a long tube and with the nonreflecting boundary conditions on the downstream boundary was carried out. The present computational grids have a singular point in the tube center, as everybody experiences, and the physical values at the singular point are averaged ones around it. However, this method is simple and can solve without numerical instability. Different initial conditions for the one-dimensional simulations are chosen in two computational domains. In the domain in the vicinity of the closed end wall, pressure was set at a high value, whereas it was set at a low value in the other domain. The results of the one-dimensional simulations, which are confirmed to be a typical ZND detonation structure, are used as initial condition for the three-dimensional simulations. The sheets of three-dimensional unburned gas mixture behind the detonation front are artificially added in order to create an initial threedimensional disturbance. These sheets are asymmetrically formed in the radial direction. The calculated

Table 1 Computational meshes and boundary conditions GS

1

2

3

4

Grid system Downstream boundary condition rθ on the wall (µm) Length of computational domain (mm) Diameter of computational domain (mm)

1201 × 41 × 213 No reflection 15 6 1

601 × 41 × 319 Fixed pressure 10 3 1

601 × 41 × 213 Fixed pressure 15 3 1

601 × 41 × 151 Fixed pressure 20 3 1

148

N. Tsuboi et al. / Combustion and Flame 149 (2007) 144–161

mean detonation velocity for the one-dimensional simulations is approximately 1980 m/s. The parallel implementation was carried out using Open MP. The present simulations require 67,000 iterations for GS1 and 83,000 for GS2 (about 300 h on 8 processors of a NEC SX-6) to obtain a steady spinning detonation.

3. Results and discussion The results for GS1 are mainly about the detonation velocity and the shock wave angle, because the calculated average detonation velocity becomes slower than the CJ value. The results for GS2 are used for the detailed flow structure and the unsteady motion behind the detonation front due to the high resolution in the present simulations.

(a)

3.1. Overall structure of spinning detonation Spinning detonation has a complicated shockwave interaction with reaction zones. Most past experimental research produced useful information about this structure on the tube wall; however, the three-dimensional unsteady detonation front structure could not be represented in detail. Fig. 2 shows the instantaneous space isosurfaces and contours in the tube and on the tube wall for GS1. The calculated average detonation velocity during a cycle reaches approximately 1980 m/s, which is the Chapman–Jouguet value in this condition after the detonation begins to be overdriven. The present calculations show that the high concentration profiles of the OH mass fraction on the tube wall resembles a “curled ribbon” wrapped in a loose spiral, as shown in Fig. 2a. The ribbon originates in the detonation spin head. A long pressure trail was experimentally observed by Schott [10]. The calculated long pressure trail on the tube wall, which is also generated from the transverse detonation, appears clearly in Fig. 2. The detail of the pressure trail is discussed in a later section. 3.2. Mach leg and whisker Voitsekhovskii and his co-workers measured the Mach configuration by examining smoked disks attached to the end plate of the detonation tube. They suggested that the Mach configuration consists of a “leg” and one or two “whiskers,” based on their experimental results [13]. They also used the term “leg,” as in “Mach leg.” Duff [11] also measured the Mach configuration by the same method. The two measured Mach configurations do not coincide completely, because the spinning detonations produced are thought

(b) Fig. 2. Instantaneous space isosurfaces and contours for GS1. The gray space isosurface is the detonation front, and the cyan one is the pressure area of 3 MPa. The white broken arrow is the propagating direction of the detonation front, and the black and pink arrows near the detonation front are the rotating direction. I—incident shock, M—Mach stem, TD—transverse detonation, and LT—long pressure trail, respectively. (a) Three-dimensional view of OH mass fraction contours on the tube wall and pressure isosurfaces in the tube. (b) Pressure contours on the wall. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

to be different states, one being stable and the other one unstable. Fig. 3 presents the typical structure of the spinning detonation in the present calculations: Fig. 3a shows the detonation front shape viewed from the front side, Fig. 3b is the pressure isosurfaces and contours behind the detonation front in the tube, and Fig. 3c is the pressure contours on the tube wall. These figures are the results obtained for GS2. At t = 23.53 µs in Fig. 3a, a triple line forming a stationary Mach configuration with a strongly developed leg is observed near the left side wall. Two whiskers, which are obviously together with the Mach leg, weaken rapidly

N. Tsuboi et al. / Combustion and Flame 149 (2007) 144–161

149

Fig. 3. Spinning detonation in the circular tube at various times for GS2. The gray space isosurface is the detonation front. The orange space isosurface is the pressure of 4 MPa. The white arrow is the propagating direction of the detonation front, and the pink arrow is the rotating direction of the transverse detonation. I—incident shock, M—Mach stem, ML—Mach leg, W—“whisker,” TD—transverse detonation, RS—reflected shock wave, and LT—long pressure trail, respectively. (a) Detonation front shape viewed from the front. (b) Pressure isosurface and contours in the tube. (c) Pressure contours on the tube wall (r = R). (d) Pressure contours on r = 0.78R. (e) Pressure contours on r = 0.57R. (f) Pressure contours on r = 0.31R.

with the increase of the distance from the wall of the tube. The Mach leg stands orthogonal to the tube wall. At t = 23.74 µs, one whisker disappears, whereas another whisker and the Mach leg merge. The inboard part of the Mach leg travels faster than its outboard

part. The former overtakes the latter at t = 23.97 µs. Then the Mach leg separates into a Mach leg and two whiskers again at t = 23.53 µs. The Mach leg rotates in the counterclockwise direction and this process repeats in the present results.

150

N. Tsuboi et al. / Combustion and Flame 149 (2007) 144–161

Fig. 4. Contours at various times on θ = 0◦ and 180◦ planes in the tube for GS2. The white arrow is the propagating direction of the detonation front. M—Mach stem, TD—transverse detonation, RS—reflected shock wave, K—keystone shape, and E—extended transverse shock wave, respectively. The bottom wall in these figures corresponds to the dashed line in Fig. 6a.

The detonation front, including the Mach leg and whiskers, is imaged by the measurement of an endplate soot pattern for a spinning detonation. Duff [11] and Voitsekhovskii et al. [13] measured a pattern produced by a single spinning detonation and presented a detonation front pattern similar to that of Fig. 3a at t = 23.53 µs. Munday et al. [44] and Schott [10] also measured this pattern and their results are similar to our simulations at t = 23.74 or 23.97 µs presented in Fig. 3a. All of their results show different front patterns; however, the present results reproduce most of them. Some work has been carried out on the extent of the Mach leg along the radius of tube. Data from this work and the stability of the spinning detona-

Table 2 The extent of the Mach leg and the stability of the spinning detonation Reference

Extent

Stability

Ul’yanitskii [45] Voitsekhovskii et al. [13] Duff [11] This work (GS2)

0.5R 0.5R 0.2R 0.2–0.5R

Unknown Uniform Irregular Irregular

tion are summarized in Table 2. The extent of the Mach leg changes between 0.2 and 0.5R as shown in Fig. 3, where R is the radius of the tube. That of the Mach leg in θ = 0◦ and θ = 180◦ planes in Fig. 4a at t = 23.87 µs is also 0.5R. The present re-

N. Tsuboi et al. / Combustion and Flame 149 (2007) 144–161

151

sults correspond to a periodically irregular spinning detonation. Ul’yanitskii [45] carried out experimental measurements and proposed an extent of the Mach leg of 0.5R. Voitsekhovskii et al. [13] reported that the extent of the Mach leg was 0.6R and his experiments show an uniform spinning detonation. The extent of the Mach leg in Duff’s experiments [11], which have an irregular spinning detonation, is 0.2R. Thus, the extent of the Mach leg appears to be dependent on the stability of the spinning mode. 3.3. Transverse detonation The shock structure of the detonation evolves continuously from a single Mach reflection to a double Mach reflection and sometimes a complex Mach reflection during a cycle in a two-dimensional simulation [46]. The complex Mach reflection, which is a characteristic of a spinning detonation, appears in the detonation limit. It consists of a transverse detonation that is locally coupled with a shock wave and a combustion wave. The transverse detonation plays an important role in the spinning detonation and its three-dimensional time-dependent shape is shown in Fig. 3b. The motion of the front Mach configuration in Fig. 3a affects the shape of the reflected shock wave behind the detonation front. While the reflected shock wave is “kneed” at t = 23.53 µs, the reflected shock wave “rise” occurs at t = 23.74 and 23.97 µs. Because of the unsteady motion of the reflected shock wave and of the Mach configuration, the transverse detonation also moves unsteadily. The extent of the transverse detonation and the reflected shock wave is 0.5–0.6R, as shown in Figs. 3d and 3e, and they coincide with that of the Mach leg. Though periodic motion of the detonation front is observed, the pressure contours on the tube wall retain an invariable shape in Fig. 3c. The transverse detonation with a high pressure region is observed to propagate from bottom to top in the figure. The detailed structure of the head of the spinning detonation on the wall is shown in Fig. 5. AB and BD are the reflected shock wave; BC is the transverse detonation. The chemical reaction is creating near the Mach stem AM, the contact surface AD, the reflected shock wave BD, the reaction front CF parallel to the incident shock AI, and the transverse detonation BC. Its location is far from the incident shock AI and the induction length becomes large. The transverse detonation BC is terminated at point C and the extended transverse shock wave CE is prolonged downstream. The interaction between AD and BD would generate another contact surface; however, it is not clearly seen in the figure. The contact surface can be observed in Figs. 6b and 6c clearly. The physical values in the detonation structure of Fig. 5a are listed in Table 3. These

Fig. 5. Detailed structure of the head of the spinning detonation on the wall. Density is represented by white contours, and the local specific energy release is also shown as a contour surface. AB, BD—reflected shock wave, BC— transverse shock wave, AD—contact surface, CF—reaction front behind the incident shock, and CE—extended transverse shock wave, respectively. The leading shock is made off the incident shock AI and Mach stem AM. Table 3 Physical values in the detonation structure of Fig. 5a Shock i Pressure (MPa) Incident Mach AB BC BD

Temperature (K)

Density (kg/m3 )

0 0.1 1 1.75

300 1106

2 3 4 5

3470 4.02 567 1865 7.62 1175 3569 15.10 451 3569–3582 10.73–15.10 451–483

5.13 5.52 19.7 14.4–19.7

0.85 3.99

Velocity 1980 1674

values are instantaneous ones at t = 33.08 µs for GS1. There are no experimental data or analytical results in the same mixture conditions though similar data are presented by Huang et al. [16] and Voitsekhovskii et al. [13] for the different mixture conditions. Then the

152

N. Tsuboi et al. / Combustion and Flame 149 (2007) 144–161

Fig. 6. Instantaneous contours on the tube wall for GS2. The white arrow is the propagating direction of the detonation front, and the pink arrow is the rotating direction of the transverse detonation. I—incident shock, M—Mach stem, TD—transverse detonation, LT—long pressure trail, and C—contact surface, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

present results will be compared with the experimental data in the same mixture conditions. The instantaneous shock pattern was experimentally observed by Schott [10] and Voitsekhovskii et al. [13] and theoretically predicted by Macpherson [47] and Voitsekhovskii et al. [13]. The present shock pattern with the transverse detonation on the wall agrees with their results. The size of the highpressure region over 10 MPa oscillates slightly near the transverse detonation in Fig. 6a.

long pressure trail is observed behind the detonation in Figs. 3b and 3c as well as in the experiments. The trail is generated from the spin head and follows after the Mach leg. The frequency of rotation of the trail coincides with that of the Mach leg. With increasing distance from the transverse detonation, the trail gradually degenerates into compression waves of finite amplitude. The length of the trail is approximately twice the tube diameter, and that of the compression waves below the trail is six times longer than the tube diameter, as shown in Fig. 2b.

3.4. Pressure trail 3.5. Keystone and unburned gas pocket The pressure trail, which extends parallel to the axis of the tube, was also discovered as a long luminous band by Campbell in 1926 [1–3]. The cause of the luminosity, whether the source of the luminescence is the pressure trail or the transverse detonation, was unclear in the past literatures; however, the existence of the long pressure trail is one of the features in spinning detonation. In the numerical simulation, the

The formation of an unburned pocket behind the detonation front has been studied numerically by Oran et al. [41,48]. Gamezo et al. [26] reported that unburned pockets are commonly observed. A recent advanced experimental technique, the method of planar laser-induced predissociated fluorescence (LIPF), has enabled Pintgen et al. to visualize the

N. Tsuboi et al. / Combustion and Flame 149 (2007) 144–161

153

Fig. 7. Spinning detonation in the circular tube at various times for GS2. The gray and green space isosurfaces in pressure are the detonation front and the pressure of 6 MPa. The green space isosurfaces in temperature is the temperature of 1400 K. The red space isosurfaces in OH and H2 mass fractions are 0.028. The white arrow is the propagating direction of the detonation front, and the pink arrow is the rotating direction of the transverse detonation. TD—transverse detonation, and LT—long pressure trail, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

OH concentration behind the detonation front [49]. They obtained images of the creation of a keystoneshaped region behind the detonation front, and they also showed that no unburned pockets could be observed. Our three-dimensional simulations in a rectangular tube [30] showed that a three-dimensional unburned gas pocket behind the front, generated by a shutter of the triple lines, takes the shape of a ring and disappears later. The shape of a peninsula was also observed in the three-dimensional OH concentrations behind the detonation front. The structure of these detonations in these experimental and nu-

merical works is cellular; however, the profile of OH concentration behind the spinning detonation has not been reported yet. Therefore we now present the OH and H2 concentrations behind the detonation front obtained in the simulations. Fig. 7c presents OH mass fraction isosurfaces and contours in the tube, while Fig. 6c shows OH mass fraction contours on the tube wall, and Fig. 4c shows OH mass fraction contours on θ = 0◦ and 180◦ planes in the tube. The high OH mass fraction area delays the appearance in Fig. 7c of the high-pressure area near the transverse detonation in Fig. 7a. This de-

154

N. Tsuboi et al. / Combustion and Flame 149 (2007) 144–161

Fig. 8. Comparison of maximum pressure histories on the wall.

lay is an induction delay. The unburned mixture is thermally activated by the high-temperature area behind the shock front, and reactions occur yielding OH radicals as intermediate products. Therefore, the time lag between OH concentration and pressure field is caused by the hydrogen/oxygen reaction system. A similar trend has already been observed in a rectangular tube [30]. High OH mass fraction areas are also observed near the transverse detonation, the Mach stem, and the extended transverse shock wave on the wall in Fig. 6. The keystone-shaped region of high OH mass fraction appears between the Mach stem and the contact surface on the wall. Fig. 4c, which resembles an experimental view of a tube cross section visualized by LIPF, shows a similar keystone-shaped region at t = 23.87 µs; however, another high OH mass fraction area is observed behind it due to the transverse detonation and the extended transverse shock wave. Fig. 7d presents H2 mass fraction isosurfaces and contours in the tube, while Fig. 6d shows H2 mass fraction contours on the tube wall, and Fig. 4d H2 mass fraction contours on θ = 0◦ and 180◦ planes in the tube. Unlike the case of the detonation with cellular structures, no unburned pocket is observed in the present results, independent of the grid resolutions. The transverse detonation can completely consume unburned H2 behind the front. The unburned pockets would be generated by shutter effects of the collision of the triple lines in the cellular detonation, but there is no such collision for the spinning detonation. Therefore no unburned pocket would be generated in the spinning detonation.

3.6. Pitch and track angle The shock structure on the wall and the pitch of the soot records, which are measured in the experiments, are presented here in order to validate our simulations. Fig. 8a shows the maximum pressure history on the tube wall. The spinning detonation pattern that is observed experimentally resembles a ribbon wrapped in a loose spiral. The ratio of the pitch to the tube diameter in the experiments was approximately 3 and its theoretical value is 3.13 [8], as mentioned in Section 1. A comparison between the experimental and numerical results is given in Table 4. The ratio in the present simulations equals 3.09 for GS1 and 3.29 for GS2, which shows that the present calculation is reliable. The track angle α, the incident shock angle Φ1 , and the Mach stem angle Φ3 in the various experimental data are listed in Table 4. Φ1 and Φ3 for the experimental results are calculated by Huang et al. [16]. Our numerical results are measured in Fig. 5. Although the experimental conditions such as compositions of the mixture, tube diameters, and initial pressures are largely different, Φ1 is 30◦ –40◦ and Φ3 is close to 90◦ . This means that the overall shock structure has two highly overdriven detonation fronts: the Mach stem AM and the transverse detonation BC in Fig. 5. The present results are not a case with uniform spin, but with a periodically irregular spinning detonation. Schott [10] measured both uniform and irregular spin and mentioned that the periodically irregular spin was commonly observed under his experimental

N. Tsuboi et al. / Combustion and Flame 149 (2007) 144–161

155

Table 4 Comparison of the track angle, the pitch angle and the flow incident angle of the incident shock Φ1 and the Mach stem Φ3 Mixture

Tube diameter (mm)

Initial pressure (Torr)

Detonation velocity (m/s)

Track angle, α (deg.)

Φ1 (deg.)

Φ3 (deg.)

P /D 

C2 H2 + 1.5O2 + 12.5Ara 2H2 + O2 a 2H2 + O2 + 3Ara 2CO + O2 + 5% H2 a 2CO + O2 + 3% H2 b 1.5H2 + 1.5O2 + 7Arc C2 H2 + 7.58O2 + 34.3Arc Corn starch particles + O2 (φ = 1)d Corn starch particles + O2 (φ = 0.45)d 2H2 + O2 e 2H2 + O2 + 3Are 2H2 + O2 + 7Are H2 + Air(Stoich.)e H2 + Air(Stoich.)f H2 + Air(Stoich.)g

21 21 21 21 27 90 90 141

45 48 40 80 76 22 30 380

1637 2688 1816 1760 1700 1325 1227 1801

49 47 46 45 44.2 46.8 48.7 45

29.2 27.8 34.2 33.6 35.6 30.0 32.0 37

89.6 89.8 87.6 87.8 87.1 89.9 88.9 85

2.73 2.93 3.03 3.14 3.23 2.95 2.76 3.14

141

380

1801

42

40 40 40 40 1 1

27 23 46 53 760 760

2350 1500 1366 1690 1980 2020

47.24 45.1 46.35 46.6 45.5 43.7

3.48

36 38

83 89

3.09 3.29

Note. D  : Circle of periphery, P : pitch, φ: equivalence ratio. Original data for Φ1 and Φ3 except the present calculations are listed in Ref. [16]. a Nikolaev et al. [53]. b Voitsekhovskii et al. [13]. c Huang et al. [16]. d Zhang and Grönig [50]. e Ul’yanitskii [45]. f This work (GS1). g This work (GS2).

Fig. 9. Instantaneous pressure isosurfaces and contours in the tube for various grid resolutions. The gray isosurfaces denote the detonation front and the green isosurface is the pressure of 6.0 MPa. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

conditions. Topchian and Ul’yanitskii [17] analyzed the instability of the spinning detonation using the shock and detonation polar; however, a more detailed analysis is necessary to find the instability region in the spinning detonation with the three-dimensional effects. As mentioned previously, the present results have a periodic fluctuation on the Mach leg, the transverse detonation, and the reflected shock wave. These fluctuations would cause the irregular spin on the maximum pressure history on the tube wall. However,

the cause of the periodic fluctuation is unclear in the present results and its identification will require future work. 3.7. Detonation velocity and angular velocity The detonation velocity was calculated for GS1 in order to confirm the CJ velocity in the simulation. Fig. 10 shows the history of the instantaneous detonation velocities along the symmetry axis in the

156

N. Tsuboi et al. / Combustion and Flame 149 (2007) 144–161

Fig. 10. Instantaneous detonation velocity at θ = 0◦ for GS1.

tube (r = 0), with θ = 0◦ along the wall (r = R), and with θ = 0◦ in the tube (r = 0.57R). The detonation velocity varies from 0.75DCJ to 1.6DCJ along the wall, but the variation of the velocity is approximately 0.9–1.1DCJ along the symmetry axis. The variation of the velocities at r = 0.57R are 0.8–1.4DCJ , which are smaller than those along the wall. Our past research for two-headed heterogeneous detonation in a tube [34] shows that the detonation velocity varies from 0.7DCJ to 1.4DCJ along the tube axis, as well as along the wall. The trend of the present results differs from that of the two-headed detonation along the tube axis. The present results indicate that there exists a zone of silence around the tube axis without the rotating Mach leg and the transverse shock wave behind the detonation front. The zone of silence around the tube axis is also observed in Fig. 4. The calculated averaged detonation velocity in a cycle is plotted in Fig. 11. The present simulations have 27 cycles. At the beginning of the spinning detonation process, the velocities are overdriven at both locations; however, the velocities decrease linearly to become the CJ detonation velocity of 1980 m/s. The calculated average velocities fall below the CJ values after 17 cycles; then they become slower than the CJ value if the simulations are extended. Zhang and Grönig [50] tried to reveal some details of the interior structure of the single spinning heterogeneous detonation by measuring pressure in the tube, and they analyzed the results using the Navier–Stokes equations with the wave-front-fixed coordinate. They assumed that the spin in the tube has a constant angular velocity Ω, where Ω is a positive value. The solutions of these equations are ur = 0,

uθ = −rΩ,

and ux = −D,

(8)

Fig. 11. Calculated average detonation velocity during one cycle for GS1.

where ur , uθ , and ux are the radial, circumferential, and axial velocity components. D and r are the axial detonation velocity and radius, respectively. The spin track angle is α = tan−1 rΩ/D.

(9)

In the present simulation r = 0.5 mm, α = 45◦ , and D = 1980 m/s for GS1; therefore Ω = −3.96 × 106 rad/s in the tube and uθ = −1980 m/s on the tube wall. Fig. 12 shows the instantaneous pressure, circumferential velocity, and angular velocity behind the detonation front in the tube. The circumferential velocity varies between negative and positive values on the five cross sections. The angular velocities near the spin head in region A in Fig. 12c are close to the value of Ω = −3.96 × 106 rad/s; however, no uniform and positive angular velocities exists in region A. The instantaneous pressure and the circumferential and angular velocities on the tube wall are shown in Fig. 13. The spirals of the circumferential and angular velocities, which are generated by the contact surface near the spin head as shown in Fig. 6, are clearly visible. The circumferential and angular velocities in the region between the Mach stem and the contact surface near the detonation front (region A in Fig. 13) are approximately −1980 m/s and −3.96 × 106 rad/s. Those in the region behind the long pressure trail (region B in Fig. 13) are half the value of those in region A in Fig. 13. Therefore both the circumferential and angular velocities on the wall also vary in the whole domain on the wall and cannot be assumed to be constant. The calculated average circumferential velocity and average angular velocity in the tube are presented

N. Tsuboi et al. / Combustion and Flame 149 (2007) 144–161

157

3.8. Detonability limit related to viscous and heat transfer effects

(a)

(b)

The present tube diameter in the simulation is 1 mm and some questions arise, such as (i) whether detonation propagates in such a small diameter and (ii) how great are the viscous and heat transfer effects near the wall. Both questions relate closely because the viscous and heat transfer effects cannot be neglected as the tube diameter becomes small. The detonation cell size λ for stoichiometric hydrogen/air at 0.1 MPa is approximately 10 mm, for example, Shepherd [51], and the empirical limit criterion λ/π is 3.2 mm. An accurate measured value for stoichiomeric hydrogen/air does not exist; however, some measured data in lean and rich fuel cases and analytical data exist [52]. They analyzed the detonation limit by one-dimensional simulation with heat and viscous losses and concluded that the minimum diameter would be about 3 mm for stoichiometric hydrogen/air at 0.1 MPa. However, the minimum diameter probably becomes smaller than 3 mm when the heat and viscous losses are neglected. In the present simulation, the largest scales generated by the shock wave might not be affected by this neglect of viscosity. We also use the results of onedimensional detonation for the start of 3D simulations and the deflagration-to-detonation (DDT) process is neglected; therefore the DDT process, which is strongly affected by the viscous and heat transfer effects, is not included in the start of 3D simulations. In the future, the viscous and heat transfer effects on the detonation limit, as well as the DDT process, will be estimated numerically. 3.9. Numerical convergence

(c) Fig. 12. (a) Instantaneous pressure, (b) circumferential velocity, and (c) angular velocity in the tube for GS1. Wall contours in all figures are pressure. TD—transverse detonation.

in Fig. 14. Both calculated average velocities are approximately 15% of the CJ value near the spin head in region A. The velocities in the other regions below the spin head are about 10% of the CJ value and the rotation below the spin head does not rotate near the tube center. These results imply that the high-pressure region behind the detonation front contributes the negative constant angular velocity and the calculated average circumferential velocities are approximately 15% of the CJ value. Therefore the assumption of the constant angular velocity proposed by Zhang and Grönig [50] is invalid.

Numerical simulation with a detailed reaction model requires high grid resolution in space. In our study, grid resolution in the circular direction as well as in the axial direction is found to be very important in simulating a spinning detonation because the spinning detonation has a transverse detonation wave rotating along the circular wall. The comparison of the instantaneous results in the tube is shown in Fig. 9. The simulation of a spinning detonation fails and a complicated detonation front is observed for GS4. However, the simulation yields a front with one transverse detonation for GS2 and GS3. The shock wave system and the unsteady motion including a Mach leg and two whiskers for GS3 are similar to that for GS2. Fig. 8 presents a comparison of maximum pressure histories on the tube. While the results for GS4 have an incomplete spinning detonation, the results for GS2 and GS3 have a spinning detonation.

158

N. Tsuboi et al. / Combustion and Flame 149 (2007) 144–161

Fig. 13. (a) Instantaneous pressure, (b) circumferential velocity, and (c) angular velocity on the wall for GS1. TD—transverse detonation.

The initial process of the spinning detonation is unclear due to the artificial conditions set at the start of calculations; however, the strength and interaction between the initial explosions due to the initial arti-

ficial perturbation would provide a rotating direction. An adequate simulation of the initiation process including the deflagration-to-detonation process for the spinning detonation should be the subject of future

N. Tsuboi et al. / Combustion and Flame 149 (2007) 144–161

(a)

(b) Fig. 14. (a) Calculated average circumferential velocity and (b) angular velocity in the tube for GS1.

research. The pitch for GS3 approximately coincides with that for GS2 and the present simulations show reasonable results. The length of the computational domain and the boundary condition on the downstream boundary should be carefully chosen, as is the grid resolution. The detonation velocities for GS2–GS4 are approximately 2020 m/s because of the short computational domain and pressure-fixed conditions on the downstream boundary. The long pressure trail extends to interact with the downstream boundary and the pressure near the boundary would be larger than the CJ pressure. The short computational domain also affects the pressure near the downstream boundary; therefore it is better if the length of the computational domain is six times the diameter and the nonreflected boundary condition is particularly suited to simulating the spinning detonation. Several unresolved phenomena remain in the spinning detonation; one of them is the structure of fine crosshatchings within the ribbon in the spinning detonation observed in the past experiments [10]. The crosshatchings would be generated by a transverse detonation with a “multiheaded” wave behind the

159

front. The present results cannot resolve such crosshatchings in the ribbon because of a lack of grid resolution. These crosshatchings depend on the stability of the spinning detonation [17]; therefore it should be studied, including stability as well as with simulations with higher grid resolution. Another phenomenon is that of detonation limits for gas composition as well as low pressure. This knowledge of the limits is, of course, of considerable practical importance for safety, many industrial processes, and combustion in pulse detonation engines. Gordon et al. [19] described the composition limit obtained by their experimental data to show that the limit for helium and argon-diluted mixture gave not only the normal spin, but also a short-wavelength spin or a longwavelength spin. The short-wavelength spin is said to be a multiple-headed detonation; however, the detailed structure of both spins is unknown. Their spins are expected to be simulated and become clear in the shock structure in further research. Thus far we have simulated a spinning detonation in a circular tube after a single spinning detonation was experimentally discovered in 1926 by Campbell and colleagues.

4. Conclusions Unsteady three-dimensional simulations were performed for hydrogen/air detonations in a circular tube and a single spinning detonation was simulated. A spinning detonation is composed of a transverse detonation rotating circumferentially around the wall, a Mach leg and whiskers consisting of triple lines, and a long pressure trail extending downstream of the tube. The present simulations revealed their threedimensional shape and their unsteady motions in the tube. The present results are an irregular mode of the spinning detonations, and the mode would depend on the compositions of the mixture, initial pressures, and tube diameters. The formation of an unburned gas pocket behind the detonation front was not observed in the present simulations because the rotating transverse detonation completely consumed the unburned gas. The calculated profiles of instantaneous OH mass fraction have a keystone shape behind the detonation front. The pitch angle, the track angle, the Mach stem angle, and the incident shock angle on the tube wall in the numerical results agree well with those in the experimental results, and they are independent of the compositions of the mixture, tube diameters, and initial pressures. Instantaneous detonation velocities in the simulations showed that the velocity varies from 0.75 DCJ to 1.6 DCJ along the wall, whereas its variation is ap-

160

N. Tsuboi et al. / Combustion and Flame 149 (2007) 144–161

proximately 0.9–1.1DCJ along the tube axis. Therefore there exists a zone of silence around the tube axis without the rotating Mach leg and the transverse detonation. The results for the calculated average angular velocity show that the assumption of a constant angular velocity is invalid except in the vicinity of the rotating head. A resolution study shows that the circumferential resolution as well as the axial resolution is very important because the rotating transverse detonation has to be resolved. Low resolution in the circumferential direction produced incomplete spinning detonation. Finally, the present simulation provided a new understanding of the physics involved in the spinning detonation, a phenomenon that was discovered in 1926.

Acknowledgments This research was carried out in collaboration with the Center for Planning and Information Systems in ISAS/JAXA. We thank Professor F. Zhang at Defence Research and Development Canada for useful and important comments.

References [1] C. Campbell, D.W. Woodhead, J. Chem. Soc. (1926) 3010–3021. [2] C. Campbell, D.W. Woodhead, J. Chem. Soc. (1972) 1572–1578. [3] C. Campbell, A.C. Finch, J. Chem. Soc. (1928) 2094– 2106. [4] W.A. Bone, R.P. Fraser, Philos. Trans. R. Soc. Sect. A 228 (1929) 197–234. [5] W.A. Bone, R.P. Fraser, Philos. Trans. R. Soc. Sect. A 230 (1932) 363–385. [6] W.A. Bone, R.P. Fraser, W.H. Wheeler, Philos. Trans. R. Soc. Sect. A 235 (1936) 29–68. [7] N. Manson, Compt. Rendus 222 (1946) 46–48. [8] J.A. Fay, J. Chem. Phys. 20 (1952) 942–950. [9] A.J. Predvoditelev, 7th Symp. (Int.) Combust. (1959) 760–765. [10] G.L. Schott, Phys. Fluids 8 (1965) 850–865. [11] R.E. Duff, Phys. Fluids 4 (1961) 1427–1433. [12] B.V. Voitsekhovskii, V.V. Mitrofanov, M.E. Topchian, 12th Symp. (Int.) Combust. (1969) 829–837. [13] B.V. Voitsekhovskii, V.V. Mitrofanov, M.E. Topchian, Izd-vo Sibirsk, Otdel. Akad. Nauk SSSR, Novosibirsk (1963), Translation 1966, Wright-Patterson Air Force Base Report, FTD-MT-64-527 (AD-633-821). [14] B.V. Voitsekhovskii, Dokl. Akad. Nauk SSSR 114 (1957) 717–720 [translation MDF-V-123]. [15] Yu.N. Denisov, Ya.K. Troshin, Proc. Acad. Sci. USSR Phys. Chem. Sect. 125 (1959) 217–221.

[16] Z.W. Huang, M.H. Lefebvre, P.J. Van Tiggelen, Shock Waves 10 (2000) 119–125. [17] M.E. Topchian, V.Y. Ul’yanitskii, Acta Astronaut. 3 (1976) 771–779. [18] F.A. Bykovskii, S.A. Zhdan, E.F. Vedernikov, Combust. Explos. Shock Waves 41 (2005) 449–459. [19] W.E. Gordon, A.J. Mooradian, S.A. Harper, 7th Symp. (Int.) Combust. (1959) 752–759. [20] H.O. Barthel, Phys. Fluids 17 (1974) 1547–1553. [21] A.R. Kasimov, D.S. Stewart, J. Fluid Mech. 466 (2002) 179–203. [22] O.V. Achasov, O.G. Penyazkov, Shock Waves 11 (2002) 297–308. [23] J.H. Lee, R.I. Soloukhin, A.K. Oppenheim, Astronaut. Acta 14 (1969) 565–584. [24] S. Taki, T. Fujiwara, AIAA J. 16 (1978) 73–77. [25] K. Kailasanath, AIAA J. 38 (2000) 1698–1708. [26] V.N. Gamezo, D. Desbordes, E.S. Oran, Shock Waves 9 (1999) 11–17. [27] D.N. Williams, L. Bauwens, E.S. Oran, Proc. Combust. Inst. 26 (1996) 2991–2998. [28] D.N. Williams, L. Bauwens, E.S. Oran, 16th International Colloquium on the Dynamics of Explosions and Reactive Systems (1997) 374–376. [29] N. Tsuboi, S. Katoh, A.K. Hayashi, Proc. Combust. Inst. 29 (2002) 2783–2788. [30] K. Eto, N. Tsuboi, A.K. Hayashi, Proc. Combust. Inst. 30 (2005) 1907–1913. [31] K. Eto, N. Tsuboi, A.K. Hayashi, AIAA Paper 20040311 (2004). [32] A.K. Hayashi, K. Eto, N. Tsuboi, 20th International Colloquium on the Dynamics of Explosions and Reactive Systems (2005) No. 85. [33] T. Washizu, T. Fujiwara, 14th ICDERS (1993) D1-9. [34] N. Tsuboi, A.K. Hayashi, Y. Matsumoto, Shock Waves 10 (2000) 277–286. [35] N. Tsuboi, K. Eto, A.K. Hayashi, 20th International Colloquium on the Dynamics of Explosions and Reactive Systems (2005) No. 71. [36] D. Stull, H. Prophet, JANAF Thermochemical Tables, second ed., NSRDS-NBS37 (1971). [37] H.C. Yee, NASA Technical Memorandum 89464 (1987). [38] Y. Liu, M. Vinokur, AIAA Paper 89-0201 (1989). [39] E.L. Petersen, R.K. Hanson, J. Propul. Power 15 (1999) 591–600. [40] K. Eto, N. Tsuboi, A.K. Hayashi, 19th International Colloquium on the Dynamics of Explosions and Reactive Systems (2003) No. 178. [41] E.S. Oran, J.W. Weber Jr., E.I. Stefaniw, M.H. Lefebvre, J.D. Anderson Jr., Combust. Flame 113 (1998) 147–163. [42] M.A. Sussman, AIAA Paper 1994-3101 (1994). [43] K. Inaba, A. Matsuo, AIAA Paper 2001-480 (2001). [44] G. Munday, A.R. Ubbelohde, I.F. Wood, Proc. R. Soc. London Ser. A 306 (1968) 171–178. [45] V.Y. Ul’yanitskii, Fiz. Goreniya Vzryva 16 (1980) 105– 111. [46] M.H. Lefebvre, E.S. Oran, Shock Waves 4 (1995) 277– 283. [47] A.K. Macpherson, Astronaut. Acta 14 (1969) 549–558.

N. Tsuboi et al. / Combustion and Flame 149 (2007) 144–161 [48] E.S. Oran, T.R. Young, J.P. Boris, J.M. Picone, D.H. Edwards, Proc. Combust. Inst. 19 (1982) 573–582. [49] F. Pintgen, C.A. Eckett, J.M. Austin, J.E. Shepherd, Combust. Flame 133 (2003) 211–229. [50] F. Zhang, H. Grönig, Phys. Fluids A 3 (8) (1991) 1983– 1990.

161

[51] J. Shepherd, http://www.galcit.caltech.edu/detn_db/ html/H2-Air1.html. [52] G.L. Agafonov, S.M. Frolov, Combust. Explos. Shock Waves (USSR) 30 (1) (1994) 91–100. [53] Y.A. Nikolaev, M.E. Topchian, V.Y. Ul’yanitskii, Fiz. Goreniya Vzryva 14 (1978) 106–109.