CN103714575A - 一种sph与动态表面网格相结合的流体仿真方法 - Google Patents

Abstract

本发明涉及一种SPH与动态表面网格相结合的流体仿真方法。该方法使显式表面追踪用粒子系统采样的流体，并在仿真过程中自动处理自交和拓扑变化，维护三角形网格的质量；该方法提出了新的基于显式表面的自适应重采样方法。SPH仿真不易展现出很细很薄的特征，尤其是在粒子分辨率不足的情况下，甚至还会形成不自然的空洞。在传统SPH基础上如果想体现出更好的细节，可以增大粒子的分辨率，但是增大分辨率会相应地增加时间和空间上的计算代价，本发明基于显式表面的自适应重采样方法既可以模仿出非常精微的流体表面细节从而提升仿真的效果，同时又可以提高仿真的效率和速度。

Description

一种SPH与动态表面网格相结合的流体仿真方法

技术领域

本发明属于计算机图形与动画技术领域，涉及一种流体仿真数据表示方法以及流体运动和动画的高效仿真计算方法。

背景技术

在自然界以及日常生活中，流体是很常见的物质形态。在很多计算机图形学应用中，流体的仿真也占据了很重要的地位，例如在动画、游戏等娱乐领域，以及在科学可视化和医疗应用领域的应用等。现在计算机图形学上的流体仿真与动画方法有很多种，主流的有基于欧拉观点的有限差分方法（通称欧拉法）和基于拉格朗日观点的粒子法（无网格法，通称拉格朗日法）。两种方法各有优缺点，各有应用。粒子法在描述流体的细节特征（波动以及涡旋）方面具有天然的优势。

欧拉法（Euler法）最早由Stam在1999年引入图形领域（Jos Stam.Stable fluids[C].Proceedings of the26th annual conference on Computer graphics and interactive techniques.ACMPress/Addison-Wesley Publishing Co.,1999,121–128）。在这篇文章中作者描述了使用有限差分法进行流体仿真的算法流程，给出了在Euler框架下求解Navier-stokes方程的方法，并给出了平流和投影的数值解法。这篇文章提供了在图形学中进行基于物理的流体仿真的基本框架。

SPH方法（光滑粒子动力学法）是一种基于Lagrange观点的典型的无网格方法。2003年M¨uller将该方法引入计算机图形学中进行流体仿真（Matthias M¨uller,David Charypar,Markus Gross.Particle-based fluid simulation for interactive applications[C].Proceedings of the2003ACM SIGGRAPH/Eurographics symposium on Computer animation.EurographicsAssociation,2003,154–159）。SPH方法避免了基于网格方法中频繁的重网格化操作，能够增加程序的鲁棒性，降低计算复杂度。SPH方法在连续介质特别是流体的仿真领域中有广泛的应用。图形学中的流体仿真有多种方法，其中比较常用的是基于Euler观点的有限差分法以及基于Lagrange方法的无网格方法（SPH）。相比于Euler方法，SPH方法能更好地展现流体的细节，例如流体的飞溅效果。另外，通过调节粒子的规模可以很容易的使用SPH方法达到实时仿真而不使效果有很大损失。

动态显式表面的主要思想是在仿真的初始阶段给定一个满足闭合流形条件的初始显式表面（通常即为三角形网格），随着仿真的不断进行，每个时间步都得到一个空间中的速度场，由该速度场驱动显式表面的顶点进行运动并保持连接关系不变。之后得到了一个可能有相交以及退化三角形的网格，对此网格上的三角形进行操作以保证三角形的质量（由面积、边长比、角度等指标来衡量），同时修正已经相交的部分使网格仍然保持闭合流形性质，在此过程中应能处理由于运动引起的拓扑变化。由于动态显式表面的输入仅需要初始网格和速度场，与底层的仿真无关，只要是能提供速度场的数值方法都可以使用。又由于动态显式表面可以很好地处理形变以及拓扑变化，很适合应用在流体仿真中作为流体的表面。传统的流体仿真使用隐式表面，即通过仿真模型生成一个距离场，再从距离场上抽取出一个等值面作为流体的表面。显式表面与隐式表面相比有如下好处：

·显式表面比隐式表面更有利于保持住表面细节，这一点上要优于动态隐式表面，后者常常会因为平流中的数值耗散而抹掉细节；

·显式表面能够很好地表示几何上退化的特征，例如三维场景下的一维或二维特征。而隐式表面则受限于空间离散化的分辨率大小，往往很难表现这类特征；

·显式表面能够很好地保持表面上的参数化特征，例如纹理坐标等等。使用动态显式表面意味着可以省去每步重建表面的过程，不仅可以节省大量的内存空间（以及通常会节省计算时间），更重要的是可以天然地保持住表面上各点在时间序列上的演变关系，从而可以使参数化特征的传递有迹可循；

·对于流体来讲，引入显式表面可以有助于以一种更直接的方式来计算表面物理，例如由于表面张力引起的表面压差和毛细波等现象。此外，有了显式表面，可以借助离散微分几何上的大量知识来理解当前的仿真系统，例如可以便捷地计算表面积以及体积等等；

·随着仿真的进行，显式表面有时间维度上的连续性。由于在每一帧中隐式表面是作为当时时间步的后处理操作来进行的，即隐式表面的形状只取决于那一时间步生成隐式表面的等值面的形状。有时在不同帧之间由于粒子分布的不同等值面的变化可能不连续，这体现为视觉效果上闪动。而显式表面则是由上一时间步的表面运动而来，保留了上一时间步的状态，很大程度上可以避免这种闪动的发生。

在使用SPH进行仿真计算时，如果粒子系统在局部受到拉伸作用，会产生拉伸不稳定（Tensile Instability）现象。这种现象的产生是由SPH的计算模型所决定的，表现出来的结果就是粒子系统在比较细、薄的特征处的粒子运动会趋向于不稳定，可能彻底分散或是聚成小团。我们希望在粒子变得过于稀疏之前可以对粒子系统进行自适应的重采样，从而使仿真可以展现出一些细薄特征。

发明内容

由上述内容可知，基于SPH思想的仿真方法不易表示和展现出很细很薄的特征，尤其是在粒子分辨率不足的情况下。在传统SPH方法基础上如果想体现出更好的细节，可以增大粒子的分辨率，但是增大分辨率会相应地增加时间和空间上的计算代价。如果能够实现自适应的重采样则可以把有限的计算资源根据细节的粗细进行合理地分配。本发明针对该问题，提出了一种显式表面与SPH流体相结合的流体仿真方法，可以使显式表面追踪用粒子系统采样的流体，并在仿真过程中自动处理自交和拓扑变化，维护三角形网格的质量。

