基于PBD框架的流體模擬

卿偉

（1.四川大學(xué)計算機學(xué)院，成都610065；2.四川大學(xué)視覺合成圖形圖像技術(shù)國防重點學(xué)科實驗室，成都610065）

0 引言

基于物理動畫的流體模擬是計算機圖形學(xué)中一個重要的研究領(lǐng)域。流體模擬應(yīng)用于自然界各種場合，如水流、泥石流、海浪等。模擬方法主要有兩大類：基于網(wǎng)格的歐拉法和基于粒子的拉格朗日法。歐拉法是在固定視點下跟蹤流體中固定的點并保存粒子的屬性；拉格朗日法是將流體視為流動的粒子，通過模擬大量的粒子運動來模擬物體的運動，每個粒子都包含相應(yīng)的物理屬性，如位置、密度、速度等。

在流體模擬中，基于粒子的拉格朗日法是常用的方法，SPH通過拉格朗日法來實現(xiàn)。將物體轉(zhuǎn)化為一群相互作用的粒子來模擬復(fù)雜場景，通過粒子的運動來模擬物體的運動。該方法的特點是簡單快速，且易于實現(xiàn)，方便應(yīng)用于各種實際場景中。但對實時性要求比較高的場景如游戲、電影等，SPH的模擬效果存在較大的局限性，如要求更小的時間步長，更多的鄰居粒子等。如果鄰居粒子數(shù)量太少，SPH將變得不穩(wěn)定，基本的解決方案是使用更小的時間步長或者增加鄰居粒子的數(shù)量，但會增加計算量，影響實時交互的效率。PBD可以彌補這些缺點，在表面增加一個力來解決鄰居粒子不足的問題，通過不可壓縮性約束來解決密度誤差問題，模擬的過程中通過雅克比迭代法逐步更改粒子的速度。

1 相關(guān)工作

傳統(tǒng)的SPH方法Muller[1]通過粒子的剛度狀態(tài)方程得出粒子間的粘度力及表面張力，計算出粒子間的速度及加速度來模擬流體的運動狀態(tài)。為了維持流體粒子間的不可壓縮性，標準的SPH方法及弱壓縮性SPH（WCSPH）Beckera and Teschner 2007[2]中所涉及到的剛度方程會導(dǎo)致張力過大，且對時間步長敏感，只在足夠小的時間步長內(nèi)模擬才會穩(wěn)定。Solenthaler and Pajarola 2009[3]提出預(yù)測-校正相結(jié)合的SPH方法（PCISPH）來緩解時間步長的限制，它使用雅可比迭代法來反復(fù)迭代得出更為精確的壓力變化值，使力收斂于某個值。

Bodin等人[3]提出將不可壓縮性作為系統(tǒng)的密度約束來處理，并以此來模擬密度均勻的流體。他指出使用線性約束函數(shù)來解決線性互補問題，線性約束可以通過Gauss-Seidel迭代來處理。但該方法無法解決非線性問題。PBD方法能很好的解決非線性問題，PBD直接在粒子上進行操作，對每個粒子通過雅可比迭代反復(fù)計算約束的差值及相應(yīng)梯度的差值，通過約束來糾正粒子的信息。

Muller等人[4]在實現(xiàn)游戲動態(tài)模擬時提出了PBD，PBD是基于Verlet積分法來實現(xiàn)的。它直接更新粒子的位置信息，通過Gauss-Seidel迭代來求解系統(tǒng)中的非線性約束，通過位置的更新來避免對力的直接處理。PBD對位置的修改是可控的，通過控制位置的改變不會因為粒子受力太大而導(dǎo)致不穩(wěn)定性；對時間步長的要求也很寬松，相比于SPH，PBD在更大的時間步長內(nèi)也能保持穩(wěn)定。

2 算法實現(xiàn)

2.1 基于PBD框架的PBF算法簡介

在PBD中，一個動態(tài)物體由N個粒子和M個約束表示，每個粒子包含各種物理屬性，如質(zhì)量、速度、位置等。每個約束由以下幾個方面組成：①n個候選粒子；②約束函數(shù)Cj；③粒子群的索引值④剛度系數(shù)⑤約束類型：相等或者不相等。若約束函數(shù)則約束類型相等；若則約束類型為不相等。剛度系數(shù)kj表示約束的強度。

PBD的基本流程如下：

