Mean-velocity profile of smooth channel flow explained by a cospectral budget model with wall-blockage

Wall-bounded turbulent flows are pervasive in a variety of domains, including the environment, biology, and industry. Yet even for one of the simplest, canonical cases—smooth, pressure-driven channel flow—their bulk properties have defied simplified theoretical explanation, relying predominantly on dimensional analysis or similarity arguments. One such property is the mean-velocity profile (MVP or U(y), where u = U + u′ is the streamwise velocity, u′ is the turbulent excursion, and y is the distance from the wall). von Karman¹ and Prandtl² first suggested the existence, at sufficiently high Reynolds number (Re), of a logarithmic region in the MVP away from the wall³ such that

where κ is the von Karman constant, B is a constant that varies with wall roughness, and the + subscript signifies normalization by wall units, i.e., the velocity and length scales u_τ and $ν 0 u τ$⁠, respectively, where u_τ is the friction velocity and ν₀ is the kinematic molecular viscosity. Closer to the wall, in the “viscous layer,” viscous stresses dominate turbulent stresses and U₊ is approximately linear in y₊, with a “buffer layer” linking the linear and logarithmic regions. However, other theories, based on “incomplete similarity” reasoning for finite Reynolds number, dispute the existence of a logarithmic layer and propose its replacement with a power-law (e.g., Ref. 4).

A framework is required to explain the full MVP that moves beyond similarity arguments. Klewicki and co-authors^5–7 have proposed an alternative, four-layer structure based on relative magnitudes of terms in the momentum budget, as observed in direct numerical simulation (DNS) and experimental data. Another potentially useful route is to consider the role of turbulent statistical fluctuations in generating macroscopic phenomena such as the MVP. L’vov et al.⁸ proposed a model based on the characteristic length-scales of energy-containing eddies and the Reynolds-averaged Navier Stokes equations for momentum, turbulent kinetic energy, and Reynolds stress. The variation of these length-scales over the flow profile was estimated from DNS data, resulting in a model with three parameters, independent of Reynolds number. Other studies have taken advantage of known scaling laws for turbulent fluctuations instead. Turbulent fluctuations follow Kolmogorov’s universal scaling at scales much smaller than the flow’s integral length scale, but much larger than its viscous dissipation scale.⁹ This scaling can be recovered by externally exciting the Navier-Stokes equations with a Gaussian random force with carefully chosen correlation structure. In previous studies, the dynamic renormalization group (RNG) method has been applied to obtain macroscopic flow properties from the randomly forced Navier-Stokes equations, directly linking the micro- and macro-scales of the flow.¹⁰ 

In a similar vein, two recent theories have proposed a “spectral link” between the MVP and the spectrum of turbulent fluctuations F, the Fourier transforms of the two-point velocity correlation functions. Gioia et al.¹¹ first proposed a mechanism for such a link, relying on a heuristic argument in which turbulent stresses were generated by the product of a longitudinal velocity excursion set in size by “dominant eddies” of radius y₊, and a vertical velocity excursion determined by the turbulent kinetic energy $TKE ( y + ) = 1 2 u i ′ 2 ¯$⁠, which was in turn obtained from the TKE spectrum F_TKE(k). This phenomenological model has been shown to yield a range of known macroscopic properties of turbulent flows. These include Manning’s formula for bulk streamwise velocity in rough channel flow¹¹ and flows over submerged aquatic canopies,¹² the Reynolds number-friction factor relation in turbulent rough-impermeable pipe flow¹³ and gravel-bed (i.e., rough-permeable) flow encountered in the hyporheic zone,¹⁴ the MVP of turbulent pipe and channel flows¹¹ and the sheared, thermally stratified atmospheric surface layer (ASL),^15,16 and the mean temperature profile of the ASL.¹⁷ 

An alternative mechanism is proposed by Katul and Manes¹⁸ (hereafter KM14): a simplified cospectral budget, with a parameterized pressure-strain correlation term, that links the MVP to F_vv(k). There are two main advantages of this approach over that of Gioia et al.¹¹ First, it includes contributions to the turbulent stress from eddies of all sizes, not just the “dominant eddies.” Second, it links the MVP to the vertical velocity variance spectrum F_vv(k) rather than F_TKE(k), consistent with the production term in the turbulent stress budget. This approach has also been used to explain macroscopic relations between the turbulent Prandtl number, Monin-Obukhov stability parameter, and Richardson number in the stratified ASL.^19–21 

