【CG物理模拟系列】流体模拟--粒子法之Position Based Fluids

版权声明:本文为博主原创文章,未经博主允许不得转载。 https://blog.csdn.net/Knight_Lyh/article/details/56282584

Position Based Fluids

我们使用粒子法进行流体模拟时,一个很重要的因素就是如何保持流体的非压缩性。非压缩性是指,随着流体的流动,密度在空间和时间上不发生变化,比如像水这样的流体,另外,与音速相比流速足够小的空气也可以当成非压缩性流体来考虑 (接近音速的时候,空气会发生压缩,音速飞机产生的冲击波即音爆便是由空气压缩产生的)。为了保持非压缩性,学者们提出了各种各样的方法,如MPS,WCSPH,PCISPH,IISPH等, 这之中,能非常高速且稳定计算的方法便是我们这次要讲的基于位置流体法(Position Based Fluids, 以下简称PBF)[1]。

PBF是基于位置法(Position-based method)[2] 中的一种。Position-based法与以前的基于力学方法不同,是根据物体位置(空间坐标值)间的制约条件直接更新物体(传统的基于力学法则是根据支配方程,按照加速度->速度->位置的顺序逐步变更),因此不会产生在传统力学方法中由逐步积分所产生的误差积累。但是,由于这种方法中,顶点及粒子的位置关系所遵循的制约条件是属于几何学的范畴,严格意义上讲不属于物理模拟。

Position-based法介绍

这里先对基本的position-based法进行说明。首先,对于一个产生位移的点,其坐标为pbf.eq1.gifpbf.eq2.gif是其制约条件。 正常状态下pbf.eq3.gif, 当粒子发生一个微小量pbf.eq4.gif的移动后,也满足这个条件,即 pbf.eq5.gif。 把这个式子泰勒展开,即

pbf.eq6.gif

忽略掉2次以上的项,我们将得到一个关于pbf.eq4.gif的方程。

pbf.eq7.gif

pbf.eq8.gif可以由制约条件的定义计算。现在,把pbf.eq1.gif当作粒子的中心坐标,在3次元空间粒子数为pbf.eq9.gif的条件下,需求得未知数的个数为pbf.eq10.gif。 如果每个粒子都分别计算一次pbf.eq11.gif的话(也就是pbf.eq12.gif), 上述方程个数为pbf.eq9.gif,未知数个数是pbf.eq10.gif,我们很难轻易求出解。因此,需要针对pbf.eq4.gif添加额外的制约条件。

假设点pbf.eq1.gif是按照某运动方程发生移动,而pbf.eq4.gif则是按照不同于运动方程的其他的制约条件pbf.eq11.gif来移动的。这样的话,就需要一个不妨碍原来的运动(平行移动和旋转)的移动方向,也就是要保留平行移动和旋转量的移动方向。 考虑到制约条件pbf.eq11.gif要对运动状态不产生影响,于是pbf.eq13.gif的方向便与原来的运动方向垂直。也就是,pbf.eq4.gif的方向与pbf.eq13.gif相同。 导入一个拉格朗日乘数pbf.eq14.gif,使得

pbf.eq15.gif

带入上面的关于pbf.eq4.gif的方程,整理得,

pbf.eq16.gif

然后把求得的pbf.eq14.gif带入pbf.eq4.gif的式子,可以得到满足条件的移动向量,进而可以更新点的位置。

Position-based法在粒子法中的应用

首先,我们要想出一个满足非压缩性条件的制约条件(Constraint)。假设流体初始密度是pbf.eq18.gif,流体不发生压缩就是指其密度pbf.eq17.gif一直保持不变,也就是要满足pbf.eq19.gif。把它转换成关于位置pbf.eq1.gif的制约pbf.eq11.gif,如下
  pbf.eq20.gif    ・・・(1)

然后,按照粒子法之一的SPH法,求解某粒子pbf.eq21.gif的密度pbf.eq22.gif

