Spatiotemporal Characterization and Suppression Mechanism of Supersonic Inlet Buzz with Proper Orthogonal Decomposition Method

Abstract

The buzz phenomenon of a typical supersonic inlet is analyzed using the unsteady Reynolds Average Navier-Stokes (RANS) simulation and proper orthogonal decomposition (POD) method. The dominant flow patterns and characteristics of the buzzed flow are obtained by decoupling the computed pressure field into spatial and temporal sub-parts based on the POD method. The supersonic inlet buzz phenomenon could be approximated as a product of decoupled temporal and spatial terms, and the one-dimensional (1D) mathematical model is therefore proposed. The standard deviations of the unsteady pressure fields from both the numerical simulation and the model prediction are compared. The limited discrepancy can be observed, and the good agreement validates the credibility of the proposed 1D model. The numerical simulation and the 1D model prediction are presented to explore the unsteady-jet control with a small perturbation. The results of the 1D model and the numerical simulation achieve good agreements with each other in terms of the overall trend. Finally, POD modal energy is employed to analyze the buzz suppression mechanism. When the jet frequency is identical to the dominant frequency of the buzz and the jet phase is opposite to the oscillation phase of the captured mass flow, the buzz suppression could be more efficient. The buzz suppression mechanism could be explained in two aspects. For one thing, the complex flow structure is suppressed and the first average modal energy in the inlet is increased. For another, the energy redistribution among each POD mode is achieved and the flow stability is gradually enhanced.

1. Introduction

The inlet buzz is an extremely unstable flow status observed in various supersonic and hypersonic air-breathing systems [1]. It usually occurs when the inlet operates in a subcritical flow regime and the entering mass flow drops under a certain threshold [2,3,4]. During inlet buzz, the inlet suffers violent and periodic shock system oscillations. As a result, the plenum chamber typically experiences severe pressure and mass flow rate oscillations, which may lead to thrust reduction and even fail of engine operation. Particularly, the accompanied periodic thermal and force loads could even damage the structure of the propulsion system [5,6].

Since the inlet buzz was first observed by Oswatitsch [7] in the experimental study of an asymmetrically configured missile inlet system, numerous remarkable experimental investigations have been carried out in the past decades. Thus, numerous cognitions concerning the mechanism of the inlet buzz have been obtained: (1) In a general sense, the inlet buzz could be classified into two types, namely the little buzz and the big buzz [3,4,8,9,10,11,12]. The former one is basically characterized by the flow oscillation with a small amplitude, while the latter features large-scale shock motion accompanied by enormous pressure oscillation and flow choke. (2) Corresponding to the two types of buzz, two prevailing triggering mechanisms, known as Ferri criterion [13] and Dailey criterion [14], have been proposed. However, recent progress by Tan et al. [11,12] revealed a possibility that both the little and the big buzz originate from the same source. (3) Although the origin of supersonic inlet buzz is still controversial, it is generally believed that the oscillatory flow is a self-excited oscillation in fluid flows and closely related to acoustic feedback mechanism [15,16].

Besides the above experimental research progress, on the theoretical side, Hankey and Shang [17] developed an acoustic resonance formula to predict the frequency of the inlet buzz and found that the result coincided with a number of inlet buzz flows. On the numerical side, it is difficult to achieve the high-precision numerical simulation for the high-frequency oscillatory flow of the big buzz [15,18]. Until 2014, Hong et al. [19] conducted the numerical simulation of the buzz evolution. However, advanced analysis methods for revealing the characteristics of the inlet buzz have not been established so far and no affirmation was made mathematically upon the flow characteristics. It impedes the in-depth understanding of the inlet buzz and the further exploration of the buzz suppression. For one thing, due to the insufficient survey points in the experiments, it is not convincing to perform a detailed description of the entire flow field. For another, it is difficult to acquire the accurate flow field of inlet buzz using computational fluid dynamics (CFD) tools. Moreover, the huge amount of data obtained from unsteady flow simulation usually makes it difficult to grasp any useful characteristics or dominant flow structures. Therefore, it is necessary to introduce an efficient method of data analysis, which aims at processing the massive transient flow data.

The proper orthogonal decomposition (POD) is a powerful and elegant method of data analysis focused on obtaining low-dimensional approximate descriptions of high-dimensional processes [20,21]. Since Lumley [22] first introduced the POD method into turbulence research, it has been successfully applied in various areas of fluid mechanics, such as coherent structures in a turbulent boundary layer [23], laminar separation boundary layer [24], flow around circular cylinder [25], Couette flow [26], wake flow downstream of a half-cylinder under the passive control [27], turbulent jet [28,29], vortex under synthetic jet control [30], pulsed jet in a blade cascade [31]. They all demonstrated that the POD method is a mature and promising tool to explore the essence of the complicated unsteady spatiotemporal flow field. In the supersonic and hypersonic field, the typical investigations include the supersonic turbulent flow over an open rectangular cavity [32] and the hypersonic separated flow over a double wedge with window proper orthogonal decomposition (WPOD) [33]. However, little work has been done in terms of applying the POD method to analyze the inlet buzz phenomena for the supersonic inlet.

Based on the above analysis, the objective of this work is to explore the spatiotemporal characterization of a representative supersonic inlet buzz phenomenon using the POD method, propose a new 1D mathematical model, and reveal the inherent pressure oscillation and buzz suppression mechanism. The rest of this paper is organized as follows: In Section 2, Nagashima et al.’s experimental work [2] was reinvestigated numerically. The numerical simulation procedure of the supersonic inlet buzz is explained, and preliminary results are validated by experimental data. Section 3 introduces the POD method and uses this method to analyze the spatiotemporal characterization of static pressure within the inlet. Section 4 discusses the 1D mathematical modeling of the static pressure oscillation of the inlet buzz. Section 5 analyzes the buzz suppression mechanism from the point of POD modal energy. Finally, Section 6 concludes this paper.

