SPH真實感流體交互模擬的改進算法

程志宇 ，徐國慶*，張嵐斌 ，許 犇

1.武漢工程大學計算機科學與工程學院，湖北 武漢 430205；

2.智能機器人湖北省重點實驗室（武漢工程大學），湖北 武漢 430205

真實感流體模擬是計算機圖像學領域中的研究熱點，不管是在電視特效和廣告、電子游戲亦或是軍用作戰(zhàn)平臺仿真、虛擬醫(yī)學等領域中都有著廣泛的應用?；谖锢淼牧黧w模擬的計算模型應該包括流體每一時刻的運動方向、受力情況以及在運動過程中流體間的相互作用力等，同時為了更好的提高真實感要增強細節(jié)模擬［1-2］，展現(xiàn)實時流動的視覺效果。

現(xiàn)階段，對于流體模擬主要采用歐拉網格法和拉格朗日粒子法兩種方法來求解計算流體動力學中的控制方程。歐拉法是模擬流體在固定網格單元上的運動，將流體離散成空間中的固定點，空間點的物理屬性包括密度、速度、壓強等，通過固定點的屬性變化，來模擬流體的運動過程，基于網格的方法可以非常有效地模擬蒸汽和煙霧這樣的體積效應［3-4］。而基于粒子的拉格朗日法則是跟蹤單個流體粒子的運動，并記錄其物理屬性變化，來完成對流體的模擬。拉格朗日粒子法在保持自由表面的質量或體積方面有明顯優(yōu)勢。

基于上述分析，針對真實感流體具有的復雜動力學特性，提出一種在基于通用光滑粒子流體動力學（smoothed particle hydrodynamics，SPH）方法上建立的流體模型，利用粒子模擬流體運動過程，SPH初始化完成粒子創(chuàng)建，相鄰粒子根據(jù)劃分的網格單元進行交互。根據(jù)粒子的空間位置構建流體密度場計算粒子密度、壓力、粘滯力、外力等，采用固定鏡像粒子方法［5］處理交互邊界問題［6-7］。通過移動立方體（marching cubes，MC）算法來重建流體網格［8-10］，實現(xiàn)對流體表面建模，光照計算實現(xiàn)反射和折射效果，完成對流體的高質量交互式渲染［11］，形成更逼真的動態(tài)演示效果。

1 流體建模算法

SPH方法是模擬流體流動的一種拉格朗日形式的無網格粒子法。它的理論基礎是粒子方法，由于粒子的運動與流體的運動相似，可以用經典牛頓流體動力學的控制方程來模擬。與傳統(tǒng)的流體模擬方法相比，采用SPH方法能夠比較真實的模擬流體運動過程，對于邊界不連續(xù)問題，大面積變形問題，模擬破裂合并以及飛濺等，SPH方法可以很好的解決。

它的核心思想是用一系列任意分布的粒子來表示問題域，粒子之間不需要任何連接，用積分形式來近似場函數(shù)，再通過局部區(qū)域（支持域）的相鄰粒子所攜帶的計算數(shù)值疊加求和來取代場函數(shù)及其導數(shù)的積分表示式，實現(xiàn)進一步近似。在每一個時間步中都要進行粒子近似，因此當前粒子的性質取決于當前局部分布的粒子。

1.1 流體運動方程

基于拉格朗日的流體運動方程可由連續(xù)性方程（1）、運動方程（2）、壓力狀態(tài)方程（3）表示［12］。

其中：v為速度，p為壓力，ρ為密度，μ為動力學黏性系數(shù)為流體高度，g為重力加速度，可引入人工黏性項∏來修正運動方程。

式（4）中：α是0到1之間的常系數(shù)，v是流體速度，r是粒子的位置矢量，η是保證分母不為零的參數(shù)，h是光滑核半徑。顯然，拉格朗日方法是針對于相對獨立的各個流體粒子來求解運動方程，通過積分計算流體粒子下一時刻的物理屬性值變化。

1.2 SPH基本方程

在SPH方法中，問題域是用具有質量和體積大小的有限粒子表示，通過對粒子物理屬性插值計算區(qū)域任一點的對應值：

其中：A（r）是位置r的物理屬性；N為粒子j在支持域中的粒子總數(shù)；mj是粒子j的質量；ρj是粒子j的密度；函數(shù)W為光滑核函數(shù)，h為光滑核半徑，光滑核函數(shù)具有2個性質：首先必須是偶函數(shù)，滿足W（-r）=W（r）；其次，必須為規(guī)整函數(shù)，滿足

在SPH方法中，導數(shù)只影響光滑核。光滑屬性函數(shù)A(r)的梯度以及拉普拉斯算子：

2 改進算法

2.1 粒子建模

流體方程計算時，首先在密度場應用式（5）計算粒子的求和密度ρ，使用Poly6核函數(shù)估計：

將式（6）和式（7）代入運動方程并進行對稱化，得到作用在粒子上的壓力和黏度力位于不同壓強區(qū)的粒子間的相互作用力不等，所以采用相互作用粒子壓強的算術平均值來優(yōu)化計算式。

