CN104574472B - 基于嵌入网格的固体碎裂模拟和动画方法 - Google Patents

Abstract

本发明涉及一种基于嵌入网格的固体碎裂模拟和动画方法，其步骤包括：将仿真空间划分成离散的三维欧拉栅格，采用欧拉场模拟固体的形变和碎裂过程；采用三角形表面网格表示待模拟物体的几何外形，并将该表面网格嵌入对应的欧拉栅格单元中；当物体断裂、碎裂时，在每一个时间步内对物体的表面和内部进行采样并计算采样点的应力，当该应力的最大特征值超过设定的最大阈值时，在采样点位置初始化裂纹并创建裂纹；对每一条裂纹，跟踪裂纹前端的生长位置，以裂纹前端节点作为采样点来确定裂纹前端的生长速度和方向，进行动态裂纹生长，当所有裂纹前端节点和表面连接起来时裂纹终止。本发明具有物体表示精确、边界清晰、能够有效减少计算复杂程度的优点。

Description

基于嵌入网格的固体碎裂模拟和动画方法

技术领域

本发明属于计算机图形学与动画技术领域，具体涉及一种基于嵌入网格的固体碎裂模拟及其动画方法。

背景技术

随着过去几十年来计算机运算速度的提升，计算机图形学的相关技术飞速发展，基于物理的动画模拟在电影、游戏、建筑等各个行业得到了广泛的应用。在这个领域中主要分成液体和固体的动画模拟，液体模拟主要有基于粒子的系统和基于流体场来两个方面。相对于液体而言，固体的模拟方式可谓花样繁多。主要有：刚体(rigid body)，非形变的固体，主要用于碰撞计算，在游戏和动画仿真中应用很多。形变体(deformation)，可以在一定范围内改变自身的形状，可以用于判断固体的内力情况。布料(cloth)和纤维(fiber)，对于自碰撞检测与处理要求很高。以及其他一些固体的模拟方式。下面首先介绍一些常见的固体碎裂模拟算法。

常见的固体模拟方法是有限元法(finite element method，FEM)和无网格法(meshless methods)。

1)有限元法：在上世纪90年代末，在碎裂模拟中对于固体的表示方式主要基于有限元法。该方法将一个固体划分成若干个小正方体或四面体，当固体因为碰撞或其他原因而导致速度改变，就会产生形变。然后可以根据固体形变的大小计算出固体的应力(stress)以及应力的方向。

当应力超出系统设定的碎裂阈值时候，就可以沿着划分好的有限元边界产生裂纹。以四面体而言，当四面体的顶点分裂开时候固体就产生了碎裂。同时，通过将四面体进行重网格化，可以使得裂纹更加精细，真实。之后Muller(M.,MCMILLAN L.,DORSEY J.,JAGNOW R.,“Real-time simulation of deformation and fracture ofstiff materials.”In Proceedings of the Euro graphics workshop on Computeranimation and simulation,2001：113–124.)将该方法进行简化实现了实时的应用，在一些高端游戏中的固体碎裂采用了该算法。

2)无网格法(meshless methods)

无网格法和有限元法不同，不需要针对形变体进行体划分，而是采用采样粒子的方式来模拟固体，从而不用维护有限元法过程中有限元之间的邻居关系。同时，有限元的裂纹需要沿着有限元的边界生成，如果需要更精细的裂纹效果则需要进行有限元的重网格化，从而造成很大的开销。无网格法通过对固体内部进行动态采样，根据采样点周围的邻居点的关系建立裂纹。例如基于MLS(moving least squares)的无网格法(STEINEMANND.M.A.,GROSS M.OTADUY.“Fast arbitrary splitting of deforming objects.”Eurographics/ACM SIGGRAPH Symposium on Computer Animation,2006:27-34.)。但是该方法对于采样点的分布以及采样点周围的邻居的数量有很多限制。由于基于粒子的边界并非显式表示的，准确找到标示固体的表面的粒子对于无网格法也是一个很大的挑战。

发明内容

本发明提出了一种新方法来模拟固体的形变和碎裂过程：欧拉框架(Eulerianmodel)表示的固体模型之上嵌入表面网格(embedded surface mesh)进行固体物理运动现象模拟计算，从而产生精妙的固体形变、产生裂纹并最终碎裂的动画现象的仿真模拟。本发明的方法避免了在传统的有限元法(finite element method，FEM)中存在的重网格化(remesh)或是无网格法(meshless method)中重采样(resample)等一些问题。

本发明采用的技术方案如下：

一种基于嵌入网格的固体碎裂模拟和动画方法，其步骤包括：

1)将仿真空间划分成离散的三维欧拉栅格，采用欧拉场模拟固体在动力学上产生的形变和碎裂过程；采用三角形表面网格表示待模拟物体的几何外形，并将该表面网格嵌入对应的欧拉栅格单元中；

2)当物体由于挤压碰撞等物理作用而导致断裂、碎裂时，在每一个时间步内对物体的表面和内部进行采样，并计算采样点的应力，当该应力的最大特征值超过设定的最大阈值时，在采样点位置初始化裂纹并创建裂纹；

3)对每一条裂纹，跟踪裂纹前端的生长位置，以裂纹前端节点作为采样点来确定裂纹前端的生长速度和方向，进行动态裂纹生长，当所有裂纹前端节点和表面连接起来时，裂纹终止。

