2.1. NACA 0012 and NACA 4412 airfoils

The NACA 0012 and NACA 4412 airfoil databases analysed in this work are described in Tanarro et al. (Reference Tanarro, Vinuesa and Schlatter2020). In this study, a well-resolved LES was performed to analyse the incompressible turbulent flow around the airfoils at a chord Reynolds number of $Re_{c} = 400\,000$ for both wing sections, with $0^{\circ }$ and $5^{\circ }$ angles of attack for the NACA 0012 and NACA 4412, respectively. The effect of the unresolved turbulent scales is modelled by a subgrid-scale model based on a relaxation-term approach (Schlatter, Stolz & Kleiser Reference Schlatter, Stolz and Kleiser2004), as validated in e.g. Eitel-Amor, Örlü & Schlatter (Reference Eitel-Amor, Örlü and Schlatter2014). An unsteady volume force was used to force transition to turbulence, at 10 % of the chord on both the top and bottom sides of the wing, for both simulations. The incompressible spectral-element Navier–Stokes solver Nek5000 (Fischer, Lottes & Kerkemeier Reference Fischer, Lottes and Kerkemeier2008) was used to carry out the simulations.

For both cases, the spatial discretisation is performed by means of Lagrange interpolants of polynomial order $N = 11$, where the points within the element are distributed as tensor products of the Gauss–Lobatto–Legendre (GLL) quadrature rule. The spatial resolution in the boundary layer is designed to reach a target resolution in wall units of ${\rm \Delta} x_{t}^{+}=18$, ${\rm \Delta} y_{n}^{+}=0.64$–$11$ and ${\rm \Delta} z^{+}=9$ in the middle of the airfoil. Here, $x_t$ and $y_n$ indicate tangential and normal directions to the wing surface, respectively, and $z$ the spanwise direction. The scaling in wall units, denoted by the $+$ superscript, is in terms of the friction velocity $u_\tau$ and kinematic viscosity $\nu$. The domain considered here is a C-mesh with streamwise length $L_x = 6c$, wall-normal length $L_y = 4c$ and periodic spanwise length $L_z = 0.1c$, where $c$ is the wing chord (see figure 1). Note that $x$ and $y$ in figure 1 denote the directions parallel and perpendicular to the chordline.

Figure 1.

Illustration of the two-dimensional domain for both airfoils, showing the spectral-element mesh in grey, without the GLL points, and the instantaneous streamwise velocity field. The airfoil NACA 0012 is at zero angle of attack, and the NACA 4412 is at $5^{\circ }$. The domain considered here is a C-mesh with streamwise length $L_x = 6c$, wall-normal length $L_y = 4c$ and periodic spanwise length $Lz = 0.1c$, where $c$ is the wing chord. (a) NACA 0012. (b) NACA 4412.

A Reynolds-averaged Navier–Stokes (RANS) simulation with a largely extended domain was used to obtain the velocity distributions used as boundary conditions in the LES. The domain size of the NACA 0012 profile is discretised using a total of $220\,000$ spectral elements, which amounts to approximately $380$ million grid points, whereas the mesh of the NACA 4412 profile has $270\,000$ spectral elements, giving a total of $466$ million grid points. Figure 1 shows the domain for both airfoils, showing the spectral-element mesh, without the GLL points, and the instantaneous streamwise velocity field. See Tanarro et al. (Reference Tanarro, Vinuesa and Schlatter2020) for more details about the airfoil databases.

We extracted from the simulations of the NACA 0012 and NACA 4412 airfoils $515$ and $801$ snapshots of the full three-dimensional flow fields, respectively, separated in time by ${\rm \Delta} t_{c} = 0.01$, a normalised time step using the free-stream velocity $U_\infty$ and the wing chord $c$, for both cases. Table 1 shows the momentum-thickness Reynolds numbers, $Re_{\theta }$, and the respective timesteps between snapshots scaled with the momentum thickness, ${\rm \Delta} t_{\theta }$, for both airfoils at stations $0.50 \leq x/c \leq 0.90$; for the NACA 4412 airfoil, the suction side is considered.

Table 1.

Momentum-thickness Reynolds number, $Re_{\theta }$, and the respective timestep, ${\rm \Delta} t_{\theta }$, for the airfoils stations $0.50 \leq x/c \leq 0.90$ on the suction side.

2.2. Zero-pressure-gradient turbulent boundary layer

Due to the large cost of the airfoil simulations, an additional set of simulations of a flat-plate boundary layer is performed to assess the influence of the spanwise domain width. The more canonical setting can then be discretised with more specialised methods. We use the pseudo-spectral solver Simson (Chevalier, Lundbladh & Henningson Reference Chevalier, Lundbladh and Henningson2007) to simulate a zero-pressure-gradient turbulent boundary layer (ZPG TBL), as in Schlatter et al. (Reference Schlatter, Li, Örlü, Hussain and Henningson2014). The horizontal and vertical lengths for all the three simulations are $L_x = 1000 \delta ^{*}$ and $L_y = 90 \delta ^{*}$, respectively, and simulations are carried out for three spanwise lengths, $L_z = 96 \delta ^{*}$, $L_z = 192 \delta ^{*}$ and $L_z = 384 \delta ^{*}$, where $\delta ^{*}$ is the displacement thickness of the laminar Blasius boundary layer at the inflow boundary, a location at which the displacement-thickness Reynolds number is set to $450$ for all cases. The inflow boundary is located at $x/\delta ^{*}=0$, as shown in figure 2(a), where the instantaneous field of the streamwise velocity component $u$ for the ZPG TBL simulation with $L_z = 192 \delta ^{*}$ is plotted. Note that $x/\delta ^{*}$ and $y/\delta ^{*}$ represent the computational coordinate, with $\delta ^{*}$ as a reference parameter.

