Smoothed Particle Hydrodynamics Simulations

Updated 13 January 2026

Smoothed Particle Hydrodynamics is a mesh-free, fully Lagrangian method that discretizes continuum equations using particle interactions.
It employs advanced numerical techniques, including adaptive kernels and artificial viscosity, to capture shocks and complex boundary dynamics.
Recent innovations, such as multiphase modeling, multiresolution strategies, and GPU acceleration, enhance simulation accuracy and performance.

Smoothed Particle Hydrodynamics (SPH) is a fully Lagrangian, meshfree method for simulating fluid and, more generally, continuum mechanics. In SPH, the continuous governing equations (e.g., the Navier–Stokes equations) are discretized on a set of particles, each carrying material properties such as mass, velocity, and internal energy. This approach facilitates large-deformation dynamics, complex boundary geometries, and multi-physics applications, at the cost of distinctive numerical challenges. The method is implemented in diverse computational frameworks for research in computational fluid dynamics, astrophysics, planetary science, engineering, and computational graphics.

1. Mathematical Foundation and Discretization Principles

SPH is derived from an integral interpolant representation of fields. For a scalar field $f(\mathbf{r})$ , the SPH kernel approximation is

$f(\mathbf{r}) \approx \int f(\mathbf{r}') W(\mathbf{r}-\mathbf{r}',h) d\mathbf{r}'$

where %%%%1%%%% is a compact-support, bell-shaped kernel of width $h$ . The particle-based discretization yields the standard interpolant

$\langle f \rangle_a = \sum_{b=1}^N \frac{m_b}{\rho_b} f_b W_{ab}$

and its gradient as

$\langle \nabla f \rangle_a = \sum_{b=1}^N \frac{m_b}{\rho_b} f_b \nabla_a W_{ab}$

with $m_b$ the particle mass, $\rho_b$ density, and $W_{ab}$ the kernel evaluated between particles $a$ and $b$ (Sutti, 2022). Conservation laws in the continuous form—mass, momentum, and energy—are recast into SPH sums. Density may be assigned either by direct summation (which exactly conserves mass) or by integrating the continuity equation.

Common choices for the smoothing kernel include the cubic spline and Wendland families. The standard cubic spline is

$W(q) = \frac{\sigma_d}{h^d} \begin{cases} 1 - \frac{3}{2}q^2 + \frac{3}{4}q^3, & 0 \le q < 1 \ \frac{1}{4}(2-q)^3, & 1 \le q < 2 \ 0 & q \ge 2 \end{cases}$

where $q = \|\mathbf{x}_a - \mathbf{x}_b\|/h$ and $\sigma_d$ normalizes the kernel by dimension (Sutti, 2022).

Artificial viscosity, most commonly the Monaghan (1992) form, is added to stabilize shocks and dissipate post-discontinuity oscillations. It introduces parameters $\alpha, \beta$ (typically $\mathcal{O}(1)$ ) and an inter-particle velocity projection $\phi_{ab}$ . Standard viscosity discretization and advanced limiters (e.g., Cullen–Dehnen) are also in use (Booth et al., 2015, Hu et al., 2014).

2. Algorithmic Workflow, Boundary Conditions, and Time Integration

A prototypical SPH simulation loop comprises:

Neighbor search: For each particle, all neighbors within $2h$ are identified, typically using a cell-linked list for $\mathcal{O}(N)$ complexity (Sutti, 2022).
Density and property update: Density is calculated either by summation or integrating the continuity form.
Pressure calculation: Pressure is closed via an equation of state (EOS), e.g., polytropic law for gases ( $p_a = (\gamma-1)\rho_a e_a$ ) or Tait EOS for weakly compressible liquids.
Force accumulation: The momentum equation is discretized, summing pressure, viscosity, external forces, and (if needed) artificial terms (e.g., for shocks or interface stabilization).
Position and velocity update: An explicit predictor–corrector or leapfrog time integration advances positions and velocities, constrained by a Courant–Friedrichs–Lewy (CFL) condition (Sutti, 2022).

Boundary conditions are handled via several strategies:

Ghost (virtual) particles—static or dynamically mirrored to enforce no-penetration or no-slip (Sutti, 2022).
Repulsive forces (Lennard–Jones or radial) for impenetrable boundaries or free surfaces.
Normal boundary force and semi-analytical corrections for accurate wall treatments (Sutti, 2022).

