Hostname: page-component-74656f9884-j7cmm Total loading time: 0 Render date: 2026-11-06T03:04:27.604Z Has data issue: false hasContentIssue false

Stochastic generation of velocity fields to reproduce the energy spectrum in wall turbulence

Published online by Cambridge University Press:  10 August 2026

Roozbeh Ehsani
Affiliation:
Saint Anthony Falls Laboratory, University of Minnesota, Minneapolis, MN 55414, USA Department of Civil, Environmental, and Geo-Engineering, University of Minnesota, Minneapolis, MN 55455, USA
Michael Heisel
Affiliation:
School of Civil Engineering, University of Sydney, Sydney, Australia
Michele Guala*
Affiliation:
Saint Anthony Falls Laboratory, University of Minnesota, Minneapolis, MN 55414, USA Department of Civil, Environmental, and Geo-Engineering, University of Minnesota, Minneapolis, MN 55455, USA
*
Corresponding author: Michele Guala, mguala@umn.edu

Abstract

Content of image described in text.

Leveraging experimental measurements in rough-wall turbulent boundary layers that provided the spatial distribution and statistics of uniform momentum zones (UMZs), we previously introduced a stochastic model to generate two-dimensional modal velocity fields representative of the roughness sublayer and logarithmic region (Ehsani et al. 2024a J. Fluid Mech., vol. 979, p. A12; J. Fluid Mech., 2024b, vol. 999, p. A56). This synthetic flow comprises a streamwise concatenation of spatially correlated step-like velocity profiles, where each step corresponds to an internal shear layer separating UMZs, and the UMZs follow scaling consistent with wall-attached eddies. Here, we further develop the model to introduce small-scale swirling motions using vortex cores generated stochastically from previously published statistics and placed opportunistically with respect to the internal shear layers. Strategies are discussed to stretch and distribute vortices and extend the spatial resolution of the synthetic field to capture the smallest flow scales. The new synthetic fields for both laboratory and atmospheric settings reproduce with reasonable accuracy the second-order turbulence statistics and the energy spectrum for most of the inertial range, capturing, for the wind tunnel datasets, the transition to exponential decay up to $\textit{k}_{1}\eta \sim 0.2$ (where $\textit{k}_{1}$ is the streamwise wave number and $\eta$ is the Kolmogorov length scale). The long-term goal is to extend and improve this bottom-up, statistically parameterised, reconstruction of rough-wall turbulence to open new wall modelling avenues for large eddy simulations.

Information

Type
JFM Papers
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited.
Copyright
© The Author(s), 2026. Published by Cambridge University Press
Figure 0

Figure 1. (a) Instantaneous streamwise velocity field u(x,z)$\textit{u}(\textit{x},\textit{z})$ for the PIV(m1) dataset. (b) The modal velocity field um(x,z)$\textit{u}_{\textit{m}}\textit{(x,z)}$ corresponding to (a). (c) Probability density function of the UMZ’s thickness at wall-normal z/δ=0.05$\textit{z}/\delta = 0.05$. (d) Generation of the step-like velocity profile by employing a random number randhmi$\textit{rand}_{\textit{h}_{\textit{m}_{i}}}$ and inverse sampling from the cumulative distribution function (c.d.f.) of the thickness. The same methodology is implemented to determine the modal and vertical velocity components of the generated zone (step). (e) Generation of a correlated modal velocity field by reorganising and spacing the stochastically generated modal velocity profiles.

Figure 1

Table 1. Experimental datasets used in the detection and statistical analysis of UMZ characteristics. The generated velocity fields are compared against these datasets. The Kolmogorov length scale η$\eta$ and Taylor micro-scale λT$\lambda _{\textit{T}}$ are estimated at z/δ=0.02$\textit{z}/\delta =0.02$. For each dataset, the aerodynamic roughness length z0$\textit{z}_{0}$ can be estimated as z0=ksexp⁡(−8.5κ)≈0.036ks$\textit{z}_{0} = \textit{k}_{\textit{s}}\exp {(-8.5\kappa )} \approx 0.036 \textit{k}_{\textit{s}}$, where κ=0.39$\kappa = 0.39$ is the Von Kármán constant.Table 1 long description.

Figure 2

Table 2. Generated low-resolution velocity fields in Ehsani et al. (2024b). The vertically averaged λT$\lambda _{\textit{T}}$ is calculated for each LR-Gen modal velocity field to generate vortex cores, and it is kept invariant across all elevations. The values of ν$\nu$, Reτ$_{\tau }$, uτ$u_{\tau }$, δ$\delta$, ks$_{\textit{s}}$ and z0$_{0}$ for the generated dataset are equal to the corresponding values in the experimental dataset. Here, λT$\lambda _{\textit{T}}$ is then recalculated after the addition of vortices (see § 4).Table 2 long description.

