1 Introduction

Transonic buffeting is usually characterised in the literature [1,2,3] as a structural response to an aerodynamic instability, which causes significant low-frequency fluctuations in the aerodynamic lift forces. This aerodynamic instability has been observed in experimental [4, 5] and numerical [6] studies of rigid airfoils as well and is known as “transonic buffet”. It is generally assumed that the structural response of the wing is triggered by resonance effects after the disturbance amplitude reaches a sufficient magnitude. It is of great interest to be able to define buffet boundaries as precisely as possible in order to fully exploit and potentially extend the safe flight envelope. However, despite large experimental efforts, the self-sustaining mechanism is still not fully understood [1, 7]. Transonic buffet is often associated with large amplitude, autonomous shock oscillations, caused by the interactions between shock waves and separated shear layers [1]. While traditional explanations (e.g. acoustic feedback and wave propagation models) have difficulties to directly couple the shock motion with the low-frequency fluctuations in the lift [7], more recent studies describe transonic buffet as a global instability [8]. It is however not quite clear whether the shock motion plays a fundamental active role or is rather a symptom of the periodically accelerating and decelerating flow over the airfoil suction side. Also, “phase-locking” between the shock motion and a global buffet mode is possible. Paladini [9] challenged the importance of the shock motion for transonic buffet, but still suggests that the shock foot plays a major role. We characterise transonic airfoil buffet phenomena here in a more general manner as low-frequency oscillations in the lift coefficient at Strouhal numbers around \(St=fc/U_{\infty } \approx 0.1\) (corresponding to typical structural resonance frequencies of aircraft wings) with amplitudes greater than 5%, instead of directly linking it to the shock-oscillation frequency. It should also be noted that buffet on swept wings occurs at higher frequencies [10,11,12] and may be a distinct phenomenon [13,14,15].

In the present contribution direct numerical simulations (DNS) are performed over wing-sections at a moderate Reynolds numbers of Re = 500, 000 (based on the chord length c) and a Mach number of M = 0.7 considering Dassault Aviation’s V2C profile [16]. In the course of the TFAST project, experimental as well as numerical analysis has been carried out on that profile under buffet conditions [17,18,19,20]. There is great interest on global stability analysis of airfoils under buffet conditions (i.e. [6, 21, 22]), but very little work has been done to date using DNS data. In particular, we wish to explore the usefulness of the dynamic mode decomposition (DMD) procedure [23] to separate and analyse interactions between complex flow phenomena, such as buffet, shock waves, acoustic waves and flow instabilities that are known to exist for this case [16]. Furthermore, it can be useful to compare DMD modes with modes obtained by linear global stability analysis [24]. Based on recent results, we want to extend the investigations for the V2C profile, exploring buffet at angles of attack α = 3∘ and α = 4∘ (denoted as ‘3∘ case’ and ‘4∘ case’, respectively), and reduced Reynolds numbers.

After outlining the well-established methodology in the following section, providing references to relevant literature for more details, we review findings and conclusions of previous work on the 4∘ case based on Fourier analysis of the pressure and lift histories [16, 25]. In Section 4, the sensitivity of convective instabilities to the low-frequency dynamics is analysed, before results from a global stability analysis, and their limitations, are discussed. DMD is applied in Section 5 to extend the study of global modes and the flow dynamics. The sensitivity of DMD modes is then examined considering also DNS results at a decreased angle of attack of α = 3∘. To establish the reduced impact of shock waves on the buffet phenomenon in the present case, an extended Fourier analysis on the 3∘ test case is presented in Section 7, before summarising the conclusions in Section 8.

2 Methodology

All direct numerical simulations in this work were carried out using the high-order fully-parallelised multi-block finite difference in-house code SBLI with details in [26, 27]. A recently published paper [16] is concerned with the set-up of airfoil simulations. The dimensionless Navier-Stokes Equations (NSE) are solved by a fourth-order finite difference scheme in space and a low-storage third-order Runge-Kutta scheme in time. The temperature dependency of the dynamic viscosity is modelled by Sutherland’s law. Zonal characteristic boundary conditions are applied at the outflow boundaries, while integral characteristic boundary conditions at the remaining domain boundaries avoid wave reflections. In the farfield, an implicit sixth-order filter increases the numerical stability of the simulation. A total variation diminishing (TVD) scheme is used to capture shock waves, but is disabled in boundary layers and near the leading edge. The computational domain is divided into three blocks consisting of one C-block around the airfoil geometry and two H-blocks enclosing the wake-region and outflow. In order to include the blunt trailing edge of the original profile, while maintaining continuous metric terms up to the second order of derivatives, an open-source grid generator was developed and released on GitHub [28]. The reference grid consists of more than one billion points, considering a spanwise domain width of 5%c. The adequacy of the grid resolution is confirmed by a grid-refinement study, based on a spectral error-indicator analysis identifying critical regions in terms of grid-to-grid point oscillations [29]. The simulation arrangement, including a grid study confirming also a sufficient domain width of the 4∘ buffet case that is further analysed in the present contribution, has already been published in [16].

In order to analyse the linear stability, the flowfield is decomposed into a mean flow with superimposed disturbances. The mean flow is obtained from time- and span-averaged DNS solutions, whereas the disturbances are modelled by a normal-mode ansatz. For local linear stability theory, the compressible Orr-Sommerfeld equations (involving parallel-flow assumptions on the linearised NSE) are solved for a 2D flowfield using the in-house code NoSTRANA [27]. Applying a temporal approach, the solution of an eigenvalue problem provides the temporal growth rate corresponding to an angular frequency for a given set of streamwise and spanwise wave numbers. More details on this methodology, including a full derivation of the equations, flow assumptions, and application to time- and space-averaged baseflows around airfoils are provided by [27] and [30]. For the global stability analysis, the linearised NSE are solved directly, applying a normal-mode ansatz with prescribed spanwise wavenumbers. The large-scale eigenvalue problem is solved in matrix-free mode [31,32,33] using the Implicitly Restarted Arnoldi Method [34] in combination with SBLI as a timestepper. More details on this approach are available in [26].

