Numerical study of three dimensional acoustic resonances in open cavities at high Reynolds numbers

JID:AESCTE AID:3323 /FLA

[m5G; v1.153; Prn:28/05/2015; 15:49] P.1 (1-11)

Aerospace Science and Technology ••• (••••) •••–•••

Contents lists available at ScienceDirect

Aerospace Science and Technology www.elsevier.com/locate/aescte

Numerical study of three dimensional acoustic resonances in open cavities at high Reynolds numbers M. Kompenhans, E. Ferrer ∗ , M. Chavez, E. Valero ETSIA-UPM – School of Aeronautics, Plaza Cardenal Cisneros 3, E-28040 Madrid, Spain

a r t i c l e

i n f o

Article history: Received 14 May 2014 Received in revised form 5 May 2015 Accepted 17 May 2015 Available online xxxx Keywords: Open cavity Compressible Navier–Stokes Unsteady RANS Implicit Large Eddy Simulation Acoustic resonance Rossiter modes

a b s t r a c t We present two and three dimensional numerical simulations over open cavities at high Reynolds (Re = 8.6 × 105 ) and various Mach numbers (0.2 < M < 0.8). Various turbulence closures are compared and the effect of the Mach number and three-dimensionality assessed. Results show that an implicit LES approach provides more accurate results that classic RANS models (e.g. k–epsilon variants) when analyzing the spectra associated to Rossiter’s acoustic resonances (regular self-sustained oscillations). In addition, we show that the simulations require relatively high Mach numbers to compare favorably to experimental data, showing the importance of compressibility effects. Finally, to capture qualitatively the dynamical content, we show that 2D simulations are not accurate enough and three dimensional simulations are necessary. © 2015 Elsevier Masson SAS. All rights reserved.

1. Introduction Predictions of noise generation (or acoustic sources) require the numerical integration of the unsteady compressible Navier–Stokes equations in complex geometries. Acoustic analogies, e.g. Lighthill [17], Ffowcs Williams–Hawkings [36], have enabled a certain degree of decoupling between the generation and propagation of noise. Following these analogies, equivalent acoustic sources can be extracted from the flow. These can be subsequently modeled and injected into a new simulation to predict the propagations of such sources. This decoupling enables fast computations, since source generation that requires highly accurate Navier–Stokes flow computations can be performed separately from the associated computations for the wave propagation. For a review of available techniques to compute aeroacoustic noise generation and propagation, the reader is referred to the monograph of Goldstein [14]. The present work is concerned with the generation of noise which requires accurate temporal and spatial numerical simulations. Indeed, to extract small acoustic scales demands fine meshes and statistically converged solution. Recently, the numerical community has drifted efforts from classical Reynolds Averaged Navier–Stokes (RANS) turbulent models, towards Large Eddy Simulation (LES) to predict unsteady aeroa-

*

Corresponding author. E-mail addresses: [email protected] (M. Kompenhans), [email protected] (E. Ferrer), [email protected] (M. Chavez), [email protected] (E. Valero). http://dx.doi.org/10.1016/j.ast.2015.05.008 1270-9638/© 2015 Elsevier Masson SAS. All rights reserved.

coustic sources (e.g. [7,33]). RANS methods consider time averaging of the NS equations where the instantaneous velocities are decomposed into their temporal mean and fluctuating components (see [35] for a detailed description), whilst LES approaches rely on spatial filtering where the large structures are resolved, reducing modeling to the small turbulent structures (i.e. small eddies), which are considered to behave in a homogeneous isotropic fashion. LES predictions show enhanced accuracy (compared to RANS methods) when computing time evolving flows such as the ones found in aeroacoustic problems [7,25]. However, it is agreed that the computational cost associated to performing LES simulations is still excessive for industrial standards (e.g. fast design cycles) and less computationally demanding solutions are needed. Towards this aim, we compare various turbulence models and asses the effect of three dimensionality when predicting tonal frequencies. Both RANS and LES methods require adequate turbulence closure models and fine meshes to capture the small fluid structures that generate sound. RANS solvers are generally used for computing steady solutions (i.e. time averaged solution) and have shown good accuracy, at a relatively low cost, when predicting aerodynamic forces for attached and mildly detached flows [35]. For unsteady, or detached flows, the steady RANS approach is no longer valid, but, the unsteady version URANS (i.e. Unsteady RANS) can be used and has shown to provide valuable results [26,35]. Although for unsteady flows, LES approaches have shown to provide better

JID:AESCTE

2

AID:3323 /FLA

[m5G; v1.153; Prn:28/05/2015; 15:49] P.2 (1-11)

M. Kompenhans et al. / Aerospace Science and Technology ••• (••••) •••–•••

