Direct numerical simulations of turbulent pipe flow up to $Re_\tau \approx 5200$

Abstract Well-resolved direct numerical simulations (DNS) have been performed of the flow in a smooth circular pipe of radius $R$ and axial length $10{\rm \pi} R$ at friction Reynolds numbers up to $Re_\tau =5200$ using the pseudo-spectral code OPENPIPEFLOW. Various turbulence statistics are documented and compared with other DNS and experimental data in pipes as well as channels. Small but distinct differences between various datasets are identified. The friction factor $\lambda$ overshoots by $2\,\%$ and undershoots by $0.6\,\%$ the Prandtl friction law at low and high $Re$ ranges, respectively. In addition, $\lambda$ in our results is slightly higher than in Pirozzoli et al. (J. Fluid Mech., vol. 926, 2021, A28), but matches well the experiments in Furuichi et al. (Phys. Fluids, vol. 27, issue 9, 2015, 095108). The log-law indicator function, which is nearly indistinguishable between pipe and channel up to $y^+=250$, has not yet developed a plateau farther away from the wall in the pipes even for the $Re_\tau =5200$ cases. The wall shear stress fluctuations and the inner peak of the axial turbulence intensity – which grow monotonically with $Re_\tau$ – are lower in the pipe than in the channel, but the difference decreases with increasing $Re_\tau$. While the wall value is slightly lower in the channel than in the pipe at the same $Re_\tau$, the inner peak of the pressure fluctuation shows negligible differences between them. The Reynolds number scaling of all these quantities agrees with both the logarithmic and defect-power laws if the coefficients are properly chosen. The one-dimensional spectrum of the axial velocity fluctuation exhibits a $k^{-1}$ dependence at an intermediate distance from the wall – also seen in the channel. In summary, these high-fidelity data enable us to provide better insights into the flow physics in the pipes as well as the similarity/difference among different types of wall turbulence.


Introduction
Turbulent flows that are constrained by a wall (referred to as 'wall turbulence') are common in nature and technology. Roughly half of the energy spent in transporting fluids through pipes or vehicles through air or water is dissipated by the turbulence near the walls (Jiménez 2012). Therefore, an improved understanding of the underlying physics of these flows is essential for modelling and control (Kim 2011;Canton et al. 2016;Yao, Chen & Hussain 2018). The spatially evolving boundary layer, the (plane) channel and the pipe are three canonical geometrical configurations of wall turbulence. Different from boundary layer and channel flows, azimuthal periodicity is inherent to pipe flows. Therefore, pipe flow is the most canonical case, being completely described by the Reynolds number (Re) and the axial length -the effect of the latter is limited if sufficiently large (El Khoury et al. 2013;Feldmann, Bauer & Wagner 2018).
Sustained interest in high-Re wall turbulence stems from numerous open questions regarding the scaling of turbulent statistics, as reviewed in Marusic et al. (2010b) and Smits, McKeon & Marusic (2011b). For example, a characteristic of high-Re wall turbulence is the logarithmic law in the mean velocity with an important parameter, namely the von Kármán constant κ, whose value and universality among different flow geometries are still highly debated (Nagib & Chauhan 2008;She, Chen & Hussain 2017). Also, there is no consensus on whether the near-wall peak of the streamwise velocity fluctuations increases continuously with Re (Marusic, Baars & Hutchins 2017) or eventually saturates at high Re (Chen & Sreenivasan 2021;Klewicki 2022). Furthermore, the existence of an outer peak in the streamwise velocity fluctuations, as indicated by experiments, is also highly debated (Hultmark et al. 2012;Willert et al. 2017). Other questions, which can be answered only with substantially higher Re values, concern the scaling and generation mechanism of the various flow structures. Large-scale motions (LSMs) and very-large-scale motions (VLSMs), with lengths 5R up to 20R, have been found experimentally in the outer region of pipe flows (Kim & Adrian 1999;Monty et al. 2007). Here, R is the radius of the pipe. Due to the increasing strength of these structures with Re, they have a footprint quite close to the wall (Monty et al. 2007), in the form of amplitude modulation as reviewed by e.g. Dogan et al. (2018).
Fundamental studies of wall turbulence require accurate representations or measurements of the flows, which were typically carried out via experiments (Zagarola & Smits 1998;McKeon, Zagarola & Smits 2005;Smits et al. 2011a;Furuichi et al. 2015Furuichi et al. , 2018Talamelli et al. 2009;Fiorini 2017). However, decades of experimental research have shown that obtaining unambiguous high-Re data, particularly near the wall, remains a challenge. This is because the smallest scales decrease with increasing Re -leading to large uncertainties in determining the probe locations and turbulence intensities. Advances in computer technology (in both hardware and software) have enshrined direct numerical simulations (DNS) as an essential tool for turbulence research. Although only moderate Re can be achieved at the current stage, DNS provide extensive, detailed data compared to experiments -even close to the walls where experimental data are very difficult to obtain. One of the earliest DNS for wall turbulence were performed by Kim, Moin & Moser (1987) for the channel flow at friction Reynolds number Re τ (≡ u τ h/ν) ≈ 180 (here, u τ is the friction velocity, h is the half-channel height, and ν is the fluid kinematic viscosity). They found good agreement between DNS and experimental data of Hussain & Reynolds (1975), except in the near-wall region. The discrepancy was speculated to be caused by the inherent near-wall hot-wire measurement errors. Eggels et al. (1994) subsequently conducted the first DNS of pipe flow at Re τ ≈ 180 to investigate the differences between channel and pipe flows.
Numerous DNS investigations have been carried out in the aftermath of these pioneering studies, with Re increasing progressively as a result of increased computational power (Moser, Kim & Mansour 1999;Wu & Moin 2008;Bernardini, Pirozzoli & Orlandi 2014). However, among them, only those with Re τ ≥ 10 3 are of particular engineering interest as this is the range of Re relevant to industrial applications. Also, it is in this range that the high-Re characteristics of wall turbulence start to manifest. One of the highest Re τ large-domain DNS was performed by Lee & Moser (2015) for channel flows at Re τ = 5186 with domain size L x × L z = 8πh × 3πh. Compared to numerous DNS for channel flows (Bernardini et al. 2014;Lozano-Durán & Jiménez 2014;Lee & Moser 2015), fewer high-Re studies have addressed pipe flow, and most of them are limited to Re τ ≈ 1000. For example,  performed DNS at Re τ ≈ 1000 with length 30R and established the existence of VLSMs of scale up to O(20R). El Khoury et al. (2013) used a spectral-element method to perform DNS for Re τ up to 1000 with length L z = 25R. Chin, Monty & Ooi (2014) found that the mean velocity profile does not exhibit a strictly logarithmic layer with Re τ up to 2000, necessitating a finite-Re correction like those introduced by Afzal (1976) and Jiménez & Moser (2007). To quantify the effects of computational length and Re, Feldmann et al. (2018) conducted DNS for 90 ≤ Re τ ≤ 1500 with L z up to 42R. They confirmed that L z = 42R is sufficiently large to capture the LSMand VLSM-relevant scales. Ahn et al. (2015) performed DNS of pipe flow at Re τ ≈ 3000 for length 30R. They claimed that the mean velocity follows a power law in the overlap region, and observed a clear scale separation between inner-and outer-scale turbulence. So far, the largest DNS of pipe flow have been done at Re τ ≈ 6000 with a relatively short length (L z = 15R) by Pirozzoli et al. (2021) based on a lower-order numerical method.
In general, one expects various simulations and experiments to agree with each other to a high degree. However, a comparison among several datasets in spatially developing turbulent boundary layers (Schlatter & Örlü 2010), channels (Lee & Moser 2015) and pipes  shows, surprisingly, considerable variations among the various DNS, even for basic measures such as the shape factor, the friction coefficient and the von Kármán constant. Accurate turbulence statistics are very much needed, both for understanding turbulence physics and for developing, adapting and validating turbulence models. Here, we present a new high-fidelity DNS dataset of turbulent pipe flow generated with a pseudo-spectral method for Re τ up to 5200 and with axial length L z /R = 10π, which is long enough to capture the LSMs and VLSMs reported in experimental studies (Guala, Hommema & Adrian 2006). The accuracy of this dataset is quantified by using the newly developed uncertainty quantification method. In addition, the dataset is compared extensively with other DNS and experimental data for turbulent pipe and channel flows.