进一步地，步骤1)中嵌入欧拉栅格单元的表面网格由顶点坐标列表和三角形列表组成，并满足两个必要条件：完整闭包；面法向量连续且指向物体外。

进一步地，步骤2)初始化裂纹时，创建裂纹的两个表面，将裂纹前端连接的裂纹的两个表面分别称为正面和负面，在初始化时正面和负面可任意标记，但保证在之后的裂纹生长过程中正面始终保存在正面列表，负面始终保存在负面列表中。

进一步地，步骤2)创建的裂纹包括表面裂纹和内部裂纹；在嵌入表面网格表面上的采样点超过阈值时开始在表面创建裂纹，并进行表面路径查找及表面网格连接工作；内部裂纹不涉及路径查找以及和表面连接的过程。

进一步地，步骤2)在表面网格新建裂纹的过程包括：

a)寻找表面路径：当准备在表面上某采样点生成新裂纹时，通过计算获得该裂纹面的法向量和采样点的位置，以确定一个裂纹面；根据表面网格和裂纹面求得相交的三角形，通过限制找到的顶点离采样点的距离，在相交三角形中上确定两个裂纹前端节点F0和F1，再寻找到F0到F1在表面网格上的路径；

b)创建新裂纹面，即插入新点和新三角形到已存在的嵌入表面网格中；

c)记录裂纹前端节点的信息，即裂纹前端之间的相邻信息和每个裂纹前端节点连接的裂纹正负面的信息；

d)确定物体内部创建好的正负裂纹和表面三角形的连接关系，对表面网格进行重连接。

进一步地，步骤3)在进行裂纹生长时，当两个相邻裂纹前端节点距离过大或过小时进行插入或合并的操作，以避免创建出的裂纹网中的三角形形状失真；当插入新节点时，采用三次样条插值(cubic spline插值)法计算新节点的位置；当合并节点时，将合并节点的相关信息保存到已经合并的顶点中去。

进一步地，步骤3)所述裂纹终止分为两个步骤：正确选择表面顶点作为裂纹前端节点的连接点，以及相邻裂纹前端节点之间的三角形面片和表面面片的连接。

本发明通过采用三角形表面网格表示物体的几何外形，并将该表面网格嵌入欧拉框架表示的物理仿真计算单元－的栅格单元(grid cell)中，这样既能够精确地描述物体表面及其外观形状变化信息，又能够使物体本质的物理运动模拟过程更加准确；当物体由于挤压碰撞等物理作用而导致其外观形态发生改变时，尤其是产生断裂、碎裂等现象时，裂纹表面会根据所受外力作用动态地创建、并随时间生长、发展或终止；并通过一系列表面网格的拓扑操作动态维护正确的物体表面网格信息及其表示的几何形状。基于本发明的方法，可以处理非连通物体处于同一个欧拉栅格单元中时造成不同运动方式的奇异情况的问题，具有物体表示精确、边界清晰，同时又有效减少了计算复杂程度的优点。

附图说明

图1是应力的表示方式示意图。

图2是物体边界处特殊处理的二维示例图。

图3是特殊栅格单元的处理示意图。

图4是物体的三角形表示的示意图。

图5是嵌入表面网格上的裂纹初始化示意图。

图6是栅格单元内部裂纹初始化示意图。

图7是表面初始化的裂纹结果示意图。

图8是表面上裂纹的生成路径示意图。

图9是最短路径示例图。

图10是根据裂纹所对应的表面路径确定候选三角形顶点的示意图。

图11是初始化的表面裂纹示意图。

图12是裂纹的生长过程示意图。

图13是相互交叉的两个表面及包含裂纹前端节点的示意图。

图14是裂纹网格以及裂纹前端节点和表面连接的示意图。

图15是裂纹的生长过程的实例效果图。

图16是球体撞墙的断裂过程示例图。

图17是复杂表面情况下多裂纹结果示例图。

图18是将图17中被裂纹分割的多块物体拖动开后的效果图。

具体实施方式

为使本发明的上述目的、特征和优点能够更加明显易懂，下面通过具体实施例和附图，对本发明做进一步说明。

本发明创新性提出基于欧拉方法来模拟固体碎裂的方式。该方法避免了在传统的有限元法中的重网格化或是无网格法中重采样等一些问题。为了解决欧拉方法下确定固体表面的问题，采用在动力学计算模拟中嵌入显式表面网格的方案。通过将一个物理仿真计算的单元“欧拉栅格单元”根据嵌入表面网格的边界划分成多个部分，判定并确定在同一欧拉栅格单元内的物体分别只受到哪些栅格单元上速度和应力的影响。这样即确保了计算固体时边界的正确性，又解决了欧拉法难以确定固体表面的问题。

本发明方法的总体步骤如下：

1.数据预处理操作，针对待动画的三维模型及其初始运动状态，创建其欧拉场表示；

2.针对仿真动画过程中的每一时间步进行动画计算

2.1首先更新欧拉栅格单元。更新的具体步骤包括首先对栅格单元进行内外检测，对于每个位于物体内部的栅格单元执行下列三个步骤的操作，1)计算该栅格的新速度，2)增加体积力(body force)，3)计算栅格单元更新之后的新应力。

2.2根据栅格单元物理数值的变化，在嵌入三维网格上计算裂纹的生长并进行裂纹和表面网格的拓扑操作。具体步骤包括对每一条裂纹进行处理，裂纹的处理又包含以下步骤：对每条裂纹中的前端节点，首先计算前端节点的生长，然后前端节点和表面节点或其他裂纹表面连接，最后分裂前端节点。