prediction capabilities, from a practical point of view URANS may (or may not) be able to capture the main acoustic tones. Original LES methods derive new filtered equations by applying a spatial filtering to the NS equations [25]. Filtering allows for scale separation: large scales (i.e. large eddies) to be resolved and small scales (i.e. subgrid scale) to be modeled. In recent years a new LES concept has emerged. The Implicit Large Eddy Simulation (ILES) approach [2,15], sometimes referred to as Monotonic Integrated Large Eddy Simulation (MILES). This technique uses the numerical dissipation inherited from the numerical scheme (e.g. from upwinding the non-linear terms) to account for subgrid scale effects. In this method it is not necessary to derive new equations to compute turbulent flows (e.g. filtered NS equations as in explicit LES), but the Direct Numerical Simulation (DNS–NS) equations can handle high Reynolds number flows provided that enough numerical dissipation is present in discretise DNS–NS equations. Indeed, this technique does not require explicit modeling of the turbulent viscosity, as in explicit LES approaches (originally proposed by Smagorinsky [28], see [25] for a review) and provides valuable results even for challenging problems such as the prediction of transitional [31] and detached [10] flows. To obtain valid implicit LES schemes, one needs to ensure a dissipative monotone character (i.e. non-oscillatory) of the discretization which can be provided for example using a Flux-limiting or Total Variation Diminishing (TVD) schemes when discretizing the non-linear terms in the NS equations [2,15]. Note that in a valid implicit LES schemes (if a finite volume discretization is used) the discretization error appears in a divergence form and hence mimics the behavior of explicit subgrid models in explicit LES methods. Further information can be found in the recent review on how to compute high speed flows [22]. A final caveat arises when considering two or three dimensional simulations when using RANS or LES models. Although the 3D nature of turbulent flows is commonly agreed, a large number of 2D computations have managed to successfully predict tonal frequencies [8] at a reduced computational cost. This results partially motivate the work presented hereafter where we compared 3D ad 2D implicit LES simulations. Finally, our work assesses the effect of compressibility effects (varying Mach number) in the prediction of tonal frequencies using Implicit LES models. To summarize, we present comparisons of two and three dimensional computations for various turbulence models including Implicit LES and URANS models (two k–epsilon variants). These models are compared for the well known but challenging case of an open cavity at Mach M = 0.8 and Reynolds number of Re = 8.6 × 105 [2,7,12,16,30]. This cavity flow was first experientially analyzed by Forestier et al. [12] and subsequently numerically reproduced using highly accurate LES computation by Larchevêque et al. [16]. These published results will be used for comparison in our work. The underlying physical mechanism responsible of the tonal noise in this configuration has been attributed [12,16] to a self sustained pressure oscillation (i.e. Rossiter modes) with a fundamental frequency of about 2000 Hz. The problem is described in details in the following section but the interested reader is referred to the more complete overview of Gloerfelt [13]. The cavity flow problem has been studied previously but a few questions are still unanswered. The following bullet points summarize the main novelties, aims and clarifications included in our work:

• We clarify the effect of turbulence modeling and show that only the ILES computations are able to predict the complex flow of this challenging case.

• We show that both 2D and 3D ILES computations are able to predict the experimental frequencies but that the tonal amplitudes are only accurately predicted using 3D computations. • We show that eddy viscosity models (also with anisotropic corrections) are unable to capture the feedback mechanism described by Rossiter. • We show that for Mach numbers below 0.5, the cavity physics are different. Hence compressibility effects are crucial to predict accurately the experimental frequencies. This paper is organized as follows, in Section 2 a short review of the physical phenomena and the experimental set-up is introduced. The computational method and turbulence modeling strategies are described in Section 3. Finally, Sections 4 and 5 present and analyze the numerical solutions, which are compared to published experiential data. In Section 5, we first discuss the influence of the turbulence model to then analyze the compressibility effects and Mach number influence in the results. 2. Problem description Open cavity flows have received a great deal of attention [8, 13], partially because the configuration is relatively simple but also because the resulting flow presents interesting and complex behaviors. In addition, this configuration presents an interesting framework for the development of active and passive flow control strategies [4,8,20]. The physical phenomenon at play is a closed loop (feedback) mechanism of self sustained oscillation. The flow leaving the leading edge (left top corner) of the cavity presents a sheared velocity profile (or a set of vortices) along the top of the cavity that impinges at the trailing edge (right top corner). When this sheared profile interacts with the corner, a hydraulic pressure wave bounces back and excites the sheared profile near the leading edge (a region of high receptivity) producing or enhancing instabilities of Kelvin–Helmholtz (K–H) type in the sheared profile. This is a two dimensional self exiting mechanism named after Rossiter [24], that leads to well defined tonal frequencies (see Eq. (4) in following sections). In addition, global cavity modes (e.g. depth-modes, spanwise modes or tunnel modes) that may be either two or three dimensional can be excited within the cavity (depending on the flow conditions) [8,38] and may interact with the previous Rossiter phenomenon. Two very recent contributions have analyzed open cavity flows in terms of linear global stability analysis [20,38]. The first has focused on the Mach number effects in laminar flow conditions, showing that K–H instabilities are weakened for high Mach numbers. This conclusion is in accordance with our non-linear results as will be shown hereafter. In [20], Mettot et al. study for the first time the linear stability associated to the high Reynolds and Mach M = 0.8 using a k–omega and Spallart–Almaras turbulence models. In their work, URANS computations compare favorably with linear stability results showing that the latter technique is appropriate to shed light in this type of complex Rossiter’s mechanism. In addition, linear stability shows that when the frequencies of the K–H phenomenon and the acoustic waves coincide (i.e. resonance or lock-on state), the amplification of the power spectra is maximum. However, in our work (see following sections) we find that URANS methods are unable to capture the dynamical content for this configuration. A sketch of the geometry selected for the simulations is shown in Fig. 1. The geometry described here is a cavity of length L = 50 mm, depth D = 120 mm and width l = 120 mm, resulting in low aspect ratios: L / D = L /l = 0.42. This cavity may be considered deep (L / D < 1) and two-dimensional (L /l < 1).

JID:AESCTE AID:3323 /FLA

[m5G; v1.153; Prn:28/05/2015; 15:49] P.3 (1-11)

M. Kompenhans et al. / Aerospace Science and Technology ••• (••••) •••–•••

3

scheme, we find necessary to provide an overview of the numerics used in the solver ISAAC. Additional details may be found in [21] and the references therein. The Navier–Stokes system is spatially discretized using a secondorder finite-volume method. Advective non-linear terms (including turbulence equations) are discretized using Roe’s approximate Riemann solver coupled with a third order MUSCL interpolation scheme. The min-mod limiter is used to avoid spurious oscillations in regions of strong gradient and ensures the monotonicity required by MILES methods [15]. Diffusive fluxes are approximated with central differences and derivatives of diffusive flux are approximated using Gauss’s divergence theorem (further details on finite volume methods can be found in [11,32]). Flow variables required at this interface are obtained from arithmetic averaging of neighboring cells. To cope with complex geometries a multiblock procedure is used. To avoid very small time steps associated to the stiffness of the turbulent viscous and compressible Navier–Stokes system, an implicit time integration is used. The implicit scheme applied to Eq. (1) yields:

