基于SPH 統(tǒng)一框架的流體細(xì)節(jié)仿真模擬方法?

華 磊 曾蘭玲 楊 洋 趙 巖

（江蘇大學(xué)計(jì)算機(jī)科學(xué)與通信工程學(xué)院 鎮(zhèn)江 212000）

1 引言

流體模擬在計(jì)算機(jī)游戲、電影特效和軍事訓(xùn)練等領(lǐng)域有廣泛的需求和應(yīng)用。經(jīng)過(guò)相關(guān)學(xué)者和專家的不斷努力，流體模擬已經(jīng)取得了很多研究成果，并應(yīng)用在不同的領(lǐng)域。但是現(xiàn)階段研究方向都是針對(duì)單一流體的仿真模擬，因此如何通過(guò)統(tǒng)一框架模擬不同流體的運(yùn)動(dòng)狀態(tài)非常具體挑戰(zhàn)性。此外，在具體細(xì)節(jié)真實(shí)感方面，如火焰、煙霧的邊緣細(xì)節(jié)優(yōu)化還有待進(jìn)一步提高。

基于粒子的方法［1~3］和基于網(wǎng)格的方法［4］是目前流體模擬最流行的離散化不同層次的典型方案。基于網(wǎng)格的方法假設(shè)在無(wú)發(fā)散條件下具有不可壓縮性，在拓?fù)滢D(zhuǎn)換上更有效率，但是無(wú)法精確有效的處理場(chǎng)的流動(dòng)，容易丟失能量和幾何合效果。SPH（Smoothed Particle Hydrodynamics）光滑顆粒的流體動(dòng)力學(xué)［5］算法就是典型的拉格朗日粒子法，其基本原理是通過(guò)粒子模擬流體的運(yùn)動(dòng)規(guī)律，然后轉(zhuǎn)換成網(wǎng)格進(jìn)行流體渲染。如何在計(jì)算成本和真實(shí)細(xì)節(jié)之間做出的合理的選擇是目前流體模擬的研究熱點(diǎn)［6~10］，但該方面的工作基于的都是單一的流體性能方面的模擬。其中文獻(xiàn)［11~14］的研究工作是采用預(yù)先存儲(chǔ)的數(shù)據(jù)，通過(guò)不同狀態(tài)的組合生成最后的效果，文獻(xiàn)［15］是通過(guò)加快鄰域粒子的搜索速度來(lái)提高計(jì)算機(jī)模擬性能。針對(duì)多樣化流體模擬研究方面的工作較少，因此，實(shí)現(xiàn)多樣化流體模擬的統(tǒng)一框架具有實(shí)際應(yīng)用意義。

由于SPH 可以有效地表示小的液體特性與低內(nèi)存成本等特點(diǎn)，本文研究如何利用SPH來(lái)實(shí)行多樣化的流體效果，并通過(guò)SPH中的光滑核函數(shù)的優(yōu)化，實(shí)現(xiàn)一個(gè)控制流體多樣化的統(tǒng)一框架。細(xì)節(jié)模擬方面，由于SPH 對(duì)粒子分布很敏感，在實(shí)現(xiàn)水波紋和煙霧波動(dòng)等細(xì)節(jié)效果，可以選擇在邊界層粒子的位置進(jìn)行細(xì)節(jié)化處理。目前的相關(guān)研究工作主要采用毛細(xì)管波和重力波模擬流體中的波紋細(xì)節(jié)效果，但毛細(xì)管波速度快、高頻率難以模擬，而重力波在細(xì)節(jié)處理方粗糙，因此，文獻(xiàn)［16~17］研究在虛擬環(huán)境中模擬聲波的傳播。由于聲波速度快，可以在流體的表面模擬聲波，且不存在數(shù)值化問(wèn)題等優(yōu)點(diǎn)，本文采用聲波，即壓縮波對(duì)流體邊界層細(xì)節(jié)進(jìn)行優(yōu)化處理。該工作背后的基本思想是基于SPH結(jié)合聲波的方法，實(shí)現(xiàn)豐富多樣化的3D 流體細(xì)節(jié)效果。

2 相關(guān)工作

流體模擬一直以來(lái)都是大規(guī)模圖形模擬的重點(diǎn)研究方向，到目前為止，已經(jīng)提出了許多流體模擬方法，并取得比較真實(shí)的結(jié)果。

SPH 是一種流體模擬算法，相對(duì)于其它的模擬算法，其特點(diǎn)就是簡(jiǎn)單快速，是典型的拉格朗日視角，它的基本原理就是通過(guò)粒子模擬流體的運(yùn)動(dòng)規(guī)律，可以有效地模擬小尺度的流體運(yùn)動(dòng)。