In order to compare the calculated stability results with the unsteady direct numerical solution of the flowfield, the streaming dynamic mode decomposition (DMD) method is applied [23, 24, 35]. In contrast to the well-known proper orthogonal decomposition (POD), DMD aims to approximate the non-linear dynamics of an unsteady flow by mapping consecutive snapshots (xt and xt+ 1) of the flowfield by a linear operator (xt+ 1 = Axt) rather than finding an orthonormal basis. The aim is to approximate the non-linear dynamics of a general problem

$$ \frac{d \boldsymbol{x}}{d t}=f(\boldsymbol{x},t,\mu) $$
(1)

by the best possible linear system

$$ \frac{d \boldsymbol{x}}{d t}\approx \boldsymbol{A} \boldsymbol{x}, $$
(2)

which has an exact solution \(\boldsymbol {x}(t)=\sum {\hat {\Phi }_{i} e^{\omega _{i} t}}\), where \(\hat {\Phi }_{i}\) and ωi are the normalised DMD modes and corresponding complex eigenvalues, respectively. In contrast to the POD method, it does not rely on a model describing the dynamics of coherent structures. DMD for flows consisting of small perturbations of a baseflow can be compared to results from global stability analysis [23]. More details on the implementation of this method are also available in [36].

3 Unsteady Flow Structures at α = 4∘

The reference simulation, at Re = 500, 000 and an angle of attack of α = 4∘, is illustrated in Fig. 1 by means of pressure gradient contours, showing a distinct supersonic region over the upper airfoil surface, bounded by the red sonic line. Two dimensional (2D) simulations already show the formation of strong Kelvin-Helmholtz (KH) vortex structures in the airfoil aft section, initialised by upstream moving pressure waves (also known as Kutta waves) that are caused (in the 2D case) by a von Karman vortex street, which first appear with Strouhal numbers of St ≈ 20 − 30. This suggests a complex cascade mechanism, also involving shock-wave/boundary-layer interaction and the Doppler effect, that allows flow structures of high frequencies to interact with flow phenomena at significantly lower frequencies. After extruding the 2D solution, a self-sustaining laminar/turbulent boundary-layer transition mechanism sets in on both sides without further artificial excitation of the flowfield. A transition mechanism similar to [37] can be observed, where vortex-stretching of near-wall rib vortices [38] that are lifted up by strong 2D vortices promotes a rapid breakdown to turbulence (inset of Fig. 1). A 2D-like silhouette of the strong vortices can still be observed in the fully turbulent section. Furthermore, those turbulent structures interact with each other as well as the potential flow.

Fig. 1
figure 1

Magnitude of pressure gradient \((|\frac {\partial p}{\partial x}|+|\frac {\partial p}{\partial y}|)\) (white = 0, black = 1), where the red curve denotes the sonic line. Sketched lines mark described flow phenomena. The insert shows Q-criteria surfaces coloured by vorticity of the upper transition region outlined by the blue dashed rectangle. The pressure gradient along the white dashed line is plotted later (Fig. 3) as a function of space and time

Full size image

In the aft section of the airfoil, strong acoustic radiation can be observed from multiple sources. Black contours in Fig. 1 indicate strong pressure gradients. Approaching the supersonic region, upstream-propagating acoustic waves can be seen to accumulate and form stronger pressure waves (labelled PW in Fig. 1). Eventually, those pressure waves turn into shock waves (labelled SW in Fig. 1) and propagate upstream. Tijdeman [39] distinguishes between three different types of shock behaviour, known as type A, type B, and type C shock motion. Type A shock motion describes a single permanent shock wave oscillating back and forth over large parts of the airfoil, while the shock wave temporarily disappears during the downstream excursion for type B shock motion. Type C shock motion is significantly different from types A and B, as there is no permanent shock wave observed. Instead, there are periodically-generated upstream-propagating shock waves leaving the airfoil at the leading edge that continue moving upstream into the oncoming free stream. All three types have been observed in simulations and experiments. Experimental studies of the same airfoil at higher Reynolds numbers [18] show type A shock motion. The difference between experiments and the present simulations is almost certainly due to Reynolds-number effects [25]. However, confirmation of this would require additional simulations of the same test case at Reynolds numbers of the order of several millions. The present simulations show how acoustic waves circumvent the supersonic region at higher speeds than the upstream-propagating shock waves and introduce additional disturbances into the supersonic region from above. A complex interaction between shock waves, pressure waves, and reflections at the boundary layer is observed. One can also observe Mach-like wave patterns (labelled MW in Fig. 1) on both sides that are caused by acoustic waves travelling upstream within the separated boundary layer (highlighted magenta in Fig. 1).

Figure 2a shows the lift coefficient CL as a function of time, highlighting low-lift phases (LLP) and high-lift phases (HLP) in blue and red respectively. The statistical sample size for the HLP and LLP is limited, as the statistics files were only written approximately every two time units (TU) of DNS run time. The time segments for the high- and low-lift phases were selected manually. A distinct low-frequency oscillation in the lift coefficient (up to 12% deviation from the mean value) is observed at a Strouhal number of St = 0.12, which is typical of buffet frequencies [7]. The low-frequency behaviour can also be clearly observed in lift spectra and pressure signals as a function of time at various locations along the surface (including the leading- and trailing-edges) and in the freestream [16]. While high-amplitude lift fluctuations associated with buffet are typically observed at the same frequency as the back- and forth-moving shock waves, there is no obvious correlation in the present case, as upstream-propagating shock waves are generated at significantly higher frequencies (St = 0.4 − 0.7) and leave the airfoil via the leading edge. While studies of the same airfoil at M = 0.7 and significantly higher Reynolds numbers close to experimental conditions (Re = 3 ⋅ 106) do not suggest transonic buffet at α = 4∘, low-frequency oscillations at St ≈ 0.1 are observed for α > 5∘ [17, 20]. The latter results, comparing delayed detached-eddy simulations (DDES), implicit large-eddy simulations (ILES), and unsteady Reynolds averaged Navier-Stokes (URANS) approaches, show high sensitivity of the low-frequency buffet phenomenon to the modelling methodology applied. This is also reported by an LES study using different grids and methods to reproduce the low-frequency behaviour of the present test case [40].

Fig. 2
figure 2

