CN109002630A - 一种超弹性材料的快速仿真方法 - Google Patents

Abstract

本发明公开一种超弹性材料的快速仿真方法。本方法为：1)对待仿真的三维模型进行均匀离散，并计算每个离散点的状态信息；2)设定离散点的邻域半径δ，构建每一离散点的邻域连接关系；3)根据待仿真的超弹性材料设置非线性函数、各向异性函数以及相关材料参数，确定三维模型的材质本构模型；4)依据材质本构模型和离散点存储的当前状态计算形变内能和对应产生的内力并施加在离散点上；5)对施加内力的离散点施加外力、摩擦力和弯曲防翻转力；然后求解下一个时刻离散点的位置和速度；6)对处理后的离散点进行碰撞检测与碰撞响应，调整离散点的位置和速度；7)重复4)～6)，直到仿真结束。本发明可有效地仿真非线性、各向异性现象。

Description

一种超弹性材料的快速仿真方法

技术领域

本发明属于计算机图形学领域，具体涉及一种基于近场动力学的超弹性材料建模与快速仿真方法。

背景技术

在计算机图形学中，复杂三维形变体仿真通常使用FEM方法。基于传统FEM方法，也已经存在了大量工作解决各向异性[Li,Y.(2014).Stable orthotropic materials.ACMSiggraph/eurographics Symposium on Computer Animation(pp.41-46).EurographicsAssociation.]、非线性材质的仿真模拟[Xu,H.,Sin,F.,&Zhu,Y.(2015).Nonlinearmaterial design using principal stretches.Acm Transactions on Graphics,34(4),75.]。但因为FEM方法在处理非连续性问题时存在天然的劣势，导致其无法很好地仿真碎裂等现象，而近年来提出的近场动力学模型[Silling,S.A.,Epton,M.,Weckner,O.,Xu,J.,&Askari,E.(2007).Peridynamic states and constitutive modeling.Journal ofElasticity,88(2),151-184.]则是利用积分形式计算巧妙处理非连续边界问题并且能够有效地做到不同维度仿真算法的统一，所以基于近场动力学的研究在工程计算和计算机图形学领域都开始逐年增加。

尽管如此，近场动力学和传统FEM方法的理论等价性还是十分模糊，而且已有的近场动力学模型只能处理简单的线性模型，无法对更为复杂的超弹性材料进行仿真。

发明内容

本发明主要针对传统近场动力学模型无法有效建模和仿真超弹性材料的问题，提出一种基于近场动力学理论的超弹性材料建模与快速仿真方法。其通过将形变内能分解成各向同性的静压力内能和剩余的各向异性内能两个独立成分并通过非线性函数和各向异性函数乘子来分别控制非线性现象和各向异性现象，并且通过增加弯曲防翻转能量有效处理了仿真过程中可能的翻转问题。

具体来说，本发明的技术方案如下：

一种基于近场动力学的超弹性材料建模与快速仿真方法，其具体步骤包括：

(1)对输入待仿真的三维模型进行均匀离散，计算每个离散点的体积、质量等基本物理量；

(2)设定离散点的邻域半径δ，构建离散点的邻域连接关系。每个离散点应与其周围任意一个距离小于邻域半径δ的离散点建立键连接。

(3)根据待仿真的超弹性材料设置非线性函数、各向异性函数以及其他材料参数如杨氏模量、泊松比确定待仿真的三维模型的材质本构模型；

(4)依据步骤(3)设置的材质本构模型，对离散点存储的当前状态(离散点的体积、位置两个物理量)计算形变内能和对应产生的内力并施加在离散点上；

(5)针对步骤(4)得到的中间结果施加外力和摩擦力；

(6)针对步骤(5)得到的中间结果施加弯曲防翻转力；

(7)利用步骤(6)得到结果，利用数值积分求解下一个时刻离散点的位置和速度；

(8)对更新后的离散点进行碰撞检测与碰撞响应，调整离散点的位置和速度；

(9)重复步骤(4)到(8)之间的操作，直到仿真结束。

进一步的，步骤(3)中的动力学建模提出了仿真方法成立的三个假设条件：一、每一个形变体微元所产生的形变内能或内力与其周围一定范围内的其余微元有关；二、形变体微元的形变内能可以分为各向同性的静压内能和剩余的各向异性内能；三、静压内能只和邻域体积变化相关，各向异性内能只和邻域内各键长有关。

进一步的，(3)中材质本构模型可以用如下能量方程进行描述：