2.3对三维网格采样计算新裂纹。具体步骤包括在物体表面和内部随机抽取采样点，针对每个采样点计算采样点特征值，如果该点特征值大于设定阈值，则创建新裂纹并增加该裂纹周围的阈值避免同一位置创建过多裂纹。

2.4更新物体的三维坐标和裂纹顶点坐标。具体步骤包括计算外插速度到物体外部栅格单元，对每一个栅格单元，计算该栅格单元的相对坐标；更新所有表面网格顶点的实际点空间坐标。

2.5三维网格表面碰撞检测与处理。

在每一时间步开始，需要对栅格单元进行内外检测，本发明采用“M.Fast and robust tracking of fluid surfaces.In Proceedings of the 2009ACMSIGGRAPH/Eurographics Symposium on Computer Animation,SCA’09,ACM,2009：237–245”中所述的内外检测方法。然后将后文所述的公式(1)分解成三步：平移物质坐标，增加体积力body force，计算应力。在平移物质坐标过程中采用“SELLE A.,FEDKIW R.,KIM B.,LIU Y.,ROSSIGNAC J.“An unconditionally stable MacCormack method.”J.Sci.Comput.35,2008：350–371.”的麦科马克格方法。之后采用后文提到的方法进行裂纹的拓扑操作，这是本发明的核心部分。在外插速度时，采用“ZHAO H.,“A fast sweepingmethod for eikonal equations.”MATHEMATICS OF COMPUTATION 74，2005：603–627”的FAST SWEEPING METHOD来计算需要外插速度的欧拉栅格单元到表面网格的有向距离。在计算平移速度时候，为了保证平移的精确性，采用四阶龙格库塔fourth order Runge Kutta(RK4)积分来获得数值解。最后处理在平移表面网格以后产生的表面网格碰撞的情况。碰撞的计算采用PQP库来处理物体的碰撞检测，PQP库可参见[LGLM99]LARSEN E.,GOTTSCHALKS.,LIN M.C.,MANOCHA D.,“Fast proximity queries with swept sphere volumes.”IEEE Intl Conf on Robotics and Automation，1999：1–32，以及GOTTSCHALK S.,LINM.C.,MANOCHA D.:Obbtree,“a hierarchical structure for rapid interferencedetection.”In Proceedings of the 23rd annual conference on Computer graphicsand interactive techniques,SIGGRAPH’96,ACM,1996：171–180。下面对各步骤进行详细说明。

1.基于嵌入表面网格的欧拉场模拟

1.1基于欧拉场的固体动力学模拟

现有技术中提出了一种基于欧拉场的固体形变模拟方法，参见KAMRINAND NAVE,J.-CK.“An eulerian approach to the simulation of deformable solids:Application to finite-strain elasticity.”ArXiv e-prints:http://arxiv.org/PS_cache/arxiv/pdf/0901/0901.3799v2.pdf.2009.，以及LEVIN D.I.W.J.,JONES G.L.,SUEDA S.,PAI D.K.LITVEN.“Eulerian solid simulation with contact.”ACMSIGGRAPH,2011：36-37.。前者将该欧拉场称作相对坐标映射(reference map)，然后利用有限差分法(finite difference method，FDM)模拟了二维空间中的固体形变。后者将该方法扩展到了三维空间，并采用该方法模拟了粘性体的接触碰撞和大形变。

本发明基于该欧拉框架来模拟物体由于受到力的作用而产生的精妙的破裂、裂纹变化以及碎裂等现象，提供了动态计算和表示的方法，并精确刻画和动画展示了上述物理现象随着时间序列而变化的过程。

一个固体的形变仿真需要计算出该物体在形变之前的相对坐标和物体产生形变后的相对坐标。然后类似于基于欧拉场的液体模拟，在物体内部设置一些粒子，通过该粒子在之后运动中坐标的偏移来描述固体的形变。为了计算物体受力之后的运动状态及其实际空间坐标位置，有如下的公式：

展开(1)式得到：

其中ρ表示物体的密度，v代表物体该点的速度，是关于v的物质导数，σ是该点的柯西应力张量(Cauchy Stress tensor)，该张量是用来计算物体受力情况的重要参数，如图1中应力的表示方式。b是该物体的体积力(body force)，例如重力g。通过最终计算出的速度可以计算物体粒子的实际坐标运动情况。

在计算σ时，需要获取任意位置x的相对坐标X，即X是关于实际空间坐标x的函数，记做X(x)，由于相对坐标在无外力时维持不变，可以通过下面的公式计算σ：

利用公式(3)计算出的结果可以计算计算形变梯度F：

为了计算形变物体的应力，首先要计算应力张量，可以通过F求得格林有限应变张量(Green finite strain tensor)E：

计算出应力张量之后，可以根据任意的本构律(constitutive law)来确定应力和应变张量之间的关系。常见的胡克定律(公式(6))是一种特殊的关系，即应力和应变张量二者成线性相关,主要适用于微小的形变状态：

σ＝εE (6)

应力和应变张量都是3×3的对称矩阵，故他们共有6个独立系数，可以展开成一个6×6的矩阵。对于各向同性的材料，胡克定律有如下形式：

其中ε代表杨氏模量，是描述物体材料抵抗形变能力的物理量，越大则越不容易产生形变。是泊松比率，用来表示物体材质参与计算量的属性。

