## Abstract

This paper presents a simplified reactive multi-gas model for the numerical simulation of detonation waves. The mathematical model is formulated based on a thermodynamically consistent and fully conservative formulation, and is extended to model reactive flow by considering the reactant and product gases as two constituents of the system and modelling the conversion between these by a simple one-step reaction mechanism. This simplified model allows simulations using more appropriate chemico-thermodynamic properties of the combustible mixture and yields close Chapman–Jouguet detonation parameters from detailed chemistry. The governing equations are approximated using a high-resolution finite volume centred scheme in an adaptive mesh refinement code, permitting high-resolution simulations to be performed at flow regions of interest. The algorithm is tested and validated by comparing results to predictions of the one-dimensional linear stability analysis of the steady detonation and through the study of the evolution of two-dimensional cellular detonation waves in gaseous hydrogen-based mixtures.

## 1. Introduction

Owing to the recent interest in the advanced propulsion concept (e.g. the pulse detonation wave engine; Kailasanath 2000), condensed phase explosives modelling (Bdzil & Stewart 2007), safety aspects and assessment of explosion hazards in the hydrogen energy economy for preventing occurrence of deflagration-to-detonation transition (Ng & Lee 2008), research in detonation physics, both experimentally and numerically, is becoming increasingly significant. Over past decades, there have been many advances in detonation modelling and experimental characterization. However, there remain many challenges to fully understand the detonation phenomena, particularly in achieving a quantitative description of the real unstable detonation structure and in predicting dynamic detonation parameters such as the initiation energy and propagation limits (Lee 2008).

A detonation wave is the most violent mode of combustion. It is a combustion-driven compression wave that propagates at supersonic speeds in an explosive medium. Typical detonation speeds in gases are of the order of 1800 m s^{−1}. For a given explosive, there corresponds a unique steady-state detonation velocity. Chapman's hypothesis of a minimum velocity (or tangency) solution, where the Rayleigh line is tangent to the equilibrium Hugoniot curve, or the equivalent sonic condition at the equilibrium plane due to Jouguet, can be used as alternate criteria for selecting the unique solution to the global conservation laws across the detonation wave. This unique solution from Chapman–Jouguet (CJ) theory is found to agree quite well with experimental observations for conditions well within the detonability limits. Although the CJ theory determines the detonation solution without the need to consider non-equilibrium, it cannot provide any information on the initiation requirement that governs the sensitivity of the explosive, nor can it describe the effect of boundary conditions on the propagation of the detonation that will lead to failure (in detonability limits). These questions and the mechanism of detonation instability can only be resolved by considering the detailed structure of the detonation wave itself. Thus, if we were to be able to predict the dynamic parameters or any departure from the CJ theory, the non-equilibrium unstable detonation structure has to be considered (Lee 2008).

Numerical simulations of a detonation wave structure, i.e. a high speed combustion wave that consists of a shock wave travelling at supersonic velocity followed by a chemical reaction zone, deal with strong nonlinear interactions between gas dynamics and chemistry. A full understanding of these coupled phenomena requires a rigorous mathematical model and highly resolved numerical simulations. One- and multi-dimensional detonation structure has been extensively studied numerically in the past by various researchers (Fickett & Wood 1966; Taki & Fujiwara 1978; Oran *et al*. 1982, 1988; Bourlioux & Majda 1991, 1992; Quirk 1994; Gamezo *et al*. 1999*a*,*b*; Sharpe & Falle 1999, 2000; Khokhlov *et al*. 2004). These numerical investigations on detonation structures are mostly carried out by assuming a single global equation of state for simplicity. One needs to select the molecular weight *M* and ratio of specific heats *γ* as average values obtained from detailed thermodynamic calculations and then adjust the energy release *Q* to match the experimentally measured detonation CJ properties of the system. This is, however, a major simplification of the genuine physical situation in most practical cases. For reactive flow, different species are being created or destroyed through chemical reactions, and the multi-gas components could have a significant effect on the thermodynamics of the flow itself. Consider, for example, a stoichiometric reacting mixture of hydrogen and oxygen; even if we assume a direct conversion from reactants to products both modelled by the ideal gas equation of state, i.e. 2H_{2}+O_{2}→2H_{2}O, then the properties of these two mixtures vary notably: for the reactants, *γ*∼1.4, and mean molecular weight *M*=12; while for the products, *γ*∼1.2 and *M*=18. This change in *γ* implies a factor of two difference in the term (*γ*−1) appearing in the equation of state and hence a marked difference in the behaviour of the two gases. In reality, the problem is even more complicated than this since the chemical reaction does not progress as a single step, but involves the conversion of the reactants to various intermediaries before the final state is reached, each of which will also have unique thermodynamic properties. For a simplified condensed-phase combustion model, the variation in the equivalent *γ* between the liquid/solid reactants and gaseous products can also be very pronounced for detonations in condensed phase explosives.