2. Numerical Simulation

2.1. Inlet Model Configuration

The axisymmetric supersonic inlet model is reproduced from Nagashima et al.’s experiments [2], as shown in Figure 1. It is composed of a cowl, a center body, and a plug. As tabulated in

Table 1. Geometrical configuration data of the inlet.

Parameter	Value
D, m	0.04
L, m	0.635
T.R.	0.67

, the total length from the cowl’s leading edge to the exit of the plenum chamber is 0.635 m. The diameter at the inner surface of the cowl wall is 0.04 m. The plug, attached at the rear part of the model, is used for controlling the inlet exit area to simulate backpressure variation from inlet downstream component. The throttling ratio (T.R.) defined by Equation (1) represents the area ratio of the exit to the minimum throat. The plug axial movement allows for adjusting of A_E to control the throttling ratio. In the subsequent study, T.R. will be kept at 0.67 constantly. Similar to Nagashima et al.’s work [2], four points on the center body and the cowl wall are selected to monitor the pressure variation during the simulation, as illustrated in Figure 1. P2 and P3 are located on the center body at the axial position of 0.044 L and 0.071 L, while P4 and P7 are located on the internal cowl surface at the axial position of 0.209 L and 0.986 L, respectively.

$$T.R.=\frac{A_E}{A^{*}}$$

(1)

2.2. Governing Equations

The unsteady Reynolds-averaged axisymmetric Navier-Stokes (URANS) equations are adopted as the governing equations. With regard to a two-dimensional axisymmetric flow, it can be written as follows:

$$\frac{\partial Q}{\partial t}+\frac{\partial F}{\partial x}+\frac{\partial G}{\partial r}+H=\frac{1}{Re}\left(\frac{\partial R}{\partial x}+\frac{\partial S}{\partial r}+E\right)$$

(2)

where the variables x and r are the axial and radial coordinate, respectively, and Re is the Reynolds number. The conservation variable Q is defined by

$$Q=\left[\begin{array}{c}\rho \\ \rho u\\ \rho v\\ e\end{array}\right]$$

(3)

in which ρ is the density, u and v are the velocity components in the axial and radial directions, and e is the internal energy. The convective flux terms, F and G, are expressed as

$$F=\left[\begin{array}{c}\rho u\\ \rho u^2+p\\ \rho uv\\ u\left(e+p\right)\end{array}\right]$$

(4)

$$G=\left[\begin{array}{c}\rho v\\ \rho uv\\ \rho v^2+p\\ v\left(e+p\right)\end{array}\right]$$

(5)

where p is the static pressure. The viscous fluxes, R and S, can be expressed as

$$R=\left[\begin{array}{c}0\\ {\tau }_{xx}\\ {\tau }_{xr}\\ u{\tau }_{xx}+v{\tau }_{xr}+q_x\end{array}\right]$$

(6)

$$S=\left[\begin{array}{c}0\\ {\tau }_{xr}\\ {\tau }_{rr}\\ u{\tau }_{xr}+v{\tau }_{rr}+q_r\end{array}\right]$$

(7)

The source terms H and E, for inviscid and viscous flow, respectively, are written as

$$H=\frac{1}{r}\left[\begin{array}{c}\rho v\\ puv\\ \rho v^2\\ v\left(e+p\right)\end{array}\right]$$

(8)

$$E=\frac{1}{r}\left[\begin{array}{c}0\\ {\tau }_{xr}\\ {\tau }_{rr}-{\tau }_{\theta \theta }\\ u{\tau }_{xr}+v{\tau }_{rr}+q_r\end{array}\right]$$

(9)

where τ is viscous stress, and q is heat flux.

We also have

$${\tau }_{xx}=\frac{2}{3}\mu \left(2\frac{\partial u}{\partial x}-\frac{\partial v}{\partial r}-\frac{v}{r}\right)$$

(10)

$${\tau }_{xr}=\mu \left(\frac{\partial u}{\partial r}+\frac{\partial v}{\partial x}\right)$$

(11)

$${\tau }_{rr}=\frac{2}{3}\mu \left(2\frac{\partial v}{\partial r}-\frac{\partial u}{\partial x}-\frac{v}{r}\right)$$

(12)

$${\tau }_{\theta \theta }=\frac{2}{3}\mu \left(2\frac{v}{r}-\frac{\partial u}{\partial x}-\frac{\partial v}{\partial r}\right)$$

(13)

$$q_x=-\kappa \frac{\partial T}{\partial x}$$

(14)

$$q_r=-\kappa \frac{\partial T}{\partial r}$$

(15)

where the viscosity coefficient µ and the thermal conductivity κ are calculated using Sutherland’s law.

In solving the flow equations, Roe’s method [34] in combination with the monotonic upwind scheme for conservation laws (MUSCL) interpolation [35] is adopted to obtain the inviscid fluxes. The turbulent flow is modeled by the Spalart–Allmaras (SA) model [36], of which the governing equations are discretized by a second-order upwind scheme. A dual-time implicit stepping method coupled with the lower-upper symmetric Gauss–Seidel scheme [37] for inner-iterations is used to simulate an unsteady flow. In addition, a no-slip adiabatic boundary condition is imposed on the solid walls. The convergence criterion for the calculation is that the residuals of continuous equations, momentum equations, energy equations, and SA equations are dropped below 10⁻³, or the residuals no longer decline as the iterations continue.

