CN104715499A - 一种层次几何细分的各向异性材料脆性破裂模拟方法 - Google Patents

Abstract

本发明公开了一种层次几何细分的各向异性材料脆性破裂模拟方法，其步骤为：1）基于弹性力学与自适应SPH的弹力分析：采用经典的弹性力学理论，自适应SPH作为数值方法，计算位移场及对应的应力和应变；2）基于能量分析的裂纹计算：通过用户预设值将碰撞点的邻居粒子分组，分别计算及裂纹方向，根据能量确定裂纹内外半径；3）根据裂纹模式，自适应的细分四面体粒子，从而得到局部高分辨率的破碎结果，并进一步渲染得到最终效果。本发明能够快速生成高细节的各向异性材料破碎细节，并提供了便于用户友好的裂纹控制模式；本发明继承了FEM方法的精确分析与SPH方法的快速简洁，具有一定的实用价值。

Description

一种层次几何细分的各向异性材料脆性破裂模拟方法

技术领域

本发明属于计算机图形学领域，具体地说是一种基于物理模型的计算机动画模拟方法，其部分技术涉及到有限元剖分算法，各向异性的自适应SPH算法，四面体几何细分算法以及弹性力学相关理论。

背景技术

脆性破碎，如破碎的玻璃、陶瓷器，是电影特效和计算机游戏中常见的现象。物理准确且细节丰富的破碎模拟在计算机图形学领域仍然是一个困难的问题。当采用弹性模型处理破碎时，其数值复杂性与不稳定性，导致了严重的计算消耗。除此之外，用户难以控制裂纹的生成与传播。作为简化的替代方法，非物理的预先设定裂纹的方法常用于游戏与电影产业，但其模拟结果欠缺真实感，且人工设计是非常耗时的。所以在基于物理的破碎动态模拟中减少计算消耗，并保持高质量的可视效果是十分必要的。

近年来，计算机图形学领域有很多脆性破碎的相关工作。有限元方法(FEM)，作为经典的网格方法，广泛应用于基于物理的固体变形模拟，但是维持高质量网络是需要极大计算开销。另一方面，无网格方法以及实现简单、交互方便的优点广受欢迎。非结构化的模式使得它适用于模拟拓扑改变现象，如流体、破碎。值得注意的是，过去的相关工作中，可变性固体的模拟或采用计算复杂的隐式方法，或选择不稳定的显式方法。

考虑到各向异性材料的特征，其脆性破裂及数值模拟更加难以处理，因为它复杂多变的破碎模式与物理特征分布。各向异性弹性矩阵包括21个独立变量，导致实现快速可控的模拟充满挑战。更糟糕的是，各向异性网格的生成和分析相对更加复杂。为解决此问题，自适应的网格重构方法被用于模拟衣料等材料，但该方法仍存在较多限制与使用约束。而各向异性的粒子方法在搜索邻居时非常耗时，同时难以保持稳定的模拟。

发明內容

本发明的目的在于继承网格方法的物理准确性与粒子方法的快速有效，并采用自适应的几何细分产生局部高分辨率细节，提供了一种层次几何细分的各向异性材料脆性破裂模拟方法，包括以下步骤：

a)基于弹性力学与自适应SPH的弹力分析，具体包括：

i)经典线性弹性力学方程：

采用线性柯西格林形变张量计算应力ε：其中J为位移场的雅可比矩阵，I为单位矩阵；进而根据胡克定律σ＝Cε分析应变σ：其中弹性矩阵C是采用受杨氏模量E和泊松比v控制的各向同性矩阵；

ii)基于自适应SPH的数值求解：

采用各向异性核的自适应SPH方法求解粒子的位移梯度应力ε、应变σ，进而根据粒子体积V得到应变能与粒子间作用力

b)可控的裂纹计算模式，具体包括：

i)分组数量：为了产生集中于一点的若干裂纹，将邻居粒子分成用户预设的N组，每组各产生一条裂纹；

ii)各向异性核函数：自适应SPH中的各向异性核可以控制裂纹的方向，将核函数主轴方向的作用改变至1/3～3倍；