Fig. 1. Sketch for the open cavity configuration.

[I + α Ai ∂i ](Un+1 − Un ) = R(Un , Un−1 ),

The cavity is immersed in a high subsonic flow characterized by a Mach number M = 0.8. The Reynolds number based on the length of the cavity L is Re = u ∞ L /ν = 8.6 × 105 , with a reference velocity u ∞ = 258 m/s and kinematic viscosity ν for air. The reference speed of sound is a∞ = 322.5 m/s. The initial momentum thickness of the boundary layer (δm ) measured 1 mm ahead of the cavity edge (leading edge or top left corner) is equal to 6.48 × 10−4 m, giving a ratio of cavity length to boundary momentum thickness of L /δm = 77.1. In addition to the reference problem with Mach number M = 0.8, we provide results for the same configuration at three different Mach numbers M = 0.2, 0.4 and 0.6. In all cases, the Reynolds number has been adjusted to Re = 8.6 × 105 and the input boundary profile has been modified to maintain the boundary layer momentum thickness defined previously (L /δm = 77.1).

where Un is the vector of conservative variables at step/time n and Ai = ∂ Fi /∂ Ui is the Jacobian matrix. The parameters α , β are det

t fined as α = (1ϑ

+ϕ ) , β = (1+ϕ ) . It can be seen that varying ϑ, ϕ , one obtains different temporal schemes. Namely, ϑ = 1, ϕ = 0 provide the first order Euler implicit scheme, ϑ = 1/2, ϕ = 0 yield the trapezoidal implicit rule and the three point backward second order implicit scheme is obtained by setting ϑ = 1, ϕ = 1/2. In this work, the latter second order implicit scheme is retained. System (2) is solved by an Alternate Direction Implicit (ADI) methodology, which gives an approximate factorization of the form:

3. Computational method

[I + α A1 ∂x ][I + α A2 ∂ y ][I + α Az ∂z ](Un+1 − Un )

In this work, the well validated compressible Navier–Stokes solver ISAAC [21] has been used for simulation. This 2D/3D block structured finite volume code includes various turbulence models and has been used in the past to asses the accuracy of turbulence models in various regimes and configurations (see for example [5,37]). The compressible Navier–Stokes equations (including RANS terms) can be written in matrix form as:

∂t U + ∂i Fi (U) = 0,

i = 1, 2, 3

(1)

where U represents the vector of conservative variables and Fi the convective, diffusive and turbulence fluxes in the three coordinate directions. Einstein’s summation convention over repeated indices has been used. The dimension of the system depends on the turbulence model. When no turbulence model is enabled (ILES or MILES case) five equations are solved with the five classic compressible unknowns (density, 3 momentum components and energy). If a two-equation turbulence model is enabled, such as the k–ε (k–eps) or the explicit algebraic stress model (k–eps ASM), then seven equation are required and two additional variables for the turbulent kinematic energy (k) and the turbulent dissipation (eps) are added to the system. 3.1. Solver implementation Since Implicit LES computations do not rely on an explicit modeled turbulent viscous terms but on the numerical viscosity of the

R(Un , Un−1 ) = −β(∂i Fi ) +

ϕ

1+ϕ

(Un − Un−1 ),

= R(Un , Un−1 ),

(2)

(3) n+1

where the term α (U − U )  O ( t ) has been neglected, which results in second order temporal accuracy. This spatial discretization, Eq. (3), results in a set of independent tridiagonal systems per block, which can be solved with a generalization of the Thomas algorithm [21]. To avoid inaccuracies if large time steps are selected in the implicit scheme, a pseudo-time iterative algorithm has been used. Thus, instead of solving system Eq. (2) directly, an additional term 2

n

3

∂Uk ∂ τ is added at each physical time step. By doing so, we ensure

that steady solutions in τ , are unsteady solutions in t. Furthermore, this unsteady treatment enables the use of geometric Vmultigrid acceleration coupled to an ADI smoother. In this work, we have calibrated the numerical scheme to efficiently and accurately solve the flow over the open cavity. The resulting set of parameters include three level multigrid with ten implicit ADI sub-iterations in pseudo-time. 3.2. Turbulence modeling Simulations using two RANS models are included in this work. Namely, the k–ε (henceforth k–eps) model for high velocity flows of Zhang et al. [39] and the explicit Algebraic Reynolds Stress Model of Abid et al. [1] (henceforth k–eps ASM). In addition, Implicit LES (ILES or MILES) simulations [2,15] have been considered. In the latter, a Baldwin–Lomax turbulence model (i.e. two-layer algebraic zero equation model [35]) is used near the wall (incoming boundary layer) whilst no explicit model is enabled

JID:AESCTE

AID:3323 /FLA

4

[m5G; v1.153; Prn:28/05/2015; 15:49] P.4 (1-11)

M. Kompenhans et al. / Aerospace Science and Technology ••• (••••) •••–•••

Fig. 2. Left: close-up of the mesh in the vicinity of the left corner. First node in the direction normal to the wall is at 3.5 × 10−5 . Right: incoming X-momentum profile obtained for different turbulence models at x/ L = −4. All profiles give the same momentum thickness at the edge of the cavity.