Specialized algorithms, such as Moving Least Squares (MLS) detection and curvature estimation, underpin advanced free-surface and interfacial modeling (Terissa et al., 2013). The implementation structure in codes such as SPHM is modular, with separate routines for each physical term and integration step.

3. Performance, Scalability, and High-Performance Computing

Efficiency of SPH codes strongly depends on neighbor search scalability, data structure layout, and parallelization strategy. For moderate problem sizes ( $N\lesssim10^5$ ), uniform cell-linked lists suffice (Sutti, 2022); for large or non-uniform cases, octrees, Barnes–Hut, or Fast Multipole (FMM) trees are deployed for both gravity and neighbor search (Meier et al., 14 Nov 2025, Guerrera et al., 2018).

GPU architectures offer substantial acceleration due to the particle-local and data-parallel structure of SPH. Modern GPU-based codes employ:

Structure-of-Arrays data layouts,
Parallel cell sorting and neighbor search via radix sort and prefix-sum,
Overlapping compute and communication,
Dedicated CUDA/OpenCL/OpenACC kernels for force and update steps. Benchmarks indicate real-time or faster performance for problems up to millions of particles on consumer GPUs, with 10–50× speedup versus CPU-only execution (Coman et al., 2019).

For exascale readiness, frameworks like SPH-EXA integrate MPI, OpenMP/threading, and accelerator offload. They utilize hybrid domain decomposition (space-filling curve, ORB), dynamic load balancing, and multi-level checkpoint/restart (Guerrera et al., 2018, Cavelan et al., 2020). At scale, bottlenecks shift from neighbor summation to global synchronizations, tree-building, and I/O; optimized designs minimize global barriers and memory footprint.

4. Applications and Benchmark Validation

SPH is validated against canonical problems:

The 1D Sod shock tube test: SPH accurately captures shock propagation, contact discontinuity, and rarefaction, with L1 errors $\lesssim10\%$ at discontinuities, $<1\%$ in smooth regions (Sutti, 2022).
Shear-driven cavity and Taylor–Couette flows: SPH recovers steady-state vortex patterns; velocity profile errors are typically $<5\%$ versus CFD/FEM benchmarks (Sutti, 2022, Kincl et al., 2022).
Free-surface and dam-break flows: Accurate interface dynamics and impact on obstacles captured, with kernel-enhanced interface tracking via gradient renormalization (Prayogo et al., 2013).
Kelvin–Helmholtz and Rayleigh–Taylor instabilities: Accurate roll-up and mixing are obtained with proper kernel, shifting and conduction strategies (Alvarado-Rodriguez et al., 2021, Chen, 2020, Hu et al., 2014).
Metal solidification: Combined flow, heat conduction, and phase change processes are captured via enthalpy formulation and high-viscosity solid phase models (S et al., 2013).

Performance benchmarks demonstrate linear or modestly superlinear scaling in strong/weak-scaling (e.g., PKDGRAV3, SPH-EXA), with high efficiency (>60–80%) to thousands of CPU cores and hundreds of GPUs (Cavelan et al., 2020, Meier et al., 14 Nov 2025).

5. Advanced Methods: Multiphysics, Multiphase, and Multiresolution SPH

Recent developments extend SPH to multiphase and multiphysics flows, solid–fluid coupling, and adaptive resolution:

Multiphase SPH incorporates variable smoothing length, shifting algorithms, and pairwise or color-gradient surface tension models. This yields sharply resolved, stable interfaces even for density ratios $\gtrsim 1000$ (Chen, 2020, Alvarado-Rodriguez et al., 2021).
Two-fluid gas–dust SPH leverages integral-gradient estimators, high neighbor counts (Wendland kernels, $N_\text{ngb} \gtrsim 400$ ), and semi-implicit drag integration to accurately model strong and weak coupling regimes (Booth et al., 2015).
Unified SPH for fluids and solids employs the SHTC (Symmetric Hyperbolic Thermodynamically Compatible) PDE framework, introducing distortion tensors and explicit relaxation for viscoelastic and plastic behavior (Kincl et al., 2022).
Multiresolution SPH uses second-order consistent operators, dynamic particle splitting/merging, and particle-shifting to maintain accuracy and regularity across resolution interfaces, with CPU time reductions up to 40–60% in 2D domains (Hu et al., 2017).