2.3. Computational Grid and Boundary Conditions

Figure 2 presents the computational domain and grid for the two-dimensional numerical simulation. The computational domain is divided into 19 blocks and composed of quadrilateral cells. In order to resolve the boundary layer accurately, the grid cells near the solid wall are denser. The non-dimensional size of the first grid point of the wall is guaranteed for y⁺ < 10. Starting from the first grid point, a constant stretching ratio of 1.2 is used in the boundary layer meshing.

The “pressure-far-field” boundary condition is given to sides AB and BC with the constant static pressure, static temperature, and Mach number. The “pressure-outlet” boundary condition is imposed on sides CD and EF, where the ambient static pressure is set. The sides AI and GH are both “axis”. The “wall” boundary condition is applied at the rest of the sides, as shown in Figure 2.

The start of the URANS simulation is based on the results of the steady-state RANS simulation at T.R. = 0.67.

Table 2. Input parameters for the proposed numerical simulation.

Parameter	Value
C0, m/s	330.2
L*, m	0.06
Ma∞	2
Re	107

lists the main parameters for the proposed numerical simulation. The stagnation speed of sound C₀ is 330.2 m/s, and the freestream Mach number is 2. The Reynolds number is 10⁷, with regard to the characteristic length (L*) of 0.06 m.

2.4. Verification of the Adopted Numerical Method

The computational error of unsteady numerical simulations highly depends on the grid resolution and the time-step size. Therefore, to solve the flow field accurately and efficiently, the grid refinement and the time-step sensitivity tests are carried out in advance. A sequence of three grids and four time-step sizes are examined, as exhibited in

Table 3. Grid refinement and time-step sensitivity test.

Grid	Number of Cells	Time-Step Size
Coarse	26,000	1.0 × 10−5
Medium	52,000	1.0 × 10−5
Fine	104,000	1.0 × 10−5
Medium	52,000	4.0 × 10−5
Medium	52,000	2.0 × 10−5
Medium	52,000	1.0 × 10−5
Medium	52,000	5.0 × 10−6

. The coarse, medium, and fine grids contain 26,000, 52,000, and 104,000 cells, respectively, whereas the grid points near the wall are identical in order to resolve the boundary layer.

By adopting the approach of Bonfiglioli and Paciorri [38], the convergence order of the numerical scheme is assessed with a posteriori error analyses. The L₁ norm of the global discretization error for density ρ is plotted against the characteristic mesh size (in 1og-1og scale), and its decay rate could be visually evaluated by the comparison with the first- and second-order slopes, as shown in Figure 3a. Actually, for a time-step size of ∆t = 1.0 × 10⁻⁵ s, the order of convergence of the medium grid and the fine grid is 1.45 and 2.04, respectively. Therefore, the numerical solution from the medium grid is in the asymptotic convergence range.

Besides the grid convergence, the time-step sensitivity is evaluated by comparing the computed buzz frequency with experimental data. Using the medium grid mentioned above, three additional time-step sizes of ∆t = 5.0 × 10⁻⁶ s, ∆t = 2.0 × 10⁻⁵ s, and ∆t = 4.0 × 10⁻⁵ s are selected to perform three computations, where 300 sub-iterations are conducted within each time-step. Examining the four computations with different time-step sizes, the numerical simulation yields a dominant buzz frequency of 360 Hz, 360 Hz, 360 Hz, and 345 Hz, respectively, as shown in Figure 3b. Since Nagashima et al.’s experimental frequency [2] is 360 Hz, the medium grid with a time-step of ∆t = 1.0 × 10⁻⁵ s is selected so as to ensure the numerical accuracy and computational efficiency.

To demonstrate the credibility of the numerical method in simulating the inlet buzz oscillatory process, qualitative comparisons are first conducted upon shock patterns on the center body and cowl region as shown in Figure 4. In Figure 4a,b, both the numerical and experimental schlieren images captured the strong bowl shock near the cowl tip region, however the cone shock induced by the center body is too weak to be visualized in experimental schlieren images but resolved in simulation. In Figure 4c, as the intense back pressure forces the bowl shock moving upstream and standing at the center body surface, both numerical and experimental results capture this flow pattern. For time-variant flow feature comparison, pressure histories of experimental and numerical results at P2, P3, P4, and P7 are plotted in Figure 5. It can be observed that numerical simulation can predict oscillatory pressure and agrees very well with experimental results, in both phase and amplitude. In general, the adopted numerical methodology is credible in predicting unsteady flow characteristics and reproducing transient pressure behavior as measured in wind tunnel experiments. Therefore, it is feasible to use the calculated results for an in-depth analysis of the inlet buzz phenomenon.

3. POD Analysis of the Inlet Buzz

3.1. Proper Orthogonal Decomposition Method

The POD algorithm was first proposed by Lumley [22] in the context of turbulent flows to identify coherent structures. A rigorous presentation of the algorithm can be found elsewhere [21,39], hence the following text briefly reviews the principal concept as applied to the analysis of supersonic inlet buzz. Although the POD method can be applied to any flow variables, this paper focuses on the pressure oscillation and the following introduction to POD is confined within the discussion of the pressure field.

The primary philosophy of POD is to determine an optimal set of orthonormal basis functions, so that a given spatiotemporal pressure field $p\left(x,t\right)$ can be decomposed as

$$p(x,t)={\displaystyle \sum _{j=1}^{\infty }{a}_j(t){\varphi }_j(x)}$$

(16)

where the functions ${\varphi }_j(x)$ are orthogonal basis functions relating to spatial characteristics, while the functions $a_j(t)$ are time-dependent functions relating to temporal characteristics. Note that functions $a_j(t)$ correspond only to ${\varphi }_j(x)$ and are independent on any other $\varphi$’s [20].

