Analysis of shock wave unsteadiness using space and time correlations applied to shadowgraph flow visualization data

Abstract

Unsteady flow characteristics of a normal shock wave, a lambda foot, and a separated turbulent boundary layer are investigated within a unique test section with supersonic inlet flow. The supersonic wind tunnel facility, containing this test section, provides a Mach number of approximately 1.54 at the test section entrance. Digitized shadowgraph flow visualization data are employed to visualize shock wave structure within the test section. These data are analyzed to determine shock wave unsteadiness characteristics, including grayscale spectral energy variations with frequency, as well as time and space correlations, which give coherence and time lag properties associated with perturbations associated with different flow regions. Results illustrate the complexity and unsteadiness of shock-wave-boundary-layer-interactions, including event frequencies from grayscale spectral energy distributions determined using a Lagrangian approach applied to shock wave location, and by grayscale spectral energy distributions determined using ensemble-averaging applied to multiple closely-located stationary pixel locations. Auto-correlation function results and two-point correlation functions (in the form of magnitude squared coherence) quantify the time-scales of periodic events, as well as the coherence of flow perturbations associated with different locations, over a range of frequencies. Associated time lag data provide information on the originating location of perturbation events, as well as the propagation direction and event sequence associated with different flow locations. Additional insight into spatial variations of time lag and flow coherence is provided by application of magnitude squared coherence analysis to multiple locations, relative to a single location associated with the normal shock wave.

1 Introduction

Shock wave boundary layer interactions occur in many aerospace applications. Some examples include ramjet isolator ducts, turbine blade tip gaps, and transonic wings. However, low frequency unsteadiness associated with shock wave boundary layer interactions is not well understood and continues to be a debated topic. This is because different studies concluded different and contradicting results. Some results indicated that shock wave unsteadiness originates from the turbulence in the upstream boundary layer. Other results indicated that a separation region formed as part of the shock wave boundary layer interaction generates the unsteadiness. Still other studies showed that both the upstream and downstream sources are responsible for shock wave unsteadiness. Clemens and Narayanaswamy [1] summarized some of the possibilities. Many methods and analysis approaches were used to collect and analyze associated data. For example, power spectral densities and cross correlations were often used, especially for wall pressure signal data. Many of these investigations also employed Schlieren and shadowgraph imaging to determine unsteady and spatially-varying shock wave structure. Erengil and Dolling [2], Handa et al. [3], Dupont et al. [4], and Bruce and Babinsky [5], among others, measured static pressure variations with time. Spatial pressure variations were then associated with shock wave position and motion. Power spectral densities were computed from time sequence pressure data at specific locations to determine the frequencies of the unsteadiness at various locations in the flow field. Bruce and Babinsky [5] also used schlieren imaging to determine the streamwise location of the shock wave. Gamba [6] used shock wave position with respect to time, in addition to the single location pressure measurements, to determine unsteadiness characteristics. Cross-spectrum and time lag were computed for pressure fluctuations and shock wave positions, such that pressure fluctuations were given relative to the position of the upstream most shock wave. Ganapathisubramani et al. [7] and Humble et al. [8] computed cross correlations between turbulent structures in the upstream boundary layer and the shock wave from particle image velocimetry data. They detected significant coherence between the turbulent structures in the upstream boundary layer and the shock wave motion. Other recent investigations which considered different methods of shock wave data analysis were described by Ganapathisubramani and Clemens [9], Wu and Martin [10], Piponniau et al. [11], Pirozzoli et al. [12], Touber and Sandham [13], and Grilli et al. [14]. Settles [15] also provided useful information in regard to schlieren and shadowgraph visualization techniques.

The present investigation differs from many previous studies by providing greater detail of the analysis techniques used, and by using digitized shadowgraph data to visualize shock wave structure. Considered are unsteady flow and structural characteristics of a normal shock wave, a lambda foot, and a separated turbulent boundary layer within a unique research test section. The supersonic wind tunnel facility, containing this test section, provides a Mach number of approximately 1.54 at the test section entrance. The digitized shadowgraph flow visualization data are analyzed to provide grayscale spectral energy variations with frequency, as well as time and space correlations, which give coherence and time lag properties associated with perturbations associated with different flow regions. As such, the current investigation is aimed at improved analysis techniques for shock-wave-boundary-layer-interactions, for better physical understanding of these flows, and for development of improved techniques for control and management of flow unsteadiness associated with particular shock wave arrangements.

2 Experimental apparatus and procedures

2.1 Supersonic wind tunnel

The supersonic wind tunnel (which is also referred to as the SS/TS/WT or SuperSonic/TranSonic/WindTunnel) is a blow down type facility, which is located within the Johnson Research Center on the University of Alabama in Huntsville campus. The low-pressure piping system used in these experiments (which is a portion of the overall piping system) consists of an air compressor, a vertical air supply tank, a series of pressure relief valves, a manual gate valve, a pneumatic flow control valve, a pressure regulating gate valve, and an air diverter plenum with 12 cubic meter volume. This plenum is then followed by the test section segment, an exhaust plenum with a volume of 2 cubic meters, and an exhaust piping system.

The volume of the vertical air supply tank is 14 cubic meters. The tank pressure rating is 2.1 MPa. As the facility is prepared for testing, one of two compressors (either a Quincy QR 350, model BM350HPDT compressor with a ZEKS 75NDQCA100 air dryer, or a Bauer BP26-E3 compressor) fills the vertical supply tank with dried air. After the tank is pressurized to the appropriate supply pressure, the compressor is shut down, and the small Ham-Let H700SSL3/4TLD manual ball valve between the compressor and the supply tank is closed.