Figure 3

Figure 2. (a) Quiver plot of the prograde Oseen vortex model. (b) Quadrant-based analysis of the Oseen vortex. Since the correlation coefficient ρu,wVorX$\rho _{\textit{u},w}^{\textit{VorX}}$ between uVorX$\textit{u}^{\textit{VorX}}$ and wVorX${w}^{\textit{VorX}}$ components of the Oseen vortex is zero, the induced Reynolds shear stress by the model is also zero. (c) Quiver plot of the transformed Oseen vortex model. (d) Quadrant-based analysis of the transformed Oseen vortex model, characterised by ρu,wVorX=−0.40$\rho _{\textit{u},w}^{\textit{VorX}}= -0.40$.

Figure 4

Figure 3. (a) The p.d.f. of the correlation coefficient between streamwise and vertical velocities within each prograde vortex core for the PIV wind tunnel datasets using Galilean invariant decomposition. The inset plot shows the profile of the mean values of the correlation coefficient. The influence of the roughness element persists up to z=2ks$ \textit{z} = 2\textit{k}_{\textit{s}}$; vortices in the roughness layer (z<2ks$\textit{z} \lt 2\textit{k}_{\textit{s}}$; white-coloured symbols with black (m1) and grey (m2) edges) contribute more to the Reynolds shear stress than vortices in the logarithmic region (z⩾2ks$\textit{z} \geqslant 2\textit{k}_{\textit{s}}$; coloured symbols). The p.d.f.s are fitted using skewed Gaussian distributions for the regions z<2ks$ \textit{z} \lt 2\textit{k}_{\textit{s}}$ (yellow solid line) and z⩾2ks$ \textit{z} \geqslant 2\textit{k}_{\textit{s}}$ (maroon dash-dotted line). (b) The p.d.f. of the radius of the prograde vortex cores at different wall-normal locations normalised by the Taylor micro-scale λT$\lambda _{\textit{T}}$, including fits using the log-normal-Pareto distribution. (c) The p.d.f. of the maximum azimuthal velocity for prograde vortex cores at different elevations normalised by the friction velocity, including fits using exponential distributions. Parameters for each fitted distribution are provided in table 3.

Figure 5

Table 3. Parameters of the p.d.f. for vortex features. Rows 1 and 2 correspond to the correlation coefficient, rows 3 and 4 to the vortex core radius and rows 5 and 6 to the maximum azimuthal velocity of the vortex core. All the parameters are empirical and derived using the MLE method.Table 3 long description.

Figure 6

Table 4. Parameters of the inverse c.d.f. for vortex features. Rows 1 and 2 correspond to the correlation coefficient, rows 3 and 4 to the vortex core radius and rows 5 and 6 to the maximum azimuthal velocity of the vortex core. All the parameters are empirical and derived using the MLE method.Table 4 long description.

Figure 7

Figure 4. Snapshot of the (a) swirling strength |λci(x,z)|>γλci(r.m.s.)(z)$\lvert \lambda _{\textit{ci}}(\textit{x,z})\rvert \gt \gamma \lambda _{\textit{ci}(\textit{r.m.s.})}(\textit{z})$ and (b) spanwise vorticity |ωy(x,z)|>γωy(r.m.s.)(z)$\lvert \omega _{\textit{y}}(\textit{x,z})\rvert \gt \gamma \omega _{\textit{y}(\textit{r.m.s.})}(\textit{z})$ fields computed on the high-resolution generated modal velocity field, where γ=0.35$\gamma = 0.35$. The sign of the values in both fields is calculated using the sign of the local spanwise vorticity. Both quantities are normalised with their respective maximum absolute values for the Gen(m1) dataset.

Figure 8

Figure 5. Distribution of the average value of vorticity in swirling strength clusters (shear layers) in the HR-SGVF in comparison with the distribution of vorticity of the generated vortices for the m1 dataset. Both p.d.f.s are normalised with the corresponding r.m.s. values. Bar plots show the deciles of the distributions. To place vortices according to the strength of the shear layer, a decile matching process is implemented.

Figure 9