具体来说，本发明采用的技术方案如下：

一种SPH与动态表面网格相结合的流体仿真方法，其步骤包括：

1）根据脚本在3D空间中对初始流体物质进行采样以生成粒子系统，计算由粒子定义的隐式函数

并由等值面

生成初始三角形网格表示流体的外表面形态，后续的流体运动仿真过程的计算将分为粒子管线和表面管线两部分组成；

2）对于粒子管线，首先进行初始粒子设置以表示流体所处的初始状态，然后根据SPH方法计算粒子所受压强力、粘滞力和表面张力，并移动粒子，其中表面张力由显式表面上的信息所提供；然后对移动后的粒子系统重新计算隐函数

并得到其等值面

3）对于表面管线，首先初始化流体的显式表面，该显式表面为具有闭合流形性质的三角形网格，然后从粒子系统得到显式表面上每一个顶点的插值速度，并相应地移动顶点至新位置；

4）将步骤3）所得的运动后的网格表面投影到步骤2）重新计算后的等值面上，在此过程中保证显式表面仍具有闭合流形性质且无明显的相交，同时处理有可能产生的拓扑变化事件，然后根据表面及两级距离场对流体进行重采样，从而完成当前时间步之内的流体仿真的所有计算，之后进入下一时间步的迭代，重复执行上述步骤2）-4）；

5）当仿真时间满足一定条件或者仿真状态满足一定条件时，整个流体仿真过程结束。

进一步的，步骤2）中通过粒子系统计算

并得到其等值面

该过程首先定义在三维空间的标量函数

对于空间中一点x，计算

其中i是粒子的下标，W是平滑函数，r是粒子i的影响域半径。在一个规则格子上对

进行采样，在进行采样之后，将粒子没有覆盖到的格点上的

都设成一个负值。这样可保证在有粒子覆盖的地方

值是正的，离粒子系统较远的地方则呈负值，等值面则可以用来定义粒子系统的表面。

进一步的，步骤3）中从粒子系统得到显式表面上每一个顶点的插值速度，该过程构造在空间任意一点x处的插值速度公式为

其中W_i为第i个采样顶点的平滑函数，V_i为其体积。对于在表面附近粒子采样比较稀少存在采样不足的现象，采用简化的插值速度，对W进行归一化，即

减少粒子稀少对插值的影响。

进一步的，步骤4）采用的投影算法包括如下步骤：

a）对于显式表面上的顶点v，其以插值速度运动之后空间位置为x，法向为n_v，在标量函数场

上计算

和其中

r是粒子系统的初始设定半径，通常是k_r倍的粒子采样间距，如果

是精确的到某个闭合曲面的距离场，则使用

可得到更快的收敛速度；

b）如果

则在从x到x+εn_v的线段上用二分法迭代查找

的0值点。

c)如果

使得表明此时顶点v的拓扑发生变化，此时如果

则将v移至x+εn_v，否则不移动v的位置。另外还要设置一个拓扑变化标志来提醒后续步骤的算法需要处理拓扑变化。

进一步的，步骤4）在处理显式表面相交和拓扑事件时，不对显式表面地整个表面进行重建，而只对发生拓扑变化或是自交的区域进行局部重建，具体处理步骤包括：

a）在空间中建立一个正则格子；

b）在表面附近的格点上计算有向距离场；

c）识别发生拓扑事件以及自交的区域，标记为重网格化区域；

d）对表及区域重网格化，并修补与外部区域网格的接缝。

进一步的，步骤4）中识别发生拓扑事件以及自交的区域，具体步骤包括：

a)依次检查一个单元是否含有复杂边、复杂面以及完全包含在单元内部的闭合曲面分支，判断是否满足复杂单元的条件；

b)距离隐式表面

超过一定距离的单元，称为深单元。从深单元开始向外扩展以包含周围所有的复杂单元，直到该区域的边界单元都是简单的，这样既能保证与拓扑变化相关的复杂单元都被标记为重网格化区域，同时尽可能保住了高分辨率的几何特征；

c)进行扩展的复杂单元的标记，即对Marching Cubes的15个模板中有4个是在一个格子面上允许有两条三角形边的情况标记为重网格化区域。

进一步的，步骤4）中重网格化标记区域的步骤是：

a）处理标记区域的边界，使内部和外部的三角形片段沿区域边界分隔开；

b）删除标记区域内部三角形，使用Marching Cubes模板在区域内部重新建立三角形；

c）处理边界，将内部外部三角形片段缝合。

进一步的，步骤4）进行重采样的步骤包括：

a)检测采样不足处的粒子，标记重采样区域；

b)在标记区域内进行粒子的快速泊松盘重采样并增加采样点；

c)重新分配重采样区域的粒子质量和动量；

d)合并符合质量条件、密度条件以及质量约束这三条合并条件的粒子。

进一步的，步骤4）中粒子的快速泊松盘重采样方法，其步骤包括：

a)在每个格子单元内部，随机生成t个探测位置。对每生成一个探测点，检测其是否离现有的粒子以及表面足够远，即探测点的位置x必须满足

且|(x-xi)|>k_pd₀且|x-v_j|>k_vd₀，其中即x必须在流体内部，i和j分别是粒子与表面顶点下标，k_p和k_v是常数，d₀是仿真开始时粒子的初始间距；

b)如探测点满足以上要求，则将探测点处插入一个新的粒子，并将粒子映射到格子上，参数k_p和k_v限制可以插入粒子的区域，即不仅在流体内部，还要离粒子和表面一定的距离；

c)对新产生的粒子进行一步松弛操作以优化新粒子的位置。对于每个新产生的粒子p_i，其以r为半径的邻域内的粒子下标集合为NP_i，邻域内的表面顶点下标集合为NV_i，计算距离p_i与其周围顶点或粒子的最小距离 $d_{\max \min }=\min \left(\underset{j\∈N{P}_i}{\min }\right|x_i-x_j|,\underset{j\∈{\mathrm{NV}}_i}{\min }|x_i-x_j\left|\right)$ 之后，在p_i的邻域范围内随机撒若干探测点。对于每个探测点如果其计算出的d_maxmin大于原来的值，则用