在仿真计算过程中，将空间沿三维坐标轴方向划分成离散的三维栅格(3D grid)，每个栅格的大小相等。对于每一个栅格单元(grid cell)根据上面推导的结果，将X和速度v记录在每个栅格单元对应面的中心点上，和图1的表示方式相同。由于应力可以通过计算得到，所以不需要保存“GOKTEKIN T.G.,BARGTEIL A.W.,O’BRIEN J.F.A method foranimating viscoelastic fluids.ACM Trans.Graph.2004：463-468.”中提到的二阶应力张量(rank-two stress tensor)。

1.2嵌入表面网格(Embedded Surface Mesh)

动态追踪被仿真物体的变化表面对于大多数物理仿真而言是一项重要而困难的工作，例如针对采用粒子法仿真液体或固体时常见的marching cube算法(Lorensen W.E.,Cline H.E.Marching cubes:a high resolution 3d surface constructionalgorithm.In Proceedings of ACM SIGGRAPH 1987：163-169.)，但是该方法重构出的表面容易产生噪声和不光滑的问题，而且仿真计算出物体表面的时序动画很可能存在前后连续帧表面不一致的视觉瑕疵。对于大多数有限元法在计算固体形变时候已经将一个固体模型剖分成多个四面体，所以可以用四面体的非连通表面当作固体的表面。

本发明选择在仿真过程中嵌入显式的多边形网格表面来表现物体的形态。由于三角形表面网格是多边形网格的一种最常见的形式，因此后续有关本发明的方法的论述以三角形表面网格/三角形网格为最典型的代表，但是本发明所述方法同样适用于其他形式的多边形网格。三角形表面网格即常见的三维模型表示形式，是由顶点列表和三角形边组成。

三角形表面网格的所有顶点在每一时间步会根据欧拉场所计算出的结果进行移动。当处理欧拉场计算时，因为处理完全在物体内部的网格单元和部分在物体内部的网格单元显然需要采用不同的计算方式，显式的表面网格也为判断网格单元的拓扑结构提供了相对简单的识别方法。当产生裂纹的时候就需要对表面网格进行划分，从而更新嵌入表面网格的拓扑和几何位置。

1.3包含物体边界欧拉栅格单元的计算

对于一个欧拉栅格单元，当物体只有部分在该栅格单元内部时才需要进行特殊处理，例如图2中V点的情况。图中，P表示待计算的物体内的物质点，P_left表示两个紧挨着的栅格单元中的左侧单元中的点，当前图中位于物体内部，P_right表示两个紧挨着的栅格单元中的右侧单元中的点，当前图中位于物体外部，P_cut表示位于物体物质边界上的点。当不处于边界的条件时，计算P点时候通过P_left和P_right插值得到。而当处于边界时，因为P_right点不属于物体内部上的点，其值变得没有意义，该插值所得到的结果也就不准确了。此时的P点的计算就变成了通过P_left和P_cut插值获得，α表示切割长度的缩小比例：

当一个栅格单元内部出现两个甚至多个不连通的物体(或者一个物体分离开形成两个或者多个部分)时，正常情况下这些对象应该采用不同的运动方式。例如将一个物体朝两边撕裂，则物体的两部分应该向相反方向运动。但是基于欧拉仿真模型对于每一个栅格单元都只有一份数据，所以需要对这种类型的栅格单元进行特殊处理。首先可以通过表面三角形在格子内部的连通关系筛选出特殊的栅格单元。然后根据每一个连通表面三角形集合可以确认栅格单元记录的信息属于哪一个单元格。图3中(a)图表示一个单元格内部出现了连个不连通的物体，(b)图表示通过筛选得到右侧的单独物体，(c)图表示通过筛选得到左侧的单独物体。(c)图中处于左侧物体使用左侧栅格上具有有效信息的两个点进行计算，而(b)图中处于右侧物体只使用栅格右上角的角点的信息进行计算。确认使用的点后，在计算时将不使用的点临时设置为外部点的值来完成计算。

2.嵌入表面网格的拓扑操作

2.1嵌入表面网格的基本介绍

嵌入表面网格采用和基本的三维模型一样的表示形式，由顶点坐标列表和三角形列表及其拓扑结构关系组成。为了能够进行有效的物理仿真模拟计算，需要附加有如下2个必要条件：1.完整闭包；2.面法向量连续且指向物体外。

本发明所涉及的模拟仿真及其动画过程，物体表面需要是完整闭包的，否则将无法判断一个点究竟是在物体内还是物体外。而物体表面的法向量则代表了表面的内外关系。面法向量连续且指向物体外，该信息是通过三角形顶点的保存顺序确定的。如图4所示，其中(a)为平面连接，(b)为夹角连接，当保存顺序为[V₀ V₁ V₂]时，按照默认的右手系法则，三角形法向垂直指向纸外；当保存顺序为[V₀ V₂ V₁]时，该则法向指向纸内。确定所有三角形面片的法向正确对于裂纹生长以及物理计算十分重要。在确定所有相邻三角形法向是否正确时，可以通过相邻边的使用情况来判断。以边<V₁，V₂>为例，当三角形列表中存在边<V₁，V₂>的时候，一定存在一个相邻三角形保存边的顺序为<V₂，V₁>。并且在闭包的三角形网内这种边的使用情况总是成对出现。这种边的对称关系和两个相邻三角形的法向夹角无关，无论是图4(a)的法向平行还是图4(b)的夹角大于90°。所以在之后的裂纹创建等一系列操作中，要保证自己创建的新裂纹表面符合该条件，并且用该方法验证表面法向是否连续。

