1 Introduction

In compressible cavity flows, self-sustained instabilities arise due to the impingement of the large-scale well-formed structures of the shear layer growing between the two corners of the cavity, resulting in a pressure feedback mechanism which induces the far field noise in the cavity. Large oscillation amplitude exceeding 160 dB may occur even in subsonic cavity flows. Hence, suppression of the pressure oscillation is quite necessary, and the closed-loop control is an effective method^[1-3]. Before applying the closed-loop control algorithm, there is a need to design a reduced-order model (ROM) of the fluid and acoustics. Moreover, cavity flow modes can be affected by the cavity geometry (e.g., the cavity length-to-depth ratio)^[4-7] and the inflow conditions (e.g., the Mach number and the thickness of the incoming boundary layer)^{[5, 8-9]}. Application of direct numerical simulation (DNS) and large eddy simulation (LES) to systematically study the long time dynamics of cavity flows often poses difficulties in terms of computational resources. Today, the ROM based on proper orthogonal decomposition (POD) and Galerkin projection has been widely used as a predictive tool to model the flow dynamics^[10].

The POD/Galerkin ROM has been applied to various problems in different fields: the study of dynamics of coherent structures in the near wall region^[11], in a mixing layer^[12], or in the wake of a cylinder^[13]. All these studies deal only with the incompressible regime, and the ROM based on incompressible Navier-Stokes (N-S) equations is quite simple. For compressible flows, the problem is a little bit more complex because one needs to consider the thermodynamic variables. Rowley et al.^[14-16] proposed the use of the simplified isentropic N-S equations to obtain a simple dynamic system for Mach 0.6 cavity flow. They used both scalar POD and vector POD to build ROMs for the purpose of bifurcation analysis and control. The results showed that ROMs including viscous terms worked better, and the dynamics can be captured only for a short time but deviated for a long time. Note that the isentropic ROM is a valid physical model only for flows at low or moderate Mach numbers. Based on this ROM, many researchers have achieved the feedback control in subsonic cavity flows^[17-20]. However, the feedback control for supersonic cavity flows is a much tougher task, and building of the ROM for supersonic cavity flows is profitable to the feedback control in experiments or numerical simulations.

The major objective of this study is to build a ROM (named the approximate full ROM) with visualized weighting inner product of the kinematic variables (the velocities $\boldsymbol{u}$) and the thermodynamic variables (the density $ρ$ and the temperature $T$) for the POD, which is able to be applied to high Mach number flows. Here, we derive new governing equations and build a new ROM, which is successfully applied to supersonic cavity flows. By comparison with the DNS results, it is shown that the proposed ROM can capture accurately the dominant flow dynamics (e.g., the Rossiter modes). In addition, the present ROM is robust.

2 Numerical methods 2.1 Proper orthogonal decomposition

POD is a mathematical technique used in many applications, including image processing, signal analysis, and data compression. The central idea consists in determining an optimal subspace of reduced dimension, in the sense that the error due to the projection on this subspace is minimal. In fluid mechanics, POD is firstly used to analyze the turbulent boundary layer^[21]. The POD modes are most easily computed using the method of snapshots, which was introduced by Sirovich^[22] for high spatially resolved data.

For a given set of data , where is the mean of the snapshots, and $H$ is a Hilbert space with inner product $\langle\cdot, \cdot\rangle$, we want to determine a subspace $S$ of dimension $M$ ($M < N$) such that the error is minimized. Here, $\|\cdot\|$ is the induced norm of the Hilbert space $H$, $P_S$ denotes the orthogonal projection onto the subspace $S$, and $E (\cdot)$ is an average operator over $k$. Note that minimizing the error is equivalent to maximizing . We write the orthogonal modes $\boldsymbol {\phi}$ on the subspace $S$ as a linear combination of the snapshots ,

(1)

Solving the optimization problem with variational method leads to an eigenvalue problem,

(2)

where $\boldsymbol {U}$ is an $N \times N$ matrix with , $\boldsymbol {c}$ is the eigenfunction, and $\boldsymbol \lambda$ is the eigenvalue. Usually, we reorder the eigenvalues as $\lambda_1>\lambda_2>\cdots>\lambda_N$, the ratio of eigenvalue $\lambda_i$, and the summation of them $\sum\limits_i \lambda_i$ denotes the normalized energy of corresponding reordered eigenfunction $\boldsymbol {c}(:, i)$. The POD modes are obtained from the eigenfunctions by

(3)

The coefficients of the POD modes are given by

(4)