The modeled MVP in KM14 matches experimental observations well in the logarithmic region (see also Ref. 22) and the viscous sublayer, where turbulent stresses are small. However, the model fails to adequately capture the strong curvature present in the buffer layer. While a number of hypotheses were proposed for the missing curvature in the buffer layer, they could not be adequately tested with the available experimental data. We use a DNS channel flow dataset at moderate Reynolds number (⁠$Re τ = u τ R ν 0 =2003$⁠, where R is the channel half-width²³) to individually test several hypotheses explaining the curvature in the buffer layer, and use the outcome to propose a new model linking the spectra F_vv to the MVP across its full range, including the buffer layer. The new physics learned from the cospectral budget model is that the distribution of turbulent vertical velocity fluctuations (the “microstate” of the flow, represented by F_vv(k)) contains sufficient information to generate the MVP (the “macrostate” of the flow). This establishes a link between two previously unrelated areas of the turbulence literature: (1) Kolmogorov’s theory of homogeneous, isotropic turbulence, including the k^−5/3 scaling of the energy spectrum and (2) the “law of the wall” in wall-bounded turbulence. The specific value of this study is in testing and refining assumptions from previous studies by using DNS data, and making modifications to the theory where conflicts arise with the DNS data.

As a starting point, the salient features of the cospectral budget model are reviewed. The once-integrated streamwise momentum budget^3,24 is given by

where F_uv+(y₊, k₊) is the cospectrum of the Reynolds stress at height y₊ and wavenumber k₊. Assuming homogeneity in the spanwise direction, the Reynolds stress budget is

where $ρ ¯$ is the mean fluid density and p′ is the turbulent excursion from the mean pressure. Under stationary conditions, the storage term is ignored. Furthermore, the turbulent transport, pressure transport, and molecular diffusion terms provide relatively minor contributions, although they become more significant near the wall;³ we neglect them for now and return to this assumption later using the DNS results.

A parameterization is required for the pressure strain rate correlation ϕ that is a function of terms in the momentum and Reynolds stress budgets, to allow closure of the equations. This term is traditionally modeled as the sum of three components: slow, rapid, and wall-blocking. The slow component, generated by turbulent fluctuations, is usually modeled using Rotta’s linear return-to-isotropy form.²⁵ The rapid part, generated by the mean rate of strain, is often modeled as proportional to the production term P_uv.²⁶ For simplicity, the wall-blocking component was neglected in KM14 and other previous studies. Combining these components, the resulting parameterization is the Launder-Reece-Rodi Isotropization of Production (LRR-IP) model,³ 

where C_R ≈ 1.8 is the Rotta constant that depends on the precise definition of T(y), C_I = 3/5 is the isotropization constant predicted by Rapid Distortion Theory³ and T(y) is a relaxation time scale. A spectral form of this parameterization is proposed in KM14,

where $P u v ( k , y ) =− F v v ( k , y ) ∂ U ∂ y ( y )$ is a spectral form of the production term, and the relaxation time scale is τ(k, y) = ϵ(y)^−1/3min(k^−2/3, k_a(y)^−2/3), where k_a(y) = 1/y is a cutoff wavenumber used in an idealized form of the F_vv spectrum,²⁷ described shortly. Noting that $u ′ v ′ ¯ =sgn ( u ′ v ′ ¯ ) ∫ 0 ∞ F u v ( k ) dk$⁠, $v ′ 2 ¯ = ∫ 0 ∞ F v v ( k ) dk$⁠, and $ϵ u v = ν 0 sgn ( u ′ v ′ ¯ ) ∫ 0 ∞ k 2 F u v ( k ) dk$⁠, and substituting these expressions into the cospectral budget corresponding to Equation (3) yields

Rearranging and nondimensionalizing to wall units yields an explicit expression for F_uv,

To close the system of equations, a simplified TKE budget is introduced, where it is assumed that TKE production balances dissipation,

This approximation is reasonable in much of the logarithmic layer, but is strongly violated closer to the wall,³ and is tested below. Combining Equations (2), (6), and (7), $∂ U + ∂ y +$ (and, therefore, the MVP) can be solved in terms of F_vv+.