2.2固体裂纹创建

2.2.1裂纹产生的计算

本发明基于Rankine–Hugoniot条件根据最大应力初始化裂纹。在每一个时间步内对物体的表面和内部进行采样，并计算采样点的应力σ，当σ的最大特征值d_max超过设定的最大阈值时在采样点位置初始化裂纹。新创建裂纹的表面法向可以通过d_max获得。

2.2.2表面裂纹和内部裂纹

初始化裂纹的时候需要创建裂纹的两个表面，同时需要跟踪裂纹前端的生长位置，用来控制裂纹的生长形式。此处引入裂纹前端节点(Crack Front Node)用来标示裂纹的前端。前端节点记录了该前端连接的裂纹的两个表面。将这两个裂纹表面分别称为正面和负面，在初始化时候正面和负面可以任意标记，只要保证在之后的裂纹生长过程中正面始终保存在正面列表，负面始终保存在负面列表中即可。同时前端节点还保存相互之间的邻接关系。

当在嵌入表面网格表面上的采样点超过阈值的时候开始在表面创建裂纹。该裂纹初始化后的形态如图5所示，其中大致水平的染色平面网格A代表初始化就存在的表面网格。半透明染色平面B是由2.2.1章节计算得到的裂纹平面所在的位置。在C区域的周围的圆点是初始化得到的裂纹前端节点。C区域标示的染色网格是最终产生并插入到嵌入表面网格的新裂纹表面，该表面有2层，在初始化时候是重叠在一起的，之后根据1.3章节中的计算方式计算表面两侧物体的运动方式，从而实现向不同方向运动。

内部裂纹初始化比较简单，不涉及和表面连接的过程。初始化的结果如图6所示。

因为裂纹往往从表面开始生成，本发明在初始化物体内部裂纹的时候采用了更高的阈值。同时当初始化新裂纹后，会增加该裂纹存在区域的阈值来保证同一个区域内同一时间不会产生过多的裂纹。

2.2.3表面网格和裂纹的连接算法

图7为表面初始化的裂纹结果。其中(a)为由裂纹初始位置确定的表面连接路径，(b)为初始化后两端的裂纹前端节点(F0，F1)之间形成的三角形分开的结果。

表面网格新建裂纹可归纳为以下4方面：

1.寻找表面路径。当准备在表面上某采样点生成新裂纹时，通过计算获得了该裂纹面的法向量，采样点的位置，这样就可以确定一个裂纹面。根据表面网格和裂纹面可以求得相交的三角形。通过限制找到的顶点离采样点的距离，可以在相交三角形中上确定两个裂纹前端节点F0和F1。再寻找到F0到F1在表面网格上的路径(F0，V0，V1，V2，F1)，该路径查找算法在2.2.4章节介绍。

2.创建图5中C区域染色部分所指示的新裂纹面，即插入新点和新三角形到已存在的嵌入表面网格中。这里要保证创建的两个表面中对应位置的三角形法向方向相反。

3.记录裂纹前端节点的信息。即裂纹前端之间的相邻信息和每个裂纹前端节点连接的裂纹正负面的信息。

4.对表面网格进行重连接。如果只是进行了上面的操作，则虽然裂纹在内部可以继续生长，但是表面网格因为在V0，V1，V2处使用的是相同的顶点而无法分开。所以这里要对V0，V1，V2连接的三角形进行分离。以分裂V0点为例，根据2.1章节中介绍的表面法向连接方法，可以确定物体内部创建好的正负裂纹(图5)和表面三角形的连接关系。这里假设正面需要连接到F0到F1下侧的三角形边界上，根据确定连接的一侧三角形，就可以确定在(F0，V0，V1)下侧的所有包括V0顶点的三角形列表。然后将这个列表中原本指向V0的顶点指向新创建的顶点V3，裂纹面上指向V0的点也指向V3，即完成了V0点的分裂。相同处理V1，V2后就得到了图7的结果。这里裂纹正面起到了寻找一侧三角形的作用，只用使用正面或是反面则均可，只要保证全部分裂顶点使用相同的一面判断即可。

对于物体内部创建的裂纹由于不需要进行表面裂纹的寻找路径以及表面连接分裂工作，所以相对简单，这里不再阐述。

2.2.4表面路径查找算法

本节讨论给定表面网格上两个顶点F₁，F₂和经过两个顶点的一个平面P的法向量N_P，求在表面网格上近似于P和表面网格交线的一条近似折线段，如图8所示，其中直线表示表面上裂纹的切线，折线表示生成的路径。该方法基于三角形表面测地线算法，为了获得较高的运算速度，认为每一条路径的权重都相等。

表面路径本质上是一个无相图中两点路径的搜索，不过搜索范围限定在了一定的顶点之间。对于无相图的搜索算法有很多，例如贪心算法：在每个点计算最短距离。广度优先搜索算法：查找使用最少顶点到达的路径。Dijkstra算法：通过动态规划求得最短路径。

