CN108986220A - 一种加速有限元求解物体网格模型弹性变形的方法 - Google Patents

Abstract

本发明公开了一种加速有限元求解物体网格模型弹性变形的方法。针对非线性弹性物体模型，给定离散为由一系列单元和顶点构成的粗网格和细网格，用分段线性矩阵形函数作为位移插值函数用有限元法求解计算获得粗网格的位移场，使用线性插值函数作为位移插值函数用有限元法求解计算获得细网格的位移场；粗网格对应的有限元法中的形函数采用特殊设计的分段线性矩阵形函数，分段线性矩阵形函数为不连续的形函数。本发明用于模拟不均匀、各项异性的非线性弹性物体的形变，在粗网格上也使用非线性本构关系对异质弹性体进行预测性模拟，获得两到三个数量级的加速，对系统整体刚度的做出良好的逼近。

Description

一种加速有限元求解物体网格模型弹性变形的方法

技术领域

本发明涉及了一种计算机物理模型和物理模拟方法，主要涉及了一种加速有限元求解物体网格模型弹性变形的方法，主要针对于模拟不均匀和非线性弹性物体的变形。

背景技术

高效的模拟复杂的变形物体是计算机动画中的长期目标。随着更复杂、更庞大的仿真的需求增加，如何使用精细的数值方法以达到即使被仿真物体的结构复杂度不断提高而计算代价却不会大幅度提高成为技术难题。当人们使用足够精细的网格来解决微小尺度的非均匀性时才可以捕捉到物理对象正确的运动，但如果几何上十分复杂或者有一个非常不同的刚度矩阵也会导致十分大的求解时间。比如器官上的静脉结构或者超材料上的细微结构。

简单的忽略小尺度的细节会导致整个物体的动力学特性。粗糙的模拟甚至无法得到最简单的变形，有时使对象比现实中的变得更加僵硬。随着视觉和制造业的保真度要求日益广泛地使用非线性几何模型和非线性本构关系，结合效率和可扩展性的数值算法设计在图形和计算物理学中越来越迫切。

发明内容

针对背景技术的需求，本发明的目的在于提出了一种加速有限元求解物体网格模型弹性变形的方法，尤其是加速比较复杂的弹性物体，用于模拟不均匀、各项异性的非线性弹性物体的形变，是一种高效的、具有良好拓展性的方法。这种方法包括优化了一组不连续的、矩阵形式的形函数，使得即使在粗粒度的网格剖分上也可以使用非线性本构关系对异质弹性体进行预测性模拟。

为实现上述目的，本发明的采用的技术方案如下：

针对非线性弹性物体模型，给定离散为由一系列单元和顶点构成的粗网格和细网格，用分段线性矩阵形函数(形函数)作为位移插值函数用有限元法求解计算获得粗网格的位移场，使用线性插值函数作为位移插值函数用有限元法求解计算获得细网格的位移场；粗网格对应的有限元法中的形函数采用特殊设计的分段线性矩阵形函数，分段线性矩阵形函数为不连续的形函数，使得粗网格所得的结果与在细网格所得的结果近似，并且求解速度更快。

所述的粗网格的网格数相比细网格的网格数更少，顶点排布更稀疏，。细网格的网格数更多，顶点数更密。

两种网格描述的物体的几何形状完全相同，细网格所占的空间领域的并集为粗网格所占的空间领域。

在有限元求解中，非线性弹性物体模型上的一点(可为顶点或者顶点以外的任意点)的位移由该点所在粗网格单元的顶点通过分段线性矩阵形函数插值得到，位移场离散在每一个顶点上，由所有位移构成位移场；并且在分段线性矩阵形函数插值过程中，采用局部标架辅助计算保证插值对坐标变换的不变性。

所述的分段线性矩阵形函数为分段线性的矩阵形函数，同时为分段线性形式和矩阵形式。

本发明中的粗网格单元和细网格单元是指三维物体模型中的一个网格面片。

所述的分段线性矩阵形函数在粗网格单元内某一点的值等于一系列采样值的线性插值，表示如下：

式中，表示粗网格单元所占的空间区域,表示粗网格单元所占的空间区域内各个细网格单元所占的空间区域；X表示非线性弹性物体模型上一点的坐标，为一个向量；是点X所在粗网格单元上顶点i的形函数，作为分段线性矩阵形函数，i为点X所在粗网格单元上的一个顶点；分段线性矩阵形函数采用矩阵形式，矩阵维度为d×d，d为空间维度，二维情况下为2，三维情况下为3；n_ij表示形函数在所在粗网格单元内细网格单元顶点j处的采样值，为一个矩阵，也是需要求解的变量；为点X所在细网格单元上顶点j的插值函数，所有细网格单元的插值函数采用线性插值函数，j为点X所在细网格单元上的一个顶点。

