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

邵緒強， 劉 艷， 王新穎

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

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

邵緒強1， 劉 艷2， 王新穎1

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

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

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

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

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

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

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

1 拉格朗日粒子流體模擬

1.1 不可壓縮粒子流體模擬

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

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

圖1 PICSPH偽代碼

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

在預測下一步的密度時，應該根據預測的位置信息更新鄰居粒子列表?？紤]到計算效率，這里不再重新搜索鄰居粒子，而是通過更新到現有鄰居粒子距離的方式來近似。雖然這樣的估計會帶來誤差，但是在最少迭代次數、時間步長、核函數影響范圍衰減等約束條件的作用下，對整個流體運動模擬結果影響是可以忽略的。采用一個最小迭代次數的方法可顯著減少誤差的局部震蕩，而且這給予粒子足夠的時間來擴散預測信息。為了提高計算效率，一般使用最少迭代3次，而3次迭代已經足夠將粒子密度誤差約束到目標范圍內。

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

1.2 流固交互模擬

對于SPH粒子流體和運動固體邊界的交互，本文結合文獻[12]的懲罰力交互方法和文獻[13]的自適應邊界粒子采樣方法實現雙向流固交互模擬，計算過程分為3步：

(1) 在每個時間步內進行流固交互計算之前，根據初始流體粒子間距r0對任一表面三角形進行自適應邊界粒子采樣，如圖 2所示。首先選取三角形頂點為邊界粒子；然后以初始流體粒子間距r0對三角形的邊進行邊界粒子采樣；確定三角形最短邊es和最長邊el，并選取向內的法線方向為掃描線方向對三角形內部進行邊界粒子采樣。自適應邊界粒子采樣可以更精確地表示固體表面的拓撲結構，從而實現更準確的固流交互。

圖2 自適應收縮壓粒子采樣

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

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

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

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

2 粒子流體表面重建

2.1 流體表面粒子檢測

流體的表面可以由位于流體表面的粒子表示和跟蹤?；谟嬎懔黧w動力學(computational fluid dynamics，CFD)領域文獻[14]方法，給出一個快速準確的流體表面粒子檢測方法。對于任一流體粒子i，首先計算一個再歸一化矩陣Bi

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

然后，計算矩陣Bi的最小特征值，并根據Bi和實驗驗證值的大小關系判定粒子 i是否屬于SPH流體表面S

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

2.2 自適應距離場重建

為了降低粒子流體表面重建過程中的內存消耗和計算復雜度，基于文獻[16]中流體表面粒子檢測和隱式表面定義方法，給出一種基于八叉樹[17]的自適應距離場構建算法，該方法只在包含流體表面粒子的葉子結點上計算距離場值。

在自適應距離場構建算法中，任一八叉樹結點包含兩個粒子集合：流體表面粒子集合和流體粒子集合。其中，流體表面粒子集合包含所有屬于該結點的流體表面粒子，流體粒子集合包含所有屬于該結點的流體粒子。本方法把八叉樹的結點分為 3類：外部結點、內部結點和表面結點。對于外部結點來說，其粒子集合流體表面粒子集合和流體粒子集合均為空；內部結點的粒子集合流體表面粒子集合為空，而粒子集合流體粒子集合非空；而對于表面結點，其粒子集合流體表面粒子集合非空。八叉樹距離場的遞歸創建過程如圖3所示。

圖3 基于八叉樹的自適應距離場的建立

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

半徑； R= 4r； f∈ [0… 1 ]并通過下式進行計算

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

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

圖4 粒子半徑r與葉子結點邊長L的大小關系

3 實例驗證

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

表1 漩渦模擬中物理參數及取值

3.1 表面重建效果分析

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

圖5 自適應流體表面重建算法

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

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

圖7給出了剛性圓柱體以10 m/s的速度在水面移動產生的漩渦細節模擬結果。從模擬結果可以看出，本文的表面重建算法可以重建出剛性圓柱體后面產生的豐富漩渦細節特征，增強了流固交互模擬的真實感。

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

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

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

3.2 表面重建性能分析

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

表2 流體表面重建算法的內存消耗對比

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

與 CAVW_2007方法相比，基于八叉數的自適應表面重建算法在大幅度降低內存消耗的同時，顯著提高了粒子流體表面重建的計算效率，適合大粒子規模下的流體表面的重建。

表3 基于八叉樹的自適應表面重建算法的內存消耗

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

表4 時間性能分析(s)

通過以上性能數據的統計分析可以得到，本文的自適應流體表面重建方法可以大幅度降低拉格朗日粒子流體表面重建過程中的內存消耗和計算時間。但本文方法在創建基于八叉樹的距離場時，對于每個葉子節點都要計算其 8個頂點的距離場值，存在計算冗余。后續工作中將研究八叉樹葉子節點的頂點共享問題，降低計算冗余。

4 結 論

本文給出的基于八叉樹的自適應流體表面重建方法可以高效地捕獲拉格朗日粒子流體的表面細節特征以進行真實感繪制，該方法通過只在流體表面建立自適應距離場，使得流體表面重建過程中的內存消耗和計算復雜度，非常適合于虛擬現實應用中大規模粒子流體場景的表面提取和繪制。

[1] 柳有權. 基于物理的計算機動畫及其加速技術的研究[D]. 北京: 中國科學院軟件研究所, 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]. 圖學學報, 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

國家自然科學基金項目(61502168)；中央高?；究蒲袠I務費專項資金(2015MS126)；河北省自然科學基金項目(F2016502069)

邵緒強(1982–)，男，山東泰安人，講師，博士研究生。主要研究方向為虛擬現實、計算機圖形學。E-mail：shaoxuqiang@163.com