The aforementioned “optimal” means that a projection of p on to the first j basis yields the smallest error with respect to a mean-square sense, which involves maximizing the expression below

$$\frac{\langle {(p,\varphi )}^2\rangle }{\left(\varphi ,\varphi \right)}=\lambda$$

(17)

where < > is the time averaging operation for spatiotemporal data, (,) is the standard inner product, and λ is the eigenvalue. Maximization of the above quantity in Equation (17) is a classical problem of variational calculus [21,30,39,40] and leads to

$${\displaystyle \int \langle p\left(x\right)p^{\ast }\left(x^{\prime }\right)\rangle }\varphi \left(x^{\prime }\right)d{x}^{\prime }=\lambda \varphi \left(x\right)$$

(18)

where the superscript * denotes the complex conjugate. In order to address the above eigenvalue issue, all fluctuating pressure components from N snapshots of m grid points can be arranged in an m × N matrix (N « m)

$$p=\left[p_1,p_2,\dots ,p_N\right]$$

(19)

Applying the method suggested by Chatterjee [20], we can compute the singular value decomposition (SVD) of the transpose of the matrix p, which is of the form

$$p^T=U\mathsf{\Sigma }{V}^T$$

(20)

where the superscript T represents matrix transpose, $U$ is an N × N orthogonal matrix, V is an m × m orthogonal matrix, and Σ is an N × m matrix with all zero elements except the diagonal elements. Then, the eigenfunction ${\varphi }_j(x)$ with corresponding eigenvalue λ_j is obtained from matrix V^T and Σ, respectively. These eigenfunctions are ordered according to a decreasing energy manner, which can be represented by λ_j with the appropriate normalization. Additionally, modal coefficient $a_j(t)$ is determined by $U\Sigma$. Thereafter, based on the results of SVD and Equation (16), we can construct a series of POD modes and perform lower-rank approximations to the original data. For any j < N, to calculate a rank j approximation to $p^T$, one needs to replace $U$ and V with the matrices of their first j column and replace Σ by its first j rows and first j columns. For more details about matrix operation in the POD method, the reader is encouraged to read Chatterjee [20].

In POD analysis of the flow field, different POD modes represent various flow structures in the flow field, and the magnitude of eigenvalues indicates the portion of the energy captured by the corresponding POD modes. This means that certain POD modes with higher eigenvalue imply large-scale coherent structures that dominate the global flow field [41]. The spatial structure of each POD mode does not vary over time, while the temporal characteristics of the instantaneous flow field are embodied by the temporal variation of modal coefficients. Figure 6 illustrates the flow chart of the general procedure of POD analysis upon the static pressure field. As shown in Figure 6, the ultimate step of POD analysis is to rebuild the flow field using a limited number of lower-order modes.

It is worth mentioning that the procedure of subtracting the mean flow field before applying POD method is skipped in this paper for the reason that the average flow field may be influenced in the case of flow control discussed in the Section 5. This is different from the common practice of POD method.

In summary, the POD is a powerful method in terms of processing and analyzing enormously complicated data. Particularly, the POD procedure could identify the most energetic contributions, and obtain the spatial and temporal characteristics of the corresponding modes from a time-varying flow field. Consequently, the fundamental and essential features of the flow field could be described and analyzed by only capturing the first few POD modes.

3.2. Analysis Based on Numerical Results and POD

With knowledge of the main principles of POD, this section focuses on POD analysis of the spatiotemporal characterization of the static pressure field during supersonic inlet buzz. Here, we employ 400 snapshots for POD analysis from the numerical simulation results, which extend over seven buzz cycles.

Considering that the data of POD analysis comes from the results of the numerical simulation, the steady and unsteady flow characteristics are analyzed firstly. Figure 7 shows the average and transient pressure distribution from the numerical results. From it we can see the average pressure inside the inlet is basically constant, yet the transient pressure distribution is uneven.

To further illustrate the unsteady flow characteristic of the pressure field, the standard deviation σ is introduced to evaluate the degree of pressure fluctuations inside the inlet, which can be calculated by

$$\sigma =\sqrt{\left[{\left(p_1-\overline{p}\right)}^2+{\left(p_2-\overline{p}\right)}^2+\cdots +{\left(p_N-\overline{p}\right)}^2\right]/N}$$

(21)

where p represents the transient pressure, $\overline{p}$ represents the average pressure, and N is the number of samples. According to this equation, the magnitude of the standard deviation σ reflects the severity of the pressure oscillations. The larger the σ is, the greater the amplitude of pressure fluctuations is.

Figure 8 shows the standard deviation distribution inside the inlet. From Figure 8 we can see that: (1) At nearly x/L = 0.07 and x/L = 0.75, the standard deviation σ approaches to zero. It means that the transient pressure does not vary with time at these positions. (2) At nearly x/L = 0.35 and x/L = 1.0, the standard deviation σ is the maximum. It means that the amplitude of pressure fluctuations is most obvious. (3) There are two stagnation points at nearly x/L = 0.07 and x/L = 0.75, and wave peak/wave trough at nearly x/L = 0.35 and x/L = 1.0.

Figure 9 illustrates the percentage of each POD mode energy to the total energy for the 400 simulation cases. A total of 400 modes were extracted and ordered according to the magnitude of the corresponding eigenvalues. As can be seen, with the increase of the mode, the corresponding proportion of modal energy decreases remarkably, and the energy percentage of higher mode quickly approaches zero. The big differences in terms of energy ratio demonstrate that there are essential discrepancies between the top n modes and higher-order modes. For example, the first-order and the second-order POD mode possess 69.1% and 12.8% of the total energy, respectively, which are much more than other higher-order modes. For higher modes, due to their little share of captured energy (e.g., the proportion of energy occupied by the 10th mode is merely 0.567% of the total energy), it is negligible when conducting lower-rank approximations to the original flow field.