The snapshots $\{\boldsymbol {q}\}$ can be reconstituted through the $M$ first POD modes,

(5)

2.2 Galerkin projection

In a dynamical system, the mean data $\overline{\boldsymbol {q}}$ and the orthonormal POD modes $\boldsymbol {\phi}_k$ are only functions of space, and the coefficients of the POD modes $\boldsymbol {a}_k$ are only functions of time. Hence, the solution $\boldsymbol {q}=\boldsymbol {q}(\boldsymbol {x}, t)$ can be expressed in term of the $M$ first POD modes as

(6)

and the time partial differential equation is

(7)

The Galerkin projection is a special case of weighted residual methods^[23]. The following dynamical system is presented:

(8)

where $\mathcal{F}$ is a space operator on $H$. Applying the inner product with the POD modes leads to

(9)

If $\mathcal{F}$ is only a linear space operator $\mathcal{F}(\boldsymbol {q})=\mathcal{L}(\boldsymbol {q})$, then

(10)

and (9) will change into

(11)

If $\mathcal{F}$ also contains a bilinear space operator $\mathcal{F}(\boldsymbol {q})=\mathcal{L}(\boldsymbol {q}) +\mathcal{Q}(\boldsymbol {q}, \boldsymbol {q})$, then

(12)

and (9) will change into

(13)

with

In one word, the Galerkin projection transforms the partial differential equation system into a system of ordinary differential equations directly, which can be solved by the explicit Runge-Kutta method and only requires very limited computational cost if the dimension $M$ is relatively small.

2.3 Isentropic ROM

In two-dimensional cavity flows, the fully compressible N-S equations for primitive variables $(ρ, u, v, p)$ or $(ρ, u, v, T)$ lead to a more complicated system. Rowley et al.^[16] suggested to use a further simplification valid for cold flows at low or moderate Mach numbers. The viscous dissipation and heat conduction in the energy equation can be neglected through the following assumptions: the dynamic and kinematic viscosities are constant; the temperature gradients are small; and the density gradients remain small and are dominated by the change of pressure. Then, the flow can be treated as isentropic, and only one thermodynamic variable (the sound speed $c^2=\gamma\mathcal{R}T$) is required, where $\gamma=1.4$ is the ratio of specific heats, $\mathcal{R}$ is the gas constant, and $T$ is the temperature.

Scaling the streamwise and crosswise components of velocities $u$, $v$ by the free stream velocity $U_{\infty}$, the local sound speed $c$ by the ambient sound speed $c_{\infty}$, the streamwise length $x$ by the cavity length $L$, the crosswise length $y$ by the initial momentum thickness $θ_0$ at the cavity leading edge, and time by $L/U_{\infty}$, the non-dimensional isentropic compressible N-S equations are given by

(14)

(15)

(16)

where ${Re}_L=U_{\infty}L/\nu$ is the Reynolds number based on the free stream velocity and cavity length, $\nu$ is the kinematic viscosity, and ${M\!a}=U_{\infty}/c_{\infty}$ is the Mach number.

$\boldsymbol {q}=(u, v, c)^\mathrm T$ is the vector of flow variables, and the system is rewritten as

(17)

with

Choosing the local sound speed $c$ as the thermodynamic variable makes it possible to introduce an energy-based inner product that includes both kinematic variables $(u, v)$ and thermodynamic variable $(c)$ in the subspace,

(18)

where $\alpha=1$ denotes using the stagnation enthalpy as the norm, and $\alpha=1/\gamma$ denotes using the stagnation energy as the norm. In this paper, we use this model for comparison and choose $\alpha=1$.

2.4 Approximate full ROM

The common POD for compressible flows is based on the weighting inner product of the primitive variables $\left (ρ, u, v, T\right)$, but the ROM constructed by this POD has not yet been reported. It is because that the fully compressible N-S equation for these primitive variables is a cubical implicit system. However, we can simplify it based on two assumptions: (i) the dynamic viscous coefficient $\mu$ is a constant; (ii) the density in the momentum and energy equations equals the free stream density after introduction of some artificial dissipation. In this way, the viscous dissipation and heat conduction in the energy equation are still maintained, in comparison with the isentropic ROM. Scaling the density $ρ$ by the free stream density $ρ_{\infty}$, the streamwise and crosswise components of velocities $u$, $v$ by the free stream velocity $U_{\infty}$, the temperature $T$ by ($\gamma-1) T_{\infty}$ with the free stream temperature equal to $T_{\infty}$, the streamwise length $x$ by the cavity length $L$, the crosswise length $y$ by the initial momentum thickness $θ_0$ at the cavity leading edge, and time by $L/U_{\infty}$, the non-dimensional governing equations are given by

