邵緒強(qiáng)， 劉 艷， 王新穎

(1. 華北電力大學(xué)控制與計(jì)算機(jī)工程學(xué)院，河北 保定 071003；2. 河北金融學(xué)院國(guó)際教育學(xué)院，河北 保定 071003)

一種拉格朗日粒子流體的高效表面重建方法

邵緒強(qiáng)1， 劉 艷2， 王新穎1

(1. 華北電力大學(xué)控制與計(jì)算機(jī)工程學(xué)院，河北 保定 071003；2. 河北金融學(xué)院國(guó)際教育學(xué)院，河北 保定 071003)

為了在幾何表示層面捕獲拉格朗日粒子流體的表面細(xì)節(jié)特征以進(jìn)行真實(shí)感繪制，提出一種高效的基于八叉樹的自適應(yīng)流體表面重建方法。首先進(jìn)行流體表面粒子的精確檢測(cè)；然后根據(jù)結(jié)點(diǎn)中是否包含流體表面粒子，對(duì)流體粒子的包圍盒進(jìn)行基于八叉樹的自適應(yīng)剖分；最后只在包含流體表面粒子的葉子結(jié)點(diǎn)里建立隱式距離場(chǎng)。該方法在流體表面重建過(guò)程中的內(nèi)存消耗和計(jì)算復(fù)雜度只取決于流體表面而不是流體體積，適合大規(guī)模粒子流體場(chǎng)景的表面提取和繪制。

拉格朗日粒子；表面重建；流體模擬；真實(shí)感繪制

波濤洶涌的大海、傾盆而下的暴雨、瞬間就要吞噬戰(zhàn)艦的巨大漩渦等流體影視特效給好萊塢特效大片的觀眾帶來(lái)了震撼的視覺享受。在欣賞這些精彩特效大片的同時(shí)，人們領(lǐng)略到了基于物理的流體模擬技術(shù)的高超魅力。由于真實(shí)物理模型的支撐，基于物理的流體模擬技術(shù)能逼真地模擬各種復(fù)雜流體現(xiàn)象，給人以震撼的視覺享受，因此吸引了國(guó)內(nèi)外圖形學(xué)研究學(xué)者們對(duì)其進(jìn)行深入研究。

基于物理的流體模擬方法[1]主要分為歐拉網(wǎng)格方法(Eulerian grid method)和拉格朗日粒子方法(Lagrangian particle method)。歐拉網(wǎng)格方法把流體占據(jù)的計(jì)算空間按一定規(guī)則剖分為網(wǎng)格，然后分析被流體所充滿的空間中每一個(gè)網(wǎng)格位置上流體的速度、壓強(qiáng)、密度等參數(shù)隨時(shí)間的變化。該方法可以較容易地保證流體的不可壓縮特性，但由于模擬過(guò)程中要同時(shí)把固體邊界網(wǎng)格離散化，不易于進(jìn)行流固交互模擬，并且會(huì)丟失精度小于網(wǎng)

格分辨率的流體細(xì)節(jié)特征。而拉格朗日粒子方法從分析流體各個(gè)微團(tuán)的運(yùn)動(dòng)著手，即研究每個(gè)流體微團(tuán)的速度、壓強(qiáng)、密度等描述流體運(yùn)動(dòng)的參數(shù)隨時(shí)間的變化狀況，這是一種以“實(shí)體”形式為流體建模的方法，建模的基本元素就是攜帶流場(chǎng)變量的粒子。與歐拉網(wǎng)格方法相比，拉格朗日粒子方法更直觀，自動(dòng)保證數(shù)值求解過(guò)程中的質(zhì)量守恒，尤其易于進(jìn)行固體邊界處理和模擬流體細(xì)節(jié)。因此，在電子游戲和電影特效等實(shí)際流體模擬應(yīng)用中，拉格朗日粒子方法變得越來(lái)越普遍。