iii)接触面积：根据固体厚度B及碰撞面积S计算裂纹半径：其中e为能量释放率，U_a，U_a+Δa表示裂纹半径为a，a+Δa时的应变能；

c)基于自适应粒子的破碎模拟，具体包括：

i)基于Delaunay四面体剖分的物理粒子初始化；

ii)基于自适应细分粒子的破碎细节生成，以完全符合裂纹计算结果；

iii)基于边界粒子的快速碰撞检测，以减少检测粒子数量至1/4～1/3。

本发明的有益效果：

现有的固体脆性破碎模拟，或采用网格方法，计算复杂耗时，或采用无网格方法，模拟难以稳定和保证细节。而本发明充分借鉴了FEM方法精确的空间剖分以初始化及重构表面，并采用快速简便的粒子方法完成了可控的模拟框架。

本发明相比传统的SPH方法，计算更加精确，使得物体的运动趋势也更加合理。充分考虑各种可控因素，保证了效果的真实感。此外，同样适用于复杂场景与混合材质物体。

总之，应用本发明可以快速有效地模拟各向异性材料的脆性破碎。在用户交互的可控性、计算的快捷性和模拟结果的逼真性上，本发明的方法都有显著提高。

附图说明

图1为各向异性SPH的核函数形态及权值分布示意图；

图2为裂纹计算过程示意图；

图3为破碎模拟过程示意图；

图4为粒子细分算法示意图；

图5为额外链接说明示意图；

图6为宠物模型表面粒子与边界粒子示意图；

图7为原始方法与细分方法对比效果图；

图8为不同各向异性核函数的破碎效果图；

图9为玻璃在不同接触面积下破碎效果图；

图10为多个模型破碎效果图；

图11为复杂的宠物模型破碎效果图；

图12为混合材质物体破碎效果图。

具体实施方式

以下结合附图对本发明进行详细描述。

本发明包括以下具体步骤：

1)基于弹性力学与自适应SPH的弹力分析：

根据弹性力学连续方程分析碰撞过程中的物理属性；采用SPH方法进行各物理属性的计算。

2)可控的裂纹计算模式：综合考虑分组数量、各向异性核函数、接触面积等因素，根据应变分析结果计算裂纹情况。

3)基于自适应粒子的破碎模拟：采用自适应细分粒子模拟破碎结果，并渲染对应结果。

本发明所述的基于弹性力学与自适应SPH的弹力分析，其具体为：

(1)当一个固体发生形变后，它从初始位置x₀移动到当前位置x_t，对应两个位置的位移场即为u＝x_t-x₀，因此该位移场的雅可比矩阵为其中I为单位矩阵。使用线性柯西-格林应变张量，其应力ε可以表达如下：

$\mathrm{\ϵ}=\frac{1}{2}(J+J^T)-I---\left(1\right)$

假设该固体由线性弹性材料组成，应变σ符合胡克定律σ＝Cε。C是一个6×6的矩阵，对于各向同性材料，C依赖于两个独立变量：杨氏模量E和泊松比v；对于各向异性材料其独立变量数目为21。

为了解决上述连续方程，必须采用一个合适的数值方法来离散化模型。出于无网格方法的有效性和灵活性，采用SPH模型，将物体离散为粒子，并计算对应的位移场梯度属性。在光滑粒子流体动力学方法SPH(Smoothed Particle Hydrodynamics)中，粒子承载着各种物理量，包括质量、速度等，通过求解粒子组的动力学方程和跟踪每个粒子的运动轨道，求得整个系统的力学行为。在SPH方法中，粒子p_i的位移场梯度可以通过一个确定邻居粒子集合与核函数W(x_ij，h)＝W(x_i-x_j，h)近似为：

$\▿u_i=\underset{j}{\mathrm{\Σ}}{V}_j{u}_{\mathrm{ji}}\▿W(x_{\mathrm{ij}},h)---\left(2\right)$

其中，u_ji＝u_j-u_i是邻居粒子p_j与p_i的距离向量，V_i是p_i的体积，h是光滑半径。