It is well recognized from early stability analyses that the ratio of specific heats *γ* is one of the fundamental parameters underlying the structure of the steady wave, parameters which play a critical role in the detonation dynamics. Hence, as a first stage to better represent the realistic structure and examine the stability properties of detonations, the simplified mathematical model should first take into account the nature of multi-component species resulting in the variation of *γ*. Although a complex set of species with chemical kinetic rate equations derived from elementary reactions could, in principle, be solved simultaneously with the reactive Euler equations within current computational capabilities (Oran *et al*. 1998; Eckett 2000; Tsuboi *et al*. 2002; Hu *et al*. 2004; Deiterding & Bader 2005), direct interpretation of the large amount of detailed information generated by such numerical simulations becomes a difficult problem, and numerical resolution issues remain a challenge. As such, the objective of this paper is to formulate a simplified rigorous multi-component numerical model for reactive flow simulation by considering the reactant and product gases as two constituents of the system, and modelling the conversion between these by a simple one-step reaction mechanism. This model provides one extra level of complexity than previous studies using one single fluid approximation, but at the same time, allows high-resolution simulations to be performed.

Perhaps the simple approach to formulate the governing mathematical model is to augment the Euler equation with the two advection equations for the ratio of specific heats and the mixture molecular weight. However, as shown by Huo & Floch (1994) and Karni (1996), among others, implementing this system within a conservative formulation commonly leads to problems in the numerical solutions, especially those associated with contact discontinuities and shock waves. The difference between the ratios of the specific heats of the two gases (reactant and products) gives rise to numerical artefacts (generation of spurious waves at their interface which propagate through the solution) (Wang *et al*. 2004; Bates *et al*. 2007). The numerical errors typically appear in both the predicted density and pressure, with those in the pressure particularly severe. These errors are often critical when considering reactive flow, since any perturbations in temperature may be amplified by the exponential temperature dependence of the reactive rate term. Various modifications to the numerical scheme have been suggested in order to correct these errors. Karni (1994) suggested a non-conservative integration scheme, while Abgrall (1996) introduced a quasi-conservative method, both of which reduced the errors. However, the non-conservative formulation causes the usual problems with predicting the locations of shocks, and the quasi-conservative approach, while accurately predicting shock locations, does not entirely remove the problems.

In this study, the reactive flow model for detonation simulations is obtained using a modified formulation, which reduces the pressure fluctuations at the material interface, by Wang *et al*. (2004). This formulation is fully conservative throughout the domain and hence shock speeds are predicted correctly. Also, since no special treatment is required at the material interface, this thermodynamically consistent and fully conservative (TCFC) model can be implemented using any available conservative solver (Bates 2006; Bates *et al*. 2007). For the non-reactive inert case, the TCFC model supplements the system of Euler equations with two advection equations given by(1.1)where *ξ*=*γ*/(*γ*−1) and *M* is the local molecular weight. In this study, the TCFC model is extended to include the source term for reactive fluid flow. To validate the present reactive multi-component model, calculations of the one-dimensional unsteady detonation propagation are carried out, and results are compared with the neutral stability boundary given by the linear stability analysis with the variation of *γ*. A series of two-dimensional simulations are also performed to model and investigate the detonation structures and evolutions in hydrogen-based mixtures.

## 2. Binary reactive gas model formulation

Neglecting the diffusive transport and viscous terms, the fluid dynamics of a reactive flow is generally given by the reactive Euler equations(2.1)where *ρ* is the mass density of the mixture; ** V** is the velocity vector;

*E*is the energy density;

*p*is the pressure; and

*λ*is the mass fraction of reactant. The parameters in the kinetic rate equation are:

*A*the pre-exponential constant;

*E*

_{a}the activation energy; and

*T*the temperature. The single-step reaction model assumes

*R*→

*P*, where

*R*is a stoichiometric mixture of all species present before reaction occurs and