Figure 2.

ZPG TBL simulation results: (a) instantaneous field of the streamwise velocity component in the computational coordinate, for the ZPG TBL simulation with $L_z = 192 \delta ^{*}$ (or $144 \theta _0$); (b) the friction coefficient as a function of $Re_{\theta }$ for the three ZPG TBL simulations compared with results from the equation deduced by Österlund et al. (Reference Österlund, Johansson, Nagib and Hites2000) (red dotted line).

The spatial resolution in the boundary layer in wall units is given by ${\rm \Delta} x^{+}=78$, ${\rm \Delta} y^{+}=0.2$–$9$ and ${\rm \Delta} z^{+}=40$, considering a reference position with $Re_{\delta ^{*}} = 450$. In addition, a fringe-region technique is used in order to employ a (periodic) Fourier discretisation in the streamwise direction: a volume force acts at the end of physical domain to reduce the boundary-layer thickness and damp out the outflowing boundary-layer disturbances. The fringe region is located at $x/\delta ^{*} \geq 900$, indicated by the right side of the red dashed line in figure 2(a).

Since we are dealing with wall-bounded turbulent flows, it is more convenient to use the momentum thickness to normalise the flow-field variables. Therefore, we choose as a reference parameter the momentum thickness where $Re_{\theta }=600$, which corresponds approximately to the same Reynolds number at the station $x/c = 0.60$ of the NACA 0012 airfoil, as shown in table 1. This reference momentum thickness will be represented here by $\theta _0$. At this boundary-layer position, the horizontal and vertical domain lengths become $L_x = 750 \theta _0$ and $L_y = 67 \theta _0$, respectively, and the three spanwise lengths become $L_z = 72 \theta _0$, $L_z = 144 \theta _0$ and $L_z = 288 \theta _0$. For these databases we extracted $1111$ snapshots from the simulations, with a non-dimensional separation of ${\rm \Delta} t_{\theta _0} = 6.7$, normalised by the free-stream velocity $U_\infty$ and $\theta _0$.

A LES with a finer mesh was also computed for the smallest spatial domain, $L_z = 72 \theta _0$, in order to show that the results throughout this paper regarding the ZPG TBL simulations do not depend on the resolution of the LES. Using the same Reynolds number as reference, $Re_{\delta ^{*}} = 450$, the spatial resolution in the boundary layer in wall units is given by ${\rm \Delta} x^{+}=26$, ${\rm \Delta} y^{+}=0.1$–$4.5$ and ${\rm \Delta} z^{+}=13$ for this simulation.

Figure 2(b) shows the friction coefficient as a function of the local momentum-thickness Reynolds number, $Re_{\theta }$, for the three ZPG TBL simulations, which is essentially identical in the three cases, as expected. These results are also in agreement with the equation deduced by Österlund et al. (Reference Österlund, Johansson, Nagib and Hites2000) (red dotted line), given by

(2.1)\begin{equation} c_f = 2 \left[ \frac{1}{\kappa} \ln(Re_{\theta}) + C \right]^{{-}2} \end{equation}

where the value of the von Kármán constant is $\kappa = 0.384$ and the additive constant $C = 4.08$. Note that the initial part of the domain, with $Re_\theta <550$, comprises not fully developed turbulence, and the sponge zone is located at $Re_{\theta } \geq 990$.

To allow comparisons between airfoils and the ZPG TBL simulations, table 2 shows the characteristics of all the numerical database analysed in this work. Note that the ZPG TBL simulation with the smallest $L_z$ has approximately the same span as the airfoils simulations.

Table 2.

Characteristics of the numerical simulations, using as reference the station $x/c = 0.60$ at the top side for the airfoil simulations and the station where $Re_{\theta }=600$ for the ZPG TBL simulations.

4.1. SPOD results

In this section we will analyse both airfoil databases using SPOD to highlight the dominant spanwise-coherent structures in the flow. Figures 3(a) and 3(b) show the square root of the first five SPOD eigenvalues as a function of the Strouhal number for the NACA 0012 and NACA 4412 airfoils, respectively. Note that SPOD modes have unit norm, and thus the square root of the eigenvalues corresponds to the amplitude that each mode has in the flow. We see from the results that the amplitudes are larger for the NACA 4412, as expected, since this is a cambered airfoil with an angle of attack of $5^{\circ }$, which induces the appearance of structures with larger amplitudes of velocity fluctuations in the suction side. We can also observe a peak in $St \approx 7$ for both airfoils, which is more pronounced for the first eigenvalue.

Figure 3.

Contribution of the leading five SPOD eigenvalues as a function of the Strouhal number for the airfoils: NACA 0012 (a,c,e) and NACA 4412 (b,d,f). The relative contributions are shown using: eigenvalues square root (a,b), ratios of the first SPOD eigenvalue and subsequent ones (c,d) and the contribution of each mode to the total TKE (e,f).