為了得到逼真的流體動(dòng)畫，在通過(guò)數(shù)值求解流體運(yùn)動(dòng)方程而用粒子描述流體運(yùn)動(dòng)信息后，拉格朗日粒子方法需要重建出光滑完整的隱式流體表面以進(jìn)行真實(shí)感繪制。拉格朗日粒子流體的隱式表面重建方法可分為密度場(chǎng)方法和距離場(chǎng)方法。對(duì)于密度場(chǎng)方法，Blinn[2]直接用一個(gè)連續(xù)性表面包裹所有流體粒子而形成密度云，并且這個(gè)表面隨流體粒子的運(yùn)動(dòng)而改變。Müller等[3]把該方法引入到光滑粒子流體動(dòng)力學(xué)(smoothed particle hydrodynamics，SPH)流體中，并利用球形核函數(shù)重新定義密度場(chǎng)，使得表面重建結(jié)果更穩(wěn)定，但存在局部區(qū)域出現(xiàn)鼓包的現(xiàn)象，從而不能重建出平滑的流體表面。為了使密度場(chǎng)方法能夠重建出足夠光滑的流體表面，Yu和Turk[4]基于粒子分布的主成分分析而用橢圓形核函數(shù)來(lái)插值計(jì)算密度場(chǎng)；為了解決密度場(chǎng)方法受限于粒子采樣的均勻性，Zhu和 Bridson[5]提出了基于距離場(chǎng)的表面重建方法，通過(guò)直接計(jì)算流體粒子到表面的距離來(lái)表示隱式流體表面，但該方法在高曲率區(qū)域存在明顯的瑕疵。Solenthaler等[6]通過(guò)動(dòng)態(tài)調(diào)整粒子的半徑使質(zhì)心移動(dòng)速度小于粒子半徑衰減速度，從而改善了距離場(chǎng)方法在高曲率區(qū)域存在的瑕疵現(xiàn)象。然而，上述流體表面重建方法都是針對(duì)于流體體積進(jìn)行的，當(dāng)對(duì)大規(guī)模拉格朗日粒子流體場(chǎng)景進(jìn)行表面重建時(shí)就會(huì)由于過(guò)高的內(nèi)存和時(shí)間消耗而失敗。針對(duì)此問(wèn)題，本文給出了一種高效的基于八叉樹的自適應(yīng)流體表面重建方法。該方法首先進(jìn)行流體表面粒子的精確檢測(cè)；然后根據(jù)結(jié)點(diǎn)中是否包含流體表面粒子，對(duì)流體粒子的包圍盒進(jìn)行基于八叉樹的自適應(yīng)剖分；最后只在包含流體表面粒子的葉子結(jié)點(diǎn)中建立隱式距離場(chǎng)。因此，流體表面重建過(guò)程中的內(nèi)存消耗和計(jì)算復(fù)雜度只取決于流體表面而不是流體體積，非常適合于大規(guī)模粒子流體場(chǎng)景的表面提取和繪制。

1 拉格朗日粒子流體模擬

1.1 不可壓縮粒子流體模擬

確保不可壓縮性的計(jì)算是SPH流體模擬最困難的部分，也是其關(guān)鍵研究的內(nèi)容。微可壓縮SPH流體[7](weakly compressible SPH，WCSPH)必須使用較小的時(shí)間步長(zhǎng)以確保不可壓縮性，嚴(yán)重影響計(jì)算效率。通過(guò)泊松方程計(jì)算壓強(qiáng)的不可壓縮方法[8-9](incompressible SPH，ISPH)可以使用較大時(shí)間步長(zhǎng)，但在每個(gè)時(shí)間步內(nèi)由于迭代求解線性方程組而導(dǎo)致計(jì)算量龐大。Ihmsen等[10]結(jié)合對(duì)稱的SPH壓力和連續(xù)方程離散方法來(lái)離散化壓強(qiáng)泊松方程，并對(duì)壓力進(jìn)行實(shí)際計(jì)算，提高了求解的收斂速度。

為了能夠在模擬過(guò)程中使用較大時(shí)間步長(zhǎng)，本文采用Solenthaler和Pajarola[11]提出的預(yù)測(cè)修正不 可 壓 縮 SPH 方 法 (predictive-corrective incompressible SPH，PCISPH)進(jìn)行不可壓縮流體模擬。PCISPH結(jié)合了WCSPH和ISPH的優(yōu)點(diǎn)，允許使用較大時(shí)間步長(zhǎng)，并簡(jiǎn)化了每個(gè)時(shí)間步內(nèi)的計(jì)算復(fù)雜度。算法的核心思路是通過(guò)預(yù)測(cè)修正密度誤差的方法迭代計(jì)算流體粒子的壓強(qiáng)和壓力，壓力的改變會(huì)影響粒子的速度和位置，然后位置的相對(duì)變化會(huì)影響粒子密度，最后粒子的相對(duì)位置和密度影響壓力大小和方向，在這樣的迭代作用的循環(huán)中找到滿足約束條件的壓強(qiáng)為止。PCISPH的偽代碼如圖1所示。

