Detonation engine fed by acetylene–oxygen mixture

Acta Astronautica 104 (2014) 134–146

Contents lists available at ScienceDirect

Acta Astronautica journal homepage: www.elsevier.com/locate/actaastro

Detonation engine fed by acetylene–oxygen mixture N.N. Smirnov a,b,n, V.B. Betelin a,b, V.F. Nikitin a,b, Yu.G. Phylippov a, Jaye Koo c a

Moscow M.V. Lomonosov State University, Moscow 119992, Russia Scientific Research Institute for System Studies of Russian Academy of Sciences, Moscow 117218, Russia c School of Aerospace and Mechanical Engineering Korea Aerospace University, Goyang-city 411-791, Republic of Korea b

a r t i c l e in f o

abstract

Article history: Received 10 July 2014 Accepted 11 July 2014 Available online 22 July 2014

The advantages of a constant volume combustion cycle as compared to constant pressure combustion in terms of thermodynamic efficiency has focused the search for advanced propulsion on detonation engines. Detonation of acetylene mixed with oxygen in various proportions is studied using mathematical modeling. Simplified kinetics of acetylene burning includes 11 reactions with 9 components. Deflagration to detonation transition (DDT) is obtained in a cylindrical tube with a section of obstacles modeling a Shchelkin spiral; the DDT takes place in this section for a wide range of initial mixture compositions. A modified ka-omega turbulence model is used to simulate flame acceleration in the Shchelkin spiral section of the system. The results of numerical simulations were compared with experiments, which had been performed in the same size detonation chamber and turbulent spiral ring section, and with theoretical data on the Chapman–Jouguet detonation parameters. & 2014 IAA. Published by Elsevier Ltd. All rights reserved.

Keywords: Detonation Engine Propulsion Thermodynamic efficiency Deflagration Transition

1. Introduction One of the schemes for producing enhanced thrust at both static and dynamic conditions is pulse detonation. The thermodynamic efficiency of the Chapman–Jouguet detonation as compared to other combustion modes is due to the minimal entropy of the exhaust jet. Efforts have been made during past several decades to show that proper utilization of the operation cycle does result in improved performance. However, there are several issues in developing this technology, which represent scientific and technological challenges. The success in resolving these problems will determine the implementation of pulse detonation propulsion. The control of detonation onset is of major importance in pulse detonating devices. The advantages of detonation

n Corresponding author at: Moscow M.V. Lomonosov State University, Moscow 119992, Russia. Tel.: þ7 495 9391190; fax: þ7 495 9393754. E-mail address: [email protected] (N.N. Smirnov).

http://dx.doi.org/10.1016/j.actaastro.2014.07.019 0094-5765/& 2014 IAA. Published by Elsevier Ltd. All rights reserved.

over constant pressure combustion stimulated developing schemes incorporating deflagration to detonation transition (DDT) as a part of working cycle and shortening the pre-detonation length. The probable application of these principles to create the new generation of engines put the problem of DDT on top of current research needs. Thus, the problem of DDT control in gaseous fuel–air mixtures became very acute. This paper contains the results of theoretical and experimental investigations of deflagration to detonation transition (DDT) processes in combustible gaseous mixtures of acetylene and oxygen in different compositions. The influence of geometrical characteristics of the confinement and fuel concentration in the unburned mixture on the onset of detonation is discussed. Detonation provides essential advantages as compared with classical combustion engine schemes from the point of view of energy conversion rate, which by three orders of magnitude surpasses the rate of energy release in classical combustion chambers. Thus the detonation regime could

N.N. Smirnov et al. / Acta Astronautica 104 (2014) 134–146

135

turn out to be much more effective in promoting fuel combustion in supersonic flows. Besides, one of the limitations in developing engines for supersonic flights is that of stabilization of combustion in supersonic flow, because flame propagation velocity is essentially subsonic, while detonation could be stabilized for the high Mach numbers. The specific impulse of detonation engine is higher as compared with specific impulse of a ram jet even operating under similar conditions and having the same fuel and oxidizer average mass flow rate. The rough estimate provides the following ratio: rffiffiffiffiffi I PDE pffiffiffi cp Z γ¼ I RAM cv Besides, the power of an engine could be increased by increasing essentially the mass flow rate of fuel and oxidizer due to the increased opportunity of the energy conversion rate. Thermodynamic efficiency of detonation mode is higher than that for deflagration, which follows from the general theory of detonation. On propagating combustion wave in a combustible mixture the final state of reaction products, as it follows from basic conservation laws, should be found on the Hugoniot curve, which has the following simplified representation: Γ ð1Þ ðp; ϑ; p0 ; ϑ0 Þ ¼ Q 0 ; Γ ð1Þ ðp; ϑ; p0 ; ϑ0 Þ ¼ eð1Þ ðp; ϑÞ eð1Þ ðp0 ; ϑ0 Þ þ 12 ðϑ ϑ0 Þðp þ p0 Þ where p is the pressure, ϑ is the specific volume, e is the specific energy, Q 0 is the specific energy released in combustion wave due to chemical reaction, superscript 1 stands for reaction products, and subscript 0 stands for unreacted mixture. At the same time, the final state of reaction products could be found on the Michelson–Rayleigh lines p1 p0 ¼  m2 ; ϑ1 ϑ0 where m2 ¼ ρ20 w20 is the second power of mass flow rate through a surface unit of combustion wave, and w0 is the combustion wave propagation velocity relatively to unburned mixture. Thus optional intersections of the Hugoniot curve and the Michelson–Rayleigh lines give us in the plane ðp; ϑÞ the geometrical place of all states of the medium, which could be obtained behind the combustion wave. Fig. 1 illustrates location of the Hugoniot curves and the Michelson–Rayleigh lines. In general, there could be a maximum of two intersections in the sector p 4p0 ; ϑ oϑ0 , which correspond to compression waves. These waves have the name of strong and weak detonation, and the limiting case of those intersections coincidence gives us a tangent corresponding to the Chapman–Jouguet detonation. The two intersections in the sector p 4 p0 ; ϑ o ϑ0 correspond to deflagration mode, which is also named flame propagation. Summarizing the advantages of detonation mode for propulsion applications as follows:

 high thermodynamic efficiency of the Chapman–Jouguet detonation as compared to other combustion modes is due to the minimal entropy of the exhaust jet,

