- Email: [email protected]
Dynamic mode decomposition of turbulent cavity flows for self-sustained oscillations
International Journal of Heat and Fluid Flow 32 (2011) 1098–1110
Contents lists available at SciVerse ScienceDirect
International Journal of Heat and Fluid Flow journal homepage: www.elsevier.com/locate/ijhff
Dynamic mode decomposition of turbulent cavity flows for self-sustained oscillations Abu Seena, Hyung Jin Sung ⇑ Department of Mechanical Engineering, KAIST, 291 Daehak-ro, Yuseong-gu, Daejeon 305-701, Republic of Korea
a r t i c l e
i n f o
Article history: Received 31 May 2011 Received in revised form 16 August 2011 Accepted 23 September 2011 Available online 15 October 2011 Keywords: Self-sustained oscillations Dynamic mode decomposition Coherent structure
a b s t r a c t Self-sustained oscillations in a cavity arise due to the unsteady separation of boundary layers at the leading edge. The dynamic mode decomposition method was employed to analyze the self-sustained oscillations. Two cavity flow data sets, with or without self-sustained oscillations and possessing thin or thick incoming boundary layers (ReD = 12,000 and 3000), were analyzed. The ratios between the cavity depth and the momentum thickness (D/h) were 40 and 4.5, respectively, and the cavity aspect ratio was L/D = 2. The dynamic modes extracted from the thick boundary layer indicated that the upcoming boundary layer structures and the shear layer structures along the cavity lip line coexisted with coincident frequency space but with different wavenumber space, whereas structures with a thin boundary layer showed complete coherence among the modes to produce self-sustained oscillations. This result suggests that the hydrodynamic resonances that gave rise to the self-sustained oscillations occurred if the upcoming boundary layer structures and the shear layer structures coincided, not only in frequencies, but also in wavenumbers. The influences of the cavity dimensions and incoming momentum thickness on the self-sustained oscillations were examined. Ó 2011 Elsevier Inc. All rights reserved.
1. Introduction Turbulent flows characterized by a wide range of scales potentially include an infinite number of degrees of freedom. However, in many situations, the flow complexity may be reduced to a set of coherent features with a few characteristic structures possessing self-excited global modes in space and time. For the past several decades, Proper Orthogonal Decomposition (POD) has been extensively used to identify the coherent structures in turbulent flows (Lumley, 1967). POD extracts a complete orthogonal set of spatial eigenfunctions corresponding to the modes of the measured second-order correlation function. POD has proven to be an effective method for identifying dominant features and events in turbulent flows based on energy rank. POD modes are computed from the two point correlation tensor which is an average quantity and results in loss of some dynamic information in the data. Dynamic information is significant for flows in which a local absolute instability prevails within a finite region. Such systems can exhibit selfsustained resonant modes at specific complex frequencies. A global stability analysis yields a large stability matrix and may be computationally expensive. The Arnoldi method requires explicit knowl-
⇑ Corresponding author. Tel.: +82 42 350 3027; fax: +82 42 350 5027. E-mail address: [email protected] (H.J. Sung). 0142-727X/$ - see front matter Ó 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.ijheatfluidflow.2011.09.008
edge of the evolution operator describing the data, which makes this method unattractive. Schmid (2010) introduced a Dynamic Mode Decomposition (DMD) approach to extract dynamic mode information from a flow field based on the Koopman analysis of dynamical system (Rowley et al., 2009). The extracted dynamic modes, which may be interpreted as a generalization of the global stability modes, may be used to describe the underlying physical mechanism that is captured within the data sequence or to project large-scale problems onto a dynamical system with a significantly large number of degrees of freedom. The decomposition algorithm relies on the extraction of a low-dimensional subspace after fitting a high-degree polynomial to the data sequence without knowledge of the evolution operator from which the data was generated. The eigenvalues and eigenvectors of the low-dimensional subspace capture the principal dynamics of the flow. Although DMD assumes a linear evolution operator, even in nonlinear flows, it extracts the structures of a linear tangent approximation to the underlying flow and describes its dynamic behavior. The mathematics underlying this decomposition is related to the Koopman operator (Lasota and Mackey, 1994; Mezic, 2005; Rowley et al., 2009), where the Koopman modes reduce to the linear global eigenmodes for the linear process while related to the Fourier modes for the nonlinear periodic problems. The Koopman operator is related to other decompositions, as described by Schmid (2010). In the present
A. Seena, H.J. Sung / International Journal of Heat and Fluid Flow 32 (2011) 1098–1110
work, the DMD algorithm is applied to turbulent flows over an open cavity. Self-sustained oscillations have been recognized as the most important characteristic behavior of separated shear layers over open cavities since Rossiter (1964) investigated the influence of cavity dimension and free stream velocity on the resonant tones. To understand the mechanism underlying self-sustained oscillations and to elucidate the large-scale vortical structures responsible for such oscillations, numerous experimental and numerical studies have been carried out on incompressible flows as well as on compressible flows. As reported in previous studies of compressible flows (Rowley et al., 2002), the mechanism underlying self-sustained oscillations has been well-characterized by extracting large-scale vortical structures over an open cavity. In incompressible turbulent flows, however, it has not been easy to identify the coherent structures responsible for oscillations due to both incompressibility and turbulence, although self-sustained oscillations are predicted. Daoud et al. (2006) found that even if there are no resonant tones in the surface pressure spectra, the acoustic feedback can exist in cavity flows. Pereira and Sousa (1995) observed periodically oscillating shear layers in the flow of a turbulent incoming boundary layer over an open cavity. In PIV studies of turbulent open cavity flows, Kang et al. (2008) observed that the incoming turbulent boundary layer gave rise to large-scale vortical structures that were responsible for self-sustained oscillations. Instantaneous fields of velocity fluctuations were used in Ashcroft and Zhang (2005) mainly to identify largescale vortical structures in a separated shear layer, although the organized nature of the coherent structures was described in the statistical analysis of the correlations. The work of Lee et al. (2008), which was the large eddy simulation (LES) of high-Reynolds number incompressible turbulent flows over an open cavity, provided instantaneous isopleths and spanwise-averaged contours for pressure fluctuations in an effort to describe the large-scale vortical structures that were responsible for the peak frequency in the energy spectra. An analysis of raw data, including instantaneous velocity or vorticity distributions, usually requires that structural analysis of separated shear layers is instantaneous and qualitative. Therefore, application of an analysis approach such as DMD is required for a quantitative and dynamic characterization of large-scale vortical structures. The objectives of the present study were to (1) identify the large-scale vortical structures responsible for hydrodynamic oscillations by applying DMD to the pressure fluctuations of incompressible turbulent flows over an open cavity; (2) obtain dynamic information about the extracted structures; (3) elucidate the quantitative characteristics of the large-scale vortical structures by comparing the oscillating behavior of a separated shear layer to the behavior of a nonoscillating system. The sequences of images over time were processed to delve into the temporal dynamics of structures in turbulent open cavity flows. To achieve these objectives, numerical data from two different flows over an open cavity, with thin or thick incoming boundary layers, were analyzed by DMD at Reynolds number based on cavity depth ReD = 12,000 and 3000, respectively. The self-resonant oscillating modes were observed for the case of a thin incoming boundary layer. Analysis of the thin boundary layer revealed a dominant peak in the DMD spectrum, whereas all dynamic modes extracted from the thick boundary layer showed the gradual decay in energy with frequency representing viscous boundary layer structures. These results permitted examination of the effects of the upcoming momentum thickness on the self-sustained oscillations in an open cavity flow. Analysis was performed over the full computational domain as well as over a subdomain. The selection of the domain significantly influenced the results, depending on the flow conditions. Furthermore, DMD was applied to velocity flow distributions
1099
in which oscillating behavior was observed in a separated shear layer over the open cavity. The vortical structures responsible for self-sustained oscillations were thereby extracted. 2. Dynamic mode decomposition Dynamic mode decomposition is a recent extension of the classical Arnoldi technique. A temporal sequence of N data fields, consisting of column vectors, vj, that are equispaced in time according to the Nyquist criteria (Schmid, 2010), can be written as
VN1 ¼ fv 1 ; v 2 ; v 3 ; . . . ; v N g:
ð1Þ
The basic premise of the method is that each snapshot in time is assumed to be generated by a linear dynamical system,
vjþ1 ¼ Avj :
ð2Þ
The eigenvalues and eigenvectors of the matrix A completely characterize the behavior of the dynamical system. DMD is a method for computing the approximate eigenvectors or Ritz vectors of a system matrix. A high-degree polynomial is fit to a Krylov sequence of flow fields. As the number of snapshots increases, the flow is assumed to approach a linear dependency after a sufficient number of snapshots. The last image is the linear combination of the previous images and may be described by
vN ¼ c1 v 1 þ c2 v 2 þ . . . þ cN1 v N1 þ r;
ð3Þ
where r is the residual vector and c = {c1, c2, c3, . . ., cN1} are the coefficients,
vN ¼ VN1 c þ r: 1
ð4Þ
Eq. (4) represents an over-determined system of equations. The coefficients may be obtained using the least squares method. The number of required snapshots, N, may increase until the residual vector r converges. Following Ruhe (1984), the least squares description of the full system matrix A may be written as
Afv 1 ; v 2 ; v 3 ; . . . ; v N1 g ¼ fv 2 ; v 3 ; v 4 ; . . . ; v N g cg þ reTN1 : ¼ fv 2 ; v 3 ; v 4 ; . . . ; VN1 1
ð5Þ
Furthermore, Eq. (5) may be written as
AVN1 1
¼ VN2 ¼ VN1 C þ reTN1 : 1
ð6Þ
The matrix C is of the companion type,
0
0
c1
1
C B c2 C B1 0 C B C¼B : : : C C; B C B 1 0 c @ N2 A 1 cN1
ð7Þ
which simply shifts the snapshot (Schmid, 2009) index from 1 to N 1. This matrix, extracted from the data sequence, represents a low-dimensional system matrix representation of the full system matrix. The characteristic solution to the matrix C approximates some of the eigenvalues (kj) of the full system matrix A, called as empirical Ritz values (Rowley et al., 2009). The logarithmic mapping of these eigenvalues provides the growth rate (rj = Re{kj}/Dt) and frequency (xj = Im{kj}/Dt). The empirical Ritz values lying on the unit circle represent the modes with zero growth rates, whereas the eigenvalues lying inside and outside the unit circle represent the damped and undamped modes, respectively. If the data sequence is extracted from a nonlinear process, then neutrally stable eigenvalues are expected as the data sequence becomes sufficiently long. This is due to the fact that instabilities have nonlinearly saturated while decaying processes have vanished from the data sequence; only neutrally stable fluid elements remain. The
1100
A. Seena, H.J. Sung / International Journal of Heat and Fluid Flow 32 (2011) 1098–1110
eigenvectors provide a linear combination of basis coefficients that may be used to extract the corresponding dynamic modes as,
Uj ¼
VN1 Tj : 1
ð8Þ
In contrast to the POD where a spatial orthogonality of the identified structures is enforced at the expense of a mixture of frequencies in each individual POD mode, the basis vector Tj obtained above aims at an orthogonality in time by identifying pure frequencies. There is no normalization step in the algorithm, and hence the modes come with the amplitudes, ||U|| that provide a way to rank their contribution to the overall energy in the set VN1 (Bagheri, 2010; Rowley et al., 2009). The companion matrix C may be a highly non-normal matrix, which yields an ill-conditioned eigenvalue decomposition problem. To improve the accuracy of the DMD method, Bagheri (2010) suggested that the matrix should first be balanced by a similarity transformation, followed by reducing to an upper Hessenberg form via a second similarity transformation, and finally, the eigenvalues should be computed using the QR algorithm. 3. Numerical simulation For an incompressible flow, the non-dimensional governing equations are
@ui @ui uj @p 1 @ 2 ui ¼ iþ ; þ @t @xj @xi ReD @x2j @ui ¼ 0; @xi
LY
i ¼ 1; 2; 3
ð9Þ
ð10Þ
where xi are the Cartesian coordinates and ui are the corresponding velocity components. The free-stream velocity U1 and the cavity depth D were used to non-dimensionalize the equations. The Reynolds number was defined as ReD = U1D/m, where m is the kinematic viscosity. The governing Eqs. (9) and (10) were integrated in time using the fully implicit decoupling method proposed by Kim et al. (2002). A schematic diagram of the computational domain is shown in Fig. 1. All terms in Eq. (9) were advanced in time using the Crank– Nicholson method and were resolved using the second-order central difference scheme in space. A block LU decomposition was used to achieve velocity–pressure decoupling as well as additional decoupling of the intermediate velocity components in conjunction with an approximate factorization. The multi-grid method was applied to obtain a solution to the pressure Poisson equation over a Tshaped domain. In the present simulations, a turbulent boundary layer with realistic velocity fluctuations, which were generated using the method of Lund et al. (1998), was provided at the inlet. A direct numerical simulation (DNS) of incompressible flows over an open cavity was performed at ReD = 3000, where the DNS data were provided at the inlet with Reh = 300, 670 and 1410. The cavity flows at high Reynolds number (ReD = 12,000) were simulated using a LES with a dynamic subgrid-scale model (Germano et al., 1991). Here, DNS data at Reh = 300 were provided at the inlet. To analyze the influence of upstream turbulence on cavity flow, four cases with D/h = 2.1, 4.5, 10 and 40 were simulated. The simulation parameters are listed in Table 1, where O and X represent the flow with and without self-sustained oscillations, respectively. Among them, we selected two typical cases with D/h = 40 and 4.5 and applied DMD to understand hydrodynamic resonance and quantitative characteristics of the large-scale vortical structures. To ascertain the reliability of the incoming turbulent boundary layer, the turbulence statistics were compared with the DNS data of Spalart (1988), and the results are presented in Fig. 2. Fig. 2a shows the profile of the mean velocity normalized by the friction velocity (u+ = u/us), with variations in y+ = yus/v, and Fig. 2b shows
LZ L
Lo D
y x Li
z
Fig. 1. Schematic diagram of the computational domain.
a comparison of the turbulence intensities. The results were in excellent agreement with those of the previous study (Spalart, 1988). Uniform grids were used in the spanwise direction. Nonuniform grids were used in the streamwise and wall-normal directions, with higher densities of grid points near the leading and þ trailing edges. The sizes of Dxþ min and Dymin were determined by obtaining convergence of the turbulent intensities near the wall and the trailing edge. The computational time steps were freed from the Courant Friedrich Levy limit by employing the fully implicit decoupling method. The time step was fixed at Dt = 0.007D/U1 for ReD = 3000. For ReD = 12,000, the time step was very small, Dt = 0.001875D/U1, due to the explicit characteristics of the dynamic subgrid-scale model. The computational domain size was determined after checking the two-point correlations of the velocity fluctuations. The domain size in the spanwise direction (Lz) was four times the cavity depth (D) for both ReD = 3000 and 12,000. Selection of these spanwise domain sizes was supported by results from previous numerical simulations of turbulent flows over a backward-facing step (Le et al., 1997). The domain size in the streamwise direction was Li = Lo = 9D for ReD = 3000, and Li = 4D, Lo = 7.5D for ReD = 12,000. Here, Li and Lo are the streamwise lengths of the development section ahead of the leading edge and behind the trailing edge. The vertical height (LY) was set to exceed 9D for ReD = 3000, and the height was 5D for ReD = 12,000. The computational details are given in Lee et al. (2010). 4. Results and discussion Before the DMD algorithm was applied to turbulent cavity flows, it was validated using a fabricated pattern. Take, for example, the pattern q, such that
q ¼ q0 þ q1 þ q2 þ q3 ;
ð11aÞ
q0 ðx; y; tÞ ¼ expðy2 =0:7Þ; qn ðx; y; tÞ ¼ an ðtÞ
m¼1 X
ð11bÞ "
ð1Þm exp
m¼1
n ¼ 1; 2; 3; . . . ;
ðx bn m cn tÞ2 dn
!
#
þ
y2 ; dn ð11cÞ
where dn = anx + bn is the diameter of the structure, an and bn are constants, bn is the wavelength factor, cn is the velocity of the pattern in the streamwise direction, and an(t) is the growth or decay factor. The pattern adopted for validation was composed of patterns of three different sizes (n = 1, 2, and 3), each traveling with a fixed velocity. The constants and factors an(t) are given in Table 2. The fabricated pattern is shown in Fig. 3a, and consisted of the superposition of the three patterns with n = 1, 2, and 3. For n = 1, the factor
1101
A. Seena, H.J. Sung / International Journal of Heat and Fluid Flow 32 (2011) 1098–1110 Table 1 Simulation conditions. D/h
Reh
ReD
Nx Ny Nz
Dxþ ; Dxþ max min
Dyþ min
Dz +
Self-sustained oscillations
2.1 4.5 10 40
1410 670 300 300
3000 3000 3000 12,000
513 383 513 513 213 257 513 133 129 897 169 257
1.00, 1.09, 1.17, 1.40,
0.196 0.188 0.196 0.560
5 5 5 10
X X X O
16.2 17.7 19.0 40.0
Fig. 2. Mean streamwise velocity profiles and velocity fluctuations in the inflow data.
an(t) was set to unity such that the pattern shown in Fig. 3b neither grew nor decayed over time. The factors an(t), for n = 2 and 3, were selected such that the structures q2 decayed over time and q3 grew in time. Next, the fabricated pattern q was processed using the DMD algorithm, as described in Section 2. The extracted Ritz values kj of the low-dimensional system matrix C are shown in Fig. 4a. The size and color of the eigenvalues indicate the amplitude of the respective mode in the data sequence. The eigenvalues kj of C are transformed via the logarithmic mapping and shown in Fig. 4b, where the unstable eigenvalues appear in the right halfplane. The eigenvalues obtained were found to be in complex conjugate pairs. As a result, the DMD spectrum was symmetric about the frequency x = 0, as shown in Fig. 4b, where r represents the growth/decay rate and x is the frequency of each mode. The global energy norms of each mode are plotted as a function of frequency in Fig. 4c. Four energetic structures, in which x = 0 contained the highest energy, were extracted from the fabricated pattern. Four distinct modes are marked 1–4 in the extracted spectra shown in Fig. 4. Fig. 5 shows their respective modes. The eigenvalue that appeared at the origin, marked 1, accounted for the steady state contained in the pattern and is shown in Fig. 5b. The second DMD mode with the second highest energy is marked 2 in the spectra and lie on the neutral line (r = 0), i.e., neither growing nor decaying over time. The modes corresponding to the real and imaginary parts, respectively, are shown in Fig. 5c and d. They clearly represented the pattern q1 with a decay factor an(t) = 1. The third highest mode, marked 3, lay in the stable region (r < 0) and decayed over time. These modes are shown in Fig. 5e and f. The extracted Table 2 Constants for the fabricated pattern. n
an
bn
bn
cn
an(t)
1 2 3
0.030 0.010 0.005
0.005 0.100 0.100
0.8 0.4 0.3
0.377 0.252 0.126
1 et/80 0.1 (1 et/15) + 0.2
mode represented the pattern q2 with a negative decay factor (an < 0). Finally, the fourth mode, marked 4, was observed in the unstable region (r > 0), as shown in Fig. 5g, and represented q3 of the fabricated pattern. The above example illustrates how the DMD method decomposes the patterns in the data and reveals important dynamic information. Next, DMD was employed to describe the open cavity data. The flows corresponding to each of the two cases considered in the present study, ReD = 3000 and 12,000, were simulated. Fig. 6 shows the mean streamlines corresponding to the two simulations. Two large recirculation eddies are observed in both simulations. Mean cavity flow is characterized by the primary clockwise rotating vortex covering more than downstream half of the cavity which elongates to the leading edge. A weaker anticlockwise rotating secondary vortex is evident near the upstream corner, induced due to the primary vortex. After a long initial transient period, a sequence of images was saved. An initial transient time of more than 10 ‘‘flow-through’’ cycles was discarded to allow the passage of unphysical fluctuations. The temporal sequences of the system snapshots were processed using the DMD algorithm to extract pertinent dynamic characteristics of the flow. One hundred fifty instantaneous snapshots of the pressure fluctuations were used for the ReD = 3000 condition, and 124 snapshots were used for the ReD = 12,000 condition. The eigenvalues of the subspace C were expected to converge toward some of the eigenvalues of the full system matrix A. The DMD algorithm may be applied in one of two possible ways. It may be applied over the entire physical domain, including both the upstream region of the leading edge and the downstream region of the trailing edge, as shown in Fig. 7a, or it may be applied over a subdomain, as shown in Fig. 7b. Both approaches were considered in the present work, and the results were compared. Analysis of the subdomain region only, shown in Fig. 7b, focused on the vortical flow inside the cavity and the shear layer that formed and detached near the right wall of the cavity. Analysis of the full domain, as shown in Fig. 7a, captured the effects of the upstream boundary layer structures in addition to the cavity flow.
1102
A. Seena, H.J. Sung / International Journal of Heat and Fluid Flow 32 (2011) 1098–1110
Fig. 3. Fabricated pattern.
Fig. 4. (a) Ritz values kj, (b) DMD spectrum and (c) energy spectrum for the fabricated pattern.
Fig. 5. Modes extracted from the fabricated pattern by the DMD algorithm.
A. Seena, H.J. Sung / International Journal of Heat and Fluid Flow 32 (2011) 1098–1110
1103
Fig. 6. Mean flow streamlines.
Fig. 7. Physical domain for application of DMD: (a) full domain and (b) subdomain.
4.1. ReD = 3000 and D/h = 4.5 DNS of incompressible cavity flow over an open cavity were performed at ReD = 3000. A subdomain analysis was performed first. One hundred fifty images of the instantaneous pressure fluctuations at a sampling rate of Dt = 0.11725 were processed using the DMD algorithm. The eigenvalues of the subspace approximated the eigenvalues of the full system matrix and provided a quantitative illustration of the underlying temporal dynamics. The eigenvalues from the subdomain analysis are shown in Fig. 8a, where the symbol color indicates the global energy norm of the associated mode. Nearly all the Ritz values are on the unit circle ||kj|| = 1, indicating that the sample points vi lie on or near an attracting set. The norm of each mode indicates the energy in the corresponding mode. The energies of the extracted modes are shown as a function of frequency in Fig. 8b, where each mode is displayed with the vertical line scaled with its magnitude at its corresponding frequency. The eigenvalue marked 1that typically appears at the origin (x = 0) with the maximum energy accounts for the steady state and represents the mean flow as such depicted by a black symbol in Fig. 8a. The corresponding dynamic mode is shown in Fig. 9a, where the low-pressure region observed in the cavity indicates the primary vortex. A high-pressure region is also observed at the trailing edge,
possibly due to the impingement of shear layer vortical structures. The other modes show the spatial structures displaying periodic behaviors in time. A few specific modes in Fig. 8b dominate the energy spectrum at x = 1.9, 2.8 and 3.9 rad/s, marked 2–4. The modes are shown in Fig. 9b–d. The solid and dotted lines represent the negative and positive pressure fluctuations, respectively. In all of the modes shown in Fig. 9b–d, the pressure fluctuations are observed above the shear layer, suggesting that they correspond to the upcoming boundary layer structures flowing over the cavity and further impinging on the trailing edge. To gain further insight into these boundary layer structures, the DMD algorithm was applied over the full cavity domain; the obtained spectrum is shown in Fig. 10a. The symbol color indicates the norm of the corresponding mode. The black symbol at zero frequency, marked 1, represents the steady state and the corresponding mode is shown in Fig. 11a. Their norms against frequency are shown in Fig. 10b. Five distinct frequency peaks are marked 2–6, and their corresponding modes are shown in Fig. 11b–f. The eigenvalue marked 2 in the spectrum indicates a very low frequency mode. The pressure fluctuations of the corresponding mode, shown in Fig. 11b, indicate the presence of a very large-scale structure relative to the size of the cavity which convects over the cavity. The mode marked 3 is shown in Fig. 11c. In addition to the upcoming boundary layer structures, the structures generated due to shear layer oscillations could be observed along the cavity lip line, which further impinged on the trailing edge. The eigenvalue marked 4 is shown in Fig. 11d. The structures that appear in this particular mode suggest that at this frequency, the wavenumber corresponding to the upcoming boundary layer structure and the perturbation generated by shear layer oscillations are equal. The modes of two other eigenvalues marked 5 and 6 are displayed in Fig. 11e and f. Note that all the modes in Fig. 11 represent the incoming viscous boundary layer structures.
Fig. 8. Ritz values and their magnitudes at ReD = 3000 for the cavity.
1104
A. Seena, H.J. Sung / International Journal of Heat and Fluid Flow 32 (2011) 1098–1110
Fig. 9. DMD modes inside the cavity at ReD = 3000.
Fig. 10. Ritz values and their magnitudes at ReD = 3000 for the full domain.
However, as compared to subdomain analysis, the energy spectrum for the full domain analysis, shown in Fig. 10b, shows no leading modes, and the magnitude gradually decays as a function of frequency, suggesting the presence of a wide range of turbulent structures. This result is consistent with the absence of self-sustained oscillations. Analysis of the subdomain and full domain suggests that the selection of only a subdomain in the DMD analysis is inappropriate for understanding the dynamics pertinent to the flow. The DMD spectrum of the full domain reveals important characteristics at the boundary layer. The dominant features illustrate the boundary layer structures over the full domain at ReD = 3000. 4.2. ReD = 12,000 and D/h = 40 LES of incompressible cavity flow over an open cavity were performed at ReD = 12,000. Realistic velocity fluctuations of Reh = 300 were provided at the inlet. One hundred twenty-four images of the instantaneous pressures at a sampling rate of Dt = 0.075 were processed using the DMD algorithm, and a low-dimensional subspace was extracted. The analysis was performed for both cavity subdomain and full domain. The extracted spectrum, i.e., the spectrum of C, is displayed in Fig. 12a and b for subdomain and Fig. 12c and d for full domain. Fig. 12a and c shows that for both cases, all the Ritz
values are on the unit circle. The symbol color indicates the global energy norms of the modes. The amplitude of modes against frequency is shown in Fig. 12b and d. Both spectra, obtained from the subdomain and the full domain, displayed dominant peaks, marked 2 and 4 with two leading modes at x = 3.5 and 4.6 rad/s, respectively. Some of the selected modes from the spectrum in Fig. 12b, marked 1–4, are shown in Fig. 13. The eigenvalue at the origin, marked 1 in the spectrum (Fig. 12b), represents the steady state component of the pressure, and its respective mode is shown in Fig. 13a. The modes marked 2 and 3 with dominant peaks and their respective modes are shown in Fig. 13b and c, respectively. One pair of negative and positive distributions represents the low-pressure fluctuations of a large-scale vortical structure and the high-pressure fluctuations of induced rotational motions. Three pairs are clearly observed between the leading and trailing edges. This was consistent with the spectral characteristics of self-sustained oscillations corresponding to NR = 3, reported by Lee et al. (2008) The streamwise length scale of the coherent structures gradually increases from 0.3D in the region immediately downstream of the leading edge to 0.8D in the region that impinges on the trailing edge. To avoid ambiguity, the length scale of the coherent structure was determined to be twice as long as the streamwise distance between the central locations of the positive and negative distributions. The transverse length scale was not sufficiently large to affect the pressure
A. Seena, H.J. Sung / International Journal of Heat and Fluid Flow 32 (2011) 1098–1110
1105
Fig. 11. DMD modes inside the cavity full domain at ReD = 3000.
Fig. 12. Ritz values and their magnitudes at ReD = 12,000 for (a) the cavity and (b) the full domain.
fluctuations on the bottom wall inside the cavity. These modes suggest the presence of self-sustained oscillations in the cavity. The symbol marked 4 represents next leading mode, and its respective mode is shown in Fig. 13d.
We performed POD to the same data set for comparison. Fig. 14 shows the first four POD eigenmodes and the solid and dotted lines represent the positive and negative distributions of eigenmodes. In the first and second modes, three pairs of well-organized
1106
A. Seena, H.J. Sung / International Journal of Heat and Fluid Flow 32 (2011) 1098–1110
Fig. 13. DMD modes inside the cavity at ReD = 12,000.
structures are observed similar to the DMD modes obtained in Fig. 13. These first two modes look like two different phase states of an identical motion as expected for the flow in experiencing global mean advection. The most energetic POD mode in Fig. 14a captured a similar spatial structure to the most energetic DMD mode shown in Fig. 13b, but the difference is shown in the time coefficients. Variations of the mode coefficient for the most energetic POD mode (Fig. 14a) and DMD mode (Fig. 13b) are shown in Fig. 15a, which are the coefficients used when reconstructing the flow field with isolated modes. Fig. 15a shows that a single DMD mode contains by construction only a single frequency component, whereas the POD modes capture the most energetic structures, resulting in the mode that contain several frequencies. This is apparent in the power spectrum of time-varying coefficients shown in Fig. 15b, where several peaks are observed in the case
of most energetic POD mode, whereas a single peak is observed in the case of most energetic DMD mode. Both methods show the dominant peak at f = 0.765 Hz, hence validating that the mode extracted in Fig. 13b is responsible for the self-sustained oscillations. To determine whether these pressure fluctuation modes arise from shear layer oscillations or from upcoming boundary layer structures, DMD analysis was performed over the full cavity domain. The spectrum for the full domain is shown in Fig. 12d. The spectrum obtained from the full domain analysis resembles the spectrum from the subdomain analysis Fig. 12b. The modes corresponding to the dominant peaks, marked 1–4 in Fig. 12d, are shown in Fig. 16. These modes resemble the modes obtained from the subdomain analysis. The size of structures along the shear layer is comparable to that of the upstream boundary layer structures.
Fig. 14. POD modes inside the cavity at ReD = 12,000.
Fig. 15. (a) Time coefficients of the DMD and POD modes and (b) power spectra of a(t).
A. Seena, H.J. Sung / International Journal of Heat and Fluid Flow 32 (2011) 1098–1110
The modes in Fig. 16b–d indicate that the structures along the shear layer coincide with the wavenumber spectra of the upstream boundary layer structures. This coincidence may have produced self-resonant modes, which will be discussed later. Although alternating pressure fluctuation patterns were regularly observed in the first and second modes, the distributions in the other modes were irregular. The other modes most likely describe the fluctuating behavior of a separated shear layer with a high wavenumber, albeit with slight variations. In previous studies of laminar cavity flows, the oscillating behavior of a separated shear layer was represented only by the first two modes, whereas the next higher modes were associated with recirculation of the primary vortex inside the cavity (Podvin et al., 2006). The simple modes of a separated shear layer produce a single peak frequency in the power spectra of laminar cavity flows. In the present study of turbulent cavity flows, however, the fluctuating behavior of the separated shear layer was represented in the higher modes as well as in the first two modes. The complex mixture of several modes produces a broad spectrum, corresponding to the pressure fluctuations in the fluctuating behavior of the separated shear layer, although the first two modes are responsible for the peak frequency corresponding to NR = 3. Fig. 17 shows the time history of the dynamic mode responsible for self-sustained oscillations. The time period T is the convection time for the vortical structure to convect along the shear layer from the leading edge to the trailing edge. The sequence of contour plots in Fig. 17a–d shows the convective patterns of large-scale vortical structures. The negative pressure fluctuations represent a clockwise rotating large-scale vortical structure (Lee et al., 2010). At time t = to, marked V1 indicates a large-scale vortical structure at the upstream of the leading edge and marked V2 represents another vortical structure located at (x/D, y/D) = (0.5, 1.0). After separation at the leading edge, the vortical structure attains energy due
1107
to the shear layer oscillations and enters the cavity shear layer. As time progresses to t = to + T/2, these structures convect along the shear layer. For t = to + 3T/4, the fringe of marked V2 begins to impinge on the trailing edge. The contour plot shows the partial clipping of marked V2 and continuous convection of marked V1. A small part of the clipped marked V2 is entrained into the cavity along the vertical wall whereas the majority is ejected out of the cavity. After escaping from the cavity, marked V2 gradually dissipates along the downstream wall. The coherent structure responsible for self-sustained oscillations was investigated through DMD analysis of the velocity data. A velocity field corresponding to each DMD mode of the pressure fluctuations was obtained. The coherent structures extracted for the leading modes at x = 4.6 and 3.5 rad/s are shown in Fig. 18a and b, respectively. This is consistent with the velocity field corresponding to the most energetic POD mode in the work of Lee et al. (2008) for the same data set. 4.3. Self-sustained oscillation modes The disturbance in the shear layer across the cavity generates coherent structures that impinge on the trailing edge. The system exhibit self-excited global modes at specific complex frequencies that are globally unstable. A key question is whether local disturbances due to incident turbulence can prevent the occurrence of global instabilities that originate from the singularity of the cavity geometry. DMD analysis suggested that the upcoming boundary layer thickness was significantly responsible for the self-sustained oscillations, and, therefore, the parameter D/h was important in a description of open cavity flows. The unsteady modes were identified at ReD = 12,000 for a thin incoming boundary layer. The dynamic modes extracted from the thick boundary layer (ReD = 3000) showed that both the boundary layer structures and
Fig. 16. DMD modes inside the cavity full domain at ReD = 12,000.
Fig. 17. Dynamic patterns of the mode marked 2R in the spectrum.
1108
A. Seena, H.J. Sung / International Journal of Heat and Fluid Flow 32 (2011) 1098–1110
Fig. 18. Streamlines showing the selected dynamic modes at ReD = 12,000 for the cavity full domain.
the shear layer structures coexist with similar frequency space but with different wavenumber space. In contrast, structures with a thin boundary layer (ReD = 12,000) showed complete coherence among the modes, which produced self-sustained oscillations. These results suggest that the hydrodynamic resonance that gives rise to self-sustained oscillations occurs when the upcoming boundary layer structures and the shear layer structures assume identical frequency space and wavenumber space. Wavenumber coincidence relies on the size of the largest structure in the upcoming boundary layer which depends on Reh, such that its size is comparable to the shear layer dimensions. Resonance may occur if the size of the largest structure of the upstream boundary layer is comparable to the wavelength of the shear layer oscillation. The cavity dimensions (length and depth) play a significant role in the shear layer dimensions. Resonance may occur when the shear layer dimensions are equal to the length scales of the two largest structures of the upcoming boundary layer. The largest structure in the upcoming boundary layer is comparable to the boundary layer thickness. In light of the above analysis, the minimum conditional requirement for restraining self-resonant modes relies on the cavity dimension. The cavity length should be less than twice the wavelength of the largest and most energetic structure in the upcoming boundary layer. The cavity depth should be larger in size than the upcoming largest boundary layer structure. This is consistent with the finding of Sarohia (1977) that no self-sustained oscillations occur below a certain critical cavity length and depth. The number of vortical structures can be evaluated using the modified form of Rossiter’s empirical relation,
L fL ¼ ¼ NR ; kx U C
ð12Þ
where kx is the streamwise wavelength of oscillation, f is the oscillating frequency, Uc is the convection velocity, and NR is the number of wavelengths of the fundamental frequency that may be contained within the cavity length in the NRth oscillation mode (Gharib and Roshko, 1987). Eq. (12) is only valid in the limit of incompressible flow because it was derived under the assumption that the lag time is zero between the impact of a hydrodynamic structure on the trailing edge and the emission of an acoustic wave. Eq. (12) also supports the present suggestion in the sense that a minimum critical cavity length of the order of wavelength of structures is required to generate self-sustained oscillations, or in other words the ratio of cavity length to wavelength should be of positive integer order describing its mode number. For each mode, the average wavelength of the structures was measured along the shear layer. To avoid ambiguity, the wavelength (kx) of the coherent structure was determined to be twice as long as the streamwise distance between the central locations of the positive and negative distributions. Fig. 19a and b plots the wavelength as a function of the frequency for both cases. The DMD algorithm extracted a wide range of structures. We found that the large structures were associated with lower frequencies, and the frequency increased as the wavelength of the structures decreased. This result agrees well with Eq. (12), which describes the wavelength as the reciprocal of the frequency of structures with the same convection velocity. Fig. 19c and d shows the frequency on a log–log scale. The slope 1 indicates that the
Fig. 19. Wavelength of the structures as a function of frequency at ReD = 3000 and 12,000.
A. Seena, H.J. Sung / International Journal of Heat and Fluid Flow 32 (2011) 1098–1110
1109
Fig. 20. Wavenumber spectra at ReD = 3000 and 12,000.
wavelength and frequency follow a reciprocal relationship. In both cases, ReD = 3000 and 12,000 (Fig 19c and d, respectively), the reciprocal trend suggests that structures convect downstream with the same velocity. All modes extracted using the DMD algorithm indicated that the flow was characterized by turbulent structures. The L2-norm of the modes yielded the energies in the associated modes. The energy (||U||) is plotted against wavenumber in Fig. 20a and b for both cases analyzed in the present work. The Kolmogorov power law (5/3) is also plotted (solid line) in the same figure. The high value of the power law prefactor at ReD = 12,000 shows that it contains higher energy modes than the cavity flow at ReD = 3000. The energies are plotted on a log–log scale in Fig. 20c and d. The region in which the data agree well with the Kolmogorov power law corresponds to the turbulence cascade region called the inertial range. The inertial range was broader for the higher Reynolds number flows, ReD = 12,000, than for the thick boundary layer flow (ReD = 3000). The modes extracted using the DMD algorithm showed a turbulent energy cascade, thus supporting the turbulent nature of the flow. 5. Conclusions In the present study, we identified large-scale vortical structures responsible for self-sustained oscillations by employing a dynamic mode decomposition of the pressure fluctuations of turbulent flows over an open cavity. DMD was applied to incompressible turbulent flows over an open cavity at ReD = 3000 and 12,000, with upstream turbulence of Reh = 670 and 300, respectively, provided at the inlet. The influence of the incoming turbulent boundary layer on the self-sustained oscillations was investigated. A wide range of structures was extracted using the DMD algorithm. The DMD algorithm identified the dynamic mode which represented the coherent structures responsible for the self-sustained oscillations. The cavity with thick incoming boundary layer represented the viscous boundary layer structures, whereas the thin incoming boundary layer showed a dominant
peak in the spectrum which represented the cavity dynamics. A mode with the dominant peak at ReD = 12,000 represented the structure responsible for self-sustained oscillations. After processing the instantaneous velocity distributions, the large-scale vortical structures responsible for self-sustained oscillations were identified. These structures were consistent with the spectral characteristics of self-sustained oscillations corresponding to NR = 3. The dynamic modes extracted from the thick boundary layer (ReD = 3000) showed that both the boundary layer structures and the shear layer structures coexist in similar frequency space but with different wavenumber space. These results suggest that the hydrodynamic resonance that produces self-sustained oscillations occurs when the upcoming boundary layer structures and the cavity shear layer structures coincide not only in frequency space, but also in wavenumber space. Hence, the upcoming boundary layer thickness was found to be a significant parameter for describing self-sustained oscillations. The parameter D/h, for flows over an open cavity, was, therefore, found to be important. The DMD energy spectra agreed well with the Kolmogorov power law, supporting the conclusion that the extracted modes represent the turbulence scales. Selection of the domain for analysis using the DMD algorithm significantly affected observation of the dynamics pertinent to the flow.
Acknowledgements This study was supported by the Creative Research Initiatives (No. 2011-0000423) program of the National Research Foundation of Korea.
References Ashcroft, C., Zhang, X., 2005. Vortical structures over rectangular cavities at low speed. Phys. Fluids 17, 015104. Bagheri, S., 2010. Analysis and control of transitional shear flows using global modes. PhD Thesis, February, Department of Mechanics, Royal Institute of Technology, Sweden.
1110
A. Seena, H.J. Sung / International Journal of Heat and Fluid Flow 32 (2011) 1098–1110
Daoud, M., Naguib, A.M., Bassioni, I., Abdelkhalek, M., Ghoneim, Z., 2006. Microphone-array measurements of the floor pressure in a low-speed cavity flow. AIAA J. 44 (9), 2018–2023. Germano, M., Piomelli, U., Moin, P., Cabot, W.H., 1991. A dynamic subrid-scale eddy viscosity model. Phys. Fluids A 3, 1760–1765. Gharib, M., Roshko, A., 1987. The effect of flow oscillations on cavity drag. J. Fluid Mech. 177, 501–503. Kang, W., Lee, S.B., Sung, H.J., 2008. Self-sustained oscillations of turbulent flows over an open cavity. Expt. Fluids 45 (4), 693–702. Kim, K., Baek, S.-J., Sung, H.J., 2002. An implicit velocity decoupling procedure for the incompressible Navier-Stokes equations. Intl. J. Numer. Methods Fluids. 38 (2), 125–138. Lasota, A., Mackey, M.C., 1994. Chaos, Fractals and Noise: Stochastic Aspects of Dynamics. Springer. Lee, S.B., Kang, W., Sung, H.J., 2008. Organized self-sustained oscillations of turbulent flows over an open cavity. AIAA J. 46 (11), 2848–2856. Le, H., Moin, P., Kim, J., 1997. Direct numerical simulation of turbulent flow over a backward-facing step. J. Fluid. Mech. 330, 349–374. Lee, S.B., Seena, A., Sung, H.J., 2010. Self-sustained oscillations of turbulent flow in an open cavity. J. Aircraft 47 (3), 820–834. Lumley, J., 1967. The structure of inhomogeneous turbulent flows. In: Yaglam, A.M., Tatarsky, V.I. (Eds.), Proceedings of the international colloquium on the fine scale structure of the atmosphere and its influence on radio wave propagation, Moscow, Nauka, pp. 166–178. Lund, T.S., Wu, X., Squires, K.D., 1998. Generation of turbulent inflow data for spatially-developing boundary layer simulation. J. Comp. Phy. 140, 233–258.
Mezic, I., 2005. Spectral properties of dynamical systems, model reduction and decompositions. Nonlinear Dyn. 41, 309–325. Pereira, J.C.F., Sousa, J.M.M., 1995. Experimental and numerical investigation of flow oscillations in a rectangular cavity. J. Fluids Eng. 117, 68–74. Podvin, P., Fraingneau, Y., Lusseyran, F., Gougat, P., 2006. A reconstruction method for the flow past an open cavity. J. Fluids Eng. 128, 531–540. Rossiter, J.E., 1964. Wind-tunnel experiments on the flow over rectangular cavities at subsonic and transonic speeds. Aeronautical Research Council Reports and Memoranda No. 3438. Rowley, C.W., Colonius, T., Basu, A.J., 2002. On self-sustained oscillations in twodimensional compressible flow over rectangular cavities. J. Fluid Mech. 455, 315–346. Rowley, C.W., Mezic, I., Bagheri, S., Schlatter, P., Henningson, D.S., 2009. Spectral analysis of nonlinear flows. J. Fluid Mech. 641, 115–127. Ruhe, A., 1984. Rational Krylov sequence methods for eigenvalue computation. Lin. Alg. Appl. 58, 279–316. Sarohia, V., 1977. Experimental investigation of oscillations in flows over shallow cavities. AIAA J. 15 (7), 984–991. Schmid, P.J., 2009. Dynamic mode decomposition of experimental data. In: 8th International Symposium on Particle Image Velocimetry. Melbourne, Victoria, Australia, August 25–28. Schmid, P.J., 2010. Dynamic mode decomposition of numerical and experimental data. J. Fluid Mech. 656, 5–25. Spalart, P.R., 1988. Direct simulation of a turbulent boundary layer up to Reh = 1410. J. Fluid Mech. 187, 61–98.
Contents lists available at SciVerse ScienceDirect
International Journal of Heat and Fluid Flow journal homepage: www.elsevier.com/locate/ijhff
Dynamic mode decomposition of turbulent cavity flows for self-sustained oscillations Abu Seena, Hyung Jin Sung ⇑ Department of Mechanical Engineering, KAIST, 291 Daehak-ro, Yuseong-gu, Daejeon 305-701, Republic of Korea
a r t i c l e
i n f o
Article history: Received 31 May 2011 Received in revised form 16 August 2011 Accepted 23 September 2011 Available online 15 October 2011 Keywords: Self-sustained oscillations Dynamic mode decomposition Coherent structure
a b s t r a c t Self-sustained oscillations in a cavity arise due to the unsteady separation of boundary layers at the leading edge. The dynamic mode decomposition method was employed to analyze the self-sustained oscillations. Two cavity flow data sets, with or without self-sustained oscillations and possessing thin or thick incoming boundary layers (ReD = 12,000 and 3000), were analyzed. The ratios between the cavity depth and the momentum thickness (D/h) were 40 and 4.5, respectively, and the cavity aspect ratio was L/D = 2. The dynamic modes extracted from the thick boundary layer indicated that the upcoming boundary layer structures and the shear layer structures along the cavity lip line coexisted with coincident frequency space but with different wavenumber space, whereas structures with a thin boundary layer showed complete coherence among the modes to produce self-sustained oscillations. This result suggests that the hydrodynamic resonances that gave rise to the self-sustained oscillations occurred if the upcoming boundary layer structures and the shear layer structures coincided, not only in frequencies, but also in wavenumbers. The influences of the cavity dimensions and incoming momentum thickness on the self-sustained oscillations were examined. Ó 2011 Elsevier Inc. All rights reserved.
1. Introduction Turbulent flows characterized by a wide range of scales potentially include an infinite number of degrees of freedom. However, in many situations, the flow complexity may be reduced to a set of coherent features with a few characteristic structures possessing self-excited global modes in space and time. For the past several decades, Proper Orthogonal Decomposition (POD) has been extensively used to identify the coherent structures in turbulent flows (Lumley, 1967). POD extracts a complete orthogonal set of spatial eigenfunctions corresponding to the modes of the measured second-order correlation function. POD has proven to be an effective method for identifying dominant features and events in turbulent flows based on energy rank. POD modes are computed from the two point correlation tensor which is an average quantity and results in loss of some dynamic information in the data. Dynamic information is significant for flows in which a local absolute instability prevails within a finite region. Such systems can exhibit selfsustained resonant modes at specific complex frequencies. A global stability analysis yields a large stability matrix and may be computationally expensive. The Arnoldi method requires explicit knowl-
⇑ Corresponding author. Tel.: +82 42 350 3027; fax: +82 42 350 5027. E-mail address: [email protected] (H.J. Sung). 0142-727X/$ - see front matter Ó 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.ijheatfluidflow.2011.09.008
edge of the evolution operator describing the data, which makes this method unattractive. Schmid (2010) introduced a Dynamic Mode Decomposition (DMD) approach to extract dynamic mode information from a flow field based on the Koopman analysis of dynamical system (Rowley et al., 2009). The extracted dynamic modes, which may be interpreted as a generalization of the global stability modes, may be used to describe the underlying physical mechanism that is captured within the data sequence or to project large-scale problems onto a dynamical system with a significantly large number of degrees of freedom. The decomposition algorithm relies on the extraction of a low-dimensional subspace after fitting a high-degree polynomial to the data sequence without knowledge of the evolution operator from which the data was generated. The eigenvalues and eigenvectors of the low-dimensional subspace capture the principal dynamics of the flow. Although DMD assumes a linear evolution operator, even in nonlinear flows, it extracts the structures of a linear tangent approximation to the underlying flow and describes its dynamic behavior. The mathematics underlying this decomposition is related to the Koopman operator (Lasota and Mackey, 1994; Mezic, 2005; Rowley et al., 2009), where the Koopman modes reduce to the linear global eigenmodes for the linear process while related to the Fourier modes for the nonlinear periodic problems. The Koopman operator is related to other decompositions, as described by Schmid (2010). In the present
A. Seena, H.J. Sung / International Journal of Heat and Fluid Flow 32 (2011) 1098–1110
work, the DMD algorithm is applied to turbulent flows over an open cavity. Self-sustained oscillations have been recognized as the most important characteristic behavior of separated shear layers over open cavities since Rossiter (1964) investigated the influence of cavity dimension and free stream velocity on the resonant tones. To understand the mechanism underlying self-sustained oscillations and to elucidate the large-scale vortical structures responsible for such oscillations, numerous experimental and numerical studies have been carried out on incompressible flows as well as on compressible flows. As reported in previous studies of compressible flows (Rowley et al., 2002), the mechanism underlying self-sustained oscillations has been well-characterized by extracting large-scale vortical structures over an open cavity. In incompressible turbulent flows, however, it has not been easy to identify the coherent structures responsible for oscillations due to both incompressibility and turbulence, although self-sustained oscillations are predicted. Daoud et al. (2006) found that even if there are no resonant tones in the surface pressure spectra, the acoustic feedback can exist in cavity flows. Pereira and Sousa (1995) observed periodically oscillating shear layers in the flow of a turbulent incoming boundary layer over an open cavity. In PIV studies of turbulent open cavity flows, Kang et al. (2008) observed that the incoming turbulent boundary layer gave rise to large-scale vortical structures that were responsible for self-sustained oscillations. Instantaneous fields of velocity fluctuations were used in Ashcroft and Zhang (2005) mainly to identify largescale vortical structures in a separated shear layer, although the organized nature of the coherent structures was described in the statistical analysis of the correlations. The work of Lee et al. (2008), which was the large eddy simulation (LES) of high-Reynolds number incompressible turbulent flows over an open cavity, provided instantaneous isopleths and spanwise-averaged contours for pressure fluctuations in an effort to describe the large-scale vortical structures that were responsible for the peak frequency in the energy spectra. An analysis of raw data, including instantaneous velocity or vorticity distributions, usually requires that structural analysis of separated shear layers is instantaneous and qualitative. Therefore, application of an analysis approach such as DMD is required for a quantitative and dynamic characterization of large-scale vortical structures. The objectives of the present study were to (1) identify the large-scale vortical structures responsible for hydrodynamic oscillations by applying DMD to the pressure fluctuations of incompressible turbulent flows over an open cavity; (2) obtain dynamic information about the extracted structures; (3) elucidate the quantitative characteristics of the large-scale vortical structures by comparing the oscillating behavior of a separated shear layer to the behavior of a nonoscillating system. The sequences of images over time were processed to delve into the temporal dynamics of structures in turbulent open cavity flows. To achieve these objectives, numerical data from two different flows over an open cavity, with thin or thick incoming boundary layers, were analyzed by DMD at Reynolds number based on cavity depth ReD = 12,000 and 3000, respectively. The self-resonant oscillating modes were observed for the case of a thin incoming boundary layer. Analysis of the thin boundary layer revealed a dominant peak in the DMD spectrum, whereas all dynamic modes extracted from the thick boundary layer showed the gradual decay in energy with frequency representing viscous boundary layer structures. These results permitted examination of the effects of the upcoming momentum thickness on the self-sustained oscillations in an open cavity flow. Analysis was performed over the full computational domain as well as over a subdomain. The selection of the domain significantly influenced the results, depending on the flow conditions. Furthermore, DMD was applied to velocity flow distributions
1099
in which oscillating behavior was observed in a separated shear layer over the open cavity. The vortical structures responsible for self-sustained oscillations were thereby extracted. 2. Dynamic mode decomposition Dynamic mode decomposition is a recent extension of the classical Arnoldi technique. A temporal sequence of N data fields, consisting of column vectors, vj, that are equispaced in time according to the Nyquist criteria (Schmid, 2010), can be written as
VN1 ¼ fv 1 ; v 2 ; v 3 ; . . . ; v N g:
ð1Þ
The basic premise of the method is that each snapshot in time is assumed to be generated by a linear dynamical system,
vjþ1 ¼ Avj :
ð2Þ
The eigenvalues and eigenvectors of the matrix A completely characterize the behavior of the dynamical system. DMD is a method for computing the approximate eigenvectors or Ritz vectors of a system matrix. A high-degree polynomial is fit to a Krylov sequence of flow fields. As the number of snapshots increases, the flow is assumed to approach a linear dependency after a sufficient number of snapshots. The last image is the linear combination of the previous images and may be described by
vN ¼ c1 v 1 þ c2 v 2 þ . . . þ cN1 v N1 þ r;
ð3Þ
where r is the residual vector and c = {c1, c2, c3, . . ., cN1} are the coefficients,
vN ¼ VN1 c þ r: 1
ð4Þ
Eq. (4) represents an over-determined system of equations. The coefficients may be obtained using the least squares method. The number of required snapshots, N, may increase until the residual vector r converges. Following Ruhe (1984), the least squares description of the full system matrix A may be written as
Afv 1 ; v 2 ; v 3 ; . . . ; v N1 g ¼ fv 2 ; v 3 ; v 4 ; . . . ; v N g cg þ reTN1 : ¼ fv 2 ; v 3 ; v 4 ; . . . ; VN1 1
ð5Þ
Furthermore, Eq. (5) may be written as
AVN1 1
¼ VN2 ¼ VN1 C þ reTN1 : 1
ð6Þ
The matrix C is of the companion type,
0
0
c1
1
C B c2 C B1 0 C B C¼B : : : C C; B C B 1 0 c @ N2 A 1 cN1
ð7Þ
which simply shifts the snapshot (Schmid, 2009) index from 1 to N 1. This matrix, extracted from the data sequence, represents a low-dimensional system matrix representation of the full system matrix. The characteristic solution to the matrix C approximates some of the eigenvalues (kj) of the full system matrix A, called as empirical Ritz values (Rowley et al., 2009). The logarithmic mapping of these eigenvalues provides the growth rate (rj = Re{kj}/Dt) and frequency (xj = Im{kj}/Dt). The empirical Ritz values lying on the unit circle represent the modes with zero growth rates, whereas the eigenvalues lying inside and outside the unit circle represent the damped and undamped modes, respectively. If the data sequence is extracted from a nonlinear process, then neutrally stable eigenvalues are expected as the data sequence becomes sufficiently long. This is due to the fact that instabilities have nonlinearly saturated while decaying processes have vanished from the data sequence; only neutrally stable fluid elements remain. The
1100
A. Seena, H.J. Sung / International Journal of Heat and Fluid Flow 32 (2011) 1098–1110
eigenvectors provide a linear combination of basis coefficients that may be used to extract the corresponding dynamic modes as,
Uj ¼
VN1 Tj : 1
ð8Þ
In contrast to the POD where a spatial orthogonality of the identified structures is enforced at the expense of a mixture of frequencies in each individual POD mode, the basis vector Tj obtained above aims at an orthogonality in time by identifying pure frequencies. There is no normalization step in the algorithm, and hence the modes come with the amplitudes, ||U|| that provide a way to rank their contribution to the overall energy in the set VN1 (Bagheri, 2010; Rowley et al., 2009). The companion matrix C may be a highly non-normal matrix, which yields an ill-conditioned eigenvalue decomposition problem. To improve the accuracy of the DMD method, Bagheri (2010) suggested that the matrix should first be balanced by a similarity transformation, followed by reducing to an upper Hessenberg form via a second similarity transformation, and finally, the eigenvalues should be computed using the QR algorithm. 3. Numerical simulation For an incompressible flow, the non-dimensional governing equations are
@ui @ui uj @p 1 @ 2 ui ¼ iþ ; þ @t @xj @xi ReD @x2j @ui ¼ 0; @xi
LY
i ¼ 1; 2; 3
ð9Þ
ð10Þ
where xi are the Cartesian coordinates and ui are the corresponding velocity components. The free-stream velocity U1 and the cavity depth D were used to non-dimensionalize the equations. The Reynolds number was defined as ReD = U1D/m, where m is the kinematic viscosity. The governing Eqs. (9) and (10) were integrated in time using the fully implicit decoupling method proposed by Kim et al. (2002). A schematic diagram of the computational domain is shown in Fig. 1. All terms in Eq. (9) were advanced in time using the Crank– Nicholson method and were resolved using the second-order central difference scheme in space. A block LU decomposition was used to achieve velocity–pressure decoupling as well as additional decoupling of the intermediate velocity components in conjunction with an approximate factorization. The multi-grid method was applied to obtain a solution to the pressure Poisson equation over a Tshaped domain. In the present simulations, a turbulent boundary layer with realistic velocity fluctuations, which were generated using the method of Lund et al. (1998), was provided at the inlet. A direct numerical simulation (DNS) of incompressible flows over an open cavity was performed at ReD = 3000, where the DNS data were provided at the inlet with Reh = 300, 670 and 1410. The cavity flows at high Reynolds number (ReD = 12,000) were simulated using a LES with a dynamic subgrid-scale model (Germano et al., 1991). Here, DNS data at Reh = 300 were provided at the inlet. To analyze the influence of upstream turbulence on cavity flow, four cases with D/h = 2.1, 4.5, 10 and 40 were simulated. The simulation parameters are listed in Table 1, where O and X represent the flow with and without self-sustained oscillations, respectively. Among them, we selected two typical cases with D/h = 40 and 4.5 and applied DMD to understand hydrodynamic resonance and quantitative characteristics of the large-scale vortical structures. To ascertain the reliability of the incoming turbulent boundary layer, the turbulence statistics were compared with the DNS data of Spalart (1988), and the results are presented in Fig. 2. Fig. 2a shows the profile of the mean velocity normalized by the friction velocity (u+ = u/us), with variations in y+ = yus/v, and Fig. 2b shows
LZ L
Lo D
y x Li
z
Fig. 1. Schematic diagram of the computational domain.
a comparison of the turbulence intensities. The results were in excellent agreement with those of the previous study (Spalart, 1988). Uniform grids were used in the spanwise direction. Nonuniform grids were used in the streamwise and wall-normal directions, with higher densities of grid points near the leading and þ trailing edges. The sizes of Dxþ min and Dymin were determined by obtaining convergence of the turbulent intensities near the wall and the trailing edge. The computational time steps were freed from the Courant Friedrich Levy limit by employing the fully implicit decoupling method. The time step was fixed at Dt = 0.007D/U1 for ReD = 3000. For ReD = 12,000, the time step was very small, Dt = 0.001875D/U1, due to the explicit characteristics of the dynamic subgrid-scale model. The computational domain size was determined after checking the two-point correlations of the velocity fluctuations. The domain size in the spanwise direction (Lz) was four times the cavity depth (D) for both ReD = 3000 and 12,000. Selection of these spanwise domain sizes was supported by results from previous numerical simulations of turbulent flows over a backward-facing step (Le et al., 1997). The domain size in the streamwise direction was Li = Lo = 9D for ReD = 3000, and Li = 4D, Lo = 7.5D for ReD = 12,000. Here, Li and Lo are the streamwise lengths of the development section ahead of the leading edge and behind the trailing edge. The vertical height (LY) was set to exceed 9D for ReD = 3000, and the height was 5D for ReD = 12,000. The computational details are given in Lee et al. (2010). 4. Results and discussion Before the DMD algorithm was applied to turbulent cavity flows, it was validated using a fabricated pattern. Take, for example, the pattern q, such that
q ¼ q0 þ q1 þ q2 þ q3 ;
ð11aÞ
q0 ðx; y; tÞ ¼ expðy2 =0:7Þ; qn ðx; y; tÞ ¼ an ðtÞ
m¼1 X
ð11bÞ "
ð1Þm exp
m¼1
n ¼ 1; 2; 3; . . . ;
ðx bn m cn tÞ2 dn
!
#
þ
y2 ; dn ð11cÞ
where dn = anx + bn is the diameter of the structure, an and bn are constants, bn is the wavelength factor, cn is the velocity of the pattern in the streamwise direction, and an(t) is the growth or decay factor. The pattern adopted for validation was composed of patterns of three different sizes (n = 1, 2, and 3), each traveling with a fixed velocity. The constants and factors an(t) are given in Table 2. The fabricated pattern is shown in Fig. 3a, and consisted of the superposition of the three patterns with n = 1, 2, and 3. For n = 1, the factor
1101
A. Seena, H.J. Sung / International Journal of Heat and Fluid Flow 32 (2011) 1098–1110 Table 1 Simulation conditions. D/h
Reh
ReD
Nx Ny Nz
Dxþ ; Dxþ max min
Dyþ min
Dz +
Self-sustained oscillations
2.1 4.5 10 40
1410 670 300 300
3000 3000 3000 12,000
513 383 513 513 213 257 513 133 129 897 169 257
1.00, 1.09, 1.17, 1.40,
0.196 0.188 0.196 0.560
5 5 5 10
X X X O
16.2 17.7 19.0 40.0
Fig. 2. Mean streamwise velocity profiles and velocity fluctuations in the inflow data.
an(t) was set to unity such that the pattern shown in Fig. 3b neither grew nor decayed over time. The factors an(t), for n = 2 and 3, were selected such that the structures q2 decayed over time and q3 grew in time. Next, the fabricated pattern q was processed using the DMD algorithm, as described in Section 2. The extracted Ritz values kj of the low-dimensional system matrix C are shown in Fig. 4a. The size and color of the eigenvalues indicate the amplitude of the respective mode in the data sequence. The eigenvalues kj of C are transformed via the logarithmic mapping and shown in Fig. 4b, where the unstable eigenvalues appear in the right halfplane. The eigenvalues obtained were found to be in complex conjugate pairs. As a result, the DMD spectrum was symmetric about the frequency x = 0, as shown in Fig. 4b, where r represents the growth/decay rate and x is the frequency of each mode. The global energy norms of each mode are plotted as a function of frequency in Fig. 4c. Four energetic structures, in which x = 0 contained the highest energy, were extracted from the fabricated pattern. Four distinct modes are marked 1–4 in the extracted spectra shown in Fig. 4. Fig. 5 shows their respective modes. The eigenvalue that appeared at the origin, marked 1, accounted for the steady state contained in the pattern and is shown in Fig. 5b. The second DMD mode with the second highest energy is marked 2 in the spectra and lie on the neutral line (r = 0), i.e., neither growing nor decaying over time. The modes corresponding to the real and imaginary parts, respectively, are shown in Fig. 5c and d. They clearly represented the pattern q1 with a decay factor an(t) = 1. The third highest mode, marked 3, lay in the stable region (r < 0) and decayed over time. These modes are shown in Fig. 5e and f. The extracted Table 2 Constants for the fabricated pattern. n
an
bn
bn
cn
an(t)
1 2 3
0.030 0.010 0.005
0.005 0.100 0.100
0.8 0.4 0.3
0.377 0.252 0.126
1 et/80 0.1 (1 et/15) + 0.2
mode represented the pattern q2 with a negative decay factor (an < 0). Finally, the fourth mode, marked 4, was observed in the unstable region (r > 0), as shown in Fig. 5g, and represented q3 of the fabricated pattern. The above example illustrates how the DMD method decomposes the patterns in the data and reveals important dynamic information. Next, DMD was employed to describe the open cavity data. The flows corresponding to each of the two cases considered in the present study, ReD = 3000 and 12,000, were simulated. Fig. 6 shows the mean streamlines corresponding to the two simulations. Two large recirculation eddies are observed in both simulations. Mean cavity flow is characterized by the primary clockwise rotating vortex covering more than downstream half of the cavity which elongates to the leading edge. A weaker anticlockwise rotating secondary vortex is evident near the upstream corner, induced due to the primary vortex. After a long initial transient period, a sequence of images was saved. An initial transient time of more than 10 ‘‘flow-through’’ cycles was discarded to allow the passage of unphysical fluctuations. The temporal sequences of the system snapshots were processed using the DMD algorithm to extract pertinent dynamic characteristics of the flow. One hundred fifty instantaneous snapshots of the pressure fluctuations were used for the ReD = 3000 condition, and 124 snapshots were used for the ReD = 12,000 condition. The eigenvalues of the subspace C were expected to converge toward some of the eigenvalues of the full system matrix A. The DMD algorithm may be applied in one of two possible ways. It may be applied over the entire physical domain, including both the upstream region of the leading edge and the downstream region of the trailing edge, as shown in Fig. 7a, or it may be applied over a subdomain, as shown in Fig. 7b. Both approaches were considered in the present work, and the results were compared. Analysis of the subdomain region only, shown in Fig. 7b, focused on the vortical flow inside the cavity and the shear layer that formed and detached near the right wall of the cavity. Analysis of the full domain, as shown in Fig. 7a, captured the effects of the upstream boundary layer structures in addition to the cavity flow.
1102
A. Seena, H.J. Sung / International Journal of Heat and Fluid Flow 32 (2011) 1098–1110
Fig. 3. Fabricated pattern.
Fig. 4. (a) Ritz values kj, (b) DMD spectrum and (c) energy spectrum for the fabricated pattern.
Fig. 5. Modes extracted from the fabricated pattern by the DMD algorithm.
A. Seena, H.J. Sung / International Journal of Heat and Fluid Flow 32 (2011) 1098–1110
1103
Fig. 6. Mean flow streamlines.
Fig. 7. Physical domain for application of DMD: (a) full domain and (b) subdomain.
4.1. ReD = 3000 and D/h = 4.5 DNS of incompressible cavity flow over an open cavity were performed at ReD = 3000. A subdomain analysis was performed first. One hundred fifty images of the instantaneous pressure fluctuations at a sampling rate of Dt = 0.11725 were processed using the DMD algorithm. The eigenvalues of the subspace approximated the eigenvalues of the full system matrix and provided a quantitative illustration of the underlying temporal dynamics. The eigenvalues from the subdomain analysis are shown in Fig. 8a, where the symbol color indicates the global energy norm of the associated mode. Nearly all the Ritz values are on the unit circle ||kj|| = 1, indicating that the sample points vi lie on or near an attracting set. The norm of each mode indicates the energy in the corresponding mode. The energies of the extracted modes are shown as a function of frequency in Fig. 8b, where each mode is displayed with the vertical line scaled with its magnitude at its corresponding frequency. The eigenvalue marked 1that typically appears at the origin (x = 0) with the maximum energy accounts for the steady state and represents the mean flow as such depicted by a black symbol in Fig. 8a. The corresponding dynamic mode is shown in Fig. 9a, where the low-pressure region observed in the cavity indicates the primary vortex. A high-pressure region is also observed at the trailing edge,
possibly due to the impingement of shear layer vortical structures. The other modes show the spatial structures displaying periodic behaviors in time. A few specific modes in Fig. 8b dominate the energy spectrum at x = 1.9, 2.8 and 3.9 rad/s, marked 2–4. The modes are shown in Fig. 9b–d. The solid and dotted lines represent the negative and positive pressure fluctuations, respectively. In all of the modes shown in Fig. 9b–d, the pressure fluctuations are observed above the shear layer, suggesting that they correspond to the upcoming boundary layer structures flowing over the cavity and further impinging on the trailing edge. To gain further insight into these boundary layer structures, the DMD algorithm was applied over the full cavity domain; the obtained spectrum is shown in Fig. 10a. The symbol color indicates the norm of the corresponding mode. The black symbol at zero frequency, marked 1, represents the steady state and the corresponding mode is shown in Fig. 11a. Their norms against frequency are shown in Fig. 10b. Five distinct frequency peaks are marked 2–6, and their corresponding modes are shown in Fig. 11b–f. The eigenvalue marked 2 in the spectrum indicates a very low frequency mode. The pressure fluctuations of the corresponding mode, shown in Fig. 11b, indicate the presence of a very large-scale structure relative to the size of the cavity which convects over the cavity. The mode marked 3 is shown in Fig. 11c. In addition to the upcoming boundary layer structures, the structures generated due to shear layer oscillations could be observed along the cavity lip line, which further impinged on the trailing edge. The eigenvalue marked 4 is shown in Fig. 11d. The structures that appear in this particular mode suggest that at this frequency, the wavenumber corresponding to the upcoming boundary layer structure and the perturbation generated by shear layer oscillations are equal. The modes of two other eigenvalues marked 5 and 6 are displayed in Fig. 11e and f. Note that all the modes in Fig. 11 represent the incoming viscous boundary layer structures.
Fig. 8. Ritz values and their magnitudes at ReD = 3000 for the cavity.
1104
A. Seena, H.J. Sung / International Journal of Heat and Fluid Flow 32 (2011) 1098–1110
Fig. 9. DMD modes inside the cavity at ReD = 3000.
Fig. 10. Ritz values and their magnitudes at ReD = 3000 for the full domain.
However, as compared to subdomain analysis, the energy spectrum for the full domain analysis, shown in Fig. 10b, shows no leading modes, and the magnitude gradually decays as a function of frequency, suggesting the presence of a wide range of turbulent structures. This result is consistent with the absence of self-sustained oscillations. Analysis of the subdomain and full domain suggests that the selection of only a subdomain in the DMD analysis is inappropriate for understanding the dynamics pertinent to the flow. The DMD spectrum of the full domain reveals important characteristics at the boundary layer. The dominant features illustrate the boundary layer structures over the full domain at ReD = 3000. 4.2. ReD = 12,000 and D/h = 40 LES of incompressible cavity flow over an open cavity were performed at ReD = 12,000. Realistic velocity fluctuations of Reh = 300 were provided at the inlet. One hundred twenty-four images of the instantaneous pressures at a sampling rate of Dt = 0.075 were processed using the DMD algorithm, and a low-dimensional subspace was extracted. The analysis was performed for both cavity subdomain and full domain. The extracted spectrum, i.e., the spectrum of C, is displayed in Fig. 12a and b for subdomain and Fig. 12c and d for full domain. Fig. 12a and c shows that for both cases, all the Ritz
values are on the unit circle. The symbol color indicates the global energy norms of the modes. The amplitude of modes against frequency is shown in Fig. 12b and d. Both spectra, obtained from the subdomain and the full domain, displayed dominant peaks, marked 2 and 4 with two leading modes at x = 3.5 and 4.6 rad/s, respectively. Some of the selected modes from the spectrum in Fig. 12b, marked 1–4, are shown in Fig. 13. The eigenvalue at the origin, marked 1 in the spectrum (Fig. 12b), represents the steady state component of the pressure, and its respective mode is shown in Fig. 13a. The modes marked 2 and 3 with dominant peaks and their respective modes are shown in Fig. 13b and c, respectively. One pair of negative and positive distributions represents the low-pressure fluctuations of a large-scale vortical structure and the high-pressure fluctuations of induced rotational motions. Three pairs are clearly observed between the leading and trailing edges. This was consistent with the spectral characteristics of self-sustained oscillations corresponding to NR = 3, reported by Lee et al. (2008) The streamwise length scale of the coherent structures gradually increases from 0.3D in the region immediately downstream of the leading edge to 0.8D in the region that impinges on the trailing edge. To avoid ambiguity, the length scale of the coherent structure was determined to be twice as long as the streamwise distance between the central locations of the positive and negative distributions. The transverse length scale was not sufficiently large to affect the pressure
A. Seena, H.J. Sung / International Journal of Heat and Fluid Flow 32 (2011) 1098–1110
1105
Fig. 11. DMD modes inside the cavity full domain at ReD = 3000.
Fig. 12. Ritz values and their magnitudes at ReD = 12,000 for (a) the cavity and (b) the full domain.
fluctuations on the bottom wall inside the cavity. These modes suggest the presence of self-sustained oscillations in the cavity. The symbol marked 4 represents next leading mode, and its respective mode is shown in Fig. 13d.
We performed POD to the same data set for comparison. Fig. 14 shows the first four POD eigenmodes and the solid and dotted lines represent the positive and negative distributions of eigenmodes. In the first and second modes, three pairs of well-organized
1106
A. Seena, H.J. Sung / International Journal of Heat and Fluid Flow 32 (2011) 1098–1110
Fig. 13. DMD modes inside the cavity at ReD = 12,000.
structures are observed similar to the DMD modes obtained in Fig. 13. These first two modes look like two different phase states of an identical motion as expected for the flow in experiencing global mean advection. The most energetic POD mode in Fig. 14a captured a similar spatial structure to the most energetic DMD mode shown in Fig. 13b, but the difference is shown in the time coefficients. Variations of the mode coefficient for the most energetic POD mode (Fig. 14a) and DMD mode (Fig. 13b) are shown in Fig. 15a, which are the coefficients used when reconstructing the flow field with isolated modes. Fig. 15a shows that a single DMD mode contains by construction only a single frequency component, whereas the POD modes capture the most energetic structures, resulting in the mode that contain several frequencies. This is apparent in the power spectrum of time-varying coefficients shown in Fig. 15b, where several peaks are observed in the case
of most energetic POD mode, whereas a single peak is observed in the case of most energetic DMD mode. Both methods show the dominant peak at f = 0.765 Hz, hence validating that the mode extracted in Fig. 13b is responsible for the self-sustained oscillations. To determine whether these pressure fluctuation modes arise from shear layer oscillations or from upcoming boundary layer structures, DMD analysis was performed over the full cavity domain. The spectrum for the full domain is shown in Fig. 12d. The spectrum obtained from the full domain analysis resembles the spectrum from the subdomain analysis Fig. 12b. The modes corresponding to the dominant peaks, marked 1–4 in Fig. 12d, are shown in Fig. 16. These modes resemble the modes obtained from the subdomain analysis. The size of structures along the shear layer is comparable to that of the upstream boundary layer structures.
Fig. 14. POD modes inside the cavity at ReD = 12,000.
Fig. 15. (a) Time coefficients of the DMD and POD modes and (b) power spectra of a(t).
A. Seena, H.J. Sung / International Journal of Heat and Fluid Flow 32 (2011) 1098–1110
The modes in Fig. 16b–d indicate that the structures along the shear layer coincide with the wavenumber spectra of the upstream boundary layer structures. This coincidence may have produced self-resonant modes, which will be discussed later. Although alternating pressure fluctuation patterns were regularly observed in the first and second modes, the distributions in the other modes were irregular. The other modes most likely describe the fluctuating behavior of a separated shear layer with a high wavenumber, albeit with slight variations. In previous studies of laminar cavity flows, the oscillating behavior of a separated shear layer was represented only by the first two modes, whereas the next higher modes were associated with recirculation of the primary vortex inside the cavity (Podvin et al., 2006). The simple modes of a separated shear layer produce a single peak frequency in the power spectra of laminar cavity flows. In the present study of turbulent cavity flows, however, the fluctuating behavior of the separated shear layer was represented in the higher modes as well as in the first two modes. The complex mixture of several modes produces a broad spectrum, corresponding to the pressure fluctuations in the fluctuating behavior of the separated shear layer, although the first two modes are responsible for the peak frequency corresponding to NR = 3. Fig. 17 shows the time history of the dynamic mode responsible for self-sustained oscillations. The time period T is the convection time for the vortical structure to convect along the shear layer from the leading edge to the trailing edge. The sequence of contour plots in Fig. 17a–d shows the convective patterns of large-scale vortical structures. The negative pressure fluctuations represent a clockwise rotating large-scale vortical structure (Lee et al., 2010). At time t = to, marked V1 indicates a large-scale vortical structure at the upstream of the leading edge and marked V2 represents another vortical structure located at (x/D, y/D) = (0.5, 1.0). After separation at the leading edge, the vortical structure attains energy due
1107
to the shear layer oscillations and enters the cavity shear layer. As time progresses to t = to + T/2, these structures convect along the shear layer. For t = to + 3T/4, the fringe of marked V2 begins to impinge on the trailing edge. The contour plot shows the partial clipping of marked V2 and continuous convection of marked V1. A small part of the clipped marked V2 is entrained into the cavity along the vertical wall whereas the majority is ejected out of the cavity. After escaping from the cavity, marked V2 gradually dissipates along the downstream wall. The coherent structure responsible for self-sustained oscillations was investigated through DMD analysis of the velocity data. A velocity field corresponding to each DMD mode of the pressure fluctuations was obtained. The coherent structures extracted for the leading modes at x = 4.6 and 3.5 rad/s are shown in Fig. 18a and b, respectively. This is consistent with the velocity field corresponding to the most energetic POD mode in the work of Lee et al. (2008) for the same data set. 4.3. Self-sustained oscillation modes The disturbance in the shear layer across the cavity generates coherent structures that impinge on the trailing edge. The system exhibit self-excited global modes at specific complex frequencies that are globally unstable. A key question is whether local disturbances due to incident turbulence can prevent the occurrence of global instabilities that originate from the singularity of the cavity geometry. DMD analysis suggested that the upcoming boundary layer thickness was significantly responsible for the self-sustained oscillations, and, therefore, the parameter D/h was important in a description of open cavity flows. The unsteady modes were identified at ReD = 12,000 for a thin incoming boundary layer. The dynamic modes extracted from the thick boundary layer (ReD = 3000) showed that both the boundary layer structures and
Fig. 16. DMD modes inside the cavity full domain at ReD = 12,000.
Fig. 17. Dynamic patterns of the mode marked 2R in the spectrum.
1108
A. Seena, H.J. Sung / International Journal of Heat and Fluid Flow 32 (2011) 1098–1110
Fig. 18. Streamlines showing the selected dynamic modes at ReD = 12,000 for the cavity full domain.
the shear layer structures coexist with similar frequency space but with different wavenumber space. In contrast, structures with a thin boundary layer (ReD = 12,000) showed complete coherence among the modes, which produced self-sustained oscillations. These results suggest that the hydrodynamic resonance that gives rise to self-sustained oscillations occurs when the upcoming boundary layer structures and the shear layer structures assume identical frequency space and wavenumber space. Wavenumber coincidence relies on the size of the largest structure in the upcoming boundary layer which depends on Reh, such that its size is comparable to the shear layer dimensions. Resonance may occur if the size of the largest structure of the upstream boundary layer is comparable to the wavelength of the shear layer oscillation. The cavity dimensions (length and depth) play a significant role in the shear layer dimensions. Resonance may occur when the shear layer dimensions are equal to the length scales of the two largest structures of the upcoming boundary layer. The largest structure in the upcoming boundary layer is comparable to the boundary layer thickness. In light of the above analysis, the minimum conditional requirement for restraining self-resonant modes relies on the cavity dimension. The cavity length should be less than twice the wavelength of the largest and most energetic structure in the upcoming boundary layer. The cavity depth should be larger in size than the upcoming largest boundary layer structure. This is consistent with the finding of Sarohia (1977) that no self-sustained oscillations occur below a certain critical cavity length and depth. The number of vortical structures can be evaluated using the modified form of Rossiter’s empirical relation,
L fL ¼ ¼ NR ; kx U C
ð12Þ
where kx is the streamwise wavelength of oscillation, f is the oscillating frequency, Uc is the convection velocity, and NR is the number of wavelengths of the fundamental frequency that may be contained within the cavity length in the NRth oscillation mode (Gharib and Roshko, 1987). Eq. (12) is only valid in the limit of incompressible flow because it was derived under the assumption that the lag time is zero between the impact of a hydrodynamic structure on the trailing edge and the emission of an acoustic wave. Eq. (12) also supports the present suggestion in the sense that a minimum critical cavity length of the order of wavelength of structures is required to generate self-sustained oscillations, or in other words the ratio of cavity length to wavelength should be of positive integer order describing its mode number. For each mode, the average wavelength of the structures was measured along the shear layer. To avoid ambiguity, the wavelength (kx) of the coherent structure was determined to be twice as long as the streamwise distance between the central locations of the positive and negative distributions. Fig. 19a and b plots the wavelength as a function of the frequency for both cases. The DMD algorithm extracted a wide range of structures. We found that the large structures were associated with lower frequencies, and the frequency increased as the wavelength of the structures decreased. This result agrees well with Eq. (12), which describes the wavelength as the reciprocal of the frequency of structures with the same convection velocity. Fig. 19c and d shows the frequency on a log–log scale. The slope 1 indicates that the
Fig. 19. Wavelength of the structures as a function of frequency at ReD = 3000 and 12,000.
A. Seena, H.J. Sung / International Journal of Heat and Fluid Flow 32 (2011) 1098–1110
1109
Fig. 20. Wavenumber spectra at ReD = 3000 and 12,000.
wavelength and frequency follow a reciprocal relationship. In both cases, ReD = 3000 and 12,000 (Fig 19c and d, respectively), the reciprocal trend suggests that structures convect downstream with the same velocity. All modes extracted using the DMD algorithm indicated that the flow was characterized by turbulent structures. The L2-norm of the modes yielded the energies in the associated modes. The energy (||U||) is plotted against wavenumber in Fig. 20a and b for both cases analyzed in the present work. The Kolmogorov power law (5/3) is also plotted (solid line) in the same figure. The high value of the power law prefactor at ReD = 12,000 shows that it contains higher energy modes than the cavity flow at ReD = 3000. The energies are plotted on a log–log scale in Fig. 20c and d. The region in which the data agree well with the Kolmogorov power law corresponds to the turbulence cascade region called the inertial range. The inertial range was broader for the higher Reynolds number flows, ReD = 12,000, than for the thick boundary layer flow (ReD = 3000). The modes extracted using the DMD algorithm showed a turbulent energy cascade, thus supporting the turbulent nature of the flow. 5. Conclusions In the present study, we identified large-scale vortical structures responsible for self-sustained oscillations by employing a dynamic mode decomposition of the pressure fluctuations of turbulent flows over an open cavity. DMD was applied to incompressible turbulent flows over an open cavity at ReD = 3000 and 12,000, with upstream turbulence of Reh = 670 and 300, respectively, provided at the inlet. The influence of the incoming turbulent boundary layer on the self-sustained oscillations was investigated. A wide range of structures was extracted using the DMD algorithm. The DMD algorithm identified the dynamic mode which represented the coherent structures responsible for the self-sustained oscillations. The cavity with thick incoming boundary layer represented the viscous boundary layer structures, whereas the thin incoming boundary layer showed a dominant
peak in the spectrum which represented the cavity dynamics. A mode with the dominant peak at ReD = 12,000 represented the structure responsible for self-sustained oscillations. After processing the instantaneous velocity distributions, the large-scale vortical structures responsible for self-sustained oscillations were identified. These structures were consistent with the spectral characteristics of self-sustained oscillations corresponding to NR = 3. The dynamic modes extracted from the thick boundary layer (ReD = 3000) showed that both the boundary layer structures and the shear layer structures coexist in similar frequency space but with different wavenumber space. These results suggest that the hydrodynamic resonance that produces self-sustained oscillations occurs when the upcoming boundary layer structures and the cavity shear layer structures coincide not only in frequency space, but also in wavenumber space. Hence, the upcoming boundary layer thickness was found to be a significant parameter for describing self-sustained oscillations. The parameter D/h, for flows over an open cavity, was, therefore, found to be important. The DMD energy spectra agreed well with the Kolmogorov power law, supporting the conclusion that the extracted modes represent the turbulence scales. Selection of the domain for analysis using the DMD algorithm significantly affected observation of the dynamics pertinent to the flow.
Acknowledgements This study was supported by the Creative Research Initiatives (No. 2011-0000423) program of the National Research Foundation of Korea.
References Ashcroft, C., Zhang, X., 2005. Vortical structures over rectangular cavities at low speed. Phys. Fluids 17, 015104. Bagheri, S., 2010. Analysis and control of transitional shear flows using global modes. PhD Thesis, February, Department of Mechanics, Royal Institute of Technology, Sweden.
1110
A. Seena, H.J. Sung / International Journal of Heat and Fluid Flow 32 (2011) 1098–1110
Daoud, M., Naguib, A.M., Bassioni, I., Abdelkhalek, M., Ghoneim, Z., 2006. Microphone-array measurements of the floor pressure in a low-speed cavity flow. AIAA J. 44 (9), 2018–2023. Germano, M., Piomelli, U., Moin, P., Cabot, W.H., 1991. A dynamic subrid-scale eddy viscosity model. Phys. Fluids A 3, 1760–1765. Gharib, M., Roshko, A., 1987. The effect of flow oscillations on cavity drag. J. Fluid Mech. 177, 501–503. Kang, W., Lee, S.B., Sung, H.J., 2008. Self-sustained oscillations of turbulent flows over an open cavity. Expt. Fluids 45 (4), 693–702. Kim, K., Baek, S.-J., Sung, H.J., 2002. An implicit velocity decoupling procedure for the incompressible Navier-Stokes equations. Intl. J. Numer. Methods Fluids. 38 (2), 125–138. Lasota, A., Mackey, M.C., 1994. Chaos, Fractals and Noise: Stochastic Aspects of Dynamics. Springer. Lee, S.B., Kang, W., Sung, H.J., 2008. Organized self-sustained oscillations of turbulent flows over an open cavity. AIAA J. 46 (11), 2848–2856. Le, H., Moin, P., Kim, J., 1997. Direct numerical simulation of turbulent flow over a backward-facing step. J. Fluid. Mech. 330, 349–374. Lee, S.B., Seena, A., Sung, H.J., 2010. Self-sustained oscillations of turbulent flow in an open cavity. J. Aircraft 47 (3), 820–834. Lumley, J., 1967. The structure of inhomogeneous turbulent flows. In: Yaglam, A.M., Tatarsky, V.I. (Eds.), Proceedings of the international colloquium on the fine scale structure of the atmosphere and its influence on radio wave propagation, Moscow, Nauka, pp. 166–178. Lund, T.S., Wu, X., Squires, K.D., 1998. Generation of turbulent inflow data for spatially-developing boundary layer simulation. J. Comp. Phy. 140, 233–258.
Mezic, I., 2005. Spectral properties of dynamical systems, model reduction and decompositions. Nonlinear Dyn. 41, 309–325. Pereira, J.C.F., Sousa, J.M.M., 1995. Experimental and numerical investigation of flow oscillations in a rectangular cavity. J. Fluids Eng. 117, 68–74. Podvin, P., Fraingneau, Y., Lusseyran, F., Gougat, P., 2006. A reconstruction method for the flow past an open cavity. J. Fluids Eng. 128, 531–540. Rossiter, J.E., 1964. Wind-tunnel experiments on the flow over rectangular cavities at subsonic and transonic speeds. Aeronautical Research Council Reports and Memoranda No. 3438. Rowley, C.W., Colonius, T., Basu, A.J., 2002. On self-sustained oscillations in twodimensional compressible flow over rectangular cavities. J. Fluid Mech. 455, 315–346. Rowley, C.W., Mezic, I., Bagheri, S., Schlatter, P., Henningson, D.S., 2009. Spectral analysis of nonlinear flows. J. Fluid Mech. 641, 115–127. Ruhe, A., 1984. Rational Krylov sequence methods for eigenvalue computation. Lin. Alg. Appl. 58, 279–316. Sarohia, V., 1977. Experimental investigation of oscillations in flows over shallow cavities. AIAA J. 15 (7), 984–991. Schmid, P.J., 2009. Dynamic mode decomposition of experimental data. In: 8th International Symposium on Particle Image Velocimetry. Melbourne, Victoria, Australia, August 25–28. Schmid, P.J., 2010. Dynamic mode decomposition of numerical and experimental data. J. Fluid Mech. 656, 5–25. Spalart, P.R., 1988. Direct simulation of a turbulent boundary layer up to Reh = 1410. J. Fluid Mech. 187, 61–98.