Simulation details
DNS of incompressible turbulent pipe flows are performed using the pseudo-spectral code 'OPENPIPEFLOW' developed by Willis (2017). The radial, axial and azimuthal directions are represented by r, z and θ, and the corresponding velocity components are u r , u z and u θ . In the periodic axial (z) and azimuthal (θ) directions, Fourier discretization is employed, and the numbers of Fourier modes in the z-and θ -directions are N z and N θ , respectively. Note that in the physical space, the numbers of grid points in the zand θ -directions increase by a factor 3/2 due to dealiasing. In the radial (r) direction, a Here, Re τ (≡ u τ R/ν) and Re b (≡ 2U b R/ν) are the frictional and bulk Reynolds numbers, respectively; N z and N θ are the number of dealiased Fourier modes in the axial and azimuthal directions, and N r is the number of grid points in the radial direction; z and (Rθ) are the grid spacings in the axial and azimuthal directions, defined in terms of the Fourier modes. In the radial direction, r + w represents the grid spacing at the wall, and r + max denotes the maximum grid spacing. Also, Tu τ /R is the total eddy-turnover time without the initial transient phase.
central finite difference scheme with a nine-point stencil is adopted. With this, the firstand second-order derivatives are, respectively, calculated to 8th and 7th order, which is reduced while approaching the boundaries. The grid points are distributed according to a hyperbolic tangent function so that high wall-normal velocity gradients in the viscous sublayer can be resolved. In addition, the first few points near r = 0 are also clustered to preserve the high order of the finite difference scheme across the pipe axis. No-slip and no-penetration boundary conditions are enforced at the wall (i.e. r/R = 1), where an actual grid point is located. To avoid numerical singularity, there is no grid point at r = 0; and the boundary conditions on the extrapolated fields are specified based on symmetry. In particular, given a Fourier mode with azimuthal index m, each mode is odd/even if m is odd/even for u z and p, and each mode is even/odd if m is odd/even for u r and u θ . The governing equations are integrated with a second-order semi-implicit time-stepping scheme. The flow is driven by a pressure gradient, which varies in time to ensure that the mass flux through the pipe remains constant. For more details about the code and the numerical methods, see Willis (2017).
Five different Reynolds numbers, Re τ ≈ 180, 550, 1000, 2000, and 5200, are considered. The detailed simulation parameters, such as domain sizes and grid sizes, are listed in table 1. The simulations are performed with resolutions comparable to those used in the prior simulations, e.g. Lozano-Durán & Jiménez (2014) and Lee & Moser (2015). In particular, for Re τ ≤ 2000, the axial and azimuthal resolutions employed here satisfy the criterion suggested by Yang et al. (2021) for capturing 99 % of the wall shear stress events. For the highest Re τ case (i.e. ≈ 5200), N z = 12 288 and N θ = 5120 Fourier modes are used in the z-and θ-directions -corresponding to an effective resolution z + = L + x /N z = 12.8 and (Rθ) + = (2πR + )/N θ = 5.1. Hereinafter, the superscript + indicates non-dimensionalization in wall units, i.e. with kinematic viscosity ν and friction velocity u τ . Along the radial direction, the minimum and maximum grid spacings in wall units are 0.2 and 8.6, respectively. Additional information on this highest Re τ simulation is provided in Appendix A.
For comparison, several DNS and experimental data from the literature are included. The details are listed in table 2. To further validate the accuracy of our simulation, an additional simulation at Re τ = 2000 is performed using NEK5000 (hereinafter, this case is denoted as NEK5000 2K). The numerical set-up and mesh generation are the same as those in El Khoury et al. (2013). The length of the pipe is chosen as L z = 35R, and the total number of spectral elements is 7 598 080. With the polynomial order set to 12, the     Re τ = 6000: Pirozzoli et al. (2021) Re τ = 10 000: Hoyas et al. (2022) Re τ = 5200: Lee & Moser (2015) y/R total number of grid points is approximately 13.1 × 10 9 . The grid spacing is comparable to that used in El Khoury et al. (2013) in all directions. The uncertainty in the flow quantities due to the finite time-averaging is estimated using the methods described in Rezaeiravesh et al. (2022) andD. Xavier et al. (personal communication 2022). For the central moments of velocity and pressure, the central limit theorem is applied to the time samples averaged over the z-and θ-directions. The associated time-averaging uncertainty is estimated using an autoregressive-based model for the autocorrelation function; see e.g. Oliver et al. (2014). For estimating the uncertainty in the combination of central moments, the method proposed by Rezaeiravesh et al. (2022) is employed. See Appendix B for further discussion on the method and estimated uncertainties in the first-and second-order velocity moments. For the Re τ ≈ 5200 case, the estimated standard deviation of the mean axial velocity (U + ) is less than 0.1 %, and the estimated standard deviation of the velocity variance (i.e. u 2 r + , u 2 θ + , u 2 z + ) and covariance ( u r u z + ) is less than 1 % in the near-wall region (y + < 100), and approximately 5 % in the core region. Hereinafter, the velocity fluctuations are denoted using the prime symbol (e.g. u r ), and the ensemble (in both time and space) averaged quantities are expressed using a capital letter or bracket (e.g. U or u r u θ ).
In addition, the mean momentum equation is employed to ensure that the simulation is statistically stationary. Due to momentum balance, the total stress, which is the sum of Reynolds shear stress u r u θ + and mean viscous stress ∂U + /∂y + , is linear in a statistically stationary turbulent pipe flow: where y = R − r. Figure 1 shows the residual in (2.1) for the Re τ ≈ 5200 case. The discrepancy between the analytic linear profile (i.e. 1 − y/R) and total stress profile (i.e. ∂U + /∂y − u r u θ + ) from the simulation is less than 0.002 in wall units and is comparable to other high-Re DNS in the literature. Note that this discrepancy is much smaller than the standard deviation of the estimated total stress (see Appendix B).