In order to evaluate the dominance of the first SPOD mode, we analysed the ratio of the first SPOD eigenvalue and subsequent ones, and the results are shown in figures 3(c) and 3(d), for the NACA 0012 and NACA 4412 airfoils, respectively. Note that the largest difference between first and second eigenvalues corresponds to the same Strouhal number of the peak, $St \approx 7$; such a difference indicates the dominance of the leading mode in the flow fluctuations, with an amplitude that is at least twice larger than the amplitudes of higher-order modes. We also evaluated the contribution of each mode to the total TKE, shown in figures 3(e) and 3(f), for both NACA 0012 and NACA 4412 airfoils, respectively. Here, we consider the TKE of spanwise-averaged fluctuations. Results show the dominance of the first SPOD mode at $St \approx 7$, which comprises approximately $60\,\%$ of the total TKE for both airfoils. The recent study by our group (Sano et al. Reference Sano, Abreu, Cavalieri and Wolf2019) using a compressible LES of a NACA 0012 airfoil showed a peak of trailing-edge noise at approximately $St \approx 7$, which is the same peak frequency found here for the velocity fields of both analysed airfoils. Therefore, we focus our analysis on $St = 7$ for both airfoils.

Figures 4 and 5 show the first SPOD mode of the velocity and pressure fields at $St = 7$ for the NACA 0012 and NACA 4412 airfoils, respectively. The real part of the modes is shown; the imaginary part forms a $90^{\circ }$ phase shifted pattern, indicating a downstream travelling wave. All the SPOD results here are multiplied by the square root of the eigenvalue, $\sqrt {\lambda _1}$, as in Sinha et al. (Reference Sinha, Rodríguez, Brès and Colonius2014), such that these correspond to fluctuation amplitudes in the flow. We can see that all the results shown in figures 4 and 5 exhibit a coherent wave through most of the airfoil, which begins just after the boundary-layer tripping, with increasing amplitude as it gets closer to the trailing edge, and peak amplitudes in the wake. Comparing the two airfoils, the leading SPOD mode of the 4412 profile has larger fluctuation amplitudes in the wake, but both airfoils have leading SPOD modes with nearly the same amplitudes inside the boundary layer considering the region $x/c \leq 0.85$, where the tangential velocity fluctuation amplitude reaches $|\widehat {u'}| \approx 0.0025$ for both airfoils. This amplitude scaled in wall units corresponds to $|\widehat {u'}^{+}| \approx 0.055$, normalised by $u_{\tau }$ and $\nu$, for the station $x/c=0.65$, where the momentum-thickness Reynolds number is approximately $690$ and $650$ for the NACA 0012 and NACA 4412 airfoils, respectively. Note that the amplitude of the SPOD modes in viscous scaling is very small with respect to dominant turbulent structures, with fluctuation amplitudes of order unity when scaled in inner units. Therefore, these leading SPOD modes show dominant structures of the TKE of spanwise-averaged fluctuations, as shown in figure 3, but these are only a small part of the bulk of turbulence. Such structures are nonetheless expected to play a significant role in trailing-edge noise, as discussed in the introduction.

Figure 4.

First SPOD mode at $St = 7$ for the NACA 0012 airfoil, representing approximately $56\,\%$ of the total kinetic energy of spanwise-averaged fluctuations; (a), (b) and (c) show fluctuations of $u$, $v$ and $p$, respectively: (a) $\sqrt {\lambda _1} \cdot \boldsymbol {\varPsi _{u1}}$; (b) $\sqrt {\lambda _1} \cdot \boldsymbol {\varPsi _{v1}}$; and (c) $\sqrt {\lambda _1} \cdot \boldsymbol {\varPsi _{p1}}$.

Figure 5.

First SPOD mode at $St = 7$ for the NACA 4412 airfoil, representing approximately $60\,\%$ of the total kinetic energy of spanwise-averaged fluctuations; (a), (b) and (c) show fluctuations of $u$, $v$ and $p$, respectively: (a) $\sqrt {\lambda _1} \cdot \boldsymbol {\varPsi _{u1}}$; (b) $\sqrt {\lambda _1} \cdot \boldsymbol {\varPsi _{v1}}$; and (c) $\sqrt {\lambda _1} \cdot \boldsymbol {\varPsi _{p1}}$.

In order to identify such coherent structures in the flow, we show in figures 6(a) and 6(b) sample instantaneous flow realisations of the spanwise-averaged fluctuations of the tangential velocity and pressure, respectively, for the NACA 0012 airfoil database. Results show that the coherent structures found in the first SPOD mode can be more clearly observed in the pressure field, as figure 6(b) displays instantaneous structures that resemble the pressure fluctuations of the leading SPOD mode shown in figure 4(c). For the velocity field, the coherent structures corresponding to the leading SPOD mode can be observed more clearly in positions further from the wall; for instance, the fluctuations in $y$ positions further from the wall have a more organised wave behaviour, which is close to the result for the leading SPOD mode of figure 4(a). The clearer observation of coherent waves further from the wall is reminiscent of wavepackets in turbulent jets (Tinney & Jordan Reference Tinney and Jordan2008; Cavalieri et al. Reference Cavalieri, Jordan and Lesshafft2019), which have a clearer structure in the near pressure field. The NACA 4412 airfoil has similar behaviour and results have not been shown here for brevity.

Figure 6.

Instantaneous flow realisations taken from the NACA 0012 airfoil database. (a) Tangential velocity fluctuations $u'$ and (b) pressure fluctuations $p'$.

Figures 7 and 8 show the second and third SPOD modes of the tangential velocity component at $St = 7$ for the NACA 0012 and NACA 4412 airfoils, respectively. Such SPOD modes also display coherent hydrodynamic waves. However, as we increase the order of the mode, the amplitudes of the overall wavy structures decrease, as expected, since their contribution to the total TKE decreases significantly. In particular, differently from the first mode, the amplitudes of velocity fluctuations for the second mode are higher in the turbulent boundary-layer region than in the wake. The third SPOD mode has higher amplitudes near the trailing-edge region, but the wavy structures are noisier due to a worse convergence of higher modes.

Figure 7.

