一种改进的定位流体方法

闫志伟， 朱晓临， 刘洋洋， 梁欣鑫（合肥工业大学数学学院，合肥230009）



一种改进的定位流体方法

闫志伟， 朱晓临， 刘洋洋， 梁欣鑫
（合肥工业大学数学学院，合肥230009）

［摘 要］针对模拟流体，提出了一种基于定位流体方法的改进模型，继承了定位动态学方法的稳定性以及光滑粒子流体动力学方法的光滑性，通过将流体粒子按速度值动态分配支集半径和给流体粒子添加密度限制来增加计算准确性并加强流体的表面细节的模拟效果．

［关键词］流体模拟；光滑粒子动力学；定位动态学；定位流体

1　引 言

流体模拟是计算机图形学的一个重要研究领域，自从20多年前粒子系统引入计算机图形学领域以来，基于粒子方法的流体模拟技术一直都在不断的发展和改善．早在1992年Monaghan就提出了SPH （Smoothed Particles Hydrodynamics）方法模拟流体［1］．谭小辉、万旺根［2］用SPH方法模拟了水流倒入水杯的动态过程．Frank Losasso［3］结合SPH和PLS（Particle Level Set）两种方法模拟了大场景的逼真海浪效果．Fengquan Zhang［4］改进了SPH方法中的压力项的求解公式，近似保证了水流模拟的不可压缩性，模拟出来的液体飞溅的效果也很好．Solenthaler［5］用预测校正粒子的压力方法来确保模拟流体的不可压缩性，这种处理巧妙地避免了求解压力项中复杂的泊松方程，同时还可以用大时间步长来模拟．Hagit Schechter［6］在粒子表面人工添加了一层“幽灵粒子”（Ghost Particle）防止粒子模拟的非常规聚集和泼溅效果．Kenneth［7］用限制流体（Constraint Fluids）的方法，不在对SPH中的Navier－Stokes方程求解，而是替换成一个较为简单的密度限制方程，虽然该方法可以加快流体模拟的速度，但是模拟出来的效果一般．Matthias Müller［8］提出了定位动态学方法（Position Based Dynamics，PBD），其直接在限制流体方法中更新粒子的位置值而不是像原来的类SPH方法中的去更新速度值．这种方法有效的避免了粒子的碰撞问题，具有很好的稳定性，模拟出来的效果也很真实．Matthias Müller［9］又提出了定位流体（Position Based Fluids，PBF）的方法；该方法基于定位动态学，用SPH方法来计算PBD方法中的密度限制的迭代过程，而且在处理压力项时，采取添加人工压力项的方式大大简化了运算过程，还改善了粒子模拟的表面张力效果，使得最终模拟出来的效果极为逼真．

在大场景的流体模拟中，定位流体方法不仅计算量小而且模拟效果真实．然而，这个计算量小只是相对于其他类SPH方法来说的，其本身的计算量还是很大；因此该模拟方法对计算机硬件的要求很高．其次，效果的逼真也是建立在大规模的粒子运动基础之上，若单就细节模拟的效果来对比，定位流体模拟是不优于SPH方法的．本文改进的目的就是为了提高粒子细节模拟效果的同时又能较好的减少粒子运动的计算量．

2　定位流体方法简介

流体模拟要求的不可压缩性体现在密度上就是要保证粒子的密度值不变，定位流体的思想就是当粒子不能保证密度不变时对粒子的位置进行修正直到该粒子所处的位置基本达到所要求其密度在限定的误差范围内．首先对粒子i，令

其中pi表示粒子i的位置，ρ0是粒子的静止密度，而粒子i的密度ρi由SPH方法中求密度公式计算

其中pj表示粒子i的支集半径内除粒子i外粒子j的位置，j=1，2，…，n;h是支集半径，表示以粒子i为中心，h为半径的范围．mj表示粒子j的质量．求密度的核函数为

其中r=pi－pj是两个粒子位置的位矢，r=|r|．

