Large Eddy Simulations of Turbulent Pipe Flows at Moderate Reynolds Numbers

Wall-bounded turbulence is relevant for many engineering and natural science applications, yet there are still aspects of its underlying physics that are not fully understood, particularly at high Reynolds numbers. In this study, we investigate fully-developed turbulent pipe flows at moderate-to-high friction velocity Reynolds numbers (361 ≤ Re τ ≤ 2 , 000), corresponding to bulk velocity-based Reynolds numbers of 11 , 700 ≤ Re b ≤ 82 , 500, using wall-modeled Large Eddy Simulations (LES) in OpenFOAM. A grid convergence study is performed for Re τ = 361, followed by an investigation of the accuracy of various subgrid-scale stress models for the same Reynolds number. Results show that the Wall-Adapting Local Eddy (WALE) model performs well compared to experiments and Direct Numerical Simulations (DNS), while One-Equation Eddy-Viscosity Model (OEEVM) and Smagorinsky (SMG) are too dissipative. LES utilizing WALE are then performed for four different Reynolds numbers with gradually refined grids, revealing excellent agreement with DNS data in the outer region. However, a significant deviation from DNS data is observed in the sub-viscous layer region, indicating the need for further mesh refinement in the wall-normal direction to accurately capture the smallest-scale motions’ behavior. Additional mesh sensitivity analysis uncovered that, as the Re τ value rises, it becomes crucial for a grid to adhere to the condition of ∆ x + ≤ 20 − 25 and ∆ z + ≤ 10 in order to precisely capture substantial large and small scale fluctuations. Overall, the WALE model enables accurate numerical simulations of high-Reynolds-number, wall-bounded flows at a fraction of the computational cost required for temporal and spatial resolution of the inner layer.