圖1 PICSPH偽代碼

圖1給出了PCISPH在每個(gè)時(shí)間步內(nèi)的迭代計(jì)算過(guò)程。在迭代計(jì)算過(guò)程中，首先對(duì)每個(gè)粒子搜索其核函數(shù)緊支域內(nèi)的鄰居粒子信息。然后計(jì)算出除壓力外的所有受力，包括粘性力、重力、表面張力、甚至漩渦力。進(jìn)入預(yù)測(cè)修正循環(huán)之前，粒子的壓強(qiáng)和壓力初始化為 0。最后便進(jìn)入PCISPH的壓強(qiáng)計(jì)算的循環(huán)中。每次循環(huán)都要檢查粒子的密度誤差的大小是否超過(guò)閾值η，其中密度誤差為。如果，則使用當(dāng)前粒子的合力來(lái)預(yù)測(cè)速度和位置，根據(jù)預(yù)測(cè)出來(lái)的粒子速度和位置重新計(jì)算粒子的密度和密度誤差，并將所有粒子的最大誤差保存在密度誤差中。在計(jì)算過(guò)程中，根據(jù)密度誤差對(duì)流體的壓強(qiáng)值進(jìn)行累加

在預(yù)測(cè)下一步的密度時(shí)，應(yīng)該根據(jù)預(yù)測(cè)的位置信息更新鄰居粒子列表?？紤]到計(jì)算效率，這里不再重新搜索鄰居粒子，而是通過(guò)更新到現(xiàn)有鄰居粒子距離的方式來(lái)近似。雖然這樣的估計(jì)會(huì)帶來(lái)誤差，但是在最少迭代次數(shù)、時(shí)間步長(zhǎng)、核函數(shù)影響范圍衰減等約束條件的作用下，對(duì)整個(gè)流體運(yùn)動(dòng)模擬結(jié)果影響是可以忽略的。采用一個(gè)最小迭代次數(shù)的方法可顯著減少誤差的局部震蕩，而且這給予粒子足夠的時(shí)間來(lái)擴(kuò)散預(yù)測(cè)信息。為了提高計(jì)算效率，一般使用最少迭代3次，而3次迭代已經(jīng)足夠?qū)⒘Ｗ用芏日`差約束到目標(biāo)范圍內(nèi)。

從PCISPH的執(zhí)行過(guò)程來(lái)看，鄰居搜索是其最大的計(jì)算瓶頸，為了提高該部分的計(jì)算性能，可建立 KD-Tree以實(shí)現(xiàn)鄰居例子的快速搜索。與WCSPH相比，PCISPH由于增加了迭代計(jì)算來(lái)計(jì)算壓力而大幅度降低了性能，但是采用的時(shí)間步長(zhǎng)卻增加了30倍左右。

1.2 流固交互模擬

對(duì)于SPH粒子流體和運(yùn)動(dòng)固體邊界的交互，本文結(jié)合文獻(xiàn)[12]的懲罰力交互方法和文獻(xiàn)[13]的自適應(yīng)邊界粒子采樣方法實(shí)現(xiàn)雙向流固交互模擬，計(jì)算過(guò)程分為3步：

(1) 在每個(gè)時(shí)間步內(nèi)進(jìn)行流固交互計(jì)算之前，根據(jù)初始流體粒子間距r0對(duì)任一表面三角形進(jìn)行自適應(yīng)邊界粒子采樣，如圖 2所示。首先選取三角形頂點(diǎn)為邊界粒子；然后以初始流體粒子間距r0對(duì)三角形的邊進(jìn)行邊界粒子采樣；確定三角形最短邊es和最長(zhǎng)邊el，并選取向內(nèi)的法線方向?yàn)閽呙杈€方向?qū)θ切蝺?nèi)部進(jìn)行邊界粒子采樣。自適應(yīng)邊界粒子采樣可以更精確地表示固體表面的拓?fù)浣Y(jié)構(gòu)，從而實(shí)現(xiàn)更準(zhǔn)確的固流交互。

圖2 自適應(yīng)收縮壓粒子采樣

(2) 考慮邊界粒子的相對(duì)密度貢獻(xiàn)[12-13]，計(jì)算任一流體粒子i的密度

其中，j為鄰居流體粒子索引；k為屬于固體邊界的鄰居粒子的索引；W為 SPH方法的核函數(shù)；為相對(duì)貢獻(xiàn)度量函數(shù)； Vbk為第k個(gè)邊界鄰居粒子當(dāng)前體積。式(2)通過(guò)φ() x度量邊界粒子的相對(duì)貢獻(xiàn)，解決了其非均勻分布特性。

(3) 對(duì)于一對(duì)相鄰的流體粒子i和固體粒子j，計(jì)算 j對(duì) i的包括法向力和切向力在內(nèi)的交互力。利用 SPH方法求解流體方程得到的壓力、粘滯力來(lái)分別近似法向力、切向力。并且，對(duì)i的作用力同樣采用φ() x進(jìn)行度量

其中，μ為物體表面磨擦系數(shù)；v為速度。

2 粒子流體表面重建

2.1 流體表面粒子檢測(cè)

流體的表面可以由位于流體表面的粒子表示和跟蹤?；谟?jì)算流體動(dòng)力學(xué)(computational fluid dynamics，CFD)領(lǐng)域文獻(xiàn)[14]方法，給出一個(gè)快速準(zhǔn)確的流體表面粒子檢測(cè)方法。對(duì)于任一流體粒子i，首先計(jì)算一個(gè)再歸一化矩陣Bi

其中，Vj是第 j個(gè)鄰居粒子的體積；W為再歸一化的高斯核函數(shù)[15]； 40h= d，d0為初始流體粒子間距。

然后，計(jì)算矩陣Bi的最小特征值，并根據(jù)Bi和實(shí)驗(yàn)驗(yàn)證值的大小關(guān)系判定粒子 i是否屬于SPH流體表面S

由式(6)可得到一個(gè)包含了SPH流體表面粒子的粒子集合。與文獻(xiàn)[14]方法相比，本文方法通過(guò)設(shè)定miniλ 的閾值為0.75，把自由流體表面附近的流體粒子都確定為流體表面粒子集合，雖然引入了一些流體內(nèi)部的粒子，但提高了檢測(cè)準(zhǔn)確度。

2.2 自適應(yīng)距離場(chǎng)重建

為了降低粒子流體表面重建過(guò)程中的內(nèi)存消耗和計(jì)算復(fù)雜度，基于文獻(xiàn)[16]中流體表面粒子檢測(cè)和隱式表面定義方法，給出一種基于八叉樹[17]的自適應(yīng)距離場(chǎng)構(gòu)建算法，該方法只在包含流體表面粒子的葉子結(jié)點(diǎn)上計(jì)算距離場(chǎng)值。

在自適應(yīng)距離場(chǎng)構(gòu)建算法中，任一八叉樹結(jié)點(diǎn)包含兩個(gè)粒子集合：流體表面粒子集合和流體粒子集合。其中，流體表面粒子集合包含所有屬于該結(jié)點(diǎn)的流體表面粒子，流體粒子集合包含所有屬于該結(jié)點(diǎn)的流體粒子。本方法把八叉樹的結(jié)點(diǎn)分為 3類：外部結(jié)點(diǎn)、內(nèi)部結(jié)點(diǎn)和表面結(jié)點(diǎn)。對(duì)于外部結(jié)點(diǎn)來(lái)說(shuō)，其粒子集合流體表面粒子集合和流體粒子集合均為空；內(nèi)部結(jié)點(diǎn)的粒子集合流體表面粒子集合為空，而粒子集合流體粒子集合非空；而對(duì)于表面結(jié)點(diǎn)，其粒子集合流體表面粒子集合非空。八叉樹距離場(chǎng)的遞歸創(chuàng)建過(guò)程如圖3所示。

圖3 基于八叉樹的自適應(yīng)距離場(chǎng)的建立

從圖 3中左圖可以看出，基于八叉樹的自適應(yīng)距離場(chǎng)的創(chuàng)建過(guò)程分為 3步：①先為整個(gè)流體模擬建立一個(gè)根結(jié)點(diǎn)。該根結(jié)點(diǎn)的物理空間為整個(gè)流體的最小立方體包圍盒，流體表面粒子集合包含所有的流體表面粒子，流體粒子集合包含所有流體粒子。②從根結(jié)點(diǎn)開始，對(duì)每個(gè)八叉樹結(jié)點(diǎn)進(jìn)行遞歸剖分，直到其粒子集合流體表面粒子集合為空或者達(dá)到最大深度。③對(duì)于每個(gè)葉子結(jié)點(diǎn)，如果其流體表面粒子集合流體表面粒子集合非空，則根據(jù)下式計(jì)算其8個(gè)頂點(diǎn)的距離場(chǎng)值

半徑； R= 4r； f∈ [0… 1 ]并通過(guò)下式進(jìn)行計(jì)算

其中，EVmax為矩陣的最大特征值；γ通過(guò)文獻(xiàn)[17]方法進(jìn)行計(jì)算。圖3右圖給出了如何確定八叉樹結(jié)點(diǎn)包含的粒子集合以避免粒子流體表面重建過(guò)程中產(chǎn)生的空洞。當(dāng)為八叉樹的任一結(jié)點(diǎn)建立粒子集合流體表面粒子集合和流體粒子集合時(shí)，本方法考慮落在比該結(jié)點(diǎn)物理空間大l的AABB包圍盒內(nèi)的所有流體表面粒子和流體粒子。在本方法中，l等于式(7)中的R。

為了保證基于八叉樹的流體表面重建算法能夠重建出僅由一層粒子構(gòu)成的薄膜結(jié)構(gòu)的表面，可給出了八叉樹葉子結(jié)點(diǎn)邊長(zhǎng)和粒子半徑的大小關(guān)系以避免表面重建失敗，如圖 4所示。在圖 4中，用二維示意圖的形式分析了八叉樹葉子結(jié)點(diǎn)邊長(zhǎng)L和粒子半徑r的大小關(guān)系，其中r值對(duì)于確定流體場(chǎng)景是固定的。圖 4(a)給出的特例表明，

圖4 粒子半徑r與葉子結(jié)點(diǎn)邊長(zhǎng)L的大小關(guān)系

3 實(shí)例驗(yàn)證

將通過(guò)一系列 3D實(shí)驗(yàn)證明本文提出的針對(duì)拉格朗日粒子流體的表面重建算法。實(shí)驗(yàn)硬件配置為：Intel Xeon E5630 CPU，6 GB RAM 和NVIDIA Geforce GTX 680 GPU。編程環(huán)境為：Windows 7，Visual Studio 2008，OpenGL 2.0。在3D場(chǎng)景中，首先利用SPH粒子方法模擬了流體運(yùn)動(dòng)，然后利用本文給出的自適應(yīng)流體表面重建方法創(chuàng)建流體的隱式表面，最后利用Marching Cubes算法[3]提取三角形網(wǎng)格面并利用開源渲染引擎Pov-Ray進(jìn)行離線繪制。實(shí)驗(yàn)中用到的物理量以及取值如表1所示。

表1 漩渦模擬中物理參數(shù)及取值

3.1 表面重建效果分析

基于八叉樹的自適應(yīng)流體表面重建算法可以高效地重建出任一粒子集合的隱式表面，并捕獲漩渦細(xì)節(jié)的幾何結(jié)構(gòu)以進(jìn)行真實(shí)感繪制。該算法獲取流體表面的過(guò)程如圖5所示。在圖5中利用基于八叉樹的自適應(yīng)流體表面重建算法對(duì)犰狳的采樣粒子集合進(jìn)行了表面重建。對(duì)于輸入的犰狳網(wǎng)格模型，首先對(duì)其所占的空間進(jìn)行規(guī)則粒子采樣得到12.1 k粒子(圖5(a))，采樣距離為0.018 m。接著利用本文給出的表面粒子檢測(cè)方法檢測(cè)出粒子集合的表面粒子；圖5(b)給出了利用本方法進(jìn)行表面粒子檢測(cè)得到的結(jié)果，表面粒子數(shù)為3.5 k。然后，根據(jù)檢測(cè)到的流體表面粒子，建立八叉樹自適應(yīng)距離場(chǎng)(圖 5(c))。最后圖 5(d)給出了Marching Cubes算法提取的三角形網(wǎng)格面的繪制結(jié)果。

圖5 自適應(yīng)流體表面重建算法

圖 6給出了初始形狀為犰狳模型的流體在重力作用下的運(yùn)動(dòng)模擬結(jié)果。通過(guò)使用本文提出的粒子流體表面重建方法建立距離場(chǎng)，并利用開源繪制引擎Pov-Ray對(duì)Marching Cubes算法提取的三角形網(wǎng)格面進(jìn)行真實(shí)感繪制。從繪制結(jié)果可以看出，本文的方法可以光滑地重建出具有復(fù)雜形狀的粒子流體表面，清晰地表現(xiàn)飛濺液滴等細(xì)小流體特征。

圖6 犰狳模型流體的表面重建

圖7給出了剛性圓柱體以10 m/s的速度在水面移動(dòng)產(chǎn)生的漩渦細(xì)節(jié)模擬結(jié)果。從模擬結(jié)果可以看出，本文的表面重建算法可以重建出剛性圓柱體后面產(chǎn)生的豐富漩渦細(xì)節(jié)特征，增強(qiáng)了流固交互模擬的真實(shí)感。

圖 8給出了多個(gè)具有不同密度的剛性物體和水壩的雙向交互模擬效果。從效果圖可以看出，本文的表面重建算法可以重建出固流交出中的飛濺液滴這一細(xì)節(jié)特征，有效提高了流固交互模擬的真實(shí)感。

圖7 剛性圓柱和流體單相交互的流體表面重建

圖8 水壩沖擊不同密度物體的流體表面重建

3.2 表面重建性能分析

使用 3種不同方法：CAVW_2007[6]、CGF_2012[18]和本文方法對(duì)不同粒子間距采樣后的犰狳模型進(jìn)行了表面重建，并比較了內(nèi)存消耗。對(duì)比結(jié)果如表 2所示。表中給出了在不同粒子半徑 r和 Cube大小下不同重建方法的內(nèi)存消耗。內(nèi)存消耗通過(guò)在重建過(guò)程中需要計(jì)算距離場(chǎng)值的立方體個(gè)數(shù)(#Cube)來(lái)進(jìn)行衡量。CAVW_2007方法對(duì)流體空間進(jìn)行規(guī)則均勻劃分，在每一個(gè)Cube里都要計(jì)算距離場(chǎng)，因此內(nèi)存消耗是最大的；CGF_2012方法則只在流體表面附近的 Cube里計(jì)算距離場(chǎng)，大幅度降低了內(nèi)存消耗，但由于不能準(zhǔn)確地檢測(cè)出流體表面粒子和設(shè)定r和L的大小關(guān)系，所以不能重建出水膜的表面；針對(duì)本文方法同樣大幅度降低了內(nèi)存消耗，雖然在八叉樹的自適應(yīng)創(chuàng)建和訪問(wèn)引入了部分計(jì)算開銷，但整體計(jì)算性能與 CGF_2012方法持平，且能重建出由單層粒子組成的水膜結(jié)構(gòu)的表面。

表2 流體表面重建算法的內(nèi)存消耗對(duì)比

本文利用基于八叉樹的自適應(yīng)表面重建方法對(duì)所有三維場(chǎng)景進(jìn)行流體表面重建，并統(tǒng)計(jì)了性能數(shù)據(jù)，如表 3所示。在表 3中，#FP、#SP和#FPsurface分別代表流體粒子數(shù)、固體粒子數(shù)以及流體表面粒子數(shù)。#Cube3和#Cube1分別代表本方法和 CAVW_2007方法在流體表面重建過(guò)程中需要計(jì)算距離場(chǎng)值的Cube數(shù)量。表示本文方法與CAVW_2007方法的運(yùn)行時(shí)間比值。#Triangles表示利用Marching Cubes算法提取三角形網(wǎng)格面后得到的三角形數(shù)量。從表 3給出的統(tǒng)計(jì)數(shù)據(jù)可以看出，在三維流體模擬過(guò)程中，流體表面粒子數(shù)量通常只占流體粒子總數(shù)的很小一部分(＜20%)。

與 CAVW_2007方法相比，基于八叉數(shù)的自適應(yīng)表面重建算法在大幅度降低內(nèi)存消耗的同時(shí)，顯著提高了粒子流體表面重建的計(jì)算效率，適合大粒子規(guī)模下的流體表面的重建。

表3 基于八叉樹的自適應(yīng)表面重建算法的內(nèi)存消耗

本文的流體表面重建算法只在整個(gè)流體模擬的計(jì)算過(guò)程中引入了較小的額外計(jì)算量。還給出了上述三維場(chǎng)景模擬過(guò)程中每個(gè)時(shí)間步內(nèi)的平均計(jì)算時(shí)間，并對(duì)本文的表面重建方法進(jìn)行時(shí)間性能分析。時(shí)間統(tǒng)計(jì)包括：物理模擬計(jì)算時(shí)間Tps和基于八叉樹的自適應(yīng)表面重建時(shí)間 Tode，統(tǒng)計(jì)結(jié)果如表4所示。從表4數(shù)據(jù)可以看出，模擬過(guò)程中主要的計(jì)算瓶頸是復(fù)雜的物理模擬部分，占據(jù)了總模擬時(shí)間的 70%以上。而自適應(yīng)表面重建部分的計(jì)算時(shí)間包括流體表面粒子檢測(cè)時(shí)間和八叉樹自適應(yīng)距離場(chǎng)創(chuàng)建時(shí)間，平均占據(jù)了總模擬時(shí)間的10%左右。

表4 時(shí)間性能分析(s)

通過(guò)以上性能數(shù)據(jù)的統(tǒng)計(jì)分析可以得到，本文的自適應(yīng)流體表面重建方法可以大幅度降低拉格朗日粒子流體表面重建過(guò)程中的內(nèi)存消耗和計(jì)算時(shí)間。但本文方法在創(chuàng)建基于八叉樹的距離場(chǎng)時(shí)，對(duì)于每個(gè)葉子節(jié)點(diǎn)都要計(jì)算其 8個(gè)頂點(diǎn)的距離場(chǎng)值，存在計(jì)算冗余。后續(xù)工作中將研究八叉樹葉子節(jié)點(diǎn)的頂點(diǎn)共享問(wèn)題，降低計(jì)算冗余。

4 結(jié) 論

本文給出的基于八叉樹的自適應(yīng)流體表面重建方法可以高效地捕獲拉格朗日粒子流體的表面細(xì)節(jié)特征以進(jìn)行真實(shí)感繪制，該方法通過(guò)只在流體表面建立自適應(yīng)距離場(chǎng)，使得流體表面重建過(guò)程中的內(nèi)存消耗和計(jì)算復(fù)雜度，非常適合于虛擬現(xiàn)實(shí)應(yīng)用中大規(guī)模粒子流體場(chǎng)景的表面提取和繪制。

[1] 柳有權(quán). 基于物理的計(jì)算機(jī)動(dòng)畫及其加速技術(shù)的研究[D]. 北京: 中國(guó)科學(xué)院軟件研究所, 2005.

[2] Blinn J F. A generalization of algebraic surface drawing [J]. ACM Transactions on Graphics (TOG), 1982, 1(3): 235-256.

[3] Müller M, Charypar D, Gross M. Particle-based fluid simulation for interactive applications [C]//In Proceedings of the 2003 ACM SIGGRAPH, Eurographics Symposium on Computer Animation. New York: CAM Press, 2003: 154-159.

[4] Yu J H, Turk G. Reconstructing surfaces of particle-based fluids using anisotropic kernels [C]//In Proceedings of the 2010 ACM SIGGRAPH, Eurographics Symposium on Computer Animation. New York: CAM Press, 2010: 217-225.

[5] Zhu Y N, Bridson R. Animating sand as a fluid [J]. ACM Transactions on Graphics (TOG), 2005, 24(3): 965-972.

[6] Solenthaler B, Schlafli J, Pajarola R. A unified particle model for fluid-solid interactions [J]. Computer Animation and Virtual Worlds, 2007, 18(1): 69-82.

[7] Becker M, Teschner M. Weakly compressible SPH for free surface flows [C]//In Proceedings of the 2007 ACM SIGGRAPH, Eurographics Symposium on Computer Animation. New York: CAM Press, 2007: 209-217.

[8] Ando R, Thurey N, Tsuruno R. Preserving fluid sheets with adaptively sampled anisotropic particles [J]. IEEE Transactions on Visualization and Computer Graphics (TVCG), 2012, 18(8): 1202-1214.

[9] He X W, Liu N, Li S, et al. Local poisson SPH for viscous incompressible fluids [J]. Computer Graphics Forum, 2012, 31(6): 1948-1958.

[10] Ihmsen M, Cornelis J, Solenthaler B, et al. Implicit incompressible SPH [J]. IEEE Transactions on Visualization and Computer Graphics (TVCG), 2014, 20(3): 426-435.

[11] Solenthaler B, Pajarola R. Predictive-corrective incompressible SPH [J]. ACM Transactions on Graphics (TOG), 2009, 28(3): 40-46.

[12] Akinci N, Ihmsen M, Akinci G, et al. Versatile rigid-fluid coupling for incompressible SPH [J]. ACM Transactions on Graphics (TOG), 2012, 31(4): 621-628.

[13] Akinci N, Cornelis J, Akinci G, et al. Coupling elastic solids with SPH fluids [J]. Computer Animation and Virtual Worlds, 2013, 24(3-4): 195-203.

[14] Marrone S, Colagrossi A, Le Touze D, et al. Fast free-surface detection and level-set function definition in SPH solvers [J]. Journal of Computational Physics, 2010, 229(10): 3652-3663.

[15] Colagrossi A, Landrini M. Numerical simulation of interfacial flows by smoothed particle hydrodynamics [J]. Journal of Computational Physics, 2003, 191(2): 448-475.

[16] Shao X Q, Zhou Z, Zhang J S, et al. Realistic and stable simulation of turbulent details behind objects in smoothed-particle hydrodynamics fluids [J]. Computer Animation and Virtual Worlds, 2015, 26(1): 79-94.

[17] 張文勝, 解 騫, 鐘 瑾, 等. 基于八叉樹鄰域分析的光線跟蹤加速算法[J]. 圖學(xué)學(xué)報(bào), 2015, 36(3): 339-344.

[18] Akinci G, Ihmsen M, Akinci N, et al. Parallel surface reconstruction particle-based fluids [J]. Computer Graphics Forum, 2012, 31(6): 1797-1809.

An Efficient Surface Reconstruction Method for Lagrangian Particle Fluid

Shao Xuqiang1, Liu Yan2, Wang Xinying1

(1. School of Control and Computer Engineering, North China Electric Power University, Baoding Hebei 071003, China; 2. International Education College, Hebei Finance University, Baoding Hebei 071003, China)

In order to simulate small-scale geometrical features of Lagrangian particle fluid surface, this paper presents an efficient octree-based adaptive fluid surface reconstruction approach. First, this approach implements the more accurate detection of fluid surface particles. Then, the approach adaptively decomposes the fluid particles’ AABB bounding boxes including fluid surface particles based on octree. Finally, the implicit distance field is only constructed in the leaf nodes including fluid surface particles. The memory consumption and computation complexity of fluid surface reconstruction depend on fluid surface not fluid volume, so our method is very suitable for large scale fluid scenes.

Lagrangian particle; surface reconstruction; fluid simulation; realistic rendering

TP 391

10.11996/JG.j.2095-302X.2016050607

A

2095-302X(2016)05-0607-07

2016-03-15；定稿日期：2016-05-03

國(guó)家自然科學(xué)基金項(xiàng)目(61502168)；中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金(2015MS126)；河北省自然科學(xué)基金項(xiàng)目(F2016502069)

邵緒強(qiáng)(1982–)，男，山東泰安人，講師，博士研究生。主要研究方向?yàn)樘摂M現(xiàn)實(shí)、計(jì)算機(jī)圖形學(xué)。E-mail：shaoxuqiang@163.com