As a blow down test is conducted, air flows through a 152.6 mm diameter pipe, which is located downstream of the vertical supply tank. A NEWCO N36726 manual gate valve, a Truline 330AITFM-CH-SQ pneumatic flow control ball valve, and a Fisher Valves Type 667-EWT-DVC6200 pressure regulating valve are located along this pipe. The manual gate valve is open fully during testing. The pneumatic flow control ball valve is also fully open during testing. The pressure regulating valve is controlled by a Fisher FIELDVUE DVC6200 Digital Valve Controller. The digital valve controller is set to maintain a constant downstream pressure of 413.7 kPa as an individual blow down test is underway. Downstream of the pressure regulating valve, the air flows into the air diverter plenum. Six Kunkle 913BFEM03 pressure relief valves provide means to vent over-pressurized air, both before and after the manual gate valve. Flow static pressure is measured before and after each control valve within the piping system, using both Ashcroft analog pressure gauges and Honeywell FPA 060-C860–15 digital pressure transducers.

Three independent wind tunnel segments are connected to the downstream parts of the air diverter plenum. Three additional Kunkle 913BFEM03 pressure relief valves and a burst disk are additionally mounted on the air diverter plenum for pressure relief in the unlikely event of over-pressurization. A spectacle blind valve or a circular sliding gate valve is located at the inlet to each of the wind tunnel segments. To conduct a test with one such segment, the associated inlet valve is open, whereas valves at the inlets of the other segments are closed. The research test section, which is employed for the present study, is used as one of the wind tunnel segments. This segment includes a spectacle blind valve, which is followed by a straight inlet duct, a converging/diverging supersonic nozzle, the test section, an exhaust diffuser and plenum, and an exhaust piping system with three independent noise baffles. Each noise baffle bolts to the associated 203 mm exhaust duct by means of a flange, and is approximately 1 m in length, with an inner diameter of 203 mm.

The converging diverging nozzle accelerates the flow to supersonic velocities. The flow then enters the test section. The test section is constructed with a flat bottom wall, a diverging top wall, and two side walls. The test section contains a shock wave holding plate and a choking flap. The angle of the choking flap, as well as the normal and streamwise positions of the shock wave holding plate, are all adjustable.

The exhaust plenum includes four air exhaust duct vents. Two 203 mm diameter exhaust ducts vent the air outside of the laboratory through 7.6 m of straight 203 mm diameter pipe. A 38 mm diameter vent vents to the air inside of the laboratory. These three vents are always open. A Kunkle 913BFEM03 pressure relief valve and a burst disk are additionally mounted on the plenum for pressure relief in the unlikely event of over-pressurization. A 102 mm diameter vent is also installed on the exhaust plenum. This particular vent includes a custom made spectacle blind valve which is adjustable. For the present tests, the 102 mm diameter exhaust passage is half-open.

2.2 Shadowgraph optical flow visualization system

During wind tunnel tests, an Edmond Optics Corporation shadowgraph system records time-varying, shock wave flow features within the test section. The system consists of a white light source and two 152.4 mm diameter focusing mirrors. The light source produces a beam of light that is reflected off the first mirror. The beam is then columnated as it passes through the test section. Because each viewed image is generally invariant in the spanwise direction, the mirrors are aligned such that the light beam is parallel to the spanwise direction and orthogonal to the side walls. This alignment provides the images with the greatest contrast and visualization of density-varying flow features. The second mirror additionally focuses the resulting shadowgraph image into the NIKKOR 200 mm lens, which is connected to a Phantom v711 camera. This high speed camera captures a time sequence of digitized flow visualization images during each wind tunnel test. Phantom Camera Control Application 2.7 software then processes the images which are captured by the camera. The Phantom v711 Camera has a 20 μm resolution at all sampling rates, and can acquire a sequence of images at rates as high as 1400 kHz. Figure 1 shows an example instantaneous shadowgraph image, obtained with this apparatus, as part of the present investigation.

Fig. 1

Shadowgraph flow visualization image for data acquired on 05 April 2018, illustrating the presence of a normal shock wave, lambda foot, and separated turbulent boundary layer within the lower flow passage, and an oblique shock wave system within the upper flow passage

2.3 Pressure measurements

Surface static pressures are measured within the test section, the exhaust plenum, and downstream of the pressure regulating valve. The static pressure tap used for this purpose within the test section is located along the centerline of the bottom wall, upstream of the shock wave. Each pressure is measured using a Honeywell FPA 060-C54985172080 digital pressure transducer. Each transducer is connected to a National Instruments Voltage Input Module, NI 9209, mounted within a NI CompactDAQ USB Chassis, NI cDAQ-9174. These devices together convert the analog voltage signal from the transducers into a digital signal which is readable using LabVIEW 2015 Full Development System version 15.0f2 software.

2.4 Experimental testing procedure

Wind tunnel tests typically last approximately eight seconds from the time the pressure regulating valve is opened to the time it begins to close. For the majority of tests, data are captured from a time before the pressure regulating valve is opened, until after the flow stops. Flow conditions are typically established and maintained for about five seconds during each test. Only this portion of the wind tunnel test is considered during the subsequent analysis. The shadowgraph images captured during start-up and shut-down, prior and after flow is established, are generally discarded.

2.5 Data analysis procedures

The digitized time sequence of flow visualization images are processed to compute frequency spectra, and different correlations using National Instruments’ LabVIEW 2015 Full Development System version 15.0f2 software, and Mathworks’ MATLAB versions R2014b through R2018a software.

2.5.1 Analysis of flow visualization time sequence data and determination of frequency spectra

