SPH: Building a Fluid Out of Particles

A fluid made of particles

There are two ways to simulate a fluid. The Eulerian way puts a grid in space and asks what flows through each cell; the Lagrangian way puts particles in the fluid and lets them carry it around. Smoothed Particle Hydrodynamics is the Lagrangian one. It was invented in 1977 by Gingold, Monaghan and (independently) Lucy for astrophysics, where there is no box to put a grid in, and it later moved into computer graphics and engineering because it handles splashes, free surfaces and separating droplets without any special-case code.

The central idea is interpolation with a smoothing kernel. Any field A can be reconstructed at any point x by summing the contributions of nearby particles, each weighted by a kernel W that falls off with distance and vanishes beyond a smoothing radius h:

A(x) = Σ_j  m_j · (A_j / ρ_j) · W(x − x_j, h)

The elegance is that derivatives of the field become derivatives of the kernel — an analytic function you differentiate once, on paper. There is no grid, no advection term, no numerical diffusion from interpolating between cells; the mass is exactly conserved because it rides on the particles.

Kernels

A kernel must be normalised (integrate to 1), have compact support (zero beyond h, so the sums are local), and be smooth enough that its gradient is well behaved. In practice several kernels are used for different terms, following Müller, Charypar and Gross's 2003 formulation, which is what most real-time SPH implementations are built on:

poly6      → density        smooth, cheap, but its gradient
                             vanishes at the centre
spiky      → pressure force  non-zero gradient at r → 0, which is
                             what stops particles clumping
viscosity  → viscous force   Laplacian is positive everywhere, so it
                             cannot inject energy

Using poly6 for the pressure gradient is a classic beginner's bug: because its gradient goes to zero as particles approach each other, the repulsive force disappears exactly when it is needed, and particles collapse into clusters. The spiky kernel exists precisely to fix that. The cubic spline kernel is the other common choice and is the standard in astrophysical SPH.

Density, pressure and the Tait equation

Every step begins with density. Each particle sums the masses of its neighbours weighted by the kernel:

ρ_i = Σ_j  m_j · W(x_i − x_j, h)      // include the particle itself

Then pressure. Truly incompressible flow would require solving a global Poisson equation each step, which is expensive; weakly compressible SPH instead uses an equation of state that punishes compression hard enough that the fluid stays almost incompressible. The Tait equation is the standard one:

p_i = B · ((ρ_i / ρ₀)^γ − 1)          // γ = 7 for water
B   = ρ₀ · c² / γ                     // c = an artificial speed of sound

The exponent γ = 7 makes the pressure rise brutally steeply once the density exceeds the rest density ρ₀, so a small density error produces a large restoring force. The stiffness constant B is set from an artificial speed of sound c, chosen not to be the real 1500 m/s of water but roughly ten times the fastest expected flow speed — that keeps density fluctuations under about 1% while allowing a time step thousands of times larger than the physical sound speed would permit. A simpler linear state equation, p = k(ρ − ρ₀), is often enough for a visual real-time demo.

The pressure force on a particle uses the symmetric form, which conserves momentum exactly because the force on i from j is equal and opposite:

f_i^pressure = − Σ_j  m_j · (p_i/ρ_i² + p_j/ρ_j²) · ∇W(x_i − x_j, h)

Viscosity and surface tension

Viscosity smooths the velocity field: each particle is dragged towards the average velocity of its neighbours, weighted by the Laplacian of the viscosity kernel. It sets how syrupy the fluid looks, and it also quietly stabilises the simulation by damping the high-frequency particle noise that the pressure term creates. Astrophysical SPH uses a different device, artificial viscosity (the Monaghan α–β term), whose job is to allow shocks to form without the particles interpenetrating.

Surface tension needs a little machinery, because a particle in the interior has no idea it is in the interior. The standard trick is a colour field: assign every particle the value 1 and smooth it with the kernel. Its gradient is nearly zero deep inside the fluid and points outward at the free surface, where the neighbour sum is one-sided. The gradient's magnitude therefore acts as a surface detector, its direction is the surface normal, and the divergence of the normalised normal is the curvature κ. The force is then f = −σ · κ · n̂, which pulls the surface flat and is what makes a free-falling blob become a sphere and small droplets stay round.

Neighbour search: the part that decides your frame rate

Every one of those sums runs over "neighbours within h". Done naively that is an O(n²) scan per step and it dominates everything else. The fix is a uniform spatial hash grid with a cell size of exactly h: then the neighbours of a particle can only be in its own cell and the 8 (2D) or 26 (3D) surrounding cells.

const c = 1 / h;                      // cell size = smoothing radius
const key = (x, y) => (Math.floor(x * c) * 92837111 ^
                       Math.floor(y * c) * 689287499) >>> 0;

// rebuild each step: counting sort into a flat bucket array
buckets.fill(0);
for (const p of particles) buckets[key(p.x, p.y) % NB]++;
// prefix-sum, scatter, then iterate only the 3×3 cell block per particle

Rebuilding the grid is O(n) with a counting sort, and the neighbour query becomes O(k) where k is the average neighbour count — typically kept around 20 in 2D and 30–40 in 3D by tuning h against the particle spacing. Too few neighbours and the density estimate is noisy; too many and every step gets quadratically more expensive in h. Storing particles sorted by cell also fixes the other half of the performance problem: memory locality, since neighbouring particles then sit next to each other in the array.

Time step and the CFL condition

SPH is explicit, so the step is bounded. A particle must not move further than a fraction of the smoothing radius in one step, and the step must also resolve the artificial sound speed and the viscous and force time scales. The standard CFL bounds are:

Δt ≤ 0.25 · h / (c + v_max)           // sound / advection
Δt ≤ 0.25 · sqrt(h / |a_max|)         // acceleration (e.g. gravity)
Δt ≤ 0.125 · h² / ν                   // viscous diffusion

Take the minimum of the three, with a safety factor. Violate them and the pressure term — which is very stiff, precisely because γ = 7 — blows up: particles overlap, density spikes, the restoring force overshoots, and the simulation explodes into a spray of NaNs within a handful of frames. This is why the artificial speed of sound is kept as low as the target compressibility allows: c appears in the denominator of the first bound, so every factor you add to c costs you directly in step size. The integrator is usually Leapfrog for the same reason as in the N-body case — it is symplectic, it costs one force evaluation, and the forces here are far too expensive to evaluate four times per step for RK4.

Frequently asked questions

Try it live

Everything above runs in your browser — open SPH Fluid and change the parameters while it is running. Nothing is installed, nothing is uploaded, the whole model lives in one tab.