Fluid Simulation Methods — SPH vs LBM vs MPM vs FVM

Every fluid simulator makes a fundamental choice: represent fluid as discrete particles (Lagrangian), as values on a fixed grid (Eulerian), or as a hybrid. Each choice determines accuracy, parallelisability, and what phenomena can be captured. This article surveys the four dominant approaches in games, VFX, and scientific computing.

1. Lagrangian vs Eulerian Framing

All fluid simulation methods can be classified by their reference frame:

Lagrangian: Track fluid parcels as they move. Each particle carries mass, velocity, pressure. The grid (if any) moves with the fluid. Advection is trivially accurate; interpolation is the hard part.
Eulerian: Fix the grid and watch fluid properties flow through it. Pressure projection and divergence-free velocity fields are solved per-cell. Good for large-scale steady flows; advection requires numerical diffusion management.
Hybrid (FLIP, APIC, MPM): Particles carry state; a background grid accumulates and projects forces. Combines advantages of both at the cost of transfer overhead.

2. SPH — Smoothed Particle Hydrodynamics

SPH represents fluid as a set of particles. Each property (density, velocity, pressure) is interpolated from neighbours using a smoothing kernel W(r, h). Originally developed for astrophysical gas dynamics (Gingold & Monaghan 1977), it is now ubiquitous in real-time water and lava simulations.

SPH density estimate at particle i: ρᵢ = Σⱼ mⱼ \cdot W(|rᵢ - rⱼ|, h) SPH pressure force: Fᵢ_pressure = -mᵢ Σⱼ mⱼ (Pᵢ/ρᵢ² + Pⱼ/ρⱼ²) \cdot \nablaW(rᵢⱼ, h) Commonly: W = Poly6 kernel for density, Spiky kernel for pressure gradient (Müller et al. 2003 "Particle-Based Fluid Simulation for Interactive Applications")

Strengths

Free surfaces handled naturally — no tracking of water surface needed
Easy to implement from scratch (weekend project)
Naturally handles large deformations — no mesh tangling
GPU-friendly with spatial hash acceleration

Weaknesses

Incompressibility problem: WCSPH uses a stiff equation of state (high sound speed), requiring very small timesteps. PCISPH / DFSPH (2019) solve pressure iteratively each step for near-incompressible results at larger dt.
Noisy density estimates at particle-sparse regions
O(N log N) per step with spatial hashing; O(N²) naive

Use SPH when: you need free surfaces, splashing, liquid-solid coupling, and can tolerate slightly compressible artefacts. Used in: games (Nvidia FleX, position-based fluids), astrophysics, ocean simulation.

3. LBM — Lattice Boltzmann Method

LBM works at the mesoscopic scale — between molecular dynamics and Navier-Stokes. Fluid is modelled as particle probability density functions f_i streaming between fixed grid cells and colliding locally. The macroscopic fluid variables emerge from moments of f.

LBM update (BGK collision): fᵢ(x + cᵢΔt, t + Δt) = fᵢ(x, t) - (1/τ)(fᵢ - fᵢᵉq) τ = relaxation time (controls viscosity: ν = (τ - 0.5)/3 in lattice units) fᵢᵉq = Maxwell-Boltzmann equilibrium distribution D2Q9 lattice: 9 velocity directions in 2D (8 diagonal/cardinal + rest) D3Q19 lattice: 19 velocity directions in 3D

Strengths

Highly parallelisable — each cell updates independently from neighbours; ideal for GPU
Handles complex geometry without mesh generation (Bounce-back BC)
Straightforward implementation of multi-phase flows, porous media
Used in industrial CFD for high-Re internal flows (PowerFLOW by Exa/3DS)

Weaknesses

Limited to incompressible, low-Mach flows (Ma < ~0.3) in standard formulation
Free surfaces require complex Volume-of-Fluid or free-surface LBM extensions
Grid stores 19–27 floats per cell (D3Q27) — memory-intensive at high resolution
Viscosity tuning via τ can become numerically unstable for τ very close to 0.5

Use LBM when: you need high-Re internal flows (ducts, around aerodynamic shapes), porous media flow, or want maximum GPU utilisation. The lattice-boltzmann simulation on this site uses D2Q9.

