流體模擬：Position Based Fluid
流體模擬：Position Based Fluid
之前部落格中SPH基礎中介紹的WCSPH或者PCISPH介紹的方法都是屬於基於力的方法。這類方法通過計算內力（比如流體內部的粘性力和壓力）和外力（比如重力和碰撞力）的合力，然後根據牛頓第二定律計算出加速度，再根據數值積分求出速度和位置資訊。
基於位置動力學（Position Based Dynamics）的方法將這些物理運動通過約束表達出來，這樣只需要一個求解器即可，更加方便地進行物理模擬。下面我們來講一下PBD。
1. 位置動力學（PBD）
我們來看以下基於力和基於位置動力學的物理碰撞更新過程的對比：
可以看到基於力的碰撞檢測首先在穿透發生時更新物體的速度，然後更新物體的位置，會有明顯的反應延遲並且需要調整剛度（stiffness）引數（比較難調）。而基於位置動力學的碰撞檢測首先只檢測是否發生穿透，然後移動位置使之不發生穿透，最後再據此更新物體的速度資訊。
1.1 演算法過程
在PBD演算法中，運動的物理由\(N\)個頂點和\(M\)
- 約束的基數為\(n_j\),即第\(j\)個約束所影響的頂點數目為\(n_j\)個
- 約束函式\(C_j\):\(R^{3nj} \rightarrow R\)
- 受約束影響的頂點索引值集合\({i_1,...,i_{nj}},i_k\in [1,...,N]\)
- 每個約束都有對應的剛性引數\(k_j\in [0,1]\),這裡我們可以理解為約束的強度；
- 約束分為兩種，一類是等式約束即\(C_j(x_{i1},x_{i2},...,x_{in_j})=0\)，另一類是不等式約束$C_j(x_{i1},x_{i2},...,x_{in_j})\ge 0 $
給定時間步長\(\Delta t\)，PBD的運動物體模擬的演算法虛擬碼如下所示:
在上面的演算法第1步到第3步中，我們首先對頂點的位置、速度和質量倒數進行初始化，其中質量的倒數\(w_i=1/m_i\),除了可以避免冗餘的除法操作外，還可以使用於靜態的物體，對於靜態的物體我們設為\(w_i=0\)(可以理解為質量無窮大),這樣在後續的更新中都不會產生位置和速度的變化量。
第5步中的\(f_{ext}\)代表不能轉換為約束形式的力（如重力），我們根據\(f_{ext}\)進行一次數值計算預測在\(f_{ext}\)的作用下的速度\(v_i\)。緊接著在第6步中我們新增阻尼的作用，阻尼可以理解為物體在運動中發生了能量耗散，從而導致速度有所衰減。
第8行主要是生成碰撞約束，物體會與周圍的環境發生碰撞，例如布料落在地板上，水碰上一面牆等，這些碰撞約束在每個時間步長都發生改變，所以每一次都需要重新生成碰撞約束。有了內部約束（如不可壓縮流體的密度約束）和外部約束（如流體與地面的碰撞約束）之後，我們需要根據這些約束做一個迭代求解，也就是上面虛擬碼中的第9行到第11行，這裡我們稱為約束投影步驟。從約束投影步驟我們得到服從給定約束的粒子位置，然後再第12行到第15行更新頂點粒子的速度和位置資訊。最後在第16行根據摩擦係數（friction）和恢復係數（restitution）更新速度，如下圖2所示。這樣，一個完整的PBD物理模擬步驟就完成了。
1.2 約束投影步驟
接下來我們就針對約束投影步驟詳細展開相關的內容，約束投影是PBD中的最難理解的核心部分，涉及的數學內容比較多一點。設有一個基數為\(n\)的約束，關聯的粒子點為\(p_i,...,p_n\),約束函式記為\(C\),剛度係數（stiffness）為\(k\),記\(p=\left[p_{1}^{T}, \ldots, p_{n}^{T}\right]^{T}\),則等式約束函式表示為：
\[C(p)=0 \tag{1} \label{1} \]我們希望計算一個位移偏移量\(\Delta p\),使得粒子頂點在\(p+\Delta p\)出約束條件依然滿足，即：
\[C(p+\Delta p)=0 \tag{2} \label{2} \]對約束條件\(C\)做一階泰勒展開，則可得：
\[C(p+\Delta p) \approx C(p)+\nabla_{p} C(p) \cdot \Delta p=0 \tag{3} \label{3} \]為了使粒子在\(p+\Delta p\)出依然滿足約束條件，我們要求方程\((3)\)得到\(\Delta p\)。PBD演算法一個絕妙之處在於他將\(\Delta p\)的方向限制在約束函式的梯度方向\(\nabla C(p)\)上。如下圖所示，約束\(C\)所涉及到了粒子位置回形成一個高維空間，下圖為該空間中滿足不同約束條件的粒子位置形成的二維等值線示意圖，其中滿足\(C\)約束條件的是黑色等值線。故當粒子處於下圖的黑色點的位置時，不滿足約束條件，如果我們沿著點所在的等值線（灰色曲線）移動，此時剛體模態（Rigid body modes）的方向與該等值線相同，新得到的位置仍然在該灰色等值線上，依然不在黑色曲線\(C\)上，即不滿足約束條件。這可以理解為，約束中存在的誤差依然沒有得到修正。以兩個粒子形成的距離約束為例，就好比同時移動了兩個粒子或者該約束繞自身旋轉，但是存在的誤差並沒有得到更正。而且這樣一來還會引入系統中不存在的一種外力，導致系統動量不守恆。所以，我們希望該點的位移方向與剛體模態方向垂直，從而保證系統動量守恆，即從黑點指向紅點的方向\(\nabla C\)
因此，令位移向量\(\nabla p\)為約束函式的梯度向量\(\nabla _pC\)再乘上一個標量縮放係數\(\lambda\)：
\[\Delta p=\lambda \nabla_{p} C(p) \tag{4} \label{4} \]其中的標量縮放係數\(\lambda\)，我們稱之為拉格朗日乘子（Lagrange multiplier）。聯立公式\((3)\)和\((4)\)我們可得：
\[\lambda=-\frac{C(p)}{\left|\nabla_{p} C(p)\right|^{2}} \tag{5} \label{5} \]再將$\lambda \(代回公式\)(4)\(我們可得\)\Delta p$的表示式：
\[\Delta p=\lambda \nabla_{p} C(p)=-\frac{C(p)}{\left|\nabla_{p} C(p)\right|^{2}} \nabla_{p} C(p) \tag{6} \label{6} \]具體到粒子\(i\)，約束投影后其對應的位移向量為：
\[\Delta p_{i}=-s \nabla_{p_{i}} C\left(p_{1}, \ldots, p_{n}\right) \tag{7} \label{7} \]其中的\(s\)如下所示，\(s\)的值對於約束函式\(C\)作用範圍內地額所有點都一樣：
\[s=\frac{C\left(p_{1}, \ldots, p_{n}\right)}{\sum_{j}\left|\nabla_{p_{j}} C\left(p_{1}, \ldots, p_{n}\right)\right|^{2}} \tag{8} \label{8} \]前面我們假定所有的粒子質量都相同，現在考慮粒子質量不同的情況。記粒子\(i\)的質量為\(m_i\),其質量的倒數為\(w_i=1/m_i\),則公式\((4)\)變為：
\[\Delta p_{i}=\lambda w_{i} \nabla_{p_{i}} C(p) \tag{9} \label{9} \]公式\((7)\)和公式\((8)\)變為：
\[\Delta p_{i}=-s w_{i} \nabla_{p_{i}} C\left(p_{1}, \ldots, p_{n}\right) \tag{10} \label{10} \] \[s=\frac{C\left(p_{1}, \ldots, p_{n}\right)}{\Sigma_{j} w_{j}\left|\nabla_{p_{j}} C\left(p_{1}, \ldots, p_{n}\right)\right|^{2}} \tag{11} \label{11} \]為了便於理解，接下來我們舉個簡單的例子應用約束投影方法。如下圖所示。
上面的約束可以表示為\(C\left(p_{1}, p_{2}\right)=\left|p_{1}-p_{2}\right|-d\),位移向量記為\(\Delta p_i\)。根據約束投影方法，我們首先約束函式\(C(p_1,p_2)\)關於\(p_1\)和\(p_2\)的梯度，也就是求偏導。注意到\(C\left(p_{1}, p_{2}\right)=\left|p_{1}-p_{2}\right|-d=\left(\sqrt{\left(p_{1}-p_{2}\right)^{2}}\right)-d\),我們可以求得以下的梯度向量表示式：
\[\begin{gathered} \nabla_{p_{1}} C\left(p_{1}, p_{2}\right)=\frac{p_{1}-p_{2}}{\left|p_{1}-p_{2}\right|} \\ \nabla_{p_{2}} C\left(p_{1}, p_{2}\right)=-\frac{p_{1}-p_{2}}{\left|p_{1}-p_{2}\right|} \end{gathered} \tag{12} \label{12} \]注意，上面求到的是一個向量，也就是我們說的梯度向量。將公式\((12)\)代入公式\((11)\)可得：
\[\begin{aligned} s &=\frac{C\left(p_{1}, \ldots, p_{n}\right)}{\Sigma_{j} w_{j}\left|\nabla_{p_{j}} C\left(p_{1}, \ldots, p_{n}\right)\right|^{2}} \\ &=\frac{\left|p_{1}-p_{2}\right|-d}{w_{1}\left|\nabla_{p_{1}} C\left(p_{1}, p_{2}\right)\right|^{2}+w_{2}\left|\nabla_{p_{2}} C\left(p_{1}, p_{2}\right)\right|^{2}} \\ &=\frac{\left|p_{1}-p_{2}\right|-d}{w_{1}+w_{2}} \end{aligned} \tag{13} \label{13} \]最後，將公式\((13)\)代入到公式\((10)\),可得約束投影計算得到的位移：
\[\begin{aligned} \Delta p_{1} &=-\frac{\left|p_{1}-p_{2}\right|-d}{w_{1}+w_{2}} w_{1} \nabla_{p_{1}} C\left(p_{1}, p_{2}\right) \\ &=-\frac{w_{1}}{w_{1}+w_{2}}\left(\left|p_{1}-p_{2}\right|-d\right) \frac{p_{1}-p_{2}}{\left|p_{1}-p_{2}\right|} \end{aligned} \]同理得\(\Delta p_2\),如下所示：
\[\Delta p_{2}=+\frac{w_{2}}{w_{1}+w_{2}}\left(\left|p_{1}-p_{2}\right|-d\right) \frac{p_{1}-p_{2}}{\left|p_{1}-p_{2}\right|} \]前面我們提到了每個約束都有對應得剛度係數\(k\),令\(k'=1-(1-5)^{1/n_s}\)去乘\(\Delta p\),這裡\(n_s\)迭代之後誤差為\(\Delta p\left(1-k^{\prime}\right)^{n_{s}}=\Delta p(1-k)\),與剛度係數成線性關係，而與迭代次數\(n_s\)無關。下一個時間步的位置如下所示：
\[\begin{aligned} &p_{1}^{t+1}=p_{1}^{t}+k^{\prime} \Delta p_{1} \\ &p_{2}^{t+1}=p_{2}^{t}+k^{\prime} \Delta p_{2} \end{aligned} \]1.3 約束投影求解器
前面的虛擬碼中我們可以看到約束投影的輸入為\(M+M_{coll}\)個約束和\(N\)個點的預測位置\(p_1,...,p_N\),所需要求解的方程組是非線性非對稱方程組或不等式組（碰撞約束產生的）。約束投影步驟的主要任務就是修正預測位置使新得到的校正位置滿足所有約束。但是一般情況下很難找到一個適當的\(\Delta p=\left[\Delta p_{1}^{T}, \ldots, \Delta p_{n}^{T}\right]^{T}\)恰好使得所有的約束都能夠同時得到滿足，故我們通常採用迭代的方法按順序依次對約束進行求解。
我們可以採用非線性高斯-賽德爾（Non-Linear Gauss-Seidel，簡稱NGS）迭代方法。高斯賽德爾（Gauss-Sedel，簡稱GS）迭代方法只能求解線性方程組，NGS在依次求解德基礎上，加入了約束投影求解這一非線性操作。與雅可比迭代方法（Jacobi method）不同，NGS求解器在一次迭代中對於頂點位置的修正立即被應用到下一個約束求解中，這樣的好處就是顯著加快了收斂速度。
但是NGS雖然穩定且容易實現，但是該方法收斂速度依然不是很快，不宜並行化。而且取決於約束求解順序，且不宜並行化。
2 基於位置動力學的流體模擬（PBF）
前面部分主要介紹了Position Based Dynamics演算法相關的內容，接下來我們就看看如何將其PBD演算法應用到流體模擬當中。
2.1 不可壓縮約束和拉格朗日乘子
在不可壓縮流體模擬中，我們希望粒子\(i\)的密度儘量與靜態密度\(\rho_0\)相同，記\(\rho_i=\rho_0\),因此，我們需要對每一個流體粒子都施加一個常量密度約束，PBF的常量密度約束如下所示：
\[C_{i}\left(p_{1}, \ldots, p_{n}\right)=\frac{\rho_{i}}{\rho_{0}}-1 \tag{14} \label{14} \]我們記粒子\(i\)的位置為\(p_i\),\(p_i,...,p_n\)使與粒子\(i\)相鄰的粒子。我們可以看到當密度約束\(=0\)時，有\(\rho_i=\rho_0\),此時流體的體積既不壓縮也不膨脹，從而保證了流體的不可壓縮條件。
流體粒子\(i\)的密度根據SPH方法的插值計算得到：
\[\rho_{i}=\sum_{j} m_{j} W\left(p_{i}-p_{j}, h\right) \tag{15} \label{15} \]在公式\((15)\)中，\(m_j\)時鄰居粒子\(j\)的質量，\(h\)是指定的光滑核半徑。將公式\((15)\)代入公式\((14)\)，我們有：
\[C_{i}\left(p_{1}, \ldots, p_{n}\right)=\frac{\Sigma_{j} m_{j} W\left(p_{i}-p_{j}, h\right)}{\rho_{0}}-1 \tag{16} \label{16} \]在公式\((15)\)的密度計算中，PBF方法採用Poly6核函式：
\[W_{p o l y 6}(r, h)=\frac{315}{64 \pi h^{3}} \begin{cases}\left(1-\frac{r^{2}}{h^{2}}\right)^{3} & 0 \leq r \leq h \\ 0 & \text { otherwise }\end{cases} \tag{17} \label{17} \]但是在計算密度的梯度時，卻又採用了Spiky核函式：
\[W_{spkiy}(r, h)=\frac{15}{\pi h^{6}} \begin{cases}\left(1-\frac{r}{h}\right)^{3} & 0 \leq r \leq h \\ 0 & \text { otherwise }\end{cases} \tag{18} \label{18} \]對公式\((18)\)求關於\(r\)的導數（注意，\(|r|=\sqrt{r^{2}}\)，不能直接對\(|r|\)求導），從而流體粒子密度的梯度如下所示：
\[\nabla W_{\text {spiky }}(r, h)=-\frac{45}{\pi h^{6}} \begin{cases}(h-|r|)^{2} \frac{r}{|r|} & 0 \leq|r| \leq h \\ 0 & \text { otherwise }\end{cases} \tag{19} \label{19} \]因此，粒子\(i\)的約束函式\((16)\)是一個關於\(p_1,...,p_n\)的非線性方程組，\(C_{i}\left(p_{1}, \ldots, p_{n}\right)\)=0,所有粒子的約束組成了一個非線性方程組。在PBF方法中，我們只考慮粒子質量相同的情況，我們省去公式\((15)\)和\((16)\)中的質量\(m_j\),即：
\[\begin{eqnarray*} \rho_{i}=\Sigma_{j} W\left(p_{i}-p_{j}, h\right) \tag{20} \label{20} \\ C_{i}\left(p_{1}, \ldots, p_{n}\right)=\frac{\Sigma_{j} W\left(p_{i}-p_{j}, h\right)}{\rho_{0}}-1 \tag{21} \label{21} \end{eqnarray*} \]然後求約束函式\(C_i\)關於\(p_k\)的梯度如下,其中\(k\in {1,2,...,n}\):
\[\nabla_{p_{k}} C_{i}=\frac{1}{\rho_{0}} \Sigma_{j} \nabla_{p_{k}} W\left(p_{i}-p_{j}, h\right) \tag{22} \label{22} \]顯然，針對\(k\)的不同，分為兩種情況。當\(k=i\)也就是粒子本身時候，連加符號中\(W\)均為關於\(p_k\)的函式；當\(k=j\)即鄰居粒子的時候，只有\(W(p_i-p_k,h)\)才有意義，其他相對於\(p_k\)來說都是常量，故導數為0（注意用到了求導的鏈式法則）：
\[\nabla_{p k} C_{i}=\frac{1}{\rho_{0}} \begin{cases}\Sigma_{j} \nabla_{p k} W\left(p_{i}-p_{j}, h\right) & \text { if } k=i \\ -\nabla_{p k} W\left(p_{i}-p_{j}, h\right) & \text { if } k=j\end{cases} \tag{23} \label{23} \]既然求出了約束函式的梯度，我們就把它應用到前面提到的拉格朗日乘子的計算公式中，聯立公式\(*(5)\)和\((23)\)，我們可得拉格朗日乘子：
\[\lambda_{i}=-\frac{C_{i}\left(p_{1}, \ldots, p_{n}\right)}{\Sigma_{k}\left|\nabla_{p k} C_{i}\right|^{2}} \tag{24} \label{24} \]2.2 軟約束（Soft constraint）與混合約束力（Constraint Force Mixing, CFM）
如果一個約束條件不能被違背，我們稱之為硬約束；而能一定程度上被違背的約束稱為軟約束。在理想的情況下，我們都希望約束始終是硬約束，但是由於誤差或者數值方法的不穩定等原因，我們有時不得不向軟約束妥協。
在\(PBF\)中，當\(|r|=h\),粒子\(i\)與粒子\(j\)之間的距離等於光滑核半徑時，粒子\(i\)和粒子\(j\)處於即將分離的狀態。注意觀察公式\((19)\)的密度梯度計算公式，此時\(\nabla W_{spiky}(r,h)=0\)。若所有的鄰居粒子與粒子\(i\)都處於這種狀態，那麼必將導致約束函式的梯度即公式\((22)\)取值為0：
\[\nabla_{p k} C_{i}=\frac{1}{\rho_{0}} \Sigma_{j} \nabla_{p k} W\left(p_{i}-p_{j}, h\right)=0 \]從而導致公式\((24)\)中的分母\(\Sigma_{k}\left|\nabla_{p k} C_{i}\right|^{2}\)為0，出現除零錯誤，這將導致PBF方法出現潛在的不穩定性。為了解決這個問題，PBF採用混合約束的方法，使密度硬約束轉變成軟約束。具體的做法就是將根據密度函式求解得到的約束力在加入到原始的約束函式中，這裡的在PBF的常量密度約束中得到的拉格朗日乘子\(\lambda\)又類似的作用，故將\(\lambda\)加入到初始的約束方程（即公式\((3)\)）:
\[C(p+\Delta p) \approx C(p)+\nabla C^{T} \nabla C \lambda+\epsilon \lambda=0 \tag{25} \label{25} \]公式\((25)\)中的$\epsilon \(是鬆弛引數，可以由使用者指定。引入公式\)(25)\(後，拉格朗日乘子的計算公式\)(24)$就變成：
\[\lambda_{i}=-\frac{C_{i}\left(p_{1}, \ldots, p_{n}\right)}{\Sigma_{k}\left|\nabla_{p k} C_{i}\right|^{2}+\epsilon} \tag{26} \label{26} \]從而可得粒子\(i\)在經過上述約束投影后對應的位移向量（包括自身密度約束以及鄰居粒子密度約束共同作用的結果。注意，這裡對應上面的公式\((4)\)，結合公式\((23)\)):
\[\begin{aligned} \Delta p_{i} &=\lambda_{i} \nabla_{p_{i}} C_{i}+\Sigma_{j} \lambda_{j} \nabla_{p_{j}} C_{i} \\ &=\frac{1}{\rho}_{0} \Sigma_{j} \lambda_{i} \nabla_{p_{i}} W(r, h)+\left(-\frac{1}{\Sigma_{j}} \lambda_{j} \nabla_{p_{j}} W(r, h)\right) \\ &=\frac{1}{\rho_{0}} \Sigma_{j} \lambda_{i} \nabla_{p_{i}} W(r, h)+\frac{1}{\rho_{0}} \Sigma_{j} \lambda_{j} \nabla_{p_{i}} W(r, h) \\ &=\frac{1}{\rho_{0}} \Sigma_{j}\left(\lambda_{i}+\lambda_{j}\right) \nabla_{p_{i}} W(r, h) \end{aligned} \tag{27} \label{27} \]2.3 拉伸不穩定性
PBF採用SPH的方法來計算流體粒子的密度，但是該方法通常需要30~40個鄰居粒子才能使密度求值結果趨於靜態密度。在鄰居粒子數量較少的情況下，通過該方法計算得到的流體密度低於靜態密度，由此會造成流體內部壓強為負數，原本粒子間的壓力變為吸引力，使得流體粒子凝聚在一起，導致流體表面的模擬效果失去真。PBF採用了一種人工排斥力的計算模型，當流體粒子距離過近時該排斥力會使它們分開，避免產生粒子聚集的現象。如下圖所示：
在公式\((24)\)的基礎上，加入一個排斥項（repulsice term）\(S_{corr}\):
\[\Delta p_{i}=\frac{1}{\rho_{0}} \Sigma_{j}\left(\lambda_{i}+\lambda_{j}+s_{\text {corr }}\right) \nabla_{p_{j}} W\left(p_{i}-p_{j}, h\right) \tag{28} \label{28} \]其中的\(S_{corr}\)計算方式如下：
\[s_{c o r r}=-k\left(\frac{W\left(p_{i}-p_{j}, h\right)}{W(\Delta q, h)}\right)^{n} \tag{29} \label{29} \]公式\((29)\)中，\(\Delta q\)表示到粒子\(i\)的一個固定距離，通常取\(|\Delta q|=0.1h,...,0.3h\),\(h\)即前面提到的光滑核半徑。此外，公式中的\(k\)可以看作表面張力引數，取值\(k=0.1\),而\(n=4\)。公式\((28)\)中的排斥項會使得流體粒子的密度稍微低於靜態密度，從而產生類似於表面張力的效果，使得流體表面的的粒子分佈均勻。通過這個排斥項，我們不再需要硬性規定流體的鄰居數量必須在30~40個，進一步提升演算法的流體模擬效率。
2.4 渦輪控制和人工粘性
由於數值耗散，PBD的方法會引入額外的阻尼，使得整個系統的能量損耗太快，導致本來應該本來一些渦流細節迅速消失。在這裡，PBF通過渦輪控制方法向整個系統重新注入能量：
\[f_{i}^{\text {vorticity }}=\epsilon\left(N \times \omega_{i}\right) \tag{30} \label{30} \]上述公式中，\(N=\frac{\eta}{|\eta|}, \eta=\nabla|\omega|_{i}\),而流體粒子的旋度\(w_i\)計算公式如下：
\[\omega_{i}=\nabla \times v=\Sigma_{j}\left(v_{j}-v_{i}\right) \times \nabla_{p_{j}} W\left(p_{i}-p_{j}, h\right) \tag{31} \label{31} \] \[\eta=\sum_{j} \frac{m_{j}}{\rho_{j}}\left\|\omega_{j}\right\| \nabla_{i} W_{i j} \]渦輪控制方法的基本思路就是：通過新增一個體積力\(f_{i}^{\text {vorticity }}\)(在演算法的第一步)，在旋度粒子（可直觀理解為比周圍粒子旋轉快的粒子，旋度\(w_i\)指向粒子\(i\)的旋轉軸）處加速例子的旋轉運動，通過這種方式來增加系統的旋度細節。公式\((30)\)中的\(\epsilon\)用於控制渦輪控制力的強度。
最後，PBF方法採用XSPH的粘度方法直接更新速度，從而產生粘性阻尼。人工粘性除了可以增加模擬的數值穩定性，還可以消除非物理的流體振盪。拉格朗日流體模擬方法中，人工粘性本質上會對流體粒子的相對運動產生阻尼作用，使流體的動能轉化為熱能：
\[v_{i}^{\text {new }}=v_{i}+c \Sigma_{j}\left(v_{i}-v_{j}\right) \cdot W\left(p_{i}-p_{j}, h\right) \tag{32} \label{32} \]在流體模擬中，我們取公式\((32)\)中的\(c=0.01\)
3. PBF演算法實現
PBF演算法的總體框架就是按照前面提到的PBD演算法，只是經典PBD演算法採用了順序高斯-賽德爾（Sequential Gauss-Seidel，SGS）迭代求解，而SGS不容易被GPU並行化，因此基於CUDA實現的PBF求解器使用了雅克比（Jacobi）迭代方法並行求解。
參考資料
【深入淺出 Nvidia FleX】(1) Position Based Dynamics