其中W是微元(即离散点)形变内能密度，θ表示微元x所处邻域体积变化，是表示微元x邻域内其它点与中心点x所连键的长度变化，a,b为常量系数其值可以由杨氏模量泊松比确定，是以微元x为中心点的邻域，w表示邻域内所有离散点的权值，G表示邻域内所有离散点的各向异性量，x代表任意一个微元(即离散点)，w、G均是定义在上的场，w<ξ>、G<ξ>代表w、G在离散点x处的场在离散点ξ处的取值，g、h为非线性函数，dV_ξ为微元x的体积。如公式所描述的，形变内能由只与邻域内体积变化有关的各向同性能量ag(θ)和与邻域内各键长独立相关的各向异性能量组成。g、h可以是任意的一元非线性函数，只要其满足一定的条件在自变量取值等于1时取到最小值0(即min(g(x))＝g(1)，min(h(x))＝h(1)；1代表未发生形变，所以在1时应该取到0值，代表形变内能为0)，二阶导数始终大于等于0等。G<ξ>可以是任意三元函数，输入是一个三元向量代表的单位向量，输出是一个正实数。它表示各向异性能量虽然是由邻域内各个键的能量累加得到的，但不是平均累加而是要乘上一个跟键方向有关的权重值进行加权累加。当然G<ξ>也得满足一定条件，比如方向对偶性即方向完全相反的单位向量应该产生同样的权重值。w<ξ>是一个只和初始键长有关的权重函数，因为每一个微元形变内能相当于是离散存储在与邻域内其他微元间组成的键当中。那么键能量除了之前所讨论的可能和键方向有关，还有可能和初始键长即微元间的初始间距有关。w<ξ>可以是任意的一元函数，只要其满足值域为正实数并且是经过归一化处理。

进一步的，步骤(4)中计算总体形变内能只需要利用步骤(3)中确定的单个微元的内能求解方式对所有微元进行累加计算即可。至于内力的计算则是通过对总内能函数对微元位置向量求导即可得到。

进一步的，步骤(6)中弯曲防翻转力的施加可以进一步详细解释为：

弯曲防翻转力的计算关键在于如何衡量邻域内某个键的弯曲状态量以及是否发生了翻转。本发明中采用的方法是利用邻域内的点的初始坐标和形变后坐标估计一个仿射变换矩阵 就代表这个邻域总体的综合形变情况。的行列式如果是负数则代表邻域发生了翻转，此时我们需要将的前两行交换位置从而使其行列式重新变为正数。作用在初始键向量上可以得到估算的形变后的未弯曲键向量，利用形变后键向量减去估算的形变后的未弯曲键向量所得的值就是键的弯曲量。然后根据弯曲量就可以构建基于其的弯曲能量，从而计算出弯曲防翻转力的大小。

进一步的，步骤(7)中数值积分方法可以采用现有的多种方法，并不受本发明中所使用进行限制。本发明中采用了拟牛顿法进行计算，将下一个时间步粒子位置向量的求解重构成一个能量最小化问题然后通过估算一个近似Hessian矩阵H并对其进行预分解，利用能量的负梯度值除以H得到下降步长Δqⁿ⁺¹。由于H是一个近似Hessian矩阵并不是准确值，所以按照此时的Δqⁿ⁺¹进行迭代并不能保证一定收敛。因此使用后向缩减线性搜索(backtracking search)方法确定一个安全的能量下降步长来更新qⁿ⁺¹←qⁿ⁺¹+αΔqⁿ⁺¹，α的值从1开始，每次减半直到某一个值使得更新后qⁿ⁺¹使得能量E(qⁿ⁺¹)下降。

至此，完成了整个方法，与现有技术相比，本发明具有如下的优点：

(1)本发明提出的基于近场动力学的材质本构模型可以有效地仿真非线性、各向异性现象；

(2)本发明中提出的弯曲防翻转力可以有效的处理翻转问题，让物体恢复原状；

(3)本发明采用的拟牛顿法求解器，通过预分解近似Hessian矩阵可以有效地加速仿真计算；

(4)形变的仿真具有可控性。

附图说明

图1为仿真算法流程图。

图2为离散化示意图；

(a)输入的三维模型，(b)为三维模型对应的离散点。

图3为邻域构建示意图。

图4为不同非线性基底函数的原函数、一阶导、二阶导函数曲线；

(a)原函数曲线、(b)一阶导函数曲线、(c)二阶导函数曲线。

图5为A₁、A₂、B₁、B₂四种非线性函数在相同外力拉伸下的仿真结果展示；

(a)为非线性函数A₁在相同外力拉伸下的仿真结果展示，

(b)为非线性函数A₂在相同外力拉伸下的仿真结果展示，

(c)为非线性函数B₁在相同外力拉伸下的仿真结果展示，

(d)为非线性函数B₂在相同外力拉伸下的仿真结果展示。

图6为各向同性与各向异性的弹簧床仿真结果比较；

(a)为各向同性的弹簧床形变，(b)各向异性版本。

图7为不同各向异性方向对挤压形变产生褶皱纹路的影响对比；