Second and third SPOD modes of the tangential velocity component at $St = 7$ for the NACA 0012 airfoil; (a) $\sqrt {\lambda _2} \cdot \boldsymbol {\varPsi _{2}}$ and (b) $\sqrt {\lambda _3} \cdot \boldsymbol {\varPsi _{3}}$.

Figure 8.

Second and third SPOD modes of the tangential velocity component at $St = 7$ for the NACA 4412 airfoil; (a) $\sqrt {\lambda _2} \cdot \boldsymbol {\varPsi _{2}}$ and (b) $\sqrt {\lambda _3} \cdot \boldsymbol {\varPsi _{3}}$.

Such coherent hydrodynamic waves appear in the first SPOD mode of the airfoils even though the boundary layer is at an already fully developed turbulent state. The trailing-edge region features a larger adverse pressure gradient in the boundary layer, which is known to enhance boundary-layer instability for laminar flows (Schmid & Henningson Reference Schmid and Henningson2001), and the present results suggest that a similar mechanism may be at work in the turbulent boundary layers near the trailing edge, especially for the cambered NACA 4412 airfoil. Furthermore, the SPOD results for the tangential velocity component in all cases (see figures 4(b) and 5(b)) show near-wall fluctuations with phase opposition to disturbances towards the boundary-layer edge, which is similar to the behaviour of Tollmien–Schlichting waves, although this is a region of turbulent boundary layer. A similar behaviour was found by Sano et al. (Reference Sano, Abreu, Cavalieri and Wolf2019) for the NACA 0012 airfoil and by Kaplan et al. (Reference Kaplan, Jordan, Cavalieri and Brès2021) for a turbulent boundary layer inside a jet nozzle. The present results indicate that such structures are not just an artefact of a particular database.

In order to further highlight the structures in the boundary-layer region, the wake contribution is neglected in SPOD by imposing zero integral weights $\boldsymbol {W}$ in the wake region, where $x/c>1$. This amounts to considering a subdomain $\varOmega$ in (3.1) that excludes the wake. We refer to these as boundary-layer-weighted SPOD modes. The results in that case are shown in figures 9(a) and 9(b) for the NACA 0012 and in figures 10(a) and 10(b) for the NACA 4412 airfoil. The SPOD results now emphasise the hydrodynamic waves near the airfoils, inside the boundary layer. Since the wake region typically has a Kelvin–Helmholtz mode which is much more energetic than the ones from the turbulent boundary layer, which has just stable modes, boundary-layer-weighted SPOD modes highlight more clearly structures within the boundary layer; inspection of the figures shows that these are the same global SPOD modes, but convergence-related noise is lower when the wake is discarded. As the wake dynamics is more thoroughly studied in the literature (see for instance Araya, Colonius & Dabiri Reference Araya, Colonius and Dabiri2017), we will thus focus for the remainder of our analysis on the boundary-layer-weighted modes.

Figure 9.

First boundary-layer-weighted SPOD mode at $St = 7$ for the NACA 0012 airfoil, representing approximately $40\,\%$ of the kinetic energy of spanwise-averaged fluctuations; (a) $\sqrt {\lambda _1} \cdot \boldsymbol {\varPsi _{1}}$ tangential velocity component $\hat {u}'$, (b) $\sqrt {\lambda _1} \cdot \boldsymbol {\varPsi _{1}}$ wall-normal velocity component $\hat {v}'$.

Figure 10.

First boundary-layer-weighted SPOD mode at $St = 7$ for the NACA 4412 airfoil, representing approximately $30\,\%$ of the kinetic energy of spanwise-averaged fluctuations; (a) $\sqrt {\lambda _1} \cdot \boldsymbol {\varPsi _{1}}$ tangential velocity component $\hat {u}'$, (b) $\sqrt {\lambda _1} \cdot \boldsymbol {\varPsi _{1}}$ wall-normal velocity component $\hat {v}'$.

4.2. Comparison between SPOD and resolvent modes

In this section we will compare the SPOD results obtained in § 4.1 for both airfoils with the most amplified responses from resolvent analysis. To do so, we will make a locally parallel assumption for the resolvent analysis, and compare the first resolvent mode with the first SPOD mode for some specific stations $x/c$ of the airfoil by taking a ‘slice’ from the boundary-layer-weighted SPOD results obtained for the entire field in figures 9 and 10.

4.2.1. The base flow

The base flow $\boldsymbol {\bar {u}}(y_n)$ was extracted from the simulations considering a time and spanwise average. Here, the mean profiles are taken in the wall-normal direction, in the overall $x$–$y$ coordinate system. For the NACA 0012 airfoil we analysed three different stations: $x/c=0.6$, $0.7$ and $0.9$. For the NACA 4412 airfoil we also analysed three different stations: $x/c=0.6$ and $0.9$ for the suction side (a) and $0.5$ on the pressure side (b), since this airfoil has no symmetry with respect to $y=0$. These choices allow consideration of boundary-layer profiles with various degrees of adverse pressure gradient, which is stronger close to the trailing edge, particularly for the NACA 4412 profile, and close to zero for the pressure side of the same airfoil. For all stations considered here, in both airfoils, the boundary-layer profiles are turbulent. Figure 11 shows the mean streamwise velocity around the airfoils, with various $x/c$ stations marked by the dashed red lines. The mean velocity profiles at these $x/c$ positions for both airfoils are shown in figure 12, scaled in outer units (wing chord $c$ and free-stream velocity $U_{\infty }$). Figure 13 shows the mean velocity profiles for the same stations in inner units, denoted with a $+$ superscript; quantities are made non-dimensional using the friction velocity $u_\tau$ and the kinematic viscosity $\nu$. The mean velocity profiles in figure 13 show the expected shape of wall-bounded turbulent flows when plotted in inner scale as a function of wall distance $y^{+}$. Adverse pressure gradients lead to departures of the typical zero-pressure-gradient profiles, as evident for the NACA 4412 results (see also Vinuesa et al. Reference Vinuesa, Bobke, Örlü and Schlatter2016).