(2)传统的SPH采用各向同性核函数，只能反映时间和空间上的内部变化，由于它不考虑不同方向的影响，会导致丢失某些方向上的邻居信息。为了解决这些问题，采用各向异性核(如椭球型核函数)的自适应SPH模型。该核函数由椭球各轴长与轴向确定(如图1)，它可以表示为一个3×3的二阶张量矩阵G。对应的各向异性核函数为：

$\begin{array}{c}W(r,G)=\frac{15}{\mathrm{\π}}\left|\right|G\left|\right|{(1-|\left|\mathrm{Gr}\right|\left|\right)}^3,\\ \▿W(r,G)=-\frac{45}{\mathrm{\π}}\left|\right|G\left|\right|\frac{{(1-|\left|\mathrm{Gr}\right|\left|\right)}^2}{\left|\right|\mathrm{Gr}\left|\right|},\end{array}---\left(3\right)$

张量G将距离向量r旋转和伸缩至Gr。值得注意的是，各向同性核函数可以视为G＝h^-1I时的一种特殊情况。

考虑到张量G的初始化，必须首先提及矩阵的SVD(Singular Value Decomposition)分解：

G＝RG_kR^T     (4)

其中R是旋转矩阵，它的列向量对应椭球的各轴方向，G_k是对角阵，G_k的对角元素表示各轴的长度。确定轴向和轴长之后即可初始化G。对于简单的物体可以直接指派，而涉及到复杂的几何体，基于加权的主成分分析方法WPCA(Weighted PCA)适用于分析局部张量。

如图1所示，当计算各向异性材料的属性时，采用各向同性核，p_i可能找不到足够的邻居，从而导致模拟崩溃或出现其它错误；反之各向异性核则可以保证模拟的稳定性。

基于自适应SPH，粒子p_i的可以重新定义为：

$\▿u_i=\underset{j}{\mathrm{\Σ}}{V}_j{u}_{\mathrm{ji}}\▿W(x_{\mathrm{ij}},G)---\left(5\right)$

弹性应变能U_i为：

$U_i=V_i\frac{1}{2}({\mathrm{\ϵ}}_i\·{\mathrm{\σ}}_i)---\left(6\right)$

假定粒子的全部体积内的应力和应变是常量，p_i受到粒子p_j的弹力f_ji为：

$f_{\mathrm{ji}}=-{\▿}_{u_j}{U}_i=-V_i(I+\▿u_i^T){\mathrm{\σ}}_i{d}_{\mathrm{ij}}---\left(7\right)$

其中d_ij可以近似计算为：

$d_{\mathrm{ij}}=V_j\▿W(x_{\mathrm{ij}},G)---\left(8\right)$

计算上述公式，即可完成碰撞过程的弹力分析，进而应用于计算裂纹。

本发明所述的可控的裂纹计算模式，其具体为：

采用郎肯的假说作为破碎的标准，即一个点的主应力如果超过材料的标准阈值，该点就会生成破碎，且破裂面方向等于主应力方向。简单的采用此模式只能生成一条裂纹，但是现实世界中的破碎通常会产生不同方向且交于一点的数条裂纹。

明显地，一个点的应力与它的所有邻居粒子的对应属性相关。如果分开计算各邻居集合中的位移场梯度及应力属性，会产生不同方向与大小的应力，即意味着不同的裂纹。这个属性可以利用于控制裂纹的数量。

裂纹的计算模式如图2所示。当碰撞发生在某点时，首先设定该点的核函数为最近粒子的各向异性核，并找到核函数范围内的所有邻居粒子。进而将所有邻居粒子分为不同组别，分组数量可以由参数设定，图2中分组数量为2。分组之后，每组的应力及应变能即可根据公式(1)及公式(6)计算，进而产生不同的裂纹。假定能量释放率e恒定，裂纹的传播半径正比于应变能：

$\mathrm{\Δa}=-\frac{U_{a+\mathrm{\Δa}}-U_a}{\mathrm{eBS}}---\left(9\right)$

