Unsteady Motion of Shock Wave for a Supersonic Compression Ramp Flow Based on Large Eddy Simulation

A large eddy simulation (LES) is conducted to investigate shock wave/turbulent boundary layer interaction in a 24° compression ramp at a high inlet Mach number of M a = 2.9 . The recycling/rescaling Method is used as the inflow turbulence generation technique. The shock wave structure and boundary layer flow in the interaction region are studied by flow visualization methods, such as vortex recognition and numerical schlieren. The distributions of turbulent kinetic energy and Reynolds normal stresses at different streamwise locations are analyzed. The results show that a strong anisotropy turbulent flow appears in the reattached boundary layer after the interaction of the shock wave. The large-scale unsteady motion of the separation shock wave is analyzed by using the intermittent factor and power spectrum. It is found that the shock moves around the averaged separation position, and the length scale is equal to 72% of the inlet boundary layer thickness. The power spectrum analysis reveals the existence of low-frequency instability in the separation region.


INTRODUCTION
The shock wave/turbulent boundary layer interaction (SWBLI) is one of the important physical phenomena in the supersonic flow, which contains complex aerodynamic and thermodynamic problems. The interaction can significantly alter the heat conduction characteristics and produce a strong pressure fluctuation. Although the 2D compression ramp is a simple geometric model in the SWBLI, however the flow phenomena contain boundary layer unsteady flow, separation, reattached flow, and turbulent fluctuation enhancement caused by a strong adverse pressure gradient. The transient flow field presents a highly three-dimensional state (Lee and Wang 1995), and the multiscale interaction produces complex flow structures (Wang 2015). Thus, the simple 2D model is the best model for the validation of the computational fluid dynamics (CFD) method.
In the past, extensive research studies have been conducted on the SWBLI by experimental and numerical methods, but there are still some issues to be explored further. Dolling (Dolling 2001) considered that the low-frequency unsteady flow should be studied for the future in detail. Settles et al. (Settles et al., 1979) conducted a surface oil flow experiment on a compression ramp with an angle of 24°. It was found that the pattern of the "node-saddle point" is alternately arranged near the reattachment line downstream of the ramp, which proves the existence of the streamwise vortex. Smits and Muck (Smits and Muck 1987) used a hot wire anemometer to measure the compression ramp flow and studied the effects of turbulent fluctuation enhancement with three different angles of the ramp. It was found that when the intensity of the shock wave is sufficiently strong, the unsteady motion of the shock wave will become important due to turbulent fluctuation enhancement. Ganapathisubramani et al. (Ganapathisubramani et al., 2007) used high-speed particle image velocimetry (PIV) technology to measure the compression ramp flow. It was found that the striplike structures with different momentum exist in the upstream boundary layer, and these strip structures cause a low-frequency pulsation of the separation bubble. Wu et al. (Wu et al., 2013) used nano-tracer planar laser scattering (NPLS) technology to study the laminar/turbulent SWBLI on the supersonic compression ramp; the overall flow field was analyzed and local fine structures were identified.
Yi Zhuang et al. (Zhuang et al., 2018a) have performed an experimental investigation on a compression ramp shock wave/ turbulent boundary layer interaction at Ma 2.83. The ice cluster-based planar laser scattering technique was applied to acquire high spatiotemporal resolution images at the center plane. Two-dimensional slices of the coherent vertical structure (CVS) were acquired and extracted from these images with a machine learning-based method. By comparing the features CVS acquired before and after the interaction, the evolution of CVSs in the SWBLI flow was analyzed.
In fact, a large number of experimental studies were carried out at higher Reynolds numbers. But Bookey et al. (Bookey et al., 2005) selected a low-density gas as the working medium to reduce the Reynolds number in the experiment and achieve the numerical simulation (DNS, LES) comparison with experiments. Adams et al. (Adams 2000) used the DNS method to study the supersonic compression ramp flow and found that the Reynolds shear stress after the SWBLI increases more greatly than the Reynolds normal stress. Loginov et al. (Loginov et al., 2006) used the LES method to study the statistical parameter distributions and fluctuation characteristics of the turbulent boundary layer in the compression ramp and compared numerical results with the experimental data. The flow field analysis showed that the large-scale threedimensional flow structure downstream of the ramp is the principal reason for the spanwise unevenness of the flow field. Wu et al. (Wu and Martin 2008) used the DNS method to study the SWBLI in the compression ramp, with the inlet flow parameters consistent with the experimental conditions given by Bookey et al. (Bookey et al., 2005). By observing the spatiotemporal evolution and using the correlation analysis of the flow field, it was found that the shock wave motion always lags behind the pulsation of the separation bubble. Therefore, a "feedback loop" model consisting of the separation bubble, shear layer, and shock wave system is proposed to explain the mechanism of the low-frequency unsteady motion of the shock wave.
Kenzo S.et al. (Sasaki et al., 2021) have investigated the mechanisms of low-frequency unsteadiness in an impinging shock wave/turbulent boundary layer interaction at a Mach number of Ma 2. The Strouhal number St fL/U ∞ was used for the space-time spectral analysis to identify the key features of the shock motion, where L is defined as the interaction length and f is the frequency of fluctuation. The dominant frequency in the vicinity of the shock was exhibited by the streamwise evolution of the pre-multiplied spectrum of pressure fluctuations. In the upstream boundary layer, the spectrum presents mainly high-frequency content (St > 1) linked to the incoming turbulent eddies. However, in the vicinity of the shock position, a low-frequency broadband range emerges and is centered approximately at St 0.03. As the author pointed out, the large gap between these frequency scales has been reported in previous investigations by Touber and Sandham (Touber and Sandham 2011) and Dupont et al. (Dupont et al., 2006).
In this study, a large eddy simulation is used to study the SWBLI phenomenon of the 24°compression ramp; detailed analysis and discussion are carried out to understand the unsteady features of the shock wave in the SWBLI. The remainder of this manuscript is organized as follows. Section 2 gives a description of the numerical methods and techniques used in this study, including the generation of the inlet turbulent boundary layer, geometric model, and mesh distribution. Section 3 compares the simulation results with experiments and numerical results in the literature to verify the reliability of the program, and the details of the flow field are discussed and analyzed. The last section gives some conclusions of this study.

