基于PBD框架的流體模擬

卿偉

（1.四川大學(xué)計(jì)算機(jī)學(xué)院，成都610065；2.四川大學(xué)視覺(jué)合成圖形圖像技術(shù)國(guó)防重點(diǎn)學(xué)科實(shí)驗(yàn)室，成都610065）

0 引言

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

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

1 相關(guān)工作

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

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

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

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

2.1 基于PBD框架的PBF算法簡(jiǎn)介

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

PBD的基本流程如下：

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

2.2 不可壓縮性

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

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

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

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

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

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

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

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

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

2.3 表面張力

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

表面張力的計(jì)算公式如下：

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

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

圖1 表面是否添加張力對(duì)比實(shí)驗(yàn)

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

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

表1

實(shí)驗(yàn)的模擬結(jié)果如圖2、圖3所示。

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

實(shí)驗(yàn)一：SPH與PBF時(shí)間步長(zhǎng)的對(duì)比實(shí)驗(yàn)

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

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

實(shí)驗(yàn)二：SPH與PBF光滑核半徑的對(duì)比

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

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

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

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

（1）PBF模擬流體運(yùn)動(dòng)時(shí)，能滿足更大的時(shí)間步長(zhǎng)；在時(shí)間步長(zhǎng)足夠大時(shí)，PBF保持相對(duì)穩(wěn)定，而SPH會(huì)出現(xiàn)粒子像四周擴(kuò)散。

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

4 結(jié)語(yǔ)

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

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

參考文獻(xiàn)：

[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.