2.1. Flow configuration

Prior to the simulation, a series of experiments have been performed. In order to be consistent with the experimental set-up, the airfoil profile used in the simulations is the same as that used in the experiments. Readers are encouraged to refer to the companion work for more information on the experimental set-up (Alva et al. Reference Alva, Yuan, Araujo, Amaral, Hanifi and Cavalieri2023). Here, a brief summary of the geometry and flow conditions is presented.

The geometry investigated here is a modified NACA 0012 airfoil. The airfoil has a zero angle of attack and the free-stream Mach number is set to be $M = U_{\infty }/a_{\infty } = 0.3$ , where $U_{\infty }$ is the free-stream velocity and $a_{\infty }$ is the free-stream speed of sound. Two values of the Reynolds number, $0.8 \times 10^{5}$ and $1.0 \times 10^{5}$ , have been chosen to be investigated for potentially different flow phenomena. The chord length of the model is $c = 100\,\textrm{mm}$ . Unlike the theoretical NACA 0012 profile, the trailing edge is truncated at around 98 % of the chord length and rounded with an arc-circle of radius $0.04c$ , in order to control the trailing-edge shape, in agreement with earlier works of Ricciardi et al. (Reference Ricciardi, Wolf and Taira2022) and Ricciardi & Wolf (Reference Ricciardi and Wolf2022). The boundary-layer displacement thickness at the trailing-edge location is expected to be much thicker than the trailing-edge thickness. Thus, trailing-edge bluntness is not expected to be dominant for trailing-edge noise (Brooks et al. Reference Brooks, Pope and Marcolini1989).

In the present study, we consider two different geometries: one with a smooth surface and one with a row of spanwise-periodic roughness elements on both sides of the airfoil. The choice of roughness elements follows the work of Fransson et al. (Reference Fransson, Brandt, Talamelli and Cossu2005). The roughness elements have a height of $0.0055c$ (around 0.3 boundary-layer thickness) and a diameter of $0.015c$ . The distance between them is $0.06c$ . These roughness elements are placed close to the mid-chord position, at $x/c=0.52$ .

The extension of the computational domain in the spanwise direction is $0.12c$ , which contains two periods of roughness elements. A periodic boundary condition is used in the spanwise direction. According to the scattering condition, in order to generate propagating acoustic waves, the spanwise wavenumber has to satisfy $k_z \lt k_0$ , where $k_0$ is the acoustic wavenumber (Nogueira, Cavalieri & Jordan Reference Nogueira, Cavalieri and Jordan2017; Sano et al. Reference Sano, Abreu, Cavalieri and Wolf2019; Demange et al. Reference Demange, Yuan, Jekosch, Hanifi, Cavalieri, Sarradj, Kaiser and Oberleithner2024b ; Yuan et al. Reference Yuan, Demange, Jekosch, Sarradj, Oberleithner, Cavalieri and Hanifi2024). In this study, spanwise-coherent two-dimensional flow structures ( $k_z = 0$ ) dominate the generation of tonal noise. Therefore, the choice of the spanwise width is expected to have a reduced effect on the noise generation. A schematic drawing of the modified airfoil profile is shown in figure 1.

Figure 1. A schematic drawing of the modified NACA 0012 airfoil with roughness elements. The roughness elements are located at $x = 0.52c$ with a diameter of $D = 0.015c$ . The distance between the elements is $H = 0.06c$ . The red dashed line indicates the original NACA 0012 airfoil profile. The modified NACA 0012 profile is truncated at $x=0.98c$ and rounded with an arc of radius $r=0.04c$ .

2.2. Numerical methods

The open-source high-order flux-reconstruction numerical framework PyFR (Witherden, Farrington & Vincent Reference Witherden, Farrington and Vincent2014) is used for wall-resolved implicit LES (iLES) to solve the compressible Navier–Stokes equations in the general Cartesian coordinates with geometries stated above. Unlike LES with subgrid-scale models, iLES introduces numerical dissipation by discretisation. In this simulation, we use an anti-aliasing technique via the approximate $L^2$ projection, which has proven to be numerically stable and efficient. The principle behind this anti-aliasing is to compute the modal expansion coefficients of the desired polynomial exactly; see Park, Witherden & Vincent (Reference Park, Witherden and Vincent2017) for details. The discrete set of equations obtained from the spatial-spectral (hp)-type discretisation is integrated in time using a total variation diminishing third-order Runge–Kutta method. The simulations are performed using fourth-order polynomials.

In this work, the length scales, velocity components, density, pressure and frequency are non-dimensionalised as $x = x^* /c$ , $u = u^* /a_{\infty }$ , $\rho = \rho ^* /\rho _{\infty }$ , $p = p^*/(\gamma a_{\infty }^2 \rho _{\infty })$ and $St = f^*c/U_{\infty }$ , respectively. Here, $\rho _{\infty }$ is the free-stream density, $\gamma$ is the specific heat ratio and the quantities with the superscript $*$ are given in dimensional units. Therefore, the frequency $St$ presented here is a Strouhal number. Furthermore, Sutherland’s law is applied for the viscosity correction.

All simulations are performed on LUMI-G AMD MI250X 128 Gb GPU nodes. A compressible Navier–Stokes solver of fourth-order polynomials with hip-backend is used. The quadrature order used in the simulation is one order higher than the solver’s polynomial order to fit the requirement of iLES (Park et al. Reference Park, Witherden and Vincent2017). Interior solution points (i.e. Gauss–Legendre points for the hexahedral elements) are used for anti-aliasing purposes. Furthermore, at least 2800 snapshots (around 84 flow overs) are taken for all cases for spectral analysis. At least 30 flow overs are discarded before collection of snapshots. Around 4000 GPU hours are consumed for each case.

2.3. Mesh presentation