現(xiàn)階段在SPH 上的主要研究工作是實(shí)現(xiàn)實(shí)時(shí)SPH 流體模擬，主要采用并行加速計(jì)算方法。SPH在模擬各種流體效果的應(yīng)用方面也取得了很多成果，如水體［18］、多相流體［19］和固液耦合［20］。Amada，et al［21］首次提出了基于圖形處理器的SPH 流體模擬，該方法在中央處理器上執(zhí)行最近鄰域搜索（Nearest Neighbor Search），并在圖形處理器上求解N-S方程，以加速流體模擬。之后，Harada，et al［22］、Zhang，et al［23］將除渲染場(chǎng)景之外的所有步驟都放在GPU上，進(jìn)一步提高了流體模擬的效率。

現(xiàn)實(shí)生活中的很多流體邊界層現(xiàn)象不能用單相流體來(lái)模擬，多流體模擬是解決這一現(xiàn)象的關(guān)鍵。Müller，et al［24］首先引入體積分?jǐn)?shù)的概念來(lái)表示不相同的空間分布，實(shí)現(xiàn)混合流體的多流體模擬系統(tǒng)。Kim［25］、Misztal，et al［26］采用網(wǎng)格的方法，并假設(shè)混合流體不可壓縮且不發(fā)散情況下，該方法在拓?fù)滢D(zhuǎn)換上獲得較好的結(jié)果。面對(duì)混合現(xiàn)象由濃度差引起的擴(kuò)散效應(yīng)，Liu，et al［27］引入了一種新的多流體模擬方法，該方法能夠真實(shí)地捕捉由于多種流體之間的相對(duì)運(yùn)動(dòng)、湍流相互作用和力變化分布而導(dǎo)致的復(fù)雜混合現(xiàn)象。湍流相互作用和力變化分布而導(dǎo)致的復(fù)雜混合現(xiàn)象。

目前，主流方法是基于SPH的表面張力來(lái)優(yōu)化流體的邊界層模擬效果。早期工作大多集中在連續(xù)體表面力CSF（Continuum Surface Force）方法上。CSF 的主要缺點(diǎn)是在自由表面情況下嚴(yán)重依賴粒子樣本。如果沒(méi)有足夠的粒子，會(huì)導(dǎo)致模擬結(jié)果數(shù)值不穩(wěn)定性和不準(zhǔn)確性的問(wèn)題。此外，He等［28］提出在微水平上將表面張力模擬為粒子之間的分子力方法，實(shí)現(xiàn)流體邊界模擬。最近，何小偉等［29］提出一種新的表面張力模型，該模型是基于擴(kuò)散界面方法處理稀疏采樣的自由表面。與其它表面張力相比，毛細(xì)管波尤其難以模擬，因?yàn)樵摲椒ㄐ枰恢潞蜏?zhǔn)確的表面張力估計(jì)。

綜上，本文提出基于SPH 算法，并結(jié)合壓縮波進(jìn)行邊界層細(xì)節(jié)優(yōu)化的方法實(shí)現(xiàn)多流體細(xì)節(jié)模擬。

3 算法實(shí)現(xiàn)

論文算法方案結(jié)構(gòu)如圖1 所示。算法流程如下。

圖1 整體算法結(jié)構(gòu)

第一步：確定所要模擬的流體屬性，選擇本文優(yōu)化的密度光滑核函數(shù)對(duì)應(yīng)的屬性因子m；

第二步：光滑核函數(shù)代入到SPH 算法，求得粒子的密度；

第三步：是否選擇優(yōu)化流體效果，如果選擇Y進(jìn)入到第四步；否則直接進(jìn)入到第六步；

第四步：通過(guò)流體的邊界層進(jìn)行邊界層粒子劃分；第五步：通過(guò)聲波公式求出代更新的邊界層密度；第六步：進(jìn)行粒子的初始化求粒子的物理屬性（例如重力、粘滯力等），得出粒子的加速度；

第七步：通過(guò)MC 算法進(jìn)行流體表面構(gòu)建得到想要的流體模擬效果。

3.1 SPH等式

SPH 是最流行的拉格朗日方法。它能有效地模擬小尺度自由表面流體。不可壓縮牛頓流體的運(yùn)動(dòng)由N-S 方程描述，而SPH 方法中選取簡(jiǎn)易的N-S 方程，分析得出作用在一個(gè)微粒上的作用力由三部分組成，如式（1）所示：