Figure 11.

Mean streamwise velocity around the airfoils, with $x/c=0.4$, $0.5$, $0.6$, $0.7$, $0.8$ and $0.9$ stations marked by the dashed red lines. (a) NACA 0012 and (b) NACA 4412.

Figure 12.

Mean velocity profiles at different positions $x/c$ for (a) NACA 0012 and NACA 4412 on (b) suction and (c) pressure sides, scaled in outer units.

Figure 13.

Mean velocity profiles at different positions $x/c$ for (a) NACA 0012 and NACA 4412 on (b) suction and (c) pressure sides, scaled in inner units. In all panels, the blue dashed line represents the linear law, $\bar {u}^{+} = y^{+}$, and the red dashed line is the log law, $\bar {u}^{+} = (1/\kappa ) \ln (y^{+}) + B$, where $\kappa = 0.41$ and $B=5.2$.

4.2.2. Resolvent parameters

Resolvent analysis was performed for $St=7$ for comparison with the SPOD results of § 4.1. The streamwise wavenumber $k_x$ is also a parameter of resolvent analysis, and is not known a priori, since different wavelengths may appear in the forcing. To study the effect of $k_x$ in the resolvent analysis, we first varied $k_x$ and studied the leading gains for the three analysed stations of each airfoil. These results are shown in figure 14 for the respective stations analysed in the NACA 0012 airfoil, and in figure 15 for the NACA 4412.

Figure 14.

First three resolvent gains as a function of the streamwise wavenumber at $St = 7$ for the NACA 0012 airfoil; (a) $x/c=0.6$, (b) $x/c=0.7$ and (c) $x/c=0.9$.

Figure 15.

First three resolvent gains as a function of the streamwise wavenumber at $St = 7$ for the NACA 4412 airfoil; (a) $x/c=0.5$, (b) $x/c=0.6$ and (c) $x/c=0.9$.

We notice a sharp peak in all results, which corresponds to disturbances outside the boundary layer, where the phase velocity is equal to the free-stream velocity. Analysis of forcing and response modes for these wavenumbers shows that they correspond to free-stream modes, similar to the observations of Dergham, Sipp & Robinet (Reference Dergham, Sipp and Robinet2013) for a backward-facing step. We have verified that forcing and response modes in these cases have spatial support mostly outside the boundary layer. Such free-stream modes were studied by Nogueira et al. (Reference Nogueira, Cavalieri, Hanifi and Henningson2020) and are not of immediate interest for the boundary-layer dynamics. Besides the sharp peak, we can also notice smooth peaks in all stations for the first gain $\sigma _{1}$. These peaks correspond to disturbances inside the boundary layer, where the phase velocity is lower than the free-stream velocity. Moreover, for the smooth peaks there is a clear dominance of the optimal forcing, which leads to a response with a gain that is at least one order of magnitude higher than the suboptimals. To perform the resolvent analysis we thus used the $k_x$ corresponding to the larger smooth peaks in each analysed section for both airfoils. For the NACA 0012 airfoil we notice that the value used in resolvent analysis, $k_x = 50$, corresponds to a wavelength $\lambda _x = 0.12$, which closely matches the wavelength of the first SPOD mode shown in figure 4.

4.2.3. Comparison between SPOD and resolvent modes

Results of the comparison between first boundary-layer-weighted SPOD mode and optimal response from resolvent analysis are shown in figures 16 and 17 for the NACA 0012 airfoil for horizontal and vertical velocity components, respectively; and in figures 18 and 19 for the NACA 4412 airfoil, also for horizontal and vertical velocity components, respectively. As commented before, the boundary-layer-weighted SPOD modes were taken to minimise statistical noise.

Figure 16.

First SPOD mode (red dashed line) compared with the optimal response from resolvent analysis (black solid line) of the streamwise velocity component $u$, for the NACA 0012 airfoil, at stations $x/c=0.6$, $0.7$ and $0.9$; (a) $\mu \approx 0.80$, (b) $\mu \approx 0.83$ and (c) $\mu \approx 0.85$.

Figure 17.

First SPOD mode (red dashed line) compared with the optimal response from resolvent analysis (black solid line) of the wall-normal velocity component $v$, for the NACA 0012 airfoil, at stations $x/c=0.6$, $0.7$ and $0.9$; (a) $\mu \approx 0.80$, (b) $\mu \approx 0.83$ and (c) $\mu \approx 0.85$.

Figure 18.

First SPOD mode (red dashed line) compared with the optimal response from resolvent analysis (black solid line) of the streamwise velocity component $u$, for the NACA 4412 airfoil, at stations $x/c=0.5$ in the pressure side, $x/c=0.6$ and $0.9$ in the suction side; (a) $\mu \approx 0.75$, (b) $\mu \approx 0.67$ and (c) $\mu \approx 0.70$.

Figure 19.

First SPOD mode (red dashed line) compared with the optimal response from resolvent analysis (black solid line) of the wall-normal velocity component $v$, for the NACA 4412 airfoil, at stations $x/c=0.5$ in the pressure side, $x/c=0.6$ and $0.9$ in the suction side; (a) $\mu \approx 0.75$, (b) $\mu \approx 0.67$ and (c) $\mu \approx 0.70$.