Flow visualization
The Re effect on the flow structure is illustrated qualitatively in figure 2, showing cross-sectional views of the instantaneous axial velocity u z . Although large scales dominate in the central region of the pipe for all Re τ cases, there is a general increase in the range of scales with increasing Re τ . The average spacing between near-wall low-speed streaks is around (Rθ) + = 100. For the lowest Re τ studied here, approximately ten evenly distributed low-speed structures are seen in figure 2(a), identified by the plume-shaped black regions ejecting from the wall. For our highest Re τ , the streak spacing is reduced to approximately 0.02R, and these fine-scale streaks can hardly be identified from the full cross-section in figure 2(e). A zoomed-in view of the near-wall region with domain size (1000, 200) in wall units in (r, θ) directions is provided to better visualize these structures, which share patterns quite similar to those in low Re τ cases.
Figures 2( f,g) show the inner-scaled instantaneous axial velocity fluctuations u z /u τ on an unrolled cylindrical surface (z − rθ) for Re τ ≈ 5200 at y + = 15 and y/R = 0.5, respectively. At y + = 15, u z /u τ shows the organization of the streaks, whose characteristics are better illustrated in a zoomed-in view with domain size (2000, 400) in wall units. Consistent with previous findings (Hellström, Sinha & Smits 2011;Bernardini et al. 2014), the flow in the outer region exhibits very long positive/negative u z /u τ patches -indicating the presence of LSMs and VLSMs.