a Lift coefficient CL as a function of time, and b wall-pressure coefficient Cp,w as a function of x. High-lift phases (HLP) and low-lift phases (LLP) are shown in red (dash-dotted) and blue (dashed), respectively. c Root-mean-square (rms) of wall-density fluctuations (\(\rho -\overline {\rho }\)) as a function of x. The black line is based on the unfiltered wall density, whereas the blue and red lines are calculated from band-pass filtered solutions of the wall density containing frequencies below and above St = 0.2, respectively

Full size image

The black line in Fig. 2b shows the suction-side wall-pressure coefficient Cp,w as a function of x, averaged over the full runtime of 25 time units. Phase-averaged Cp,w corresponding to low- and high-lift phases (as shown in Fig. 2a) are denoted by the blue dashed and red dash-dotted lines in Fig. 2b, respectively. The main differences are observed at x ≈ 0.7 (region of turbulent breakdown), as the flow is decelerated further downstream by a steeper increase in the surface pressure during HLP compared to LLP. Furthermore, a small plateau is observed at x ≈ 0.4 during HLP, which does not appear in the overall- and LLP-averaged Cp,w. The black line in Fig. 2c shows the root-mean-square (rms) of wall-density fluctuations (\(\rho ^{\prime }\)) as a function of x, which is calculated at each chord position according to

$$ \begin{array}{@{}rcl@{}} \rho_{w,\text{rms}}(x)=\frac{1}{T_{tot}}{\int}_{t=0}^{T_{tot}} \sqrt{(\rho^{\prime}(x,t))^{2}} dt, \end{array} $$
(3)
$$ \begin{array}{@{}rcl@{}} \text{with} \ \rho^{\prime}(x,t) = \rho(x,t) - \overline{\rho(x)}=\rho(x,t) -\frac{1}{T_{tot}}{\int}_{t=0}^{T_{tot}} \rho(x,t) dt, \end{array} $$
(4)

where the total runtime for the present case denotes Ttot = 25. For isothermal surfaces, the wall density is directly proportional to the wall pressure so that

$$ C_{p,w}=\frac{2(\rho_{w}-1)}{\gamma M^{2}_{\infty}}. $$
(5)

As already suggested in Fig. 2b, we observe strong fluctuations for 0.2 < x < 0.4 and 0.6 < x < 0.8. To better understand the contribution of the low-frequency phenomenon to the total wall-density fluctuations, the instantaneous wall-density signal is filtered by applying a Fourier band-pass filter. After transforming one-dimensional time signals of \(\rho ^{\prime }\) at each x-position along the airfoil surface into Fourier space, the power-spectral density is set to zero beyond selected frequency ranges and then transformed back to physical time space. The root-mean-square is again calculated according to Eq. 3 using filtered wall-density fluctuations instead of \(\rho ^{\prime }\). The fluctuations corresponding to the low-frequency phenomenon with St < 0.2 (denoted by the blue line in Fig. 2c) are dominant in the fore part of the airfoil (x < 0.4) and at x ≈ 0.7. The contribution of the low-frequency phenomenon to the overall ρw,rms drops significantly between these regions, where upstream-propagating shock waves are observed in DNS visualisations generated at frequencies in the range 0.4 < St < 0.7. Also in the turbulent part for x > 0.8, fluctuations with St > 0.2 (corresponding to the red line in Fig. 2c) are significantly stronger. While the blue line corresponding to the low-frequency fluctuations exhibits a distinct peak at x ≈ 0.7, the red line shows a plateau for 0.65 < x < 0.75, where upstream propagating shock waves are formed. Based on Fig. 2c, no direct correlation between fluctuations at low-frequencies for St < 0.2 (corresponding to transonic buffet) and higher frequencies (corresponding to convective structures and shock waves) can be established. It is, however, not possible to rule out a potential interaction between the buffet phenomenon (St = 0.12) and phenomena at significantly higher frequencies (particularly near x = 0.7). Therefore, it is necessary to study the physical phenomena, such as shock waves, acoustic waves, and convective structures more precisely by other means.

The full spatio-temporal behaviour of shock- and pressure-waves along the white-dashed line over the airfoil suction side in Fig. 1 is provided by Fig. 3. Figure 3a shows the probability of compression- (red bars) and expansion-waves (blue bars) occurring as a function of the x-location. Cells cut by the sonic line are assigned a value of unity, while remaining cells are denoted by zero. An integration over time suggests how much time shocks spend at a given location. The probability is obtained by dividing the result by the total runtime of 25 time units. A value of unity would imply a steady compression- or expansion-wave occurring at that fixed position for the full time span. As mentioned before, there are no steady or permanently visible shock waves observed at this flow condition. Even though the behaviour of the shock waves seems rather chaotic, there are regions on the airfoil where shock waves tend to spend more time. Considering the red bars, corresponding to compression waves, one can observe increased probability of shock waves occurring at x ≈ 0.42, x ≈ 0.37, just before half chord at x ≈ 0.48, and within the transition region 0.6 < x < 0.775. It is notable that the red line in Fig. 2 showing ρw,rms for St > 0.2 is rather smooth and does not exhibit local peaks as suggested by Fig. 3a.

Fig. 3
figure 3

a Probability of compression- (red) and expansion-waves (blue) occurring as a function of the x coordinate. b Streamwise pressure gradient on the suction side along a iso-curve at η = 200 as a function of space and time. The dotted line sketches the transition region moving at the buffet frequency of St = 0.12. c Probability of the compression-wave occurrence (red bars) and lift coefficient CL (black line) as a function of time

Full size image

The pressure gradient along the gridline at η = 200 (defined by the dashed white line in Fig. 1) is plotted as contours in Fig. 3b as a function of the chordwise location (x) and simulation time (t). Green curves indicate the sonic line (iso-curves for M = 1) and red contours indicate strong adverse pressure gradient (∂p/∂x > 0), while blue contours indicate favourable pressure gradients. Figure 3b also shows strong downstream-convecting high-frequency pressure fluctuations at the right hand side of the plot (x > 0.7). These high-frequency oscillations mainly reflect the acoustic nearfield generated by large downstream-propagating turbulent vortex structures. Three particularly noisy patches near the trailing edge appear at the main buffet frequency (St = 0.12), caused by turbulent structures in the thick shear layer reaching up to the monitored η-gridline.