接下来要找出一个粒子的位置修正值Δp，使之满足等式

为了更快速地得到粒子位置校正值Δp，可以限定Δp与函数梯度C(p)的方向一致，这样就可以确定一个标量λ满足

再由(3)式可以进一步得出

即

此处的梯度用下面的核函数来近似估计

综上可得

公式(8)在粒子支集半径内的粒子很少的时候会导致数值的不稳定，因此不加选择的使用外力混合的方法

ε是一个调节参数;此时公式(8)变为

为了加速粒子位置校正值的更新，把支集半径内的粒子的密度限制系数λj也考虑到粒子位置的变化因素中，此时Δp的值就变为

在SPH方法中，还有一个常见的问题就是粒子总是会有聚成一团的情形，这主要是因为粒子支集半径内的粒子数较少，某些粒子的压力项为负值形成的．如果给粒子人工添加一个非负的压力项，这个问题就可以避免．这里加入的压力项也是一个平滑核函数

其中当

时，计算出的数据稳定性就比较理想．此时公式(11)变为

这样通过数次迭代后，进行位置校正过的粒子密度会控制在我们可以接受的误差范围内，并且粒子会将支集半径内的粒子拉向自己，产生模拟的流体表面张力效果．

后面为了保持模拟过程的能量守恒，还需要给粒子加入涡流限制(vorticityconfinement)．添加每个粒子的旋转力形式如下

其中ωi表示旋量，v表示速度表示粘性力，ε*是一个调节参数．此处的核函数计算估计形式如下

最后为了保证每个粒子的连续运动的实时效果，给粒子人工添加了适当的粘性力

3　定位流体的改进

(i)在计算公式(1)时，当粒子出现在边界或与固体碰撞时，其支集半径内的粒子数过少，会导致计算密度限制方程精度很低．因此我们将边界处的粒子的密度限制方程改写成如下形式

(ii)虽然PBF方法稳定性很好，进行大规模场景的流体模拟时可以降低对粒子支集半径的要求，但是PBF方法的计算量还是很大．因此，我们提出如下改进．

在常规流体模拟的场景中，除边界外流体都是连续的．从粒子的角度看，一个内部粒子与其邻近粒子的属性值的差别很小，而正是利用这种微小的差别，我们可以在计算相对“安静”的粒子属性值时降低计算量．根据粒子速度值的大小，我们将所有粒子分为三种:活跃粒子(activeparticle)、普通粒子(ordinaryparticle)、惰性粒子(inactiveparticle)．计算活跃粒子时，用较大的迭代次数和支集半径;相对地，计算惰性粒子时，则用较小的迭代次数和支集半径．

上述改进的第1步虽然增加了少许的计算量，但却可以加强模拟的细节效果;改进的第2步则是在减少计算量的同时，增强粒子边界、表面处模拟的细节效果．

4　程序步骤

步骤4 计算支集半径内每个粒子的λi；

步骤5 根据公式（13）计算支集半径内每个粒子的Δpi；

步骤6 通过迭代进一步的校正Δpi；

步骤9 对更新的速度添加涡流限制和粘性力．

5　效果对比

实验平台是Windows 8系统，程序开发的软件选用Microsoft Visual Studio 2008，三维图形的编程接口为OpenGL．算法实现的硬件环境是Intel（R）Core（TM）i3－4130的CPU、4G内存的PC机．下图给出了我们的实验结果．

图1　模拟了粒子从左侧向右流动趋近静止的运动状态

对比图分别为第4，15，33帧图像．左：定位流体方法，帧率为2．0080，；右：改进后的定位流体方法，帧率为2．1276．从上面的对比图中，我们可以很明显的看出改进后的方法粒子的分布位置更自然，针对5500粒子能够计算效率提升近7%．

最后，使用Marching Cubes方法进行表面重建，其效果如图2所示．

6　总 结

