2.1. Flow facility and rough wall

The Virginia Tech Stability Wind Tunnel (VTSWT) was used as the flow facility for all measurements conducted in this study. The VTSWT is a closed-loop, low-speed wind tunnel with a square test section of dimensions 1.85 $\times$ 1.85 $\times$ 7.32 m. A diagram of the VTSWT test section for this particular study is given in figure 2. A trip strip was placed upstream of the test section within the contraction to induce the formation of a research-quality turbulent boundary layer along all walls of the test section. A 0.914 m chord length NACA 0012 aerofoil is positioned in the centre of the test section and extends from the floor to the ceiling to generate different pressure gradient conditions along the port-side wall. The aerofoil quarter chord was located 3.45 m from the origin at the beginning of the test section, such that the measurement area has a near-constant pressure gradient. The most upstream sensor in the sub-convective pressure sensing array is positioned 2.46 m from the origin and extends for 443 mm in the streamwise direction. Measurements of the boundary layer parameters were made with a Pitot-static boundary layer rake positioned 2.64 m from the origin, wherein the base of the boundary layer rake was flush with the base of the roughness.

Figure 2.

Schematic of the experimental set-up within the test section of the Virginia Tech Stability Wind Tunnel. The port wall is fully configured as a rough wall and the starboard wall is configured as a smooth wall.

The port-side wall of the VTSWT test section was fully outfitted with 0.61 $\times$ 0.61 m panels covered in staggered cylindrical roughness elements beginning at the origin marked by the $x_{1}$ and $x_{2}$ axes downstream of the trip strip. The roughness was homogeneous and $k$ -type, and for the baseline flow case, this corresponded to a roughness Reynolds number of $k_{s}^{+}\approx 352$ , indicating a fully rough flow regime. Each roughness element had a diameter of 3.14 mm and was 2 mm in height (figure 3). Using a least squares regression of mean velocity profiles, the effective sand grain roughness was found to be $k_{s} \approx 1.6k_{g} = 3.2$ mm. This value for the sand grain roughness is in agreement with independent measurements made by Vishwanathan et al. (Reference Vishwanathan, Fritsch, Lowe and Devenport2023). The roughness elements were staggered, the pitch between elements was approximately 6.93 mm and was chosen such that the roughness panels tile exactly across the entire port wall of the test section, forming a fetch 1.85 m wide and 5.5 m in length. Figure 3 displays a schematic of the specific dimensions of this rough surface as well as a view of the rough wall as it was installed in the VTSWT. This rough-wall geometry was also selected because it has been studied extensively in the past, both experimentally and computationally (Mulchandani et al. Reference Mulchandani, Adams, Garcia-Mayoral, Vishwanathan, Fritsch, Lowe and Devenport2021; Fritsch et al. Reference Fritsch, Vishwanathan, Todd Lowe and Devenport2022; Fritsch et al. Reference Fritsch, Vishwanathan, Roy, Todd Lowe and Devenport2023; Vishwanathan et al. Reference Vishwanathan, Fritsch, Lowe and Devenport2023). This allows for the direct comparison between the data measured in this study and data measured previously. Outside of the array itself, the roughness along the port wall was moulded to aluminium panels using epoxy resin and CNC-machined moulds.

Figure 3.

(a) Zoom-in view of the rough wall showing arrangement of the staggered cylinder rough wall ( $s=6.93$ mm, $k_{g}=2$ mm, $d=3.14$ mm) and (b) view of the rough wall installed in the Virginia Tech Stability Wind Tunnel looking upstream.

The primary test condition for this study is a near-zero pressure gradient flow of $ \textit{Re}_{\tau }=10\,247$ and is nominally similar to that boundary layer studied by Fritsch et al. (Reference Fritsch, Vishwanathan, Todd Lowe and Devenport2022), Fritsch et al. (Reference Fritsch, Vishwanathan, Roy, Todd Lowe and Devenport2023) and Vishwanathan et al. (Reference Vishwanathan, Fritsch, Lowe and Devenport2023). Both these previous studies and this current study investigated a Reynolds number based on an aerofoil chord of $ \textit{Re}_{c} = 2.0\times 10^{6}$ . These previous studies used an identical set-up in the VTSWT test section, but made velocity, turbulence and point-wise pressure fluctuation measurements aimed at assessing boundary layer similarity for smooth- and rough-wall flows. They found that indirect methods for determining skin friction, as is done in this study, are highly sensitive to the selection of the logarithmic region. This has implications for comparisons based on wall-similarity, and for this study, is relevant for comparisons across Reynolds numbers and pressure gradients due to normalisation of the spectrum based on wall shear stress. They have also shown that when compared with similar smooth-wall flows, rough-wall flows under a pressure gradient show a better degree of wall similarity when using Reynolds number scaling, which bodes well for unifying wavenumber–frequency spectrum results across pressure gradients at constant Reynolds number. Additionally, they have shown that for single-point pressure spectra, scaling on outer variables is successful in unifying results across conditions. This approach is done throughout this study for both the single-point and wavenumber–frequency spectra with similar success, as is further discussed in §§ 4.2 and 4.5. Direct comparisons between the present study and the data collected by Fritsch et al. (Reference Fritsch, Vishwanathan, Roy, Todd Lowe and Devenport2023) are made to demonstrate the similarity between these boundary layer flows, with § 4.1 discussing the mean boundary layer properties and § 4.2 discussing the single-point wall-pressure spectra. Comparisons across studies are limited to boundary layer profiles and single-point pressure spectra, because the sub-convective pressure sensing array described in § 2.2 had not yet been developed.