（1）-（3）初始化粒子屬性。該算法的主要思想是（9）-（11）、（12）-（14）及（18）-（21），第（10）對每個粒子運用歐拉積分預(yù)估位置信息；第（12）-（14）通過迭代計算約束函數(shù)并糾正粒子的位置，使其滿足所有的約束。第（19）-（20）通過粒子的位置信息計算粒子的速度，并更新粒子的速度和位置信息。迭代中，每次循環(huán)只在一個時間步長內(nèi)更新粒子的速度和位置信息，這保證了PBD模擬的穩(wěn)定性。第（6）行通過一個外力來影響粒子的運動。第（8）行通過外力或者其他方法來降低速度，使其模擬效果更接近現(xiàn)實。第（22）行根據(jù)摩擦力和恢復(fù)系數(shù)處理粒子碰撞后的速度信息。

2.2 不可壓縮性

綜上可知PBD通過計算粒子的約束來模擬物理。PBF是基于PBD的基本思想，流體的約束主要考慮流體的不可壓縮性，即粒子的密度約束。密度約束指使每個粒子的密度保持在常量密度ρ0下。每個粒子的約束函數(shù)由粒子的位置及其鄰居粒子的位置信息組成。[Bondin2012]中對第i個粒子密度約束的定義：

其中，ρ0為靜態(tài)密度。ρi根據(jù)標準SHP得出的密度預(yù)估值。即：

在模擬中，假定所有流體粒子的質(zhì)量相同。對光滑核函數(shù)的計算中，密度的預(yù)估使用Poly6核系數(shù)，密度梯度的預(yù)估使用Spiky核系數(shù)。

PBD的目的是找出粒子位置的糾正值?p；因此把密度約束具體化為如下形式：

在梯度方向上，粒子信息的改變量是最大的，因此可以沿著梯度的方向來預(yù)估粒子的位置信息，通過牛頓迭代算法，得出如下公式：

在SPH上定義了粒子梯度的處理函數(shù)，將其運用到每個粒子k，可得到梯度的約束函數(shù)：

將其代入公式（5）,得出約束解的因子值λ，如下所示：

在相同的條件下，同性質(zhì)的粒子具有相同的約束，因此通過上式（8）可求出每個粒子的約束解。由于密度約束函數(shù)（1）是非線性函數(shù)，光滑核的邊界處會出現(xiàn)梯度消失的情況，因此光滑核邊界處的粒子是離散分布的，方程（8）的分母（梯度）趨于0，λ的值將趨于無窮大，由此將導(dǎo)致模擬的不穩(wěn)定。為了解決此問題，PCI?SPH給定一個校正值來處理，校正值的大小與鄰居粒子數(shù)量有關(guān)。在CFM中通過增加一個約束力來調(diào)整約束，其基本思想在約束函數(shù)中添加一個約束力。PBD也采用這種方法來解決核邊界處粒子不穩(wěn)定的情況。公式如下：

再根據(jù)自身的密度約束（λi）及鄰居粒子的密度約束（λj）來更新糾正位置Δpi：

2.3 表面張力

SPH模擬時粒子會發(fā)生聚簇現(xiàn)象，這是由于當一個粒子的鄰居粒子不足時將無法滿足密度約束條件，將產(chǎn)生一個負壓力使粒子成簇。為了解決此問題，通常的情況是增加一個非負的力來抵消負壓力，但會降低粒子間的凝聚而使粒子分散。Alduan,Otaduy[2011]提出使用有限元力（DEM）法使粒子間的距離為核半徑的一半；Schechter,Bridson[2012]提出在粒子的自由表面生成一個幽靈力。PBF采用[Monaghan 2000]提出的方法，在光滑核范圍內(nèi)增加一個拉力，將所有離散的粒子拉向光滑核范圍內(nèi)。PBF加入表面張力與無表面張力的對比效果圖如圖1所示。

表面張力的計算公式如下：

其中?q是光滑核內(nèi)粒子到核中心的固定距離，k為正常數(shù)。當取時模擬效果最好。將該公式代入方程（11）可得：

流體模擬時還需考慮流體本身存在一定的粘性，因此我們應(yīng)用XSPH粘度來更新粒子的速度信息。XSPH的計算公式如下：參數(shù)c可取0.01。

圖1 表面是否添加張力對比實驗

3 實驗結(jié)果