COMPRESSION RAMP AND NUMERICAL METHOD
The work in this study is based on the NUAA-Turbo CFD solver developed by our research group. In the large eddy simulation (LES), the finite volume method is used and the ROE scheme is used for the evaluation of convective fluxes, the WENO_ZQ scheme (Zhu and Qiu 2017) with fifth-order precision for the interface reconstruction, and the sixth-order central difference scheme for the spatial discretization of viscous fluxes. Time discretization uses the Runge-Kutta method with a total variation reduction property of third-order accuracy (Shu and Osher 1989). The dynamic sub-grid scale model is considered, and the sub-grid scale viscosity coefficient is determined with the two consecutive filtering by the method of Germano et al. (Germano 1991). The inlet boundary layer thickness δ of the compression ramp is used to make the length scale dimensionless. The computation domain consists of two parts: the flat plane computation domain (auxiliary computation domain) and the 24°compression ramp computation domain (primary computation domain). The "recycling/rescaling" method was used to generate the turbulent boundary layer in the auxiliary computation domain and the computed turbulence information as the inlet boundary condition of the primary computation domain. In the auxiliary computation, the distance from the inlet plane "7.3δ" is set as the recycled plane. The schematic of the computation domain is shown in Figure 1. The coordinate system origin is located at the ramp, and the coordinate axes "x, y, and z" indicate the streamwise, spanwise, and normal directions respectively. The upstream and downstream lengths of the ramp are both 7.73δ, the spanwise width is 2.15δ, and the normal height along the wall is 5.23δ. The number of grid points in the three directions is "505 × 89 × 112". The grid is evenly distributed in the spanwise and refined along the streamwise direction at the ramp and in the normal direction to guarantee z + ≈ 1 in the first layer of the grid near the wall. The flow field of the recycled plane in the auxiliary computation domain is extracted as the inlet boundary condition for the primary computation. It is specially noted that when the primary/auxiliary computation domain uses grids of different resolutions, the process of flow field extraction needs to interpolate the variables. For details, this method can refer to Zhong et al., (2021).
For the primary computation domain, the upper boundary of the computational domain and the outlet are set to the subsonic outlet boundary condition. The wall condition of the non-slip isothermal is adopted to the wall, and the wall isothermal temperature is 307 K. The spanwise boundary adopts periodic boundary condition, and the inlet turbulent boundary conditions are dynamically given by the auxiliary computation. For the auxiliary computation, the boundary layer conditions are consistent with the main computation, and the flow Mach number is Ma 2.9, the flow static temperature is 108.1 K, and the flow density is 0.074kg/m 3 .
In order to verify the reliability of the LES software, the numerical results of the compression ramp will be compared with the experimental results under the same inlet flow conditions.    (Tong et al., 2017) and the experimental data of Bookey et al. (Bookey et al., 2005) are also presented for the purpose of comparison. By solving the position where the averaged skin friction coefficient C f is zero, the averaged separation point x sep −2.7δ and the re-attachment point x rea 0.8δ in the ramp flow are obtained. Overall, the computed pressure distribution and friction coefficient distribution are in good agreement with the experimental values. At upstream of the ramp, it agrees well with the experimental values and DNS results. But it is slightly higher downstream of the ramp. The reason is that the predicted separation region size is slightly smaller and the separated shear layer completes the reattachment process in advance downstream of the ramp; thus, the wall pressure increases rapidly and the range of the pressure platform is shortened. Figure 3 illustrates the instantaneous flow structure in the compression ramp computation domain. Among them, the quasiorder vortex structure of the turbulent boundary layer is displayed using the Q criterion, and the streamwise direction velocity is used for coloring. The translucent gray surface in the figure is the dimensionless pressure isosurface p/(ρ ∞ U 2 ∞ ) 1.7, which is used to represent the three-dimensional structure of the separation shock wave. As can be seen in the figure, the flow field in the interaction region has significant three-dimensional characteristics. When a large vortex structure passes through the root of the shock wave, it breaks under the effect of turbulent fluctuation enhancement at the ramp. At the same time, the shock wave deforms at the root, and the shock surface wrinkles along the spanwise direction. While away from the interaction region, the shock wave keeps still the typical two-dimensional structure characteristic.
In order to further illustrate the interaction process between the shock wave and turbulent boundary layer, a slice contour in three directions is used to show the flow field details inside the interaction region, as shown in Figure 4. The streamwise and spanwise directions show the numerical schlieren diagram, which is used to show the shock wave, turbulent boundary layer, and flow structure in the separation bubble. It can be seen that the large-scale structure in the interaction region breaks into smallscale structures, which is well agreed with Zhuang Y.et al. (Zhuang et al., 2018b) experimental observation. But, on the slope downstream of the ramp, a large-scale quasi-order structure in the boundary layer is re-established. The instantaneous streamwise direction velocity contour (y/δ 0.02, approach laminar sublayer) is shown on the slice parallel to the wall. It can be seen that a velocity strip structure alternately arranged along the spanwise direction is observed both upstream and downstream of the ramp in this region. In addition, the flow scale in the separation bubble along the spanwise direction is very different, indicating that the separation bubble is a multiscale flow structure.

