- Email: [email protected]
Numerical study of cellular detonation structures of methane mixtures
Journal of Loss Prevention in the Process Industries 36 (2015) 394e403
Contents lists available at ScienceDirect
Journal of Loss Prevention in the Process Industries journal homepage: www.elsevier.com/locate/jlp
Numerical study of cellular detonation structures of methane mixtures Anatoliy V. Trotsyuk*, Pavel A. Fomin, Anatoly A. Vasil'ev Lavrentyev Institute of Hydrodynamics Siberian Branch of the Russian Academy of Science (LIH SB RAS), 15 Ac. Lavrentyev Ave., Novosibirsk 630090, Russia
a r t i c l e i n f o
a b s t r a c t
Article history: Received 15 September 2014 Received in revised form 10 March 2015 Accepted 11 March 2015 Available online 12 March 2015
A two-step model of the kinetics of detonation combustion of methane in mixtures with an oxygen and air has been developed. The proposed model of the kinetics of detonation combustion of methane appears to be reasonably accurate within the confines of the assumptions and is consistent with the second law of thermodynamics. Constants of the model have a clear physical meaning. The model is useful for multi-dimensional numerical simulations of detonation processes. A numerical simulation of a twodimensional (2D) structure of the detonation wave (DW) in a stoichiometric methaneeair mixture in a wide range of channel height has been performed. The changes in the 2D structure of the selfsustaining DW on variations of the height of the channel have been studied. From the analysis of the flow structure, the dominant size of the detonation cell for stoichiometric methaneeair mixture determines to be 34 ± 1 cm, and 0.3÷0.35 cm for methaneeoxygen mixture. These values are in good agreement with all known experimental data. The transverse size of the detonation cell is a characteristic size of the DW front structure. Based on this value, it is possible to determine such parameters of detonation as critical conditions of detonation combustion, the critical energy of direct initiation, detonation propagation limits, etc., i.e. to estimate the detonation hazard of gaseous mixture and its use in various schemes of detonation propulsion systems. The 2D simulations were reproduced the DW irregular cellular structure with all its main features observed in the experiment: a chaotic uncoordinated movement of the main transverse waves; the presence of a fine (secondary) cell structure at the transverse waves themselves; the existence of numerous secondary transverse waves at the leading shock front, forming a hierarchy of the decreasing the size DW front perturbations; and a significant number of pockets of the unburned mixture at a considerable distance behind the DW front, etc. Numerical simulation of irregular structure of the multi-front DW has been conducted for a real hydrocarboneair mixture with an adequate description of all its thermal and chemical properties. © 2015 Elsevier Ltd. All rights reserved.
Keywords: Gaseous detonation Transverse wave Detonation cell Critical tube diameter Critical initiation energy Detonation velocity Methane Chemical kinetics Irregular multifront structure
1. Introduction Extensive experimental and theoretical studies of gas detonations (Voitsekhovskii et al., 1966; Fickett and Davis, 1979) have shown that the actual structure of the detonation front is very different from that predicted by the one-dimensional Zeldovich theory (Zeldovich et al., 1955). The flow behind a plane detonation wave front is unstable. Due to the growth of the initial disturbances,
Abbreviations: DW, Detonation wave; TW, Transverse wave; CJ, ChapmaneJouguet; TVD, Total variation diminishing; MUSCL, Monotonic UpstreamCentered Scheme for Conservation Laws; ASIRK, Additive semi-implicit RungeKutta; HLLC, Harten, Lax, van Leer, Contact (surface); CPU, Central Processing Unit. * Corresponding author. E-mail address: [email protected] (A.V. Trotsyuk). http://dx.doi.org/10.1016/j.jlp.2015.03.012 0950-4230/© 2015 Elsevier Ltd. All rights reserved.
the leading shock front is no longer smooth and exhibits a complex time-periodic structure. The main element of this structure is a triple configuration consisting of a Mach stem (overdriven DW), an incident wave (decaying DW), and a reflected (transverse) wave adjacent to the first two waves at the triple point (Voitsekhovskii et al., 1966; Fickett et al., 1979). Collisions of transverse waves moving over the DW front in opposing directions lead to the reproduction of the front structure in time. The trajectories of triple points are two intersecting families of lines that form a network of diamond cells. These detonation structures are fixed in the experiment on the soot imprints (Voitsekhovskii et al., 1966; Fickett et al., 1979; Nettleton, 1987). The transverse size of the cell, a0, determined as the average distance between the trajectories of one family measured in the direction normal to DW propagation, is a characteristic size of the DW front
A.V. Trotsyuk et al. / Journal of Loss Prevention in the Process Industries 36 (2015) 394e403
structure. One of the most important problems in the theory of dynamic systems is the formation and destruction of ordered gas-dynamic structures in a reactive medium. An example of such a selforganizing system is a multifront (cellular) DW steadily propagating in a reacting gas mixture. The complex three-dimensional and time-dependent structure of the front of these waves propagating in constant cross-section channels shows some order, whose geometrical parameter is the size of the elementary cell of the DW a0. Based on this value, it is possible to determine such parameters of detonation as critical conditions of detonation combustion, the critical energy of initiation (Vasil'ev and Grigor'ev, 1980), detonation propagation limits, etc., i.e. to estimate the detonation hazard of gaseous mixture. Based on the currently available experimental data, gas mixtures can be divided into several groups (Voitsekhovskii et al., 1966; Fickett et al., 1979; Nettleton, 1987) according to the degree of regularity of cells: 1) mixtures with a very regular detonation-cell structure (all cells at soot imprints are strikingly similar in size and shape); 2) mixtures with a slightly irregular structure (cells are slightly different in size and shape); 3) slightly regular mixtures with an appreciable variation in cell size; 4) mixtures with a highly irregular structure (cells vary greatly in both shape and size, and the characteristic size of the structure can be determined quite subjectively). The first group includes, for example, hydrogeneoxygen mixtures highly diluted with argon. The last group includes hydrocarboneair mixtures. A particularly high degree of irregularity is observed fir methaneeair mixtures, for which it is very difficult to identify a characteristic ordered spatial structure and the size of the detonation cell a0 can only be interpreted as a quantity averaged over the soot imprints in different parts of the detonation tube and over many experiments. In our previous numerical simulation (Trotsyuk, 1999) of the two-dimensional cellular DW structure in hydrogeneoxygen mixtures, good agreement was obtained between the numerical results and experimental data on the size a0 over a wide range of initial pressures and degrees of dilution with argon. A two-scale (bifurcation) cellular detonation structure was first found in numerical studies of the DW structure in hydrogeneoxygen mixtures with the addition of hydrogen peroxide (Fomin et al., 2006; Vasil'ev et al., 2010). Note that in the vast majority of studies, regular cellular structures in model gas mixtures have been simulateded; see the literature review in papers by Fomin et al. (2006) and Vasil'ev et al. (2010). Real mixtures are studied using detailed kinetic mechanisms of gas reactions. The use of this method in multi-dimensional numerical simulations is hampered by its cumbersomeness, programming difficulties, and high demands on computing power. Therefore, such studies are relatively few, and they also model DWs with regular or slightly irregular cell structures. As regards studies of irregular detonation structures, an irregular cellular structure was obtained in calculations of Gamezo et al. (1999, 2000) for a model mixture with thermochemical parameters very different from those of a real gas mixture resulted. Use was made of a simple chemical kinetics in the form of a single Arrhenius-type equation which described an irreversible exothermic reaction with high activation energy. Two-dimensional irregular structure of detonation in methaneeair mixture is numerically obtained firstly by Kessler et al. (2010, 2011), using a one-step irreversible Arrhenius equation also. The objectives of this work are to develop an approximate model for the kinetics of chemical reactions to describe the detonation combustion of methane based gaseous mixtures and to perform a numerical simulation of the irregular multifront (cellular) structure of a two-dimensional detonation wave in a
395
stoichiometric methaneeair mixture under the standard initial conditions. 1.1. Governing equations and model of chemical kinetics The dynamics of the compressible chemically reactive medium is described by the two-dimensional Euler equations. The chemical reaction in the DW is described according to our model for the combustion of methane using the two-step model of the detonation kinetics (induction period and main heat release stage step) first proposed by Korobeinikov et al., 1972). The induction period for methane mixtures is calculated by the empirical formula of Soloukhin (1970) suggested in paper (Vasil'ev, 1998). As shown by Vasil'ev (1998), the induction kinetics of Soloukhin (1970) provides the best fit to the available experimental data on the size of the detonation cell a0 in methane based mixtures when using an analytical model of the detonation cell (Vasil'ev and Nikolaev, 1976). A two-step model for the kinetics of chemical reactions in hydrogeneair mixtures (NFZT model) was developed in Nikolaev (1978); Nikolaev and Fomin (1982); Nikolaev and Zak (1988); Fomin and Trotsyuk (1995) As in (Korobeinikov et al., 1972), the reaction is divided into two steps: the induction period and the heat-release step. The model includes one ordinary differential equation for the molar mass of the mixture after the induction period and explicit algebraic formulas for the internal energy and the adiabatic index as functions of the pressure p and temperature T. The model adequately describes chemical equilibrium in hydrogeneair and hydrocarboneair mixtures, allowing calculations of the molar mass, internal energy, and isentropic index of the mixture as functions of T and p using explicit algebraic formulas. The model is highly accurate and consistent with the second law of thermodynamics, and the constants of the model have a clear physical meaning. In this paper, the two-step kinetic model NFZT is modified to describe the chemical reaction in detonation waves in methaneeair mixtures: CH4 þ a1H2 þ a2O2 þ a3N2 þ a4Ar þ a5H2O, a2 > 0.5. We replace the real multistep processes occurring during the induction period by some bruttoereaction chosen from the following considerations. At the end of the induction period, the temperature increase (and, hence, the total heat release of the chemical reactions) is small. Therefore, we choose the bruttoereaction so that its heat Qt is much lower than the maximum possible heat release Qmax from the total recombination of the reaction products with the formation of CO2 and H2O molecules. Considering that the induction period involves chemical reactions associated with the formation and increase in number of active centers and the disintegration of methane molecules, we assume that during the induction period, all methane molecules decomposed to form CO molecules. In accordance with the foregoing, during the induction period, all methane molecules will sooner or later undergo the following chemical reactions: CH4 þ O2 / CO þ OH þ 1.5H2, if a2 1,
(1)
CH4 þ a2O2 / CO þ (2a2 1)OH þ (2.5 a2)H2, if 1 > a2 0.5.(2) If the initial mixture contains molecules of H2O, H2, N2, and Ar, we assume that they are not involved in chemical reactions during the induction period.
396
A.V. Trotsyuk et al. / Journal of Loss Prevention in the Process Industries 36 (2015) 394e403
It is easy to show that for reaction (1) Qt << Qmax, and for reaction (2), Qt0.1Qmax for a3 0.8. Thus, the assumption of the model that the heat release during the induction period is small is valid for a3 0.8. Let b be the fraction of methane molecules that have not undergone chemical reactions by a certain time t. The quantity b decreases monotonically as the mixture moves in the induction zone: b ¼ 1 at t ¼ 0 and b ¼ 0 at the end of the induction period t*. If the quantity b is known, then, according to the foregoing, the chemical composition of the gas during the induction period is also known, and its thermodynamic parameters can easily be calculated using the corresponding tabulated polynomials or the approximate formulas of the NFZT model. Note the quantity b affects only the wave parameter profiles in the induction zone. The wave velocity, the flow parameters at the ChapmaneJouguet (CJ) point, and the main heat release zone do not depend on b. Therefore, once the above conditions on the quantity b are satisfied, the specific form of the equation for calculating this quantity affects the wave parameter profiles in the induction zone only quantitatively and has no fundamental significance. As a rule, the chemical reaction rate increases at the end of the induction period. In this connection, the formula for calculating b should be chosen so as to satisfy the boundary conditions and so that the rate of decrease of b increases as the mixture moves in the induction zone. For example, we can use the formula
b ¼ 1 Y n ; n 1: Here Y is the fraction of the induction period:
Zt YðtÞ ¼ 0
dt ; ti
ti is the induction period at constant parameters of the gas. At t ¼ 0, we have Y ¼ 0; at t ¼ t*, Y ¼ 1. The value of ti is calculated by an experimental algebraic formula by Soloukhin (1970): E ti ¼ Aind ,exp ind ,½O2 1 eff ; RT where ½O2 eff ½mol=cm3 is the oxygen concentration, Aind ¼ 6.0$1015 mol s/cm3, Eind ¼ 33,200 cal/mol. In the main heat release zone, the molar mass m of the gas is calculated by integrating the differential equation (Nikolaev and Zak, 1988):
2 dm r2 m ¼ 4Kþ 1 dt mmax m 3=2 m AT 3=4 1 eq=T r 1 eE=RT ; mmin
q, mmin, mmax, Kþ, and A are constants of the model NFZT: E ¼ 4.7235,1012 erg/mol, q ¼ 3000 K, Kþ ¼ 0.6,1015 cm6/(mol2,s), A ¼ 6.045,1014 cm3/(mol,s,K3/4), mmin ¼ 18.817 g/mol, mmax ¼ 27.739 g/mol. The internal energy and the adiabatic index of the mixture are calculated by algebraic equations (Nikolaev and Zak, 1988). Table 1 shows the results of calculation of CJ parameters and Von Neumann spike parameters obtained by the proposed kinetics for ideal plane detonation wave. The corresponding results of the calculations, presented by Nikolaev and Topchiyan (1977) are shown. In this paper detailed equation of chemical equilibrium are used for calculation chemical composition and thermodynamic parameters of the gas in detonation wave. The induction period ti and induction zone length xind in Table 1 were estimated using formula by Soloukhin (1970). Similarly to the NFZT model, the proposed kinetic model for the detonation combustion of methane is consistent with the second law of thermodynamics and its constants have a clear physical meaning and are calculated from the tabulated thermochemical parameters of the mixture before the two-dimensional numerical calculation of the DW structure. Determinations of the CJ wave parameters and the DW cell size showed high accuracy of the model. An advantage of this kinetics is its simplicity and ease of use in multi-dimensional gas dynamics codes. This kinetic model is described in more detail in our papers (Fomin et al., 2013, 2014). Numerical calculation of Zel'dovich-von Neumann-Doering structure of one-dimensional stationary DW in methaneeair mixture based on this model of chemical kinetics has been performed by (Fedorov et al., 2014). The system of governing equations was closed by the wellknown thermal equation of state for an ideal gas. Numerical simulation was performed for a stoichiometric methaneeair mixture at normal initial conditions. DW propagation in a two-dimensional straightline channel was studied. The slip conditions (impermeable solid wall) were imposed on the upper, lower, and left boundaries of the channel according to Godynov (1976), and the right boundary was subjected to the condition of an undisturbed initial state of the gas. The detonation wave was initiated near the left closed end of the channel and propagated from left to right along the x axis. The size and energy of the initiation source (Trotsyuk, 1999) were chosen sufficiently large to ensure a supercritical regime of DW direct initiation in methaneeair mixture. The problem of determining the minimum critical initiation energy was not considered in this study. 2. Numerical method
where r is the density of the mixture, R is universal gas constant, E,
The resultant systems of equations were solved numerically. The following discretization in space and time was performed. A uniform motionless grid with the total number of cells Ny and cell size hy ¼ H/Ny was used in the y direction, where H is a height of
Table 1 Velocity and parameters of plane ideal DW in stoichiometric CH4eair mixture at T0 ¼ 298.15 K, p0 ¼ 1 atm. ChapmaneJouguet point parameters
The present approximate model Nikolaev and Topchiyan (1977)
D, m/s
TCJ, K
pCJ, atm
rCJ, g/cm3
uCJ, m/s
cCJ,e, m/s
mCJ, g/mol
1809 1801
2801 2780
17.37 17.19
2.051$103 2.048$103
809.0 803.8
999.7 997.2
27.13 27.18
Von Neumann spike parameters T10, K The present approximate model Nikolaev and Topchiyan (1977)
1553 1531
p10, atm 31.57 31.33
r10, g/cm3 3
6.874$10 6.913$103
u10, m/s
ti, ms
xind, cm
1510 1506
6.02 6.98
0.1796 0.2058
A.V. Trotsyuk et al. / Journal of Loss Prevention in the Process Industries 36 (2015) 394e403
channel. A moving grid with the total number of cells Nx was used in the x direction, Nx1 of these cells formed a uniform grid with the cell size hx. Almost in all calculations hx ¼ hy. The uniform grid moved together with the leading shock front of the DW in such a way that it covered flow regions with large gradients of parameters and/or shock waves, and thus the fine DW structure was captured. The remaining Nx2 ¼ NxNx1 cells formed a nonuniform moving grid that occupied the region between the closed left end of the channel and the beginning of the uniform grid. The use of adaptive moving grids with local refinement allowed us to have a fine resolution where necessary, using a significantly smaller number of grid cells, as compared with a uniform splitting of the computational domain. A typical of numbers of grid cells for high resolution computations were Nx ¼ 2000, Nx1 ¼ 1875, Ny ¼ 4000, and 256 coprocessors were used. The finite-volume Godunov-type scheme (Godynov, 1976) has been used with the MUSCL TVD reconstruction procedures. The fourth-order reconstruction by Yamamoto and Daiguji (1993), Daiguji et al. (1997) are used to calculate the fluxes along every axes in a case of the uniform grid zone, and a third-order scheme Chakravarthy and Osher (1985), Lin and Chin (1995) is used in the non-uniform grid zone along x axes. The advanced HLLC solver by Batten et al. (1997) is used for an approximate solution of the Riemann problem on the control volume boundaries. To implement this algorithm for obtaining numerical fluxes for the case of a chemically reacting mixture, the ‘‘energy relaxation method’’ by Coquel and Perthame (1998) is applied. This method eliminates the problems of numerical solution of the Riemann problem for a medium with a complicated nonlinear equation of state (including that with a variable ratio of specific heats for gaseous mixtures or any form of equation of state for condensed materials). A special work has been done to fit these algorithms to simulations of reacting detonations flows. The system of conservation equations for simulating the detonation process in a reacting medium is a hyperbolic system of Euler equations with stiff relaxation term in the right-hand sides. The stiffness of this system during its integration in time is caused by the source terms, which model chemical kinetics in detonation wave. In our numerical code we use the ASIRK method (Shen and Zhong, 1996) of the second order accuracy for integration on time. The semi-implicit method means that the stiff relaxation terms in this method are treated by an implicit manner, and the non-stiff convection terms by an explicit calculation. The time step is determined at each time layer of the solution from the stability condition, see (Godynov, 1976). The typical values of the Courant number in detonation problems simulations are CFL ¼ 0.3e0.4. The code is parallelized with MPI library using the domain decomposition technique. Two approaches were used. For a comparatively small number of CPUs (up to 20), the entire computational domain is decomposed into equal subdomains along one of the axes according to the number of CPUs used in the computation. The data exchange is performed only between the adjacent subdomains. The amount of the exchanged information depends on the type of used scheme. For our case it is three halo cells (three layers of boundary cells) from each side of virtual subdomain's boundary. But an amount of exchanged information also depends on strategy of decomposition of numerical computation domain to subdomains. The main idea is to minimize the volume of information transferred between subdomains for accelerating computations. For example, in the case of 2D simulations for a very big number of CPUs the decomposition only along one of the axes would result in the subdomains in form of very thin strips with a great number of boundary cells. It could lead to a prohibitive increase of volume of exchanged information. The optimal strategy is
397
to divide the entire computational domain into equal subdomains along the axes x and y according to the number of CPUs used in computation. The decomposition is performed in such a way, that subdomains have a square form in computational space Nx Ny. Although simple, this parallelization technique yields good results in terms of parallel efficiency and scalability.
Fig. 1. Detonation wave structure in a stoichiometric methaneeair mixture in channel with H ¼ 35 cm at xfront ¼ 1500 cm, pair of primary transverse wave AA and BB: (a) numerical Schlieren-visualization of all flowfield; (b) magnified fragment of the flowfield near an upper wall of a channel.
398
A.V. Trotsyuk et al. / Journal of Loss Prevention in the Process Industries 36 (2015) 394e403
Fig. 3. Detonation wave structure in stoichiometric methaneeair mixture in channel with H ¼ 35 cm at xfront ¼ 1500 cm “Low” resolution: hx ¼ hy z 160 mm. Numerical Schlieren visualization of all flowfield. AA and BB e primary transverse waves.
Fig. 2. Detonation wave structure in a stoichiometric methaneeair mixture in channel with H ¼ 35 cm at xfront ¼ 1500 cm: (a) normalized density flowfield; (b) temperature (K) flowfield.
3. Results of calculations The results of the simulation of the stoichiometric methaneeair mixture are presented in Figs. 1e3, 5e9, which show the fields of different parameters in the DW, including a numerical Schlieren visualization of the flowfield. The use of this parameter to visualize the structure of the detonation wave is due to its advantages over the other parameter, for example, pressure. The density gradient field accurately represents the flow pattern. Density changes
Fig. 4. Front-propagation velocity (solid line) as a function of front position in channel with H ¼ 35 cm, average velocity Daver ¼ 1840 m/s (dashed line), and ideal CJ velocity DCJ ¼ 1809 m/s (dash-dot line).
visualize not only the shape of the leading shock front of the DW and the structure of transverse waves, but also contact discontinuities. Significant changes in density at the end of the induction period, which are associated with the beginning of the main heat release stage, clearly show the two-dimensional structure of the DW induction region. In all the figures, the x and y coordinates are given in centimeters. The detonation cell size for the methaneeair mixture was determined in a procedure similar to that used in the numerical study of regular and double-scale cellular structures, see (Trotsyuk,
A.V. Trotsyuk et al. / Journal of Loss Prevention in the Process Industries 36 (2015) 394e403
Fig. 5. Detonation wave structure in a channel with H ¼ 35 cm at xfront ¼ 650 cm, Daver ¼ 1669 m/s, DCJ ¼ 1809 m/s, one primary transverse wave AA: (a) numerical Schlieren-visualization flowfield; (b) normalized density flowfield.
1999; Fomin et al., 2006; Vasil'ev et al., 2010). The initial height of the channel H was chosen rather arbitrarily; usually, it was assumed to be slightly higher than the presumed size of the cell a0. This simulation was started with H ¼ 15 cm, which was found later to be considerably smaller than the size of the cell. Further variation of H in the calculations affected the number of primary transverse waves remaining at the DW front after its initiation and the establishment of a self-sustained CJ detonation regime, thus, the cell size was determined. Figures 1 and 2 show the DW structure in a channel of height
399
Fig. 6. Detonation wave structure in a channel with H ¼ 25 cm at xfront ¼ 1140 cm, one primary transverse wave AA: (a) normalized density flowfield; (b) temperature (K) flowfield.
H ¼ 35 with the wave front located at xfr ¼ 1500 cm. Analysis of the change in the average velocity of the DW as it propagates away from the initiation source shows that a steady-state CJ detonation had already been established by this time. The figure shows a system of the two main (primary) transverse waves (TWs), AA and BB, which are located approximately at y z 12 cm and y z 19 cm. At this time, transverse waves move in the opposite direction to each other, in antiphase, as in mixtures with regular cellular
400
A.V. Trotsyuk et al. / Journal of Loss Prevention in the Process Industries 36 (2015) 394e403
Fig. 7. Detonation wave structure in a channel with H ¼ 40 cm at xfront ¼ 1170 cm, three primary transverse wave AA, BB and CC: (a) normalized density flowfield; (b) temperature (K) flowfield.
structure (compare with the figures of our papers (Trotsyuk, 1999; Vasil'ev & Trotsyuk, 2003). Moreover, the DW front has a system of secondary (smaller) transverse waves of different size, intensity, and direction of motion. One of these secondary TW at y z 30 cm is denoted as ‘aa’. The movement of these perturbations in the leading shock front of the DW is quite random. To identify the primary, most intense and longest, transverse waves (forming the
Fig. 8. Detonation wave structure in a channel with H ¼ 33 cm at xfront ¼ 600 cm, one primary transverse wave AA: (a) normalized density flowfield; (b) temperature (K) flowfield.
detonation cell) among these perturbations, it is necessary to analyze the flowfield of other parameters, in particular, the density and temperature, see Fig. 2. In all figures, the main transverse waves are denoted as AA, BB, etc. Fig. 1b presents an enlarged fragment of the DW front in the region y z 26 cm, which shows more details of the numerous secondary transverse waves forming the hierarchy of perturbations on the leading DW front decreasing in size. The secondary transverse waves that form the fine structure of the leading DW front,
A.V. Trotsyuk et al. / Journal of Loss Prevention in the Process Industries 36 (2015) 394e403
Fig. 9. Detonation wave structure in a channel with H ¼ 33 cm at xfront ¼ 1500 cm, pair of primary transverse wave AA and BB: (a) normalized density flowfield; (b) temperature (K) flowfield.
had been observed in the experiments of Manzhalei (1977). Fig. 1b shows the fine cellular (multifront) structure of the front one of the TW ‘aa’, located at y z 30 cm. The TWs on the front of wave ‘aa’ are denoted as 1,2,3,4. Thus, this simulation reproduces a characteristic feature of the detonation of hydrocarbon mixtures e instability of the front of the transverse waves itself at the leading shock front of the DW. This feature has been observed in soot imprints of typical trajectories of TWs in numerous experiments of Manzhalei & Mitrofanov (1973). In general, the shock-wave structure of the
401
DW front is extremely irregular (compare with the figures from Trotsyuk (1999); Vasil'ev & Trotsyuk (2003)). Furthermore, in Fig. 1b, one can see a significant amount of unburned pieces of gas at a considerable distance behind the leading shock front of the DW. The formation of these packages is in agreement with experimental observations (Subbotin, 1975; Kiyanda and Higgins, 2013). This phenomenon was not observed in simulations of regular mixtures based on hydrogen (Trotsyuk, 1999; Vasil'ev & Trotsyuk, 2003). The presence of these small-scale transverse waves on DW front and the unburned pockets of gas in the mixture under study leads to the need to maintain an extremely high spatial resolution of the computational grid in the induction zone and in the reaction zone behind the leading shock front. Most calculations were carried out with a characteristic mesh size hy ¼ hx z 80 mm, see Figs. 1, 2, 4e9. To study grid convergence of the numerical solutions the calculations have been performed with different sizes of the computational mesh. Numerical simulations with twice greater grid size hy ¼ hx z 160 mm in channel with H ¼ 35 cm show the same number of primary transverse waves, see Fig. 3. Of course, the fine structures of the flow were more blurred and smoothed as compared with fine grid simulations. So, the grid convergence of the numerical solutions was obtained. High spatial resolution led to a small allowed step in time and a large number of computational meshes. This, in turn, required parallelization of the numerical code and the use of a multi-processor supercomputer for numerical simulation of the detonation of the methane mixture. Fig. 4 shows a front-propagation velocity (solid line) as a function of DW front position after the initiation in a channel H ¼ 35 cm, parameters of calculations correspond to Figs. 1 and 2. Denotation has reached a quasi-steady propagation regime at the last 5 m. The average DW velocity obtained on this base is Daver ¼ 1840 m/s (dashed line). This value is slightly higher than the ChapmaneJouguet velocity (dash-dot line) DCJ ¼ 1809 m/s derived from thermodynamic calculations. Approximately the same absolute difference in the velocity of the multifront DW in a self-sustained CJ regime with DCJ was observed for hydrogeneoxygen mixtures with regular cell structure (Trotsyuk, 1999). The calculations were performed up to distances of more than ten meters, traveled by DW, because methane mixtures have extremely high critical initiation energy (Vasil'ev, 1998; Vasil'ev et al., 2000). Therefore, in the calculations, we had to use a highenergy initiator, which led to unsteadiness of the detonation process up to large distances. As an example, Fig. 5 shows the DW structure in the same channel with the front located at xfr ¼ 650 cm. In contrast to the pattern in Figs. 1 and 2, here we see only one main TW at y z 5 cm, and strong secondary TW ‘aa’ at y z 26 cm. The system of secondary transverse waves of different size at the leading front and the presence of regions of unburned gas are similar to the flow in s 1,2. At this point, the DW velocity averaged over two-dimensional pulsations, Daver ¼ 1669 m/s is much lower than the CJ velocity DCJ ¼ 1809 m/s. This suggests that at this distance, a self-sustained CJ detonation has not yet been established. Calculations with varying H showed that the cellular structure of the DW is significantly different from the structures observed with the channel sizes in Figs. 1 and 2. Thus, at H ¼ 25 cm, detonation flow with a single transverse wave at the DW front is formed. The normalized density field and the temperature field for the flow are shown in Fig. 6. Here we see one intense and extended TW AA and a no less extended, but much less intense degenerate (weak) TW ‘aa’ near the upper wall of the channel. The induction zone behind the leading shock front has a very irregular shape, and there are many large regions of unburned mixture in the region y z 2.5 cm. If H ¼ 40 cm, we have at least three main transverse waves, AA, BB, and CC (see Fig. 7); otherwise the structure of the front is similar to that shown in Figs. 1 and 2.
402
A.V. Trotsyuk et al. / Journal of Loss Prevention in the Process Industries 36 (2015) 394e403
To investigate the influence the numerical grid on the results of simulation of the DW structure, we performed additional calculations for H ¼ 33 cm with an increased number of cells along the x axis: Nx ¼ 4000 and Nx1 ¼ 3500. The number of cells along the y axis remained unchanged, Ny ¼ 4000, and hy ¼ hx in the zone with a fine mesh in the vicinity of the leading shock front. Figs. 8 and 9 show the simulation results for a channel with H ¼ 33 cm. As can be seen from the flowfields, the boundary between fine uniform grid and the nonuniform grid along the x axis is much less pronounced in these calculations than in those shown in the previous figures. At xfr ¼ 600 cm in Fig. 8, we have structure of an unsteady DW in the transient state, which is largely similar to the DW structure in Fig. 5. The multifront DW structure of steady self-sustaining CJ detonation is shown at xfr ¼ 1500 cm in Fig. 9. As in the case of detonation in the channel with H ¼ 35 cm (see Figs. 1 and 2), at the leading shock front of the DW there are two main transverse waves, AA and BB, and numerous secondary transverse waves of different lengths and amplitudes. In addition, the flow pattern has features that distinguish it from that in Figs. 1 and 2. It is seen that in the channel of this size, the movement of the transverse waves is less consistent; at the time presented in the figure, the transverse waves move in the same direction toward the bottom wall of the channel. Obviously, after the reflection of the BB wave from the bottom wall, it will collide with the AA wave not exactly on the axis of the channel. This will result in a distorted shape of the detonation cell, which is typical of all hydrocarboneair mixtures. In the case of mixtures with regular cellular structure, it could be argued that with this behavior of transverse waves, the channel width is slightly smaller than the transverse dimension of the cell a0. However, for the type of cellular structure considered here, this conclusion would be incorrect. From an analysis of the numerical results, the dominant transverse dimension of the irregular detonation cell formed by the main (primary) transverse waves AA and BB for stoichiometric methaneeair mixture can be determined to be a0 ¼ 34 ± 1 cm This value is in good agreement with the experimental value of the cell for this mixture (see (Nettleton, 1987; Vasil'ev, 1998; Vasil'ev et al., 2000; Bull et al., 1982; Knystautas et al., 1984; Moen et al., 1984)) and lies in the middle of the range of the experimental data. A numerical simulation of a two-dimensional structure of the DW in a stoichiometric methaneeoxygen mixture has been carried out. The fine multifront structure of the primary TWs is not clearly revealed in the present simulations. Calculations have been performed with a mesh size hy ¼ hx z 1.5 mm. Perhaps this resolution was not enough to resolve the fine structure of the transverse wave, requires further numerical studies. The size of the detonation cell for stoichiometric methaneeoxygen mixture from our numerical computations determines to be a0 ¼ 0.3÷0.35 cm. This value also is in good agreement with all known experimental data (Nettleton, 1987; Vasil'ev, 1998; Vasil'ev et al., 2000). 4. Conclusions An approximate two-step model (induction parameter model) for the kinetics of detonation combustion of methane was developed. Because of the simplicity and high accuracy (within the confines of the assumptions), the model can be used in multidimensional numerical simulations of detonation flows of methaneeoxygen and methaneeair mixtures. From the analysis of two-dimensional flow structures, the dominant size of the detonation cell for stoichiometric methaneeair mixture determines to be 34 ± 1 cm, and 0.3÷0.35 cm for methaneeoxygen mixture. These data reveal good agreement with available experimental data. This simulation takes into account the real thermophysical and
chemical properties of the mixture, and reproduces an irregular cell structure with all of its main features: random inconsistent movement of the main transverse waves; numerous secondary transverse waves forming the hierarchy of perturbations of decreasing size at the leading shock front of the DW; regions of unburned mixture at a considerable distance behind the DW front; fine (cellular) structure of the transverse wave. These calculations are the first in which the irregular structure of detonation wave is obtained using a two-step kinetic model. In experiments the fine structure of TW is manifested in the typical shape of transversewave trajectories recorded in soot imprints. The secondary waves at the leading shock front are responsible for the small-cell network inside the main cell on soot imprints in experiments with methane mixtures. Acknowledgments The authors gratefully acknowledge the financial contribution from the Russian Foundation for Basic Research (Grant No. 14-0100406 and Grant No. 14-03-00838), the Program of Basic Research No. 26.2 entitled “Combustion and explosion” of the Presidium of the Russian Academy of Sciences. The computations were performed on an MVS-100K supercomputer at the Joint SuperComputer Center (JSCC) of the Russian Academy of Sciences (Moscow). References Batten, P., Leschziner, M.A., Goldberg, U.C., 1997. Average-state Jacobians and implicit methods for compressible viscous and turbulent flows. J. Comput. Phys. 137, 38e78. Bull, D.C., Elsworth, J.E., Shuff, P.J., 1982. Detonation cell structures in fuel/air mixtures. Combust. Flame 45, 7e22. Chakravarthy, S.R., Osher, S., 1985. A new class of high accuracy TVD schemes for hyperbolic conservation laws. In: Proceedings of the AIAA 23rd Aerospace Sciences Meeting. Reno, Nevada, pp. 85e0363. AIAA Paper. Coquel, F., Perthame, B., 1998. Relaxation of energy and approximate Riemann solvers for general pressure laws in fluid dynamics. SIAM J. Numer. Anal. 35, 2223e2249. Daiguji, H., Yuan, X., Yamamoto, S., 1997. Stabilization of higher-order high resolution schemes for the compressible Navier-Stokes equation. Int. J. Numer. Methods Heat Fluid Flow 7, 250e274. Fedorov, A.V., Fomin, P.A., Tropin, D.A., 2014. Simple kinetics and detonation wave structure in a methane-air mixture. Combust. Explos. Shock Wave 50, 87e96. Fickett, W., Davis, W.C., 1979. Detonation. University of California Press, Berkeley, CA. Fomin, P.A., Trotsyuk, A.V., 1995. An approximate calculation of the isentrope of a gas in chemical equilibrium. Combust. Explos. Shock Waves 31, 455e458. Fomin, P.A., Trotsyuk, A.V., Vasil’ev, A.A., Mitropetros, K., Hieronymus, H., Roekaerts, D., 2006. Model of chemical reaction kinetics for calculating detonation processes in gas and heterogeneous mixtures containing hydrogen peroxide. Combust. Sci. Technol. 178, 895e919. Fomin, P.A., Trotsyuk, A.V., Vasil’ev, A.A., 2013. Approximate models of chemical kinetics for detonation processes in mixtures of CH4, H2O2 and O3 with air. In: Proceedings of the 24th ICDERS. July 28eAugust 2, 2013, Taipei, Taiwan. Paper 0161. Fomin, P.A., Trotsyuk, A.V., Vasil’ev, A.A., 2014. Approximate model of chemical reaction kinetics for detonation processes in mixture of CH4 with air. Combust. Sci. Technol. 186, 1716e1735. Gamezo, V.N., Desbordes, D., Oran, E.S., 1999. Formation and evolution of twodimensional cellular detonations. Combust. Flame 116, 154e165. Gamezo, V.N., Vasil’ev, A.A., Khokhlov, A.M., Oran, E.S., 2000. Fine cellular structures produced by marginal detonations. In: Proceedings of the 28th Symposium (International) on Combustion. The Combustion Institute, Pittsburgh, pp. 611e617. Godynov, S.K. (Ed.), 1976. Numerical Solution of the Multi-dimensional Problems of Gaseous Dynamics. Nauka Press, Moscow (in Russian). Kessler, D.A., Gamezo, V.N., Oran, E.S., 2010. Simulations of flame acceleration and deflagration-to-detonation transitions in methane-air systems. Combust. Flame 157, 2063e2077. Kessler, D.A., Gamezo, V.N., Oran, E.S., 2011. Multilevel detonation cell structures in methane-air mixtures. Proc. Combust. Inst. 33, 2211e2218. Kiyanda, C.B., Higgins, A.J., 2013. Photographic investigation into the mechanism of combustion in irregular detonation waves. Shock Waves 23, 115e130. Knystautas, R., Guirao, C., Lee, J.H., Sulmistras, A., 1984. Measurements of cell size in hydrocarbon-air mixtures and prediction of critical tube diameter, critical initiation energy, and detonability limits. Prog. Astronaut. Aeronaut. 94, 23e37.
A.V. Trotsyuk et al. / Journal of Loss Prevention in the Process Industries 36 (2015) 394e403 Korobeinikov, V.P., Levin, V.A., Markov, V.V., Chernyi, G.G., 1972. Propagation of blast waves in a combustible gas. Astronaut. Acta 17, 529e537. Lin, S.-Y., Chin, Y.-S., 1995. Comparison of higher resolution Euler schemes for aeroacoustic computations. AIAA J. 33, 237e245. Manzhalei, V.I., Mitrofanov, V.V., 1973. The stability of detonation shock waves with a spinning configuration. Combust. Explos. Shock Waves 9, 614e620. Manzhalei, V.I., 1977. Fine structure of the leading front of a gas detonation. Combust. Explos. Shock Waves 13, 402e404. Moen, I.O., Funk, J.W., Ward, S.A., Rude, G.M., Thibault, P.A., 1984. Detonation length scales for fuel-air explosives. Prog. Astronaut. Aeronaut. 94, 55e79. Nettleton, M.A., 1987. Gaseous Detonations: Their Nature, Effects and Control. Chapman and Hall, London, New York, ISBN 0-412-27040-4. Nikolaev, YuA., 1978. Model of the kinetics of chemical reactions at high temperatures. Combust. Explos. Shock Waves 14, 468e471. Nikolaev, YuA., Fomin, P.A., 1982. Analysis of equilibrium flows of chemically reacting gases. Combust. Explos. Shock Waves 18, 53e58. Nikolaev, YuA., Topchiyan, M.E., 1977. Analysis of equilibrium flows in detonation waves in gases. Combust. Explos. Shock Waves 13, 327e338. Nikolaev, YuA., Zak, D.V., 1988. Agreement of models of chemical reactions in gases with the second law of thermodynamics. Combust. Explos. Shock Waves 24, 461e464. Shen, J.W., Zhong, X., 1996. Semi-implicit Runge-Kutta schemes for nonautonomous differential equations in reactive flow computations. In: Proceedings of the AIAA 27th Fluid Dynamics Conference. New Orleans, LA, June 17-20, 1996. AIAA Paper 961969. Soloukhin, R.I., 1970. Measurement methods and basic results of shock-tube experiments. In: Proceedings of the 7th Int. Shock Tubes Symposium. 23e25 June
403
1969, Toronto, Canada. University of Toronto Press, Toronto and Buffalo, ISBN 08020-1729-0, pp. 662e706. Subbotin, V.A., 1975. Collision of transverse detonation waves in gases. Combust. Explos. Shock Waves 11, 411e414. Trotsyuk, A.V., 1999. Numerical simulation of the structure of two-dimensional gaseous detonation of an H2eO2eAr mixture. Combust. Explos. Shock Waves 35, 549e558. Vasil'ev, A.A., 1998. Effect of nitrogen on multifront detonation parameters. Combust. Explos. Shock Waves 34, 72e76. Vasil'ev, A.A., Nikolaev, YuA., 1976. Model of a multifront gas detonation cell. Combust. Explos. Shock Waves 12, 744e754. Vasil’ev, A.A., Grigor’ev, V.V., 1980. Critical conditions for gas detonation in sharply expanding channels. Combust. Explos. Shock Waves 16, 579e585. Vasil’ev, A.A., Valishev, A.I., Vasil’ev, V.A., Panfilova, L.V., Topchian, M.E., 2000. Detonation hazards of methane mixtures. Arch. Combust. 20, 31e48. Vasil'ev, A.A., Vasil'ev, V.A., Trotsyuk, A.V., 2010. Bifurcation structures in gas detonation. Combust. Explos. Shock Waves 46, 196e206. Vasil'ev, A.A., Trotsyuk, A.V., 2003. Experimental investigation and numerical simulation of an expanding multifront detonation wave. Combust. Explos. Shock Waves 39, 80e90. Voitsekhovskii, B.V., Mitrofanov, V.V., Topchian, M.E., 1966. The Structure of a Detonation Front in Gases. WrightPatterson (AFB AD633821). Yamamoto, S., Daiguji, H., 1993. Higher-order-accurate upwind schemes for solving the compressible Euler and Navier-Stokes equations. Comput. Fluids 22, 259e270. Zel'dovich, YaB., Kompaneets, A.S., 1955. The Theory of Detonation. Gostekhteorizdat, Moscow (in Russian).
Contents lists available at ScienceDirect
Journal of Loss Prevention in the Process Industries journal homepage: www.elsevier.com/locate/jlp
Numerical study of cellular detonation structures of methane mixtures Anatoliy V. Trotsyuk*, Pavel A. Fomin, Anatoly A. Vasil'ev Lavrentyev Institute of Hydrodynamics Siberian Branch of the Russian Academy of Science (LIH SB RAS), 15 Ac. Lavrentyev Ave., Novosibirsk 630090, Russia
a r t i c l e i n f o
a b s t r a c t
Article history: Received 15 September 2014 Received in revised form 10 March 2015 Accepted 11 March 2015 Available online 12 March 2015
A two-step model of the kinetics of detonation combustion of methane in mixtures with an oxygen and air has been developed. The proposed model of the kinetics of detonation combustion of methane appears to be reasonably accurate within the confines of the assumptions and is consistent with the second law of thermodynamics. Constants of the model have a clear physical meaning. The model is useful for multi-dimensional numerical simulations of detonation processes. A numerical simulation of a twodimensional (2D) structure of the detonation wave (DW) in a stoichiometric methaneeair mixture in a wide range of channel height has been performed. The changes in the 2D structure of the selfsustaining DW on variations of the height of the channel have been studied. From the analysis of the flow structure, the dominant size of the detonation cell for stoichiometric methaneeair mixture determines to be 34 ± 1 cm, and 0.3÷0.35 cm for methaneeoxygen mixture. These values are in good agreement with all known experimental data. The transverse size of the detonation cell is a characteristic size of the DW front structure. Based on this value, it is possible to determine such parameters of detonation as critical conditions of detonation combustion, the critical energy of direct initiation, detonation propagation limits, etc., i.e. to estimate the detonation hazard of gaseous mixture and its use in various schemes of detonation propulsion systems. The 2D simulations were reproduced the DW irregular cellular structure with all its main features observed in the experiment: a chaotic uncoordinated movement of the main transverse waves; the presence of a fine (secondary) cell structure at the transverse waves themselves; the existence of numerous secondary transverse waves at the leading shock front, forming a hierarchy of the decreasing the size DW front perturbations; and a significant number of pockets of the unburned mixture at a considerable distance behind the DW front, etc. Numerical simulation of irregular structure of the multi-front DW has been conducted for a real hydrocarboneair mixture with an adequate description of all its thermal and chemical properties. © 2015 Elsevier Ltd. All rights reserved.
Keywords: Gaseous detonation Transverse wave Detonation cell Critical tube diameter Critical initiation energy Detonation velocity Methane Chemical kinetics Irregular multifront structure
1. Introduction Extensive experimental and theoretical studies of gas detonations (Voitsekhovskii et al., 1966; Fickett and Davis, 1979) have shown that the actual structure of the detonation front is very different from that predicted by the one-dimensional Zeldovich theory (Zeldovich et al., 1955). The flow behind a plane detonation wave front is unstable. Due to the growth of the initial disturbances,
Abbreviations: DW, Detonation wave; TW, Transverse wave; CJ, ChapmaneJouguet; TVD, Total variation diminishing; MUSCL, Monotonic UpstreamCentered Scheme for Conservation Laws; ASIRK, Additive semi-implicit RungeKutta; HLLC, Harten, Lax, van Leer, Contact (surface); CPU, Central Processing Unit. * Corresponding author. E-mail address: [email protected] (A.V. Trotsyuk). http://dx.doi.org/10.1016/j.jlp.2015.03.012 0950-4230/© 2015 Elsevier Ltd. All rights reserved.
the leading shock front is no longer smooth and exhibits a complex time-periodic structure. The main element of this structure is a triple configuration consisting of a Mach stem (overdriven DW), an incident wave (decaying DW), and a reflected (transverse) wave adjacent to the first two waves at the triple point (Voitsekhovskii et al., 1966; Fickett et al., 1979). Collisions of transverse waves moving over the DW front in opposing directions lead to the reproduction of the front structure in time. The trajectories of triple points are two intersecting families of lines that form a network of diamond cells. These detonation structures are fixed in the experiment on the soot imprints (Voitsekhovskii et al., 1966; Fickett et al., 1979; Nettleton, 1987). The transverse size of the cell, a0, determined as the average distance between the trajectories of one family measured in the direction normal to DW propagation, is a characteristic size of the DW front
A.V. Trotsyuk et al. / Journal of Loss Prevention in the Process Industries 36 (2015) 394e403
structure. One of the most important problems in the theory of dynamic systems is the formation and destruction of ordered gas-dynamic structures in a reactive medium. An example of such a selforganizing system is a multifront (cellular) DW steadily propagating in a reacting gas mixture. The complex three-dimensional and time-dependent structure of the front of these waves propagating in constant cross-section channels shows some order, whose geometrical parameter is the size of the elementary cell of the DW a0. Based on this value, it is possible to determine such parameters of detonation as critical conditions of detonation combustion, the critical energy of initiation (Vasil'ev and Grigor'ev, 1980), detonation propagation limits, etc., i.e. to estimate the detonation hazard of gaseous mixture. Based on the currently available experimental data, gas mixtures can be divided into several groups (Voitsekhovskii et al., 1966; Fickett et al., 1979; Nettleton, 1987) according to the degree of regularity of cells: 1) mixtures with a very regular detonation-cell structure (all cells at soot imprints are strikingly similar in size and shape); 2) mixtures with a slightly irregular structure (cells are slightly different in size and shape); 3) slightly regular mixtures with an appreciable variation in cell size; 4) mixtures with a highly irregular structure (cells vary greatly in both shape and size, and the characteristic size of the structure can be determined quite subjectively). The first group includes, for example, hydrogeneoxygen mixtures highly diluted with argon. The last group includes hydrocarboneair mixtures. A particularly high degree of irregularity is observed fir methaneeair mixtures, for which it is very difficult to identify a characteristic ordered spatial structure and the size of the detonation cell a0 can only be interpreted as a quantity averaged over the soot imprints in different parts of the detonation tube and over many experiments. In our previous numerical simulation (Trotsyuk, 1999) of the two-dimensional cellular DW structure in hydrogeneoxygen mixtures, good agreement was obtained between the numerical results and experimental data on the size a0 over a wide range of initial pressures and degrees of dilution with argon. A two-scale (bifurcation) cellular detonation structure was first found in numerical studies of the DW structure in hydrogeneoxygen mixtures with the addition of hydrogen peroxide (Fomin et al., 2006; Vasil'ev et al., 2010). Note that in the vast majority of studies, regular cellular structures in model gas mixtures have been simulateded; see the literature review in papers by Fomin et al. (2006) and Vasil'ev et al. (2010). Real mixtures are studied using detailed kinetic mechanisms of gas reactions. The use of this method in multi-dimensional numerical simulations is hampered by its cumbersomeness, programming difficulties, and high demands on computing power. Therefore, such studies are relatively few, and they also model DWs with regular or slightly irregular cell structures. As regards studies of irregular detonation structures, an irregular cellular structure was obtained in calculations of Gamezo et al. (1999, 2000) for a model mixture with thermochemical parameters very different from those of a real gas mixture resulted. Use was made of a simple chemical kinetics in the form of a single Arrhenius-type equation which described an irreversible exothermic reaction with high activation energy. Two-dimensional irregular structure of detonation in methaneeair mixture is numerically obtained firstly by Kessler et al. (2010, 2011), using a one-step irreversible Arrhenius equation also. The objectives of this work are to develop an approximate model for the kinetics of chemical reactions to describe the detonation combustion of methane based gaseous mixtures and to perform a numerical simulation of the irregular multifront (cellular) structure of a two-dimensional detonation wave in a
395
stoichiometric methaneeair mixture under the standard initial conditions. 1.1. Governing equations and model of chemical kinetics The dynamics of the compressible chemically reactive medium is described by the two-dimensional Euler equations. The chemical reaction in the DW is described according to our model for the combustion of methane using the two-step model of the detonation kinetics (induction period and main heat release stage step) first proposed by Korobeinikov et al., 1972). The induction period for methane mixtures is calculated by the empirical formula of Soloukhin (1970) suggested in paper (Vasil'ev, 1998). As shown by Vasil'ev (1998), the induction kinetics of Soloukhin (1970) provides the best fit to the available experimental data on the size of the detonation cell a0 in methane based mixtures when using an analytical model of the detonation cell (Vasil'ev and Nikolaev, 1976). A two-step model for the kinetics of chemical reactions in hydrogeneair mixtures (NFZT model) was developed in Nikolaev (1978); Nikolaev and Fomin (1982); Nikolaev and Zak (1988); Fomin and Trotsyuk (1995) As in (Korobeinikov et al., 1972), the reaction is divided into two steps: the induction period and the heat-release step. The model includes one ordinary differential equation for the molar mass of the mixture after the induction period and explicit algebraic formulas for the internal energy and the adiabatic index as functions of the pressure p and temperature T. The model adequately describes chemical equilibrium in hydrogeneair and hydrocarboneair mixtures, allowing calculations of the molar mass, internal energy, and isentropic index of the mixture as functions of T and p using explicit algebraic formulas. The model is highly accurate and consistent with the second law of thermodynamics, and the constants of the model have a clear physical meaning. In this paper, the two-step kinetic model NFZT is modified to describe the chemical reaction in detonation waves in methaneeair mixtures: CH4 þ a1H2 þ a2O2 þ a3N2 þ a4Ar þ a5H2O, a2 > 0.5. We replace the real multistep processes occurring during the induction period by some bruttoereaction chosen from the following considerations. At the end of the induction period, the temperature increase (and, hence, the total heat release of the chemical reactions) is small. Therefore, we choose the bruttoereaction so that its heat Qt is much lower than the maximum possible heat release Qmax from the total recombination of the reaction products with the formation of CO2 and H2O molecules. Considering that the induction period involves chemical reactions associated with the formation and increase in number of active centers and the disintegration of methane molecules, we assume that during the induction period, all methane molecules decomposed to form CO molecules. In accordance with the foregoing, during the induction period, all methane molecules will sooner or later undergo the following chemical reactions: CH4 þ O2 / CO þ OH þ 1.5H2, if a2 1,
(1)
CH4 þ a2O2 / CO þ (2a2 1)OH þ (2.5 a2)H2, if 1 > a2 0.5.(2) If the initial mixture contains molecules of H2O, H2, N2, and Ar, we assume that they are not involved in chemical reactions during the induction period.
396
A.V. Trotsyuk et al. / Journal of Loss Prevention in the Process Industries 36 (2015) 394e403
It is easy to show that for reaction (1) Qt << Qmax, and for reaction (2), Qt0.1Qmax for a3 0.8. Thus, the assumption of the model that the heat release during the induction period is small is valid for a3 0.8. Let b be the fraction of methane molecules that have not undergone chemical reactions by a certain time t. The quantity b decreases monotonically as the mixture moves in the induction zone: b ¼ 1 at t ¼ 0 and b ¼ 0 at the end of the induction period t*. If the quantity b is known, then, according to the foregoing, the chemical composition of the gas during the induction period is also known, and its thermodynamic parameters can easily be calculated using the corresponding tabulated polynomials or the approximate formulas of the NFZT model. Note the quantity b affects only the wave parameter profiles in the induction zone. The wave velocity, the flow parameters at the ChapmaneJouguet (CJ) point, and the main heat release zone do not depend on b. Therefore, once the above conditions on the quantity b are satisfied, the specific form of the equation for calculating this quantity affects the wave parameter profiles in the induction zone only quantitatively and has no fundamental significance. As a rule, the chemical reaction rate increases at the end of the induction period. In this connection, the formula for calculating b should be chosen so as to satisfy the boundary conditions and so that the rate of decrease of b increases as the mixture moves in the induction zone. For example, we can use the formula
b ¼ 1 Y n ; n 1: Here Y is the fraction of the induction period:
Zt YðtÞ ¼ 0
dt ; ti
ti is the induction period at constant parameters of the gas. At t ¼ 0, we have Y ¼ 0; at t ¼ t*, Y ¼ 1. The value of ti is calculated by an experimental algebraic formula by Soloukhin (1970): E ti ¼ Aind ,exp ind ,½O2 1 eff ; RT where ½O2 eff ½mol=cm3 is the oxygen concentration, Aind ¼ 6.0$1015 mol s/cm3, Eind ¼ 33,200 cal/mol. In the main heat release zone, the molar mass m of the gas is calculated by integrating the differential equation (Nikolaev and Zak, 1988):
2 dm r2 m ¼ 4Kþ 1 dt mmax m 3=2 m AT 3=4 1 eq=T r 1 eE=RT ; mmin
q, mmin, mmax, Kþ, and A are constants of the model NFZT: E ¼ 4.7235,1012 erg/mol, q ¼ 3000 K, Kþ ¼ 0.6,1015 cm6/(mol2,s), A ¼ 6.045,1014 cm3/(mol,s,K3/4), mmin ¼ 18.817 g/mol, mmax ¼ 27.739 g/mol. The internal energy and the adiabatic index of the mixture are calculated by algebraic equations (Nikolaev and Zak, 1988). Table 1 shows the results of calculation of CJ parameters and Von Neumann spike parameters obtained by the proposed kinetics for ideal plane detonation wave. The corresponding results of the calculations, presented by Nikolaev and Topchiyan (1977) are shown. In this paper detailed equation of chemical equilibrium are used for calculation chemical composition and thermodynamic parameters of the gas in detonation wave. The induction period ti and induction zone length xind in Table 1 were estimated using formula by Soloukhin (1970). Similarly to the NFZT model, the proposed kinetic model for the detonation combustion of methane is consistent with the second law of thermodynamics and its constants have a clear physical meaning and are calculated from the tabulated thermochemical parameters of the mixture before the two-dimensional numerical calculation of the DW structure. Determinations of the CJ wave parameters and the DW cell size showed high accuracy of the model. An advantage of this kinetics is its simplicity and ease of use in multi-dimensional gas dynamics codes. This kinetic model is described in more detail in our papers (Fomin et al., 2013, 2014). Numerical calculation of Zel'dovich-von Neumann-Doering structure of one-dimensional stationary DW in methaneeair mixture based on this model of chemical kinetics has been performed by (Fedorov et al., 2014). The system of governing equations was closed by the wellknown thermal equation of state for an ideal gas. Numerical simulation was performed for a stoichiometric methaneeair mixture at normal initial conditions. DW propagation in a two-dimensional straightline channel was studied. The slip conditions (impermeable solid wall) were imposed on the upper, lower, and left boundaries of the channel according to Godynov (1976), and the right boundary was subjected to the condition of an undisturbed initial state of the gas. The detonation wave was initiated near the left closed end of the channel and propagated from left to right along the x axis. The size and energy of the initiation source (Trotsyuk, 1999) were chosen sufficiently large to ensure a supercritical regime of DW direct initiation in methaneeair mixture. The problem of determining the minimum critical initiation energy was not considered in this study. 2. Numerical method
where r is the density of the mixture, R is universal gas constant, E,
The resultant systems of equations were solved numerically. The following discretization in space and time was performed. A uniform motionless grid with the total number of cells Ny and cell size hy ¼ H/Ny was used in the y direction, where H is a height of
Table 1 Velocity and parameters of plane ideal DW in stoichiometric CH4eair mixture at T0 ¼ 298.15 K, p0 ¼ 1 atm. ChapmaneJouguet point parameters
The present approximate model Nikolaev and Topchiyan (1977)
D, m/s
TCJ, K
pCJ, atm
rCJ, g/cm3
uCJ, m/s
cCJ,e, m/s
mCJ, g/mol
1809 1801
2801 2780
17.37 17.19
2.051$103 2.048$103
809.0 803.8
999.7 997.2
27.13 27.18
Von Neumann spike parameters T10, K The present approximate model Nikolaev and Topchiyan (1977)
1553 1531
p10, atm 31.57 31.33
r10, g/cm3 3
6.874$10 6.913$103
u10, m/s
ti, ms
xind, cm
1510 1506
6.02 6.98
0.1796 0.2058
A.V. Trotsyuk et al. / Journal of Loss Prevention in the Process Industries 36 (2015) 394e403
channel. A moving grid with the total number of cells Nx was used in the x direction, Nx1 of these cells formed a uniform grid with the cell size hx. Almost in all calculations hx ¼ hy. The uniform grid moved together with the leading shock front of the DW in such a way that it covered flow regions with large gradients of parameters and/or shock waves, and thus the fine DW structure was captured. The remaining Nx2 ¼ NxNx1 cells formed a nonuniform moving grid that occupied the region between the closed left end of the channel and the beginning of the uniform grid. The use of adaptive moving grids with local refinement allowed us to have a fine resolution where necessary, using a significantly smaller number of grid cells, as compared with a uniform splitting of the computational domain. A typical of numbers of grid cells for high resolution computations were Nx ¼ 2000, Nx1 ¼ 1875, Ny ¼ 4000, and 256 coprocessors were used. The finite-volume Godunov-type scheme (Godynov, 1976) has been used with the MUSCL TVD reconstruction procedures. The fourth-order reconstruction by Yamamoto and Daiguji (1993), Daiguji et al. (1997) are used to calculate the fluxes along every axes in a case of the uniform grid zone, and a third-order scheme Chakravarthy and Osher (1985), Lin and Chin (1995) is used in the non-uniform grid zone along x axes. The advanced HLLC solver by Batten et al. (1997) is used for an approximate solution of the Riemann problem on the control volume boundaries. To implement this algorithm for obtaining numerical fluxes for the case of a chemically reacting mixture, the ‘‘energy relaxation method’’ by Coquel and Perthame (1998) is applied. This method eliminates the problems of numerical solution of the Riemann problem for a medium with a complicated nonlinear equation of state (including that with a variable ratio of specific heats for gaseous mixtures or any form of equation of state for condensed materials). A special work has been done to fit these algorithms to simulations of reacting detonations flows. The system of conservation equations for simulating the detonation process in a reacting medium is a hyperbolic system of Euler equations with stiff relaxation term in the right-hand sides. The stiffness of this system during its integration in time is caused by the source terms, which model chemical kinetics in detonation wave. In our numerical code we use the ASIRK method (Shen and Zhong, 1996) of the second order accuracy for integration on time. The semi-implicit method means that the stiff relaxation terms in this method are treated by an implicit manner, and the non-stiff convection terms by an explicit calculation. The time step is determined at each time layer of the solution from the stability condition, see (Godynov, 1976). The typical values of the Courant number in detonation problems simulations are CFL ¼ 0.3e0.4. The code is parallelized with MPI library using the domain decomposition technique. Two approaches were used. For a comparatively small number of CPUs (up to 20), the entire computational domain is decomposed into equal subdomains along one of the axes according to the number of CPUs used in the computation. The data exchange is performed only between the adjacent subdomains. The amount of the exchanged information depends on the type of used scheme. For our case it is three halo cells (three layers of boundary cells) from each side of virtual subdomain's boundary. But an amount of exchanged information also depends on strategy of decomposition of numerical computation domain to subdomains. The main idea is to minimize the volume of information transferred between subdomains for accelerating computations. For example, in the case of 2D simulations for a very big number of CPUs the decomposition only along one of the axes would result in the subdomains in form of very thin strips with a great number of boundary cells. It could lead to a prohibitive increase of volume of exchanged information. The optimal strategy is
397
to divide the entire computational domain into equal subdomains along the axes x and y according to the number of CPUs used in computation. The decomposition is performed in such a way, that subdomains have a square form in computational space Nx Ny. Although simple, this parallelization technique yields good results in terms of parallel efficiency and scalability.
Fig. 1. Detonation wave structure in a stoichiometric methaneeair mixture in channel with H ¼ 35 cm at xfront ¼ 1500 cm, pair of primary transverse wave AA and BB: (a) numerical Schlieren-visualization of all flowfield; (b) magnified fragment of the flowfield near an upper wall of a channel.
398
A.V. Trotsyuk et al. / Journal of Loss Prevention in the Process Industries 36 (2015) 394e403
Fig. 3. Detonation wave structure in stoichiometric methaneeair mixture in channel with H ¼ 35 cm at xfront ¼ 1500 cm “Low” resolution: hx ¼ hy z 160 mm. Numerical Schlieren visualization of all flowfield. AA and BB e primary transverse waves.
Fig. 2. Detonation wave structure in a stoichiometric methaneeair mixture in channel with H ¼ 35 cm at xfront ¼ 1500 cm: (a) normalized density flowfield; (b) temperature (K) flowfield.
3. Results of calculations The results of the simulation of the stoichiometric methaneeair mixture are presented in Figs. 1e3, 5e9, which show the fields of different parameters in the DW, including a numerical Schlieren visualization of the flowfield. The use of this parameter to visualize the structure of the detonation wave is due to its advantages over the other parameter, for example, pressure. The density gradient field accurately represents the flow pattern. Density changes
Fig. 4. Front-propagation velocity (solid line) as a function of front position in channel with H ¼ 35 cm, average velocity Daver ¼ 1840 m/s (dashed line), and ideal CJ velocity DCJ ¼ 1809 m/s (dash-dot line).
visualize not only the shape of the leading shock front of the DW and the structure of transverse waves, but also contact discontinuities. Significant changes in density at the end of the induction period, which are associated with the beginning of the main heat release stage, clearly show the two-dimensional structure of the DW induction region. In all the figures, the x and y coordinates are given in centimeters. The detonation cell size for the methaneeair mixture was determined in a procedure similar to that used in the numerical study of regular and double-scale cellular structures, see (Trotsyuk,
A.V. Trotsyuk et al. / Journal of Loss Prevention in the Process Industries 36 (2015) 394e403
Fig. 5. Detonation wave structure in a channel with H ¼ 35 cm at xfront ¼ 650 cm, Daver ¼ 1669 m/s, DCJ ¼ 1809 m/s, one primary transverse wave AA: (a) numerical Schlieren-visualization flowfield; (b) normalized density flowfield.
1999; Fomin et al., 2006; Vasil'ev et al., 2010). The initial height of the channel H was chosen rather arbitrarily; usually, it was assumed to be slightly higher than the presumed size of the cell a0. This simulation was started with H ¼ 15 cm, which was found later to be considerably smaller than the size of the cell. Further variation of H in the calculations affected the number of primary transverse waves remaining at the DW front after its initiation and the establishment of a self-sustained CJ detonation regime, thus, the cell size was determined. Figures 1 and 2 show the DW structure in a channel of height
399
Fig. 6. Detonation wave structure in a channel with H ¼ 25 cm at xfront ¼ 1140 cm, one primary transverse wave AA: (a) normalized density flowfield; (b) temperature (K) flowfield.
H ¼ 35 with the wave front located at xfr ¼ 1500 cm. Analysis of the change in the average velocity of the DW as it propagates away from the initiation source shows that a steady-state CJ detonation had already been established by this time. The figure shows a system of the two main (primary) transverse waves (TWs), AA and BB, which are located approximately at y z 12 cm and y z 19 cm. At this time, transverse waves move in the opposite direction to each other, in antiphase, as in mixtures with regular cellular
400
A.V. Trotsyuk et al. / Journal of Loss Prevention in the Process Industries 36 (2015) 394e403
Fig. 7. Detonation wave structure in a channel with H ¼ 40 cm at xfront ¼ 1170 cm, three primary transverse wave AA, BB and CC: (a) normalized density flowfield; (b) temperature (K) flowfield.
structure (compare with the figures of our papers (Trotsyuk, 1999; Vasil'ev & Trotsyuk, 2003). Moreover, the DW front has a system of secondary (smaller) transverse waves of different size, intensity, and direction of motion. One of these secondary TW at y z 30 cm is denoted as ‘aa’. The movement of these perturbations in the leading shock front of the DW is quite random. To identify the primary, most intense and longest, transverse waves (forming the
Fig. 8. Detonation wave structure in a channel with H ¼ 33 cm at xfront ¼ 600 cm, one primary transverse wave AA: (a) normalized density flowfield; (b) temperature (K) flowfield.
detonation cell) among these perturbations, it is necessary to analyze the flowfield of other parameters, in particular, the density and temperature, see Fig. 2. In all figures, the main transverse waves are denoted as AA, BB, etc. Fig. 1b presents an enlarged fragment of the DW front in the region y z 26 cm, which shows more details of the numerous secondary transverse waves forming the hierarchy of perturbations on the leading DW front decreasing in size. The secondary transverse waves that form the fine structure of the leading DW front,
A.V. Trotsyuk et al. / Journal of Loss Prevention in the Process Industries 36 (2015) 394e403
Fig. 9. Detonation wave structure in a channel with H ¼ 33 cm at xfront ¼ 1500 cm, pair of primary transverse wave AA and BB: (a) normalized density flowfield; (b) temperature (K) flowfield.
had been observed in the experiments of Manzhalei (1977). Fig. 1b shows the fine cellular (multifront) structure of the front one of the TW ‘aa’, located at y z 30 cm. The TWs on the front of wave ‘aa’ are denoted as 1,2,3,4. Thus, this simulation reproduces a characteristic feature of the detonation of hydrocarbon mixtures e instability of the front of the transverse waves itself at the leading shock front of the DW. This feature has been observed in soot imprints of typical trajectories of TWs in numerous experiments of Manzhalei & Mitrofanov (1973). In general, the shock-wave structure of the
401
DW front is extremely irregular (compare with the figures from Trotsyuk (1999); Vasil'ev & Trotsyuk (2003)). Furthermore, in Fig. 1b, one can see a significant amount of unburned pieces of gas at a considerable distance behind the leading shock front of the DW. The formation of these packages is in agreement with experimental observations (Subbotin, 1975; Kiyanda and Higgins, 2013). This phenomenon was not observed in simulations of regular mixtures based on hydrogen (Trotsyuk, 1999; Vasil'ev & Trotsyuk, 2003). The presence of these small-scale transverse waves on DW front and the unburned pockets of gas in the mixture under study leads to the need to maintain an extremely high spatial resolution of the computational grid in the induction zone and in the reaction zone behind the leading shock front. Most calculations were carried out with a characteristic mesh size hy ¼ hx z 80 mm, see Figs. 1, 2, 4e9. To study grid convergence of the numerical solutions the calculations have been performed with different sizes of the computational mesh. Numerical simulations with twice greater grid size hy ¼ hx z 160 mm in channel with H ¼ 35 cm show the same number of primary transverse waves, see Fig. 3. Of course, the fine structures of the flow were more blurred and smoothed as compared with fine grid simulations. So, the grid convergence of the numerical solutions was obtained. High spatial resolution led to a small allowed step in time and a large number of computational meshes. This, in turn, required parallelization of the numerical code and the use of a multi-processor supercomputer for numerical simulation of the detonation of the methane mixture. Fig. 4 shows a front-propagation velocity (solid line) as a function of DW front position after the initiation in a channel H ¼ 35 cm, parameters of calculations correspond to Figs. 1 and 2. Denotation has reached a quasi-steady propagation regime at the last 5 m. The average DW velocity obtained on this base is Daver ¼ 1840 m/s (dashed line). This value is slightly higher than the ChapmaneJouguet velocity (dash-dot line) DCJ ¼ 1809 m/s derived from thermodynamic calculations. Approximately the same absolute difference in the velocity of the multifront DW in a self-sustained CJ regime with DCJ was observed for hydrogeneoxygen mixtures with regular cell structure (Trotsyuk, 1999). The calculations were performed up to distances of more than ten meters, traveled by DW, because methane mixtures have extremely high critical initiation energy (Vasil'ev, 1998; Vasil'ev et al., 2000). Therefore, in the calculations, we had to use a highenergy initiator, which led to unsteadiness of the detonation process up to large distances. As an example, Fig. 5 shows the DW structure in the same channel with the front located at xfr ¼ 650 cm. In contrast to the pattern in Figs. 1 and 2, here we see only one main TW at y z 5 cm, and strong secondary TW ‘aa’ at y z 26 cm. The system of secondary transverse waves of different size at the leading front and the presence of regions of unburned gas are similar to the flow in s 1,2. At this point, the DW velocity averaged over two-dimensional pulsations, Daver ¼ 1669 m/s is much lower than the CJ velocity DCJ ¼ 1809 m/s. This suggests that at this distance, a self-sustained CJ detonation has not yet been established. Calculations with varying H showed that the cellular structure of the DW is significantly different from the structures observed with the channel sizes in Figs. 1 and 2. Thus, at H ¼ 25 cm, detonation flow with a single transverse wave at the DW front is formed. The normalized density field and the temperature field for the flow are shown in Fig. 6. Here we see one intense and extended TW AA and a no less extended, but much less intense degenerate (weak) TW ‘aa’ near the upper wall of the channel. The induction zone behind the leading shock front has a very irregular shape, and there are many large regions of unburned mixture in the region y z 2.5 cm. If H ¼ 40 cm, we have at least three main transverse waves, AA, BB, and CC (see Fig. 7); otherwise the structure of the front is similar to that shown in Figs. 1 and 2.
402
A.V. Trotsyuk et al. / Journal of Loss Prevention in the Process Industries 36 (2015) 394e403
To investigate the influence the numerical grid on the results of simulation of the DW structure, we performed additional calculations for H ¼ 33 cm with an increased number of cells along the x axis: Nx ¼ 4000 and Nx1 ¼ 3500. The number of cells along the y axis remained unchanged, Ny ¼ 4000, and hy ¼ hx in the zone with a fine mesh in the vicinity of the leading shock front. Figs. 8 and 9 show the simulation results for a channel with H ¼ 33 cm. As can be seen from the flowfields, the boundary between fine uniform grid and the nonuniform grid along the x axis is much less pronounced in these calculations than in those shown in the previous figures. At xfr ¼ 600 cm in Fig. 8, we have structure of an unsteady DW in the transient state, which is largely similar to the DW structure in Fig. 5. The multifront DW structure of steady self-sustaining CJ detonation is shown at xfr ¼ 1500 cm in Fig. 9. As in the case of detonation in the channel with H ¼ 35 cm (see Figs. 1 and 2), at the leading shock front of the DW there are two main transverse waves, AA and BB, and numerous secondary transverse waves of different lengths and amplitudes. In addition, the flow pattern has features that distinguish it from that in Figs. 1 and 2. It is seen that in the channel of this size, the movement of the transverse waves is less consistent; at the time presented in the figure, the transverse waves move in the same direction toward the bottom wall of the channel. Obviously, after the reflection of the BB wave from the bottom wall, it will collide with the AA wave not exactly on the axis of the channel. This will result in a distorted shape of the detonation cell, which is typical of all hydrocarboneair mixtures. In the case of mixtures with regular cellular structure, it could be argued that with this behavior of transverse waves, the channel width is slightly smaller than the transverse dimension of the cell a0. However, for the type of cellular structure considered here, this conclusion would be incorrect. From an analysis of the numerical results, the dominant transverse dimension of the irregular detonation cell formed by the main (primary) transverse waves AA and BB for stoichiometric methaneeair mixture can be determined to be a0 ¼ 34 ± 1 cm This value is in good agreement with the experimental value of the cell for this mixture (see (Nettleton, 1987; Vasil'ev, 1998; Vasil'ev et al., 2000; Bull et al., 1982; Knystautas et al., 1984; Moen et al., 1984)) and lies in the middle of the range of the experimental data. A numerical simulation of a two-dimensional structure of the DW in a stoichiometric methaneeoxygen mixture has been carried out. The fine multifront structure of the primary TWs is not clearly revealed in the present simulations. Calculations have been performed with a mesh size hy ¼ hx z 1.5 mm. Perhaps this resolution was not enough to resolve the fine structure of the transverse wave, requires further numerical studies. The size of the detonation cell for stoichiometric methaneeoxygen mixture from our numerical computations determines to be a0 ¼ 0.3÷0.35 cm. This value also is in good agreement with all known experimental data (Nettleton, 1987; Vasil'ev, 1998; Vasil'ev et al., 2000). 4. Conclusions An approximate two-step model (induction parameter model) for the kinetics of detonation combustion of methane was developed. Because of the simplicity and high accuracy (within the confines of the assumptions), the model can be used in multidimensional numerical simulations of detonation flows of methaneeoxygen and methaneeair mixtures. From the analysis of two-dimensional flow structures, the dominant size of the detonation cell for stoichiometric methaneeair mixture determines to be 34 ± 1 cm, and 0.3÷0.35 cm for methaneeoxygen mixture. These data reveal good agreement with available experimental data. This simulation takes into account the real thermophysical and
chemical properties of the mixture, and reproduces an irregular cell structure with all of its main features: random inconsistent movement of the main transverse waves; numerous secondary transverse waves forming the hierarchy of perturbations of decreasing size at the leading shock front of the DW; regions of unburned mixture at a considerable distance behind the DW front; fine (cellular) structure of the transverse wave. These calculations are the first in which the irregular structure of detonation wave is obtained using a two-step kinetic model. In experiments the fine structure of TW is manifested in the typical shape of transversewave trajectories recorded in soot imprints. The secondary waves at the leading shock front are responsible for the small-cell network inside the main cell on soot imprints in experiments with methane mixtures. Acknowledgments The authors gratefully acknowledge the financial contribution from the Russian Foundation for Basic Research (Grant No. 14-0100406 and Grant No. 14-03-00838), the Program of Basic Research No. 26.2 entitled “Combustion and explosion” of the Presidium of the Russian Academy of Sciences. The computations were performed on an MVS-100K supercomputer at the Joint SuperComputer Center (JSCC) of the Russian Academy of Sciences (Moscow). References Batten, P., Leschziner, M.A., Goldberg, U.C., 1997. Average-state Jacobians and implicit methods for compressible viscous and turbulent flows. J. Comput. Phys. 137, 38e78. Bull, D.C., Elsworth, J.E., Shuff, P.J., 1982. Detonation cell structures in fuel/air mixtures. Combust. Flame 45, 7e22. Chakravarthy, S.R., Osher, S., 1985. A new class of high accuracy TVD schemes for hyperbolic conservation laws. In: Proceedings of the AIAA 23rd Aerospace Sciences Meeting. Reno, Nevada, pp. 85e0363. AIAA Paper. Coquel, F., Perthame, B., 1998. Relaxation of energy and approximate Riemann solvers for general pressure laws in fluid dynamics. SIAM J. Numer. Anal. 35, 2223e2249. Daiguji, H., Yuan, X., Yamamoto, S., 1997. Stabilization of higher-order high resolution schemes for the compressible Navier-Stokes equation. Int. J. Numer. Methods Heat Fluid Flow 7, 250e274. Fedorov, A.V., Fomin, P.A., Tropin, D.A., 2014. Simple kinetics and detonation wave structure in a methane-air mixture. Combust. Explos. Shock Wave 50, 87e96. Fickett, W., Davis, W.C., 1979. Detonation. University of California Press, Berkeley, CA. Fomin, P.A., Trotsyuk, A.V., 1995. An approximate calculation of the isentrope of a gas in chemical equilibrium. Combust. Explos. Shock Waves 31, 455e458. Fomin, P.A., Trotsyuk, A.V., Vasil’ev, A.A., Mitropetros, K., Hieronymus, H., Roekaerts, D., 2006. Model of chemical reaction kinetics for calculating detonation processes in gas and heterogeneous mixtures containing hydrogen peroxide. Combust. Sci. Technol. 178, 895e919. Fomin, P.A., Trotsyuk, A.V., Vasil’ev, A.A., 2013. Approximate models of chemical kinetics for detonation processes in mixtures of CH4, H2O2 and O3 with air. In: Proceedings of the 24th ICDERS. July 28eAugust 2, 2013, Taipei, Taiwan. Paper 0161. Fomin, P.A., Trotsyuk, A.V., Vasil’ev, A.A., 2014. Approximate model of chemical reaction kinetics for detonation processes in mixture of CH4 with air. Combust. Sci. Technol. 186, 1716e1735. Gamezo, V.N., Desbordes, D., Oran, E.S., 1999. Formation and evolution of twodimensional cellular detonations. Combust. Flame 116, 154e165. Gamezo, V.N., Vasil’ev, A.A., Khokhlov, A.M., Oran, E.S., 2000. Fine cellular structures produced by marginal detonations. In: Proceedings of the 28th Symposium (International) on Combustion. The Combustion Institute, Pittsburgh, pp. 611e617. Godynov, S.K. (Ed.), 1976. Numerical Solution of the Multi-dimensional Problems of Gaseous Dynamics. Nauka Press, Moscow (in Russian). Kessler, D.A., Gamezo, V.N., Oran, E.S., 2010. Simulations of flame acceleration and deflagration-to-detonation transitions in methane-air systems. Combust. Flame 157, 2063e2077. Kessler, D.A., Gamezo, V.N., Oran, E.S., 2011. Multilevel detonation cell structures in methane-air mixtures. Proc. Combust. Inst. 33, 2211e2218. Kiyanda, C.B., Higgins, A.J., 2013. Photographic investigation into the mechanism of combustion in irregular detonation waves. Shock Waves 23, 115e130. Knystautas, R., Guirao, C., Lee, J.H., Sulmistras, A., 1984. Measurements of cell size in hydrocarbon-air mixtures and prediction of critical tube diameter, critical initiation energy, and detonability limits. Prog. Astronaut. Aeronaut. 94, 23e37.
A.V. Trotsyuk et al. / Journal of Loss Prevention in the Process Industries 36 (2015) 394e403 Korobeinikov, V.P., Levin, V.A., Markov, V.V., Chernyi, G.G., 1972. Propagation of blast waves in a combustible gas. Astronaut. Acta 17, 529e537. Lin, S.-Y., Chin, Y.-S., 1995. Comparison of higher resolution Euler schemes for aeroacoustic computations. AIAA J. 33, 237e245. Manzhalei, V.I., Mitrofanov, V.V., 1973. The stability of detonation shock waves with a spinning configuration. Combust. Explos. Shock Waves 9, 614e620. Manzhalei, V.I., 1977. Fine structure of the leading front of a gas detonation. Combust. Explos. Shock Waves 13, 402e404. Moen, I.O., Funk, J.W., Ward, S.A., Rude, G.M., Thibault, P.A., 1984. Detonation length scales for fuel-air explosives. Prog. Astronaut. Aeronaut. 94, 55e79. Nettleton, M.A., 1987. Gaseous Detonations: Their Nature, Effects and Control. Chapman and Hall, London, New York, ISBN 0-412-27040-4. Nikolaev, YuA., 1978. Model of the kinetics of chemical reactions at high temperatures. Combust. Explos. Shock Waves 14, 468e471. Nikolaev, YuA., Fomin, P.A., 1982. Analysis of equilibrium flows of chemically reacting gases. Combust. Explos. Shock Waves 18, 53e58. Nikolaev, YuA., Topchiyan, M.E., 1977. Analysis of equilibrium flows in detonation waves in gases. Combust. Explos. Shock Waves 13, 327e338. Nikolaev, YuA., Zak, D.V., 1988. Agreement of models of chemical reactions in gases with the second law of thermodynamics. Combust. Explos. Shock Waves 24, 461e464. Shen, J.W., Zhong, X., 1996. Semi-implicit Runge-Kutta schemes for nonautonomous differential equations in reactive flow computations. In: Proceedings of the AIAA 27th Fluid Dynamics Conference. New Orleans, LA, June 17-20, 1996. AIAA Paper 961969. Soloukhin, R.I., 1970. Measurement methods and basic results of shock-tube experiments. In: Proceedings of the 7th Int. Shock Tubes Symposium. 23e25 June
403
1969, Toronto, Canada. University of Toronto Press, Toronto and Buffalo, ISBN 08020-1729-0, pp. 662e706. Subbotin, V.A., 1975. Collision of transverse detonation waves in gases. Combust. Explos. Shock Waves 11, 411e414. Trotsyuk, A.V., 1999. Numerical simulation of the structure of two-dimensional gaseous detonation of an H2eO2eAr mixture. Combust. Explos. Shock Waves 35, 549e558. Vasil'ev, A.A., 1998. Effect of nitrogen on multifront detonation parameters. Combust. Explos. Shock Waves 34, 72e76. Vasil'ev, A.A., Nikolaev, YuA., 1976. Model of a multifront gas detonation cell. Combust. Explos. Shock Waves 12, 744e754. Vasil’ev, A.A., Grigor’ev, V.V., 1980. Critical conditions for gas detonation in sharply expanding channels. Combust. Explos. Shock Waves 16, 579e585. Vasil’ev, A.A., Valishev, A.I., Vasil’ev, V.A., Panfilova, L.V., Topchian, M.E., 2000. Detonation hazards of methane mixtures. Arch. Combust. 20, 31e48. Vasil'ev, A.A., Vasil'ev, V.A., Trotsyuk, A.V., 2010. Bifurcation structures in gas detonation. Combust. Explos. Shock Waves 46, 196e206. Vasil'ev, A.A., Trotsyuk, A.V., 2003. Experimental investigation and numerical simulation of an expanding multifront detonation wave. Combust. Explos. Shock Waves 39, 80e90. Voitsekhovskii, B.V., Mitrofanov, V.V., Topchian, M.E., 1966. The Structure of a Detonation Front in Gases. WrightPatterson (AFB AD633821). Yamamoto, S., Daiguji, H., 1993. Higher-order-accurate upwind schemes for solving the compressible Euler and Navier-Stokes equations. Comput. Fluids 22, 259e270. Zel'dovich, YaB., Kompaneets, A.S., 1955. The Theory of Detonation. Gostekhteorizdat, Moscow (in Russian).