Near the transition region, a narrow band (centred on the dotted curve in Fig. 3b) of upstream-propagating high-frequency acoustic waves can be observed, which is hidden further downstream by the noise of the large turbulent vortices. These upstream-propagating acoustic waves slow down as they approach the supersonic region, where shock waves are formed. The location of this band is unsteady and varies mainly with the buffet frequency of St = 0.12. Upstream of this band, the flow is mainly supersonic (bounded by the green sonic lines) and shielded from high-frequency oscillations. In the left half of Fig. 3b, pairs of compression- (red) and expansion-waves (blue) are observed, which is also shown by the V-shape of the shocks in the snapshot of Fig. 1. As the waves propagate upstream, the subsonic regions bounded by the green sonic lines between compression- and expansion-waves disappear near the airfoil surface.

The bars in Fig. 3c on the right of b correspond to the number of shock waves counted at each time step. Similarly to Fig. 3a, cells cut by the sonic line are assigned a value of unity, while remaining cells are denoted by zeros. An integration over x gives the values for the red bars in Fig. 3c. The black curve denotes the lift coefficient, importantly showing no direct correlation between the low-frequency buffet phenomenon and the number of apparent shock waves. At least for the present configuration it seems that buffet and shock waves are independent phenomena. The low-lift phases correlate well with the patches of high pressure gradient, at x ≈ 0.9, that were mentioned in connection with Fig. 3b. It can also be observed that strong pressure waves are halted during low-lift phases at x ≈ 0.35, before they continue moving upstream. This continuation of the upstream motion of strong pressure waves sets in during phases where the lift is recovering after a minimum and the flow over the suction side is again accelerated. Further discussion of the unsteady flow phenomena is given in [16, 25], and [36].

4 Local and Global Linear Stability

This section aims to identify linear instabilities and describe them by means of frequencies, growth rates, and spatial coherence. Before attempting global linear stability analysis, a temporal local linear stability theory (LST) approach is applied to time- and span-averaged flowfields of the suction side of a direct numerical simulation of an airfoil at α = 4∘ with a total run time of 25 time units. As a consequence of the observed buffet phenomenon, the flowfield and in particular the boundary- and shear-layers change significantly. The contour-plots of Fig. 4a and b show the z-vorticity component (ωz) of phase-averaged flowfield considering only high-lift phases (HLP) or low-lift phases (LLP) according to Fig. 2a, respectively. Due to the limited number of simulated low-frequency cycles, there are still traces of instantaneous flow features (especially in the freestream) that are not completely averaged out. During low-lift phases, the separation bubble moves upstream and the flow separation becomes more pronounced. The local peak of the blue line in Fig. 2c associated with low-frequency ρw,rms at x ≈ 0.3 seems to correspond to the intermittent flow separation in that region. Iso-curves in Fig. 4c show a direct comparison of the shear layers for HLP (red) and LLP (blue). For HLP, the flow follows the contour longer, whereas the shear at the lower corner of the blunt trailing edge is significantly reduced. The shear layer along the suction side surface does not change significantly.

Fig. 4
figure 4

a Phase-averaged vorticity fields for high-lift phase and b low-lift phase. c Iso-curves for ωz = 50 of high-lift phase and low-lift phase in red and blue, respectively

Full size image

Using a mean flow that is averaged over the total run time of the simulation fails to take the periodic variations in lift into account, and it is clear that shear- and boundary-layer characteristics are significantly influenced by the low-frequency oscillations. Therefore the phase averaged flows shown in Fig. 4 are also analysed with respect to linear instabilities. This is justified on the ground that there still exists a wide frequency separation between the buffet mode and the KH modes. Figure 5 shows the temporal growth rate (ωi) as a function of surface distance s for an angular wave number of ωr = 125 (St ≈ 20), considering only 2D modes (spanwise wavenumber β = 0). The frequency range agrees with Kelvin-Helmholtz instabilities reported by [16] for the present test case. Similar flow structures (St ≈ 25) in a simulation of a high-pressure turbine vane could also be associated with linear instabilities by [30]. The wavy pattern is likely to be due to upstream-moving pressure waves interacting with the boundary layer combined with short time averaging. We can clearly see that the total time average (black curve) underestimates linear instabilities in comparison to high- (red curve) and low-lift (blue curve) phases. Compared to LLP, the boundary layer at HLP shows higher temporal growth rates peaking closer to the trailing edge. Instantaneous snapshots confirm that KH roll-ups form further upstream at LLP. Particularly the LST results at LLP show increased growth rates locally corresponding to regions with high shock probability in Fig. 3a.

Fig. 5
figure 5

Temporal growth rate ωi as a function of the chord position for an angular frequency of ωr = 2π ⋅ St = 125. The line colour corresponds to time- and span-averaged mean flows over high-lift phases (red), low-lift phases (blue) and the full run time (black)

Full size image

We next analyse the global instability of the 2D mean flow with spanwise wavenumber β = 0. Considering the flow around a NACA 0012 airfoil at Re = 200, 000 and M = 0.4, Fosas De Pando et al. [41] report a region in the spectrum that is dominated by distinct equally-spaced frequencies (tonal noise) around a maximum peak at St ≈ 7 that corresponds to a stable global mode. An impulse response analysis by [42] showed the vivid interaction between suction and pressure side at that frequency and suggested a feedback mechanism due to pressure waves that are scattered at the trailing edge and form upstream moving acoustic waves. Analysing the time- and span-averaged mean flow over 24 time units, a similar stable mode can be observed in global stability results at St = 5.89 suggesting a growth rate of ωi = − 0.019 (negative growth rates denote damping). The divergence field of that global mode in Fig. 6a shows regions of high growth rates in the separation regions on both sides. This global mode also involves upstream-travelling acoustic waves originating at the trailing edge. While those waves can travel along the pressure side without any restrictions, they are slowed down on the suction side approaching the supersonic region. Near the shock wave that bounds the supersonic region in the downstream direction, those waves are compressed as the phase speed decreases to zero. The pressure waves circumventing the supersonic region slide along the sonic line and introduce disturbances that are reflected at the airfoil surface. This mode thus contains a dynamic coupling between separation regions and the trailing edge, via upstream moving pressure waves. The z-vorticity of this eigenmode is shown in Fig. 6b.

Fig. 6
figure 6

Divergence field (top) and vorticity field (bottom) of a global mode at St = 5.89

Full size image

