📚 Learning Series June 2026 ~9 min read

How SPH Fluid Simulation Works: The Math Behind Fluid Particles

Real fluids are continuous — infinitely many molecules interacting simultaneously. SPH turns that into something a computer can handle: a finite set of particles that each carry a piece of the fluid's mass, momentum, and energy. Here is how the mathematics works.

The Problem with Simulating Fluids

Fluid behaviour is governed by the Navier-Stokes equations — a set of partial differential equations relating velocity, pressure, density, and viscosity across a continuous medium. The problem is the word continuous. Real water contains roughly 3 × 1025 molecules per litre. Simulating each one is impossible.

Two broad strategies exist for handling this. Grid-based methods (Eulerian approaches) divide space into fixed cells and track how fluid flows between them. Particle-based methods (Lagrangian approaches) move with the fluid, tracking parcels of mass as they travel through space.

Smoothed Particle Hydrodynamics (SPH) is the most widely used particle method. It was originally invented in 1977 by Lucy and Gingold & Monaghan for astrophysics simulations, and later adapted for fluid dynamics. Its key advantage is the natural handling of free surfaces — boundaries between fluid and empty space — which are notoriously difficult for grid-based methods.

The SPH Kernel: Spreading Point Masses into Volumes

The core mathematical trick in SPH is the kernel function W(r, h), also called the smoothing kernel. Every particle has a smoothing length h — a radius of influence. The kernel determines how a property carried at a point is distributed over the surrounding volume.

A kernel must satisfy several properties:

The cubic spline kernel, introduced by Monaghan, is the most common choice. It is smooth, computationally cheap, and well-behaved numerically. The Poly6 and Spiky kernels, introduced by Müller et al. in their landmark 2003 paper on interactive SPH, are used specifically for density and pressure calculations respectively — a key insight that dramatically improves simulation stability.

Estimating Quantities from Particles

Once you have a kernel, you can estimate any scalar or vector quantity A at a point r by summing contributions from all nearby particles j: 

A(r) = Σⱼ mⱼ (Aⱼ / ρⱼ) W(r − rⱼ, h)

Where mⱼ is the mass of particle j, Aⱼ is the value of A carried by that particle, ρⱼ is its density, and W is the kernel. The density itself is estimated the same way: 

ρ(r) = Σⱼ mⱼ W(r − rⱼ, h)

Derivatives are computed analytically by differentiating the kernel — no finite differences needed. The gradient of A at particle i is: 

∇A(rᵢ) = Σⱼ mⱼ (Aⱼ / ρⱼ) ∇W(rᵢ − rⱼ, h)

This is the mathematical machinery that translates the continuous Navier-Stokes equations into a discrete particle system.

Pressure Force

Fluid resists compression. When particles get too close together, density rises above a rest value ρ₀, and a restoring pressure force pushes them apart. This is computed via the equation of state — the simplest being the ideal gas law: 

p = k(ρ − ρ₀)

Where k is a stiffness constant. Higher k makes the fluid more incompressible but can cause numerical instability (particles bouncing violently). Real incompressible SPH uses more sophisticated pressure solvers (PCISPH or DFSPH) that iteratively enforce zero divergence, at higher computational cost.

The pressure force on particle i from its neighbours j is: 

fᵢ_pressure = −mᵢ Σⱼ mⱼ ((pᵢ + pⱼ) / (2ρⱼ)) ∇W(rᵢ − rⱼ, h)

The symmetrised pressure term (pᵢ + pⱼ)/2 is important — it ensures Newton's third law holds and the total momentum of the system is conserved.

Viscosity Force

Real fluids resist shearing — this is viscosity. Without it, fluid particles that pass near each other transfer no energy, and the simulation looks like bubbles rather than water. The SPH viscosity term adds a damping force proportional to the velocity difference between neighbouring particles: 

fᵢ_viscosity = μ Σⱼ mⱼ ((vⱼ − vᵢ) / ρⱼ) ∇²W(rᵢ − rⱼ, h)

Where μ is the dynamic viscosity coefficient and ∇²W is the Laplacian of the kernel. Müller et al. proposed a special kernel for this term whose Laplacian is always positive — avoiding the sign-change instabilities that plagued earlier implementations.

Surface Tension

Water forms droplets because interior molecules are attracted equally in all directions, while surface molecules feel a net inward pull. In SPH this is modelled with a colour field — a scalar that is 1 inside the fluid and 0 outside. The gradient of this field points toward the fluid surface; its curvature determines the surface tension force.