Figure 6. Schematic of different types of added vortices to HR-SGVF. (a–c) Are samples of vortices placed at shear layers with |λci(x,z)|>γλci(r.m.s.)(z)$\lvert \lambda _{\textit{ci}}(\textit{x,z})\rvert \gt \gamma \lambda _{\textit{ci}(\textit{r.m.s.})}(\textit{z})$, while (d) shows samples of near-wall vortices. The radii rω$\textit{r}_{\omega }$, azimuthal velocity uω$\textit{u}_{\omega }$, and the correlation coefficient between streamwise and vertical velocity components of the vortices ρu,wVorX$\rho _{\textit{u},w}^{\textit{VorX}}$ are generated stochastically, while the direction of rotation is given by the sign of the local vorticity. (a) A high-λci$\lambda _{\textit{ci}}$ cluster with sgn(ωy)=+1$\textit{sgn}(\omega _{\textit{y}}) = +1$: a primary retrograde vortex is placed at the centroid of the cluster, where the vortex footprint is given by 2rω$2 \: \textit{r}_{\omega }$ as shown in the figure 2. Secondary vortices are then placed at the same wall-normal location to fill the streamwise extent of the cluster. Secondary vortices whose xcentroid$\textit{x}_{\textit{centroid}}$ extends beyond the cluster margins are not added. (b) A high-λci$\lambda _{\textit{ci}}$ cluster with sgn(ωy)=−1$\textit{sgn}(\omega _{\textit{y}}) = -1$ and prograde vortices. A secondary vortex was not added to the right side of the primary vortex because the generated vortex extended beyond the right margin. (c) The primary vortex exceeds the high-λci$\lambda _{\textit{ci}}$ cluster margins, leaving no space for secondary vortices. (d) Square-shaped near-wall vortices with an edge length of amax=2rω$a_{\textit{max}} = 2\:\textit{r}_{\omega }$. A randomly generated near-wall vortex is placed if uHR-SGVF′>0$\textit{u}^{\prime }_{\textit{HR-SGVF}}\gt 0$ at the centre of the vortex, and the randomly generated binary number is randbinary∼B(0.5)=1$\textit{rand}_{\textit{binary}} \sim \mathcal{B}(0.5) = 1$. A quiver plot of the generated near-wall vortex (2) is included in the inset. The primary effect of the near-wall vortices is to reduce the mean streamwise velocity and increase the Reynolds shear stress along the bottom of the vortex, with minimal impact on the field away from the bottom.

Figure 10

Figure 7. (a) Streamwise velocity contour plot of the HR-SGVF after adding vortex cores (HR + VorX-Gen(m1)). (b) Experimental wind tunnel (m1) dataset. The velocity components are normalised with the inner velocity scale uτ$u_{\tau }$. The convection velocity of uconv+=8$\textit{u}^{+}_{\textit{con}v} = 8$ and uconv+=11$\textit{u}^{+}_{\textit{con}v} = 11$ is subtracted from the generated and experimental datasets, respectively, to show the fluctuating velocity field. The top left green arrows, which correspond to 7.8 u+$\textit{u}^{+}$, are provided for reference.

Figure 11

Figure 8. The mean streamwise velocity profile for different datasets (m1, m2 and ASL) is shown. The yellow solid line is the theoretical log law U¯/uτ=1/κln⁡(z/z0)$\overline {\textit{U}}/u_{\tau } = 1/\kappa \ln (\textit{z}/\textit{z}_{0})$ where κ=0.39$\kappa = 0.39$. Due to the limitations of the histogram-based approach for detecting UMZs near the wall, we were unable to accurately reproduce the mean streamwise velocity profiles close to the wall in the LR-Gen dataset (Ehsani et al.2024a). The effect of the near-wall vortices on the correction of the mean streamwise velocity profile is evident in the HR + VorX-Gen fields (i.e. 15⪅z/z0⪅30$15\lessapprox \textit{z}/\textit{z}_{0}\lessapprox 30$ for HR + VorX-Gen(m1,m2) and 350⪅z/z0⪅600$350\lessapprox \textit{z}/\textit{z}_{0}\lessapprox 600$ for HR + VorX-Gen(ASL)).

Figure 12

Figure 9. Profiles of the second moment of statistics for different datasets. Dataset HR + VorX-Gen is compared with LR-Gen (Ehsani et al.2024a) and the actual dataset. Increasing the resolution of the generated field from LR-Gen to HR-SGVF attenuates the second moment of statistics while preserving the trend observed in LR-Gen. The effect of different types of vortices at different wall-normal elevations on the statistics can be studied separately. (a) Streamwise velocity variance profile. The inset plot shows the wall-normal elevation non-dimensionalised by δ$\delta$. The theoretical profile (yellow line), as proposed by Townsend (1976), follows: (u′u′¯/uτ2)=−Aln⁡(z/δ)+B$({\overline {u^{\prime }u^{\prime }}}/{u_{\tau }^{2}}) = -A\ln {(\textit{z}/\delta )}+B$, where A=1.45$A = 1.45$, and B=1.57$B = 1.57$. (b) Vertical velocity variance profile. (c) Reynolds shear stress (covariance) profile.