其中稱為外部力，一般就是重力，如式（2）所示：

是由流體內(nèi)部的壓力差產(chǎn)生的作用力，如式（3）所示：

是由粒子間的速度差引起的，跟流體的黏度系數(shù)以及速度差有關(guān)，如式（4）所示：

根據(jù)式（2）~式（4）帶入到式（1），并帶入公式，可以得到：

加速度則為

以上就是粒子運(yùn)動(dòng)學(xué)的計(jì)算方法，其中u 表示速度，ρ表示密度，p 表示壓力，g 為重力表示外力，μ表示粘度系數(shù)。

SPH 算法與其他的流體力學(xué)中的數(shù)學(xué)方法類似，同樣涉及到“光滑核”的概念，因?yàn)榱黧w中每個(gè)位置參與運(yùn)算的值都是由周圍一組粒子累加起來(lái)的，則該處某項(xiàng)屬性Attribute 的累加公式，如式（7）所示：

其中Attributej是要累加的某種屬性，mj和ρj是周圍粒子j的質(zhì)量和密度，和是粒子i 和j的位置，h是光滑核半徑，而W就是光滑函數(shù)。

求解N-S 方程需要滿足動(dòng)量守恒和質(zhì)量守恒。由于SPH粒子的數(shù)量在產(chǎn)生后不會(huì)發(fā)生變化，所以SPH中的粒子i在其運(yùn)動(dòng)過(guò)程中可以滿足質(zhì)量守恒。SPH粒子i的加速度可以由式（5）計(jì)算得出：

其中Fi是等式（1）外力之和。ai為本文更新粒子i的位置的加速度，推算得到的加速度公式如式（6）所示。

3.2 基于SPH的優(yōu)化

由于基于SPH 的原始方法以及擴(kuò)展方法都無(wú)法實(shí)現(xiàn)一個(gè)統(tǒng)一的流體控制框架，針對(duì)這個(gè)問(wèn)題，本文對(duì)SPH的光滑核函數(shù)進(jìn)行優(yōu)化，計(jì)算粒子密度的光滑核函數(shù)可以改寫(xiě)成式（9）所示：

接下來(lái)根據(jù)式（7），用ρ代替Attribute 可以得到式（10）：

計(jì)算密度中將不再使用之前光滑核函數(shù)稱為Poly6 的光滑核函數(shù)，而是使用式（9）。現(xiàn)將其帶入到式（10）中，則處的密度計(jì)算公式最終為式（11）：

根據(jù)式（7），壓力計(jì)算中使用的光滑核函數(shù)稱為Spiky 函數(shù)，因此在ri處由壓力產(chǎn)生的加速度部分為式（12）所示：

根據(jù)式（7），粘滯力計(jì)算中使用的光滑核函數(shù)稱為Wviscosity，在ri處由粘滯力產(chǎn)生的加速度部分為式（13）所示：

將式（12）、式（13）代入到公式中得到式（14）所示：

通過(guò)修改本文采用的密度光滑核函數(shù)參數(shù)，對(duì)基本的流體仿真模型進(jìn)行復(fù)雜形變，包括流體（火、煙）抖動(dòng)形變和流體濺落形變。對(duì)基本流體仿真模型進(jìn)行的復(fù)雜形變通過(guò)修改密度光滑核函數(shù)中的屬性因子m、分母系數(shù)k 實(shí)現(xiàn)，在可控區(qū)域內(nèi)，m=1、0.5、h（核半徑）是最接近現(xiàn)實(shí)流體的仿真構(gòu)建，粒子顯示結(jié)果如圖2所示。

圖2 統(tǒng)一框架呈現(xiàn)的粒子效果雛形

3.3 基于聲波的細(xì)節(jié)處理

本文在流體仿真細(xì)節(jié)優(yōu)化處理中，使用聲波的波動(dòng)方程來(lái)更新密度變化，如式（15）所示：

其中c為波速，ρouter表示流體表面邊界層中的粒子關(guān)于位置x和時(shí)間t的波強(qiáng)度的一個(gè)測(cè)量。

式（16）所示，其中γ是表面張力，k 是波數(shù)，ρ0的流體密度。對(duì)于水，論文使用ρ0=1.0*103kg/m3，對(duì)于氣體相對(duì)應(yīng)的可以直接使用空氣的密度ρ0=1.293 kg/m3。

