Prediction of the nonlinear response of an acoustically forced flow withlinearized Navier-Stokes equations
Abel FAURE BEAULIEU(1)† , Claire BOURQUARD(2) , Nicolas NOIRAY(3)(1) CAPS Lab, Mechanical and Process Engineering Department, ETH Zürich, Zürich 8092, Switzerland, abelf@ethz.ch(2) CAPS Lab, Mechanical and Process Engineering Department, ETH Zürich, Zürich 8092, Switzerland, clairebo@ethz.ch(3) CAPS Lab, Mechanical and Process Engineering Department, ETH Zürich, Zürich 8092, Switzerland, noirayn@ethz.ch

Abstract
The prediction of self-sustained aeroacoustic oscillations has been the topic of numerous of theoretical studies,but obtaining quantitative predictions on practical configurations is still challenging. In this study we use particle
image velocimetry to record flow data at an acoustically excited T-Junction subject to low-Mach grazing flow.
Data is recorded for several acoustic excitation amplitudes and frequencies. This allows to formulate a model
predicting the non-linear response of the system with Navier-Stokes equations linearized around different mean
flows measured with different excitation amplitudes and frequencies. This linearized perturbation approachgives a good prediction of the spatial repartition of the velocity fluctuations and predicts their amplitude to areasonable degree. The test rig serves as a surrogate for the aeroacoustic behavior of deep closed cavities.
Keywords: Aeroacoustics, Linearized Navier-Stokes equations

1

INTRODUCTION

Certain geometrical configurations placed in a moving flow can lead to self-sustained pressure oscillations due
to a constructive feedback loop between the unsteady motion of the flow and the acoustics. This leads to aloud tonal sound commonly called whistling. A typical example of configuration which can produce whistlingis the rectangular closed cavity subject to a grazing flow. The velocity of the air in the cavity is low comparedto the velocity of the grazing flow, so that a shear layer is present at the mouth of the cavity. If the velocityof the flow is large enough, this shear layer is unsteady and vortical structures are convected in it, as well inlaminar regime as in turbulent regime. Since rectangular side cavities are present in many industrial domains(e.g aircraft wheels wells or bomb bays, pantographs cavities for trains, side branches in pipe networks), a
lot of literature already exists on the topic. The effect of the cavity aspect ratio and the grazing flow’s Mach
number was investigated through the early wind tunnel experiments of Krishnamurthy (1) and Plumblee et al.
(2). (3, 4, 5) propose models to predict the aeroacoustic resonance frequencies and the conditions under which
they happen. In particular, East (4) studies specifically the aeroacoustic resonance of deep cavities, whose depth
D is greater than the length L of the cavity’s mouth in the flow direction. He shows the resonance can betriggered when the characteristic convective time of the vortical disturbances in the shear layer across the cavity
width matches the acoustic frequency of a depth acoustic mode of the cavity. Under these conditions, twomechanisms are interacting. The first one is the excitation of a cavity acoustic mode by the unsteady vorticity
in the shear layer, which is referred to as ’vortex sound mechanism’. The second mechanism is the forcingof the shear layer by the acoustic oscillations: The shear layer being located at an acoustic velocity antinode,
it undergoes vertical acoustic velocity fluctuations triggering the periodic formation of vortices. When the two
mechanisms are constructively interacting, it can lead to large amplitudes oscillations. Nonlinearities come theninto play and lead to a saturation of the response (6, 7). Predicting the final amplitude of the self-sustained
oscillations is still a challenge; It requires quantitative models both for the vortex sound mechanism (see theworks of Howe (8, 9)) and for the nonlinear response of the flow to acoustic oscillations. The linearized NavierStokes equations (LNSE) around the time-averaged unsteady flow appears to be an interesting tool to predict
the flow’s response to an external forcing: the linearized Navier-Stokes operator (LNS) can be used to describe5299

