Grid-Based Fluid Methods

Marker-and-Cell (MAC) Grid

We want to be able to simulate Continuum Mechanics, and more specifically Fluid Mechanics.

A crude and naive way to do is to split the world into a 3D grid of cells, labeled whether it is an obstacle, fluid, or empty space. The fluid velocity is defined on the cell wall.

$"\\usepackage{tikz}\n\n\\begin{document}\n\\begin{tikzpicture}[scale=0.8, every node/.style={scale=1.0}]\n\n$

source code

The fluid velocity has 3 components

$$\mathbf{u}=(u^1,u^2,u^3)$$

like 3 scalar fields, with $u^1$ defined on the faces normal to the $x-$direction (and similarly with the other two). Velocity on a wall next to a solid cell must be $0$.

$"\\usepackage{amsmath}\n\\usepackage{amssymb}\n\\usepackage{amsfonts}\n\\usepackage{xcolor}\n\n\\begin{document}\n\\begin{tikzpicture}[>=stealth, line join=round, line cap=round, \n x={(1.3cm,0cm)}, y={(0cm,1.3cm)}, z={(0.6cm,0.45cm)}]\n\n$

source code

• Let $\mathcal{F}$ denote the set of faces incident to both fluid cells.
• Let $\mathcal{C}$ denote the set of fluid cells.
• Velocity field is a data $\mathbf{u}\in {\mathbb{R}}^{\mathcal{F}}$.
• Scalar quantities like pressure and density are stored at cell center.

The divergence of the velocity field is defined on fluid cells

$$q=\text{div }\mathbf{u}\in {\mathbb{R}}^{\mathcal{C}}$$

where $\text{div}:{\mathbb{R}}^{\mathcal{F}}\to {\mathbb{R}}^{\mathcal{C}}$ (a $|\mathcal{C}|\times |\mathcal{F}|$ sparse matrix) is defined by aggregating the total flux from each cell wall.

Pressure Projection

Given non divergence free $\tilde{\mathbf{u}}\in {\mathbb{R}}^{\mathcal{F}}$, (some velocity field generated without respect to pressure), we need to try and project it to $\text{ker}(\text{div})\subset {\mathbb{R}}^{\mathcal{F}}$. The idea here is that $\text{ker}(\text{div})$ consists of all velocity fields with exactly zero divergence. The goal of the algorithm is to find the closest valid vector in this subspace to our invalid input vector $\tilde{\mathbf{u}}$.

Idea

This is mathematically the least-squares projection problem. Again, treat $\text{div}$ as a linear operator (matrix).

We want to solve for pressure $p\in {\mathbb{R}}^{\mathcal{C}}$ where

$$\mathbf{u}=\tilde{\mathbf{u}}-{\text{div}}^{\top }p$$

and $\text{div }\mathbf{u}=0$ (to satisfy incompressibility). Mathematically this means we need to find a $p$ where

$$\begin{array}{rl}\mathbf{u}& =\tilde{\mathbf{u}}-{\text{div}}^{\top }p\\ \text{div }\mathbf{u}& =\text{div}(\tilde{\mathbf{u}}-{\text{div}}^{\top }p)\\ 0& =\text{div }\tilde{\mathbf{u}}-\text{div}{\text{div}}^{\top }p\\ \text{div}{\text{div}}^{\top }p& =\text{div }\tilde{\mathbf{u}}\end{array}$$

Note that ${\text{div}}^{\top }:{\mathbb{R}}^{\mathcal{C}}\to {\mathbb{R}}^{\mathcal{F}}$ is a discrete (negative) gradient operator. So, solving this system is the same as finding $Ax=b$ where

$$\underset{A}{\underbrace{\text{div}{\text{div}}^{\top }}}p=\underset{b}{\underbrace{\text{div }\tilde{\mathbf{u}}}}$$

and $x:=p$. We can also see that the matrix $\text{div}{\text{div}}^{\top }\in {\mathbb{R}}^{|\mathcal{C}|\times |\mathcal{C}|}$ is a discrete Laplacian. Fortunately, this means we can solve this by any linear algebra routine such as conjugate gradient.