PyFR supports both structured and unstructured mesh topology. To take advantage of this, we used the following meshing strategy. A wall-resolved O-grid with maximum resolution at around mid-chord is given by $\Delta x^+ \lt 20$ , $\Delta z^+ \lt 10$ , $\Delta y^+ \lt 0.9$ , in a region where the flow is basically laminar. The $+$ superscript denotes viscous units. Considering that the flow on the airfoil is still transitional or fully developed to turbulence in the high-Re cases, at the trailing edge the resolution in the streamwise and wall-normal directions is kept as $\Delta x^+ \lt 4$ , $\Delta y^+ \lt 0.7$ . Near-wall mesh elements are lifted to second order to accommodate geometry curvatures. The structured mesh is also applied in the wake-resolved region until $1.5c$ after the trailing edge. The first layer of elements from the airfoil surface has an initial height of $\Delta y = 0.00055c$ and the growth rate of the structured layers is 1.08. For the cases with roughness elements, a prism type of mesh is applied in the near-roughness region to have a better fit of the roughness shapes and to avoid stretching of the mesh. The prism layer has much smaller elements with $\Delta s = 0.002c$ in the streamwise direction, which can reduce mesh-to-mesh discontinuity due to the discontinuous nature of the flux-reconstruction method.

In the acoustic region, unstructured grids are applied to reduce the large aspect ratio and significantly reduce the computational cost. Analysis of the simulation data shows that acoustic waves are still resolved up to Strouhal number $St = 15$ . The element aspect ratio in this region is kept below 5. Meshing details are shown in figure 2. The total grid number is around 70 million (polynomial order of 4). To allow a proper representation of the acoustic field and to have less interference from boundary conditions to the airfoil, the computational domain has been extended to 15 chords away from the airfoil surface. A designed sponge zone is added around the outer boundaries to absorb hydrodynamic wake and acoustic waves (Freund Reference Freund1997; Bodony Reference Bodony2006). The sponge is added as a forcing term $-\sigma (q-q_0)$ to the solver. Here, $q_0$ denotes free-stream quantities in the conservative format and the sponge strength parameter $\sigma = 0.008$ is carefully chosen such that flow solutions are damped to free-stream values. At the outflow boundaries we apply the non-reflective Riemann-invariant boundary conditions, which creates the outflow boundary conditions with free-stream values, together with the sponge zones, minimising reflections at the boundaries. At the inflow boundary, the Riemann-invariant condition is also used to provide a homogeneous potential flow with $M$ = 0.3. A periodic boundary condition is used in the spanwise direction. A non-slip adiabatic wall condition is applied to the airfoil surface.

Figure 2. High-order computational mesh elements (polynomial order of 4) near the airfoil (a) without roughness elements (clean case) and (b) with roughness elements (rough case). (c) Three-dimensional view of the wall surface mesh with roughness elements. Grid points inside the elements are skipped for a better visualisation. Structured grids are applied in the near-wall region, while unstructured grids are applied in the acoustic region.

2.4. Mesh validation

The mesh of the current simulation is designed targeting at the higher Reynolds number, i.e. Re = 1 00 000. The mesh has a distribution of points in the streamwise and spanwise directions on the airfoil surface given by $n_x = 830$ and $n_z = 130$ , respectively. From the airfoil surface to the height of the roughness elements ( $\triangle h = 0.0055c$ ), the distribution of points in the wall-normal direction is given by $n_y = 45$ . For clarification, we name this mesh as $M0$ which has a total of 4.86 million grid points in this near-wall region. A refined mesh $M1$ is designed as $n_x = 1185$ , $n_y = 50$ and $n_z = 170$ which leads to 10.07 million grid points in the near-wall region. The refinement is more concentrated from the roughness elements to the trailing-edge region where the flow gradients introduced by flow transition are much larger. Moreover, the aspect ratio is maintained. In the end, $M0$ leads to 70 million and $M1$ results in 80 million grid points in the computational domain.

Figure 3 illustrates a comparison between two mesh designs. Velocity profiles are acquired at the streamwise locations $x/c = 0.85$ and $x/c = 0.95$ . The velocity profiles show great convergence as overlapping for both mesh designs. The spectrum of the wall pressure signal is presented in figures 3(c) and 3(d) at the same locations and for the rough cases of both Reynolds numbers. Both mesh set-ups are able to capture the tonal peaks and secondary tones at Re = 80k and give the same broadband trend as for the high-Reynolds-number case. There are slight differences due to the fact that the time series acquired for the $M1$ mesh is shorter than that for the $M0$ mesh. For the case without roughness elements, the flow around the airfoil is expected to display weaker gradients, such that our mesh resolution $M0$ is more than sufficient. Note that both mesh set-ups use fourth-order polynomials and $h$ -type refinement is carried out. Since our cases are performed with iLES with anti-aliasing via projection acting on the surface flux (Park et al. Reference Park, Witherden and Vincent2017), $p$ -type refinement, i.e. increasing the polynomial order, will lead to a change in the filter as well as a change in the effect of the filter caused by the time scheme. Thus, $p$ -type refinement with the LES filter for grid convergence study is not considered in this study. The data computed from $M0$ are employed in the following data analysis.

Figure 3. Mesh refinement study. Mean flow profiles measured for roughness cases at $x/c = 0.85$ and $x/c = 0.95$ with (a) Re = 80k and (b) Re = 100k. The wall pressure spectrum obtained at $x/c = 0.85$ and $x/c = 0.95$ for (c) Re = 80k case and (d) Re = 100k case.

3.1. Acoustics

In order to compare with the experimental results, the power spectrum density (PSD) of the pressure fluctuations at the location $x/c = y/c = 1$ is illustrated in figure 4. The origin is taken at the leading edge of the airfoil, and this location is thus at 90° from the trailing edge. In the experiments, a beamforming technique is applied to compute the PSD of the pressure fluctuation, and a detailed description can be found in the experimental part of this work (Alva et al. Reference Alva, Yuan, Araujo, Amaral, Hanifi and Cavalieri2023). A Welch method (Welch Reference Welch1967) is used to calculate the PSD. For the LES data, the segment length of the Fourier transform is 512 and the overlap percentage is set to 75 %, which results in 22 blocks. The length of the time signal is 90 flow overs ( $c/U_{\infty }$ ) with sampling frequency $St = 33.3$ to ensure a well converged dataset. Results are presented here in terms of the sound pressure level, which is defined in decibels ( $\rm dB$ ) as