4. MPM — Material Point Method

MPM, developed by Sulsky et al. (1994) and popularised by Disney Research for Frozen (2013), is a hybrid Lagrangian-Eulerian method. Particles carry mass and deformation gradient; a background grid accumulates forces, solves momentum, and transfers velocities back to particles.

MPM loop per timestep: 1. P2G: Scatter particle mass/momentum to grid nodes (bilinear interpolation) 2. Grid update: Apply forces, boundary conditions, integrate velocities 3. G2P: Gather new velocities back to particles (APIC / FLIP transfer) 4. Advect: Move particles: x_p += v_p \cdot dt 5. Update deformation gradient F via velocity gradient from grid

Strengths

Handles snow, sand, slime, paste, cloth, and fluid with a single codebase — just change the constitutive model (elastoplasticity, viscoplasticity, neo-Hookean)
No mesh tangling — grid is fixed background structure
Used in major VFX for snow (Frozen), mud (Brave), sand (Moana)

Weaknesses

P2G/G2P transfers introduce numerical diffusion; APIC (Jiang 2015) and PolyPIC improve this significantly
Grid cell size limits resolution; adaptive/nested grids are complex
More implementation complexity than SPH

5. FVM / FEM — Grid-Based (Navier-Stokes direct)

Classic Eulerian CFD discretises the Navier-Stokes equations directly on a mesh using Finite Volume (FVM) or Finite Element (FEM) methods. The velocity-pressure coupling is enforced via a Poisson equation solved each timestep (Chorin's projection method).

Incompressible Navier-Stokes: \partialu/\partialt + (u\cdot\nabla)u = -(1/ρ)\nablaP + ν\nabla²u + f (momentum) \nabla\cdotu = 0 (incompressibility) Chorin projection: 1. Advect and diffuse: u* = u + dt(-(u\cdot\nabla)u + ν\nabla²u + f) 2. Solve pressure Poisson: \nabla²P = (ρ/dt)\nabla\cdotu* 3. Project: u^(n+1) = u* - (dt/ρ)\nablaP

Strengths

Most physically accurate for well-resolved flows
Rich literature, mature solvers (OpenFOAM, Fluent, SU2)
Conserves mass, momentum, energy exactly (when formulated correctly)

Weaknesses

Free surfaces require level set or VOF tracking — significant extra coding
Mesh generation for complex geometry is a major engineering effort
Pressure solve scales poorly with resolution (sparse linear system solver)

6. Comparison Table

Method	Family	Free Surfaces	Parallelism	Incompressibility	Best For
SPH / DFSPH	Lagrangian	✅ Natural	Good (GPU hash)	⚠️ Iterative	Splashing water, lava, astrophysics
LBM	Mesoscopic grid	⚠️ Extension needed	✅ Excellent (1 cell = 1 thread)	✅ Built-in	Aerodynamics, porous media, industrial CFD
MPM / APIC	Hybrid	✅ Natural	Good	Depends on constitutive law	Snow, sand, mud, multi-material VFX
FVM (projection)	Eulerian grid	⚠️ Level set needed	Moderate (sparse solve)	✅ Projection enforced	Engineering CFD, steady flows
FLIP / PIC	Hybrid	✅ Natural	Good	Grid projection	Large-scale liquids (Houdini FLIP)

7. How to Choose

Real-time game fluid (water, ~1K–100K particles): SPH or Position-Based Fluids (PBF). Fast to implement, GPU-friendly, visually convincing free surfaces.
High-Re aerodynamics simulation: LBM. PowerFLOW/XFlow use LBM for car/aircraft simulations up to Re = 10⁷ on GPU clusters.
Snow/sand/paste in VFX or games: MPM. Disney's Snow model (Stomakhin 2013) is open-sourced as reference.
Scientific accuracy, turbulence: FVM with RANS/LES (OpenFOAM). High-Re flows require turbulence models that SPH and LBM can approximate but not match.
Smoke/fire in real-time: Grid-based Eulerian (Jos Stam 1999 "Stable Fluids"). Semi-Lagrangian advection avoids CFL timestep restriction; industry-standard approach.