One-point statistics for turbulent pipe flow up to $Re_{\tau} \approx 6000$

We study turbulent flows in a smooth straight pipe of circular cross--section up to $Re_{\tau} \approx 6000$ using direct--numerical-simulation (DNS) of the Navier--Stokes equations. The DNS results highlight systematic deviations from Prandtl friction law, amounting to about $2\%$, which would extrapolate to about $4\%$ at extreme Reynolds numbers. Data fitting of the DNS friction coefficient yields an estimated von K\'arm\'an constant $k \approx 0.387$, which nicely fits the mean velocity profile, and which supports universality of canonical wall-bounded flows. The same constant also applies to the pipe centerline velocity, thus providing support for the claim that the asymptotic state of pipe flow at extreme Reynolds numbers should be plug flow. At the Reynolds numbers under scrutiny, no evidence for saturation of the logarithmic growth of the inner peak of the axial velocity variance is found. Although no outer peak of the velocity variance directly emerges in our DNS, we provide strong evidence that it should appear at $Re_{\tau} \gtrsim 10^4$, as a result of turbulence production exceeding dissipation over a large part of the outer wall layer, thus invalidating the classical equilibrium hypothesis.


Introduction
Turbulent flow in circular pipes has always attracted the interest of scientists, owing to its prominent importance in the engineering practice and because of the beautiful simplicity of the setup. In this respect, the pioneering flow visualizations of Reynolds (1883) may be regarded as a milestone for the understanding of turbulent and transitional flows. The most extensive experimental measurements of high-Reynolds-number pipe flows have been carried out in modern times in the Princeton Superpipe pressurized facility (Zagarola & Smits 1998;McKeon et al. 2005;Hultmark et al. 2010). Those investigations have allowed scientists to measure the main flow features as friction and mean velocity profiles with high precision, and they currently constitute the most comprehensive database for the study of pipe turbulence. However, even the use of specialized micro-fabricated hot-wire probes could not provide fully reliable information about the viscous and buffer layers at high Reynolds numbers (Hultmark et al. 2012). Additional experimental studies of pipe turbulence have been carried out in the high-Reynolds-number actual flow facility (Hi-Reff), a water tunnel with relatively large diameter, which allows for accurate estimation of friction (Furuichi et al. 2015(Furuichi et al. , 2018. Recently, the CICLoPE facility of the University of Bologna (Fiorini 2017;Willert et al. 2017) has been set up, whose large diameter (about 1m) offers a wellestablished turbulent flow with relatively large viscous scales, thus granting higher spatial resolution. Flows in different facilities seem to have sensibly different properties in terms of friction and mean velocity profiles, which we will comment on.
Numerical simulation of pipe turbulence flow has received less interest than other canonical flows, the plane channel in particular, because of additional difficulties involved with discrete solution of the Navier-Stokes equations in cylindrical coordinates, with special reference to treatment of the geometrical singularity at the pipe axis. Early numerical simulations of turbulent pipe flow were carried out by Eggels et al. (1994), at friction Reynolds number Re = 180 (Re = / , with = ( / ) 1/2 the friction velocity, the pipe radius, and the fluid kinematic viscosity). Effects of drag reduction associated with pipe rotation were later studied by Orlandi & Fatica (1997). Higher Reynolds numbers (up to Re ≈ 1140) were reached by Wu & Moin (2008), which first allowed to observe a near logarithmic layer in the mean velocity profile. Flow visualizations and two-point correlation statistics pointed to the existence of high-speed wavy structures in the pipe core region which are elongated in the axial direction, and whose streamwise and azimuthal dimensions do not change substantially with the Reynolds number, when normalized in outer units. Further follow-up DNS studies have been carried out by El Khoury et al. (2013); Chin et al. (2014); Ahn et al. (2013). At present, the highest Reynolds number in pipe flow (Re ≈ 3000) has been reached in the study of Ahn et al. (2015). Although no sizeable logarithmic layer is present yet at those conditions, some effects associated with significant scale separation between innerand outer-scale turbulence were observed, as the presence of a −1 ( being the wavenumber in any wall-parallel direction) power-law ranges in the velocity spectra.
Despite inherent limitations in the Reynolds numbers which can be attained, DNS has the advantage over experiments of yielding immediate access to the near-wall region, and of allowing scientists to measure some flow properties, e.g. the turbulence dissipation rate, which can hardly be measured in experiments. Hence, it is generally claimed that DNS data at increasing Reynolds numbers are needed to prove or disprove theoretical claims related to departure (or not) of the statistical properties of wall-bounded turbulence from the universal wall scaling (Cantwell 2019; Monkewitz 2021; Chen & Sreenivasan 2021). In this paper we thus present DNS data of turbulent flow in a smooth circular pipe at Re ≈ 6000, which is two times higher that the previous state of art. Relying on the DNS data, we revisit current theoretical inferences and discuss implications about possible trends in the extreme Reynolds number regime.

The numerical dataset
The code used for DNS is the spin-off of an existing solver previously used to study Rayleigh-Bénard convection in cylindrical containers at extreme Rayleigh numbers (Stevens et al. 2013). That code is in turn the evolution of the solver originally developed by Verzicco & Orlandi (1996), and used for DNS of pipe flow by Orlandi & Fatica (1997). A second-order finite-difference discretization of the incompressible Navier-Stokes equations in cylindrical coordinates is used, based on the classical marker-and-cell method (Harlow & Welch 1965), with staggered arrangement of the flow variables to remove odd-even decoupling phenomena and guarantee discrete conservation of the total kinetic energy in the inviscid flow limit. Uniform volumetric forcing is applied to the axial momentum equation to maintain constant , and are the number of grid points in the azimuthal, radial and axial directions, respectively, Re = 2 / is the bulk Reynolds number, = 8 /( 2 ) is the friction factor, Re = / is the friction Reynolds number, is the time interval used to collect the flow statistics, and = / is the eddy turnover time. mass flow rate in time. The Poisson equation resulting from enforcement of the divergencefree condition is efficiently solved by double trigonometric expansion in the periodic axial and azimuthal directions, and inversion of tridiagonal matrices in the radial direction (Kim & Moin 1985). An extensive series of previous studies about wall-bounded flows from this group proved that second-order finite-difference discretization yields in practical cases of wall-bounded turbulence results which are by no means inferior in quality to those of pseudospectral methods (e.g. Pirozzoli et al. 2016;Moin & Verzicco 2016). A crucial issue is the proper treatment of the polar singularity at the pipe axis. A detailed description of the subject is reported in Verzicco & Orlandi (1996), but basically, the radial velocity in the governing equations is replaced by = ( is the radial space coordinate), which by construction vanishes at the axis. The governing equations are advanced in time by means of a hybrid third-order low-storage Runge-Kutta algorithm, whereby the diffusive terms are handled implicitly, and convective terms in the axial and radial direction explicitly. An important issue in this respect is the convective time step limitation in the azimuthal direction, due to intrinsic shrinking of the cells size toward the pipe axis. To alleviate this limitation we rely on implicit treatment of the convective terms in the azimuthal direction (Akselvoll & Moin 1996;Wu & Moin 2008), which enables marching in time with similar time step as in planar domains flow in practical computations. In order to minimize numerical errors associated with implicit time stepping, in the present code explicit and explicit discretizations of the azimuthal convective terms are linearly blended with the radial coordinate, in such a way that near the pipe wall the treatment is fully explicit, and near the pipe axis it is fully implicit. The code was adapted to run on clusters of graphic accelerators (GPUs), using a combination of CUDA Fortran and OpenACC directives, and relying on the CUFFT libraries for efficient execution of FFTs (Ruetsch & Fatica 2014). The DNS were carried out on the Marconi-100 machine based at CINECA (Italy), relying on NVIDIA Volta V100 graphic cards. Specifically, 1024 GPUs were used for DNS-F.
Numerical simulations are carried out with periodic boundary conditions in the axial ( ) and azimuthal ( ) directions. The velocity field is then controlled by two parameters, namely the bulk Reynolds number (Re = 2 / , with the pipe radius, the fluid bulk 0.02136 ± 0.13% 24.07 ± 0.18% 7.995 ± 0.29% 14.66 ± 0.073% 0.1697 ± 0.37% DNS-C-SH 0.02164 ± 0.14% 24.09 ± 0.20% 8.071 ± 0.44% 14.37 ± 0.11% 0.1952 ± 0.54% DNS-C-LO 0.02128 ± 0.16% 24.17 ± 0.11% 7.965 ± 0.29% 14.62 ± 0.058% 0.1704 ± 0.40% DNS-C-FT 0.02114 ± 0.12% 24.28 ± 0.14% 7.948 ± 0.27% 14.66 ± 0.078% 0.1691 ± 0.34% DNS-C-FR 0.02132 ± 0.25% 24.10 ± 0.12% 7.886 ± 0.31% 14.41 ± 0.096% 0.1741 ± 0.60% DNS-C-FZ 0.02132 ± 0.21% 24.07 ± 0.26% 8.168 ± 0.38% 14.89 ± 0.14% 0.1727 ± 0.44% DNS-D 0.01839 ± 0.25% 25.56 ± 0.34% 8.397 ± 0.43% 14.79 ± 0.098% 0.1822 ± 0.57% DNS-E 0.01658 ± 0.26% 26.47 ± 0.27% 8.681 ± 0.69% 14.87 ± 0.13% 0.1903 ± 0.93% DNS-F 0.01428 ± 0.36% 28.05 ± 0.35% 9.108 ± 0.72% 15.14 ± 0.20% 0.1993 ± 1.10% Table 2: Uncertainty estimation study: mean values of representative quantities and standard deviation of their estimates. is the friction factor, + is the mean pipe centerline velocity, < 2 > + is the peak axial velocity variance and + is its distance from the wall, and + 11 is the dissipation rate of < 2 > at the wall.  (2015) DNS (channel) 180-5200 ⚪ Table 3: List of other references for data used in the paper velocity, and its kinematic viscosity), and the relative pipe length, / . A list of the main simulations that we have carried out is given in table 1. The mesh resolution is designed based on well-established criteria in the wall turbulence community. In particular, the collocation points are distributed in the wall-normal direction so that about thirty points are placed within + 40 ( = − is the wall distance, and the + superscript is used to denote normalization with respect to and ), with the first grid point at + ≈ 0.05. The mesh is progressively stretched in the outer wall layer in such a way that the mesh spacing is proportional to the local Kolmogorov length scale, which there varies as + ≈ 0.8 +1/4 (Jiménez 2018), and the radial spacing at the pipe axis is Δ + ≈ 8.8. Additional details are provided in a specifically focused publication (Pirozzoli & Orlandi 2021). Regarding the axial and azimuthal directions, finite-difference simulations of wall-bounded flows yield grid-independent results as long as Δ + ≈ 10, + Δ ≈ 4.5 (Pirozzoli et al. 2016), hence the associated number of grid points scales as ≈ / × Re /10, ∼ 2 × Re /4.5. All DNS have been carried out at CFL number close to unity, based on the radial convective time step limitation. The CFL number along the axial direction is typically smaller by a factor two. The time step expressed in wall units ( / 2 ) ranges from Δ + = 0.55 in DNS-A to Δ + = 0.15 in DNS-F. According The sampling errors for some key properties discussed in this paper have been estimated using the method of Russo & Luchini (2017), based on extension of the classical batch means approach. The results of the uncertainty estimation analysis are listed in table 2, where we provide expected values and associated standard deviation for the friction factor ( ), mean centerline velocity ( ), peak axial velocity variance and its position ( 2 and , respectively), and dissipation rate of streamwise velocity variance ( 11 ). Here and elsewhere, capital letters are used to denote flow properties averaged in the homogeneous spatial directions and in time, brackets denote the averaging operator, and lower-case letters to denote fluctuations from the mean. We find that the sampling error is generally quite limited, being larger in the largest DNS, which have been run for shorter time. In particular, in DNS-F the expected sampling error in friction, centerline velocity and peak velocity variance is about 0.5%, whereas it is about 1% for the wall dissipation. Additional tests aimed at establishing the effect of axial domain length and grid size have been carried out for the DNS-C flow case, whose results are also reported in table 2. We find that doubling the pipe length yields a change in the basic flow properties of about 0.2 − 0.3%, whereas halving it yields changes of about 1% in friction and peak velocity variance, and up to 10% in the wall dissipation. Hence, consistent with previous studies (Chin et al. 2010), we believe that the selected pipe length ( / = 15) is representative of an infinitely long pipe, at least for the purposes of the present study. In order to quantify uncertainties associated with numerical discretization, additional simulations have been carried out by doubling the grid points in the azimuthal, radial and axial directions, respectively. Based on the data reported in the table, after discarding the short pipe case, we can thus quantify the uncertainty due to numerical discretization and limited pipe length to be about 0.3% for the friction coefficient and pipe centerline velocity, 0.6% for the peak velocity variance, and 0.9% for the wall dissipation.

Results
Qualitative information about the structure of the flow field is provided by instantaneous perspective views of the axial velocity field, provided in figure 1. Although finer-scale details are visible at the higher Re, the flow in the cross-stream planes is always characterized by a limited number of bulges distributed along the azimuthal direction, which closely recall the POD modes identified by Hellström & Smits (2014), and which correspond to alternating intrusions of high-speed fluid from the pipe core and ejections of low-speed fluid from the wall. Streaks are visible in the near-wall cylindrical shells, whose organization has clear association with the cross-stream pattern. Specifically, regardless of the Reynolds number, -sized low-streaks are observed in association with large-scale ejections, whereas -sized high-speed streaks occur in the presence of large-scale inrush from the core flow. At the same time, smaller streaks scaling in wall units appear, corresponding to buffer-layer ejections/sweeps. Hence, organization of the flow on at least two length scales is apparent here, whose separation increases with Re . Mean friction is obviously a parameter of paramount importance as it is related to power expenditure to sustain the flow. In figure 2, we show the friction factor, namely  A correlation generally used for smooth pipes is the Prandtl friction law, We believe that this difference may be related to different grid resolution in the azimuthal direction, which was + Δ = 7 − 8 in those previous studies, and 4 − 5 in our DNS. Our data in fact show minimal overshoot at low Reynolds number, and consistent undershoot from Prandtl law by about 2%. Regarding experiments, Superpipe data typically tend to lie above the theoretical curve by about 2%, whereas the CICLoPE and Hi-Reff data tend to fall short of it. Although the range of data overlap is not extensive, it appears that DNS data tend to be more consistent with the CICLoPE and Hi-Reff data than with other datasets. Fitting the current DNS data with a functional relationship as (3.2), yields ≈ 2.102, ≈ 1.148, with an inferred value of the von Kármán constant of = 0.387 ± 0.004, with uncertainty estimates based on 95% confidence bounds from the curve-fitting procedure. This value is extremely close to that suggested by Furuichi et al. (2018), who reported = 0.386 as an average value over a very wide range of Reynolds numbers, and also very close to values reported in boundary layers (Nagib & Chauhan 2009) and channels (Lee & Moser 2015). If this trend is extrapolated, deviations of about 4% from the standard Prandtl law would result at Re = 10 7 .
The mean velocity profile in turbulent pipes has received extensive attention from theoretical studies, much of the early debate being focused on whether a log law or a power law better fits the experimental results (Barenblatt et al. 1997), mainly carried out in Overall, good agreement is observed across various sources as far as the inner and the overlap regions are concerned, with data gradually approaching a logarithmic distribution, here identified by visual fitting as + = 1/ log + + 4.53, using the value of = 0.387 determined from friction data. This is quite close to estimates based on direct fitting of the mean velocity profile in pipe flow (Marusic et al. 2013), which yielded + = 1/0.391 log + + 4.34. The DNS velocity profiles for Re 10 3 follow this distribution with deviations of no more than 0.1 wall units from + ≈ 30 to / ≈ 0.15, whence the core region develops. Differences with respect to previous DNSs are concentrated in the core region, which seemingly stronger wake in some datasets, including our own, Wu & Moin (2008) and Ahn et al. (2013), and weaker in others (El Khoury et al. 2013;Chin et al. 2014), reflecting previously noted differences in the friction coefficient. Especially satisfactory is the excellent agreement between our DNS-E velocity profile and the data of Ahn et al. (2015) at Re ≈ 3000. Comparison of our DNS dataset with experimental data also shows overall good agreement, although some differences are quite clear in the core region, in which Superpipe experiments consistently yield lower + , which translates into lower friction.
More refined information on the behaviour of the mean velocity profile can be gained from inspection of the log-law diagnostic function which is shown in figure 4, and whose constancy would imply the presence of a genuine logarithmic layer in the mean velocity profile. The figure supports universality of the innerscaled axial velocity for Re 10 3 , up to + ≈ 100, where Ξ attains a minimum, and the presence of an outer maximum at / ≈ 0.6. Between these two sites the distribution is roughly linear, as can be better appreciated in panel (b), with nearly constant slope when expressed in outer coordinates. Approximate linear variation of the diagnostic function in channel flow was observed by Jiménez & Moser (2007), who, based on refined overlap arguments expressed by Afzal & Yajnik (1973), proposed the following fit where , are adjustable constants, and is the von Kármán constant. Here we find that the set of constants = 0.387, = 2.0, = 0, yields overall good approximation of the pipe DNS data. The consequence is that a genuine logarithmic layer would only be attained at infinite Reynolds number. In this respect, Superpipe data seem to suggest the formation of a plateau at Re 10 4 , although the scatter of points is quite significant. Hence, DNS at higher Reynolds number would be most welcome to confirm or refute our findings, and possibly determine more accurate values of the extended log-law constants in (3.4).
Comparison with Superpipe data is presented in outer units in figure 5, limited to the higher Re cases. Despite differences in the Reynolds number, the velocity profiles now agree very well, throughout the outer layer. This observation would suggest problems with correct estimation of the friction velocity, which however seems unlikely both in DNS, in which we independently evaluate friction velocity by computing the wall derivative of the velocity profile and from momentum balance, and in experiments, as measurements of the pressure drop are regarded to have low uncertainty. Hence, reasons for this discrepancy are not known, and additional experiments as those currently carried out in the large CICLoPE facility would be especially useful and welcome. Unfortunately, velocity profiles along the full radial span are not available at the moment for that facility.
The structure of the core region is examined in detail in figure 6, where the mean velocity profiles are shown in defect form. Although full outer-layer similarity is not reached at the conditions of our DNS study (also see the inset of figure 3(a)), scatter across the Reynolds number range and with respect to Superpipe profiles for / 0.2 is no larger than 5%. As suggested by Pirozzoli (2014), the core velocity profiles can be closely approximated with a simple quadratic function, reflecting near constancy of the eddy viscosity. In particular, we find that the formula  suggest logarithmic increase with Re according to where we find = = 0.387 as for the friction law, and = 5.85. For convenience, the trend of / is also presented, having in fact the same logarithmic growth with Re . With some previously noted differences, all pipe flow DNSs seem to exhibit a consistent trend in the accessible range. While the trend is very similar at low Reynolds number, experimental data yield consistently lower values of + , especially those from the Superpipe. At Reynolds numbers higher than about Re = 10 4 , experiments seem to suggest milder growth rate, although significant differences emerge between the Superpipe and the Hi-Reff datasets. Hence, whether this is the result of a change of behaviour at high Reynolds number, or some form of shortcoming of experiments is difficult to say. As a result of the observed identity (or very close vicinity) of the von Kármán constant for the centerline and for the bulk velocity, figure 7(b) highlights that their ratio approaches unity at large Re, supporting the inference that pipe flow asymptotes to plug flow in the infinite-Reynolds-number limit (Pullin et al. 2013). Regarding that study, it is worthwile noticing that one of the assumptions made in the analysis is that the wall-normal location of the onset of the logarithmic region is either finite, or increases no faster than Re . Interpreting the near-wall minimum of the diagnostic function in figure 4 as the root of the (near) logarithmic layer, our data well support that assumption. Whereas the curvature of the core velocity profile is not changing substantially when expressed in wall units (see figure 6), it would become vanishingly small when expressed in outer units. However, as figure 7(b) suggests, this trend is extremely slow. Interestingly, again despite some scatter, DNS and experiments here seem to indicate a common trend with overall monotonic decrease, perhaps with a 'bump' in the range of Reynolds numbers in the low thousands. The DNS data points at the highest Reynolds numbers (DNS-D,E,F) now appear to be in good agreement with Superpipe experiments, which is in line with the previously noted agreement of the outer-scaled mean velocity profiles.
The shows overall good agreement of all turbulence intensities. Comparison with Superpipe data at Re = 3000 is also very good, with exception of the near-wall peak which is likely to be over-estimated in experiments. DNS-F data seem to be well bracketed by Superpipe and CICLoPE measurements at nearby Reynolds numbers, and also compare very well with experimental data for plane channel flow (Schultz & Flack 2013).
Distributions of the turbulent shear stress are shown in figure 9. As is well established (e.g. Lee & Moser 2015), the shear stress profiles tend to become flatter at higher Re , the peak value rises towards unity, and its position moves farther from the wall, in inner units. In particular, exploiting mean momentum balance and assuming the presence of a logarithmic layer in the mean axial velocity, the following prediction follows for the position of the turbulent shear stress peak (Afzal 1982) which is intermediate between inner and outer scaling. This observation has led some authors to argue about the relevance of a 'mesolayer' (Long & Chen 1981;Wei et al. 2005). The asymptotic relationship (3.8) is satisfied with good accuracy starting at Re ≈ 10 3 , reflecting the onset of a near logarithmic layer. The behaviour of the Reynolds stresses when expressed as a function of the outer-scaled wall distance, which is shown in figure 10 is also of great theoretical interest. In fact, according to the attached-eddy model (Townsend 1976; Marusic & Monty 2019), the wallparallel velocity variances are expected to decline logarithmically with the wall distance in the outer layer, hence where , are universal constants. Regarding the axial stress, Marusic et al. (2013) argued that Superpipe data at the highest available Reynolds number are best fit with 1 = 1.23, 1 = 1.56, with a sensible logarithmic layer only emerging at Re > 10 4 , in the range of wall distances 3Re 1/2 + 0.15Re . DNS data only show the formation of a near logarithmic layer farther away from the wall, which is not where it is expected from theoretical arguments. Hence, little can be said in this respect. The azimuthal velocity variance, shown in figure 10(b), has a more benign behaviour, and it features clear logarithmic layers even at modest Re . Fitting the DNS data yields 3 = 0.40, 3 = 1.0, which is very close to what found in channels (Bernardini et al. 2014;Lee & Moser 2015). Measurements of pipe flow carried out in the CICLoPE facility (Örlü et al. 2017) yielded 3 = 0.63, 3 = 1.21, hence much larger values than in DNS. Possible overestimation of the wall-normal and azimuthal Reynolds stresses was in fact acknowledged by the authors of that paper.
Quantitative insight into Reynolds number effects is provided by inspection of the amplitude of the inner peak of the axial velocity variance, which we show in figure 11. The general theoretical expectation is that the peak grows logarithmically with Re owing to the increasing influence of distant, inactive eddies (Marusic & Monty 2019). However, some recent experimental results (Willert et al. 2017), and theoretical arguments (Chen & Sreenivasan 2021) suggest that such growth should eventually saturate. Although difference between slow logarithmic growth and attainment of an asymptotic value is quite subtle in practice, the theoretical interest is high as in the latter case universality of wall scaling would be eventually restored. Within the investigated range of Reynolds numbers, our DNS data in fact support continuing logarithmic increase. Comparison with channel data (Lee & Moser 2015) shows some difference, which might result from stronger geometrical confinement of distant eddies in the pipe geometry. However, differences tend to becomes smaller at higher Re . In quantitative terms, we find the slope of logarithmic increase to be about 0.67, a bit steeper than found in channel flow DNS ( While our DNS data cannot be used to directly evaluate the theoretical predictions owing to limited achievable Reynolds number, they can be used to better scrutinize the foundations of the theoretical arguments. The main argument made by Chen & Sreenivasan (2021), although not thoroughly justified in our opinion, was that since turbulence kinetic energy production is bounded, the wall dissipation must also stay bounded. Hence, let = − d /d be the turbulence kinetic energy production rate, and 11 = |∇ | 2 be the dissipation rate of the axial velocity variance, those authors first argue that the wall limiting value of 11 should scale as 11 + = 1/4 − /Re 1/4 , (3.10) with a suitable constant. In figure 11 we explore deviations of and of the peak turbulence kinetic energy production, say , from their asymptotic value, namely 1/4. According to analytical constraints (Pope 2000), we find that production tends to its asymptotic value quite rapidly, as about 1/Re . Consistent with equation (3.10), the wall dissipation also tends to 1/4, more or less at the predicted rate, thus empirically validating the first assumption. The next argument advocated by Chen & Sreenivasan (2021) is that wall balance between viscous diffusion and dissipation and Taylor series expansion of the axial velocity variance near the wall yields 2 + ∼ + 11 +2 , (3.11) whence, from assumed invariance of the peak location of 2 (say, + ), saturation of growth of the peak velocity variance would follow. Table 2 suggests that this second assumption is in fact violated, as the position of the peak slightly increases with Re , with non-negligible effect on the peak variance as it appears in squared form in equation (3.11). As a consequence, logarithmic growth of the peak velocity variance still holds, at least in the range of Reynolds numbers currently accessible to DNS. Although no distinct outer peak of the axial velocity variance is found at the Reynolds numbers under scrutiny, it is nevertheless instructive to explore the scaling of the velocity fluctuations in the range of wall distances where the peak is expected to occur. For that purpose, we consider the outer position where the second logarithmic derivative of the velocity variance vanishes, which in the present DNS ranges from + ≈ 115 for DNS-A, to + ≈ 140 for DNS-F. Weak dependence of the inner-scaled outer peak position on Re , although at much higher Reynolds number, was also noticed by Hultmark et al. (2012). The resulting distribution is shown in figure 12. All DNS data fall nicely on a logarithmic fit, and they seem to connect smoothly to the experimental results, whose scatter and uncertainty is expected to be much less than for the inner peak. Experiments indicate a change of behaviour to a shallower logarithmic dependence with slope of about 0.63 (Fiorini 2017;Pullin et al. 2013), which would be very close to the growth rate of the inner peak (see figure 11). The figure suggests that verification of this effect would require Re of about 10 4 .
As pointed out by Hultmark et al. (2012), the formation and growth of an outer peak of the axial velocity variance has important theoretical and practical implications. From the modeling standpoint, no current RANS model is capable of predicting non-monotonic behaviour of Reynolds stresses outside the buffer layer. From the fundamental physics standpoint, the presence of an outer peak is suggestive of violation of equilibrium between turbulence production and dissipation, thus invalidating one of the possible arguments advocated for the formation of a logarithmic mean velocity profile (Townsend 1976). DNS allows to substantiate this scenario, and for that purpose in figure 13, we show the relative excess of turbulent kinetic energy production ( ) over its total dissipation rate, here defined as = ∇ 2 , which lumps together dissipation rate and viscous diffusion. Data confirm the presence of a near-universal region confined to the buffer layer (say, 8 + 35), in which production exceeds dissipation by up to 40%. Data also show the onset, starting from DNS-B, of another region farther from the wall with positive unbalance, whose inner limit is constant in inner units, at + = 100, and whose outer limit tends to become constant at high Re in outer units, at / ≈ 0.4. The peak unbalance at high Reynolds number is about 17%, and its position seems to scale more in inner than in outer units. Turbulence kinetic energy production excess in the presence of a (near) logarithmic mean velocity profile can be interpreted by recalling that only part of the kinetic energy which is being generated is converted to active motions which carry turbulent shear stress, and the rest is used to feed inactive motions. This finding clearly indicates that at high enough Reynolds number the outer wall layer becomes a dynamically active part of the flow, having the potential to transfer energy both to the core flow, and towards the wall, in the form of imprinting on the near-wall layer (Marusic & Monty 2019).

Concluding comments
Although DNS of wall turbulence is still confined to a moderate range of Reynolds numbers, it is beginning to approach a state in which some typical phenomena of the asymptotically high-Re emerge. Given its ability to resolve the near-wall layer, DNS lends itself to testing theories of wall turbulence and to in-depth scrutiny of experimental data. In this work, DNS of flow in a smooth pipe has been carried out up to Re ≈ 6000, which, although still far from what achievable in experimental tests, allows to uncover a number of interesting issues, in our opinion. First, we have noted that DNS data fall systematically short of the classical Prandtl friction law, by as much as 2%. This evidence is not consistent with data from the Superpipe facility, although other recent data from the CICLoPE and Hi-Reff facilities seem to yield similar trends. DNS data fitting suggests that a logarithmic law as (3.2) still holds, however with a von Kármán constant ≈ 0.387, which matches extremely well the value quoted by Furuichi et al. (2018), and which would reconcile pipe flow with plane channel and boundary layer flows, thus corroborating claims made by Marusic et al. (2013). A logarithmic profile with ≈ 0.387 well fit the mean axial velocity distributions for 30 + 0.15Re , although linear deviations are clearly visible, as argued by Afzal & Yajnik (1973); Luchini (2017), which is taken into account yield excellent representation of the velocity profiles up to / ≈ 0.5. It is remarkable that the same value of the von Kármán constant also well fits the mean centerline velocity distribution (see figure 7), which is found to grow logarithmically throughout the range of Re under investigation. This finding is quite reasonable as it corroborates that the eventual state of turbulent flow in a pipe should be plug flow, as argued by Pullin et al. (2013), hence → as Re → ∞. This would however seemingly contrast recent measurements made in the CICLoPE facility (Nagib et al. 2017), which rather suggest a different von Kármán constant for the bulk and the centerline velocity. Experimental data at Re 10 4 in fact suggest deviations of + from the logarithmic trend found DNS, however this effect requires further confirmation, as data are quite scattered. The core velocity profile is found to be to a good approximation parabolic, with curvature which is nearly constant in wall units, and decreasing in outer units.
Regarding the velocity fluctuations, we find evidence for continuing logarithmic increase of the inner-peak magnitude with Re . Some experiments and theoretical arguments would indicate that beyond Re ≈ 10 4 a change of behaviour might occur, which however is very difficult to quantify. DNS is probably of little use in this respect, as in order to clearly discern among the various trends, Re in excess of 10 5 are likely to be needed. As predicted by the attached-eddy hypothesis, the wall-parallel velocity variances in the outer layer tend to form logarithmic layers, which are especially evident in the azimuthal velocity. Although we do not find direct evidence for the existence of an outer peak of the axial velocity variance, our results highlight the occurrence of an outer site with substantial turbulence production excess over dissipation, thus contradicting the classical equilibrium hypothesis and likely to yield a distinct peak at Re ≈ 10 4 , as Superpipe and CICLoPE experiments suggest. Investigating these and other violations of universality of wall turbulence to extrapolate asymptotic behaviours is a formidable challenge for theoreticians in years to come. D