Fig. 1. (a) Hugoniot curves and Michelson–Rayleigh lines; (b) Huginiot function variation along Michelson–Rayleigh lines of different inclinations; (c) entropy variation along Hugoniot curves and Michelson– Rayleigh lines.

 CO emission reduction,  high rate of energy conversion (103 times), and  specific impulse increase. The first studies on pulse detonation engines as reported in [1] were performed by Hoffmann in 1940, wherein mixtures of acetylene, benzyl and liquid oxygen were used. More than 10 years later Nicholls et al. investigated pulse detonation devices fed by mixtures of hydrogen and/or acetylene with oxygen [2]. The first investigations of pulse detonating devices using mixtures of hydrocarbon fuels with atmospheric air were performed in the former Soviet Union in the beginning of the 1980's [3], and demonstrated the principle possibility of creating pulse detonation engines fed by available liquid hydrocarbon fuels and atmospheric air as an oxidizer. Investigations of pulse detonation engines were continued in the 1990's of the XX-th century and beginning of the XXI-st century [4–8]. Most of the effective schemes were based on deflagration to detonation transition principles [1–4,8]. Schemes based on detonation transmission from a narrow

136

N.N. Smirnov et al. / Acta Astronautica 104 (2014) 134–146

tube into a chamber of large diameter were also investigated [9], as well as schemes using rotating detonation waves in a supersonic flow [10–12], and schemes of high frequency resonators [13]. Advantages and limitations of pulse detonation engines were discussed in [14]. The most important issue for developing pulse detonation engines is the guarantee of a periodical onset of stable detonation waves [15–17]. Due to this reason fundamental investigations of deflagration to detonation transition constitute an important part of pulse detonation engine investigations. A new scheme of transmission of detonation from an outer narrow gap into a wider tube by means of converging shock formation for hydrogen–oxygen mixtures was studied in [29]. In 1881 Mallard and Le Chatelier [18], and Bertelot and Vieille [19] were the first to confront deflagration to detonation transition process. On investigating flame propagation in homogeneous gaseous mixtures they unexpectedly detected the onset of a supersonic mode of combustion wave traveling at velocities of thousands of meters per second. This combustion mode got the name of “false” or “out of tone” combustion, which sounds in French like “detonation” being originated from the French verb “détonner”. The existence of two alternative velocities for the combustion process needed theoretical explanation, which was given in 1893 by an Associate Professor of Moscow University V.A. Mikhelson in his paper “On normal combustion velocity of explosive gaseous mixtures” [20]. Based on papers, having been published by that time by Rankine [21] and Hugoniot [22], Mikhelson was the first to explain that the mechanism of flame propagation in detonation was not the heat conductivity, but “adiabatic heating up to the ignition point in shock waves”. This was the birth of the classical theory of detonation, which found its further development in the papers by Chapman [23] and Jouguet [24]. Investigations of deflagration to detonation transition (DDT) in hydrogen–oxygen mixtures [25–28] and later in hydrocarbons–air mixtures [3,31] showed the multiplicity of the transition process scenarios. The various modes of the detonation onset were shown to depend on the particular flow pattern created by the accelerating flame, thus making the transition process non-reproducible in its detailed sequence of events. By now, there exist different points of view on the DDT mechanism: the “explosion in the explosion” mechanism by Oppenheim [26,28] and the gradient mechanism of “spontaneous flame” by Zeldovich [32]. The later theoretical analysis showed that micro-scale non-uniformities (temperature and concentration gradients) arising in local exothermic centers (“hot spots”) ahead of the flame zone could be sufficient for the onset of detonation or normal deflagration [33–39]. Analysis and comparison of theoretical and experimental results showed that self-ignition in one or in a number of hot spots ahead of the accelerating flame followed by the onset of either detonation or deflagration waves brings to a multiplicity of the transition scenarios [40]. The common feature of many scenarios is the formation of local exothermic centers according to the stochastic Oppenheim mechanism followed by the onset of detonation at a