的位置替换p_i。最后，p_i将被移动到其邻域内与周围顶点和粒子距离近似最大的位置。

进一步的，步骤4）中重采样区域内的粒子的质量和动量的重新分配，其步骤包括：

a)将重采样区域进行聚类，并按照聚类结果进行分配。聚类的方法是设置一个队列并将一个重采样单元放入队尾。之后进行迭代循环，将队首单元弹出，加入集合S，并检查这个单元的邻接单元，如果也为重采样单元则加入队尾，直到队列为空，此时S中的单元聚为一类。如所有重采样单元都已被检查过，则聚类完成，否则继续将一个未检查的单元放入队尾进行以上操作;

b)在质量的分配上，采用的方法是平分一类内部的粒子质量，设粒子数为n_c，原粒子的下标集为P_c，粒子p_i的质量为m_i，则重采样后粒子上的采样质量为

计算完m′之后将其与一个质量阈值m_min比较，如有m′<m_min则舍弃该聚类内部的所有重采样点并停止重采样过程，维持聚类粒子状态不变。

c)在动量的分配上，在新插入的粒子处使用插值估算出速度

其中v_j是原粒子的速度。之后计算在一个类内部原粒子的合动量

然后通过调整粒子的速度对合动量进行归一化。

本发明提出的显式表面与SPH流体相结合的流体仿真方法，可以使显式表面追踪用粒子系统采样的流体，并在仿真过程中自动处理自交和拓扑变化，维护三角形网格的质量。同时基于显式表面的自适应重采样方法提高了仿真的效率和速度。与现有技术相比，本发明的优点如下：

1）可以模仿非常精微的流体表面细节。流体表面的网格质量比传统的表面重建法要好，主要是细节处效果的提升。而且在时间的连续性也要好得多，不会带来表面重建法中表面上细碎的噪声特征。

2）自适应重采样可以给仿真带来显著的效果提升。这一算法可以使表面变得更加连续，去除表面上因粒子局部稀疏采样不足而造成的空洞。能够自动识别流体表面内部粒子采样比较稀少的区域并以合理的密度进行重采样，在采样之后可以对采样的物理量进行合理的重新分配。应用该算法可以使粒子法流体仿真也展现出一些薄片或是细线的细节特征。

附图说明

图1是本发明的流体仿真方法的步骤流程图。

图2是表面附近采样不足的示意图。

图3是投影过程的示意图。

图4是复杂单元识别并进行融合的示意图。

图5是复杂单元识别并进行分离的示意图。

图6是处理边界上的格子边示意图。

图7是处理边界上的格子面示意图。

图8是重新生成三角形并进行缝合的示意图。

图9是标记重采样区域的示意图。

图10是可行的采样区域示意图。

图11是流体表面仿真效果及其局部放大图。

图12是另一个流体表面仿真效果图。

具体实施方式

下面通过具体实施例和附图，对本发明做详细的说明。

本发明提出了一种显式表面与SPH流体相结合的流体仿真方法。该方法可以使显式表面追踪用粒子系统采样的流体，并在仿真过程中自动处理自交和拓扑变化，维护三角形网格的质量。在该方法中提出了一种新的基于显式表面的自适应重采样方法。SPH仿真不易展现出很细很薄的特征，尤其是在粒子分辨率不足的情况下。在传统SPH基础上如果想体现出更好的细节，可以增大粒子的分辨率，但是增大分辨率会相应地增加时间和空间上的计算代价。如果能够实现自适应的重采样则可以把有限的计算资源根据细节的粗细进行合理地分配。

本方法是基于SPH方法的流体仿真框架，在算法中使用了动态显式表面和自适应重采样。算法初始化生成粒子系统和显式表面，之后两者分别随着时间步的推进进行演化。算法的流程如图1所示。在图1中，各个阶段流程如下：

a)根据脚本在3D初始流体物质进行采样以生成粒子系统。计算由粒子定义的隐式函数

并由其定义的等值面生成初始三角形网格。图1中实心圆形表示流体粒子，空心圆形表示重采样之后新加入系统中的粒子，实线边界表示流体的显式表面形式（三角形网格表示）而虚线边界表示隐式表面（如图1d）。接下来的过程分为两个管线，分别为粒子管线和表面管线，其中粒子管线包括状态b),c),d)，表面管线包括状态e),f)。之后两条管线的结果相结合称为状态g)，进行重采样后进入状态h)；

b)初始粒子设置，表示流体所处的初始状态；

c)根据SPH方法计算粒子所受压强力、粘滞力、表面张力，并移动粒子。其中表面张力是由显式表面上的信息所提供的；

d)对移动后的粒子系统重新计算

图中虚线表示等值面

e)初始显式表面，该显式表面为具有闭合流形性质的三角形网格；

f)然后从粒子系统得到显式表面上每一个顶点的插值速度，并相应地移动顶点至新位置；

g)将运动后的网格表面投影到d)中的等值面上。在此过程中要做一系列处理，保证显式表面仍具有闭合流形性质且无明显的相交，同时能够处理分离和融合等拓扑事件；

h)根据表面及两级距离场对流体进行重采样，图中空心粒子为重采样所得到的粒子。在此状态之后算法进入下一时间步的迭代，重新执行c)-h)的所有步骤；

i)当仿真时间满足一定条件或者仿真状态满足一定条件时，整个流体仿真过程结束

下面具体说明一些关键步骤的实施细节：

1.显式表面的运动（步骤f）

1.1显式表面的定义与运动原理

在离散微分几何上曲面的参数化方法有很多种形式，本发明所提到的显式表面都是离散化的三角形网格。本发明的三角形网格使用半边结构进行存储。符号定义如下：

·v_i――三角形网格的第i个顶点

·e_i――三角形网格的第i条边

·n_vi――第i个顶点的法向

在半边结构中，边（两个顶点编号构成）的信息和三角形（三个顶点编号）的信息只定义了网格的拓扑连接信息，其空间位置信息则是由顶点列表来定义的，顶点列表中的每个分量是一个顶点在三维空间中的坐标。故只要通过速度来更新每个顶点的坐标则就能使表面运动起来。

对于顶点v_i，在时间步k结束后，设其位置为

速度为

则其在时间步k+1结束后的位置为： $x_i^{k+1}=x_i^k+v_i^k\&CenterDot;\mathrm{\&Delta;t}---\left(1\right),$