(3.1) \begin{equation} L_p = 10\log _{10}\frac {P_{xx}}{p^2_0} \end{equation}

with a reference sound pressure $p_0 = 2 \times 10^{-5}$ $\rm Pa$ ; we denote the PSD of the pressure fluctuations by $P_{xx}$ . Furthermore, the acoustic Mach number scaling (Williams & Hall Reference Ffowcs Williams and Hall1970; Brooks & Marcolini Reference Brooks and Marcolini1985), i.e. $M^5$ , is applied to the numerical dataset. This scaling helps to match the tonal acoustic amplitude from both the numerical and experimental datasets.

Figure 4. Acoustic spectra for (a) numerical simulations and (b) experiments measured in $\text{dB}/St$ . In (b), results of Re = 80k cases are shifted 55 dB lower for better visualisation. Tonal frequencies are labelled. Grey line: wind tunnel background noise.

For all numerical clean cases, the spectrum is dominated by a main frequency of $St = 5.5$ for the Re = 80k clean case and $St = 5.2$ for the Re = 100k clean case. As the Reynolds number increases, the main tone shifts to a slightly lower frequency but it has a stronger tonal noise amplitude. It is evident that secondary tones emerge in the spectrum at frequencies slightly lower and higher than the primary tone. As suggested by Tam (Reference Tam1974), the combination of the primary tone and the secondary tones is also referred to as ladder-like structures. With roughness elements being added to the airfoil surface, the main tone in the Re = 80k case has been broken down into multiple tones with much lower amplitude and wider frequency bands. Those tones have more or less similar amplitudes, and it becomes harder to distinguish a primary tone in the spectrum. In the Re = 100k rough case, the spectrum becomes broadband. Both the main tone and secondary tones are eliminated. The overall sound pressure level (OASPL) is calculated by integrating the spectrum in figure 4(a). The OASPL decreases from 89.1 to 85.6 $\rm dB$ in the case of Re = 80k and from 94.1 to 72.8 $\rm dB$ in the case of Re = 100k, from the clean to rough geometry, respectively. The OASPL results show that in the case of Re = 80k, the acoustic energy is redistributed to tones with lower amplitude, but with an overall decrease of about 3 $\rm dB$ . In the case of Re = 100k, the decrease of the main tones leads to a large overall reduction in noise.

Figure 4(b) presents the experimental measurement of the acoustic spectrum. A direct comparison between simulation and experiment is not possible, as is discussed later, but similar trends appear. For visualisation purposes, the results for Re = 80k cases are shifted 55 $\rm dB$ lower. For Re = 80k cases, both clean and rough surfaces can generate tonal noise. For the clean geometry the main tones have frequency $St = 3.2$ and $St = 4.1$ . The noise level for the rough case is lower and the peak frequency is shifted away from that of the clean case. More tones can be observed at higher frequencies, i.e. $St= 3.1$ , $3.9$ and $4.8$ . For the Re = 100k cases, very strong tonal noise can be observed at $St = 3.5$ , $4.0$ , $7.0$ and $8.0$ for the clean geometry, while most of the tones are eliminated from the spectrum with the presence of roughness elements, with the exception of tones at $St = 5.6$ and $6.2$ with much lower amplitudes.

Accounting for a Mach number scaling following $M^5$ (Williams & Hall Reference Ffowcs Williams and Hall1970; Brooks & Marcolini Reference Brooks and Marcolini1985), the tonal amplitudes between the numerical and experimental results are comparable, namely around 100 $\rm dB$ for case Re = 100k clean, 97 $\rm dB$ for case Re = 80k clean and 91 $\rm dB$ for case Re = 80k rough. The difference between the simulation and experimental results is within 3 $\rm dB$ . However, the tonal frequencies in the two datasets differ. This is likely for a number of reasons. First, in the simulation, free-stream turbulence (FST) level is assumed to be zero. However, this FST can lead to early transition to turbulence which eventually can affect the noise generation. Second, the experiments are carried out in a closed-section wind tunnel, unlike the free-field conditions used in the simulation. Finally, and perhaps more significantly, there is a Mach number difference between experiments and simulations (Mach = 0.3 in numerical simulations and Mach = 0.037 for experiments). In the experimental work of Pröbsting (Reference Pröbsting2022), the authors found that at zero incidence the tonal noise regime in the low-to-moderate Reynolds-number regime is sensitive to the variation in Mach number. The authors used the theoretical model from Arbey & Bataille (Reference Arbey and Bataille1983) and successfully predicted the tonal noise frequency intervals experimentally. In Appendix A, we employ the same analysis and find that both numerical and experimental datasets follow Arbey and Bataille’s model. Furthermore, we also perform the simulation with $M = 0.1$ and Re = 80k. The results are presented in Appendix B and the results show the relation of dominant tonal frequency and the Mach number, which follows the findings of Ricciardi, Arias-Ramirez & Wolf (Reference Ricciardi, Arias-Ramirez and Wolf2020) and Pröbsting (Reference Pröbsting2022). The influence of the FST is discussed in Appendix C which shows the limited influence on the results with low turbulence intensity (TI = 0.16 %) as in the experiments. In this work, we did not intend to match the experiments exactly due to the cost of the simulations. However, it is already encouraging that the same qualitative trends related to the use of roughness elements are observed in both simulation and experiment.

