Characterization of Supersonic Compressible Fluid Flow Using High-Speed Interferometry

This paper presents a very effective interference technique for the sensing and researching of compressible fluid flow in a wind tunnel facility. The developed technique is very sensitive and accurate, yet easy to use under conditions typical for aerodynamic labs, and will be used for the nonintrusive investigation of flutter in blade cascades. The interferometer employs a high-speed camera, fiber optics, and available “of-the-shelf” optics and optomechanics. The construction of the interferometer together with the fiber optics ensures the high compactness and portability of the system. Moreover, single-shot quantitative data processing based on introducing a spatial carrier frequency and Fourier analysis allows for almost real-time quantitative processing. As a validation case, the interferometric system was successfully applied in the research of supersonic compressible fluid discharge from a narrow channel in a wind tunnel. Density distributions were quantitatively analyzed with the spatial resolution of about 50 μm. The results of the measurement revealed important features of the flow pattern. Moreover, the measurement results were compared with Computational Fluid Dynamics (CFD) simulations with a good agreement.

However, the field of the transonic or even supersonic flow generated in wind tunnel facilities is very specific, and the current state of the art of interferometric techniques in this field allows only their limited use. Environmental harshness, particularly vibrations, the lack of an optical table, temporal flow instability, and high fluid density gradients leading to complex interference patterns often make the quantitative evaluation of interferometric data very difficult and cumbersome. This is why the use of brute force interferogram analysis [43] or capturing interferograms for qualitative analysis and flow visualization in wind tunnels is still common [22,[44][45][46][47][48][49].
In this paper, we introduce a robust, portable, high-speed and high-resolution interferometric method that can be easily applied in a wind tunnel facility, providing very accurate results with high spatial resolution in almost real time. The method utilizes a high-speed camera allowing for short exposures to avoid the blur of interference patterns typical for unstable flow, as well as to minimize environmental disturbances. Fiber optics ensure the high compactness and portability of the interferometric system, and fast single-shot quantitative data processing based on spatial carrier frequency and Fourier analysis is employed.
The method is intended for the investigation of high-speed flow in planar blade cascades, namely, for studying dynamic events during transonic and supersonic blade flutter. As a first step, it was successfully tested in the case of a supersonic highly underexpanded flow at the outlet of a narrow channel by measuring the density distribution of the flow. The present interferometric technique has proven to be an excellent tool to meet these objectives. In addition to a description of the principles, experimental arrangement and results, numerical simulations showing the high reliability of the measured data are presented.

Interferometry and Fluid Flow
Interferometric techniques generally provide whole-field, non-invasive and highly accurate measurements. These properties can be advantageously used for measuring and visualizing such flow phenomena that affect light wave [25]. The amplitude of the propagating light wave through the transparent/semi-transparent medium under investigation is not significantly affected, while the optical phase is [50]. These phenomena are called phasesensitive, implying that when a change takes place in the medium, it affects the refractive index of the medium ∆n, which as a consequence introduces phase deviation ∆ϕ: with dl denoting the differential distance along the optical path between points L 1 and L 2 , and λ stands for the wavelength of the laser light. The x,y coordinates define the position within the Cartesian coordinates system. The general solution of (1) leads to tomographic techniques [51] that are difficult or cumbersome to apply in real-time, dynamic measurements. However, there are two assumptions that allow a significant simplification of (1): (i) wind tunnel experiments are designed in order to keep the refractive index along the z-axis constant, i.e., the optical axis, and (ii) the bending of the rays due to a spatial variation in the refractive index can be neglected. Assuming (i) and (ii), Equation (1) can be simplified to: where L denotes the length of the object in the z-axis direction (i.e., wind tunnel testing section width). The refractive index change distribution is usually not the main concern. An important equation in the optical study of fluid flow is the Gladstone-Dale equation, which links the change of refractive index ∆n to the density ρ as: where K represents the Gladstone-Dale constant and ρ re f is the reference density. The Gladstone-Dale equation is almost independent of pressure or temperature under moderate physical conditions [50]. For a case of an isentropic compressible flow of ideal gas, the measured outlet density ρ can be further linked with the isentropic Mach number by the following equation: where ρ in is the inlet air density and κ = 1. 4 denotes the heat capacity ratio of an ideal diatomic gas. In addition, other quantities, such as velocity distribution, pressure distribution, temperature distribution, etc., can be analyzed quantitatively by interferometry. The investigation within this paper leads to the characterization of the air flow in/behind a narrow channel (see Figure 1). Air enters the test section at ambient conditions, and thus the total pressure p in is equal to the barometric pressure. The static pressure at the outlet of the narrow channel p out is measured by a pressure tap on a wall in the settling chamber using a pressure transducer (from Huba Control AG, Würenlos, Switzerland). For the case of an isentropic compressible flow of ideal gas, the density distribution ρ M in the settling chamber can be calculated as: where In (6), T 0 denotes the ambient air temperature and R = 287. 1 J/(kg·K) is the specific gas constant. The isentropic Mach number (4), however, can provide relevant information only in flows where the drops in total pressure and total density throughout the flow field are relatively small, e.g., the flow in blade cascades. This does not apply to the investigated phenomenon, and therefore the isentropic Mach number was not evaluated.

Spatial Carrier Interferometry
Equations (2)-(6) provide the relation between the phase change ∆ϕ and other quantities, such as density ρ M ; however, it is not clear yet how the phase change ∆ϕ is obtained.
The interference phase ϕ can be retrieved from an interference pattern created on a charge-coupled device (CCD), or a complementary metal-oxide semiconductor (CMOS) chip, by the superposition of a reference and an object beam. In (7), A is the additive component and B is the multiplicative component. The interference phase ϕ is coded in a cosine modulate fringe pattern. There are several ways to retrieve the interference phase [52]; however, for the characterization of fluid flow in wind tunnels, where real-time and high-resolution measurements are required, most of the interferogram processing techniques fail. An option to tackle the issues is spatial carrier interferometry [52], which introduces small angles of incidence θ x , θ y in the x, y directions on the camera sensor in the reference beam. This results in further modulation of the interference pattern with cosine fringes pattern of constant spatial frequencies k x0 = 2π λ sinθ x and k y0 = 2π λ sinθ y , respectively. Let us rewrite (7) using a complex exponential including the carrier frequencies k x0 and k y0 : where with the spatial frequency coordinates k x , k y is composed of a central DC (direct current) term and two conjugated components located symmetrically from the center of the spectrum. The roof symbolˆdenotes Fourier spectrum and the superscript symbol * denotes complex conjugation. The introduced angles θ x , θ y play a major role in the clear separation of all spectral componentsÂ,Ĉ andĈ * in the Fourier domain and must be sufficiently high to avoid any overlap. On the other hand, angles cannot be too high, in order to meet the Nyquist criterion.
Once the spectral components are separated, using a bandpass filter around the specific spatial frequencies, it is possible to filter out the desired spectral componentĈ containing the interference phase ϕ. By filtering only the desired component, the filtered spectrumĈ F (subscript F denotes filtered component) is no longer Hermitian, so its inverse Fourier transform will contain non-zero both real and imaginary parts, and hence the interference phase can be calculated as: The phase along with other quantities is stored in a computer's memory as a matrix of numbers. Let us assume a measurement at a steady/reference state without the presence of fluid flow. Quantities measured in the reference state can be denoted using a subscript 0, e.g., ϕ 0 for the interference phase. Such a measurement carries information about interferometer optical aberrations, including aberrations introduced by, e.g., optical windows of the wind tunnel experiment. Other measurements performed with the presence of a phenomenon (denoted with subscript 1) can be related to the reference state measurement in order to suppress undesired optical aberrations. The phase change ∆ϕ occurring in (1) can be calculated as: It is worth noting that spatial carrier interferometry could also be called off-axis image plane digital holographic interferometry. As phenomena in fluid flow do not lead to speckles, and there is also no need for the numerical propagation/focusing/reconstruction of wave fields typical for digital holography, we tend to call it interferometry.

Materials and Methods
Interferometry-based research of a fluid flow in a narrow channel including a jet flow at its outlet was carried out at the Laboratory of Internal Flows of the Institute of Thermomechanics of the Czech Academy of Sciences. The laboratory is equipped with a modular in-draft wind tunnel facility that enables low mass flow rates. The vacuum tank is connected to the test section through a pipe with a quick-acting valve and a nozzle that offers the possibility of changing the pressure inside the testing section from 0. 01 MPa to values close to those of the ambient atmospheric pressure. The three-dimensional computer-aided design (3D CAD) assembly of the narrow channel is shown in Figure 2. The experimental arrangement was designed in order to maintain a constant refractive index along the z-axis, i.e., the optical axis. Two parallel walls made of stainless steel constitute the flow channel. The test section has optical windows on the walls, allowing the passage of light beams for optical measurements. The channel height (along the y-axis) was set to 3 mm. The channel is 100 mm long (x-axis) and the width (z-axis) of the channel is L = 100 mm to achieve a sufficient aspect ratio to eliminate secondary flow effects and to provide the sufficiently long optical path necessary for the adequate sensitivity of the interferometry. In order to obtain quantitative values of measurements, the ambient atmospheric pressure (P in ) as well as the outlet static pressure in the settling chamber of the wind tunnel (P out ) were measured.
The experimental setup, consisting of the narrow channel experimental section and a Mach-Zehnder interferometer, is illustrated in Figure 3. A laser beam was generated by a pigtailed single-frequency distributed feedback laser (LAS) of wavelength λ = 773 nm and fiber output power of 38 mW. The laser beam was split using a fibersplitter (FS) into the reference wave and the object wave. The reference wave was guided by an optical fiber to the non-polarizing beamsplitter (NBS), while the free-space object wave passed through a collimating imaging lens L0, which sent light through the measuring area (MA) within the wind tunnel (WT), before it passed through the lens L1. Lenses L1 and L2 acted as a beam expander with magnification of Mag ≈ 0. 6 as well as an imaging system. Both reference and object waves were recombined by the NBS and collimated by L2. The object wave impinged the digital camera sensor (CAM) normally, while the reference wave angles were θ x = 0. 34 • , θ y = 0. 30 • , respectively. The angles of the reference wave were adjusted by a slight translation of the reference beam's fiber ferrule (FF).
The intensity interference pattern (interferogram) of both superposed beams was captured by a Phototron FASTCAM Mini WX100 high-speed camera with a resolution of 2048 × 2048 pixels (10 µm × 10 µm pixel size) and a frame rate of 1080 fps. The exposure time of 10 µs was set in order to avoid the blurring of the interference pattern while maintaining sufficient brightness (signal). In the first step, a reference interferogram at a steady state (P out = P in ) was captured and processed, and the filtered complex field C F0 was stored in a computer. The same procedure was repeated after a change in the pressure inside the testing section (P out = P in ), resulting in air flow. In order to obtain enough data for statistical analysis, a sequence of M = 50 interferograms was captured, of which one is presented in Figure 4a. The amplitude Fourier spectrum of the interferogram in the logarithmic scale is shown in Figure 4b. The termsÂ,Ĉ,Ĉ * are well separated by the spatial carrier frequencies f x0 = k x0 2π = sinθ x λM = 12.3 mm −1 and f y0 = 10.8 mm −1 , respectively. The black circle represents the filtering bandpass Hanning window, with bandwidth f B = 20 mm −1 , applied in order to retrieve the complex field C Fm , where m = 1, 2 . . . M is an integer denoting the frame number and M is the number of captured interferograms. Such bandwidth results in a spatial (lateral) resolution of about 50 µm. It is important to note that the applied filtering must remain unchanged from that in the reference interferogram processing. The optical phase change computed by (11) is wrapped within the [−π, π] radians interval due to the harmonic nature of the light wave-see Figure 4c. Therefore, a spatial unwrapping was applied with the starting point in the position of the Pitot probe, i.e., far from the channel without the influence of the flow, obtaining a sequence of unwrapped phase fields ∆ϕ m -see Figure 4d. In order to suppress the random time fluctuations of the flow, the mean of the unwrapped phase was computed: The mean unwrapped phase change distribution ∆ϕ was used to calculate the density distribution as described in Section 2.1.
In addition to the interferometric measurements, computational fluid dynamics (CFD) simulation software was used to cross-check both the measurement and the simulation. The numerical computations were performed using Ansys Fluent 2021R1 commercial code. This software package uses the finite volume method to solve the governing Navier-Stokes equations for compressible fluid flow. The problem was simulated in two dimensions (2D) using the Reynolds-averaged Navier-Stokes equations (RANS) approach to close the system of equations. The flow of air as an ideal gas was simulated as fully turbulent using the two-equation k-ω shear stress transport (SST) turbulence model. The computational mesh consisted of approximately 150, 000 cells. The mesh sensitivity analysis was performed using two finer meshes, yielding a relative error of the mass flow rate of 0. 3 % for the coarse mesh, which is very acceptable. The computational domain was assumed to be symmetrical around the channel axis with the no-slip boundary condition at the channel walls. The static pressure according to the desired pressure drop across the channel was prescribed as the outlet boundary condition. The following boundary conditions were used at the inlet: (a) constant total temperature T 0 = 292. 67 K, (b) uniform total pressure distribution p 0 = 98. 05 kPa, (c) constant turbulence intensity Tu = 1. 2 % and (d) constant eddy viscosity ratio µ t µ = 10.

Results and Discussion
The results of the interferometric measurement of air flowing out of the narrow channel, and their comparison to the results of CFD simulation, are presented in this section.
The fluid flow in the narrow channel was investigated in five different regimes defined by different pressure gradients through the channel. The magnitude of the pressure gradient was controlled by a valve that determined the outlet pressure in the settling chamber downstream of the narrow channel. Both inlet and outlet pressure were measured, yielding the pressure ratio η = P out P in . The inlet pressure (the ambient atmospheric pressure) was constant during all measurements: P in = 98. 050 kPa. Parameters under different regimes are summarized in Table 1. Other constants considered for further flow analysis and CFD simulations are the ambient temperature T 0 = 295. 9 K, the Gladstone-Dale constant K = 0. 000225 m 3 kg −1 , and the specific gas constant of dry air R = 287. 1 J/(kg·K).
The density distributions computed using (5), together with the results of the CFD simulation, are presented in Figure 5. It can be seen that the interferometric measurements and CFD simulations are in a very good agreement. Small discrepancies can be attributed to minor asymmetries in the geometry of the experimental setup. The boundary conditions of the CFD simulation could be modified in order to derive an even better agreement; however, this is not this paper's objective.
The measured distributions in Figure 5 have significantly dynamic ranges of values. Therefore, some faint details of the flow might be difficult to observe. In order to visualize the high-frequency spatial variations in the flow, the slopes of the measured density distribution in the horizontal (x) and vertical (y) directions ∂ρ M (x, y)/∂x, ∂ρ M (x, y)/∂y were calculated, and these, together with the magnitude of the slope (∂ρ M /∂x) 2 + (∂ρ M /∂y) 2 , are shown in Figure 6. Such results could be directly compared to the Schlieren imaging techniques, which are sensitive to the ray deflections caused by refractive index gradients ∂n(x, y)/∂x and ∂n(x, y)/∂y, respectively.
Five different regimes with a gradually decreasing ratio η were measured. In the first regime, #1 (η = 0. 560), the supercritical pressure ratio across the channel results in the subsonic discharge of air into the settling chamber. Thus, a subsonic exit jet can be seen in Figure 6. There is also visible compression just before the exit of the channel. This indicates that the regime is close to aerodynamic choking, and is approaching supersonic velocity within the channel.
In regime #2 (η = 0. 481), the exit jet is supersonic with a shock wave present. Thus, the flow is already aerodynamically choked. In this regime, the flow is most likely overexpanded (i.e., the pressure at the exit cross-section of the channel is lower than the outlet pressure in the settling chamber), since two oblique shocks are formed at the exit of the channel first.  In the supersonic regime #3, as the ratio further decreases (η = 0.400), the jet flow's inner pattern starts to take an "X" or a "diamond" shape, which is characteristic of a moderately underexpanded jet flow. Here, expansion waves form at the edge of the channel, causing the flow to turn outward in a plume, creating a characteristic underexpanded jet flow pattern. Expansion waves "reflect" from the jet's constant pressure boundary as compression waves, which eventually form oblique shock waves. A periodical series of compression and expansion waves appears, until the jet is dissipated and static pressure within the jet flow reaches equilibrium with the surrounding pressure (backpressure).
In the supersonic regimes #4 (η = 0.287) and #5 (η = 0.243), the flow first expands through expansion waves, before turning into oblique shocks and Mach disks. Such highly underexpanded flow regimes exhibit a "barrel-shaped" pattern that differs in the shapes and lengths of the shock cells with different η. The flow patterns repeat until mixing dissipates the flow.
The typical resolution of the optical phase measured by interferometry is better than 1/100 of 2π, leading to a refractive index resolution of about 4 × 10 −7 in this particular case. However, interferometric measurements in a wind tunnel are burdened with errors. The first group of errors is directly due to interferometric sensing-the measurement of the optical phase. When retrieving ∆n from (2), there is an error in laser wavelength u λ , thickness u L , and phase noise u ∆ϕ . Let us define the refractive index combined error as: The phase noise for our measurement is estimated to be u ∆ϕ = 0. 07 rad. This value was derived from experimental data as the standard deviation of an unwrapped phase map (in different regions) after low-order polynomial subtraction (i.e., high-spatial frequency components were considered as noise). The wavelength of the used laser source depends on the temperature and current of the laser diode. The temperature and the current are measured/controlled; however, there is not a specific relation between these and the output wavelength, thus there is a wavelength error, which can be estimated to be u λ = 0. 5 nm. The error of thickness, u L = 0. 1 mm, comprises, on one hand, the limited accuracy of channel manufacturing, and on the other hand the fact that rays are deflected due to refractive index non-homogeneities. Such ray deflection is not considered in (2). Moreover, in locations with a large refractive index gradient, the ray deflection misplaces information geometrically (up to 200 µm), leading to so-called mapping errors. Such distortion is below ten pixels in the areas with the largest gradients, otherwise it is negligible. By substituting in (13), the maximal error is u ∆n ≈ 1 × 10 −6 and the relative maximal error results in 100 u ∆n ∆n ≈ 2%. Assuming the refractive index error u ∆n , along with the 2% relative error of the pressure/temperature sensors and 1% error of the Gladstone-Dale constant in (5), the density error yields u ρ M ≈ 0. 015 kg/m 3 , or relatively 100 Another source of error is usually related to the repeatability of measurements. However, the repeatability of measurements in this particular case can only be determined under a steady state (without flow) due to the unstable behavior of the flow. Steady state repeatability has a negligible effect when compared to the stability of the flow. Therefore, time stability was examined. As aforementioned, we captured 50 frames during one measurement sequence, taking 50 ms, and used averaged phase maps ∆ϕ for further processing-see (12). However, the standard deviation, for each pixel of (x, y) coordinates can also be computed from the data set. The standard deviation can be used to calculate a time stability map τ(x, y) == 1 − 3σ(x, y) that (assuming the normal probability distribution) equals one in pixels of completely stable flow, while it approaches zero at locations of random signal. The time stability maps of flow in different flow regimes are shown in Figure 7. It is clear that in the narrow channel, the flow is very stable, but becomes much less stable in shear layers and regions of flow dissipation. This effect is more apparent at high flow velocities with a more turbulent character.

Conclusions
This paper presents a very effective technique for sensing and investigating compressible fluid flow. The developed interferometric technique is very sensitive and accurate, yet easy to use, even under pretty harsh environmental conditions (vibrations, flow instability, unclean room) typical for aerodynamic labs. The interferometer employs a high-speed camera, fiber optics, and readily available "of-the-shelf" optics and optomechanics. The interferometer consists of an illumination and a sensing unit, each weighing less than 5 kg. An optical fiber delivers the reference wave from the illumination to the sensing unit, while the object wave emitted from the illumination unit propagates through a measured volume into the sensing unit. Both waves are superimposed, and an interference pattern is captured by a high-speed camera. The construction of the interferometer together with the fiber optics ensures the high compactness and portability of the system. Moreover, single-shot quantitative data processing based on introducing a spatial carrier frequency and processing in the Fourier domain allows for almost real-time quantitative processing, which is also suitable for very fast-evolving phenomena.
The interferometer was successfully applied for the study of a fluid flow in a narrow channel including a free jet at the channel exit. The experimental section was built in a wind tunnel, allowing for a high pressure difference and thus high flow velocity through the channel. The flow was investigated under five different regimes, ranging from subsonic to supersonic flow with a highly underexpanded jet. The density distributions were quantitatively analyzed with a spatial resolution of about 50µm. The results of the measurement were able to reveal important features of the flow pattern, such as shock waves and slip lines.
The measurement results were compared with CFD simulations, showing good agreement in terms of the flow pattern similarity (e.g., size of the shock cells) as well as the absolute values of density. An error analysis revealed the error range in units of percent. It can be concluded that the method constitutes a promising tool for studying highly dynamic events, even during transonic and supersonic blade fluttering.