outside the near wall region. Far from the wall only the numerical viscosity contributes to the physical dissipation, which acts as a subgrid scale model. 3.3. Numerical simulation details The computational domain is divided in three blocks (multiblock mesh topology). The first one discretizing the area upstream the cavity, the second resolves the cavity and the region above it, and a third region to discretise the trailing edge of the cavity (see Fig. 3). Each block is solved independently by means of an ADI implicit solver. The interface connection with the other blocks is performed in an explicit fashion. When using the implicit ADI formulation, it is important to include the sheared layer region (top of the cavity) in a unique block to enable larger time steps. If the large gradients present in this region were to be resolved using an explicit scheme the required time steps (imposed by a CFL restriction) would have need to be of the order of 100 times smaller than selected in our computations. However, let us recall that our implicit time step is fine enough as it enables 300 points per cycle of the fundamental frequency of 2000 Hz. For numerical stability reasons this area must be solved implicitly, because of the high gradients that appear in the shear layer. The mesh has been designed to ensure that the y + values are less than one in the incoming boundary layer. Mesh stretching close to the wall is obtained by imposing a potential law with initial size y 0 = 10−5 . All the turbulence models selected are able to resolve the boundary layer (without needing wall functions) provided that enough grid points are located near the wall. To save computational resources, wall viscous boundary condition are imposed only in walls outside the cavity whilst inviscid boundary conditions are used inside the cavity. In the 3D cases, the sides (front and back) of the open cavity are modeled using periodic boundary conditions. The final mesh has 33 282 nodes in 2D (129 × 129 + 129 × 65 + 129 × 65) and 1 597 536 nodes in 3D (the 2D mesh repeated 48 times in the spanwise direction). The design of the computational mesh is different to the one proposed by Larcheveque et al. [16]. These authors propose a lower definition at the boundary layer with a y +  15 which gives between 5–8 cells inside the boundary layer. Their boundary layer is solved by using a wall function, whilst the mesh is refined to get a

definition of around x+  300 in the streamwise direction. In the present analysis the boundary layer is solved without wall functions, requiring almost twice the computational points close to the wall. To limit the total number of nodes, the streamwise density distribution of nodes is increased such that x+ grows to values close to 1000 upstream of the cavity and around 200 in the shear layer area. Fig. 2(left) shows a detail of the computational mesh. The figure illustrates the mentioned stretching in the streamwise (x-direction) together with the fine wall resolving grid used in the wall normal direction. The incoming boundary layer for the different turbulence models is computed using a separate flat boundary layer computation. In all cases, the incoming profile at the inflow boundary x/ L = −4 (see Fig. 2(right)) is modified to keep the momentum thickness at the cavity edge equal to 6.48 × 10−4 m which is the value measured in the experiments [12]. Computations are performed in two steps. Firstly, a quasi-steady solution is computed and will serve as an “initial solution” for the time-resolved final simulation. The “initial solution” is computed setting the velocities inside the cavity to zero and using the input boundary layer profile (previously computed in auxiliary computations) in the remaining of the computational domain. This “initial condition” is advanced in time using a non-dimensional time step of 0.3, large in absolute terms but small enough to capture the fundamental frequency of f = 2000 Hz observed in the experiments. More precisely, the selected characteristic time step is L /a∞  2 × 10−4 whilst the non-dimensional time corresponding to the fundamental frequency is t  2.5. With these setting, the calculation is advanced for 10 000 time units, which are enough to reach an “initial solution” for the accurate simulation. Once this “initial solution” has been obtained, the simulation is evolved for 1000 non-dimensional time units, with a time step of 0.01. This gives around 300 points per cycle for the fundamental frequency of 2000 Hz. Note that this time step is an order of magnitude higher than the time step used previously [16]. 3.4. Validation with DLR-TAU solver In addition to the main computations presented using the ISAAC solver, we provide results issued form the DLR-TAU solver [27]. This code is commonly used within the European aeronautical industry and has been widely validated. In our calculations, the DLR-TAU code solves the Detached Eddy Simulation (DES)

JID:AESCTE AID:3323 /FLA

[m5G; v1.153; Prn:28/05/2015; 15:49] P.5 (1-11)

M. Kompenhans et al. / Aerospace Science and Technology ••• (••••) •••–•••

5

Fig. 3. Left: computational cavity and measurement points: showing streamtraces and contours of vorticity in the shear layer. Right: implicit LES 3D computation, vortical structures depicted using iso-surface of vorticity and slice showing Shadowgraph type contours. (For interpretation of the colors in this figure, the reader is referred to the web version of this article.)

Fig. 4. Mean horizontal velocity profile at three locations x/ L = 0.2, 0.4 and 0.8.

equations on unstructured hybrid grids employing a second-order finite-volume discretization. Convective terms are discretized with central fluxes and a scalar Jameson–Schmidt–Turkel (JST) numerical dissipation model. For time-accurate computations, the solver uses a dual time stepping methodology, where a series of fictitious pseudo-time steady problems are integrated at each time step. A backward implicit scheme is used in the real time whilst Runge–Kutta and Multigrid acceleration are utilized in the dual time steady problem. The Detached Eddy Simulation [29] approach implemented in the DLR-TAU solver, blends the one equation Spalart–Allmaras RANS model when solving boundary layers and the Large Eddy Simulation approach away from walls. We perform 2D and 3D simulations. The 2D mesh has 25 150 elements and is extruded in the third dimension to obtain a 3D mesh of size 1 207 200 elements. In the next section we compare Implicit LES results using ISAAC to DES DLR-TAU simulations. 4. Results In this section we present the main results that will be used for the discussion. We first include the mean velocity profiles at various points located at the cavity shear layer (top of the cavity). Secondly, we compare the frequency spectra for various calculations: implicit LES, DES using DLR-TAU and various eddy viscosity models.

4.1. Mean velocity gradients Time history of flow quantities and associated power spectra are computed at a particular point located at the upstreamed wall with coordinates (x/ L , z/ L ) = (0.0, −0.7) and is depicted as a red dot in Fig. 3(left). In addition, we compare numerical and experimental time averaged velocity profiles at three shear layer location (x/ L , z/ L ) = (0.2, 0.0), (x/ L , z/ L ) = (0.4, 0.0) and (x/ L , z/ L ) = (0.8, 0.0), green, blue and orange points in Fig. 3(left). In addition and for illustrative purposes, we show in Fig. 3(right) a snapshot of the 3D Implicit LES computation, where vortical structures and acoustic waves are shown using vorticity iso-surfaces and Shadowgraph contours. We start by presenting a comparison, for the various simulations, of the averaged streamwise velocity in the three location near the sheared layer region. Note that the average velocity is obtained in time and space (using the average of various spanwise sections). The resulting sheared profiles are shown in Fig. 4. The figure shows that the results provided by the 3D ILES simulation agree well with the experimental data. In addition, it can be seen that the 2D ILES provides a similar profile but over-estimates the velocity inside the upper limit of the cavity, leading to the wrong sheared profile. Finally, all URANS simulations (that use eddy viscosity assumptions) show similar behavior, and although a good agreement with experiments in the upper shear layer region is observed (outside the cavity), the sheared profile is not well captured due to the over-predicted gradients near the limit

