Fluid Dynamics
July 2026 · 15 min read · SPH · Free Surface · Surface Tension · Last updated: 3 July 2026

SPH with Free Surface — Simulating Splashes, Droplets & Waves

Written by MySimulator Team · Reviewed by MySimulator Editorial Review

A wave breaking on a beach, a droplet detaching from a tap, water splashing out of a bucket — every one of these events involves a free surface: a boundary between fluid and empty space (or air) whose shape is not known in advance and must emerge from the simulation itself. Grid-based solvers struggle with this because they need to track where the fluid boundary sits inside a fixed mesh. Smoothed Particle Hydrodynamics sidesteps the problem almost for free — the surface is simply wherever the particles stop.

1. Why Particle Methods Suit Free Surfaces

In a grid-based (Eulerian) solver, every cell in the domain holds fluid-state values whether or not fluid actually occupies that region. Tracking the boundary between fluid and air requires an auxiliary representation — a level set, a volume-of-fluid (VOF) fraction field, or marker particles seeded purely for interface tracking. When a splash breaks into droplets, or two streams merge, that auxiliary field must handle topology changes correctly, which is a genuinely hard numerical problem (thin films disappear, droplets fail to pinch off cleanly, mass leaks away).

SPH is Lagrangian: the fluid is the particles. There is no separate interface field to maintain — the free surface is implicitly wherever particle density trails off. A droplet separating from a jet is just a group of particles whose neighbour connectivity changes; no special-case code is required. This is why SPH (and its cousins, MPM and PBF) dominate splash-heavy VFX and many research CFD codes for dam breaks, wave impact, and coastal engineering.

Trade-off: what free-surface tracking gives away for free, SPH pays for elsewhere — pressure noise, weaker incompressibility enforcement than grid pressure projection, and O(N) to O(N log N) neighbour search cost that grid methods do not incur.

2. The Density Deficiency Problem

SPH estimates density at a particle i by summing a kernel-weighted contribution from every neighbour j within the smoothing radius h:

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

Deep inside the bulk of the fluid, a particle is surrounded on all sides, so this sum converges to the correct rest density. Near a free surface, however, the kernel support extends into empty space where no particles exist — the sum under-counts, and computed density comes out artificially low, even though the true physical density at the surface is exactly the same as in the bulk.

This is the density deficiency problem. Left uncorrected, the equation of state (used to convert density into pressure in weakly-compressible SPH) reads this deficiency as negative pressure, which incorrectly pulls surface particles inward — the simulated fluid appears to have an artificial "skin" that resists the natural formation of thin sheets and fine splashes.

P = k(ρ − ρ₀) // Tait / weakly-compressible EOS ↑ near the surface ρ < ρ₀ even for physically correct fluid → P < 0 → spurious inward pull

Common mitigations include clamping negative pressure to zero for free-surface particles, using a corrected (renormalised) density kernel that accounts for the missing kernel mass, or switching to a divergence-free / incompressible SPH (IISPH, DFSPH) formulation where pressure is solved implicitly rather than read directly off a possibly-deficient density estimate.

3. Detecting the Free Surface

Many parts of the simulation — surface tension, rendering, and adaptive resolution — need to know which particles actually sit on the boundary. The classic technique is the colour field method (Morris, 2000): treat each particle as carrying a scalar "colour" c = 1 for fluid and interpolate it like any other SPH field:

c(r) = Σⱼ (mⱼ/ρⱼ) · W(r − rⱼ, h) n = ∇c // surface normal points from fluid toward air

In the bulk, the gradient ∇c is close to zero because the colour field is roughly uniform in every direction. Near the surface, ∇c has significant magnitude and points outward — this both flags surface particles (|∇c| above a threshold) and supplies the surface normal needed for rendering and surface-tension forces.

A cheaper alternative counts neighbour deficiency directly: a particle whose neighbour count, or whose kernel-weighted density, falls below a fraction of the expected bulk value (typically 80–95%) is classified as a surface particle. This avoids computing a gradient at all and is fast enough to run every substep.