Turbulent Fluctuation
When the flow passes through the interaction region, the turbulent fluctuation is significantly enhanced due to the strong inverse gradient pressure in this region. Figure 5 is a contour of the time-space averaged turbulent kinetic energy. The definition of turbulent kinetic energy is as follows: k (ρu′u′ + ρv′v′ + ρw′w′)/2ρ ∞ U 2 ∞ ). It can be seen from the figure that the intensity of the turbulent kinetic energy after the interaction region rapidly increases and reaches its peak value downstream of the ramp. The increase of the turbulent kinetic energy means significantly that a large number of small-scale structures are generated after the SWBLI. Figure 6 presents the variation of Reynolds normal stress R ii at different locations in the interaction region. It can be seen that near the ramp (x 0.66δ), the three Reynolds stress components produce a sharp peak value. As the separated boundary layer reattaches, the Reynolds stress continues to decrease and tends to return to the pre-interaction state at (x 2.68δ). It is noted that the distribution of Reynolds stress R 22 in the boundary layer is relatively smooth and appears as a bimodal distribution near the ramp.

Unsteady Motion of Shock Waves
In order to study the shock wave unsteady motion, intermittent factors are usually used to measure the range of shock wave motion along the streamwise direction. The intermittent factor is defined as the time when the instantaneous wall pressure somewhere in the streamwise direction is greater than the given threshold takes up the proportion of the total flow field time. The computation formula is as follows (Dolling and Or 1985): where p w represents the instantaneous wall pressure; p w and σ(p w ) represent the average wall pressure of the inlet boundary layer and the standard deviation of the wall pressure, respectively, and the sum of the two is set as the threshold for the computation of the intermittent factor.  Figure 7 exhibits the local distribution of the intermittent factor in the compression ramp. It is known from the definition of the intermittent factor that in the undisturbed upstream boundary layer and downstream region of the ramp, the intermittent factor is always equal to 0 (or 1). However, in the vicinity of the average separation point, the separation shock wave shows a strong unsteady flow characteristic, which is represented by the intermittent factor: 0 < λ < 1. Wu et al. (Wu and Martin 2008) considered that the motion of the shock wave is closely related to the pulsation of the separation bubble. For comparison, the upper and lower limitations of the given intermittent factor are 98 and 4%, respectively, to estimate the streamwise direction range of the shock wave motion as follows: −3.2δ ≤ x ≤ − 2.48δ . It can be seen that the flow movement of the separation shock wave in the intermittent region is carried out around the average separation point(x −2.7δ). By computation, the length of the intermittent region in this simulation is L i 0.72δ which is slightly smaller than that in the DNS result of the study by Tong et al. (Tong et al., 2017L i

Görtler Flow Vortex
From the abovementioned analysis, the supersonic compression ramp flow has significant three-dimensional characteristics. Especially on the slope downstream of the ramp, the flow field shows a strong unevenness in the spanwise direction. This phenomenon is closely related to the presence of Görtler flow vortex (Dolling and Or 1985;Grilli et al., 2013). Figure 8 shows a time-averaged skin friction coefficient contour. The white dashed line represents the location of the ramp. The solid black line in the figure is used to show the separation and reattachment positions of the boundary layer near the ramp defined by the time-averaged friction coefficient equal to 0. It can be observed that the skin friction coefficient of the inlet boundary layer represents a strip pattern alternately arranged which is consistent with the spanwise distribution of the time-averaged separation line. This indicates that the distribution of the separation line is, in a great measure, determined by the skin friction coefficient of the inlet boundary layer. Due to the lowfrequency pulsation of the separation bubble, the average reattachment still exhibits strong transient flow characteristics. When the separated boundary layer reattaches, the skin friction coefficient distribution appears as an obvious "V" structure, which is caused by the presence of Görtler vortex pairs. The existence of this "V" structure was also confirmed by Fu-Lin et al., (2016).
The formation of Görtler vortices is related to the curvature of the concave surface. In order to describe the spatial development of the Görtler vortex in detail, Figure 9 shows the flow field details on six different sections downstream of the average separation    Figure 9A, in which the background is spacetime-averaged numerical schlieren. Figure 9B shows the streamwise direction vorticity ω x contour, and streamline distribution on six sections, where the vorticity contour has a range −0.5 ≤ ω x ≤ 0.5. In sections A and B, the streamlines in the near wall region are gradually bent, and a smaller streamwise vortex is formed. When the flow develops to section C, it can be clearly seen in this figure that there is a large pair of vortexes in the streamwise direction (B) and a small pair of vortexes in the streamwise direction (A), and the entire boundary layer has a large amount of vortexes in this direction. With the separated boundary layer reattach in the sections D, E and F, the larger pair of the streamwise vortexes gradually occupies almost the entire span width (about 2δ) and appears as a classical kidney eddy structure. Figure 9B illustrates that it can be estimated that the center of the Görtler vortex pair is located at y 1.2δ. In general, the position of the Görtler vortex in the spanwise direction changes constantly with time due to the turbulent flow (Floryan 1991), but Figure 10 shows the time-averaged skin friction coefficient along the span in different sections, including the inlet section of the ramp computation domain and the A, C, D, E, and F sections in Figure 9 so that the variation of the spanwise distribution of the skin friction coefficient with the development of the flow field can be presented. The dotted line in the figure indicates the timespan-averaged result at this section. At the inlet section, it can be found that the skin friction coefficient along the spanwise direction is not evenly distributed, and a low skin friction coefficient region appears near the spanwise center. Sections A, C, and D are all located in the separation region, so the friction coefficients of these three locations are close to 0. In addition, due to the shock wave/ boundary layer interaction and Görtler instability, the skin friction coefficient of these three sections along the spanwise direction shows strong fluctuation. Especially, an obvious low friction coefficient region is formed between y/δ 1~1.5. In sections E and F, the amplitude of this low-friction coefficient region is further increased, and a V-shaped distribution is formed. The fluid near the wall converges between the pair of Görtler vortexes and moves toward the outer layer of the boundary layer, locally resulting in a lower friction coefficient. In section F, the skin friction coefficient at the spanwise position y 1.2δ reaches a minimum, which is exactly the same as the center position of the Görtler vortex pair. In fact, by observing the friction coefficient distribution at the six different sections, it can be found that the low-friction coefficient region of the inlet section has a high degree of coincidence with the downstream section in the spanwise position. This phenomenon indicates that the spanwise position of the Görtler vortex downstream of the ramp is likely to be affected by the spanwise friction coefficient distribution at the inlet section.

Low-Frequency Instability
In order to study the low-frequency motion of the separated shock wave, it is necessary to collect the signal of numerical pressure for spectrum analysis. Therefore, 528 wall pressure probes are arranged along the flow direction in the middle section of the calculation domain. The arrangement of these probes is consistent with the grid distribution, that is, one probe is arranged in the center of each wall grid cell. It is worth noting that the signal acquisition process is carried out after the SWBLI flow field is fully established. In order to ensure sufficient time resolution, data recording is carried out every ten time steps in the large eddy simulation, and the physical time interval between two adjacent recording points is Δt 0.01δ/U ∞ . According to the Nyquist sampling theorem, the highest physical frequency that can be captured is f s 1/(2Δt) 50U ∞ /δ, which is far greater than the characteristic frequency U ∞ /δ of the turbulent boundary layer. Because the low-frequency instability of the SWBLI has a wide frequency range, in order to accurately capture the lowfrequency characteristics, the signal acquisition time should cover at least one low-frequency instability period. Table 1 shows the relevant information of pressure signal acquisition in this large eddy simulation. Figure 11 presents the variation of wall pressure with time at two different positions upstream of the corner (the corner position is x 0 ) in the center plane, and the wall pressure is nondimensionlized by using the inlet pressure P ∞ . The blue solid line represents the pressure signal graph collected at the inlet Compared with the inlet, the pressure signal at the average separation point has larger amplitude of fluctuation due to the interaction between the shock wave and turbulent boundary layer. In addition, from the red solid line, the pressure signal at the separation point not only includes the high-frequency fluctuation but also has the lowfrequency fluctuation component as shown in Figure 11. Next, we will make a strict quantitative analysis of the wall pressure fluctuation signal from the perspective of spectral analysis. The Welch method (Barbe et al., 2010) is used to segment the discrete pressure signal, aiming to obtain a smoother and less variance power spectrum density (PSD) curve. In this study, the collected pressure signals are divided into three segments, and the coincidence rate between the segments is 50%. The Hanning window is used to add windows for each segment to improve the variance performance. It should be noted that in the Welch method, the more segmented the data, the smoother the power spectrum curve and the smaller the noise, but at the same time, the resolution of the power spectrum will be affected. Therefore, in the process of data segmentation, we must consider the balance of noise and resolution in the power spectrum curve. Figure 11B shows the PSD distribution curve corresponding to the abovementioned two wall pressure signals, where the abscissa is the dimensionless frequency proposed by Dussauge et al. (Dussauge et al., 2006), which is defined as St fL sep /U e . It can be seen from the figure that the peak value of pressure fluctuation energy in the upstream boundary layer is located in the high-frequency region, while the peak value of pressure fluctuation energy at the average separation point is located in the low-frequency region. Figure 12 shows the weighted power spectrum density (WPSD) expressed as WPSD(f) fpPSD(f)df. From this figure, we can see more clearly the frequency components of pressure signals at different flow direction positions. For the upstream boundary layer, the fluctuating energy is mainly distributed between 10 0~1 0 1 , and the peak frequency is about St 3, which is consistent with the characteristic frequency of the fully developed turbulent boundary layer. For the signal at the average separation point, there is also a wide small peak area in the high-frequency region, which is similar to the WPSD curve obtained in the upstream turbulent boundary layer. At the same time, the red curve shows a more obvious energy peak appearing in the low-frequency region, and the corresponding characteristic frequency is St 0.04, which is completely consistent with the result (0.02-0.05) summarized by Dussauge et al. (Dussauge et al., 2006). In addition, Wang Bo et al. (Wang 2015), Kenzo S. et al. (Sasaki et al., 2021), Tong et al. (Tong et al., 2017), Touber and Sanham (Touber and Sandham 2009), and Pasquariello et al. (Pasquariello et al., 2017) obtained similar results in their respective numerical simulations.
In order to study the distribution of low-frequency instability characteristics in the whole field of the compression corner, Figure 13 shows the distribution of WPSD on all pressure probes. The white solid line is used to represent the position of the average separation point x sep and the average reattachment point x rea and the white dotted line is used to represent the position of the corner. Once again, it can be clearly seen that the turbulence generation technology in this study does not introduce any lowfrequency energy at the inlet. By observing the whole flow field, it is found that this low-frequency instability mainly exists near the average separation point, while in other locations of the flow field, the energy of pressure fluctuation almost exists in the highfrequency range. Especially when the boundary layer is reattached in the downstream, its power spectrum distribution is almost the same as that in the upstream, and the reattached boundary layer shows a process of flow recovery. This study considers that this lowfrequency motion is caused by the unsteady motion of large-scale separation bubbles. Obviously, the scale of the separation bubble is larger than that of any flow in the boundary layer.

CONCLUSION
In this study, a large eddy simulation (LES) is conducted to investigate the shock wave and turbulent boundary layer interaction in a 24°c ompression ramp with an inflow Mach number of Ma 2.9 and compared with experimental and numerical simulation results in the literature. The shock wave/turbulent boundary layer interaction (SWBLI) was analyzed from two aspects: time domain and frequency domain. The main conclusions are as follows: The shock wave/turbulent boundary layer interaction represents some flow structure characteristics. When the large scale vortex passes the root of the shock wave, the shock wave surface wrinkles due to the intermittence of the large-scale vortex in the spanwise direction. Under the strong inverse gradient pressure of the compression ramp, the large-scale vortex breaks into a small-scale vortex, the turbulent boundary layer manifests a strong anisotropy characteristic, and the turbulent kinetic energy rapidly increases in the outer layer of the boundary layer downstream of the ramp.
In the region of interaction, the separated shock wave represents an unsteady motion along the streamwise direction. This phenomenon is successfully captured by the intermittent factor. The range of the shock wave motion is −3.2δ ≤ x ≤ − 2.48δ, indicating that the unsteady shock wave motion is around the average separation point.
The LES simulation represents a pair of streamwise vortices, called Görtler vortex, occupying almost the entire span width about 2δ on the compression ramp. It is found that the spanwise skin friction coefficient distribution at the inlet of the turbulent boundary layer not only determines the spanwise separation location of the turbulent boundary but also makes the downstream Görtler vortex to be fixed at a certain spanwise position, thus forming a time stable flow structure.
The low-frequency instability in the SWBLI is successfully captured by using the power spectrum analysis method. The corresponding characteristic frequency is St 0.04 near the average separation point. The spectrum analysis indicates that the shock wave motion is related to the pulsation of the separation bubble, while in other locations of flow field, the energy of pressure fluctuation represents almost exclusively the high-frequency characteristics of the fully developed turbulent boundary layer.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.