Applying FFT (Spectral Method)

If the domain has a simple periodic boundary condition, then we can use the Fast Fourier Transform (FFT) to turn this system into a diagonalized equation with direct solution.

For a fluid to have a simple periodic boundary condition, it means a fluid flows out one side of some volume (like a cube) and flows in seamlessly on the left¹.

Let $(L_1,L_2,L_3)$ be the tuple that describes the “size” of the periodic domain. The goal is to solve the Poisson problem $\Delta \mathbf{u}=f$² such that

$$\frac{{\partial }^2\mathbf{u}}{\partial x^2}+\frac{{\partial }^2\mathbf{u}}{\partial y^2}+\frac{{\partial }^2\mathbf{u}}{\partial z^2}=f$$

² $\Delta$ is actually the Laplacian operator, the continuous equivalent of the discrete $\text{div}{\text{div}}^{\top }$ matrix!

We will discretize both $\mathbf{u}$ and $f$ on a regular grid with resolution $(N_1,N_2,N_3)$. Then apply FFT:

$$\begin{array}{rl}u(j_1,j_2,j_3)& =\frac{1}{N_1{N}_2{N}_3}\sum _{k_1=0}^{N_1-1}\sum _{k_2=0}^{N_2-1}\sum _{k_3=0}^{N_3-1}\hat{u}(k_1,k_2,k_3)\exp \left(2\pi i\left(\frac{j_1{k}_1}{N_1}+\frac{j_2{k}_2}{N_2}+\frac{j_3{k}_3}{N_3}\right)\right)\\ \hat{u}(k_1,k_2,k_3)& =\sum _{j_1=0}^{N_1-1}\sum _{j_2=0}^{N_2-1}\sum _{j_3=0}^{N_3-1}u(j_1,j_2,j_3)\exp \left(-2\pi i\left(\frac{j_1{k}_1}{N_1}+\frac{j_2{k}_2}{N_2}+\frac{j_3{k}_3}{N_3}\right)\right)\end{array}$$

where $j_1,j_2,j_3$ represent the individual grid cells and $k_1,k_2,k_3$ represent the individual frequencies. Recall that $\mathbf{u},f$ represent variables in the spatial domain and $\hat{\mathbf{u}},\hat{f}$ represent variables in the frequency domain.

This is fast because we can represent the pressure field as a linear combination of Fourier basis functions. More simply, image $\mathbf{u}$ is the sum of one wave, $\mathbf{u}=e^{ikx}$ where $x$ is the position in space and $k$ is the frequency of the wave. The derivative of $\mathbf{u}$ is $\frac{\partial }{\partial x}\mathbf{u}=ik{e}^{ikx}$, making this operation on a computer much easier.

Since the simulation exists on a discrete grid of fixed pixels, we need to prevent aliasing² via the following definition of $p_1$:

$$\widehat{\frac{\partial }{\partial x}\mathbf{u}(k)}=p_1\hat{\mathbf{u}}(k),p_1=\{\begin{array}{ll}2\pi i\frac{k_1}{L_1}& 0\le k_1<N_1/2\\ 0& k_1=N_1/2\\ 2\pi i\frac{k_1-N_1}{L_1}& N_1/2<k_1<N_1\end{array}$$

There are two distinct subspaces that our final velocity field $\mathbf{u}$ can exist in

1. The subspace of all possible velocity fields where the fluid is perfectly incompressible everywhere (i.e. $\text{div }\mathbf{u}=0$).
2. The subspace of all possible velocity fields where the fluid inside and immediately on the boundary of the solid obstacle is completely stationary (i.e. $\mathbf{u}=0$).

The solver finds the intersection of both. The alternating projection (pressure project $\to$ set velocity to zero by enforcing obstacles) will converge to their intersection.

Obstacles

If we want to add obstacles, we can just do Pressure Projection and set velocity in the obstacle to zero.

Time Splitting

Consider the incompressible fluid Euler equation:

$$\{\begin{array}{rl}\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}& =-\nabla p\\ \nabla \cdot \mathbf{u}& =0\end{array}$$

Projection Method (Chorin)

For each time step, we will

1. solve the advection term (for half timestep): $$\frac{\partial \mathbf{u}}{\partial t}+{\nabla }_{\mathbf{u}}\mathbf{u}=0$$
2. pressure “reflect” $$\mathbf{u}\leftarrow \mathbf{u}-\nabla {\Delta }^{-1}\nabla \cdot \mathbf{u}$$

Because we advected the fluid while ignoring pressure, the fluid followed the wrong trajectory for that entire $\Delta t$. Projecting it straight back down does not return it to where it would have been if pressure had been applied continuously.

Reflection Method (Strang)

We can solve this by: for each time step, we will

1. solve the advection term (for half timestep): $$\frac{\partial \mathbf{u}}{\partial t}+{\nabla }_{\mathbf{u}}\mathbf{u}=0$$
2. pressure “reflect” $$\begin{array}{rl}{\mathbf{u}}^{\ast }& =\mathbf{u}-\nabla {\Delta }^{-1}\nabla \cdot \mathbf{u}\\ \mathbf{u}& \leftarrow 2{\mathbf{u}}^{\ast }-\mathbf{u}\end{array}$$ This reflection across the div-free subspace allows you to push the fluid further along the correct trajectory (via the advection step) before projecting it back down to the div-free subspace.
3. solve the advection equation (for half timestep): $$\frac{\partial \mathbf{u}}{\partial t}+{\nabla }_{\mathbf{u}}\mathbf{u}=0$$
4. pressure project

$"\\usepackage{amsmath}\n\\usepackage{amssymb}\n\\usepackage{xcolor}\n\n\\begin{document}\n\\begin{tikzpicture}[>=stealth, scale=1.3]\n\n$

source code

Advection

Let $\mathbf{v}$ be a time-independent vector field over the world coordinate $W$. A time dependent scalar function $f$ is advected by $\mathbf{v}$ if

$$\frac{\partial }{\partial t}f+\mathbf{v}\cdot \nabla f=0\text{ or }\frac{\partial }{\partial t}f+df[[\mathbf{v}]]=0\text{ or }\frac{\partial }{\partial t}f+{\mathcal{L}}_{\mathbf{v}}f=0$$

where ${\mathcal{L}}_{\mathbf{v}}$ is the Lie Derivative. Equivalently, if ${\phi }_t:W\to W$ is the flow map generated by $\mathbf{v}$ ($\dot{\phi }=\mathbf{v}\circ \phi$ and ${\phi }_0={\text{id}}_W$) and let ${\psi }_t:={\phi }_t^{-1}$ (inverse flow map), then an advected function satisfies

$$f_t({\phi }_t(X))=f_0(X)\text{ or }{f}_t(x)=f_0({\psi }_t(x))$$

In other words, tracing the flow with ${\psi }_t$ allows us to see where the fluid originally came from. We can convert this to the language of pullbacks.

$$\begin{array}{rl}{f}_t({\phi }_t(X))=f_0(X)& ⟶f_0={\phi }_t^{\ast }{f}_t\\ f_t(x)=f_0({\psi }_t(x))& ⟶f_t={\psi }_t^{\ast }{f}_0\end{array}$$

In fact,

$${\phi }_t^{\ast }=\exp (t{\mathcal{L}}_{\mathbf{v}})\text{ and }{\psi }_t^{\ast }=\exp (-t{\mathcal{L}}_{\mathbf{v}})$$

connecting continuous geometry (the pullback) to discrete numerical methods (the $t$).

We have two numerical methods for advection.

1. Eulerian method: approximate the $\exp (\text{derivative})$ or the pullback operator
2. Lagrangian method: set quantites on particles (or other geometry) and let the particle move by $\mathbf{v}$.

1D Advection of a Function

Suppose we have an advection equation with flow velocity $a=a(x)$.

$$\frac{\partial }{\partial t}u(t,x)+a\frac{\partial }{\partial x}u(t,x)=0\text{ or }\dot{u}+a{u}^{\prime }=0$$

Let $u_i^n$ denote the approximate solution at the $i-$th gridpoint and $n-$th timestep.

Aside: Von Neumann Stability Analysis

We can analyze the stability of a numerical method by looking at how it amplifies or dampens perfect continuous sine waves $r^n\exp (ikx)$. The scaling factor is $r\in \mathbb{C}$. If $|r|>1$, the wave amplifies and the method is unstable. If $|r|<1$, the wave dampens and the method is stable. If $|r|=1$, the wave is preserved and the method is neutrally stable.

Central Difference Method

A straightforward approach is to look at the neighboring cells to find the slope. We then step forward in time by $\Delta t$ using the formula

$$u_i^{n+1}=u_i^n-\frac{a\Delta t}{\Delta x}(u_{i+1}^n-u_{i-1}^n)$$

We can apply Von Neumann Stability Analysis. Let $u_i^n=r^n\exp (\mathbf{i}kx)$. Then we get that

$$r=1-\mathbf{i}\frac{a\Delta t}{\Delta x}\sin (k)$$

Because its magnitude is greater than 1, this method is unconditionally unstable, meaning no matter how small we make $\Delta t$, the method will always explode.

Using ODE Solvers

We can rewrite the advection equation as an ODE by treating $x$ as a parameter.

$$\frac{d}{dt}{u}_i=-\frac{a}{2\Delta x}(u_{i+1}-u_{i-1})$$

This equation can still be stable using ODE solvers like Runge-Kutta Method (RK4) or Forward Euler Method or Backward Euler Method.

But this leaves us with some Taylor expansion errors since we will try to approximate a continuos function on discrete grid points. We must analyze the modified euquation to see what errors we are making. For example, the central difference method generates

$$\dot{u}=-a{u}^{\prime }-\underset{\begin{array}{c}\text{numerical}\\ \text{dispersion}\end{array}}{\underbrace{\frac{\Delta x^2}{12}a{u}^{\prime \prime \prime }}}$$

The numerical dispersion term means waves of different frequencies will travel at different speeds, and cause rippling near sharp corners (like in a square wave).

Lax-Friedrichs Method

One quick fix to the central difference method is to stop using the current cell value, and instead take the average of its left and right neighbors. This gives us the following update rule:

$$u_i^{n+1}=\frac{1}{2}\left(u_{i-1}^n+u_{i+1}^n\right)-\frac{a\Delta t}{2\Delta x}\left(u_{i+1}^n-u_{i-1}^n\right)$$

Using Von Neumann Stability Analysis, we can find that

$$r=\cos (k)-\mathbf{i}\frac{a\Delta t}{\Delta x}\sin (k)$$

It is stable iff

$$\left|\frac{a\Delta t}{\Delta x}\right|\le 1$$

and so it is conditionally stable. We get the modified equation:

$$\dot{u}=-a{u}^{\prime }+\underset{\begin{array}{c}\text{numerical}\\ \text{diffusion}\end{array}}{\underbrace{\frac{\Delta x^2}{2}{u}^{\prime \prime }}}$$

The average injects an artificial diffusion term, which causes the solution to blur over time.

Upwind Method

Recall the $\dot{u}+a{u}^{\prime }=0$ equation. The Upwind Method examines how $a$, the transport coefficient, changes sign.

$$u_i^{n+1}=\{\begin{array}{rl}{u}_i^n-\frac{a\Delta t}{\Delta x}(u_i^n-u_{i-1}^n)& a>0\\ u_i^n-\frac{a\Delta t}{\Delta x}(u_{i+1}^n-u_i^n)& a<0\end{array}$$

This is stable iff

$$\left|\frac{a\Delta t}{\Delta x}\right|\le 1$$

It is conditionally stable. The modified equation is

$$\dot{u}=-a{u}^{\prime }+\underset{\begin{array}{c}\text{numerical}\\ \text{diffusion}\end{array}}{\underbrace{\frac{\Delta x}{2}{u}^{\prime \prime }}}$$

The numerical diffusion is better than Lax-Friedrichs.

Courant-Friedrichs-Lewy (CFL) Condition

The CFL condition ensures Lax-Friedrichs Method and Upwind Method are conditionally stable. They only work if

$$\left|\frac{a\Delta t}{\Delta x}\right|\le 1$$

To understand this physically,

1. $a\Delta t$ is the actual physical distance the fluid travels in one time step.
2. $\Delta x$ is the physical width of one grid cell.

The CFL condition is a strict speed limit: the fluid cannot travel more than one grid cell per time step.

Semi-Lagrangian Scheme

Consider again the advection equation $\dot{u}+a{u}^{\prime }=0$. Sometimes we want a larger time step. We can do this by tracing the flow backwards in time.

$$u_i^{n+1}=u(x_i-a\Delta t)$$

How it works:

1. Stand at the current grid cell center $x_i$.
2. Trace a trajectory backwards in time by $\Delta t$ seconds.
3. Arrive at the physical location ${\psi }_{\Delta t}(x_i)=x_i-a\Delta t$.
4. The starting location will almost certainly not land perfectly on a grid cell center, so interpolate the values of the surrounding grid cells to estimate the fluid data at that exact continuous coordinate.
5. Copy that value into the current grid cell.

This method is unconditionally stable. The modified equation is the same as the Upwind Method.

nD Advection of a Function

The advection equation is now

$$\dot{f}+{\mathcal{L}}_{\mathbf{v}}f=0$$

We want to approximate

$$f^{n+1}={\psi }_{\Delta t}^{\ast }{f}^n$$

The Semi-Lagrangian scheme is somewhat the same. For scalar functions, we use

$$f_i^{n+1}=\text{Interpolate}(f^n,x_i-\Delta t{\mathbf{v}}_i)$$

In 1D, interpolation is simple between two points. In 2D, we need interpolate the surronding 4 grid cells (pixels). In 3D, we need to interpolate the surrounding 8 grid cells (voxels).

Let us call one semi-Lagrangian step as $f^{n+1}-{\mathcal{A}}_{\mathbf{v},\Delta t}^{\text{sL}}(f^n)$. By combining multiple semi-Lagrangian steps, we can bootstrap a higher order method.

Back-and-Forth Error Correction and Compensation (BFECC)

How can we remove the blur caused by the $u^{\prime \prime }$ term if the grid fundamentally requires interpolation? The idea is to trick the grid into calculating its own error and then compensate (subtract) for it.

$${\mathcal{A}}_{\mathbf{v},\Delta t}^{\text{BFECC}}:=\underset{f^{\ast }}{\underbrace{{\mathcal{A}}_{\mathbf{v},\Delta t}^{\text{sL}}}}-\left({\mathcal{A}}_{\mathbf{v},\Delta t}^{\text{sL}}\right)\frac{1}{2}\left(\underset{f^{\ast \ast }}{\underbrace{{\mathcal{A}}_{-\mathbf{v},\Delta t}^{\text{sL}}\circ {\mathcal{A}}_{\mathbf{v},\Delta t}^{\text{sL}}}}-1\right)$$

Here, we are calculating the forward step via $f^{\ast }$, and then finding the backward step by applying the opposite velocity $-\mathbf{v}$ to $f^{\ast }$ giving us $f^{\ast \ast }$. The remaining term is the error, which we take the semi-Lagrangian step of and subtract from $f^{\ast }$. Indeed,

$$f^{\text{final}}=f^{\ast }-{\mathcal{A}}_{\mathbf{v},\Delta t}^{\text{sL}}\underset{\text{error}}{\underbrace{\left(\frac{f^{\ast \ast }-f}{2}\right)}}$$

The tradeoff is that we remove numerical dissipation (blurring) but we get some numerical dispersion (spatial oscillation). We can detect and cutoff overshoots using a limiter.

Limiter

If the value of the field after BFECC is lying outside some convex hull of the neighborhood values of stable semi-Lagrangian, blend with or use the stbale semi-Lagrangian result.

$"\\usepackage{tikz}\n\n\\begin{document}\n\\begin{tikzpicture}[>=stealth, scale=1.1]\n\n$

source code

Footnotes

1. Like in PacMan! ↩

2. This is a term from computer graphics. It refers to the phenomenon where high-frequency signals are misinterpreted as low-frequency signals when sampled at a discrete rate. ↩ ↩² ↩³