物理模擬——PBF流體模擬（1）

本文是有關我大學畢設的一個總結，畢設題目為：基於粒子法流體動力學的物理仿真引擎開發，實際工作為基於圖形API，搭建起一套簡單的渲染框架，並且基於此框架實現CPU端流體模擬算法Position Based Fluid，倉庫地址為：https://gitee.com/FlyingZiming/fluid-simulation-engine

目錄

• 介紹
• 基於物理的動畫
  • 基礎概念
  • 流體動畫
• 流體模擬基礎
  • 梯度
  • 散度
  • 旋度
  • 拉普拉斯算子
  • 粘性
  • 納維-斯托克斯偏微分方程組
  • 拉格朗日視角和歐拉視角
• SPH
  • 原理及公式
  • 過程
  • 加速
• 引用

介紹

PBF全稱Position Based Fluid，意為基於位置的流體，是基於Position Based Dynamics實現的流體模擬，本文會講解物理動畫概念、流體模擬基礎、SPH、PBD、PBF以及流體渲染。

筆者水平有限，有說的不對的地方歡迎指出。

基於物理的動畫

基礎概念

基於物理的動畫意味着利用物理原理，通常是經典力學，去為復雜的現象建模，模擬時間推進，實現成模擬算法。

他們的目標是使用一系列的建模方法，隨着時間推演，生成一系列的系統狀態，再把這些狀態渲染為圖像。這種動畫方法通常是作用於復雜的動畫，其表現通常是不被預設好的，根據環境變化有所變化的，因此無法按照傳統的制作關鍵幀方式來進行動畫制作，比如一個彈性球落地彈起、布料隨風飄動、雨水落在地面。

創建物理模擬動畫有兩個元素：模型和模擬。模型是一系列定律或規則，通常是物理定律。模擬是將模型封裝進一個框架之中從而預測該模型在時間推演下的演化。基於物理的動畫的模型就是物理的，是模擬現實世界的行為，但這種模型其實也是被簡化過的，影響較小的部分通常會被忽略，模型注重宏觀行為，在具體細節上不是特別在意。

具體的基於物理的動畫在使用中通常要考慮很多要素，動畫希望動畫具有實時性好、相對穩定的以及可以方便控制等等的性質。為了達到這些效果，動畫師有機會重新發明物理定律，以不同的方式描述我們的世界，只要他們模擬的時候看起來像。

流體動畫

3d游戲中的流體動畫一般的實現比較復雜，不同規模的流體有不同的做法，甚至不同的材料都會使用不同的模擬方法。

小規模的流體，比如一杯水，可以使用SPH等粒子法去模擬，游戲中典型的實時粒子法就是Position Based Fluid，NVIDIA 的Flex便是使用此算法模擬的流體，跑在GPU上的；

中等規模的流體，通常使用網格法，歐拉視角求解NS方程，網格法適合處理時常充滿空間的流體，比如煙霧等；

大規模的流體，一般就不能使用物理的方法了，而是使用波函數疊加的方式，去模擬一個場，用在海平面的高度場模擬，具體如Simulating Ocean Water。

流體模擬基礎

接下來簡單講解一下針對模擬領域需要使用的基礎知識、概念等，不需要完全理解。

梯度

梯度是一個矢量，他的每一個分量是一個關於空間一個軸的偏導數，2維的梯度見式：

\[\nabla f(x,y) =(\frac {\partial f}{\partial x},\frac {\partial f}{\partial y}) \]

3維的梯度見式：

\[\nabla f(x,y,z) =(\frac {\partial f}{\partial x},\frac {\partial f}{\partial y},\frac {\partial f}{\partial z}) \]

形象的理解，梯度就是函數在空間變化最快的方向。梯度可以用來估算函數的值，具體表現為同一階泰勒泰勒展開一樣的表現形式。

散度

散度面向的對象是一個向量場，用於衡量某一點的矢量的聚集或者發散程度，結果是一個標量可以將其理解為一個矢量和梯度算子的點乘，值為該矢量的三個空間偏導數的和，如果矢量場的散度為0，稱該矢量場無散度。

二階的散度表示見式：

\[\nabla \cdot \vec{u} = \nabla \cdot (u,v) = \frac{\partial u}{\partial x}+ \frac {\partial v}{\partial y} \]

旋度

旋度表示圍繞某一點的旋轉速度，三維形式旋度是一個向量，二維形式的旋度是一個標量，旋度不能簡單的進行維度推廣，旋度符號可以理解為是梯度算子和一個矢量的叉乘。二維的旋度表示見式：