2.2. Sub-convective pressure sensing array

Measurements of the fluctuating wall pressure were made using a linear array of sub-resonant sensors. This array of sensors was designed, manufactured and implemented with the aim of eliminating spatial aliasing in the sub-convective regime of the streamwise wavenumber–frequency spectrum. Each sub-resonant sensor was designed to communicate the pressures at the wall through a cavity to a microphone in the base of the sensor through the use of many small-diameter pores, thus allowing for the system to be modelled as a Helmholtz resonator with multiple necks. The sensor is operated at frequencies below that at which sound waves can form inside the cavity, since then the pressure at the base of the cavity is a coherent sum of the pressure at each of the pinholes. Figures 4(a) and 4(c) display schematics of the individual sensor profile, which features a tapered internal cavity to space sensors more tightly than the diameter of the microphones located within their base. Each sensor used a rectangular cross-section at the flow interface with dimensions of 50 mm in the spanwise direction and 6.43 mm in the streamwise direction. Each sensor was outfitted with 116 0.4 mm diameter pores evenly distributed in a $29\times 4$ grid over the rectangular cross-section through which the wall pressures were communicated to the sub-resonant sensor below. Sensor response is sensitive to three primary variables: the number of pores, the depth of each pore and the volume of the cavity. In the construction of these sensors, the goal was to maximise the resonant frequency of each sensor, which is done by increasing the number of pores, decreasing the depth of each pore and decreasing the volume of the cavity, resulting in a resonant frequency of approximately 2.8 kHz. Sensing areas were spaced every 6.93 mm in the streamwise direction (allowing for a 0.5 mm gap between adjacent areas), which equated to a maximum separation distance of 436.4 mm, or approximately $4.6\delta$ for the baseline flow condition. The roughness elements along the top of each sensor were hollowed such that the depth of each pore was the same for pores appearing in the roughness elements and on the substrate. This ensured the local geometry of, and around, each of the pores was identical. Each sensor cross-section was tapered such that a G.R.A.S Type 40PH-S5 1/4ʹʹ microphone could be inserted into the bottom of each cavity despite the separation between cavities being smaller than the microphone diameter. This had the added benefit of reducing the volume of each sensor; therefore, increasing the resonant frequency and expanding the frequency range of this sensing system. Adjacent sensors have mirrored sensor contours to allow for the system to be as compact as possible; therefore, reducing the separation distance between measurement locations. Due to the small size of each of the pores and an open area ratio of 4.4 %, this sensing strategy is effectively non-intrusive, as supported by previous analysis done by Damani et al. (Reference Damani, Devenport, Alexander, Szőke, Balantrapu and Starkey2025b ) and Szoke, Glegg & Devenport (Reference Szoke, Glegg and Devenport2021), and as displayed in the unperturbed velocity contours over the top of the sensing array displayed by Butt et al. (Reference Butt, Sharma, Damani, Totten, Devenport and Lowe2024).

Figure 4.

(a) Design and dimensions of a single-sensor profile. (b) View of manufactured array from the flow side highlighting the presence of the pores both through the roughness elements and the surrounding surface. (c) Side view of a single sensor highlighting hollow roughness elements. (d) View of the manufactured array from behind the port wall highlighting the alternating orientation of the sensor profile.

The advantage of this design is that the large sampling area above each sub-resonant cavity effectively area-averages the pressure. The effect of this spatial averaging is to filter much of the convective pressure fluctuations, which allows for an improved signal-to-noise ratio in the measurement of the sub-convective pressure fluctuations. Damani et al. (Reference Damani, Butt, Totten, Devenport and Lowe2025a ) used a similarly designed sensor array to measure the sub-convective pressure fluctuations beneath smooth-wall boundary layers. Sections 3.1 and 3.2 of Damani et al. (Reference Damani, Butt, Totten, Devenport and Lowe2025a ) derive the approach for modelling the wavenumber–frequency spectrum measured by an array of sensors of this type. They begin by modelling a single sensor, for which the spatial sensitivity can be determined according to (2.1). Here, $\boldsymbol{x}$ refers to the sensor coordinate space, $\boldsymbol{y}$ is the location of the centre of the sensor, $\boldsymbol{z}_{j}$ is the location of each pore and $s_{\kern-1pt j}$ is the relative sensitivity of the sensor to that pore. It should be noted that $s_{\kern-1pt j}$ is a complex value, but is subject to the restriction given in (2.2).

(2.1) \begin{align} \mathbb{A}\left (\boldsymbol{x}-\boldsymbol{y}\right )=\sum _{j} s_j\delta (\boldsymbol{x}-\boldsymbol{y}-\boldsymbol{z}_j), \\[-28pt] \nonumber \end{align}

(2.2) \begin{align} \sum _{j} s_j=1. \\[0pt] \nonumber \end{align}

Damani et al. (Reference Damani, Butt, Totten, Devenport and Lowe2025a ) use this relationship to determine what an array of sensors of this nature will measure. By first finding the wavenumber transform of the sensitivity function according to (2.3), one can derive an estimate of the cross-spectral density function that will be measured by the array according to (2.4), which can be evaluated directly to get the relation given in (2.5). This allows for the estimation of what an array of sensors of this type will measure the cross-spectral density function to be, $S_{qq}$ , given some model of the cross-spectral density function, $S_{pp}$ . Here, $n$ and $m$ refer to different sensors, and $N$ and $M$ are the number of pores on each of those sensors, respectively. For the array in this study, the number of pores per sensor is 116.