其中U_a，U_a+Δa表示裂纹半径为a，a+Δa时的应变能，B是模型的厚度，S是接触面积。通常情况下，裂纹数量等于分组数量。假设应变能呈径向分布，引入内部半径r_inner，其r_inner范围内的点应变能大于阈值U_limit，使该范围内的点全部粉碎，以逼近真实场景。

本发明所述的基于自适应粒子的破碎模拟，其具体为：

(1)物理粒子初始化：采用基于FEM域分解，Delaunay驱动的四面体剖分来近似表达模型，以初始化物理粒子。首先将模型预处理为层状结构，进而使用Delaunay算法生成模型的四面体表达。在每个四面体的中心位置分布一个粒子，每个粒子的质量等于四面体的质量，粒子半径r取决于四面体体积V：每个粒子还存储了另外两种信息：初始的四面体顶点位置，与邻接粒子之间的链接关系。存储初始的四面体顶点用于保持与生成表面，以及计算内力。固体的表面可以通过邻接四面体之间的链接约束进行重建，以解决几何拓扑改变的问题。

(2)自适应细分粒子：需要一定的几何处理以得到前文中计算所得的裂纹模式。图3(b-c)解释了裂纹的基本生成方法。找到所有与分析结果相交的链接，移除这些链接，即可生成破碎片段。但是，在这种方法中模拟的细节与碎片的形状会受限于初始四面体的情况。为了保证物理真实性与可视效果，四面体的体积必须足够小，但这意味着极大的粒子数量，将会严重影响计算效率。所以采用自适应细分的思想，在裂纹处通过粒子细分产生局部高分辨率细节的同时保持快速的模拟框架。

粒子细分模式如图3(d)所示。首先标记接近裂纹表面的所有粒子，对每个标记粒子对应的四面体进行细分，生成新的小四面体(粒子)以逼近裂纹计算的结果。然后生成新粒子的链接，以保证正确的表面重建。在标准的数据结构中，一个粒子最多只能拥有4个邻居。此时需要一个列表来存储不同层次粒子间的额外链接。最后，移除对应的链接即可得到局部高细节的破碎片段，如图3(e)所示。

对于每一个粒子的细分，考虑四面体与平面相交的情况，四面体的4个顶点与平面有两种相对位置：有1个顶点在一侧，另外3个顶点在另一侧；各有2个顶点在两侧。为行文方便，在后文分别简称其为‘31’与‘2-2’模式。

算法过程如图4所示。借鉴于Marching Cubes算法，用符号‘+’与‘-’标记四面体各顶点与裂纹平面的相对位置，并在一个预先存储了16种可能性的数组中查询相交情况。如果所有的顶点都是同一符号，该粒子未与断裂平面相交。反之，则依据其相交情况进行如图4所示的细分。

对于‘3-1’模式，首先将三个交点与底面中心互连，这样形成了一个顶部及一个中心位置的四面体。接下来，将底面中心与底面各顶点相连，可以将剩余几何体分割为3个四棱锥。对每个四棱锥，取底面中心与各点相连，即可分为4个四面体。依上述方式，一个‘3-1’模式的粒子就被分成了14个新的粒子。

对于‘2-2’模式，平面将四面体分割为两个三棱柱。对于每个三棱柱，我们选择切面的中心，将它与其他顶点相连，这样就可以将之分割为2个四面体和2个三棱锥。三棱锥可以按‘3-1’模式的方法分解。依此，一个‘2-2’模式的粒子就被分成了20个新的粒子。

当‘3-1’粒子分解时，新生成的粒子会与原始的底面邻居产生不平衡的链接，如图5(a)所示，这种链接会被存储到额外链接列表中，并在移除链接时进行额外的判断。5(b)是两个‘3-1’粒子底面相连的情况，不会产生额外的链接。5(c)是底面邻居在后续的分解中属于‘2-2’模式的情况，该情况仍需进行额外的存储与判断。对于‘2-2’粒子，其所有邻居肯定与平面相交，不会产生额外的链接。

