To simulate fluids, we can simulate a velocity field which is applied to liquid or suspended particles.

There are two types of simulation:

  • Lagrangian, where material is discretized as in smoothed-particle hydrodynamics (SPH).
  • Eulerian, where space is discretized as in stable fluids.

Smoothed-Particle Hydrodynamics§

See: review of Smoothed particle hydrodynamics (Monaghan, 2005).

  • Use particles to represent fluid.
    • Particle motion is the same as material transport.
  • The particles hold information.
    • The information is accurate at the particle’s position, however, it falls off with distance
  • The fall-off can be given by a weight function $W(r) = W(\|\vec{x}_i - \vec{x}_j\|),$ typically called a kernel for SPH.
    • Characteristics of $W:$
      • Symmetric, radially
      • Finite support
      • Integrated volume
      • Monotonically decreasing, with max at $r = 0$
      • Smooth gradient at boundary
      • See: poly6 kernel

If particles have some characteristic that they carry (e.g. velocity), we can compute a weighted average of influence to determine the characteristic at any point in space.

Let $\varphi_i$ be some characteristic measured at particle $i.$ Then,

$$ \varphi(\vec{x}) = \frac{\sum_{i} \varphi_i m_i W\left(\|\vec{x} - \vec{x}_i\|\right)}{\sum_{i} m_i W\left(\|\vec{x} - \vec{x}_i\|\right)} $$

where $m_i$ is the mass of the particle, typically $m_i = 1$ everywhere.

The denominator of the above is the particle density,

$$ \rho_i = \sum_j m_j W(\|\vec{x}_i - \vec{x}_j\|). $$

Then, a characteristic can be computed for any $\vec{x} \in \mathbb{R}^3$, not just particle position $\vec{x}_k:$

$$ \varphi(\vec{x}) = \sum_i \varphi_i \frac{m_i}{\rho(\vec{x})} W(\|\vec{x} - \vec{x}_i\|) $$

This also makes it easy to compute the derivatives of characteristics,

$$ \begin{align*} \nabla\varphi(\vec{x}) &= \nabla \sum_i \varphi_i \frac{m_i}{\rho(\vec{x})} W(\|\vec{x} - \vec{x}_i\|)\\ &= \sum_i \varphi_i \frac{m_i}{\rho(\vec{x})} \nabla W(\|\vec{x} - \vec{x}_i\|)\\ \nabla^2 \varphi(\vec{x}) &= \sum_i \varphi_i \frac{m_i}{\rho(\vec{x})} \nabla^2 W(\|\vec{x} - \vec{x}_i\|) \end{align*} $$

Fluid Forces§

Force is defined by $$ F = m\vec{a}. $$

In fluids, instead of mass we have density $$ F = \rho\vec{a} \iff \vec{a} = \frac{1}{\rho} F. $$

Then, we can define the pressure force as $-\nabla p$:

$$ \frac{\nabla p_j}{\rho p_j} = \sum_i m_i \left(\frac{p_i}{\rho^2_i} + \frac{p_j}{\rho^2_j}\right) \nabla W(\|\vec{x}_i - \vec{x}_j\|). $$

$$ \begin{align*} p &= \frac{k \rho_0}{\gamma} \left[\left(\frac{\rho}{\rho_0}\right)^{\gamma} - 1\right]\\ p_i &= k (\rho_i - \rho_0) \end{align*} $$

where $\rho_0$ is rest density, $k$ is a stiffness constant, and usually $\gamma = 1.$

We can compute $p_i$ for each particle, then compute the pressure force from that.

Diffusion and Viscosity§

$$ \begin{align*} F^{\textrm{viscosity}} &= \nu \nabla^2 \vec{v} = \nu \nabla \cdot \nabla \vec{v} = \nu\,\textrm{div}(\nabla \vec{v})\\ &= \sum_j m_i \frac{\vec{v}_i - \vec{v}_j}{\rho_j} \nabla^2 W(\|\vec{x}_i - \vec{x}_j\|). \end{align*} $$

where $\nu$ is the kinematic viscosity.

Other forces can also be applied such as gravity or surface tension.

Stable Fluids (1999)§

See: Stable Fluids (Stam, 1999).

The the algorithm is as follows:

Stable Fluids

while simulation loop do:   $\vec{u}^{(n)} \leftarrow$ velocity at time $n.$   $\vec{w}_0 = \vec{u}^{(n)}$   $\vec{w}_1 = \vec{w}_0 + \Delta t \frac{\vec{F}}{\rho}$   $\vec{w}_2 = \verb|backtrace|(\vec{w}_1, \Delta t)$   $\hat{w}_2 = \mathrm{FFT}(\vec{w}_2)$   $\hat{w}_3 = \hat{w}_2 \left(\mathbf{I} - \nu \Delta t \nabla^2\right)^{-1}$ [a]   $\hat{w}_4 = \hat{w}_3 - \nabla \varphi$ [b]   $\vec{w}_4 = \mathrm{FFT}^{-1}(\hat{w}_4)$

Backtrace

In each step, want to know what the new velocity is at a point based on the prior velocity.

In other words, we want to advect the fluid’s velocity.

If we simply move the velocity forward along its direction, we may end up with more velocity (energy) than we started with.

A simple solution is to use a ‘semi-Lagrangian’ method trace a “particle” backward from the point that we want to know.

[a]: This step applies the fluid’s viscosity, as noted in the paper it is equivalent to the following diffusion equation:

$$ \frac{\partial \vec{w}_2}{\partial t} = \nu \nabla^2 \vec{w}_2. $$

By recognizing the fact that we can discretize the Laplacian $\nabla^2$ (because we have discretized the space), we can apply the Fourier transform to convolve the space with the Laplacian kernel then do the inverse Fourier to recover the field.

[b]: This step is about projecting the field onto a divergence-free field.

The scalar field $\varphi$ is a result of the Helmholtz-Hodge decomposition.

Helmholtz-Hodge Decomposition

Any vector field $\vec{w}$ can be split intop the sum of a divergence-free vector field $\vec{u}$ and a curl-free vector field $\vec{v}.$

$$ \vec{w} = \vec{u} + \vec{v}, \nabla \cdot \vec{u} = \mathbf{0}, \nabla \times \vec{v} = \mathbf{0}. $$

Divergence free: $\vec{u}$ is incompressible (no sources/sinks)

Curl free $\iff \vec{v}$ is the gradient of some scalar field $\varphi,$ i.e. $\vec{v} = \nabla \varphi.$

Following Helmholtz-Hodge:

$$ \begin{align*} \vec{u} &= \vec{w} - \vec{v}\\ &= \vec{w} - \nabla \varphi\\ \nabla \cdot \vec{w} &= \nabla \cdot \vec{u} + \nabla \cdot \vec{v}\\ &= \nabla \vec{v}\\ &= \nabla^2 \varphi \end{align*} $$

Then, plugging in $\vec{w}_3$, we get that $\nabla \cdot \vec{w}_3 = \nabla^2 \varphi$ in which we solve for $\varphi$.