在式(1)的数值积分中使用了前向Euler迭代法，其中Δt是时间步的大小。

1.2在SPH流体中驱动显式表面的运动

根据式(1)，只要计算出三角形网格上每个顶点v上的速度，即可使顶点进行运动。根据SPH的插值方法，构造在空间任意一点x处的插值速度为

其中W_i为第i个采样顶点的平滑函数（又称为核函数或者权重函数），V_i为其体积。由于在表面附近粒子采样比较稀少（通常少于在液体内部的），如图2所示，其中虚线圆圈是表示液体表面某处的平滑区域范围，存在采样不足的现象，实线圆圈是表示流体内部某处的平滑区域，故有

所以使用简化的插值速度，对W进行归一化，即减少粒子稀少对插值的影响。

1.3投影（步骤g）

当使用SPH插值速度移动表面一步之后，数值误差可能使表面与粒子系统趋于分离。在经过若干时间步之后，如果这种分离得不到纠正，那么数值误差的积累将使表面与粒子系统愈来愈分离。为了确保表面能够追踪粒子系统，我们使用粒子系统定义一个标量函数

并将运动后的显式表面投影到隐式表面

上。图3是投影过程示意图，其中黑色弧线表示的是粒子系统定义的标量函数折线代表的是使用插值速度运动后的显式表面。

对标量函数

习惯上一般当x在粒子系统所影响的区域内部时有

在外部时则该区域的表面即为等值面

本发明中因为实现时使用了与习惯标法相反的标法，故后文提到

均以

表示外部，以表示内部。下面重点介绍投影算法。对于显式表面上的顶点v，其空间位置为x（该位置为以插值速度运动之后），法向为n_v。在标量函数场

上计算

和

其中

r是粒子系统的初始设定半径，通常是k_r倍的粒子采样间距，实现时可以设置k_r=2.5。接下来，如果

则在从x到x+εn_v的线段上用二分法迭代查找

的0值点。通常最多使用4次迭代并以误差0.001为终止条件，可以取得较好的效率及精度。另外，如果

是精确的到某个闭合曲面的距离场，则使用

可得到更快的收敛速度。这种投影方法有效的一个前提条件就是在每个时间步中，使用插值速度进行运动后显式表面不会偏离粒子系统太多，这个偏移的门限值是εnv，一旦显式表面与粒子系统定义的等值面

在其法向上偏离超过ε则投影算法无法正确地将显式表面投影到

等值面。

在三角形网格运动过程中，如无拓扑变化，偏离超过门限值得情况是很少见的，但是一旦拓扑变化发生则会引起

使得

故

可以作为在v附近可能有拓扑变化产生的一个指示。以拓扑融合为例，当粒子系统的两个不同分量进行融合时，函数的拓扑发生变化，显式表面必须也相应进行更改。此时在

进行融合的区域附近，将有部分三角形网格的顶点出现此时如果

则将v移至x+εn_v，否则不移动v的位置。另外还要设置一个标志来提醒后续步骤的算法需要处理拓扑变化。

2.显式表面相交和拓扑事件的处理（步骤g中的处理要点）

针对上面所述显式表面的运动计算方法，表面不可避免地会发生自交、退化等情况，既影响计算的鲁棒性，也影响渲染模型的美观。另外，如果不能动态处理表面的拓扑，当粒子系统拓扑发生变化时，表面就不能与粒子系统始终保持一致。所以处理拓扑事件和网格自交是使用动态显式表面时的一个重要问题。

显式表面在运动中大部分情况下并不需要对网格进行重网格化，但在少部分情况下，网格需要进行拓扑变化（网格的融合和分离等）或是网格发生自交，在这种情况下需要修改这些区域的三角形片，恢复网格的性质。Glimm等人曾经尝试过在发生类似状况时就全局重建表面（James Glimm,John W Grove,XL Li,DC Tan.Robust computational algorithms for dynamicinterface tracking in three dimensions[J].SIAM Journal on Scientific Computing.2000,21(6):2240–2256)，这个方法避免了每个时间步都使用Marching Cubes建立新的表面(William ELorensen,Harvey E Cline.Marching cubes:A high resolution3D surface constructionalgorithm[C].ACM Siggraph Computer Graphics.ACM,1987,vol.21,163–169)一定程度上保持了表面在时间上的连续性，但是在需要重建的时间步，表面的时间连续性还是会间断。尤其是拓扑变化通常需要若干个连续时间步来完成，故在这一时间段内表面也有可能会出现由于每步重建带来的抖动。

本发明的创新点之一就是不针对显式表面地整个表面进行重建，而只对发生拓扑变化或是自交的区域进行局部重建。遵循Marching Cubes方法的框架，只要能够检测到被拓扑变化与自交的三角形片段所覆盖的格子单元并使用Marching Cubes的模板在这些格子里面重建出三角形片段，并和周围的网格缝合起来，即可得到新的网格。

2.1算法流程

处理表面拓扑与相交的算法和计算模型（如SPH）无关，任何计算模型只要能够为表面上的每个顶点提供速度即可驱动表面运动。

本拓扑处理算法的步骤如下：

1）在空间中建立一个正则格子（grid）；

2）在表面附近的格点上计算有向距离场；

3）识别发生拓扑事件以及自交的区域，标记为重网格化区域；

4）对表及区域重网格化，并修补与外部区域网格的接缝。

在经过以上步骤之后，三角形网格已经被修正为一个符合闭合流形性质且无明显自交的网格，且在拓扑上与标量函数

大致相同。接下来将详细介绍各个步骤。

●建立正则格子

格子（grid）是一种常用的数据结构，将空间划分为划分为k×m×n个立方体单元（cell）。使用格子结构可以对网格进行空间分割（例如将网格的顶点位置存放在格子单元中），能够加快对网格元素的搜索速度。这一过程称为空间Hash化。除了空间Hash化之外，正则格子还有很多作用，例如标量函数的值就是存储在格子上的。在有限差分法仿真中格子的格点是物理量的采样点。另外，Marching Cubes方法建立三角形网格的方法就是使用格子在每个立方体单元（cube）中根据标量函数建立三角形。

●计算有向距离场

在格子上为离显式表面较近的格点计算到显式表面的有向距离。设空间中任一点x，其到闭合曲面S的有向距离的定义为：

$d_{\mathrm{sign}}\left(x\right)=s\underset{x^{\&prime;}\&Element;S}{\min }|x-x^{\&prime;}|$