Figure 13

Figure 10. Second-order structure function for the generated m1 dataset. The yellow line (slope = 2.9) represents the logarithmic D11$\textit{D}_{11}$ law, signature of the attached eddies, where the slope equals 2A$2A$ (A$A$ being the Townsend Perry constant). Inset: power spectrum displaying the equivalent k1−1$\textit{k}_{1}^{-1}$ scaling for attached-eddy structures.

Figure 14

Figure 11. The streamwise velocity spectra from the wind tunnel datasets are compared with those computed on the increased resolution field (HR-Gen) and HR + VorX-Gen to show the effect of near-wall and shear layer vortices. (a) For the m1 dataset at z/δ=0.02$\textit{z}/\delta =0.02$ and (b) for m2 dataset at z/δ=0.04$\textit{z}/\delta =0.04$ (relatively close to the wall).

Figure 15

Figure 12. The streamwise velocity spectra from the wind tunnel datasets are compared with those computed on the increased resolution field (HR-Gen) and HR + VorX-Gen to show the effect of added shear layer vortices. The reduced number of shear layers at higher elevations results in a slight underestimate of turbulent kinetic energy at the higher end of the inertial subrange. (a) For the m1 dataset at z/δ=0.12$\textit{z}/\delta =0.12$ and (b) for m2 dataset at z/δ=0.13$\textit{z}/\delta =0.13$ (relatively far from the wall).

Figure 16

Figure 13. The streamwise velocity power spectral density for ASL datasets. (a) The near-wall and shear layer added vortices effectively reproduce energy at z/δ=0.01$\textit{z}/\delta =0.01$ for HR + VorX-Gen(ASL). (b) The energy spectrum of HR + VorX-Gen at z/δ=0.03$\textit{z}/\delta =0.03$ moderately underestimates the inertial subrange decay. This is due to the use of the vortex radius distribution derived from the wind tunnel dataset, which contains relatively small vortex sizes. Moreover, the insufficient number of shear layers at higher elevations limits the addition of vortices, leading to inadequate resolution of the energy spectrum at the higher elevations.

Figure 17

Figure 14. Profiles of added vortex densities per λT2$\lambda _{\textit{T}}^{2}$ for the generated wind tunnel (m1) and ASL datasets compared with the experimental PIV(m1) dataset. For the experimental dataset, the density of swirling strength clusters per λT2$\lambda _{\textit{T}}^{2}$ is estimated. (a) Number of prograde shear layer vortices, (b) retrograde shear layer vortices and (c) prograde, square-shape, near-wall vortices.

Figure 18

Figure 15. (a) Streamwise fluctuating velocity signal for different velocity fields close to the wall. The velocity generated in Ehsani et al. (2024b) is denoted as LR-Gen. Here, HR-Gen is the high-resolution velocity field using the MAkima method without adding any types of vortices. The resolution of LR-Gen(m1) is increased 32 times. Dataset HR+VorX-Gen is the high-resolution velocity field after adding vortices without applying the Gaussian filter and HR+VorX+Filter-Gen is the high-resolution generated velocity field after adding vortices and passing the Gaussian filter. The filter just smooths out the sharp edges, which are due to vortex addition. (b) Power spectral density for different velocity fields. Passing the filter dampens the energy at the high end of the wavenumber (dissipation range) and has a viscous effect.

Figure 19

Figure 16. Distribution of prograde and retrograde vortex parameters at different wall-normal elevations. Curves are fitted to the distribution of prograde vortex parameters, which are then used in the vortex generation. The details about the p.d.f. fitting equations are summarised in (3.4), (3.5) and table 3. (a) The p.d.f. of the radius normalised with the Taylor micro-scale, with log-normal-Pareto distributions employed for fitting. (b) The p.d.f. of the maximum azimuthal velocity normalised with the friction velocity, with exponential functions fitted to the distribution.

Supplementary material: File

Ehsani et al. supplementary movie

Visualization of the high-resolution, stochastically generated, instantaneous velocity field for the wind tunnel m1 dataset. The contour plot represent the streamwise velocity component normalized with the shear velocity scale u*= 0.39 [m/s]. Both spatial coordinates are normalized with boundary layer thickness δ= 0.4 [m]. A convection velocity u+= 8 has been subtracted from the streamwise flow field to depict the velocity fluctuating components in a vectorial, quiver plot.
Download Ehsani et al. supplementary movie(File)
File 45.2 MB