Despite the success in extracting the above global mode, in general we found global analysis to be of limited value for the current case. This is partly because the flow is already undergoing moderate buffet and the linearised approach is no longer relevant. In addition, the convergence of global modes was limited by under-resolution of the acoustic field far from the airfoil. As an alternative, we next consider the DMD approach.

5 Dynamic Mode Decomposition

In this section, a dynamic mode decomposition (DMD) of the test case with 4∘ angle of attack is performed considering two sets of instantaneous snapshots over 25 time units, which are in both cases a combination of 2D xy- and xz-planes (located within the shear layers) to capture 3D characteristics of the flow, while keeping the required memory low. The following section will consider the sensitivity of the results to the angle of attack. As the modes discussed in this contribution are mainly 2D, visualisations will focus on the xy-plane. After discussing the observed DMD modes, a sensitivity analysis using additional data sets with different sample rates and time segments is presented in the next section, so that the reader can better interpret the observations. Figure 7 shows the growth rates (denoted by blue circles corresponding to the left-hand scale) as a function of St for DMD modes corresponding to one set of 250 snapshots with a step size of 0.1 time units (top plot) and another set of 1249 snapshots with a step size of 0.02 time units (bottom plot), where negative growth rates correspond to damping. The green symbols (+ ) in the top plot of Fig. 7 denote the amplitudes normalised by the mean flow mode corresponding to the right-hand scale and for the current case suggest a correlation between normalised amplitudes and growth rates of DMD modes. Due to the limitations in random-access memory, it is necessary to either keep the sampling rate low to compute a larger number of modes (as a POD projection is applied to approximate the mapping matrix A in Eq. 2), or vice-versa. A lower sampling rate, however, is not able to reconstruct high-frequency modes, as shown in Fig. 7. Consistency between data sets with different sampling rates is established in [36] comparing a mode that appears in both sets.

Fig. 7
figure 7

Growth rate of DMD modes (denoted by blue circles corresponding to the left-hand-scale) as a function of their Strouhal number considering 250 snapshots with a step size of 0.1 time units (top plot) and another set of 1249 snapshots with a step size of 0.02 time units (bottom plot). The top plot shows also normalised amplitude denoted by green + -symbols corresponding to the right-hand scale

Full size image

A DMD mode with a similar frequency (St = 6.39) to the global mode (St = 5.89) and decreased damping at St = 6.39 (marked green in the bottom plot of Fig. 7) is shown in Fig. 8. The divergence field of the DMD mode (Fig. 8a) shows acoustic waves originating from the trailing edge on both sides of the airfoil, which are also observed in the global mode in Fig. 6. Both the DMD mode and the global mode show traces of shock waves near the half chord position. The DMD mode in Fig. 8a, however, is very noisy within the shear layers and wake, so that the structures shown by the global mode at x ≈ 0.8 on the suction side and at x ≈ 0.9 on the pressure side are hidden. These structures can be better seen in the velocity and density shown in Fig. 8b and c.

Fig. 8
figure 8

DMD eigenmode at St = 6.39 corresponding to a global mode at St = 5.89 showing divergence field (top), vorticity field (middle), and density field

Full size image

After having identified a DMD mode, which is similar to a global mode, we extend the study to other flow features. Selected eigenfunctions of the modes marked red in Fig. 7 are shown in Fig. 9, plotting contours of density (left column) and streamwise velocity component (right column). The top row shows a mode at St = 19.3, picking up the KH roll-ups, which are associated with convective linear instabilities. The shape of the eigenmode is reminiscent of structures that are observed in movie snapshots like Fig. 1. Figure 7 shows several modes at frequencies between 15 < St < 25 with similar or even lower damping rates with similar shapes of eigenmodes. Considering the results of the local linear stability analysis, the frequency of the KH instabilities is expected to vary significantly due to the low-frequency dynamics of the flowfield.

Fig. 9
figure 9

Selected DMD eigenmodes at different frequencies. The left and right columns show the density (ρ) and streamwise velocity (u) field, respectively

Full size image

There are big turbulent vortices observed downstream of the laminar/turbulent transition region with Strouhal numbers around St = 1.8. These modes can also be found in the DMD spectrum. The density field in the second row of Fig. 9 shows strong oscillations in the aft section of the airfoil and in the wake. In addition, upstream-moving pressure waves, originating at the trailing edge can be observed. A phase shift of the oscillations within the upper-side shear layer and the freestream can be seen in the velocity field of that mode.

Besides the low-frequency peak at St = 0.12 (the last row of Fig. 9) corresponding to the buffet phenomenon, distinct low-frequency modes at St = 0.6 and St = 0.45 are shown in the third and fourth rows of Fig. 9, respectively. The variation of the spacing between green sonic lines at x ≈ 0.6 in Fig. 3b in the range of Δt = 1.5 − 2.5 time units suggests a production rate of upstream-propagating shock waves in the range of St ≈ 0.4 − 0.7. The modes at St = 0.6 and St = 0.45 in Fig. 9 are chosen as representative modes for the shock motion, as they have high normalised amplitudes and low damping rates in the DMD spectrum (marked red on Fig. 7). Furthermore, high power-spectral density is also observed in the Fourier spectra of surface pressure probes at these frequencies [25]. Both eigenmodes look similar, but have a clear phase shift. These modes seem to be strongly coupled with the shock dynamics and fluctuations in the wake. Their eigenmodes are mainly active on the suction side, but also include acoustic waves moving over the pressure side and also comprise acoustic waves circumventing the supersonic region. These modes can be interpreted in the light of other observations of higher frequencies on airfoils. Performing a large-eddy simulation of a supercritical laminar airfoil (OALT25) at α = 4∘ and M = 0.735, but at a higher Reynolds number of Re = 3 ⋅ 106, Dandois et al. [43] reported shock motion (a permanent back and forth moving shock wave) at significantly higher Strouhal numbers (St = 1.2) compared to typical buffet frequencies, and linked it to a breathing phenomenon of the separation bubble associated with downstream convecting vortices. Similar phenomena were reported by Memmolo et al. [17] analysing the V2C airfoil at α = 7∘ and high Reynolds numbers using large-eddy simulation. In a URANS parameter study exploring the buffet domain varying the angle of attack (3 ≤ α ≤ 9.5) and Mach number (0.55 < M < 0.75), an interesting coexistence between type A and type C shock motion was reported by Giannelis et al. [44]. The disturbances associated with the type C shock motion emerged from downstream-convecting recirculation pockets that produced oscillations at the trailing edge at significantly higher frequencies compared to the low-frequency lift oscillation associated with buffet. In the present case, an acoustic propagation model, originally proposed by Lee [7] to predict transonic buffet, is able to approximate the frequency range associated with the periodic shock motion, but not the low-frequency buffet phenomenon [36].

