[Note] Position Based Fluids

原論文：M.Macklin and M. Müller. Position Based Fluids. ACM Transactions on Graphics, 2013.

1. Introduction

對於互動的環境來說，流體模擬的穩定性非常重要。Smoothed Particle Hydrodynamics (SPH) 是一個在互動環境中被廣泛使用的 Particle-Based 方法，然而他卻有個致命的缺點：當粒子周圍的相鄰粒子過少時會很難維持住這個算法的穩定性，尤其以在流體表面、或是在邊緣的流體粒子更常發生。將每個 Time step 調到足夠小，或是增加足夠的粒子數量即使可避免掉這項問題，卻會大幅度增加計算成本。

Position Based Dynamics (PBD) 在遊戲開發或是影片製作中都是很受歡迎的一套物理模擬方式，作者選擇使用 PBD 正是因為其具有 Unconditionally Stable 以及穩定性等性質。這篇 Paper 將介紹如何使用 PBD Framework 來模擬 Incompressible flow，並克服上述在 Free Surfaces 發生 Particle Deficiency 的問題。

Muller [2003] 等人在 Particle-Based Fluid Simulation for Intereactive Applicatoins 中提出使用 Smoothed Particle Hydrodynamics (SPH) 模擬具有 Viscosity 跟 Surface Tension 的流體
為了維持住 Incompressibility，Weakly Compressible SPH (WCSPH) [2007] 與標準 SPH 所使用的 Stiff Equation 會大大限制住 Time step 的長度。
Predictive-Corrective Incompressible SPH [2009] 利用 Iterative Jacobi-style 方法，藉由不斷迭代、累積流體壓力並逐步施力的方式，來確保流體能夠在較長的Time step 上也能夠穩定存在，而不需要額外設置 Stiffness value 且可以分散掉不斷矯正相鄰粒子密度的計算成本。
Bodin [2012] 將 Incompressibility 組成一個 Velocity Constraints 的線性方程，並使用 Gauss-Seidel Iteration 解出線性方程來確保流體密度的一致。相反的，Position-Based 跟 PCISPH 則是使用類似 Jacobi Iteraion 的方法，解出非線性方程，並不斷重新估計誤差與梯度。
Hybrid Method，像是 Fluid Implicit-Particle (FLIP) 則是使用 Grid 來解出流體壓力並將流體速度擴散到粒子上模擬。Zhu 與 Bridson [2005] 結合 PIC/FLIP 方法。Raveendran [2011] 使用 Coarse grid 方法混合 SPH 趨近 divergence free velocity field。
Clavet [2005] 使用 Position Based 方法模擬 Viscoelastic Fluids。但是他們的方法只是 Conditionally Stable。
Position Based Dynamics [2007] 基於 Verlet integration 提供一個在遊戲中模擬物理的方法。

3. Enforcing Incompressibility

對每一個粒子使用 Density Constraint 來確保流體密度不變。Constraint 是每個粒子位置與鄰近粒子位置的函數，位置寫作 $\mathbf{p}_1,\cdots,\mathbf{p}_n$，以下是對第$i$個粒子的 Constraint Function：

$C_i\big(\mathbf{p}_1,\cdots,\mathbf{p}_n\big)=\frac{\rho_i}{\rho_0}-1\tag{1}$

其中 $\rho_0$是靜止密度 (Rest Density)，$\rho_i$ 是由 SPH density estimator 計算出來的粒子密度:

$\rho_i=\sum_jm_jW\big(\mathbf{p}_i-\mathbf{p}_j, h\big)\tag{2}$

使用 Poly6 作為 Density Estimator 的 Kernel。並使用 Spiky Kernel 計算 gradient

註: 由於原論文中考慮粒子質量 $m_j$ 為定值且相同，因此在後續推倒公式中皆省略不寫，但是在本篇 Note 中會將他歸至原位。

註: Poly 6 為 $6^{th}$ degree polynomial kernel 的簡稱，以下為其函數:

$W_{\text{poly6}}\big(\mathbf{r},h\big)=\frac{315}{64\pi h^9}\begin{cases} \big(h^2-\lVert\mathbf{r}\rVert^2\big)^3, &0\le\lVert\mathbf{r}\rVert\le h\\ 0, &\lVert\mathbf{r}\rVert\gt h\end{cases}$

Poly 6 的 Gradient:

$\nabla W_{\text{poly6}}\big(\mathbf{r},h\big)=-\frac{945}{32\pi h^9}\mathbf{r} \big(h^2-\lVert\mathbf{r}\rVert^2\big)^2$

Poly 6 的 Laplacian:

$\nabla^2 W_{\text{poly6}}\big(\mathbf{r},h\big)=-\frac{945}{32\pi h^9}\big(h^2-\lVert\mathbf{r}\rVert^2\big)\big(3h^2-7\lVert\mathbf{r}\rVert^2\big)$

註: Spiky Kernel 函數:

$W_{\text{spiky}}\big(\mathbf{r},h\big)=\frac{15}{\pi h^6}\begin{cases} \big(h-\lVert\mathbf{r}\rVert\big)^3, &0\le\lVert\mathbf{r}\rVert\le h\\ 0, &\lVert\mathbf{r}\rVert\gt h\end{cases}$

Spiky Kernel 的 Gradient:

$\begin{gathered}\nabla W_{\text{spiky}}\big(\mathbf{r},h\big)=-\frac{45}{\pi h^6}\frac{\mathbf{r}}{\lVert\mathbf{r}\rVert}\big(h-\lVert\mathbf{r}\rVert\big)^2,\\\lim_{r\rightarrow0^-}\nabla W_{\text{spiky}}\big(\mathbf{r},h\big)=\frac{45}{\pi h^6},\\\lim_{r\rightarrow0^+}\nabla W_{\text{spiky}}\big(\mathbf{r},h\big)=-\frac{45}{\pi h^6},\\\end{gathered}$

Spiky Kernel 的 Laplacian:

$\begin{gathered}\nabla^2 W_{\text{spiky}}\big(\mathbf{r},h\big)=-\frac{90}{\pi h^6}\frac{1}{\lVert\mathbf{r}\rVert}\big(h-\lVert\mathbf{r}\rVert\big)\big(h-2\lVert\mathbf{r}\rVert\big),\\\lim_{r\rightarrow0}\nabla^2 W_{\text{spiky}}\big(\mathbf{r},h\big)=-\infty\end{gathered}$

PBD 的目標在於找出粒子位置的修正項 $\Delta\mathbf{p}$來滿足Constraint:

$C\big(\mathbf{p}+\Delta\mathbf{p}\big)=0\tag{3}$

使用牛頓法來求解:

$\Delta\mathbf{p}\approx\nabla C(\mathbf{p})\lambda \tag{4}$ $\begin{align} C\big(\mathbf{p}+\Delta\mathbf{p}\big) &\approx C\big(\mathbf{p}\big)+\nabla C^T\Delta\mathbf{p}=0\tag{5}\\ &\approx C\big(\mathbf{p}\big)+\nabla C^T\nabla C\lambda=0\tag{6} \end{align}$

註: $\lambda$ 是一個變量。這個思路其實很簡單，我們要追求每一個 Time step 都要滿足約束函數 $C$，但是粒子的位置 $\mathbf{p}$ 不一定會滿足，因此我們必須求出修正項 $\Delta \mathbf{p}$ 來修正粒子位置直到滿足約束條件。而修正項可以定義為往 $\nabla C(\mathbf{p})$ 方向乘上一個微小變量 $\lambda$。根據 Taylor 一階展開式，$C\big(\mathbf{p}+\Delta\mathbf{p}\big)\approx C\big(\mathbf{p}\big)+\nabla C^T\Delta\mathbf{p}$ 最後再將 $\Delta \mathbf{p}$ 帶入，得到式(6)。

註: $\nabla C^T$ 中的 $T$ 是 Transpose，由於 $\nabla C$ 是向量，$\nabla C^T\Delta\mathbf{p}$ 表示的是向量$\nabla C$與向量$\Delta\mathbf{p}$ 的內積，因此 $\nabla C$ 會需要 Transpose。

根據 SPH 方法，我們可以將粒子 $i$ 對粒子 $k$ 的約束函數 $C$ 梯度定義為：

$\nabla_{\mathbf{p}_k}C_i=\frac{1}{\rho_0}\sum_j m_j\nabla_{\mathbf{p}_k}W\big(\mathbf{p}_i-\mathbf{p}_j, h\big)\tag{7}$

根據 $i$ 與 $k$ 的關係分為:

$\nabla_{\mathbf{p}_k}C_i=\frac{1}{\rho_0} \begin{cases} \displaystyle \sum_j m_j\nabla_{\mathbf{p}_k}W\big(\mathbf{p}_i-\mathbf{p}_j, h\big) &\text{if $k=i$}\\ \displaystyle -m_j\nabla_{\mathbf{p}_k}W\big(\mathbf{p}_i-\mathbf{p}_j, h\big) &\text{if $k=j$} \end{cases} \tag{8}$