micro-scale in one of the exothermic centers according to Zeldovich mechanism [40–42]. Experimental and theoretical investigations [41,42] of the reflected shock–laminar flame interactions bringing to the onset of detonation showed, as well, that the transition to detonation in a hot spot takes place through the gradient mechanism, while the shock and flame interactions were important for creating the proper conditions for the hot spots to occur. In [48] a mechanism of detonation onset was described different from gradient mechanism. The deflagration propagating accelerates due to the flame front interaction with the flow ahead of the front, which on the final stage brings to a pressure peak and consequently shock wave formation inside the reaction zone. There exist two potential methods of adjusting detonation for propulsion applications: pulse detonation engines, and rotating detonation engines. Numerical simulation of both are present in [29,30]. In this paper we will concentrate on pulse detonation devices base on DDT principle. To promote DDT in tubes effective measures were suggested introducing the Shchelkin spiral in the ignition section [43], incorporating wider turbulizing chambers in the ignition section (Smirnov's cavities) [29,40], blocking the initial part of the tube with orifice plates. To enhance thrust nozzles are used, which bring peculiarities for outflow conditions [44,45]. Wider turbulizing chambers were discovered to provide for DDT both a promoting effect and an inhibiting effect depending on their number and location [46]. Insights on peculiarities of detonation onset in DDT are provided in [47–50]. 2. Problem statement A scheme of the detonation chamber used for experimental studies and numerical simulations is shown in Fig. 2. A more detailed description of the experimental procedure can be found in [51,52]. In our simulations, we did not calculate the process of filling the gas chamber with acetylene and oxygen, neither we did consider nitrogen diluting the mixture. Instead, we regarded a fully premixed non-moving gas in the chamber, which consisted of the ignition section, the spiral section, and the tube section. The whole system was closed, and had smooth internal walls, except for the spiral section. We used the following sizes in our calculations: D ¼24.4 mm – diameter of the detonation chamber; L ¼1240 mm – total length of the chamber; Lign ¼100 mm – length of the ignition section; Lspir ¼300 mm – length of the spiral section. The turbulent spiral ring had the following features: thickness 3.5 mm and pitch 7.1 mm, thus giving the blockage ratio of 49.8%, as in the experimental setup. The ignition was simulated by introducing energy into a small section of the system during a relatively short period; the location of the ignitor was at 70 mm from the left closed edge of the system (where the nitrogen valve is depicted in Fig. 2). Radius of the section where the energy was introduced was 2 mm, and the energy flux was

N.N. Smirnov et al. / Acta Astronautica 104 (2014) 134–146

137

Fig. 2. Detonation chamber used in experiments and numerical simulations.

2  1011 W/m3 during 10  5 s, resulting in about 64 mJ of ignition energy. The energy release was enough to ignite the mixture, but the corresponding shock was not intensive to lead to instant detonation formation. Therefore, the detonation was obtained each time via DDT process in the Shchelkin spiral section of the chamber.

The turbulence model includes 2 equations of the Wilcox ka-omega model [53], accomplished with the third equation modeling turbulent viscosity lag [54]. The model applies both for fully developed turbulence, and for the turbulence onset. The viscosity lag equation enables better modeling of the transient processes.

3. Mathematical model

  ∂ρK þ ∇ðρKuÞ ¼ ∇ ðμ þσ n μT Þ∇K þ P T  βn ρKω; ∂t

3.1. Equations and algebraic relationships The gas dynamic model consists of the equations of mass balance of components, momentum balance, and energy balance. It also incorporates 3 partial differential equations for turbulent parameters, and multiple algebraic equations of state, transport, turbulence model, and the kinetic model. The last describes mass transition between components in chemical processes, and the corresponding energy release and consumption. The gas dynamic model is as follows: ∂ρk _ k; þ ∇ðρk uÞ þ ∇Ik ¼ ω ∂t

k ¼ 1; …; N;

ð1Þ

∂ρu þ ∇ðρuu þ pÞ ¼ ∇τ; ∂t

ð2Þ

∂ρE _ þ∇ðρEu þ puÞ þ∇J ¼ ∇ðτuÞ þ q: ∂t

ð3Þ

Here ρk is the partial density of k-th species (total N components in gas mixture), u is the mean velocity vector, _ k is the rate of the species Ik is the diffusion vector, ω formation in chemical processes (mass/volume unit), ρ is the density of the mixture, p is the pressure, τ is the stresses tensor deviator, J is the conductive energy flux, and E is the total energy including internal, kinetic and chemical parts. The last term q_ is the intensity of external energy modeling the ignitor. The gas dynamic models (1)–(3) includes the features from the turbulence model, in particular, the spherical part of the turbulence stresses tensor, which is equal to ð 2ρK=3ÞU by definition, is being added to pressure (with an opposite sign). Also, the total energy includes an additional term turbulent kinetic energy, which is equal to K.

  ∂ρω ω þ ∇ðρωuÞ ¼ ∇ ðμ þσμT Þ∇ω þ α P T  βρω2 ; ∂t K   ∂μT K þ∇ðμT uÞ ¼ A ρ  μT : ω ∂t

ð4Þ ð5Þ ð6Þ

Here K is the kinetic energy of turbulence, μ is the dynamic molecular viscosity of the mixture, μT is the turbulent viscosity, P T is the turbulence production term, ω is the turbulence dissipation rate; σ, σ n , α, β, βn and A are the parameters of the model that either depend algebraically on flow parameters and turbulent Reynolds ReT and Mach M T numbers [53,55], or remain constant. Those turbulence model relations look as follows:   9 1 þ 680χ 2k 5=18 þðReT =8Þ4 3 1 þ FðM T Þ ; βn ¼ 2 4 100 1 þ 400χ k 1 þðReT =8Þ 2 ! 9 9 5=18 þ ðReT =8Þ4 3 1 þ70χ ω  FðMT Þ ; β¼ 125 100 1 þ ðReT =8Þ4 2 1 þ80χ ω 13 0:1 þ ReT =2:7 1 þ ReT =6 ; 25 1 þ ReT =2:7 β=3 þ ReT =6 35 0:01 þ ReT ; AT ¼ 100 1 þ ReT

α¼

1 ρK 2K σ ¼ σ n ¼ ; ReT ¼ ; M2T ¼ 2 ; FðMT Þ ¼ max f0; ðM2T  0:0625Þg; 2 ωμ a

  Ω Ω S 1 ∂K ∂ω ij jk ki χ k ¼ max 0; 3 ; χω ¼ : ð0:09ωÞ3 ω ∂xj ∂xj

ð7Þ

Here Sij are the components of the deformation velocity tensor, and Ωij are the components of the rotation tensor; both tensors are symmetric and antisymmetric parts of the velocity gradient 1 1 S ¼ ð∇u þ∇uT Þ; Ω ¼ ð∇u  ∇uT Þ: 2 2

ð8Þ

138

N.N. Smirnov et al. / Acta Astronautica 104 (2014) 134–146

turbulence model). They are as follows:

Turbulence production is expressed in terms of S P T ¼ μT S: S:

ð9Þ

In accordance with a Boussinesq turbulent model for fluxes, diffusion and energy conduction are modeled as follows:     μ μ ρ þ T ∇ k ; ð10Þ Ik ¼  Sc ScT ρ 

 μ μT þ ∇h ðμ þμT Þ∇K: Pr PrT

J¼ 

ð11Þ

Here Sc and ScT are the constant Schmidt numbers, molecular and turbulent, correspondingly; Pr and PrT are the constant Prandtl numbers, and h is the mixture enthalpy per mass unit. We use the constant Schmidt number model for molecular diffusion because the turbulent term is usually prevailing over the molecular one. We use enthalpy instead of temperature because the total energy includes the chemical part, as the enthalpy does. Also, the energy flux (11) contains a turbulence energy term. The full stress tensor contains spherical and deviator parts, and is modeled as  

 2 2  pU þτ ¼  p þ ρK Uþ ðμþ μT Þ 2S ð∇uÞU : ð12Þ 3 3 Here U is the unity tensor and p is thermodynamic pressure. Thermodynamic pressure is obtained from the state equation of a perfect gases mixture N

N Yk ¼ RG T ∑ X k ; W k k¼1 k¼1

p ¼ RG Tρ ∑

ð13Þ

where RG is the universal gas constant, T is the absolute temperature, Y k ¼ ρk =ρ is the mass share of a component, W k is the molar mass of the component, and X k is the partial molar density usually referred in literature as “molar concentration”. The system of concentrations X k is preferred in thermodynamics and chemical kinetics due to operations with molar mass values being omitted. Total energy E is expressed as follows: E ¼ eðY j ; TÞ þ

N uu þ K; e ¼ ∑ Y k ek ðTÞ: 2 k¼1

Density of the mixture is a sum of partial densities: N

K ¼ K 0 ; μT ¼ ReT0 μðX j ; TÞ; ω ¼ μT 1 ρK:

ð16Þ

k¼1

