  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 \leq Re_{\tau} \leq 2,000$), corresponding to
bulk velocity-based Reynolds numbers of $11,700 \leq Re_{b} \leq 82,500$, using
wall-modeled Large Eddy Simulations (LES) in OpenFOAM. A grid convergence study
is performed for $Re_{\tau} = 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_{\tau}$ value rises, it becomes crucial for a grid
to adhere to the condition of $\Delta x^{+} \leq 20 - 25$ and $\Delta z^{+}
\leq 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.