pbf.eq23.gif

这里的pbf.eq24.gif是粒子pbf.eq21.gif的临近粒子的集合,pbf.eq25.gif是粒子质量,pbf.eq26.gif是有效半径。在临近粒子数为pbf.eq9.gif的时候,与粒子pbf.eq21.gif的密度相关的制约,是一个与临近粒子坐标pbf.eq27.gif相关的函数,pbf.eq14.gif的计算公式如下:
  pbf.eq28.gif    ・・・(2)

pbf.eq29.gif可由把密度pbf.eq22.gif的公式代入式(1)求得。 pbf.eq30.gif如下,

pbf.eq32.gif

这里要注意的是,由于pbf.eq31.gif是与pbf.eq33.gifpbf.eq34.gif(和pbf.eq26.gif)相关的函数,满足pbf.eq35.gif的微分为0的条件只有pbf.eq36.gifpbf.eq37.gif。于是上式可以改写成下面这样:
  pbf.eq38.gif    ・・・(3)

pbf.eq36.gif的时候由于pbf.eq40.gif包含了pbf.eq35.gif所以直接计算便可。 pbf.eq41.gif时,由于只有pbf.eq42.gif的时候包含pbf.eq40.gif,所以只计算这一种情况即可 (这里的负号是由于pbf.eq43.gif)。另外,实际计算的时候一般都会把所有的粒子质量pbf.eq25.gif设为定值,简化运算。

由式(2)和式(3)可以计算出各粒子的拉格朗日乘数pbf.eq45.gif,下面我们来求各粒子的位置修正量pbf.eq46.gif。为使得对临近粒子的影响以及受临近粒子的影响保持对称(和SPH法的圧力项同理),我们用如下公式计算:
  pbf.eq47.gif    ・・・(4)

 原本这里的pbf.eq45.gifpbf.eq49.gif的平均要乘以pbf.eq50.gif的, 但是与SOR法的加速缓和系数同理,不乘的话收敛会更加迅速。

张力的稳定化

SPH法中,当临近粒子数很少的时候,会发生由负压力所产生的粒子过度集中的现象(称作particle clustering或particle stacking)。这里参照Monaghan[3] 的解决方案,把下式添加到式(4)里,将pbf.eq55.gif的部分改成pbf.eq56.gif

pbf.eq51.gif

这里的参数值pbf.eq52.gif,pbf.eq53.gif,pbf.eq54.gif。 

边界粒子

前面所说的particle stacking现象在边界附近也会发生。粒子与固体边界间的碰撞处理,由于固体内部不存在粒子,因此粒子会在边界附近聚集。一个比较好的解决方法便是在固体边界内设置位置固定的粒子(边界粒子)。通常,为了符合边界粒子的有效半径会使用复数层的粒子配置,这里我们参照Akinci的方法[4],使用提前计算好的幻想体积,只设置一层粒子来计算。把粒子沿边界设置好之后,只对边界粒子按照[4]中的下式计算幻想体积。

pbf.eq59.gif

此时,在计算PBF中的边界粒子质量时,替换成pbf.eq60.gif即可 。

参考文献 

[1] M. Macklin and M. Muller, "Position based fluids" ACM Trans. Graph., 32, pp.104:1-104:12, 2013.

[2] M. Muller, B. Heidelberger, M. Hennix and J. Ratcliff, "Position based dynamics", J. Vis. Comun. Image Represent., 18, pp.109-118, 2007.

[3] J. Monaghan, "SPH without a Tensile Instability", Journal of Computational Physics, 159, pp.290-311, 2000.

[4] N. Akinci, M. Ihmsen, G. Akinci, B. Solenthaler and M. Teschner, "Versatile rigid-fluid coupling for incompressible SPH", ACM Trans. Graph., ACM, 31, pp.62:1-62:8, 2012.


展开阅读全文

没有更多推荐了,返回首页