Fluid Simulation

8 Nov 2022

Notes Computer Animation Fluid Simulation

This was a guest lecture by Dr. John Keyser

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(∥x⃗i−x⃗j∥),W(r) = W(\|\vec{x}_i - \vec{x}_j\|), typically called a kernel for SPH.
    • Characteristics of W:W:
      • Radial symmetry
      • Finite support
      • Integrated volume
      • Monotonically decreasing, with max at r=0r = 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 φi\varphi_i be some characteristic measured at particle i.i. Then,

φ(x⃗)=∑iφimiW(∥x⃗−x⃗i∥)∑imiW(∥x⃗−x⃗i∥)\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 mim_i is the mass of the particle, typically mi=1m_i = 1 everywhere.

The denominator of the above is the particle density,

ρi=∑jmjW(∥x⃗i−x⃗j∥).\rho_i = \sum_j m_j W(\|\vec{x}_i - \vec{x}_j\|).

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

φ(x⃗)=∑iφimiρ(x⃗)W(∥x⃗−x⃗i∥)\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,

∇φ(x⃗)=∇∑iφimiρ(x⃗)W(∥x⃗−x⃗i∥)=∑iφimiρ(x⃗)∇W(∥x⃗−x⃗i∥)∇2φ(x⃗)=∑iφimiρ(x⃗)∇2W(∥x⃗−x⃗i∥)\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=ma⃗.F = m\vec{a}.

In fluids, instead of mass we have density

F=ρa⃗  ⟺  a⃗=1ρF.F = \rho\vec{a} \iff \vec{a} = \frac{1}{\rho} F.

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

∇pjρpj=∑imi(piρi2+pjρj2)∇W(∥x⃗i−x⃗j∥).\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\|).

p=kρ0γ[(ρρ0)γ−1]pi=k(ρi−ρ0)\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 ρ0\rho_0 is rest density, kk is a stiffness constant, and usually γ=1.\gamma = 1.

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

§ Diffusion and Viscosity

Fviscosity=ν∇2v⃗=ν∇⋅∇v⃗=ν div(∇v⃗)=∑jmiv⃗i−v⃗jρj∇2W(∥x⃗i−x⃗j∥).\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
    • u⃗(n)←\vec{u}^{(n)} \leftarrow velocity at time n.n.
    • w⃗0=u⃗(n)\vec{w}_0 = \vec{u}^{(n)}
    •  
    • w⃗1=w⃗0+ΔtF⃗ρ\vec{w}_1 = \vec{w}_0 + \Delta t \frac{\vec{F}}{\rho}
    •  
    • w⃗2=backtrace(w⃗1,Δt)\vec{w}_2 = \verb|backtrace|(\vec{w}_1, \Delta t)
    • w^2=FFT(w⃗2)\hat{w}_2 = \mathrm{FFT}(\vec{w}_2)
    •  
    • w^3=w^2(I−νΔt∇2)−1\hat{w}_3 = \hat{w}_2 \left(\mathbf{I} - \nu \Delta t \nabla^2\right)^{-1} [a]
    •  
    • w^4=w^3−∇φ\hat{w}_4 = \hat{w}_3 - \nabla \varphi
    • w⃗4=FFT−1(w^4)\vec{w}_4 = \mathrm{FFT}^{-1}(\hat{w}_4) [b]

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:

∂w⃗2∂t=ν∇2w⃗2.\frac{\partial \vec{w}_2}{\partial t} = \nu \nabla^2 \vec{w}_2.

By recognizing the fact that we can discretize the Laplacian ∇2\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 w⃗\vec{w} can be split intop the sum of a divergence-free vector field u⃗\vec{u} and a curl-free vector field v⃗.\vec{v}.

w⃗=u⃗+v⃗,∇⋅u⃗=0,∇×v⃗=0.\vec{w} = \vec{u} + \vec{v}, \nabla \cdot \vec{u} = \mathbf{0}, \nabla \times \vec{v} = \mathbf{0}.

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

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

Following Helmholtz-Hodge:

u⃗=w⃗−v⃗=w⃗−∇φ∇⋅w⃗=∇⋅u⃗+∇⋅v⃗=∇v⃗=∇2φ\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 w⃗3\vec{w}_3, we get that ∇⋅w⃗3=∇2φ\nabla \cdot \vec{w}_3 = \nabla^2 \varphi in which we solve for φ\varphi.