在格点上计算对显式表面的d_sign时，只要寻找离格点最近的三角形并计算其距离就行了。还有另外一个问题是如何判断一个格点x是在显式表面的外面还是里面。首先初始化一个相交次数计数器为0，然后从x向某一方向发射一条射线（通常是沿坐标轴的正向或负向），计算射线与显式表面的相交次数。如果射线在一次相交中射线进入显式表面内部，计数器-1；如果射线离开显式表面，则计数器+1。如果最终计数器为+1，那说明点x在显式表面内部，如果计数器为0，则说明点x在表面外部。如果计数器出现其它值，说明该格点处于明显的网格相交之中。

●识别并标记拓扑变化与自交区域

能够识别拓扑变化所发生的区域是处理拓扑变化的先决条件。在Brochu的文章中（TysonBrochu,Robert Bridson.Robust topological operations for dynamic explicit surfaces[J].SIAMJournal on Scientific Computing.2009,31(4):2472–2493），作者使用距离来判断合并事件，即当两个曲面分支距离小于一定阈值时则将两个分支合并。而曲面分支的距离则是通过由距离最近的网格边―边对的距离来定义的。这个方法有几个问题：

1）曲面的距离判断十分耗时，尤其是在曲面比较复杂的情况下，距离相近的边―边对可能特别多；

2）通过距离判断只能筛查出距离较近的两条边，而处理也只能在这两条边上进行，即在合并的双方分别删除这两条边，在曲面上留下两个空洞，然后使用三角形模板将空洞联通。曲面合并的时候往往在一个区域内会有诸多距离近的边―边对，这样的联通方式很可能在合并区域内形成许多联通而使整个联通分支变成一个高度复联通的分支，这既影响视觉效果也影响程序的鲁棒性。

本发明采用局部的Marching Cubes方法识别拓扑变化的方法：

首先定义一个概念―复杂单元。复杂单元是正则格子中的一个单元，这个单元中所包含的显式表面网格片段与使用Marching Cubes模板生成的三角形片段形状显著不同。所谓显著不同，是指该单元含有复杂边、复杂面或为复杂单元。以上提到的边、面均为格子单元立方体的边和面，以后提到的时候会称之为格子边、格子面，与三角形网格的边和面相区分。

关于复杂边、复杂面和复杂单元的定义：

复杂边：如果一条格子边与网格相交超过一次则认为是复杂边。Marching Cubes最多只能在一条格子边上识别一个交点；

复杂面：如果一个格子面的四条边中有复杂边，或者网格穿过格子面的交线为一条闭合曲线（折线），则认为是复杂面。Marching Cubes不能构建出只穿越格子面而不于其四条边相交的三角形。

复杂单元：含有复杂边或者复杂面的单元即为复杂单元。另外，如果在单元内部有一个闭合曲面分支，则该单元也为复杂单元，因为该分支对分辨率的要求已经高于Marching Cubes的分辨率；

使用复杂单元可以帮助确定拓扑变化发生的位置。上文为显式表面附近的格子单元的格点计算了到显式表面的有向距离，而离表面较远处的格子则没有考虑。针对这些计算过距离场的单元，我们按照上面提到的方法分别判断是不是复杂单元。判断的方法则是依次检查一个单元是否含有复杂边、复杂面以及完全包含在单元内部的闭合曲面分支。复杂单元将被标记，并在后面的步骤中对内部的三角形进行重网格化。

图4是复杂单元识别并进行融合的示意图，图5是复杂单元识别并进行分离的示意图。在图4和图5左侧图中，白色格子单元即为识别出的复杂单元。从图中可以看出，该方法每次可处理的拓扑变化区域大小与处理后该处网格的精细程度与使用的格子分辨率有关。

使用复杂单元作为判断拓扑变化标准的两条准则：

深单元的标记：并非所有的复杂单元处理后都会形成拓扑变化。在流体表面处可能有分辨率大于格子的细节，而很多情况下这些细节并不会引起拓扑变化，如果这种情况也作为拓扑变化处理则可能使显式表面在演进中丢失很多几何细节。故我们在标记拓扑变化的时候只标记距离隐式表面

超过一定距离的单元，这些单元则叫做深单元。之后从深单元开始向外扩展以包含周围所有的复杂单元，直到该区域的边界单元都是简单的，这样既能保证与拓扑变化相关的复杂单元都被标记上，同时尽可能保住了高分辨率的几何特征；

扩展的复杂单元的标记：以上识别复杂单元的方案是根据Marching Cubes模板的性质制定的。Marching Cubes模板中有的情况在局部Marching Cubes中是有歧义的。Marching Cubes的15个模板中有4个是允许在用一个格子面上允许有两条三角形边的，这样的模板对是否在这个格子面上的两条边之间是否生成三角形有歧义，如果避免不了这种情况则有可能在重网格化后的表面上留下空洞。故我们在实现过程中将符合这四个模板的格子单元也标记为复杂，即扩展的复杂单元。

●重网格化标记区域

重网格化标记区域需要将标记区域内部的网格删除，并在区域内部重新建立三角形片段。重网格化的步骤是：

1）处理标记区域的边界，使内部和外部的三角形片段沿区域边界分隔开；

2）删除标记区域内部三角形，使用Marching Cubes模板在区域内部重新建立三角形；

3）处理边界，将内部外部三角形片段缝合。

处理区域边界的方法是：

1）首先对于每个与标记区域边界面的格子边相交的三角形，在交点处插入一个顶点（三角形网格顶点），并将该三角形细分。这一过程见图6，其中黑色边框为标记区域边界的一个格子面，灰色三角形为显式表面的一个片段。

2）之后对于显式表面上每一条与格子面相交的边，在交点处插入一个顶点，并将边细分。这一过程见图7，其中浅灰色三角形为显式表面在标记区域外部的三角形片段，深灰色三角形为显式表面在标记区域内部的三角形片段，浅灰色三角形片段和深灰色片段之间的线为沿着区域边界格子面的内、外部分界线。