(a)高强度方向与水平方向夹角为0°，(a)高强度方向与水平方向夹角为45°，

(c)高强度方向与水平方向夹角为90°。

图8为armadillo模型下坠至地面过程，

(a)为线性材料发生了明显翻转，(b)为非线性材料有效防止了翻转。

图9为不同泊松比材质在同样拉伸条件下的形变结果；

(a)泊松比0.0，(b)泊松比0.3，(c)泊松比0.49。

图10为扭曲形变示意图。

图11为使用不同弯曲防翻转力系数对布料的褶皱行为影响；

(a)较小弯曲控制系数，(b)较大弯曲控制系数。

图12为兔子耳朵被固定，受重力拉伸的形变示意图；

(a)为FEM方法StVK材质，(b)为本发明方法和对应的各向异性非线性材质。

图13为一盆植物的仿真示意图；

(a)受重力下垂，(b)牵拉叶子，(c)牵拉茎。

具体实施方式

为使本发明的上述目的、特征和优点能够更加明显易懂，下面通过具体实现和附图对本发明做详细说明,但不构成对本发明的限制。

本发明方法的流程图如图1所示，整个流程分为两部分。第一部分是预处理部分，需要对输入的三维模型进行体积、质量的离散化，然后根据邻域半径值构建邻域关系，最后设定本构模型中自定义非线性基底函数、各向异性函数和杨氏模量、泊松比等常用材质参数确定材质本构模型。第二部分为仿真主循环部分，由计算内力，施加外力及摩擦力，施加弯曲防翻转力，数值积分求解下一时间步位置，碰撞检测与响应。

1、三维模型离散化。离散化方法有多种，可以采取体素化成六面体网格提取网、格顶点，也可以采用利用有向距离场提取三维模型内部网格顶点。由于采用的离散化方法是均匀离散，所以每个离散点的体积也是相同的，均设为离散分辨率d的3次方。质量的离散则是用体积乘以该点处的密度。如图2所示，左边是一个三维模型，右边是通过体素化方法提取的离散点数据。

2、邻域关系构建。通过设定邻域半径值δ(一般为d的1～3倍之间)，采用一般的空间hash方法，计算每个离散点的邻域点信息，并将其保存下来。由于本方法采用了GPU加速，所以邻域包含的粒子数一般会设置上限(30～50)。如果有粒子的邻域内存在超过上限的邻居粒子，则只取距离中心点最近的若干个点直至达到上限。如图3所示，离散点x_i,x_j拥有各自的邻居粒子，被包括在他们各自的球形区域内。在构建邻域的同时，一般会同时计算邻域内的权重值w_ij。w_ij可以多种设置模式，比如等。设置完之后还需要对权重进行归一化保证，∑_jw_ijV_j＝1，其中V_j为离散点j的体积。

3、设定本构模型中自定义非线性函数、各向异性函数、杨氏模量、泊松比。虽然g、h等非线性函数可以是满足前述条件的任意函数，但本方法仍然给出了一组非线性基底函数，用户只需要将g、h函数定义成这组基底函数的线性组合即可。

如图4所示从左到右分别是A1、A2、B1、B2四种基底函数的原函数、一阶导、二阶导图像。图5是其他所有参数相同，仅非线性函数采用不同基底在相同外力下产生不同的非线性形变现象。

各向异性函数G虽然可以是满足前述条件的任意函数，但本发明中采用了如下形式：

其中Σ为3x3的对角矩阵，对角元代表各自各向异性方向上的强度修正值。R是一个旋转矩阵，三个列向量分别代表3个互相垂直的各向异性方向，x_i代表未形变时的离散点x的位置，y_i代表离散点x的当前位置。如图6所示，左侧是各向同性的弹簧床形变，右侧是各向异性版本。又如图7所示，Σ不变，从左到右R矩阵旋转0°、45°、90°，三个被挤压的圆柱形薄壳会产生不同数量的皱纹而且会产生不同的纹路。

4、通过形变内能对粒子位置进行微分计算内力。本发明采用的实现方法如下

x_i代表未形变时的离散点x的位置，y_i代表离散点x的当前位置，其他符号含义均与之前表述相符。上述表示步骤3的内能表达式通过对离散点位置求导可以得到离散点的内力表达式。

5、施加外力与摩擦力。本发明中采用的摩擦力就是简单的和质量速度成正比的摩擦力可表示为f＝-αMv，其中α为阻尼系数。

6、施加弯曲防翻转力。利用邻域内粒子的初始位置信息和当前位置信息，可以估算出邻域的形变梯度 的估算方法如下公式描述：

在得到形变梯度估计量之后，再通过下计算弯曲翻转力

其中是一个翻转力系数，它可以是任意一个合理的正实数，但在本发明中使用键在未形变时的刚度值替代。