(2.3) \begin{align} \mathbb{A}\left (\boldsymbol{k}\right )=\frac {1}{\left (2\pi \right )^2}\ \sum _{j} s_je^{-i\boldsymbol{k}.\boldsymbol{z}_j}, \\[-28pt] \nonumber \end{align}

(2.4) \begin{align} S_{qq}(\Delta \boldsymbol{y}, \omega ) = 16\pi ^4 \iint _{-\infty }^{\infty } \varPhi _{pp}(\boldsymbol{k},\omega )|\mathbb{A}(\boldsymbol{k})|^2 e^{i \boldsymbol{k}.\Delta \boldsymbol{y}}\, \text{d}k_1 \,\text{d}k_3, \\[-28pt] \nonumber \end{align}

(2.5) \begin{align} S_{qq}\big (\Delta \boldsymbol{y}^{(m,n)},\omega \big )=\sum _{N}\sum _{M}{{s_N^{(n)}}^\ast s_M^{(m)}}S_{pp}\big (\boldsymbol{y}^{(n)}+\boldsymbol{z}_N^{(n)}-\boldsymbol{y}^{(m)}-\textbf {z}_M^{(m)},\omega \big ). \\[-12pt] \nonumber \end{align}

An estimate of the measured wavenumber–frequency spectrum can then be determined according to (2.6) by taking the Fourier transform of the cross-spectral density function in space, where $\overline {w^{2}}$ is the windowing energy,

(2.6) \begin{align} \varPhi _{qq}(k_1,\omega ) = \frac {\Delta y_1}{2\pi \sqrt {\overline {w^2}}} \textit{DFT}\left (S_{qq}^{(w)}(\Delta y_1^{(m,n)},\omega )\right ). \end{align}

To assess the error associated with a measurement of such an array, an estimate of its measurement of the wavenumber–frequency spectrum, $\varPhi _{qq}$ , must be compared with the true spanwise averaged wavenumber–frequency spectrum, $\varPhi _{AA}$ . Damani et al. (Reference Damani, Butt, Totten, Devenport and Lowe2025a ) shows that $\varPhi _{AA}$ can be found through the numerical integration of (2.7), where $N$ is a discrete number of subdivisions of the spanwise length of each sensor, $L_{3}$ , and these are related through $N \Delta x_{3} = L_{3}$ ,

(2.7) \begin{align} \varPhi _{AA} (k_1,\omega )=\frac {1}{N^2} \sum _{q=1-N}^{N-1}(N-|q|)\varPhi _{pp} (k_1,q\Delta x_3,\omega ) .\end{align}

To find both $\varPhi _{qq}$ and $\varPhi _{AA}$ , an input spectrum must be selected. For this analysis, the Hwang, Bonness & Hambric (Reference Hwang, Bonness and Hambric2003) modified Corcos model, as given in (2.8), is a good selection because its constants may be tailored to more closely approximate the true pressure spectrum. Note that there must be a selection for an input pressure spectrum $S_{pp}(\omega )$ , as well as the constants $\alpha _{1}$ and $\alpha _{3}$ , which control the decay rate of the spectrum in the streamwise and spanwise directions, respectively. A modified version of the Goody model was employed, and the two constants were selected as $\alpha _{1}=0.3275$ and $\alpha _{3}=0.77$ . This model construction was identical to that of Damani et al. (Reference Damani, Butt, Totten, Devenport and Lowe2025a ). As will be further discussed in §§ 3.2 and 4.5.1, there currently do not exist models that can approximate the wavenumber–frequency spectrum beneath a rough wall boundary layer flow accurately, so using the modified Corcos model was simply a best approximation.

(2.8) \begin{align} \varPhi _{pp} (k_1,\Delta x_3,\omega ) = \frac {S_{pp} (\omega ) U_c}{\pi \omega } \frac {2\alpha _1^3 e^{-|{\Delta x_3 \omega \alpha _3}/{U_c}|}}{ \big[\alpha _1^2+(U_c k_1/\omega -1)^2 \big]^2}. \end{align}

The final, and perhaps most important, input for modelling the error are the sensitivities associated with each pore $s_{\kern-1pt j}$ . A simple approach is to assume the sensitivity is uniform across the measurement area of the sensor, but both experiment and simulation have confirmed this to not be the case. By treating the sensor as an invertible system, the sensitivity at each pore can be estimated by exciting the cavity at its base and measuring the response at each of the pore locations. Figure 5 displays the relative sound pressure level measured 1 mm above the surface of a single sensor at a frequency of 1000 Hz. The locations of the roughness elements (large diameter circles) and each of the pores (small circles) have been indicated above the contour. It is clear that assuming a uniform spatial sensitivity is not sufficient, so further modelling was performed using COMSOL to determine the appropriate sensitivity for each pore. Figure 6 displays the corresponding contour of the relative sound pressure level throughout the domain as modelled using COMSOL. We observe decent agreement between the measured sensitivity profile in figure 5 and the computed sensitivity profile in figure 6, with the measured results displaying a slightly more uniform sensitivity distribution. This may be due to contamination from background noise, which the model eliminates completely.

Figure 5.

Experimental measurements of relative sound pressure level at 1000 Hz and 1 mm above a single sensor for an excitation signal provided through its base.

Figure 6.

COMSOL model of the relative sound pressure level at 1000 Hz and 1 mm above a single sensor for an excitation signal provided through its base.