Friction factor
The mean friction (or wall shear stress), which is proportional to the pressure drop or the amount of energy required to sustain the flow, is an important parameter and has been studied extensively (Blasius 1913;McKeon et al. 2005;Furuichi et al. 2015;Pirozzoli et al. 2021). A semi-empirical relation between the friction factor λ = 8τ z,w /(ρU 2 b ) and Re is given as (known as the Prandtl friction law) where the constant A is related to the von Kármán constant as A = 1/(2κ √ 2 log 10 (e)). Curve-fitting the experimental data over 3.1 × 10 3 < Re b < 3.2 × 10 6 by Nikuradse (1933) yields A = 2.0 and B = 0.8, which corresponds to κ = 0.407. However, notable deviations were observed when comparing the Prandtl friction law with other experimental data. For example, McKeon et al. (2005) showed that for the Princeton Superpipe data (in the range 3.1 × 10 4 ≤ Re b ≤ 3.5 × 10 7 ), the constants of the Prandtl law work only over a limited range of Re b . New constants (i.e. A = 1.920 and B = 0.475) and additional Re-dependent corrections are needed to better fit the data in the entire Re b range. For the 'Hi-Reff' data, Furuichi et al. (2015) found that λ deviates from the Prandtl law by approximately 2.5 % in the lower Re b region, and −3 % in the high Re b region. In addition, λ, although agreeing with the Superpipe data in the low Re b range, deviates for Re b > 2 × 10 5 .  its effect on one-point statistics appears limited when L z /R ≥ 7 (Feldmann et al. 2018;Pirozzoli et al. 2021). All DNS and experimental data seem to follow λ p . However, the scatter is better highlighted by examining the relative error with respect to the Prandtl law (i.e. λ/λ p − 1). As depicted in figure 3(b), all DNS data overshoot λ p at the low Re b (i.e. ≤ 4 × 10 4 ). Our data, which agree with Wu & Moin (2008), El Khoury et al. (2013) and Chin et al. (2014) and the new simulation at Re = 2000 using NEK5000, exceed λ p by approximately 2 %, while the results of Ahn et al. (2013) and Pirozzoli et al. (2021) are closer to λ p (within 1 % for Re b < 5 × 10 4 ). Pirozzoli et al. (2021) attributed this discrepancy to the different grid resolutions employed in the θ-direction, which is (Rθ) + = 4-5 in theirs and Ahn et al. (2013), but 7-8 for Wu & Moin (2008) and Chin et al. (2014). However, the data from El Khoury et al. (2013) and our DNS have azimuthal resolutions comparable to those in Ahn et al. (2013) and Pirozzoli et al. (2021), but produce results similar to Wu & Moin (2008) and Chin et al. (2014) -suggesting that the azimuthal resolution is not the main reason for such discrepancy. Interestingly, the data of Ahn et al. (2013) and Pirozzoli et al. (2021) are consistently lower than our and other results in the whole Re b range. In particular, at the highest Re b , the Pirozzoli et al. (2021) data undershoot by 2 % from λ p , but our data are lower by only approximately 0.6 %. Table 3 asserts that such differences in λ are beyond the uncertainty limit. Our DNS and the experimental data of Furuichi et al. (2015) agree well for the Re b range studied. Fitting our DNS data with (3.1) yields A = 2.039 ± 0.083, B = 0.948 ± 0.364 with uncertainty estimates based on 95 % confidence bounds, giving κ = 0.399 ± 0.015. The value reported by Pirozzoli et al. (2021) (i.e. A = 2.102, B = 1.148) are slightly larger than ours but still within the uncertainty range. However, large uncertainty is present in the fitted values due to the limited data points in Re. In addition, as Re b is still relatively low, the reported value of κ = 0.399 ± 0.015 should not be used outside of the given Re range. As will be shown in § 3.4, even for the highest Re b case, a distinct logarithmic region does not manifest itself in the mean velocity profile U( y). Higher Re data (e.g. Re τ ≥ 10 4 ) are required to better estimate the constants in (3.1) and the associated κ values.

Wall shear stress fluctuations
The Re-dependence of axial wall shear stress fluctuation τ 2 z,w + is one of the highly debated issues in wall turbulence. Note that τ 2 z,w + is also equivalent to the wall dissipation of the axial Reynolds stress component + z,w , the azimuthal vorticity variance at the wall ω 2 θ + or the limiting value of u 2 z + /U 2 at the wall (Örlü & Schlatter 2011). Previous DNS studies of channel, pipe and turbulent boundary layer observed an increase in τ 2 z,w + with Re τ , which reflects the increased contribution of large-scale motions to wall shear stress at high Re values (Marusic, Mathis & Hutchins 2010a). However, the exact dependence of τ 2 z,w + on Re τ is not well established. For example, Örlü & Schlatter (2011) suggested that the r.m.s. τ 2 where the two constants C and D are chosen as 0.298 and 0.018 based on the DNS of turbulent boundary layer data. Some works (Yang & Lozano-Durán 2017;Smits et al. 2021) also suggested that  (2015) NEK5000 2K where G is the asymptotic value at infinite Re, and H is a coefficient. The assumption for (3.4) is that the dissipation of the axial Reynolds stress component ( + z = ν (∂u + z /∂x + j )(∂u + z /∂x + j ) ) balances the turbulent kinetic energy production (P k = − u r u z + (∂U + /∂y + )) near the location of peak production. The fact that P k is bounded by 1/4 (Sreenivasan 1989;Pope 2000) implies that + z,w may also stay bounded, which is further assumed by Chen & Sreenivasan (2021) to be the same bound of P k , i.e. G = 1/4 (Chen & Sreenivasan 2021). This argument was later criticized by Smits et al. (2021) for the following two reasons. First, based on the DNS data, the location of peak production is actually the place where the largest imbalance of production and dissipation occurs. Second, as the balance between different terms in the Reynolds stress transport equation changes rapidly near the wall, it is unclear how the balance between P k and + z,w can be extended up to the wall, where P k = 0, and + z,w equals the viscous diffusion. The axial wall shear stress fluctuation τ 2 z,w + as a function of Re τ is depicted in figure 4(a). Akin to the observation in El Khoury et al. (2013), the values for the pipe are slightly lower than those for the channel, but the difference decreases with increasing Re τ . Note that the data of Pirozzoli et al. (2021) are slightly lower than others, which is consistent with the lower λ in figure 3. The logarithmic law (3.2) proposed by Örlü & Schlatter (2011) is higher than the DNS data. This is somehow expected as the fitting coefficients were obtained based on the DNS of turbulent boundary layer data, which are higher than in pipe and channel. In the high Re τ range, both the logarithmic (3.3) and defect-power-law scalings (3.4) agree well with the data. However, both scalings exhibit notable disagreements in the low Re τ range. The discrepancy seems not particularly surprising, given that the parameters in these equations are obtained from different datasets and different Re τ ranges. As a reference, the inset in figure 4(a) shows the fitting results For the Re τ range studied, the data seem to match better the defect-power law but with a slightly higher asymptotic value than suggested by Chen & Sreenivasan (2021). Additional data at higher Re τ are needed to confirm this finding. . Again, the defect-power law seems to match better the DNS data, but the agreement is not as good as for τ 2 z,w + .