I. INTRODUCTION
Two common features of engineering and environmental fluid flows are wall-boundedness [1, 2] and high Reynolds number (Re) turbulence [3][4][5].The study of turbulent flow traces back over a century, and its widespread application has kept this field under constant development [6].The spatially evolving boundary layer, the channel, and the pipe are three canonical building block flows for wall bounded turbulent flows.Consequently, pipe flow is fully described by Re and the length of the pipe.The effect of the pipe length is insignificant only if considered sufficiently large [7].During the transport of fluids through pipes, roughly half of the energy is dissipated by turbulence near the walls [8].The need to investigate fundamental wall-bounded turbulence at high Re is widely agreed in literature.Clauser [9] and Coles & Hirst [10] provided a comprehensive review of classical wall-bounded turbulent scaling and structure for the case of a smooth wall.However, one of the most recent revealing research result comes from efforts seeking to achieve an everhigher Re number and, in doing so, shedding light on some of the longstanding questions about the very nature of wallbounded turbulence, including how well our present analytical and computational models match its complex behavior [3].
To date, accurate high-Re-number data are still primarily achievable only via well-designed experiments [11].With the progress of computational technology, simulations have emerged to play a vital role in research on turbulent flows and wall-bounded turbulence.One of the most cited articles on turbulent pipe flow is by Eggels et al. [12], showing compar-isons between Direct Numerical Simulation (DNS) and experimental results of Hot Wire Anemometry, Particle Image Velocimetry, and Laser Doppler Anemometry at nominally similar friction-velocity-based turbulent Re τ = u τ h/ν ≃ 180, where u τ the friction velocity, h the characteristic length, and ν the kinematic viscosity.Numerous studies (e.g., [13][14][15][16]) have extensively documented the effects at low Re τ .These studies illustrate the augmentation in the maxima of streamwise turbulence intensity and Reynolds shear stress, along with inner-scaled turbulence statistics correlated with Re τ .These effects are predominantly observed in the logarithmic law of the mean velocity profiles, attributed to higher pressure gradients.Nevertheless, decades of experimental research have shown that obtaining detailed high-Re τ data, notably near the wall, remains challenging.That is because as Re τ increases, the smallest length scales decrease, and as a result lead to significant uncertainties in determining the probe locations and turbulence intensities.Recently, DNS of turbulent pipe flows with Re τ ranging from 500 to 6,000 have been presented [7,[17][18][19][20][21][22][23][24].These Re τ are considered modest compared to turbulent pipe flow experiments carried out with a turbulent Re range of Re τ ≃ O(10 3 − 10 5 ) [11,[25][26][27][28].
As Re τ increases, the range of turbulent scales extends, demanding higher resolutions for DNS.However, the number of required grid points for DNS scales roughly as O(Re 9/4 τ ) [29], presently falling short of computing capability at engineering Re τ values.Capturing the turbulent flow and accounting for additional computing overhead requires computer resources that scale nominally as Re 4  τ [30].Considering that many real-world flows have Re τ > 10 6 , the computing costs associated with DNS quickly becomes prohibitive.
On the other hand Large Eddy Simulations (LES) takes a different approach, by filtering of the smallest scales and replacing their effects on the large resolved scales by subgrid models [31].LES of unbounded homogeneous and shear turbulence is well established [32][33][34][35], whereas LES of wallbounded flows has remained challenging, firstly, because the presence of a wall causes the length scales to progressively decrease toward the wall, and secondly, owing to the nearwall anisotropy.
When a wall is present the LES grid is typically constrained by the near-wall region, where turbulent and viscous momentum transport occur on small scales.Piomelli & Balars [36] showed that for near-wall-resolved LES, around 99% of grid points are used to resolve the wall layers in fullydeveloped channel flow.As a result, for wall-bounded flows the computational cost of near-wall-resolved LES to resolve the viscous sublayer is estimated to scale as Re 1. 8  τ , which demonstrates the potential savings over DNS [36].Nevertheless, the development of wall-modeled LES [36][37][38][39][40] has widened the range of applications of LES for wall-bounded flows.In near-wall-modeled LES, the effects of the nearwall anisotropic fine scales are not resolved but modeled by wall models, approximating the near-wall physics.Nearwall-modeled LES is computationally more efficient than near-wall-resolved LES since only the outer layer is explicitly computed, with the required grid size determined by the outer flow eddies [36,[41][42][43][44]. Thus in wall-modeled LES, the grid requirements could be relaxed; it is estimated that the computational costs would scale as Re 0.2 τ [45].Sagaut [31] gives a thorough description of the (still on-going) efforts in the formulation of more accurate LES models.Yet few simulations have been performed at high Re.
To the best of our knowledge, the highest turbulent Re τ used in LES of turbulent pipe flows were Re τ ≃ 2,000 and 180,000, by Berrouk et al. [46] and Ferrer et al. [47].The authors present explicit and implicit LES of fully-developed turbulent pipe flows using a continuous Galerkin spectral element solver, motivated by the work of Chung & Pullin [40].For explicit LES they utilized the explicit stretchedvortex model of Chung & Pullin [40].From the literature review, it is evident that only a limited number of studies have been conducted with wall-modeled LES focusing on the moderate-to-high Re τ wall-bounded turbulent flow.The objective of the present study is to carry out wall-modeled LES for moderate to high Re τ and compare the results to DNS at matched Re τ for internal wall-bounded flows.To perform wall-modeled LES, different meshes have been generated, with y + values from 4 to 11.We study the performance of different subgrid models, investing their accuracy at predicting the statistical behavior of moderate Re fully-developed turbulent pipe flows.Five subgrid models are considered: the Wall-Adapting Local Eddy viscosity model (WALE) [48], the Smagorinsky model (SMG) [49], the One-Equation Eddy Viscosity Model (OEEVM) [50], the Localized Dynamic Kinetic energy Model (LDKM) [51,52], and the Differential Stress Equation Model (DSEM) [53].Only the WALE model is selected to investigate the dependence on grid resolution, as a trade-off between computational cost and numerical accuracy.For detailed investigation, four friction-velocity-based turbulent Re are considered: Re τ = 361, 550, 1,000, and 2,000.The simulations are evaluated based on predictions of zeroth-, first-, and second-order turbulence statistics and flow features such as friction coefficient, mean velocity profiles, velocity correlations, mean turbulent shear stress profiles, quadrant analysis, and anisotropy invariant maps.Numerical results for each model are compared with DNS data [7,23].

A. Mathematical formulation
In LES, the governing equations are the low pass filtered conservation equations of incompressible flow mass and momentum as represented by, where u i is the filtered velocity, p, the filtered pressure, ν, the kinematic viscosity, and τ sgs i j = u i u j −u i u j , the subgrid stress tensor.Mathematical analysis and physical arguments, e.g., [54][55][56], suggest that: (i) τ sgs i j is invariant under a change of frame; (ii) τ sgs i j is positive-definite symmetric, provided that the filter function, G ∆ , is symmetric; and (iii) the inequalities k = 1 2 tr(τ , and det(τ sgs i j ) ≥ 0 must be satisfied for τ sgs i j to be positive-definite.

B. Subgrid modeling
Two methodologies are usually adopted for subgrid modeling: functional modeling, which consists in modeling the action of the subgrid scales on the resolved scales, and structural modeling, whose principle is to model the subgrid stresses without incorporating any knowledge about the interactions between the subgrid and resolved scales [31].However, for the purpose of this paper, we instead prefer to focus on isotropic and anisotropic subgrid models since complex flows with high Re numbers are often characterized by anisotropy on a wide range of scales -typically reaching the range of scales that require modeling, e.g., in the nearwall region.The most frequently used subgrid models are the eddy-viscosity models.This type of models computes the deviatoric part of the subgrid stress tensor by relating it, in a linear correspondence, to the rate-of-strain tensor, S i j = 1 2 (∂ u i /∂ x j + ∂ u j /∂ x i ) [57].By doing so, the subgrid stress tensor may be thought of being decomposed into an isotropic part and an anisotropic part.The former is also known as the normal stresses and is directly related to the turbulent kinetic energy while the latter measures the deviations from isotropy.Hence, the subgrid stress tensor, based on the Boussinesq hypothesis, is expressed as where ν sgs is the subgrid eddy viscosity and δ i j is the Kronecker delta.Models for the eddy viscosity, ν sgs , and the subgrid turbulent kinetic energy, k sgs = 1 2 tr(τ sgs ii ), are needed to close Eq.(3).For this purpose, we assume the existence of characteristic length and velocity scales and we infer separation between the resolved and subgrid scales.
In this work, we limit ourselves to five well-known LES subgrid models: Firstly, the One-Equation Eddy Viscosity Model (OEEVM) [50], in which a modeled transport equation for k sgs is solved, and ν sgs = c k ∆k 1/2 sgs is modeled, where c k = 0.094 and c ε = 1.048.
The second model is the Smagorinsky (SMG) model [49] which assumes a local equilibrium between production and dissipation of k sgs .As a result, the subgrid viscosity and kinetic energy becomes, ν sgs = c 2 s ∆ 2 (2S i j S i j ) 1/2 and k sgs = c I ∆ 2 (2S i j S i j ), where c k = 0.094, c ε = 1.048, c I = c k /c ε = 0.089, and c s = c k c k /c ε are the model coefficients [49,58].
The third model is the Localized Dynamic Kinetic energy Model (LDKM) based on the work of Kim and Menon [51,52].In the LDKM, scale similarity between the subgridscale stress tensor and the test-scale Leonard stress tensor is assumed to evaluate the model coefficients, for ν sgs and the subgrid dissipation, ε.The transport equation used in this model is similar to that in OEEVM model.The expression for ν sgs is similar to the one of the OEEVM model: ν sgs = c k ∆k 1/2 sgs , except that the model coefficient, c k , is now evaluated dynamically as c k = 1 2 , where M i j = − ∆K 1/2 S i j , L i j = u i u j − u i u j , K = 1 2 L ii , and the tilde symbol, (.), denotes the explicit filtering operation, usually defined as twice the original filter, e.g., ∆ = 2∆.The subgrid dissipation is modeled as ε ≈ c ε K 3/2 / ∆ in which c ε is evaluated dynamically as c ε = g i j g i j − g i j g i j [(ν+νsgs) ∆] −1 K 3/2 , with g i j = ∂ u i /∂ x j being the filtered velocity gradient.
The fourth model is the Wall-Adapting Local Eddy viscosity (WALE) model [48].This method allows to obtain a subgrid viscosity proportional to the wall-normal distance cubed (ν sgs ∝ y 3 ), and additionally allows the subgrid viscosity to vanish near the solid walls.The WALE model can be defined via the following relation for subgrid viscosity, ν sgs = c w ∆ 2 (S d i j S d i j ) 3/2 (S i j S i j ) 5/2 +(S d i j S d i j ) 5/4 , where S d i j = 1 2 (g 2 i j + g 2 ji ) − 1 3 δ i j g 2 i j is the traceless symmetric part of the square of the velocity gradient tensor.
Lastly, we consider the Differential Stress Equation Model (DSEM) by Deardorff [53], which does not rely on the local isotropy assumption within the subgrid scales.The DSEM employed here uses a modeled transport equation for the sub-grid stress tensor, τ i j , following [53], The fourth, fifth, and sixth terms on the right-hand side of Eq. 4, which represent triple correlation, pressure-strain correlation, and the dissipation term, are modeled, and details can be found in [31,59].
The van Driest wall-damping function is used in the wall treatment for the SMG and OEEVM subgrid models.Whereas, a wall-model is instead used for the WALE, DSEM, and LDKM models.This wall model is based on the Spalding's law of the wall from which the friction velocity can be estimated based on which a wall viscosity can be derived at the first off-wall grid point, following the ideas of [41].

III. SIMULATION DETAILS
The LES equations are solved numerically using the finite volume method with a second-order cell-centered discretization scheme using OpenFOAM.Time stepping is performed using the second-order accurate implicit Adams-Bashforth method [60] with a maximum Courant-Friedrichs-Lewy (CFL) number of 0.5 for numerical stability.The PISO algorithm for LES is used to decouple the pressuremomentum equations [61].The resulting pressure equation is solved using a Generalized Geometric-Algebraic Multigrid (GAMG) method with a Diagonal Incomplete Cholesky (DIC) smoother [62].Velocity, turbulent kinetic energy, and dissipation rate are solved using a Preconditioned BiConjugate Gradient (PBiCG) solver [63] with a Diagonal Incomplete Lower-Upper (DILU) preconditioning method [62].
The turbulent bulk Re b number based on the pipe diameter, D, and bulk velocity, U b , is given by yielding the four values Re b = 11,700, 19,000, 37,700, and 82,500 used in this study.The motivation behind selecting these Re b values is the well-established DNS results previously reported [7,23].To provide a direct comparison of LES results with DNS, we use the same values of Re b as reported in these papers.In addition to the bulk velocity, another characteristic velocity scale is commonly introduced in connection with wall-bounded flows.The friction velocity, u τ , is defined in terms of the wall shear stress, τ w , and the fluid density, ρ, as The turbulent friction Re τ number, based on the pipe radius, R, and friction velocity, u τ , is given by yielding the four values Re τ = 361, 550, 1,000, and 2,000 used in this study.Hexahedral butterfly grids, have been used for all simulations.These grids are composed of a Cartesian part at the center of the pipe surrounded by four transitional regions from the inner square to the outer circle of all cross-sections of the pipe.Such grids thus require the definition of a partition of the pipe into five blocks prior to the grid generation, but generally provide the best grid quality, i.e., high orthogonality and low grid density.The generation of these grids is automated by means of an m4 macro [64] in OpenFOAM.The length of the computational domain is L. The radial, axial, and azimuthal directions are represented by r, x, and θ , and the corresponding velocity components, by u, v, and w.The wall-normal distance, i.e., along the radial direction, is denoted by y = R − r, where R = D/2 is the pipe radius.Here, u ′ , v ′ , and w ′ are the velocity fluctuations in the streamwise, radial, and azimuthal directions, respectively.
Since the fully-developed turbulent pipe flow considered here is homogeneous in the streamwise direction, periodic boundary conditions are imposed in this direction.No-slip boundary conditions are imposed at the pipe wall for all velocity components, whereas zero-Neumann boundary condition is used for the pressure.An additional external force term is introduced into the NSE to mimic the effect of the pressure gradient in the context of periodic boundary conditions in the streamwise direction.The magnitude of this force is determined by the bulk velocity.For the present study we chose the pipe length to be five times the pipe diameter, i.e., L = 10R.Indeed, to ensure that statistics are not affected by the streamwise periodic boundary conditions, we have investigated one-dimensional temporal and spatial autocorrelation functions of velocity fluctuations.The results are presented in Appendix B and show that the chosen domain length is sufficiently large for velocity fluctuations to uncorrelate.All flow statistics are averaged over approximately 100 to 120 flow-through times.

IV. RESULTS AND DISCUSSIONS
This discussion is mainly subdivided into three main sections.In Sec.IV A simulations for statistically fullydeveloped turbulent pipe flow at Re τ = 361 are presented.These simulations correspond to four different mesh resolutions.The results are grouped according to the computed quantites, e.g., mean flow profiles, Reynolds shear stresses, and friction factor.In Sec.IV B fully-developed turbulent pipe flow simulations at fixed Re τ = 361, for five LES models are performed.In Sec.IV C LES-WALE simulations for four moderate-to-high turbulent Re τ are carried out.Again, in both sections results are grouped based on the computed quantities.
A. Mesh independence study of fully-developed turbulent pipe flows LES-WALE simulations for pipe flows are perfromed at Re τ = 361 used four different meshes, named M1, M2, M3, and M4, with cell counts varying from 0.13 to 5 million.Table I provides details about these meshes, such as total cell count, cells along the flow direction, maximum skewness and aspect ratio, and ∆y + ranges.
Figure 1 shows the resulting flow using isosurfaces of the second invariant of the velocity gradient tensor, λ 2 [65], obtained through the LES-WALE simulations at different mesh resolutions.Vortices are identified using the λ 2 -criterion and presented for λ 2 /λ 2 max = −0.01,where more negative values indicate stronger vortices.The pipe shows a dense collection of turbulent eddies from the wall up to y/R ≃ 0.5.Coarser grids (M1 and M2) exhibit fewer vortical structures in the boundary layers, forming elongated axial vortices with length scales around 0.1R.Higher mesh resolution helps capture finer structures near the wall, reducing numerical dissipation.Despite differences in mesh, fluctuations appear similar from M1 to M4. Meshes M3 and M4 show similar turbulence statistics, suggesting convergence and independence from mesh sensitivity.
In Fig. 2, (a) shows mean velocity profiles in inner variables for four grid resolutions and includes DNS data from El Khoury et al. [7] in a pipe flow.The profiles are computed along a vertical line through the pipe center.Higher mesh resolution causes a downward shift in the overlap region compared to DNS data.Although there is a slight difference in the inner layer for M1, the profiles maintain a similar shape across different mesh resolutions.Unlike turbulent channel flows in experiments or numerical simulations [66][67][68][69], the mean velocity profiles do not align with the logarithmic velocity distribution, ⟨u⟩ + = (1/κ) log(y + ) + B + for y + > 200, where κ = 0.4 is the Von Kármán constant and B + = 5.2.This discrepancy is likely due to a wake region near the pipe centerline [7].Additionally, (b) tests the law of the wall by plotting mean profiles in velocity-deficit form in Fig. 2. Differences in mesh resolution become more evident.Transitioning from M1 to M4, profiles shift toward DNS data, and excellent overall agreement with DNS is achieved for M3 and M4.
Figure 3 assesses mesh sensitivity by presenting streamwise, wall-normal, and azimuthal root-mean-squared velocity profiles against the wall distance with inner scaling.As mesh resolution increases, results progressively align with DNS data (dashed black line).Notably, M1 and M2 show different flow behavior in the outer layer (y + > 50) compared to M3 and M4.The finest grid yields excellent agreement with DNS across the entire region, while M3 and M4 produce similar results despite their resolution difference, indicating statistical independence for these two grids.The positive impact of grid refinement is evident in the gradual increase in result accuracy.Additionally, root-meansquared velocity fluctuations, u + rms , v + rms , and w + rms , converge toward DNS results with grid refinement.The peak value of u + rms consistently occurs at y + ≃ 15, observed in various studies across wall-bounded flows.In Fig. 3 (d), pro-  τ .As mesh resolution increases, results increasingly agree with DNS data (dashed black line).M3 and M4 stand out, providing exceptionally good results, while M1, being a very coarse mesh, leads to the underprediction of most turbulent properties.
To conclude the mesh sensitivity study, friction factor ( f D ) values are calculated using f D = 8τ w /(ρ ⟨u⟩ 2 ) [71] for each simulation and compared with values from the Colebrook correlation [72] in Table II.In all cases, f D values are slightly underpredicted compared to empirical formulations, as expected due to the nature of turbulence models used in LES, which do not fully resolve wall dynamics.However, with a gradual increase in grid resolution and refinement near the wall, the results improve noticeably.The relative error in the friction coefficient, To evaluate the predictive capabilities of the subgrid models, numerical simulations of a turbulent pipe flow at a fixed turbulent Re τ = 361, are performed, comparing results against previously reported DNS findings [7].This comparison involves five different subgrid models (SMG, OEEVM, WALE, LDKM, and DSEM).As we refine the mesh resolu- tion, the results from the LES gradually approach the DNS results, as shown in Figs. 2 and 3 of Sec.IV A. The reason for this is that as ∆ decreases k sgs and ν sgs also decrease approximately as ∆ 2 .Since the objective of this sensitivity study is to examine the influence of the subgrid model, it is advisable to use the coarser meshes.The coarser meshes will emphasize the differences between the subgrid models more effectively.On finer meshes, we will not be able to detect any noticeable differences between the models.Therefore, to facilitate an effective comparison, mesh M2 has been se- lected.
Table III presents a priori and a posteriori estimates of Re τ , estimates of f D num , f D exp , and f D relErr for each subgrid model.It is noted that, across all cases, the simulated f D values are slightly underpredicted compared to targeted values.However, no discernible trend linking f D with subgrid models is evident, with negligible differences in f D values among the various models.Although additional refinement near the wall could potentially enhance results, the impact on different models remains uncertain.Nevertheless, in line with the LES philosophy, such improvement would come at the expense of increased computational requirements -an outcome undesirable in the context of this subgrid model comparison.
To assess predictive abilities in both near-wall and outer regions, we calculate mean velocity and mean velocitydefect profiles normalized in wall units.Fig. 4 (a) displays mean velocity profiles in wall units, while Fig. 4 (b) shows mean velocity-defect profiles in the outer region for different subgrid models.Analyzing Fig. 4, all subgrid models exhibit distinct behaviors, clustering based on their predictive capabilities.Near the wall, LES-WALE and LES-LDKM accurately predict mean velocity profiles, while LES-DSEM is also accurate but slightly less so.In contrast, LES-SMG and LES-OEEVM slightly underpredict near-wall mean flow but show reasonable agreement with DNS data and other models in the core flow region.
Our computational models yield slightly lower Re τ than El Khoury et al. [7].This results in slight overestimates in mean velocity profiles when normalized using inner units compared to DNS data [7].In Fig. 4 (b), showing mean velocitydefect profiles, there is a slight underestimation compared to DNS values.This occurs as we used the theoretical value of u τ as a normalization parameter instead of the actual value measured in each simulation, see Fig. 14 in Appendix A. For y/R > 0.65, all curves collapse onto DNS results.Additionally, all mean velocity profiles exhibit a log-law behavior above the buffer layer, i.e., y + ≥ 30 (see the orange solid line in Fig. 4 (a)).In our results, the logarithmic representation of turbulent profiles shows a universal behavior across subgrid models, with values 1/κ = 2.5 and B + = 5.2, proposed and accepted for various canonical flows.
RMS velocity fluctuation values are shown in Figs. 5  (a-c).LES-SMG and LES-OEEVM underpredict peak values and overpredict their locations compared to LES-WALE, LES-LDKM, LES-DSEM, and DNS.LES-WALE, LES-LDKM, and LES-DSEM accurately predict peak locations but overestimate the values.This aligns with observations in other channel flow studies where u + rms tends to be overpredicted in LES compared to DNS [68].For v + rms and w + rms , all subgrid models exhibit excellent qualitative behavior compared to DNS, as shown in Figs. 5 (b-c).However, quantitatively, all subgrid models underpredict the fluctuation values.LES-OEEVM and LES-SMG are highly diffusive near the wall and fail to predict peak values accurately.In contrast, LES-LDKM and LES-WALE behave closely to DNS but with some underprediction.Surprisingly, LES-DSEM emerges as the most accurate candidate when compared to DNS.Both profiles of velocity fluctuations (cyan data in Figs. 5 (b-c)) show good agreement with DNS due to their anisotropic nature, albeit with slight underprediction.This may be related to fixed values of coefficients c m and c e in the kinetic energy distribution equation (Eq.( 4)) in the DSEM subgrid model.While dynamically computing these coefficients could improve predictions, such enhancements are beyond the scope of this paper.
R + u ′ v ′ profiles in Fig. 5 (d) show similar behavior to DNS for LES-WALE, LES-LDKM, and LES-DSEM models.In the viscous wall region (y + < 15) and outer region (y + > 50), R + u ′ v ′ becomes linear.Quantitatively, LES-WALE and LES-LDKM reproduce DNS tendencies, including a peak at y + ≃ 30 and semi-linear behavior in the outer region.Both subgrid models accurately predict the peak location, albeit with a slight underestimation of its value.Surprisingly, LES-DSEM closely agrees with DNS compared to LES-WALE and LES-LDKM.In contrast, LES-SMG and LES-OEEVM tend to overpredict the peak location while underpredicting the peak value, consistent with rms values of velocity fluctuations in Figs. 5 (a-c).
C. Investigation of fully-developed turbulent flows for Re τ ranging from 361 to 2, 000 The analysis thus far suggests that the DSEM, WALE, and LDKM subgrid models accurately predict flow behavior.Therefore, we proceed with the WALE model for further exploration at high Re τ values.We choose not to use LDKM and DSEM, as our earlier observation suggests similar pre-  diction capabilities among WALE, LDKM, and DSEM, simplifying the computational process.This simplifies the computational process by selecting one model.Results are compared against previously reported DNS results [7,23] for four different Re τ values: 361, 550, 1,000, and 2,000.In all simulations, Re τ values are initially undefined and computed a posteriori, with Re b values taken from El Khoury et al. [7] and Pirozzoli et al. [23].To explore near-wallmodeled LES behavior under coarse grid conditions and high Re τ , we choose the coarse M2 grid as the base grid for Re τ = 361, where approximately 85% of turbulent Reynolds shear stresses and turbulent kinetic energy are resolved.Following Pope [73], this implies using a grid marginally acceptable for LES in the high Re study, expecting errors to be at least as high as for Re τ = 361 or possibly higher.

Mesh refinement based on Taylor scaling
To resolve a relatively similar amount of turbulent Reynolds shear stresses and turbulent kinetic energy for higher Re τ , we estimate the integral, Taylor, and Kolmogorov scales using conventional approximations.The integral length scale, l I , is set as the pipe diameter, i.e., l I = D.The Taylor scale, λ , an intermediate length scale in the inertial subrange, is estimated as λ = (15ν/ε)u rms [74].Here, u rms = (u ′2 + v ′2 + w ′2 ) /3 is the root mean square of the velocity fluctuations.The Kolmogorov scale, η, dissipation subrange, is estimated as η = (ν 3 /ε) 1/4 [75].Combining the expression for the Kolmogorov scale with the dissipation rate, ε ≃ k 3/2 /l I , where k = 3 2 u 2 rms , yields η/l I = C(u rms l I /ν) −3/4 = CRe −3/4 l , with C ≃ 0.86 as a constant and Re l as the integral scale Re number.The value of u rms , can be approximated using the turbulent intensity, I, as u rms = IU b , where I = 10% for high turbulence intensity applications.The grid resolution needed to resolve a relatively similar amount of the turbulent Reynolds shear stresses and the turbulent kinetic energy as at Re τ = 361 for higher Re τ values is computed from the expression of λ /l I and provided in Table IV.
Table IV presents a priori and a posteriori Re τ estimations.The chosen Re b value is deemed adequate.As expected, Re τ is underpredicted in all simulations (Table IV).However, the earlier grid resolution analysis suggests significant improvement with refinement, indicating eventual convergence  In Fig. 6 (a), mean velocity profiles from LES-WALE are compared to DNS data for various Re τ values.A generally good agreement is observed, especially in the inner layer and overlap regions, with data gradually converging to a logarithmic distribution.Fig. 6 (b) shows profiles of streamwise turbulent velocity fluctuations normalized by u 2 τ in relation to Re τ .For low Re τ , results match DNS profiles, but as Re τ increases, differences between LES-WALE and DNS become more pronounced.Specifically, LES-WALE consistently overpredicts the peak value of u + rms for each Re τ , intensifying with higher Re τ .Conversely, in the outer region, LES-WALE tends to underestimate velocity fluctuations, more noticeable with increasing Re τ .Meshes based on Taylor scaling do not capture fluctuations' overshoot in the outer region observed in DNS.Results from Berrouk et al. [46] and Ferrer et al. [47] also show discrepancies with DNS data.This could be due to our grid resolution and seems inherent in LES subgrid models, as seen in our previous grid resolution study, even for grid M3.A more detailed investigation with advanced subgrid models like DSEM is recommended, although beyond the scope of this study.

Optimizing Mesh: Directional Refinement for Improved Results
To address the mesh resolution issue related to largescale fluctuations, we revisited our meshing strategy, initially based on Taylor's scale.We conducted simulations at Re τ = 1000 for several new generated grids (G31, G32, G33, G34) with decreasing mesh spacing in all directions.Similar meshes were generated for Re τ = 361, 550, 2,000.Additional details on the grids are available in Table .V. Results for M1-M4 are already presented in Sec.IV A.
Figure 7 displays mean turbulent shear stress and rms streamwise velocities for various meshes at Re τ = 1000, along with DNS data for comparison.G32, G33, and G34 simulations show collapsed profiles for turbulent shear stress, while G3 and G31, using coarser resolution based on Taylor's scale, underestimate the results.This suggests grid convergence is achieved for ∆y + wall ∼ 3, ∆z + ∼ 15, ∆x + ∼ 30.However, for similar resolution, Fig. 7 (b) shows that rms velocity fluctuations are not collapsed.Good agreement with DNS data is only observed with G34.Comparing profiles for G3, G31, G32, and G33, it is observed that to accurately capture small-scale fluctuations, especially overshoots in the outer region, in the presence of wall models, ∆y + values can be relaxed, but ∆x + and ∆z + values play a significant role and must be strictly constrained.A trend can even be established between the peak rms velocity and the grid resolution: the former improves when the mesh resolution is respected as follows: ∆z + ≤ 10, ∆x + = 20 − 25.Moreover, the G34 mesh resolution is determined by the Kolmogorov scales, specifically the grid points utilized in LES-WALE, denoted as N LES ≈ (2.25 • Re τ ) 9/4 .It is worth noting that the factor 2.25 differs from the one employed in DNS studies, which is 3 [76].
Before doing further quantitative analysis, we show the evolution of the vortical structures for different Re τ values in Fig. 8. Vortices, identified using the λ 2 -criterion and plotted for λ 2 /λ 2 max = −0.01,appear densely along the pipe from the wall upto the wall-normal location y/R ≃ 0.5.Notably, with increasing Re τ , fluctuations and resolved vortical structures become finer and more amplified.Streaks become discernible in the near-wall cylindrical shells, arranged in a cross-stream pattern.Regardless of turbulent Re, lowspeed streaks (deep blue color, approximately the size of the pipe radius) align with large-scale ejections, while highspeed streaks (red-yellow color, of same size as low-speed streaks) correspond to substantial inflow from the core flow.Smaller streaks appear, corresponding to ejection and sweep events in the buffer layer.The flow organization on at least two length scales is apparent, with the separation increasing as Re τ values escalate.
In Fig. 9 (a), we present mean velocity profiles computed with LES-WALE for Re τ = 361, 550, 1,000, and 2,000 using new meshes M4, G21, G34, and G41.A good agreement is observed between LES-WALE and DNS, particularly in the inner layer and overlap regions, with the data gradually converging to a logarithmic distribution.This log-law is visually fitted as ⟨u⟩ + = (1/κ)log(y + ) + B, with coefficients κ = 0.387 and B = 4.53 from Pirozzoli et al. [23].The results closely aligns with the estimates from direct fitting of mean velocity profiles in pipe flows [77], yielding ⟨u⟩ + = (1/κ)log(y + ) + B, with κ = 0.391 and B = 4.34.Not much difference is obtained compared to results from coarser meshes, shown in Fig. 6 (a).
Figure 10 (a) provides insights into the turbulent shear stress distribution for varying Re τ values.As Re τ increases, the shear stress profiles flatten, with the peak value rising toward unity, and its position shifting further from the wall.This well-established behavior [69] is characterized by the asymptotic relationship that accounts for the position of the turbulent shear stress peak, y + p [78], The precision of Eq. ( 8), with κ = 0.387, is notably accurate starting from Re τ ≃ 1,000 onwards, indicating the onset TABLE V: Parameters of the new grid used for the grid convergence study; mesh spacings are given in wall units; ∆x + , ∆z + : mesh spacings in the streamwise and azimuthal directions; ∆y + wall , ∆y + center : mesh spacings in the wall normal direction at the wall and at the center of the channel.  of a near-logarithmic layer [79].In our simulations, the predicted peak position, y + p , closely aligns with expected values for Re τ = 1,000 and 2,000, with deviations of approximately 5%.However, for Re τ < 1,000, the predictions deviate significantly from the asymptotic relation, showing deviations of roughly 27%.
Based on the proposed criteria we have further looked upon the velocity fluctuations results for Re τ = 361, 550, 1000 and 2,000 in Fig. 10 (b-d).Compared with the results obtained using Taylor scaling (Fig. 6 (b)), the present results matches with the DNS, as shown in Fig. 10.Slight discrepancies between DNS results and LES at Re τ = 2,000 can still be observed; nevertheless, both the primary peaks and the secondary overshoot region are well captured in comparison to existing literature [46,47].Unlike previous results, the mesh requirement for Re τ = 2000 is ∼ 220M cells and increases tremendously with Re τ .Therefore, we restrict this study to Re τ = 2000.

Single point statistics
To quantitatively depict near-wall regions at different Re τ , we utilize quadrant analysis for streamwise and wall-normal velocity fluctuations.Quadrant analysis divides the (u ′ , v ′ ) plane into four quadrants (Q 1 , Q 2 , Q 3 , Q 4 ) based on sign combinations of velocity fluctuations.Specifically, Q 1 corresponds to positive values of both u ′ and v ′ , Q 2 represents u ′ < 0, v ′ > 0, Q 3 indicates u ′ < 0, v ′ < 0, and Q 4 corresponds to u ′ > 0, v ′ < 0. Ejection-like events occur in Q 2 , while sweep-like events are in Q 4 .Quadrant analysis provides insights into the dynamics of near-wall streaks, initially observed by Kline et al. [80] and further explored by Corino & Brodkey [81] and Wallace et al. [82].Willmarth & Lu [83] computed quadrant contributions to the total shear stress, highlighting the importance of ejection and sweep events in the production of turbulent Reynolds shear stress.The shear stress plays a crucial role in turbulence dynamics, appearing in the LES equation and the turbulent kinetic energy equation's production term, P i j , multiplied by the mean shear stress.
Fig. 11 shows contour plots of the joint probability density function (JPDF) of velocity fluctuations, P(u ′ ,v ′ ), normalized by their local rms values, σ u ′ and σ v ′ .These plots are at y + = 3 √ Re τ for four Re τ values.The JPDF iso-contours have an elliptical shape, mainly spanning quadrants Q 2 and Q 4 , with the maximum value in Q 4 , indicating a higher frequency of {u ′ v ′ (t) < 0} events, particularly with more extreme values in Q 2 (ejections).Normalized PDFs are also displayed on the top (u ′ ) and right (v ′ ) corners in each plot.Fig. 12 presents PDFs of streamwise and wall-normal velocity fluctuations, P(u ′ ) and P(v ′ ), at various distances from the wall for different Re τ .These PDFs concern fluctuations relative to the mean and are normalized by the root-mean-square value.In the near-wall region (Fig. 12 (a)), the distribution of streamwise velocity fluctuations, u ′ , deviates significantly from a Gaussian shape, indicating intermittency with more intense and intermittent high-speed streaks close to the wall.This behavior contributes to the positive skewness observed in experiments [84] and DNS [85].The PDFs of wall-normal fluctuations, v ′ , (Fig. 12 (c)), are weakly asymmetric, favoring descents over ascents.In the outer region of the pipe (y + = 150), the PDF shapes approach a Gaussian distribution (Fig. 12 The state of anisotropy can then be characterized with the two variables η and ξ defined as and Reynolds stress tensor states are confined within the Lumley triangle, defined in the (ξ ,η) plane, with different realizable states characterized by the invariants of the normalized anisotropy tensor, b i j .Details on these states can be found in [86].Anisotropy maps in Fig. 13 depict variations for four Re τ values across y + .These maps utilize all somain cells, and post-processed values of η and ζ are represented with a color map of y + , indicating the likelihood of specific turbulent states.Near the wall, turbulence is highly twodimensional (dominated by u ′ u ′ and w ′ w ′ fluctuations), transitioning toward a one-dimensional, one-component state as u ′ u ′ peaks around y + ≃ 8, as highlighted in the inset of Fig. 13.As y + increases, data points move toward the origin, aligning along the axisymmetric state, signifying an approach to isotropy.These findings align well with Pope's analysis [75] of DNS data from Kim et al. [66], where twocomponent turbulence was identified near the wall (y + < 5), transitioning to axisymmetric expansion in the channel's remaining portion.The anisotropy maps exhibit a universal behavior across different Re τ values.

V. CONCLUSIONS
The present work assesses the prediction capabilities of wall-modeled Large Eddy Simulations (LES) at moderateto-high Reynolds numbers and in the limit of coarse meshes, in the OpenFOAM framework.Initially, results for a turbulent pipe at Re τ = 361 were compared against DNS data available in the literature.The simulation results obtained using the WALE subgrid model showed excellent qualitative and quantitative prediction accuracy when compared to DNS data; nonetheless, with slight underprediction of all turbulence characteristics.We further assessed the accuracy of the WALE subgrid model by performing a grid resolution study on four different grids (M1, M2, M3, and M4, where M1 is the coarsest grid, and M4 is the finest grid).The results showed a gradual increase in accuracy of zeroth-, first-, and second-order critical turbulence statistics with grid resolution.However, as a trade-off between numerical accuracy and computational wall clock time, we decided to proceed with a rather coarse M2 grid for further detailed investigations.
Based on the grid resolution study, we performed LES for five different subgrid models: the Wall-Adapting Local Eddy viscosity model (WALE), the Smagorinsky model (SMG), the One-Equation Eddy Viscosity Model (OEEVM), the Localized Dynamic Kinetic energy Model (LDKM), and the Differential Stress Equation Model (DSEM) on grid M2.The results obtained demonstrated that WALE, LDKM, and DSEM were able to reproduce the main and most critical characteristics of turbulent flows, featuring regions of high strain rates and large-scale vortices, with higher accuracy than OEEVM and SMG, both quantitatively and qualitatively, in comparison to DNS data from the literature.Observing the flow statistics, it could be evidenced that the estimates of the mean flow profiles and turbulent Reynolds shear stresses were consistent for the three subgrid models further investigation of grid convergence using Kolmogorov scaling, it was found that to accurately capture second-order statistics, particularly secondary overshoots in fluctuations, it is necessary to utilize ∆x + ≤ 20 − 25 and ∆z + ≤ 10.This leads to an actual number of grid points required for wellresolved LES of pipe flow, denoted as N LES , at approximately (2.25 • Re τ ) 9/4 .The coefficient differs from DNS, which is 3.While the computational savings compared to DNS are extremely significant, they amount to at least 30%.
To conclude, our findings suggest that the scaling strategy we have explored in this paper is not effective, and we advise against others attempting to replicate this approach.Rather, we recommend a different methodology that involves a direct measurement of the resolution using appropriate criteria such as modeled kinetic energy and grid variations that ensure comparable resolutions across different Re τ .Although this approach can be time-consuming and challenging, it is essential for achieving relatively similar resolutions across different Reynolds numbers.Therefore, we suggest either the use of grid spacing constraint proposed above or the use of simulation methods that replace the filter width with a modeled length scale and include the modeled kinetic en-ergy for resolution calculations.By adopting these methods, one can ensure a more accurate and reliable resolution for the simulations in various physical and engineering applications.TABLE VI: Correlation times and lengths in the streamwise direction of axial, u ′ , wall-normal, v ′ , and azimuthal, w ′ , velocity fluctuation components.The correlation times are computed as t i ′ = T max 0 C i ′ (T )dT , while the correlation lengths are deduced from the zero-crossing lengths of the corresponding autocorrelation functions, C i ′ (X).
100 %, tends towards zero.Once again, the measurement of f D shows that M3 and M4 exhibit similar behavior and outperform other grids in terms of accuracy.B. Different subgrid model behavior for Re τ = 361

FIG. 2 :FIG. 3 :
FIG. 2: Comparison of LES-WALE and DNS data for mean velocity profiles (a) in inner variables and (b) in velocity-defect form.U + cl is the normalized centerline velocity.DNS data at Re τ = 361 is shown in dashed black line.The dashed-dotted pink and solid orange lines in Fig. 2 (a) represent the universal behavior of the turbulent velocity profile.

FIG. 4 :
FIG. 4: (a) Mean velocity profiles and (b) mean velocity-defect profiles for Re τ = 361 and five different subgrid models, normalized by u 0 τ which is a priori estimate of the friction velocity.U + cl is the normalized centerline velocity.The dashed black line presents DNS data from El Khoury et al. [7] and the dashed pink and solid orange lines in Fig. 4 (a) represent the universal behavior of the turbulent velocity profile.

FIG. 5 :
FIG. 5: Comparison of (a) axial, u + rms , (b) wall-normal, v + rms , and (c) azimuthal, w + rms , turbulence intensities; (d) turbulent Reynolds shear stress, R + u ′ v ′ , with DNS data (dashed lines) as functions of y + , at Re τ = 361, for different subgrid models.The thin diagonal line represents the total shear stress, with the difference between this line and the data indicating the viscous shear stress attributed to the mean flow.

u= 2 FIG. 6 :
FIG. 6: Comparison of (a) ⟨u⟩ + and (b) u + rms , as functions of y + for different values of Re τ and the WALE subgrid model.u + rms profiles are plotted with an offset of ∆ ⟨.⟩ + = 1 for clarity.In all figures, the DNS data are represented in dashed black lines.

FIG. 7 :
FIG. 7: (a) Turbulent shear stresses (R + u ′ v ′ ) and (b) axial turbulent velocity fluctuations (u + rms ) on a semi-logarithmic scale, are shown as a function of y + at Re τ = 1,000 for the WALE subgrid model.In all subfigures, the DNS data are represented in dashed black lines.
Fig. 11 shows contour plots of the joint probability density function (JPDF) of velocity fluctuations, P(u ′ ,v ′ ), normalized by their local rms values, σ u ′ and σ v ′ .These plots are at y + = 3 √ Re τ for four Re τ values.The JPDF iso-contours have an elliptical shape, mainly spanning quadrants Q 2 and Q 4 , with the maximum value in Q 4 , indicating a higher frequency of {u ′ v ′ (t) < 0} events, particularly with more extreme values in Q 2 (ejections).Normalized PDFs are also displayed on the top (u ′ ) and right (v ′ ) corners in each plot.Fig.12presents PDFs of streamwise and wall-normal velocity fluctuations, P(u ′ ) and P(v ′ ), at various distances from the wall for different Re τ .These PDFs concern fluctuations relative to the mean and are normalized by the root-mean-square value.In the near-wall region (Fig.12 (a)), the distribution of streamwise velocity fluctuations, u ′ , deviates significantly from a Gaussian shape, indicating intermittency with more intense and intermittent high-speed streaks close to the wall.This behavior contributes to the positive skewness observed in experiments[84] and DNS[85].The PDFs of wall-normal fluctuations, v ′ , (Fig.12 (c)), are weakly asymmetric, favoring descents over ascents.In the outer region of the pipe (y + = 150), the PDF shapes approach a Gaussian distribution (Fig.12(b), (d)), indicating that Re τ has minimal impact on the overall dominance of ejection-and sweep-like events.

( a )
FIG. 13: Anisotropy-invariant mapping of turbulence from present LES-WALE data at (a) Re τ = 361, (b) Re τ = 550, (c) Re τ = 1,000, and (d) Re τ = 2,000.Data points for each subgrid model are based on all cells in the domain and are colored with normalized wall distance values, y + .The color map varies from pink to yellow.Insets in all subfigures show the same zoom of the one-component turbulence area.

TABLE I :
Mesh properties.Mesh Ncells along x-direction Total size ∆y + Maximum skewness Maximum aspect ratio

TABLE II :
Computed global flow parameters for different grid resolutions.Re τ = u 0 τ D/(2ν), and Re ap τ = u τ D/(2ν) are the a priori and a posteriori estimates of Re τ , respectively.u 0τ and u τ are the a priori and a posteriori estimates of the friction velocity, respectively.fD num and f D exp are numerical estimation and target value of the Darcy friction coefficient, respectively.fD relErr the relative error in the numerical computation of f D .

TABLE III :
Computed global flow parameters for different subgrid models.Re τ = u 0 τ D/(2ν), and Re ap τ = u τ D/(2ν) are the a priori and a posteriori estimates of Re τ , respectively.u 0τ and u τ are the a priori and a posteriori estimates of the friction velocity, respectively.fD num and f D exp are numerical estimation and target value of the Darcy friction coefficient, respectively.fD relErr the relative error in the numerical computation of f D .

TABLE IV :
Computed global flow parameters for different turbulent Re τ .Re τ = u 0 τ D/(2ν), and Re ap τ = u τ D/(2ν) are the a priori and a posteriori estimates of Re τ , respectively.u 0 τ and u τ are the a priori and a posteriori estimates of the friction velocity, respectively.f D num and f D exp are numerical estimation and target value of the Darcy friction coefficient, respectively.f D relErr the relative error in the numerical computation of f D .Re τ .Numerical estimations of f D num for various turbulent Re τ , along with f D exp values, are in Table IV.Deviations of approximately 5-20% from expected results are observed.
Re τ T u ′ /t T v ′ /t T w ′ /t l u ′ /L l v ′ /L l w ′ /L