電影工業中的流體模擬（六）--- 平滑粒子動力學（1）

編輯：@碩鼠醬，@Raymond

排版：@碩鼠醬

審校：@Raymond @李旻辰

版權宣告：

本文采用CC BY-NC-ND 4.0協議許可，轉載需在顯要位置註明原文出處並附連結，任何組織和個人均不得將本材料用於商業用途。作者保留對侵權行為追究責任的權利。

這一節我們來介紹平滑粒子動力學（Smoothed Particle Hydrodynamics， SPH）的最初的形式。本文的撰寫主要參考［Kelager2006］，@Raymond 的本科論文［Fei2013］以及SPH的原始文獻［Monaghan1977］。

SPH本質上是一種核密度估計（Kernel Density Estimation，KDE）。把空間中的物理量用它周圍一個範圍內的相同物理量透過逼近Delta函式的核函式來進行插值。

聽起來有點繞？我們先來解釋一下核密度估計。

KDE

$\hline$

KDE應用在天文學中計算氣體的物理量是由Lucy等人在1977年的一篇論文［Lucy1977］中提出的。這裡我們需要引入一個特殊的函式，Dirac Delta函式

$\delta(x)$

：

$\delta (x)={\begin{cases}+\infty ,&x=0\\0,&x\neq 0\end{cases}}$

並且有

$\int _{-\infty }^{\infty }\delta (x)\,dx=1.$

（注意，這裡並不是嚴格的定義，詳細的定義請參考［

Dirac delta function

］）

如果我們記空間中處在位置

$\mathbf{r}$

的粒子的物理量為

$A(\mathbf{r})$

，那麼

$A(\mathbf{r})$

在全空間上可以近似表示為

