Numerical Simulation of Pulse Detonation Engine ... - Springer Link
Journal of Scientific Computing, Vol. 19, Nos. 1–3, December 2003 (© 2003)
Numerical Simulation of Pulse Detonation Engine Phenomena X. He 1 and A. R. Karagozian 2 Received June 17, 2002; accepted (in revised form) October 25, 2002 This computational study examines transient, reactive compressible flow phenomena associated with the pulse detonation wave engine. The PDWE is an intermittent combustion engine that relies on unsteady detonation wave propagation for combustion and compression elements of the propulsive cycle. The present computations focus on high order numerical simulations of the generic PDWE configuration with simplified reaction kinetics, so that rapid, straightforward estimates of engine performance may be made. Both one- and two-dimensional simulations of the high speed reactive flow phenomena are performed and compared to determine the applicability of 1D simulations for performance characterization. Examination of the effects of the combustion reaction mechanism and the use of a pressure relaxation length for 1D simulations is made. Characteristic engine performance parameters, in addition to engine noise estimates within and external to the detonation tube, are presented. KEY WORDS: Detonations; gasdynamics; reactive Euler equations; shock capturing schemes.
1. INTRODUCTION AND BACKGROUND The Pulse Detonation Wave Engine (also called the Pulse Detonation Engine) or PDWE is an intermittent combustion engine that relies on unsteady (pulsating) detonation wave propagation for combustion and compression elements of the propulsive cycle. While specific configurations vary, in general the PDWE consists of propagating detonation waves generated periodically within an engine tube, with associated reflected expansion and compression waves which can act in periodic fashion to produce high forward thrust [1–3]. The generic reactive and non-reactive wave phenomena occurring in the PDWE cycle are described in Fig. 1. The combustion process (as a detonation or a shock/deflagration wave which rapidly 1
2
Graduate Student, Department of Mechanical and Aerospace Engineering, University of California, Los Angeles, California 90095-1597. E-mail: [email protected] Professor, Department of Mechanical and Aerospace Engineering, University of California, Los Angeles, California 90095-1597. E-mail: [email protected]
201 0885-7474/03/1200-0201/0 © 2003 Plenum Publishing Corporation
202
He and Karagozian reflected expansion wave
detonation front
reactants reactants
a)
f)
expansion wave
propagating detonation
b)
products
reactants
g)
products
h)
products
i)
products
reactants
shock/detonation reflection
expansion wave
e)
reactants
compression wave
reflected expansion wave
d)
reactants
reflected compression wave
detonation
c)
products
enter
j)
reactants
Fig. 1. The generic Pulse Detonation Wave Engine (PDWE) cycle.
transitions to a detonation) is initiated at the inlet of the thrust or detonation tube, propagating into a combustible mixture of reactants (Figs. 1a and 1b), assumed in the figure to initially fill the entire tube. The propagating detonation leaves a mixture of combustion products behind the wave. The detonation reaches the end of the tube, propagates outward as a shock into the surrounding air, and an expansion wave simultaneously reflects back into the tube (Figs. 1c and 1d). The expansion wave produces lower pressure gas in its downstream region, further expelling the combustion products situated behind the wave. When the expansion reaches the thrust wall or plate at the forward end of the thrust tube, it is reflected as an expansion wave (Fig. 1e). The pressure behind the reflected expansion is further reduced, below ambient, and air and fuel flow into the tube through their respective inlets during this ‘‘aspiration phase’’, as the expansion wave propagates toward the exit (Fig. 1f ). When the expansion wave reaches the exit (Fig. 1g), it is reflected as a compression wave, coalescing to a shock, from the open end and the wave propagates upstream toward the thrust wall (Fig. 1h). The reflection of the wave at the thrust wall delivers thrust to the vehicle and is reflected as a shock, allowing the combustible mixture of reactants to ignite, with transition to a propagating detonation wave (Figs. 1i and 1j) and a return to the initial condition of the cycle (Fig. 1a).
Numerical Simulation of Pulse Detonation Engine Phenomena
203
During propagation of the detonation, combustion of reactants takes place so rapidly that the process effectively takes place at a constant volume, which is theoretically more efficient than a constant pressure combustion process as essentially occurs in conventional airbreathing or rocket engines. The absence of rotating machinery suggests that the PDWE concept holds real promise for high thrust density, low fuel consumption propulsion applications, either in an airbreathing or rocket engine mode. A number of alternative configurations for the PDWE have been proposed and tested [4], and numerical simulations of PDWE phenomena are being actively pursued. Helman et al. [5] were the first to demonstrate a successful air-breathing detonation engine, with intermittent detonation frequencies of 25 Hz. Since then, extensive testing has been and continues to be pursued by a variety of groups [3, 4]. An excellent historical overview of the development of the PDWE as a constant volume combustion engine concept is provided in Eidelman et al. [1]. Overviews of past and ongoing numerical simulations of PDWEs are described in recent articles by Kailasanath [2, 6–8]. Kailasanath and Patnaik [9] explore the effect of employing a variable ‘‘pressure relaxation length’’ external to the end of the PDWE’s detonation tube. The relaxation length is used in order to mimic multidimensional wave effects in the pressure near the tube exit when using a 1D simulation, in particular when the local exit flow is subsonic over a portion of the cycle. Kailasanath and Patnaik find that the relaxation length can have a significant effect on engine performance, with increases in relaxation length leading to increases in engine impulse. The impulse I is typically defined as .
I — A F Dptw (t) dt
(1)
0
where A is the area of the thrust wall and Dptw (t) is the time-dependent pressure differential (compared with initial tube conditions) at the thrust wall. The engine’s specific impulse Isp may also be quantified as a performance parameter: Isp —
I rVg
(2)
where r is the initial mass of the gas mixture in the tube, V is the detonation tube volume, and g is the earth’s gravitational acceleration. An alternative performance parameter, the fuel-based specific impulse, Isp, f , is defined as Isp, f —
Isp Yf
(3)
where Yf is the fuel mass fraction present within the initial premixed reactants in the tube. The effect of partially filling the detonation tube with reactants, in the vicinity of the thrust wall, and having the remainder of the tube filled with air, is explored by Cambier and Tegner [10] and by Li et al. [11, 12], in both 1D and 2D simulations. Reducing the region of the tube occupied by reactants results in a higher Isp, f for a given engine cycle, clearly, due to the reduction in the effective Yf . This occurs
204
He and Karagozian
despite the fact that when the detonation propagates through the interface between reactants and air, reflected expansion waves return to the thrust wall and reduce the pressure at the wall at an earlier time than the reduction that occurs when the expansion waves reflected from the tube exit encounter the wall. An analytical model for the impulse of the single cycle PDWE has been developed by Wintenberger et al. [13]. Employing a similarity solution for the temporally constant pressure region within the tube (as occurs during the portion of the cycle shown in Figs. 1a–1e), with the subsequent ‘‘blowdown’’ process modeled using dimensional analysis and empirical relations, the model is able to impressively capture the essential physical phenomena present in the PDWE. Correspondence of the model results to performance parameters quantified in experiments and in numerical simulations is quite good. The present study focuses on both one- and two-dimensional simulations of the pulse detonation engine using high order numerical schemes, developed over the years by Stanley Osher and collaborators. We focus not only on exploring the PDWE’s reactive flow processes and performance parameters, but also on examining noise generation by the PDWE. NASA’s interest in the PDWE for advanced vehicle propulsion [14] is incumbent upon the ability of the engine to efficiently generate thrust without having to pay a significant penalty in engine noise. 2. PROBLEM FORMULATION AND NUMERICAL METHODOLOGY 2.1. Governing Equations The governing equations which describe the flow and reaction evolution in the PDWE with a single step, irreversible chemical reaction consist of the conservation of mass, momentum, energy, and species. The following representation was used, for example, in simulations of the two-dimensional planar PDWE configurations under investigation: F (U F )x +G F (U F )y =S F (U F) F t +F U
(4)
F , the flux vectors, F F and G F , and where the vector containing conserved variables, U F , are, respectively, the vector containing source terms, S
RS R S R S R S r
ru F = rv U
ru
ru +p
F (U F )= F
ruv
0
rv
2
0
ruv
F (U F )= G
rv 2+p
E
(E+p) u
(E+p) v
rY
ruY
rvY
F )= S(U
0 0
− KrYe −
12 Ti T
(5) and where E may be written r(u 2+v 2) p + +rqY E= c−1 2
(6)
Numerical Simulation of Pulse Detonation Engine Phenomena
205
Here r represents density, p is the static pressure, u and v are the x- and y-components of the velocity vector, respectively, and c is the ratio of specific heats. q is a heat release parameter which characterizes the amount of energy released during the reaction, and Ti represents the activation temperature. Y is the reactant mass fraction, which varies from 0 to 1, while K is the reaction-rate multiplier for the reaction, which also sets the spatial and temporal scales in the problem. One-step chemical reactions were considered in the present study, mostly for the methane-oxygen reaction, but also for the hydrogen-oxygen reaction. These were sufficient to enable reasonable estimates for performance parameters to be derived from the computations. The reaction terms K and Ti were determined for the methane-oxygen reaction from data provided in [15]. For the hydrogen-oxygen reaction, the STANJAN package was used to compute the equilibrium conditions downstream of an H 2 –O2 detonation, and K and Ti were fit to these data. We found that the PDWE cycle results were not sensitive to our selections of K and Ti , within reasonable limits, for this reaction, and hence these values were chosen to be the same as those for the methane-oxygen reaction. The values of these reaction quantities for the two reactions considered are shown in Table I. Two methods for initiation of the detonation reaction were explored. One method employed a computational ‘‘spark’’ adjacent to the thrust wall, by which a narrow high pressure, high temperature region (8 grid cells in width) initiated the propagating shock and ignited the reactants. The second method initiated the computation via an incoming shock moving through the reactants toward the thrust wall, which reflected from the wall and ignited the detonation. These initiation mechanisms produced virtually the same overall performance parameters, since the incoming shock strength in the second method could be chosen so that the pressure at the thrust wall rapidly rose to the same extent using either ignition method. As such, only results for the ‘‘spark’’ ignition of reactants are shown here. In most results shown, except as otherwise stated, the computational ‘‘spark’’ was represented by a region with an initial pressure of 20 atm and an initial temperature of 6000 K. At much lower spark region initial temperatures and pressures, of course, there would be a much more gradual transition from the deflagration (subsonic combustion) to the detonation process [16]. Both one- and two-dimensional numerical simulations were performed in this study. In addition to the 2D planar problem described in Eq. (4) above, corresponding equations for the axisymmetric PDWE configuration were solved, as were the equations for a 1D PDWE configuration. The percentage ‘‘fill’’ of the reactants
Table I. Values of Constant Terms in Relations (5) and (6) for Approximate Single Step Reactions Reaction
q (J/kg)
K
Ti (K)
CH 4 –O2 H 2 –O2
1.0 × 10 7 1.344 × 10 7
3.224 × 10 8 3.224 × 10 8
24358 24358
206
He and Karagozian
within the tube was varied in some of our computations, but results for the influence of this parameter are not shown in this paper since they basically replicated the trends described in Li and Kailasanath [12]. In the case of the 1D simulations, the use and effect of employing a pressure relaxation length exterior to the end of the detonation tube, after Kailasanath and Patnaik [9], was also explored. Incorporating a relaxation length pre-defines a spatial decay in the pressure immediately external to the PDWE tube after the exit of the detonation/shock (see Fig. 1c). When the relaxation length l was set to be greater than zero, the external pressure was allowed to spatially decay in a linear fashion over the length l to atmospheric pressure, rather than having the exit plane pressure exposed immediately to atmospheric pressure. The relaxation length assumption is suggested by Kailasanath and Patnaik [9] to enable the flow interior to the PDWE tube (i.e., the nature of the reflection of an expansion wave back into the detonation tube) to be more closely representative of the phenomena that actually take place in a multi-dimensional PDWE, particularly when exit conditions from the tube become locally subsonic and the interior flow can respond to changes outside of the tube. In the present 1D simulations, however, the computational domain consisted primarily of the detonation tube (containing 600 grid points), since the effect of a pressure relaxation length on the wave reflected back into the detonation tube could be represented using only a few grid points beyond the tube end in order to capture the spatial slope in pressure. In the 2D simulations, on the other hand, the air external to the detonation tube was assumed to be uniformly at atmospheric pressure, and the computational domain extended well downstream of the end of the tube, in general at least one and one half tube lengths downstream. The 2D domain extended at least two tube diameters away from the detonation tube in the dimension perpendicular to the axial dimension. Hence for the 2D simulations, no relaxation length approximation was necessary, and outflow boundary conditions were employed at the edges of the computational domain. Comparisons were thus possible between the 2D simulations, which captured flow characteristics external and internal to the PDWE, and 1D simulations with and without a pressure relaxation length, in which only the interior of the PDWE and a few exterior grid cells were resolved.
2.2. Numerical Methodology and Code Validation The present study used the Weighted Essentially Non-Oscillatory (WENO) method [17], a derivative of the Essentially Non-Oscillatory (ENO) method, introduced and developed by Stanley Osher and his colleagues [18–21] for spatial interpolation of the system of governing equations. ENO methods constitute a class of high accuracy, shock capturing numerical schemes for hyperbolic systems of conservation laws, based on upwind biased differencing in local characteristic fields. It has high accuracy (third order or higher) in smooth regions of the flow, and captures the motion of unresolved steep gradients in general without introducing spurious oscillations. Hence ENO is particularly well suited for resolution of
Numerical Simulation of Pulse Detonation Engine Phenomena
207
flowfields in which there are shocks or flame fronts. ENO uses an adaptive polynomial interpolation constructed on the basis of decisions to avoid steep gradients in the data; these decisions involve selection of the smoothest computational stencil in the vicinity of the gradient to approximate fluxes. The polynomial is also biased to extrapolate data from the direction of information propagation (i.e., upwind) for physical consistency and stability. WENO schemes, developed by Jiang and Shu [17], are based on ENO and incorporate all possible stencils in approximating fluxes, so that higher order accuracy may be achieved with a smaller number of points in a given stencil [22]. For example, with 3 points per stencil, an ENO scheme becomes third order in accuracy, but a WENO scheme becomes fifth order accurate in smooth regions and remains third order accurate in the vicinity of discontinuities. To avoid entropy-violating expansion shocks near sonic points, where characteristic velocities change sign, high order dissipation was added in the present study via the Local Lax Friedrichs (LLF) scheme [23], which adds extra numerical viscosity throughout the computational domain at each time step. Once the numerical approximation of the spatial terms was completed using WENO, the conservation equations (4) were written in the form F (U F) F t =f U
(7)
A third order TVD Runge–Kutta method was used for the time discretization [23], since the method has high accuracy and, more importantly, a large time step stability region which includes a segment of purely imaginary linear growth rates. This feature ensures that for a sufficiently small time step (according to the Courant–Friedrichs–Levy or CFL condition), the time discretization will not introduce spatial oscillations to the result. Since the current simulation only involved a single step reaction, use of the TVD Runge–Kutta method was preferred. The CFL condition in this case incorporated the source term to insure stability, since reaction rates often impose the greatest constraints on time integration. The ENO/WENO schemes were tested on a variety of problems, including that of the classical one-dimensional, overdriven, pulsating detonation [24]. It is found that the degree of spatial resolution of the reaction zone can have a profound effect on the peak pressure amplitude and period of oscillation of the overdriven pulsating detonation wave. A minimum required resolution of 20 grid points per detonation reaction zone half-length is suggested in [24] for the canonical pulsating detonation problem with an overdrive of 1.6; much finer grids may be required for other types of spatial integration schemes [25, 26]. Grid resolution tests in the present PDWE study, however, suggested that even 5 grid points per reaction zone half-length were sufficient for the cases under consideration. The present WENO scheme was also tested on a non-reactive problem with an analytic solution, that of the one-dimensional propagation of a normal shock wave out of a straight tube into air at atmospheric conditions. This problem yields reflection of an expansion fan back into the tube, mimicing the processes taking place when a detonation leaves the PDWE tube (see Figs. 1c and 1d). As with the one-dimensional PDWE simulations, the computational grid consisted of only the
208
He and Karagozian 250000
250000
Exact solution dx = 0.10 m dx = 0.02 m dx = 0.01 m
225000
200000
Pressure (pa)
Pressure (pa)
200000
175000
175000
150000
150000
125000
125000 100000
100000
75000
75000
50000
0
2
4
6
8
50000
10
0
2
4
6
X (m)
X (m)
(a)
(b)
250000
8
10
250000
Exact solution dx = 0.10 m dx = 0.02 m dx = 0.01 m
225000
225000
200000
Exact solution dx = 0.10 m dx = 0.02 m dx = 0.01 m
Pressure (pa)
Pressure (pa)
200000
175000
175000
150000
150000
125000
125000
100000
100000
75000
75000
50000
Exact solution dx = 0.10 m dx = 0.02 m dx = 0.01 m
225000
0
2
4
6
8
10
50000
0
2
4
6
X (m)
X (m)
(c)
(d)
8
10
Fig. 2. Validation of the 3rd order WENO scheme on a one-dimensional, non-reactive problem with an analytic solution. At time t=0, a diaphragm separating high pressure gas (phigh =2.0 × 10 5 Pa, rhigh =1.78 kg/m 3, x < 7 m) from lower pressure gas (plow =1.0 × 10 5 Pa, rlow =1.0 kg/m 3, x > 7 m) is burst. On either side of the diaphragm, c=1.2 and initially, T=300K. The tube end is situated at x=10 m, and for x > 10 m at all times, the pressure is 1 × 10 5 Pa. Shown are spatial distributions of pressure in the vicinity of the open tube end (within the tube) at various times: (a) t=0.8 ms, (b) t=4.0 ms, (c) t=12.0 ms, and (d) t=16.0 ms. Comparisons are made among simulation results for various spatial grid sizes (symbols) and the analytic or exact solution (lines). Propagation of the shock wave to the right ((a), (b)), out of the tube, results in reflection of an expansion wave back into the tube ((c), (d)). Here the flow throughout the domain shown remains subsonic (in the laboratory reference frame) for all times shown.
tube and a few grid cells exterior to or downstream of the tube. Figure 2 shows the computed pressure distribution in this non-reactive test problem at four different times. There was virtually no difference between the 3rd order WENO-computed solution and the analytic solution when the spatial grid spacing was less than about 2 cm in this test case. While ENO schemes can become expensive due to the required decisions on the stencils, WENO schemes were found to be less costly, and moreover, the ability to use relatively coarse grids, as seen in [24], helps to reduce
Numerical Simulation of Pulse Detonation Engine Phenomena
209
the overall computational time. Both the shock propagating out of the tube (right side of Figs. 2a and 2b) and the expansion wave reflected back into the tube (right side of Figs. 2c and 2d) were captured very well, even though only a few grid points outside of the end of the tube at atmospheric pressure were employed in the computation. The flow throughout the tube remained subsonic (or, prior to passage of the shock, quiescent) during the computation. This comparison lends confidence that the present scheme and boundary conditions, even in 1D, were capable of capturing important compressible wave dynamics at the end of the detonation tube. 2.3. Noise Estimates As described above and in Wintenberger et al. [13], for example, a variety of performance parameters may be derived from the spatial and temporal evolution of pressure and other flow properties in the PDWE. In addition to the impulse I, specific impulse Isp , and fuel specific impulse, Isp, f , the sound pressure level (SPL) at various locations within and external to the detonation tube were computed in this study. These noise levels were roughly estimated by examining the Fourier transform of the time-dependent pressure measured at various locations within the computational domain. A peak in the spectrum often appeared at roughly 330 Hz, which corresponded to the inverse of the period of the PDWE cycle for a tube length of one meter and a CH 4 –O2 reaction. Smaller peaks at higher harmonics were usually observed. The maximum pressure peak was used to estimate the sound pressure level or SPL at a given physical location, where the SPL is defined as SPL=20 log10
p pref
(8)
and where the value of pref is 2 × 10 −5 N/m 2 or Pascals (Pa). Hence the present computations allow rough estimates of the distribution of sound pressure level, or noise, emanating from pulse detonation engine to be made. 3. RESULTS The temporal evolution of flow and reaction structures within the pulse detonation engine were computed over a single cycle. An example of the temporal evolution of the computed 2D planar pressure field associated with the PDWE, for the CH 4 -O2 reaction, is shown in Fig. 3. A straight, planar PDWE tube one meter in length and 0.2 m in width was simulated. The tube was assumed here to be completely filled with a stoichiometric mixture of methane and oxygen at atmospheric pressure. In the 2D simulations, the grid spacing in the downstream dimension was Dx=1 cm or less, while the grid size in the transverse direction was Dy=0.3 cm. A CFL number of 0.1 was used in the time discretization. The propagation of the detonation to the right, within the tube, is readily apparent; its exit from the tube end resulted in the propagation of a vortical structure, coincident with the leading shock, into the surrounding air. Simultaneous to the exit of the shock and initiation of the vortex, a reflected expansion fan entered the detonation
210
He and Karagozian
Fig. 3. Temporal evolution of the 2D planar pressure field within and external to the PDWE over one cycle, with data shown at times corresponding to (a) 0.0 ms, (b) 0.06 ms, (c) 0.15 ms, (d) 0.26 ms, (e) 0.49 ms, (f ) 0.76 ms, (g) 1.01 ms, and (h) 2.89 ms. A CH 4 –O2 reaction was simulated here. The color map for the pressure distribution is shown, with units of Pascals.
Numerical Simulation of Pulse Detonation Engine Phenomena
211
Fig. 4. Temporal evolution of the centerline pressure within and external to the PDWE over one cycle, extracted from a 2D planar simulation of the flowfield, as in Fig. 3, with data shown at times corresponding to (a) 0.0 ms, (b) 0.06 ms, (c) 0.15 ms, (d) 0.26 ms, (e) 0.49 ms, (f ) 0.76 ms, (g) 1.01 ms, and (h) 2.89 ms.
212
He and Karagozian (e)
(a) .
(b) . Detonation tube
.
(f)
(c) .
.
(d) . Patm
Fig. 5. Locations within and exterior to the detonation tube (shown here without a nozzle) at which pressure data and SPL are computed. Here, along the centerline of the tube, (a) refers to the thrust wall, (b) refers to the tube midpoint, (c) refers to the exit plane of the tube, and (d) refers to a point one tube length external to the tube. Points (e) and (f ) lie 0.4 tube lengths above the tube centerline, aligned with the tube exit plane and at a location one tube length beyond the exit plane, respectively.
tube, whose subsequent reflection from the thrust wall resulted in nonlinear wave interactions evolving within the tube. The temporal evolution of the pressure distribution along the centerline of the tube for the same 2D simulation as in Fig. 3 is shown in Fig. 4. Here the initiation (Fig. 4a) and propagation of the detonation wave through the tube (Figs. 4b–4d), the exit of the shock from the tube (Figs. 4d and 4e), and the reflection of the expansion fan from the exit plane back into the tube (Figs. 4e–4g) are more obvious. The temporal evolution of pressure at various locations within and external to the PDWE may be derived from results such as those shown in Fig. 3. Data are presented here at a variety of points in the flowfield, labeled as shown in Fig. 5. For the 2D planar geometry, corresponding to the conditions in Fig. 3, the time-dependent pressures at the various locations (a)–(f ) identified in Fig. 5 are plotted in Figs. 6a–6f, respectively. A comparable computation but for the 2D axisymmetric detonation tube with a diameter of 0.2 m produces the results shown in Fig. 7. Interestingly, while the pressure evolution at the thrust wall and within the tube was not very different between the planar and axisymmetric configurations, the propagation and evolution of the shock structure exterior to the tube broadened over a longer time period for the axisymmetric configuration as compared with the planar configuration. In both planar and axisymmetric cases, the local pressure at point (c), external to the PDWE, was observed to drop to lowered density conditions for a short time (see Figs. 6e and 7e). This corresponded to the time at which the low pressure region associated with the inviscid, compressible vortex centers passed this point. The influence of these external flow differences between 2D and axisymmetric cases on performance parameters for the PDWE were very small, however, as seen by the plots of I, Isp , and Isp, f as a function of time shown in Fig. 8 for both planar and axisymmetric configurations. It is further noted that the time-dependent pressure at the thrust wall is of the same functional form as that observed in recent PDWE experiments performed at Stanford University [27] for a stoichiometric ethylene-oxygen mixture but with a 1.35 m long detonation tube.
Numerical Simulation of Pulse Detonation Engine Phenomena
213
Fig. 6. Temporal evolution of the pressure over one PDWE cycle, for points located at (a)–(f ) as labeled in Fig. 5. The 2D planar configuration with a CH 4 –O2 reaction was simulated here.
214
He and Karagozian
Fig. 7. Temporal evolution of the pressure over one PDWE cycle, for points located at (a)–(f ) as labeled in Fig. 5. The 2D axisymmetric computation with a CH 4 –O2 reaction was simulated here.
Numerical Simulation of Pulse Detonation Engine Phenomena
215
Fig. 8. Performance parameters for the 2D planar and 2D axisymmetric simulations of the detonation tube: (a) Impulse I per unit area, (b) specific impulse Isp , and (c) fuel specific impulse Isp, f as a function of time. A CH 4 –O2 reaction was simulated here.
Time-series pressure data such as those in Figs. 6 or 7 may be used to estimate the noise generated at various points in the flowfield over a single PDWE cycle. Pressure spectra for the various locations within and external to an axisymmetric PDWE tube are shown in Fig. 9, with correspondence to the points labeled in Fig. 5 and the data in Fig. 7. In all locations we observe the highest peak in pressure to appear close to the frequency associated with the period of the PDWE cycle, roughly 330 Hz for these conditions. Clearly, the numerical shocks, with an artificial thickness, contributed at least some of the higher frequency peaks in these spectra, but as a rough approximation, the SPL may be computed from the lowest frequency’s peak pressure, at 330 Hz. Hence an approximation for the noise levels within and about the detonation tube may be quantified, shown in Table II. While clearly the highest noise levels may be found within the detonation tube, it is interesting to note that the noise level even one tube length downstream of the
216
He and Karagozian
Fig. 9. Sample pressure spectra at points labeled as shown in Fig. 5 for the 2D axisymmetric PDWE tube with a CH 4 –O2 reaction.
Numerical Simulation of Pulse Detonation Engine Phenomena
217
Table II. Computed Sound Pressure Level at Flowfield Locations Labeled in Fig. 5 Location
a
b
c
d
e
f
SPL (dB)
212
211
202
173
174
173
exit plane of the tube may be as high as 173 dB. These data were generally consistent with noise measurements in PDWE tests performed at NASA Glenn Research Center [28], although since the present simulation was inviscid, dissipative phenomena which could reduce external noise levels were not represented here. The decay in SPL from 202 dB at the tube exit (c) to 173 dB one tube length downstream (d) reflects the spatial dispersion of the exiting shock wave. Reduction in the shock/detonation strength within the tube could significantly reduce the noise generation by the PDWE, but of course this has the disadvantage of reducing engine performance. This tradeoff between performance and noise, and the effects of the higher frequencies shown in Fig. 7, will be explored in future studies. To this point, all results shown have been for 2D simulations, either planar or axisymmetric. A significant savings in computational time may be derived if performance parameters and noise estimates are computed from 1D simulations. In performing the 1D simulations, we first explored the effect of employing a pressure relaxation length l, as utilized in the computations of Kailasanath and Patnaik [9]. As soon as the detonation/shock proceeded from the PDWE tube exit in the present computations, the pressure immediately external to the tube was assumed to ‘‘relax’’ linearly between its time-dependent value at the exit to atmospheric pressure at a distance I from the exit. Only a few ‘‘ghost points’’ downstream of the exit were employed in the 1D simulation, so that specification of l and thus the spatial pressure gradient were sufficient to impose this exit condition. Increasing l, from zero up to half a PDWE tube length L, was observed to increase the PDWE’s impulse, shown in Fig. 10 and also as seen in [9]. The reason for this increase in PDWE impulse may be discerned by examining the time-dependent pressure at the thrust wall and at the tube exit, shown in Figs. 11a and 11b, respectively, for cases where the pressure relaxation lengths are l=0 and l=0.5L. While the pressure at the thrust wall was the same for both cases until about 1 msec (the time at which the reflected expansion wave reached the thrust wall), after this time the thrust wall pressure decayed more rapidly when l=0 than when l=0.5L. A similar phenomenon was observed at the tube exit. The pressure relaxation length forced a more gradual spatial drop in the pressure outside of the tube, causing the pressure at the exit plane to exceed that of the case where l=0. In fact, in some cases, incorporating a non-zero pressure relaxation length caused the flow at the tube exit to remain subsonic when it would have become supersonic if l had been set to zero, as will be seen below. To understand if similar PDWE data may be derived from 1D simulations as for 2D simulations, comparisons for the 2D axisymmetric detonation tube with a
218
He and Karagozian
Impulse per area (pa.s)
2000
1500
relaxation length = 0 relaxation length = 0.01 relaxation length = 0.02 relaxation length = 0.10 relaxation length = 0.30 relaxation length = 0.50
1000
L L L L L
500
0
0
0.001
0.002
0.003
0.004
Time (sec.)
Fig. 10. Comparisons of time-dependent PDWE impulse per unit area from 1D simulations using different pressure relaxation lengths, given in terms of the detonation tube length L (1 m). A methane-oxygen reaction is simulated.
methane-oxygen reaction are made in Fig. 12. Here the 1D simulations were conducted for l=0 and l=0.5L, and comparisons were made with the 2D axisymmetric computations by altering the initial temperature and pressure of the computational ‘‘spark’’ in the 1D simulation so that equivalent initial conditions were attained in all three cases. The spatial pressure evolution shows that, prior to the detonation wave’s leaving the tube (Figs. 12a–12d), all three sets of simulations yielded nearly identical results, which is appropriate since the flow within the tube
(a)
40
(b)
35
35 30
relaxation length = 0 relaxation length = 0.50 L
25
Pressure (atm)
Pressure (atm)
30
20
15
20 15
10
10
5
5
0
0
0.001
0.002
Time (sec.)
0.003
0.004
relaxation length = 0 relaxation length = 0.50 L
25
0
0
0.001
0.002
0.003
0.004
Time (sec.)
Fig. 11. Temporal evolution of the pressure over one PDWE cycle, for points located at (a) the thrust wall and (b) the PDWE tube exit. Here the 1D simulation is employed, for the cases of zero relaxation length (l=0) and a relaxation length l corresponding to half a PDWE tube length. A CH 4 –O2 reaction was simulated here.
Numerical Simulation of Pulse Detonation Engine Phenomena
219
was essentially one-dimensional for a single-step reaction, and the boundary conditions at and just beyond the tube exit were initially identical. But when the detonation left the tube, becoming a shock, the expansion fan reflected back into the tube had characteristics that are different between the two different 1D simulations, depending on the boundary specifications at the tube exit and hence on the assumption of a pressure relaxation length. If the pressure just downstream of the exit was forced to spatially decay slowly away from the exit for subsonic outflow, as with the l=0.5L case, rather than abruptly to atmospheric pressure, as in the 1D simulation without a relaxation length, the pressure immediately external to the
Fig. 12. Spatial distribution of the centerline pressure at different times over one PDWE cycle, with comparisons among the 2D axisymmetric simulations (solid line), 1D simulations without the inclusion of a pressure relaxation length (dashed line), and 1D simulations with a pressure relaxation length of half the PDWE tube length (l=0.5L) (dashed-dotted line). A methane-oxygen reaction was simulated. Distributions shown are for approximately the same times: (a) 0.00 ms, (b) 0.06 ms, (c) 0.150 ms, (d) 0.260 ms, (e) 0.370 ms, ( f ) 0.490 ms, (g) 0.760 ms, (h) 1.01 ms, and (i) 2.88 ms.
220
He and Karagozian
tube was forced to become higher for l=0.5L than for l=0. This caused the reflected expansion fan to become relatively weak for l=0.5L as compared with the 2D simulation. The abrupt pressure change at the exit for the 1D simulation with l=0 was actually more severe than for the 2D simulation (see Fig. 12e), but was closer to the 2D results over time (see Figs. 12g–12i). As noted above, the imposition of the pressure relaxation length when the flow at the exit plane was locally subsonic caused the flow at the tube exit to remain locally subsonic after the detonation/shock left the tube, while the case where l=0 in fact produced slightly supersonic exit flow in a short time (see Fig. 13), as did the 2D simulation. While neither 1D simulation precisely replicated the 2D axisymmetric flow’s Mach number solution at the tube exit resulting from the propagating shock/vortical structure, for the operating conditions examined, the l=0 case was able to capture the flow somewhat better. As a consequence, performance parameters were predicted by the 1D simulation with l=0 that were closer to the 2D results, as seen in Fig. 14. This type of comparison was found for many relaxation lengths considered here, however, different initial conditions, tube lengths, and reaction mechanisms could render a non-zero relaxation length more realistic, as seen by Kailasanath and Patnaik [9]. Hence a great deal of care is required in employing the pressure relaxation length concept in order to replicate multidimensional wave phenomena.
2
2D simulation 1D without pressure relaxation 1D with pressure relaxation
Mach Number
1.5
1
0.5
0
0
0.001
0.002
0.003
Time (sec.)
Fig. 13. Comparisons of the Mach number at the PDWE tube exit plane for a methane-oxygen reaction. Results are shown for the 2D axisymmetric PDWE simulations, and for the 1D simulations with a pressure relaxation length of half a tube length (l=0.5L) and with no pressure relaxation length (l=0). Here the 1D simulations incorporated a computational ‘‘spark’’ consisting of a pressure of 10 atm and a temperature of 3000 K in order to match the initial conditions for the 2D simulation.
Numerical Simulation of Pulse Detonation Engine Phenomena
221
Fig. 14. Comparisons of (a) time-dependent pressure at the thrust wall, (b) time-dependent impulse per unit area, and (c) time-dependent specific impulse for both 1D and 2D axisymmetric simulations of the PDWE tube with a CH 4 –O2 reaction. 1D simulations explored the use of a pressure relaxation length l=0.5L. Here the 1D simulations incorporated a computational ‘‘spark’’ consisting of a pressure of 10 atm and a temperature of 3000 K in order to match the initial conditions for the 2D simulation.
Finally, it is noted that when an alternative combustion reaction was simulated, that of hydrogen-oxygen, there was a quantifiable increase in specific impulse as compared with methane-oxygen. This faster H 2 –O2 reaction had a higher q value in the reaction source term (Eq. (6)) and as such generated higher temperatures and pressures behind the detonation, but the H 2 –O2 reaction also had a lower mixture mass than did the CH 4 –O2 reaction. There was thus a correspondingly quantifiable increase in specific impulse, as compared with methane-oxygen, as shown in Fig. 15 for the 1D computation without employment of a relaxation length.
222
He and Karagozian
Specific impulse (sec.)
200
150
CH 4 + O 2 H 2 + O2
100
50
0
0
0.001
0.002
0.003
Time (sec.)
Fig. 15. Comparisons of time-dependent PDWE specific impulse computed from 1D simulations using different single step reaction mechanisms, that for hydrogen and oxygen, and that for methane and oxygen. A zero pressure relaxation length was employed. Here the 1D simulations each incorporated a computational ‘‘spark’’ consisting of a pressure of 20 atm and a temperature of 6000 K.
4. CONCLUSIONS The present simulations suggest that a great deal of useful performance and noise related data may be obtained even from one-dimensional computations of the pulse detonation wave engine with simplified reaction kinetics. Third order WENO schemes were observed to provide robust, accurate simulations of the reactive and non-reactive flow processes occurring within the PDWE. 1D and 2D simulations were observed to produce similar pressure profiles throughout the PDWE tube, especially in the vicinity of the thrust wall, and hence similar engine impulse values, without the inherent need for inclusion of a pressure relaxation length for the conditions examined here. Alternative reaction mechanisms, tube lengths, and initial conditions may require such relaxation, however. Straightforward estimates of sound pressure levels within and external to the engine tube may also be made, from both 2D and 1D simulations, although the generation and effects of higher frequency overtones requires further examination. Future studies will explore these effects on the acoustic field generated by the pulse detonation engine, in addition to the effects of the presence of a nozzle on noise generation. ACKNOWLEDGMENTS The authors wish to extend sincere congratulations to Professor Stanley Osher on the occasion of his 60th Birthday and thanks for a number of years of collegial collaboration at UCLA. The authors also wish to thank Professor Chi-Wang Shu of Brown University for his helpful advice in speeding the WENO computations. This work has been supported at UCLA by NASA Dryden Flight Research Center
Numerical Simulation of Pulse Detonation Engine Phenomena
223
under Grant NCC4-153, with Dr. Trong Bui and Dave Lux as technical monitors, and by the Office of Naval Research under Grant ONR N00014-97-1-0027, with Dr. Wen Masters as technical monitor.
REFERENCES 1. Eidelman, S., Grossmann, W., and Lottati, I. (1991). Review of propulsion applications and numerical simulations of the pulsed detonation engine concept. J. Propulsion and Power 7(6), 857–865. 2. Kailasanath, K. (2002). Recent Developments in the Research on Pulse Detonation Engines, AIAA Paper 2002-0470 (Invited), AIAA 40th Aerospace Sciences Meeting. 3. Roy, G. D. (2001). Practical Pulse Detonation Engines–How Far Are They? ISABE Paper 20011170. 4. Pulse-Detonation Engine (PDE) Workshop, sponsored by the Office of Naval Research, Naval Postgraduate School, Monterey, CA, 1997. 5. Helman, D., Shreeve, R. P., and Eidelman, S. (1986). Detonation Pulse Engine, AIAA Paper 86-1683. 6. Kailasanath, K., Patnaik, G., and Li, C. (2002). The Flowfield and Performance of Pulse Detonation Engines, to appear in the Proc. of the 29th Symposium (Intl.) on Combustion. 7. Kailasanath, K. (2001). A Review of PDE Research – Performance Estimates, AIAA Paper 20010474, AIAA 39th Aerospace Sciences Meeting. 8. Kailasanath, K. (2000). Review of propulsion applications of detonation waves. AIAA Journal 38(9), 1698–1708. 9. Kailasanath, K., and Patnaik, G. (2000). Performance Estimates of Pulsed Detonation Engines, Proc. of the 28th Symposium (Intl.) on Combustion. 10. Cambier, J.-L., and Tegner, J. K. (1998). Strategies for pulsed detonation engine performance optimization. Journal of Propulsion and Power 14(4), 489–498. 11. Li, C., Kailasanath, K., and Patnaik, G. (2000). A Numerical Study of Flow Field Evolution in a Pulsed Detonation Engine, AIAA Paper 2000-0314, AIAA 38th Aerospace Sciences Meeting. 12. Li, C., and Kailasanath, K. (2001). A Numerical Study of Reactive Flows in Pulsed Detonation Engines, AIAA Paper 2001-3933, 37th AIAA/ASME/SAE/ASEE Joint Propulsion Conference. 13. Wintenberger, E., Austin, J. M., Cooper, M., Jackson, S., and Shepherd, J. E. (2001). An Analytical Model for the Impulse of a Single Cycle Pulse Detonation Engine, AIAA Paper 2001-3811, 37th AIAA/ASME/SAE/ASEE Joint Propulsion Conference. 14. Three Pillars for Success: NASA’s Response to Achieve the National Priorities in Aeronautics and Space Transportation, NASA Office of Aeronautics and Space Transportation Technology, 1997. 15. Turns, S. (2000). An Introduction to Combustion, 2nd Ed., McGraw–Hill. 16. Kuznetsov, M. S., Alekseev, V. I., and Dorofeev, S. B. (2000). Comparison of critical conditions for DDT in regular and irregular cellular detonation systems. Shock Waves 10(3), 217–223. 17. Jiang, G. S., and Shu, C. W. (1996). Efficient implementation of weighted ENO schemes. J. Comput. Phys. 126, 202–228. 18. Harten, A., Osher S. J., Engquist, B. E., and Chakravarthy, S. R. (1986). Some results on uniformly high-order accurate essentially nonoscillatory schemes. J. Appl. Numer. Math. 2, 347–377. 19. Shu, C. W., and Osher, S. (1989). Efficient implementation of essentially non-oscillatory shock capturing schemes II. J. Comput. Phys. 83, 32–78. 20. Fedkiw, R. P., Merriman, B., Donat, R., and Osher, S. (1998). The penultimate scheme for systems of conservation laws: Finite difference ENO with marquina’s flux splitting. In Hafez, M. (ed.), Progress in Numerical Solutions of Partial Differential Equations, Arachon, France. 21. Fedkiw, R. P., Merriman, B., and Osher, S. (1997). High accuracy numerical methods for thermally perfect gas flows with chemistry. J. Comput. Phys. 132(2), 175–190.
224
He and Karagozian
22. Lee, T. K., and Zhong, X.-L. (1999). Spurious numerical oscillations in simulation of supersonic flows using shock-capturing schemes. AIAA Journal 37(3), 313–319. 23. Lahey, C. B. (1998). Computational Gasdynamics, Cambridge University Press. 24. Hwang, P., Fedkiw, R., Merriman, B., Karagozian, A. R., and Osher, S. J. (2000). Numerical resolution of pulsating detonation waves. Combust. Theory Model. 4(3), 217–240. 25. Bourlioux, A., Majda, A. J., and Roytburd, V. (1991). Theoretical and numerical structure for unstable one-dimensional detonations. SIAM J. Appl. Math. 51, 303–343. 26. Papalexandris, M. V., Leonard, A., and Dimotakis, P. E. (1997). Unsplit schemes for hyperbolic conservation laws with source terms in one space dimension, J. Comp. Phys. 134(1), 31–61. 27. Jenkins, T. P., Sanders, S. T., Kailasanath, K., Li, C., and Hanson, R. D. (2000). Diode-Laser Based Measurements for Model Validation in Pulse Detonation Flows, Proc. of the 25th JANNAF Airbreathing Propulsion Meeting, Monterey, CA. 28. Perkins, H. D., private communication.
Numerical Simulation of Pulse Detonation Engine Phenomena X. He 1 and A. R. Karagozian 2 Received June 17, 2002; accepted (in revised form) October 25, 2002 This computational study examines transient, reactive compressible flow phenomena associated with the pulse detonation wave engine. The PDWE is an intermittent combustion engine that relies on unsteady detonation wave propagation for combustion and compression elements of the propulsive cycle. The present computations focus on high order numerical simulations of the generic PDWE configuration with simplified reaction kinetics, so that rapid, straightforward estimates of engine performance may be made. Both one- and two-dimensional simulations of the high speed reactive flow phenomena are performed and compared to determine the applicability of 1D simulations for performance characterization. Examination of the effects of the combustion reaction mechanism and the use of a pressure relaxation length for 1D simulations is made. Characteristic engine performance parameters, in addition to engine noise estimates within and external to the detonation tube, are presented. KEY WORDS: Detonations; gasdynamics; reactive Euler equations; shock capturing schemes.
1. INTRODUCTION AND BACKGROUND The Pulse Detonation Wave Engine (also called the Pulse Detonation Engine) or PDWE is an intermittent combustion engine that relies on unsteady (pulsating) detonation wave propagation for combustion and compression elements of the propulsive cycle. While specific configurations vary, in general the PDWE consists of propagating detonation waves generated periodically within an engine tube, with associated reflected expansion and compression waves which can act in periodic fashion to produce high forward thrust [1–3]. The generic reactive and non-reactive wave phenomena occurring in the PDWE cycle are described in Fig. 1. The combustion process (as a detonation or a shock/deflagration wave which rapidly 1
2
Graduate Student, Department of Mechanical and Aerospace Engineering, University of California, Los Angeles, California 90095-1597. E-mail: [email protected] Professor, Department of Mechanical and Aerospace Engineering, University of California, Los Angeles, California 90095-1597. E-mail: [email protected]
201 0885-7474/03/1200-0201/0 © 2003 Plenum Publishing Corporation
202
He and Karagozian reflected expansion wave
detonation front
reactants reactants
a)
f)
expansion wave
propagating detonation
b)
products
reactants
g)
products
h)
products
i)
products
reactants
shock/detonation reflection
expansion wave
e)
reactants
compression wave
reflected expansion wave
d)
reactants
reflected compression wave
detonation
c)
products
enter
j)
reactants
Fig. 1. The generic Pulse Detonation Wave Engine (PDWE) cycle.
transitions to a detonation) is initiated at the inlet of the thrust or detonation tube, propagating into a combustible mixture of reactants (Figs. 1a and 1b), assumed in the figure to initially fill the entire tube. The propagating detonation leaves a mixture of combustion products behind the wave. The detonation reaches the end of the tube, propagates outward as a shock into the surrounding air, and an expansion wave simultaneously reflects back into the tube (Figs. 1c and 1d). The expansion wave produces lower pressure gas in its downstream region, further expelling the combustion products situated behind the wave. When the expansion reaches the thrust wall or plate at the forward end of the thrust tube, it is reflected as an expansion wave (Fig. 1e). The pressure behind the reflected expansion is further reduced, below ambient, and air and fuel flow into the tube through their respective inlets during this ‘‘aspiration phase’’, as the expansion wave propagates toward the exit (Fig. 1f ). When the expansion wave reaches the exit (Fig. 1g), it is reflected as a compression wave, coalescing to a shock, from the open end and the wave propagates upstream toward the thrust wall (Fig. 1h). The reflection of the wave at the thrust wall delivers thrust to the vehicle and is reflected as a shock, allowing the combustible mixture of reactants to ignite, with transition to a propagating detonation wave (Figs. 1i and 1j) and a return to the initial condition of the cycle (Fig. 1a).
Numerical Simulation of Pulse Detonation Engine Phenomena
203
During propagation of the detonation, combustion of reactants takes place so rapidly that the process effectively takes place at a constant volume, which is theoretically more efficient than a constant pressure combustion process as essentially occurs in conventional airbreathing or rocket engines. The absence of rotating machinery suggests that the PDWE concept holds real promise for high thrust density, low fuel consumption propulsion applications, either in an airbreathing or rocket engine mode. A number of alternative configurations for the PDWE have been proposed and tested [4], and numerical simulations of PDWE phenomena are being actively pursued. Helman et al. [5] were the first to demonstrate a successful air-breathing detonation engine, with intermittent detonation frequencies of 25 Hz. Since then, extensive testing has been and continues to be pursued by a variety of groups [3, 4]. An excellent historical overview of the development of the PDWE as a constant volume combustion engine concept is provided in Eidelman et al. [1]. Overviews of past and ongoing numerical simulations of PDWEs are described in recent articles by Kailasanath [2, 6–8]. Kailasanath and Patnaik [9] explore the effect of employing a variable ‘‘pressure relaxation length’’ external to the end of the PDWE’s detonation tube. The relaxation length is used in order to mimic multidimensional wave effects in the pressure near the tube exit when using a 1D simulation, in particular when the local exit flow is subsonic over a portion of the cycle. Kailasanath and Patnaik find that the relaxation length can have a significant effect on engine performance, with increases in relaxation length leading to increases in engine impulse. The impulse I is typically defined as .
I — A F Dptw (t) dt
(1)
0
where A is the area of the thrust wall and Dptw (t) is the time-dependent pressure differential (compared with initial tube conditions) at the thrust wall. The engine’s specific impulse Isp may also be quantified as a performance parameter: Isp —
I rVg
(2)
where r is the initial mass of the gas mixture in the tube, V is the detonation tube volume, and g is the earth’s gravitational acceleration. An alternative performance parameter, the fuel-based specific impulse, Isp, f , is defined as Isp, f —
Isp Yf
(3)
where Yf is the fuel mass fraction present within the initial premixed reactants in the tube. The effect of partially filling the detonation tube with reactants, in the vicinity of the thrust wall, and having the remainder of the tube filled with air, is explored by Cambier and Tegner [10] and by Li et al. [11, 12], in both 1D and 2D simulations. Reducing the region of the tube occupied by reactants results in a higher Isp, f for a given engine cycle, clearly, due to the reduction in the effective Yf . This occurs
204
He and Karagozian
despite the fact that when the detonation propagates through the interface between reactants and air, reflected expansion waves return to the thrust wall and reduce the pressure at the wall at an earlier time than the reduction that occurs when the expansion waves reflected from the tube exit encounter the wall. An analytical model for the impulse of the single cycle PDWE has been developed by Wintenberger et al. [13]. Employing a similarity solution for the temporally constant pressure region within the tube (as occurs during the portion of the cycle shown in Figs. 1a–1e), with the subsequent ‘‘blowdown’’ process modeled using dimensional analysis and empirical relations, the model is able to impressively capture the essential physical phenomena present in the PDWE. Correspondence of the model results to performance parameters quantified in experiments and in numerical simulations is quite good. The present study focuses on both one- and two-dimensional simulations of the pulse detonation engine using high order numerical schemes, developed over the years by Stanley Osher and collaborators. We focus not only on exploring the PDWE’s reactive flow processes and performance parameters, but also on examining noise generation by the PDWE. NASA’s interest in the PDWE for advanced vehicle propulsion [14] is incumbent upon the ability of the engine to efficiently generate thrust without having to pay a significant penalty in engine noise. 2. PROBLEM FORMULATION AND NUMERICAL METHODOLOGY 2.1. Governing Equations The governing equations which describe the flow and reaction evolution in the PDWE with a single step, irreversible chemical reaction consist of the conservation of mass, momentum, energy, and species. The following representation was used, for example, in simulations of the two-dimensional planar PDWE configurations under investigation: F (U F )x +G F (U F )y =S F (U F) F t +F U
(4)
F , the flux vectors, F F and G F , and where the vector containing conserved variables, U F , are, respectively, the vector containing source terms, S
RS R S R S R S r
ru F = rv U
ru
ru +p
F (U F )= F
ruv
0
rv
2
0
ruv
F (U F )= G
rv 2+p
E
(E+p) u
(E+p) v
rY
ruY
rvY
F )= S(U
0 0
− KrYe −
12 Ti T
(5) and where E may be written r(u 2+v 2) p + +rqY E= c−1 2
(6)
Numerical Simulation of Pulse Detonation Engine Phenomena
205
Here r represents density, p is the static pressure, u and v are the x- and y-components of the velocity vector, respectively, and c is the ratio of specific heats. q is a heat release parameter which characterizes the amount of energy released during the reaction, and Ti represents the activation temperature. Y is the reactant mass fraction, which varies from 0 to 1, while K is the reaction-rate multiplier for the reaction, which also sets the spatial and temporal scales in the problem. One-step chemical reactions were considered in the present study, mostly for the methane-oxygen reaction, but also for the hydrogen-oxygen reaction. These were sufficient to enable reasonable estimates for performance parameters to be derived from the computations. The reaction terms K and Ti were determined for the methane-oxygen reaction from data provided in [15]. For the hydrogen-oxygen reaction, the STANJAN package was used to compute the equilibrium conditions downstream of an H 2 –O2 detonation, and K and Ti were fit to these data. We found that the PDWE cycle results were not sensitive to our selections of K and Ti , within reasonable limits, for this reaction, and hence these values were chosen to be the same as those for the methane-oxygen reaction. The values of these reaction quantities for the two reactions considered are shown in Table I. Two methods for initiation of the detonation reaction were explored. One method employed a computational ‘‘spark’’ adjacent to the thrust wall, by which a narrow high pressure, high temperature region (8 grid cells in width) initiated the propagating shock and ignited the reactants. The second method initiated the computation via an incoming shock moving through the reactants toward the thrust wall, which reflected from the wall and ignited the detonation. These initiation mechanisms produced virtually the same overall performance parameters, since the incoming shock strength in the second method could be chosen so that the pressure at the thrust wall rapidly rose to the same extent using either ignition method. As such, only results for the ‘‘spark’’ ignition of reactants are shown here. In most results shown, except as otherwise stated, the computational ‘‘spark’’ was represented by a region with an initial pressure of 20 atm and an initial temperature of 6000 K. At much lower spark region initial temperatures and pressures, of course, there would be a much more gradual transition from the deflagration (subsonic combustion) to the detonation process [16]. Both one- and two-dimensional numerical simulations were performed in this study. In addition to the 2D planar problem described in Eq. (4) above, corresponding equations for the axisymmetric PDWE configuration were solved, as were the equations for a 1D PDWE configuration. The percentage ‘‘fill’’ of the reactants
Table I. Values of Constant Terms in Relations (5) and (6) for Approximate Single Step Reactions Reaction
q (J/kg)
K
Ti (K)
CH 4 –O2 H 2 –O2
1.0 × 10 7 1.344 × 10 7
3.224 × 10 8 3.224 × 10 8
24358 24358
206
He and Karagozian
within the tube was varied in some of our computations, but results for the influence of this parameter are not shown in this paper since they basically replicated the trends described in Li and Kailasanath [12]. In the case of the 1D simulations, the use and effect of employing a pressure relaxation length exterior to the end of the detonation tube, after Kailasanath and Patnaik [9], was also explored. Incorporating a relaxation length pre-defines a spatial decay in the pressure immediately external to the PDWE tube after the exit of the detonation/shock (see Fig. 1c). When the relaxation length l was set to be greater than zero, the external pressure was allowed to spatially decay in a linear fashion over the length l to atmospheric pressure, rather than having the exit plane pressure exposed immediately to atmospheric pressure. The relaxation length assumption is suggested by Kailasanath and Patnaik [9] to enable the flow interior to the PDWE tube (i.e., the nature of the reflection of an expansion wave back into the detonation tube) to be more closely representative of the phenomena that actually take place in a multi-dimensional PDWE, particularly when exit conditions from the tube become locally subsonic and the interior flow can respond to changes outside of the tube. In the present 1D simulations, however, the computational domain consisted primarily of the detonation tube (containing 600 grid points), since the effect of a pressure relaxation length on the wave reflected back into the detonation tube could be represented using only a few grid points beyond the tube end in order to capture the spatial slope in pressure. In the 2D simulations, on the other hand, the air external to the detonation tube was assumed to be uniformly at atmospheric pressure, and the computational domain extended well downstream of the end of the tube, in general at least one and one half tube lengths downstream. The 2D domain extended at least two tube diameters away from the detonation tube in the dimension perpendicular to the axial dimension. Hence for the 2D simulations, no relaxation length approximation was necessary, and outflow boundary conditions were employed at the edges of the computational domain. Comparisons were thus possible between the 2D simulations, which captured flow characteristics external and internal to the PDWE, and 1D simulations with and without a pressure relaxation length, in which only the interior of the PDWE and a few exterior grid cells were resolved.
2.2. Numerical Methodology and Code Validation The present study used the Weighted Essentially Non-Oscillatory (WENO) method [17], a derivative of the Essentially Non-Oscillatory (ENO) method, introduced and developed by Stanley Osher and his colleagues [18–21] for spatial interpolation of the system of governing equations. ENO methods constitute a class of high accuracy, shock capturing numerical schemes for hyperbolic systems of conservation laws, based on upwind biased differencing in local characteristic fields. It has high accuracy (third order or higher) in smooth regions of the flow, and captures the motion of unresolved steep gradients in general without introducing spurious oscillations. Hence ENO is particularly well suited for resolution of
Numerical Simulation of Pulse Detonation Engine Phenomena
207
flowfields in which there are shocks or flame fronts. ENO uses an adaptive polynomial interpolation constructed on the basis of decisions to avoid steep gradients in the data; these decisions involve selection of the smoothest computational stencil in the vicinity of the gradient to approximate fluxes. The polynomial is also biased to extrapolate data from the direction of information propagation (i.e., upwind) for physical consistency and stability. WENO schemes, developed by Jiang and Shu [17], are based on ENO and incorporate all possible stencils in approximating fluxes, so that higher order accuracy may be achieved with a smaller number of points in a given stencil [22]. For example, with 3 points per stencil, an ENO scheme becomes third order in accuracy, but a WENO scheme becomes fifth order accurate in smooth regions and remains third order accurate in the vicinity of discontinuities. To avoid entropy-violating expansion shocks near sonic points, where characteristic velocities change sign, high order dissipation was added in the present study via the Local Lax Friedrichs (LLF) scheme [23], which adds extra numerical viscosity throughout the computational domain at each time step. Once the numerical approximation of the spatial terms was completed using WENO, the conservation equations (4) were written in the form F (U F) F t =f U
(7)
A third order TVD Runge–Kutta method was used for the time discretization [23], since the method has high accuracy and, more importantly, a large time step stability region which includes a segment of purely imaginary linear growth rates. This feature ensures that for a sufficiently small time step (according to the Courant–Friedrichs–Levy or CFL condition), the time discretization will not introduce spatial oscillations to the result. Since the current simulation only involved a single step reaction, use of the TVD Runge–Kutta method was preferred. The CFL condition in this case incorporated the source term to insure stability, since reaction rates often impose the greatest constraints on time integration. The ENO/WENO schemes were tested on a variety of problems, including that of the classical one-dimensional, overdriven, pulsating detonation [24]. It is found that the degree of spatial resolution of the reaction zone can have a profound effect on the peak pressure amplitude and period of oscillation of the overdriven pulsating detonation wave. A minimum required resolution of 20 grid points per detonation reaction zone half-length is suggested in [24] for the canonical pulsating detonation problem with an overdrive of 1.6; much finer grids may be required for other types of spatial integration schemes [25, 26]. Grid resolution tests in the present PDWE study, however, suggested that even 5 grid points per reaction zone half-length were sufficient for the cases under consideration. The present WENO scheme was also tested on a non-reactive problem with an analytic solution, that of the one-dimensional propagation of a normal shock wave out of a straight tube into air at atmospheric conditions. This problem yields reflection of an expansion fan back into the tube, mimicing the processes taking place when a detonation leaves the PDWE tube (see Figs. 1c and 1d). As with the one-dimensional PDWE simulations, the computational grid consisted of only the
208
He and Karagozian 250000
250000
Exact solution dx = 0.10 m dx = 0.02 m dx = 0.01 m
225000
200000
Pressure (pa)
Pressure (pa)
200000
175000
175000
150000
150000
125000
125000 100000
100000
75000
75000
50000
0
2
4
6
8
50000
10
0
2
4
6
X (m)
X (m)
(a)
(b)
250000
8
10
250000
Exact solution dx = 0.10 m dx = 0.02 m dx = 0.01 m
225000
225000
200000
Exact solution dx = 0.10 m dx = 0.02 m dx = 0.01 m
Pressure (pa)
Pressure (pa)
200000
175000
175000
150000
150000
125000
125000
100000
100000
75000
75000
50000
Exact solution dx = 0.10 m dx = 0.02 m dx = 0.01 m
225000
0
2
4
6
8
10
50000
0
2
4
6
X (m)
X (m)
(c)
(d)
8
10
Fig. 2. Validation of the 3rd order WENO scheme on a one-dimensional, non-reactive problem with an analytic solution. At time t=0, a diaphragm separating high pressure gas (phigh =2.0 × 10 5 Pa, rhigh =1.78 kg/m 3, x < 7 m) from lower pressure gas (plow =1.0 × 10 5 Pa, rlow =1.0 kg/m 3, x > 7 m) is burst. On either side of the diaphragm, c=1.2 and initially, T=300K. The tube end is situated at x=10 m, and for x > 10 m at all times, the pressure is 1 × 10 5 Pa. Shown are spatial distributions of pressure in the vicinity of the open tube end (within the tube) at various times: (a) t=0.8 ms, (b) t=4.0 ms, (c) t=12.0 ms, and (d) t=16.0 ms. Comparisons are made among simulation results for various spatial grid sizes (symbols) and the analytic or exact solution (lines). Propagation of the shock wave to the right ((a), (b)), out of the tube, results in reflection of an expansion wave back into the tube ((c), (d)). Here the flow throughout the domain shown remains subsonic (in the laboratory reference frame) for all times shown.
tube and a few grid cells exterior to or downstream of the tube. Figure 2 shows the computed pressure distribution in this non-reactive test problem at four different times. There was virtually no difference between the 3rd order WENO-computed solution and the analytic solution when the spatial grid spacing was less than about 2 cm in this test case. While ENO schemes can become expensive due to the required decisions on the stencils, WENO schemes were found to be less costly, and moreover, the ability to use relatively coarse grids, as seen in [24], helps to reduce
Numerical Simulation of Pulse Detonation Engine Phenomena
209
the overall computational time. Both the shock propagating out of the tube (right side of Figs. 2a and 2b) and the expansion wave reflected back into the tube (right side of Figs. 2c and 2d) were captured very well, even though only a few grid points outside of the end of the tube at atmospheric pressure were employed in the computation. The flow throughout the tube remained subsonic (or, prior to passage of the shock, quiescent) during the computation. This comparison lends confidence that the present scheme and boundary conditions, even in 1D, were capable of capturing important compressible wave dynamics at the end of the detonation tube. 2.3. Noise Estimates As described above and in Wintenberger et al. [13], for example, a variety of performance parameters may be derived from the spatial and temporal evolution of pressure and other flow properties in the PDWE. In addition to the impulse I, specific impulse Isp , and fuel specific impulse, Isp, f , the sound pressure level (SPL) at various locations within and external to the detonation tube were computed in this study. These noise levels were roughly estimated by examining the Fourier transform of the time-dependent pressure measured at various locations within the computational domain. A peak in the spectrum often appeared at roughly 330 Hz, which corresponded to the inverse of the period of the PDWE cycle for a tube length of one meter and a CH 4 –O2 reaction. Smaller peaks at higher harmonics were usually observed. The maximum pressure peak was used to estimate the sound pressure level or SPL at a given physical location, where the SPL is defined as SPL=20 log10
p pref
(8)
and where the value of pref is 2 × 10 −5 N/m 2 or Pascals (Pa). Hence the present computations allow rough estimates of the distribution of sound pressure level, or noise, emanating from pulse detonation engine to be made. 3. RESULTS The temporal evolution of flow and reaction structures within the pulse detonation engine were computed over a single cycle. An example of the temporal evolution of the computed 2D planar pressure field associated with the PDWE, for the CH 4 -O2 reaction, is shown in Fig. 3. A straight, planar PDWE tube one meter in length and 0.2 m in width was simulated. The tube was assumed here to be completely filled with a stoichiometric mixture of methane and oxygen at atmospheric pressure. In the 2D simulations, the grid spacing in the downstream dimension was Dx=1 cm or less, while the grid size in the transverse direction was Dy=0.3 cm. A CFL number of 0.1 was used in the time discretization. The propagation of the detonation to the right, within the tube, is readily apparent; its exit from the tube end resulted in the propagation of a vortical structure, coincident with the leading shock, into the surrounding air. Simultaneous to the exit of the shock and initiation of the vortex, a reflected expansion fan entered the detonation
210
He and Karagozian
Fig. 3. Temporal evolution of the 2D planar pressure field within and external to the PDWE over one cycle, with data shown at times corresponding to (a) 0.0 ms, (b) 0.06 ms, (c) 0.15 ms, (d) 0.26 ms, (e) 0.49 ms, (f ) 0.76 ms, (g) 1.01 ms, and (h) 2.89 ms. A CH 4 –O2 reaction was simulated here. The color map for the pressure distribution is shown, with units of Pascals.
Numerical Simulation of Pulse Detonation Engine Phenomena
211
Fig. 4. Temporal evolution of the centerline pressure within and external to the PDWE over one cycle, extracted from a 2D planar simulation of the flowfield, as in Fig. 3, with data shown at times corresponding to (a) 0.0 ms, (b) 0.06 ms, (c) 0.15 ms, (d) 0.26 ms, (e) 0.49 ms, (f ) 0.76 ms, (g) 1.01 ms, and (h) 2.89 ms.
212
He and Karagozian (e)
(a) .
(b) . Detonation tube
.
(f)
(c) .
.
(d) . Patm
Fig. 5. Locations within and exterior to the detonation tube (shown here without a nozzle) at which pressure data and SPL are computed. Here, along the centerline of the tube, (a) refers to the thrust wall, (b) refers to the tube midpoint, (c) refers to the exit plane of the tube, and (d) refers to a point one tube length external to the tube. Points (e) and (f ) lie 0.4 tube lengths above the tube centerline, aligned with the tube exit plane and at a location one tube length beyond the exit plane, respectively.
tube, whose subsequent reflection from the thrust wall resulted in nonlinear wave interactions evolving within the tube. The temporal evolution of the pressure distribution along the centerline of the tube for the same 2D simulation as in Fig. 3 is shown in Fig. 4. Here the initiation (Fig. 4a) and propagation of the detonation wave through the tube (Figs. 4b–4d), the exit of the shock from the tube (Figs. 4d and 4e), and the reflection of the expansion fan from the exit plane back into the tube (Figs. 4e–4g) are more obvious. The temporal evolution of pressure at various locations within and external to the PDWE may be derived from results such as those shown in Fig. 3. Data are presented here at a variety of points in the flowfield, labeled as shown in Fig. 5. For the 2D planar geometry, corresponding to the conditions in Fig. 3, the time-dependent pressures at the various locations (a)–(f ) identified in Fig. 5 are plotted in Figs. 6a–6f, respectively. A comparable computation but for the 2D axisymmetric detonation tube with a diameter of 0.2 m produces the results shown in Fig. 7. Interestingly, while the pressure evolution at the thrust wall and within the tube was not very different between the planar and axisymmetric configurations, the propagation and evolution of the shock structure exterior to the tube broadened over a longer time period for the axisymmetric configuration as compared with the planar configuration. In both planar and axisymmetric cases, the local pressure at point (c), external to the PDWE, was observed to drop to lowered density conditions for a short time (see Figs. 6e and 7e). This corresponded to the time at which the low pressure region associated with the inviscid, compressible vortex centers passed this point. The influence of these external flow differences between 2D and axisymmetric cases on performance parameters for the PDWE were very small, however, as seen by the plots of I, Isp , and Isp, f as a function of time shown in Fig. 8 for both planar and axisymmetric configurations. It is further noted that the time-dependent pressure at the thrust wall is of the same functional form as that observed in recent PDWE experiments performed at Stanford University [27] for a stoichiometric ethylene-oxygen mixture but with a 1.35 m long detonation tube.
Numerical Simulation of Pulse Detonation Engine Phenomena
213
Fig. 6. Temporal evolution of the pressure over one PDWE cycle, for points located at (a)–(f ) as labeled in Fig. 5. The 2D planar configuration with a CH 4 –O2 reaction was simulated here.
214
He and Karagozian
Fig. 7. Temporal evolution of the pressure over one PDWE cycle, for points located at (a)–(f ) as labeled in Fig. 5. The 2D axisymmetric computation with a CH 4 –O2 reaction was simulated here.
Numerical Simulation of Pulse Detonation Engine Phenomena
215
Fig. 8. Performance parameters for the 2D planar and 2D axisymmetric simulations of the detonation tube: (a) Impulse I per unit area, (b) specific impulse Isp , and (c) fuel specific impulse Isp, f as a function of time. A CH 4 –O2 reaction was simulated here.
Time-series pressure data such as those in Figs. 6 or 7 may be used to estimate the noise generated at various points in the flowfield over a single PDWE cycle. Pressure spectra for the various locations within and external to an axisymmetric PDWE tube are shown in Fig. 9, with correspondence to the points labeled in Fig. 5 and the data in Fig. 7. In all locations we observe the highest peak in pressure to appear close to the frequency associated with the period of the PDWE cycle, roughly 330 Hz for these conditions. Clearly, the numerical shocks, with an artificial thickness, contributed at least some of the higher frequency peaks in these spectra, but as a rough approximation, the SPL may be computed from the lowest frequency’s peak pressure, at 330 Hz. Hence an approximation for the noise levels within and about the detonation tube may be quantified, shown in Table II. While clearly the highest noise levels may be found within the detonation tube, it is interesting to note that the noise level even one tube length downstream of the
216
He and Karagozian
Fig. 9. Sample pressure spectra at points labeled as shown in Fig. 5 for the 2D axisymmetric PDWE tube with a CH 4 –O2 reaction.
Numerical Simulation of Pulse Detonation Engine Phenomena
217
Table II. Computed Sound Pressure Level at Flowfield Locations Labeled in Fig. 5 Location
a
b
c
d
e
f
SPL (dB)
212
211
202
173
174
173
exit plane of the tube may be as high as 173 dB. These data were generally consistent with noise measurements in PDWE tests performed at NASA Glenn Research Center [28], although since the present simulation was inviscid, dissipative phenomena which could reduce external noise levels were not represented here. The decay in SPL from 202 dB at the tube exit (c) to 173 dB one tube length downstream (d) reflects the spatial dispersion of the exiting shock wave. Reduction in the shock/detonation strength within the tube could significantly reduce the noise generation by the PDWE, but of course this has the disadvantage of reducing engine performance. This tradeoff between performance and noise, and the effects of the higher frequencies shown in Fig. 7, will be explored in future studies. To this point, all results shown have been for 2D simulations, either planar or axisymmetric. A significant savings in computational time may be derived if performance parameters and noise estimates are computed from 1D simulations. In performing the 1D simulations, we first explored the effect of employing a pressure relaxation length l, as utilized in the computations of Kailasanath and Patnaik [9]. As soon as the detonation/shock proceeded from the PDWE tube exit in the present computations, the pressure immediately external to the tube was assumed to ‘‘relax’’ linearly between its time-dependent value at the exit to atmospheric pressure at a distance I from the exit. Only a few ‘‘ghost points’’ downstream of the exit were employed in the 1D simulation, so that specification of l and thus the spatial pressure gradient were sufficient to impose this exit condition. Increasing l, from zero up to half a PDWE tube length L, was observed to increase the PDWE’s impulse, shown in Fig. 10 and also as seen in [9]. The reason for this increase in PDWE impulse may be discerned by examining the time-dependent pressure at the thrust wall and at the tube exit, shown in Figs. 11a and 11b, respectively, for cases where the pressure relaxation lengths are l=0 and l=0.5L. While the pressure at the thrust wall was the same for both cases until about 1 msec (the time at which the reflected expansion wave reached the thrust wall), after this time the thrust wall pressure decayed more rapidly when l=0 than when l=0.5L. A similar phenomenon was observed at the tube exit. The pressure relaxation length forced a more gradual spatial drop in the pressure outside of the tube, causing the pressure at the exit plane to exceed that of the case where l=0. In fact, in some cases, incorporating a non-zero pressure relaxation length caused the flow at the tube exit to remain subsonic when it would have become supersonic if l had been set to zero, as will be seen below. To understand if similar PDWE data may be derived from 1D simulations as for 2D simulations, comparisons for the 2D axisymmetric detonation tube with a
218
He and Karagozian
Impulse per area (pa.s)
2000
1500
relaxation length = 0 relaxation length = 0.01 relaxation length = 0.02 relaxation length = 0.10 relaxation length = 0.30 relaxation length = 0.50
1000
L L L L L
500
0
0
0.001
0.002
0.003
0.004
Time (sec.)
Fig. 10. Comparisons of time-dependent PDWE impulse per unit area from 1D simulations using different pressure relaxation lengths, given in terms of the detonation tube length L (1 m). A methane-oxygen reaction is simulated.
methane-oxygen reaction are made in Fig. 12. Here the 1D simulations were conducted for l=0 and l=0.5L, and comparisons were made with the 2D axisymmetric computations by altering the initial temperature and pressure of the computational ‘‘spark’’ in the 1D simulation so that equivalent initial conditions were attained in all three cases. The spatial pressure evolution shows that, prior to the detonation wave’s leaving the tube (Figs. 12a–12d), all three sets of simulations yielded nearly identical results, which is appropriate since the flow within the tube
(a)
40
(b)
35
35 30
relaxation length = 0 relaxation length = 0.50 L
25
Pressure (atm)
Pressure (atm)
30
20
15
20 15
10
10
5
5
0
0
0.001
0.002
Time (sec.)
0.003
0.004
relaxation length = 0 relaxation length = 0.50 L
25
0
0
0.001
0.002
0.003
0.004
Time (sec.)
Fig. 11. Temporal evolution of the pressure over one PDWE cycle, for points located at (a) the thrust wall and (b) the PDWE tube exit. Here the 1D simulation is employed, for the cases of zero relaxation length (l=0) and a relaxation length l corresponding to half a PDWE tube length. A CH 4 –O2 reaction was simulated here.
Numerical Simulation of Pulse Detonation Engine Phenomena
219
was essentially one-dimensional for a single-step reaction, and the boundary conditions at and just beyond the tube exit were initially identical. But when the detonation left the tube, becoming a shock, the expansion fan reflected back into the tube had characteristics that are different between the two different 1D simulations, depending on the boundary specifications at the tube exit and hence on the assumption of a pressure relaxation length. If the pressure just downstream of the exit was forced to spatially decay slowly away from the exit for subsonic outflow, as with the l=0.5L case, rather than abruptly to atmospheric pressure, as in the 1D simulation without a relaxation length, the pressure immediately external to the
Fig. 12. Spatial distribution of the centerline pressure at different times over one PDWE cycle, with comparisons among the 2D axisymmetric simulations (solid line), 1D simulations without the inclusion of a pressure relaxation length (dashed line), and 1D simulations with a pressure relaxation length of half the PDWE tube length (l=0.5L) (dashed-dotted line). A methane-oxygen reaction was simulated. Distributions shown are for approximately the same times: (a) 0.00 ms, (b) 0.06 ms, (c) 0.150 ms, (d) 0.260 ms, (e) 0.370 ms, ( f ) 0.490 ms, (g) 0.760 ms, (h) 1.01 ms, and (i) 2.88 ms.
220
He and Karagozian
tube was forced to become higher for l=0.5L than for l=0. This caused the reflected expansion fan to become relatively weak for l=0.5L as compared with the 2D simulation. The abrupt pressure change at the exit for the 1D simulation with l=0 was actually more severe than for the 2D simulation (see Fig. 12e), but was closer to the 2D results over time (see Figs. 12g–12i). As noted above, the imposition of the pressure relaxation length when the flow at the exit plane was locally subsonic caused the flow at the tube exit to remain locally subsonic after the detonation/shock left the tube, while the case where l=0 in fact produced slightly supersonic exit flow in a short time (see Fig. 13), as did the 2D simulation. While neither 1D simulation precisely replicated the 2D axisymmetric flow’s Mach number solution at the tube exit resulting from the propagating shock/vortical structure, for the operating conditions examined, the l=0 case was able to capture the flow somewhat better. As a consequence, performance parameters were predicted by the 1D simulation with l=0 that were closer to the 2D results, as seen in Fig. 14. This type of comparison was found for many relaxation lengths considered here, however, different initial conditions, tube lengths, and reaction mechanisms could render a non-zero relaxation length more realistic, as seen by Kailasanath and Patnaik [9]. Hence a great deal of care is required in employing the pressure relaxation length concept in order to replicate multidimensional wave phenomena.
2
2D simulation 1D without pressure relaxation 1D with pressure relaxation
Mach Number
1.5
1
0.5
0
0
0.001
0.002
0.003
Time (sec.)
Fig. 13. Comparisons of the Mach number at the PDWE tube exit plane for a methane-oxygen reaction. Results are shown for the 2D axisymmetric PDWE simulations, and for the 1D simulations with a pressure relaxation length of half a tube length (l=0.5L) and with no pressure relaxation length (l=0). Here the 1D simulations incorporated a computational ‘‘spark’’ consisting of a pressure of 10 atm and a temperature of 3000 K in order to match the initial conditions for the 2D simulation.
Numerical Simulation of Pulse Detonation Engine Phenomena
221
Fig. 14. Comparisons of (a) time-dependent pressure at the thrust wall, (b) time-dependent impulse per unit area, and (c) time-dependent specific impulse for both 1D and 2D axisymmetric simulations of the PDWE tube with a CH 4 –O2 reaction. 1D simulations explored the use of a pressure relaxation length l=0.5L. Here the 1D simulations incorporated a computational ‘‘spark’’ consisting of a pressure of 10 atm and a temperature of 3000 K in order to match the initial conditions for the 2D simulation.
Finally, it is noted that when an alternative combustion reaction was simulated, that of hydrogen-oxygen, there was a quantifiable increase in specific impulse as compared with methane-oxygen. This faster H 2 –O2 reaction had a higher q value in the reaction source term (Eq. (6)) and as such generated higher temperatures and pressures behind the detonation, but the H 2 –O2 reaction also had a lower mixture mass than did the CH 4 –O2 reaction. There was thus a correspondingly quantifiable increase in specific impulse, as compared with methane-oxygen, as shown in Fig. 15 for the 1D computation without employment of a relaxation length.
222
He and Karagozian
Specific impulse (sec.)
200
150
CH 4 + O 2 H 2 + O2
100
50
0
0
0.001
0.002
0.003
Time (sec.)
Fig. 15. Comparisons of time-dependent PDWE specific impulse computed from 1D simulations using different single step reaction mechanisms, that for hydrogen and oxygen, and that for methane and oxygen. A zero pressure relaxation length was employed. Here the 1D simulations each incorporated a computational ‘‘spark’’ consisting of a pressure of 20 atm and a temperature of 6000 K.
4. CONCLUSIONS The present simulations suggest that a great deal of useful performance and noise related data may be obtained even from one-dimensional computations of the pulse detonation wave engine with simplified reaction kinetics. Third order WENO schemes were observed to provide robust, accurate simulations of the reactive and non-reactive flow processes occurring within the PDWE. 1D and 2D simulations were observed to produce similar pressure profiles throughout the PDWE tube, especially in the vicinity of the thrust wall, and hence similar engine impulse values, without the inherent need for inclusion of a pressure relaxation length for the conditions examined here. Alternative reaction mechanisms, tube lengths, and initial conditions may require such relaxation, however. Straightforward estimates of sound pressure levels within and external to the engine tube may also be made, from both 2D and 1D simulations, although the generation and effects of higher frequency overtones requires further examination. Future studies will explore these effects on the acoustic field generated by the pulse detonation engine, in addition to the effects of the presence of a nozzle on noise generation. ACKNOWLEDGMENTS The authors wish to extend sincere congratulations to Professor Stanley Osher on the occasion of his 60th Birthday and thanks for a number of years of collegial collaboration at UCLA. The authors also wish to thank Professor Chi-Wang Shu of Brown University for his helpful advice in speeding the WENO computations. This work has been supported at UCLA by NASA Dryden Flight Research Center
Numerical Simulation of Pulse Detonation Engine Phenomena
223
under Grant NCC4-153, with Dr. Trong Bui and Dave Lux as technical monitors, and by the Office of Naval Research under Grant ONR N00014-97-1-0027, with Dr. Wen Masters as technical monitor.
REFERENCES 1. Eidelman, S., Grossmann, W., and Lottati, I. (1991). Review of propulsion applications and numerical simulations of the pulsed detonation engine concept. J. Propulsion and Power 7(6), 857–865. 2. Kailasanath, K. (2002). Recent Developments in the Research on Pulse Detonation Engines, AIAA Paper 2002-0470 (Invited), AIAA 40th Aerospace Sciences Meeting. 3. Roy, G. D. (2001). Practical Pulse Detonation Engines–How Far Are They? ISABE Paper 20011170. 4. Pulse-Detonation Engine (PDE) Workshop, sponsored by the Office of Naval Research, Naval Postgraduate School, Monterey, CA, 1997. 5. Helman, D., Shreeve, R. P., and Eidelman, S. (1986). Detonation Pulse Engine, AIAA Paper 86-1683. 6. Kailasanath, K., Patnaik, G., and Li, C. (2002). The Flowfield and Performance of Pulse Detonation Engines, to appear in the Proc. of the 29th Symposium (Intl.) on Combustion. 7. Kailasanath, K. (2001). A Review of PDE Research – Performance Estimates, AIAA Paper 20010474, AIAA 39th Aerospace Sciences Meeting. 8. Kailasanath, K. (2000). Review of propulsion applications of detonation waves. AIAA Journal 38(9), 1698–1708. 9. Kailasanath, K., and Patnaik, G. (2000). Performance Estimates of Pulsed Detonation Engines, Proc. of the 28th Symposium (Intl.) on Combustion. 10. Cambier, J.-L., and Tegner, J. K. (1998). Strategies for pulsed detonation engine performance optimization. Journal of Propulsion and Power 14(4), 489–498. 11. Li, C., Kailasanath, K., and Patnaik, G. (2000). A Numerical Study of Flow Field Evolution in a Pulsed Detonation Engine, AIAA Paper 2000-0314, AIAA 38th Aerospace Sciences Meeting. 12. Li, C., and Kailasanath, K. (2001). A Numerical Study of Reactive Flows in Pulsed Detonation Engines, AIAA Paper 2001-3933, 37th AIAA/ASME/SAE/ASEE Joint Propulsion Conference. 13. Wintenberger, E., Austin, J. M., Cooper, M., Jackson, S., and Shepherd, J. E. (2001). An Analytical Model for the Impulse of a Single Cycle Pulse Detonation Engine, AIAA Paper 2001-3811, 37th AIAA/ASME/SAE/ASEE Joint Propulsion Conference. 14. Three Pillars for Success: NASA’s Response to Achieve the National Priorities in Aeronautics and Space Transportation, NASA Office of Aeronautics and Space Transportation Technology, 1997. 15. Turns, S. (2000). An Introduction to Combustion, 2nd Ed., McGraw–Hill. 16. Kuznetsov, M. S., Alekseev, V. I., and Dorofeev, S. B. (2000). Comparison of critical conditions for DDT in regular and irregular cellular detonation systems. Shock Waves 10(3), 217–223. 17. Jiang, G. S., and Shu, C. W. (1996). Efficient implementation of weighted ENO schemes. J. Comput. Phys. 126, 202–228. 18. Harten, A., Osher S. J., Engquist, B. E., and Chakravarthy, S. R. (1986). Some results on uniformly high-order accurate essentially nonoscillatory schemes. J. Appl. Numer. Math. 2, 347–377. 19. Shu, C. W., and Osher, S. (1989). Efficient implementation of essentially non-oscillatory shock capturing schemes II. J. Comput. Phys. 83, 32–78. 20. Fedkiw, R. P., Merriman, B., Donat, R., and Osher, S. (1998). The penultimate scheme for systems of conservation laws: Finite difference ENO with marquina’s flux splitting. In Hafez, M. (ed.), Progress in Numerical Solutions of Partial Differential Equations, Arachon, France. 21. Fedkiw, R. P., Merriman, B., and Osher, S. (1997). High accuracy numerical methods for thermally perfect gas flows with chemistry. J. Comput. Phys. 132(2), 175–190.
224
He and Karagozian
22. Lee, T. K., and Zhong, X.-L. (1999). Spurious numerical oscillations in simulation of supersonic flows using shock-capturing schemes. AIAA Journal 37(3), 313–319. 23. Lahey, C. B. (1998). Computational Gasdynamics, Cambridge University Press. 24. Hwang, P., Fedkiw, R., Merriman, B., Karagozian, A. R., and Osher, S. J. (2000). Numerical resolution of pulsating detonation waves. Combust. Theory Model. 4(3), 217–240. 25. Bourlioux, A., Majda, A. J., and Roytburd, V. (1991). Theoretical and numerical structure for unstable one-dimensional detonations. SIAM J. Appl. Math. 51, 303–343. 26. Papalexandris, M. V., Leonard, A., and Dimotakis, P. E. (1997). Unsplit schemes for hyperbolic conservation laws with source terms in one space dimension, J. Comp. Phys. 134(1), 31–61. 27. Jenkins, T. P., Sanders, S. T., Kailasanath, K., Li, C., and Hanson, R. D. (2000). Diode-Laser Based Measurements for Model Validation in Pulse Detonation Flows, Proc. of the 25th JANNAF Airbreathing Propulsion Meeting, Monterey, CA. 28. Perkins, H. D., private communication.