其中，基于計算精度和效率，估計計算選擇核函數(shù)形式如下：

最后計算得出流體的加速度a，進行速度場的更新。

模擬流體動態(tài)效果的SPH流體粒子模型算法的基本流程如圖1所示。

圖1 SPH粒子系統(tǒng)建模流程Fig.1 Modeling flowchart of SPH particle system

此SPH模型適合用來求解高速碰撞過程中的流體大變形的問題，通過不斷循環(huán)迭代修正密度差［13-14］，對流體的形狀和運動狀態(tài)都能模擬較真實效果，模擬具有有效性和高效性。

2.2 自由表面建模和渲染

流體表面重建是流體模擬的重要內容。在未實現(xiàn)自由界面的重構和渲染時，用SPH方法得到的是由一堆離散粒子組成的流體系統(tǒng)，流體的真實感以及運動連貫性欠缺。采用一種基于移動立方體（MC）算法來提取流體運動表面的真實效果，實現(xiàn)流體表面建模。MC算法的核心思想是在三維離散數(shù)據(jù)場中按一定規(guī)則劃分體元，在每個體元內部根據(jù)設定值對每條邊進行插值，通過線性插值得到的三角面片來逼近等值面，完成等值面提取。在動畫模擬過程中圖形處理器（graphics processing unit，GPU）加速進行實時渲染［15］，將流體表面渲染成透明的，并進行光照計算，實現(xiàn)了反射和折射效果。同時流體表面的光照和紋理渲染提高了流體場景的真實感。

3 結果與討論

3.1 實驗結果

實驗平臺是 Intel（R）Core（TM）i7-8550U CPU 1.8 GHZ，主存 16 GB，圖形卡為 NVIDIA GeForce MX150，顯卡內存為2 GB，開發(fā)環(huán)境Visual Studio 2017，圖形的編程接口由OpenGL提供。在圖形硬件上，實驗完成了對復雜流體方程的求解，同時針對于交互結果進行了實時渲染展示。實驗結果顯示了本文方法模擬的有效性和高效性，實時場景渲染的真實感有所提高。

圖2（a）顯示了粒子系統(tǒng)初始形態(tài)，粒子數(shù)約為4 000，室外場景為背景，圖2（b）是渲染后的流體系統(tǒng)，流體采用透明渲染，用MC算法構造了自由移動的流體表面的效果。其中每個粒子包含重力、壓力、粘滯力、交互粒子間的引力與斥力作用。場景包括粒子的物理屬性值計算、自由表面建模、光照計算和材質渲染等。

圖2 SPH模型初始狀態(tài)：（a）粒子模型，（b）流體模型Fig.2 Initial state of SPH models：（a）particle model，（b）fluid model

圖3（a）展示了球體渲染模型在外部輸入力的作用下運動的效果，粒子數(shù)約為4 000。圖3（b）是透明渲染流體在外力作用下的運動效果。在模擬過程中，用戶可以通過鼠標拖拽來調整流體的受力，其中每個粒子除了受到基礎作用力外還受到外部輸入的力的作用。在外力作用下流體表面發(fā)生劇烈變形，形成飛濺效果。

圖3 SPH運動模型：（a）球體模型，（b）流體模型Fig.3 Models of SPH motion：（a）spherical model，（b）fluid model

實驗將木塊運動與實時流體模擬相結合，通過計算木塊下落位置和速度來觸發(fā)與流體的動態(tài)交互，使流體產生相應形變。圖4是流體和木塊交互的動態(tài)演示效果，包括木塊下沉和上浮過程、場景包括相鄰網格粒子搜索和復雜的邊界交互計算。

圖4 流體和木塊交互模型：（a）木塊下沉，（b）木塊上浮Fig.4 Models of fluid and wood interaction：（a）wood sinking，（b）wood floating

3.2 對比分析

小尺度的流體模型渲染速度會受到粒子數(shù)和光滑半徑的影響，實驗過程中通過設置對照組進行對比實驗，不同粒子數(shù)和光滑半徑下的渲染速度關系如表1所示。實驗結果表明光滑半徑相同，隨著粒子數(shù)增加，渲染速度變慢；粒子數(shù)不變，隨著光滑半徑的增加，渲染速度變慢。

表1 不同粒子數(shù)和光滑半徑下的渲染幀率Tab.1 Rendering frame rates of different particle numbers and smoothing radii frames/s

4 結 語

本文主要討論真實感流體的模擬問題，利用基于SPH的數(shù)值計算法進行粒子建模，再采用MC算法對粒子表面建模。實驗結果表明，本文提出的對于真實感流體的模擬改進算法，綜合利用體繪制和面繪制法實現(xiàn)流體三維模擬，真實感提高，優(yōu)化效果明顯。但對于數(shù)值方法的固有計算損耗、細節(jié)特征模擬以及硬件加速方面要做進一步改進，提高流體模擬的效率和精度。