Surface tension is what makes SPH fluid simulations beautiful: the water beads up, droplets merge, and thin sheets naturally break. Without it, fluid particles disperse into a featureless cloud.

Why SPH Is Good for Free-Surface Flows

Grid-based methods track which grid cells contain fluid and which are empty. As the surface moves, cells must be classified, filled, and emptied — a bookkeeping challenge that requires special methods like the Volume of Fluid (VOF) or Level Set approaches.

SPH has no grid. The particles are the fluid. Wherever particles go, the fluid surface naturally follows. Splashing, merging, breakup, droplet formation — all emerge automatically from the particle dynamics without any explicit surface tracking.

This makes SPH the method of choice for:

See SPH in action in the browser — the interactive simulation at /fluid/ runs thousands of particles in real time. Try pouring fluid from different angles, or adding obstacles to see the flow split and reform.

The Performance Challenge

The naive SPH algorithm checks every particle against every other particle to find neighbours within radius 2h. That is O(n²) work — doubling the particle count quadruples the computation.

In practice, all real-time SPH implementations use a spatial hash: divide space into cells of size 2h and assign each particle to a cell. Neighbour queries only need to check the 27 surrounding cells (in 3D). This reduces complexity to O(n) on average for uniformly distributed particles.

Even so, thousands of particles per frame is close to the CPU limit for real-time simulation. Modern high-fidelity SPH uses GPU compute shaders to run the neighbour search and force calculation in parallel across thousands of GPU threads simultaneously, enabling millions of particles at interactive rates.

SPH Fluid Dynamics Mathematics Physics Learning Series

Frequently Asked Questions

What is Smoothed Particle Hydrodynamics (SPH)?

Smoothed Particle Hydrodynamics (SPH) is a computational fluid dynamics method that represents fluids as a collection of discrete particles rather than a fixed grid. Each particle carries properties like mass, density, pressure, and velocity, and interactions are computed by smoothing these properties over a kernel function.

How does pressure force work in SPH?

In SPH, pressure at each particle is computed from its local density using an equation of state. The pressure gradient force pushes particles from high-pressure to low-pressure regions, preventing compression. This is approximated by summing weighted contributions from neighboring particles within the smoothing kernel radius.

What is the smoothing kernel in SPH?

The smoothing kernel (also called the weighting function or kernel function) defines how much influence neighboring particles have based on distance. Popular choices include the cubic spline kernel and the Wendland kernel. The kernel must be symmetric, normalized, and approach a Dirac delta as the smoothing length approaches zero.

How is viscosity modeled in SPH?

Viscosity in SPH is typically modeled as an artificial viscosity term added to the force calculation between particle pairs. The viscous force resists relative motion between neighboring particles, smoothing out velocity differences. This term prevents particle interpenetration and stabilizes the simulation.

What determines the number of particles in an SPH simulation?

Particle count is a balance between accuracy and computational cost. More particles give higher resolution and smoother results but scale as O(N log N) for neighbor searches. Typical interactive simulations use 1,000–100,000 particles, while production simulations for films use millions.

What is the continuity equation in fluid simulation?

The continuity equation states that mass must be conserved — what flows in must flow out. In SPH, it governs how density changes over time based on particle velocities and their neighbors. It ensures the simulation maintains physical mass conservation throughout the computation.

How does surface tension work in SPH?

Surface tension in SPH is modeled using cohesive forces between particles at the fluid surface. Particles at the boundary have fewer neighbors than interior particles, creating an imbalance that pulls surface particles inward, mimicking the molecular forces responsible for real surface tension.

What are the limitations of SPH for fluid simulation?

SPH struggles with tensile instability (particles clumping artificially), maintaining incompressibility (often requiring small time steps), boundary conditions (walls are tricky to implement), and computational cost for high-resolution scenes. Techniques like position-based fluids or grid-based methods often complement SPH for different scenarios.

What is the difference between SPH and grid-based fluid simulations?

Grid-based methods (like finite differences or finite volumes) solve fluid equations on fixed spatial grids, making them efficient for large homogeneous fluid regions. SPH is Lagrangian — it moves with the fluid — making it natural for free surfaces, splashing, and large deformations, but less efficient for bulk fluid regions.

How is neighbor search optimized in SPH?

Naive neighbor search is O(N²) since every particle must check every other particle. SPH simulations use spatial hashing, uniform grids, or k-d trees to reduce this to O(N log N) or even O(N) in practice. Each particle only needs to check neighbors within the smoothing kernel radius h.