An idealized functional form for F_vv+ is proposed in KM14, and is shown in Fig. 1(a) (inset). This form consists of (i) the well-known $C 0 ϵ + 2 / 3 k + − 5 / 3$ Kolmogorov scaling⁹ in the inertial range (where C₀ is the Kolmogorov constant), (ii) a high-wavenumber dissipative range correction, (iii) a breakpoint at k_a+ = 1/y₊ and transition to $k + 0$ scaling, consistent with Townsend’s attached-eddy hypothesis,²⁸ and (iv) a low-wavenumber correction (Fig. 1(a), inset). The bottleneck at the crossover from inertial to viscous ranges is ignored and its effects were shown to be minor by KM14 when using the Meyers and Meneveau²⁹ formulation.

The missing curvature in the modeled MVP buffer layer (evident in Fig. 1(a)) may be due to the violation of one or more of the following assumptions in KM14: (1) balance between production and dissipation in the TKE budget, (2) the idealized shape of the F_vv+ spectrum used, and/or (3) use of the spectral LRR-IP parameterization for the pressure-strain term in the cospectral budget. These assumptions are addressed sequentially.

First, the assumption of a local balance between TKE production and dissipation is tested by forcing the model with DNS-reported dissipation rather than estimating it from Equation (7). The resulting profile is given in Fig. 1(b). The only discernible effect of this change is in the buffer layer, where the mean velocity is slightly increased. The reason for this change is evident from the profile of the production-to-dissipation ratio obtained from the DNS (inset of Fig. 1(b)). Around y₊ = 12, the production term is approximately 1.75 times the dissipation rate. However, as y₊ increases, the imbalance quickly fades to tolerably small values. Therefore, while this assumption introduces a small bias in the buffer layer, it quickly fades and does not play a significant role in the missing curvature. A possible explanation for the lack of sensitivity of the MVP to the imbalance between production and dissipation in TKE is due to F_uv+ scaling as $ϵ + 1 / 3$⁠, a sub-unity exponent.

Second, in addition to forcing the cospectral budget model with the DNS dissipation rate, we test the effect of the assumed idealized form of the F_vv+ spectra by forcing the model with the DNS-reported F_vv+ spectra. We approximate the three-dimensional F_vv+ spectra with their streamwise one-dimensional equivalents F_vvx+, effectively assuming isotropy. The resulting profile (Fig. 1(c)) overestimates the DNS-observed MVP. Directly comparing the DNS spectra with the idealized spectra (Fig. 2) reveals that the idealized form consistently overestimates F_vv+ for most wavenumbers. This results in a consistently overestimated vertical velocity variance throughout the profile (Fig. 3), which was also observed in KM14 when comparing with several experiments.

Substituting Equation (6) into Equation (2) and rearranging yields

It can be seen that correcting the overestimate (i.e., systematically reducing F_vv+) results in an increase in $∂ U + ∂ y +$ and the MVP. The mismatch between the DNS and modeled spectra is due to both the scaling and shape of the idealized form in Fig. 1(a) (inset). Inertial-range Kolmogorov scaling, as assumed in the idealized form for k₊ > k_a+, only applies when the spectral separation between y₊ and the Kolmogorov microscale is sufficiently large. Given that the DNS is at a moderate Reynolds number (Re_τ = 2003), it is not surprising that the Kolmogorov inertial-range scaling tends to result in overestimated F_vv+. Furthermore, the DNS spectra exhibit additional regimes not captured by the current idealized form. For y₊ < ∼ 120, the cutoff wavenumber k_a+ appears to remain roughly constant rather than varying with height as $k a + ∼ 1 y +$ assumed in KM14. It is only for y₊ > ∼ 120 that the cutoff wavenumber appears to begin varying with height. The invariance of k_a+ near the wall has been noted in previous studies (e.g., Ref. 30) and is possibly the result of streaks within the buffer region that are known to not scale with distance from the wall.³¹ The four-layer model proposed by Klewicki and co-authors^5–7 proposes a self-similar region governed by a hierarchy of length-scales linearly related to y₊, spanning $y + =2.6 Re τ$ to y = 0.5R. For this DNS dataset, the self-similar region (and associated scaling with y₊) is predicted to begin at $y + =2.6× 2003 =116$⁠, in good agreement with the observed height at which k_a+ begins varying with y₊ in the DNS dataset. In addition, there appear to be non-negligible contributions from a low-wavenumber mode (around k_x+ = 2 × 10⁻⁴ for y₊ > ∼ 120) in the DNS spectra that are not accounted for in the idealized form. Given the sensitivity of low-wavenumber components of calculated spectra to biases caused by smoothing and averaging, it is unclear if this is an artifact or a genuine source of variability. Since the Reynolds number of the DNS is still relatively low, we do not attempt to devise a parameterization for the DNS F_vvx+ spectra to avoid overfitting to potentially Re-dependent features. We leave this to future studies.