As discussed above, modal coefficients ($a_j(t)$) present alteration of mode energy intensity, which can be derived from the projection of the instantaneous flow field on each POD mode. Each modal coefficient indicates that the magnitude of energy captured by each mode in the transient flow field. Figure 10 depicts the evolution of several modal coefficients during a few buzz cycles. It shows that different from other higher-order modes, the first-order modal coefficient fluctuates around the value of—$22\times {10}^6$. With the largest absolute value, the first-order modal coefficient (a₁(t)) again implies the dominant role of the first-order mode in terms of energy. For the second-order mode, its modal coefficient fluctuation curve is similar to a sinusoidal curve, and the oscillatory amplitude is obviously greater than that of the first-order mode and other higher-order modes. Moreover, using a fast Fourier transform (FFT) analysis, the dominant frequency of the second-order modal coefficient (a₂(t)) oscillation is about 360 Hz, which is identical to the buzz frequency obtained from the experiment [2] and numerical simulation. Therefore, it is perceived that the second-order POD mode could reflect the transient pressure fluctuation component and buzz dominant frequency characteristics. Furthermore, the third/fourth/fifth modal coefficients with a relatively small fluctuation amplitude have a frequency-doubling characteristic of 360 Hz and, even for the 10th modal coefficient, the frequency-doubling characteristic is still existing. However, when it comes to higher-order modes, such as 20th/50th mode, their modal coefficients present obvious high-frequency and low-amplitude noise characteristics which might keep multiple frequencies. In general, the top 20 modes and other higher modes could be considered as a frequency-doubling mode zone and a noise mode zone, respectively.

In order to explore the spatial characteristics, it is feasible to approximately reconstruct the entire flow field by employing the POD modes. Figure 11a,b displays the instantaneous static pressure snapshots of the first and the second-order POD mode. The snapshot of the first-order mode is almost the same as the average pressure distribution shown in Figure 7. Considering the above description of the first-order modal coefficient and the spatial distribution of the first-order mode, it is considered that the first-order mode could reflect the average flow phenomenon and dominate the global flow field. Meanwhile, for the snapshot of the second-order mode, the spatial pressure distribution in the supersonic inlet is observed, which is similar to the transient pressure distribution shown in Figure 7. Correspondingly, Figure 11c shows the static pressure fluctuation curve of the second-order mode along the internal flow passage. There are a peak and valley in the curve, both of which represent the most violent pressure fluctuations inside the internal flowpath. Besides, as a whole, the curve is similar to a 3/4 wavelength sine curve. Comparing the spatiotemporal characteristics of the first-order and the second-order modes, the second-order mode could be more suitable to analyze the pressure oscillations of the supersonic inlet buzz.

4. Mathematical Modeling

Theoretical Analysis

The introduction of the POD method could obtain the dominant features in the flow field. Based on these dominant features, it is possible to simplify the flow phenomenon and get a simple mathematical description, which could promote the comprehension of the spatiotemporal characteristics of the inlet buzz. According to the above POD analysis results, the first-order mode reflects the average flow phenomenon while the second-order POD mode could reflect the transient pressure fluctuation component and buzz dominant frequency characteristics. Therefore, we model the transient spatiotemporal characteristics of the buzz phenomenon with particular attention given to the second mode.

Figure 12 shows the variation of static pressure fluctuation contour of the second-order POD mode in one buzz period (t from 3.5 T to 4.5 T) and the corresponding modal coefficient curve. At t = 3.5 T, the static pressure of upstream air is higher than the static pressure of downstream air inside the inlet, and at the same time, the modal coefficient corresponds to a maximum value. After that, the upstream high-pressure air propagates downstream, making the upstream pressure gradually decrease and the downstream pressure continually increase. Accordingly, the value of the modal coefficient decreases from the maximum to the minimum, as displayed at t = 3.75 T and t = 4 T. Then, due to the exit blockage and the reduction of upstream stored mass, the downstream high-pressure air in the plenum chamber propagates upstream. Consequently, the downstream pressure continually decreases, and the upstream pressure progressively increases, as presented on t = 4.25 T and t = 4.5 T. Correspondingly, the value of the modal coefficient raised from valley to peak. From t = 3.5 T to t = 4.5 T, an inlet buzz period is completed, and the flow field at t = 4.5 T recovers to be the same as t = 3.5 T. As time goes, the inlet buzz turns to the next period.

Referring to some interpolation method mentioned in the reference [42,43], such as radial basis functions (RBF) and two-level discretization method, the variation pattern of the second-order modal coefficient a₂(t) is processed and sinusoidal function is determined to described the a₂(t). Thereby, the second-order modal coefficient a₂(t) can be approximately written as

$$a_2(t)\approx A_1\sin (2\pi \cdot \frac{t}{T}-\frac{\pi }{2})$$

(22)

where A₁ is the oscillation amplitude of the second-order modal coefficient and A₁ is about $5.65\times {10}^6$.

Besides the 2D contour of the static pressure fluctuation as shown in Figure 12 and Figure 13, Figure 11c correspondingly presents the transient static pressure curve along the internal flowpath. As can be seen, the pressure at each location inside the inlet oscillates between its own maximum and minimum value over time within one buzz cycle. At t = 3.5 T, points B and D reach a peak and valley on the pressure curve, respectively. Then, the pressure at position B gradually decreases, while the pressure at location D continually increases. As a result, at t = 4 T, the pressure at points B and D reach a minimum and maximum value, respectively. Afterward, the pressure at point B continuously raises to the peak, whereas the pressure at point D progressively decreases to the valley, as shown at t = 4.5 T. Different from the periodical pressure oscillation at most locations inside the inlet, it is interesting that the pressure fluctuation at points A and C stay always zero.