在对边界上所有格子面做完这一步之后，可保证标记区域内部与外部的三角形完全分隔开，而图7中的分界线在整个区域范围内将成为闭合曲线。此时检查整个显式表面，如果某个三角形全部包含在标记区域内部，则将三角形删除，如图8(b)所示。在实施例中判断三角形是否完全在标记区域内的方法是：对于三角形的三个顶点，如果都在标记区域内或属于闭合分界线上的顶点，则判断三角形的中心，如果在标记区域内则该三角形可判定为全部在标记区域内部。内部三角形删除之后，显式表面上就会留有以闭合分界线为边界的空洞。用Marching Cubes模板在标记区域内生成三角形，如图8(c)所示。由于Marching Cubes生成的三角形在格子面上都是直线，故在生成三角形之后将标记区域内部的三角形和外部的三角形缝合起来。具体做法是：将内部与表面相邻的三角形删除，该三角形的三个顶点有两个在边界格子面的格子边上（顶点v₁v₂），另一个在标记区域内部（v₀）。将v₀与v₁、v₂之间的折线中的每一线段相连，如图8(d)所示。最后为了不使连接处三角形分布过于复杂，我们将在边界格子面上的折线边坍缩掉，形成比较整洁的边界，如图8(e)所示。在这一过程中，在执行坍缩之前首先会检查是否会使格子边上的点过早相连，即两点之间尚有折线段未坍缩之前便因为其它边的坍缩而相连。如果是则放弃此次操作。

●显式表面网格质量的保持与鲁棒性

为了保证显式表面运动和相交检测的鲁棒性，表面上的三角形需要保持良好的形状。我们自适应地细分过大的三角形、合并较小的三角形与接近退化的三角形。在实施例中，进行以下操作：边细分，边坍缩，边翻转，顶点偏移。

边细分：当显式表面上有边大于平均边长的1.5倍时，就将这条边细分。具体的方法是在边的中点上插入一个新顶点，并把与这条边相邻接的两个三角形细分为四个三角形。

边坍缩：当显式表面上有边小于平均边长的0.5倍时，就坍缩这条边。方法是将这条边的两个顶点合并，同时将邻接的边合并。

边翻转：当三角形过于狭长的时候我们执行边翻转。检查每条边是否比其邻接三角形的两个顶点的距离长，如果长则翻转。在实施例中，禁止一切引起相交与引起退化三角形的翻转操作。因为边翻转可能在相邻边上引起反复，对整个显式表面遍历5次进行边翻转操作。

顶点偏移：在标记复杂单元和自交区域与重网格化过程中需要执行大量的格子边与显式表面的求交操作，求交的鲁棒性极其重要。某些退化的情况可能引起求交结果的错误，比如格子边正好穿过显式表面某一三角形的边或点，故应该避免这些情况的发生。解决方法是检测每个网格顶点的坐标，如果与某一格子边的坐标很相近的话就随机移动一个较小的距离，使之偏离格子边一点。

3.等值面生成（步骤d）

为了实现显式表面追踪粒子系统，本发明将显式表面投影到由粒子定义的标量函数的等值面（隐式表面）上。在进行这一操作时，等值面的形状就直接决定了投影后显式表面网格的形状，故如何从粒子生成一个较好的隐式表面，即一个较好的标量函数，对最终显式表面的质量有很大影响。

首先定义在三维空间的标量函数

对于空间中一点x，

其中i是粒子的下标，W是平滑函数，可以使SPH风格核函数，也可以是其它核函数（例如高斯函数），r是粒子i的影响域半径，例如

一般情况下我们在一个规则格子上对进行采样，在进行采样之后，将粒子没有覆盖到的格点上的都设成一个负值。这样可保证在有粒子覆盖的地方

值是正的，离粒子系统较远的地方则呈负值，而

等值面则可以用来定义粒子系统的表面。

4.自适应局部重采样（步骤h）

基于Lagrange观点的数值方法与基于Euler观点的数值方法的一个区别就是前者的采样点会随着物体的运动而运动，而后者的采样点则始终保持自己的空间位置。传统SPH的粒子（采样点）的位置随着仿真对象的运动而不停改变，粒子之间的相对位置也不断地变化。在某些体现仿真对象细节特征的地方，如薄片、尖角、细线等，使用SPH方法进行计算可能会出现粒子采样不足的情况。粒子采样不足会导致很多问题：

·粒子的邻居数显著减少，插值精度降低，破坏场的连续性；

·粒子密度减小，导致产生负压强；

·根据粒子建立或者追踪产生的表面会产生很多空洞，影响视觉效果。

本发明通过对粒子采样不足的区域进行自适应的局部重采样，使采样点的分布密集程度达到初始化水平，这样可以改善细节特征处的动力学特性与视觉效果。在增加采样点后，为了保持系统的动量与质量守恒，将系统的质量和动量重新分配。另外，增加采样会使粒子数量增加，从而增加计算代价，所以在适当的条件下将质量较小的粒子合并。该过程对应图1中的步骤g)→h)。自适应重采样的算法流程是：

1.检测采样不足处的粒子，标记重采样区域；

2.在标记区域内增加采样点；

3.重新分配重采样区域的粒子质量和动量；

4.合并符合合并条件的粒子。

接下来详细介绍各个步骤的方法：

4.1标记重采样区域

首先将整个空间离散化，建立正则格子。在这里使用的格子与前面步骤所述的格子相重叠，具有相同的空间位置与分辨率。该格子既用于标示重采样区域，也用于加速显式表面顶点和粒子的邻域搜索。

粒子p_i的密度不仅与其邻域N_i中粒子的质量有关系，而且与其周围粒子的分布有关。如果N_i中粒子过少，则p_i的密度通常会较小。通过检查粒子的密度来搜寻那些可能处于采样不足区域的粒子，下面将这些粒子简称为稀疏粒子。设置一个阈值k_dmin，当p_i的密度ρ_i<k_dmin·ρ₀时则将p_i标记为稀疏粒子，其中ρ₀为初始密度。

标记过稀疏粒子后，将粒子p_i半径r_i内覆盖的格子区域均标记为重采样区域。为了减少重采样时的计算量，标记重采样区域时会检测每个在粒子覆盖范围内的单元的八个顶点的

值，当对于一个单元的八个顶点存在一个顶点上的

时，该单元则被标记为重采样单元。这一过程见图9，该图表示的是在一个流体薄片细节附近进行标记重采样区域的过程。图中黑色粒子为稀疏粒子，灰色粒子为分布合理（不稀疏）的粒子，灰色曲线所示的边界为流体边界，阴影区域为重采样区域。左图中表示的重采样区域是粒子p_i所覆盖的区域，r为p_i的半径。图9右图中的重采样区域是由稀疏粒子最终确定的重采样区域。