每一个粗网格单元上的每一个顶点的形函数有一系列的采样值，采样值的位置为细网格单元上的各个顶点位置，且采样值的数量等于该粗网格单元细分后的细网格单元中的顶点总数。

针对所述采样值进行求解，具体是：使形函数所对应的位移场转化为线性二次规划问题，使满足几何条件、物理条件和正则化，对线性二次规划问题进行求解获得采样值，从而得到分段线性矩阵形函数；对于每一个粗网格单元上的形函数当满足条件时代入得到一个线形二次规划问题，线形二次规划问题公式如下：

式中，Ι为d×d的单位矩阵，表示分段线性矩阵形函数的导数，M表示单位四阶张量，为点X所在粗网格单元上顶点i的坐标，为点X所在粗网格单元上顶点k的坐标，表示点X所在细网格单元上顶点j的坐标，a、b为三维的笛卡尔坐标系下x、y、z三轴方向中的任意两个坐标方向，且可以取相同值，为坐标处ab方向的典型位移；tr()为迹函数，δ_ij为克罗内克函数，[·]_×来表示矩阵外积，当i＝k时，δ_ik为1，当i≠k时，δ_ik为0；

所述的典型位移计算是求解以下公式表示的带诺伊曼边界条件的线性化弹性静力学问题：

式中，表示散度计算，σ为应力张量，e_a为第a个坐标方向的单位向量，e_b为第b个坐标方向的单位向量，n为边界上的法向量，T表示转置，σ(h_ab)表示典型位移下的应力张量，h_ab表示ab方向的典型位移场，insideΩ表示在物体内，on表示在物体边界上；

求解上述偏微分微分方程组的公式时使用有限元法求解，并且将应变能的海森矩阵设为一个常量矩阵，等于弹性体处于静止状态时的海森矩阵，形函数采用线性插值函数，并且附加以下的第零阶和第一阶惯性矩约束：

式中，为变形后粗网格上的顶点参数坐标，为初始状态粗网格顶点坐标的平均值。

本发明中满足几何条件、物理条件和正则化具体如下：

A、所述的几何条件分为以下三条约束：

(1)当对一个弹性体作用一个平移时，该弹性体内部的应变场不会变化。为了满足平移不变性，任一粗网格单元上所有顶点的形函数在任意位置的和都为一个单位阵。用公式描述如下：

式中，Ι为d×d的单位矩阵。

(2)当对一个弹性体作用一个旋转时，该弹性体内部的应变场不会变化。为了满足旋转不变性，若弹性体以角速度ω旋转，任一粗网格单元所占空间中的一点在旋转后的位置应当等于粗网格单元上的顶点在旋转后的位置分别乘上其对应的形函数用公式描述如下：

为了方便便于描述，使用矩阵记号[·]×来表示矩阵外积，比如外积a×b表示为[b]×a。上述约束写为：

式中，为某一粗网格单元上第i个顶点的参数坐标。

(3)为了满足拉格朗日插值条件，任一粗网格单元上任一顶点上的形函数应当在该顶点的位置等于一个单位阵，在其他顶点上等于零矩阵。用公式描述如下：

式中，为点X所在粗网格单元上顶点k的坐标，δ_ik为克罗内克函数；

B、分段线性的矩阵形函数能够精确的重建一组典型位移，也就是某一个粗网格单元内的形变场等于该粗网格单元顶点上的典型位移取值分别乘上对应顶点的形函数

所述的物理条件公式描述如下：

式中，a、b为三维的笛卡尔坐标系下x、y、z三轴方向中的任意两个坐标方向，且可以取相同值，为坐标处ab方向的典型位移；

所述的典型位移的计算方法如下：

求解以下带诺伊曼边界条件的线性化的弹性静力学问题：

式中，表示散度计算，σ为应力张量，e_a为第a个坐标方向的单位向量，e_b为第b个坐标方向的单位向量，n为边界上的法向量，T表示转置。

求解上述偏微分微分方程组时使用现有有限元法求解，并且将应变能的海森矩阵设为一个常量矩阵，等于弹性体处于静止状态时的海森矩阵。形函数采用线性插值函数，并且附加第零阶和第一阶惯性矩约束：