如果不对粒子的细分进行限制，粒子数目会产生爆炸式的增长，因为接近碰撞点的粒子会不断的分解，所有只对初始化的物理粒子进行上述处理。新粒子的体积和质量会被重新计算，速度则按照所属破碎片段设定。

(3)边界粒子：邻居数量小于4的粒子即为表面粒子，它们会有某些面作为物体的表面。通常直接采用表面粒子检测碰撞，但这种方法效率较低，因此引入边界粒子进行碰撞检测。

对于所有表面粒子，采用n个最近的邻接粒子，产生一个边界粒子，n是用户设定的采样比例，在实验中为4。边界粒子的质量是所有粒子的和，位置是平均位置。不可忽视的是，表面粒子数目不会总是n的整倍数，而且某些边角处会有一些孤立粒子，对于这些粒子，直接设定其为边界粒子。图6是一个宠物模型的边界粒子示意图。

本发明涉及的层次几何细分的各向异性材料脆性破裂模拟方法，效果展示说明如下：

图7中展示了瓷碗掉落的场景，左图采用基本方法模拟，由于粗糙的初始化，而产生了人工折线的裂纹痕迹；右图则为采用本发明的结果，可以明显看到该方法带来了更准确的裂纹细节和后续运动情况。

图8-10则展示了裂纹的相关控制手段。图8中两个瓷碗的核函数经过调整，分别将沿碗面的作用调整到到1/2和2倍，可以看到不同的核函数趋向影响了裂纹的延展方向趋势。图9中采用了形状一致、半径分别为3和5的球体与玻璃碰撞，导致了不同的破碎半径。值得注意的是，该示例中本发明在保证玻璃整体结构完整的前提下生成极微小的裂纹。图10是一个多个瓷碗破碎的复杂场景，其中每个物体的分组数量(1～3)，核函数影响(1/3～3)都不相同。

图11是陶瓷宠物模型被球击碎的模拟，其中的宠物模型初始含23664个粒子，并在后续的模拟中生成了34218个细分粒子。图12是一个混合材料的太极饰品模型，上方的材料为纵向易碎材料，下方的材料为横向易碎材料，经撞击后上方材料破碎，并传播至下方材料。该方法亦可扩展至多种材料的物体。

以上列举的仅是本发明的具体实施例。显然，本发明不限于以上实施例，还可以有许多变形。本领域的普通技术人员能从本发明公开的内容直接导出或联想到的所有变形，均应认为是本发明的保护范围。

Claims (1)

1.一种层次几何细分的各向异性材料脆性破裂模拟方法，其特征在于包括以下步骤：

a)基于弹性力学与自适应SPH的弹力分析，具体包括：

i)经典线性弹性力学方程：

采用线性柯西格林形变张量计算应力ε：其中J为位移场的雅可比矩阵，I为单位矩阵；进而根据胡克定律σ＝Cε分析应变σ：其中弹性矩阵C是采用受杨氏模量E和泊松比v控制的各向同性矩阵；

ii)基于自适应SPH的数值求解：

采用各向异性核的自适应SPH方法求解粒子的位移梯度应力ε、应变σ，进而根据粒子体积V得到应变能与粒子间作用力

b)可控的裂纹计算模式，具体包括：

i)分组数量：为了产生集中于一点的若干裂纹，将邻居粒子分成用户预设的N组，每组各产生一条裂纹；

ii)各向异性核函数：自适应SPH中的各向异性核能够控制裂纹的方向，将核函数主轴方向的作用改变至1/3～3倍；

iii)接触面积：根据固体厚度B及碰撞面积S计算裂纹半径：其中e为能量释放率，U_a，U_a+Δa表示裂纹半径为a，a+Δa时的应变能；

c)基于自适应粒子的破碎模拟，具体包括：

i)基于Delaunay四面体剖分的物理粒子初始化；

ii)基于自适应细分粒子的破碎细节生成，以完全符合裂纹计算结果；

iii)基于边界粒子的快速碰撞检测，以减少检测粒子数量至1/4～1/3。