- Email: [email protected]
A numerical study of wedge-induced detonations
A Numerical Study of Wedge-Induced Detonations MILTIADIS V. PAPALEXANDRIS
Graduate Aeronautical Laboratories, California Institute of Technology, Mail Code 301-46, Pasadena, CA 91125 The structure of detonations generated by the flow of a combustible mixture over a wedge is examined numerically. A recently designed, unsplit algorithm which integrates the convective and reaction source terms simultaneously, is employed for the simulations. The chemical kinetics is assumed to follow a simple one-step Arrhenius law. Two different geometrical configurations are considered. In the first one, the wedge is assumed to be long enough so that its top corner can not affect the structure of the reaction zone. In the second one, the wedge is short so that this corner can influence the reaction zone. It is observed that for moderate wedge angles the leading shock curves smoothly until it reaches a final angle. A stable detonation and a steady flow-field is established near the wedge. For higher wedge angles, an explosion occurs at the leading shock and the detonation becomes unstable. Grid-convergence studies for both stable and unstable detonations are presented and comparisons with earlier results are made. © 2000 by The Combustion Institute
INTRODUCTION The subject of the present article is the numerical study of oblique detonations. These detonations can be generated, for example, in supersonic flow of combustible mixtures past projectiles. The particular problem examined here is the flow past a wedge. The oblique shock that is formed at the tip of the wedge curves upward due to dilatation of the gas caused by the reaction that takes place behind the shock. Oblique detonations have attracted considerable attention because of potential applications to hypersonic propulsion. A very popular device that can be associated with oblique detonations is the ram accelerator; see, for example Brackett and Bogdanoff [1], Grismer and Powers [2], Lefebvre and Fujiwara [3], and Yungster et al. [4]. Another relevant concept that makes use of these waves for propulsion purposes is the oblique detonation wave engine (ODWE). The idea is to use the thrust from the wedge-induced detonation of a fuel–air mixture. The problem of premixing the fuel in a safe manner has hampered, so far, the construction of such an engine. There is an important question associated with this type of flows, and it has to do with the stability of the oblique wave. If the detonation behind the wave is stable, one should expect that the asymptotic structure of the shock is a Corresponding author. E-mail: [email protected] 0010-2180/00/$–see front matter PII S0010-2180(99)00113-3
standing oblique ZND wave. Nonetheless, this structure might not be realizable for certain wedge-angles and detonation parameters. For a detailed discussion of the possible configurations see, among others, Buckmaster [5] and Shepherd [6]. An asymptotic treatment of this problem in the hypersonic limit has been presented by Powers and Stewart [7]. Work on the numerical simulation of wedgeinduced detonations was done by various authors in the recent past [1– 4, 8 –15]. Results for standing detonations induced by long wedges (wedges that are too long for their top corner to affect the structure of the reaction zone) have been presented by Glenn and Pratt [8], Cambier et al. [9], Li et al. [10], Gerlinger et al. [11], and recently by Thaker and Chelliah [12]. Li et al. [10], in particular, provided numerical evidence for the stability of standing detonations even in the presence of very strong perturbations and they have included results for unsteady detonations. All the above works considered a multispecies reaction mechanism, either detailed or reduced. The recent paper of Thaker and Chelliah [12] provides comparisons between different mechanisms. Simulations on supersonic ramjet configurations were presented, among others, by Brackett and Bogdanoff [1], Grismer and Powers [2], Yungster et al. [4], and Yungster and Radhakrishnan [13]. These works, which also considered multispecies reaction mechanisms, provided parametric studies and numerical verCOMBUSTION AND FLAME 120:526 –538 (2000) © 2000 by The Combustion Institute Published by Elsevier Science Inc.
WEDGE-INDUCED DETONATIONS ification of the various propulsion concepts that rely on stable detonations. Other geometrical configurations have been studied too. Lefebvre and Fujiwara [3] and Fujiwara et al. [14] considered detonations induced by blunt bodies. More recently, Grismer [15] studied numerically a different wedge problem. His wedge was curved in such a way that the leading shock front was straight. The present article deals with two generic geometrical configurations: long wedges, and short wedges (wedges of relatively small length so that their top corner can possibly influence the structure of the reaction zone). Emphasis is given to wedge angles that result in unstable detonations. The main contribution of this work is the detailed numerical representation of many of the fine-length structures that such flows generate close to the wedge. The study of such structures is particularly interesting in the case of short wedges because the location of the top corner of the wedge plays a significant role in the far-field solution. Detailed numerical results that resolve fine structures generated by short wedges are not, as yet, available in the literature. A major difference from the works mentioned above is that in the present article a simplified, one-step reaction model has been assumed, instead of a multispecies model. This simplified model has been employed by many authors for the numerical study of planar detonations; see, among others, Oran et al. [16], Bourlioux and Majda [17], Cai [18], and Papalexandris et al. [19]. Despite its simplicity, it can reproduce many of the important structural features of detonations. This main reason for choosing this model is to save computing time. The other reason is to compare the flow-characteristics of wedge-induced detonations with those of planar detonations that have been studied numerically in the past. One should be aware, however, of the restrictions and limitations of this particular model. For example, the number of the unstable modes of the detonations predicted by this model always increases as the activation energy of the reaction increases, but this is not always the case in real-life detonations. For a discussion of the drawbacks of the model, see Clavin et al. [20]. The numerical simulations have been per-
527 formed on SGI-Onyx and DEC-Alpha workstations with an unsplit algorithm recently developed by Papalexandris et al. [19, 21]. The grid size used is at the level of 1200 ⫻ 250 cells (it varies from case to case).
FORMULATION OF THE PROBLEM Consider a simple model of chemical interaction of two calorically perfect gases, A 3 B, assuming one-step, irreversible, Arrhenius kinetics, and absence of dissipation mechanisms. The conservation equations of this system are given by: ⭸ ⫹ ⵜ 䡠 共 u兲 ⫽ 0, ⭸t
(1a)
1 ⭸u ⫹ u 䡠 ⵜu ⫹ ⵜp ⫽ 0, ⭸t
(1b)
⭸p ⫹ u 䡠 ⵜp ⫹ ␥ pⵜ 䡠 u ⫽ Kq 0共 ␥ ⫺ 1兲 ze ⫺Ea/T, ⭸t (1c) ⭸z ⫹ u 䡠 ⵜz ⫽ ⫺Kze ⫺Ea/T. ⭸t
(1d)
The equation of state of the reacting system satisfies the perfect-gas equation: T⫽
p ,
(2)
where the temperature, T, has been normalized by the gas constant. In the equations above, u ⫽ (u, v) is the velocity vector and z is the reactant mass fraction, satisfying 0 ⱕ z ⱕ 1. The parameters of the system are:
␥,
the specific heat ratio, assumed common for both species; q 0 , the heat-release parameter; E a, the activation-energy parameter; and K, an amplitude parameter that sets the spatial and temporal scales. For every detonation admitted by this system of equations, the shock velocity, D, has to be greater than D CJ, the velocity of the corresponding Chapman-Jouguet detonation. The parameter f, defined as
528
f⬅
M. V. PAPALEXANDRIS
冉 冊 D D CJ
2
,
is the overdrive factor of the detonation. The half-reaction length, L1/2, i.e., the distance between the shock wave and the point where z ⫽ 0.5, has been used as unit length throughout. The half-reaction length divided by the sound speed upstream, ahead of the shock, gives the unit time. The value of the amplitude parameter K is chosen in each detonation such that the half-reaction length, L 1/ 2 equals to unity. The evaluation of K requires the calculation of the 1-D, steady-state solution of (1), i.e., the ZND profile of the detonation, through numerical integration. As mentioned above, the algorithm that is used to solve these equations numerically is described in [19]. This algorithm is a generalization of the numerical scheme introduced by Lappas et al. [21] for the compressible Euler equations. The algorithm yields second-order accurate results in both space and time at smooth parts of the flow. It belongs to the family of MUSCL-type, shock-capturing schemes. It is an unsplit algorithm, in the sense that it integrates all the terms of the conservation laws simultaneously, in a fully coupled time-step, thus avoiding dimensional or time-splitting. Appropriate families of space–time manifolds are introduced, along with the conservation equations decouple to the characteristic equations of the corresponding 1-D homogeneous system. This system is then solved numerically along paths that lie on the manifolds. In all the tests performed the algorithm turned out to be stable even in the presence of very strong shocks. Therefore, there is no need to employ standard numerical techniques for stabilization of strong shocks, such as explicit artificial-viscosity terms or flux-splitting. Further, there is no need to supplement the algorithm with extra terms that are often used to avoid entropy-violating shocks. The extension of the algorithm to treat multicomponent reacting systems is straightforward.
DETONATIONS INDUCED BY LONG WEDGES In this section the wedge is assumed to be long enough so that its geometry downstream can
Fig. 1. Schematic of the computational domain for the problem of a detonation induced by a long wedge.
not influence the flow-field near the reaction zone. A schematic of the geometrical configuration and the computational domain is given in Fig. 1. The wedge is placed instantly in uniform flow of a reactive gas. A shock is immediately formed at the wedge. If the surrounding gas was inert, this shock would be oblique, at a prescribed angle with respect to the centerline of the wedge. Since the gas is assumed to be reactive, the shock is curved due to the dilatation of the reacting material behind the shock. A detonation might be established downstream if the shock temperature is sufficiently high. The wedge angle, , is an important parameter of this problem. It is expected that for small wedge angles the shock turns smoothly and the flow far downstream consists of a standing oblique ZND wave, i.e., a ZND wave with a nonzero transversal velocity component. Pratt et al. [23] have studied the parametric values of and M that result in such standing waves. For small wedge angles, both the inert and the equilibrium shock-polars admit solutions for the shock angle, . The shock near the tip is essentially inert and its angle can be computed from the inert shock-polar. The angle of the ZND wave, far downstream, can be computed from the equilibrium shock-polar. Given the state ahead of the shock, denoted by the subscript
WEDGE-INDUCED DETONATIONS
529
“1”, one can determine the flow variables behind the shock, denoted by the subscript “s”, by employing the standard kinematic relations for oblique shocks: F⬅
2 2 ⫺1兲2 ⫺2共␥ ⫹1兲 M1n q0/共cpT1兲 1⫹␥M1n ⫾ 冑共M1n , 2 共␥ ⫹1兲 M1n
(3a) tan共  ⫺ 兲 1 ⫽F⫽ , s tan
(3b)
p s ⫽ p 1 ⫹ u 21共1 ⫺ F兲 sin2 ,
(3c)
where M 1n⬘ is the normal Mach number ahead of the shock: M 1n ⬅
u 1 sin 
冑␥ T
.
(4)
For wedge angles larger than a certain value (but small enough so that the equilibrium shock-polar admits a solution for ), the shock can not turn smoothly, and a strong explosion is expected to occur on the front. This explosion is caused by the interaction of pressure waves inside the reaction zone. These pressure waves are emitted from the points near the wedge at which the material burns rapidly. Results for three different sets of detonation parameters are presented. The parameters are chosen so that they match the parameters of some of the planar detonations that were studied numerically in [17] and [19]. Inflow conditions have been assigned at the left boundary and at the first 7 cells of the bottom boundary. Reflecting conditions have been assigned at the rest of the bottom boundary. Finally outflow conditions (zero normal derivatives) have been assumed at the top and right boundaries. The flow at these boundaries is supersonic, and one must ensure that all information for the evaluation of the boundary fluxes comes from inside the computational domain. This is performed by copying the values from the boundary cells to their corresponding dummy cells, which also ensures zero normal derivatives at the boundary. The values of the specific heat ratio, ␥, and the heat release coefficient, q 0 , are set as follows:
␥ ⫽ 1.2,
q 0 ⫽ 50.
Upstream, the pressure and density of the gas are equal to unity. The thermodynamic variables of the system of equations have been made dimensionless with respect to this upstream state. Contour plots of the pressure, temperature, vorticity, and reactant mass fraction are presented. Each plot consists of 30 contours, equally distributed between the extremal values, except for the plots of the reactant mass fraction which contain 11 contours, at the levels z ⫽ 0.01, 0.1, 0.2, . . . , 0.9, 0.99. The CFL number for these simulations is set to 0.70. Case A This is a low activation energy case: E a ⫽ 10,
K ⫽ 3.1245.
As a first test, the wedge angle is set to ⫽ 20°. The upstream velocity of the gas is u 1 ⫽ 12.171. Given these parameters and upstream conditions, the theoretical prediction is that far downstream this detonation reduces to an oblique ZND wave of overdrive factor f ⫽ 1.2 and at an angle  ⫽ 34.02°. The equivalent 1-D detonation is linearly stable (see Erpenbeck [24], Fickett and Davis [25], and Lee and Stewart [26]). The equivalent 2-D planar detonation problem, however, has three unstable transversal modes (linear stability analysis of 2-D one-step detonations was initiated by Erpenbeck [24], and was later improved by Majda and Roytburd [27], Short and Stewart [28], Short [29], and others). Numerical simulations for this 2-D planar problem were performed, among others, in [17], [18], and [19]. Those simulations demonstrated that the cellular patterns of this detonation are formed in a periodic fashion, resulting in the formation of regular vortical structures in the wake. In the present study, the computational domain consists of 960 ⫻ 400 cells. This corresponds to a nominal resolution of 8 points per half-reaction length of the steady (ZND) solution. In this simulation the shock wave turns smoothly until it reaches the asymptotic angle, . Contour plots for this case are presented in Fig. 2. These plots are taken at t ⫽ 24.0. No change on the flow variables is observed at later
530
Fig. 2. Case A, ⫽ 20°: contour plots of the flow variables at t ⫽ 24.0.
times, which implies that the part of the flowfield that is covered by the computational domain reaches a steady state. In the area near the tip of the wedge, the shock is essentially inert. The pressure and particle velocity in this region are almost constant along a streamline. The temperature increase across the oblique shock is small and, consequently, the source term in the species equation (1d) remains very small. This source term, however, is small but not zero, so the material reacts, albeit very slowly. At the end of this slow-burning region the temperature has risen high enough to initiate and sustain fast
M. V. PAPALEXANDRIS burning of the gas. As a result of this fast burning close to the wedge, pressure waves are transmitted to the shock front. These waves interact with the slow-burning region behind the shock, causing the fluid particles to burn faster and dilate and the shock front to turn. Simulations were also carried for smaller wedge angles and the resulting flow-fields were found to have the same characteristics as in this case. Next, the wedge angle is increased to ⫽ 27°. This is the maximum angle for which the equilibrium shock-polar admits a solution for . For a detonation of overdrive factor f ⫽ 1.2 to occur, the upstream velocity has to be u 1 ⫽ 9.255. The asymptotic limit of the shock angle is readily found to be  ⫽ 53.7°. The computational domain for this simulation consists of 1020 ⫻ 402 cells, with a nominal resolution of 6 points per half-reaction length. Contour plots of the variables for this problem, at t ⫽ 50.0, are given in Fig. 3. Several triple points have been formed along the leading front. These points can be observed, for example, in the vorticity plot of Fig. 3, since the contact discontinuities that are attached to the leading front emanate from the triple points. All the triple-point structures are convected by the flow but have different velocities. Therefore, collisions between the triple points are expected to occur after enough time, most likely outside the area covered by the computational mesh, thus forming the cellular patterns that are typically encountered in detonations. According to the linear stability studies mentioned above, formation of triple points far downstream should have been observed in the previous simulation (when ⫽ 20°), had the computational domain been large enough. Case B The activation energy and the stiffness coefficients are now set at: E a ⫽ 50,
K ⫽ 99.762.
The wedge angle is ⫽ 20°. The upstream velocity is u 1 ⫽ 20.58. Theoretically, the shock angle tends asymptotically to a shock angle  ⫽ 27.9°. At this limit the detonation is an oblique ZND wave of overdrive factor f ⫽ 2.0.
WEDGE-INDUCED DETONATIONS
Fig. 3. Case A, ⫽ 27°: contour plots of the flow variables at t ⫽ 50.0.
Three different grid sizes have been used for this problem. In the first test the computational domain consists of 285 ⫻ 60 cells (Fig. 4a). In the second test it consists of 570 ⫻ 120 cells (Fig. 4b). In the final test it consists of 1140 ⫻ 240 cells (Fig. 4c), corresponding to a nominal resolution of 6 points per half-reaction length of the one-dimensional, steady-state solution. The shock turns smoothly until it reaches a steady state, very close to the asymptotic limit. The results presented in Fig. 4 are taken at t ⫽ 18.0. By this time, the solution has already reached the steady state. Planar detonations with these parameters, however, are unstable. It is therefore expected that formation of triple points will eventually occur further downstream, outside the computational domain. It can be
531 observed that the same steady-state solution has been computed with all three meshes. As expected, the shock profiles on the coarse meshes are more smeared than the ones on the fine mesh. Next, the wedge angle is increased to ⫽ 35°, which is near the maximum angle for which the equilibrium shock-polar admits a solution. As before, three different mesh sizes have been used for this problem, consisting of 560 ⫻ 64 cells, 560 ⫻ 128 cells, and 1140 ⫻ 240 cells, respectively. The finest computational domain corresponds to a nominal resolution of 8 points per half-reaction length. The upstream velocity is u 1 ⫽ 11.509. The asymptotic limit is a detonation with overdrive factor f ⫽ 1.2 and at an angle  ⫽ 56.8°. Results for this simulation, taken at t ⫽ 50.0, are given in Figs. 5a, 5b, and 5c. The flow can not turn smoothly in this case, because of the high value of the wedge angle. As a result, a strong explosion takes place at the front. The center of this explosion is a triple point. The incident shock and the Mach stem are the two parts of the main front, below and above the triple point respectively. Another shock emanates from the triple point, the reflected shock, which hits the wedge and reflects back. Additionally, a contact discontinuity (shear layer) is formed between the Mach stem and the reflected shock. The material behind the incident shock burns very slowly because the shock is relatively weak and the temperature increase is small. The material behind the Mach stem, however, burns fast because of the high temperature rise. Consequently, there is a strong density and temperature gradient across the shear layer. The shear layer becomes unstable quickly and generates strong vortical structures that are convected downstream (Fig. 5c). Diffusion in these simulations is introduced by the discretization of the equations and the truncation error of the scheme. Consequently, the size and the velocity of the vortical structures are determined by the size of the grid and the implicit artificial viscosity of the algorithm. The convective Mach number for the shear layer is between 0.5 ⬃ 0.7, therefore the layer can be considered slightly compressible. It is worth mentioning that the flow between the incident shock and the wedge is steady. Furthermore, the region between the shear
532
M. V. PAPALEXANDRIS
Fig. 4. (a) Case B, ⫽ 20°: contour plots of the flow variables at t ⫽ 18.0. Grid size, 285 ⫻ 60 cells. (b) Case B, ⫽ 20°: contour plots of the flow variables at t ⫽ 18.0. Grid size, 570 ⫻ 120 cells. (c) Case B, ⫽ 20°: contour plots of the flow variables at t ⫽ 18.0. Grid size, 1140 ⫻ 240 cells.
layer and the Mach stem is subsonic in the vicinity of the triple point. The boundary of this pocket is a sonic line which emanates from the triple point and lies in the shear layer. After some distance the sonic line curves upward and ends up in the Mach stem. Beyond this subsonic pocket, the flow has all the typical characteristics of planar detonations, such as formation of colliding triple points and transverse waves, and
vorticity generation at the main front. Qualitatively, the results on the three different meshes are the same. Interestingly, the location of the explosion at the front is the same regardless of mesh size. The results on the coarse meshes, however, are more diffusive than than the ones on the fine mesh: the shock profiles are more smeared, and the structures of the shear layer and the reaction zone are more diffuse.
WEDGE-INDUCED DETONATIONS
533
Fig. 5. (a) Case B, ⫽ 35°: contour plots of the flow variables at t ⫽ 50.0. Grid size, 280 ⫻ 64 cells. (b) Case B, ⫽ 35°: contour plots of the flow variables at t ⫽ 50.0. Grid size, 560 ⫻ 128 cells. (c) Case B, ⫽ 35°: contour plots of the flow variables at t ⫽ 50.0. Grid size, 1120 ⫻ 256 cells.
Since the equilibrium shock-polar admits a solution for the shock angle , it is expected that far downstream the Mach stem reduces to an oblique ZND wave, at an angle  ⫽ 56.8°. Linear stability analysis, however, reveals that two-dimensional detonations are intrinsically unstable. Consequently, the instability mechanisms that are encountered in detonations propagating in channels and lead to the formation of triple points (see, for example, [16 –19]) should also appear on the flow-field of the oblique
detonations under study. This region, however, is too far downstream to be included in the computational domain. For higher wedge angles, the equilibrium shock polar can not give a solution for . In such situations, a strong triple point is also formed at the main front. But since there is no solution to the equilibrium shock-polar, the Mach stem can not reduce to an oblique ZND wave downstream, and it is everywhere curved. For even higher wedge angles, none of the shock-polars
534
M. V. PAPALEXANDRIS
Fig. 6. Case C, ⫽ 20°: contour plots of the flow variables at t ⫽ 18.0.
admits a solution, and the main front detaches from the wedge. Such high wedge angles have not been considered in the present work. Case C The activation energy is the same as in case B, but the stiffness coefficient has now been increased: E a ⫽ 50,
K ⫽ 230.75.
The wedge angle is ⫽ 20°. The upstream velocity is u 1 ⫽ 18.051. The theoretical prediction is that far downstream the detonation will be a ZND wave of overdrive factor f ⫽ 1.60 and at an angle  ⫽ 28.5°. The computational domain consists of 1200 ⫻ 200 cells, corresponding to a nominal resolution of 8 points per half-reaction length of the one-dimensional, steady-state solution. As in case B, the shock turns smoothly until it reaches a steady state. Contour plots of the flow variables are given in Fig. 6. These results are taken at time t ⫽ 18.0. No change in the flow variables could be observed after that time. The shock angle at the right boundary is very close to the asymptotic limit  ⫽ 28.5°. The case of a higher wedge angle, namely ⫽ 30°, has also been considered. The computational domain consists of 1440 ⫻ 240 cells, corresponding to a nominal resolution of 8 points per half-reaction length. The upstream velocity is now set at u 1 ⫽ 12.035. Under this
Fig. 7. Case C, ⫽ 35°: contour plots of the flow variables at t ⫽ 36.0.
initial condition, the theoretical prediction is a detonation of overdrive factor f ⫽ 1.6, at a shock angle  ⫽ 45.7°. Results for this simulation, taken at time t ⫽ 36.0, are given in Fig. 7. These results are similar to the ones obtained in the high-angle simulation of case B. The basic features that were encountered in the results of case B (such as the explosion at the front, the development of an unsteady shear layer, and the formation of a subsonic pocket behind the Mach stem in the vicinity of the triple point) can also be observed in this simulation. It is worth mentioning that Li et al. [10] also presented simulations for wedge-induced detonations using the FCT algorithm. Their numerical results agree qualitatively with the results obtained with the proposed unsplit algorithm. In particular, they also observed that for small wedge angles the main front turns smoothly, whereas an explosion occurs at the front if large wedge angles are considered. Quantitative comparisons between the two studies can not be made because Li et al. considered a multispecies combustion model and did not use dimensionless quantities for their simulations. DETONATIONS INDUCED BY SHORT WEDGES In the cases examined so far, it was assumed that the wedge is long enough for the explosion
WEDGE-INDUCED DETONATIONS
535
Fig. 8. Flow past a wedge of weight of height h w ⫽ 70 and angle ⫽ 35°. Contour plots of the flow variables at t ⫽ 46.0.
to occur upstream of the corner of the wedge. In such cases, the location of the triple point and the flow-field in that neighborhood are determined completely by the kinetics of the reaction. If, however, the wedge is not long enough, then the explosion will take place near the corner. The effect of the corner in such situations has also been studied numerically, and the results are described below. The chemical reaction increases the temperature of the fluid, while the expansion at the corner decreases it. These two mechanisms “compete” against each other. Furthermore, by including the corner in the computational domain, a second characteristic length is introduced. This is the height of the wedge, h w , the first one being the half-reaction length. It is the combination of these two length-scales with the wedge angle that ultimately determines which of the two mechanisms will dominate. In the numerical simulations, the activation
energy and the stiffness coefficient are set at E a ⫽ 50 and K ⫽ 99.762, respectively. The wedge angle is ⫽ 35°. The upstream state of the fluid is p 1 ⫽ 1.0,
1 ⫽ 1.0,
u 1 ⫽ 11.509.
The flow that is produced by these parameters and upstream conditions in the case of a long wedge was examined in the previous section (see Fig. 5c). Three different wedge heights are considered. First, the wedge height is set at h w ⫽ 70.0. The computational domain consists of 930 ⫻ 330 cells. The length of the domain is 155.0, and its width from the upper wall of the wedge is 55.0. Results for this simulation, taken at t ⫽ 46.0, are shown in Fig. 8. In this case the explosion occurs upstream with respect to the corner. The expansion at the corner does not affect the explosion because the flow in the
536
M. V. PAPALEXANDRIS
Fig. 9. Flow past a wedge of weight of height h w ⫽ 50 and angle ⫽ 35°. Contour plots of the flow variables at t ⫽ 46.0.
corner is supersonic. The leading shock is expected to reduce far downstream to a Chapman-Jouguet (CJ) wave, just like the one-dimensional detonation initiated by a moving piston that comes suddenly to a rest. The flow-field in the neighborhood of the inert shock and the explosion is the same as in the corresponding case with a long wedge (see Fig. 5c). The interaction of the reflected shock with the shear layer causes the fluid to decelerate, thus forming another subsonic region in the flow-field, besides the one in the vicinity of the triple point. It is also worth mentioning that the expansion that takes place at the corner affects the evolution of the shear layer. As the fluid below the shear layer expands, the pressure and density drop, generating even higher entropy gradients across the layer. Therefore, the expansion at the corner leads to an increased amount of vorticity generation across the shear layer.
Subsequently, the wedge height is reduced to h w ⫽ 50.0. The computational domain consists of 810 ⫻ 420 cells. The length of the domain is 135.0, and its width, as measured from the upper wall of the wedge, is 70.0. Results for this simulation are shown in Fig. 9. They are taken at time t ⫽ 46.0. In this case, the explosion takes place downstream of the corner. The flow in the subsonic area in the vicinity of the triple point (between the Mach stem and the shear layer) is influenced by the expansion at the corner. Additionally, the expansion at the corner affects the curvature of the leading front and the reflected shock. A fluid element moving parallel to the wedge does not have time to increase its temperature substantially, because of the small length of the wedge. Consequently, it remains almost unreacted when it reaches the head of the expansion. The expansion produces a further decrease of
WEDGE-INDUCED DETONATIONS
537
Fig. 10. Flow past a wedge of weight of height h w ⫽ 25.0 and angle ⫽ 35°. Contour plots of the flow variables at t ⫽ 10.0.
the temperature which delays the initiation of the reaction even more. It is observed that the material on the upper wall of the wedge has remained only partially reacted (about 20%), as a result of the temperature decrease caused by the expansion at the corner. As a final test, the wedge height is lowered further to h w ⫽ 25.0. The simulation is performed on a domain of 1266 ⫻ 270 cells. The length of the domain is 210.0 and its width is 45.0. Results from this simulation, taken at t ⫽ 10.0, are presented in Fig. 10. It is observed that the expansion at the corner reduces the temperature so much that the gas near the wedge remains almost unreacted because the time needed for the temperature to increase high enough to establish rapid reaction becomes very large. As a result, a detonation can not be established and the shock wave is expected to reduce to a Mach wave downstream. It is also observed that there is a small pocket of slightly reacted material near the wedge. It can be verified by looking at results taken at early times that the formation of this pocket is a transient phenomenon caused by the interaction of the shock wave and the expansion at the corner, in the beginning of the simulation. It is convected downstream with the fluid velocity. The minimum value of the reactant mass fraction inside the pocket is z ⯝ 0.9. If the material reacted completely, pressure waves would be transmitted to the shock wave and a CJ detonation could be established. But it is observed that the reaction rate, ˙z ⫽ ⫺Kz exp (⫺E a/T), inside the pocket decays with time, which suggests that
the reaction process will not be completed inside the pocket. After some time the pocket exits the computational domain, and the shock front eventually assumes a fixed position. No change in the flow variables can be observed after that. The material all along the wedge will remain only partially reacted, and a detonation will not be established. CONCLUDING REMARKS A numerical study of wedge-induced detonations has been presented. The numerical scheme that was used in the simulations integrates all the terms of the governing equations simultaneously, without the use of flux-splitting or additional artificial viscosity. It was confirmed that, for small wedge angles, the shock that is attached to the wedge turns smoothly to an oblique ZND wave. For high wedge angles, however, such a smooth turn is not possible and an explosion takes place at the front. The center of the explosion is a triple point that assumes a fixed position after some time. A shear layer emanates from the triple point that eventually becomes unstable. The effect of the corner of the wedge was also studied numerically. It was shown that when the explosion of the leading front takes place upstream of the corner, the expansion at the corner does not affect the evolution of the front, which reduces to a CJ wave. When the explosion occurs downstream of the corner, the curvature of the front and the reaction process depend on
538
M. V. PAPALEXANDRIS
the expansion at the corner. It appears that for wedge heights small enough, a detonation can not be established downstream, and the front decays downstream to a Mach wave. Future directions include the numerical simulation of multicomponent reacting systems. This will help the study of various problems arising in supersonic combustion, such as the effect of chain-branching. The better understanding of these mechanisms can be helpful in the design of the next-generation, high-speed propulsion systems. The author would like to thank Prof. Leonard for his support and many fruitful discussions. Thanks are also due to Prof. Dimotakis and Prof. Shepherd for their suggestions. This work is sponsored by the Air Force Office of Scientific Research Grant Nos. F49620-94-1-0353 and F49626-93-1-00338, whose support is gratefully acknowledged. REFERENCES 1. 2. 3. 4. 5. 6.
7. 8. 9. 10.
Brackett, D. C., Bogdanoff, D. W., J. Propul. Power 5:276 –281 (1989). Grismer, M. J., Powers, J. M., J. Propul. Power 11:105– 111 (1995). Lefebvre, M. H. and Fujiwara, T., Combust. Flame 100:85–93 (1995). Yungster, S., Radhakrishnan, K., Rabinowitz, M. J. (1995). AIAA paper 95–2489. Buckmaster, J. D., Combust. Sci. Technol. 72:283–296 (1990). Shepherd, J. E., in Combustion in High-Speed Flows (J. Buckmaster et al., Eds.), Kluwer Academic Publishers, 1994, pp. 373– 420. Powers, J. M., Stewart, D. S., AIAA J. 30:726 –736 (1992). Glenn, D. E., Pratt, D. T. (1998). AIAA paper 88-0440. Cambier, J. L., Adelman, H., Menees, G. P., J. Propul. Power 5:482– 491 (1989). Li, C., Kailasanath, K., Oran, E. S., Phys. Fluids 6:1600 –1611 (1994).
11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33.
Gerlinger, P., Algermissen, J., Bruggemann, D., AIAA J. 33:1865–1870 (1995). Thaker, A. A., Chelliah, H. K., Combust. Theory Modelling 1:347–376 (1997). Yungster, S., and Radhakrishnan, K. (1995). SA Report 95-19. Fujiwara, T., Matsuo, A., Nomoto, H. (1988). AIAA paper 88-0098. Grismer, M. J. (1994). Ph.D. thesis, University of Notre Dame. Oran, E. S., Kailasanath, K., Guirguis, R. H., Prog. Aeronaut. Astronaut. 114:155 (1988). Bourlioux, A., Majda, A. J., Combust. Flame 90:211– 229 (1992). Cai, W., AIAA J. 33:1248 –1255 (1995). Papalexandris, M. V., Leonard, A., Dimotakis, P. E., GALCIT Report 97-14, submitted to J. Computa. Phys. Clavin, P., He, L., Williams, F. A., Phys. Fluids 9:3764 – 3785 (1997). Lappas, T., Leonard, A., Dimotakis, P. E., SIAM J. Sci. Comp. 20:1481–1512 (1999). Papalexandris, M. V., Leonard, A., Dimotakis, P. E., J. Comp. Phys. 134:31– 61 (1997). Pratt, D. T., Humphrey, J. W., Glenn, D. E., J. Propul. Power 7:837– 845 (1991). Erpenbeck, J. J., Phys. Fluids 7:684 – 696 (1964). Fickett, W., Davis, W. C., Detonation, University of California Berkeley Press, Berkeley, CA, 1979. Lee, H. I., Stewart, D. S., J. Fluid Mech. 206:103–132 (1990). Majda, A., Roytburd, V., Studies Appl. Math. 87:135– 174 (1992). Short, M., Stewart, D. S., J. Fluid Mech. 340:249 –295 (1997). Short, M., SIAM J. Appl. Math. 57:307–326 (1997). Klein, R., Stewart D. S., SIAM J. Appl. Math. 53:1401– 1435 (1993). Lasseigne, D. G., Hussaini, M. Y., Phys. Fluids A 5:1047–1058 (1993). Menikoff, R., Lackner, K. S., Bukiet, B. G., Combust. Flame 104:219 –240 (1996). Stewart, D. S., Bdzil, J. B., Combust. Flame 72:311–323 (1988).
Received 26 January 1998; revised 20 May 1999; accepted 21 July 1999
Graduate Aeronautical Laboratories, California Institute of Technology, Mail Code 301-46, Pasadena, CA 91125 The structure of detonations generated by the flow of a combustible mixture over a wedge is examined numerically. A recently designed, unsplit algorithm which integrates the convective and reaction source terms simultaneously, is employed for the simulations. The chemical kinetics is assumed to follow a simple one-step Arrhenius law. Two different geometrical configurations are considered. In the first one, the wedge is assumed to be long enough so that its top corner can not affect the structure of the reaction zone. In the second one, the wedge is short so that this corner can influence the reaction zone. It is observed that for moderate wedge angles the leading shock curves smoothly until it reaches a final angle. A stable detonation and a steady flow-field is established near the wedge. For higher wedge angles, an explosion occurs at the leading shock and the detonation becomes unstable. Grid-convergence studies for both stable and unstable detonations are presented and comparisons with earlier results are made. © 2000 by The Combustion Institute
INTRODUCTION The subject of the present article is the numerical study of oblique detonations. These detonations can be generated, for example, in supersonic flow of combustible mixtures past projectiles. The particular problem examined here is the flow past a wedge. The oblique shock that is formed at the tip of the wedge curves upward due to dilatation of the gas caused by the reaction that takes place behind the shock. Oblique detonations have attracted considerable attention because of potential applications to hypersonic propulsion. A very popular device that can be associated with oblique detonations is the ram accelerator; see, for example Brackett and Bogdanoff [1], Grismer and Powers [2], Lefebvre and Fujiwara [3], and Yungster et al. [4]. Another relevant concept that makes use of these waves for propulsion purposes is the oblique detonation wave engine (ODWE). The idea is to use the thrust from the wedge-induced detonation of a fuel–air mixture. The problem of premixing the fuel in a safe manner has hampered, so far, the construction of such an engine. There is an important question associated with this type of flows, and it has to do with the stability of the oblique wave. If the detonation behind the wave is stable, one should expect that the asymptotic structure of the shock is a Corresponding author. E-mail: [email protected] 0010-2180/00/$–see front matter PII S0010-2180(99)00113-3
standing oblique ZND wave. Nonetheless, this structure might not be realizable for certain wedge-angles and detonation parameters. For a detailed discussion of the possible configurations see, among others, Buckmaster [5] and Shepherd [6]. An asymptotic treatment of this problem in the hypersonic limit has been presented by Powers and Stewart [7]. Work on the numerical simulation of wedgeinduced detonations was done by various authors in the recent past [1– 4, 8 –15]. Results for standing detonations induced by long wedges (wedges that are too long for their top corner to affect the structure of the reaction zone) have been presented by Glenn and Pratt [8], Cambier et al. [9], Li et al. [10], Gerlinger et al. [11], and recently by Thaker and Chelliah [12]. Li et al. [10], in particular, provided numerical evidence for the stability of standing detonations even in the presence of very strong perturbations and they have included results for unsteady detonations. All the above works considered a multispecies reaction mechanism, either detailed or reduced. The recent paper of Thaker and Chelliah [12] provides comparisons between different mechanisms. Simulations on supersonic ramjet configurations were presented, among others, by Brackett and Bogdanoff [1], Grismer and Powers [2], Yungster et al. [4], and Yungster and Radhakrishnan [13]. These works, which also considered multispecies reaction mechanisms, provided parametric studies and numerical verCOMBUSTION AND FLAME 120:526 –538 (2000) © 2000 by The Combustion Institute Published by Elsevier Science Inc.
WEDGE-INDUCED DETONATIONS ification of the various propulsion concepts that rely on stable detonations. Other geometrical configurations have been studied too. Lefebvre and Fujiwara [3] and Fujiwara et al. [14] considered detonations induced by blunt bodies. More recently, Grismer [15] studied numerically a different wedge problem. His wedge was curved in such a way that the leading shock front was straight. The present article deals with two generic geometrical configurations: long wedges, and short wedges (wedges of relatively small length so that their top corner can possibly influence the structure of the reaction zone). Emphasis is given to wedge angles that result in unstable detonations. The main contribution of this work is the detailed numerical representation of many of the fine-length structures that such flows generate close to the wedge. The study of such structures is particularly interesting in the case of short wedges because the location of the top corner of the wedge plays a significant role in the far-field solution. Detailed numerical results that resolve fine structures generated by short wedges are not, as yet, available in the literature. A major difference from the works mentioned above is that in the present article a simplified, one-step reaction model has been assumed, instead of a multispecies model. This simplified model has been employed by many authors for the numerical study of planar detonations; see, among others, Oran et al. [16], Bourlioux and Majda [17], Cai [18], and Papalexandris et al. [19]. Despite its simplicity, it can reproduce many of the important structural features of detonations. This main reason for choosing this model is to save computing time. The other reason is to compare the flow-characteristics of wedge-induced detonations with those of planar detonations that have been studied numerically in the past. One should be aware, however, of the restrictions and limitations of this particular model. For example, the number of the unstable modes of the detonations predicted by this model always increases as the activation energy of the reaction increases, but this is not always the case in real-life detonations. For a discussion of the drawbacks of the model, see Clavin et al. [20]. The numerical simulations have been per-
527 formed on SGI-Onyx and DEC-Alpha workstations with an unsplit algorithm recently developed by Papalexandris et al. [19, 21]. The grid size used is at the level of 1200 ⫻ 250 cells (it varies from case to case).
FORMULATION OF THE PROBLEM Consider a simple model of chemical interaction of two calorically perfect gases, A 3 B, assuming one-step, irreversible, Arrhenius kinetics, and absence of dissipation mechanisms. The conservation equations of this system are given by: ⭸ ⫹ ⵜ 䡠 共 u兲 ⫽ 0, ⭸t
(1a)
1 ⭸u ⫹ u 䡠 ⵜu ⫹ ⵜp ⫽ 0, ⭸t
(1b)
⭸p ⫹ u 䡠 ⵜp ⫹ ␥ pⵜ 䡠 u ⫽ Kq 0共 ␥ ⫺ 1兲 ze ⫺Ea/T, ⭸t (1c) ⭸z ⫹ u 䡠 ⵜz ⫽ ⫺Kze ⫺Ea/T. ⭸t
(1d)
The equation of state of the reacting system satisfies the perfect-gas equation: T⫽
p ,
(2)
where the temperature, T, has been normalized by the gas constant. In the equations above, u ⫽ (u, v) is the velocity vector and z is the reactant mass fraction, satisfying 0 ⱕ z ⱕ 1. The parameters of the system are:
␥,
the specific heat ratio, assumed common for both species; q 0 , the heat-release parameter; E a, the activation-energy parameter; and K, an amplitude parameter that sets the spatial and temporal scales. For every detonation admitted by this system of equations, the shock velocity, D, has to be greater than D CJ, the velocity of the corresponding Chapman-Jouguet detonation. The parameter f, defined as
528
f⬅
M. V. PAPALEXANDRIS
冉 冊 D D CJ
2
,
is the overdrive factor of the detonation. The half-reaction length, L1/2, i.e., the distance between the shock wave and the point where z ⫽ 0.5, has been used as unit length throughout. The half-reaction length divided by the sound speed upstream, ahead of the shock, gives the unit time. The value of the amplitude parameter K is chosen in each detonation such that the half-reaction length, L 1/ 2 equals to unity. The evaluation of K requires the calculation of the 1-D, steady-state solution of (1), i.e., the ZND profile of the detonation, through numerical integration. As mentioned above, the algorithm that is used to solve these equations numerically is described in [19]. This algorithm is a generalization of the numerical scheme introduced by Lappas et al. [21] for the compressible Euler equations. The algorithm yields second-order accurate results in both space and time at smooth parts of the flow. It belongs to the family of MUSCL-type, shock-capturing schemes. It is an unsplit algorithm, in the sense that it integrates all the terms of the conservation laws simultaneously, in a fully coupled time-step, thus avoiding dimensional or time-splitting. Appropriate families of space–time manifolds are introduced, along with the conservation equations decouple to the characteristic equations of the corresponding 1-D homogeneous system. This system is then solved numerically along paths that lie on the manifolds. In all the tests performed the algorithm turned out to be stable even in the presence of very strong shocks. Therefore, there is no need to employ standard numerical techniques for stabilization of strong shocks, such as explicit artificial-viscosity terms or flux-splitting. Further, there is no need to supplement the algorithm with extra terms that are often used to avoid entropy-violating shocks. The extension of the algorithm to treat multicomponent reacting systems is straightforward.
DETONATIONS INDUCED BY LONG WEDGES In this section the wedge is assumed to be long enough so that its geometry downstream can
Fig. 1. Schematic of the computational domain for the problem of a detonation induced by a long wedge.
not influence the flow-field near the reaction zone. A schematic of the geometrical configuration and the computational domain is given in Fig. 1. The wedge is placed instantly in uniform flow of a reactive gas. A shock is immediately formed at the wedge. If the surrounding gas was inert, this shock would be oblique, at a prescribed angle with respect to the centerline of the wedge. Since the gas is assumed to be reactive, the shock is curved due to the dilatation of the reacting material behind the shock. A detonation might be established downstream if the shock temperature is sufficiently high. The wedge angle, , is an important parameter of this problem. It is expected that for small wedge angles the shock turns smoothly and the flow far downstream consists of a standing oblique ZND wave, i.e., a ZND wave with a nonzero transversal velocity component. Pratt et al. [23] have studied the parametric values of and M that result in such standing waves. For small wedge angles, both the inert and the equilibrium shock-polars admit solutions for the shock angle, . The shock near the tip is essentially inert and its angle can be computed from the inert shock-polar. The angle of the ZND wave, far downstream, can be computed from the equilibrium shock-polar. Given the state ahead of the shock, denoted by the subscript
WEDGE-INDUCED DETONATIONS
529
“1”, one can determine the flow variables behind the shock, denoted by the subscript “s”, by employing the standard kinematic relations for oblique shocks: F⬅
2 2 ⫺1兲2 ⫺2共␥ ⫹1兲 M1n q0/共cpT1兲 1⫹␥M1n ⫾ 冑共M1n , 2 共␥ ⫹1兲 M1n
(3a) tan共  ⫺ 兲 1 ⫽F⫽ , s tan
(3b)
p s ⫽ p 1 ⫹ u 21共1 ⫺ F兲 sin2 ,
(3c)
where M 1n⬘ is the normal Mach number ahead of the shock: M 1n ⬅
u 1 sin 
冑␥ T
.
(4)
For wedge angles larger than a certain value (but small enough so that the equilibrium shock-polar admits a solution for ), the shock can not turn smoothly, and a strong explosion is expected to occur on the front. This explosion is caused by the interaction of pressure waves inside the reaction zone. These pressure waves are emitted from the points near the wedge at which the material burns rapidly. Results for three different sets of detonation parameters are presented. The parameters are chosen so that they match the parameters of some of the planar detonations that were studied numerically in [17] and [19]. Inflow conditions have been assigned at the left boundary and at the first 7 cells of the bottom boundary. Reflecting conditions have been assigned at the rest of the bottom boundary. Finally outflow conditions (zero normal derivatives) have been assumed at the top and right boundaries. The flow at these boundaries is supersonic, and one must ensure that all information for the evaluation of the boundary fluxes comes from inside the computational domain. This is performed by copying the values from the boundary cells to their corresponding dummy cells, which also ensures zero normal derivatives at the boundary. The values of the specific heat ratio, ␥, and the heat release coefficient, q 0 , are set as follows:
␥ ⫽ 1.2,
q 0 ⫽ 50.
Upstream, the pressure and density of the gas are equal to unity. The thermodynamic variables of the system of equations have been made dimensionless with respect to this upstream state. Contour plots of the pressure, temperature, vorticity, and reactant mass fraction are presented. Each plot consists of 30 contours, equally distributed between the extremal values, except for the plots of the reactant mass fraction which contain 11 contours, at the levels z ⫽ 0.01, 0.1, 0.2, . . . , 0.9, 0.99. The CFL number for these simulations is set to 0.70. Case A This is a low activation energy case: E a ⫽ 10,
K ⫽ 3.1245.
As a first test, the wedge angle is set to ⫽ 20°. The upstream velocity of the gas is u 1 ⫽ 12.171. Given these parameters and upstream conditions, the theoretical prediction is that far downstream this detonation reduces to an oblique ZND wave of overdrive factor f ⫽ 1.2 and at an angle  ⫽ 34.02°. The equivalent 1-D detonation is linearly stable (see Erpenbeck [24], Fickett and Davis [25], and Lee and Stewart [26]). The equivalent 2-D planar detonation problem, however, has three unstable transversal modes (linear stability analysis of 2-D one-step detonations was initiated by Erpenbeck [24], and was later improved by Majda and Roytburd [27], Short and Stewart [28], Short [29], and others). Numerical simulations for this 2-D planar problem were performed, among others, in [17], [18], and [19]. Those simulations demonstrated that the cellular patterns of this detonation are formed in a periodic fashion, resulting in the formation of regular vortical structures in the wake. In the present study, the computational domain consists of 960 ⫻ 400 cells. This corresponds to a nominal resolution of 8 points per half-reaction length of the steady (ZND) solution. In this simulation the shock wave turns smoothly until it reaches the asymptotic angle, . Contour plots for this case are presented in Fig. 2. These plots are taken at t ⫽ 24.0. No change on the flow variables is observed at later
530
Fig. 2. Case A, ⫽ 20°: contour plots of the flow variables at t ⫽ 24.0.
times, which implies that the part of the flowfield that is covered by the computational domain reaches a steady state. In the area near the tip of the wedge, the shock is essentially inert. The pressure and particle velocity in this region are almost constant along a streamline. The temperature increase across the oblique shock is small and, consequently, the source term in the species equation (1d) remains very small. This source term, however, is small but not zero, so the material reacts, albeit very slowly. At the end of this slow-burning region the temperature has risen high enough to initiate and sustain fast
M. V. PAPALEXANDRIS burning of the gas. As a result of this fast burning close to the wedge, pressure waves are transmitted to the shock front. These waves interact with the slow-burning region behind the shock, causing the fluid particles to burn faster and dilate and the shock front to turn. Simulations were also carried for smaller wedge angles and the resulting flow-fields were found to have the same characteristics as in this case. Next, the wedge angle is increased to ⫽ 27°. This is the maximum angle for which the equilibrium shock-polar admits a solution for . For a detonation of overdrive factor f ⫽ 1.2 to occur, the upstream velocity has to be u 1 ⫽ 9.255. The asymptotic limit of the shock angle is readily found to be  ⫽ 53.7°. The computational domain for this simulation consists of 1020 ⫻ 402 cells, with a nominal resolution of 6 points per half-reaction length. Contour plots of the variables for this problem, at t ⫽ 50.0, are given in Fig. 3. Several triple points have been formed along the leading front. These points can be observed, for example, in the vorticity plot of Fig. 3, since the contact discontinuities that are attached to the leading front emanate from the triple points. All the triple-point structures are convected by the flow but have different velocities. Therefore, collisions between the triple points are expected to occur after enough time, most likely outside the area covered by the computational mesh, thus forming the cellular patterns that are typically encountered in detonations. According to the linear stability studies mentioned above, formation of triple points far downstream should have been observed in the previous simulation (when ⫽ 20°), had the computational domain been large enough. Case B The activation energy and the stiffness coefficients are now set at: E a ⫽ 50,
K ⫽ 99.762.
The wedge angle is ⫽ 20°. The upstream velocity is u 1 ⫽ 20.58. Theoretically, the shock angle tends asymptotically to a shock angle  ⫽ 27.9°. At this limit the detonation is an oblique ZND wave of overdrive factor f ⫽ 2.0.
WEDGE-INDUCED DETONATIONS
Fig. 3. Case A, ⫽ 27°: contour plots of the flow variables at t ⫽ 50.0.
Three different grid sizes have been used for this problem. In the first test the computational domain consists of 285 ⫻ 60 cells (Fig. 4a). In the second test it consists of 570 ⫻ 120 cells (Fig. 4b). In the final test it consists of 1140 ⫻ 240 cells (Fig. 4c), corresponding to a nominal resolution of 6 points per half-reaction length of the one-dimensional, steady-state solution. The shock turns smoothly until it reaches a steady state, very close to the asymptotic limit. The results presented in Fig. 4 are taken at t ⫽ 18.0. By this time, the solution has already reached the steady state. Planar detonations with these parameters, however, are unstable. It is therefore expected that formation of triple points will eventually occur further downstream, outside the computational domain. It can be
531 observed that the same steady-state solution has been computed with all three meshes. As expected, the shock profiles on the coarse meshes are more smeared than the ones on the fine mesh. Next, the wedge angle is increased to ⫽ 35°, which is near the maximum angle for which the equilibrium shock-polar admits a solution. As before, three different mesh sizes have been used for this problem, consisting of 560 ⫻ 64 cells, 560 ⫻ 128 cells, and 1140 ⫻ 240 cells, respectively. The finest computational domain corresponds to a nominal resolution of 8 points per half-reaction length. The upstream velocity is u 1 ⫽ 11.509. The asymptotic limit is a detonation with overdrive factor f ⫽ 1.2 and at an angle  ⫽ 56.8°. Results for this simulation, taken at t ⫽ 50.0, are given in Figs. 5a, 5b, and 5c. The flow can not turn smoothly in this case, because of the high value of the wedge angle. As a result, a strong explosion takes place at the front. The center of this explosion is a triple point. The incident shock and the Mach stem are the two parts of the main front, below and above the triple point respectively. Another shock emanates from the triple point, the reflected shock, which hits the wedge and reflects back. Additionally, a contact discontinuity (shear layer) is formed between the Mach stem and the reflected shock. The material behind the incident shock burns very slowly because the shock is relatively weak and the temperature increase is small. The material behind the Mach stem, however, burns fast because of the high temperature rise. Consequently, there is a strong density and temperature gradient across the shear layer. The shear layer becomes unstable quickly and generates strong vortical structures that are convected downstream (Fig. 5c). Diffusion in these simulations is introduced by the discretization of the equations and the truncation error of the scheme. Consequently, the size and the velocity of the vortical structures are determined by the size of the grid and the implicit artificial viscosity of the algorithm. The convective Mach number for the shear layer is between 0.5 ⬃ 0.7, therefore the layer can be considered slightly compressible. It is worth mentioning that the flow between the incident shock and the wedge is steady. Furthermore, the region between the shear
532
M. V. PAPALEXANDRIS
Fig. 4. (a) Case B, ⫽ 20°: contour plots of the flow variables at t ⫽ 18.0. Grid size, 285 ⫻ 60 cells. (b) Case B, ⫽ 20°: contour plots of the flow variables at t ⫽ 18.0. Grid size, 570 ⫻ 120 cells. (c) Case B, ⫽ 20°: contour plots of the flow variables at t ⫽ 18.0. Grid size, 1140 ⫻ 240 cells.
layer and the Mach stem is subsonic in the vicinity of the triple point. The boundary of this pocket is a sonic line which emanates from the triple point and lies in the shear layer. After some distance the sonic line curves upward and ends up in the Mach stem. Beyond this subsonic pocket, the flow has all the typical characteristics of planar detonations, such as formation of colliding triple points and transverse waves, and
vorticity generation at the main front. Qualitatively, the results on the three different meshes are the same. Interestingly, the location of the explosion at the front is the same regardless of mesh size. The results on the coarse meshes, however, are more diffusive than than the ones on the fine mesh: the shock profiles are more smeared, and the structures of the shear layer and the reaction zone are more diffuse.
WEDGE-INDUCED DETONATIONS
533
Fig. 5. (a) Case B, ⫽ 35°: contour plots of the flow variables at t ⫽ 50.0. Grid size, 280 ⫻ 64 cells. (b) Case B, ⫽ 35°: contour plots of the flow variables at t ⫽ 50.0. Grid size, 560 ⫻ 128 cells. (c) Case B, ⫽ 35°: contour plots of the flow variables at t ⫽ 50.0. Grid size, 1120 ⫻ 256 cells.
Since the equilibrium shock-polar admits a solution for the shock angle , it is expected that far downstream the Mach stem reduces to an oblique ZND wave, at an angle  ⫽ 56.8°. Linear stability analysis, however, reveals that two-dimensional detonations are intrinsically unstable. Consequently, the instability mechanisms that are encountered in detonations propagating in channels and lead to the formation of triple points (see, for example, [16 –19]) should also appear on the flow-field of the oblique
detonations under study. This region, however, is too far downstream to be included in the computational domain. For higher wedge angles, the equilibrium shock polar can not give a solution for . In such situations, a strong triple point is also formed at the main front. But since there is no solution to the equilibrium shock-polar, the Mach stem can not reduce to an oblique ZND wave downstream, and it is everywhere curved. For even higher wedge angles, none of the shock-polars
534
M. V. PAPALEXANDRIS
Fig. 6. Case C, ⫽ 20°: contour plots of the flow variables at t ⫽ 18.0.
admits a solution, and the main front detaches from the wedge. Such high wedge angles have not been considered in the present work. Case C The activation energy is the same as in case B, but the stiffness coefficient has now been increased: E a ⫽ 50,
K ⫽ 230.75.
The wedge angle is ⫽ 20°. The upstream velocity is u 1 ⫽ 18.051. The theoretical prediction is that far downstream the detonation will be a ZND wave of overdrive factor f ⫽ 1.60 and at an angle  ⫽ 28.5°. The computational domain consists of 1200 ⫻ 200 cells, corresponding to a nominal resolution of 8 points per half-reaction length of the one-dimensional, steady-state solution. As in case B, the shock turns smoothly until it reaches a steady state. Contour plots of the flow variables are given in Fig. 6. These results are taken at time t ⫽ 18.0. No change in the flow variables could be observed after that time. The shock angle at the right boundary is very close to the asymptotic limit  ⫽ 28.5°. The case of a higher wedge angle, namely ⫽ 30°, has also been considered. The computational domain consists of 1440 ⫻ 240 cells, corresponding to a nominal resolution of 8 points per half-reaction length. The upstream velocity is now set at u 1 ⫽ 12.035. Under this
Fig. 7. Case C, ⫽ 35°: contour plots of the flow variables at t ⫽ 36.0.
initial condition, the theoretical prediction is a detonation of overdrive factor f ⫽ 1.6, at a shock angle  ⫽ 45.7°. Results for this simulation, taken at time t ⫽ 36.0, are given in Fig. 7. These results are similar to the ones obtained in the high-angle simulation of case B. The basic features that were encountered in the results of case B (such as the explosion at the front, the development of an unsteady shear layer, and the formation of a subsonic pocket behind the Mach stem in the vicinity of the triple point) can also be observed in this simulation. It is worth mentioning that Li et al. [10] also presented simulations for wedge-induced detonations using the FCT algorithm. Their numerical results agree qualitatively with the results obtained with the proposed unsplit algorithm. In particular, they also observed that for small wedge angles the main front turns smoothly, whereas an explosion occurs at the front if large wedge angles are considered. Quantitative comparisons between the two studies can not be made because Li et al. considered a multispecies combustion model and did not use dimensionless quantities for their simulations. DETONATIONS INDUCED BY SHORT WEDGES In the cases examined so far, it was assumed that the wedge is long enough for the explosion
WEDGE-INDUCED DETONATIONS
535
Fig. 8. Flow past a wedge of weight of height h w ⫽ 70 and angle ⫽ 35°. Contour plots of the flow variables at t ⫽ 46.0.
to occur upstream of the corner of the wedge. In such cases, the location of the triple point and the flow-field in that neighborhood are determined completely by the kinetics of the reaction. If, however, the wedge is not long enough, then the explosion will take place near the corner. The effect of the corner in such situations has also been studied numerically, and the results are described below. The chemical reaction increases the temperature of the fluid, while the expansion at the corner decreases it. These two mechanisms “compete” against each other. Furthermore, by including the corner in the computational domain, a second characteristic length is introduced. This is the height of the wedge, h w , the first one being the half-reaction length. It is the combination of these two length-scales with the wedge angle that ultimately determines which of the two mechanisms will dominate. In the numerical simulations, the activation
energy and the stiffness coefficient are set at E a ⫽ 50 and K ⫽ 99.762, respectively. The wedge angle is ⫽ 35°. The upstream state of the fluid is p 1 ⫽ 1.0,
1 ⫽ 1.0,
u 1 ⫽ 11.509.
The flow that is produced by these parameters and upstream conditions in the case of a long wedge was examined in the previous section (see Fig. 5c). Three different wedge heights are considered. First, the wedge height is set at h w ⫽ 70.0. The computational domain consists of 930 ⫻ 330 cells. The length of the domain is 155.0, and its width from the upper wall of the wedge is 55.0. Results for this simulation, taken at t ⫽ 46.0, are shown in Fig. 8. In this case the explosion occurs upstream with respect to the corner. The expansion at the corner does not affect the explosion because the flow in the
536
M. V. PAPALEXANDRIS
Fig. 9. Flow past a wedge of weight of height h w ⫽ 50 and angle ⫽ 35°. Contour plots of the flow variables at t ⫽ 46.0.
corner is supersonic. The leading shock is expected to reduce far downstream to a Chapman-Jouguet (CJ) wave, just like the one-dimensional detonation initiated by a moving piston that comes suddenly to a rest. The flow-field in the neighborhood of the inert shock and the explosion is the same as in the corresponding case with a long wedge (see Fig. 5c). The interaction of the reflected shock with the shear layer causes the fluid to decelerate, thus forming another subsonic region in the flow-field, besides the one in the vicinity of the triple point. It is also worth mentioning that the expansion that takes place at the corner affects the evolution of the shear layer. As the fluid below the shear layer expands, the pressure and density drop, generating even higher entropy gradients across the layer. Therefore, the expansion at the corner leads to an increased amount of vorticity generation across the shear layer.
Subsequently, the wedge height is reduced to h w ⫽ 50.0. The computational domain consists of 810 ⫻ 420 cells. The length of the domain is 135.0, and its width, as measured from the upper wall of the wedge, is 70.0. Results for this simulation are shown in Fig. 9. They are taken at time t ⫽ 46.0. In this case, the explosion takes place downstream of the corner. The flow in the subsonic area in the vicinity of the triple point (between the Mach stem and the shear layer) is influenced by the expansion at the corner. Additionally, the expansion at the corner affects the curvature of the leading front and the reflected shock. A fluid element moving parallel to the wedge does not have time to increase its temperature substantially, because of the small length of the wedge. Consequently, it remains almost unreacted when it reaches the head of the expansion. The expansion produces a further decrease of
WEDGE-INDUCED DETONATIONS
537
Fig. 10. Flow past a wedge of weight of height h w ⫽ 25.0 and angle ⫽ 35°. Contour plots of the flow variables at t ⫽ 10.0.
the temperature which delays the initiation of the reaction even more. It is observed that the material on the upper wall of the wedge has remained only partially reacted (about 20%), as a result of the temperature decrease caused by the expansion at the corner. As a final test, the wedge height is lowered further to h w ⫽ 25.0. The simulation is performed on a domain of 1266 ⫻ 270 cells. The length of the domain is 210.0 and its width is 45.0. Results from this simulation, taken at t ⫽ 10.0, are presented in Fig. 10. It is observed that the expansion at the corner reduces the temperature so much that the gas near the wedge remains almost unreacted because the time needed for the temperature to increase high enough to establish rapid reaction becomes very large. As a result, a detonation can not be established and the shock wave is expected to reduce to a Mach wave downstream. It is also observed that there is a small pocket of slightly reacted material near the wedge. It can be verified by looking at results taken at early times that the formation of this pocket is a transient phenomenon caused by the interaction of the shock wave and the expansion at the corner, in the beginning of the simulation. It is convected downstream with the fluid velocity. The minimum value of the reactant mass fraction inside the pocket is z ⯝ 0.9. If the material reacted completely, pressure waves would be transmitted to the shock wave and a CJ detonation could be established. But it is observed that the reaction rate, ˙z ⫽ ⫺Kz exp (⫺E a/T), inside the pocket decays with time, which suggests that
the reaction process will not be completed inside the pocket. After some time the pocket exits the computational domain, and the shock front eventually assumes a fixed position. No change in the flow variables can be observed after that. The material all along the wedge will remain only partially reacted, and a detonation will not be established. CONCLUDING REMARKS A numerical study of wedge-induced detonations has been presented. The numerical scheme that was used in the simulations integrates all the terms of the governing equations simultaneously, without the use of flux-splitting or additional artificial viscosity. It was confirmed that, for small wedge angles, the shock that is attached to the wedge turns smoothly to an oblique ZND wave. For high wedge angles, however, such a smooth turn is not possible and an explosion takes place at the front. The center of the explosion is a triple point that assumes a fixed position after some time. A shear layer emanates from the triple point that eventually becomes unstable. The effect of the corner of the wedge was also studied numerically. It was shown that when the explosion of the leading front takes place upstream of the corner, the expansion at the corner does not affect the evolution of the front, which reduces to a CJ wave. When the explosion occurs downstream of the corner, the curvature of the front and the reaction process depend on
538
M. V. PAPALEXANDRIS
the expansion at the corner. It appears that for wedge heights small enough, a detonation can not be established downstream, and the front decays downstream to a Mach wave. Future directions include the numerical simulation of multicomponent reacting systems. This will help the study of various problems arising in supersonic combustion, such as the effect of chain-branching. The better understanding of these mechanisms can be helpful in the design of the next-generation, high-speed propulsion systems. The author would like to thank Prof. Leonard for his support and many fruitful discussions. Thanks are also due to Prof. Dimotakis and Prof. Shepherd for their suggestions. This work is sponsored by the Air Force Office of Scientific Research Grant Nos. F49620-94-1-0353 and F49626-93-1-00338, whose support is gratefully acknowledged. REFERENCES 1. 2. 3. 4. 5. 6.
7. 8. 9. 10.
Brackett, D. C., Bogdanoff, D. W., J. Propul. Power 5:276 –281 (1989). Grismer, M. J., Powers, J. M., J. Propul. Power 11:105– 111 (1995). Lefebvre, M. H. and Fujiwara, T., Combust. Flame 100:85–93 (1995). Yungster, S., Radhakrishnan, K., Rabinowitz, M. J. (1995). AIAA paper 95–2489. Buckmaster, J. D., Combust. Sci. Technol. 72:283–296 (1990). Shepherd, J. E., in Combustion in High-Speed Flows (J. Buckmaster et al., Eds.), Kluwer Academic Publishers, 1994, pp. 373– 420. Powers, J. M., Stewart, D. S., AIAA J. 30:726 –736 (1992). Glenn, D. E., Pratt, D. T. (1998). AIAA paper 88-0440. Cambier, J. L., Adelman, H., Menees, G. P., J. Propul. Power 5:482– 491 (1989). Li, C., Kailasanath, K., Oran, E. S., Phys. Fluids 6:1600 –1611 (1994).
11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33.
Gerlinger, P., Algermissen, J., Bruggemann, D., AIAA J. 33:1865–1870 (1995). Thaker, A. A., Chelliah, H. K., Combust. Theory Modelling 1:347–376 (1997). Yungster, S., and Radhakrishnan, K. (1995). SA Report 95-19. Fujiwara, T., Matsuo, A., Nomoto, H. (1988). AIAA paper 88-0098. Grismer, M. J. (1994). Ph.D. thesis, University of Notre Dame. Oran, E. S., Kailasanath, K., Guirguis, R. H., Prog. Aeronaut. Astronaut. 114:155 (1988). Bourlioux, A., Majda, A. J., Combust. Flame 90:211– 229 (1992). Cai, W., AIAA J. 33:1248 –1255 (1995). Papalexandris, M. V., Leonard, A., Dimotakis, P. E., GALCIT Report 97-14, submitted to J. Computa. Phys. Clavin, P., He, L., Williams, F. A., Phys. Fluids 9:3764 – 3785 (1997). Lappas, T., Leonard, A., Dimotakis, P. E., SIAM J. Sci. Comp. 20:1481–1512 (1999). Papalexandris, M. V., Leonard, A., Dimotakis, P. E., J. Comp. Phys. 134:31– 61 (1997). Pratt, D. T., Humphrey, J. W., Glenn, D. E., J. Propul. Power 7:837– 845 (1991). Erpenbeck, J. J., Phys. Fluids 7:684 – 696 (1964). Fickett, W., Davis, W. C., Detonation, University of California Berkeley Press, Berkeley, CA, 1979. Lee, H. I., Stewart, D. S., J. Fluid Mech. 206:103–132 (1990). Majda, A., Roytburd, V., Studies Appl. Math. 87:135– 174 (1992). Short, M., Stewart, D. S., J. Fluid Mech. 340:249 –295 (1997). Short, M., SIAM J. Appl. Math. 57:307–326 (1997). Klein, R., Stewart D. S., SIAM J. Appl. Math. 53:1401– 1435 (1993). Lasseigne, D. G., Hussaini, M. Y., Phys. Fluids A 5:1047–1058 (1993). Menikoff, R., Lackner, K. S., Bukiet, B. G., Combust. Flame 104:219 –240 (1996). Stewart, D. S., Bdzil, J. B., Combust. Flame 72:311–323 (1988).
Received 26 January 1998; revised 20 May 1999; accepted 21 July 1999