基于有限元分析的非局部超弹性与接触建模方法及系统

本发明属于计算力学、近场动力学与物理仿真，具体涉及一种基于有限元分析的非局部超弹性与接触建模方法。本发明以近场动力学的非局部点对相互作用为建模框架，有限元分析用于获得四面体离散下的几何量及超弹性力学量。

背景技术：

1、模拟超弹性形变体（如橡胶、生物软组织、布料等）广泛应用于电影特效、视频游戏、虚拟手术及机器人仿真中。传统的有限元方法（finite element method, fem）是描述此类连续介质形变的标准方法。在经典连续介质力学中，本构关系通常表示为位置及其空间导数的点函数，并通过求解偏微分方程（pde）来计算材料响应。

2、然而，传统fem在处理大变形及实时仿真时面临巨大挑战：

3、一、计算复杂度高：为了求解非线性系统，通常采用牛顿-拉夫逊（newton-raphson）方法，这需要计算能量密度的二阶导数（即海森矩阵，hessian matrix）。海森矩阵的组装和求解（通常涉及稀疏线性系统求解）极其耗时，且占用大量内存。

4、二、gpu并行困难：传统隐式求解器涉及全局矩阵的构建与求解，难以充分利用现代gpu的众核并行架构。虽然已有基于对角化海森矩阵或局部求解（如projectivedynamics）的加速方法，但它们往往牺牲了物理准确性或收敛速度。

5、三、不连续性处理困难：基于局部微分算子的fem在处理断裂、接触等不连续现象时，空间导数定义变得困难，通常需要复杂的额外约束处理。

6、近场动力学（peridynamics）属于非局部连续介质力学方法，其以一定邻域内材料点之间的积分型点对相互作用替代局部空间微分算子，因而适于描述断裂、接触等不连续现象。现有近场动力学模型在连续超弹性响应的离散精度、几何权重构造以及接触与摩擦的统一求解方面仍有改进空间。有限元分析能够提供成熟的四面体离散、形函数梯度及超弹性应力计算手段，但若直接沿用传统有限元的局部单元组装及海森矩阵求解方式，仍存在计算与存储开销较大、gpu并行效率受限等问题。因此，亟需一种以近场动力学为框架、借助有限元分析获得离散几何及一阶力学量，并将其重构为非局部点对作用的高效建模方法。

技术实现思路

1、本发明的目的在于提出一种基于有限元分析的非局部超弹性与接触建模方法。该方法以近场动力学的非局部点对相互作用为主体框架，借助有限元分析手段获得四面体离散下的几何量、变形梯度及超弹性应力，并将其重构为非局部点对作用，以解决非局部超弹性与接触统一建模中几何权重构造、海森矩阵计算复杂及并行求解效率受限的技术问题。

2、为实现上述目的，本发明采用以下技术方案。

3、一种基于有限元分析的非局部超弹性与接触建模方法，包括以下步骤：

4、采用有限元分析手段对仿真对象进行四面体离散，并根据各四面体单元的节点位置和变形状态确定单元几何信息、变形梯度及超弹性节点力，作为非局部建模的输入；

5、在近场动力学非局部框架下，根据各四面体单元的边向量和形状张量建立形函数梯度与所述边向量的对应关系，并基于格拉姆矩阵确定各节点对的几何权重；

6、根据所述对应关系、所述几何权重和所述变形梯度，将采用有限元分析获得的所述超弹性节点力重构为沿所述边向量作用的非局部点对作用力，得到各节点对之间的点对刚度张量；

7、根据所述非局部点对作用力和所述点对刚度张量，组装具有拉普拉斯矩阵结构的全局运动方程；

8、将所述点对刚度张量分解为正定分量和负定分量，根据所述正定分量构造迭代系数矩阵，并对所述全局运动方程进行半隐式迭代求解，得到各节点的更新位置。

9、进一步地，采用有限元分析手段对仿真对象进行四面体离散，并根据各四面体单元的节点位置和变形状态确定单元几何信息、变形梯度及超弹性节点力，作为非局部建模的输入，包括：

10、采用有限元分析手段对所述仿真对象进行四面体离散，得到多个四面体单元；

11、根据所述多个四面体单元的节点位置和变形状态，确定所述单元几何信息、所述变形梯度及所述超弹性节点力，作为非局部建模的输入。

12、进一步地，根据所述多个四面体单元的节点位置和变形状态，确定所述单元几何信息、所述变形梯度及所述超弹性节点力，作为非局部建模的输入，包括：