Sample snapshots of the two-dimensional acoustic fields of all four cases are illustrated in figure 5 with pressure fluctuation $p^{\prime}$ . In the clean cases, strong acoustic radiation can be observed from the trailing-edge region. From lower to higher Reynolds number, the amplitude of the acoustic waves becomes stronger. With the presence of roughness elements, the amplitude of acoustic waves is lower and the wavelength is changed compared with the clean cases. For Re = 100k rough case, the acoustic radiation is much lower than for other cases and it is hard to identify under the same level range. These observations are consistent with the results of the acoustic spectra in figure 4.

Figure 5. Acoustic field contour plots for pressure fluctuation $p^{\prime}$ . The airfoil leading-edge point is at $x/c = y/c = 0$ . (a) Re = 80k clean, (b) Re = 100k clean, (c) Re = 80k rough and (d) Re = 100k rough.

The directivity of cases Re = 80k clean/rough and Re = 100k clean are presented in figure 6. The directivity is obtained along an arc-circle with the origin at the trailing edge and radius of $2c$ , and wake region is excluded. The directivity corresponding to the tonal frequencies is presented. Similar to Demange et al. (Reference Demange, Yuan, Jekosch, Hanifi, Cavalieri, Sarradj, Kaiser and Oberleithner2024a ), the directivity of the tonal frequencies shows as a modified cardioid pattern, arising from interferences between trailing- and leading-edge scattering (Howe Reference Howe2001; Roger & Moreau Reference Roger and Moreau2005).

Figure 6. Directivity plots of the acoustic field for (a) Re = 80k clean, $St = 5.3$ (), (b) Re = 100k clean, $St = 5.1$ () and (c) Re = 80k rough, $St = 3.4$ ( ), $St = 4.5$ (), $St = 5.5$ (), $St = 6.5$ ().

3.2. Flow structures

As an overview of the flow structures, the visualisation of the flow field is presented in figure 7 in terms of the Q-criterion (Jeong & Hussain Reference Jeong and Hussain1995) coloured by the streamwise velocity component $u$ . The Q-criterion is defined as $1/2 ( ||R||^2 - ||S||^2)$ , where $R$ and $S$ are the rotation and strain rate tensors, respectively. To enhance visualisation and conserve computational resources, a boxed region covering the upper side of the airfoil is utilised for calculating the Q-criterion. The Q-criterion can help us to identify the flow states and dominant structures on the airfoil. The Q-criterion values used in the visualisations for the clean cases are one order of magnitude higher than those for the rough cases. Videos corresponding to four cases can be seen as supplementary movies available at https://doi.org/10.1017/jfm.2025.10321.

Figure 7. Iso-surface of the Q-criterion coloured by streamwise velocity $u$ .

For clean cases, the dominant structures are K–H-type rollers. These spanwise-coherent structures are responsible for the tonal noise generation at the trailing edge as also identified in the works of Pröbsting et al. (Reference Pröbsting, Scarano and Morris2015) and Ricciardi et al. (Reference Ricciardi, Wolf and Taira2022). With roughness elements being added to the airfoil surface for both Reynolds numbers, streaky structures are generated by these roughness elements. These streaks will develop into lambda structures (also called hairpins) while propagating downstream and modulating K–H rollers. At the lower Re, K–H-type rollers are modulated but maintain a certain level of spanwise coherence. An earlier breakdown of the K–H rollers to three-dimensional structures is spotted close to the trailing edge. With an increase of Re, the flow structure on the airfoil surface becomes dominated by three-dimensional lambda structures. The trace of spanwise-coherent structures is very weak. From these observations, we can infer that these structures may link the acoustic field presented in figure 4 and the hydrodynamic field. Tonal noise is generated by a feedback mechanism involving the spanwise-coherent K–H rollers. The streak-modulated K–H rollers have lower spanwise coherence, which ultimately decreases the acoustic amplitude and shifts the tonal frequencies. Stronger streaks at higher Re could lead to an earlier breakdown to three-dimensional structures which will generate broadband type of noise rather than tonal noise. Further analysis of these aspects is presented in § 4.

Figure 8. Skin friction coefficient ( $C_f$ ) of span-averaged field of the clean cases at both Reynolds numbers: (a) Re = 80k; (b) Re = 100k. Solid red line and dashed black line represent upper and lower sides of the airfoil, respectively. Blue dashed line indicates zero values.

In the numerical work of Ricciardi et al. (Reference Ricciardi, Wolf and Taira2022), separation bubbles play an important role, as they lead to strong K–H instabilities. In order to reveal the changes of the separation bubble induced by streaks, the span- and time-averaged skin friction coefficients ( $C_f$ ) for clean cases are shown in figure 8. Here the skin friction coefficient is defined as

(3.2) \begin{equation} C_f = \frac{\tau_w}{0.5\rho_{\infty} U_{\infty}^2}, \end{equation}

where $\tau _w$ represents the wall shear stress and subscript $\infty$ refers to the free-stream quantities. For the clean cases, a long separation bubble can be spotted on the both sides of airfoil. At Re = 80k, the separation bubble initiates at $x/c = 0.65$ and at $x/c = 0.67$ for the higher-Re case, extending up to the trailing edge. A significant difference of $C_f$ at the two sides of airfoil close to the trailing edge can be seen for the higher-Re case. This phenomenon has also been documented in the work of Ricciardi et al. (Reference Ricciardi, Arias-Ramirez and Wolf2020) where a series of two-dimensional simulations were conducted at low Reynolds numbers and Mach numbers. At Re = 100k and Mach = 0.3, the dynamics displays symmetry breaking, with a stronger separation bubble at one of the airfoil sides which occasionally switches to the other side. Apparently, our dataset does not cover such a long time evolution, which could be the reason for the difference of $C_f$ between the two sides of the airfoil. Though this could slightly affect the tonal noise characteristics, it does not change the fact that streaks can attenuate the tonal noise level. And as discussed in the following analysis, with the presence of the streaks, the $C_f$ distributions of both sides of the airfoil are identical. So in order to avoid long and expensive simulations to have fully converged statistics, we keep the current dataset for further analysis.