*P*are the combustion products. The presence of all the intermediate species formed during the reaction progress is neglected.

To formulate a viable model for reaction within a multi-gas system, the above reactive Euler equations with the single-step reaction mechanism is combined with the TCFC multi-gas formulation in such a way as to account for a conversion between a mixture of products and the mixture of reactants as the reaction progresses. We also assume that each of these two mixtures obey the ideal gas equation of state, albeit with different gas parameters,(2.2)and hence are parametrized by the values of *γ*, the ratio of specific heats, and *M*, the molecular weight. The mixture in any volume of gas consists of a spatially varying mixture of these two components; this may be expressed by the scalar field of the local mole fraction of reaction *X*. Equivalently, the mole fraction of the product is determined as (1−*X*).

To extend the TCFC multi-gas formulation for reactive flow, it is required to determine the reactive source term within the TCFC model. We can express the mole fraction *X* of reactants from the mixture rules by(2.3)where *ξ*=*γ*/(*γ*−1). The source terms for the variables 1/*M* and *ξ*/*M* due to the chemical reaction within the TCFC model can be obtained using the application of the chain rule as(2.4)Instead of expressing the above source terms in mole fraction *X*, it is more useful to obtain an expression in term of the mass fraction *λ* of the reactant where(2.5)Taking the derivative of equations (2.3) and (2.5) yields(2.6)Using equations (2.3)–(2.6) and after some simple mathematical manipulation, the reactive source terms can be expressed as(2.7)Hence, for completeness, the full three-dimensional system of governing equations is obtained by appending these source terms to the right-hand side of the TCFC system and combining with the reactive Euler equations. This can be expressed in vector form as(2.8)

## 3. Numerical method

Because the above extended system of governing equations is written in conservative form with source terms, its solution is approximated in the present study using a conservative finite volume centred scheme, namely the slope limiter centred (SLIC) scheme (Toro 1999; Toro & Billett 2000), in an adaptive mesh refinement code (CAMR), which allows high-resolution simulations to be performed on desktop computers. The mesh is dynamically refined around shocks, flame fronts and in regions of large gradient of density (Berger & Colella 1989). The SLIC method belongs to the high-resolution class of methods, which is conservative, explicit and second-order accurate in space and time. The SLIC method employs slope-limiting in the region of solution discontinuities to prevent the formation of spurious oscillations, which occur typically for unlimited second-order schemes. Centred schemes do not require information about the characteristic structure of the hyperbolic equation system to be provided; therefore, they generally have a lower computational cost and simple structure compared with most Riemann solver-based methods, although the latter typically have the advantage when it comes to accuracy. Therefore, they can be applied very easily to any hyperbolic system which is in flux-conservative form (Anile *et al*. 2000). The algorithm is very convenient when computing the numerical solution of reactive flow with a system of hyperbolic conservation laws containing advection equations of multi-component mixture gas. Strang's fractional step operator splitting approach is employed in the process of solving the full multi-dimensional system of equations with reactive source terms. The original numerical CAMR code has been used extensively and is well-validated from our previous detonation simulation studies (Ng *et al*. 2005*a*,*b*).

## 4. Results and discussion

### (a) One-dimensional detonation instability and neutral stability boundary

Since the pioneering work of Erpenbeck (1964), Fickett & Wood (1966) and Lee & Stewart (1990), many investigations have been carried out to study the one-dimensional pulsating detonation as a first stage in understanding detonation dynamics, both theoretically using formal linear stability analysis and computationally with high-resolution numerical simulations. In addition, the neutral stability boundary given by the results of the linear stability analysis for longitudinal perturbations has long been used to test numerical schemes (e.g. Bourlioux & Majda 1991). In general, on the stable side of the stability boundary, initial perturbations caused by the numerics are quickly damped out and the wave reverts to the Zeldovich–Neumann–Döring (ZND) solution. For conditions above the stability limit, oscillations initially grow and saturate to a limit-cycle behaviour. Near the boundary, the pulsations are regular and harmonic, but become nonlinear and chaotic further away from the boundary. Richer nonlinear dynamics have also been recently observed through numerical simulations when more complex chemical kinetic mechanisms or other equations of state for condensed phase detonations are considered (Gorchkov *et al*. 2007; Short *et al*. 2008).