The comparisons show that the resolvent modes have overall features that correspond to the leading SPOD modes from the databases. In order to perform a quantitative comparison between the first SPOD mode and the optimal response from resolvent analysis, we quantified the agreement between the SPOD eigenfunction $\varPsi _{1}$ and optimal response from the resolvent analysis $U_1$ using the normalised projection given by (3.5), as in Abreu et al. (Reference Abreu, Cavalieri, Schlatter, Vinuesa and Henningson2020a,Reference Abreu, Cavalieri, Schlatter, Vinuesa and Henningsonb). Values of $\mu$ close to 1 are indicative of agreement of both amplitudes and phases for both velocity components. The agreement between resolvent and SPOD is closer for the 0012 airfoil, with normalised projections higher than $0.8$ for all stations. The comparison for the 4412 airfoil has the best agreement at the pressure side, with a normalised projection of $\mu =0.75$. The suction side has agreement metrics close to $0.7$. These results can be explained by the fact that the locally parallel resolvent analysis is based on properties in the near wall, which are simpler in the lower adverse pressure gradient from the NACA 0012 than in the complex adverse pressure gradient of the NACA 4412 (interaction of a wider range of scales with stronger adverse pressure gradients), which leads to a stronger inhomogeneity in $x$, such that the accuracy of the locally parallel assumption becomes questionable. However, the general features of the SPOD mode are captured by the optimal response for all cases presented in figures 16–19. The double peak structure of the streamwise velocity is related to a phase jump in $y$, which may be observed in figure 4(a). Such phase jump is also obtained in the leading resolvent mode, as exemplified in figure 20, which shows phase comparisons for a station of the NACA 0012 airfoil. Other stations show similar results and are not shown here for brevity.

Figure 20.

Phase of the first SPOD mode (red dashed line) compared with the phase of the optimal response from resolvent analysis (black solid line) for the streamwise (a) and wall-normal velocity components (b), at the station $x/c=0.7$ of the NACA 0012 airfoil. (a) Phase of $u$ (rad) and (b) Phase of $v$ (rad).

Since we deal with a turbulent boundary layer, an exact match is not expected, but the comparisons presented here show that resolvent analysis provides a good reduced-order model for the spanwise-coherent disturbances considered here, which play an important role in the radiation of trailing-edge noise (Sano et al. Reference Sano, Abreu, Cavalieri and Wolf2019). There is close agreement with streamwise velocity fluctuations of the leading SPOD mode, which has an amplitude distribution similar to that of Tollmien–Schlichting waves, with a phase jump at a position with vanishing amplitude, also visible in the SPOD fields in figures 4–10. The relatively low amplitude of the two-dimensional disturbances within the boundary layer tends to ease the application of linearised models.

Resolvent analysis also shows consistent trends when the two airfoils are compared. The comparison of figures 14 and 15 shows that resolvent gains are larger for boundary-layer stations with adverse pressure gradient. This is observed in the SPOD modes of figures 4–10, which display larger amplitudes at the suction side of the NACA 4412 airfoil, particularly close to the trailing edge. The NACA 0012 airfoil fluctuation amplitudes are less significant, and the lowest levels are seen at the pressure side of the NACA 4412 profile. These trends are also observed in the resolvent gains.

5.1. Turbulence statistics and spectral analysis

Using standard Reynolds decomposition $u=\bar {u} + u'$, the mean horizontal velocity profiles $\bar {u}^{+}$ and the variance profiles of the horizontal velocity fluctuations $(\overline { u'}^{+})^{2}$ in inner scaling are shown in figures 21(a) and 21(b), respectively, for the position with $Re_{\theta }=990$, which corresponds to approximately the station $x/c \approx 0.80$ of the NACA 0012 airfoil and $x/c \approx 0.75$ of the NACA 4412 profile. The mean velocity profiles for all the three simulations collapse as expected, and show the typical shape of wall-bounded turbulent flows when plotted as a function of wall distance $y^{+}$, following linear and log laws. Variance profiles also collapse for the three simulations, and have the expected pattern characteristic of wall-bounded turbulent flows, with a near-wall peak in the buffer layer at $y^{+} \approx 15$ (Eitel-Amor et al. Reference Eitel-Amor, Örlü and Schlatter2014). Note that, despite the moderate resolution of the ZPG TBL simulation, the results presented below do not depend on this choice of lower resolution, which was due to the computational efficiency. In order to verify this statement, results for a finer LES for the smallest spatial domain are shown in Appendix.

Figure 21.

(a) The mean flow and (b) horizontal velocity fluctuations at $Re_{\theta }=990$, scaled in inner units. In both panels, the blue dashed line represents the linear law, $\bar {u}^{+} = y^{+}$; and the red dashed line is the log law, $\bar {u}^{+} = (1/\kappa ) \ln (y^{+}) + B$, where $\kappa = 0.41$ and $B=5.2$.