Mean velocity profile
In an overlap region between the inner and outer flows, there is a logarithmic variation of the mean axial velocity U + profile, which is given as In a true log layer, the indicator function β = y + (∂U + /∂y + ) is constant and equals 1/κ. The U + profiles at different Re τ are compared to previous DNS data in figure 5(a). First, as expected, the wake in U + for the pipe is stronger than in the channel by Lee & Moser (2015). Second, our data in the outer region agree well with those of El Khoury et al. (2013), but not with Pirozzoli et al. (2021). This discrepancy was also noted by Pirozzoli et al. (2021), who found that their data, along with those of Wu & Moin (2008) and Ahn et al. (2013), differ from those of El Khoury et al. (2013) and Chin et al. (2014). This disparity is probably due to the numerical methods, where all of the former used low-order finite difference methods, while the latter two used high-order spectral-element methods.
The log-law diagnostics function β is shown in figure 5(b) for all high DNS data (i.e. Re τ > 5000). Interestingly, β for our Re τ ≈ 5200 agrees with the channel data of Lee & Moser (2015) and Hoyas et al. (2022) up to y + ≈ 250. Note that a relatively small domain size was used by Hoyas et al. (2022), but it is assumed large enough to yield accurate flow statistics (Lozano-Durán & Jiménez 2014). This suggests a near-wall universality of the inner scaled mean velocity -similar to that observed by Monty et al. (2009). For all three cases, the trough is β ≈ 2.30, located at y + ≈ 70. However, the data of Pirozzoli et al. (2021) deviate from others for y + > 40 and have a larger magnitude of the trough, which is somewhat consistent with the slight upward shift of the U + observed in figure 5(a). This discrepancy is significantly larger than the statistical uncertainty (see Appendix B). Unlike channel flow, where a plateau starts to develop for Re τ ≥ 5200, there is no plateau for pipe flow -suggesting that the minimum Re τ for U + to develop a logarithmic region should be higher in the pipe than in the channel. The β in the wake region of the pipe is distinctly larger than the channel, implying a notable difference in flow structures in the core region between these two flows (Chin et al. 2014).
High-order corrections to the log-law relation (3.5) were sometimes introduced to better describe the mean velocity profile in the overlap region (Buschmann & Gad-el Hak 2003;Luchini 2017;Cantwell 2019). For example, based on refined overlap arguments expressed by Afzal & Yajnik (1973), Jiménez & Moser (2007) proposed the indicator function where α 1 and α 2 are adjustable constants, and κ ∞ is the asymptotic von Kármán constant. Equation (3.6) allows for a Re-dependence of κ = κ ∞ + (α 1 /Re τ ) −1 and introduces a linear dependence on y. By fitting our Re τ ≈ 5200 data in the region between y + = 300 and y = 0.16, we obtain κ = 0.401 and α 2 = 2.0. This κ is very close to 0.399 estimated from the friction factor relation (3.1) and 0.402 reported by Jiménez & Moser (2007) using channel data of Re τ = 1000 from Del Alamo et al. (2004), and Re τ = 2000 from Hoyas & Jiménez (2006). It is slightly larger than 0.387 by Pirozzoli et al. (2021) and 0.384 by Lee & Moser (2015). In addition, α 2 is generally much larger in the pipe than in the channel -suggesting a strong geometry effect on β. The value of α 2 = 2 is consistent with the finding by Luchini (2017), who suggested that the logarithmic law of the velocity profile is universal across different geometries of wall turbulence, provided that the perturbative effect of the pressure gradient is taken into consideration. Furthermore, a good collapse between our data and the analytical prediction by Luchini (2017) is seen in figure 5(b).

