Resolving high Reynolds numbers in SPH simulations of subsonic turbulence

Accounting for the Reynolds number is critical in numerical simulations of turbulence, particularly for subsonic flow. For Smoothed Particle Hydrodynamics (SPH) with constant artificial viscosity coefficient alpha, it is shown that the effective Reynolds number in the absence of explicit physical viscosity terms scales linearly with the Mach number - compared to mesh schemes, where the effective Reynolds number is largely independent of the flow velocity. As a result, SPH simulations with alpha=1 will have low Reynolds numbers in the subsonic regime compared to mesh codes, which may be insufficient to resolve turbulent flow. This explains the failure of Bauer and Springel (2011, arXiv:1109.4413v1) to find agreement between the moving-mesh code AREPO and the GADGET SPH code on simulations of driven, subsonic (v ~ 0.3 c_s) turbulence appropriate to the intergalactic/intracluster medium, where it was alleged that SPH is somehow fundamentally incapable of producing a Kolmogorov-like turbulent cascade. We show that turbulent flow with a Kolmogorov spectrum can be easily recovered by employing standard methods for reducing alpha away from shocks.

💡 Research Summary

The paper addresses a critical issue in numerical simulations of subsonic turbulence: the ability of Smoothed Particle Hydrodynamics (SPH) to reproduce a Kolmogorov‑5/3 inertial cascade. The authors begin by recalling that turbulence is characterised by two dimensionless numbers – the Mach number (M = V/cₛ) and the Reynolds number (Re = VL/ν). While the Mach number determines the relative importance of inertial versus pressure forces, the Reynolds number controls the balance between inertial forces and viscous dissipation, and therefore dictates whether a flow becomes turbulent and how large the inertial range can be.

In SPH, dissipation is introduced only through the artificial viscosity term. By comparing the SPH formulation with the Navier–Stokes equations, the authors show that the effective shear viscosity is ν ≈ (1/10) α cₛ h (and bulk viscosity ζ ≈ (1/6) α cₛ h), where α is the artificial‑viscosity coefficient and h the smoothing length. For low‑Mach flows the signal speed v_sig≈cₛ, so the Reynolds number can be expressed simply as

 Re = 10 α M (L/h).

Thus, at fixed Mach number the Reynolds number depends only on α and the numerical resolution L/h. Using the parameters of Bauer & Springel (2011) – a subsonic drive with M≈0.3 and particle numbers 64³, 128³, 256³ – and assuming a constant α = 1, the resulting Reynolds numbers are only ≈150, 300 and 600 respectively. These values are well below the ≈10³ threshold required for a fully developed turbulent cascade, explaining why their SPH runs failed to show a Kolmogorov spectrum. By contrast, mesh‑based codes (e.g., AREPO) have numerical dissipation that scales with the flow velocity, making their effective Reynolds numbers roughly independent of Mach number; consequently they achieve much higher Re at the same resolution.

The key remedy proposed is to reduce the artificial viscosity away from shocks. The classic Morris & Monaghan (1997) switch makes α a time‑dependent quantity that rises in convergent flows (∇·v < 0) and decays to a floor value (often α_min ≈ 0). Implementing this switch with α_min = 0 and a typical post‑shock peak α≈0.1 raises the effective Reynolds numbers by roughly an order of magnitude: Re≈1500, 3000 and 6000 for the three resolutions. The authors performed new simulations with the PHANTOM SPH code, using the same driving scheme as Bauer & Springel, and obtained power spectra that clearly exhibit a k⁻⁵/³ inertial range at the highest resolution (256³ particles, Re≈6000). The spectra agree well with those from AREPO, demonstrating that SPH can indeed reproduce Kolmogorov turbulence once the viscosity is properly controlled.

The paper also discusses more advanced viscosity switches, notably Cullen & Dehnen (2010), which combine shock detection based on higher‑order derivatives with rapid decay of α in smooth regions. These switches can reduce α to ≪0.01 in quiescent flow, potentially pushing SPH Reynolds numbers into the 10⁵–10⁶ regime, comparable to the highest‑resolution mesh calculations.

In the discussion, the authors refute the claim by Bauer & Springel that SPH’s gradient errors are responsible for the missing cascade. Their own results, using identical kernel functions and neighbour numbers, show that the gradient estimator is not the limiting factor; rather, insufficient Reynolds numbers are. They also caution against arbitrarily setting α to very low fixed values (e.g., α = 0.1 or α = 0) because artificial viscosity is still needed to capture physical shock dissipation, especially at Mach ≈ 0.3 where non‑linear steepening is non‑negligible.

Finally, the paper places the achievable numerical Reynolds numbers in astrophysical context. For the cold interstellar medium, physical Re≈10⁵–10⁷, which is within reach of current SPH resolutions if viscosity is minimised. For the intracluster/intragalactic medium, estimates range from Re≈50 (laminar) up to >10¹⁰ when magnetic fields suppress viscosity. While the latter extreme remains out of reach, the authors argue that SPH, equipped with modern viscosity switches, can attain sufficiently high Re to study the transition to turbulence and the resulting cascade in realistic astrophysical settings.

In summary, the study demonstrates that SPH is not fundamentally incapable of modelling subsonic turbulence. The apparent failure in earlier work stemmed from using a constant, large artificial viscosity, which kept the effective Reynolds number too low. By employing time‑dependent viscosity switches, SPH can achieve Reynolds numbers comparable to mesh codes and recover the expected Kolmogorov spectrum, thereby restoring confidence in SPH for a wide range of astrophysical turbulence problems.