13、根据所述多个四面体单元中各四面体单元的节点位置，确定包括边向量、单元体积和形函数梯度的所述单元几何信息；

14、根据所述多个四面体单元中各四面体单元的变形状态，确定所述变形梯度和第一皮奥拉-克希霍夫应力张量；

15、根据所述单元体积、所述形函数梯度和所述第一皮奥拉-克希霍夫应力张量，确定所述超弹性节点力，并将所述单元几何信息、所述变形梯度及所述超弹性节点力作为非局部建模的输入。

16、进一步地，在近场动力学非局部框架下，根据各四面体单元的边向量和形状张量建立形函数梯度与所述边向量的对应关系，并基于格拉姆矩阵确定各节点对的几何权重，包括：

17、根据各四面体单元的所述边向量确定对应的形状张量；

18、根据所述形状张量建立所述形函数梯度与所述边向量的对应关系；

19、根据所述边向量的内积构建格拉姆矩阵，并基于所述格拉姆矩阵确定各节点对的几何权重。

20、进一步地，根据所述边向量的内积构建格拉姆矩阵，并基于所述格拉姆矩阵确定各节点对的几何权重，包括：

21、从各四面体单元的顶点中分别确定参考顶点，并根据由所述参考顶点指向其余顶点的所述边向量构建参考格拉姆矩阵；

22、对所述参考格拉姆矩阵求逆，并基于得到的逆矩阵确定所述参考顶点与所述其余顶点之间的几何权重；

23、根据所述逆矩阵的矩阵元素和对应四面体单元的顶点仿射关系，确定其余节点对的几何权重，得到各节点对的所述几何权重。

24、进一步地，根据所述对应关系、所述几何权重和所述变形梯度，将采用有限元分析获得的所述超弹性节点力重构为沿所述边向量作用的非局部点对作用力，得到各节点对之间的点对刚度张量，包括：

25、根据所述对应关系将所述超弹性节点力分解至各节点对，并根据各四面体单元的单元体积、所述变形梯度、第一皮奥拉-克希霍夫应力张量和所述几何权重，确定各四面体单元对应的单元点对刚度张量；

26、汇总包含同一节点对的各四面体单元对应的所述单元点对刚度张量，得到各节点对之间的所述点对刚度张量；

27、根据所述点对刚度张量和变形后的节点位置差，确定沿所述边向量作用的所述非局部点对作用力。

28、进一步地，根据所述对应关系将所述超弹性节点力分解至各节点对，并根据各四面体单元的单元体积、所述变形梯度、第一皮奥拉-克希霍夫应力张量和所述几何权重，确定各四面体单元对应的单元点对刚度张量，包括：

29、将各四面体单元的所述第一皮奥拉-克希霍夫应力张量与对应的所述变形梯度的逆张量进行组合，得到应力映射张量；

30、根据所述对应关系将所述超弹性节点力分解至各节点对，并根据所述应力映射张量、所述单元体积和所述几何权重，确定所述单元点对刚度张量。

31、进一步地，根据所述非局部点对作用力和所述点对刚度张量，组装具有拉普拉斯矩阵结构的全局运动方程，包括：

32、对所述仿真对象进行碰撞检测，确定碰撞点对，并根据所述碰撞点对构建虚拟边界单元；

33、根据所述虚拟边界单元的变形状态确定各向异性势能，并根据所述各向异性势能确定碰撞点对作用力和碰撞点对刚度张量；

34、将所述碰撞点对作用力和所述碰撞点对刚度张量分别与所述非局部点对作用力和所述点对刚度张量进行组装，得到具有拉普拉斯矩阵结构的所述全局运动方程。

35、进一步地，根据所述虚拟边界单元的变形状态确定各向异性势能，并根据所述各向异性势能确定碰撞点对作用力和碰撞点对刚度张量，包括：

36、根据所述虚拟边界单元的变形状态和所述碰撞点对的接触法向，确定法向变形分量、切向变形分量和共旋变形状态；

37、根据所述法向变形分量确定法向势能，并根据所述切向变形分量和所述共旋变形状态确定切向势能；

38、将所述法向势能和所述切向势能组合为所述各向异性势能，并根据所述各向异性势能确定所述碰撞点对作用力和所述碰撞点对刚度张量。

39、进一步地，将所述碰撞点对作用力和所述碰撞点对刚度张量分别与所述非局部点对作用力和所述点对刚度张量进行组装，得到具有拉普拉斯矩阵结构的所述全局运动方程，包括：