JID:AESCTE

AID:3323 /FLA

6

[m5G; v1.153; Prn:28/05/2015; 15:49] P.6 (1-11)

M. Kompenhans et al. / Aerospace Science and Technology ••• (••••) •••–•••

Fig. 5. Comparison of LES and DES solvers. Power Spectral Density (PSD) of the time evolving density at monitoring point z/ L = −0.7, Left: 3D simulations and Right: 2D simulations.

of the cavity. The recovery of the velocity profiles (as we progress in the streamwise x-direction) for the ILES computation agrees better with the experimental results than when using RANS models, showing that these models are not adequate to accurately capture sheared layer flows and Kelvin–Helmholtz instabilities. Finally, 2D and 3D ILES simulations show very similar profiles at x/ L = 0.8 even when a different behavior is observed near the cavity leading edge (x/ L = 0.8). 4.2. Pressure spectra at the upstream wall We first show the power spectra (using a long time histories) for the pressure signal monitored at the upstream wall of the cavity at z/ L = −0.7 (red point in Fig. 3(left)). In Fig. 5 we compare frequencies issued from ILES simulations to DES DLR-TAU computations. All spectra have been smoothed using Welch method [34], which is often used when post processing numerical and experimental data [6,23]. 3D results are shown in Fig. 5(left) and 2D in Fig. 5(right). In addition, we include published LES results from Larcheveque et al. [16] and experimental data from Forestier et al. [12]. It can be seen that both 3D and 2D ILES computations capture the two most important peaks at 2000 and 4000 Hz and also the smaller peak at 3500 Hz. However, DES computations although predicting accurately the main peaks, fail in capturing the mid frequency at 3500 Hz, which is associated to a longitudinal mode in the cavity [16]. Let us note that neither Larcheveque et al. in the numerical work [16], nor Forestier et al. in the experimental report [12], detail the filtering (or smoothing) procedure used to produce the data shown in the papers. However, we have shown that good agreement is obtained when comparing our results (using Welch method) with previous published data. We now detail the results for the ILES calculations and depict in Fig. 6(right) a small portion of the time history for the density. In addition, the power spectra (using a long time history) for the pressure signal is presented in Fig. 6(left). The top figures show results for 2D and 3D Implicit LES results, that compare favorably to experimental data. The bottom figures compare the various URANS turbulent approaches, showing discrepancies with experimental data. In particular, it becomes apparent that URANS k–epsilon solutions over-damp high frequency fluctuations

leading to a degrease in the background PSD level after the first peak. However, the k–eps anisotropic model (k–eps ASM) shows an accurate background PSD level for all frequencies. Despite this better behavior, this enhance k–epsilon model fails in predicting the main frequency peaks. All figures illustrate the dependency of the spectra (amplitude and frequency) with the turbulence model selected for the simulation. Indeed, it will be shown that only 3D computations provide accurate solutions for the two fundamental frequencies. In addition to the computations that compare turbulent closures, we also quantify the effect of varying the Mach number in various Implicit LES simulations. A summary for the first two frequency peaks of the power spectra density of the pressure is provided in Table 1. The table includes the frequencies determined by the empirical correlation for the Rossiter feedback mechanism [24]:

freqRoss =

u ∞ n −α L 1 +M κ

n = 1, 2, . . . ,

(4)

where α = 0 for this particular aspect ratio [16], κ = 0.57 (as shown in [24]) and M = 0.8 is the Mach number. In addition, we include the frequency for the fundamental resonant mode as defined in [9,38]

freqres =

u∞ 0.25 L M {1+0.65( L / D )0.75 }

(5)

with L / D = 0.42 as defined previously. As shown in the table, the frequency for the first Rossiter mode (freqRoss = 2020 Hz) is remarkably close to the fundamental resonant mode (freqres = 2159 Hz) and hence only one unique peak appears in the spectra (which is a combination of the two phenomenon). Note that these empirical formulas only provide the frequency for the phenomenon and do not provide pressure levels. The following section will discuss the presented results and will attempt to clarify some of the observed results. 5. Discussion and analysis 5.1. Turbulence modeling and 2D/3D comparisons at Mach M = 0.8 A detail of the two first peaks in the range of 1500 to 4500 Hz is shown in Fig. 6. The fundamental frequency of 1975 Hz, the first

JID:AESCTE AID:3323 /FLA

[m5G; v1.153; Prn:28/05/2015; 15:49] P.7 (1-11)

M. Kompenhans et al. / Aerospace Science and Technology ••• (••••) •••–•••

7

Fig. 6. Comparisons using ISAAC solver. Left: power Spectral Density of the time evolving density at monitoring point z/ L = −0.7. Right: partial plot of the unsteady density as a function of the nondimensional time at z/ L = −0.7. Top: Implicit LES calculations. Bottom: URANS models.