Reynolds stresses
The non-zero components of the Reynolds stress tensor (or the velocity variances and covariance) are examined in this subsection (figures 6-10). For all datasets, the inner-scaled velocity variances and covariance increase with Re τ in the whole wall-normal range. In terms of the axial velocity variance u 2 z + (figure 6a), our data agree well with El Khoury et al. (2013) but differ from Pirozzoli et al. (2021), which is notably smaller in the near-wall region, particularly at low Re τ . For the highest Re τ cases, the agreement is reasonably good near the wall. However, note that the Pirozzoli et al. (2021) simulation is at a slightly higher Re τ (i.e. ∼ 6000). The differences between our case and that of Pirozzoli et al. (2021) can be better highlighted in the diagnostic plot (figure 6b), where the r.m.s. axial velocity fluctuation u z,rms is plotted against the mean velocity U + . The diagnostic plot was introduced by Alfredsson & Örlü (2010) to assess if the mean velocity and velocity fluctuation profiles behave correctly without the need to determine the friction velocity or the wall position. Consistent with the observation in Alfredsson & Örlü (2010), the diagnostic plot collapses in the outer region for Re τ > 180, and has a clear Re trend around the peak value. Most importantly, the data by Pirozzoli et al. (2021) are consistently lower than ours, particularly in the near-wall region. Such inconsistency is also observed for u 2 θ + (figure 10). The agreement for u 2 r + (figure 7a) is reasonably good among different pipe flow datasets, which is slightly larger than the channel, particularly in the outer region. The Reynolds shear stress u r u z + shows excellent agreement among different datasets, even including the channel.
Let us focus now on the inner peak of the axial velocity variance u 2 z + . The inner peak is assumed to increase logarithmically with Re τ -similar to the wall shear stress fluctuations due to the increased modulation effect of the large-scale structures in the logarithmic layer (Marusic & Monty 2019). Chen & Sreenivasan (2021) suggested recently that the growth of u 2 z + would eventually saturate. The argument is based on the balance between the viscous diffusion and dissipation at the wall, and the Taylor series expansion of the axial velocity variance near the wall, given as (3.7)  independent of Re τ , then (3.7) suggests that the peak of axial velocity variance should also be bounded in a defect-power form similar to the wall shear stress fluctuation: where M is the asymptotic value and N is the coefficient. This validity of (3.8) was challenged recently by Pirozzoli et al. (2021), who, based on their data, observed a slight increase of y + z,p with Re τ , from y + z,p = 14.28 for Re τ ≈ 500 to 15.14 for Re τ ≈ 6000. We emphasize that such variation of y + z,p with Re τ is not observed in our case, where much finer near-wall resolutions than in Pirozzoli et al. (2021) are used. The value of y + z,p is approximately 15 for all Re τ (e.g. y + z,p = 15.07, 15.03, 15.50 for Re τ = 180, 2000, 5000, respectively) -akin to the findings by many others (Moser et al. 1999;Jiménez et al. 2010;Chin et al. 2014;Smits et al. 2021).
Figure 6(c) shows u 2 z + p for all the DNS data listed in table 2, along with the logarithmic law u 2 z + p = 3.8 + 0.64 ln(Re τ ) by Marusic et al. (2017) and the power law u 2 z + p = 11.5 − 19.32Re −1/4 τ by Chen & Sreenivasan (2021). The difference between different DNS datasets is relatively small, except for those from Pirozzoli et al. (2021), which are consistently lower than others for all Re τ . Note that this discrepancy is much larger than the uncertainty (standard deviation), which is less than 0.5 % (see table 3). Both the logarithmic and defect-power laws fit well with the data in the high Re τ range but have certain discrepancies at low Re τ . This suggests that there might exist a transitional scaling -similar to that found for the Reynolds shear stress (Chen, Hussain & She 2019). The parameters in these two scaling laws can be adjusted to better fit our dataset. The inset shows the fitting results for our data without the one at Re τ = 180: u 2 z + p = 3.251 + 0.687 ln(Re τ ) and u 2 z + p = 11.132 − 17.402 Re −1/4 τ . With these new constants, the agreement is improved for both scaling laws. In summary, for the Re τ range studied, both scaling laws provide a good match with u 2 z + p data when the fitting parameters are properly adjusted. Data at even higher Re τ are required to determine which law is more consistent with the data.
According to Townsend's attached eddy hypothesis, at sufficiently high Re, the Reynolds stress components in a certain y range satisfy where A i and B i are universal constants. Consistent with these relations, the radial velocity variance u 2 r + slowly develops a flat region as Re τ increases. In addition, the Reynolds shear stress −u r u z + profiles also tend to become flattened at higher Re τ . As noted by Afzal (1982), the peak Reynolds shear stress at high Re τ follows −u r u z κ Re τ , and the corresponding position y + m shifts away from the wall following y + m ∼ √ Re τ /κ. Chen et al. (2019) suggested that there is a non-universal scaling transition, where the peaks at low Re τ scales as u r u z + p ≈ Re −2/3 τ and their locations scales as y + m ∼ Re 1/3 τ . Figure 8 shows −u r u z + p and the corresponding y + m as a function of Re τ . For Re τ > 1000, Re −1/2 for −u r u z + p and Re −1/2 for y + m are satisfied with good accuracy, and at the low Re τ range, Re ) -suggesting that no clear logarithmic region develops for the Re τ range considered. As discussed in Lee & Moser (2015), Re τ = 5200 is not quite high enough to exhibit such a region. Based on the Superpipe data, Marusic et al. (2013) suggested that a sensible logarithmic layer emerges only for Re τ > 10 4 . Consistent with the findings in Lee & Moser (2015) and Pirozzoli et al. (2021), the azimuthal velocity variance u 2 θ + develops the logarithmic layer at lower Re τ . The indicator function of u 2 θ + shows a distinct plateau (figure 9b). Interestingly, when compared with other cases, the range of the logarithmic layer is wider in our Re τ ≈ 5200 case. Fitting the data in the range 120 ≤ y + ≤ 800 yields A 3 = 0.921, B 3 = 0.420, which are close to A 3 = 1.0, B 3 = 0.40 by Pirozzoli et al. (2021), and between A 3 = 1.08, B 3 = 0.387 obtained by Lee & Moser (2015) for the turbulent channel and A 3 = 0.8, B 3 = −0.45 by Sillero, Jiménez & Moser (2013) in the turbulent boundary layer. Figure 11(a) shows the production P + k and dissipation + k (= ν (∂u + i /∂x + j )(∂u + i /∂x + j ) ) of the turbulent kinetic energy (i.e. k + = ( u 2 r + + u 2 θ + + u 2 z + )/2). Other terms in the transport equations of the turbulent kinetic energy and individual Reynolds stress components are available at https://dataverse.tdl.org/dataverse/turbpipe. The production P + k has a peak at around y + ≈ 11, and the magnitude approaches the asymptotic value 1/4 as Re τ increases. Despite notable differences in the mean axial/streamwise velocity profile observed in the outer region, P + k is quite similar between pipes and channels. This explains why the higher velocity gradient of the pipe does not contribute an effect to the turbulence intensities. The magnitude of dissipation + k increases continuously with Re τ , 956 A18-17  Figure 11. (a) Production P + k and dissipation + k of the turbulent kinetic energy, and (b) balance of production and dissipation P + k / + k − 1, as functions of y + . and the difference between the pipe and channel is mainly in the near-wall region and decreases with increasing Re τ . At sufficiently high Re τ , there is an intermediate region where production balances dissipation. Recent numerical results (Lee & Moser 2015;Pirozzoli et al. 2021) suggest that such equilibrium between production and dissipation is violated due to the presence of LSMs and VLSMs. Figure 11(b) shows the relative excess of production over dissipation (P + k / + k − 1). First, there is a near-wall region (8 ≤ y + ≤ 35) where P + k distinctly exceeds + k . At high Re τ , another region of P + k / + k > 1 develops, and the magnitude increases with Re τ . For the Re τ ≈ 5200 case, the peak imbalance is approximately 11 % (located at y + ≈ 330), which is slightly larger than in the channel (i.e. 8 %).