Using the COMSOL simulation, a spatial sensitivity distribution was generated for every 50 Hz from 100 Hz to 2000 Hz. Figure 7 displays the estimated error between the measured wavenumber–frequency spectrum of the sensor array and the true spanwise-averaged wavenumber–frequency spectrum for both uniform sensitivities and the sensitivities generated using COMSOL. Assuming a uniform sensitivity distribution yields precisely what the array was designed to accomplish, which is close to zero error ( $\lt$ 1 dB) in the sub-convective regime, with larger errors ( ${\sim}4$ dB) at higher wavenumbers. When modelling the error using a non-uniform sensitivity distribution, the estimated error map remains relatively unchanged, except for the presence of a region of increased error ( $\sim$ 5–7 dB) in the acoustic cone at high frequencies. It should be noted that investigations into the spatial sensitivity were not performed until after data collection in the VTSWT had occurred. Future measurements using this array concept should improve upon this sensitivity distribution to ensure the minimum amount of error is present in the measurement. As previously mentioned, the sensors in this study aimed to minimise volume, which was done through a tapering of the internal geometry. Tapering can cause non-uniform sensitivity at frequencies below the cross-mode of the sensors, so future sensors should be designed with minimal tapering to ensure more uniform sensitivity and subsequently less error in the wavenumber–frequency spectrum.

Figure 7.

Estimated error between model estimate of the wavenumber–frequency spectrum measured by the sub-convective pressure sensing array and the modified Corcos model for (a) uniform sensitivities and (b) non-uniform sensitivities. The dashed lines represent the mean convective line, $k_{c}=2\pi f/U_{c}$ .

The final, manufactured array was built in two parts using stereolithography, which allowed for proper inspection and cleaning of the interior of each cavity. The two sections of the array were adhered together and sealed using grease. Each sub-resonant cavity was checked individually to ensure it was sealed from each other cavity. With this array design, the wavenumber resolution was 7.09 rad m⁻¹, and the Nyquist wavenumber limit was 453.5 rad m⁻¹. This sensor array was calibrated using typical microphone calibration procedures in an anechoic chamber to determine the dynamic response of each sensor individually. A B&K Type 4929-L OmniPower Sound Source was used to produce an omnidirectional sound source that was first measured with a B&K Type 4138 1/8ʹʹ microphone covered with a pinhole cap. This measurement acted as a reference with which to calculate the transfer function associated with each individual sensor. With appropriate time delay corrections to account for the variations in the sensor locations, both amplitude and phase corrections were determined for each sub-resonant sensor in the array. These corrections were smoothed using a polynomial fitting algorithm to account for uncertainties present in the measurements and experimental set-up. These corrections were applied to all data measured in the VTSWT to more accurately capture the wall-pressure field measured by the sensor array. Above 2 kHz, the calibration for each sensor began to diverge, which introduced uncertainty into the measurements, necessitating a rstriction to data below 2 kHz for all final data products. Notably, a cutoff frequency of 2 kHz is sufficient for capturing the range of influence of the pressures generated by the large-scale motions (LSMs) and very-large-scale motions (VLSMs) within the boundary layer for the flow conditions examined in this study.

Once the sensor array was installed in the VTSWT, measurements were made using a 16-bit data acquisition system and simultaneous sampling of all 64 channels. Data were acquired using a sampling rate of 25.6 kHz over a total sampling time of 256 s. To determine the auto-spectral density and cross-spectral density associated with each sensor and each sensor pair, Welch’s method was implemented for a record length of 4096 samples, a 50 % overlap and windowed using a Hanning window. This corresponded to a total of 3199 records per measurement, or 521 349 $\delta ^{*}/U_{\infty }$ for the baseline flow case. Additionally, each measurement was repeated six times to further improve averaging and aid in data convergence. Combining these six measurements yielded a total of 19 199 records per flow condition, or 3 128 091 $\delta ^{*}/U_{\infty }$ for the baseline flow case. Further discussion of the advantages of combining datasets during post-processing is given in § 4.3. In addition to the sub-convective pressure sensing array, measurements of the single-point wall-pressure spectrum were made using B&K Type 4138 1/8ʹʹ microphones with 0.5 mm diameter pinhole caps. Data were acquired using a sampling rate of 65.536 kHz for a sampling duration of 600 s. The auto-spectral density was calculated using Welch’s method for a record length of 4096 samples, a 50 % overlap and windowed using a Hanning window. This corresponded to a total of 19 199 records per measurement, just as with the measurements made with the sub-convective pressure sensing array.

3.1. Single-point wall-pressure spectrum models

This study will examine four models of the single-point pressure spectrum, of which three incorporate the effects of roughness. For a more in-depth comparison and discussion of these models, refer to Fritsch et al. (Reference Fritsch, Vishwanathan, Roy, Todd Lowe and Devenport2023).

3.1.1. Goody model

The Goody (Reference Goody2004) model is a widely used empirical model and forms the basis for the development of single-point models incorporating the effects of roughness. The original form documented for smooth-wall boundary layers is given in (3.1), where $R_{T}$ is the ratio between the outer layer and inner layer time scales defined as $ (u_{\tau }\delta /\nu )\sqrt {C_{f}/2}$ . This parameter is generally responsible for the roll off in the spectrum at high frequencies, where a larger value for $R_{T}$ corresponds to a slower roll off because of the larger range of spatial and temporal scales present in the boundary layer.

(3.1) \begin{align} \frac {G_{pp}(\omega )U_{e}}{\tau _{w}^{2}\delta } = \frac {3.0\left ( \omega \delta /U_{e}\right )^{2}}{\left [\left (\omega \delta /U_{e}\right )^{0.75}+0.5\right ]^{3.7} + \left [\left (1.1R_{T}^{-0.57}\right )\left (\omega \delta /U_{e}\right )\right ]^{7}} .\end{align}