A formal normal-mode linear stability analysis of planar detonations with arbitrary equations of state, such as the consideration of multi-component fluid mixtures with temperature-dependent thermo-chemical properties, has recently been formulated by Gorchkov *et al*. (2007). Their results, thus, can provide some canonical values to test the present numerical formulation. Figure 1 shows predictions of the linear stability analysis obtained with the *E*(*v*, *p*) used in the present study as given in equation (2.2). This result is obtained for a particular case in Gorchkov *et al*. (2007), i.e. with the parameter in their formulation, for a two-component mixture with non-dimensional heat release *Q*=50, specific heat ratio *γ*_{reactant}=1.4 and varying *γ*_{product} and non-dimensional activation energies *E*_{a}. These parameters have been made dimensionless by reference to the uniform unburned state ahead of the detonation front. To compare with the numerical results, we choose a particular result of *γ*_{product}=1.25. From the linear stability analysis, the one-dimensional neutral boundary is found to be at .

One-dimensional numerical simulations using the present multi-gas reactive flow model and numerical methods are carried out to compare with this prediction of the linear stability analysis of the steady detonation by Gorchkov *et al*. (2007). To begin the computation and look at the transient development of the detonation wave, an arbitrary initial condition is considered using a strong blast initiation method. The initial transient is highly overdriven to suppress the instability and this indeed allows easier identification of the stability boundary, above which instability will manifest asymptotically. For all simulations, an effective numerical resolution of 128 points per half-reaction zone length of the steady ZND detonation *L*_{1/2} is used to ensure the detailed features of the pulsating front are properly resolved.

Figure 2 shows the von Neumann pressure normalized by the steady ZND solution versus shock position normalized by the half-reaction zone length for different activation energies. From these plots, there appears to be good agreement between the numerical simulations and the linear stability predictions. Below the neutral stability boundary obtained from the linear stability analysis (i.e. *E*_{a}=52.0), a steady-state detonation is established asymptotically from the strong blast wave, as shown in figure 2*a*. Any initial perturbation caused by the initiation would indeed be damped out subsequently, and a steady-state detonation would again be established. As shown in figure 2*b*, for *E*_{a}=52.1, just slightly above the neutral stability limit, the instability would start to develop when the blast wave decays near the CJ conditions. Oscillations continue to grow and it eventually reaches a saturated mode. Further above the stability limit (i.e. *E*_{a}=52.5 in figure 2*c*), the amplitude of the oscillation becomes larger.

The nonlinear evolution away from the neutral stability boundary is illustrated in figure 3. It is shown that the results are similar to those previously observed for one-dimensional pulsating detonations in the idealized polytropic model with single adiabatic index *γ*, i.e. the transition from harmonic to nonlinear chaotic behaviour is through a series of bifurcations (Ng *et al*. 2005*a*). From these results, the two-component numerical model correctly elucidates both the linear and nonlinear behaviours of the one-dimensional unsteady detonation propagation.

### (b) Detonation structure in H_{2}–air and H_{2}–O_{2}–Ar mixtures

It is well established that multi-dimensional detonation waves generally exhibit a complex and unsteady reaction zone structure. Two-dimensional waves are characterized by an ensemble of interacting transverse waves sweeping laterally across the leading shock front of the detonation wave. The interaction of incident shocks, Mach stems and transverse waves form a characteristic cellular pattern, producing so-called detonation cells. In the present study, we conducted two-dimensional computations by solving the two-dimensional augmented Euler equations to validate the present multi-fluid algorithm. We are particularly interested in the gaseous hydrogen-based mixtures.

The computational geometry and details are shown in figure 4. The adaptive grid method is used to improve the accuracy of numerical solutions behind the detonation shock front and inside the reaction zone structure with steep pressure and density gradients or regions with strong gas dynamic fluctuations. Computations were performed in a moving domain with a base grid of 125×75. Continuous flow is coming from the right-hand boundary at a velocity *D*_{CJ} and leaving at the other end. Following Gamezo *et al*. (1999*a*), the boundary condition on the left is an extrapolated outflow defined for every boundary value *Y*_{b}(*Y*=*ρ*, ** V**,

*p*,

*E*,

*γ*,

*M*) as

*Y*

_{b}=

*Y*

_{1}(1−

*r*)+

*Y*

_{e}

*r*, where

*r*is the relaxation rate coefficient having a value of 0.05.

*Y*

_{1}is the current value in the first cell and

*Y*