Mean pressure and r.m.s. pressure fluctuation
The mean pressure and r.m.s. pressure fluctuations are displayed in figure 12. First, the mean pressure P + has different behaviour in the outer region between pipe and channel flows, with P + being substantially lower in the wake of the pipe. As discussed in El Khoury et al. (2013), this difference is related to the mean radial momentum equation, which, in pipe flow, is given as (Hinze 1975) By changing variable (i.e. r = R − y) and then integrating the above equation, the mean pressure for pipe flow with the wall value set to zero can be expressed as (3.14) In channel flow, the last term on the left-hand side of (3.13) is absent, and the mean pressure is solely balanced by the wall-normal velocity fluctuation, i.e. P + ( y) = − v 2 + . From figure 7(a), it is clear that the wall-normal velocity fluctuation is comparable between pipe and channel flows. However, as u 2 r + < u 2 θ + in pipe flow, the extra term in (3.14) is zero at the wall and decreases with increasing y -resulting in a lower pressure in pipes than in channels. Similar to that observed by El Khoury et al. (2013), the r.m.s. pressure fluctuation p + rms exhibits similar behaviour between the pipe and channel, except for slightly lower values for the latter. Minor differences are observed between our data and those of Pirozzoli et al. (2021), particularly near the peak value. The difference between our data and the channel data of Lee & Moser (2015) in the near-wall region decreases with increasing Re τ . Figure 13 further shows the peak and wall values of p + rms , which has similar Re-dependence as for other measures, such as wall shear stress fluctuations and axial velocity variances. Again, for the Re τ studied, both the logarithmic and defect-power laws fit the data well.
peaks coincide with those found by Lee & Moser (2015) in the channel at the same Re τ , but the physical scales are smaller than in the channel. It is well known that the inner peak at y + = 13 is associated with the streaks that are generated through the near-wall self-sustaining cycle (Waleffe 1997;Schoppa & Hussain 2002). As seen frequently in experiments (Hutchins & Marusic 2007;Monty et al. 2009;Rosenberg et al. 2013), the outer peak results from VLSMs. The outer peak in the θ-direction (at y = 0.192) is located farther away than the streamwise one (at y = 0.077), which according to Wu, Baltzer & Adrian (2012) suggests that the VLSMs in the outer region maintain their energy in the θ-direction more strongly than in the z-direction. The pre-multiplied energy spectra of the Reynolds shear stress in axial (k z E rz /u 2 τ ) and azimuthal (k θ E rz /u 2 τ ) directions as functions of y + are shown in figures 14(c,d), respectively. The inner peak is located at y + = 30 with k z R = 49.2 (λ + θ = 664) for k z E rz /u 2 τ , and k θ = 268 (λ + θ = 120) for k θ E rz /u 2 τ . Compared with the axial velocity spectra, although the wavelength of the outer peak remains identical, the magnitude is much weaker and farther away from the wall. Figure 15(a) shows the one-dimensional pre-multiplied energy spectra k z E zz /u 2 τ at different y locations. For comparison, the channel data of Lee & Moser (2015) at the same Re τ are also included. First, good agreement is observed at y + = 15 between channel and pipe, particularly at higher wavenumbers -suggesting insignificance of pipe curvature on fine-scale near-wall structures. The scaling analysis of Perry, Henbest & Chong (1986) suggests that the energy spectral density of the axial velocity fluctuations k z E zz /u 2 τ should 956 A18-20  Figure 15. Comparison of the pre-multiplied energy spectra between pipe (solid) and channel (dashed) for Re τ = 5200: (a) k z E zz /u 2 τ , and (b) k θ E zz /u 2 τ .
vary as k −1 z in the overlap region. The k −1 z region has previously been observed in the high-Re experiments (Nickels et al. 2005;Rosenberg et al. 2013). Recently, such k −1 z has also been discovered in DNS of turbulent channel flow at Re τ = 5200 (Lee & Moser 2015). Similarly, a plateau in the region 6 ≤ k z R ≤ 10 is observed for 90 ≤ y + ≤ 170, and the magnitude 0.8 agrees with experiments (Nickels et al. 2005;Rosenberg et al. 2013). A bimodal is observed for y + = 90, with the peak magnitude at low wavenumbers (k z R = 1) being smaller than at high wavenumbers (k z R = 30). Interestingly, k z E zz /u 2 τ values at low wavenumbers are slightly smaller in the pipe than in the channel. Figure 15(b) further shows the one-dimensional pre-multiplied energy spectra k θ E zz /u 2 τ at different y locations. Again, k θ E zz /u 2 τ agrees well between the pipe and the channel at high wavenumbers. Consistent with those in the channel, a plateau appears for 5 ≤ k θ ≤ 30 in the overlap region, with the magnitude increasing with y + . Such a plateau is present even in the viscous sublayer, which is the footprint of LSMs and VLSMs near the wall (Mathis, Hutchins & Marusic 2009;Hwang et al. 2016).