Equation (3.1) can be generalised according to (3.2), where the constants $a{-}h$ are free to be altered, and the parameters $SS$ and $FS$ correspond to the spectral and frequency scaling, respectively.

(3.2) \begin{align} G_{pp}(\omega )\textit{SS} = \frac {a\left ( \omega \textit{FS}\right )^{b}}{\left [\left (\omega \textit{FS}\right )^{c}+d\right ]^{e} + \left [\left (fR_{T}^{g}\right )\left (\omega \textit{FS}\right )\right ]^{h}}. \end{align}

This form has been adapted in the models that follow to more accurately account for the effects of roughness in the boundary layer by controlling the magnitude and the roll-off at higher frequencies.

3.1.2. Rough-wall models

The general form of the Goody (Reference Goody2004) model has been adapted by Catlett et al. (Reference Catlett, Bryan, Chang, Hemingway and Anderson2021), Joseph, Devenport & Glegg (Reference Joseph, Devenport and Glegg2022) and Fritsch et al. (Reference Fritsch, Vishwanathan, Roy, Todd Lowe and Devenport2023) to incorporate the effects of roughness into its construction. Table 1 lists the values of the parameters for each of these models. Here, $k_{s}^{+}$ is the roughness Reynolds number, $U_{e}^{+}$ is the viscous-scaled edge velocity, $u_{\tau }$ is the friction velocity, $u_{\nu }$ is the shear friction velocity and $\overline {U}$ is the mean velocity in the boundary layer as described by Zagarola & Smits (Reference Zagarola and Smits1998). The Catlett et al. (Reference Catlett, Bryan, Chang, Hemingway and Anderson2021) model improves upon the general form of the Goody (Reference Goody2004) model by increasing the magnitude of the spectrum at its peak and steepening the roll-off at higher frequencies, which is more representative of the single-point spectrum in rough-wall boundary layers. However, as will be documented more thoroughly in § 4.2, the magnitude increase accounted for in this model is not sufficiently large for the surface in this study. The Joseph et al. (Reference Joseph, Devenport and Glegg2022) model further increases the peak magnitude and steepens the roll-off in the mid-frequency and high-frequency ranges, and the Fritsch et al. (Reference Fritsch, Vishwanathan, Roy, Todd Lowe and Devenport2023) model has the largest increase in the peak magnitude, with a similar roll-off in the mid-frequency and high-frequency ranges to that of the Joseph et al. (Reference Joseph, Devenport and Glegg2022) model. All three of these roughness focused models yield significantly different forms, indicating the continued need for additional high-quality measurements of the pointwise wall-pressure spectrum to be adequately verified. Section 4.5.1 will compare these models with the data collected in this study as well as to the data collected by Fritsch et al. (Reference Fritsch, Vishwanathan, Roy, Todd Lowe and Devenport2023).

Table 1.

Comparison between parameter values in the generalised Goody model for rough-wall boundary layer flows. Adapted from Fritsch et al. (Reference Fritsch, Vishwanathan, Roy, Todd Lowe and Devenport2023).

3.2. Wavenumber–frequency spectrum modelling

Currently, there does not exist an empirical model of the three-dimensional wavenumber–frequency spectrum of pressure fluctuations in rough-wall turbulent boundary layers. Certain models such as the Corcos (Reference Corcos1967), Chase (Reference Chase1980), Chase (Reference Chase1987), Hwang et al. (Reference Hwang, Bonness and Hambric2003) modified Corcos and Smol’yakov (Reference Smol’yakov2006) models have been developed for smooth wall turbulent boundary layers. However, as documented by Damani et al. (Reference Damani, Butt, Totten, Devenport and Lowe2025a ), most of these models do not accurately capture the magnitude of the spectrum at the convective peak, in the sub-convective range, or the frequency dependence of the convective ridge, although the Chase model has shown promise. It is therefore necessary to improve upon these models to more accurately capture the nature of the wavenumber–frequency spectrum. Additionally, these models do not account for the acoustic pressure fluctuations generated via the scattering of the hydrodynamic pressure fluctuations from the rough surface or the effects due to the roughness sublayer. However, this does not mean these models do not have any utility when beginning to build models of the wavenumber–frequency spectrum beneath rough-wall boundary layers.

Glegg (Reference Glegg2023) proposed a method of determining the wavenumber–frequency transform of the pressure sources resulting from roughness-induced scattering from any rigid rough surface, as given by (3.3). This equation holds true for rough surfaces of any arbitrary geometry and material, assuming the height of the roughness is much smaller than the acoustic wavelength. The terms $\tilde {\tilde {p}}_{h}(k_{1},k_{3},\omega )$ and $\tilde {\tilde {p}}_{s}(k_{1},k_{3},\omega )$ represent the wavenumber transform of the hydrodynamic pressure fluctuations at the surface and the wavenumber transform of the pressure fluctuations generated due to scattering, respectively. The term $H(\kappa _{1},\kappa _{3},k_{1},k_{3},\omega )$ controls how the pressure fluctuations generated at the surface defined in terms of wavenumbers $\kappa _{1}$ and $\kappa _{3}$ are scattered into acoustic and evanescent pressure waves at wavenumbers $k_{1}$ and $k_{3}$ . This function depends on the specifics of the roughness distribution responsible for inducing the scattering and takes the form given in (3.4). The term $\mu$ is given in (3.5) and the term $\tilde {\tilde {\xi }}(\kappa _{1}-k_{1},\kappa _{3}-k_{3})$ represents the wavenumber-shifted wavenumber transform of the surface, where $\xi (x_{1},x_{3})$ represents the local height of the surface.

