SPH Particle Fluids

Splash a bucket of water across a floor and it does something no polygon mesh ever wants to do: it tears apart, throws off droplets, sheets into a thin film, then puddles and merges again. The topology changes every frame. One way to tame this chaos for a film shot is to stop thinking of the water as a shape at all, and think of it as a crowd of tiny particles — each carrying a little mass, each nudging its neighbours — that together behave like a fluid.

That is Smoothed Particle Hydrodynamics (SPH): the Lagrangian (particle-following) way to simulate liquids, as opposed to the Eulerian (fixed-grid) way of Stable Fluids. Born in 1977 for simulating stars in astrophysics, SPH turned out to be a wonderfully natural fit for the splashy, merging, tearing liquids that grids find awkward. This page builds the SPH idea from its one central trick — the kernel-weighted sum — up to a working pressure force, and is honest about the two things it makes you pay for.

The one idea: estimate any field by a weighted sum over neighbours

In SPH the fluid is the particles. Each particle j carries a mass m_j, a position \mathbf{x}_j, a velocity, and whatever field values it samples (density, pressure, colour…). The magic is a recipe for reading off the value of any field A at any point in space, even where no particle sits, by a smoothed average of the nearby particles:

A(\mathbf{x}) \;=\; \sum_{j} m_j \, \frac{A_j}{\rho_j} \, W(\mathbf{x} - \mathbf{x}_j,\, h).

Here \rho_j is particle j's density and W is the smoothing kernel: a smooth, bump-shaped weight that is largest at the centre and falls to exactly zero beyond a support radius h. So the sum is really only over the handful of particles within h of \mathbf{x} — near neighbours count a lot, far ones not at all, and everything past h contributes nothing. The factor m_j/\rho_j is just the little volume that particle j represents, so the formula is a volume-weighted, blurred sample of the field.

The fluid is discretised into particles carrying mass, not into grid cells.
Any field value at a point is a kernel-weighted sum over the neighbours within the support radius h.
The kernel W(\mathbf{r}, h) is smooth, normalised (\int W \, dV = 1), peaks at \mathbf{r}=0, and is exactly zero for |\mathbf{r}| > h.
Because W is a known smooth function, derivatives are free: \nabla A(\mathbf{x}) = \sum_j m_j \tfrac{A_j}{\rho_j}\, \nabla W(\mathbf{x}-\mathbf{x}_j, h) — you differentiate the kernel, not the messy data.

That last point is the quiet superpower. In a grid method you approximate gradients with finite differences between cells; in SPH the gradient of a field is a sum against the analytic gradient of a smooth kernel you chose. Pressure gradients, viscosity Laplacians — all become sums over neighbours of \nabla W or \nabla^2 W.

The neighbourhood picture

Everything in SPH is local. To evaluate any quantity at particle i you look only inside a disc (a sphere in 3-D) of radius h centred on it. The dark points below sit inside that radius and contribute to i's sums; the faint points are farther than h away and the kernel zeroes them out entirely.

Choosing h is a real trade-off: too small and a particle has too few neighbours, so its sums are noisy and lumpy; too large and every particle sees hundreds of others and the simulation crawls. A common target is a few dozen neighbours inside the radius. And because the sums only ever touch nearby particles, the whole cost hinges on being able to find those neighbours quickly — which is why SPH is glued to a neighbour search structure (more below).

From particles to a force: density → pressure → motion

The loop that makes a puddle move has three beats each step. First, set A = \rho in the master formula. The \rho_j and the A_j = \rho_j cancel, leaving the beautifully simple density estimate:

\rho_i \;=\; \sum_{j} m_j \, W(\mathbf{x}_i - \mathbf{x}_j,\, h).

A particle is dense where its neighbours crowd in close (the kernel is large) and light where they spread out. Second, turn that density into a pressure with a stiff equation of state — the water pushes back when squeezed past its rest density \rho_0:

p_i \;=\; k\,(\rho_i - \rho_0),

where k is a stiffness constant. Third, drive the particles with the pressure force, which pushes fluid from high pressure to low — the negative gradient of pressure, evaluated as a kernel sum (a symmetrised form keeps the forces equal-and-opposite):

\mathbf{f}^{\text{pressure}}_i \;=\; -\sum_{j} m_j \, \frac{p_i + p_j}{2\,\rho_j}\, \nabla W(\mathbf{x}_i - \mathbf{x}_j,\, h).