The pressure curve at t = 3.5 T in Figure 13 is quite similar to a sine wave, while the temporal variation of the pressure curve is similar to a standing wave. Moreover, the fluctuation characteristics of the two are similar as well. Particularly points A and C correspond to the nodes on a standing wave, while points B and D correspond to the antinodes on a standing wave. Therefore, the oscillatory flow ought to link with the standing wave formed inside the inlet. This conclusion is also in agreement with previous speculations [2,17,19]. Inspired by the standing wave expression, the fluctuation component of the pressure inside the inlet should take a formulation as:

$$p^{\prime }\left(x,t\right)=p\left(x,t\right)-\overline{p}\left(x,t\right)\approx a_2\left(t\right){\varphi }_2\left(x\right)$$

(23)

where p(x,t) is the given spatiotemporal pressure, $\overline{p}\left(x,t\right)$ is the time-averaged pressure obtained by Reynolds decomposition. Correspondingly, the spatial variation of the pressure should also be a sine function. In addition, the wavelength λ and the length of inlet model L satisfy

$$L=\frac{3}{4}\lambda$$

(24)

Similarly, we have

$${\varphi }_2\left(x\right)\approx A_2\sin \left(2\pi \cdot \frac{x}{\lambda }\right)$$

(25)

where $A_2$ is the amplitude of the spatial fluctuation.

Therefore, according to Equations (22)–(25), the second-order POD mode approximation to the original static pressure field inside the chamber during inlet buzz can be manifested by the following mathematic formula

$$p^{\prime }\left(x,t\right)\approx a_2\left(t\right){\varphi }_2\left(x\right)=A_1{A}_2\sin \left(2\pi \cdot \frac{t}{T}-\frac{\pi }{2}\right)\sin \left(2\pi \cdot 0.75\frac{x}{L}\right)$$

(26)

where the product of the two coefficients $A_1{A}_2=p_{2nd-max}$ is the maximum fluctuation of the antinode B show in the Figure 13. The value of the $p_{2nd-max}$ is 2.27 times the $p_{\infty }$ and it could be obtained by reconstructing the second-order modal flow field. Its physical meaning reflects the oscillation amplitude of the inlet buzz.

According to Equation (26), we successfully model the complicated inlet buzz phenomena using a simple mathematical formulation. It should be noted that Equation (26) is an approximation to the original pressure field during the inlet buzz. In more detail, we use the second-order POD mode to approximate the original pressure field. With this proposed model, the inlet buzz phenomena can be easily analyzed.

In order to validate the proposed model in Equation (26), Figure 14 shows the comparison of σ of the static pressure between the original numerical simulation and the second-order POD approximation (Equation (26)). As illustrated in Figure 14, there are two valleys and peaks on both two curves, corresponding to the locations of the nodes and antinodes, respectively. Furthermore, the overall trend of the two curves is basically consistent. In general, the mathematical model developed with the POD method can describe the periodical oscillation of the static pressure field during the inlet buzz, and the discrepancy is attributed to the mere utilization of the second-order mode to approximate the original flow field.

5. Analysis of Buzz Suppression Mechanism

For one thing, the development of the POD mathematical model is to deeply comprehend the inlet buzz phenomenon by analyzing the characteristics which dominate the oscillation flow field. For another, it aims to provide recommendations for the control of supersonic inlet buzz. Therefore, based on the POD method, mechanism analysis for the preliminary exploration of buzz suppression is performed.

As a result of the interference of the shock wave/boundary-layer near the lip, the self-excited energy produced by moving shock wave allows some parameters, like mass flow rate and pressure, to vary periodically. Considering the above description, a small perturbation jet with different control parameters is introduced into the original flow field. Due to the difference of phase and frequency between the small perturbation signal and the moving shock wave, it is possible to achieve the supplement and absorption of the air when the mass flow rate captured by inlet decreases and increases.

In order to obtain the mathematical model for exploring the buzz suppression, a numerical simulation of a small perturbation jet is conducted firstly. The location of the jet slot is shown in Figure 15. The unsteady “mass-flow-inlet” boundary condition is imposed on the jet slot JK and could be expressed as:

$$m_{jet}=\psi \cdot {\overline{m}}_{out}\cos (2\pi \cdot f^{\prime }\cdot t+\mathsf{\Delta }\phi )$$

(27)

where $m_{jet}$ represents the mass flow rate of the unsteady jet, ${\overline{m}}_{out}$ represents the average mass flow rate at the inlet exit under the uncontrolled state, $f^{\prime }$ represents the jet frequency, and $\psi$ is defined as jet intensity coefficient and could be evaluated by the percentage of the unsteady-jet amplitude to the average mass flow rate at the inlet exit.

It is assumed that the introduced small perturbation signal and the original oscillation flow field conform to the same mathematical model. The perturbation pressure may be described as follow:

$${p^{\prime }}_{\mathrm{s}}\left(x,t\right)=\psi p_{2nd-\max }\sin \left(2\pi \cdot \frac{t}{T^{\prime }}-\frac{\pi }{2}+\mathsf{\Delta }\phi \right)\sin \left(2\pi \cdot 0.75\frac{T}{T^{\prime }}\frac{x}{L}\right)$$

(28)

