Particle Fluid Methods

Particle in Cell Method

Particle in Cell (PIC) is a numerical method that moves away from the grid-based method into a hybrid of Lagrangian and Eulerian solvers.

• Grids (Eulerian): Excellent for solving the global pressure projection constraint because the math sits on a rigid, predictable matrix.
  • However, they are terrible at advection because interpolating data across the grid causes massive numerical diffusion (blurring).
• Particles (Lagrangian): Perfect for advection. Simply move the particle through space ($\frac{\partial \varphi }{\partial t}=\mathbf{v}$). There is exactly zero numerical diffusion.
  • However, calculating the global pressure constraint is a nightmare because particles are disorganized and require expensive nearest-neighbor searches to figure out who is pushing whom.

First, discretize the material space $M$ as particles $M=\{1,\dots ,\text{\# particles}\}$. Each particle $p\in M$ has a mass $m_p$ that is invariant over time. The particle location is our flow map $\varphi :M\to W={\mathbb{R}}^3$. Computing pressure projection is difficult using this particle method.

Denote an Eulerian function and its Lagrangian twin:

$$\begin{array}{rl}f& :W\to \mathbb{R}\\ \hat{f}& :M\to \mathbb{R}\end{array}$$

where

$$\hat{f}={\varphi }^{\ast }f=f\circ \varphi$$

To do advection, we need to move the particles according to the fluid velocity field. This is done by doing it in the Lagrangian coordinate:

$$\begin{array}{rl}\frac{\partial }{\partial t}\hat{f}& =0\\ \frac{\partial }{\partial t}\varphi & =\mathbf{u}\circ \varphi \end{array}$$

This is much easier as there is no diffusion term. This is unlike the Eulerian method where

$$\frac{\partial }{\partial t}f+{\mathcal{L}}_{\mathbf{u}}f=0$$

which requires heavy interpolation when discretized. This causes numerical diffusion and blurring.

Steps

At the beginning of each iteration, each particle $p\in M$ has a velocity ${\mathbf{v}}_p$ and ${\varphi }_p$.

1. Advection step $\frac{\partial }{\partial t}{\varphi }_p={\mathbf{v}}_p$ or $\mathbf{u}=\text{P2G}(\mathbf{v})$ and $\frac{\partial }{\partial t}\varphi =\mathbf{u}\circ \varphi$.
2. Transfer particle data to grid $\mathbf{u}=\text{P2G}(\mathbf{v})$.
3. Pressure projection on grid $\mathbf{u}=\text{PressureSolve}(\mathbf{u})$. Recall we do this via Applying FFT (Spectral Method).
4. Transfer grid data back to particles $\mathbf{v}=\text{G2P}(\mathbf{u})$.
5. End iteration.

Fluid Implicit Particle Solver (FLIP)

At the beginning of each iteration, each particle $p\in M$ has a velocity ${\mathbf{v}}_p$ and ${\varphi }_p$.

1. Advection step $\frac{\partial }{\partial t}{\varphi }_p={\mathbf{v}}_p$ or $\mathbf{u}=\text{P2G}(\mathbf{v})$ and $\frac{\partial }{\partial t}\varphi =\mathbf{u}\circ \varphi$.
2. Transfer particle data to grid $\mathbf{u}=\text{P2G}(\mathbf{v})$.
3. Pressure projection on grid ${\mathbf{u}}^{\prime }=\text{PressureSolve}(\mathbf{u})$. Recall we do this via Applying FFT (Spectral Method).
4. Transfer the incremental effect of pressure back to particle. $\mathbf{v}\leftarrow \mathbf{v}+\text{G2P}({\mathbf{u}}^{\prime }-\mathbf{u})$.
5. End iteration.

This works because the particle keeps its own original, perfectly advected velocity, and only adds the acceleration caused by the pressure gradient.

Particle-Grid Transfer

How do we calculate the transfer of data between continuous floating particles and a rigid, discrete grid? In particular, how do we calculate $\text{P2G}(\cdot )$ and $\text{G2P}(\cdot )$?

Before we can transfer data, we need a rule for how much a particle affects a grid cell. This is the basis function $N_i(\mathbf{x})$ for $i\in \text{Grid}$ and $\mathbf{x}\in {\mathbb{R}}^3$. It has the following properties:

• Local: $N_i(\mathbf{x})=0$ when $\mathbf{x}$ is far away from the grid point $i$.
• Partition of unity: ${\sum }_{i\in \text{Grid}}{N}_i(\mathbf{x})=1$ for all $\mathbf{x}$.
• Nonnegativity: $N_i(\mathbf{x})\ge 0$ for all $i,\mathbf{x}$.

Idea

Imagine a particle glowing like a lightbulb. $N_i(\mathbf{x})$ describes how bright that light is at grid cell $i$ when the particle is at position $\mathbf{x}$.

• Local: The light fades to exactly $0$ if the particle is too far away.
• Partition of unity: If you add up the light hitting all the surrounding grid cells, it must equal exactly $100\%$. You cannot artificially create or destroy mass/momentum during the transfer.
  • Usually, this is calculated using a B-spline, which is just a mathematically smooth bell curve.

Thus, we can calculate Grid to Particle ($\text{G2P}(\cdot )$) as

$${\hat{f}}_p\approx \sum _{i\in \text{Grid}}{f}_i{N}_i({\varphi }_p)$$

To find the particle’s new velocity, we look at the surround grid cells and multiply each grid velocity by the basis weight, and add them up (kind of like a 2D convolution!).

And for Particle to Grid ($\text{P2G}(\cdot )$):

$$f_i\approx \frac{{\sum }_{p\in M}{\hat{f}}_p{m}_p{N}_i({\varphi }_p)}{{\sum }_{p\in M}{m}_p{N}_i({\varphi }_p)}$$

This is just the weighted average of the particles, which causes the smoothing effect.

Affine Particle-Grid Transfer

We can upgrade both $\text{G2P}(\cdot )$ and $\text{P2G}(\cdot )$ to be Affine and not just linear. The idea is that each particle stores more information than just velocity.

Affine Grid to Particle ($\text{AffineG2P}(\cdot )$):

$$\begin{array}{rl}{\hat{f}}_p& \approx \sum _{i\in \text{Grid}}{f}_i{N}_i({\varphi }_p)\\ \widehat{\nabla f_p}& \approx \underset{\mathbf{A}}{\mathrm{argmin}}\sum _{i\in \text{Grid}}{\left|\frac{{\sum }_{p\in M}({\hat{f}}_p+\mathbf{A}[[{\mathbf{x}}_i-{\varphi }_p]])m_p{N}_i({\varphi }_p)}{{\sum }_{p\in M}{m}_p{N}_i({\varphi }_p)}-f_i\right|}^2\end{array}$$

Affine Particle to Grid ($\text{AffineP2G}(\cdot )$):

$$f_i\approx \frac{{\sum }_{p\in M}({\hat{f}}_p+\widehat{\nabla f_p}[[{\mathbf{x}}_i-{\varphi }_p]])m_p{N}_i({\varphi }_p)}{{\sum }_{p\in M}{m}_p{N}_i({\varphi }_p)}$$

The primary difference here is that the numerator now includes a term for the gradient of the particle’s properties. This is merely a first-order Taylor expansion of the particle’s influence on the grid.

Idea

Instead of a particle being a single point moving through space with velocity $\mathbf{v}$, we can imagine the particle is a tiny, squishy rubber ball. It can also stretch, shear, and spin. We store this inside $\widehat{\nabla f_p}$.

So how do we calculate the optimization problem? The goal is to find a specific $\mathbf{A}=\widehat{\nabla f_p}$ gradient matrix that captures how the particle’s properties are changing in space.

• Let ${\mathbf{r}}_i=({\mathbf{x}}_i-{\varphi }_p)$ be the distance vector from the particle to grid node $i$.
• Our approximation of the velocity at node $i$ is ${\hat{f}}_p+\mathbf{A}{\mathbf{r}}_i$. The actual velocity at node $i$ is $f_i$.
• We also only care about the nodes that are close to the particle.

Thus, we get an error function for a given $\mathbf{A}$:

$$E(\mathbf{A})=\sum _{i\in \text{Grid}}{N}_i({\varphi }_p){\Vert ({\hat{f}}_p+\mathbf{A}{\mathbf{r}}_i)-f_i\Vert }^2$$

We want to minimize this error. So we need to take the derivative and set it to $0$ to find the Local Minima.

$$\frac{\partial E}{\partial \mathbf{A}}=2\sum _{i\in \text{Grid}}{N}_i({\varphi }_p)\left({\hat{f}}_p+\mathbf{A}{\mathbf{r}}_i-f_i\right){\mathbf{r}}_i^{\top }=0$$

With some rearranging, we get

$$\sum _i{N}_i({\varphi }_p)\mathbf{A}{\mathbf{r}}_i{\mathbf{r}}_i^{\top }=\sum _i{N}_i({\varphi }_p)(f_i-{\hat{f}}_p){\mathbf{r}}_i^{\top }$$

For this specific particle, $\mathbf{A}$ is a constant matrix, so we can pull it out of the summation.

$$\mathbf{A}\left[\sum _i{N}_i({\varphi }_p){\mathbf{r}}_i{\mathbf{r}}_i^{\top }\right]=\sum _i{N}_i{f}_i{\mathbf{r}}_i^{\top }-{\hat{f}}_p\left[\sum _i{N}_i({\varphi }_p){\mathbf{r}}_i^{\top }\right]$$

The very last term is important. Because standard B-spline basis functions are perfectly symmetrical, if you sum up all the distances to the surrounding grid nodes weighted by the spline, they perfectly balance out to exactly $0$. This causes ${\hat{f}}_p$ to completely vanish! We are left with

$$\mathbf{A}\left[\sum _i{N}_i({\varphi }_p){\mathbf{r}}_i{\mathbf{r}}_i^{\top }\right]=\sum _i{N}_i({\varphi }_p)f_i{\mathbf{r}}_i^{\top }$$

Now, define

$$\begin{array}{rl}\mathbf{B}& :=\sum _{i\in \text{Grid}}{N}_i({\varphi }_p)f_i({\mathbf{x}}_i-{\varphi }_p{)}^{\top }\\ \mathbf{D}& :=\sum _{i\in \text{Grid}}{N}_i({\varphi }_p)({\mathbf{x}}_i-{\varphi }_p)({\mathbf{x}}_i-{\varphi }_p{)}^{\top }\end{array}$$

where $\mathbf{D}$ is the “Inertia” tensor (note how $N_i$ is like a mass distribution and ${\mathbf{r}}_i{\mathbf{r}}_i^{\top }$ is like the moment of inertia) and $\mathbf{B}$ is the “Angular Momentum” tensor (note how $N_i$ is like a mass distribution, $f_i$ is like the velocity, and ${\mathbf{r}}_i$ is like the distance from the center of mass).

Thus, $\mathbf{AD}=\mathbf{B}$ and $\widehat{\nabla f_p}=\mathbf{B}{\mathbf{D}}^{-1}$.

Lattice Boltzmann Method (LBM)

Previously, we have been solving the Eulerian fluid equations by directly discretizing the PDEs. The Lattice Boltzmann Method (LBM) simulates a mesoscopic¹ statistical view of the fluid.

1. Continuum Mechanics is only a macroscopic model where the fluid is a continuous blob. Every point in space has one single, averaged velocity $\mathbf{u}$.
2. Molecular dynamics simulates every single fluid molecule bouncing off each other. One $\text{mol}$ of water is $10^{23}$ particles (computationally expensive!)
3. Mesoscopic (statistical) dynamics tracks the probability of particles moving in certain directions at a given point in space instead of tracking one average velocity or many.

The state of the system is a Probability Measure $f(\mathbf{x},\mathbf{v},t)$ over the tangent bundle $TW$. In particular, $f(\mathbf{x},\mathbf{v},t)$ is the probability density of a particle being at position $\mathbf{x}$ and velocity $\mathbf{v}$ at time $t$. With this, we can calculate:

• The mass density: $${\rho }_{(\mathbf{x},t)}={\int }_{v\in T_{\mathbf{x}}M}f(\mathbf{x},\mathbf{v},t)$$
• The momentum density: $$(\rho \mathbf{u}{)}_{(\mathbf{x},t)}={\int }_{v\in T_{\mathbf{x}}W}\mathbf{v}f(\mathbf{x},\mathbf{v},t)$$
• Internal energy: $$e_{(\mathbf{x},t)}={\int }_{v\in T_{\mathbf{x}}M}\frac{1}{2}\Vert \mathbf{v}-\mathbf{u}{\Vert }^{2}f(\mathbf{x},\mathbf{v},t)=\frac{3k{T}_{\mathbf{x},t}}{2}$$ where $T$ is the temperature and $k$ is the Boltzmann constant.
• Stress-energy tensor: $${\Sigma }_{(\mathbf{x},t)}={\int }_{v\in T_{\mathbf{x}}M}(\mathbf{v}-\mathbf{u})(\mathbf{v}-\mathbf{u}{)}^{\top }f(\mathbf{x},\mathbf{v},t)$$

Boltzmann Equation

The Boltzmann equation describes how the probability density $f$ evolves over time.

$$\underset{\text{advection}}{\underbrace{\frac{\partial }{\partial t}f(\mathbf{x},\mathbf{v},t)+\mathbf{v}\cdot f(\mathbf{x},\mathbf{v},t)+{\mathbf{g}}^{\text{ext}}\cdot \frac{\partial }{\partial \mathbf{v}}f(\mathbf{x},\mathbf{v},t)}}=Q(f{)}_{(\mathbf{x},\mathbf{v},t)}$$

The probability distribution is advected by its own velocity and external forces ${\mathbf{g}}^{\text{ext}}$. The right hand side “can be”² a pointwise operator $Q$ modeling the internal collision that mutates the distribution (obeying some conservation law).

Instead of particles moving in any direction, they are only allowed to move in a handful of fixed directions. This dramatically simplifies the computation required.

• $f$ becomes discrete
• Advection becomes a simple shift of the probability distribution along the lattice directions.
• Local collision runs $Q(f)$ on each node.

$"\\usepackage{tikz}\n\n\\begin{document}\n\\begin{tikzpicture}\n\n$

source code

Other Methods

• Vorticity is curl of velocity. Vortex lines are curves parallel to vorticity.
• Vortex lines are transported in the Euler equation.
• One can reconstruct divergence-free velocity from the vorticity field (on simply-connected domains).
• Simulations of fluid using this principle are called vortex methods.
• When velocity is treated as 1-form (covector), it is also known as “impulse”. They are Lie-transported in Euler fluids.
• They correspond to covector methods/impulse methods.

Footnotes

1. Mesoscopic means it is between the microscopic particle level and the macroscopic fluid level. ↩

2. The reason why it is “can be” is because it depends on the kind of fluid. For example, we could have a dilute fluid (like upper atmosphere) where the particles are so far apart that they almost never collide. Since collisions are statistically negligible, we can set $Q=0$. Dense fluids will have a more complicated $Q$ that models the collisions. ↩