盡管聲波和水體的毛細(xì)管波雖然很多的行為相似，但是還有一定的區(qū)別：液體毛細(xì)管波只在表面?zhèn)鞑ィ暡梢栽谝后w的體積內(nèi)傳播。如果只是簡(jiǎn)單地求解式（16），模擬流體會(huì)得到不一樣的模擬效果，尤其在模擬流體的粒子數(shù)目過(guò)多的情況下。為了實(shí)現(xiàn)更加逼真效果，針對(duì)上述問(wèn)題，論文采用的解決方案是檢測(cè)流體表面邊界層的粒子，并對(duì)粒子進(jìn)行范圍的劃分，讓聲波在不同邊界層范圍內(nèi)進(jìn)行傳播，該方法可以實(shí)現(xiàn)較好的流體細(xì)節(jié)效果。因此，本文引入了流體的邊界層劃分算法，如式（17）所示：

其中式（17）表示布拉修斯對(duì)平板層流邊界層求得的精確解，其中v 表示運(yùn)動(dòng)粘度，v∞表示的是主體流速，δ表示的是邊界層厚度，邊界層厚度指的就是壁面到流速與主體流速v∞相差不到1%的流體層厚度，如式（18）所示：

其中μ表示流體的動(dòng)力粘度（粘性的度量，也稱粘性系數(shù)），ρ表示液體密度。

根據(jù)式（14），可以求得每個(gè)粒子的加速度情況，因?yàn)樗俣仁鞘噶浚敲粗髁黧w的速度，如式（19）所示：

其中公式t表示時(shí)間，n表示粒子數(shù)。

如圖3 所示，接下來(lái)就是將式（19）中得到的主流體流速代入到式（17）中，得到了流體的邊界層，將式（16）帶入到式（15）可以得到ρouter，更新邊界層中粒子的密度為ρouter。就可以將聲波形成的細(xì)節(jié)效果帶入到流體的表面。以水流為例，實(shí)現(xiàn)效果與沒(méi)有聲波豐富細(xì)節(jié)明顯的區(qū)別如圖4 所示，對(duì)比效果，圖（a）細(xì)節(jié)化程度明顯弱于圖（b）采用的本文方法。且圖（a）中文獻(xiàn)［30］的實(shí)現(xiàn)方法需要構(gòu)造一個(gè)密集的點(diǎn)集的拓?fù)渥兓砑恿怂惴ǖ膹?fù)雜度，因此本文中的方法更加容易實(shí)現(xiàn)。

圖3 邊界層粒子劃分

圖4 （a）（b）表面波紋對(duì)比

為了進(jìn)一步驗(yàn)證本文細(xì)節(jié)優(yōu)化處理的有效性，通過(guò)粒子實(shí)驗(yàn)來(lái)進(jìn)行驗(yàn)證，計(jì)算參數(shù)有：流體初始密度為ρ=1000kg/m3，流體模擬粒子數(shù)為2200 個(gè)，邊界粒子數(shù)為120 個(gè)，粒子影響半徑r=0.02m。從圖5 中可以看出，在其他條件相同的情況下，圖（a）相對(duì)圖（b）的流體表面均勻，圖（b）的凹凸效果更加明顯，因此本文方法可以提高流體表面細(xì)節(jié)仿真效果。

圖5 粒子方法對(duì)比

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

本文實(shí)驗(yàn)是在一臺(tái)配備英特爾酷睿i7-8700K 3.70GHz、英偉達(dá)GTX1080 顯卡、16G 內(nèi)存，操作系統(tǒng)為Windows 10 的電腦上測(cè)試文中的算法。整個(gè)實(shí)現(xiàn)基于OpenGL，使用的是MC 算法（Marching Cube Algorithm）進(jìn)行流體液面繪制。這里給出了一些實(shí)驗(yàn)結(jié)果，這些結(jié)果顯示了本文所提出的方法的有效性。經(jīng)過(guò)實(shí)驗(yàn)不僅完成了整個(gè)基于物理的流體方程的推導(dǎo)，同時(shí)也將計(jì)算的結(jié)果實(shí)時(shí)繪制出來(lái)。