和这些路径搜索算法一样,为了计算表面路径，首先必须初始化保存点和点之间的邻接关系。同时为了确定三角形的相连关系，还要保存每个顶点连接的三角形列表。为了适应固体碎裂裂纹是沿给定裂纹表面方向产生的，不能简单的采用最短路径算法。以图9为例，假设计算出的裂纹生长平面理想情况是垂直于该椭球体，但是对于椭球体两侧的点，如果简单地根据网格表面距离计算最短路径，则该路径会错误地计算出如图9中所示的水平平面和椭球体的交线。这显然不是我们期望的结果。

通过观察可以发现，实际上是在全部顶点集V的一个特定子集V_p中寻找最短路径。顶点子集V_p具有条件：如果该三角形有一条经过顶点v的边和给定裂纹平面P存在交点，则v∈V_p。给定平面上一点p₁，法向量N_p，一个三角形的两个顶点V₁，V₂，可以计算顶点到平面的有向距离d₁，d₂：

d＝(V-p₁)·N_p (9)

如果当d₁,d₂之积小于等于0时，即两个顶点在平面上或分别在该平面的两侧，则认为该边和平面有交点。为了避免计算的误差，当二者任意点到平面距离较近时同样将该点加入到顶点子集中。由于三角形的边长可以近似认为相等，就可以在通过公式(9)确定的候选顶点子集(例如图10)中进行广度优先搜索来确定最短路径，再通过2.2.2章节介绍的方式建立裂纹。

最后建立的裂纹结果如图11所示。其中下半部分为原先表面顶点。中间区域的深色染色部分采用的新建顶点将表面和裂纹连接起来。在某些极特殊情况下该方法将无法查找到路径，在该方法查找路径失败的情况下再使用无限制的BFS获取路径。

2.3裂纹的动态生长

本发明采用与“PAULY M.,KEISER R.,ADAMS B.,P.,GROSS M.,GUIBASL.J.:Meshless animation of fracturing solids.ACM SIGGRAPH,2005：957-964.”中相同的显式表示裂纹前端的表示方式，将裂纹信息通过裂纹前端节点记录了下来。在处理裂纹的生长时只要以裂纹前端节点作为采样点来确定裂纹前端的生长速度和方向。对于任意的裂纹前端节点C_i，可以计算出它的生长方向向量为：

d_i＝α_iλ_i(v_i×t_i) (10)

其中α_i是根据物体材质属性来控制裂纹生长速度的参数，λ_i是该点应力张量的最大特征值，v_i是根据应力张量计算出的该点裂纹面的法向量。t_i是根据该节点相邻两个节点计算出的向量。

图12左侧图(a)显示了一个待生长的裂纹面，其中箭头所在位置为裂纹前端节点，箭头方向是通过公式(10)计算得出的生长方向。注意到这时裂纹面是有正反两层表面，两层表面的法向相反，且只有在裂纹前端(即图中边界位置)的顶点处相交。右侧图(b)是生长之后的新裂纹。裂纹前端移动到了新的边界位置，并且将旧位置顶点分为两份，将正面和负面连接到对应的顶点上完成裂纹正负面的分裂。

当两个相邻裂纹前端节点距离过大或过小时，需要插入或合并的操作。否则创建出的裂纹网中的三角形形状会失真。当插入新节点时，采用三次样条插值(cubic spline插值)法计算新节点的位置。该算法根据给定的节点列表确定一条闭合曲线，根据插值的需要获取到一个近似在曲线上的插值点。

当合并节点时候，需要将合并节点的相关信息保存到已经合并的顶点中去。图12右侧裂纹可以看到上面的两个节点因为距离过大，而插入了一个新前端节点。而下面的两个节点因为距离过小，合并成了一个节点。

2.4裂纹生长终止

裂纹的终止对于碎裂是十分重要的一步，只有裂纹的两端完全分离物体才能产生实际的碎裂效果。此处阐述在一个固体内部单一生长一个裂纹直到最后完全分开的过程。裂纹的终止在本发明的实现方法中即裂纹前端节点和已经存在的表面网格的连接过程。当所有裂纹前端节点都和表面连接起来时，裂纹完全终止。此时物体的表面网格会分为两块不连通的新表面网格。

对于裂纹终止而言，实际上分为两个步骤：正确选择表面顶点作为裂纹前端节点的连接点，以及相邻裂纹前端节点之间的三角形面片和表面面片的连接。

2.4.1表面顶点选择

对于一个裂纹前端节点来说，如何判断该节点需要连接到表面是本节最重要的问题。当对裂纹前端节点进行表面连接时，实际上是将裂纹前端节点和表面上一点连接起来的过程。需要避免在连接中发生以下几个问题：

1.连接到自身已经使用的顶点，例如节点自己之前时间步使用过的顶点。

2.连接到非连通的顶点，例如一个节点连接到另外一个裂纹上，因为另外一个裂纹还没分开，所以是有2个很接近的平行表面，即上面所说的正负面。这个节点本来处于裂纹的正面一侧，如果连接到了另外一个裂纹的反面，则造成了连接错误。

3.生长过程中插入和合并节点的新位置问题。该新位置有可能直接被插入到了物体的外侧，当两个物体表面连接很近的时候甚至可能被插入到另外一个物体内。

问题1在布料等可产生特大形变的材料模拟中可能会考虑处理，准备模拟的材料并不具有那么大形变的可能。通过记录自身使用过的节点，在寻找要连接的表面顶点时过滤掉这些本裂纹使用过的节点。

