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. Four equations, and the traps that sit next to them.

Switch the kernel among Gaussian, cubic spline, poly6, and spiky. The bell and the slope at the center change together.

Poly6 goes flat near the center. Spiky stays pointed there. Which kernel belongs in the pressure term is decided by that difference.

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. Another weight that keeps a gradient on a skewed stencil is moving least squares.

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. Gaussian, cubic spline, poly6, and spiky differ in support width and in the slope at the center.

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.

Walk the same formulas with the script's numbers. $N=80$ over $y$ from 0.5 to 1.0, so the spacing is $0.5/79\approx0.00633$. Particle mass is $m=0.5\times1000/80=6.25$. Then $m/\rho_0=6.25/1000=0.00625$, almost the spacing: at rest density the particle "volume" matches the gap. With $h=0.018$, $2h=0.036$, and one side of the support holds $0.036/0.00633\approx5.7$ neighbors, five or six. Dropping $h$ to 0.03 still leaves particles inside the support. Go shorter than that and the sum gets sparse.

Put $\gamma=7$ and $K=200$ into Tait and a 1% density error is no longer small. At $\rho/\rho_0=1.01$, $(1.01)^7\approx1.072$, so $P=K\times0.072=14.4$. One percent of density becomes a pressure of 14. The same $K$, $\gamma$, and $\rho_0$ give a sound speed $c_0=\sqrt{7\times200/1000}=\sqrt{1.4}\approx1.18$. The artificial viscosity in the script does not use that value. It uses $c_0=20$. The viscous term is about an order of magnitude stiffer than the equation-of-state sound speed. $\alpha=0.5$, and the $\beta$ term is absent from the code. A gentle collapse survives on $-\alpha c_0\mu/\bar\rho$ alone. A hard shock, with the $\mu^2$ piece missing, lets particles overlap again.

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#

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#

Slide $h$ up to 0.07 and each particle sees more neighbors, so the colors even out. Drop it to 0.03 and gaps open between particles while the pressure starts to chatter.

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.

When h is shorter than a few spacings#

$h=0.018$ is about three times the spacing 0.00633. One side of the support radius $2h$ holds more than five neighbors. The density sum stays smooth at that width. Cut $h$ down toward one spacing and $W$ mostly sees the particle itself, so $\rho_i$ chatters. $\gamma=7$ turns that chatter into pressure. A step of $5\times10^{-5}$ run for 8000 steps is 0.4 s of time. Over that interval the top particle is said to fall from 1.0 to about 0.55. The pressure term is what caught the weight.

Drop the symmetric pair and the forces no longer cancel. Skip $\nabla(P/\rho)$ and the total momentum drifts every step. Artificial viscosity has a switch too. $\Pi_{ij}$ turns on only while the pair approaches. The script keeps only the linear piece at $\alpha=0.5$. On a free surface the density count falls, pressure goes negative, and particles pull on each other until the surface jitters. That pairing is the subject of tensile instability with a fixed kernel.

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.