實(shí)驗(yàn)結(jié)果驗(yàn)證了本文算法所對(duì)產(chǎn)生豐富的細(xì)節(jié)多樣化流體效果，論文中的實(shí)驗(yàn)例子包括多樣化的流體雛形效果（如圖2），圖2 的實(shí)驗(yàn)得到了一個(gè)基于SPH 不局限單一流體的框架。在圖2 多樣化的流體雛形的框架基礎(chǔ)上，添加粒子數(shù)，為粒子添加對(duì)應(yīng)流體的顏色，結(jié)合MC 的流體表面繪制，呈現(xiàn)小型流體波紋效果（如圖3），火焰燃燒（圖6）的流體效果，以及煙霧吹動(dòng)的效果（圖7）。這些例子說(shuō)明了本文算法的正確性，以及類似現(xiàn)實(shí)世界中的流體的波動(dòng)效果。在圖5 中，文中方法（圖（b）、（c））比較Müller's（圖（a））的方法在火焰上實(shí)現(xiàn)的對(duì)比，可以看出文中的方法優(yōu)化流體細(xì)節(jié)方向這一塊，使得流體顯示的效果更加鱗次櫛比。經(jīng)過(guò)圖7的實(shí)驗(yàn)，圖（a）是將何曉偉等［30］實(shí)現(xiàn)毛細(xì)管波的方法來(lái)豐富煙霧的細(xì)節(jié)化，效果是近似現(xiàn)實(shí)效果，但計(jì)算時(shí)間長(zhǎng)，圖（b）、（c）實(shí)驗(yàn)展示本文的方法應(yīng)用到煙霧模擬中的，圖（b）中粒子細(xì)節(jié)效果更加豐富，再進(jìn)行表面繪制得到的圖（c）的結(jié)果，經(jīng)過(guò)實(shí)現(xiàn)得到近似現(xiàn)實(shí)世界煙霧效果，在煙霧剛開(kāi)始的起始階段基于文中算法得到的細(xì)節(jié)效果好于文獻(xiàn)［30］的方法。

圖6 豐富的細(xì)節(jié)化火焰

圖7 豐富的細(xì)節(jié)化煙霧

由表1~2 中的數(shù)據(jù)結(jié)合圖5 和圖7 的模擬效果可以得出在相同的條件下，隨著時(shí)間的增加，本文的最高速度與最低速度差明顯大于Müller's的方法和文獻(xiàn)［30］的方法，進(jìn)一步得出本文方法速度波動(dòng)大，流體模擬細(xì)節(jié)越顯著。

表1 N=45000數(shù)據(jù)

表2 N=70000數(shù)據(jù)

本文采用的算法實(shí)現(xiàn)多樣化的流體效果，是局限小型流體，大型流體的效果會(huì)出現(xiàn)過(guò)度偏差細(xì)節(jié)化優(yōu)化也會(huì)出現(xiàn)失真的效果，并且流體運(yùn)動(dòng)過(guò)長(zhǎng)，細(xì)節(jié)方面也會(huì)出現(xiàn)精度問(wèn)題。算法在細(xì)節(jié)優(yōu)化步驟上是基于壓縮波的，相對(duì)于液體的毛細(xì)管波等現(xiàn)實(shí)波紋細(xì)節(jié)方面，并不意味著物理上的精確。文中算法不指定流體的具體密度，是根據(jù)單一的密度來(lái)實(shí)現(xiàn)多樣化的流體細(xì)節(jié)優(yōu)化的效果。最后，基于文中算法的基礎(chǔ)得到的結(jié)果在某種程度上是取決于表面的粒子顆粒上的。

5 結(jié)語(yǔ)

論文提出了一種框架，基于SPH實(shí)現(xiàn)豐富細(xì)節(jié)的多樣化的流體模擬。基于SPH 的基礎(chǔ)下優(yōu)化光滑核函數(shù)實(shí)現(xiàn)了多樣化的流體框架不再局限于單一的流體模擬，并且結(jié)合聲波實(shí)現(xiàn)細(xì)節(jié)優(yōu)化的方法。經(jīng)過(guò)實(shí)驗(yàn)，表明了這是一個(gè)合理的假設(shè)。經(jīng)過(guò)實(shí)驗(yàn)研究進(jìn)一步表明了，在SPH框架聲波更新的表面粒子可以為優(yōu)化流體的細(xì)節(jié)提供很好的方向。

本文中的方法可以實(shí)現(xiàn)多樣化的流體模擬，但是將該方法應(yīng)用于大型流體模擬是我進(jìn)一步研究的工作。另外，在未來(lái)，希望基于以上算法的基礎(chǔ)，進(jìn)行理論的研究，實(shí)現(xiàn)多種多樣流體控制。并且想進(jìn)一步研究氣液耦合方面的壓強(qiáng)問(wèn)題，以獲得更精確的細(xì)節(jié)優(yōu)化效果。更進(jìn)一步的話，將通過(guò)機(jī)器學(xué)習(xí)來(lái)提高文中算法的性能，提高整體算法的模擬速度。