对于问题2，要处理的情况相对复杂。首先在计算过程中，可能无法避免两个表面出现相互交叉的情况，如图13的(a)图所示：左侧和右侧分别是两块不同的区域，他们在中间相互交叉。而如果此时在左侧区域有一个裂纹，如图13的(b)图所示，则当简单地按照距离或是射线方式判断节点该连到的表面顶点时，该裂纹前端上面的2个节点会正常连接到左侧表面。中间一个节点该连接到哪个表面就不太容易确定，而下面两个节点显然会被连接到右侧表面。所以简单地通过裂纹前端节点和表面上顶点的距离来判断要连接的顶点是不合适的。而通过裂纹前端节点和前端节点生长速度方向确定一条射线，通过和射线相交的第一个表面连接也无法应对图13中下侧两个节点的问题。为了解决这个问题，本发明采取了多种策略组合的判断方式，通过多步筛选去掉不合适的顶点，最终在剩余顶点中选取合适节点进行连接。

首先根据空间划分的欧拉栅格单元将顶点和裂纹前端节点放入对应单元中去。然后选取裂纹前端节点周围的栅格单元，计算这些栅格单元中顶点到裂纹前端节点之间的距离，将距离低于设定阈值的顶点作为候选，同时根据问题1筛除掉自身裂纹的顶点。

然后可以根据节点的位置和节点的速度得到一个向量N_c，用该向量计算相交的表面上的三角形面片。通过计算三角形面片的法向量N_p和N_c的夹角大小和节点位置在该三角形面片所在平面的里侧或是外侧来判断是否选取该三角形上的顶点作为候选点。

最后通过有步数限制的广度优先搜索判断候选顶点是否和周围已连接节点可达，从可达的候选顶点中选择距离最近的顶点作为该裂纹前端节点连接的表面顶点。

2.4.2裂纹网和表面网格的连接

裂纹网和表面网格在数据形式上是相同的顶点三角形模式。不同的是裂纹网的前端顶点记录在了裂纹前端节点中。根据裂纹插入和合并阈值的不同设置，裂纹网和表面网格在密度上会存在差异，图14的(a)图为一个从球体左端开始生长的裂纹面。当相邻两个裂纹前端节点都连接到表面时就需要将这两个相邻节点之间的裂纹三角形和表面三角形进行连接，图14的(b)图显示了裂纹前端节点和表面连接的实例。可以看到上面两个节点已经和表面连接，因为该示例中表面网格比裂纹网要密，所以原本两个前端之间只有一组(正负两面)的裂纹面被重新划分成了多组。裂纹前端节点在表面上的路径搜索和连接和创表面建裂纹的操作非常类似，也是先寻找路径，然后再创建细化后的三角形同表面连接。详细路径寻找和连接算法见本文2.2章节。

当一个裂纹前端节点连接到表面上时候，他所连接的这个节点还暂时不能分裂成两份。图14的(b)图中从上往下第一个节点已经分成了两份，即该节点位置上有两个顶点，一个连接表面上一侧顶点和裂纹的正面，一个连接表面另一侧顶点和裂纹负面。第二个节点因为他的邻居(第三个节点)还没有连接到表面，故暂时还不能将之分裂。当裂纹前端节点的左右邻居都连接到表面以后才能进行分裂。分裂方法和2.2.2章节中创建裂纹时分裂路径经过顶点的方式一致。最终当一个裂纹的全部裂纹前端顶点都和表面或是其他裂纹的正负面连接完毕，该裂纹完全终止，物体被分为两半。当然也可以设定裂纹前端处的应力阈值使得裂纹可以生长到一半停止，这样物体上就可以产生未完全裂开的裂纹。

2.5结果示例

2.5.1单一裂纹示例

一个球体运动撞墙后产生裂纹到裂纹完全生长结束，将物体分为两半。裂纹完全终止后的球断裂成的两部分。图15为裂纹的生长过程示意图，图16为球体撞墙的断裂过程示意图，图16中(a)为球体撞墙产生断裂，(b)为断裂球体下侧，(c)为断裂球体上侧。

2.5.2多裂纹示例

Armadillo模型(国际通用测试三维模型)撞墙后由背部一点产生多条裂纹，如图17所示。图18为将被裂纹分割的多块物体拖动开后的效果。未分开的为裂纹停止生长的情况。

2.6结果分析

第一个示例为一个球体垂直撞上墙后产生形变，当应力的最大值超过阈值之后产生了一条裂纹。随着碰撞的继续和应力的传播，裂纹逐步生长，最终将球体分裂成两个不连通区域。实现碎裂的过程。该结果反映了固体的形变以及裂纹的生长。

第二个示例为一个复杂固体撞上墙后产生多条裂纹的过程，该示例和第一个示例在物理计算上完全相同。在裂纹扩展上取消了只产生一个裂纹的限制，实现了多个裂纹生长，并相互交叉连接，符合实际裂纹生长于分裂的需求。

本发明的仿真模拟方法的由于采用嵌入式表面网格，具有了如下特殊的优点：1)不需要对物体进行四面体划分，2)也不需要通过粒子重建表面，提供了一种新型的物体碎裂的表面的表示方式。该方式实时速度较快，并且可以提供比有限元法更丰富的表面变化，对于实时性的碎裂模拟来说是也一种新的思路。

上述步骤中涉及的常微分方程组数值解的求解方法，除了龙格库塔法之外的其他可行的求解方法也适用于本发明。

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

Claims (10)

1.一种基于嵌入网格的固体碎裂模拟和动画方法，其步骤包括：

