基于PBD框架的流体模拟

卿伟

（1.四川大学计算机学院，成都610065；2.四川大学视觉合成图形图像技术国防重点学科实验室，成都610065）

0 引言

基于物理动画的流体模拟是计算机图形学中一个重要的研究领域。流体模拟应用于自然界各种场合，如水流、泥石流、海浪等。模拟方法主要有两大类：基于网格的欧拉法和基于粒子的拉格朗日法。欧拉法是在固定视点下跟踪流体中固定的点并保存粒子的属性；拉格朗日法是将流体视为流动的粒子，通过模拟大量的粒子运动来模拟物体的运动，每个粒子都包含相应的物理属性，如位置、密度、速度等。

在流体模拟中，基于粒子的拉格朗日法是常用的方法，SPH通过拉格朗日法来实现。将物体转化为一群相互作用的粒子来模拟复杂场景，通过粒子的运动来模拟物体的运动。该方法的特点是简单快速，且易于实现，方便应用于各种实际场景中。但对实时性要求比较高的场景如游戏、电影等，SPH的模拟效果存在较大的局限性，如要求更小的时间步长，更多的邻居粒子等。如果邻居粒子数量太少，SPH将变得不稳定，基本的解决方案是使用更小的时间步长或者增加邻居粒子的数量，但会增加计算量，影响实时交互的效率。PBD可以弥补这些缺点，在表面增加一个力来解决邻居粒子不足的问题，通过不可压缩性约束来解决密度误差问题，模拟的过程中通过雅克比迭代法逐步更改粒子的速度。

1 相关工作

传统的SPH方法Muller[1]通过粒子的刚度状态方程得出粒子间的粘度力及表面张力，计算出粒子间的速度及加速度来模拟流体的运动状态。为了维持流体粒子间的不可压缩性，标准的SPH方法及弱压缩性SPH（WCSPH）Beckera and Teschner 2007[2]中所涉及到的刚度方程会导致张力过大，且对时间步长敏感，只在足够小的时间步长内模拟才会稳定。Solenthaler and Pajarola 2009[3]提出预测-校正相结合的SPH方法（PCISPH）来缓解时间步长的限制，它使用雅可比迭代法来反复迭代得出更为精确的压力变化值，使力收敛于某个值。

Bodin等人[3]提出将不可压缩性作为系统的密度约束来处理，并以此来模拟密度均匀的流体。他指出使用线性约束函数来解决线性互补问题，线性约束可以通过Gauss-Seidel迭代来处理。但该方法无法解决非线性问题。PBD方法能很好的解决非线性问题，PBD直接在粒子上进行操作，对每个粒子通过雅可比迭代反复计算约束的差值及相应梯度的差值，通过约束来纠正粒子的信息。

Muller等人[4]在实现游戏动态模拟时提出了PBD，PBD是基于Verlet积分法来实现的。它直接更新粒子的位置信息，通过Gauss-Seidel迭代来求解系统中的非线性约束，通过位置的更新来避免对力的直接处理。PBD对位置的修改是可控的，通过控制位置的改变不会因为粒子受力太大而导致不稳定性；对时间步长的要求也很宽松，相比于SPH，PBD在更大的时间步长内也能保持稳定。

2 算法实现

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

在PBD中，一个动态物体由N个粒子和M个约束表示，每个粒子包含各种物理属性，如质量、速度、位置等。每个约束由以下几个方面组成：①n个候选粒子；②约束函数Cj；③粒子群的索引值④刚度系数⑤约束类型：相等或者不相等。若约束函数则约束类型相等；若则约束类型为不相等。刚度系数kj表示约束的强度。

PBD的基本流程如下：

（1）-（3）初始化粒子属性。该算法的主要思想是（9）-（11）、（12）-（14）及（18）-（21），第（10）对每个粒子运用欧拉积分预估位置信息；第（12）-（14）通过迭代计算约束函数并纠正粒子的位置，使其满足所有的约束。第（19）-（20）通过粒子的位置信息计算粒子的速度，并更新粒子的速度和位置信息。迭代中，每次循环只在一个时间步长内更新粒子的速度和位置信息，这保证了PBD模拟的稳定性。第（6）行通过一个外力来影响粒子的运动。第（8）行通过外力或者其他方法来降低速度，使其模拟效果更接近现实。第（22）行根据摩擦力和恢复系数处理粒子碰撞后的速度信息。