the evolution of small perturbations around the mean flow. Barkley (10) used it on a unstable laminar cylinder
wake and managed to predict the frequency of the vortex shedding in the cylinder’s wake at different Reynolds
number only knowing the mean flow in each case. In (11), the studied flow is a linearly stable laminar jet
and the LNSE allows to predict the response of the flow to an external harmonic forcing. In (7), the studied
configuration is a deep rectangular cavity whose top wall is considered as anechoic, which is different fromthe closed cavities studied by East (4) since the absence of reflection on the top wall opens the aeroacoustic
feedback loop: the sound emitted by the vorticity is evacuated through the anechoic wall instead of coming backto the shear layer. The configuration of (7) allows therefore to study specifically how the shear layer responds
to an acoustic excitation. LES simulations were conducted on turbulent flows undergoing different amplitudes
of forcing. For each simulation, the time-averaged flow is taken and the LNS operator is computed. Although
only linearized perturbations are considered, the nonlinearities affect the structure of the mean flows obtainedfor different forcing amplitudes. The hydrodynamic responses predicted by the LNSE study are shown to be inagreement with the ones observed in the LES simulations and the nonlinear saturation is well predicted. The
present paper shows an experimental study on a cavity similar to the one studied in (7). The incoming flow ishighly turbulent, and the instantaneous flow is resolved using fast particle image velocimetry (PIV). In 2, theexperimental setup and the operating conditions are presented. In 3, the LNSE around the forced mean flow arederived. In 4, the coherent velocity fluctuations extracted from the PIV measurements are compared with theones predicted by the LNSE.

2EXPERIMENTAL CONFIGURATION

2.1 Geometry and experimental conditions
The experimental configuration used for this study is a channel of square cross-section with a perpendicularside branch of rectangular cross-section. A compressor sends a turbulent flow of mean velocity U at low Mach
number in the main channel. The x axis is defined as the direction of the main channel and is oriented so thatthe mean flow’s velocity is positive. The y axis is oriented towards the side branch and the z axis corresponds
to the spanwise direction. The edge of the main channel’s section is D=6.2 cm, the width of the side branchin the x direction is W=3 cm. Figure (1) shows the main features of the experimental setup. An unsteadyshear layer is present at the junction. To prevent acoustic reflections, the end of the side branch is terminatedwith a horn opening to a box lined with an acoustic absorbing foam. Both of the main channel’s ends are also
horn-shaped anechoic ends. The shear layer is acoustically forced with two loudspeakers placed on the side
branch. Several microphones are placed on each branch to reconstruct spatially the acoustic field. To get thevelocity field in real time, particle image velocimetry (PIV) is used. To do so, a laser sheet is used to enlighten
the particles at the middle of the main channel’s span. This means that only a 2D slice of the velocity field canbe measured. In this study, 3D effects due to the finite spanwise extension of the channel are not considered.
The experimental setup allows to impose different flow velocities, acoustic forcing levels and frequencies. The
range of investigated mean flow velocities goes between U=41 m/s and U=76 m/s. We define a Reynolds
number ReD = UD/ν based on the channel’s hydraulic diameter, where ν is the kinematic viscosity of air at
30◦ C. ReD is between 1.5 · 105 and 2.7 · 105 . The incoming pipe length between the horn contraction and the
T-junction is about 1 m, which is not sufficient to get a fully developed velocity profile (12) but enough to getstrong turbulent fluctuations in the whole channel’s cross section. Each measurement gives a set of 2000 images
acquired at a frequency of 5kHz. For each measurement, a mean flow is computed by taking the average of all
the images.
2.2 Compacity assumption
Looking at the experimental results, the zone reacting the most to the acoustic forcing is the shear layer, which
remains confined in a horizontal band of 1 cm of thickness. The speed of sound in the experimental conditions
is c = 349 m/s. The tested frequencies for the acoustic forcing are between 800 Hz and 1800 Hz. For thehighest frequency, the wavelength is 0.2 m, which is great compared to the shear layer’s dimensions. We can

5300

Figure 1. Simplified drawing of the experimental test rig. Left: whole setup. The horn-shaped expansions areused to prevent acoustic reflections. Right: zoom on the T-junction. W = 3 cm, D = 6.2 cm. The position of
the shear layer is shown in green. In blue, the zone observed by the PIV camera.
then assume that the shear layer is acoustically compact, so that compressibility effects can be neglected.

3

LINEARIZED NAVIER-STOKES EQUATIONS

3.1 Linearization
The acoustic forcing is modelled as a horizontal oscillating wall Γa in the cavity imposing harmonic oscillationsin the vertical direction : fa = ua cos(ωat)ex . We start from the incompressible Navier-Stokes equations for a
newtonian fluid of kinematic viscosity ν and density ρ:∇·u = 0∂u1
+ u · ∇u = − ∇p + ν∇2 u
∂t
ρ

