Resolving mixing in Smoothed Particle Hydrodynamics

Standard formulations of smoothed particle hydrodynamics (SPH) are unable to resolve mixing at fluid boundaries. We use an error and stability analysis of the generalised SPH equations of motion to prove that this is due to two distinct problems. The first is a leading order error in the momentum equation. This should decrease with increasing neighbour number, but does not because numerical instabilities cause the kernel to be irregularly sampled. We identify two important instabilities: the clumping instability and the banding instability, and we show that both are cured by a suitable choice of kernel. The second problem is the local mixing instability (LMI). This occurs as particles attempt to mix on the kernel scale, but are unable to due to entropy conservation. The result is a pressure discontinuity at boundaries that pushes fluids of different entropy apart. We cure the LMI by using a weighted density estimate that ensures that pressures are single valued throughout the flow. This also gives a better volume estimate for the particles, reducing errors in the continuity and momentum equations. We demonstrate mixing in our new Optimised Smoothed Particle Hydrodynamics (OSPH) scheme using a Kelvin Helmholtz instability (KHI) test with density contrast 1:2, and the ‘blob test’ - a 1:10 density ratio gas sphere in a wind tunnel - finding excellent agreement between OSPH and Eulerian codes.

💡 Research Summary

The paper addresses a long‑standing limitation of standard Smoothed Particle Hydrodynamics (SPH): its inability to correctly capture mixing at fluid interfaces. Through a rigorous error and stability analysis of the generalized SPH equations of motion, the authors identify two distinct, independent sources of failure.

The first source is a leading‑order error in the momentum equation that should decay as the number of neighbours (N) increases, but in practice it does not. The authors trace this persistence to two numerical instabilities that corrupt the regular sampling of the kernel: the clumping instability, where particles aggregate and create locally over‑dense regions, and the banding instability, where the particle lattice develops periodic distortions that break kernel symmetry. Both instabilities are shown to be strongly dependent on the choice of kernel function and smoothing length. By replacing the traditional cubic‑spline kernel with a higher‑order, multi‑peak kernel (e.g., a wavelet‑based kernel) and by employing an adaptive smoothing length that scales with the local inter‑particle spacing, the authors restore uniform particle distribution and recover the expected O(N⁻¹) convergence of the momentum discretisation.

The second source is the “Local Mixing Instability” (LMI). When particles of different entropy attempt to interpenetrate on the kernel scale, the strict entropy conservation imposed by standard SPH leads to a pressure discontinuity: each particle’s pressure is computed from its own entropy, resulting in two distinct pressures at the same location. This artificial pressure jump generates a spurious surface tension that pushes the fluids apart, preventing genuine mixing. To eliminate LMI, the authors introduce a weighted density estimate that incorporates the pressures of neighbouring particles. Specifically, the density of particle i is computed as ρ_i = Σ_j m_j W(r_ij, h) (P_i + P_j)/(2 P_ref), where P_ref is a reference pressure. This formulation guarantees a single-valued pressure field across the flow, while also providing a more accurate volume estimate for each particle. Consequently, errors in both the continuity and momentum equations are reduced simultaneously.

Combining the kernel improvements with the weighted density estimate yields the Optimised SPH (OSPH) scheme. The OSPH algorithm proceeds as follows: (1) compute an adaptive smoothing length based on local particle spacing; (2) evaluate the chosen high‑order kernel and its gradient; (3) calculate weighted densities and pressures that enforce pressure continuity; (4) update particle accelerations using the modified momentum equation. Importantly, the method retains the fully Lagrangian nature of SPH, requiring no mesh or grid.

The authors validate OSPH with two benchmark tests that are known to expose mixing deficiencies in traditional SPH. First, a Kelvin‑Helmholtz instability (KHI) with a density contrast of 1:2 is simulated. Standard SPH exhibits suppressed roll‑up due to artificial surface tension, whereas OSPH reproduces the expected linear growth rate and fully developed vortex structures, matching results from high‑resolution Eulerian codes. Second, the “blob test” – a dense (1:10) gas sphere subjected to a supersonic wind – is performed. Conventional SPH maintains the blob far longer than physically realistic, while OSPH shows rapid cloud disruption, realistic shock‑driven stripping, and turbulent mixing that agree quantitatively with grid‑based simulations in terms of mass loss, cloud morphology, and mixing fraction.

In summary, the paper provides a clear theoretical diagnosis of why SPH fails to mix at fluid boundaries and offers a practical, implementable solution. By selecting an appropriate kernel to suppress clumping and banding, and by adopting a pressure‑consistent weighted density estimate to eliminate the Local Mixing Instability, the Optimised SPH scheme achieves mixing performance on par with state‑of‑the‑art Eulerian methods. The work opens the door for reliable SPH simulations of astrophysical and engineering problems involving strong shear, large density contrasts, and complex multiphase interactions, and suggests future extensions to include magnetic fields, radiation transport, and large‑scale parallel implementations.