一种sph与动态表面网格相结合的流体仿真方法【专利摘要】本发明涉及一种SPH与动态表面网格相结合的流体仿真方法。该方法使显式表面追踪用粒子系统采样的流体，并在仿真过程中自动处理自交和拓扑变化，维护三角形网格的质量；该方法提出了新的基于显式表面的自适应重采样方法。SPH仿真不易展现出很细很薄的特征，尤其是在粒子分辨率不足的情况下，甚至还会形成不自然的空洞。在传统SPH基础上如果想体现出更好的细节，可以增大粒子的分辨率，但是增大分辨率会相应地增加时间和空间上的计算代价，本发明基于显式表面的自适应重采样方法既可以模仿出非常精微的流体表面细节从而提升仿真的效果，同时又可以提高仿真的效率和速度。【专利说明】一种SPH与动态表面网格相结合的流体仿真方法【
技术领域：
】[0001]本发明属于计算机图形与动画【
技术领域：
】，涉及一种流体仿真数据表示方法以及流体运动和动画的高效仿真计算方法。【
背景技术：
】[0002]在自然界以及日常生活中，流体是很常见的物质形态。在很多计算机图形学应用中，流体的仿真也占据了很重要的地位，例如在动画、游戏等娱乐领域，以及在科学可视化和医疗应用领域的应用等。现在计算机图形学上的流体仿真与动画方法有很多种，主流的有基于欧拉观点的有限差分方法(通称欧拉法)和基于拉格朗日观点的粒子法(无网格法，通称拉格朗日法)。两种方法各有优缺点，各有应用。粒子法在描述流体的细节特征(波动以及涡旋)方面具有天然的优势。[0003]欧拉法(Euler法)最早由Stam在1999年引入图形领域(JosStam.Stablefluids[C].Proceedingsofthe26thannualconferenceonComputergraphicsandinteractivetechniques.ACMPress/Addison-ffesleyPublishingC0.,1999,121-128)。在这篇文章中作者描述了使用有限差分法进行流体仿真的算法流程，给出了在Euler框架下求解Navier-stokes方程的方法，并给出了平流和投影的数值解法。这篇文章提供了在图形学中进行基于物理的流体仿真的基本框架。[0004]SPH方法(光滑粒子动力学法)是一种基于Lagrange观点的典型的无网格方法。2003年M--ulIer将该方法引入计算机图形学中进行流体仿真(MatthiasM"uller,DavidCharypar,MarkusGross.Particle-basedfluidsimulationforinteractiveapplications[C]?Proceedingsofthe2003ACMSIGGRAPH/EurographicssymposiumonComputeranimation.EurographicsAssociation,2003，154-159)。SPH方法避免了基于网格方法中频繁的重网格化操作，能够增加程序的鲁棒性，降低计算复杂度。SPH方法在连续介质特别是流体的仿真领域中有广泛的应用。图形学中的流体仿真有多种方法，其中比较常用的是基于Euler观点的有限差分法以及基于Lagrange方法的无网格方法(SPH)。相比于Euler方法，SPH方法能更好地展现流体的细节,例如流体的飞派效果。另外,通过调节粒子的规模可以很容易的使用SPH方法达到实时仿真而不使效果有很大损失。[0005]动态显式表面的主要思想是在仿真的初始阶段给定一个满足闭合流形条件的初始显式表面(通常即为三角形网格)，随着仿真的不断进行，每个时间步都得到一个空间中的速度场，由该速度场驱动显式表面的顶点进行运动并保持连接关系不变。之后得到了一个可能有相交以及退化三角形的网格，对此网格上的三角形进行操作以保证三角形的质量(由面积、边长比、角度等指标来衡量)，同时修正已经相交的部分使网格仍然保持闭合流形性质，在此过程中应能处理由于运动引起的拓扑变化。由于动态显式表面的输入仅需要初始网格和速度场，与底层的仿真无关，只要是能提供速度场的数值方法都可以使用。又由于动态显式表面可以很好地处理形变以及拓扑变化，很适合应用在流体仿真中作为流体的表面。传统的流体仿真使用隐式表面，即通过仿真模型生成一个距离场，再从距离场上抽取出一个等值面作为流体的表面。显式表面与隐式表面相比有如下好处:[0006]?显式表面比隐式表面更有利于保持住表面细节，这一点上要优于动态隐式表面，后者常常会因为平流中的数值耗散而抹掉细节；[0007]?显式表面能够很好地表示几何上退化的特征，例如三维场景下的一维或二维特征。而隐式表面则受限于空间离散化的分辨率大小，往往很难表现这类特征；[0008]?显式表面能够很好地保持表面上的参数化特征，例如纹理坐标等等。使用动态显式表面意味着可以省去每步重建表面的过程，不仅可以节省大量的内存空间(以及通常会节省计算时间)，更重要的是可以天然地保持住表面上各点在时间序列上的演变关系，从而可以使参数化特征的传递有迹可循；[0009]?对于流体来讲，引入显式表面可以有助于以一种更直接的方式来计算表面物理，例如由于表面张力引起的表面压差和毛细波等现象。此外，有了显式表面，可以借助离散微分几何上的大量知识来理解当前的仿真系统，例如可以便捷地计算表面积以及体积等等；[0010]?随着仿真的进行，显式表面有时间维度上的连续性。由于在每一帧中隐式表面是作为当时时间步的后处理操作来进行的，即隐式表面的形状只取决于那一时间步生成隐式表面的等值面的形状。有时在不同帧之间由于粒子分布的不同等值面的变化可能不连续，这体现为视觉效果上闪动。而显式表面则是由上一时间步的表面运动而来，保留了上一时间步的状态，很大程度上可以避免这种闪动的发生。[0011]在使用SPH进行仿真计算时，如果粒子系统在局部受到拉伸作用，会产生拉伸不稳定(TensileInstability)现象。这种现象的产生是由SPH的计算模型所决定的,表现出来的结果就是粒子系统在比较细、薄的特征处的粒子运动会趋向于不稳定，可能彻底分散或是聚成小团。我们希望在粒子变得过于稀疏之前可以对粒子系统进行自适应的重采样，从而使仿真可以展现出一些细薄特征。【
发明内容】[0012]由上述内容可知，基于SPH思想的仿真方法不易表示和展现出很细很薄的特征，尤其是在粒子分辨率不足的情况下。在传统SPH方法基础上如果想体现出更好的细节，可以增大粒子的分辨率，但是增大分辨率会相应地增加时间和空间上的计算代价。如果能够实现自适应的重采样则可以把有限的计算资源根据细节的粗细进行合理地分配。本发明针对该问题，提出了一种显式表面与SPH流体相结合的流体仿真方法，可以使显式表面追踪用粒子系统采样的流体，并在仿真过程中自动处理自交和拓扑变化，维护三角形网格的质量。[0013]具体来说，本发明采用的技术方案如下:[0014]一种SPH与动态表面网格相结合的流体仿真方法，其步骤包括:[0015]I)根据脚本在3D空间中对初始流体物质进行采样以生成粒子系统，计算由粒子定义的隐式函数V并由等值面穸=0生成初始三角形网格表示流体的外表面形态，后续的流体运动仿真过程的计算将分为粒子管线和表面管线两部分组成；[0016]2)对于粒子管线，首先进行初始粒子设置以表示流体所处的初始状态，然后根据SPH方法计算粒子所受压强力、粘滞力和表面张力，并移动粒子，其中表面张力由显式表面上的信息所提供；然后对移动后的粒子系统重新计算隐函数V并得到其等值面-=Oi[0017]3)对于表面管线，首先初始化流体的显式表面，该显式表面为具有闭合流形性质的三角形网格，然后从粒子系统得到显式表面上每一个顶点的插值速度，并相应地移动顶点至新位置；[0018]4)将步骤3)所得的运动后的网格表面投影到步骤2)重新计算后的等值面0上，在此过程中保证显式表面仍具有闭合流形性质且无明显的相交，同时处理有可能产生的拓扑变化事件，然后根据表面及两级距离场对流体进行重采样，从而完成当前时间步之内的流体仿真的所有计算，之后进入下一时间步的迭代，重复执行上述步骤2)-4)；[0019]5)当仿真时间满足一定条件或者仿真状态满足一定条件时，整个流体仿真过程结束。[0020]进一步的，步骤2)中通过粒子系统计算-并得到其等值面?=0,该过程首先定义在三维空间的标量函数勢对于空间中一点X，计算【权利要求】1.一种SPH与动态表面网格相结合的流体仿真方法，其步骤包括:1)根据脚本在3D空间中对初始流体物质进行采样以生成粒子系统，计算由粒子定义的隐式函数妒并由等值面黌=O生成初始三角形网格表示流体的外表面形态，后续的流体运动仿真过程的计算分为粒子管线和表面管线两部分；2)对于粒子管线，首先进行初始粒子设置以表示流体所处的初始状态，然后根据SPH方法计算粒子所受压强力、粘滞力和表面张力，并移动粒子，其中表面张力由显式表面上的信息所提供；然后对移动后的粒子系统重新计算隐函数f并得到其等值面-=Oi3)对于表面管线，首先初始化流体的显式表面，该显式表面为具有闭合流形性质的三角形网格，然后从粒子系统得到显式表面上每一个顶点的插值速度，并相应地移动顶点至新位置；4)将步骤3)所得的运动后的网格表面投影到步骤2)重新计算后的等值面-=0上，在此过程中保证显式表面仍具有闭合流形性质且无明显的相交，同时处理有可能产生的拓扑变化事件，然后根据表面及两级距离场对流体进行重采样，从而完成当前时间步之内的流体仿真的所有计算，之后进入下一时间步的迭代，重复执行上述步骤2)-4)；5)当仿真时间满足一定条件或者仿真状态满足一定条件时，整个流体仿真过程结束。2.如权利要求1所述的方法，其特征在于:步骤2)通过粒子系统计算-并得到其等值面黌=0,该过程首先定义在三维空间的标量函数料对于空间中一点X，计算其中i是粒子的下标，W是平滑函数，r是粒子i的影响域半径；在一个规则格子上对史进行采样，在进行采样之后，将粒子没有覆盖到的格点上的炉都设成一个负值，以保证在有粒子覆盖的地方妒值是正的，离粒子系统较远的地方则呈负值，V=O等值面则可以用来定义粒子系统的表面。3.如权利要求1所述的方法，其特征在于:步骤3)从粒子系统得到显式表面上每一个顶点的插值速度，该过程构造在空间任意一点X处的插值速度公式为4.如权利要求1所述的方法，其特征在于，步骤4)采用的投影算法包括如下步骤:a)对于显式表面上的顶点V，其以插值速度运动之后空间位置为X，法向为nv，在标量函数场w上计算树-r)和斜X—--,.)，其中C=Algfl(树r是粒子系统的初始设定半径，通常是b倍的粒子采样间距，如果妒是精确的到某个闭合曲面的距离场，则使用e=sign{(p{x))?可得到更快的收敛速度；b)如果则在从X到X+env的线段上用二分法迭代查找黌的0值点；c)如果3x,使得钱X)嗔x+env)>0，表明此时顶点v的拓扑发生变化，此时如果I史(x+env)j<I斜x)|，则将v移至X+env,否则不移动v的位置，另外设置一个拓扑变化标志来提醒后续步骤的算法需要处理拓扑变化。5.如权利要求1所述的方法，其特征在于:步骤4)在处理显式表面相交和拓扑事件时，不对显式表面地整个表面进行重建，而只对发生拓扑变化或是自交的区域进行局部重建，具体处理步骤包括:a)在空间中建立一个正则格子；b)在表面附近的格点上计算有向距离场；c)识别发生拓扑事件以及自交的区域，标记为重网格化区域；d)对表及区域重网格化，并修补与外部区域网格的接缝。6.如权利要求5所述的方法，其特征在于，所述识别发生拓扑事件以及自交的区域的具体步骤包括:a)依次检查一个单元是否含有复杂边、复杂面以及完全包含在单元内部的闭合曲面分支，判断是否满足复杂单元的条件；b)距离隐式表面炉=O超过一定距离的单元，称为深单元，从深单元开始向外扩展以包含周围所有的复杂单元，直到该区域的边界单元都是简单的，这样既能保证与拓扑变化相关的复杂单元都被标记为重网格化区域，同时尽可能保住了高分辨率的几何特征；c)进行扩展的复杂单元的标记，即对MarchingCubes的15个模板中有4个是在一个格子面上允许有两条三角形边的情况标记为重网格化区域。7.如权利要求5所述的方法，其特征在于，所述重网格化标记区域的步骤是:a)处理标记区域的边界，使内部和外部的三角形片段沿区域边界分隔开；b)删除标记区域内部三角形,使用MarchingCubes模板在区域内部重新建立三角形；c)处理边界，将内部外部三角形片段缝合。8.如权利要求1所述的方法，其特征在于，步骤4)进行重采样的步骤包括:a)检测采样不足处的粒子，标记重采样区域；b)在标记区域内进行粒子的快速泊松盘重采样并增加采样点；c)重新分配重采样区域的粒子质量和动量；d)合并符合质量条件、密度条件以及质量约束这三条合并条件的粒子。9.如权利要求8所述的方法，其特征在于，所述粒子的快速泊松盘重采样方法的步骤包括:a)在每个格子单元内部，随机生成t个探测位置，对每生成一个探测点，检测其是否离现有的粒子以及表面足够远，即探测点的位置X必须满足树x}>0且I(X-Xi)|>kpdQ且Ix-VjlSkvCltl,其中树4>0即X必须在流体内部，i和j分别是粒子与表面顶点下标，kp和kv是常数，d0是仿真开始时粒子的初始间距；b)如探测点满足以上要求，则将探测点处插入一个新的粒子，并将粒子映射到格子上，参数kp和kv限制可以插入粒子的区域，即不仅在流体内部，还要离粒子和表面一定的距离；c)对新产生的粒子进行一步松弛操作以优化新粒子的位置，对于每个新产生的粒子Pi,其以r为半径的邻域内的粒子下标集合为NPi,邻域内的表面顶点下标集合为NVi,计算距离Pi与其周围顶点或粒子的最小距离⑶:.r,.,v；之后，在Pi的邻域范围内随机撒若干探测点，对于每个探测点‘，如果其计算出的Clnlaxnlin大于原来的值，则用ft的位置替换Pi;最后，Pi将被移动到其邻域内与周围顶点和粒子距离近似最大的位置。10.如权利要求8所述的方法，其特征在于，重采样区域内的粒子的质量和动量的重新分配，其步骤包括:a)将重采样区域进行聚类，并按照聚类结果进行分配，聚类的方法是设置一个队列并将一个重采样单元放入队尾，之后进行迭代循环，将队首单元弹出，加入集合S，并检查这个单元的邻接单元，如果也为重采样单元则加入队尾，直到队列为空，此时S中的单元聚为一类，如所有重采样单元都已被检查过，则聚类完成，否则继续将一个未检查的单元放入队尾进行以上操作；b)在质量的分配上，采用的方法是平分一类内部的粒子质量，设粒子数为n。，原粒子的下标集为P。，粒子Pi的质量为IV则重采样后粒子上的采样质量为【文档编号】G06T17/00GK103714575SQ201310744592【公开日】2014年4月9日申请日期:2013年12月30日优先权日:2013年12月30日【发明者】李胜,任毅,王衡,汪国平申请人:北京大学