_{e}is the extrapolation limit, which is chosen as the ambient fluid parameters (except

*λ*=1,

*M*=

*M*

_{product}and

*γ*=

*γ*

_{product}). Reflective boundary conditions are used along the top and bottom boundaries. In all cases, the simulations were initiated by placing a shocked region of pressure

*P*

_{CJ}in the domain. To accelerate the formation of cellular structure of the detonation, a small curvature is introduced at the contact surface. For the simulations, three levels of grid refinement are used (2, 2, 4), and the numerical resolution in the highest level contains approximately 24 grid points per steady ZND half-reaction zone length

*L*

_{1/2}(i.e. where the mass fraction of reactant equals 0.5). As pointed out by Sharpe (2001), it is important to have resolutions of at least 20 points per

*L*

_{1/2}to capture qualitatively the physical features of the detonation reactive flow. An even higher resolution may be needed to capture correctly the fine details of transverse wave structures. Nevertheless, a resolution study was carried out in this study and results indicate that simulations using a grid level higher than 24 points per

*L*

_{1/2}(for example, 32/

*L*

_{1/2}) does not change significantly the salient detonation flow features for the chosen mixtures in this study.

Input physical parameters of the combustible system for the simulations are given in table 1. The material and chemistry parameters of the model were obtained using a detailed chemistry calculation for a real stoichiometric H_{2}–O_{2} mixture with and without argon dilution, and a stoichiometric H_{2}–air mixture (Schultz & Shepherd 2000). Using the binary mixture formulation, the CJ detonation velocity for the multi-gas case can be readily determined as (Lee & Guirao 1981)(4.1)where(4.2)and the CJ properties can be obtained using the reactive Rankine–Hugoniot relationships(4.3)The pre-exponential constant *A* in the single-step Arrhenius model is chosen to approximate the steady ZND reaction zone structure obtained using detailed chemistry and such that the half-reaction zone length approximates roughly the ZND induction zone length (table 2).

Before carrying out the unsteady simulation, it is of significance to compute first the ideal one-dimensional ZND structure using the present multi-component model to look at the characteristics of the steady reaction zone. In fact, it is very important for any simplified model to describe correctly the ZND structure because it forms the mathematical basis for any linear and nonlinear stability analysis and it provides important information with regards to the reaction zone structure and any chemical length scales. Properly modelling the ZND structure represents the first stage in understanding the complex dynamics of real detonations. A comparison of temperature profiles calculated using the simple binary reactive gas model with single-step Arrhenius kinetics and the full detailed chemistry model is given in figure 5. One can see that the von Neumann (VN) and CJ states are well approximated with the present multi-gas model. Without using a multi-component model formulation, one must estimate average values for both the heat release *Q* and the specific heat ratio *γ* and it is clear that the average values cannot give correctly both the solutions of the VN and CJ states. The present multi-component formulation which takes into account the thermodynamic properties (i.e. the specific heat ratio and molecular weight) of both reactant and products mixtures allows one to model better these two states in the detonation structure, as illustrated in figure 5. The post-shock conditions and the equilibrium CJ states are important because of their influence on the stability of the detonation structure (Lee & Stewart 1990).

As a first step in developing a TCFC reactive flow model for multi-component computations in this work, only a single-step Arrhenius reaction is considered; this allows simulations with high numerical resolution to be performed. It becomes clear that the induction zone region cannot be approximated using single-step Arrhenius kinetics (Ng *et al*. 2005*b*). Significant improvements can be obtained by using a two-step or three-step chain-branching kinetic model (Short & Quirk 1997; Clifford *et al*. 1998; Sichel *et al*. 2002; Ng & Lee 2003; Short & Sharpe 2003; Liang & Bauwens 2005). This is outside the scope of the present study. It will be considered in a future paper to include two- or three-step chemical reaction kinetics to improve the present TCFC multi-component reactive flow model.

The numerical results from the unsteady simulations showing the unstable detonation structure in different gaseous mixtures are now discussed. The ‘numerical smoked foils’ showing the time-integrated maximum pressure contour from the numerical simulations, which corresponds to the trajectories of the triple shock interactions (i.e. triple points), are shown in figure 6 for the three cases considered (table 2). Early in the cell formation process, the cell patterns are rather regular in all three cases. However, as the detonation wave continue to propagates, the results show a transition from weakly unstable (displaying regular detonation cells) in the stoichiometric H_{2}–O_{2} mixture with and without argon dilution, to a strongly unstable or irregular cellular structure in the stoichiometric H_{2}–air mixture. These results are in accord with previous experimental and numerical studies (Gamezo *et al*. 1999*a*,*b*).