式中为变形后粗网格上的顶点参数坐标，为初始状态粗网格顶点坐标的平均值。该物理条件在三维情况下提供6条约束方程(二维情况下3条约束方程)。满足此物理条件的形函数将精确代表典型位移的任意线性组合。

C、为了使得位移更平滑，任一粗网格单元内的所有顶点上的形函数的散度的Frobenius范数的和应当取最小值。所述的正则化使用以下公式，使得形函数应使下述泛函数取得最小值：

式中，tr()为迹函数，M为d×d×d×d的四阶单位张量，为形函数的导数，为一个d×d×d的三阶张量。

使用以下公式表示的带有局部标架的分段线性矩阵形函数插值获得位移场时，从而保证插值对坐标变换的不变性：

式中，u(X)表示物体上X处的位移场；R_e为局部标架；R_e等于X所在粗网格的变形梯度的极分解中的旋转矩阵，变性梯度中的x′为初始时刻坐标X在物体变形后的坐标，用线性插值函数得到，变性梯度等于坐标x′对坐标X求雅可比矩阵。

模拟一个物体在外力场作用和位置约束下的变形的大部分过程与常规的有限元分析求解方法一样，除了使用上述的分段线性的矩阵形函数替换线形插值函数计算变形梯度。

计算新的位移场下的新变形梯度时将每个粗网格单元都映射到一个等参元去计算：

式中，x表示物体初始时刻X处变形后的坐标，使用分选线形矩阵形函数插值，表示使用新的位移场下的新变形梯度，是变形后的位置x的雅可比矩阵，是位移场的雅可比矩阵；表示张量积：表示张量缩并；表示分段线性矩阵行函数对X求偏导数；为X所在粗网格单元的顶点j对应的线性插值函数，ξ为物体上X处映射到等参元后的坐标，为对ξ偏导数；-1表示矩阵求逆；R_e和的张量积会生成一个三阶张量，再与缩并。

本发明是在粗网格上某个单元上某个顶点的形函数在其所在粗网格单元上为一个分段的函数矩阵。单元的领域可以由一系列的细网格领域划分。划分后该粗网格单元内部(包括该粗网格单元上)的一系列细网格单元的顶点上形函数的采样值n_ij为待求解值(自由度)。为了求解每个粗网格单元上每个顶点的一系列采样值n_ij，使该形函数所刻画的位移场满足几何条件、物理条件和正则化，从而转化为线性二次规划问题，可有下述方法求解。

本发明提出的分段线性矩阵形函数在粗网格单元上将通过求解一个线性二次规划问题得到，其中一些识别出来的关键条件能够很好地平衡单元间连续性和单元内刚度，从而对系统整体刚度的做出良好的逼近。

本发明与背景技术相比具有的有益效果是：

1)形函数的矩阵形式。本发明中的形函数使用矩阵的形式而不是标量，从而可以非常好的利用复杂物体的不可避免的各向异性行为。这种插值方式可提供更多的自由度以耦合不同维度的线性变化并且自然地处理各项异性。

2)形函数的分段线性形式。本发明中的形函数采用分段线性函数，从而可以更好的近似单元内部的位移场。并且在约束中加入了调和位移，使得形函数可以表达调和位移的任何线性组合，从而抑制了使用局部标架Re可能造成的几何不连续性。

3)求解时间。由于本发明中在粗网格上的求解大大减少了顶点的数量从而降低了自由度，虽然在求解调和位移和形函数时会有额外的时间开销，但是本发明中的方法依然可以大量减少求解时间并且保证结果的误差不大，特别是在动力学问题中效果显著。和现有有限元计算相比，本发明方法可获得两到三个数量级的加速。

4)材料非线性与不均匀性和边界非线性。在有限元求解中使用本发明中的形函数可以更好的处理非线性本构关系的材料，不均质的材料和碰撞等动力学问题。

附图说明

图1是本发明分段线性形函数在二维平面四边形单元上的部分函数图像。

图2是本发明求解有限元问题的算法流程图。

图3是本发明中求解得到的典型位移的一个算例结果图。

图4是本发明中的方法求解一个动力学问题的算例与现有方法的对比图。

图5是本发明中的方法与Nesme方法求解非均质的一根棒和一个柱子在不同的受力状态下的变形的一个对比图。

图6是本发明中的方法与Kharevych方法求解非均质的一根棒和一个柱子在不同的受力状态下的变形的一个对比图。

具体实施方式

下面结合附图和实施例对本发明作进一步说明。

如图1所示，按照本发明完整方法实施的实施例如下：