harmonic at 3950 Hz and the small peak at approximately 3500 Hz are clearly captured by the 3D ILES computations and reasonably well by the 2D ILES. Nevertheless, the 3D predicted decibels are around 10–15 dB higher than experimentally predicted. The reasons for this over-prediction may be related to the inviscid wall used as boundary conditions inside the cavity. This non-viscous boundary can reduce the damping in the frequency peak and could also be responsible of the increase of the background level in the 3D simulation. The 2D ILES computations provides accurate results, the main peaks are well observed in the spectra and the 3500 Hz frequency appears slightly shifted with respect to the 3D case. Fig. 6 shows that the main frequency (i.e. with highest PSD level) is in the 2D case around 4000 Hz and not at 2000 Hz as in the 3D case and the experiments. Let us finally note that the background level is always lower in the 2D case. We have not found any explanation to the latter observation and further analysis for the differences of between 2D and 3D are required. However it is clear that the

2D ILES simulation provides overall a good estimation of the tonal frequencies in this challenging case. Regarding the unsteady RANS simulations, neither the 3D nor 2D computations provide sufficient accuracy as they are unable to predict the experimental tonal frequencies. As shown in Fig. 6(right), the URANS simulations provide smoother unsteady (more damped) solutions than the corresponding ILES runs. This is not surprising since the eddy viscosity assumptions used in RANS modeling are often over dissipative, resulting in very low values of decibels of the pressure spectra (damped time histories). To further illustrate this behavior, we provide the phase portrait of the velocity vector monitored at z/ L = −0.7. Fig. 7(left) compares ILES and k–eps models. Since flow velocities at the monitoring point are constrained in the x-direction by the presence of the wall, we find that the only relevant and quantitatively important velocity fluctuations occur in the depicted y–z plane. Note that 2D simulation do not present fluctuations in the y-direction (spanwise direction) and hence appear as a line in the figure. Fig. 7(left)

JID:AESCTE

AID:3323 /FLA

[m5G; v1.153; Prn:28/05/2015; 15:49] P.8 (1-11)

M. Kompenhans et al. / Aerospace Science and Technology ••• (••••) •••–•••

8

Table 1 Frequency spectra: only the two larger peaks (higher pressure levels) are shown for different turbulent models and Mach numbers. 2D = two dimensional computations. 3D = three dimensional computations. ILES = Implicit LES, k–eps = k–epsilon turbulence model and k–eps-ASM = k–epsilon Anisotropic Stress Model. Method

Mach

Dimension

1st peak pressure

2nd peak pressure

Hz

dB

Hz

dB

Experimental [12] Rossiter (freqRoss ) [24] Resonance (freqres ) [9,38] ILES ILES k–eps k–eps k–eps-ASM

0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8

3D 2D 2D 3D 2D 3D 2D 2D

1975 2020 2159 1971 1974 1793 1830 1849

155 – – 152 136 160 150 150

3950 4040 4318 3952 3969 3587 3647 3326

153 – – 142 152 135 133 144

DES DLR-TAU DES DLR-TAU

0.8 0.8

3D 2D

2050 2050

157 157

4102 4101

155 155

ILES ILES

0.6 0.6

3D 2D

340 328

160 165

704 489

159 149

ILES ILES

0.4 0.4

3D 2D

431 431

169 167

869 863

154 159

ILES ILES

0.2 0.2

3D 2D

424 419

161 154

467 470

167 170

Fig. 7. Left: velocity phase portrait comparison of various turbulence models. Right: velocity phase portrait for ILES when varying the Mach number. Monitoring point is z/ L = −0.7.

shows that URANS modeling provide limited velocity fluctuations when compared to Implicit LES methods. The k–eps models shows a well organized periodic path. The velocity path is constrained by one unique physical attractor but it is not clear that this attractor is located where the real flow solution lives. When observing the velocity profiles predicted by the ILES models, it can be seen a more disorganized and less periodic behavior, showing the existence of various attractors that interplay. Fig. 7(left) depicts the phase portrait of the velocities for various ILES computations, where the Mach number has been varied. It can be seen that the higher the Mach number, the more spatially constrained appears the velocity path. This is in good agreement with previous observations that have shown that K–H instabilities are attenuated as the compressibility effects increase [18,19, 38]. Note also that in the context of flow stability [38], the effect of compressibility has shown to decrease Rossiter type of interactions in open cavities at low Reynolds numbers (laminar flow regimes). Having seen that the ILES computations provide the most accurate results, we now try to shed light on the differences observed between the 2D and 3D simulations. Fig. 6(left) showed that the pressure spectra for 2D and 3D computations are similar

but that some differences could be observed. The 2D simulation has a dominant fundamental frequency larger that in the 3D case (3969 Hz with 152 dB), which agrees with the time pressure evolution shown in Fig. 6(right). Additionally, the 2D spectra presents a peak at 2000 Hz (Fig. 6(left)). A possible explanation is provided by Fig. 8, where a snapshot of the vorticity (a) and normal velocity (b) are presented. It can be observed that the 3D computations produce less but larger vortices in the shear layer. These vortices move slower and grow longer provoking more energetic structures which impact against the edge of the cavity. The impact of these structures induces energetic acoustic waves which travel upstream and close Rossiter’s loop mechanism. This wave that travels upstream excites a resonance inside the cavity. These two effects can be observed in Fig. 9, where a comparison of the magnitude of the density gradient ||∇ ρ || for the 3D (left) and 2D (right) simulations shows the existence of these waves for the 3D cases. However, it can be seen that in the 2D case only the wave outside the cavity is observed. This difference explains the existence of the peaks (the K–H frequencies are exited) and the disparity in the absolute magnitude of these peaks since the wave energy differs between 2D and 3D simulations. Finally, these snap-

JID:AESCTE AID:3323 /FLA

[m5G; v1.153; Prn:28/05/2015; 15:49] P.9 (1-11)

M. Kompenhans et al. / Aerospace Science and Technology ••• (••••) •••–•••

9

Fig. 8. Snapshot of ILES simulations. (a) vorticity map and (b) normal velocity contour. Comparison between the 3D (left) and 2D (right) computations at M = 0.8.

is substantially different in 2D and 3D simulations, which has a noticeable effect in the frequency spectra. 5.2. Mach number effects

Fig. 9. Snapshot of ILES simulations displaying ||∇ ρ ||. Comparison between the 3D (left) and 2D (right) computations at M = 0.8.