Godunov SPH variants (GSPH) replace artificial viscosity with Riemann-solver–based fluxes between particle pairs, greatly improving shock resolution and capturing mixing at contact discontinuities (Murante et al., 2011).

6. Best Practices, Limitations, and Directions for Improvement

Practical guidance for SPH simulation setup includes:

Particle spacing $\Delta x$ and smoothing length $h$ control resolution and stability ( $h \approx 1.2\Delta x$ – $1.5\Delta x$ ) (Sutti, 2022).
Time step $\Delta t$ must satisfy all stability constraints (CFL, viscous, and force-based), using safety factors $<0.3$ (Coman et al., 2019).
Choice of kernel impacts accuracy and stability: Wendland kernels are preferred to suppress tensile instability with large neighbor numbers (Booth et al., 2015, Chen, 2020, Meier et al., 14 Nov 2025).
Artificial viscosity should be tuned to the presence of shocks; switches (e.g., Cullen–Dehnen) are essential to avoid smearing in shear flows (Booth et al., 2015, Hu et al., 2014).
Initial conditions: Weighted Voronoi Tessellation (WVT) and similar methods enable low-noise particle distributions matching arbitrary density fields, achieving interpolation errors $\lesssim1\%$ in both uniform and adaptive settings (Diehl et al., 2012).
Visualization: Kernel-based rendering (e.g., SPLASH) is necessary for accurate interpretation of SPH data; scatter plots are insufficient (0709.0832).

Limitations include: lack of inherent turbulence modeling and explicit free-surface tension (unless extended), single-threaded performance in simple codes, boundary deficiency for mass summation at boundaries, and the need for calibration of artificial terms. Modern research targets improved wall treatments, adaptive time stepping, robust multiphase and solid interface tracking, and communication-minimized exascale parallelism (Sutti, 2022, Cavelan et al., 2020, Meier et al., 14 Nov 2025).

References:

(Sutti, 2022) SPHM: a MATLAB package for Smoothed Particle Hydrodynamics simulations
(Coman et al., 2019) Hydrodynamic Simulations using GPGPU Architectures
(Guerrera et al., 2018) Towards a Mini-App for Smoothed Particle Hydrodynamics at Exascale
(Terissa et al., 2013) Three-Dimensional Smoothed Particle Hydrodynamics Simulation for Liquid Droplet with Surface Tension
(Kincl et al., 2022) Unified description of fluids and solids in Smoothed Particle Hydrodynamics
(S et al., 2013) Three-Dimensional Smoothed Particle Hydrodynamics Simulation for Liquid Metal Solidification Process
(Chen, 2020) Meshfree simulation of multiphase flows with SPH family methods
(Booth et al., 2015) Smoothed particle hydrodynamics simulations of gas and dust mixtures
(Cavelan et al., 2020) A Smoothed Particle Hydrodynamics Mini-App for Exascale
(Diehl et al., 2012) Generating Optimal Initial Conditions for Smoothed Particle Hydrodynamics Simulations
(Koschier et al., 2020) Smoothed Particle Hydrodynamics Techniques for the Physics Based Simulation of Fluids and Solids
(Prayogo et al., 2013) Three-Dimensional Smoothed Particle Hydrodynamics Method for Simulating Free Surface Flows
(Alvarado-Rodriguez et al., 2021) Multiphase flows simulation with the Smoothed Particle Hydrodynamics method
(Murante et al., 2011) Hydrodynamic simulations with the Godunov SPH
(Hu et al., 2017) A Consistent Multi-Resolution Smoothed Particle Hydrodynamics Method
(Meier et al., 14 Nov 2025) Smoothed Particle Hydrodynamics in pkdgrav3 for Shock Physics Simulations I: Hydrodynamics
(0709.0832) SPLASH: An interactive visualisation tool for Smoothed Particle Hydrodynamics simulations
(Hu et al., 2014) SPHGal: Smoothed Particle Hydrodynamics with improved accuracy for Galaxy simulations
(Borrow et al., 2020) Sphenix: Smoothed Particle Hydrodynamics for the next generation of galaxy formation simulations