The inner-scaled premultiplied power-spectral density of horizontal velocity fluctuations $k_{z}^{+} E_{uu}^{+}$ as a function of spanwise wavelength $\lambda _z^{+}=2{\rm \pi} / k_z^{+}$ is shown in figure 22(a) at position $Re_{\theta }=990$, for all the three simulations. Similarly to the mean profile and the variance, the premultiplied spectra for the three simulations collapse, showing further indication of independence of domain size. Results in figure 22(a) show a peak in the spanwise wavelength corresponding to $\lambda _{z}^{+} \approx 100$, which is representative of the signature of the near-wall cycle of streaks and quasi-streamwise vortices. This has been documented in many studies across a range of Reynolds numbers and flow types (Jiménez Reference Jiménez2013). The same inner-scaled power-spectral density of horizontal velocity fluctuations, without premultiplication $E_{uu}^{+}$, is shown in figure 22(b) at the same position, $Re_{\theta }=990$, for all the three simulations. Figure 22(c) shows the power-spectral density of figure 22(b) for selected wall-normal positions, highlighting the low-wavenumber behaviour. Results show that the power-spectral density contours do not close when premultiplication is not considered, even with the increasing spanwise length $L_z$, indicating that there is non-zero power-spectral density for very low wavenumbers $k_z^{+} \to 0$, or, equivalently, very large values of $\lambda _z^{+}$. We can see that there are structures with $k^{+}_{z}=0$, which correspond to two-dimensional perturbations in the computational domain, corresponding to the structures studied for the airfoils for $k^{+}_{z}=0$.

Figure 22.

Spanwise power-spectral density of the streamwise velocity fluctuation $u'$ at location $Re_{\theta }=990$ for the three ZPG TBL simulations: (a) with premultiplication $k_{z}^{+} E_{uu}^{+}$ and $x$-axis as $\lambda _{z}^{+}$ shown in log scale; and (b) without premultiplication $E_{uu}^{+}$ and $k_{z}^{+}$ shown in linear scale; and (c) without premultiplication $E_{uu}^{+}$ and $k_{z}^{+}$ at specific $y^{+}$ positions.

The spectrum converges to a finite value at $k^{+}_{z}=0$ and there is no evidence in the current results, which already include much wider domains than usual, that the power of $k^{+}_{z}=0$ will decay to zero. The wavenumber $k^{+}_{z} \to 0$ corresponds to spanwise-coherent velocity fluctuations, which have not been studied in detail in the literature and are the main focus of this study. This limit is thus independent of domain size for the computations shown here. As discussed by Lozano-Durán & Jiménez (Reference Lozano-Durán and Jiménez2014), two-dimensional fluctuations in the simulation, for $k_z^{+}=0$, represent disturbances whose spanwise length scale is larger than the size of the computational box; they should be thus understood as fluctuations with $\lambda _z \gg \delta$, which are responsible for the low-wavenumber behaviour of power-spectral density shown in figure 22(b).

5.2. SPOD and resolvent results

As discussed previously, in order to compare the results obtained with the ZPG TBL simulations with those from the airfoils, we normalised the parameters using the momentum thickness of reference, $\theta _0$, corresponding to the one taken at $60\,\%$ of the NACA 0012 airfoil chord, where the momentum-thickness Reynolds number is approximately $600$. The region of interest in the ZPG TBL was considered to be the one corresponding to the trailing-edge region, i.e. between $60\,\%$ and $80\,\%$ of the airfoil chord length, where the hydrodynamic waves have larger amplitudes. Moreover, the frequency used to perform the SPOD in the airfoils, $St \approx 7$, corresponds to $St_{\theta } \approx 0.010$ at $60\,\%$ of the NACA 0012 airfoil chord. We will then perform the SPOD in the region of interest of the ZPG TBL databases, by proper selection of a domain of interest $\varOmega$, considering $St_{\theta _0} = 0.010$.

Similarly to the analysis of the airfoil databases, we have taken a spanwise average of flow fluctuations before performing the SPOD, to focus our analysis on spanwise-coherent structures. SPOD for the three ZPG TBL databases was performed in a fully turbulent region delimited by $600 \leq Re_{\theta } \leq 990$, which corresponds approximately to the NACA 0012 airfoil stations $0.60 \leq x/c \leq 0.80$ (see table 1 for the correspondent NACA 4412 airfoil stations). The focus on the region of interest ensures that only fluctuations in the developed turbulent region are considered, with minimal influence from the upstream tripping. Furthermore, the size of the region ensures that disturbances are able to travel across the entire region within the duration of each time segment used in Welch's method.

Figure 23 shows the SPOD eigenvalues for the three ZPG TBL simulations for $k_z = 0$, which show nearly no change among the three domain sizes. We thus conclude that the coherent structures in all cases have nearly the same amplitude. Figures 24(a–c) and 25(a–c) show the first SPOD mode for horizontal and vertical velocity components, respectively, for all the three ZPG TBL simulations. We can observe the presence of coherent waves for all the cases, with minimal changes as the spanwise length of the domain is varied. The SPOD results for the ZPG TBL databases in these figures display nearly the same mode shapes and amplitudes for the three considered domain sizes, and the structures are very similar to those around the airfoils, shown in figures 9 and 10. For the ZPG TBL modes, the amplitude of the horizontal velocity fluctuation reaches $|\widehat {u'}| \approx 0.0028$, or $|\widehat {u'}^{+}| \approx 0.061$ in wall units, normalised by $u_{\tau }$ and $\nu$ from the beginning of the domain, $x/\theta _0=0$, where $Re_{\theta }=600$. Note that the amplitudes of the SPOD mode of the ZPG TBL are very close to those from the airfoils, and they also make a small contribution to the overall turbulence, but have a potential effect in sound radiation. The present SPOD results indicate that the coherent hydrodynamic waves found in the present study, corresponding to spanwise-coherent disturbances in turbulent boundary layers, have a physical meaning and are not an artefact from the simulations.

Figure 23.

SPOD eigenvalues for the three ZPG TBL simulations.

Figure 24.

First SPOD mode of the horizontal velocity fluctuations (a,b) compared with the optimal response from resolvent analysis (c,d) for the three analysed ZPG TBL simulations at $k_z = 0$; (a) $\varPsi _{u1}$ $L_z = 72 \theta _0$, (b) $\varPsi _{u1}$ $L_z = 144 \theta _0$, (c) $\varPsi _{u1}$ $L_z = 288 \theta _0$ and (d) optimal response.