Frequency spectra are computed in six steps: extract data for a particular pixel location, filter the data, transform and normalize the data into the frequency domain, smooth the frequency domain results, remove the white noise, and ensemble average the data associated with multiple pixel locations.

2.5.2 Extract pixel data

For a selected collection of images, a MATLAB program extracts the time sequence of grayscale values at any specified pixel. During a wind tunnel test, flow conditions are established and maintained for a time interval. Only this portion of the wind tunnel test is considered during the analysis. The MATLAB code used to extract grayscale value data from specific pixels in the images has three necessary inputs. The user indicates which pixel to analyze, the file path to the folder of bitmap images during established flow conditions, and the output file name. MATLAB’s “imread” function determines the grayscale value of the specified pixel for each image in the folder. The output is a time sequenced array of grayscale values. MATLAB’s “xlswrite” function saves the associated array in a Microsoft Excel spreadsheet.

2.5.3 Filter the time sequence

The time sequence of grayscale values must be read into LabVIEW software in order to be filtered. The file path to the saved text file is specified within a subroutine in LabVIEW called the “Read from Measurements File” Express VI. The subroutine imports the time sequence of grayscale values into LabVIEW from the tab delimited text file. Once in LabVIEW, the data is converted from the dynamic data datatype to an array of double precision floating point values.

A low-pass Butterworth filter is used to filter the time sequence data. A Butterworth filter is used, as opposed to another type of filter, because of its relatively uniform transfer function scaling over the frequencies of interest. The transfer function is the ratio of the filtered value to the unfiltered value as a function of frequency. Uniform scaling over the transfer function is important because it does not add extraneous frequency content to the data. This ensures that any peaks in the frequency domain data are from real events, not distortion from the filter.

Figure 2 shows the magnitude of the transfer function with respect to frequency ratio for the first through fifth order Butterworth filters. The first through fifth order filters are labeled as A through E, respectively, in the figure. The frequency ratio on the abscissa is the frequency divided by the cut-off frequency. The figure shows how the transfer function is generally uniform but drops rapidly as the frequency approaches the cut-off frequency. The cut-off frequency is defined as the frequency where the filtered frequency content is 3 dB lower than the unfiltered frequency content. The drop-off of the transfer function’s magnitude near the cut-off frequency is more abrupt for higher order filters than lower order filters. The fifth order filter is used in this analysis because of the abrupt decrease in magnitude of the transfer function near the cut-off frequency.

Fig. 2

The transform function of a Low Pass Butterworth Filter, where lines A through E represent orders 1 through 5 [17]

The time sequence array of grayscale values is filtered with a fifth-order, low-pass Butterworth filter, using a cut-off frequency of 1 Hz less than the Nyquist folding frequency. The filtering is done with LabVIEW’s “Butterworth Filter” VI subroutine. The data, the cut-off frequency, the sampling rate, the type of filter (low pass, for example), and the order of the filter are inputs to this subroutine. This output is the filtered data as an array of double precision floating point values. After being filtered, the data are converted to a dynamic data datatype. The resulting filtered data are saved to a tab delimited spreadsheet.

2.5.4 Transformation to the frequency domain

The transformation from the time domain to the frequency domain is accomplished using MATLAB software. This transformation gives the square root of the power spectra of grayscale values as it varies with frequency, also referred to as the grayscale spectral energy distribution. The filtered data is read into MATLAB from the spreadsheet using the “xlsread” function. This function’s inputs are the file path and the range of cells in the spreadsheet that contain the data. The function’s output is an array of the filtered data. The array of the filtered time sequence data undergoes a single sided discrete fast Fourier transform. This transform is done using the Fast Fourier Transform —“fft”— function in MATLAB software. The filtered time sequenced array of grayscale values is input into the Fast Fourier Transform function. The function gives the double sided frequency transform of that time sequenced array. Next, the second half of the frequency transform array is discarded. The resulting array is the single sided frequency transform.

This transform is normalized such that it becomes the square root of the power spectra of the time domain function. With the exception of the first term in the frequency spectra array, this relationship is described by the equation given by

$$ Y(f)=\frac{2}{N}\cdotp \left|\mathcal{F}\left(y(t)\right)\right| $$

(1)

Equation (1) indicates that the magnitude of each term in the single sided frequency transform is multiplied by 2 and divided by the number of indices in the time sequence array. The first term in this transform represents the time averaged grayscale value and corresponds to a frequency of zero. The value of this term is set equal to half of the value attained by using Eq. (1). Each term in the frequency transform corresponds to a unique frequency. The frequencies steadily increase from zero to the Nyquist folding frequency. The frequency resolution is the difference in frequency between any two consecutive terms in the frequency transform. As such, the frequency resolution is determined using

$$ \varDelta f=\frac{Fs}{N} $$

(2)

The frequency resolution is defined as the sampling frequency divided by the number of points in the time sequence array. The frequency resolution is therefore a function of the sampling frequency and the amount of time that useful data are collected.

The transformation preserves spectral energy content such that the variance of the time domain function is equal to half of the sum of the squares of the frequency domain function, excluding the steady state term. This relationship is described using the equation given by

$$ \frac{1}{2}\cdotp \overline{Y{(f)}^2}=\overline{y(t){\prime}^2} $$

(3)

The variance formula used on the time domain data to attain the energy is expressed using

$$ \overline{y(t){\prime}^2}=\frac{\sum_{k=1}^N{\left(y{(t)}_k-{\sum}_{j=1}^N\frac{y{(t)}_j}{N}\right)}^2}{N}=\overline{{\left(y{(t)}_k-\overline{y{(t)}_k}\right)}^2} $$