定位流体方法是目前适应大场景流体模拟最好的方法，同时也是最耗时的方法．我们的实验结果表明文章中的改进确实可以较大幅度的减少定位流体方法的计算量，同时还能够较好的提高模拟的表面细节效果．

流体模拟的不可压缩性是很难保证的［10］，除去方法计算的误差，模拟过程的数值耗散和能量守恒问题大多由此而来．定位流体方法也只能做到模拟过程的近似不可压缩．如何更好的去解决流体模拟的不可压缩性问题也一直是诸多流体模拟方法不断发展与创新的关键所在［11－12］．本文虽然给出了提高计算速度的改进定位流体方法，但是对于定位流体本身巨大的计算量来说还是略显不足．我们未来的工作方向将主要放在具体的真实感绘制和更有效的解决定位流体的数值耗散问题，这方面文献［13］对自由表面的再构建的创新方法将对我们有所帮助．

图2　水流的仿真图

［参 考 文 献］

［1］ Monaghan J J．Smoothed particle hydrodynamics［J］．Annual review of astronomy and astrophysics，1992，30：543－574．

［2］ 谭小辉，万旺根，黄炳，崔滨．基于SPH的三维流体模拟［J］．计算机应用于软件，2009，26（12）：222－258．

［3］ Losasso F，Talton J O，Kwatra N，Fedkiw R．Two－Way Coupled SPH and Particle Level Set Fluid Simulation ［J］．IEEE Translations on visualization and computer graphics，2008，14（4）：797－804．

［4］ Zhang F，Shen X，Long X，et al．Real－time particle fluid simulation with wcsph［J］．Pacific Graphics Short Papers，2012：29－34．

［5］ Solenthaler B，Pajarola R．Predictive－Corrective Incompressible SPH［J］．ACMSIGGRAPH，2009，28（3）：401－406．

［6］ Schechte H，Bridson R．Ghost SPH for Animating Water［J］．ACM Transactions on Graphics（TOG），2012，31（4）：61－69．

［7］ Bodin K，Lacoursiere C，Servin M．Constraint Fluids［J］．IEEE Transactions on visualization and computer graphics，2012，18（3）：516－527．

［8］ Muller M，Charypar D．Gross M．Particle－Based Fluid Simulation for Interactive Applications［J］．Eurographics／SIGGRAPH Symposium on Computer Animation，2003，Proceeding SCA’03：154－159．

［9］ Macklin M，Müller M．Position based fluids［J］．ACM Transactions on Graphics（TOG），2013，32（4）：104－109．

［10］ Khayyer A，Gotoh H，Shao S D．Development of CISPH method for accurate water－surface tracking in breaking waves［J］．Coastal Engineering，2008，55（3）：236－273．

［11］ Solenthaler B，Gross M．Two－scale particle simulation［C］．ACM Transactions on Graphics（TOG）．ACM，2011，30（4）：81－88．

［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）：62－70．

［13］ Yu J，Turk G．Reconstructing surfaces of particle－based fluids using anisotropic kernels［J］．ACM Transactions on Graphics（TOG），2013，32（1）：5－17．

An Improved Method of Position Based Fluids

YAN Zhi－wei， ZHU Xiao－lin， LIU Yang－yang， LIANG Xin－xin
（School of Mathematics，Hefei University of Technology，Hefei 230009，China）

Abstract：This paper puts forward an improved model for fluid simulation according to Position Based Fluids，which inherits the stability of Position Based Dynamics（PBD）method and the smoothness of Smoothed Particles Hydrodynamics（SPH）method．In this method，we distrabuted different supported radius for different particles according to velocity and added density constraints for paticles，which improved calculation accuracy and simulation effects．

Key words：fluid simulation；smoothed particles hydrodynamics（SPH）；position based dynamics（PBD）；position based fluids（PBF）

［基金项目］国家自然科学基金资助项目（61272024）

［收稿日期］2015－03－15

［中图分类号］TP391

［文献标识码］A

［文章编号］1672－1454（2016）01－0038－06