(3.3) \begin{align} \tilde {\tilde {p}}_{s}(k_{1},k_{3},\omega ) = \int _{-\infty }^{\infty }\int _{-\infty }^{\infty }\tilde {\tilde {p}}_{h}(k_{1},k_{3},\omega )H(\kappa _{1},\kappa _{3},k_{1},k_{3},\omega )\,\text{d}\kappa _{1}\,\text{d}\kappa _{3}, \\[-28pt] \nonumber \end{align}

(3.4) \begin{align} H(\kappa _{1},\kappa _{3},k_{1},k_{3},\omega )=\left (\frac {i}{\mu }\right )\left (k_{1}\left (\kappa _{1}-k_{1}\right )+k_{3}\left (\kappa _{3}-k_{3}\right )\right )\tilde {\tilde {\xi }}(\kappa _{1}-k_{1},\kappa _{3}-k_{3}) , \\[-28pt] \nonumber \end{align}

(3.5) \begin{align} \mu = \sqrt {k_{1}^{2} + k_{3}^{2}-\left (\frac {\omega }{c_{\infty }}\right )^{2}}. \\[-12pt] \nonumber \end{align}

For a specific surface description, it is then possible to determine the wavenumber–frequency spectrum of the pressure fluctuations due to scattering. For this study, the rough surface that has been tested is constructed from a staggered pattern of cylindrical elements, as seen in figure 3. To begin, we can define a method for analytically describing the surface height from a single element, $\xi _{0}(x_{1},x_{3})$ , which can be done through the use of a cylinder function, $C(x_{1},x_{3})$ . Here, $k_{g}$ represents the height of the element, $r$ is the radius of the element and $(a_{0},b_{0})$ is the location of the centre of the element in the $x_{1}{-}x_{3}$ plane.

(3.6) \begin{align} \xi _{0}(x_{1},x_{3}) = k_{g}C(x_{1},x_{3},r,a_{0},b_{0}) = \begin{cases} k_{g} & \text{for} \def\lumina\hspace {0.5cm} \sqrt {(x_{1} - a_{0})^{2} + (x_{3} - b_{0})^{2}} \leqslant r, \\[5pt] 0 & \text{for} \def\lumina\hspace {0.5cm} \sqrt {(x_{1} - a_{0})^{2} + (x_{3} - b_{0})^{2}} \gt r .\\ \end{cases} \end{align}

By setting $(a_{0},b_{0})$ = (0, 0), we may find the shifted wavenumber transform of this function analytically, using the approach outlined by Weisstein (Reference Weisstein2024), where $J_{1}(\kappa _{1} - k_{1},\kappa _{3} - k_{3})$ is a Bessel function of the first kind.

(3.7) \begin{align} \tilde {\tilde {\xi }}_{0}(\kappa _{1} - k_{1},\kappa _{3} - k_{3}) = \frac {k_{g}}{4 \pi ^{2}}\int _{-\infty }^{\infty }\int _{-\infty }^{\infty } C \left (x_{1}, x_{3} \right )e^{i \left ( \kappa _{1} - k_{1} \right )x_{1} + i \left ( \kappa _{3} - k_{3} \right )x_{3}}\, \text{d}x_{1}\, \text{d}x_{3}, \\[-28pt] \nonumber \end{align}

(3.8) \begin{align} \tilde {\tilde {\xi }}_{0}(\kappa _{1} - k_{1},\kappa _{3} - k_{3}) = \frac {k_{g}rJ_{1}(r\sqrt {(\kappa _{1}-k_{1})^{2}+(\kappa _{3}-k_{3})^{2}})}{2 \pi \sqrt {(\kappa _{1}-k_{1})^{2}+(\kappa _{3}-k_{3})^{2}}}. \\[0pt] \nonumber \end{align}

The form given in (3.8) only describes the wavenumber transform of an element positioned at the origin. However, we may use this to obtain the result for an infinite collection of elements. By recognising that each element will have a transform that is shifted by the position of its centre, we may construct the wavenumber transform for the entire surface according to (3.9), where $(a_{j},b_{j})$ refers to the centre of each element in the rough surface.

(3.9) \begin{align} \tilde {\tilde {\xi }}(\kappa _{1} - k_{1},\kappa _{3} - k_{3}) = \left ( \frac {k_{g}rJ_{1} \left ( r \sqrt {(\kappa _{1}-k_{1})^{2}+(\kappa _{3}-k_{3})^{2}} \right )}{2 \pi \sqrt {(\kappa _{1}-k_{1})^{2}+(\kappa _{3}-k_{3})^{2}}} \right ) \sum _{j=0}^{\infty }e^{i \left ( \kappa _{1} - k_{1} \right )a_{j} + i \left ( \kappa _{3} - k_{3} \right )b_{j}}. \end{align}

For this surface, we may rewrite the infinite sum of complex exponential terms by using the relationship given in (3.10) and further simplifying it to the form given in (3.11). Here, $s$ refers to the pitch between elements, as shown in figure 3. Note that the expression given in (3.11) holds regardless of where the element corresponding to $(l,m)=(0,0)$ is in relation to the origin in the $x_{1}{-}x_{3}$ plane.

