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

01-26

編輯：@碩鼠醬，@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)={egin{cases}+infty ,&x=0\0,&x eq 0end{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}$ 來代替。其中 $h$ 是這個函數支撐集（support，也就是函數值不為0的點所組成的集合）的大小。也就是集合

$X = operatorname {supp} (W)={xin X,|,W(x) eq 0}$ 的大小。

當然，我們需要對 $W$ 做一些限制：

$int_{Omega} W(mathbf{r} , h) mathrm{d}mathbf{r} = 1$ ，（注意到這是delta函數的性質）
$lim_{h o 0^+ } W(mathbf{r},h) = delta(mathbf{r})$ ，（注意到 $W$ 是delta函數的一個近似）
$W(mathbf{r}, h) geq 0$ ，（這是自然的）

以上過程詳細的推導可以參考[Lucy1977]。

SPH

$hline$

如果我們把連續的空間劃分成體積元 $V_j$ ，每個體積元包含一個粒子，該粒子的物理量為 $A_j$ ，由KDE， $A(mathbf{r})$ 可以近似表示成：

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

注意到 $V = frac{m}{ ho}$ ，上式又可以寫成

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

這便是SPH的標準形式。有了對量 $A$ 的表示，我們自然會關心 $abla A$ 和 $abla^2 A$ ：

對於 $abla A$ ，考慮空間中的一個方向 $x$ ，我們有 $frac{partial}{partial x} left( A_j frac{m_j}{ ho_j} W(mathbf{r} - mathbf{r}$

這裡使用了求導法則，並且注意到，對於空間中的一個粒子，由於我們是使用拉氏方法的角度看待問題，那麼 $A_j, m_j, V_j$ 與空間位置是無關的。

所以我們有

$abla A_s(mathbf{r}) = sum_j A_j frac{m_j}{ ho_j} abla W(mathbf{r} - mathbf{r}$

同樣地，我們還有

$abla^2 A_s(mathbf{r}) = sum_j A_j frac{m_j}{ ho_j} abla^2 W(mathbf{r} - mathbf{r}$

SPH與NS方程

$hline$

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

回顧前文提到的動量方程，把等式左右兩邊同時乘 $ho$ 我們得到（ $f$ 代表外力，我們把重力算作外力）： $ho frac{D mathbf{u}}{D t} = - abla p + mu abla cdot abla mathbf{u} + f$

其中 $mu = ho u$ 為動力粘性係數（dynamic viscosity）。（感謝評論區王路飛指正）

容易看出，我們需要求解的為密度 $ho$ ，流體壓強 $p$ ，還有粘滯力以及外力。

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

$ho_i = sum_j m_j W(mathbf{r}_i - mathbf{r}_j, h), ext{with}~i eq j$

對於壓強，流體的壓強保證流體不會像氣體一樣耗散，需要對流體施加一個標準密度 $ho_0$ ，在計算流體力學（Computational Fluid Dynamics, CFD）中，我們有 $p = K(left( frac{ ho}{ ho_0} ight)^gamma -1)$ 。動機來自下圖[Kelager2006]，為了粒子之間在分離的時候產生吸引力，集中的時候產生互斥力：

不同的 $k$ 以及 $gamma$ 的選取會帶來不同的視覺效果。這裡我們取 $p = k ( ho - ho_0)$ 。

對於SPH方程，有

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

這個方程是不對稱的，違反牛頓第三定律。一個對稱形式為

$abla p_i = - ho_i sum_{j eq i} left(frac{p_i}{ ho_i^2} + frac{p_j}{ ho_j^2} ight) m _j abla W(mathbf{r}_i - mathbf{r}_j, h)$

上式的推導來自於，對於任意量 $A$ ，由求導法則:

$frac{ abla A}{ ho} = abla left( frac{A}{ ho} ight) + frac{A}{ ho^2} abla ho$

所以

$abla A = ho left( abla left( frac{A}{ ho} ight) + frac{A}{ ho^2} abla p ight)$

進而

$abla A(mathbf{r}) = ho sum_{j} left(frac{A_j}{ ho^2_j} + frac{A}{ ho^2} ight) m_j abla W(mathbf{r} - mathbf{r}_j, h)$

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

$mu abla cdot abla v = mu sum_{j eq i} (mathbf{u}_j - mathbf{u}_i) frac{m_j}{ ho_j} abla^2 W(mathbf{r} - mathbf{r}_j, h)$

最後放一張總結圖（來自[Kelager2006]）

核函數選取

$hline$

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

其中第一個稱作Poly6的平滑核用於插值密度 $ho$ （如下圖左）：

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

其梯度和Laplacian分別為

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

$abla^2 W(mathbf{r},h)=egin{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 & ext{otherwise}end{cases}$

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

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

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

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

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

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

$W(mathbf{r},h)=egin{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 & ext{otherwise}end{cases}$

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

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

SPH的求解

$hline$

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

詳細的演算法過程如下：

對於全部的粒子建立空間索引；
對於全部的粒子計算粘性力、合外力，壓強，密度，並根據這些力的作用計算粒子的加速度；
根據粒子的加速度與時間步長得到下一時刻的速度和位置（時間積分, Time integration），並處理必要的碰撞。

需要注意的是：

對於步驟1，由於SPH是考察一個粒子被周圍粒子影響效應的疊加，我們需要在整個粒子空間中遍歷多次，時間複雜度為 $O(n)$ 。加速方法可以採用快速最近鄰查找演算法（原始文獻參見[Teschner2003]）。[Kelager2006] 給這個演算法做了一個簡單的概括；
對於步驟2，對於某些確定的物理量，我們需要決定哪些粒子對我們考察的粒子產生影響。這些粒子的範圍是由核函數的支撐集長度決定的；流體的不可壓縮性是由我們人為引入的壓強所保證的；
對於步驟3，我們可以使用很多時間積分方法，比如前向歐拉（Forward Euler），蛙跳（Leap Frog）等。蛙跳法的具體形式如下

$egin{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}$

其中

$egin{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長期接受投稿，如果你想投稿給我們可以通過上面的方式聯繫我們！__
推薦閱讀：

※用Python做數值計算-給定任意公式生成波數法公式的程序
※很想掌握如何編程求解一維水動力模型，具備基本理論基礎，請問如何一步步實現這一目標？
※數值實驗怎麼進行參數率定？
※單晶高溫合金行業現狀及其數值模擬有什麼應用？

TAG:计算流体力学CFD | 计算机图形学 | 数值模拟 |