4.2重采样

使用重采样算法插入的新采样点需要满足两个要求，一是必须保证插入的采样点必须在流体内部；二是采样点之间，采样点和表面之间必须满足合理的间距，既令采样点的分布能够充满整个流体，又能保证计算的稳定性和表面的稳定性。本发明设计并提出面向SPH粒子的快速泊松盘重采样方法。

在确定了需要进行重采样的区域之后，采样都在重采样区域内部进行。在每个格子单元内部，随机生成t个探测位置。对每生成一个探测点，检测其是否离现有的粒子以及表面足够远，即探测点的位置x必须满足

且|(x-x_i)|>k_pd₀且|x-v_j|>k_vd₀，其中

是前文所提到的标量函数，即x必须在流体内部。i和j分别是粒子与表面顶点下标，k_p和k_v是常数，d₀是仿真开始时粒子的初始间距。如探测点满足以上要求，则将探测点处插入一个新的粒子，并将粒子映射到格子上。参数k_p和k_v限制了可以插入粒子的区域，即不仅在流体内部，还要离粒子和表面一定的距离，如图10。图中被实线边界包围的阴影区域是被已有的粒子和表面覆盖的区域，落在该区域内部的探测点不能生成采样点，而重采样单元内的白色区域范围内可以进行采样。

在重采样之后，对新产生的粒子进行一步松弛操作以优化新粒子的位置。对于每个新产生的粒子p_i，其以r为半径的邻域内的粒子下标集合为NP_i，邻域内的表面顶点下标集合为NV_i，计算距离p_i与其周围顶点或粒子的最小距离 $d_{\max \min }=\min \left(\underset{j\&Element;N{P}_i}{\min }\right|x_i-x_j|,\underset{j\&Element;{\mathrm{NV}}_i}{\min }|x_i-x_j\left|\right)$ 之后，在p_i的邻域范围内随机撒若干探测点。对于每个探测点

如果其计算出的d_maxmin大于原来的值，则用

的位置替换p_i。最后，p_i将被移动到其邻域内与周围顶点和粒子距离近似最大的位置。

4.3质量和动量的重新分配

若想使经重采样新插入的粒子参与到计算中，必须对其赋予相应的物理量。由于重采样是在一个区域之内进行，我们希望重采样后的物理量能够反映区域的特征，故将本区域内的所有粒子，包括原有的粒子与新产生的粒子，全部进行重新分配。主要考虑质量和动量的重新分配。

聚类：在每个时间步里发生重采样的区域可能在空间上并不一定相连，显然不能将多个不连通区域中的粒子放在一起进行质量和动量的重新分配。故在重新分配质量和动量之前将重采样区域进行聚类，并按照聚类结果进行分配。聚类的方法是设置一个队列并将一个重采样单元放入队尾。之后进行迭代循环，将队首单元弹出，加入集合S，并检查这个单元的邻接单元，如果也为重采样单元则加入队尾，直到队列为空，此时S中的单元聚为一类。如所有重采样单元都已被检查过，则聚类完成，否则继续将一个未检查的单元放入队尾进行以上操作。

质量和动量的重新分配：重采样区域的一个聚类标记的是仿真域中较独立一块采样较为稀疏的区域，重采样实际上是加密了该区域的采样点，故质量和密度的重新分配应在一个类内部进行。在质量的分配上，采用的方法是平分一类内部的粒子质量，设粒子数为n_c，原粒子的下标集为P_c，粒子p_i的质量为m_i，则重采样后粒子上的采样质量为计算完m′之后将其与一个质量阈值m_min比较，如有m′<m_min则舍弃该聚类内部的所有重采样点并停止重采样过程，维持聚类粒子状态不变。这样做是为了防止粒子的无限分割致使流体无法分离，以及产生质量极小的粒子导致计算不稳定。

除了质量之外还有动量也需要重新分配。先在新插入的粒子处使用SPH插值估算出速度

其中v_j是原粒子的速度。之后计算在一个类内部原粒子的合动量

然后通过调整粒子的速度对合动量进行归一化。设

为重采样之后的下标集，对于粒子p_i，调整之后的速度是其中的除法与乘法运算均对向量的每个分量分别进行，而不是向量点乘操作。

4.4合并

当一个粒子p_i满足如下条件时：

就将粒子p_i和p_j合并。该式中ρ₀和m₀是初始密度和初始质量，k_mmin，k_dmax是常数。第一条是质量条件，当粒子质量比较小的时候才可以触发合并。第二条是密度条件，必须在密度较大的地方可以合并，这一条是避免在细节处粒子反复地合并与重采样，因为细节处通常有退化情形，密度通常低于流体内部。第三条是质量约束，如果p_i的邻域之内有粒子质量也较小则可以进行合并。设置质量约束是为了控制合并后粒子的质量的分布，防止产生质量很大的粒子。

图11和图12即为采用本发明方法的流体仿真效果以及局部微观放大的效果图，证明了本发明方法的有效性。

以上实施例仅用以说明本发明的技术方案而非对其进行限制，本领域的普通技术人员可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明的精神和范围，本发明的保护范围应以权利要求所述为准。

Claims (10)

1.一种SPH与动态表面网格相结合的流体仿真方法，其步骤包括：

1）根据脚本在3D空间中对初始流体物质进行采样以生成粒子系统，计算由粒子定义的隐式函数

并由等值面

生成初始三角形网格表示流体的外表面形态，后续的流体运动仿真过程的计算分为粒子管线和表面管线两部分；

2）对于粒子管线，首先进行初始粒子设置以表示流体所处的初始状态，然后根据SPH方法计算粒子所受压强力、粘滞力和表面张力，并移动粒子，其中表面张力由显式表面上的信息所提供；然后对移动后的粒子系统重新计算隐函数

并得到其等值面

3）对于表面管线，首先初始化流体的显式表面，该显式表面为具有闭合流形性质的三角形网格，然后从粒子系统得到显式表面上每一个顶点的插值速度，并相应地移动顶点至新位置；

4）将步骤3）所得的运动后的网格表面投影到步骤2）重新计算后的等值面

上，在此过程中保证显式表面仍具有闭合流形性质且无明显的相交，同时处理有可能产生的拓扑变化事件，然后根据表面及两级距离场对流体进行重采样，从而完成当前时间步之内的流体仿真的所有计算，之后进入下一时间步的迭代，重复执行上述步骤2）-4）；