(1)

with the boundary conditions:
u=0

on Γw

(walls)u = faonΓawhere u is the velocity field and p the pressure field. As in (13), the fluctuating variables are decomposed as
˜ corresponding to the response to the harmonic
the sum of a time-averaged value (·), a coherent fluctuation (·)
forcing, and a non-coherent fluctuation (·)0 coming from the turbulent noise. We therefore have u(x,t) =u(x) + ũ(x,t) + u0 (x,t). The time averaging of eq. (1) gives an equation for the mean flow:∇·u = 0
1u · ∇u = − ∇p + ν∇2 u − ∇ · ũũ + u0 u0ρ(2)with the boundary conditions u = 0 on Γa ∪ Γw . The last term of the equation shows that the velocity fluctuations have an impact on the mean flow. Taking the phase-average of eq. (1) and subtracting eq. (2), the5301

following equations are obtained for the coherent fluctuations:∇ · ũ = 0


1
∂ ũ
0 u0
+ u · ∇ũ + ũ · ∇u = − ∇ p̃ + ν∇2 ũ − ∇ · f̃uũ + ug∂t
ρ(3)with the boundary conditions ũ = fa on Γa and ũ = 0 on Γw .
0 u0 needs a model. Following (14) and (7), we use the
The coherent part of the turbulent Reynolds tensor ∇ · ug
Boussinesq assumption and we introduce a turbulent viscosity νt so that:2−u0 u0 + k0 = νt (∇ + ∇T )u3g−u0 u0 = νt (∇ + ∇T )ũ

(4)
(5)

where k is the turbulent kinetic energy. The Boussinesq assumption is not valid in the general case, but was
found to be reasonable for simple shear flows (15). Using eq. (4), the turbulent viscosity field νt (x) is obtainedby extracting the random fluctuations from the PIV measurements.
Injecting eq. (5) in eq. (3), the turbulent viscosity and the molecular kinematic viscosity can be combined inan effective viscosity νeff = νt + ν. Eq. (3) becomes:∂ ũ1+ u · ∇ũ + ũ · ∇u = − ∇ p̃ + νeff ∇2 ũ + (∇u)T ∇νt − ∇ · f̃
uũ
∂t
ρ

(6)

In the experimentals results, the term (∇u)T ∇νt was found to be small compared to νeff ∇2 ũ and will therefore
be omitted in the following.
3.2 Equation in the frequency domain
Assuming that the coherent response contains only the forcing frequency ωa and its harmonics nωa , n ∈ N∗ , ũwrites as:ũ(x,t) = ∑ ûn (x)einωa t

(7)

n6=0

where û−n = û∗n . The Fourier component of the equation (6) at the forcing frequency is:1iωa û1 + u · ∇û1 + û1 · ∇u = − ∇ p̂1 + νeff ∇2 û1 − ∑ ∇ · ûn û1−nρn6=0,n6=1(8)with the boundary conditions û1 = ua on Γa and û1 = 0 on Γw . If we assume that the response at the higher
harmonics is negligible compared to the response at the forcing frequency, this reduces to:1iωa û1 + u · ∇û1 + û1 · ∇u = − ∇ p̂1 + νeff ∇2 û1ρ(9)

Now that we consider only the response at the fundamental frequency, the indices (·)1 will be dropped.
We define the state vector q = {û, p̂}. Eq. (9) can be written as:
R(ωa , u)q ≡ iωa Iq + L(u)q = P(fa )(10)1∇ p̂ − νeff ∇2 û. I is the identity operator and P isρa projection operator from the forced boundary Γa to the whole domain and its boundaries. R = iωa I + L(u) iscalled the resolvent operator.
where L(u) is the LNS operator: L(u)q = u · ∇û + û · ∇u +

5302

(a)