7、数值积分求解下一时间时刻位置。

时间积分的数值解法有很多可选方案，最简单的是欧拉显式积分。但是出于仿真效率和稳定性考虑，本方法采用隐式欧拉法并利用拟牛顿法进行求解加速。隐式方法求解下一个时间步形变体位置信息，等同于求解如下能量最小化问题

记

其中qⁿ⁺¹表示要求的下一个时间步的位置向量，表示初始估计位置，M是质量对角矩阵，Δt表示时间步长，W(θ_i,τ_i)表示步骤(3)中计算的微元形变内能密度，V_i表示微元体积。

优化此能量函数时，先计算其近似Hessian矩阵其中

也是采用了键的未形变时的刚度值，A_ij是一个向量其i分量为1，j分量为-1，其余分量为0，I₃为3x3的单位矩阵，为kronecker积。反复利用近似Hessian矩阵和梯度负方向计算前进步长，并加上线性搜索(linear search)操作保证收敛性，直至qⁿ⁺¹收敛到误差范围内。

8、碰撞检测与响应。碰撞检测与响应采用了基于位置的方法(position basedmethod)方法直接在每个时间步计算的结尾利用碰撞约束调整粒子位置和速度值。具体可以参考公知的文献如：[Müller,M.,Heidelberger,B.,Hennix,M.,&Ratcliff,J.(2007).Position based dynamics.Journal of Visual Communication&ImageRepresentation,18(2),109-118.]

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

Claims (10)

1.一种超弹性材料的快速仿真方法，其步骤为：

1)对待仿真的三维模型进行均匀离散，并计算每个离散点的状态信息；

2)设定离散点的邻域半径δ，构建每一离散点的邻域连接关系；

3)根据待仿真的超弹性材料设置非线性函数、各向异性函数以及相关材料参数，确定待仿真的三维模型的材质本构模型；

4)依据所述材质本构模型和离散点存储的当前状态计算形变内能和对应产生的内力并施加在离散点上；

5)对施加内力的离散点施加外力、摩擦力和弯曲防翻转力；然后求解下一个时刻离散点的位置和速度；

6)对步骤5)处理后的离散点进行碰撞检测与碰撞响应，调整离散点的位置和速度；

7)重复步骤4)～6)，直到仿真结束。

2.如权利要求1所述的方法，其特征在于，所述材质本构模型中每一个形变体离散点所产生的形变内能或内力与其周围设定范围内的其余离散点有关；所述形变体离散点的形变内能包括各向同性的静压内能和剩余的各向异性内能；所述静压内能只和邻域体积变化相关，各向异性内能只和邻域内各键长有关。

3.如权利要求1或2所述的方法，其特征在于，所述材质本构模型为

其中，W是离散点x处的形变内能密度，θ表示离散点x所处邻域体积变化，τ表示离散点x邻域内其它点与中心点x所连键的长度变化，a,b为常量系数，是以离散点x为中心点的邻域，w表示邻域内所有离散点的权值，G表示邻域内所有离散点的各向异性量，g、h为非线性函数，dV_ξ为离散点x的体积，τ<ξ>、w<ξ>、G<ξ>代表τ、w、G在离散点x处的场在离散点ξ处的取值。

4.如权利要求3所述的方法，其特征在于，非线性函数g、h分别为一组基底函数的线性组合；其中该组基底函数为：

5.如权利要求3所述的方法，其特征在于，G<ξ>是任意三元函数，输入是一个三元向量代表的单位向量，输出是一个正实数，且满足方向对偶性，即方向完全相反的单位向量产生同样的权重值。

6.如权利要求3所述的方法，其特征在于，w<ξ>是一个只和初始键长有关的权重函数，其可以是任意的一元函数，其满足值域为正实数并且是经过归一化处理。

7.如权利要求1所述的方法，其特征在于，所述步骤5)中，利用数值积分求解下一个时刻离散点的位置和速度。

8.如权利要求1所述的方法，其特征在于，计算施加的所述弯曲防翻转力的方法为：利用邻域内的点的初始坐标和形变后坐标估计一个仿射变换矩阵 代表该邻域的综合形变情况；的行列式如果是负数，即邻域发生了翻转，则将的前两行交换位置从而使其行列式重新变为正数；通过在初始键向量上得到估算的形变后的未弯曲键向量，利用形变后键向量减去估算的形变后的未弯曲键向量所得的值作为键的弯曲量，然后构建基于弯曲量的弯曲能量，从而计算出弯曲防翻转力的大小。

9.如权利要求1所述的方法，其特征在于，所述状态信息包括体积、质量；所述相关材料参数包括杨氏模量、泊松比。

10.如权利要求1所述的方法，其特征在于，所述形变内能包括静压力内能和各向异性内能。