40、分别汇总作用于同一节点对的所述非局部点对作用力与所述碰撞点对作用力，以及对应的所述点对刚度张量与所述碰撞点对刚度张量，得到总点对作用力和总点对刚度张量；

41、根据所述总点对作用力对应的节点对关系，将所述总点对刚度张量确定为全局刚度矩阵的非对角元素，并将每一节点对应的非对角元素之和的负值确定为对应的对角元素；

42、根据各节点的当前位置和当前速度确定预测位置；

43、根据节点质量确定质量矩阵，并根据所述总点对作用力、所述质量矩阵、所述全局刚度矩阵、时间步长和所述预测位置构建所述全局运动方程。

44、进一步地，将所述点对刚度张量分解为正定分量和负定分量，根据所述正定分量构造迭代系数矩阵，并对所述全局运动方程进行半隐式迭代求解，得到各节点的更新位置，包括：

45、将各节点对的所述点对刚度张量分解为正定分量和负定分量；

46、根据所述正定分量汇总各节点对应的对角系数，并结合节点质量和时间步长构造所述迭代系数矩阵；

47、根据所述迭代系数矩阵、所述负定分量、所述全局运动方程中的预测位置和上一迭代的节点位置，对所述全局运动方程进行半隐式不动点迭代，得到各节点的所述更新位置。

48、进一步地，根据所述迭代系数矩阵、所述负定分量、所述全局运动方程中的预测位置和上一迭代的节点位置，对所述全局运动方程进行半隐式不动点迭代，得到各节点的所述更新位置，包括：

49、根据所述负定分量、所述预测位置和所述上一迭代的节点位置，确定各节点的迭代右端项；

50、采用雅可比迭代方式，根据所述迭代系数矩阵和所述迭代右端项分别求解各节点的当前迭代位置；

51、将所述当前迭代位置作为下一迭代的节点位置进行重复求解，直至满足收敛条件，得到各节点的所述更新位置。

52、一种基于有限元分析的非局部超弹性与接触建模系统，包括：

53、有限元分析模块，用于采用有限元分析手段对仿真对象进行四面体离散，并根据各四面体单元的节点位置和变形状态确定单元几何信息、变形梯度及超弹性节点力，作为非局部建模的输入；

54、几何权重确定模块，用于在近场动力学非局部框架下，根据各四面体单元的边向量和形状张量建立形函数梯度与所述边向量的对应关系，并基于格拉姆矩阵确定各节点对的几何权重；

55、点对力转换模块，用于根据所述对应关系、所述几何权重和所述变形梯度，将采用有限元分析获得的所述超弹性节点力重构为沿所述边向量作用的非局部点对作用力，得到各节点对之间的点对刚度张量；

56、运动方程组装模块，用于根据所述非局部点对作用力和所述点对刚度张量，组装具有拉普拉斯矩阵结构的全局运动方程；

57、迭代求解模块，用于将所述点对刚度张量分解为正定分量和负定分量，根据所述正定分量构造迭代系数矩阵，并对所述全局运动方程进行半隐式迭代求解，得到各节点的更新位置。

58、本发明取得了如下有益效果。

59、1.本发明以近场动力学的非局部点对相互作用为主体框架，通过有限元分析获得四面体离散下的形函数梯度及超弹性节点力，再利用形状张量和格拉姆矩阵将其重构为沿边方向作用的非局部点对作用力。

60、2.本发明利用四面体单元的顶点仿射关系，根据一个参考顶点对应的格拉姆矩阵逆矩阵确定各节点对的几何权重，每个四面体单元仅需执行一次矩阵求逆，减少了几何权重预计算过程中的重复求逆运算，提高了初始化效率。

61、3.本发明根据第一皮奥拉-克希霍夫应力张量、变形梯度、单元体积和几何权重构造点对刚度张量，使点对作用力仅依赖一阶力学量，无需计算应力张量相对于变形梯度的导数，避免了海森矩阵及四阶张量计算，降低了单步迭代的计算量、显存占用和寄存器需求，并提高了gpu并行度和计算吞吐量。

62、4.本发明将点对刚度张量组装为具有拉普拉斯矩阵结构的全局运动方程，并通过正定分量与负定分量分离，以正定分量构造半隐式迭代系数矩阵，能够提高非线性运动方程迭代求解的收敛性和数值稳定性。

63、5.本发明针对碰撞点对构建虚拟边界单元，基于法向变形分量和共旋变形状态的切向分量构建各向异性势能，将接触作用力和摩擦作用力分别限制于接触法向和接触切向，实现碰撞力与弹性力的统一点对建模，并避免传统惩罚方式产生的伪体积膨胀。