where $\mathsf{\Delta }\phi$ is the phase difference between the small perturbation signal and the original flow field. $T^{\prime }$ is the cycle time of the small perturbation signal. $\psi$ represents the intensity coefficient of the small perturbation signal. Therefore, after introducing a small perturbation signal, the mathematical model of the buzz suppression could be obtained by adding a control term of the small perturbation signal and could be expressed as

$$\begin{array}{ll}{p^{\prime }}_t\left(x,t\right)\approx p^{\prime }+{p^{\prime }}_{\mathrm{s}}& =p_{2nd-\max }\sin \left(2\pi \cdot \frac{t}{T}-\frac{\pi }{2}\right)\sin \left(2\pi \cdot 0.75\frac{x}{L}\right)\\ & +\psi p_{2nd-\max }\sin \left(2\pi \cdot \frac{t}{T^{\prime }}-\frac{\pi }{2}+\mathsf{\Delta }\phi \right)\sin \left(2\pi \cdot 0.75\frac{T}{T^{\prime }}\frac{x}{L}\right)\end{array}$$

(29)

Figure 16 shows the comparison of control effects between the mathematical model and the CFD numerical simulation under different phases and amplitudes but the same frequency (360 Hz). In order to facilitate quantitative comparison, the overall control effect is evaluated by ${\sum }_{x=0}^{x=L}\sigma$, which represents the sum of the standard deviations of the entire duct and can measure the total intensity of the oscillation in the duct. The sum of the standard deviations along the flow path in different control states is normalized by that of the uncontrolled state. According to Equation (29), the control effect curves of the mathematical model and the numerical simulation are displayed in Figure 16. Under two control amplitudes ($\psi$ = 2%, 1%), the control effects decrease when the phase difference ($\mathsf{\Delta }\phi$) varies from 0 to π and then increase when phase difference varies from π to 2π. It is worth mentioning that the inlet buzz is suppressed when the phase difference is in the range of π/2 to 3/2π. Although the trends of the model and the numerical simulation are consistent, wide discrepancies still exist at the peak and the trough of the control effect curves and it may attribute to the simply linear superposition.

In order to correct the model of buzz suppression, the effect of the perturbation pressure should be taken into consideration. Referring to the harmonically excited vibration, the small perturbation signal could be regarded as excitation and might cause the resonance phenomenon under the frequency of 360 Hz. The significant difference in the fluctuation amplitude may result from the resonance phenomenon. Therefore, the resonance amplification factor M is introduced to describe the resonance amplitude and the modified model could be written as:

$$\begin{array}{ll}{p^{\prime }}_t\left(x,t\right)\approx p^{\prime }+& M\cdot {p^{\prime }}_{\mathrm{s}}=p_{2nd-\max }\sin \left(2\pi \cdot \frac{t}{T}-\frac{\pi }{2}\right)\sin \left(2\pi \cdot 0.75\frac{x}{L}\right)\\ & +M\cdot \psi p_{2nd-\max }\sin \left(2\pi \cdot \frac{t}{T^{\prime }}-\frac{\pi }{2}+\mathsf{\Delta }\phi \right)\sin \left(2\pi \cdot 0.75\frac{T}{T^{\prime }}\frac{x}{L}\right)\end{array}$$

(30)

The resonance amplification factor M could be written as:

$$M=\frac{1}{\sqrt{{(1-r^2)}^2+{(2\zeta r)}^2}}$$

(31)

where r is the frequency ratio of the small perturbation signal to the original flow field. $\zeta$ is the damping ratio of this flow field. An appropriate value is required to be given to the resonance amplification factor M. Under the condition of $\psi$ = 2%, the value of M is adjusted so that the peak point of the modified model on the vertical axis coincides with the peak of the CFD simulation, thereby determining the resonance amplification factor M as 2.78. Thus, the damping ratio is approximately determined as 0.18. The results of the modified mathematical model are shown in Figure 17. The two curves of $\psi$ = 2% achieve a good agreement with each other in the term of the overall trend and fluctuation amplitude. As for the case of $\psi$ = 1%, the coincidence between the modified model and CFD simulation is not as well as that of $\psi$ = 2%. The limited error could be owing to the rough solution method of the damping ratio. The damping ratio of the jet flow field typically is difficult to be measured for the reason that it might be affected by the intensity of the small perturbation signal.

In order to analyze the control mechanism, Figure 18 shows the comparison in the phase relationship among the total energy, first/second/third modal energy. The total energy and each modal energy in different control states are presented as relative energy, using that of the uncontrolled state as the denominator. From the perspective of the POD modal energy, when the phase difference is in the range of π/2 to 3/2π, the first average modal energy increases. But for the second dominant frequency mode (360 Hz) and the third double-frequency mode (720 Hz), their energy is reduced, which indicates that the buzz phenomenon and the corresponding complex flow structure is suppressed. In addition, the reason that the energy of the first mode is increased might be that the energy of the second/third mode transfers to the first average mode.

The different curves of control effects are obtained under the same conditions of jet intensity ($\psi$ = 2%) but different jet frequencies. As shown in Figure 19, the control effects in the cases of 180 Hz and 540 Hz are both extremely insensitive to the variation of phase and close to the uncontrolled state. The CFD result of the frequency characteristic basically coincides with what is predicted by Equation (31) and further verifies the reliability of the modified mathematical model. The frequency characteristics imply that the unsteady jet with different frequencies might affect different modes. The jet with a frequency of 180 Hz (or 540 Hz) and the jet of 360 Hz might affect the noise modes and the frequency-doubling modes, respectively.