Concluding remarks
A new direct numerical simulation providing reliable high-fidelity data of turbulent pipe flow for Re τ up to 5200 is presented. Particular focus has been put on providing data that are as accurate as possible, by using a high-order numerical method, large domains and sufficient integration time with quantified uncertainty. The DNS are performed with a pseudo-spectral code OPENPIPEFLOW, and the axial extent of the domain is 10πR (where R is the pipe radius), which can be considered sufficiently long to capture all the relevant structures. A wealth of statistical data with uncertainty, including mean velocity, Reynolds stress and their budgets, pressure and its fluctuations, and energy spectra, are gathered (available online at https://dataverse.tdl.org/dataverse/turbpipe). An extensive comparison between our new pipe data and other simulation and experimental data is made, and small but still substantial and systematic differences between the various datasets are identified. For example, consistent lower values of the friction factor, wall shear stress fluctuations, and the inner peak of the axial velocity variance are observed for data generated using low-order methods, such as Ahn et al. (2013) and Pirozzoli et al. (2021). In pipe flow simulation, the only parameter apart from Re is the length of the pipe. Once the latter is chosen large enough, all data should, in principle, be the same.
Such discrepancies between different simulations thus highlight the need for high-order accurate methods for this particular flow case. This argument is further complemented by performing additional DNS at Re τ = 2000 with a spectral-element code NEK5000, where all the statistical data generated are found to match well the results obtained using OPENPIPEFLOW.
Different from turbulent channel flow, the mean velocity has not yet developed a logarithmic region at Re τ = 5200, yet the diagnostic function collapses well between our Re τ ≈ 5200 and the channel data of Lee & Moser (2015) and Hoyas et al. (2022) up to y + ≈ 250 -suggesting a near-wall universality of the inner scaled mean velocity. The wall shear stress fluctuations, the inner peak of axial velocity variance, and the wall and peak of r.m.s. pressure fluctuations continuously increase with Re τ , and their difference between pipe and channel decreases with increasing Re τ . In addition, at the Re τ range considered, the Re dependence of these quantities agrees with both the logarithmic and defect-power scaling laws (Chen & Sreenivasan 2021). Consistent with observations in the channel, the one-dimensional spectrum of the axial velocity exhibits a k −1 dependence at intermediate distances from the wall.

Appendix B. Estimated uncertainties for one-point statistics
We explain briefly the approach employed to estimate the uncertainties in the mean velocity and Reynolds stress components of the DNS of the pipe flow. During the simulations, the time samples of the quantities contributing to statistical terms are averaged over the azimuthal (θ) and axial (z) directions. To compute the central statistical moments, the temporal correlations between the spatially averaged quantities are preserved by, for instance, writing a Reynolds stress component as u r u z = u r u z −ū rūz , where overbar means averaging over θ and z. In practice, the sample-mean estimator (SME) is used to estimate a from a finite number of time series samples {a i } n i=1 , where a i = a(t i ) are equispaced time samples. The SME for a is defined aŝ whereÊ[a] is the estimated expectation of a. Based on the central limit theorem, for a sufficiently large number of samples, the SME converges to the true expectation via a Gaussian distribution,μ a ∼ N (μ a , σ 2 (μ a )).
To estimate σ (μ a ) and hence quantify the time-averaging uncertainty inμ a , an analytical expression can be derived that depends on the autocorrelation of time series a at different lags; see e.g. Oliver et al. (2014) and Rezaeiravesh et al. (2022). To avoid inaccuracy in σ (μ a ) due to the oscillations in the sample-estimated autocorrelations, especially at higher lags, an autoregressive model is first fitted to the samples {a i } n i=1 , which is then used to construct a smooth model for the autocorrelations. The details of the approach can be found in Oliver et al. (2014) andD. Xavier et al. (personal communication 2022). The main hyperparameters of this uncertainty quantification (UQ) approach are the order of the autoregressive model and the number of training lags when modelling the autocorrelation function (ACF). In the present study, the optimal values of these hyperparameters are chosen based on the sampling frequency of the data at each Re. All estimations of uncertainties have been performed using UQit (Rezaeiravesh, Vinuesa & Schlatter 2021).
Following the above approach, the uncertainty in the statistical moments of any order can be estimated accurately. However, there are various turbulence statistics that are defined as a combination of the exponents of various moments; for instance, consider the turbulence intensity, r.m.s. fluctuations, turbulent kinetic energy, and various terms in the transport equations of the Reynolds stress components. The uncertainty in such terms can be estimated by applying the approach described in Rezaeiravesh et al. (2022). The main idea for estimating the uncertainty in a compound statistical term is to estimate the uncertainty in its constitutive statistical moments and also estimate the cross-covariance between the SMEs corresponding to them. Following this procedure, in the DNS database reported online in connection with the present study, all statistics are accompanied by an accurate estimation of the corresponding time-averaging uncertainty. An important aspect of this procedure is that for the statistics expressed in wall units, the uncertainty of the wall friction velocity is also taken into account. This means, for instance, that for u r u z + , the uncertainty of both u r u z and u τ 2 are considered applying a Monte-Carlo-based UQ forward problem, which does not require any linearization. Table 4 summarizes the sampling interval and the total number of samples used for UQ for simulations at different Re τ . Our investigation showed that for the collected samples, an autoregressive model of order 20 along with the sample-estimated ACF at the first 20 lags, for Re τ = 180, 550 and 1000, and 40 lags, for Re τ = 2000 and 5200, leads to accurate models for autocorrelation of various quantities. For low-order moments, using sample-estimated ACF at a higher number of lags, especially near the centre of the pipe, could lead to slightly more accurate models for ACF. However, the difference in the resulting estimated uncertainty is below 1 %. Figure 17 shows the standard deviation σ (see (B2)) of the sample estimation of different inner-scaled statistical terms. Clearly, the estimated uncertainties vary between the moments and also in the wall-normal direction. However, for all quantities, the lowest uncertainty (corresponding to highest certainty) is observed near the wall. Moreover, the estimated uncertainty for each quantity exhibits a similar variation in the wall-normal direction for different Re τ .