\[\nabla \times \vec{u} = \nabla \times (u,v) = \frac{\partial v}{\partial x}- \frac {\partial u}{\partial y} \]

拉普拉斯算子

拉普拉斯算子定義為梯度的散度，符號表示為 ，拉普拉斯算子是一個二階的微分算子，偏微分方程 被稱為拉普拉斯方程，右邊如果是非零項則被稱為泊松方程。拉普拉算子2維的表示見式：

\[\nabla \cdot \nabla f =\frac{\partial^2 f}{\partial x^2}+ \frac {\partial^2 f}{\partial y^2} \]

粘性

粘性是流體持續剪切變形時內部產生剪切力的性質，他宏觀上表現為流體之間的動摩擦力，本質上是分子間的作用力，不同速度分子的化學鍵的斷裂和重新生成，兩層液體分子之間的牽扯和擠壓。

描述粘性使用兩種粘性系數，動力粘性系數和運動粘性系數。

$\mu $ 動力粘性系數表示相同變形率是的粘性大小， $ \upsilon $ 運動粘性系數表示的是相同加速度時的粘性大小，在Navier-Stokes方程中出現的即是運動粘性系數。

納維-斯托克斯偏微分方程組

流體模擬面臨的最重要的方程就是Navier-Stokes方程，不論以何種模擬方法，都是求解這個流體運動方程，見式：

\[\frac{\partial \vec{u}}{\partial {t}} + \vec{u}\cdot\nabla\vec{u}+\frac{1}{\rho}\nabla p = \vec{g}+\upsilon \nabla \cdot \nabla \vec{u} \]

\[\nabla \cdot \nabla \vec{u} = 0 \]

該方程表示的是流體的動量方程，本質上是一個牛頓第二定律的形式，描述流體上的力是如何影響流體的運動的，流體粒子所受的力包括粒子的重力、粒子之間所給與的壓力、流體粒子之間相互運動產生的粘滯阻力。

Navier-Stokes方程是一個二階非線性偏微分方程，求解起來非常復雜，在模擬領域使用計算機對NS方程進行求解，通常將方程分解為幾個部分，然后按順序推進。

拉格朗日視角和歐拉視角

拉格朗日視角和歐拉視角是研究流體和可變形固體運動的兩大常用方法，其中拉格朗日視角一般對應粒子法，而歐拉視角一般對應網格法，近年來也出現了混合歐拉拉格朗日視角的方法。

拉格朗日視角基於粒子，將流體的數據，比如密度、質量、速度全部都記錄在一個具體的粒子上，粒子會隨着模擬的推進而改變位置，這個粒子就相當於一個數據采樣點，更加具體一點理解，可以將每個粒子當成一個流體微團，拉格朗日方法模擬的是流體微團之間的相互作用得到的運動。這是一種對物質本身進行離散的NS方程求解方法，這種方法有很多的好處，比如：理解起來直觀、模擬中的Advection部分是易於守恆的，但是拉格朗日視角下有不擅長的點，其中最大的就是粒子的鄰居查找，這塊需要使用空間的數據結構，每次執行算法前都需要更新粒子的鄰居結構，這塊查找鄰居的時間一般是算法的瓶頸。

歐拉視角基於網格，將流體的數據記錄在一組固定的空間點上面，這個空間點上的屬性的含義表示通過這個空間點的屬性是多少，空間點位置是不會隨着模擬的推進而改變的。這是一種對空間進行離散的NS方程求解方法，這種做法的好處就是他在處理整個場的數據是自然離散到了空間點上，在做鄰居查找時也十分的遍歷，可以通過固定的點加上偏移得到周圍空間點的數據，但是歐拉視角在處理粒子運動時經常會出現能量耗散現象，導致整個流體看起來粘性很大，他並不想拉格朗日視角處理Advection那樣的自然的守恆。

兩視角的主要區別就在於他們的采樣點是否隨着模擬推進而運動，胡淵鳴說：拉格朗日視角是隨波逐流，歐拉視角是巋然不動。

近年來最火的是MPM物質點法，他屬於混合歐拉拉格朗日方法，但是筆者沒有具體了解過所以就不多說了。

SPH

原理及公式

Smoothed-Particle Hydrodynamics光滑粒子流體動力學是一種用來模擬連續介質動力學的離散方法，他屬於拉格朗日法，最初用於天體物理領域，近年來在計算機圖形學領域熱度逐漸提高，常用來模擬流體運動。由於其無網格的特性，SPH在模擬復雜邊界、扭曲拉伸、湍流上具有天然的優勢。