1)给定离散化一个弹性体后的粗、细网格(包括单元顶点列连接关系和初始时刻的顶点坐标)，并已知每个细网格单元的材质特性(一般包括泊松比和杨氏模量，也可为拉梅系数)、约束、受力。

2)在细网格上求解典型位移h_ab。求解以下带诺伊曼边界条件的线性化的弹性静力学问题：

式中，表示散度计算，σ为应力张量，e_a为第a个坐标方向的单位向量，e_b为第b个坐标方向的单位向量，n为边界上的法向量，T表示转置。

求解上述偏微分微分方程组时使用有限元法求解，使用现有线性插值函数插值位移场，并且附加第零阶和第一阶惯性矩约束：

式中为变形后粗网格上的顶点参数坐标，为初始状态粗网格顶点坐标的平均值。该物理条件在三维情况下提供6条约束方程(二维情况下3条约束方程)。满足此物理条件的形函数将精确代表典型位移的任意线性组合。

实施例求解得到的典型位移例子结果如图3所示，对一个三维立方体求它的典型位移，得到如图的六组结果，a、b等于x、y、z中的任意两个。第一排从左到右为典型位移h_xx,h_xy,h_xz；第二排从左到右为h_yy,h_yz,h_zz。

3)求粗网格上每个单元每个顶点的形函数。求解以下带约束的最优化问题：

求解这个最优化问题时先将一给定单元所有的采样值集成为一个向量n_e后，这些约束可以表达成一个线性系统：

Cn_e＝y

其中，C为一个表达上述约束的矩阵，n_e为求解的采样值，y为上述约束方程中右手项集成的向量。求解该线性系统时，将n_e分解为：

式中，U为矩阵C的零空间，q为零空间的坐标值，为一个任意解。所述零空间和一个任意解使用Suitesparse的QR分解得到。

4)由此形成分段线性形函数，在二维平面四边形单元上如图1所示，左图表示一个粗网格单元并且标注了其中的一个顶点。中间图为该粗网格单元上所对应的一系列细网格单元并且标注了一个细网格单元的顶点。右图为中间图所标注的细网格单元顶点的线性插值函数，高度为函数值。该线性插值函数乘上求得的采样值n_ij，对所有该粗网格内所有细网格结点重复上述步骤并且求和，就得到了粗网格单元上顶点i的分段线性形函数。

用分段线性矩阵形函数作为位移插值函数用有限元法集成能量求解计算获得粗网格的位移场，得到某一时刻下弹性物体的位置x_i。即根据形函数计算新的位移场下的新变形梯度，然后计算弹性势能和其他势能，附加约束，通过求导得到能量梯度和海森矩阵。最后使用隐式欧拉等方法计算每一个时间步长后的粗网格上每个顶点的位置x_i。

根据得到的每一时刻的位置x_i进行渲染呈现结果。

具体实施中按照本发明方法求解一个动力学问题的算例与现有方法对比如图4所示。该问题为模拟一个空腔心脏模型从空中落到地上的动态过程。将细网格结果作为真实状况作为参照。左边一列为使用线性插值函数在细网格上的求解结果，中间一列为使用本发明的方法在粗网格上的求解结果，右边一列为使用线性插值函数在粗网格上的求解结果。第一排为初始时刻，第二排为第20帧时的状态，第三排为第27帧时的状态。可以看到，第二列的状态基本与第一列相同，并且比第三列相比求解更快。

具体实施中按照本发明方法与现有Kharevych方法均求解非均质的一根棒和一个柱子在不同的受力状态下的变形情况，对比如图5所示。将细网格结果作为真实状况作为参照。左边一列为使用线性插值函数在细网格上的求解结果，中间一列为使用本发明的方法在粗网格上的求解结果，右边一列为使用Nesme的方法在粗网格上的求解结果。第一排为一根棒左端固定，右端在重力的作用下下坠，第二排为在一根棒的右端固定，在左端施加一个扭矩，棒子扭转的结果。第三排为一个柱子上端固定，在重力作用下下坠。第四排为一个柱子上端固定，施加一个弯矩，柱子向上弯的结果。可以看到本发明的方法的结果与参照的真实状况基本一致。

具体实施中按照本发明方法与现有Kharevych方法均求解非均质的一根棒和一个柱子在不同的受力状态下的变形情况，对比如图6所示。将细网格结果作为真实状况作为参照。左边三列中，从左到右分别为使用线性插值函数在细网格上的求解结果、使用本发明的方法在粗网格上的求解结果和使用Kharevych的方法在粗网格上的求解结果。其中第一排为一根棒左端固定，右端在重力的作用下下坠，第二排为在一根棒的右端固定，在左端施加一个扭矩，棒子扭转的结果。右边三列同理。其中第一排为一根棒左端固定，右端在重力的作用下下坠，第二排为在一根棒的右端固定，在左端施加一个扭矩，棒子扭转的结果。可以看到，方法与参照的真实状况基本一致，能够在细节上表现的更好，并且求解速度更快。

