SPH Without Tears — Solving Fluids with Particles, from Kernel to Artificial Viscosity

In 1992 Joe Monaghan made a provocation in an astrophysics review article. Stop laying grids around stars, he argued. Instead, scatter the gas as thousands of overlapping blobs and let each blob feel its neighbors through a smooth pressure field. Thirty years later that idea — Smoothed Particle Hydrodynamics, or SPH — owns a quiet but durable corner of astrophysics, free-surface flow, and impact simulation. This post tracks the four equations that make SPH work, the traps people fall into, and lets you push particles around in the browser.

A fluid solver without a mesh#

Eulerian methods nail down the coordinate system and watch the flow stream past. SPH does the opposite. Each particle carries its own velocity and density and rides along with the flow. In Lagrangian form the continuity and momentum equations are

Here $\rho$ is density, $\mathbf{u}$ velocity, $P$ pressure, $\mathbf{g}$ gravity. The time derivative on the left becomes a plain $dt$ along the particle path. The catch is the spatial derivative on the right. Without a regular grid of neighbors, how do we evaluate $\nabla P$? That is the question SPH was built to answer.

The kernel W(r, h) — a virtual volume around each particle#

The trick is convolution. Drop a bell-shaped function $W(r, h)$ around every particle and approximate any field $A(\mathbf{x})$ as a weighted sum over neighbors:

$m_j$ is the particle mass, $\rho_j$ its density, $h$ the smoothing length. The ratio $m_j/\rho_j$ behaves like the volume that particle $j$ occupies. Differentiation just moves onto the kernel — $\nabla A \approx \sum_j (m_j/\rho_j) A_j \nabla W$. The grid derivative becomes a kernel derivative.

A valid kernel must integrate to one and converge to a Dirac delta as $h \to 0$. The most popular choices have compact support — they vanish outside $r < h$ or $r < 2h$, so each particle only talks to a handful of neighbors. Compare the four classics below.

The cubic spline drops smoothly to zero at $r = 2h$. The poly6 kernel is great for density but its gradient flattens near $r = 0$, so graphics SPH uses a separate spiky kernel for pressure. Gaussians have no compact support, which kills performance. Astrophysics SPH almost always uses cubic spline.

Density and pressure as sums#

Density at particle $i$ is just a kernel sum over neighbors:

Pressure follows from an equation of state. For weakly compressible free-surface flow the Tait EOS is the workhorse:

$\rho_0$ is the rest density, $K$ the stiffness, usually $\gamma = 7$. A small density swing produces a strong pressure response, mimicking incompressibility. Picking $K$ so the sound speed $c_0 = \sqrt{\gamma K / \rho_0}$ is roughly ten times the maximum flow speed keeps the Mach number low.

A symmetry trick that saves momentum#

The naive SPH form of $\nabla P / \rho$ does not conserve momentum. The fix is the identity

Plugging it in produces a pressure term that is symmetric in pairs:

The force on $i$ from $j$ is exactly the negative of the force on $j$ from $i$. Newton's third law holds on every pair, total momentum and angular momentum drift only with floating-point noise. This symmetrization is half the secret of SPH stability.

Artificial viscosity for shocks#

The bare equations let particles tunnel through each other when a shock arrives. Monaghan and Gingold (1983) added an artificial viscosity that only switches on for approaching pairs:

Typical values are $\alpha \approx 0.5$–$1.0$ and $\beta = 2\alpha$. Because the viscosity vanishes in tension, it does not glue particles in expansion regions. The downside is that it also fires inside shear layers, which inspired switches like Balsara's later on.

A 1D column collapse in Python#

Let us watch eighty particles, stacked vertically, collapse under gravity.

import numpy as np
 
def cubic_spline_1d(r, h):
    """1D cubic spline kernel W(r, h)."""
    q = abs(r) / h
    c = 2.0 / (3.0 * h)
    if q < 1:
        return c * (1 - 1.5 * q * q + 0.75 * q ** 3)
    if q < 2:
        return c * 0.25 * (2 - q) ** 3
    return 0.0
 
def cubic_spline_dW(r, h):
    """First derivative dW/dr; sign follows r."""
    q = abs(r) / h
    c = 2.0 / (3.0 * h * h)
    s = 1.0 if r > 0 else (-1.0 if r < 0 else 0.0)
    if q < 1:
        return c * (-3 * q + 2.25 * q * q) * s
    if q < 2:
        return c * (-0.75 * (2 - q) ** 2) * s
    return 0.0
 