(4)

The first step in determining the variance is to calculate the average of the grayscale values in the time sequence. Next, the difference between each of the grayscale values in the time sequence and that average value is computed. Each of these differences is then squared. Finally, all of those squared differences are averaged together. That average gives a single value that is the variance of the time sequence. It is also equal to the total spectral energy.

The overall energy from the frequency domain is evaluated using an equation of the form

$$ \frac{1}{2}\cdotp \overline{Y{(f)}^2}=\frac{1}{2}\ast {\sum}_{k=2}^{\frac{N}{2}-1}Y{(f)}_k^2 $$

(5)

Equation (5) states that the total spectral energy is equal to half of the sum of the squared grayscale values in the single sided frequency spectra. The summation does not include the value of the first term where the frequency is zero. Associated results may be expressed with respect to frequency or with respect to the Strouhal Number. The Strouhal Number is determined using the equation given by

$$ {St}_r=\frac{2\pi \cdotp f\cdotp {\delta}_o}{u_{\infty }} $$

(6)

The Strouhal Number is the frequency scaled by the boundary layer thickness and the incoming freestream velocity at the entrance to the test section. Similar to frequency, the Strouhal number has a resolution. The resolution of the Strouhal Number is determined using an equation given by

$$ \varDelta {St}_r=\frac{2\pi \cdotp \varDelta f\cdotp {\delta}_o}{u_{\infty }} $$

(7)

The resolution of the Strouhal number is the dimensionless form of the frequency resolution.

2.5.5 Smoothing the frequency transform

The transformed arrays have large variations of grayscale spectral energy distribution values, given as Y(f), over small ranges of frequency. As mentioned, $1/2\ast \overline{Y{(f)}^2}$ is the spectral energy content. The variation of the grayscale spectral energy distribution varies with frequency. Figure 3 shows a typical example of the frequency transform result with respect to frequency, obtained using Eq. (1). Note that data smoothing is required to determine detailed frequency content.

Fig. 3

Grayscale spectral energy distribution before smoothing, for the data associated with pixel [736, 329] on the shock wave of the test on 13 November 2017. The spectrum is representative of unsmoothed grayscale spectral energy distribution data

Smoothing is imposed using a simple running average which is applied to the normalized transform array data, such as are shown in Fig. 3. For a simple running average, each data point in an array is replaced by the average of the data within a specified interval centered on that data point. Table 1 summarizes the smoothing process used in the analysis. With this approach, the first 20 points in each array, corresponding to the points at the 20 lowest frequencies, are not averaged; they retain initial values. For data between the 21st point and the data point corresponding to 20 Hz, each data point is replaced by the value which is equal to the average of 3 data points comprised of the considered data point and the adjacent point on each side of the data point considered. A similar scheme is done for data corresponding to frequencies greater than 20 Hz, but the number of points per side employed to determine the smoothed value is given by

$$ 0.05\cdotp \frac{f}{\varDelta f} $$

(8)

Table 1 Spectral local data averaging procedures over different frequency ranges

This expression is equal to 5% of the array index of the data point considered. The array index is the frequency at the considered point divided by the frequency resolution. This expression is a function of the data point’s frequency. Note that the expression 0.05 · f/Δf is inapplicable for very high frequency data because there are not 0.05 · f/Δf points at higher frequencies remaining in the array. The expression becomes invalid for frequencies that are greater than approximately 95% of the Nyquist folding frequency which equals Fs/2. This limit is given by the expression

$$ f\le 0.95\cdotp \frac{Fs}{2} $$

(9)

As a result, averaging is not applied to the data points at frequencies greater than this limit. For the data corresponding to frequencies greater than approximately 95% of the Nyquist folding frequency, each value is replaced by a single value that is determined by averaging of all the non-smoothed high frequency data. This is the average of the frequency content over the high frequency range for a specific pixel during a specific test. Each of the points in this range is subsequently given that value. The data smoothing procedure is typically implemented one time for each spectral data set. Figure 4 is representative of the smoothed data. It shows the same frequency content as the result shown in Fig. 3, but after the smoothing process is complete. Note that spectral peaks in Fig. 4 are much more evident than in Fig. 3.

Fig. 4

Grayscale spectral energy distribution after smoothing, using the data from Fig. 3 for data acquired on 13 November 2017

2.5.6 Remove white noise

The white noise is associated with the background electronic noise within the time series signal. The next step in the analysis is to subtract the white noise from the frequency domain data. The average value that replaces the grayscale spectral energy distribution data corresponding to frequencies greater than 95% of the Nyquist folding frequency during the smoothing process is considered the white noise energy content level. The same value of the white noise energy content level is subtracted from every spectral data point at each frequency within the smoothed frequency array. Note that each image pixel is associated with a unique frequency transform and therefore a unique white noise energy content level. Removing the white noise gives greater definition of the spectral variations on a log-log plot, especially at higher frequencies. Figure 5 shows a typical example of the smoothed data before and after the white noise is removed.

Fig. 5

Grayscale spectral energy distribution with and without removing the white noise. The data is from the pixel [736, 329] from data acquired on 13 November 2017. The spectrum is representative of data before and after removal of white noise

2.5.7 Ensemble averaging the frequency data

Grayscale spectral energy distribution results are also ensemble-averaged. Five pixels located within the same flow structure are selected for this purpose at specific positions. Grayscale spectral energy distribution data associated with the five pixel locations on the shock wave (from the test on 13 November 2017) are shown with respect to frequency in Fig. 6. Grayscale spectral energy distribution data for the five pixel locations are averaged together at each distinct frequency. The ensemble-averaged grayscale spectral energy distribution with respect to frequency for this data is then shown in Fig. 7.

