- Email: [email protected]
Detection of shock and detonation wave propagation by cross correlation
ARTICLE IN PRESS Mechanical Systems and Signal Processing 23 (2009) 1098–1111
Contents lists available at ScienceDirect
Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/jnlabr/ymssp
Detection of shock and detonation wave propagation by cross correlation F.K. Lu a,, A.A. Ortiz a, J.-M. Li b,1, C.H. Kim c, K.-M. Chung d a
Aerodynamics Research Center, Mechanical and Aerospace Engineering Department, University of Texas at Arlington, 214B Woolf Hall, 500 West 1st Street, Arlington, TX 76019-0018, USA Institute of Aeronautics and Astronautics, National Cheng Kung University, Tainan, Taiwan, Republic of China c Mechanical Engineering Department, Purdue University, West Lafayette, IN 47907, USA d Aerospace Science and Technology Research Center, National Cheng Kung University, Tainan, Taiwan, Republic of China b
a r t i c l e i n f o
abstract
Article history: Received 28 August 2007 Received in revised form 3 November 2008 Accepted 10 November 2008 Available online 14 November 2008
The propagation speed of a shock or detonation wave in a shock tube is usually determined by a time-of-flight method by dividing the distance between two transducers with the propagation time of the disturbance signal. Some arbitrariness is inherent in determining the propagation time by this method. An improved method for objectively determining the propagation time using a nonstationary cross-correlation technique is described. The method requires the choice of an integration window that includes the nonstationary event. The method was first tested against a number of model functions with different noise levels. It was then applied to propagating and reflected shock and detonation waves, including an example of a transitioning detonation wave propagating past six transducers. In addition, the nonstationary CCF technique was also applied to evaluate the uncertainty in estimating the deflagrationto-detonation transition run-up distance. In all cases, the time delay and its standard deviation could be obtained. & 2008 Elsevier Ltd. All rights reserved.
Keywords: Detonation Shock Shock tube Wave propagation
1. Introduction The determination of the propagation speed of a shock or detonation wave as in shock or detonation tube experiments typically makes use of the time-of-flight (TOF) method. In this method, the speed of a propagating shock, say, is obtained by determining the passage time of the wave DT between two transducers with a known separation distance DL, to yield u^ s ¼ DL=DT .
(1)
Some recent examples of the TOF method for shock and detonation wave studies include [1–3] and for other propagating signals include [4,5]. The TOF method appears to be suitable for strong, distinct transient signals such as shock and detonation waves. In situations where the signal is weak, the TOF method requires further treatment or may even have to be replaced by more powerful methods. Corresponding author. Tel.: +1 817 272 2083; fax: +1 817 272 5010.
E-mail addresses: [email protected] (F.K. Lu), [email protected] (A.A. Ortiz), [email protected] (J.-M. Li), [email protected] (C.H. Kim), [email protected] (K.-M. Chung). 1 Present address: Temasek Laboratories, National University of Singapore, Singapore. 0888-3270/$ - see front matter & 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.ymssp.2008.11.001
ARTICLE IN PRESS F.K. Lu et al. / Mechanical Systems and Signal Processing 23 (2009) 1098–1111
1099
Even though the TOF method is valid in principle and is simple to apply to propagating shock or detonation waves, its practical implementation can be subjective. This is because the shock front is not measured as a step rise, as in theory, but is spread over a small time interval. This interval depends on a number of factors, such as the size of the transducer, viscous effects and the data-sampling rate. Different arbitrary criteria can, thus, be formulated for determining the shock arrival over a transducer. For example, the shock may be deemed to have arrived at the transducer when the signal rises above an arbitrary threshold, or when the signal is midway between arbitrarily defined upstream and downstream values. These criteria, which yield signal arrival times close to one another, are subjective. They may be unreliable if the error of the time estimate is comparable to DT. Such a situation may occur if the transducers are closely spaced or if the sampling rate is low. For the same record length, increasing the sampling rate does not fundamentally overcome this difficulty although it may reduce the uncertainty in the DT estimate. To illustrate this arbitrariness, consider the propagation of a fully-developed Chapman–Jouguet (CJ) detonation wave past two transducers in a stoichiometric oxyhydrogen mixture initially at 2.01 atm, Fig. 1, the data being reported previously in [6]. (Fig. 1(a) includes other features, including a subsequent reflection that will be described later.) The region where the pressure initially rises is enlarged and shown in Fig. 1b. The symbols represent actual data, sampled at 100 kHz per channel with a simultaneous sample-and-hold system. Fig. 1b shows that the two wave fronts X(t) and Y(t) are spread over a few data points. The double-headed arrows indicate possible pairs of data points for the TOF method, resulting in estimates of DT from 0.21 to 0.23 ms. The 0.02 ms uncertainty represents a 10% error on the estimate of the velocity. A cross-correlation approach has been developed to improve the estimate of DT. This method also ensures an objective ˆ s. Moreover, the method is accurate even for noisy signals. However, since the wave method for estimating the error on u signals are not random, stationary processes [7], an extension of the usual correlation analysis, which assumes data stationarity, was needed. This nonstationary cross-correlation approach was developed and tested against a number of model functions before being applied to field signals. It can be noted that the cross-correlation technique has been widely used in various disciplines for detecting propagating signals, especially if they are weak, under the general method of timedelay estimation [8–14].
Fig. 1. Arbitrariness of the TOF method: (a) Detonation wave propagating past two transducers and (b) enlargement of wave fronts.
ARTICLE IN PRESS 1100
F.K. Lu et al. / Mechanical Systems and Signal Processing 23 (2009) 1098–1111
1.1. The cross-correlation function (CCF) For an infinitely long record of a random, stationary process, the mean value is constant, that is, it is independent of t [7]. Consider two random, stationary processes x(t) and Z(t), both with zero means. The CCF between these two processes is defined as Z 1 T xðtÞZðt þ tÞ dt. (2) RxZ ðtÞ ¼ lim T!1 T T If x(t) and Z(t) are totally uncorrelated, then RxZ(t) will be zero for all values of t. If x(t) and Z(t) are correlated, RxZ(t) will display a peak at some value of t. If the peak occurs at a positive value of t, the correlated event in Z(t) occurs after an event in x(t), and vice versa if the peak occurs at a negative value. In practice, data cannot be obtained over an infinite time. Eq. (2) is then modified for two data records of the same finite record length T to yield an estimate of the CCF Z T=2 1 xðtÞZðt þ tÞ dt. (3) R^ xZ ðtÞ ¼
T
T=2
The (4) will be dropped in the subsequent discussion. For digital implementations, the above equation is discretized to yield RxZ ðtj Þ ¼
N 1X xðtk ÞZðtk þ tj Þ. N k¼1
(4)
If there is a correlation between the two signals x(t) and Z(t), then the peak value of RxZ(t) occurs at tM, which represents the best correlation between x(t) and Z(t). If x(t) and Z(t) are from a propagating signal, tM would be the time for the signal to propagate past the two respective measurement locations. Together with the displacement between the two locations, the propagation speed can be deduced from Eq. (1), the same equation used for the TOF method. Thus, for the present application, the cross-correlation approach can be regarded as a sophisticated extension of the TOF method. An important advantage of the cross-correlation method over the TOF method is that it can provide an objective means of determining the propagation time. The resulting error in estimating the propagation speed is also quantifiable. 1.2. Nonstationary CCF The above discussion holds for random data which can be assumed to be stationary although there are many classes of data where stationarity does not hold. However, the cross-correlation technique may be extended to nonstationary data. Let the record from a transient process be expressed as XðtÞ ¼ xðtÞ þ nðtÞ,
(5)
where x(t) is a deterministic function and n(t) is a random ‘‘noise’’ with a stationary mean value of zero. Before applying a nonstationary cross-correlation, it is necessary to remove the noise as much as possible. If the variation of x(t) is slow compared to the lowest frequency in n(t), then x(t) can be separated from n(t) by low-pass filtering operations on X(t), However, low-pass filtering cannot be applied to shock wave data. The technique suggested in [7] is to use segmented mean value estimates by short-time-averaging. In this approach, the original record is divided into segments and an average value of each segment is used. The effect is a localized low-pass filtering at a frequency equal to the reciprocal of the averaging time t. The selection of t is arrived by trial-and-error. As t becomes short, the random error from n(t) increases ˆ(t) departs from the ‘‘true’’ x(t). Thus, the selection of an appropriate averaging time but as t becomes long, the estimated x involves a compromise between random and bias errors. A suitable value of t for the present study is three times the sampling time. As can be seen from Fig. 1, the shock is detected within 3–5 times of the sampling interval. A different value of t may be more appropriate for other situations. The possibility for such a modification is incorporated within the algorithm (block diagram in Appendix A). The procedure for performing the nonstationary CCF computation requires one more step than that for stationary data. Eq. (4) indicates that the CCF for a stationary process is a function of t only. However, for a nonstationary process, the CCF also depends on the time window, that is, it is also a function of t. This is expressed as a dependence on ti in general, as Rxy ðti ; tj Þ ¼
N i X 1 xðtiþk Þyðtiþk þ tj Þ. N i þ 1 k¼1
Moreover, there is a number of restrictions on ti since (1) x(ti+k) and y(ti+k+tj) must contain the nonstationarity event within k ¼ 0-Ni, (2) ti4t1, and (3) ti+kotN.
(6)
ARTICLE IN PRESS F.K. Lu et al. / Mechanical Systems and Signal Processing 23 (2009) 1098–1111
1101
Condition (1) means that the nonstationarity of interest (specifically a shock or detonation wave) must exist between ti and tN. If N datapoints can be selected from the entire range of each record, the value N constitutes another parameter as required in Eq. (6), since tN itself can take on a range of values as long as conditions (1) and (3) are satisfied. To satisfy the above conditions, therefore, Eq. (6) can be formulated as Rxy ðT L ; T U ; tj Þ ¼
U X 1 xðtk Þyðtk þ tj Þ, U L þ 1 k¼L
(7)
where L and U are the lower and upper limits (LoU) of the CCF calculation. Thus, the nonstationary event (of both records) must exist between TL and TU. Finally, since the interest here is the peak correlation time tM (corresponding to the propagation time between two points), the correlation function can be normalized against the peak value Rxy(TL,TU,tM) rxy ðT L ; T U ; tÞ ¼
1 Rxy ðT L ; T U ; tM ÞT
Z
TL
xðtÞyðt þ tÞ dt,
(8)
TL
where T ¼ T U T L . 2. Algorithm Algorithms were developed in FORTRAN, MatlabTM and LabVIEW TM to extract the propagation time of a shock or detonation wave past pairs of transducers.2 The FORTRAN program, known as XCOR, will be described. There is practically no difference in functionality between the three implementations. As seen in Eq. (8), the nonstationary CCF is a function of the time delay, and of the lower and upper limits of the sample window. The time delay t can be varied from T to þTbut the lower and upper limits of the integration window T must be carefully determined according to the characteristics of the given data. Also, the averaging time t should be determined to appropriately extract the nonstationary deterministic component from the given nonstationary random data. XCOR has the routines to perform these functions and it is described in more detail in the Appendix. It can also be noted here that, although the examples involve nonstationary pressures, the algorithm is a strictly numerical procedure and can be applied to pairs of nonstationary data in general. 3. Model calculation Before applying XCOR on experimental data, it was tested against a number of model functions with pseudo-random noise. These are
unit step functions X(t) and Y(t), both with SNR ¼ 2–32 and with Y(t) delayed by tM ¼ 50; a unit step function X(t) at t ¼ 0 with SNR ¼ 16, and a step function Y(t) at tM ¼ 50 with an amplitude of A ¼ 22–23 and SNR ¼ 16;
exponentially decaying functions X(t) ¼ ebt/100 and Y ðtÞ ¼ ebðtt
M Þ=100
, where b ¼ 22–22 and SNR ¼ 16 for both
functions, and tM ¼ 50 and
linearly decaying functions X(t) ¼ 1ct and X(t) ¼ 1c(ttM), where c ¼ 22–22 and SNR ¼ 16 for both functions, and tM ¼ 50. The SNR is the inverse of the rms noise level with respect to unity. Examples of these functions are shown in Fig. 2, while their cross correlations obtained using XCOR are shown in Fig. 3. For all models, the number of data points N is 1001 and the interval between the data points Dt is unity. The time delay tM between X(t) and Y(t) was fixed at 50 for all these cases. Thus, this exercise checks the ability of XCOR to compute the maximum of rxy(t) at this delay. Since these tests were to examine the dependence of the CCF upon data characteristics, the routine to vary the limits of the integration window [TL,TU] was not included for expediency. The data displayed in Fig. 3(a) indicate that the time delay estimated from the cross correlations for model 1 was not seriously affected by noise if the SNR was larger than 4. Thus, the SNR was fixed as 16 for the other 3 model functions. The amplitude ratio A of the second step function was varied from 0.25 to 8 in model 2 to further test the algorithm’s ability to detect the propagation of a signal subjected to attenuation or amplification. For the model profiles shown in Fig. 2(b), the computations obtained correlation peaks at t ¼ 50 as seen in Fig. 3(b). Therefore, the ability of the routine for detecting peak correlation is unaffected by variations in the amplitude. Based on the results of models 1 and 2, the SNR and amplitude were set as 16 and 1, respectively, for models 3 and 4. Model 3 simulates a detonation wave with Taylor rarefaction. The value of b for the exponential decay varied from 0.25 to 4. For all the cases in this model, the peak of the cross correlation at t ¼ tM ¼ 50 was well captured as indicated in Fig. 3(c). (For comparison, the CCF of model 1 is also plotted.) It can be seen that as the profile decays faster, the CCF becomes 2
These algorithms are available from the corresponding author or from the website http://arc.uta.edu/.
ARTICLE IN PRESS 1102
F.K. Lu et al. / Mechanical Systems and Signal Processing 23 (2009) 1098–1111
Fig. 2. Model functions. (a) Model 1: unit step function with SNR ¼ 16. (b) Model 2: variable step function with SNR ¼ 16. (c) Model 3: exponential decay; X(t) only, SNR ¼ 16. (d) Model 4: linear decay; X(t) only, SNR ¼ 16.
sharper. Finally, linear decaying profiles are shown in Fig. 2(d), with c varying from 0.25 to 4. The CCF has a peak at a time delay at t ¼ tM ¼ 50 as well, as can be seen in Fig. 3(d). From the model calculations, it can be concluded that the computed value of tM ¼ 50 is not a function of profile shape or amplitude of the transient signal. However, it might be affected by excessive noise. These tests also indicate that a sharp signal such as that from a detonation wave with Taylor rarefaction enables tM to be more distinctly resolved. Moreover, the wave system in a shock tube may involve the passage of two nearby events, such as an incident and a reflected shock. Therefore, the determination of an integration window [TL,TU] is even more important in separating out the two events for further processing. The functionality of XCOR was thus further tested using a model with two peaks to simulate the propagation and reflection of a detonation wave. Fig. 4 shows the profiles of model 5, which has a maximum pressure level of unity and an exponential decay with b ¼ 0.01 for both the incident and reflected waves. The propagation time of the incident and reflected waves were set to 20 and 40, respectively, and the SNRs were N, 16 and 8. From these profiles, tM was derived using XCOR. The absolute values of tM as a function of TL and TU are shown in Fig. 5. For the incident wave, tM is positive, but for the reflected wave, tM is negative, since Rxy(t) ¼ Ryx(tM). When SNR ¼ N, tM occurs on the tM ¼ 20 plane for the incident wave and on the tM ¼ 40 plane for the reflected wave. (Such flat surfaces are hereafter known as the ideal tM-surface.) However, when there is noise, tM deviates from the ideal tM-surface as can be seen in Fig. 5(b) and (c). Such departures were eliminated by the algorithm using pre-determined criteria in order to compute the average and standard deviation of tM. These computed values are listed in Table 1. It can be seen that noise caused a poorer estimate of the propagation time by both the peak-to-peak TOF method and the CCF method. While the TOF method apparently gave a closer estimate to the true propagation delay, the arbitrariness of the method remains and an estimate of the standard deviation may be difficult. The table shows that the propagation time for the incident wave was well captured when the SNR is larger than 8. However, for the reflected wave, the averaged value of tM was more difficult to obtain since the incident wave may not be clearly excluded from the integration window. Nevertheless, the nonstationary CCF technique can provide an objective method for obtaining the standard deviation in the wave propagation time. The common TOF method does not appear capable of achieving this. The nonstationary CCF technique is next demonstrated by applying it to 4 sets of experimental data.
ARTICLE IN PRESS F.K. Lu et al. / Mechanical Systems and Signal Processing 23 (2009) 1098–1111
1103
Fig. 3. Cross correlations of model functions. (a) Model 1: unit step function with SNR ¼ 16. (b) Model 2: variable step function with SNR ¼ 16. (c) Model 3: exponential decay; X(t) only, SNR ¼ 16. (d) Model 4: linear decay; X(t) only, SNR ¼ 16.
4. Application of nonstationary CCF 4.1. Shock wave propagation Fig. 6a shows voltages from two pressure transducers, where the step rise in the signal indicates the passage of a shock wave [15]. XCOR was applied to the wave and the resulting tM(TL,TU) is shown in Fig. 6(b). It can be seen that the tM surface is relatively flat except for the steep rise when the upper limits of the integration window approach the time of shock arrival. The average and the standard deviation of tM were computed based on the data in the relatively flat region of the tM surface, excluding the steep outlying data in Fig. 6(b). These respective values were found to be 0.252 and 0.013 ms. A TOF peak-to-peak estimate of the propagation time tpp was 0.244 ms, which lies within the uncertainty range estimated from the cross-correlation technique. 4.2. Detonation wave propagation The pressures from a propagating detonation wave in a stoichiometric oxyhydrogen mixture initially at 2 atm and 298 K are displayed in Fig. 7(a) [16]. Six pressure transducers were mounted 76.2 mm apart consecutively from an arc igniter at the downstream face of a detonation tube. These transducers are numbered 1–6 from the igniter. The figure shows that the wave is accelerating, indicating the possibility of a transition from deflagration to detonation (DDT). However, the existence of the transition is not easily determined by the TOF method as it gave an abrupt change in the propagation time between transducers from 0.08 to 0.03 ms, see Table 2. XCOR was applied to each pair of consecutive pressure data and the results are shown in Fig. 7(b). The flat tM-surface exists at lower values for transducer pairs further from the ignition point, indicating that the wave sped up. The wave speeds were calculated between each pair of transducers using the resultant E(tM) values as listed in Table 2. An estimate of the error in the wave speed was made based on the standard deviation of the time delays s(tM) and the inherent error er. Here s(tM) is the error resulting from the actual wave front spreading over a small time interval instead of a step rise. The inherent error er is the error associated with the inherent rise time of the sensors. In practice, er is twice the rise time.
ARTICLE IN PRESS 1104
F.K. Lu et al. / Mechanical Systems and Signal Processing 23 (2009) 1098–1111
Fig. 4. Model 5. Propagating and reflecting detonation pressure models. (a) SNR ¼ N. (b) SNR ¼ 16. (c) SNR ¼ 8.
Fig. 5. Cross-correlation surface obtained for model 5. Reflected wave values are absolute for clarity. (a) SNR ¼ N. (b) SNR ¼ 16. (c) SNR ¼ 8.
The total error of the propagating time et can be expressed as
et ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sðtM Þ2 þ e2r
(9)
ARTICLE IN PRESS F.K. Lu et al. / Mechanical Systems and Signal Processing 23 (2009) 1098–1111
1105
Table 1 Time delay results for waves in Fig. 4. SNR
Wave
tpp
E(tM)Xa
s(tM)Xb
N
Incident Reflected
20 40
20.0 140.0
0.0 0.0
16
Incident Reflected
21 44
21.3 46.5
1.2 1.2
8
Incident Reflected
22 44
21.7 46.5
2.9 1.1
a b
Expected or average value computed from XCOR. Standard deviation computed from XCOR.
Fig. 6. Example: nonstationary cross correlation of shock propagation: (a) Shock propagation past two transducers and (b) the non-stationary crosscorrection function.
Fig. 8 shows uX, the wave speed derived by XCOR as circles, with 7Du as the error bars, which can be expressed as uX Du ¼
DL EðtM Þ et
(10)
where the rise time of the pressure transducers in the experimental setup of Fig. 8 is 1 ms, thus, er ¼ 2 ms. The abscissa is the distance from the igniter. The error bar in the distance axis is simply the distance between the transducers whose locations are marked as triangles at the bottom of the figure. The speeds determined using the peak-to-peak time upp are displayed as squares. The detonation wave accelerated down the tube. It can be seen that upp is within the error bars derived by XCOR. Moreover, uX clearly shows the possible transition up to the theoretical CJ speed uCJ. This indication could not be obtained through a peak-to-peak TOF approach.
ARTICLE IN PRESS 1106
F.K. Lu et al. / Mechanical Systems and Signal Processing 23 (2009) 1098–1111
Fig. 7. Example: nonstationary cross correlation of detonation wave propagation: (a) Detonation wave propagating past six transducers and (b) the nonstationary cross-correction function.
Table 2 Propagation time by TOF and XCOR for shock, PDE and detonation data.
tpp, ms
E(tM)X, ms
s(tM)X, ms
Fig.
Resolution, ms
Wave
Fig. 6
0.002
– 1–2 2–3
0.244 0.08 0.08
0.252 0.081 0.099
0.013 0.007 0.022
Fig. 7(b)
0.01
3–4 4–5 5–6
0.03 0.03 0.03
0.048 0.031 0.028
0.021 0.004 0.005
Figs. 1 and 9
0.01
Incident Reflected
0.22 0.39
0.208 0.385
0.006 0.005
Fig. 8. Detonation wave speed past six transducers.
4.3. Detonation wave reflection XCOR was applied to the pressure data shown in Fig. 1. The values of tM for the incident and reflected waves were calculated and shown in Fig. 9. The tM-surface is well established for the incident detonation wave. From
ARTICLE IN PRESS F.K. Lu et al. / Mechanical Systems and Signal Processing 23 (2009) 1098–1111
1107
Fig. 9. tM for detonation wave.
the tM-surface, E(tM)X and s(tM)X were found to be 0.208 and 0.006 ms, respectively, as listed in Table 2. The propagation time using the TOF method was 0.22 ms [6], which is within two standard deviations of the averaged tM obtained by XCOR. 4.4. DDT run-up distance DDT is the process wherein a slow subsonic deflagration becomes a fast supersonic detonation. In this process, the required transition distance before the onset of detonation occurs, defined as the DDT run-up distance normalized by the internal diameter of the tube Xddt, is of great interest, especially in the development of pulse detonation engines (PDEs) [17]. There have been many experimental studies on DDT since the 1960s. However, there is no standard definition and apparatus to measure the onset of detonation and Xddt. Results of the DDT process for a stoichiometric propane/oxygen mixture initially at ambient conditions are shown in Fig. 10 [18]. Six pressure transducers and six photodiodes were mounted along the streamwise direction from an arc igniter at the closed end to measure the trajectories of the pressure waves and the combustion waves, respectively. The signals from the pressure transducers and photodiodes are shown as solid and dashed lines, respectively. This figure shows the creation of localized explosion, forward propagating overdriven detonation immediately after transition, and backward propagating retonation in the DDT process, which are the common features for the onset of detonation based on the photographic observation [19]. Most of the previous studies use these common features to determine Xddt. These are Method 1: TOF. Xddt is defined as the distance from the ignition source to the middle position of the two sensors, where the pressure or combustion wave propagation speed by the TOF method is initially greater than the CJ detonation speed. Method 2: Visible and near ultraviolet emission. Since the onset of detonation shows a very different emission from a deflagration, the variations in the intensity and spectrum of the flame emission are used as a criterion for the onset of detonation, shown in Fig. 10(b). Knowing the trajectory of the pressure or combustion waves allows Xddt to be determined. Method 3: Onset of retonation. The intersection of the combustion and retonation trajectories allows Xddt to be determined. Among these methods, XCOR was applied to each pair of waveforms to compute the propagation time E(tM) and its error et. The propagation time was used to determine the speeds and trajectories of the pressure wave, combustion wave, and retonation wave. These results were used to estimate Xddt. Further, et was used to estimate the uncertainty of Xddt in methods 2 and 3. For method 1, the uncertainty of Xddt is only related to the spacing between the sensors. Since et comes from statistical information, namely, Eq. (9), the uncertainties in methods 2 and 3 are referred to the overall error in the root-sum-square (RSS) formula as [20] sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qf 2 qf 2 qf 2 E arss ¼ Dz1 þ Dz2 þ þ Dzn (11) qz1 qz2 qzn where Xddt ¼ f(z1, z2, y, zn). The n independent variables zi, i ¼ 1, y, n are the estimated quantities (e.g. the propagating time E(tM) and the wave speeds uX) and have uncertainties of error by 7Dzi, (e.g. et and Du). The normalized run-up distance Xddt and its uncertainty estimated by the three methods for different equivalence ratios are displayed in Fig. 11. The
ARTICLE IN PRESS 1108
F.K. Lu et al. / Mechanical Systems and Signal Processing 23 (2009) 1098–1111
Fig. 10. Displacement–time diagram of a DDT process.
Fig. 11. Xddt by methods 1, 2 and 3.
figure shows that the uncertainties of Xddt in methods 2 and 3 are clearly defined rather than being only limited to the spacing distance of the sensors. This also indicates that the statistical information of the wave propagation time provided by the nonstationary CCF technique enables Xddt and its uncertainty to be systematically estimated.
ARTICLE IN PRESS F.K. Lu et al. / Mechanical Systems and Signal Processing 23 (2009) 1098–1111
Fig. 12. XCOR flowchart.
1109
ARTICLE IN PRESS 1110
F.K. Lu et al. / Mechanical Systems and Signal Processing 23 (2009) 1098–1111
5. Discussion and conclusions A nonstationary cross-correlation technique was developed to provide an objective estimate of the propagation time of shock and detonation waves past a pair of transducers. The nonstationary CCF depends on the time delay and on an integration window. Thus, a time delay surface is obtained. Objective estimates of the time delay and its standard deviation could be obtained by limiting the data to a certain band close to the horizontal portion of this surface. The technique was tested against a number of model functions with different amounts of noise. The cross-correlation peak was found to be accurately obtained for a SNR of more than 16. The algorithm was able to reject spurious signal pairs, thereby allowing it to be used for detecting a set of propagating and reflecting wave past a pair of transducers. The propagation times of a detonation wave past an array of six transducers were also successfully obtained. The DDT transition run-up distance and its uncertainty can be clearly defined. When compared to the more conventional TOF technique in estimating the wave speeds and DDT transition run-up distance, the present technique was objective and yielded useful statistical information. Appendix A. Numerical algorithm of XCOR A block diagram of XCOR is shown in Fig. 12. Input data are first read to determine the short-averaging time t. A trialand-error procedure was used with a value of t given by c (signal rise time), where c is constant. The finite rise time is estimated by XCOR from the measured shock wave profile. It was found that if the SNR is larger than 4, values of c from 2.5 to 3.5 did not produce much difference in determining the peak correlation time and a value of 3 was used in this study. By averaging over the time interval t for all data, the nonstationary component x(t) was extracted from the original data. Next, XCOR requires that a choice be made of the signal of interest in cases that a data record consists of more than one wave, for example, an incident and a reflected shock. This is determined by input of the initial selection of the integration window bT Lo ; T U o c as schematically shown in Fig. 1. As long as the initial window includes the nonstationarity of interest, the initial selection is arbitrary. However, TL and TU are separately expanded by the algorithm, as indicated by the doubleheaded arrows in Fig. 1, to find reasonable ranges according to three criteria, which are that Tshould (1) not include another wave, (2) exclude data that are more than three standard deviations from the mean, thereby excluding spurious signals which may not be the random noise but which may be the electrical disturbances, and (3) not include the range in which x(t) and y(t) cross each other. These criteria were coded into XCOR, so that it automatically finds the possible ranges of T L ¼ T L1 ! T L2 and T U ¼ T U 1 ! T U 2 , as shown in arrows in Fig. 1. Once bT L1 ; T L2 c and bT U 1 ; T U 2 c are determined, Eq. (8) is repeatedly applied over all pairs of TL and TU to compute tM(TL,TM). The tM surface will still include spurious results due to noise. Therefore, a criterion that the variation of tM must be the time interval of the data that was included. Finally, the average and standard deviation on tM are calculated. References [1] K. Ishii, H. Kataoka, T. Kojima, Initiation and propagation of detonation waves in combustible high speed flows, in: Proceedings of the Combustion Institute, 2009, doi:10.1016/j.proci.2008.05.029. [2] K. Kontis, R. An, D. Kounadis, H. Zare-Behtash, Head-on collision of shock wave induced vortices with a cylinder and a sphere, International Journal of Heat and Fluid Flow 29 (5) (2008) 1380–1392. [3] A.E. Rakitin, A.Y. Starikovskii, Mechanisms of deflagration-to-detonation transition under initiation by high-voltage nanosecond discharges, Combustion and Flame 155 (1–2) (2008) 343–355. [4] G. Borgioli, L. Capineri, P. Falorni, S. Matucci, C.G. Windsor, The detection of buried pipes from time-of-flight radar data, IEEE Transactions on Geoscience and Remote Sensing 46 (8) (2008) 2254–2266. [5] F. Amorini, A. Anzalone, R. Bassini, C. Boiano, G. Cardella, S. Cavallaro, E. De Filippo, P. Guazzoni, E. La Guidara, G. Lanzano`, A. Pagano, M. Papa, S. Pirrone, G. Politi, F. Porto, F. Riccio, F. Rizzo, S. Russo, P. Russotto, L. Zetta, Digital signal processing for mass identification in a 4p-detector, using time of flight measurement, IEEE Transactions on Nuclear Science 55 (2) (2008) 717–722. [6] F.K. Lu, H.-C. Liu, D.R. Wilson, Electrical conductivity channel for a shock tube, Measurement Science and Technology 16 (9) (2005) 1730–1740. [7] J.S. Bendat, A.G. Piersol, Random Data: Analysis and Measurement Procedures, second ed., Wiley-Interscience, New York, 1986. [8] C.H. Knapp, G.C. Carter, The generalized correlation method for estimation of time delay, IEEE Transactions on Acoustics, Speech, and Signal Processing 24 (4) (1976) 320–327. [9] M. Azaria, D. Hertz, Time delay estimation by generalized cross correlation methods, IEEE Transactions on Acoustics, Speech, and Signal Processing 32 (2) (1984) 280–285. [10] D. Avitzour, Time delay estimation at high signal-to-noise ratio, IEEE Transactions on Aerospace and Electronic Systems 27 (2) (1991) 234–237. [11] A. Kumar, Y. Bar-Shalom, Time-domain analysis of cross correlation for time delay estimation with an autocorrelated signal, IEEE Transactions on Signal Processing 41 (4) (1993) 1664–1668. [12] J. Benesty, J. Chen, Y. Huang, Time-delay estimation via linear interpolation and cross correlation, IEEE Transactions on Speech and Audio Processing 12 (5) (2004) 509–519. [13] S. S- is- bot, A cross-correlation technique as a system evaluation tool; application to blood flow measurement in extra-corporeal circuits, Flow Measurement and Instrumentation 16 (1) (2005) 27–34. [14] A. Hribernik, G. Bombek, Improved method for shot particle velocity measurement within a shotblasting chamber, Flow Measurement and Instrumentation 17 (2) (2006) 99–105. [15] F.K. Lu, K.M. Kinnear, Characterization of thin-film heat flux gauges, AIAA Journal of Thermophysics and Heat Transfer 13 (4) (1999) 548–550.
ARTICLE IN PRESS F.K. Lu et al. / Mechanical Systems and Signal Processing 23 (2009) 1098–1111
[16] [17] [18] [19]
1111
S.B. Stanley, W.S. Stuessy, D.R. Wilson, Experimental investigation of pulse detonation wave phenomenon, AIAA Paper 95–2197, 1995. K. Kailasanath, Recent developments in the research on pulse detonation engines, AIAA Journal 41 (2) (2003) 145–159. J. Li, W.H. Lai, K. Chung, F.K. Lu, Uncertainty analysis of deflagration-to-detonation run-up distance, Shock Waves 14 (5–6) (2005) 413–420. P.A. Urtiew, A.K. Oppenheim, Experimental observations of the transition to detonation in an explosive gas, Proceedings of the Royal Society of London Series A: Mathematical and Physical Sciences 295 (1440) (1966) 13–28. [20] E.O. Doebelin, Measurement System Application and Design, fourth ed., McGraw-Hill, New York, 1989.
Contents lists available at ScienceDirect
Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/jnlabr/ymssp
Detection of shock and detonation wave propagation by cross correlation F.K. Lu a,, A.A. Ortiz a, J.-M. Li b,1, C.H. Kim c, K.-M. Chung d a
Aerodynamics Research Center, Mechanical and Aerospace Engineering Department, University of Texas at Arlington, 214B Woolf Hall, 500 West 1st Street, Arlington, TX 76019-0018, USA Institute of Aeronautics and Astronautics, National Cheng Kung University, Tainan, Taiwan, Republic of China c Mechanical Engineering Department, Purdue University, West Lafayette, IN 47907, USA d Aerospace Science and Technology Research Center, National Cheng Kung University, Tainan, Taiwan, Republic of China b
a r t i c l e i n f o
abstract
Article history: Received 28 August 2007 Received in revised form 3 November 2008 Accepted 10 November 2008 Available online 14 November 2008
The propagation speed of a shock or detonation wave in a shock tube is usually determined by a time-of-flight method by dividing the distance between two transducers with the propagation time of the disturbance signal. Some arbitrariness is inherent in determining the propagation time by this method. An improved method for objectively determining the propagation time using a nonstationary cross-correlation technique is described. The method requires the choice of an integration window that includes the nonstationary event. The method was first tested against a number of model functions with different noise levels. It was then applied to propagating and reflected shock and detonation waves, including an example of a transitioning detonation wave propagating past six transducers. In addition, the nonstationary CCF technique was also applied to evaluate the uncertainty in estimating the deflagrationto-detonation transition run-up distance. In all cases, the time delay and its standard deviation could be obtained. & 2008 Elsevier Ltd. All rights reserved.
Keywords: Detonation Shock Shock tube Wave propagation
1. Introduction The determination of the propagation speed of a shock or detonation wave as in shock or detonation tube experiments typically makes use of the time-of-flight (TOF) method. In this method, the speed of a propagating shock, say, is obtained by determining the passage time of the wave DT between two transducers with a known separation distance DL, to yield u^ s ¼ DL=DT .
(1)
Some recent examples of the TOF method for shock and detonation wave studies include [1–3] and for other propagating signals include [4,5]. The TOF method appears to be suitable for strong, distinct transient signals such as shock and detonation waves. In situations where the signal is weak, the TOF method requires further treatment or may even have to be replaced by more powerful methods. Corresponding author. Tel.: +1 817 272 2083; fax: +1 817 272 5010.
E-mail addresses: [email protected] (F.K. Lu), [email protected] (A.A. Ortiz), [email protected] (J.-M. Li), [email protected] (C.H. Kim), [email protected] (K.-M. Chung). 1 Present address: Temasek Laboratories, National University of Singapore, Singapore. 0888-3270/$ - see front matter & 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.ymssp.2008.11.001
ARTICLE IN PRESS F.K. Lu et al. / Mechanical Systems and Signal Processing 23 (2009) 1098–1111
1099
Even though the TOF method is valid in principle and is simple to apply to propagating shock or detonation waves, its practical implementation can be subjective. This is because the shock front is not measured as a step rise, as in theory, but is spread over a small time interval. This interval depends on a number of factors, such as the size of the transducer, viscous effects and the data-sampling rate. Different arbitrary criteria can, thus, be formulated for determining the shock arrival over a transducer. For example, the shock may be deemed to have arrived at the transducer when the signal rises above an arbitrary threshold, or when the signal is midway between arbitrarily defined upstream and downstream values. These criteria, which yield signal arrival times close to one another, are subjective. They may be unreliable if the error of the time estimate is comparable to DT. Such a situation may occur if the transducers are closely spaced or if the sampling rate is low. For the same record length, increasing the sampling rate does not fundamentally overcome this difficulty although it may reduce the uncertainty in the DT estimate. To illustrate this arbitrariness, consider the propagation of a fully-developed Chapman–Jouguet (CJ) detonation wave past two transducers in a stoichiometric oxyhydrogen mixture initially at 2.01 atm, Fig. 1, the data being reported previously in [6]. (Fig. 1(a) includes other features, including a subsequent reflection that will be described later.) The region where the pressure initially rises is enlarged and shown in Fig. 1b. The symbols represent actual data, sampled at 100 kHz per channel with a simultaneous sample-and-hold system. Fig. 1b shows that the two wave fronts X(t) and Y(t) are spread over a few data points. The double-headed arrows indicate possible pairs of data points for the TOF method, resulting in estimates of DT from 0.21 to 0.23 ms. The 0.02 ms uncertainty represents a 10% error on the estimate of the velocity. A cross-correlation approach has been developed to improve the estimate of DT. This method also ensures an objective ˆ s. Moreover, the method is accurate even for noisy signals. However, since the wave method for estimating the error on u signals are not random, stationary processes [7], an extension of the usual correlation analysis, which assumes data stationarity, was needed. This nonstationary cross-correlation approach was developed and tested against a number of model functions before being applied to field signals. It can be noted that the cross-correlation technique has been widely used in various disciplines for detecting propagating signals, especially if they are weak, under the general method of timedelay estimation [8–14].
Fig. 1. Arbitrariness of the TOF method: (a) Detonation wave propagating past two transducers and (b) enlargement of wave fronts.
ARTICLE IN PRESS 1100
F.K. Lu et al. / Mechanical Systems and Signal Processing 23 (2009) 1098–1111
1.1. The cross-correlation function (CCF) For an infinitely long record of a random, stationary process, the mean value is constant, that is, it is independent of t [7]. Consider two random, stationary processes x(t) and Z(t), both with zero means. The CCF between these two processes is defined as Z 1 T xðtÞZðt þ tÞ dt. (2) RxZ ðtÞ ¼ lim T!1 T T If x(t) and Z(t) are totally uncorrelated, then RxZ(t) will be zero for all values of t. If x(t) and Z(t) are correlated, RxZ(t) will display a peak at some value of t. If the peak occurs at a positive value of t, the correlated event in Z(t) occurs after an event in x(t), and vice versa if the peak occurs at a negative value. In practice, data cannot be obtained over an infinite time. Eq. (2) is then modified for two data records of the same finite record length T to yield an estimate of the CCF Z T=2 1 xðtÞZðt þ tÞ dt. (3) R^ xZ ðtÞ ¼
T
T=2
The (4) will be dropped in the subsequent discussion. For digital implementations, the above equation is discretized to yield RxZ ðtj Þ ¼
N 1X xðtk ÞZðtk þ tj Þ. N k¼1
(4)
If there is a correlation between the two signals x(t) and Z(t), then the peak value of RxZ(t) occurs at tM, which represents the best correlation between x(t) and Z(t). If x(t) and Z(t) are from a propagating signal, tM would be the time for the signal to propagate past the two respective measurement locations. Together with the displacement between the two locations, the propagation speed can be deduced from Eq. (1), the same equation used for the TOF method. Thus, for the present application, the cross-correlation approach can be regarded as a sophisticated extension of the TOF method. An important advantage of the cross-correlation method over the TOF method is that it can provide an objective means of determining the propagation time. The resulting error in estimating the propagation speed is also quantifiable. 1.2. Nonstationary CCF The above discussion holds for random data which can be assumed to be stationary although there are many classes of data where stationarity does not hold. However, the cross-correlation technique may be extended to nonstationary data. Let the record from a transient process be expressed as XðtÞ ¼ xðtÞ þ nðtÞ,
(5)
where x(t) is a deterministic function and n(t) is a random ‘‘noise’’ with a stationary mean value of zero. Before applying a nonstationary cross-correlation, it is necessary to remove the noise as much as possible. If the variation of x(t) is slow compared to the lowest frequency in n(t), then x(t) can be separated from n(t) by low-pass filtering operations on X(t), However, low-pass filtering cannot be applied to shock wave data. The technique suggested in [7] is to use segmented mean value estimates by short-time-averaging. In this approach, the original record is divided into segments and an average value of each segment is used. The effect is a localized low-pass filtering at a frequency equal to the reciprocal of the averaging time t. The selection of t is arrived by trial-and-error. As t becomes short, the random error from n(t) increases ˆ(t) departs from the ‘‘true’’ x(t). Thus, the selection of an appropriate averaging time but as t becomes long, the estimated x involves a compromise between random and bias errors. A suitable value of t for the present study is three times the sampling time. As can be seen from Fig. 1, the shock is detected within 3–5 times of the sampling interval. A different value of t may be more appropriate for other situations. The possibility for such a modification is incorporated within the algorithm (block diagram in Appendix A). The procedure for performing the nonstationary CCF computation requires one more step than that for stationary data. Eq. (4) indicates that the CCF for a stationary process is a function of t only. However, for a nonstationary process, the CCF also depends on the time window, that is, it is also a function of t. This is expressed as a dependence on ti in general, as Rxy ðti ; tj Þ ¼
N i X 1 xðtiþk Þyðtiþk þ tj Þ. N i þ 1 k¼1
Moreover, there is a number of restrictions on ti since (1) x(ti+k) and y(ti+k+tj) must contain the nonstationarity event within k ¼ 0-Ni, (2) ti4t1, and (3) ti+kotN.
(6)
ARTICLE IN PRESS F.K. Lu et al. / Mechanical Systems and Signal Processing 23 (2009) 1098–1111
1101
Condition (1) means that the nonstationarity of interest (specifically a shock or detonation wave) must exist between ti and tN. If N datapoints can be selected from the entire range of each record, the value N constitutes another parameter as required in Eq. (6), since tN itself can take on a range of values as long as conditions (1) and (3) are satisfied. To satisfy the above conditions, therefore, Eq. (6) can be formulated as Rxy ðT L ; T U ; tj Þ ¼
U X 1 xðtk Þyðtk þ tj Þ, U L þ 1 k¼L
(7)
where L and U are the lower and upper limits (LoU) of the CCF calculation. Thus, the nonstationary event (of both records) must exist between TL and TU. Finally, since the interest here is the peak correlation time tM (corresponding to the propagation time between two points), the correlation function can be normalized against the peak value Rxy(TL,TU,tM) rxy ðT L ; T U ; tÞ ¼
1 Rxy ðT L ; T U ; tM ÞT
Z
TL
xðtÞyðt þ tÞ dt,
(8)
TL
where T ¼ T U T L . 2. Algorithm Algorithms were developed in FORTRAN, MatlabTM and LabVIEW TM to extract the propagation time of a shock or detonation wave past pairs of transducers.2 The FORTRAN program, known as XCOR, will be described. There is practically no difference in functionality between the three implementations. As seen in Eq. (8), the nonstationary CCF is a function of the time delay, and of the lower and upper limits of the sample window. The time delay t can be varied from T to þTbut the lower and upper limits of the integration window T must be carefully determined according to the characteristics of the given data. Also, the averaging time t should be determined to appropriately extract the nonstationary deterministic component from the given nonstationary random data. XCOR has the routines to perform these functions and it is described in more detail in the Appendix. It can also be noted here that, although the examples involve nonstationary pressures, the algorithm is a strictly numerical procedure and can be applied to pairs of nonstationary data in general. 3. Model calculation Before applying XCOR on experimental data, it was tested against a number of model functions with pseudo-random noise. These are
unit step functions X(t) and Y(t), both with SNR ¼ 2–32 and with Y(t) delayed by tM ¼ 50; a unit step function X(t) at t ¼ 0 with SNR ¼ 16, and a step function Y(t) at tM ¼ 50 with an amplitude of A ¼ 22–23 and SNR ¼ 16;
exponentially decaying functions X(t) ¼ ebt/100 and Y ðtÞ ¼ ebðtt
M Þ=100
, where b ¼ 22–22 and SNR ¼ 16 for both
functions, and tM ¼ 50 and
linearly decaying functions X(t) ¼ 1ct and X(t) ¼ 1c(ttM), where c ¼ 22–22 and SNR ¼ 16 for both functions, and tM ¼ 50. The SNR is the inverse of the rms noise level with respect to unity. Examples of these functions are shown in Fig. 2, while their cross correlations obtained using XCOR are shown in Fig. 3. For all models, the number of data points N is 1001 and the interval between the data points Dt is unity. The time delay tM between X(t) and Y(t) was fixed at 50 for all these cases. Thus, this exercise checks the ability of XCOR to compute the maximum of rxy(t) at this delay. Since these tests were to examine the dependence of the CCF upon data characteristics, the routine to vary the limits of the integration window [TL,TU] was not included for expediency. The data displayed in Fig. 3(a) indicate that the time delay estimated from the cross correlations for model 1 was not seriously affected by noise if the SNR was larger than 4. Thus, the SNR was fixed as 16 for the other 3 model functions. The amplitude ratio A of the second step function was varied from 0.25 to 8 in model 2 to further test the algorithm’s ability to detect the propagation of a signal subjected to attenuation or amplification. For the model profiles shown in Fig. 2(b), the computations obtained correlation peaks at t ¼ 50 as seen in Fig. 3(b). Therefore, the ability of the routine for detecting peak correlation is unaffected by variations in the amplitude. Based on the results of models 1 and 2, the SNR and amplitude were set as 16 and 1, respectively, for models 3 and 4. Model 3 simulates a detonation wave with Taylor rarefaction. The value of b for the exponential decay varied from 0.25 to 4. For all the cases in this model, the peak of the cross correlation at t ¼ tM ¼ 50 was well captured as indicated in Fig. 3(c). (For comparison, the CCF of model 1 is also plotted.) It can be seen that as the profile decays faster, the CCF becomes 2
These algorithms are available from the corresponding author or from the website http://arc.uta.edu/.
ARTICLE IN PRESS 1102
F.K. Lu et al. / Mechanical Systems and Signal Processing 23 (2009) 1098–1111
Fig. 2. Model functions. (a) Model 1: unit step function with SNR ¼ 16. (b) Model 2: variable step function with SNR ¼ 16. (c) Model 3: exponential decay; X(t) only, SNR ¼ 16. (d) Model 4: linear decay; X(t) only, SNR ¼ 16.
sharper. Finally, linear decaying profiles are shown in Fig. 2(d), with c varying from 0.25 to 4. The CCF has a peak at a time delay at t ¼ tM ¼ 50 as well, as can be seen in Fig. 3(d). From the model calculations, it can be concluded that the computed value of tM ¼ 50 is not a function of profile shape or amplitude of the transient signal. However, it might be affected by excessive noise. These tests also indicate that a sharp signal such as that from a detonation wave with Taylor rarefaction enables tM to be more distinctly resolved. Moreover, the wave system in a shock tube may involve the passage of two nearby events, such as an incident and a reflected shock. Therefore, the determination of an integration window [TL,TU] is even more important in separating out the two events for further processing. The functionality of XCOR was thus further tested using a model with two peaks to simulate the propagation and reflection of a detonation wave. Fig. 4 shows the profiles of model 5, which has a maximum pressure level of unity and an exponential decay with b ¼ 0.01 for both the incident and reflected waves. The propagation time of the incident and reflected waves were set to 20 and 40, respectively, and the SNRs were N, 16 and 8. From these profiles, tM was derived using XCOR. The absolute values of tM as a function of TL and TU are shown in Fig. 5. For the incident wave, tM is positive, but for the reflected wave, tM is negative, since Rxy(t) ¼ Ryx(tM). When SNR ¼ N, tM occurs on the tM ¼ 20 plane for the incident wave and on the tM ¼ 40 plane for the reflected wave. (Such flat surfaces are hereafter known as the ideal tM-surface.) However, when there is noise, tM deviates from the ideal tM-surface as can be seen in Fig. 5(b) and (c). Such departures were eliminated by the algorithm using pre-determined criteria in order to compute the average and standard deviation of tM. These computed values are listed in Table 1. It can be seen that noise caused a poorer estimate of the propagation time by both the peak-to-peak TOF method and the CCF method. While the TOF method apparently gave a closer estimate to the true propagation delay, the arbitrariness of the method remains and an estimate of the standard deviation may be difficult. The table shows that the propagation time for the incident wave was well captured when the SNR is larger than 8. However, for the reflected wave, the averaged value of tM was more difficult to obtain since the incident wave may not be clearly excluded from the integration window. Nevertheless, the nonstationary CCF technique can provide an objective method for obtaining the standard deviation in the wave propagation time. The common TOF method does not appear capable of achieving this. The nonstationary CCF technique is next demonstrated by applying it to 4 sets of experimental data.
ARTICLE IN PRESS F.K. Lu et al. / Mechanical Systems and Signal Processing 23 (2009) 1098–1111
1103
Fig. 3. Cross correlations of model functions. (a) Model 1: unit step function with SNR ¼ 16. (b) Model 2: variable step function with SNR ¼ 16. (c) Model 3: exponential decay; X(t) only, SNR ¼ 16. (d) Model 4: linear decay; X(t) only, SNR ¼ 16.
4. Application of nonstationary CCF 4.1. Shock wave propagation Fig. 6a shows voltages from two pressure transducers, where the step rise in the signal indicates the passage of a shock wave [15]. XCOR was applied to the wave and the resulting tM(TL,TU) is shown in Fig. 6(b). It can be seen that the tM surface is relatively flat except for the steep rise when the upper limits of the integration window approach the time of shock arrival. The average and the standard deviation of tM were computed based on the data in the relatively flat region of the tM surface, excluding the steep outlying data in Fig. 6(b). These respective values were found to be 0.252 and 0.013 ms. A TOF peak-to-peak estimate of the propagation time tpp was 0.244 ms, which lies within the uncertainty range estimated from the cross-correlation technique. 4.2. Detonation wave propagation The pressures from a propagating detonation wave in a stoichiometric oxyhydrogen mixture initially at 2 atm and 298 K are displayed in Fig. 7(a) [16]. Six pressure transducers were mounted 76.2 mm apart consecutively from an arc igniter at the downstream face of a detonation tube. These transducers are numbered 1–6 from the igniter. The figure shows that the wave is accelerating, indicating the possibility of a transition from deflagration to detonation (DDT). However, the existence of the transition is not easily determined by the TOF method as it gave an abrupt change in the propagation time between transducers from 0.08 to 0.03 ms, see Table 2. XCOR was applied to each pair of consecutive pressure data and the results are shown in Fig. 7(b). The flat tM-surface exists at lower values for transducer pairs further from the ignition point, indicating that the wave sped up. The wave speeds were calculated between each pair of transducers using the resultant E(tM) values as listed in Table 2. An estimate of the error in the wave speed was made based on the standard deviation of the time delays s(tM) and the inherent error er. Here s(tM) is the error resulting from the actual wave front spreading over a small time interval instead of a step rise. The inherent error er is the error associated with the inherent rise time of the sensors. In practice, er is twice the rise time.
ARTICLE IN PRESS 1104
F.K. Lu et al. / Mechanical Systems and Signal Processing 23 (2009) 1098–1111
Fig. 4. Model 5. Propagating and reflecting detonation pressure models. (a) SNR ¼ N. (b) SNR ¼ 16. (c) SNR ¼ 8.
Fig. 5. Cross-correlation surface obtained for model 5. Reflected wave values are absolute for clarity. (a) SNR ¼ N. (b) SNR ¼ 16. (c) SNR ¼ 8.
The total error of the propagating time et can be expressed as
et ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sðtM Þ2 þ e2r
(9)
ARTICLE IN PRESS F.K. Lu et al. / Mechanical Systems and Signal Processing 23 (2009) 1098–1111
1105
Table 1 Time delay results for waves in Fig. 4. SNR
Wave
tpp
E(tM)Xa
s(tM)Xb
N
Incident Reflected
20 40
20.0 140.0
0.0 0.0
16
Incident Reflected
21 44
21.3 46.5
1.2 1.2
8
Incident Reflected
22 44
21.7 46.5
2.9 1.1
a b
Expected or average value computed from XCOR. Standard deviation computed from XCOR.
Fig. 6. Example: nonstationary cross correlation of shock propagation: (a) Shock propagation past two transducers and (b) the non-stationary crosscorrection function.
Fig. 8 shows uX, the wave speed derived by XCOR as circles, with 7Du as the error bars, which can be expressed as uX Du ¼
DL EðtM Þ et
(10)
where the rise time of the pressure transducers in the experimental setup of Fig. 8 is 1 ms, thus, er ¼ 2 ms. The abscissa is the distance from the igniter. The error bar in the distance axis is simply the distance between the transducers whose locations are marked as triangles at the bottom of the figure. The speeds determined using the peak-to-peak time upp are displayed as squares. The detonation wave accelerated down the tube. It can be seen that upp is within the error bars derived by XCOR. Moreover, uX clearly shows the possible transition up to the theoretical CJ speed uCJ. This indication could not be obtained through a peak-to-peak TOF approach.
ARTICLE IN PRESS 1106
F.K. Lu et al. / Mechanical Systems and Signal Processing 23 (2009) 1098–1111
Fig. 7. Example: nonstationary cross correlation of detonation wave propagation: (a) Detonation wave propagating past six transducers and (b) the nonstationary cross-correction function.
Table 2 Propagation time by TOF and XCOR for shock, PDE and detonation data.
tpp, ms
E(tM)X, ms
s(tM)X, ms
Fig.
Resolution, ms
Wave
Fig. 6
0.002
– 1–2 2–3
0.244 0.08 0.08
0.252 0.081 0.099
0.013 0.007 0.022
Fig. 7(b)
0.01
3–4 4–5 5–6
0.03 0.03 0.03
0.048 0.031 0.028
0.021 0.004 0.005
Figs. 1 and 9
0.01
Incident Reflected
0.22 0.39
0.208 0.385
0.006 0.005
Fig. 8. Detonation wave speed past six transducers.
4.3. Detonation wave reflection XCOR was applied to the pressure data shown in Fig. 1. The values of tM for the incident and reflected waves were calculated and shown in Fig. 9. The tM-surface is well established for the incident detonation wave. From
ARTICLE IN PRESS F.K. Lu et al. / Mechanical Systems and Signal Processing 23 (2009) 1098–1111
1107
Fig. 9. tM for detonation wave.
the tM-surface, E(tM)X and s(tM)X were found to be 0.208 and 0.006 ms, respectively, as listed in Table 2. The propagation time using the TOF method was 0.22 ms [6], which is within two standard deviations of the averaged tM obtained by XCOR. 4.4. DDT run-up distance DDT is the process wherein a slow subsonic deflagration becomes a fast supersonic detonation. In this process, the required transition distance before the onset of detonation occurs, defined as the DDT run-up distance normalized by the internal diameter of the tube Xddt, is of great interest, especially in the development of pulse detonation engines (PDEs) [17]. There have been many experimental studies on DDT since the 1960s. However, there is no standard definition and apparatus to measure the onset of detonation and Xddt. Results of the DDT process for a stoichiometric propane/oxygen mixture initially at ambient conditions are shown in Fig. 10 [18]. Six pressure transducers and six photodiodes were mounted along the streamwise direction from an arc igniter at the closed end to measure the trajectories of the pressure waves and the combustion waves, respectively. The signals from the pressure transducers and photodiodes are shown as solid and dashed lines, respectively. This figure shows the creation of localized explosion, forward propagating overdriven detonation immediately after transition, and backward propagating retonation in the DDT process, which are the common features for the onset of detonation based on the photographic observation [19]. Most of the previous studies use these common features to determine Xddt. These are Method 1: TOF. Xddt is defined as the distance from the ignition source to the middle position of the two sensors, where the pressure or combustion wave propagation speed by the TOF method is initially greater than the CJ detonation speed. Method 2: Visible and near ultraviolet emission. Since the onset of detonation shows a very different emission from a deflagration, the variations in the intensity and spectrum of the flame emission are used as a criterion for the onset of detonation, shown in Fig. 10(b). Knowing the trajectory of the pressure or combustion waves allows Xddt to be determined. Method 3: Onset of retonation. The intersection of the combustion and retonation trajectories allows Xddt to be determined. Among these methods, XCOR was applied to each pair of waveforms to compute the propagation time E(tM) and its error et. The propagation time was used to determine the speeds and trajectories of the pressure wave, combustion wave, and retonation wave. These results were used to estimate Xddt. Further, et was used to estimate the uncertainty of Xddt in methods 2 and 3. For method 1, the uncertainty of Xddt is only related to the spacing between the sensors. Since et comes from statistical information, namely, Eq. (9), the uncertainties in methods 2 and 3 are referred to the overall error in the root-sum-square (RSS) formula as [20] sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qf 2 qf 2 qf 2 E arss ¼ Dz1 þ Dz2 þ þ Dzn (11) qz1 qz2 qzn where Xddt ¼ f(z1, z2, y, zn). The n independent variables zi, i ¼ 1, y, n are the estimated quantities (e.g. the propagating time E(tM) and the wave speeds uX) and have uncertainties of error by 7Dzi, (e.g. et and Du). The normalized run-up distance Xddt and its uncertainty estimated by the three methods for different equivalence ratios are displayed in Fig. 11. The
ARTICLE IN PRESS 1108
F.K. Lu et al. / Mechanical Systems and Signal Processing 23 (2009) 1098–1111
Fig. 10. Displacement–time diagram of a DDT process.
Fig. 11. Xddt by methods 1, 2 and 3.
figure shows that the uncertainties of Xddt in methods 2 and 3 are clearly defined rather than being only limited to the spacing distance of the sensors. This also indicates that the statistical information of the wave propagation time provided by the nonstationary CCF technique enables Xddt and its uncertainty to be systematically estimated.
ARTICLE IN PRESS F.K. Lu et al. / Mechanical Systems and Signal Processing 23 (2009) 1098–1111
Fig. 12. XCOR flowchart.
1109
ARTICLE IN PRESS 1110
F.K. Lu et al. / Mechanical Systems and Signal Processing 23 (2009) 1098–1111
5. Discussion and conclusions A nonstationary cross-correlation technique was developed to provide an objective estimate of the propagation time of shock and detonation waves past a pair of transducers. The nonstationary CCF depends on the time delay and on an integration window. Thus, a time delay surface is obtained. Objective estimates of the time delay and its standard deviation could be obtained by limiting the data to a certain band close to the horizontal portion of this surface. The technique was tested against a number of model functions with different amounts of noise. The cross-correlation peak was found to be accurately obtained for a SNR of more than 16. The algorithm was able to reject spurious signal pairs, thereby allowing it to be used for detecting a set of propagating and reflecting wave past a pair of transducers. The propagation times of a detonation wave past an array of six transducers were also successfully obtained. The DDT transition run-up distance and its uncertainty can be clearly defined. When compared to the more conventional TOF technique in estimating the wave speeds and DDT transition run-up distance, the present technique was objective and yielded useful statistical information. Appendix A. Numerical algorithm of XCOR A block diagram of XCOR is shown in Fig. 12. Input data are first read to determine the short-averaging time t. A trialand-error procedure was used with a value of t given by c (signal rise time), where c is constant. The finite rise time is estimated by XCOR from the measured shock wave profile. It was found that if the SNR is larger than 4, values of c from 2.5 to 3.5 did not produce much difference in determining the peak correlation time and a value of 3 was used in this study. By averaging over the time interval t for all data, the nonstationary component x(t) was extracted from the original data. Next, XCOR requires that a choice be made of the signal of interest in cases that a data record consists of more than one wave, for example, an incident and a reflected shock. This is determined by input of the initial selection of the integration window bT Lo ; T U o c as schematically shown in Fig. 1. As long as the initial window includes the nonstationarity of interest, the initial selection is arbitrary. However, TL and TU are separately expanded by the algorithm, as indicated by the doubleheaded arrows in Fig. 1, to find reasonable ranges according to three criteria, which are that Tshould (1) not include another wave, (2) exclude data that are more than three standard deviations from the mean, thereby excluding spurious signals which may not be the random noise but which may be the electrical disturbances, and (3) not include the range in which x(t) and y(t) cross each other. These criteria were coded into XCOR, so that it automatically finds the possible ranges of T L ¼ T L1 ! T L2 and T U ¼ T U 1 ! T U 2 , as shown in arrows in Fig. 1. Once bT L1 ; T L2 c and bT U 1 ; T U 2 c are determined, Eq. (8) is repeatedly applied over all pairs of TL and TU to compute tM(TL,TM). The tM surface will still include spurious results due to noise. Therefore, a criterion that the variation of tM must be the time interval of the data that was included. Finally, the average and standard deviation on tM are calculated. References [1] K. Ishii, H. Kataoka, T. Kojima, Initiation and propagation of detonation waves in combustible high speed flows, in: Proceedings of the Combustion Institute, 2009, doi:10.1016/j.proci.2008.05.029. [2] K. Kontis, R. An, D. Kounadis, H. Zare-Behtash, Head-on collision of shock wave induced vortices with a cylinder and a sphere, International Journal of Heat and Fluid Flow 29 (5) (2008) 1380–1392. [3] A.E. Rakitin, A.Y. Starikovskii, Mechanisms of deflagration-to-detonation transition under initiation by high-voltage nanosecond discharges, Combustion and Flame 155 (1–2) (2008) 343–355. [4] G. Borgioli, L. Capineri, P. Falorni, S. Matucci, C.G. Windsor, The detection of buried pipes from time-of-flight radar data, IEEE Transactions on Geoscience and Remote Sensing 46 (8) (2008) 2254–2266. [5] F. Amorini, A. Anzalone, R. Bassini, C. Boiano, G. Cardella, S. Cavallaro, E. De Filippo, P. Guazzoni, E. La Guidara, G. Lanzano`, A. Pagano, M. Papa, S. Pirrone, G. Politi, F. Porto, F. Riccio, F. Rizzo, S. Russo, P. Russotto, L. Zetta, Digital signal processing for mass identification in a 4p-detector, using time of flight measurement, IEEE Transactions on Nuclear Science 55 (2) (2008) 717–722. [6] F.K. Lu, H.-C. Liu, D.R. Wilson, Electrical conductivity channel for a shock tube, Measurement Science and Technology 16 (9) (2005) 1730–1740. [7] J.S. Bendat, A.G. Piersol, Random Data: Analysis and Measurement Procedures, second ed., Wiley-Interscience, New York, 1986. [8] C.H. Knapp, G.C. Carter, The generalized correlation method for estimation of time delay, IEEE Transactions on Acoustics, Speech, and Signal Processing 24 (4) (1976) 320–327. [9] M. Azaria, D. Hertz, Time delay estimation by generalized cross correlation methods, IEEE Transactions on Acoustics, Speech, and Signal Processing 32 (2) (1984) 280–285. [10] D. Avitzour, Time delay estimation at high signal-to-noise ratio, IEEE Transactions on Aerospace and Electronic Systems 27 (2) (1991) 234–237. [11] A. Kumar, Y. Bar-Shalom, Time-domain analysis of cross correlation for time delay estimation with an autocorrelated signal, IEEE Transactions on Signal Processing 41 (4) (1993) 1664–1668. [12] J. Benesty, J. Chen, Y. Huang, Time-delay estimation via linear interpolation and cross correlation, IEEE Transactions on Speech and Audio Processing 12 (5) (2004) 509–519. [13] S. S- is- bot, A cross-correlation technique as a system evaluation tool; application to blood flow measurement in extra-corporeal circuits, Flow Measurement and Instrumentation 16 (1) (2005) 27–34. [14] A. Hribernik, G. Bombek, Improved method for shot particle velocity measurement within a shotblasting chamber, Flow Measurement and Instrumentation 17 (2) (2006) 99–105. [15] F.K. Lu, K.M. Kinnear, Characterization of thin-film heat flux gauges, AIAA Journal of Thermophysics and Heat Transfer 13 (4) (1999) 548–550.
ARTICLE IN PRESS F.K. Lu et al. / Mechanical Systems and Signal Processing 23 (2009) 1098–1111
[16] [17] [18] [19]
1111
S.B. Stanley, W.S. Stuessy, D.R. Wilson, Experimental investigation of pulse detonation wave phenomenon, AIAA Paper 95–2197, 1995. K. Kailasanath, Recent developments in the research on pulse detonation engines, AIAA Journal 41 (2) (2003) 145–159. J. Li, W.H. Lai, K. Chung, F.K. Lu, Uncertainty analysis of deflagration-to-detonation run-up distance, Shock Waves 14 (5–6) (2005) 413–420. P.A. Urtiew, A.K. Oppenheim, Experimental observations of the transition to detonation in an explosive gas, Proceedings of the Royal Society of London Series A: Mathematical and Physical Sciences 295 (1440) (1966) 13–28. [20] E.O. Doebelin, Measurement System Application and Design, fourth ed., McGraw-Hill, New York, 1989.