(b) ux(c) uyFigure 2. a): 2D domain used for the discretization. Γa : , Γw : rigid walls, Γ f : stress-free boundaries. b,c): x
and y components of the mean flow velocity for an incoming mean velocity U=76 m/s and an acoustic forcing
imposing harmonic acoustic oscillations of amplitude 0.5m/s at 1000Hz in the shear layer.
3.3 Numerical method
For each experiment, the 2D fields u and νt are extracted. The equivalent variational forms of the linear
operators R(ωa , u) and P are discretized on a 2D grid, using finite elements. The degrees of freedom for eachvelocity component are represented by P2 finite elements, the pressure and turbulent viscosity, by P1 elements.
The 2D domain has the same shape as the PIV window and is shown in Fig. (2a). For a given forcing fa , thelinearized perturbations are obtained by inverting the finite element matrix of the resolvent operator. The 2D
mesh has 79508 triangular cells. Computing the perturbations for a given forcing on a given mean flow can beachieved in less than one minute on a desktop computer.

4

COMPARISON BETWEEN EXPERIMENTAL RESULTS AND LNSE PREDICTIONS

4.1 Harmonic response
Experiments were conducted on a flow of mean velocity U=76 m/s. Different acoustic frequencies and amplitudes are applied. The harmonic component of the velocity fluctuations at the excitation frequency is extracted
for each experiment, by computing the Fourier transform of the velocity at each spatial location. At the sametime, the LNSE is used to predict the response of the mean flow to a forcing fa of the same acoustic amplitudeand frequency as in the experiment. The experimental value of fa is based on the acoustic velocity at the shear
layer, obtained with a multi-microphone method (16). Figure 3 shows the spatial distribution of the velocity
fluctuations from the LNSE and the experiment for a given forcing amplitude and different frequencies. The
spatial shapes are well reproduced, but the LNSE overestimates the amplitudes. Figure 4 shows the |ûy | profiles
along the shear layer. For forcing frequencies around 1400 Hz, the LNSE predicts a peak of response which
is not present in the experimental results. This peak was also observed in (7) and corresponds to a matching
between the forcing frequency and the characteristic convection time of a vortical structure through the cavity’s
mouth.
4.2 Spectral Proper Orthogonal Decomposition
In (7), the LNSE was able to predict the right amplitudes for the coherent response. In that study however,
the turbulent fluctuations were weak compared to the ones encountered in the present study. The strong noncoherent velocity fluctuations could therefore have a non-negligible impact on the coherent structures and might
affect their frequency through non-linear phenomena which are not taken into consideration by the LNSE model.
In order to capture coherent energetic structures whose frequencies are close to the forcing frequency, we makeuse of the spectral proper orthogonal decomposition (SPOD) introduced in (17). This improved POD (proper

5303

(a) ûx . left: LNSE. right: Harmonic componentof the experimental flow(b) ûy . left: LNSE. right: Harmonic componentof the experimental flow

Figure 3. Comparison of the velocity fluctuations predicted with the LNSE and the ones extracted from the PIV
measurements. Experimental conditions: U=76m/s, ua = 0.5 m/s. Colorbar values in m/s. The colorbars limit
are based on the maximum absolute value of the displayed variable.

Figure 4. |ûy | profiles on vertical slices across the cavity’s mouth. Blue plain lines: fourier component at theforcing frequency, extracted from the experimental results. Red dotted-dashed lines: LNSE. ua =0.2 m/s, U=76
m/s.

5304

Figure 5. |ûy | profiles on vertical slices across the cavity’s mouth. Blue plain lines: most energetic SPOD modeof the experimental flow in the vicinity of the forcing frequency. Red dotted-dashed lines: LNSE. ua =0.2 m/s,
U=76 m/s.
orthogonal decomposition) relies on a diagonal filtering of the autocorrelation matrix to extract spatially andspectrally coherent structures. SPOD decomposes the coherent fluctuations of a flow into the sum a few ’mostenergetic’ modes having distinct spectral contents. SPOD allows to extract a coherent velocity structure whose
spectrum is located in a band around the forcing frequency. Figure 5 shows a comparison between the LNSE
prediction and the SPOD mode. The amplitude of the coherent response obtained from the SPOD is comparable
to the LNSE response. The maximum amplitude of the SPOD mode is obtained at 1200Hz, while the LNSE
predicts a peak of response at 1400Hz. This difference could be due to the fact that the experimental velocity
profile is not constant in the spanwise direction, but its values decreases in the boundaries layer attached to
the lateral walls. The 2D slice observed by the PIV is located at a half-span position, where the velocity is
maximal, leading to an overestimation of the average convective velocity in the spanwise direction, which could
be an explanation for the overestimation of the maximal response frequency.