Fig. 6

Grayscale spectral energy distributions with respect to frequency for the five pixel locations analyzed near the shock wave for data acquired on 13 November 2017

Fig. 7

Ensemble-averaged grayscale spectral energy distribution with respect to frequency, determined from Fig. 6 data, for pixel locations near the shock wave from data acquired on 13 November 2017

2.5.8 Shock wave position tracking

The MATLAB R2013a code is used to track the position of a shock wave in a time sequence of shadowgraph images. Shock waves generally appear as a dark line next to a bright line within shadowgraph images. Typically, the dark portion has the most contrast relative to the background. Hence, streamwise shock wave positions are represented by locations of the darkest pixels in a particular portion of the shadowgraph visualization images.

In order to track the shock wave position, the user first specifies the file path containing the time sequence of shadowgraph images. The file path and name of the first image in the time sequence is passed into the “imshow” function in MATLAB. This function displays the image in a MATLAB figure. The user must select a location near the shock wave in a particular shadowgraph visualization figure for this purpose. The Get Points function, “getpts,” determines the pixel coordinates at the location clicked by the user. Only the vertical coordinate is used. It is paired with all horizontal, i, pixel location coordinates, creating a horizontal line of pixels. Figure 8 shows a white horizontal line at the selected vertical pixel coordinate from data obtained on 05 April 2018.

Fig. 8

Shadowgraph flow visualization image, where horizontal white line indicates pixel locations where the grayscale data are analyzed to obtain results given in Figs. 9, 10, 11, and 13 for data acquired on 05 April 2018

The coordinates of each pixel along this line are input into the “imread” function. This “imread” function obtains the grayscale pixel value at every pixel along this line. The program determines the grayscale pixel value along the same horizontal line for each image in the time sequence array. An example of the variation of grayscale pixel values with respect to horizontal pixel coordinate, i, at one instant of time for one image, is shown in Fig. 9. This figure shows the grayscale pixel values for the last image in the time sequence associated with the 05 April 2018 test. The user must specify a region of streamwise, i, pixel coordinates where the shock wave is expected to appear. Doing this minimizes the possibility that the shock wave finding algorithm mistakes an image effect (not related to flow structure) for the shock wave. As an example, scratches on the side wall of the wind tunnel appear as dark spots within the shadowgraph images. A scratch is identified in Fig. 9, as are the edges of the viewing window. Microsoft Paint is used to determine the pixel range in which the shock wave is expected. The Minimum function, “min,” in MATLAB then determines the pixel with the lowest grayscale value within the specified region which is present along the horizontal line. This is done for each image within the time sequence. This pixel is the darkest pixel. As such, it represents the location of the shock wave. The horizontal pixel number, i, is saved to an array. The darkest pixel in the range between i = 500 and 800 is identified as the shock wave in Fig. 9.

Fig. 9

Grayscale value variation with horizontal pixel location of the last frame in the time sequence for vertical pixel 277 for data acquired on 05 April 2018

With the shock wave location array in hand, the average pixel number representing the average shock wave position is then determined. This average value is subtracted from each value in the time sequence so that a value of 0 corresponds to the average shock wave location. Negative values indicate the shock wave is downstream of its average position. Positive values indicate the shock wave is upstream of its average position. The physical location corresponding to the pixel location is then determined. A scaling factor is used to convert the pixel coordinate of the shock wave into a physical distance in milllimeters. The scaling factor is determined by relating the physical length of a piece of tape in shadowgraph images to the number of pixels along its length. The number of pixels is found using Microsoft Paint. The length is measured using Mitutoyo 500–196-30 Absolute Solar digital calipers. The streamwise shock wave position, relative to its average position, is determined as it varies with time. Time is the index of the time sequence array associated with each shock wave position divided by the sampling frequency. Figure 10 shows an example of the resulting streamwise shock wave position with respect to time. Figure 11 shows a subset of these data that further illustrates the motion of the shock wave. Figure 13 then shows the associated grayscale spectral energy of shock wave position with respect to frequency.

Fig. 10

Normal shock wave location (with respect to the time-averaged shock wave position) as a function of time for vertical pixel number 277 for data acquired on 05 April 2018

Fig. 11

A subset of the data from Fig. 10 showing the normal shock wave position (with respect to the time-averaged shock wave position) from 1.40 to 1.45 s, as a function of time for vertical pixel number 277 for data acquired on 05 April 2018

2.5.9 Correlations

The MATLAB code is also used to determine different spatial and time correlations as they vary with location, frequency, phase lag magnitude, and time lag magnitude. Filtered grayscale value time-varying data, as described earlier, are used as inputs for correlation determination.

2.5.10 Single point auto-correlation

An auto-correlation of grayscale data from a single pixel is computed as it varies with time. The auto-correlation function is defined by

$$ {ACF}_k=\frac{c_k}{c_o} $$

(10)

This gives the normalized correlation between the grayscale value time sequence data and the same data sequence at a time lag of k. The correlation is expressed by the equation

$$ {c}_k=\frac{1}{N}{\sum}_{n=1}^{N-k}\left({y}_n-\overline{y}\right)\left({y}_{n+k}-\overline{y}\right) $$

(11)

The magnitude of ACF_k given by Eq. (10) is normalized by the variance of the time sequence, which is determined using the equation given by

$$ {c}_o=\frac{1}{N}{\sum}_{n=1}^N{\left({y}_n-\overline{y}\right)}^2 $$

(12)