Claims (6)

1.一种加速有限元求解物体网格模型弹性变形的方法，其特征在于：针对非线性弹性物体模型，给定离散为由一系列单元和顶点构成的粗网格和细网格，用分段线性矩阵形函数作为位移插值函数用有限元法求解计算获得粗网格变形后的位移场，使用线性插值函数作为位移插值函数用有限元法求解计算获得细网格变形后的位移场；粗网格对应的有限元法中的形函数采用特殊设计的分段线性矩阵形函数，分段线性矩阵形函数为不连续的形函数。

2.根据权利要求1所述的一种加速有限元求解物体网格模型弹性变形的方法，其特征在于：在有限元求解中，非线性弹性物体模型上的一点的位移由该点所在粗网格单元的顶点通过分段线性矩阵形函数插值得到，位移场离散在每一个顶点上，由所有位移构成位移场；并且在分段线性矩阵形函数插值过程中，采用局部标架辅助计算保证插值对坐标变换的不变性。

3.根据权利要求1所述的一种加速有限元求解物体网格模型弹性变形的方法，其特征在于：所述的分段线性矩阵形函数在粗网格单元内某一点的值等于一系列采样值的线性插值，表示如下：

式中，表示粗网格单元所占的空间区域,表示粗网格单元所占的空间区域内各个细网格单元所占的空间区域；X表示非线性弹性物体模型上一点的坐标，为一个向量；是点X所在粗网格单元上顶点i的形函数，作为分段线性矩阵形函数，i为点X所在粗网格单元上的一个顶点；分段线性矩阵形函数采用矩阵形式，矩阵维度为d×d；n_ij表示形函数在所在粗网格单元内细网格单元顶点j处的采样值，为一个矩阵；为点X所在细网格单元上顶点j的插值函数，所有细网格单元的插值函数采用线性插值函数，j为点X所在细网格单元上的一个顶点。

4.根据权利要求3所述的一种加速有限元求解物体网格模型弹性变形的方法，其特征在于：采样值的位置为细网格单元上的各个顶点位置，且采样值的数量等于该粗网格单元细分后的细网格单元中的顶点总数。

5.根据权利要求3所述的一种加速有限元求解物体网格模型弹性变形的方法，其特征在于：针对所述采样值进行求解，具体是：使形函数所对应的位移场转化为线性二次规划问题，使满足几何条件、物理条件和正则化，对线性二次规划问题进行求解获得采样值，从而得到分段线性矩阵形函数；线形二次规划问题公式如下：

式中，Ι为d×d的单位矩阵，表示分段线性矩阵形函数的导数，M表示单位四阶张量，为点X所在粗网格单元上顶点i的坐标，为点X所在粗网格单元上顶点k的坐标，表示点X所在细网格单元上顶点j的坐标，a、b为三维的笛卡尔坐标系下x、y、z三轴方向中的任意两个坐标方向，且可以取相同值，为坐标处ab方向的典型位移；tr()为迹函数，[·]_×来表示矩阵外积，当i＝k时，δ_ik为1，当i≠k时，δ_ik为0；

上述公式中的典型位移计算是求解以下公式表示的带诺伊曼边界条件的线性化弹性静力学问题：

式中，表示散度计算，σ为应力张量，e_a为第a个坐标方向的单位向量，e_b为第b个坐标方向的单位向量，n为边界上的法向量，T表示转置，σ(h_ab)表示典型位移下的应力张量，h_ab表示ab方向的典型位移场，insideΩ表示在物体内，表示在物体边界上；

求解上述公式时使用有限元法求解，并且将应变能的海森矩阵设为一个常量矩阵，等于弹性体处于静止状态时的海森矩阵，形函数采用线性插值函数，并且附加以下的第零阶和第一阶惯性矩约束：

式中，为变形后粗网格上的顶点参数坐标，为初始状态粗网格顶点坐标的平均值。

6.根据权利要求3所述的一种加速有限元求解物体网格模型弹性变形的方法，其特征在于：使用以下公式表示的带有局部标架的分段线性矩阵形函数插值获得位移场时：

式中，u(X)表示物体上X处的位移场；R_e为局部标架；R_e等于X所在粗网格的变形梯度的极分解中的旋转矩阵。