Add a viscosity force (a kernel Laplacian of velocity that smooths neighbouring velocities together — this is the fluid's internal friction) and a surface-tension force (which pulls the ragged boundary inward, letting droplets ball up), throw in gravity, and you have the full acceleration \mathbf{a}_i = (\mathbf{f}^{\text{pressure}}_i + \mathbf{f}^{\text{visc}}_i + \mathbf{f}^{\text{tens}}_i + \mathbf{f}^{\text{ext}}_i)/\rho_i. Integrate that forward in time and the crowd of particles sloshes like water.

Worked example: one particle's density and pressure

Let's do the first two beats by hand. Take a particle i with exactly two neighbours inside its radius, and use the classic poly6 density kernel in 2-D-ish constant form. To keep the arithmetic clean we'll use a kernel whose value at distance r is

W(r, h) \;=\; \frac{4}{\pi h^{2}}\left(1 - \frac{r^{2}}{h^{2}}\right)^{3} \quad \text{for } r \le h, \; 0 \text{ otherwise.}

Say every particle has mass m = 1, the support radius is h = 2, and i's two neighbours sit at distances r_1 = 1 and r_2 = 1.5. Include the particle's own contribution at r_0 = 0. The kernel prefactor is 4/(\pi \cdot 4) = 1/\pi \approx 0.318.

Self term (r_0=0): (1 - 0)^3 = 1, so W_0 = 0.318.
Neighbour 1 (r_1=1): 1 - 1/4 = 0.75, cubed is 0.4219, so W_1 = 0.318 \times 0.4219 \approx 0.1342.
Neighbour 2 (r_2=1.5): 1 - 2.25/4 = 0.4375, cubed is 0.0837, so W_2 = 0.318 \times 0.0837 \approx 0.0266.

Sum with the masses for the density:

\rho_i \;=\; 1\cdot 0.318 + 1\cdot 0.1342 + 1\cdot 0.0266 \;\approx\; 0.479.

Now the pressure. Suppose the rest density is \rho_0 = 0.4 and the stiffness is k = 100. Then

p_i \;=\; k\,(\rho_i - \rho_0) \;=\; 100\,(0.479 - 0.4) \;=\; 7.9.

Because i is a touch denser than rest, its pressure is positive, so the pressure force will push its neighbours away — the fluid resisting compression, exactly as it should. If those neighbours had been farther out, \rho_i would drop below \rho_0, the pressure would go negative, and the force would pull inward instead.

Finding the neighbours: a spatial grid

Every sum is over particles within h. Checking all N particles against all others is O(N^2) — hopeless for the hundreds of thousands of particles a splash needs. The standard fix is a uniform spatial grid (or spatial hash) with cell size h: bin every particle into its cell, and to find a particle's neighbours you only scan its own cell and the eight (in 3-D, twenty-six) touching cells. Every neighbour within h is guaranteed to live in that 3\times3\times3 block, so the search drops to roughly O(N). Rebuilding this grid each step is often the single biggest chunk of an SPH frame's time budget.

Because there is no shape to keep consistent. In a mesh-based or grid-based method, water tearing into droplets is a topology change — the surface has to split, and tracking that split is fiddly and error-prone. In SPH there is nothing to split: a droplet flying off is just a few particles that happened to drift beyond kernel range of the rest, so they stop feeling the main body's pressure and sail away on their own. Two streams merging is just their particles coming within h of each other and starting to share forces. Splitting and merging aren't special cases you code — they fall out of the neighbourhood sums for free. That is exactly why SPH is the go-to for the wild, foamy, topology-changing water that grids dread.

Two traps bite everyone who writes their first SPH solver.

1. Weakly-compressible SPH needs a stiff pressure → tiny timesteps. The simple equation of state p = k(\rho - \rho_0) only fights compression if k is large — otherwise the "water" squashes like a sponge and looks springy and wrong. But a stiff pressure means huge, sudden restoring forces, and an explicit integrator can only follow those with a very small timestep. Push the step too far and the fluid jitters, boils, and then explodes — particles flung to infinity. That is the weakly-compressible (WCSPH) bargain: stiff enough to look incompressible, and therefore slow. This is why the field invented PCISPH, IISPH and DFSPH, which enforce (near-)incompressibility with predictor–corrector or implicit pressure solves and so allow far larger, stable steps than raw WCSPH.

2. There is no explicit surface. You have a cloud of points, not a water boundary. Render the particles raw and you get a lumpy field of blobs, not a smooth liquid. To get a renderable surface you must reconstruct one from the particles each frame — typically build a density (or signed-distance) field on a grid and run marching cubes (or a screen-space filter) to extract the isosurface. Forgetting this step is why a first SPH render looks like a swarm of beads instead of water.

Contrast: particles (Lagrangian) vs grid (Eulerian)

SPH's opposite number is the grid approach of Stable Fluids, where the fluid quantities live on a fixed mesh of cells and the fluid "flows past" the grid. Each style is strong exactly where the other is weak.

Question	SPH (Lagrangian, particles)	Grid (Eulerian, cells)
Where does data live?	on moving particles	at fixed grid points
Mass conservation	exact — particles never vanish	needs care at cell boundaries
Splash / droplets / merging	free — no topology to track	hard — must track the surface
Incompressibility	awkward — stiff EOS or PCISPH/IISPH/DFSPH	natural — a global pressure Poisson solve
Empty space cost	free — no particles, no work	you still store empty cells
Renderable surface	must be reconstructed from points	often a level set already present

The rough rule of thumb: reach for particles when the shot is splashy, thin, and topology-changing (a crashing wave, a thrown bucket, spray), and for a grid when you want smooth, strongly incompressible bulk motion (a calm pool, smoke, a filling tank). Modern production solvers often blend the two — FLIP/PIC methods carry particles and a grid, taking incompressibility from the grid and detail from the particles.