本實驗?zāi)M環(huán)境如表1：

表1

實驗的模擬結(jié)果如圖2、圖3所示。

本實驗的數(shù)據(jù)：光滑核半徑：0.004-0.02；時間間隔?t：0.0001-0.005。

實驗一：SPH與PBF時間步長的對比實驗

在實驗?zāi)M的過程中，當SPH和PBF的光滑核半徑h=0.016，?t在0.0044以后，SPH將逐漸變得不穩(wěn)定，最后所有的粒子會向四周散開，而PBF能一直保持在穩(wěn)定狀態(tài)。

從實驗數(shù)據(jù)中可以得出，SPH在h=2.5r時，即h=0.01時模擬將不穩(wěn)定；而PBF在h=2.2r時，即h=0.0088時仍然保持穩(wěn)定，直到h=2.1r時，即h=0.0084時將不穩(wěn)定。

實驗二：SPH與PBF光滑核半徑的對比

圖2 SPH模擬隨?t的變化情況

圖3 PBF模擬隨?t的變化情況

圖4 SPH隨h的變化情況，粒子半徑為r=0.004，h為核半徑

在相同的實驗環(huán)境下通過SPH和PBF兩種方式模擬流體，我們可以得出幾個結(jié)論：

（1）PBF模擬流體運動時，能滿足更大的時間步長；在時間步長足夠大時，PBF保持相對穩(wěn)定，而SPH會出現(xiàn)粒子像四周擴散。

（2）在其他相同條件下，PBF比SPH滿足更小光滑核函數(shù)。在足夠小的光滑核半徑內(nèi)，PBF模擬的穩(wěn)定性更好。

4 結(jié)語

綜上所述，基于PBD的方法比SPH能更好地模擬流體運動，其時間步長的要求更為寬松，對模擬的光滑核半徑也更小，所以在鄰居粒子很少的情況下，PBF也能很好地模擬流體。本文中只對流體進行模擬，沒有處理流體對剛體的交互效果，后期的工作主要處理流體對剛體的碰撞效果以及剛體之間的碰撞檢測效果，同時對PBD的一些限制條件也需要及時驗證。

圖5 PFB隨h的變化情況，粒子半徑為r=0.004，h為核半徑

參考文獻：

[1]Ivan Alduan，Miguel A Otaduy．SPH Granular Flow with Friction and Cohesion.In Proceedings of the 2011 ACM SIGGRAPH/Eurographics Symposiumon Computer Animation,pages 25-32.ACM,2011.

[2]Markus Becker，Matthias Teschner．Weakly Compressible SPH for Free Surface Flows.In Proceedings of the 2007 ACM SIGGRAPH/Eurographics Symposium on Computeranimation,pages209-217.Eurographics Association,2007.

[3]Kenneth Bodin,Claude Lacoursiere,Martin Servin.Constraint Fluids.IEEE Transactions on Visualization and Computer Graphics,18（3）:516-526,2012.

[4]Simon Clavet,Philippe Beaudoin,Pierre Poulin.Particle-based Viscoelastic Fluid Simulation.In Proceedings of the 2005 ACM SIGGRAPH/Eurographics Symposium on Computer Animation,pages 219-228.ACM,2005.

[5]Jeong-Mo Hong,Ho-Young Lee,Jong-Chul Yoon,Chang-Hun Kim.Bubbles alive.In ACM Transactions on Graphics（TOG）,volume27,page48.ACM,2008.

[6]Joseph JMonaghan.Sphwithouta Tensile Instability.Journal of Computational Physics,159（2）:290-311,2000.

[7]Matthias Muller,David Charypar,Markus Gross.Particle-based Fluid Simulation for Interactive Applications.In Proceedings of the 2003 ACM SIGGRAPH/Eurographics Symposium on Computer Animation,pages 154-159.Eurographics Association,2003.

[8]Matthias Muller,Bruno Heidel berger,Marcus Hennix,John Ratcliff.Position Based Dynamics.Journal of Visual Communication and Image Representation,18（2）:109-118,2007.

[9]Hagit Schechter and Robert Bridson.GhostSPH for Animating Water.ACM Transactions on Graphics（TOG）,31（4）:61,2012.

[10]Barbara Solenthaler，Renato Pajarola.Predictive-Corrective Incompressible SPH.In ACM Transactions on Graphics（TOG）,volume28,page40.ACM,2009.