The fact that the modeled MVP degrades when forced with the DNS F_vvx+ and dissipation rates suggests errors caused by the third assumption (the parameterization of the pressure-strain term) must be compensating. To test this, pressure-strain profiles from the DNS are compared to those from the LRR-IP model (Fig. 4(a)). The model significantly overestimates the DNS-reported pressure-strain profile. In the simplified cospectral budget (Equation (5)), the production term effectively balances the pressure-strain term (since the dissipation term is small). Therefore, overestimating the pressure-strain term results in a negatively biased production term. To achieve this, since F_vv+ is prescribed, the model overestimates $∂ U + ∂ y +$ and the MVP.

Clearly, the pressure-strain correlation parameterization needs revising. The missing wall-blocking component of the pressure-strain rate correlation can be included using a standard parameterization^26,32

where L(y) = TKE(y)^3/2/ϵ(y) is a turbulent length scale. In this formulation, near the wall, the wall-blocking term increases the decorrelation effect of the slow pressure-strain component, and decreases the effect of the rapid component. A spectral version of Equation (9) may be proposed, given as

where $λ ( k , y ) =min ( 1 k , 1 k a ( y ) )$ is a spectral length-scale. Choosing $C R ′ =0$ and $C I ′ =0.6$ for the new parameters, the modeled pressure-strain profile $∫ 0 ∞ π ( k , y ) dk$ fits the DNS pressure-strain profile reasonably outside the viscous region (Fig. 4(b)). Furthermore, it fits the sum of the DNS pressure-strain and excluded transport/diffusion terms even better, perhaps justifying the neglect of these terms in Equation (3). A desirable feature of this pressure-strain model is that, in the intermediate region (η ≪ k⁻¹ ≪ y, where η = (ν³/ϵ)^1/4 is the Kolmogorov microscale), $λ ( k , y ) y = k y − 1 ≪1$⁠, so it reduces to the LRR-IP model. The LRR-IP model has been shown to generate good estimates of the inertial range cospectrum $F u v = C u v ∂ U + ∂ y + ϵ + 1 / 3 k − 7 / 3$⁠, including the constant C_uv ≈ 0.15.^18,22

The resulting modeled MVP (Fig. 1(d)) fits the DNS MVP well across the entire boundary-layer depth. The fit is particularly good given that errors in estimating $∂ U + ∂ y +$ accumulate in the estimated MVP as y⁺ increases. The only “tuning” that occurs is in developing the spectral pressure-strain model. Here, the parameters $C R ′$ and $C I ′$ are obtained by manually minimizing two objective functions: (1) the difference between $∫ 0 ∞ π ( k , y ) dk$ and the DNS pressure strain profile, using the DNS ϵ, F_uv and $∂ U ∂ y$ fields as inputs to Equation (11), (2) the difference between $∫ 0 ∞ F u v ( k , y ) dk$ and the DNS $u ′ v ′ ¯$ profile, using the DNS ϵ, F_vv and $∂ U ∂ y$ fields as inputs. No tuning is performed to match the modelled MVP to the DNS MVP, further emphasizing the impact of the form of the turbulent spectrum on the MVP.