The mode with St = 0.12 (in the last row of Fig. 9) agrees with typical transonic buffet frequencies in the literature and can also be extracted from the lift coefficient over time. In particular, the density fluctuations in the rear part of the airfoil agree well with observations by Fukushima and Kawai [45]. Strong fluctuations of the modes associated with the shock motion at x ≈ 0.6 are located slightly upstream of high amplitudes in the buffet mode (St = 0.12) at x ≈ 0.7. High amplitudes are also observed in supersonic regions around x = 0.4, were Sartor et al. [46] reported high receptivity of the global buffet mode along characteristic lines for a OAT15A airfoil at buffet conditions. In the velocity field, the suction and pressure sides are clearly separated by a phase shift, indicating an opposed streamwise oscillation of the shear layers. Furthermore, there seems to be a phase shift between the separation regions and the shear layers. In particular, the velocity field of the DMD buffet mode is qualitatively very similar to the unstable global buffet mode reported by [47]. Density fluctuations of the modes with St < 1 are not only observed in the streamwise direction, but also in the wall-normal direction.

Fluctuations near the trailing edge in DMD modes corresponding to shock motion (St = 0.45 and St = 0.6) are reminiscent of the buffet mode at St = 0.12. The difficulties in fully separating shock motion from the buffet phenomenon using DMD can be either due to the sensitivity to the selected data sets, or to large lift oscillations at well-established buffet conditions. Therefore, the robustness of DMD to sample rates and sampling time needs to be studied, before analysing a test case closer to onset conditions at α = 3∘ in order to establish the (minor) impact of shock motion on transonic buffet in the following sections.

6 Sensitivity Study of DMD Results

Dynamic mode decomposition has been widely used to study the flow dynamics over airfoils (e.g. [48,49,50]). Comparing DMD with POD methods, [51] found DMD favourable to study transonic buffet. In this section additional data sets of the simulation at α = 4∘ are analysed using DMD. Furthermore, a simulation is performed at a decreased angle of α = 3∘, using the same DNS set-up as described above to study the sensitivity of low-frequency phenomena (observed at St < 1) to the angle of attack. Figure 10 shows the lift coefficient as a function of time for the 4∘ case (grey and black lines) and 3∘ case (red lines). The 3∘ case was restarted from a solution of the 4∘ case at t = 13 and underwent a transient between t = 13 − 32 (denoted by the light red-coloured line), before reaching the developed quasi-periodic flow (denoted by the dark red line). For the mode analysis it is useful to perform a sensitivity analysis of the dominant DMD modes to the sampling period. In order to save random-access memory, it would be also interesting to find out whether 2D xy-slices alone are sufficient to accurately capture the DMD modes corresponding to low-frequency phenomena. Therefore, three data sets of the 4∘ case are analysed and summarised in the first three rows of Table 1, before comparing the DMD modes of data sets at different angles of attack (last two rows in Table 1).

Fig. 10
figure 10

Lift coefficient of test cases at 4∘ and 3∘ angles of attack, corresponding to black curves and red curves, respectively. The lighter coloured lines were disregarded from the statistics

Full size image
Table 1 Summary of data sets used for dynamic mode decomposition for test cases at angles of attack α = 3∘ and α = 4∘
Full size table

Data set 4a was discussed in the previous section in connection with the top plot of Fig. 7 and the bottom three DMD modes in Fig. 9. Data set 4b contains the same number of snapshots (250) and sampling rate (sampling every 0.1 time units), but uses only xy-slices. Both data sets 4a and 4b consider the full runtime of 24 time units (corresponding to the grey and black lines in Fig. 10) sampling 10 snapshots per time unit. Data set 4c has a slightly lower sampling rate of 8.33 snapshots per time unit (sampling every 0.12 time units) using 140 xy-slices over 16 time units (corresponding to the black line in Fig. 10). Data set 3a comprises data from the test case at a lower angle of attack and can be compared to data set 4c.

Figure 11a shows the normalised amplitude of DMD modes for the case at α = 4∘, where data set 4a is represented by black rectangles, while blue and red triangular symbols correspond to data sets 4b and 4c, respectively. Perfect agreement between all three data sets can be observed for the buffet mode at St = 0.12 (outlined red). Also the DMD mode at St = 0.6 seems to be consistent, while the amplitude slightly varies (modes are outlined green). For data set 4a, we find a DMD mode with a high amplitude at St = 0.45. Looking at the DMD modes within the magenta ellipse in Fig. 11a, data set 4b again shows good agreement, whereas data set 4c with a shorter sampling time suggests modes at slightly higher and lower frequencies (St = 0.42 and St = 0.61, respectively). This is due to the variation of time scales corresponding to the shock motion, which was also observed in Fig. 3b. DMD modes found in the range of St = 0.4 − 0.6 are, however, qualitatively similar and can be clearly linked to the shock behaviour. Data set 4c gives a DMD mode with increased amplitude at St = 0.24 (outlined by the blue ellipse in Fig. 11a). This DMD mode is thus a super-harmonic of the buffet mode (St = 0.12). Varying the sampling rate, [51] also reported low sensitivity of low-frequency DMD modes (St < 1) for sampling rates greater than 3 snapshots per time unit (corresponding to 24 snapshots per buffet cycle).

Fig. 11
figure 11

Normalised amplitude of low-frequency DMD modes (St < 1) considering a three data sets of the 4∘ case (first three rows in table 1), and b data sets of the 4∘ and 3∘ test cases (last two rows in table 1). The modes corresponding to the filled symbols are shown in Fig. 12

Full size image