def sph_density_pressure(x, m, h, rho0, K, gamma=7.0):
    n = len(x)
    rho = np.zeros(n)
    for i in range(n):
        for j in range(n):
            rho[i] += m * cubic_spline_1d(x[i] - x[j], h)
    p = K * ((rho / rho0) ** gamma - 1.0)
    return rho, p
 
def sph_acceleration(x, v, rho, p, m, h, alpha, c0, g):
    n = len(x)
    a = np.full(n, g)
    for i in range(n):
        for j in range(n):
            if i == j:
                continue
            r = x[i] - x[j]
            if abs(r) >= 2 * h:
                continue
            grad = cubic_spline_dW(r, h)
            a[i] -= m * (p[i] / rho[i] ** 2 + p[j] / rho[j] ** 2) * grad
            v_ij = v[i] - v[j]
            if v_ij * r < 0:                       # approaching pair
                mu = h * v_ij * r / (r * r + 0.01 * h * h)
                rho_bar = 0.5 * (rho[i] + rho[j])
                Pi = -alpha * c0 * mu / rho_bar
                a[i] -= m * Pi * grad
    return a
 
# 1D column: 80 particles spread evenly over y in [0.5, 1.0]
N = 80
x = np.linspace(0.5, 1.0, N)
v = np.zeros(N)
m = 0.5 * 1000.0 / N        # column mass / particle count
h = 0.018
rho0, K = 1000.0, 200.0
alpha, c0 = 0.5, 20.0
g, dt = -9.81, 5e-5
 
for step in range(8000):
    rho, p = sph_density_pressure(x, m, h, rho0, K)
    a = sph_acceleration(x, v, rho, p, m, h, alpha, c0, g)
    v += a * dt
    x += v * dt
    below = x < 0.0                                 # reflect at bottom wall
    x[below] = -x[below]
    v[below] = -0.5 * v[below]
 
print(f"final mean rho = {rho.mean():.1f}, top particle y = {x.max():.3f}")

After 8,000 steps the mean density oscillates around $\rho_0$ and the highest particle settles near $y \approx 0.55$ on top of a denser bed. That tiny script already contains the entire SPH skeleton — sums for density, EOS for pressure, symmetric pairs for acceleration, artificial viscosity for stability.

Push the particles around yourself#

Try the simulation below directly. Pushing the smoothing length $h$ up to 0.07 gives each particle more neighbors and the colors smooth out — density estimates get less noisy. Drop $h$ back to 0.03 and gaps appear between particles, with pressure oscillations rising sharply.

Cranking up the stiffness $K$ approximates incompressibility better but raises the sound speed, so the time step has to shrink. Setting viscosity to zero produces interpenetration streaks where particles slip past each other — a quick visual answer to why artificial viscosity exists.

Where SPH stumbles#

Every bright side has a shadow. SPH bleeds numerical viscosity in shear-dominated flows, often swamping the physical value. Cluster patterns break isotropy, which hurts surface-tension-driven phenomena. Variable smoothing lengths $h_i$ buy adaptivity but force corrective terms in momentum (Springel and Hernquist 2002 introduced the $f_i$ correction). Near a free surface the density count drops abruptly, pressure can go negative, and tensile instability ripples the surface — XSPH and $\delta$-SPH formulations were designed to suppress that.

What to take away#

Three takeaways. SPH replaces grid derivatives with kernel derivatives over a Lagrangian sum. Without symmetrizing the pressure term you lose momentum conservation, and without artificial viscosity shocks make particles tunnel. Compact-support kernels plus a sensible $h$ govern both cost and accuracy. Next time a free-surface or explosion simulation lands on your desk, consider 10,000 floating particles instead of yet another adaptive mesh.

References#

• Monaghan, J. J. (1992). Smoothed Particle Hydrodynamics. Annual Review of Astronomy and Astrophysics, 30, 543.
• Monaghan, J. J., and Gingold, R. A. (1983). Shock simulation by the particle method SPH. Journal of Computational Physics, 52(2), 374.
• Müller, M., Charypar, D., and Gross, M. (2003). Particle-based fluid simulation for interactive applications. SCA '03.
• Springel, V. (2005). The cosmological simulation code GADGET-2. MNRAS, 364, 1105.