(19)

(20)

(21)

(22)

where ${Pr}=0.7$ is the Prandtl number, ${Re}_L=ρ_{\infty}U_{\infty}L/\mu$ is the Reynolds number based on the free stream velocity and cavity length, and ${M\!a}=U_{\infty}/c_{\infty}$ is the Mach number. The ambient sound speed $c_{\infty}$ gets from $c^2_{\infty} = \gamma \mathcal{R} T_{\infty}$. We scale the pressure by $ρ c^2_{\infty}$. Then, the non-dimensional pressure is obtained by $p=(\gamma-1)ρ T/ \gamma$.

The vector of flow variables becomes $\boldsymbol {q}=(ρ, u, v, T)^\mathrm T$, and then the system is expressed as

(23)

with

For the primitive variables $(ρ, u, v, T)$, the energy-based inner product cannot be obtained. According to the work of Lumley and Poje^[24], the weighting inner product is defined by

(24)

where $\alpha$ is a weighting number, and the value of $\alpha = 0.5$ is chosen in this work.

3 Results and discussion 3.1 DNS results

The cavity in the DNS is shown in Fig. 1, which has a length-to-depth ratio of $L/D=4$, the free stream Mach number is $M\!a=1.8$, and the initial incoming momentum thickness gives $L/{θ_0} \approx 68.2$. The Reynolds number based on the cavity length and the free stream velocity is $Re_L=18000$. The initial condition is a polynomial expression of the laminar Blasius boundary layer profile. A monitor is located in the shear layer near the cavity trailing edge $(x/L = 0.99$ and $ y/θ_0 = 0)$.

The DNS is based on the OpenFOAM C$^{++}$ Library, with spatial discretization by central upwind schemes of Kurganov and Tadmor^[25] and time integration by a second-order implicit backward approach. A small time-step of 4.0$\times{10^{-8}}$ s is employed, and the flow fields are recorded every 250 iterations. A non-uniform structured grid is used with the mesh points clustered near walls of the cavity. The length from the inlet to the cavity leading edge is equal to $6.82D$, and $12D$ is extended from the cavity trailing edge to the outlet. The upper cavity region is built up by $656\times150$ rectangular cells, while $200\times152$ rectangular cells are used inside the cavity region. The minimum cell length equals 1.0$\times{10^{-4}}$ m. Details about validation of the numerical methods and grid convergence study are available in the previous work^[9].

The time trace of the streamwise velocity at the monitor in the DNS is shown in Fig. 2. One can see that the self-sustained oscillations are periodic. The spectrum in Fig. 3 shows that the oscillations are dominated by a single frequency. The value of the Strouhal number $St=fL/U_{\infty}=0.464$ is in reasonable agreement with the value of 0.492 determined by the Rossiter formula^[26] $St=(n-0.25)/({M\!a}+1/0.57)$ for $n=2$, which means that the supersonic cavity flow oscillates in the single Rossiter Ⅱ mode. The flow structure can be seen in the previous work^[9], by the dynamic mode decomposition (DMD) for the B2 case.

3.2 Application of reduced-order models

When the self-sustained oscillations are well established (after almost 9000000 iterations of the DNS), the sub-zones $(x/L, y/θ_0)\in[-1, 3]\times [-15, 50]$ of 200 recorded flow fields corresponding to nearly three periods of oscillations are chosen as snapshots. Following the numerical methods given above, the isentropic ROM and the approximate full ROM can be obtained.

The distribution of the eigenvalues $\lambda_k/(\sum\limits_k\lambda_k)$ as a function of the mode number $k$ for both ROMs is shown in Fig. 4. One can see that most of the low-order modes for the approximate full ROM can capture more fluctuation energy than the isentropic ROM, and the high-order modes of the isentropic ROM contain more fluctuation energy than the other one. Hence, it is possible that one can use less POD modes in the approximate full ROM to model the supersonic cavity flow dynamics, which is verified below by numerical tests.