Figure 9. (a,b) Skin friction coefficient ( $C_f$ ) heat maps of mean flow field of roughness cases. Red is positive and blue is negative. Dashed grey line is zero-contour line. (c) Oil visualisation from the experimental campaign at Re = 1 50 000. Red dashed line: the zero-contour line from Re = 100k simulation projected on the oil visualisation.

For the cases with roughness elements, the structures on the airfoil are three-dimensional. Hence, spanwise-averaged statistics are no longer appropriate. Thus, to understand how the existence of structures generated by the roughness elements modifies the surface flow field, the time-averaged surface skin friction coefficient heat maps are presented in figure 9(a,b). The colour map is saturated such that blue indicates negative and red represents positive values. The grey dashed line indicates the zero-contour line. A re-circulation region can be spotted prior to and after the cylinder roughness in both cases. Downstream of the location of the roughness elements, a long separation bubble (characterised by negative $C_f$ regions) can be spotted in both cases. However, the separation bubble has been sliced/modulated into subdomains by the streamwise-elongated structures highlighted by the positive $C_f$ values. Fast streaks (indicated by regions of positive $C_f$ ) are induced downstream of roughness elements, in agreement with earlier studies (Fransson et al. Reference Fransson, Brandt, Talamelli and Cossu2004, Reference Fransson, Brandt, Talamelli and Cossu2005). At lower Reynolds numbers, streaks have significantly lower amplitudes. The initial locations of separation bubbles are roughly equivalent to those in the clean case, yet the re-attachment location is further upstream compared with the clean case. For the high-Re case, streaks have higher amplitude immediately after the roughness elements, maintaining it until the trailing-edge region. The end of the separation bubble moves upstream, to around $x/c = 0.86$ .

Figure 9(c) illustrates flow visualisation from the experiment at Re = 1 50 000 in which the surface oil flow visualisation technique is employed to reveal flow patterns, with technical details provided in the companion paper. It should be noted that a higher-Reynolds-number case is presented to enhance the visualisation effect, as lower-speed flows did not lead to clear visualisations in the wind tunnel. The overall mechanism remains the same: this visualisation method captures the fast and slow streaks directly behind the roughness elements, with streak traces extending up to around $x/c = 0.9$ . This finding is highly consistent with the $C_f$ distribution from the simulation depicted in the same figure. The zero-velocity contour line from the simulation Re = 100k case is projected to the oil visualisation and, clearly, the location and width of the streaks are very similar in both cases.

4.1. Spectral proper orthogonal decomposition

The previous visualisations are of steady structures induced by the roughness elements. We now characterise fluctuating fields, which are expected to be directly related to acoustic radiation, using modal analysis. The POD technique is a data-driven method which extracts a set of orthogonal basis functions from flow realisations (Lumley Reference Lumley1970; Berkooz, Holmes & Lumley Reference Berkooz, Holmes and Lumley1993) that maximise the mean square energy, considering an appropriate inner product. In SPOD (Picard & Delville Reference Picard, Delville, Rodi and Laurence1999; Towne, Schmidt & Colonius Reference Towne, Schmidt and Colonius2018), this basis is defined by the eigenvalue decomposition of the cross spectral density matrix of the Fourier-transformed realisations at each frequency, ensuring spatio-temporal coherence of the structures. We consider the standard Reynolds decomposition of the flow variables, i.e. $q = \bar {q}+q^{\prime}$ , where $\bar {q} = [\overline {\rho },\overline {u},\overline {v},\overline {w},\overline {T}]^{\text{T}}$ stands for the time average and $q^{\prime}$ is the fluctuation component. Here in this work, we define the fluid component as $q^{\prime} = [\rho ^{\prime},u^{\prime},v^{\prime},w^{\prime},T^{\prime}]^{\text{T}}$ as the density, tangential, wall-normal and spanwise velocity components and temperature fluctuations, respectively. Since SPOD is applied to compressible flow components at the streamwise cut locations, the norm used here is the compressible energy (Chu Reference Chu1965; Mack Reference Mack1984; Hanifi, Schmid & Henningson Reference Hanifi, Schmid and Henningson1996). Following the so-called snapshot method (Schmidt & Colonius Reference Schmidt and Colonius2020), firstly a windowed fast Fourier transform is performed on the fluid field data in time to obtain the flow field $\hat {q}$ (where a hat denotes a Fourier-transformed quantity) at each discrete frequency. The decomposition is performed at each frequency based on the inner product between different flow realisations $\hat {q}_{{i}}$ , given by

(4.1) \begin{equation} \langle \hat {q}_{{i}}, \hat {q}_{{j}}\rangle = \int _{\Omega }\hat {q}_{{j}}^H\mathcal{W}\hat {q}_{{i}}\, \text{d}\textbf{x}= \hat {q}_{{j}}^H\textbf {W}\hat {q}_{{i}}{,} \end{equation}

where $\Omega$ is the region of interest, the superscript $H$ denotes the Hermitian transpose and the discretised weighting operator, $\textbf {W}$ , is chosen such that the sum of the eigenvalues defines a compressible energy (related to the variance of the cross spectral density). Following the original derivation from Chu (Reference Chu1965), $\textbf {W}$ is a positive-definite weighting tensor:

(4.2) \begin{equation} \textbf {W} =\int _{\Omega }\begin{bmatrix} \frac {\overline {T}}{\gamma \overline {\rho }M^2} & & & &\\ & \overline {\rho } & & &\\ & &\overline {\rho } & &\\ & & &\overline {\rho } &\\ & & & & \frac {\overline {\rho }}{\gamma (\gamma -1)\overline {T}M^2} \end{bmatrix}\text{d}\textbf{x}, \end{equation}

where $\gamma$ is the specific heat ratio and $M$ is the Mach number.