shots suggest that the jitter phenomenon [13] (i.e. intermittency between partial clipping and partial scape of the impinging vortices)

Let us examine the role of the Mach number and its relation to 2D/3D dimensionality of the simulation for Implicit LES methods. Table 1 showed that the first two main frequencies do not coincide when varying the Mach number. Indeed, the first and second peaks for Mach numbers 0.2 to 0.6 appear at lower frequencies than the fundamental frequency of 2000 Hz. If the Mach number used for simulation is 0.8 then we retrieve the experimental frequencies. Let us now study the effect of compressibility together with the simulation dimensionality (2D vs 3D ILES). Fig. 10 shows frequency results at the monitoring point z/ L = −0.7 for 2D and 3D simulations and two Mach numbers (M = 0.4 and 0.6). At these Mach numbers, the frequencies are lower than in the results presented in Fig. 6 for M = 0.8. It is remarkable the similarity between the spectra obtained at Mach M = 0.4 for 2D and 3D computations, Fig. 10(left), where a very accused fundamental frequency of 430 Hz and its harmonics appear. Fig. 9(right) shows that this

Fig. 10. Comparison of the PSD of 3D (left) and 2D (right) computations at two Mach numbers 0.4 (left) and 0.6 (right). Experimental data is for Mach 0.8.

JID:AESCTE

AID:3323 /FLA

[m5G; v1.153; Prn:28/05/2015; 15:49] P.10 (1-11)

M. Kompenhans et al. / Aerospace Science and Technology ••• (••••) •••–•••

10

∂p

Fig. 11. Plan view of the spanwise pressure gradient ( ∂ y ) at section z/ L = 0.1. Top: M = 0.8 and Bottom: M = 0.4.

agreement is lost when compressibility effects become important (e.g. M = 0.6). In the latter case, the 3D solution behaves differently than the 2D version. It also interesting to note that the background noise level are higher for in 3D computations when compared to the 2D results for a fixed Mach number. This behavior can also be observed in the computations that use the original Mach number M = 0.8, in Fig. 6. The importance of 3D effects can also be seen in Fig. 11, where plane cuts at z/ L = 0.1 for the spanwise gradient of the pressure ∂p ∂ y are shown for M = 0.4 (bottom) and M = 0.8 (top). These figures confirm that 3D effects at M = 0.4 are negligible when compared to results at M = 0.8. The importance of 3D structures as a function of the Mach number has been previously reported [3] for lower Reynolds numbers (i.e. no turbulence modeling required) using linear stability analysis. Our analysis extends these results to high Reynolds numbers and, in addition, quantifies the effects of compressibility and dimensionality in the frequency spectra. 6. Conclusions In this work, the capabilities of various turbulence models to predict the dynamical content of open cavity flows have been studied and compared. Unsteady simulations have been performed using four turbulence models: k–epsilon, k–epsilon with anisotropic Reynolds stresses, a detached eddy simulation technique and an implicit LES approach. A cavity of length to depth ratio of L / D = 0.24, at Reynolds of 8.6 × 105 and various Mach numbers, has been used in the calculations. In this analysis, we extended previous studies and observations by considering different Mach numbers and 2D/3D effects. The results have shown that only Implicit LES simulations provide good comparisons, both qualitatively and quantitatively, with experimental data. The unsteady turbulence URANS models have shown to provide inaccurate results, regardless of the isotropic or anisotropic treatment of the explicit turbulent viscosity. 2D ILES computations show qualitatively good agreement with experiments and main frequencies, but fail in predicting the absolute values for the power spectra. The general flow features appear to be strongly two-dimensional but our results suggest that the three dimensional structures have a strong influence in the flow behavior. Finally, when studying compressibility effects, it is observed that there is a limiting Mach number (M ≈ 0.6) below which three dimensional effects are unimportant. For Mach numbers above this limit 3D simulations become compulsory. Conflict of interest statement None declared.