5）当仿真时间满足一定条件或者仿真状态满足一定条件时，整个流体仿真过程结束。

2.如权利要求1所述的方法，其特征在于：步骤2）通过粒子系统计算并得到其等值面

该过程首先定义在三维空间的标量函数

对于空间中一点x，计算其中i是粒子的下标，W是平滑函数，r是粒子i的影响域半径；在一个规则格子上对

进行采样，在进行采样之后，将粒子没有覆盖到的格点上的

都设成一个负值，以保证在有粒子覆盖的地方

值是正的，离粒子系统较远的地方则呈负值，等值面则可以用来定义粒子系统的表面。

3.如权利要求1所述的方法，其特征在于：步骤3）从粒子系统得到显式表面上每一个顶点的插值速度，该过程构造在空间任意一点x处的插值速度公式为其中W_i为第i个采样顶点的平滑函数，V_i为其体积，对于在表面附近粒子采样比较稀少存在采样不足的现象，采用简化的插值速度，对W进行归一化，即

减少粒子稀少对插值的影响。

4.如权利要求1所述的方法，其特征在于，步骤4）采用的投影算法包括如下步骤：

a）对于显式表面上的顶点v，其以插值速度运动之后空间位置为x，法向为n_v，在标量函数场上计算

和其中

r是粒子系统的初始设定半径，通常是k_r倍的粒子采样间距，如果

是精确的到某个闭合曲面的距离场，则使用

可得到更快的收敛速度；

b)如果

则在从x到x+εn_v的线段上用二分法迭代查找

的0值点；

c)如果

使得

表明此时顶点v的拓扑发生变化，此时如果

则将v移至x+εn_v，否则不移动v的位置，另外设置一个拓扑变化标志来提醒后续步骤的算法需要处理拓扑变化。

5.如权利要求1所述的方法，其特征在于：步骤4）在处理显式表面相交和拓扑事件时，不对显式表面地整个表面进行重建，而只对发生拓扑变化或是自交的区域进行局部重建，具体处理步骤包括：

a）在空间中建立一个正则格子；

b）在表面附近的格点上计算有向距离场；

c）识别发生拓扑事件以及自交的区域，标记为重网格化区域；

d）对表及区域重网格化，并修补与外部区域网格的接缝。

6.如权利要求5所述的方法，其特征在于，所述识别发生拓扑事件以及自交的区域的具体步骤包括：

a)依次检查一个单元是否含有复杂边、复杂面以及完全包含在单元内部的闭合曲面分支，判断是否满足复杂单元的条件；

b)距离隐式表面

超过一定距离的单元，称为深单元，从深单元开始向外扩展以包含周围所有的复杂单元，直到该区域的边界单元都是简单的，这样既能保证与拓扑变化相关的复杂单元都被标记为重网格化区域，同时尽可能保住了高分辨率的几何特征；

c)进行扩展的复杂单元的标记，即对Marching Cubes的15个模板中有4个是在一个格子面上允许有两条三角形边的情况标记为重网格化区域。

7.如权利要求5所述的方法，其特征在于，所述重网格化标记区域的步骤是：

a）处理标记区域的边界，使内部和外部的三角形片段沿区域边界分隔开；

b）删除标记区域内部三角形，使用Marching Cubes模板在区域内部重新建立三角形；

c）处理边界，将内部外部三角形片段缝合。

8.如权利要求1所述的方法，其特征在于，步骤4）进行重采样的步骤包括：

a)检测采样不足处的粒子，标记重采样区域；

b)在标记区域内进行粒子的快速泊松盘重采样并增加采样点；

c)重新分配重采样区域的粒子质量和动量；

d)合并符合质量条件、密度条件以及质量约束这三条合并条件的粒子。

9.如权利要求8所述的方法，其特征在于，所述粒子的快速泊松盘重采样方法的步骤包括：

a)在每个格子单元内部，随机生成t个探测位置，对每生成一个探测点，检测其是否离现有的粒子以及表面足够远，即探测点的位置x必须满足且|(x-x_i)|>k_pd₀且|x-v_j|>k_vd₀，其中

即x必须在流体内部，i和j分别是粒子与表面顶点下标，k_p和k_v是常数，d₀是仿真开始时粒子的初始间距；

b)如探测点满足以上要求，则将探测点处插入一个新的粒子，并将粒子映射到格子上，参数k_p和k_v限制可以插入粒子的区域，即不仅在流体内部，还要离粒子和表面一定的距离；

c)对新产生的粒子进行一步松弛操作以优化新粒子的位置，对于每个新产生的粒子p_i，其以r为半径的邻域内的粒子下标集合为NP_i，邻域内的表面顶点下标集合为NV_i，计算距离p_i与其周围顶点或粒子的最小距离 $d_{\max \min }=\min \left(\underset{j\&Element;N{P}_i}{\min }\right|x_i-x_j|,\underset{j\&Element;{\mathrm{NV}}_i}{\min }|x_i-x_j\left|\right)$ 之后，在p_i的邻域范围内随机撒若干探测点，对于每个探测点

如果其计算出的d_maxmin大于原来的值，则用的位置替换p_i；最后，p_i将被移动到其邻域内与周围顶点和粒子距离近似最大的位置。

10.如权利要求8所述的方法，其特征在于，重采样区域内的粒子的质量和动量的重新分配，其步骤包括：

a)将重采样区域进行聚类，并按照聚类结果进行分配，聚类的方法是设置一个队列并将一个重采样单元放入队尾，之后进行迭代循环，将队首单元弹出，加入集合S，并检查这个单元的邻接单元，如果也为重采样单元则加入队尾，直到队列为空，此时S中的单元聚为一类，如所有重采样单元都已被检查过，则聚类完成，否则继续将一个未检查的单元放入队尾进行以上操作；

b)在质量的分配上，采用的方法是平分一类内部的粒子质量，设粒子数为n_c，原粒子的下标集为P_c，粒子p_i的质量为m_i，则重采样后粒子上的采样质量为

计算完m′之后将其与一个质量阈值m_min比较，如有m′<m_min则舍弃该聚类内部的所有重采样点并停止重采样过程，维持聚类粒子状态不变；

c)在动量的分配上，在新插入的粒子处使用插值估算出速度

其中v_j是原粒子的速度，之后计算在一个类内部原粒子的合动量

然后通过调整粒子的速度对合动量进行归一化。

Images