Figure 5 shows the first four POD modes of the streamwise velocity and percent energy captured for the isentropic ROM, and the corresponding results for the approximate full ROM are shown in Fig. 6. One can see that the flow structures of the first two POD modes for both ROMs are almost the same, and the difference exists in the percent of captured energy. The first POD mode of the approximate full ROM contains $58.20\%$ of the fluctuation energy, which is $5.12\%$ larger than the isentropic ROM. This mode mainly describes the distribution of the streamwise velocity in the shear layer. The second POD mode of the approximate full ROM contains $36.17\%$ of the fluctuation energy, which is $5.63\%$ smaller than the isentropic ROM. This mode mainly describes the distribution of the streamwise velocity in the cavity. The third POD mode of the isentropic ROM is similar to the fourth POD mode of the approximate full ROM, which focuses on the flow structures in the cavity. Meanwhile, the fourth POD mode of the isentropic ROM is similar to the third POD mode of the approximate full ROM, with flow structures mainly describing the dynamics of the shear layer. Overall, compared with the isentropic ROM, the new model captures more fluctuation energy linked to the dynamics of the shear layer.

The reconstruction of the streamwise velocity at the monitor with both ROMs is shown in Figs. 7 and 8, and the influence of the number of POD modes retained for the projection is included. At the initial stage, the amplitude of the model is in good agreement with that of the DNS. After a long time integration, the amplitudes increase to fixed larger values, which may be due to a lack of dissipation, to the differences between the original and projected equations, or to the rounding errors. For the isentropic ROM, this deviation is still large when 14 POD modes are retained for reconstruction, but there is no further modification by using 16, 18, or 20 POD modes. For the approximate full ROM, the deviation is large when 10 POD modes are retained for reconstruction, but there is no further modification by using 12, 14, or 16 POD modes. Hence, the isentropic ROM needs 16 POD modes in the Galerkin projection, while the approximate full ROM only needs 12 POD modes. Time trace for the coefficient of the first POD mode is shown in Fig. 9. The Runge-Kutta integration for both ROMs can stably proceed, although the amplitudes are overestimated slightly finally.

To check the ability of the ROMs in prediction of the flow dynamics, we compare the spectra of the streamwise velocity at the monitor based on the ROMs with the DNS results, which are shown in Fig. 10. The dominant frequencies and amplitudes given by two ROMs are almost the same, which in general agree with those of the DNS. The main differences between the two ROMs and the DNS results appear in high frequency components, in which the two ROMs have relatively larger amplitudes. We also compare the transient fields of the streamwise velocity for two models with the DNS results, which are shown in Fig. 11. The flow structures are quite similar, indicating that the ROMs can predict the flow dynamics reasonably.

We also apply both ROMs to supersonic cavity flows with $L/{θ_0} \approx 113.67$. The DNS results in the previous work^[9] show that the Rossiter Ⅱ mode is the dominant flow mode, and the Rossiter Ⅲ mode and the low frequency mode also appear, which means that the flow is more complex. It has been suggested that a $99\%$ of the energy used as cutoff can represent the flow field accurately^[22]. When 60 POD modes are used to construct both ROMs, $99.94\%$ of the fluctuation energy can be captured. Using the fourth-order Runge-Kutta method, the coefficients of the POD modes are integrated for the two ROMs, and the results are shown in Fig. 12. As we can see, the coefficients given by the isentropic ROM diverge finally, while the coefficients given by the present ROM converge to a quasi-periodic state. Apparently, the performance of the present ROM is much more stable than the isentropic ROM.

4 Conclusions

For incompressible flows, the method of POD only needs to consider the kinematic variables, and the energy inner product is simply used. In compressible flows, the POD is a little bit more complex because one also needs to consider the thermodynamic variables. Previously, the isentropic assumption is widely used to obtain an energy inner product for compressible flows^[16], which is just valid for flows at low or moderate Mach numbers. A more suitable weighting inner product for compressible flows in a wide range of Mach numbers was given by Lumley and Poje^[24]. Considering this point, a new POD/Galerkin reduced-order model based on this inner product about the primitive variables $(ρ, u, v, T)$ for fully compressible N-S equations is built, which is theoretically valid for supersonic flows. The results show that the presented ROM is successfully applied in supersonic cavity flows, and the dominant flow dynamics can be captured accurately. For instance, the dominant frequency and amplitude for the Rossiter Ⅱ mode given by the spectra of the proposed model are in agreement with the DNS results performed based on the OpenFOAM. Comparison between the transient flow fields and that of the DNS demonstrates that the proposed ROM can reproduce the flow field reasonably. Furthermore, the performance of the present ROM is shown to be more stable than that of the isentropic ROM. It is worth mentioning that the model reduction method can also be easily extended to other supersonic flows.

Acknowledgements In the DNS, the authors especially thank the OpenFOAM Foundation for the free, open source computational fluid dynamics software package OpenFOAM. The author Zhenhua WAN especially thanks the Fundamental Research Funds for the Central Universities.