Large positive values of the ACF_k auto-correlation function indicate that the data are similar at both instances of time. Negative values outside of a confidence interval indicate that the data are strongly different, with less positive correlation. The confidence interval employed is 95%, as applied to the correlation to show time lags at which the auto-correlation is significant.

MATLAB 2013a code is used to determine the auto-correlation of the grayscale time sequence data associated with a pixel located on the shock wave. The time sequence data is an input into the “autocorr” function. The number of lags to compute is the second input into the “autocorr” function. For this input, a number one less than the number of data points in the time sequence is used. Therefore, the auto-correlation is calculated over a range of possible time lag values. The auto-correlation, the time lag, and the confidence interval bounds are outputs to the “autocorr” function. The “autocorr” function assumes a sampling frequency of 1 Hz for determination of time lag values, which is not correct for the collected data. The time lag that is output by the “autocorr” function is then divided by the sampling frequency to determine actual time lag values in seconds.

2.5.11 Squared coherence magnitude

Data associated with two separate pixel locations, each from a different region, are inputs to a magnitude squared coherence calculation. The magnitude squared coherence is the ratio of the cross power spectral density, P_y1y2, to the product of the power spectral densities of the two functions, P_y1y1 and P_y2y2, expressed as

$$ {C}_{y_1{y}_2}=\frac{{\left|{P}_{y_1{y}_2}\right|}^2}{P_{y_1{y}_1}\ast {P}_{y_2{y}_2}} $$

(13)

The power spectral density is represented by the equation

$$ {P}_{y_1{y}_1}=Y{(f)}^2 $$

(14)

The cross power spectral density is then given by

$$ {P}_{y_1{y}_2}={\sum}_{m=-\infty}^{\infty }E\left\{{\left({y}_1\right)}_{n+m}{\left({y}_2\right)}_n^{\ast}\right\}\cdotp {e}^{-j2\pi f\cdotp m} $$

(15)

Equation (15) represents the Fourier transform of the cross correlation of the time sequences of grayscale data associated with two different pixel locations, designated y₁ and y₂. Larger values of magnitude squared coherence indicate greater coherence between the grayscale data associated with the two pixel locations. Magnitude squared coherence, power spectral density, and cross power spectral density all vary with frequency. This is because unsteadiness is more prominent at certain frequencies than at others, and prominent frequencies change with spatial position within the flow field.

The magnitude squared coherence is determined using the “mscohere” subroutine in MATLAB. Filtered time sequence data from two separate pixel locations are inputs to the subroutine. The sampling frequency is also an input to the “mscohere” subroutine. This subroutine uses Welch’s Overlapped Segment Averaging procedure. Use of this procedure requires input information related to window type, window size, and size of overlapped region. The averaging scheme employs Hanning windowing. The window size is equivalent to one quarter of the length of the filtered time sequence array of grayscale data. This is the same as a quarter of the number of shadowgraph flow visualization images acquired and analyzed. For example, there are 21,620 images captured during the test on 05 April 2018, so each window consists of 5405 data points (21,620 divided by 4). The default size of the overlapped region is also required. This default size is half of the size of the window length, rounded down. Hence, the windows overlap by 2702 data points.

Magnitude squared coherence for five pairs of pixel locations within the same regions are ensemble-averaged. For example, data from five pixel locations in the downstream boundary layer are correlated, respectively, with data from five pixels locations on the shock wave. Those five coherence results are ensemble-averaged. Next, the resulting ensemble-averaged data are smoothed. However, smoothing begins at the sixth data point, not the twentieth. Table 2 summarizes this procedure.

Table 2 Magnitude squared coherence local data averaging procedures over different frequency ranges

2.5.12 Cross power spectral density, phase lag, and time lag

The cross power spectral density, which is given by Eq. (15), is used to calculate the phase lag and the time lag between time sequences of grayscale data associated with two different pixel locations. The time lag is a function of frequency. Perturbations occur at a location and propagate to other locations. Time is required for flow perturbations to travel, so the signals detected at two locations are similar, but differ by some time lag. All the perturbations do not originate in the same location, nor do they propagate along the same path at the same rate. Hence, the perturbations at unique frequencies have unique time and phase lags when considered at the same two pixel locations.

Time-sequences of grayscale values associated with individual pixels in different regions are inputs to the cross power spectral density, “cpsd,” function in MATLAB. For example, one way used for such correlation determination is the variable containing the data associated with the boundary layers is entered as the first parameter, and the variable containing the data associated with the shock wave or lambda foot is entered as the second parameter. The “cpsd” function also uses Hanning windowing. The window size is equivalent to one quarter of the length of the filtered time sequence array of grayscale data. The default overlap is used, which is half of the window length. The sampling frequency is the last input parameter in the “cpsd” function.

In regard to cross power spectral density results, the angle of the phase lag is calculated based on the real and imaginary parts of the complex result at each frequency. This is done using the negative of the “angle” function in MATLAB. The resulting phase angles are smoothed using the smoothing process given in Table 2. Phase angle results at different frequencies associated with five pairs of pixel locations from the same regions are ensemble-averaged.

Time lag magnitudes are determined by dividing phase lag values by associated frequency values, as given by

$$ \varDelta t=\frac{\phi }{2\pi f} $$

(16)

Although the phase and time lags are calculated at each discrete frequency between 0 and the Nyquist Folding frequency, the only meaningful phase lag results are at frequencies associated with strong coherence. Because the boundary layer signal is represented by the first term in the cross-power spectral density function, positive times and positive phase lag values indicate that the perturbation in the boundary layer occurs prior to the perturbation in the shock wave. Negative values indicate that the boundary layer disturbance occurs after the shock wave disturbance at the same frequencies. Note that the range of phase lag in this analysis is limited to values between -π and π. Spurious results may arise if the actual time lag corresponds to a phase lag that is outside of these limits (at a value of 3π/2, for example).