Figure 20 displays the relative energy of all modes and first/second/third mode in the cases of 180 Hz, 360 Hz, and 540 Hz with the intention of comparing the frequency characteristic of energy. The different dot marks represent different phase states at the corresponding frequencies. Intuitively, in the 180 Hz and 540 Hz control states, the distributions of the total energy and each modal energy are similar and significantly concentrated, especially for the first-order modal energy (concentrated on about 100%) and the total energy (concentrated on about 103%). However, the distribution of relative energy under the jet frequency of 360 Hz is obviously dispersed. The total energy and the first/second/third modal energy are scattered above or under the uncontrolled state line, which behaves the phase characteristic of energy. In addition, in the 180 Hz/540 Hz control state, the relative total energy is about 103% while the energy of the first-order mode is equivalent to the uncontrolled. There is almost no control effect in this case. It indicates that evaluation of the flow control should focus on the first-order modal energy.

Figure 21 shows the cumulative energy ratio of the top n modes in the typical control states to that of the uncontrolled state. The cumulative energy ratio could be written as ${\displaystyle {\sum }_n\lambda {}_{\left(\mathrm{control}\right)}}/{\displaystyle {\sum }_n\lambda {}_{\left(\mathrm{no}\mathrm{control}\right)}}$. Under the same phase conditions ($\mathsf{\Delta }\phi$ = 1.00π), it is found that: (1) The modal energy of the top modes in the cases of 180 Hz and 540 Hz jet frequencies have almost no change, so the energy ratio curves overlap the uncontrolled state line in the several low-order modes. (2) The rapid-change range of modal energy extends to about the 80-order mode in the cases of 180 Hz and 540 Hz, and the growth slope is significantly smaller than the case of 360 Hz. Generally speaking, each POD mode might keep multiple frequencies. If a mode contains the same frequency component as the jet frequency, it would be affected by the perturbation jet. There is a little component of jet frequency in the frequency-doubling modes but a relatively large part of jet frequency in the noise modes. Owing to this fact, under the condition of non-dominant frequency, the jet produces an enhanced effect directly on these higher-order modes causing the energy of corresponding modes and the total energy both to increase.

When the jet frequency is identical to the dominant frequency of 360 Hz, as shown in Figure 21, it is observed that: (1) When the phase difference $\mathsf{\Delta }\phi$ = 0.00π, the cumulative energy of the top n modes increases sharply as the increasing of n from 1 to 20, and then grows extremely gentle over the high-order noise mode zone of approximately 20 to 400. (2) In the no control effect state ($\mathsf{\Delta }\phi$ = 0.50π), the energy ratios show a trend of decreasing at first and then increasing. The slight energy might be shifted from the frequency-doubling modes to the high-order noise modes under the influence of the unsteady jet. The total energy does not change much. (3) In the buzz suppression state ($\mathsf{\Delta }\phi$ = 1.00π), the declining extent of the frequency-doubling modal energy is much larger than that of the no control effect state; and the rising extent of the noise modal energy of the buzz suppression state is much less than that of the no control effect state. The effect of the unsteady jet on the dominant-frequency mode and the frequency-doubling modes determines the final control effect (enhanced, suppressed, or invalid), and the effect on the noise modes just causes a slightly increasing for the total energy.

In summary, the buzz suppression mechanism could be explained in two aspects. (1) The energy of the second/third mode is decreased, and the weakening of the dominant frequency and the frequency-doubling modal energy illustrates that the corresponding flow structure and oscillation phenomenon are suppressed. (2) The energy redistribution among each POD mode is achieved. The energy of the second/third mode and the other frequency-doubling modes tends to transfer to the first average mode, and the energy of the first mode is increased. Therefore, flow stability and buzz-suppressed effect are gradually enhanced.

6. Conclusions

The spatiotemporal characteristic of the supersonic inlet buzz was discussed using numerical and modeling methods. The static pressure field was used to study the oscillation pattern of the inlet buzz. The POD method was introduced to spatiotemporally decouple the static pressure oscillation characterization. The numerical simulation and the model prediction are presented for exploring the small perturbation unsteady-jet control. Finally, POD modal energy is employed to analyze the flow control mechanism. According to the results, some conclusions could be summarized as follows:

• The dominant flow patterns and characteristics of the buzzed flow could be obtained by the POD method from the energy ratio. Compared with higher POD modes, the first mode reflects the average flow phenomenon and dominates the global flow field. The second mode reflects the buzz dominant frequency characteristics. The buzz phenomenon could be modeled based on the spatiotemporal behaviors obtained by the second POD mode. The supersonic inlet buzz phenomenon could be approximated as a product of the decoupled temporal and spatial terms. The standard deviations of the unsteady pressure field from both the numerical simulation and the model prediction are compared. The limited discrepancy can be observed, and the good agreement validates the credibility of the proposed 1D model. The tiny error between them is due to the utilization of the second POD mode, which only approximates the pressure fluctuation.
• The model results and the numerical results for exploring the small perturbation unsteady-jet control achieve good agreements with each other in terms of the overall trend. When the jet frequency is identical to the dominant frequency of the buzz, the unsteady jet will have a significant suppression or enhancement effect on the dominant frequency and the frequency-doubling modes but little impact on the high-order noise mode. Through supplying and absorbing the air when the mass flow captured by the inlet decrease and increase, the buzz suppression could be more efficient. The buzz suppression mechanism could be explained in two aspects. (1) The weakening of the dominant frequency and the frequency-doubling modal energy illustrates that the corresponding flow structure and oscillation phenomenon are suppressed. (2) The energy redistribution among each POD mode is achieved. The energy of high-order frequency-doubling modes tends to transfer to the first average mode, and the energy of the first mode is increased. Therefore, flow stability and buzz-suppressed effect are gradually enhanced.