Acknowledgement Support of the Marie Curie Grant PIRSES-GA-2009-247651 FP7-PEOPLE IRSES: ICOMASEF Instability and Control of Massively Separated Flows is gratefully acknowledged. In addition, MK, EF and EV would like to thank the European Commission for the financial support of the ANADE project (Advances in Numerical and Analytical tools for DEtached flow prediction) under grant contract PITN-GA-289428. References [1] R. Abid, J.H. Morrison, T.B. Gatski, C.G. Speziale, Prediction of complex aerodynamic flows with explicit algebraic stress models, AIAA Paper 96-0565, 1996. [2] J.P. Boris, F.F. Grinstein, E.S. Oran, R.L. Kolbe, New insights into large eddy simulation, Fluid Dyn. Res. 10 (4–6) (1992) 199. [3] G.A. Bres, T. Colonius, Three-dimensional instabilities in compressible flow over open cavities, J. Fluid Mech. 509 (2009) 309–339. [4] L. Cattafesta, F. Alvi, D.R. Williams, C.W. Rowley, Review of active control of flow-induced cavity oscillations, in: 33rd AIAA Fluid Dynamics Conference and Exhibit, 2003. [5] C.F. Chenault, P.S. Beran, R.D.W. Bowersox, Numerical investigation of supersonic injection using a Reynolds-stress turbulence model, AIAA Pap. 37 (10) (1999). [6] S.L. Clainche, J. Li, V. Theofilis, J. Soria, Flow around a hemisphere-cylinder at high angle of attack and low Reynolds number. Part I: Experimental and numerical investigation, Aerosp. Sci. Technol. (2014), http://dx.doi.org/10.1016/ j.ast.2014.03.017. [7] T. Colonius, S. Lele, Computational aeroacoustics: progress in nonlinear problems of sound generation, Prog. Aerosp. Sci. 40 (6) (2004) 345–416. [8] T. Colonius, S. Lele, P. Moin, Sound generated in a mixing layer, J. Fluid Mech. 330 (1997) 375–409. [9] L. East, Aerodynamically induced resonance in rectangular cavities, J. Sound Vib. 3 (3) (1966) 277–287. [10] E. Ferrer, A high order discontinuous Galerkin–Fourier incompressible 3d Navier–Stokes solver with rotating sliding meshes for simulating cross-flow turbines, Ph.D. thesis, University of Oxford, 2012. [11] J. Ferziger, M. Peric, Computational Methods for Fluid Dynamics, Springer, 1997. [12] N. Forestier, L. Jacquin, P. Geffroy, The mixing layer over a deep cavity at highsubsonic speed, J. Fluid Mech. 471 (2003) 101–145. [13] X. Gloerfelt, Cavity Noise, in: VKI Lectures: Aerodynamic Noise from WallBounded Flows, von Karman Institute, 2009. [14] M.E. Goldstein, Aeroacoustics, McGraw–Hill International Book Co., New York, 1976. [15] F. Grinstein, L. Margolin, W. Rider, Implicit Large Eddy Simulation: Computing Turbulent Fluid Dynamics, Cambridge Univ. Press, Leiden, 2007. [16] L. Larcheveque, P. Sagaut, I. Mary, O. Labbe, Large-eddy simulation of a compressible flow past a deep cavity, Phys. Fluids 15 (1) (2003) 193–210. [17] M.J. Lighthill, On sound generated aerodynamically. I: General theory, Proc. R. Soc. 221 (1952) 564. [18] P. Meliga, D. Sipp, J. Chomaz, Absolute instability in axisymmetric wakes: compressible and density variation effects, J. Fluid Mech. 600 (2008) 373–401. [19] P. Meliga, D. Sipp, J. Chomaz, Effect of compressibility on the global stability of axisymmetric wake flows, J. Fluid Mech. 660 (2010) 499–526. [20] C. Mettot, F. Renac, D. Sipp, Computation of eigenvalue sensitivity to base flow modifications in a discrete framework: application to open-loop control, J. Comput. Phys. 269 (0) (2014) 234–258. [21] J.H. Morrison, A compressible Navier–Stokes solver with two-equations and Reynolds stress turbulence closure models, NASA CR-4440, 1992.

JID:AESCTE AID:3323 /FLA

[m5G; v1.153; Prn:28/05/2015; 15:49] P.11 (1-11)

M. Kompenhans et al. / Aerospace Science and Technology ••• (••••) •••–•••

[22] S. Pirozzoli, Numerical methods for high-speed flows, Annu. Rev. Fluid Mech. 43 (1) (2011) 163–194. [23] W. Press, S. Teukolsky, W. Vetterling, B. Flannery, Numerical Recipes in C, 2nd edn., The Art of Scientific Computing, Cambridge University Press, New York, NY, USA, 1992. [24] J. Rossiter, Wind tunnel experiments on the flow over rectangular cavities at subsonic and transonic speeds, Tech. rep., Royal Aircraft Establishment RAE, Ministry of Aviation, Farnborough, 1964. [25] P. Sagaut, Large Eddy Simulation for Incompressible Flows: An Introduction, Scientific Computation, Springer, 2001. [26] P. Sagaut, S. Deck, M. Terracol, Multiscale and Multiresolution Approaches in Turbulence, Imperial College, 2006. [27] D. Schwamborn, T. Gerhold, R. Heinrich, The DLR tau-code: recent applications in research and industry, in: E.O.P. Wesseling, J. Priaux (Eds.), ECCOMAS CFD European Conference on Computational Fluid Dynamics, Delft University of Technology, The Netherlands, 2006. [28] J. Smagorinsky, General circulation experiments with the primitive equations, Mon. Weather Rev. 91 (3) (1963) 99–164. [29] P.R. Spalart, W. Jou, M. Strelets, S.R. Allmaras, Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach, in: C. Liu, Z. Liu (Eds.), First AFOSR International Conference on DNS/LES, in: Advances in DNS/LES, Greyden, Columbus, OH, Ruston, LA, 1997. [30] B. Thornber, D. Drikakis, Implicit large-eddy simulation of a deep cavity using high-resolution methods, AIAA J. 46 (10) (2008).

11

[31] A. Uranga, P. Persson, M. Drela, J. Peraire, Implicit large eddy simulation of transition to turbulence at low Reynolds numbers using a discontinuous Galerkin method, Int. J. Numer. Methods Eng. 87 (1–5) (2011). [32] H. Versteeg, W. Malalasekera, An Introduction to Computational Fluid Dynamics: The Finite Volume Method, Longman Scientific & Technical, 1995. [33] M. Wang, J. Freund, S. Lele, Computational prediction of flow-generated sound, Annu. Rev. Fluid Mech. 38 (1) (2006) 483–512. [34] P. Welch, The use of fast Fourier transform for the estimation of power spectra: a method based on time averaging over short, modified periodograms, IEEE Trans. Audio Electroacoust. 15 (1967) 70–73. [35] D. Wilcox, Turbulence Modeling for CFD, DCW Industries Inc., La Cañada, CA, 1993. [36] J.E.F. Williams, D.L. Hawkings, Sound generation by turbulence and surfaces in arbitrary motion, Philos. Trans. R. Soc. A 264 (1151) (1969) 321–342. [37] S.L. Woodruff, J.M. Seiner, M.Y. Hussaini, G. Erlebacher, Evaluation of turbulence-model performance in jet flows, AIAA J. 39 (12) (2001). [38] S. Yamouni, D. Sipp, L. Jacquin, Interaction between feedback aeroacoustic and acoustic resonance mechanisms in a cavity flow: a global stability analysis, J. Fluid Mech. 717 (2013) 134–165. [39] H.S. Zhang, R.M.C. So, T.B. Gatski, C.G. Speziale, A near wall second order closure for compressible turbulent flows, in: Near-Wall Turbulent Flows, Elsevier Science Publishers, 1993.