(3.10) \begin{align} & \sum _{j=0}^{\infty }e^{i \left ( \kappa _{1} - k_{1} \right )a_{j} + i \left ( \kappa _{3} - k_{3} \right )b_{j}} =\left ( \sum _{l=-\infty }^{\infty } e^{i\left ( \kappa _{1} - k_{1} \right )2ls} \right ) \left ( \sum _{m=-\infty }^{\infty } e^{i\left ( \kappa _{3} - k_{3} \right )ms} \right ) \nonumber \\&\quad + \left ( \sum _{n=-\infty }^{\infty } e^{i\left ( \kappa _{1} - k_{1} \right )\left (2ns + s\right )} \right ) \left ( \sum _{q=-\infty }^{\infty } e^{i\left ( \kappa _{3} - k_{3} \right )\left (qs + \tfrac {s}{2}\right )} \right ), \\[-12pt] \nonumber \end{align}

(3.11) \begin{align} &\sum _{j=0}^{\infty }e^{i \left ( \kappa _{1} - k_{1} \right )a_{j} + i \left ( \kappa _{3} - k_{3} \right )b_{j}} \nonumber\\&\quad =\left ( 1 + e^{i\left ( \kappa _{1} - k_{1} \right )s + i\left ( \kappa _{3} - k_{3} \right )\left (\tfrac {s}{2}\right )}\right )\left ( \sum _{l=-\infty }^{\infty } e^{i\left ( \kappa _{1} - k_{1} \right )2ls} \right ) \left ( \sum _{m=-\infty }^{\infty } e^{i\left ( \kappa _{3} - k_{3} \right )ms} \right ). \\[10pt] \nonumber \end{align}

Using the relation given in (3.12), we may recast this infinite sum of exponential terms as an infinite sum of delta functions, also referred to as a Dirac comb, to yield the relationship given in (3.13).

(3.12) \begin{align} \sum _{n=-\infty }^{\infty } e^{i\omega n \Delta t } = \frac {2 \pi }{\Delta t} \sum _{n = -\infty }^{\infty }\delta \left ( \omega -\frac {2\pi n}{\Delta t}\right ), \\[-35pt] \nonumber \end{align}

(3.13) \begin{align} &\sum _{j=0}^{\infty }e^{i \left ( \kappa _{1} - k_{1} \right )a_{j} + i \left ( \kappa _{3} - k_{3} \right )b_{j}} =\left ( 1 + e^{i\left ( \kappa _{1} - k_{1} \right )s + i\left ( \kappa _{3} - k_{3} \right ) \left (\tfrac {s}{2}\right )}\right ) \nonumber \\&\quad \left (\frac {4\pi ^{2}}{s^2}\right )\left ( \sum _{l=-\infty }^{\infty } \delta \left (\left ( \kappa _{1} - k_{1} \right )-\frac {\pi l}{s}\right ) \right ) \left ( \sum _{m=-\infty }^{\infty } \delta \left (\left ( \kappa _{3} - k_{3} \right )-\frac {2\pi m}{s}\right ) \right ). \\[-12pt] \nonumber \end{align}

We may now substitute the form given in (3.13) into (3.9) to obtain the final functional form of the wavenumber transform of the rough surface being studied.

(3.14) \begin{align} \tilde {\tilde {\xi }}(\kappa _{1} - k_{1},\kappa _{3} - k_{3}) & = \left ( \frac {2 \pi k_{g}rJ_{1} \left ( r \sqrt {(\kappa _{1}-k_{1})^{2}+(\kappa _{3}-k_{3})^{2}} \right )}{s^{2} \sqrt {(\kappa _{1}-k_{1})^{2}+(\kappa _{3}-k_{3})^{2}}} \right ) \nonumber \\[5pt] &\quad \left ( 1 + e^{i\left ( \kappa _{1} - k_{1} \right )s + i\left ( \kappa _{3} - k_{3} \right )\left (\tfrac {s}{2}\right )}\right ) \left ( \sum _{l=-\infty }^{\infty } \delta \left (\left ( \kappa _{1} - k_{1} \right )-\frac {\pi l}{s}\right ) \right ) \nonumber \\[5pt]& \quad \left ( \sum _{m=-\infty }^{\infty } \delta \left (\left ( \kappa _{3} - k_{3} \right )-\frac {2\pi m}{s}\right ) \right ). \end{align}

We may now use the result for the wavenumber transform of the rough surface to determine the relationship between $\tilde {\tilde {p}}_{s}(k_{1},k_{3},\omega )$ and $\tilde {\tilde {p}}_{h}(\kappa _{1},\kappa _{3},\omega )$ . First, we construct $H(\kappa _{1},\kappa _{3},k_{1},k_{3})$ .

(3.15) \begin{align} & H(\kappa _{1},\kappa _{3},k_{1},k_{3}) \nonumber \\& = \left ( \frac {2 i \pi k_{g}r \left ( k_{1} \left ( \kappa _{1} - k_{1} \right ) + k_{3} \left ( \kappa _{3} - k_{3} \right )\right ) J_{1} \left ( r \sqrt {(\kappa _{1}-k_{1})^{2}+(\kappa _{3}-k_{3})^{2}} \right )}{s^{2} \sqrt {(\kappa _{1}-k_{1})^{2}+(\kappa _{3}-k_{3})^{2}}\left ( \sqrt {k_{1}^{2} + k_{3}^{2} - \left ( {\omega }/{c_{\infty }} \right )^{2}} \right )} \right )\nonumber \\ & \quad \big ( 1 + e^{i\left ( \kappa _{1} - k_{1} \right )s + i\left ( \kappa _{3} - k_{3} \right )\left ({s}/{2}\right )}\big ) \left ( \sum _{l=-\infty }^{\infty } \delta \left (\left ( \kappa _{1} - k_{1} \right )-\frac {\pi l}{s}\right ) \right ) \nonumber \\& \quad \left ( \sum _{m=-\infty }^{\infty } \delta \left (\left ( \kappa _{3} - k_{3} \right )-\frac {2\pi m}{s}\right ) \right ). \end{align}