The present simulation correctly modelled the chemical and thermodynamics properties (i.e. specific heat ratio and molecular weight) of the reactants and products as used in detailed chemistry calculations. As suggested from previous studies, the effect of the equation of state and the influence of argon dilution on the specific heat ratio account mainly for, and could be the main mechanism explaining, the unstable detonation front behaviours (Khokhlov *et al*. 2004). For stoichiometric H_{2}–O_{2} without dilution (see table 1), a smaller value of *γ*_{react} leads to a smaller shock temperature *T*_{s}, and thus a larger reduced activation energy *E*_{a}/*RT*_{s}. In such cases, instabilities will manifest due to the high temperature sensitivity and an irregular cellular structure will result. With increasing argon dilution, a higher *γ*_{react} leads to a higher shock temperature and a weaker coupling between chemical reactions and fluid dynamics, hence leading to a more stable or regular detonation front. The difference in shock temperature for H_{2}–O_{2} with and without argon dilution can be clearly seen in figure 5 (i.e. the results of the solid and dotted curves at *x*=0). For fuel–air mixtures such as the stoichiometric H_{2}–air mixture, the activation energy *E*_{a} is higher compared with the fuel–oxygen case (see table 1) and hence results in a more unstable detonation front. These observations are equivalent to those from one-dimensional simulations with varying *E*_{a} shown earlier in this section (see figures 2 and 3).

To look at the flow field behind the detonation front in more detail, figure 7 shows snapshots of the Schlieren-type density gradient plots and the reactant mass fraction plots. In all three cases, the contact surfaces and reaction fronts are properly simulated without evidence of spurious oscillations or numerical noise. These results, therefore, indicate the capability of the present numerical multi-fluid algorithm for producing highly complex simulations of multi-dimensional detonation cell phenomena.

For unstable detonations, as in the case of the stoichiometric H_{2}–air mixture, the appearance of unburned pockets of reactants and the compressible turbulent nature of the flow are quite apparent in figure 7*c*. Recent experimental investigation also indicates that a unique characteristic of a cellular detonation structure is the keystone feature, which arises from the reaction rate sensitivity to shock conditions (Pintgen *et al*. 2003). Such a keystone feature is also clearly and sharply revealed within the framework of the multi-component reactive flow model with single-step Arrhenius chemistry, as shown in the mass fraction plots of figure 7.

## 5. Conclusion

In this work, an improved reactive binary-mixture model for combustion is considered. It is based on the TCFC formulation, which is extended in the present study to include reactive source terms to model unsteady detonation wave propagations. The benefit of using the binary-mixture formulation is that the individual thermodynamic properties of the combustion reactants and products can be considered without performing any averaging and matching for the computation. Here, the augmented Euler equations are solved using a total variation diminishing centred scheme with adaptive mesh refinement for high-resolution simulations. To assess the validity of the present numerical formulation, it is shown that the one-dimensional neutral stability boundary obtained from the numerical simulation agrees well with that predicted by the normal-mode linear stability analysis. Results showing the nonlinear behaviour away from the stability boundary are consistent with those previously reported in the literature. From the numerical examples of cellular detonation propagation in H_{2}–air and H_{2}–O_{2}–Ar mixtures, the present two-component model correctly simulates all the salient feature of unstable cellular detonation structure. The model prevents spurious waves at material interfaces, which have a particular detrimental effect for reactive flow simulation. Therefore, the present reactive TCFC model is proven to be useful to explore, in more detail, the effect of the variation of the specific heat ratio *γ* and molecular weight *M* within the reaction zone on the dynamic structure of unstable cellular detonations. This model formulation can also be further modified and extended in the future to include multi-step reaction kinetic mechanisms. Although in this paper, only simulation of the propagation of a gaseous detonation was focused on, the present multi-component formulation can also be applied to other high-speed reactive flow studies, ranging from the fundamentals of compressible turbulent flow and chemical processes interaction to applied explosion safety, propulsion or propellant design.

## Acknowledgments

K.R.B. was funded by an EPSRC Doctoral Training Award. H.D.N. would like to acknowledge the support by the Natural Sciences and Engineering Research Council of Canada. G.C. was supported by the CESAer internship programme.

## Footnotes

- Received September 10, 2008.
- Accepted March 17, 2009.

- © 2009 The Royal Society