基於SPH的流體模擬過程
轉載自：https://blog.csdn.net/changbaolong/article/details/13172079
http://blog.sina.com.cn/s/blog_6f638fb60100shw0.html
SPH的流體模擬是目前大多數遊戲所採用的模擬流體方法，特點是簡單，十分容易實現，相比與基於Grid的Eulerian方法更加簡單和高速，本文主要介紹一下使用SPH的流體模擬中一些常用的技巧和資料結構。
目前流體模擬中常用的2類方法, 分別代表了從2種不同的方面來解釋Navier-Stokes的流體方程：
1、Eulerian方法
2、Lagrangian方法則將液體看作是跟隨著流動的Particle。
Eulerian方法比較複雜，常用作離線模擬，可以產生非常逼真的流體效果，具體可以參考SIGGRAPH 2007的《FLUID SIMULATION》，作者是Bridson 和 Muller。
現在遊戲中常用的實時流體模擬通常使用Lagrangian的方法來模擬流體，也就是基於粒子的方法，而大多數採用的都是基於SPH(Smooth Particle hydrodynamic)的方法。主要參考文章為Muller的《Particle-Based Fluid Simulation for Interactive Applications
SPH方法非常簡單和容易實現，從我自己做過的剛體和柔體模擬的經驗來看，可以說基於SPH的流體比剛體和柔體的模擬簡單和容易實現的多。SPH流體模擬主要包括3部分：
1、水分子模擬
2、水面資訊抽取
3、水體渲染
其中第一部分，第三部分比較簡單，但是好的快速的效果需要有良好、真實地水面資訊抽取演算法支援。Marching Cube演算法可以通過加大網格精度來不斷提高視覺效果，速度卻得不到保障，所以需要一個更加出色的演算法來解決，因此如何開發一個高效、真實的適用於流體的Surface Extraction演算法會是個很大的挑戰。
（一） SPH方法概括的說：
1：對每個粒子，通過它附近周圍的所有和它的距離半徑小於r的所有粒子來計算出該粒子點的密度r被稱為Smooth Length，通俗的說就是每個粒子的最大影響半徑。
2：對於每一個粒子，通過理想氣體狀態方程根據上一步計算出的密度從而算出它的壓強。
3：對於每一個粒子，根據它附近周圍距離半徑小於r的所有粒子的壓強（上一步已經計算出）來計算出該粒子由於附近壓強差而導致受到的壓力差。
4：對於每一個粒子根據它附近周圍距離半徑小於r的所有粒子的速度差來計算出該粒子由於附近速度差而導致受到的粘滯力。
注： (1)對於第一步求解密度來說，其中r代表粒子i和粒子j之間的距離，h代表前面說的Smooth Length
(2)對於第二步來說使用公式 其中 是粒子的密度 是流體的穩定密度比如水的話一般取1000kg/立方米。 (3) 第3，4步也可以依次類推。具體請參考Muller的文章《Particle-Based Fluid Simulation for Interactive Applications》。實現中的常用技巧
1、推導Smooth Kernel函式的梯度和拉普拉辛運算元
注意，根據上面一篇文章的SPH方法，需要推導Smooth Kernel函式的梯度和拉普拉辛運算元。但Muller在文章中只給出了3個Kernel(一個計算密度，一個計算壓強差力，一個計算粘滯力)函式，並沒有幫我們推匯出後2個函式的梯度公式和拉普拉辛運算元公式。一開始我還花時間自己推導了下，還反覆檢查生怕推導錯了(以前大學學數學我的粗心導致1，2個符號錯誤是常有的事情)。後來發現《SPH survival kit》這片小文章裡已經幫你都推導好了。。。汗。所以你只要照抄就可以了。
2、空間Grid的資料結構
為了快速定位一個粒子周圍的例子，我使用了空間Grid的資料結構，通常的說就是將空間等分成尺寸為h的立方格，這樣對於一個粒子來說只要查詢周圍3X3X3的立方格中所包含的粒子就可以了。
為了實現無限大範圍的Grid系統，我使用了空間稀疏雜湊的方法來儲存Grid,這是因為空間大部分Grid都是空著的所以使用Hash表儲存和快速定位那些包含粒子的Grid，具體可以參考《Optimized Spatial Hashing for Collision Detection of Deformable Objects》。
每對粒子之間的力只需要計算一次就可以了（這樣提高2倍的速度），也就是說對於粒子A,B來說。A粒子由於受到B粒子作用的壓強力分量，和B粒子由於受到 A粒子作用的壓強力分量是相同的只差一個正負號而已，為了避免計算A和B所受的作用力時重複計算，我在每一楨開始模擬前，首先構造每個粒子的鄰接粒子列表：同時要保證A和B的鄰接關係要麼只出現在A的鄰接列表中，要麼只出現在B的鄰接列表中。具體做法很簡單，首先清空Hash Grid中每個Grid中的粒子資訊，然後對於每一個粒子首先查詢周圍27個鄰接Grid中的粒子，將這些粒子中和自己的距離長度小於Kernel Smooth Length的那些例子加入自己的鄰接列表，最後將自己加入到自己所在的Grid中，這樣就能保證A和B的鄰接關係出現並且只出現一次（將會出現在後面的那個粒子的鄰接列表中，因為對於前一個粒子來說，在統計它的鄰接粒子的時候，後一個粒子還沒有被加入到Hash Grid中~~）。
實際上Muller給出的那些公式中有不少可以提取公因式的部分，你可以仔細去研究下會發現很多可以優化速度的地方。
（二）水面資訊抽取
上面的步驟只是完成了一個粒子系統的模擬，這樣你看到的還只是一堆水分子，在正式渲染的時候需要通過第一步計算出的密度場資訊來構造出水錶面。我這裡也是採用Muller論文裡的Marching Cube演算法，這個演算法簡單容易實現，但是效率極低，我模擬1500個粒子在一臺Core2 4400的Intel cpu上可以達到150fps，一旦經過Marching Cube演算法後一下子掉到52fps.所以不大推薦這種方法，我個人看到過一片比較好的文章是《Interactive Water Streams with Sphere Scan Conversion》,感覺挺不錯，作者號稱速度比Marching Cube快了很多倍。不過我自己也沒有時間去做。
Marching Cube演算法的論文<<Marching Cubes:A High Resolutin 3D Surface Construction Algorithm>>,注意論文中提到的構造一個頂點資訊構造表，這個表你可以自己寫程式構造。。但是更加方便的技巧是你去下載一個叫做 “PolyVox”的開源體素構造庫，它裡面已經幫你把表格生成好了。
（三）水體渲染
基本上需要模擬的包括水對於背景的折射效果，水面的Specular ,以及菲涅爾效果。其中前2項是最重要的，折射我使用的方法是先渲染水體以外的所以其他場景物體，然後將該Back Buffer拷到一張渲染紋理上，然後在渲染水的時候在PS裡面對這張紋理進行一些扭曲（根據法向量）。最後將畫素點寫回back buffer,同時根據計算出的菲涅爾係數將畫素點的顏色值和反射顏色值進行融合。（這個折射方法是同事教的，正好他在做畫面扭曲效果。。）Specular的話和一般的Specular計算方法相同。