Abstract
The work is dedicated to the numerical study of detonation wave initiation and propagation in the variable cross-section axisymmetric channel filled with the model hydrogen-air mixture. The channel models the large-scale device for the utilization of worn-out tires. Mathematical model is based on two-dimensional axisymmetric Euler equations supplemented by global chemical kinetics model. The finite volume computational algorithm of the second approximation order for the calculation of two-dimensional flows with detonation waves on fully unstructured grids with triangular cells is developed. Three geometrical configurations of the channel are investigated, each with its own degree of the divergence of the conical part of the channel from the point of view of the pressure from the detonation wave on the end wall of the channel. The problem in consideration relates to the problem of waste recycling in the devices based on the detonation combustion of the fuel.
1. Introduction
Mathematical modeling of two-dimensional flows with detonation waves (DW) as a result of solving of the Euler equations for the inviscid gas supplemented by the chemical reaction kinetics model originates in the late 1970s [1, 2]. Since then with the development of the computational methods for the solution of gas dynamics problems, the improvement of the kinetic schemes, and the growth of available computational resources, the continuous clarification of qualitative and quantitative characteristics of the process has occurred (for example, obtaining the three-dimensional spin mode of propagation, the thin structure of cellular detonation, and the detonation limits [3]). One can agree with [4] in which the current understanding of the mechanisms of DW propagation is put into direct dependence on the maximum computer performance achieved at the moment. At the same time, the growth of computational powers and as a result the possibility of carrying out the calculations with more and more detailed spatial-temporal resolution revealed a number of unusual effects of the detonation modeling that have not had full explanation yet. The effects include, for example, the possible detonation decay for the kinetics parameters close to the hydrocarbon fuels in two-dimensional calculations of the long time propagation of DW in the plane channel [5]. Thus, despite almost forty-year history of computational works in the field of detonation, a number of fundamental questions of DW mathematical modeling remain.
Numerical research of problems of DW initiation and propagation in technical systems and facilities involves the necessity of consideration of computational domains of a complex shape. This fact leads to the set of questions related to the construction of computational grids in such domains, construction of monotone in some sense schemes of high approximation order on the selected class of grids, and also infrastructure problems with large amounts of poorly structured data including visualization issues. Despite all the achievements in the field of multiprocessor computing, the numerical studies of problems of physical and chemical hydrodynamics with a volume exceeding 108 computational cells are at the limit of the capabilities of the researchers due to the factors listed above. Thus, the development of computational technologies for the modeling of high speed flows with chemical reactions in areas of complex shape remains the actual issue.
In [5, 6] for carrying out three-dimensional computations of the DW initiation and propagation in complex-shaped domains, the block-structured grids and the classical first approximation order scheme of S. K. Godunov [6] or its modification based on some interpolation schemes for the approximation order increase [5] are used. At the same time, in relation to gas dynamics problems of chemically inert media the apparatus of high approximation order schemes on completely unstructured computational grids including shock wave problems has been developed; see the most known paper [7]. The analysis of publications about the usage of such approaches to the study of high speed flows with chemical reactions gives a few works [8–10]. In [9, 10], the finite element approach is applied. The most probable reason is the fact that although the integration of the gas dynamics equations is a key element of the computational algorithm for the simulation of high speed flows with chemical reactions, the presence of strongly nonlinear sources in the right-hand sides of the equations significantly reduces the possibilities of applying many numerical methods that successfully deal with the gas dynamics problems of inert media.
2. Statement of the Problem
The axisymmetric channel of a variable cross-section filled with quiescent model hydrogen-air mixture under normal conditions (see Figure 1) that simulates the facility for utilization of used automobile tires [11, 12] is considered. The geometry of the problem corresponds to that considered in [12]. The channel consists of a short narrow segment for detonation initiation (preliminary chamber), a conical part for DW passing into a segment of a larger diameter (detonation chamber), and a working chamber. The cases of different expansion angles of the conical part corresponding to the values cm, 30 cm, and 50 cm are investigated. For the detonation initiation, the pressure 40 atm and temperature 1500 K are set in a region with the length of 4 cm in the preliminary chamber (colored gray in Figure 1). The boundary conditions of impermeability are set on all boundaries of the computational area. The pressure is recorded by the sensors and . The objectives are the description of the mechanism of the development of detonation process in the channel of such geometry, the comparative analysis of the pressure curves at the sensors under variation of the geometrical parameter , comparison of the obtained results with the data from [12], and research of the applicability of the developed computational algorithm of second approximation order on completely unstructured computational grids for solving large-scale practical tasks. Note that the problem in consideration relates to the fundamental problem of the investigation of mechanism of DW passing from a narrow channel into the wide, which was discussed in theoretical and experimental works of many authors; see, for example, [13, 14]. The theory of this phenomenon which would predict the realized regime depending on the geometrical characteristics of the channel, mixture properties, and degree of DW overdriving has not been built.
3. Mathematical Model
The mathematical model is based on two-dimensional system of Euler equations written in the cylindrical frame under the assumption of the axial symmetry of the flow supplemented by one-stage kinetics of hydrogen combustion in air [15]:Here is the time, is the gas density, and are the radial and axial velocity components, is the pressure, is the specific internal energy of the gas, is the total energy per the unit of volume, is the heat release of the chemical reactions, is the mass fraction of the reactive component of the mixture, is the rate of chemical reactions, is the preexponent factor, is the activation energy, is the molar mass, is the universal gas constant, and is the temperature. The gas is considered to be ideal with the specific heat ratio . The following parameters of the mixture are used [15]:The used kinetic model gives the following values of Chapman-Jouguet (CJ) parameters:and the half-reaction length in the stationary Zeldovich-von Neumann-Doering solution equals 0.2 mm.
4. Computational Algorithm
The main feature of the computational technique is the use of completely unstructured computational grids with triangular cells. A Delaunay triangulation based on the arrangements of the computational mesh nodes to satisfy the Delaunay condition is carried out to construct the grid.
The computational algorithm is based on the Strang splitting principle in terms of physical processes [16]. When passing from one time layer to another one, one first integrates the gas dynamics equations without considering the chemical reactions ( in (1)) and thereby performs the first stage of the splitting procedure. On the second stage, the chemical reactions are taken into account without considering the convection (the second stage of splitting).
The spatial part of (1) is discretized using the finite volume method:Here, the spatial index corresponds to the number of calculated cells of area , and is the radius of the center mass point of the triangular. The summation is conducted over all edges of the cell , is the length of the edge with the index , and is the vector of numerical flux through the edge in the direction of the unit normal that is external to the edge. To calculate the fluxwe pass from the initial global frame to the local one that is associated with the outer normal to the edge and the tangent vector . Denote the matrix of transformation as . Vector denotes the conservative variables vector in the triangular that has the common edge with the considered triangular . The flux in the local frame is calculated using AUSM scheme [17] extended for the case of a two-component mixture.
For the approximation order increase, the following reconstruction of the grid functions is applied [18]. The numerical flux in the local frame is computed using the following parameters from the left and the right sides of the edge (see Figure 2):Here denotes the set of triangles that have a common node with the triangle where the node is an opposite to the edge . Similar denotes the set of triangles that have a common node with the triangle where the node is an opposite to the edge . The minmod limiter applied component-wise is considered as a limiter function :
Time integration is carried out using explicit Runge-Kutta method of the second approximation order:The upper tilde denotes that the obtained solution is the results of the first stage of the splitting procedure. The time step is chosen dynamically from the stability condition.
On the second stage of the algorithm, the system of ordinary differential equations of chemical kinetics for variable and temperature in each computational cell of the grid is solved.
The estimation of the practical approximation order of the algorithm on the problem of isentropic vortex evolution [7] gave the value near 2. Note that the computations for few individual triangles near the places of sharp changes of the channel geometry were carried out using the first approximation order scheme to ensure stability.
The computations were carried out on computational grids with cells number between 3 and 4 million depending on the geometry. Because of the large scale of the considered problem, the grid resolution was too coarse with respect to the half-reaction length for the considered kinetics model in order to obtain effects associated with the DW front instability, but sufficient for estimating the integral gas dynamics parameters of the problem.
5. Results of the Numerical Experiments
Figure 3 illustrates the main stages of the process. The direct initiation of detonation in the preliminary chamber leads to the overdriven regime of DW propagation (see Figures 3(a) and 4). The pressure in the rarefaction wave when DW enters the conical part is reduced from 30 atm to the pressures close to (see Figure 3(b)). Long away behind the DW front the relatively weak transverse waves are formed. Much more strong transverse wave is formed when DW passes from the conical part of the detonation chamber to a segment of a constant cross-section (see Figure 3(c)).
(a) 0.06 ms
(b) 0.18 ms
(c) 0.30 ms
(d) 0.38 ms
(e) 0.43 ms
(f) 0.54 ms
In the case of cm, the DW first reaches the side wall of the working chamber, then the interaction of the DW with the end wall occurs (see Figure 5(c)).
(a) сm, 0.36 ms
(b) сm, 0.47 ms
(c) сm, 0.56 ms
In the case of cm, these events occur at one time moment (see Figure 5(b)). In the case of сm, the DW first reaches the end wall of the channel (see Figure 5(a)). In this case, the wave pattern after the cumulation of the reflected wave at the symmetry axis is the most complicated since the wave interacts with the set of transverse waves which were formed in the conic expansion part (see Figures 3(d) and 3(e)). The cumulation occurs at different time moments with respect to the moment of the DW arrival to the end wall of the working chamber for different channel configurations. The following stages of the process include the interference of the waves reflected from the end and the side walls of the working chamber. For example, for the case of cm, the highest pressure after the interference is obtained in the “corner” point (0.4, 0.8) (see Figure 3(e), red colored area). The next noticeable event is the interaction of the jet from the reflected wave cumulation at the symmetry axis with the end wall (see Figure 3(f), red colored area near the point (0.0, 0.8)).
Figure 6 illustrates the comparison of the pressure curves at the sensors and for the case of сm with the data from [12]. The curves are aligned at the moment of the arrival of the DW at the sensors since there is no exact information about the initiation parameters in [12]. The peak pressure at the sensor is approximately 2 times higher than the analogous value at the sensor , because of the interference of the waves reflected from the side and end walls of the working chamber (see Figure 3(e)). For each of the three considered configurations, the pressure impulses applied to the end wall of the working chamber were measured by calculating the area under the pressure curve on the sensor . The greatest impulse is realized in the case of cm. In the case of the smallest angle of the conical expansion part, the DW undergoes the least attenuation when the DW passes into the working chamber in comparison with other configurations. In the case of sm, the impulse is reduced by 6%. In the case of sm, the impulse is reduced by 13%.
(a)
(b)
6. Concluding Remarks
The work demonstrates the possibility of the numerical modeling of initiation and propagation of gaseous detonation waves in the two-dimensional axisymmetric statement on completely unstructured computational mesh with triangular cells. The computational algorithm based on the finite volume method of the second approximation order is presented and described in detail.
The work also demonstrates the efficiency of the suggested computational algorithm for the numerical simulation of the detonation waves propagation in the large-scale channel of the variable cross-section which models the facility for crushing tires. The correctness of the obtained results is confirmed by the comparison of the parameters of the self-sustained detonation wave with the analytic estimates of the Chapmen-Jouguet parameters and the comparison with the computational results of other authors. The mechanism of the process development for three device configurations is described. The pressure impulse acting at the end wall of the working chamber is maximal at the minimum angle of the conical expansion solution of the detonation chamber.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper or the received funds.
Acknowledgments
The work is partially funded by the Federal Agency for Scientific organizations and the Russian Foundation for Basic Research (Grant no. 15-31-70004).