5

CONCLUSIONSThe linearized Navier-Stokes operator around a time-averaged acoustically forced flow was shown to be able
to predict the spatial distribution of the acoustic fluctuation. Amplitude prediction is still challenging due tothe strong turbulent fluctuations, but kinetic energy based data processing methods like SPOD allow to extract
coherent structures whose amplitude is in agreement with the LNSE predictions. Both the experiments and the
LNSE predicts a peak of response at a given frequency, but in the present configuration, a 2D analysis seems
not to be sufficient to predict accurately the frequency at which the peak of response occurs.

ACKNOWLEDGEMENTS
This project has received funding from the European Union’s Horizon 2020 research and innovation programme
under Grant Agreement No 765998.

REFERENCES[1] Krishnamurthy, K., Acoustic radiation from two-dimensional rectangular cutouts in aerodynamic surfaces,
Tech. Rep. 3487, California Institute of Technology - NACA, 1955.
[2] Plumblee, H.E.; Gibson, J.S.; Lassiter, L.W., A theoretical and experimental investigation of the acoustic
response of cavities in an aerodynamic flow, Tech. Rep. WADD-61-75, U.S. Air Force, 1962.
[3] Rossiter, J.E., Wind tunnel experiments on the flow over rectangular cavities at subsonic and transonic
speeds, Tech. Rep. 3438, Aeronautical Research Council - Reports and memoranda, 1964.
[4] East, L.F., Aerodynamically induced resonance in rectangular cavities, J. Sound Vib, Vol 3 (3), 1966, pp
277–287.

5305

[5] Tam, C.K.; Block, P.J.W., On the tones and pressure oscillations induced by flow over rectangular cavities,
J. Fluid Mech., Vol 89 (2), 1978, pp 373–399.
[6] Mantič-Lugo, V.; Gallaire, F., Self-consistent model for the saturation mechanism of the response toharmonic forcing in the backward-facing step flow, J. Fluid Mech., Vol 793, 2016, pp 777–797.
[7] Boujo, E.; Bauerheim, M.; Noiray, N., Saturation of a turbulent mixing layer over a cavity: response toharmonic forcing around mean flows, J. Fluid Mech., Vol 853, 2018, pp 386–418.
[8] Howe, M.S., Contributions to the theory of aerodynamic sound, with application to excess jet noise andthe theory of the flute, J. Fluid Mech., Vol 71 (4), 1975, pp 625–673.
[9] Howe, M.S., Theory of vortex sound, Cambridge University Press, 2002.
[10] Barkley, D., Linear analysis of the cylinder wake mean flow, Europhys. Lett., Vol 75 (5), 2006, pp
750–756.
[11] Garnaud, X.; Lesshaft, L.; Schmid, P.J.; Huerre, P., The preferred mode of incompressible jets: linear
frequency response analysis, J. Fluid Mech., Vol 716, 2012, pp 189–202.
[12] Lien, K.; Monty, J.P.; Chong, M.S.; Ooi, A., The entrance length for fully developed turbulent channel
flow, In 15th Australasian Fluid Mechanics Conference, Sydney, 2004.
[13] Hussain, A.K.M.F.; Reynolds, W.C., The mechanics of an organised wave in turbulent shear flow, J. Fluid
Mech., Vol 41 (2), 1970, pp 241–258.
[14] Viola, F.; Iungo, G.; Camarri, S.; Porté-Agel, F.; Gallaire, F., Prediction of the hub vortex instability in awind turbine wake: stability analysis with eddy-viscosity models calibrated on wind tunnel data, J. Fluid
Mech., Vol 750, 2014.
[15] Pope, S.B., Turbulent flows, Cambridge University Press, 2000.
[16] Seung-Ho, J.; Jeong-Guon, I., On the multiple microphone method for measuring in-duct acoustic properties in the presence of mean flow, The journal of the acoustical society of america, Vol 103 (1520),
1998.
[17] Sieber, M.; Paschereit, C.O.; Oberleithner, K., Spectral proper orthogonal decomposition, J. Fluid Mech.,
Vol 792, 2016, pp 798–828.

5306