Divergence-based detection (comparing Σⱼ (rⱼ − rᵢ)·∇W against a threshold, per Marrone et al.'s CSPM criterion) is more robust in thin sheets and highly curved splashes than a simple neighbour count, at a modest extra cost.

4. Surface Tension — the CSF Model

Surface tension arises from an imbalance of intermolecular cohesive forces at an interface: molecules in the bulk are pulled equally in all directions, but molecules at the surface feel a net inward pull, which macroscopically manifests as a force proportional to interface curvature (the Young–Laplace relation, ΔP = σκ). The Continuum Surface Force (CSF) model of Brackbill, Kothe & Zemach (1992) converts this into a per-particle body force using the colour field from the previous section:

κ = −∇²c / |∇c| // surface curvature Fₜₑₙₛᵢₒₙ = −σ · κ · n · δₛ // n = ∇c / |∇c|, δₛ = |∇c| (surface delta)

Here σ is the surface tension coefficient (≈ 0.072 N/m for water at 20 °C), κ is the local curvature computed from the Laplacian of the colour field, and δₛ localises the force to a thin band near the actual interface. In practice ∇²c is noisy at typical SPH particle spacing, so production implementations often use the simpler and more robust pairwise cohesion model of Akinci et al. (2013), which applies a short-range attractive spline force directly between neighbouring particles instead of differentiating a field twice:

Fᵢⱼ = −σ · mᵢmⱼ · C(|rᵢⱼ|) · r̂ᵢⱼ C(r) = (H−r)³r³ for 0 < r ≤ H/2 2(H−r)³r³ − H⁶/64 for H/2 < r ≤ H 0 otherwise

This spline C(r) is carefully normalised so integrating it over a full neighbourhood gives zero net force in the bulk (cohesion cancels out) while producing a net inward pull exactly at curved surfaces — reproducing droplet rounding, capillary rise in thin tubes, and the beading of water on non-wetting surfaces without ever computing a second derivative.

5. Tensile Instability and Particle Clumping

Standard SPH pressure forces use the gradient of a smoothing kernel (typically the Spiky kernel) which becomes attractive rather than repulsive under negative pressure — exactly the regime that free surfaces and surface tension push particles into. The result is the well-known tensile instability: particles clump into an unphysical lattice or cluster instead of forming a smooth sheet, most visible in thin splashes and droplets where negative pressure is common.

Monaghan's artificial pressure correction adds a small repulsive term whenever the kernel value grows relative to its value at the initial particle spacing, directly counteracting the attractive branch:

f(rᵢⱼ) = ( W(rᵢⱼ) / W(Δp) )ⁿ Fₜₑₙₛᵢₗₑ = ε · f(rᵢⱼ)ⁿ · ∇W(rᵢⱼ) // ε ≈ 0.01–0.2, n ≈ 4

An alternative that avoids tuning constants is to switch kernel families: the Wendland C2 kernel has a gradient that stays finite and well-behaved at small particle separations and is the de facto standard in modern SPH codes (SPlisHSPlasH, DualSPHysics) precisely because it resists clumping without an ad hoc patch term.

6. Pressure Solvers Near the Surface — WCSPH vs PCISPH

Weakly Compressible SPH (WCSPH) reads pressure directly from the Tait equation of state shown above. It is simple and fully explicit, but requires a stiff equation of state (large k) to keep density errors small, which in turn forces a tiny timestep via the Courant–Friedrichs–Lewy condition — the numerical sound speed must stay comfortably above the fastest fluid velocity in the scene, and splashes routinely produce fast-moving droplets that make this expensive.

Predictive-Corrective Incompressible SPH (PCISPH) and its relatives (IISPH, DFSPH) instead solve for pressure implicitly, iterating a predict-correct loop each step until the density deviates from rest density by less than a tolerance (often <1%):

repeat: predict velocities & positions using current pressure recompute density ρ* from predicted positions Δpᵢ += δ · (ρ* − ρ₀) // pressure update, δ from a scaling factor until max |ρ* − ρ₀| < tolerance

Near a free surface this iterative correction must be paired with the density-deficiency fix from Section 2 — otherwise the solver "corrects" the artificially low surface density by inflating pressure, producing an over-stiff, trampoline-like surface. Most modern implementations detect surface particles first and either clamp their target density or skip the incompressibility constraint for them entirely, letting surface tension and gravity govern their motion instead.

💧 Interactive Fluid Simulation Watch particle-based fluid form waves, splashes, and droplets in real time

7. Solid Boundaries and Ghost Particles

A free surface is only half the boundary problem — SPH fluid must also interact correctly with solid walls (a glass, a beach, a boat hull). The simplest approach lines solid surfaces with static boundary particles that participate in density and pressure summation exactly like fluid particles, but never move. Their presence fills in the kernel support at the wall, preventing fluid particles from reporting a spurious density deficiency there (the same issue as Section 2, but at a rigid rather than free interface).

Naively-placed boundary particles still let fluid penetrate corners or thin walls at high speed, so refinements are common: ghost particles mirrored across the boundary at simulation time, pressure-extrapolated boundary handling (Akinci 2012) that computes a physically consistent pressure for boundary particles instead of a fixed value, and packed multi-layer boundary shells (2–3 particles deep) that give the kernel full support even at grazing angles.

A common bug: forgetting to include boundary particles' mass in the fluid particle's density sum. The fluid then reports low density near walls exactly as it does near the free surface, and the pressure solver pushes fluid away from solid boundaries as if they were vacuum — an easy diagnostic if particles appear to "levitate" a few particle radii off the container floor.

8. Rendering the Surface — Marching Cubes and Metaballs

Particles alone render as a cloud of spheres, which reads as foam rather than a coherent liquid surface. Production renderers instead build a smooth implicit surface from the particle positions and extract a polygon mesh from it every frame:

For real-time web demos, a screen-space variant is far cheaper: render particles as depth-only spheres, blur the resulting depth buffer to smooth out particle graininess, reconstruct normals from the blurred depth via a screen-space gradient, and shade with refraction and Fresnel reflectance — no 3D voxel grid or mesh extraction required, which is what most WebGL/Three.js real-time SPH demos use to get a convincing water look at interactive frame rates.