SPH的高層思想就是使用一堆隨着模擬運動而更新位置的粒子去攜帶物理量，使用核函數來近似一個連續的場，某一個位置的物理量的值等於周圍粒子的物理量的加權平均，核函數就是用來決定權值，距離該位置越近的權值越大。

變種非常多，比如：PBF、PCI-SPH、WCSPH等等，這里以WCSPH的算法作引入：

其核心的表述求一點的屬性的公式為，A表示一個抽象屬性，可以是密度、壓力、速度等物理量：

\[A(x) = \sum _i A_i \frac{m_i}{\rho_i}W(||x-x_j||_2,h) \]

除了屬性的值還有計算屬性的梯度的方法，這一點上他的梯度不是由上面的公式推導而來的，而是自己做的一個近似公式，但這條公式是對稱的，使用時可以保證動量守恆，所以大家還是接受這條求梯度的公式：

\[\nabla A_i =\rho_i \sum _j {(\frac{A_i}{\rho_i^2}+\frac{A_j}{\rho_j^2}) \nabla_{x_i} W(||x_i-x_j||_2,h))} \]

有了這兩條公式還需要流體相關的公式，需要計算流體粒子的加速度，還有模擬受到的各種粘性力、表面張力等物理定律描述的力，粒子的加速度使用材料導數，材料導數包括了兩項，粒子速度隨時間的變化和粒子速度在空間上的變化，在流體力學領域速度的材料導數表示為：

\[\frac{Dv}{Dt} = -\frac{1}{\rho}\nabla p +g \]

加速度公式涉及到壓強梯度、密度以及外力g。

計算密度就是用上面的抽象屬性公式，將密度代入A：

\[\rho(x) = \sum _i {m_i}W(||x-x_j||_2,h) \]

計算壓強的公式使用了Tait方程，是一個理想氣體狀態方程的變種，不過實際使用還是用的SPH梯度算法來求的

\[p = B((\frac{\rho}{\rho_0})^\gamma - 1) \]

B：bulk modulus，$ \rho_0$ 是一個理想密度，密度很大壓強就很大，壓強很高就會嘗試推着周圍的粒子離開自己。

過程

有了這些公式，現在就可按照一個模擬的積分公式來確定整個SPH的算法流程：

1、對每個粒子計算它所在位置的密度；

2、對每個粒子計算它的壓強梯度，進而得到它的Dv/Dt（其實就相當於求一個力）；

3、Symplectic Euler Step半隱式歐拉方法推進模擬，先根據算出來的力更新速度，然后用更新過的速度來更新位移。

如此一來SPH的流程遍完成了，這是WCSPH的做法，半隱式的時間積分經常需要考慮CFL條件問題，從材料剛度和粒子運動速度方面考慮有時間步長的限制，CFL條件是一個簡單而直觀的判斷計算是否收斂的必要條件。它的比較直觀的物理解釋就是時間推進求解的速度必須大於物理擾動傳播的速度，因為只有這樣才能將物理上所有的擾動俘獲到。不滿足的話整個模擬容易發生數值爆炸等問題，滿足CFL條件意味着當單位時間位移值和時間步長趨於取極限時，數值計算所求的解就會收斂到原微分方程的解。（有關CFL條件，筆者現在也不是理解很透徹，在代碼里，他主要是使得我們的代碼是變步長時間積分的，而不是定步長的）

加速

SPH算法在使用過程中一般最大瓶頸就在粒子的鄰居搜索，因為粒子隨模擬運動而運動，鄰居在每次算法執行前都需要進行更新，要精確的找到一個粒子的周圍粒子這個直接暴力的時間復雜度是O（n^2），為了降低這個瓶頸帶來的影響加速SPH算法，一般會使用特殊的數據結構來存儲粒子鄰居，比如一個將位置信息通過哈希算法投射到一個數組中，一般的步驟是先將連續的粒子位置信息離散為幾個大塊，分成格子，然后做格子位置到粒子數組的映射，格子之間訪問鄰居非常簡單只需要加上偏移量。這樣一來一個粒子周圍的粒子就等於他所在的格子的所有粒子加上周圍一圈格子的所有粒子，每次算法結束步驟時進行維護，根據位置再對格子里的粒子做一個更新。

SPH鄰居搜索的加速方法是很多，比如使用CPU並行進行鄰居搜索、使用CUDA等GPU計算來加速。

引用

https://yangwc.com/2019/05/01/fluidSimulation/

https://diglib.eg.org/bitstream/handle/10.2312/PE.PG.PG2012short.029-034/029-034.pdf

https://book.douban.com/subject/35287308/

2021.8.25