2.2 不可压缩性

综上可知PBD通过计算粒子的约束来模拟物理。PBF是基于PBD的基本思想，流体的约束主要考虑流体的不可压缩性，即粒子的密度约束。密度约束指使每个粒子的密度保持在常量密度ρ0下。每个粒子的约束函数由粒子的位置及其邻居粒子的位置信息组成。[Bondin2012]中对第i个粒子密度约束的定义：

其中，ρ0为静态密度。ρi根据标准SHP得出的密度预估值。即：

在模拟中，假定所有流体粒子的质量相同。对光滑核函数的计算中，密度的预估使用Poly6核系数，密度梯度的预估使用Spiky核系数。

PBD的目的是找出粒子位置的纠正值∆p；因此把密度约束具体化为如下形式：

在梯度方向上，粒子信息的改变量是最大的，因此可以沿着梯度的方向来预估粒子的位置信息，通过牛顿迭代算法，得出如下公式：

在SPH上定义了粒子梯度的处理函数，将其运用到每个粒子k，可得到梯度的约束函数：

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

在相同的条件下，同性质的粒子具有相同的约束，因此通过上式（8）可求出每个粒子的约束解。由于密度约束函数（1）是非线性函数，光滑核的边界处会出现梯度消失的情况，因此光滑核边界处的粒子是离散分布的，方程（8）的分母（梯度）趋于0，λ的值将趋于无穷大，由此将导致模拟的不稳定。为了解决此问题，PCI⁃SPH给定一个校正值来处理，校正值的大小与邻居粒子数量有关。在CFM中通过增加一个约束力来调整约束，其基本思想在约束函数中添加一个约束力。PBD也采用这种方法来解决核边界处粒子不稳定的情况。公式如下：

再根据自身的密度约束（λi）及邻居粒子的密度约束（λj）来更新纠正位置Δpi：

2.3 表面张力

SPH模拟时粒子会发生聚簇现象，这是由于当一个粒子的邻居粒子不足时将无法满足密度约束条件，将产生一个负压力使粒子成簇。为了解决此问题，通常的情况是增加一个非负的力来抵消负压力，但会降低粒子间的凝聚而使粒子分散。Alduan,Otaduy[2011]提出使用有限元力（DEM）法使粒子间的距离为核半径的一半；Schechter,Bridson[2012]提出在粒子的自由表面生成一个幽灵力。PBF采用[Monaghan 2000]提出的方法，在光滑核范围内增加一个拉力，将所有离散的粒子拉向光滑核范围内。PBF加入表面张力与无表面张力的对比效果图如图1所示。

表面张力的计算公式如下：

其中∆q是光滑核内粒子到核中心的固定距离，k为正常数。当取时模拟效果最好。将该公式代入方程（11）可得：

流体模拟时还需考虑流体本身存在一定的粘性，因此我们应用XSPH粘度来更新粒子的速度信息。XSPH的计算公式如下：参数c可取0.01。

图1 表面是否添加张力对比实验

3 实验结果

本实验模拟环境如表1：

表1

实验的模拟结果如图2、图3所示。

本实验的数据：光滑核半径：0.004-0.02；时间间隔∆t：0.0001-0.005。

实验一：SPH与PBF时间步长的对比实验

在实验模拟的过程中，当SPH和PBF的光滑核半径h=0.016，∆t在0.0044以后，SPH将逐渐变得不稳定，最后所有的粒子会向四周散开，而PBF能一直保持在稳定状态。

从实验数据中可以得出，SPH在h=2.5r时，即h=0.01时模拟将不稳定；而PBF在h=2.2r时，即h=0.0088时仍然保持稳定，直到h=2.1r时，即h=0.0084时将不稳定。

实验二：SPH与PBF光滑核半径的对比

图2 SPH模拟随∆t的变化情况

图3 PBF模拟随∆t的变化情况

图4 SPH随h的变化情况，粒子半径为r=0.004，h为核半径

在相同的实验环境下通过SPH和PBF两种方式模拟流体，我们可以得出几个结论：

（1）PBF模拟流体运动时，能满足更大的时间步长；在时间步长足够大时，PBF保持相对稳定，而SPH会出现粒子像四周扩散。

（2）在其他相同条件下，PBF比SPH满足更小光滑核函数。在足够小的光滑核半径内，PBF模拟的稳定性更好。

4 结语

综上所述，基于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.