1)将仿真空间划分成离散的三维欧拉栅格，采用欧拉场模拟固体在动力学上产生的形变和碎裂过程；采用三角形表面网格表示待模拟物体的几何外形，并将该表面网格嵌入对应的欧拉栅格单元中；

2)在仿真模拟和动画过程中，当物体受到物理作用而导致断裂、碎裂时，在每一个时间步内对物体的表面和内部进行采样，并根据欧拉场固体动力学原理计算采样点的应力，当该应力的最大特征值超过设定的最大阈值时，在采样点位置初始化裂纹并创建裂纹；

3)对每一条裂纹，跟踪裂纹前端的生长位置，以裂纹前端节点作为采样点来确定裂纹前端的生长速度和方向，进行动态裂纹生长，当所有裂纹前端节点和表面连接起来时，裂纹终止。

2.如权利要求1所述的方法，其特征在于，步骤1)中嵌入欧拉栅格单元的表面网格由顶点坐标列表和三角形列表组成，并满足两个必要条件：完整闭包；面法向量连续且指向物体外。

3.如权利要求1所述的方法，其特征在于：步骤2)初始化裂纹时，创建裂纹的两个表面，将裂纹前端连接的裂纹的两个表面分别称为正面和负面，在初始化时正面和负面可任意标记，但保证在之后的裂纹生长过程中正面始终保存在正面列表，负面始终保存在负面列表中。

4.如权利要求1所述的方法，其特征在于：步骤2)创建的裂纹包括表面裂纹和内部裂纹；在嵌入表面网格表面上的采样点超过阈值时开始在表面创建裂纹，并进行表面路径查找及表面网格连接工作；内部裂纹不涉及路径查找以及和表面连接的过程。

5.如权利要求1所述的方法，其特征在于，步骤2)在表面网格新建裂纹的过程包括：

a)寻找表面路径：当准备在表面上某采样点生成新裂纹时，通过计算获得该裂纹面的法向量和采样点的位置，以确定一个裂纹面；根据表面网格和裂纹面求得相交的三角形，通过限制找到的顶点离采样点的距离，在相交三角形中确定两个裂纹前端节点F0和F1，再寻找到F0到F1在表面网格上的路径；

b)创建新裂纹面，即插入新点和新三角形到已存在的嵌入表面网格中；

c)记录裂纹前端节点的信息，即裂纹前端之间的相邻信息和每个裂纹前端节点连接的裂纹正负面的信息；

d)确定物体内部创建好的正负裂纹和表面三角形的连接关系，对表面网格进行重连接。

6.如权利要求5所述的方法，其特征在于，所述步骤a)在全部顶点集V的一个特定子集V_p中寻找最短路径；顶点子集V_p具有条件：如果该三角形有一条经过顶点v的边和给定裂纹平面P存在交点，则v∈V_p；给定平面上一点P₁，法向量N_p，一个三角形的两个顶点V₁,V₂，计算顶点到平面的有向距离d₁,d₂：

d＝(V-P₁)·N_p，

如果当d₁,d₂之积小于等于0时，即两个顶点在平面上或分别在该平面的两侧，则认为该边和平面有交点，当二者任意点到平面距离较近时同样将该点加入到顶点子集中；三角形的边长近似认为相等，在通过上述公式确定的候选顶点子集中进行广度优先搜索来确定最短路径。

7.如权利要求1所述的方法，其特征在于：步骤3)在进行裂纹动态生长时，对于任意的裂纹前端节点C_i，其生长方向向量为

d_i＝α_iλ_i(v_i×t_i)，

其中α_i是根据物体材质属性来控制裂纹生长速度的参数，λ_i是该点应力张量的最大特征值，v_i是根据应力张量计算出的该点裂纹面的法向量，t_i是根据该节点相邻两个节点计算出的向量。

8.如权利要求1所述的方法，其特征在于：步骤3)在进行裂纹生长时，当两个相邻裂纹前端节点距离过大或过小时进行插入或合并的操作，以避免创建出的裂纹网中的三角形形状失真；当插入新节点时，采用三次样条插值法计算新节点的位置；当合并节点时，将合并节点的相关信息保存到已经合并的顶点中去。

9.如权利要求1所述的方法，其特征在于：步骤3)所述裂纹终止分为两个步骤：正确选择表面顶点作为裂纹前端节点的连接点，以及相邻裂纹前端节点之间的三角形面片和表面面片的连接。

10.如权利要求9所述的方法，其特征在于：采取多种策略结合的判断方式选择表面顶点，通过多步筛选去掉不合适的顶点，最终在剩余顶点中选取合适节点进行连接，具体方法是：

首先根据空间划分的欧拉栅格单元将顶点和裂纹前端节点放入对应单元中去，然后选取裂纹前端节点周围的栅格单元，计算这些栅格单元中顶点到裂纹前端节点之间的距离，将距离低于设定阈值的顶点作为候选，同时筛除掉自身裂纹的顶点；

然后根据节点的位置和节点的速度得到一个向量N_c，用该向量计算相交的表面上的三角形面片，通过计算三角形面片的法向量N_p和N_c的夹角大小和节点位置在该三角形面片所在平面的里侧或是外侧来判断是否选取该三角形上的顶点作为候选点；

最后通过有步数限制的广度优先搜索判断候选顶点是否和周围已连接节点可达，从可达的候选顶点中选择距离最近的顶点作为该裂纹前端节点连接的表面顶点。