Lagrangian Fluids

Smoothed-particle hydrodynamics from first principles. From my Computational Physics final project.

Please Scroll

move your cursor to stir the fluid

The Lagrangian View

The Navier–Stokes momentum equation, written at fixed points in space (the Eulerian frame):

\[ \frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v}\cdot\nabla)\mathbf{v} = -\frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{v} + \mathbf{g} \]

The advection term \((\mathbf{v}\cdot\nabla)\mathbf{v}\) is nonlinear and is the main source of difficulty. The Lagrangian formulation discretizes the fluid into parcels and follows them. In a parcel's frame the material derivative absorbs advection:

\[ \frac{D\mathbf{v}}{Dt} = -\frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{v} + \mathbf{g} \]

Advection is handled by the particles moving. Mass conservation is automatic, since the particles carry the mass. The faint grid is the Eulerian mesh this method does not use. Each streak is one parcel's velocity.

Smoothed Particles

Each particle is a sample point carrying mass and velocity. Continuous fields are reconstructed by kernel-weighted sums over neighbors:

\[ A(\mathbf{r}) \;\approx\; \sum_j \frac{m_j}{\rho_j}\, A_j\, W(\lvert\mathbf{r}-\mathbf{r}_j\rvert,\, h) \]

Setting \(A = \rho\) gives the density estimate, computed first every step:

\[ \rho_i = \sum_j m_j\, W(\lvert\mathbf{r}_i-\mathbf{r}_j\rvert,\, h), \qquad W_{\text{poly6}} = \frac{4}{\pi h^8}\,(h^2 - r^2)^3 \]

The smoothing length \(h\) sets the interaction range. The rings show it on three tracked particles: neighbors inside the ring contribute to that particle's density, neighbors outside do not. Color encodes density.

Pressure

Exact incompressibility requires a Poisson solve every step. Weakly compressible SPH (Müller et al., 2003) uses a stiff equation of state instead:

\[ p = k\,(\rho - \rho_0), \qquad p \ge 0 \]

Compressed regions develop pressure and push back. The clamp at zero avoids the tensile instability. The force is symmetrized so momentum is conserved pairwise:

\[ \mathbf{f}_i^{\,\text{press}} = -\sum_j m_j\, \frac{p_i + p_j}{2\rho_j}\, \nabla W_{\text{spiky}} \]

The gradient uses the spiky kernel rather than poly6: the poly6 gradient vanishes as \(r \to 0\), so close pairs would feel no repulsion and clump. Color encodes pressure. It peaks where the fluid meets the walls.

Viscosity

Viscosity diffuses momentum between neighboring parcels:

\[ \mathbf{f}_i^{\,\text{visc}} = \mu \sum_j m_j\, \frac{\mathbf{v}_j - \mathbf{v}_i}{\rho_j}\, \nabla^2 W_{\text{visc}}, \qquad \nabla^2 W = \frac{40}{\pi h^5}(h - r) \]

The Laplacian of this kernel is positive everywhere, so the force only damps relative motion and never adds energy (Müller, 2003). It also stabilizes the explicit integration.

Here \(\mu\) is raised three-fold. The same dam break becomes a slow collapse.

Finding Neighbors

The kernels have compact support, so every SPH sum runs over neighbors within \(h\). Testing all pairs costs \(\mathcal{O}(N^2)\). A uniform grid with cell size \(h\) reduces it:

\[ \text{bin all particles} \;\rightarrow\; \text{scan } 3\times3 \text{ cells} \;\Rightarrow\; \mathcal{O}(N) \]

Particles are binned with a counting sort each substep. The outlined cells are the occupied bins; the highlighted square is one particle's \(3\times3\) search neighborhood, with its kernel radius drawn inside.

The same loop ported to the GPU (Taichi, one thread per particle) ran 10,000+ particles at 60 fps, since each particle only reads its own neighborhood. Stability requires \(\Delta t \cdot v_{\max} < h\), so a particle cannot cross a cell in one step.

Zero Gravity

Gravity off. A weak cohesive attraction between neighbors stands in for surface tension, and the fluid is given an initial spin. Pressure pushes out, cohesion pulls in, and the balance is a rotating drop.

The full project ran these experiments in 2D and 3D with Taichi: spinning drops with spin-down measured against viscosity, drop oscillation under a tension force, and collisions of free-floating blobs.

This page runs the same algorithm in plain JavaScript: density → pressure → forces → integration, four substeps per frame, about two thousand particles.