In all the analyses studied here, we have used Hamming windows for the Welch method, since they are commonly used in narrowband applications. The analysis is performed with the numerical code from Schmidt & Colonius (Reference Schmidt and Colonius2020), using a short-time fast Fourier transform considering $N_{{FFT}}=128$ frequency bins per block with an overlap of $75\,\%$ , resulting in a number of blocks $N_{{blocks}}=172$ and $\Delta T = 0.03c/U_{\infty }$ . In order to increase the resolution of SPOD analysis, data from the two periods of roughness elements on one side of the airfoil are treated as two independent realisations. As the case with roughness elements does not allow a Fourier transform along the spanwise direction, SPOD modes have a non-trivial z dependency. To limit the data requirements, we apply SPOD separately at the streamwise positions illustrated in figure 10.

The energy spectrum ( $\gamma _i$ ) calculated at the station $x/c = 0.75$ is presented in figure 11 for both Reynolds numbers. At lower Re, a group of noticeable peaks can be observed in the spectrum at $St = 5.5$ , $4.4$ and $6.5$ . Those peaks in the spectrum are related to the tonal noise illustrated in the figure 4. As expected, for the higher Re, there are no clear peaks in the spectrum and no dominant frequency can be detected. In order to investigate the effects of roughness elements, we choose the frequency $St = 5.2$ in the Re = 100k case which corresponds to the tonal noise in the clean case. It is worthwhile to note that from the spectrum, the energy of the leading mode ( $\gamma _1$ ) of chosen frequency is around one order of magnitude higher than that of the second mode ( $\gamma _2$ ) at the frequencies of interest, which indicates that the system dynamics is of low rank, i.e. the flow is dominated by the first mode. The leading SPOD mode can thus be used as an effective representation for flow structures.

Figure 11. The SPOD energy ( $\gamma _i$ ) calculated at the station $x/c = 0.75$ at Re = 80k (a) and Re = 100k (b). Data correspond to the rough cases.

Figure 12 shows the absolute values of the leading SPOD mode of the wall-tangential velocities $\hat {u}$ at the selected streamwise stations. Note that SPOD modes have a unit norm, so to see the growth of structures at streamwise stations, SPOD modes are multiplied with the square root of the eigenvalue to show the actual fluctuation amplitudes in the flow. From the first and second stations ( $x/c = 0.55$ and $0.65$ ), we see the fast growth of the amplitude of unsteady streak-type structures with a change of almost two orders of magnitude for higher-Re case. At the station $x/c = 0.75$ , for the lower-Re case, the growth of spanwise-modulated K–H rollers can be identified. Those structures sit at the location $|z-z^{\prime}|/c \geqslant 0.015$ . For the higher Re, similar patterns are observed. Although the spanwise coherence is lower for Re = 100k, the structures display patterns similar to those of the Re = 80k case, suggesting that, while of difficult identification in the Q-criterion snapshots, strongly modulated K–H rollers remain even for higher Reynolds numbers. From station $x/c = 0.85$ onward ( $x/c \geqslant 0.85$ ), the amplitude of structures begins to decrease for the higher-Re case and continues to increase for the lower-Re case. With the growth of the boundary layer, flow structures start to spread out, marking the ongoing process of the transition to turbulence at these stations. Therefore, in the subsequent discussions, we primarily focus on the lower-Re scenario at around station $x/c = 0.75$ .

Figure 12. (a,b) Absolute value of the leading SPOD mode of streamwise velocity $\tilde {u}$ . From left to right: stations $x/c = 0.55$ , $0.65$ , $0.75$ , $0.85$ and $0.95$ .

Figure 13 depicts the absolute values of the leading and secondary SPOD modes at $x/c = 0.75$ for Re = 80k. At this station, streaks are fully developed and the amplitudes of both streaks and K–H rolls have reached the same order of magnitude. From the figure, the leading SPOD mode describes a streak-modulated K–H roll, while the second mode contains only thinner localised structure akin to streak fluctuations. However, at this frequency, the flow field is dominated by the leading SPOD mode according to the energy ranking presented in figure 11. Hence, the leading SPOD mode is a good representation of the flow state. The most energetic structures are located $0.01c$ away from the airfoil surface.

Figure 13. The absolute value of the leading (a) and secondary (b) SPOD modes calculated at the station $x/c = 0.75$ for Re = 80k. From left to right: velocity components $\hat {u}$ , $\hat {v}$ and $\hat {w}$ .

4.2. Local spatial stability analysis

In order to investigate how the streaks attenuate growth of the K–H rolls, a two-dimensional spatial–local stability analysis in the $z{-}y$ plane is performed at the streamwise station $x/c = 0.75$ . Here, we are interested in the analysis of the hydrodynamic part of the flow field, and because of the low Mach number of the flow, the incompressible linearised Navier–Stokes formulation is sufficient for the analysis:

(4.3) \begin{equation} \begin{gathered} \frac {\partial \textbf{u}^{\prime}}{\partial t} + (\bar {\textbf{U}}\cdot \nabla)\textbf{u}^{\prime} + (\textbf{u}^{\prime}\cdot \nabla)\bar {\textbf{U}} = -\nabla p^{\prime} + \frac {1}{Re}\nabla ^2\textbf{u}^{\prime}, \\ \nabla \cdot \textbf{u}^{\prime} = 0, \end{gathered} \end{equation}

where $\bar {\textbf{U}}$ indicates the mean flow, upon which linearisation is performed and variables with primes indicate perturbation around the mean flow $\bar {\textbf{U}}$ .

We carry out a stability analysis in the local setting, neglecting x-variation of the mean flow. Hence, the ansatz for the perturbation has the form

(4.4) \begin{equation} \{u^{\prime},v^{\prime},w^{\prime},p^{\prime}\} = \{\hat {u},\hat {v},\hat {w},\hat {p}\}(y,z)\text{e}^{\text{i}(\alpha x - \omega t)}. \end{equation}

Substituting the ansatz (4.4) into the linearised Navier–Stokes equation (4.3) with assumption of slow variation in the streamwise direction, the following stability problem is established:

(4.5) \begin{equation} \begin{gathered} i\alpha \hat {u} + \frac {\partial \hat {v}}{\partial y} + \frac {\partial \hat {w}}{\partial z} = 0, \\ -\bar {U}\alpha \hat {u} + \xi \hat {u} + i\hat {v}\frac {\partial \bar {U}}{\partial y} + i\hat {w}\frac {\partial \bar {U}}{\partial z} - \alpha \hat {p} + \alpha ^2\frac {i}{Re}\hat {u} = -\omega \hat {u}, \\ -\bar {U}\alpha \hat {v} + \xi \hat {v} + i\hat {v}\frac {\partial \bar {V}}{\partial y} + i\hat {w}\frac {\partial \bar {V}}{\partial z} - \frac {\hat {p}}{y} + \alpha ^2\frac {i}{Re}\hat {v} = -\omega \hat {v}, \\ -\bar {U}\alpha \hat {w} + \xi \hat {w} + i\hat {v}\frac {\partial \bar {W}}{\partial y} + i\hat {w}\frac {\partial \bar {W}}{\partial z} - \frac {\hat {p}}{z} + \alpha ^2\frac {i}{Re}\hat {w} = -\omega \hat {w}, \\ \xi = i\bar {V}\frac {\partial }{\partial y} + i\bar {W}\frac {\partial }{\partial z} - \frac {i}{Re}\left(\frac {\partial ^2}{\partial y^2} + \frac {\partial ^2}{\partial z^2}\right). \end{gathered} \end{equation}

Letting $\textbf{q} = \{\hat {u},\hat {v},\hat {w},\hat {p}\}^T$ , the equation system can be rearranged to

(4.6) \begin{equation} \left (\alpha \mathcal{A}_1 + \alpha ^2\mathcal{A}_2 \right)\textbf{q} = \left (\omega \mathcal{B} - \mathcal{A}_0 \right)\textbf{q}, \end{equation}

where $\mathcal{A}_0,\mathcal{A}_1,\mathcal{A}_2$ and $\mathcal{B}$ are functions of the mean flow and flow parameters. Following the spatial stability ansatz, we take $\omega \in \mathcal{R}$ to be the temporal frequency input while $\alpha \in \mathcal{C}$ is the complex output eigenvalue. Notice that a nonlinear term $\alpha ^2$ appears in the formulation. Thus, a linear eigenvalue problem, considering an eigenfunction given by $[\textbf{q}, \alpha \textbf{q}]$ , is solved instead (Bridges & Morris Reference Bridges and Morris1984; Piot, Casalis & Rist Reference Piot, Casalis and Rist2008).

The effect of roughness-generated streaks is investigated through tracking the unstable modes from the mean state without roughness elements to that with roughness elements (Lajús et al. Reference Lajús, Sinha, Cavalieri, Deschamps and Colonius2019) by a continuous variation of the mean flow ( $\bar {\textbf{U}}$ ) as

(4.7) \begin{align} \bar {\textbf{U}} = (1 - \sigma) {\bar{\textbf{U}}_{clean}} + \sigma {\bar{\textbf{U}}_{rough}}, \end{align}

where $0\leqslant \sigma \leqslant 1$ . The tracking mode technique aids in identifying the growth rate changes between two distinct states. Ideally, a continuous change of the parameter $\sigma$ allows us to trace the progression of an unstable mode through these states. For the sake of reducing computational effort, we utilise five discrete values of $\sigma$ , which are sufficient to demonstrate the evolution of the eigenvalues. The spectrum is illustrated in figure 14. As can be seen there, there are three unstable modes in the clean case. The most unstable mode stands for two-dimensional K–H instability, with the mode shape depicted in the right-hand subplot of figure 15. The other two modes have identical low growth rate and streamwise wavenumber, corresponding to degenerate oblique K–H modes with positive and negative spanwise wavenumbers, as in Lajús et al. (Reference Lajús, Sinha, Cavalieri, Deschamps and Colonius2019). With an increase in the value of the parameter $\sigma$ , the growth rate of the most unstable mode is decreasing from $\alpha _i = 0.1699$ to $\alpha _i = 0.1231$ . However, the other two unstable modes in the clean state diverge and become more unstable with increasing $\sigma$ . Two other modes appear from the stable regime, becoming unstable as $\sigma$ approaches 1.

Figure 14. The spectrum traced with a gradual increase of the mean state modulation parameter $\sigma$ . Open circles stand for the unstable modes and open squares represent stable modes.

The shapes of the most unstable mode for various values of $\sigma$ are presented in figure 15. Note that the eigenvectors are normalised with respect to the maximum value of streamwise velocity component $\tilde {u}$ at each $\sigma$ value. From the clean state ( $\sigma = 0$ ) to the rough state ( $\sigma = 1$ ), a progressive spanwise modulation appears at the middle of domain. The K–H mode is heavily modulated in this region. The black dashed line denotes the critical layer, where $U=\omega /\alpha$ , which is the location where the modulated modes peak. The appearance of streaks increases the amplitude of the spanwise velocity component $\tilde {w}$ while it decreases the strength of the wall-normal velocity component $\tilde {v}$ . Comparing the eigenmode for $\sigma =1$ with the leading SPOD mode at the same station, the modes show many similarities, indicating that the SPOD mode is indeed representative of a spanwise-modulated K–H roller. We also extend the stability analysis to the upstream station $x/c = 0.70$ for both Reynolds numbers, as presented in Appendix D. The results suggest the same trend that the growth rate of the modulated K–H mode decreases with the increasing modulation of the mean flow made by streaks.

Figure 15. The leading eigenmodes (K–H modes) traced with gradual increase of mean state modulation parameter $\sigma$ . Here $\sigma$ changes from 1 to 0 from left to right. (a), (b) and (c) rows present the velocity components $\tilde {u}$ , $\tilde {v}$ and $\tilde {w}$ , respectively. Black dashed line indicates the critical layer.