註: $i$ 是自身粒子，$j$ 是鄰近粒子。由於 Kernel $W$ 是 $\mathbf{p}_i-\mathbf{p}_j$ 的函數，因此對於其他 $k\ne i$ 與 $k\ne j$ 的粒子來說其梯度皆為 $0$，只有 $k=i$ 與 $k=j$ 會滿足條件。至於對 $\mathbf{p}_i$ 或 $\mathbf{p}_j$ 的梯度差別只在於方向: 由於對象都是座標 $\mathbf{p}$，因此梯度皆為針對座標每個維度作微分，只是由於 $W$ 是 $\mathbf{p}_i-\mathbf{p}_j$ 的函數，如果是對 $\mathbf{p}_i$ 取 $W$ 梯度則根據 Chain Rule 得到 $\nabla_{\mathbf{p}_i}W=W’\nabla_{\mathbf{p}_i}(\mathbf{p}_i-\mathbf{p}_j)=W’$，而對 $\mathbf{p}_j$ 取 $W$ 梯度則得到 $\nabla_{\mathbf{p}_j}W=W’\nabla_{\mathbf{p}_j}(\mathbf{p}_i-\mathbf{p}_j)=-W’$，因此式(8)在 $k=j$的情況下加了一個負號，將方向導正。

將式(8)代入式(6)求得變量 $\lambda$:

$\lambda_i=-\frac{C_i\big(\mathbf{p}_i,\cdots,\mathbf{p}_j\big)}{\sum_k\big|\nabla_{\mathbf{p}_k}C_i\big|^2}\tag{9}$

註: $|\nabla_{\mathbf{p}_k}C_i|^2$ 是因為式(6)中的 $\nabla C^T\nabla C$

由於約束函數(1)是非線性，當粒子逐漸分離時，Kernel的邊緣值會趨近0，導致式(9)的梯度分母逐漸趨近0，進而造成整體模擬的不穩定性。使用 Constraint Force Mixing (CFM) 則可以避免這個情況，將式(6)改寫為:

$C\big(\mathbf{p}+\Delta\mathbf{p}\big) \approx C\big(\mathbf{p}\big)+\nabla C^T\nabla C\lambda+\varepsilon\lambda=0\tag{10}$

其中 $\varepsilon$ 是使用者自訂的 Relaxation Parameter，$\lambda_i$變成:

$\lambda_i=-\frac{C_i\big(\mathbf{p}_i,\cdots\mathbf{p}_j\big)}{\sum_k\big|\nabla_{\mathbf{p}_k}C_i\big|^2+\varepsilon}\tag{11}$

註: 簡而言之就是在分母的部份加上一個不為 $0$ 的微小常數來確保分母不會為 $0$

最後修正項 $\Delta \mathbf{p}_i$ 的函數變成:

$\Delta \mathbf{p}_i=\frac{1}{\rho_0}\sum_j\big(\lambda_i+\lambda_j\big)m_j\nabla W\big(\mathbf{p}_i-\mathbf{p}_j, h\big)\tag{12}$

註:

$\begin{align}\Delta \mathbf{p}_i &= \lambda_i\nabla_{\mathbf{p}_i}C_i+\sum_j\lambda_j\nabla_{\mathbf{p}_j}C_i\\ &= \frac{1}{\rho_0}\sum_j \lambda_i m_j\nabla_{\mathbf{p}_i}W\big(\mathbf{p}_i-\mathbf{p}_j, h\big) + \Big(-\frac{1}{\rho_0}\sum_j\lambda_j m_j\nabla_{\mathbf{p}_j}W\big(\mathbf{p}_i-\mathbf{p}_j, h\big)\Big)\\ &= \frac{1}{\rho_0}\sum_j \lambda_i m_j\nabla_{\mathbf{p}_i}W\big(\mathbf{p}_i-\mathbf{p}_j, h\big) + \frac{1}{\rho_0}\sum_j\lambda_j m_j\nabla_{\mathbf{p}_i}W\big(\mathbf{p}_i-\mathbf{p}_j, h\big) \\ &= \frac{1}{\rho_0}\sum_j\big(\lambda_i+\lambda_j\big)m_j\nabla_{\mathbf{p}_i}W\big(\mathbf{p}_i-\mathbf{p}_j, h\big)\end{align}$

其中第2行到第3行是因為:

$\because\nabla_{\mathbf{p}_j}W\big(\mathbf{p}_i-\mathbf{p}_j, h\big)=-\nabla_{\mathbf{p}_i}W\big(\mathbf{p}_i-\mathbf{p}_j, h\big)$

4. Tensile Instability