We may now find the relationship between the hydrodynamic pressure fluctuations and the pressure fluctuations due to scattering according to (3.16).

(3.16) \begin{align} &\tilde {\tilde {p}}_{s}(k_{1},k_{3},\omega ) =\int _{-\infty }^{\infty } \int _{-\infty }^{\infty }\tilde {\tilde {p}}_{h}(\kappa _{1},\kappa _{3},\omega ) \nonumber\\ &\quad \left ( \frac {2 i \pi k_{g}r \left ( k_{1} \left ( \kappa _{1} - k_{1} \right ) + k_{3} \left ( \kappa _{3} - k_{3} \right )\right ) J_{1} \left ( r \sqrt {(\kappa _{1}-k_{1})^{2}+(\kappa _{3}-k_{3})^{2}} \right )}{s^{2} \sqrt {(\kappa _{1}-k_{1})^{2}+(\kappa _{3}-k_{3})^{2}}\left ( \sqrt {k_{1}^{2} + k_{3}^{2} - \left ({\omega }/{c_{\infty }} \right )^{2}} \right )} \right ) \nonumber\\&\qquad\left ( 1 + e^{i\left ( \kappa _{1} - k_{1} \right )s + i\left ( \kappa _{3} - k_{3} \right )\left ({s}/{2}\right )}\right ) \left ( \sum _{l=-\infty }^{\infty } \delta \left (\left ( \kappa _{1} - k_{1} \right )-\frac {\pi l}{s}\right ) \right ) \nonumber\\&\quad \left ( \sum _{m=-\infty }^{\infty } \delta \left (\left ( \kappa _{3} - k_{3} \right )-\frac {2\pi m}{s}\right ) \right )\,\text{d}\kappa _{1} \,\text{d}\kappa _{3}. \end{align}

Evaluating this integral analytically allows for the determination of the wavenumber–frequency transform of the scattered pressures, $\tilde {\tilde {p}}_{s}(k_{1},k_{3},\omega )$ , to be determined exactly.

(3.17) \begin{align} &\tilde {\tilde {p}}_{s}(k_{1},k_{3},\omega ) =\sum _{l=-\infty }^{\infty }\sum _{m=-\infty }^{\infty }\tilde {\tilde {p}}_{h} \left(k_{1}+\frac {\pi l}{s},k_{3}+\frac {2 \pi m}{s},\omega \right) \nonumber \\ & \quad \left ( \frac {2 i \pi k_{g}r \left ( k_{1}l + 2k_{3}m\right ) J_{1}\left (\frac {r\pi }{s}\sqrt {l^{2}+4m^{2}}\right )}{s^{2} \sqrt {l^{2}+4m^{2}}\left ( \sqrt {k_{1}^{2} + k_{3}^{2} - \left ( {\omega }/{c_{\infty }} \right )^{2}} \right )} \right ) \big ( 1 + e^{i\pi \left (l+m \right )}\big ). \end{align}

The wavenumber–frequency spectrum is related to the wavenumber–frequency transform of any arbitrary variable according to (3.18), where $R_{\infty }$ and $T$ refer to the limits of integration in the integral definition of the spectrum. This allows for the direct determination of the wavenumber–frequency spectrum of the scattered surface pressure generated by the surface given in figure 3, as given in (3.19).

(3.18) \begin{align} \varPhi _{aa}\left ( k_{1},k_{3},\omega \right ) = \frac {\pi ^{3}}{R_{\infty }^{2}T} E \big [ \tilde {\tilde {a}}^{*} ( k_{1},k_{3},\omega )\tilde {\tilde {a}} ( k_{1},k_{3},\omega )\big ] , \\[-28pt] \nonumber \end{align}

(3.19) \begin{align} \varPhi _{ss}\left ( k_{1},k_{3},\omega \right ) = \left (\frac {8\pi ^{2}k_{g}^{2}r^{2}}{s^{4}}\right )\sum _{l = -\infty }^{\infty }\sum _{m = -\infty }^{\infty }\varPhi _{hh}\left (k_{1}+\frac {\pi l}{s},k_{3}+\frac {2\pi m}{s},\omega \right ) \\ \nonumber \left (\frac { \left (k_{1}l+2k_{3}m\right )J_{1}\left (\frac {r\pi }{s}\sqrt {l^{2}+4m^{2}}\right )}{\left (\sqrt {l^{2}+4m^{2}}\right )\left |\sqrt {k_{1}^{2} + k_{3}^{2} - \left ( {\omega }/{c_{\infty }} \right )^{2}}\right |}\right )^{2}\left (1+\cos {\left (\pi \left (l+m\right )\right )}\right ). \\[10pt] \nonumber \end{align}

It should be noted that the form given in (3.19) is undefined for $ (l,m ) = (0,0 )$ . However, the limit exists and can be shown to equal zero. Therefore, for the implementation of this model, all instances of $(l,m)=(0,0)$ can be ignored. While (3.19) is a closed, analytical form of the spectrum due to scattering from the surface given in figure 3, it does not come without drawbacks. Namely, calculating an infinite sum is infeasible, so this must be reduced to a more reasonable approximation. Further discussion of this is given in § 4.5.1.