Abstract
Numerical simulation of an oblique detonation induced by a wedge is performed to investigate the formation and evolution of the oblique cellular detonation structure and the quantitative comparison of the cellular structure of a normal and an oblique detonation. The compressible reactive Euler equations are solved using the seventh-order WENO scheme on an adaptive mesh based on the open source program Adaptive Mesh Refinement in Object-oriented C++ (AMROC). The numerical results show that there are two sets of transverse waves, the left-running transverse waves (LRTW) and the right-running transverse waves (RRTW), which form the oblique cellular detonation structure. Although both sets of transverse waves are convected downstream, they propagate at almost the same relative velocity in the opposite direction. The LRTWs start in the transition zone because of the detonation instability. However, the RRTWs form due to the interaction between the deformed detonation front and the reflected shock wave from the wedge. For the same degree of overdrive, numerical simulations reveal that the characteristic cell size of an oblique detonation is almost the same as that of a normal detonation.
1. Introduction
An oblique detonation engine is a potential propulsion device based on the stationary oblique detonation theory [1–6]. The oblique detonation in the combustion chamber is induced by supersonic reactive inflow over a wedge, which is relevant to practical configuration for hypersonic propulsion operation [5]. The cellular structure of oblique detonation is particularly important in terms of studying the self-sustaining and stable mechanism of the oblique detonation on the wedge.
A wealth of analytic [1, 3, 4, 7], experimental [2, 8–10], and numerical studies on the fundamental issues about the ODW induced by the wedge can be found in the literatures. Due to the complex physical phenomena of oblique detonation and the very demanding requirements on experimental conditions, high-resolution numerical simulation is still a more effective research method on oblique detonation initiation mechanism than an analytic and experimental method. In the aspect of numerical simulations, there are a lot of authors investigating the formation of a stable detonation induced by a wedge. Li et al. [11] first addressed that the two-dimensional oblique detonation wave (ODW) structure consists of a thermally neutral induction zone, a transition zone, and an oblique detonation surface. Another characteristic of the ODW is the instability of the formed detonation. Choi et al. [12, 13] numerically showed that the activation energy has a strong effect on the ODW instability, and these unstable oblique detonation structures cannot be captured without sufficient numerical resolution. In their simulations, only LRTWs were observed. By using a higher-order numerical scheme, Gui et al. [14] revealed two families of transverse waves, i.e., LRTWs and RRTWs. Both of the transverse waves are convected downstream along the ODW. Teng et al. [15–17] conducted simulations of oblique detonation structures with a semi-infinite wedge and cones using with a one-step hydrogen-air chemistry model. They studied the effect of the activation energies, wedge angles, and conical shocks on the ODW structures and instabilities. In their simulations, they found that high activation energy and a low wedge angle are beneficial to the LRTW formation, and the RRTW formation is more affected by the activation energy and the existence of the LRTW. Verreault et al. [18] used high-resolution grids with an adaptive mesh refinement (AMR) method to conduct the simulations of the ODW. They observed that the formation of RRTW comes from the compressive wave, which may be amplified along the oblique surface. The criteria and location of transverse wave formation are numerically studied by Wang et al. [19] and Miao et al. [20], and both studies show that transverse waves at the region of the primary triple exist in the abrupt oblique shock wave (OSW) to ODW transition but cannot stand in the smooth transition. Recently, most of numerical simulation have focused on the criteria of detonation initiation and propagation modes on the different geometric wedge [21–23]. Huang et al. [24] simulated the wedge-induced oblique detonation transitions to investigate the nondetonation reaction zone structures that include induction and deflagration regions and found that the hypothesis of constant pressure combustion along the inclined wedge surface is more appropriate than the hypothesis of constant volume combustion.
Although these two families of transverse waves were observed and studied, their formation and evolution are not yet fully understood. For instance, is the characteristic cell size of an ODW structure the same as that of a normal detonation wave for the same degree of overdrive? Are the velocities of the LRTW and RRTW the same? What caused the LRTW to be initiated earlier than the RRTW? What are the origin and formation mechanism of these two transverse waves? In this paper, numerical investigations on the cellular structures of ODW are conducted to address these issues.
As some authors [4, 25] pointed out, the difficulty in investigations of this problem is that a long computational domain and high resolution grids are required to allow the cellular structures to be fully developed in order to reveal the final transverse wave structure of the ODW. In this study, numerical simulations are conducted with an AMR method, and using a WENO scheme for the convective flux discretization that is of order seven in space and thereby high resolution grids can cover the large enough computational domain to investigate the formation and evolution of these two families of transverse waves.
The paper is organized as follows. Section 2 summarizes the governing equations for the simulations and the numerical methods used to solve these equations. In Section 3, the method is validated using a lower activation energy case. The formation and evolution of the two families of the transverse waves are investigated in detail with numerical smoke foil record. Furthermore, the cellular structure of an ODW and a normal detonation wave is compared quantitatively. Section 4 contains our conclusions.
2. Numerical Methods
2.1. Governing Equations
The reactive Euler equations can be written as where denotes the density, denotes the component of the velocity in the direction, p is pressure, denotes the Kronecker symbol ( if , and otherwise). denotes the mass fraction of the product, and is the total energy is given by where is the amount of heat release, is the universal gas constant, is the temperature, and is the ratio of specific heats at a constant pressure and volume. The chemical reaction is modeled by a single step irreversible reaction as follows: where is a preexponential factor (frequency factor) and is the activation energy. The equation of state for the ideal gas mixture has the following form: where is the gas constant.
The thermodynamic and hydrodynamic quantities is nondimensionalized using the free stream values for density, pressure, and temperature, so resulting in the following normalizations for the thermodynamic and hydrodynamic quantities:
The distances and times were further normalized in the problem as follows: where is the frequency factor and is the half reaction length. Hence, the nondimensional transport equation for the product mass fraction becomes where the nondimensionalized value of the preexponential factor is chosen so that the half-reaction length of a planar CJ detonation is unity.
2.2. Numerical Scheme
The numerical simulations are carried out using the open-source code AMROC [26, 27] and based on a structured adaptive mesh refinement (SAMR) framework [28]. An operator splitting technique (or the method of fractional steps) for the computation of the time-dependent reactive flow (see [29]) is used to solve Equation (1). This technique allows a decoupled treatment for the time-implicit discretization of the local source term and the time-explicit discretization of the hydrodynamic transport term. The convection terms in Equation (1) are discretized with the seventh-order WENO-symmetric-order optimized (WENO-SYMOO) scheme as shown by Martín et al. [30]. The optimal third-order strong stability preserving (SSP) Runge-Kutta scheme [31] is used for time integration in combination with time-splitting and the fourth-order accurate semi-implicit GRK4A method [25] for source term integration.
2.3. Simulation Settings and Computational Domain
Figure 1 shows the two-dimensional computational domain and boundary conditions for the wedge-induced oblique detonation simulations. The horizontal wedge surface is employed, and the Cartesian grids are aligned with the wedge surface. The free stream conditions are imposed in both left and upper boundaries of the domain as follows: where and denote the components of the inflow velocity in the and directions, respectively. Slip boundary conditions are employed at the wedge surface. The outflow conditions are imposed on the right boundary. The wedge angle is . The chemistry parameters are set as , , and .
(a)
(b)
3. Results and Discussion
3.1. Numerical Validation and Grid Independence Study
To validate the numerical scheme, the numerical results are compared with those from the method of characteristics (MoC) [32]. The flow behind an attached oblique detonation is supersonic, and the Method of Characteristics (MoC) is a common technique to solve supersonic flow. Verreault et al. [32] developed a MoC program for simulation of supersonic uniform reactive flow passing through a semi-infinite wedges, which is to validate the CFD codes, explore the formation process of an ODW, and determine the steady-state oblique detonation angle.
In the simulations, because the oblique detonation remains stable for a low activation energy, the activation energy is decreased to . The computational domain, boundary conditions, and grid resolutions are the same as those in the higher activation energy () case. Figures 2(a) and 2(b) show the shock trajectory and the profile of the reaction progress variable along the wedge surface. It is seen that the AMROC solutions are in excellent agreement with the MoC results [32].
(a)
(b)
To test the mesh convergence, three different refinement levels of grid resolution with a ratio of 2 (i.e., 4, 5, and 6 refinement level mesh) are considered for the higher activation energy case (). Figure 3 shows the density and temperature profiles with the three grid systems. It is found that the density and temperature profiles are almost the same in all three grid levels. Figure 4 shows the density gradient contour with 4, 5, and 6 refinement level grid in the range of . The transverse waves can be observed clearly in the last two grid systems. To balance the calculation resolution and cost, the 5 refinement level grid system is selected for all the following numerical simulations. The 5 levels of adaptive grids in the range of are shown in Figure 1(b).
(a)
(b)
(a)
(b)
(c)
3.2. Analysis on the Oblique Cellular Detonation Structures
When a reactive supersonic inflow passes over a wedge, under the following assumption: the is constant before and after the detonation wave; the gas before and after the detonation wave is an ideal gas; and the combustion surface of the chemical reaction is infinitely thin and the heat is released instantaneously, the induced oblique detonation angle is given by the following equation [4]: where , where is the inflow Mach number. The CJ detonation speed can be calculated by the following equation [26]:
By solving Equations (10) and (11) based on the given inflow and chemical parameters, the theoretical oblique detonation angle, which is 40.0°, and the degree of overdrive, which is , where is the velocity normal to the ODW, could be obtained.
The computational domain is (-0.5, 30) and (0, 10) in the and directions, respectively. The tip of the wedge is located at (0, 0). The adaptive grid refinement algorithm is employed. The adaptive grid has () cells at the coarsest level and four additional levels of refinement with a ratio of 2. This provides a resolution of 640 points per half-reaction length (pts/hrl) of a CJ detonation, or 55 pts/hrl of an overdriven ODW with .
Figure 5 shows the pressure contour with at time . A leading oblique shock forms when the supersonic inflow impinges the wedge surface. The flow field downstream is not stable, so a weak compression wave is formed between the oblique shock and the flow field downstream to balance the pressure in the flow field. A triple point is formed on the shock front and moves downstream, and a thin chemical reaction layer is formed on the wedge surface. As time goes on, the oblique shock is strengthened, which make the reaction gradually accelerate and shorten the reaction zone behind the oblique shock. As a result, the oblique detonation forms eventually. The region ahead of the oblique detonation is known as the induction zone and transition zone [11, 33].
Figure 6 shows the leading shock trajectory with at time . The red line represents the leading shock trajectory, and the black line represents a straight line with a theoretical oblique detonation angle of 40°. Point A is the onset of the oblique detonation, which means that the region ahead of point A is composed of the induction zone and the transition zone. From the figure, it shows that the overdriven oblique detonation is not a straight line but an oscillating curve. The oscillation starts at point B in the transition zone due to the detonation instability.
In order to study the evolution of the instabilities in time from the numerical simulations, a series of individual pressure contours at the leading shock taken from different times are to be superimposed and translated. The translation direction is opposite to that of the free stream, and the distance is obtained by the velocity of the freestream multiplied by the time interval between the two figures. Figure 7 shows the resulting cellular oblique detonation structure, constructed with 100 frames. The time starts at , and the time interval is .
From the trajectories of the triple points in Figure 7, there are four zones on the leading oblique shock. Zone A in the range of is the induction zone. Zone B in the range of is located in the transition zone. Both of zones A and B are in a steady state. In zone C in the range of , detonation instabilities are formed and only the LRTWs are recognized. In the range of (or zone D), both of the transverse waves (LRTWs and RRTWs) are present. These two sets of transverse waves form the cellular structure. Also, from Figure 7, the triple points of the same family on the ODW move along parallel straight lines at a constant speed. Based on the traveling distance and time of a triple point in the figure, the relative speed of the triple point on the LRTW, which is 2.79, can be obtained. Since the inflow velocity that is tangential to the ODW is 6.13, the triple points on the LRTW move downstream with an absolute speed of 3.34 along the front of the ODW. Similarly, the triple points on the RRTW move downstream as well with a relative speed of 2.77 and an absolute speed of 8.90. That means that the two sets of transverse waves propagate at almost the same relative velocity in the opposite direction along the detonation front even though both are convected downstream.
Figure 8 shows the pressure distributions along the leading shock at . From Figure 8(a), it can be seen that the oscillations start at around and are amplified along the direction. Figure 8(b) displays the pressure distributions along the leading shock at the range of . In this figure, three triple points (LTW1, LTW2, and LTW3) can be clearly identified. The facing upstream pressure distributions also infer the transverse wave is the LRTWs.
(a)
(b)
Figure 9 shows the temporal evolution of the pressure density contour and the LRTWs at different times. In these figures, the line refers to the reaction front location where the mass fraction of the product . In Figure 10(a), due to the detonation instability in the transition zone, the oscillations of the ODW at points LTW1, LTW2, and LTW3 are visible. As a result, these oscillations cause the shock angle to fluctuate. Because the oscillations are still small at the current stage, the triple points have not been generated yet. As time evolves, the oscillations are amplified, as shown in Figures 10(b) and 10(c). The triple points are generated where the shock angles are larger and their intensities are stronger. The transverse waves and slip lines attached to the triple points are formed. Because the leading shock is strong near the triple point, the pressure at the triple point is higher than the theoretical overdrive oblique detonation pressure. It leads to that the detonation is formed and the reaction front is coupled with the leading shock near the triple point. After the detonation, the heat release is unable to maintain the high pressure at the triple point, which will lead to the oblique detonation angle to gradually decrease downstream until the next triple point. At this region, the nonreactive leading shock is decoupled with the reaction front.
(a)
(b)
(c)
(a)
(b)
Figure 10 shows the flow structure with an LRTW at . In Figure 10(b), it shows second triple point B along the transverse wave ABD behind primary triple point A. B is caused by the interaction of the strong shock wave BC with the transverse wave. ACE is the slip line emanating from primary triple point A. Figure 10(a) displays the transverse wave BD is stronger than the transverse wave AB. The combustible gas flowing through the leading shock is ignited because they are compressed by the transverse wave AB.
In Figure 11, at location , the interaction between the leading shock front and the reflection waves from the wedge surface causes the leading shock surface to deform and generates the inflection point as a result, which is the starting position of the RRTW (also see Figure 7).
Figure 12 shows the temporal evolution of the density gradient contour and the RRTW at different times. Figure 13 presents the evolution of the RRTW in the numerical smoke foil. Figure 12(a) displays the starting position (RTW) of RRTW at . As time elapses, the RRTW is strengthened after every collision with different LRTWs (see Figures 12(b)–12(d)). As a result, the point RTW is gradually evolved to the triple point from the inflection point. At , the primary and second triple points of the RRTW are clearly visible. AE is the slip line of the RRTW. Hence, similar to the LRTW, the triple point structure of the RRTW is a double Mach reflection structure.
(a)
(b)
(c)
(d)
Teng et al. [15] revealed the existence of the unreacted region behind the cellular oblique detonation. However, they did not analyze the reasoning behind it. Figure 14 displays the formation of an unreacted region behind the cellular oblique detonation. After the cellular oblique detonation is formed, Figure 14(a) shows the flow structure when the two triple points (LTW and RTW) hit in one single point. After the collision, both triple points break down and the old slip lines are torn off immediately, as shown in Figure 14(b). It is similar to the reinitiation procedure in a normal detonation [26]. An unreacted region forms when the combustion zone is separated behind the leading shock wave from the reaction regions. It is bounded to the down by the old slip lines. The unreacted region is burned quickly because it is surrounded by the high temperature region. Figure 14(c) shows the formation of the new triple points, transverse waves, and slip lines after the collision.
(a)
(b)
(c)
3.3. Comparison of the Cellular Structure of a Normal and an Oblique Detonation
In order to compare the cellular structure of a normal and an oblique detonation, two-dimensional numerical simulations of a normal detonation were conducted. The computational domain is a rectangle of the size . The chemistry parameters and the degree of overdrive are the same as those in the simulations of the ODW. The exact ZND solution [26] is used as an initial condition, but all the velocities are shifted by (where is the dimensionless detonation speed). Due to this transformation, the detonation front remains quasistationary at its initial position (). Inflow and outflow boundary conditions are applied to the right and left side and slip boundary conditions are used in the direction. A transverse perturbation of the planar detonation is initiated by increasing the pressure by 15% inside the pocket . The adaptive grid has () at the coarsest level as well as four additional levels of refinement with a ratio of 2. As a result, this provides a resolution of 160 pts/hrl of a CJ detonation.
Figure 15 shows the numerical smoke foil of the normal detonation for the time interval (0, 30). The normal detonation propagates from left to right. As shown in Figure 15, the transverse waves are induced by the pressure pulse inside the pocket just behind the detonation front. At the earlier stage, the transverse waves consist of the up-running transverse waves (where ) and the down-running transverse waves (where ). The propagation directions of the two sets of transverse waves are changed after their collisions with the upper and lower boundaries, respectively. Later on, these two sets of transverse waves form the cellular structure after colliding with each other. From the figure, it is obtained that these two sets of transverse waves propagate with almost the same speed of 2.75. The length and width of a regular normal detonation cell are 3.7 and 2.2, and the ratio of the width and length is around 0.6.
In order to compare the characteristic cell size of a normal and an oblique detonation, a larger computational domain was employed to conduct another simulation of the oblique detonation. Other computational parameters and grid resolution are kept the same as those that are employed with a smaller computational domain. Figure 16 shows the local numerical smoke foil of the oblique detonation in the region of with 200 frames used to construct the figure. The time starts at , and the interval is 0.008. From this figure, it is obtained that the length and width of a characteristic oblique detonation cell are 3.5 and 2.1, and the ratio of the width and length is about 0.6. Therefore, it can be found that the characteristic cell size of a normal and an oblique detonation is almost the same.
4. Conclusions
In this study, two-dimensional simulations were carried out in order to investigate the formation and evolution of the cellular structure in the overdrive oblique detonations and quantitatively compare the characteristic cell size of an oblique and a normal detonation. The numerical results show that the LRTW starts at the transition zone because of the detonation instability. The pressure at the triple points is higher than the theoretical overdrive oblique detonation pressure. Because the heat release from the reaction zone is unable to maintain the high pressure at the triple points, the oblique detonation angle gradually decreases until the next triple point. To add on, such pressure could make the front of the detonation deform. The RRTW forms after the deformed detonation front interacts with the reflection waves from the wedge. The numerical smoke foil shows that the RRTW is enforced after the collision with the LRTW. The two sets of transverse waves propagate at almost the same relative velocity in the opposite direction along the detonation front even though both of them are convected downstream. After the quantitative comparison of the oblique and normal detonation, it is found that for the same degree of overdrive, the characteristic cell size of a normal and an oblique detonation is almost the same.
Data Availability
Data are available on request from the authors.
Conflicts of Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Acknowledgments
This work was supported by the General Program of National Natural Science Foundation of China (Nos. 51876182, 11972331, and 91941103) and the Fundamental Research Funds for the Central Universities (No. 20720180058).