SPH 中常見的問題就是因鄰近粒子數量不足無法達到 Rest Density，而產生負壓力所導致的粒子的聚集效應 (粒子間的排斥力變成吸力)。

上圖為聚集效應產生的不自然現象 / 下圖為施加人工壓力項的結果

一種解決方案是對壓力做 Clamping 讓壓力不為負值，但是這會讓粒子的內聚力衰弱。

註:

1
2
3

clamp(x) {
    return x>0 ? x:0;
}

其他解決方案:

Clavet [2005] 使用 second near pressure term
Alduan and Otaduy [2011] 使用 discrete element (DEM) forces
Schechter and Bridson [2012] 使用 ghost particles

此篇論文使用 Monaghan [2000] 的方法，加上一個人工壓力項:

$s_{corr}=-k\bigg(\frac{W(\mathbf{p}_i-\mathbf{p}_j, h)}{W(\Delta \mathbf{q}, h)}\bigg)^n\tag{13}$

其中 $\Delta \mathbf{q}$ 是固定的距離，$k$ 是微小的正常數。通常取 $|\Delta \mathbf{q}|=0.1h, \cdots,0.3h$，$k=0.1$，$n=4$。將此項納入修正項:

$\Delta \mathbf{p}_i=\frac{1}{\rho_0}\sum_j\big(\lambda_i+\lambda_j+s_{corr}\big)m_j\nabla W\big(\mathbf{p}_i-\mathbf{p}_j, h\big)\tag{14}$

這個醜陋的項會保持粒子略低於 Rest Density，最終產生出類似 Surface Tension 的效果。

5. Vorticity Confinement and Viscosity

PBD 方法通常會產生額外的 Damping 造成渦流快速消散。Fedkiw [2001] 引進 Vorticity Confinement 方法來克服模擬煙霧時的消散問題(Numerical Dissipation)，後來由 Lentine [2011] 當作 Energy Conserving 引入流體模擬。Hong [2008] 展示如何將 Vorticity Confinement 導入 Hybrid 的模擬方法。

左邊沒有加 Vorticity Confinement / 右邊有加 Vorticity Confinement

此篇文章使用 Vorticity Confinement 來補充流失的能量。此方法需要先計算粒子位置的 Vorticity，Monaghan [1992]:

$\omega_i=\nabla\times\mathbf{v}=\sum_j m_j\mathbf{v}_{ij}\times\nabla\mathbf{p}_j W\big(\mathbf{p}_i-\mathbf{p}_j, h\big)\tag{15}$ $\mathbf{f}^{\text{vorticity}}_i=\varepsilon\big(\mathbf{N}\times\omega_i\big)\tag{16}$

其中 $\mathbf{N}=\frac{\eta}{|\eta|}$，$\eta=\nabla|\omega|_i$， $\mathbf{v}_{ij}=\mathbf{v}_{j}-\mathbf{v}_{i}$ 。

註: 這邊有點小 Tricky，原論文裏面也沒有寫清楚，首先這邊目標是加速粒子的渦流速度，因此他先算出旋度向量 $\omega_i$，也就是遵守右手定則的旋轉軸心方向 (大姆指指的方向)。接著用 $\omega_i$ 的 Gradient L2-norm 算出往旋轉軸心方向的向量 $\eta$，算出單位向量 $\mathbf{N}$，最後再用外積算出旋轉方向並乘上 $\varepsilon$ 算出 Vorticity forces

註:

$\begin{align}\eta &= \nabla|\omega|=\nabla\bigg(\sum_{k=1}^n\omega_k^2\bigg)^{\frac{1}{2}}\\ &= \sum_{j=1}^n\frac{\partial}{\partial \omega_j}\bigg(\sum_{k=1}^n\omega_k^2\bigg)^{\frac{1}{2}}\hat{\omega_j}\\ &= \frac{1}{2}\sum_{j=1}^n\frac{2\omega_j}{\Big(\displaystyle\sum_{k=1}^n\omega_k^2\Big)^{\frac{1}{2}}}\hat{\omega_j}\\ &=\sum_{j=1}^n\frac{\omega_j}{|\omega|}\hat{\omega_j}\end{align}$

此外，此篇論文也使用 XSPH viscosity 來模擬流體的黏滯力:

$\mathbf{v}_i^{new}=\mathbf{v}_i+c\sum_jm_j\mathbf{v}_{ij}\cdot W\big(\mathbf{p}_i-\mathbf{p}_j, h\big)\tag{17}$

其中，參數 $c$ 通常設為 $0.1$。

6. Algorithm

7. Rendering

使用 GPU based ellipsoid splatting technique 詳情見液体渲染：一种屏幕空间方法