Figure 25.

First SPOD mode of the vertical velocity fluctuations (a,b) compared with the optimal response from resolvent analysis (c,d) for the three analysed ZPG TBL simulations at $k_z = 0$; (a) $\varPsi _{v1}$ $L_z = 72 \theta _0$, (b) $\varPsi _{v1}$ $L_z = 144 \theta _0$, (c) $\varPsi _{v1}$ $L_z = 288 \theta _0$ and (d) optimal response.

Figures 24(d) and 25(d) show the optimal response from resolvent analysis for horizontal and vertical velocity components, respectively, for $k_z = 0$. In this case, a global resolvent analysis was carried out, allowing a better quantitative comparison with the SPOD modes, since the streamwise wavenumber $k_x$ is no longer specified, and appears naturally in the structure of the optimal flow response. As the mean flows of the three simulations are nearly identical, resolvent analysis was computed using the mean velocity for the smallest computational domain. The frequency used to compute the resolvent analysis here was $St_{\theta _0}=0.010$. We observe a very close agreement between the SPOD modes and the optimal response from resolvent analysis, for both velocity components. The observed agreement further confirms that spanwise-coherent structures in turbulent boundary layers can be modelled as the most amplified flow response, obtained using the linearised N–S operator.

We now examine non-zero spanwise wavenumbers in order to evaluate if the spanwise-averaged fluctuations of $k_z=0$ are representative of disturbances for other low spanwise wavenumbers. Figures 26, 27 and 28 show the first SPOD mode and the optimal response from resolvent analysis of the ZPG TBL simulations for the horizontal, vertical and spanwise velocity components, respectively, for $k_z = 2{\rm \pi} /L_z$, with $L_z$ taken for the smallest domain size, such that this non-zero wavenumber is obtained for the three simulations. We can observe that these results are very similar to the results for $k_z = 0$, also showing a coherent wave with similar features. As these are no longer two-dimensional modes, they display a spanwise velocity component, which also has a wave structure. All three velocity components are similarly matched by the leading resolvent mode, which confirms the pertinence of resolvent analysis to model such structures. The observed similarity of $u$ and $v$ fluctuations of non-zero $k_z$ to the $k_z=0$ counterparts strongly suggests that two-dimensional fluctuations in the simulations are not spurious modes, and are actually representative of a range of low wavenumbers, or, alternatively, of disturbances with large spanwise coherence.

Figure 26.

First SPOD mode of the horizontal velocity fluctuations (a,b) compared with the optimal response from resolvent analysis (c,d) for the three analysed ZPG TBL simulations at $k_z = 2{\rm \pi} /L_z$; (a) $\varPsi _{u1}$ $L_z = 72 \theta _0$, (b) $\varPsi _{u1}$ $L_z = 144 \theta _0$, (c) $\varPsi _{u1}$ $L_z = 288 \theta _0$ and (d) optimal response.

Figure 27.

First SPOD mode of the vertical velocity fluctuations (a,b) compared with the optimal response from resolvent analysis (c,d) for the three analysed ZPG TBL simulations at $k_z = 2{\rm \pi} /L_z$; (a) $\varPsi _{v1}$ $L_z = 72 \theta _0$, (b) $\varPsi _{v1}$ $L_z = 144 \theta _0$, (c) $\varPsi _{v1}$ $L_z = 288 \theta _0$ and (d) optimal response.

Figure 28.

First SPOD mode of the spanwise velocity fluctuations (a,b) compared with the optimal response from resolvent analysis (c,d) for the three analysed ZPG TBL simulations at $k_z = 2{\rm \pi} /L_z$; (a) $\varPsi _{w1}$ $L_z = 72 \theta _0$, (b) $\varPsi _{w1}$ $L_z = 144 \theta _0$, (c) $\varPsi _{w1}$ $L_z = 288 \theta _0$ and (d) optimal response.

Figures 29(a) and 29(c) show the corresponding forcing from the optimal response presented in figures 24(d) and 25(d), considering $k_z = 0$, for the horizontal and vertical velocity components, respectively. Figures 29(b), 29(d) and 29(e) show the corresponding forcing from the optimal response presented in figures 26(d) and 27(d), considering $k_z = 2{\rm \pi} /L_z$, for the horizontal, vertical and spanwise velocity components, respectively. Notice that, whereas the domain of interest for the flow response was taken in the downstream region of the computational domain, corresponding to the region analysed using SPOD, the optimal forcing was determined throughout the physical domain, as it is known that the optimal forcing typically appears upstream of the flow response. The forcing modes also display a wave structure, with fluctuations inclined against the shear, which suggests that the forcing exploits the Orr mechanism, similar to the optimal excitation of Tollmien–Schlichting waves (Åkervik et al. Reference Åkervik, Ehrenstein, Gallaire and Henningson2008), and to coherent boundary-layer disturbances inside a nozzle (Garnaud et al. Reference Garnaud, Lesshafft, Schmid and Huerre2013; Kaplan et al. Reference Kaplan, Jordan, Cavalieri and Brès2021).

Figure 29.

Optimal forcing from the resolvent analysis for the ZPG TBL simulations; (a) $f_x$ for $k_z = 0$, (b) $f_x$ for $k_z = 2{\rm \pi} /L_z$, (c) $f_y$ for $k_z = 0$, (d) $f_y$ for $k_z = 2{\rm \pi} /L_z$ and (e) $f_z$ for $k_z = 2{\rm \pi} /L_z$.