2.5.13 Spatial variations

The magnitude squared coherence and the time lag are determined for different pixel locations for certain frequencies. The magnitude squared coherence and the time lag are determined between pixels along a line, relative to a single pixel location. In this process, the magnitude squared coherence and time lag are calculated and smoothed, but are not ensemble-averaged. As such, the magnitude squared coherence and time lag results at a frequency of 40 Hz are considered with respect to spatial position. With this approach, the magnitude squared coherence and time lag at this frequency are determined for a single location associated with the normal shock wave, relative to other locations near and adjacent to the shock wave.

3 Experimental results and discussion

3.1 Test section flow characteristics

The main flow Mach number at the test section inlet is approximately 1.54. This value is determined from the inlet nozzle design, and is confirmed by measurements of local static pressure and local stagnation pressure. The associated test section mass flow rate is approximately 12.5 kg/s. Local boundary layer thickness is estimated to be approximately 10.5 mm on the top and bottom walls of the test section at a streamwise location which is just upstream of the leading edge of the shock wave holding plate. When steady flow conditions are established within the test section, the absolute pressure downstream of the pressure regulator is approximately 380 kPa. Static absolute pressure upstream of the shock wave in the test section is approximately 97 to 98 kPa, and the absolute pressure within the exhaust plenum is 130 kPa..

3.2 Overall flow structure

Figure 1 presents an instantaneous shadowgraph flow visualization image for the 05 April 2018 experimental result, illustrating the presence of a normal shock wave, a lambda foot, and a separated turbulent boundary layer located near the entrance of the lower flow passage, and an oblique shock wave system within the upper flow passage. Flow within this image is from right to left, with the shock wave holding plate evident within the left-hand side of the image. As such, the present experimental arrangement and resulting flow structure are similar to the ones described by Ogawa and Babinsky [16]. Flow visualization data are captured at an acquisition rate of 10.0 kHz during the test on 05 April 2018. A total of 21,619 images are collected, as steady flow conditions are maintained during this test. The frequency resolution is 0.4625 Hz. Each image size is 1024 × 512 pixels. To achieve the instantaneous result presented in Fig. 1, flow at the test section inlet must be uniform and steady, with relatively low freestream turbulence intensity, and no significant flow disturbances or non-uniformities. Note that a few oblique lines are evident within Fig. 1 on the right-hand side of the image, which are located just downstream of the test section inlet. These lines are due to the presence of Mach waves. According to Ogawa and Babinsky [16], “faint oblique lines represent Mach waves stemming from joints of parts or small disturbances on the tunnel walls, which have negligible effects on the flow.” Also evident in Fig. 1 are a few dark spots, which are due to imperfections on the outside surfaces of the side walls of the wind tunnel. These spots have no effect on the flow.

3.3 Shock wave Streamwise position grayscale spectral energy result

A Lagrangian approach is used to determine the spectral energy distribution associated with the streamwise location of the normal shock wave, as discussed earlier. At a vertical location marked by the white line in Fig. 12, the shock wave position changes approximately 6 mm in either direction, relative to its average position, due to shock wave unsteadiness. Note that the shock wave is positioned slightly downstream of its average position in this shadowgraph image. The resulting grayscale energy spectrum is shown in Fig. 13. This is compared with the ensemble-averaged grayscale spectral energy distribution (determined for five closely-located pixel locations) near the shock wave, shown in Fig. 14. The grayscale spectral energy distributions at five pixel locations are ensemble-averaged to generate these results. The region in which these pixels are located is denoted by a white rectangle in Fig. 12. Note that Figs. 13 and 14 show peaks at similar frequencies. For example, both data sets show spectral peaks at approximately 40 Hz and between 2 Hz and 9 Hz. Figures 13 and 14 also show that the spectral energy decreases significantly for frequencies greater than 200 Hz.

Fig. 12

Shadowgraph flow visualization image, where short horizontal white line indicates pixel locations where the grayscale data are analyzed to obtain results given in Fig. 13, and where white rectangle indicates pixel locations where grayscale data are analyzed to obtain the ensemble-averaged results given in Fig. 14, for data acquired on 05 April 2018

Fig. 13

Spectral energy distribution associated with the shock wave position (relative to average shock wave position) as a function of frequency for data acquired on 05 April 2018

Fig. 14

Ensemble-averaged grayscale spectral energy distribution for locations near the shock wave for data acquired on 05 April 2018

3.4 Auto-correlation function result

The auto-correlation is calculated for data associated with a specific pixel location. For the present result, a pixel location associated with the shock wave is considered. The auto-correlation of this time sequence data is shown in Fig. 15. Significant positive values are evident at 0.0 s, 1.4 s, and 0.75 s. Values for approximately 0.75 s indicate that the associated correlation is relatively weaker, since the auto-correlation magnitude barely exceeds the 95% confidence interval. Negative values of the auto-correlation indicate that the data are strongly dissimilar at time lags of 0.4 s and 0.9 s. Overall, the data in Fig. 15 show that gray scale time sequence data are roughly cyclical, with the most pronounced oscillation time period equal to approximately 1.4 s.

Fig. 15

Single point auto-correlation variation with time, determined from data from a single pixel near the shock wave. Horizontal blue lines indicate the range of a 95% confidence level relative to the zero value, for data acquired on 05 April 2018

3.5 Correlation function variations between shock wave and downstream boundary layer regions