$A_I(\mathbf{r}) = \int_{\Omega} A(\mathbf{r}$

。

這裡我們使用一個近似

$\delta(x)$

的函式

$W(\mathbf{r} - \mathbf{r}$

來代替。其中

是這個函式支撐集（support，也就是函式值不為0的點所組成的集合）的大小。也就是集合

$X = \operatorname {supp} (W)=\{x\in X\,|\,W(x)\neq 0\}$

的大小。

當然，我們需要對

做一些限制：

$\int_{\Omega} W(\mathbf{r} , h) \mathrm{d}\mathbf{r} = 1$

，（注意到這是

delta函式

的性質）

$\lim_{h\to 0^+ } W(\mathbf{r},h) = \delta(\mathbf{r})$

，（注意到

是delta函式的一個近似）

$W(\mathbf{r}, h) \geq 0$

，（這是自然的）

以上過程詳細的推導可以參考［Lucy1977］。

SPH

$\hline$

如果我們把連續的空間劃分成體積元

，每個體積元包含一個粒子，該粒子的物理量為

，由KDE，

$A(\mathbf{r})$

可以近似表示成：

$A_S(\mathbf{r}) = \sum_j A_j V_j W(\mathbf{r} - \mathbf{r}_j, h)$

注意到

$V = \frac{m}{\rho}$

，上式又可以寫成

$A_S(\mathbf{r}) = \sum_j A_j \frac{m_j}{\rho_j}W(\mathbf{r} - \mathbf{r}_j, h)$

這便是SPH的標準形式。有了對量

的表示，我們自然會關心

$\nabla A$

和

$\nabla^2 A$

：

對於

$\nabla A$

，考慮空間中的一個方向

，我們有

$\frac{\partial}{\partial x} \left( A_j \frac{m_j}{\rho_j} W(\mathbf{r} - \mathbf{r}$

這裡使用了求導法則，並且注意到，對於空間中的一個粒子，由於我們是使用拉氏方法的角度看待問題，那麼

與空間位置是無關的。

所以我們有

$\nabla A_s(\mathbf{r}) = \sum_j A_j \frac{m_j}{\rho_j} \nabla W(\mathbf{r} - \mathbf{r}$

同樣地，我們還有

$\nabla^2 A_s(\mathbf{r}) = \sum_j A_j \frac{m_j}{\rho_j} \nabla^2 W(\mathbf{r} - \mathbf{r}$

SPH與

NS方程

$\hline$

下面我們來使用SPH計算NS方程中出現的各種物理量。

回顧前文提到的

動量方程

，把等式左右兩邊同時乘

$\rho$

我們得到（

代表外力，我們把重力算作外力）：

$\rho \frac{D \mathbf{u}}{D t} = - \nabla p + \mu\nabla \cdot \nabla \mathbf{u} + f$

其中

$\mu = \rho \nu$

為

動力粘性係數

（dynamic viscosity）。（感謝評論區

王路飛

指正）

容易看出，我們需要求解的為密度

$\rho$

，

流體壓強

，還有粘滯力以及外力。

對於密度，代入SPH方程我們有

$\rho_i = \sum_j m_j W(\mathbf{r}_i - \mathbf{r}_j, h), \text{with}~i \neq j$

對於壓強，流體的壓強保證流體不會像氣體一樣耗散，需要對流體施加一個

標準密度

$\rho_0$

，在計算流體力學（Computational Fluid Dynamics， CFD）中，我們有

$p = K(\left( \frac{\rho}{\rho_0}\right)^\gamma -1)$

。動機來自下圖［Kelager2006］，為了粒子之間在分離的時候產生吸引力，集中的時候產生互斥力：

不同的

以及

$\gamma$

的選取會帶來不同的視覺效果。這裡我們取

$p = k (\rho - \rho_0)$

。

對於SPH方程，有

$\nabla p_i = -\sum_{j \neq i} p \frac{m_j}{p_j} \nabla W(\mathbf{r}_i - \mathbf{r}_j, h)$

這個方程是不對稱的，違反

牛頓第三定律

。一個對稱形式為

$\nabla p_i = -\rho_i \sum_{j \neq i} \left(\frac{p_i}{\rho_i^2} + \frac{p_j}{\rho_j^2} \right) m _j \nabla W(\mathbf{r}_i - \mathbf{r}_j, h)$

上式的推導來自於，對於任意量

，由求導法則：

$\frac{\nabla A}{\rho} = \nabla \left( \frac{A}{\rho} \right) + \frac{A}{\rho^2} \nabla \rho$

所以

$\nabla A = \rho \left(\nabla \left( \frac{A}{\rho}\right) + \frac{A}{\rho^2} \nabla p\right)$

進而

$\nabla A(\mathbf{r}) = \rho \sum_{j} \left(\frac{A_j}{\rho^2_j} + \frac{A}{\rho^2} \right) m_j \nabla W(\mathbf{r} - \mathbf{r}_j, h)$

對於粘滯力，代入SPH方程一樣可得

$\mu \nabla \cdot \nabla v = \mu \sum_{j \neq i} (\mathbf{u}_j - \mathbf{u}_i) \frac{m_j}{\rho_j} \nabla^2 W(\mathbf{r} - \mathbf{r}_j, h)$

最後放一張總結圖（來自［Kelager2006］）

核函式選取

$\hline$

平滑核的選擇對於模擬的速度、穩定性以及質量影響都非常之大。對於這一點，Müller 等人［Muller2003］提供了一套非常有效的平滑核，分別用於不同的物理量的計算，他們的方法也被目前大部分關於SPH的工作沿用。

其中第一個稱作Poly6的平滑核用於插值密度

$\rho$

（如下圖左）：

$W(\mathbf{r},h)=\begin{cases} \frac{315}{64\pi h^9} (h^2 - r^2)^3 & 0 \leq r \leq h\\ 0 & \text{otherwise} \end{cases}$

其梯度和Laplacian分別為

$\nabla W(\mathbf{r},h)=\begin{cases} -\frac{945}{32 \pi h^9} (h^2 - r^2)^2 \mathbf{r} & 0 \leq r \leq h\\ 0 & \text{otherwise} \end{cases}$

$\nabla^2 W(\mathbf{r},h)=\begin{cases} \frac{945}{8 \pi h^9} (h^2 - r^2)^2 (r^2 - \frac{3}{4}(h^2 - r^2)) & 0 \leq r \leq h\\ 0 & \text{otherwise} \end{cases}$

由於Poly6核的梯度在中心變為0，因此它不適用於插值壓力，為了使粒子

接近時具有較大的壓力，必須使用另一種梯度在0 點有較大取值的平滑核來插

值壓力，稱為Spiky核（如下圖中）：

$W(\mathbf{r},h)=\begin{cases} \frac{15}{ \pi h^6} (h - r)^3 (r^2 - \frac{3}{4}(h^2 - r^2)) & 0 \leq r \leq h\\ 0 & \text{otherwise} \end{cases}$

$\nabla W(\mathbf{r},h)=\begin{cases} -\frac{45}{ \pi h^6 r} (h - r)^2\mathbf{r} & 0 \leq r \leq h\\ 0 & \text{otherwise} \end{cases}$

由於Poly6核的拉氏量在接近0點時很快變為負值，使用它對粘性力插值將會對粒子起到加速的作用，這與粘性力降低粒子速度的特性不符，因此需要使用另一種平滑核插值粘性力，稱為粘度核（如下圖右）：

$W(\mathbf{r},h)=\begin{cases} \frac{15}{2 \pi h^3}(-\frac{r^3}{2h^3} + \frac{r^2}{h^2} + \frac{h}{2r} -1) & 0 \leq r \leq h\\ 0 & \text{otherwise} \end{cases}$

$\nabla^2 W(\mathbf{r},h)=\begin{cases} \frac{45}{\pi h^5}(1-\frac{r}{h}) & 0 \leq r \leq h\\ 0 & \text{otherwise} \end{cases}$

下圖展示上面提到的平滑核和他們的梯度，Laplacian（圖片出處［Fei2013］）

SPH的求解

$\hline$

標準SPH的求解過程可以由下圖概括（［Kelager2006］）

詳細的演算法過程如下：

對於全部的粒子建立空間索引；

對於全部的粒子計算粘性力、合外力，壓強，密度，並根據這些力的作用計算粒子的加速度；

根據粒子的加速度與時間步長得到下一時刻的速度和位置（

時間積分

， Time integration），並處理必要的碰撞。

需要注意的是：

對於步驟1，由於SPH是考察一個粒子被周圍粒子影響效應的疊加，我們需要在整個粒子空間中遍歷多次，時間複雜度為

。加速方法可以採用快速最近鄰查詢演算法（原始文獻參見［Teschner2003］）。［Kelager2006］給這個演算法做了一個簡單的概括；

對於步驟2，對於某些確定的物理量，我們需要決定哪些粒子對我們考察的粒子產生影響。這些粒子的範圍是由核函式的支撐集長度決定的；流體的不可壓縮性是由我們人為引入的壓強所保證的；

對於步驟3，我們可以使用很多時間積分方法，比如前向尤拉（Forward Euler），蛙跳（Leap Frog）等。蛙跳法的具體形式如下

$\begin{eqnarray} \mathbf{u}_{t+\frac{1}{2}\Delta t} &=& \mathbf{u}_{t- \frac{1}{2}\Delta t} + \Delta t \mathbf{a}_t, \\ \mathbf{r}_{t+ \Delta t} &=& \mathbf{r}_t + \Delta t \mathbf{u}_{t + \frac{1}{2}\Delta t} \end{eqnarray}$

其中

$\begin{eqnarray} \mathbf{u}_{-\frac{1}{2}\Delta t} &=& \mathbf{u}_0 - \frac{1}{2} \Delta t \mathbf{a}_0, \\ \mathbf{u}_t &\approx &\frac{\mathbf{u}_{t+\frac{1}{2}\Delta t} + \mathbf{u}_{t- \frac{1}{2}\Delta t} }{2} \end{eqnarray}$

實現例子

$\hline$

SPH的實現例子（此處是作者的一個小失誤，這是SPH的另一個變種，而不是我們以後要介紹的PCISPH，程式碼的實現實際上是參考了［Clavet2005］）：

我們來結合這份程式碼看一下具體的實現過程（PS，執行這個程式碼需要安裝OpenGL/GLUT以及Eigen3。OpenGL一個比較不錯的參考教程在這裡。

程式顯示的原理是GLUT不斷重新整理當前的顯示視窗（glutDisplayFunc（Render）），並且在後臺處理下一幀的顯示（glutIdleFunc（Update）），兩個佈景（前景/後景）交替顯示達到動畫的效果（glutSwapBuffers（））。程式使用Render 繪製效果圖，並使用Update 更新各個粒子的狀態以供下一步顯示；

void InitGL（void）和 void InitSPH（void），設定外觀引數；

void InitSPH（void）中的 void GridInsert（void）對空間建立粒子索引，給每個粒子編號，並用連結串列的形式儲存方便遍歷粒子列表；

來看Update 。（不同於標準的SPH，這部分的過程對應於［Clavet2005］的Algorithm 1）：

for（int i = 0； i < SOLVER_STEPS； i++）{

ApplyExternalForces（）；

// 施加外力，可以看到這裡僅考慮了重力；

Integrate（）； //時間積分，這裡使用了前向尤拉方法；

GridInsert（）； //更新時間積分以後的粒子狀態；

PressureStep（）； //計算壓強和密度；

Project（）； //施加壓強，表面張力，粘性力，並計算下一步粒子的位置；

Correct（）；//計算速度；

GridInsert（）；//回寫粒子資訊，準備下一輪使用；

EnforceBoundary（）； //施加邊界條件，防止粒子溢位容器。

}

// 詳細的演算法解釋請參考［Clavet2005］。

標準形式的SPH就介紹完了。下一章我們將講述預校不可壓SPH（PCISPH），敬請期待！

P。S。來自評論區網友@葉言由

現在暫時沒有成熟的商業軟體。國內一般翻譯成光滑粒子動力學比較多。北京大學

劉謀斌

老師（曾任中科院力學所研究員）曾經出過第一本SPH的專著，叫“一種無網格化方法-光滑粒子動力學”。初學的可以用這本教材。該方法在爆炸，多相流等方面均有應用。仍處在發展階段。計算量是一個缺點，但是在模擬大密度變化，自由液麵，流固耦合，多相流有天然優勢。克服缺點的辦法就是大規模並行或者GPU運算。有一個開原始碼叫dualSPHysics，大家想學習可以去下載，歐洲幾個大學的學者共同開發的，有比較詳細的技術手冊供大家學習。

REFERENCES

［Clavet2005］Clavet， Simon， Philippe Beaudoin， and Pierre Poulin。 “Particle-based viscoelastic fluid simulation。” Proceedings of the 2005 ACM SIGGRAPH/Eurographics symposium on Computer animation。 ACM， 2005。

［Dirac delta function］ Dirac

delta function

， wikipedia

［Interpolation］ Interpolation， wikipedia

［Jin2005］Hongbin， Jin， and Ding Xin。 “On criterions for smoothed particle hydrodynamics kernels in stable field。” Journal of Computational Physics 202。2 （2005）： 699-709。

［Kelager2006］ Kelager， Micky。 “Lagrangian fluid dynamics using smoothed particle hydrodynamics。” University of Copenhagen： Department of Computer Science （2006）。

［Lucy1977］ Lucy， Leon B。 “A

numerical approach

to the testing of the fission hypothesis。” The

astronomical journal

82 （1977）： 1013-1024。

［Monaghan1977］ Monaghan， Joe J。 “Smoothed particle hydrodynamics。” Annual review of astronomy and astrophysics 30。1 （1992）： 543-574。

［Muller2003］ Müller， Matthias， David Charypar， and Markus Gross。 “Particle-based fluid simulation for interactive applications。” Proceedings of the 2003 ACM SIGGRAPH/Eurographics symposium on Computer animation。 Eurographics Association， 2003。

［Fei2013］

費昀

。不可壓平滑粒子流體動力學演算法 GPU 並行加速及其應用研究。 Diss。 2013。

［Teschner2003］ Teschner， Matthias， et al。 “Optimized Spatial Hashing for Collision Detection of Deformable Objects。” Vmv。 Vol。 3。 2003。

\_（：3」∠）\_ \_（・ω・”∠）\_ \_（：з）∠）\_ ∠（ ᐛ 」∠）＿ \_（：зゝ∠）\_

請毫不猶豫地關注我們：

我們的網站：GraphiCon

知乎專欄：GraphiCon圖形控 - 知乎專欄

公眾號：GraphiCon

如果你有什麼想法，建議，或者想加入我們，你可以：

給我們發郵件：hi@graphicon。io

加入我們的QQ群：SIQGRAPH（342086343）

加入我們的slack群 Join us on Slack！

__GraphiCon長期接受投稿，如果你想投稿給我們可以透過上面的方式聯絡我們！__

小蜜蜂問答

小蜜蜂問答

電影工業中的流體模擬（六）--- 平滑粒子動力學（1）

推薦文章

小蜜蜂問答

小蜜蜂問答

電影工業中的流體模擬（六）--- 平滑粒子動力學（1）

相關文章

在職場中年危機不可避免，你知道如何應對嗎？

僅次於光速的速度？

使用者自定義函式包含的兩個部分是什麼？

Java高階特性——註解，這也許是最簡單易懂的文章了

推薦文章