Inside SPH: The Algorithm That Makes Our Fluid Simulation Tick

When you drag your cursor through our fluid simulation and watch liquid ripple, slosh, and splash, you're seeing Smoothed Particle Hydrodynamics (SPH) running in your browser. SPH was invented in 1977 independently by Lucy and by Gingold and Monaghan for simulating astrophysical flows. Today it underpins real-time fluid effects in games, visual effects pipelines, and engineering simulations. Here's the mathematics behind it.

The Core Idea: Particles as Fluid Samples

Classical computational fluid dynamics solves the Navier-Stokes equations on a fixed grid. Every cell of the grid has a density, velocity, and pressure, and the equations describe how fluid flows between adjacent cells. This works well for flows that stay within a fixed domain — airflow around a wing, for instance — but becomes complex when the fluid surface itself is moving: waves, splashes, drops.

SPH takes a different approach. Instead of a fixed grid, the fluid is represented by a large number of particles, each carrying mass, position, velocity, and other quantities. The continuous fluid field is reconstructed from the particles at any point in space by weighted averaging over nearby particles. The particles move with the fluid, automatically handling free surfaces — where the fluid ends and air begins — without any special treatment.

This Lagrangian (particle-following) approach is the key advantage of SPH over Eulerian (grid-fixed) methods for free-surface flows: the surface is tracked implicitly by where the particles are.

Kernel Functions: The Mathematics of Smoothing

The "smoothed" in SPH comes from the kernel function W(r, h) — a bell-shaped function that assigns a weight to each neighbouring particle based on distance r and smoothing length h. The kernel must be normalised (integrates to 1), symmetric (depends only on |r|), and vanish beyond some cutoff distance (typically 2h).

The most commonly used kernel in real-time SPH is the Poly6 kernel for density estimation and the Spiky kernel for pressure forces. The Spiky kernel, proposed by Müller et al. in their 2003 SIGGRAPH paper that brought SPH to games and computer graphics, has a non-zero gradient at the particle itself — crucial for computing repulsive pressure forces that prevent particles from collapsing together.

The smoothing length h is typically set to cover 20–50 particles in 2D. Larger h means smoother but less detailed fluid; smaller h preserves fine features but requires more particles for stability.

Density and Pressure: Step by Step

Each simulation step follows a sequence of computations for every particle:

• 1 Find neighbours. For each particle i, collect all particles j within distance 2h. A spatial hash grid makes this O(N) rather than O(N²).
• 2 Compute density. ρ_i = Σ_j m_j W(|r_i − r_j|, h) — a weighted sum of neighbouring particle masses.
• 3 Compute pressure. Using an equation of state: P_i = k(ρ_i − ρ₀), where k is the stiffness constant and ρ₀ is the rest density.
• 4 Compute pressure force. F_pressure,i = −Σ_j m_j (P_i + P_j) / (2ρ_j) ∇W(|r_i − r_j|, h)
• 5 Compute viscosity force. F_visc,i = μ Σ_j m_j (v_j − v_i) / ρ_j ∇²W(|r_i − r_j|, h)
• 6 Integrate. Add gravity, apply total force, update velocity and position using a leapfrog or symplectic Euler integrator.

Pressure Forces and the Equation of State

Real water is nearly incompressible — its density barely changes under pressure. SPH with a simple equation of state P = k(ρ − ρ₀) is weakly compressible: it allows small density fluctuations in exchange for avoiding the expensive solve that strict incompressibility requires. The stiffness parameter k must be high enough to keep density variation below a few percent, but not so high that the timestep must be made impractically small to keep the simulation stable.

For more accurate incompressible SPH — used in engineering applications — algorithms like IISPH (Implicit Incompressible SPH) or DFSPH (Divergence-Free SPH) solve for pressure implicitly at each timestep. These converge to near-zero density deviation but require iterative solvers that are too slow for real-time rendering.

Surface Tension and Viscosity

Surface tension in SPH is modelled using either the colour field method (compute the interface normal from density gradients and apply a curvature-based force) or cohesion forces between particles. Our simulation uses a simplified cohesion model that creates the characteristic tendency of water drops to pull back into spheres.

Viscosity damps relative velocity between neighbouring particles — the viscosity force term makes fast-moving particles slow down neighbours and vice versa. High viscosity produces honey-like flow; zero viscosity (inviscid flow) produces sharp, unrealistic velocity discontinuities. The Laplacian kernel used for viscosity — with its positive second derivative — ensures the viscosity force is always dissipative, preventing energy injection.

Why SPH Beats Grid Methods for Free Surfaces

The decisive advantage of SPH for free-surface flows is the absence of any interface tracking. In grid-based methods, the free surface must be explicitly tracked using methods like VOF (Volume of Fluid) or Level Set. These add significant complexity and can produce numerical diffusion that smears sharp interfaces over time.

SPH's surface is wherever the particles are. Drops naturally separate, merge, and bounce because the particles simply follow the forces. This makes SPH the dominant method in real-time fluid simulation for games and interactive applications, and a strong competitor in offline visual effects — the ocean simulation in many modern films is SPH or a hybrid SPH-grid method.

The trade-off is that SPH struggles with thin fluid sheets and low-viscosity turbulent flows, where grid methods excel. Real production fluid simulators often use hybrid approaches — SPH near the surface where free-surface tracking matters, grid methods in the interior where incompressibility must be enforced precisely.