4.3. Spectral proper orthogonal decomposition versus local stability analysis

From the spatial stability analysis, it can be concluded that streaks can be used to decrease the growth of K–H instabilities. However, we are using the mean flow rather than the base state of the flow due to the difficulty of getting a laminar solution in the present setting. In order to verify this result, one can calculate the growth of hydrodynamic waves from the leading SPOD mode. Nonetheless, due to the existence of strong acoustic waves in the field, SPOD modes contain energy of the hydrodynamic and acoustic perturbations. In order to extract hydrodynamic energy from the total energy, we utilise the bi-orthogonality of state and adjoint eigenfunctions (see e.g. Tumin Reference Tumin2011; Shahriari, Hanifi & Henningson Reference Shahriari, Hanifi and Henningson2016). In order to use our existing tool, we consider the problem in the framework of temporal stability analysis. Consider the direct and adjoint problem in a compact form as

(4.8) \begin{equation} \textbf{N} \textbf{q} = \omega \textbf{M} \textbf{q}, \end{equation}

(4.9) \begin{equation} \textbf{N}^+ \textbf{q}^+ = \omega ^+ \textbf{M}^+ \textbf{q}^+,\end{equation}

where superscript $+$ indicates the adjoint and $\omega ^+=\text{conjg}(\omega)$ . Let us define the weighted inner product with operator M as the weight operator as

(4.10) \begin{equation} \big \langle \textbf{a},\textbf{b}\big \rangle _{\textbf{M}} = \big \langle \textbf{M}\textbf{a},\textbf{b}\big \rangle = \int _{z_{min}}^{z_{max}}\int _0^{\infty } \textbf{b}^H\textbf{M}\textbf{a}\ \text{d}y\ \text{d}z, \end{equation}

where $H$ denotes complex conjugate transpose. Then the bi-orthogonality yields

(4.11) \begin{equation} \big \langle \textbf{q}_i^{\ },\textbf{q}^+_j\big \rangle _{\textbf{M}}=\mathcal{E}\delta _{ij}, \end{equation}

where $\delta _{ij}$ is the Kronecker delta function and $\mathcal{E}$ is an arbitrary constant. Therefore, by projecting the LES solution onto the adjoint of hydrodynamic modes, one extracts their amplitude (Rodríguez et al. Reference Rodríguez, Cavalieri, Colonius and Jordan2015). The amplitude of the hydrodynamic wave, $C_{h}$ , at each xlocation can be evaluated through

(4.12) \begin{equation} C_{h} = \frac {\big \langle \hat {\boldsymbol{\Phi }}, \textbf{q}_{h}^+\big \rangle _{\textbf{M}}}{\big \langle \textbf{q}_{h}^{\ }, \textbf{q}_{h}^+\big \rangle _{\textbf{M}}}.\end{equation}

Here, $\textbf{q}_{h}$ is the eigenvector corresponding to the hydrodynamic mode and $\hat {\boldsymbol{\Phi }}$ = $\langle \hat {\rho }_s,\hat {u}_s,\hat {v}_s,\hat {w}_s,\hat {T}_s\rangle$ is the vector containing the SPOD mode. The denominator in (4.12) can be considered as a normalisation constant for the adjoint eigenfunction at each station. It is also important to note that the computed hydrodynamic mode amplitudes depend on the normalisation of the direct hydrodynamic eigenmode. Here, we have normalised hydrodynamic eigenmode, $\textbf{q}_{h}$ , with maximum streamwise velocity.

The eigenvalues of the direct and adjoint analysis are illustrated in figure 16 at $x/c = 0.75$ . The complex-conjugated adjoint eigenvalues overlap with the direct eigenvalues which indicates the convergence of mesh and resolution. The frequencies of the unstable modes correspond to the tonal frequency and a modulated K–H mode, $St = 5.5$ , and possibly to a secondary instability of streaks, $St = 3.9$ . The mode corresponding to the modulated K–H mode is shown in figure 17. The direct mode is closely aligned with the spatial and SPOD analyses. Furthermore, the adjoint mode reveals the sensitivity of flow, and the least unstable mode indicates sensitivity of the spanwise-modulated K–H instability.

Figure 16. Eigenvalues of direct and adjoint modes calculated at the station $x/c = 0.75$ . Note that adjoint eigenvalues are complex-conjugated.

Figure 17. The absolute value of the modulated K–H mode calculated at the station $x/c = 0.75$ . From (a) to (c) and (d) to (f): velocity components $\tilde {u}$ , $\tilde {v}$ and $\tilde {w}$ . (a–c): direct mode; (d–f): adjoint mode.

Five close stations $x/c = 0.75, 0.75 \pm 4\times 10^{-4}, 0.75 \pm 8\times 10^{-4}$ are used to calculate the spatial growth of K–H modes around the station $x/c = 0.75$ . For this purpose, first, the amplitude of the dominant K–H mode is extracted at these stations using the projection given in (4.12). Then the maximum absolute value of the spanwise-averaged velocity field (to extract the two-dimensional K–H mode) is used to calculate the growth rate as

(4.13) \begin{equation} \text{Growth rate}=\frac {1}{|u|_{max}}\frac {\text{d} |u|_{max}}{\text{d} x}. \end{equation}

The results are listed in table 1. As can be seen there, the projected growth rates of the K–H mode are very close to the prediction given by the spatial stability analysis. This confirms the proposed mechanism for reduction of tonal noise: the introduction of roughness elements induces steady streaks, which render the laminar separation bubble three-dimensional, creating a spanwise modulation of the K–H instability of the shear layer, with a significant reduction of the growth rate. This weakens one of the components of the feedback loop leading to tonal airfoil noise, reducing tones emitted to the far field.

Table 1. Growth rate at station $x/c = 0.75$ from spatial stability analysis and projected SPOD modes via (4.12) for both clean ( $\sigma = 0$ ) and rough ( $\sigma = 1$ ) states.