3.2. Initial and boundary conditions The initial conditions for gas dynamic equations stand for non-moving uniform gas mixtures with low level of turbulence (zero turbulence is invalid in ka-omega

ð17Þ

The turbulence parameters initial state is determined by initial kinetic turbulent energy K 0 , and initial turbulent Reynolds number, which is in fact relation between steady-state eddy viscosity and molecular viscosity. Given initial pressure and temperature, the species initial partial densities are determined with a single parameter x; they are revealed as follows: ρ0C 2 H2 ¼ W C 2 H2

p0 p0 x ; ρ0 ¼ W O2 ; RG T 0 ð1 þ xÞ RG T 0 ð1 þxÞ O2

ð18Þ

and zero initial partial density (and concentration) for all other species. The boundary conditions stand for conditions on a smooth impermeable adiabatic wall with no adsorption of components and no catalytic reactions. This results in zero normal velocity, and zero normal gradient of temperature and mass flux at the wall. un ¼ 0;

  ∂T ∂ ρk ¼ 0; ¼ 0: ∂n ∂n ρ

ð19Þ

The conditions on the tangential velocity component and on turbulent parameters are based on the theoretical profiles of those parameters in the vicinity of a smooth wall (so called wall laws) [56]   ρu2 u 2 n ∇u þ ∇uT  ∇uU ¼  T ; ð20Þ 3 μ juj n∇K ¼ 0; n∇ω ¼ 

ð21Þ ∂ω uT ¼ pffiffiffiffiffi : ∂δ κ βn δ2

ð22Þ

Here δ is a small distance from the wall. Parameter uT can be determined from the logarithmic velocity profile in the vicinity of the wall.

ð14Þ

Enthalpy in the expression (8) is composed of enthalpies of species as   N N p RG : ð15Þ h ¼ ∑ Y k hk ðTÞ ¼ eþ ¼ ∑ Y k ek ðTÞ þ ρ k¼1 Wk k¼1

ρ ¼ ∑ ρk :

u ¼ 0; p ¼ P 0 ; T ¼ T 0 ; X C 2 H2 : X O2 ¼ 1: x;

juj ¼ uT

  1 ρuT δ ln þ BT ; κ μ

κ ¼ 0:41; BT ¼ 5:2:

ð23Þ

3.3. Thermodynamic and transport models The heat capacity at constant pressure, enthalpy, and entropy of the gas mixture are calculated on the basis of their values for components dependent on temperature T, and partial molar density of each component X k , as follows: N N ρ ρcP ¼ RG ∑ C~ k ðTÞX k ; ρh ¼ RG T ∑ H~ k ðTÞX k ; X k ¼ k : W k k¼1 k¼1

ð24Þ Here C~ k and H~ k are the dimensionless heat capacity and the enthalpy of species per mole. Internal energy and heat

N.N. Smirnov et al. / Acta Astronautica 104 (2014) 134–146

capacity per unit volume are derived from those values as N

N

p ρcV ¼ ρcP  ¼ RG ∑ ðC~ k  1ÞX k ; ρe ¼ ρh  p ¼ RG T ∑ ðH~ k  1ÞX k : T k¼1 k¼1

ð25Þ

The dimensionless values C~ k and H~ k dependence on temperature is assumed polynomial; coefficients are computed from tabular data in order to fit heat capacity with least squares, and to fit enthalpy at the reference temperature. Those dependencies are as follows: 3 3 α α k;j j T þ k;4 : C~ k ¼ ∑ αk;j T j ; H~ k ¼ ∑ T j¼0 j ¼ 0 jþ1

ð26Þ

Thus, 5 constants for each of species αk;j determine the thermodynamic properties of the mixture, including the chemical energy which is “hidden” mainly in the constant αk;4 . The transport model uses the constant Prandtl and Schmidt assumptions in order to simplify details of diffusion and thermal conduction modeling; this does not simplify the molecular viscosity modeling, which is significant in determining transport, both laminar and turbulent. The turbulence is also governed by molecular viscosity due to its role in wall laws, which determine boundary conditions for velocity and turbulent parameters. The viscosity mixture depends on species concentrations X k and components viscosity μk ðTÞ non-linearly [58] μk X k : m k ¼ 1 ∑j ¼ 1 ϕkj X j m

μ¼ ∑

ð27Þ

The matrix ϕkj can be calculated on the basis of μk as follows: 2 32 !1=2    1=2 W j 1=4 5 2W j 14 μk ϕkj ¼ 1 þ : ð28Þ 4 μj Wk Wk þWj The dependence of μk upon T (rather complicated, if found by molecular statistical methods) was fit in [58] to depend on 4 constants for each of the species, so that

 B C μk ¼ 10  7 exp Ak ln T þ k þ 2k þDk : ð29Þ T T

3.4. Kinetic model The full kinetic model of acetylene consists of hundreds of reactions with hundreds of components. On the other hand, acetylene is a high-energetic fuel, and the mixture must contain radical components. The reduced kinetic mechanism incorporates the following parts. 1. Acetylene oxidation: the reaction governing the decay of acetylene is set as [57,59] proposed as C2 H2 þ O2 2CO þ 2H, but its kinetics is changed as to fit with the experimental results [60], which read as C2 H2 þ O2 products; with overall intensity coefficient k ¼ 4:6 1015 T  0:54 expð  45; 000=RTÞ½C 2 H2 ½O2 , 2. CO oxidizing as in the Watanabe–Otaka model [61] as described in [62], with 2 changes: reversibility of CO

139

oxidation, and no hydrogen oxidation (because it is taken into account further) CO þ0:5O2 -CO2 ; k ¼ 3:98  1014 expð  40; 000=RTÞ½CO½O2 0:25 ½H 2 O0:5 CO þH 2 OC-O2 þ H 2 ; k ¼ 2:75  1012 expð  20; 000=RTÞ½CO½H 2 O; 3. Hydrogen, and short radicals exchange and recombination kinetics is taken from Maas and Pope's work [63]. Exchange H þ O2 -OHþ O; k ¼ 2:00  1014 expð  8455=TÞ½H½O2 

O þH 2 -OH þH; k ¼ 5:06  104 T 2:67 expð  3163=TÞ½O½H 2  OH þ H 2 -H2 O þH; k ¼ 1:00  108 T 1:6 expð 1660=TÞ½OH½H 2  2OH-H2 Oþ O; k ¼ 1:50  109 T 1:14 expð  48:1=TÞ½OH2 Recombination/dissociation 2H þ M-H2 þ M; k ¼ 1:80  1018 T  1 ½H2 ½M 2O þM-O2 þM; k ¼ 2:90  1017 T  1 ½O2 ½M H þOH þ M-H2 O þM; k ¼ 2:20  1022 T  1 ½H½OH½M ðcm; mol; KÞ: 4. Another recombination/dissociation reaction, which accelerates water vapor release, is taken from Ref. [64] O þH þM-OH þM; k ¼ 4:71  1018 T  1 ½H½O½M ðcm; mol; KÞ: The third-body component effect in recombination and dissociation reactions is based on the Chaperon coefficients from Ref. [65] ½M ¼ 6:5½H 2 O þ 0:4½O2  þ 0:4½N 2  þ 0:75½CO þ 1:5½CO2  þ 3:0½C 2 H2  þ 1:0 ½others

4. Results and discussion Numerical simulations and experimental investigations were performed for two geometrical characteristics of the apparatus shown in Fig. 2, a long tube (1240 mm) with a long spiral (300 mm) and a short tube (910 mm) with a short spiral (150 mm). Fig. 3 shows successive stages of detonation onset in the spiral region of the long tube in the case of the equivalence ratio of fuel to oxidizer being equal to unity. The spiral occupies the zone from 100 to 400 mm. The section of spiral zone from 200 to 300 mm from the left

140

N.N. Smirnov et al. / Acta Astronautica 104 (2014) 134–146

Fig. 3. Successive stages of detonation onset in the obstacle region of a long tube for the case equivalence ratio of fuel to oxidizer is equal to unity. Upper part – pressure and lower part – CO2 concentration.

edge is shown in Fig. 3. Upper part of each color map illustrates the pressure field, and lower part illustrates concentration of CO2. The last parameter traces the flame front fairly well because its concentration grows higher as the cold mixture ignites but then carbon monoxide decomposes into COþO2 due to high temperature behind the flame front.

As the flame propagates, the pressure grows in the system, and the flame accelerates (compare times 0.808, 0.823 and 0.838 ms). The pressure field shows numerous shocks produced by the spiral elements. Transition to detonation takes place between 0.846 and 0.848 ms in the vicinity of the axis of the system, after the flame ignites the unburned pocket of gas surrounded by hot

N.N. Smirnov et al. / Acta Astronautica 104 (2014) 134–146

reaction products as seen in this location at 0.838 and 844 ms. This coincides with conclusions of [48] on detonation onset. First the shape of the detonation wave

141

is almost hemi-spherical (0.848 ms), and then it flattens (0.857 ms). The level of pressure behind the wave is above 20 bar.

Fig.4. Typical pressure–time diagrams for various cases of detonation onset within the spiral part of the system. First curve computed at 600 mm and second at 700 mm: (a) long tube, equivalence ratio 1.0; (b) long tube, equivalence ratio 2.5; (c) short tube, equivalence ratio 1.0; (d) short tube, equivalence ratio 2.5.

Fig. 5. Pressure history at 600 and 700 mm markers, deflagration case: (a) long tube, equivalence ratio 4.17; (b) short tube, equivalence ratio 4.17.

142

N.N. Smirnov et al. / Acta Astronautica 104 (2014) 134–146

Typical pressure–time diagrams for cross-section averaged pressure are shown in Fig. 4 for the case in which detonation onset took place. Time difference between two

curves calculated at distances of 600 mm and 700 mm from the beginning of the tube was used to estimate detonation velocity both in numerical simulations and in experiments.

Fig. 6. Pressure history at 600 and 700 mm markers, late onset of detonation: (a) long tube, equivalence ratio 0.625; (b) short tube, equivalence ratio 0.625.

Fig. 7. Flame velocity evolution versus tube length: (a) long tube, equivalence ratio 1.0; (b) long tube, equivalence ratio 2.5; (c) short tube, equivalence ratio 1.0; (d) short tube, equivalence ratio 2.5.

N.N. Smirnov et al. / Acta Astronautica 104 (2014) 134–146

143

Fig. 8. Flame velocity evolution versus tube length: (a) long tube, equivalence ratio 0.625; (b) short tube, equivalence ratio 0.625.

Fig. 9. Flame velocity evolution versus tube length: (a) long tube, equivalence ratio 4.17; (b) short tube, equivalence ratio 4.17.

Typical curves in Fig. 5 illustrate the deflagration case (rich mixture, equivalence ratio 2.5/0.6¼4.17). Under these conditions onset of detonation does not take place neither in the spiral section, nor after it in tube being too short for it. Late onset of detonation (beyond 600 and 700 mm markers) in lean mixture (equivalence ratio 2.5/4.0¼ 0.625) is shown in Fig. 6. One can see low initial pressure growth followed by rapid growth due to retonation and reflection waves passing near those locations. In the case of the shorter tube, there was a galloping defagration mode, which could bring to detonation onset finally but the length of the system prevented it. Evidently pressure history at those locations 600 and 700 mm cannot help estimating flame velocity in the case of deflagration or late detonation onset. We used the other method – flame position was estimated as the rightmost location of reaction products; then this curve was smoothed using the robust loss (locally weighted scatter plot smooth) method using MATLAB, and differentiated. Typical results for velocity evolution versus tube length are

shown in Figs. 7–9. Fig. 7 illustrates typical detonation cases when the detonation onset is either in the spiral section, or a small instance beyond it. Then the typical deflagration cases (lean mixture) are shown in Fig. 8. One can see from Fig. 8 that flame is accelerating; features of galloping regime are seen on the shorter tube. Velocity increases by the end of the spiral section in both cases, then decreases, and rises again, for longer tube more constantly and for shorter tube with some jumps (due to pressure oscillations invoked in the spiral region, shorter for the 2nd case). Flame velocity is somewhat less than that of the leading shock. Typical cases of later detonation onset (or approaching onset in the shorter tube case) are shown in Fig. 9. Galloping acceleration of flame in rich mixture (equivalence ratio 4.17, Fig. 9) finally leads to detonation onset in the longer tube (left) of to approaching detonation onset (right, shorter tube). The results on velocity variation of flame front (detonation wave, or turbulent flame) for different mixture equivalence ratios are present in Fig. 10. Solid curve illustrates

144

N.N. Smirnov et al. / Acta Astronautica 104 (2014) 134–146

Fig. 10. Combustion wave velocity variation versus mixture equivalence ratio: simulations, experiments and CJ detonation.

CJ detonation velocity; calculations are based on the assumption of chemical equilibrium being established instantly behind the leading shock; the method is described in [66]. The NASA CEA program was used to compute chemical equilibrium [67]. The species used to compute the composition were C2H2, CO, CO2, O2, H2, H, O, OH, and H2O. Dashed curve with square markers – results for long tube (1240 mm) with long spiral (300 mm). Markers show the calculated cases. Dashed curve with diamond markers – results of calculation for short tube (910 mm) with short spiral (150 mm). Triangles turned up and down correspond to experimental data measured in long and short tubes respectively. It is seen from Fig. 10 that for lean mixtures (lower equivalence ratio) the flame velocity in short tube is higher in cases when the calculated curve does not correspond to detonation but traces accelerated deflagration. Flame velocity here is taken as average on the base of 500–800 mm from the system left edge. The flame velocity shown on the figure in the case of detonation is calculated as in the experiments; by time difference between rapid pressure increase at 600 and 700 mm locations from the left edge of the system. Numerical simulations show maximal velocity in the range of equivalence ratios from 1.0 up to 3.0. Fig. 11 illustrates peak pressure at location 600 mm as a function of equivalence ratio. Solid curve illustrates pressure behind CJ detonation wave calculated in accordance with [66,67], curve with diamond markers – shorter tube with short spiral, and curve with square markers – longer tube with longer spiral. It is seen that near stoichiometry the calculated curves lay near the theoretical one (typically higher, however, due to slightly overdriven detonation; CJ detonation being theoretical limit for an overdriven one). Closer to concentration limits, peak pressure is usually much higher than the theoretical curve, because the detonation onset takes place nearer to the location 600 mm from the system left edge, and the pressure peak is much higher near the spontaneous detonation onset location than the peak predicted by CJ wave calculations based on thermodynamics only. Beyond the limits onset of detonation does not take place within the

Fig. 11. Maximal pressure at location 600 mm calculated in long (squares) and short (diamonds) tubes for different equivalence ratios. CJ detonation pressure is illustrated by solid curve.

section of the tube under investigation and combustion velocity, as well as shock wave pressure, are much lower than that provided by CJ calculations. To provide comparison with experimental data average pressure was calculated at the location 600 mm from the system left edge. The average pressure in calculations was obtained at the wall. To obtain the average data, we integrated it upon a set of time intervals where p1 4 p0 , where p0 is the initial pressure in the system, and p1 ðtÞ is the pressure in the given location as a function of time. The integration could be processed as follows: ( ) R t calc p1 ðtÞ; p1 4p0 dt 0 0; p1 rp0 ( ) pave ¼ : ð30Þ R t calc 1; p1 4p0 dt 0 0; p1 r p0

It is seen from formula (30) that the detonation Von Neuman pressure peak despite of being very narrow also contributes to the average pressure. However, inertial pressure transducers could hardly detect it. Due to that reason for comparison of numerical and experimental results we excluded Von Neuman peaks out of consideration by using the following averaging strategy: ( ) R t calc p1 ðtÞ; p0 op1 o pave þ pd dt 0 0; p1 rp0 or p1 Z pave þ pd ð1Þ ( ) pave ¼ ; ð31Þ R t calc 1; p0 op1 o pave þ pd dt 0 0; p1 rp0 or p1 Z pave þ pd where pd is a model parameter characterizing height of cutting peaks above average values (1 bar for present averaging). The averaged pressure provided by formula (31) was compared with averaged pressure measured in experiments. Results of comparison being function of equivalence ratio are shown in Fig. 12. The results show that in general the qualitative and quantitative agreement of experimental and simulation

N.N. Smirnov et al. / Acta Astronautica 104 (2014) 134–146

145

Acknowledgements The authors gratefully acknowledge the financial support from the Russian Foundation for Basic Research (Grant 13-03-00003) and the National Research Foundation of Korea (NRF)-RFBR Joint Research Program (Grant 12-08-91702). References

Fig. 12. Average pressure rise obtained at 600 mm location in numerical simulations (dashed lines) and experimental measurements (triangles).

data is rather good. However, close to stoichiometry numerical simulations provide a higher value of pressure than experimental data, while for equivalence ratio 2.5, which corresponds to a highest CJ pressure (Fig. 11), the coincidence of numerical and experimental data in terms of average pressures is very good. Numerical simulations provide higher values for pressure rise at the spot 600 mm for long tube with long spiral than that obtained in short tube with short spiral for the same spot. This difference manifests in the equivalence ratios range: 1:0 o ϕo 2:5. This could be explained by a faster detonation onset in tube with a longer spiral. In experimental data one cannot detect such a big difference in pressures measured in long and short tubes.

5. Conclusions Deflagration to detonation transition in acetylene mixed with oxygen in various proportions is obtained in a cylindrical tube with a section of obstacles modeling Shchelkin spiral both experimentally and in numerical simulations. Detailed numerical simulation showed that detonation onset takes place in the unburned pocket of gas surrounded by hot reaction products as soon as ignition occurs in this pocket. Depending on the mixture equivalence ratio different deflagrations to detonation transition scenarios happened. For mixture equivalence ratios 1:0 o ϕ o2:5, onset of detonation took place in the zone of obstacles, and then the detonation wave propagated along the tube. For lean mixtures accelerated turbulent flame was observed. For rich mixtures ðϕ ffi4:0Þ a galloping combustion regime was observed, which resulted the onset of detonation in longer tube, while in a shorter one it did not take place due to predetonation distance limitation. Main detonation characteristics, measured experimentally, developed in numerical simulations, or calculated directly based on thermodynamic analysis are in good agreement for the equivalence ratio interval 1:0 o ϕo 3:5.

[1] D.F. Helman, R.P. Shreeve, S. Eidelman, Pulsed Detonation Engine, AIAA Paper No. 86-1683,1986. [2] J.A. Nicholls, H.P. Wilkinson, R.B. Morrison, Intermittent detonation as a thrust producing mechanism, Jet Propuls. 27 (5) (1957) 534–541. [3] N.N. Smirnov, A.P. Boichenko, Deflagration to detonation transition in gasoline–air mixtures, Combust. Explos. Shock Waves 22 (2) (1986) 65–67. [4] S. Eidelman, W. Grossmann, Pulsed Detonation Engine – Experimental and Theoretical Review, AIAA Paper No. 92-3168, 1992. [5] G. Roy, S.M. Frolov, K. Kailasanath, N.N. Smirnov (Eds.), Advances in Experimentation and Computation of Detonations, Moscow ENAS Publ., 1998. [6] G.D. Roy, S. Frolov, K. Kailasanath, N. Smirnov (Eds.), Gaseous and Heterogeneous Detonations: Science to Applications, ENAS Publ., Moscow, 1999. [7] N.N. Smirnov, V.F. Nikitin, Modeling and simulation of hydrogen combustion in engines, Int. J. Hydrog. Energy 39 (2) (2014) 1122–1136. [8] F. Schauer, J. Stutrud, R. Bradley, V. Katta, J. Hoke, Detonation studies and performance results for a research pulse detonation engine, in: G. Roy, et al., (Eds.), Confined Detonations and Pulse Detonation Engines, Torus Press, Moscow, 2003, pp. 287–302. [9] C.M. Brophy, J.O. Sinibaldi, D. Netzer, K. Kailasanath, Initiator diffraction limits in a pulse detonation engine, in: G. Roy, et al., (Eds.), Confined Detonations and Pulse Detonation Engines, Torus Press, Moscow, 2003, pp. 59–72. [10] V.G. Alexandrov, A.N. Kraiko, K.S. Reent, Mathematical model of a supersonic pulsed detonation Ramjet engine, in: G. Roy, et al., (Eds.), High Speed Deflagration and Detonation: Fundamentals and Control, ELEX-KM Publ., Moscow, 2001, pp. 333–344. [11] A.A. Vasil'ev, About a detonation engine with external combustion, in: G. Roy, et al., (Eds.), High Speed Deflagration and Detonation: Fundamentals and Control, ELEX-KM Publ., Moscow, 2001, pp. 303–314. [12] F.A. Bykovskii, S.A. Zhdan, E.F. Vedernikov, Continuous detonation in the regime of self-oscillating oxidant supply, Combust. Explos. Shock Waves 47 (2) (2011) 101–111. [13] V.A. Levin, J.N. Nechaev, A.I. Tarasov, A new approach to organizing operating cycles in pulsed detonation engines, in: G. Roy, et al., (Eds.), High Speed Deflagration and Detonation: Fundamentals and Control, ELEX-KM Publ., Moscow, 2001, pp. 223–238. [14] N.N. Smirnov, Pulse detonation engines: advantages and limitations, in: N. Syred, A. Khalatov (Eds.), Advanced Combustion and Aerothermal Technologies, Springer, 2007. [15] V.P. Korobeinikov, V.V. Markov, I.V. Semenov, D.D. Pedrow, S. Wojcicki, On modeling of a pulse detonation engine, in: G.D. Roy, et al., (Eds.), High Speed Deflagration and Detonation, ELEKS-KM Publishers Moscow, 2001, p. 289. [16] P. Wolanski, Detonative propulsion, Proc. Combust. Inst. 34 (2013) 125–158. [17] Yu.G. Phylippov, V.R. Dushin, V.F. Nikitin, V.A. Nerchenko, N.V. Korolkova, V.M. Guendugov, Fluid mechanics of pulse detonation thrusters, Acta Astronaut. 76 (2012) 115–126. [18] E. Mallard, H.L. Le Chatelier, Comptes Rendus Acad. Sci. Paris 93 (1883) 145. [19] M. Berthélot, P. Vieille, Comptes Rendus Acad. Sci. Paris 93 (1883) 18. [20] V.A. Mikhelson, On normal combustion velocity of explosive gaseous mixtures, Imp. Mosc. Univ. Sci. Bull.: Phys. Math. Ser. 10 (1893) 1–92. (in Russian). [21] W.J.M. Rankine, Philos. Trans. (1870) 277–288. [22] H. Hugoniot, J. Liouville 3 (1887) 477–492; H. Hugoniot, J. Liouville 4 (1888) 153–167. [23] D.L. Chapman, Philos. Mag. 47 (1899) 90. [24] E.J. Jouguet, Mathematics (1905) 347. [25] G.D. Salamandra, T.V. Bazhenova, I.M. Naboko, J. Theor. Phys. 29 (11) (1959) 1354.

146

N.N. Smirnov et al. / Acta Astronautica 104 (2014) 134–146

[26] A.K. Oppenheim, P.A. Urtiew, Experimental observations of the transition to detonation in an explosive gas, Proc. R. Soc. A 295 (1966) 13. [27] R.I. Soloukhin, Methods of Measure and Main Results of Experiments in Shock Tubes, Novosibirsk State University Publ., Novosibirsk, 1969. [28] A.K. Oppenheim, R.I. Soloukhin, Annu. Rev. Fluid Mech. 5 (1973) 31. [29] V.F. Nikitin, V.R. Dushin, Y.G. Phylippov, J.C. Legros, Pulse detonation engines: technical approaches, Acta Astronaut. 64 (2009) 281–287. [30] Y. Wang, J. Wang, Y. Li, Y. Li, Induction for multiple rotating detonation waves in the hydrogen–oxygen mixture with tangential flow, Int. J. Hydrog. Energy (2014), http://dx.doi.org/10.1016/j. ijhydene.2014.05.162. [31] N.N. Smirnov, M.V. Tyurnikov, Experimental Investigation of deflagration to detonation transition hydrocarbon–air gaseous mixtures, Combust. Flame 100 (1995) 661–668. [32] Ya.B. Zeldovich, V.B. Librovich, G.M. Makhviladze, G.I. Sivashinsky, On the onset of detonation in a non-uniformly pre-heated gas, Sov. J. Appl. Mech. Tech. Phys. 2 (1970) 76. [33] A.G. Merzhanov, On critical conditions for thermal explosion of a hot spot, Combust. Flame 10 (1966) 341–348. [34] A.A. Borisov, On the origin of exothermic centers in gaseous mixtures, Acta Astronaut. 1 (1974) 909–920. [35] K. Kailasanath, E.S. Oran, Ignition of flamelets behind incident shock waves and the transition to detonation, Combust. Sci. Technol. 34 (1983) 345–362. [36] Ya.B. Zeldovich, B.E. Gelfand, S.A. Tsyganov, S.M. Frolov, A.N. Polenov, Concentration and temperature non-uniformities of combustible mixture as a reason of pressure waves generation, in: A. Kuhl, et al., (Eds.), Dynamics of Explosions, vol. 114, AIAA Inc., New York, 1988, p. 99. [37] N.N. Smirnov, M.V. Tyurnikov, A study of deflagration and detonation in multiphase hydrocarbon–air mixtures, Combust. Flame 96 (1994) 130–140. [38] P. Wolanski, Arch. Combust. (3–4) (1991) 143–149. [39] N.N. Smirnov, I.I. Panfilov, Deflagration to detonation transition in combustible gas mixtures, Combust. Flame 101 (1995) 91–100. [40] N.N. Smirnov, V.F. Nikitin, A.P. Boichenko, M.V. Tyurnikov, V.V. Baskakov, Deflagration to detonation transition in gases and its application to pulse detonation devices, in: G.D. Roy, et al., (Eds.), Gaseous and Heterogeneous Detonations: Science to Applications, ENAS Publ., Moscow, 1999, pp. 65–94. [41] C.J. Brown, G.O. Thomas, Experimental studies of shock-induced ignition and transition to detonation in ethylene and propane mixtures, Combust. Flame 117 (1999) 861–870. [42] A.M. Khohlov, E.S. Oran, Numerical simulation of detonation initiation in a flame brush: the role of hot spots, Combust. Flame 119 (1999) 400–416. [43] K.I. Shchelkin, Y.a.K. Troshin, Gas Dynamics of Combustion, USSR Academy Sci. Publ., Moscow, 1963. [44] M.V. Silnikov, M.V. Chernyshov, V.N. Uskov, Two-dimensional overexpanded jet flow parameters in supersonic nozzle lip vicinity, Acta Astronaut. 97 (2014) 38–41. [45] M.V. Silnikov, M.V. Chernyshov, V.N. Uskov, Analytical solutions for Prandtl–Meyer wave–oblique shock overtaking interaction, Acta Astronaut. 99 (2014) 175–183. [46] N.N. Smirnov, V.F. Nikitin, The influence of confinement geometry on deflagration to detonation transition in gases, J. Phys. IV Fr. 12 (Pr7) (2002) 341–351. [47] N.N. Smirnov, V.F. Nikitin, S. Alyari Shurekhdeli, Investigation of selfsustaining waves in metastable systems: deflagration-to-detonation transition, J. Propuls. Power 25 (3) (2009) 593–608.

[48] M.A. Liberman, M.F. Ivanov, A.D. Kiverin, M.S. Kuznetsov, Deflagration-to-detonation transition in highly reactive combustible mixtures, Acta Astronaut. 67 (7–8) (2010) 688–701. [49] N.N. Smirnov, V.F. Nikitin, M.V. Tyurnikov, A.P. Boichenko, J.C. Legros, V.M. Shevtsova, Control of detonation onset in combustible gases, in: G.D. Roy, et al., (Eds.), High Speed Deflagration and Detonation, Elex-KM Publ., Moscow, 2001, pp. 3–30. [50] N.N. Smirnov, V.F. Nikitin, Yu.G. Phylippov, Deflagration to detonation transition in gases in tubes with cavities, J. Eng. Phys. Thermophys. 83 (6) (2010) 1287–1316. [51] M. Son, C. Seo, K.W. Lee, J. Koo, N.N. Smirnov, Effect of spiral turbulent ring on detonation performances of acetylene–oxygen mixture, J. Korean Soc. Propuls. Eng. 17 (2) (2013) 9–15. [52] Min Son, Chanwoo Seo, N.N..Smirnov, Jaye Koo, Detonation Characteristics of Acetylene–Oxygen with Additional Liquid Fuel, AIAA, ISABE-2013-1669, 2013, 6pp. [53] D.C. Wilcox, Turbulence Modeling for CFD, DCW Industries, La Cañada, California, 1994. [54] Q. Xiao, H.M. Tsai, F. Liu, Computation of transonic diffuser flows by a lagged k–ω turbulence model, J. Propuls. Power 19 (3) (2003) 168–177. [55] D.C. Wilcox, Turbulence Modeling for CFD, 2nd ed. DCW Industries, La Cañada, California, 1998. [56] H. Grotjans, F. Menter, Wall functions for general application CFD codes, in: Proceedings of the 4th European Computational Fluid Dynamics Conference, ECCOMAS 98, John Wiley & Sons, 1998, pp. 1112–1117. [57] N.M. Marinov, W.J. Pitz, C.K. Westbrook, M. Hori, N. Matsunaga, An experimental and kinetic calculation of the promotion effect of hydrocarbons on the NO–NO2 conversion in a flow reactor, Proc. Combust. Inst. 27 (1998) 389–396. [58] R. Svehla, Transport Coefficients for the NASA Lewis Chemical Equilibrium Program, NASA Technical Memorandum 4647, April 1995. [59] Y. Hidaka, K. Hattori, T. Okuno, K. Inami, T. Abe, T. Koike, Combust. Flame 107 (1996) 401. [60] D.L. Baulch, C.J. Cobos, R.A. Cox, P. Frank, G. Hayman, T.h. Just, J.A. Kerr, T. Murrells, M.J. Pilling, J. Troe, R.W. Walker, J. Warnatz, Evaluated kinetic data for combustion modeling, J. Phys. Chem. Ref. Data 23 (1994) 847. [61] H. Watanabe, M. Otaka, Numerical simulation of coal gasification in entrained flow coal gasifier, Fuel 85 (2006) 1935–1943. [62] O.A. Marzouk, E.D. Huckaby, A comparative study of eight finite-rate chemistry kinetics for CO/H2 combustion, Eng. Appl. Comput. Fluid Mech. 4 (3) (2010) 331–356. [63] U. Maas, S. Pope, Simplifying chemical kinetics: intrinsic lowdimensional manifolds in composition space, Combust. Flame 88 (1992) 239–264. [64] Z. Hong., An improved hydrogen/oxygen mechanism based on shock tube/laser absorption measurements, A Dissertation Submitted to the Department of Mechanical Engineering and the Committee on Graduate Studies of Stanford University, November 2010. [65] V.V. Lissianski, V.M. Zamansky Jr., W.C. Gardiner, Gas-phase combustion chemistry, in: W.C. Gardiner Jr. (Ed.), Chapter Combustion Chemistry Modeling, Springer-Verlag, New York, 2000. [66] S. Browne, J. Ziegler, J.E. Shepherd, Numerical solution methods for shock and detonation jump conditions, 91125 GALCIT Report FM2006.006, Aeronautics and Mechanical Engineering, California Institute of Technology, Pasadena, CA USA, July 2004 (revised 29.08.08). [67] S. Gordon, B.J. McBride, Computer program for calculation of complex chemical equilibrium compositions and applications I, Analysis, NASA RP-1311, October 1994.