流体力学学习笔记3--Adaptively sampling particle fluids

1 简介：

这篇文章提出了基于粒子的流体仿真的自适应采样算法。介绍了基于局部几何特征尺寸的采样条件，这使计算资源能集中于几何复杂的区域，同时减少了流体内部和厚平面粒子数。而且根据视觉重要性改变采样密度能进一步的提高性能。另外，我们提出了一种新的基于粒子到表面距离的流体表面判定方法，每个粒子都储存到表面的距离信息并被适时更新。这种表面重建方法和现有方法相比有几点优势，包括粒子重采样下的稳定性和更适合于构建光滑平平面。这篇文章证明了自适应采样算法和基于距离的表面重建算法在内存和时间占用上的巨大改进。

2 物理框架：

仿真框架是基于SPH方法的，每个粒子有位置Xi，支持域Ri，和质量Mi。在仿真过程中，粒子的采样密度被动态更新，粒子被分成等级 l = 0，1 ,  2.....。等级0的粒子质量为m，支持域为r，其他等级粒子的质量和支持域分别为：粒子间的间距被初始化为ri/h，我们取h=2.5。如果，那么粒子pj被认为粒子pi的邻域粒子。邻域的搜索方法采用k-d树等范围查询算法。

3 特征自适应采样：

我们的目标是运用自适应粒子采样算法来动态的分配可利用的计算资源到我们感兴趣的区域。为了真实的重构流体表面，我们必须保证在细流和小水滴区域有足够的粒子密度。然后，对于流体内部和厚平表面区域的粒子，需要用来表面重建的粒子数会减少。下面我们会定义一个扩展局部特征尺寸来捕获粒子的这个特征。我们可以证明，虽然我们的方法是基于几何的，但是这并不会影响流动特性，特别是在感兴趣区域。

3.1 扩展局部特征尺寸

我们的几何采样条件是基于局部特征尺寸描述符的。封闭的区域表面M上的一个点y的局部特征尺寸被定义为y到中区域中轴线上最近点的距离。我们定义这个特征尺寸为这个局部特征尺寸可扩展到被区域包围的体积域v中，.即点x的扩展局部特征尺寸被定义为到点y的距离加上y的局部特征尺寸和的最小值。图2展示了一个二维区域的局部特征尺寸定义。

假定对于每个粒子pi我们都知道离它最近的在表面的y点，我们将高效的算出方程3.在4.2中我们将介绍怎么找到y点。我们通过计算流体面的中轴线粒子来估算扩展局部特征尺寸，这用计算表面附近粒子的扩张局部特征尺寸。运用快速步进法，内部粒子的局部特征尺寸也将很快得到。

近似中轴构造：

考虑2个位置分别为xi和xj的邻近粒子pi和pj，和各自对应的点yi和yj。如果满足下式，则两个粒子在中轴附近。

我们取区分角度为60度。第一个条件表示两个粒子应该在中轴的两侧，第二个条件将会裁剪掉一些分支，尤其是在表面附近。只考虑那些邻域中至少有50%的粒子在中轴线附近的粒子，中轴粒子将进一步减少。最后那些孤立的粒子也将被算作中轴粒子。点y的局部特征尺寸被定义为它到中轴粒子的最小距离，从点y的局部特征尺寸出发，我们将很快计算出粒子位置的扩展局部特征尺寸。

扩张局部特征尺寸的计算：

方程3的扩展局部特征尺寸可以通过快速行进算法来计算。初始化时，对于在表面附近的粒子，如，.这些粒子将按照扩展局部特征尺寸从小到大的顺序被插入到一个优先级队列中来。其他的粒子的扩展局部特征尺寸都被初始化为无穷大。接下来，每次从队列顶部取出粒子pi，算法将从外到里计算粒子的局部特征尺寸。pi对应的表面点yi‘用来计算pi的所有邻域粒子，例如：将被设置为如果比之前的值小。接着，pj将被加入到队列中，yj’被替换为yi‘。更新完所有邻域后，算法继续从队列中取出下一个粒子。

3.2自适应采样

如果粒子pi将被分成两个粒子pj和pk，如果，这个粒子将可能和其邻域粒子合并。为了防止震荡上下抽样，α应该比β小，我们这里取α=2，β=3.

分离：

粒子pi分离成两个粒子pj和pk，lj = lk = li-1.两个新粒子pj和pk分别被放置在距离pi的位置上。pj和pk的位置被精确选择以保证pj和pk在流体域里而且pj和pk和任何其他粒子的距离都不能超过d来避免大压力。最后一个条件的限制如图4c所示。所有和pi的距离在2d之内的粒子定义了半径为d，中心为xi的无效球区域。这些区域在图4c中以红弧表示。剩下的区域如果pj和pk能相对pi对称的放置在无效区域之外，这些区域则是有效区域。图4c中有效区域以蓝弧表示。通过在球坐标系上的球-球布尔运算，算法将高效的实现。注意到可能有多个有效位置，我们将从中随机选择一个。在极少数的情况下，会出现没有有效区域的情况，我们将迭代的减小d直到找到有效区域。我们将看到当d取一定值时，各个球区域相切，球体内区域都成为有效区域，红弧将消失。这也就意味着新粒子不能放置在距离存在粒子比粒子间距的二分之一还小的位置。最后，新粒子的速度将继承pi的速度。

合并：

一个粒子将被合并如果它和它的邻域粒子都符合合并条件，而且等级相同。合并产生的新的粒子pk，lk=li+1=lj+1，xk = (xi+xj)/2 （需满足xk在流体域里且没有粒子和pk的距离小于rk/(2h)(为了防止大压力））。如果这些条件不能满足的话，粒子将取消合并，并进入下一个时间步。新粒子的速度置为老粒子的平均速度。

未完待续。。。。。