Magnitude squared coherence and time lag variations are determined between data associated with the shock wave and the downstream boundary layer. The regions containing the locations of the pixels that are analyzed are shown by white rectangles in Fig. 16. Figure 17 shows the magnitude squared coherence of the time sequence grayscale signals in these two regions with respect to frequency. Associated data values evidence significant coherence between the shock wave, and downstream boundary layer regions, at frequencies of approximately 6 Hz, 20 Hz, 40 Hz, and 100 Hz, which correspond to respective Strouhal numbers of 0.00079, 0.00263, 0.00526, and 0.0132. Note that the Hanning windowing creates some fluctuations in the results at high frequencies, which are not representative of true coherence.

Fig. 16

Shadowgraph flow visualization image, where rectangles mark pixel locations in the downstream boundary layer and near the normal shock wave, where time-varying data are analyzed to obtain the results shown in Figs. 17 and 18, for data acquired on 05 April 2018

Fig. 17

Magnitude squared coherence variation with frequency of data associated with locations near the shock wave and in the downstream boundary layer, for data acquired on 05 April 2018

Time lag values from grayscale flow visualization results are shown in Fig. 18, also for perturbations between the downstream boundary layer and the shock wave. Perturbations of approximately 20 Hz, 30 Hz, and 70 to 100 Hz occur in the downstream boundary layer prior to the same frequency events in the shock wave. This is determined because values of time lag are positive for these different frequencies. At 100 Hz, perturbations in the downstream boundary layer occur approximately 1 ms before they do in the shock wave. However, for events at frequencies of approximately 5 Hz, 10 Hz, 40 Hz, and 60 Hz, time lag values are negative, indicating that the perturbations in the shock wave occur prior to the ones in the downstream boundary layer. Note that all of these time lags are approximately an order of magnitude longer than time lags associated with the estimated advection speed of the fluid.

Fig. 18

Time lag magnitude variation with frequency of data associated with locations near the shock wave and in the downstream boundary layer, for data acquired on 05 April 2018. Note that positive time lag values indicate the signal in the boundary layer (for a particular frequency) occurs prior to the signal (at the same frequency) which is associated with the normal shock wave

3.6 Correlation function and time lag results to illustrate spatial variations

The magnitude squared coherence and time lag at certain frequencies are determined for one location, relative to a range of other spatial locations. The locations are represented by the vertical white line within the flow visualization image of Fig. 19. Note that coordinate scale locations along this line are included. The magnitude squared coherence and time lag are calculated between these locations and a location on the shock wave, which is indicated by a white dot in Fig. 19. The origin is chosen to be the location of the pixel on the shock wave. The other pixel locations are measured relative to that location.

Fig. 19

Shadowgraph flow visualization image, where white dot and white line show locations where magnitude squared coherence and time lag are determined and analyzed at a frequency of 40 Hz to obtain the results shown in Fig. 20, for data acquired on 05 April 2018

Magnitude squared coherence values and magnitudes of time lag between the same shock wave pixel location and the pixels along the vertical line on the shock wave are given in Fig. 20 for a frequency of 40 Hz. This figure shows a magnitude squared coherence of 1.0, and a time lag of 0 s, at the location where y₁ is the same as y₂. This is because any signal is completely coherent with itself. For a frequency of 40 Hz, a peak in the coherence is evident approximately 12.7 mm beneath the location of y₂. Note that this location approximately corresponds to the intersection point between the normal shock wave and the two oblique shock wave legs of the lambda foot. Associated time lag results in Fig. 20 show that little coherence is generally present for locations where time lag values are notable. The most significant time lag value is negative and is evident at a location about 7 mm below the single shock wave pixel location. This means that events for the single shock wave position y₂ occur prior to events associated with the location 7 mm below, which is associated with y₁. With larger values of coordinate z, time lag values in Fig. 20 are then either near zero or show large amounts of data scatter.

Fig. 20

Magnitude squared coherence and time lag value, as they vary with location associated with y₁, at frequency of 40 Hz, where y₁ is located on a line and y₂ is positioned at one shock wave location, for data acquired on 05 April 2018. Positive values of time lag indicate that perturbations which originate at line locations occur prior to events which occur at the single shock wave location

4 Summary and conclusions

Within the present study, investigated are unsteady flow characteristics of a normal shock wave, a lambda foot, and a separated turbulent boundary layer within a unique research test section. The supersonic wind tunnel facility, containing this test section, provides a Mach number of approximately 1.54 at the test section entrance. A shadowgraph optical system is employed to visualize shock wave structure within the test section. Digitized shadowgraph flow visualization data are analyzed to determine shock wave structure and unsteadiness characteristics, including grayscale spectral energy variations with frequency, as well as time and space correlations, which give coherence and time lag properties associated with perturbations associated with different flow regions.

A Lagrangian approach is used to determine the spectral energy distribution associated with the streamwise location of the normal shock wave, The resulting grayscale energy spectrum is compared with an ensemble-averaged grayscale spectral energy distribution (determined for five closely-located pixel locations), which shows that both results evidence peaks at similar frequencies at approximately 40 Hz and between 2 Hz and 9 Hz. Auto-correlation function results show that the time sequence of grayscale values, for particular shock wave locations, is approximately cyclical, with a period of approximately 1.4 s. Two-point correlation functions, as they vary with frequency, indicate that perturbations at a frequency of 100 Hz (St_r = 0.0130) originate downstream of the shock wave and propagate upstream. Magnitudes of time lag between one stationary shock wave pixel location, and the pixels along the vertical line along and near to the shock wave, show that the most significant time lag value is negative and is evident at a location about 7 mm below the single shock wave pixel location. This means that events for the single shock wave position occur prior to events associated with the location 7 mm below that position.