The modeled and DNS cospectra F_uv+(k₊, y₊) are compared in Fig. 5. While the scale of the modeled F_uv+ is reasonable across the profile, the shapes vary significantly. For (y₊ > ∼400), the fit in the inertial range steadily improves and matches the DNS F_uvx+ well, consistent with other arguments.^18,22 However, a low wavenumber mode (around k₊ = 5 × 10⁻⁴) dominates the DNS F_uvx+ cospectrum but is not present in the modeled cospectrum. The low-wavenumber mode may be an important source of missed variability. It may also be a purely numerical artifact of smoothing the estimated spectra. To further examine these different explanations, the one-dimensional streamwise and spanwise vertical velocity variance, horizontal velocity variance, and associated cospectra are shown in Figs. 6 and 7, respectively. The horizontal velocity variance spectra F_uux+ are greater than the corresponding vertical velocity variance spectra F_vvx+ across most wavenumbers, consistent with the observed profiles of horizontal and vertical velocity variance (Fig. 3). An upper bound on the cospectrum can be obtained from the Cauchy-Schwarz inequality $| F u v + ( k + ) |≤ F u u + ( k + ) F v v + ( k + )$ (see, e.g., Ref. 33). This serves as a useful check in determining if the low wavenumber mode in the DNS F_uvx+ cospectrum is anomalous. As shown in Fig. 6 (bottom row), the low wavenumber mode in the compensated DNS cospectrum exceeds the upper bound given by the Cauchy-Schwarz inequality. Therefore, the low wavenumber mode appears to be an artifact of the DNS data, possibly due to the windowing and tapering process, known to particularly impact estimates of the lowest wavenumbers of the estimated (co)spectrum. This is further confirmed by examining the spanwise spectra (Fig. 7): the spanwise DNS F_uvz+ satisfies the Cauchy-Schwarz inequality and fits the modelled cospectra reasonably well. This is expected, since low wavenumber components of the spanwise spectra are naturally dampened by the presence of side-walls in the channel and are therefore less susceptible to artifacts from smoothing and averaging. The symmetry in the spanwise direction also forces some degree of homogeneity on the flow, while the streamwise direction is highly inhomogeneous in the buffer layer due to the presence of streaks.³⁴ 

In summary, a cospectral budget model is shown to be capable of reproducing the full law-of-the-wall MVP with no tunable parameters, when forced with the right vertical-velocity variance spectra F_vv+ and TKE dissipation rates. Modeling the dissipation rates by assuming a perfect balance with TKE production (and, hence, varying with $∂ U + ∂ y +$⁠) introduces only small biases to the predicted MVP. Parameterizations used for F_vv+ used in previous studies overestimate the true spectra, at least at the relatively low Reynolds number considered in this study. A major departure between the idealized F_vv+ previously employed and those reported from DNS is the apparent invariance with distance from the wall, in the viscous- and lower buffer-layer, of the cross-over scale 1/k_a in F_vv+. This invariance may be due to the presence of streaks within the buffer region whose dimensions do not scale with distance from the wall.³¹ The modeled cospectra match the DNS cospectra well in their inertial ranges in regions of the flow where the turbulence is closest to fully developed. Low-wavenumber differences between modelled and DNS-observed cospectra are likely due to numerical artifacts in the DNS-observed cospectra introduced through smoothing and averaging. This study further establishes the “spectral link” between the MVP and F_vv+ across the full boundary-layer depth. Recent studies have highlighted the significant role of very-large scale structures on near-wall turbulent stresses.^35,36 The cospectral budget model opens up the possibility of exploring the link between these very-large scale motions (which can be parsimoniously represented by low-wavenumber peaks in the F_vv+ spectra used as inputs to the model) and the MVP. Furthermore, relating the micro- and macro-states (i.e., the distribution of turbulent fluctuations F_vv+, and the MVP, respectively) through a cospectral budget may prove fruitful in explaining macroscopic properties of other wall-bounded turbulent flows, such as those subject to heating or surface roughness effects, and will be the subject of future studies.

The authors thank Professor J. Jimenez, Professor R. D. Moser, Dr. J.-C. del Alamo, Dr. S. Hoyas, and Dr. P. S. Zandonade for making their DNS datasets publicly available; and two anonymous reviewers for providing critical feedback on earlier versions of this manuscript. K.A.M. is funded by the National Science Foundation’s Graduate Research Fellowship Program. G.G.K. acknowledges support from the US Department of Energy through the Office of Biological and Environmental Research Terrestrial Carbon Processes program (Nos. DE-SC0006967 and DE-SC0011461) and No. NSF-EAR-1344703. P.G. is funded by the NASA New Investigator Program NNX14AI36G, U.S. Department of Energy Early Career Award No. DE-SC0014203, and the National Science Foundation CAREER award.