程志宇 ，徐國(guó)慶*，張嵐斌 ，許 犇

1.武漢工程大學(xué)計(jì)算機(jī)科學(xué)與工程學(xué)院，湖北 武漢 430205；

2.智能機(jī)器人湖北省重點(diǎn)實(shí)驗(yàn)室（武漢工程大學(xué)），湖北 武漢 430205

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

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

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

1 流體建模算法

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

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

1.1 流體運(yùn)動(dòng)方程

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

其中：v為速度，p為壓力，ρ為密度，μ為動(dòng)力學(xué)黏性系數(shù)為流體高度，g為重力加速度，可引入人工黏性項(xiàng)∏來(lái)修正運(yùn)動(dòng)方程。

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

1.2 SPH基本方程

在SPH方法中，問(wèn)題域是用具有質(zhì)量和體積大小的有限粒子表示，通過(guò)對(duì)粒子物理屬性插值計(jì)算區(qū)域任一點(diǎn)的對(duì)應(yīng)值：

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

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

2 改進(jìn)算法

2.1 粒子建模

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

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

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

最后計(jì)算得出流體的加速度a，進(jìn)行速度場(chǎng)的更新。

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

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

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

2.2 自由表面建模和渲染

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

3 結(jié)果與討論

3.1 實(shí)驗(yàn)結(jié)果

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

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

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

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

圖3 SPH運(yùn)動(dòng)模型：（a）球體模型，（b）流體模型Fig.3 Models of SPH motion：（a）spherical model，（b）fluid model

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

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

3.2 對(duì)比分析

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

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

4 結(jié) 語(yǔ)

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