After having shown the robustness of DMD modes to the selection of data sets, Fig. 11b shows the normalised amplitudes comparing DMD modes at different angles of attack considering a sampling time of 16 time units. Similar to Fig. 11a, red symbols correspond to α = 4∘ and black circles to α = 3∘. The filled symbols are the modes of interest. Again, we can see good agreement for the buffet mode St ≈ 0.12. At lower angles of attack, the buffet frequency seems to slightly increase, while the amplitude decreases. Also the lift fluctuations in Fig. 10 for the 3∘ case are smaller than for the 4∘ case. Even though the frequency of well-established buffet oscillations often increases with the angle of attack [1], some experimental studies of laminar-flow airfoils near buffet-onset conditions [52] show a similar trend as for the present test case. For the α = 4∘ simulation, this can also be due to the short sampling time. A significant decrease in normalised amplitude is observed for the DMD modes associated with the shock motion and there is no significant peak observed in the DMD spectrum at α = 3∘ in the range of St = 0.4 − 1.

Figure 12 shows density fluctuations of the DMD modes corresponding to the black filled circles and red filled triangles of Fig. 11b on the right hand side and left hand side, respectively. Modes at comparable frequencies, but different angles of attack, are qualitatively similar and animations, which are available in the corresponding online repository (https://doi.org/10.5258/SOTON/D1004), help to clarify the dynamics of the discussed modes. In particular, the DMD modes corresponding to the low-frequency phenomena (the bottom plots of Fig. 12) have a very similar structure, albeit with a 180∘ phase shift. Towards the front of the airfoil, the St = 0.12 mode shows oblique structures, which can also be observed in the modes at higher frequencies at α = 4∘. These structures may be reflections of strong pressure waves near the leading edge. These are weaker at α = 3∘ and we do not observe oblique structures for the St = 0.15 mode. The modes at St = 0.61 and St = 0.42 in the 4∘ case show very strong fluctuations at the trailing edge reminiscent to those in the buffet mode, whereas the corresponding fluctuations at St = 0.55 and St = 0.44 for the 3∘ case are significantly weaker. For the 3∘ case, fluctuations in the shear layers and wake are more clearly visible, which are linked to downstream-convecting turbulent vortex structures. The DMD mode at St = 0.44 for the 3∘ case shows fluctuations on the suction side at the half-chord position and represents mainly the upstream motion of the shock waves, while the DMD mode at St = 0.55 also involves acoustic waves circumventing the supersonic region and downstream excursions of shock waves. In both the 3∘ and 4∘ cases, strong streamwise fluctuations of modes in the range St = 0.4 − 0.6 seem to be in between regions where the low-frequency buffet mode shows high amplitudes.

Fig. 12
figure 12

Contour plots showing density fluctuations of DMD eigenmodes corresponding to the filled symbols in Fig. 11b. The left hand side and right hand side plots correspond an angles of attack of 4∘ and 3∘. Sketches in the bottom-row plots indicate suction-side main features of the buffet mode

Full size image

A comparison of DMD modes at two different angles of attack (α = 3∘ and α = 4∘) suggests that modes corresponding to shock motion (St = 0.4 − 0.7) have increased amplitudes for higher angles of attack. For the α = 4∘ case, these modes contain features reminiscent of the buffet mode at St = 0.12 (marked by x-symbols in Fig. 12). These strong fluctuations near the trailing edge are significantly less pronounced at α = 3∘. This suggests a potential coupling or phase-locking of shock motion and transonic buffet at higher angles of attack, which would lead to the traditional observations of a single shock wave oscillating back and forth at a distinct buffet frequency corresponding to large-amplitude lift oscillations [18,19,20, 53, 54]. In the present cases (in particular the 3∘ case), there is no clear evidence on an interaction between shock waves and low-frequency lift oscillations. Instead, the shock motion appears as a consequence to acceleration and deceleration of the flow due to the buffet mode, as regions of high amplitudes in the DMD modes at St = 0.4 − 0.7 are located between regions of strong fluctuation in the buffet modes (marked by dotted lines and x-symbols in the bottom plots of Fig. 12). The following section aims to confirm the minor effect of the shock motion on the surface pressure.

7 Fourier Analysis of the Wall Density

After having gained a good impression of the global dynamics in the xy-plane, we now perform Fourier analysis along representative gridlines with z = 0 to study acoustic phenomena outside the boundary layers and, in particular, their foot-print on the wall density (which is proportional to the wall pressure for isothermal boundary conditions). We focus on the 3∘ case, as the amplitudes of large-scale fluctuations in the lift coefficient in the range of 27 < t < 52 (denoted by the dark red line in Fig, 10) are more regular compared to the 4∘ case. An instantaneous snapshot in Fig. 13a at t = 37.9 shows contours of the streamwise density gradient during the high-lift phase with the sonic line in green. Flow quantities are monitored on the airfoil surface, and along the dashed and dash-dotted lines, to generate similar space/time diagrams as in Fig. 3b for the 4∘ case. Contours of the streamwise density gradient are shown in Fig. 13b as a function of x and time, where the horizontal dashed black line corresponds to the dashed black curve in Fig. 13a. Similar plots are shown in Fig. 13c and d showing the density gradient as a function of x and time along the surface and along the dash-dotted line outside the supersonic region in Fig. 13a, respectively. The green lines in Fig. 13b–d correspond to the sonic lines as a function of x and time along the dashed black line in Fig. 13a and represent the shock dynamics. While the wall-density gradient in Fig. 13c represents a superposition of the buffet phenomenon, downstream-convecting turbulent structures, and traces of acoustic phenomena, Fig. 13d comprises mainly acoustic pressure waves circumventing the supersonic region (sketched by the black-dotted curves), which propagate upstream significantly faster than the shock waves indicated by the green lines.

Fig. 13
figure 13

Space/time diagrams corresponding to grid lines marked in a a 2D snapshot, which are located b just above the shear layer (dashed line in (a)), c on the surface, and d in the freestream (dash-dotted line in (a)). All plots show contours of streamwise pressure gradient. The green lines denote the sonic lines along a grid line corresponding to the dashed curve in Fig. 13a

Full size image

To get a clearer picture of the dynamics of the wall density, the same Fourier low-pass filter technique from Section 3 is applied to wall-density fluctuations for the 4∘ case. Filtered space/time diagrams at the surface are shown in Fig. 14 for a St < 0.2 and b 0.2 < St < 1. In Fig. 14a, large amplitudes in the low-frequency wall-density fluctuations on the suction side are observed at x ≈ 0.4 and at x ≈ 0.8, which agree with observations from the DMD analysis. We observe significantly lower amplitudes in regions near flow separation (0.4 < x < 0.7). Before and after that region, we can observe respectively the upstream- and downstream-propagation of peaks and troughs. The propagation speeds, however, are significantly slower than either the convection or acoustic velocities, which can be observed in Fig. 13b and d (sketched by dotted curves). While the fluctuations on both sides are in phase near the trailing edge in Fig. 14a, a clear phase shift can be observed in transition regions on the pressure and suction sides at x ≈ 0.8 and x ≈ 0.7, respectively. The surface-density fluctuations upstream of transition regions are in phase with the shear-layer dynamics, whereas the strong fluctuations near the trailing edge are 180∘ phase shifted, as shown in the DMD mode at St = 0.15 in Fig. 12. While the low-frequency fluctuations for x > 0.5 on the pressure side seem to be stationary, for x < 0.5 they propagate upstream at similar speeds as acoustic waves. A comparison of low-frequency oscillations on the suction-side surface with the shock motion, illustrated by the green lines in Fig. 14a, suggests that the global buffet mode sets the speed of upstream-propagating shock waves. During high-density phases (red contours), upstream-propagating shock waves are able to maintain subsonic regions between shock- and reflection-waves further upstream than during low-density phases (blue contours). The dotted curve in Fig. 14a sketches the potential path of a single shock wave oscillating between strong fluctuations at x ≈ 0.4 and x ≈ 0.8, always moving from alternating high- to low-density peaks. Whether this appears at higher Reynolds numbers requires further investigations.

Fig. 14
figure 14

Fourier band-filtered space/time diagrams of the surface density for frequency ranges between aSt < 0.2 and b 0.2 < St < 1. The green lines denote the sonic lines along a grid line corresponding to the dashed curve in Fig. 13a. The dotted curve in (a) suggests the potential path of a single shock wave

Full size image

Figure 14b, showing phenomena in the range of 0.2 < St < 1, suggests that the surface density is not significantly affected by shock waves. Even though upstream-propagating shock waves are observed in this frequency range, the green lines associated with the shock motion do not align with the contours of the wall-density gradient. Instead, we can observe traces of upstream-propagating acoustic waves circumventing the supersonic region at significantly higher speeds compared to the upstream-moving shock waves (as was shown in Fig. 13d). Towards the rear of the suction side, downstream-convecting structures appear in this frequency range, which are associated with large turbulent structures. These structures generate new acoustic waves when they arrive at the trailing edge.

The Fourier analysis in this section suggests that the shock motion is modified by the global buffet mode, as the irregular shock pattern (corresponding to St = 0.4 − 0.7) follows large-scale surface-density fluctuations associated with the St = 0.12 phenomenon.

8 Conclusion

DNS data has firstly been analysed in terms of local and global linear stability in order to investigate the transonic buffet mechanism of a narrow wing section at α = 4∘ and Mach and Reynolds numbers of M = 0.7 and Re = 500, 000, respectively. Local linear stability theory associates KH roll-ups seen in the DNS with convective linear instabilities, as has been previously observed for a high-pressure turbine vane at similar freestream conditions [30]. The shear layer on the suction side shows significantly different characteristics during high- and low-lift phases. The analysis of the time- and span-averaged flowfield underestimates the growth rates of those instabilities compared to phase-averaged flowfields for high-lift and low-lift phases. During the high-lift phase, the unstable region in the boundary layer is further upstream, with higher growth rates at higher frequencies, compared to the low-lift phase.

Global stability analysis captures a tonal mode at St = 5.89 that has been previously reported in literature with respect to the coupled dynamics of separation bubbles [42], but other modes fail to converge. Dynamic mode decomposition is able to capture the KH instabilities as well as the global mode at St = 5.89. Furthermore, the DMD shows an eigenmode at St = 1.87 that is associated with large turbulent vortices in the suction-side aft section of the airfoil. Those roll-ups seem to be coupled with the shock region and large-scale structures in the wake that are also observed in instantaneous visualisations. A phase shift in the streamwise velocity field of that mode occurs on the upper side between the wake and freestream. Modes at lower frequencies with 0.3 < St < 1 have similar structures with high amplitudes in similar regions and correspond to flapping of the separated shear layers and involve the shock motion as well as acoustic waves circumventing the supersonic region.

At a lower angle of attack, these spectral intermediate peaks are less pronounced. A DMD mode at St = 0.44 mainly describes the shock motion, whereas DMD modes at higher frequencies again comprise shock motion as well as acoustic waves. The existence of a bifurcation in the shock dynamics was already suggested by Giannelis et al. [44], where at certain flow conditions upstream-propagating pressure waves are periodically generated at the trailing edge at significantly higher frequencies than the co-existing permanent shock wave oscillating at a typical frequency associated with buffet (St ≈ 0.1). In the present case, no permanent shock wave is apparent so that only type C shock motion is observed, similar to experiments by McDevitt at al. [55]. The DMD modes associated with buffet at St = 0.12 − 0.15 show high fluctuations in the rear part of the airfoil and around x ≈ 0.4, but do not seem to be directly coupled with the shock motion. In contrast to the modes at high frequencies, the buffet DMD mode does not change significantly when decreasing the angle of attack.

The DMD buffet modes are reminiscent of the global modes reported by Crouch et al. [8] and Sartor et al. [46]. Propagation speeds of low-frequency structures do not correspond to acoustic speeds, suggesting that Lee’s acoustic feedback mechanism of buffet [7] is not active here. Instead it seems that the shock motion is modulated by the low-frequency buffet. The current observations suggest that transonic buffet does not rely on an interaction with shock waves, but does not exclude interactions and phase-locking between low-frequency phenomena and shock waves. Simulations at higher Reynolds numbers or higher angles of attack would be useful to see if there is a change to a single shock wave (visible at all times) oscillating between high- and low-pressure regions as sketched by the dotted curve in Fig. 14a.

Based on the conclusions of this work, it may be necessary to revisit the definition of transonic buffet to distinguish between transonic “shock buffet”, where a single shock wave oscillates at a frequency that agrees with low-frequency oscillations in the lift signal, and a form of “incipient buffet”, where the shock motion is not necessarily correlated with the dominant oscillations in the lift signal.