CN113806951B - 一种基于半边数据结构的自然邻近点搜索的弹性仿真方法 - Google Patents

Abstract

本发明提供了一种基于半边数据结构的自然邻近点搜索的弹性仿真方法，具体为：建立计算对象模型，对模型进行节点离散，读入计算对象的材料数据和位移约束信息；根据离散节点划分Voronoi图与Delaunay三角形，其中涉及搜索每个积分点的自然邻近点；求解离散***的平衡方程，计算节点位移、应变、应力等物理量，具体为：通过高斯消元法求解平衡方程，得到节点位移，通过物理方程与几何方程求得计算模型的应力、应变；本发明提供的方法明显减少了三角面片外接圆检验计算量，减少了大量无效的外接圆检验计算，大大提高了自然邻近点的搜索效率，对于固定形状的弹性问题求解域有相当高的计算效率，适合弹性静力学求解问题。

Description

一种基于半边数据结构的自然邻近点搜索的弹性仿真方法

技术领域

本发明涉及弹性仿真领域，特别是指一种基于半边数据结构的自然邻近点搜索的弹性仿真方法。

背景技术

自然单元法有效地解决了有限元法的前处理网格划分复杂与网格畸变问题，与其他无网格方法无法直接施加本质边界条件的问题；但自然单元法的计算量远大于有限单元法，因而自然单元法求解问题的计算时间复杂度通常也大于有限单元法。造成自然单元法计算效率低的主要因素有：1)前处理离散节点的Voronoi图的形成；2)与其他无网格法一样，形函数构造均有一个搜索计算点影响域内节点的过程，自然单元法则是搜索计算点的自然邻近点以构造计算点的二次Voronoi图；3)形函数计算过程中需要的影响域内的节点数较多且计算过程较为复杂，但是在整个形函数求解过程中占用的时间比例较少，且离散节点的Voronoi图与自然单元法形函数的构造都是基于自然邻近点。

发明内容

本发明的主要目的在于克服现有技术中的上述缺陷，提出一种基于半边数据结构的自然邻近点搜索的弹性仿真方法，能够提高自然邻近点的搜索效率与形函数的计算效率，从而降低求解弹性问题的计算代价。

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

一种基于半边数据结构的自然邻近点搜索的弹性仿真方法，包括如下步骤：

步骤一：建立计算对象模型，读入计算对象的材料和位移约束信息；

步骤二：根据离散节点划分Voronoi图与Delaunay三角网格，具体步骤如下：

1)采用Delaunay三角剖分库对离散节点进行剖分；

2)将Delaunay三角网格的每一条边分解为两条反向的有向矢量边，构建Delaunay三角网格的半边数据结构；

步骤三：搜索每个积分点的自然邻近点，具体步骤如下：

3)在每个Delaunay三角面片中分布高斯积分点；

4)循环所有Delaunay三角面片；

3)搜索某个Delaunay三角面片内高斯积分点的自然邻近点；具体为：根据外接圆准则，该三角面片的三个顶点为该面片内的高斯点的自然邻近点，将其存入高斯点的自然邻近点数组中；根据半边数据结构的邻接关系，用外接圆准则检验相邻的第二个三角面片，如果第二个三角面片的外接圆包含积分点则将其三个顶点无重复地存入自然邻近点的数组中；将搜索移动到与第二个三角面片相邻的第三个三角面片，重复上述步骤，若第三个三角面片的外接圆不包含高斯点，则将当前搜索面片移回到上一个三角面片；当第二个三角面片的所有相邻三角面片都搜索完毕，将当前搜索移回到该三角面片并从与该三角面片另一条边相邻的三角面片出发搜索；重复上述步骤，直到与该面片的三条边均相邻的面片搜索完毕。

4)根据Sibson插值或Non-Sibson插值算法计算形函数；

5)组装总体刚度方程；

6)施加载荷与边界条件；

步骤四：求解离散***的平衡方程，计算节点位移、应变、应力等物理量，具体为：通过高斯消元法求解平衡方程，得到节点位移，通过物理方程与几何方程求得计算模型的应力、应变。

具体地，所述根据Sibson插值或Non-Sibson插值算法计算形函数，具体步骤为：

根据离散节点生成Voronoi图与二次Voronoi图；

采用Sibson插值构造形函数Φ_I(x)，即

进行形函数的构造，A_I(x)为计算点x的二次Voronoi图与其第I个自然邻近点的Voronoi图的相交面的面积，n为计算点x的自然邻近点数目；

采用Non-Sibson插值算法进行形函数的构造，即

式中s_J(x)为与节点J所对应的计算点x的二次Voronoi结构的边长，h_J(x)为计算点x与其自然邻近点J之间的距离。

具体地，所述组装单元刚度方程和总体刚度方程，具体为：

根据问题域Ω^h内的控制方程：

式中：σ为应力张量，b表示体力向量，Γ_u和Γ_t分别为问题域Ω^h内位移和力边界，u为位移向量，

为Γ_u上的位移向量，/>

为Γ_t上的面力向量，n为法向向量，/>

为微分算子；

令

和/>

分别为试探函数与检验函数的有限维子空间，则控制方程的Galerkin法变分弱形式为求/>

从而得：

式中，σ(u^h)和ε(v^h)分别为计算点x的应力与应变，计算点x处的位移试探函数u^h和检验函数v^h由相同的形函数线性插值而成，其插值公式为：

式中n为计算点x的自然邻近点个数。

得到矩阵形式的离散格式：

Ku＝f；

K＝∑K_e；

K_e＝∑K_IJ

K_IJ＝∫_ΩB_I ^TDB_JdΩ

式中，f为外力向量；K为总体刚度矩阵，K_e为自然相邻单元的单元刚度矩阵，I、J为自然相邻单元的节点；B_I为形函数偏导数矩阵，D为弹性矩阵；

其中，Φ_I,x为节点I处形函数对x的偏导数，Φ_I,y为节点I处形函数对y的偏导数，μ为泊松比。

具体地，所述施加载荷与边界条件，所述载荷的矩阵的离散格式为：

其中，f_I为载荷相对应的矩阵。

具体地，所述通过物理方程与几何方程求得应力应变，其具体过程为：

根据几何方程，求得节点应变：

ε＝Ba_e

式中，B为应变矩阵，B＝[B₀B₁…B_I…B_m]，

其中m＝0，1，2，…为积分点的自然邻近点数，u_m和v_m为自然邻近点m在x，y方向的速度；

根据物理方程，求得节点应力：

σ＝Sa_e

式中S＝[S₀S₁…S_I…S_m]，

/>

由上述对本发明的描述可知，与现有技术相比，本发明具有如下有益效果：

(1)本发明提供了一种基于半边数据结构的自然邻近点搜索的弹性仿真方法，对于自然邻近点的搜索，整个过程都是以高斯点(计算点)所在三角面片为中心的，搜索过程是一种局部搜索，相比全局三角面片循环搜索，本发明提出的方法明显减少了三角面片外接圆检验计算量，减少了大量无效的外接圆检验计算，大大提高了自然邻近点的搜索效率，对于固定形状的弹性问题求解域有相当高的计算效率，适合弹性静力学求解问题。

附图说明

图1为：Delaunay三角网格的半边数据结构示意图；

图2为：基于半边数据结构的自然邻近点搜索方法步骤示意图；其中图2(a)为搜索过程的示意图1，其中图2(b)为搜索过程的示意图2；

图3为：右端受集中力载荷的悬臂梁结构示意图；

图4为：悬臂梁问题域的离散节点及其节点的Delaunay三角化剖分；其中图4(a)悬臂梁问题域的节点离散；图4(b)悬臂梁的离散节点Delaunay三角化剖分；

图5为：采用Half-Edge NEM、三节点单元的FEM计算悬臂梁在y＝0处y方向的位移数值解与精确解对比图。

具体实施方式

以下通过具体实施方式对本发明作进一步的描述。

步骤一：建立计算对象模型，读入计算对象的材料、位移约束信息等。

步骤二：根据离散节点划分Voronoi图与Delaunay三角形，具体步骤如下：

1、使用Delaunay三角剖分库对离散节点进行剖分；

2、将Delaunay三角网格的每一条边分解为两条反向的有向矢量边，构建Delaunay三角网格的半边数据结构；

步骤三：搜索每个积分点的自然邻近点，具体步骤如下：

1)在每个Delaunay三角面片中分布高斯积分点；

2)循环所有Delaunay三角面片；

3)搜索某个Delaunay三角面片内高斯积分点的自然邻近点；具体为：根据外接圆准则，该三角面片的三个顶点为该面片内的高斯点的自然邻近点，将其存入高斯点的自然邻近点数组中；根据半边数据结构的邻接关系，用外接圆准则检验相邻的第二个三角面片，如果第二个三角面片的外接圆包含积分点则将其三个顶点无重复地存入自然邻近点的数组中；将搜索移动到与第二个三角面片相邻的第三个三角面片，重复上述步骤，若第三个三角面片的外接圆不包含高斯点，则将当前搜索面片移回到上一个三角面片；当第二个三角面片的所有相邻三角面片都搜索完毕，将当前搜索移回到该三角面片并从与该三角面片另一条边相邻的三角面片出发搜索；重复上述步骤，直到与该面片的三条边均相邻的面片搜索完毕；

具体结构如图1所示，半边数据结构的核心思想是一种基于边的分解假想，即将每一条边假想地分解为两条反向的有向矢量边，这两条半边组合在一起称为一条边，也就是一条边等于一对半边。每一条半边存储共这条边的其中一个三角面片、该半边的反向半边、前一条半边、后一条半边以及半边的出发节点；每一个顶点存储着以此顶点为出发点的半边以及点的坐标信息；每一个Delaunay三角面片存储一条半边。这样就形成了Delaunay三角网格的半边数据结构。半边数据结构(Half-Edge data structure)是一种以半边为中心的描述多边形几何体中点、边和面片之间关系的一种高效的数据结构。由于在这个结构中并不会存储网格的边的信息，存储的是半边，因此称为半边数据结构。

本发明提出的关于基于Delaunay三角化网格半边数据结构的自然邻近点局部快速搜索方法适用于以Delaunay三角形为背景网格的多点高斯背景网格积分。如图2所示，以11个节点情况为例，离散节点剖分后得到12个Delaunay三角面片；

面片1内采用三点高斯积分，高斯点位置如图2(a)所示。从面片1出发，显然面片1的外接圆包含高斯点，根据定义，该三角面片的三个顶点为高斯点(计算点)的自然邻近点，将其存入高斯点的自然邻近点数组中，并标记面片1为已检验(checked)。

(2)如图2(b)所示，由半边数据结构的邻接关系，用外接圆准则检验与面片1相邻的面片2，若面片2满足条件则将其三个顶点无重复的存入自然邻近点的数组中，并标记为已检验(checked)；

(3)将搜索三角面片移动到与面片2相邻的面片3，重复步骤(2)，若三角面片3的外接圆不包含高斯点(计算点)，那么将当前搜索面片移回到面片2。

(4)与三角面片2相邻的所有面片都已经搜索，此时将当前搜索面片移动回到面片1并从与面片1相邻的另一条三角面片出发搜索，重复步骤(1)-(3)直到面片1的三个相邻面片均搜索完毕，最后回到面片1，算法完毕。

(5)、根据Sibson插值或Non-Sibson插值算法计算每个积分点的形函数，其具体步骤为：

1)根据离散节点生成Voronoi图与二次Voronoi图；

2)采用Sibson插值，即

进行形函数的构造，A_I(x)为计算点x的二次Voronoi图与其第I个自然邻近点的Voronoi图的相交面的面积，n为计算点x的自然邻近点数目,或采用Non-Sibson插值算法进行形函数的构造，即

式中s_J(x)为与节点J所对应的计算点x的二次Voronoi结构的边长，h_J(x)为计算点x与其自然邻近点J之间的距离。

(6)、组装单元刚度方程和总体刚度方程，具体步骤如下：

1)根据问题域Ω^h内的控制方程：

式中：σ为应力张量，b表示体力向量，Γ_u和Γ_t分别为问题域Ω^h内位移和力边界，

为Γ_u上的位移向量，/>

为Γ_t上的面力向量，n为法向向量；

令

和/>

分别为试探函数与检验函数的有限维子空间，则式(1)的Galerkin法变分弱形式为求/>

从而得：

式中，σ(u^h)和ε(v^h)分别为计算点x的应力与应变。

2)计算点x处的位移试探函数u^h和检验函数v^h由相同的形函数线性插值而成，其插值公为：

式中n为计算点x的自然邻近点个数。

而后得到矩阵形式的离散格式：

Ku＝f

式中K为总体刚度矩阵，f为外力向量。

K＝∑K_e (4)

K_e为自然相邻单元的单元刚度矩阵

K_e＝∑K_IJ

I、J为自然相邻单元的节点

K_IJ＝∫_ΩB_I ^TDB_JdΩ (5)

式中，B_I为形函数偏导数矩阵，D为弹性矩阵。

其中，Φ_I,x为节点I处形函数对x的偏导数，Φ_I,y为节点I处形函数对y的偏导数，μ为泊松比。

6)、施加载荷与边界条件；

其载荷矩阵的离散格式为：

步骤三：求解离散***的平衡方程，得到节点位移u，通过物理方程与几何方程求得计算模型的应力应变，具体步骤如下：

通过高斯消元法求解平衡方程Ku＝f，得到节点位移。

根据几何方程，求得节点应变：

ε＝Ba_e

式中，B＝[B₀B₁…B_I…B_m]，

(m＝0，1..为积分点的自然邻近点数)

根据物理方程，求得节点应力：

应力：σ＝Sa_e

式中S＝[S₀S₁…S_I…S_m]，

u_m和v_m为节点m在x，y方向的速度

本发明提出的方法通过悬臂梁算例验证上述两种搜索算法在弹性问题中的正确性，并与传统自然邻近点搜索算法对比，证明基于半边数据结构的自然邻近点局部快速搜索算法在时间效率上的优越性，与有限元法对比证明自然单元法精度上的优势。

考虑如图3所示的悬臂梁问题，数值模拟中定义悬臂梁末端受竖直向下的集中力P＝1000N，长L＝480mm，宽D＝120mm，弹性模量E＝10⁷Pa，泊松比μ＝0.3，梁的厚度为1mm，不计自重，按平面应力分析。

悬臂梁在y方向的位移解析式为

本发明采用面向对象的程序语言C++编写了三节点单元的有限元方法程序(FEM)、r半径范围的局部搜索算法的自然单元法(MNEM)(蔡永昌,朱合华.基于局部搜索算法的自然邻接点方法[J].力学学报,2004,36(5):623-628.)、基于半边数据结构的自然邻近点局部快速搜索方法(Half-Edge NEM)和全局Delaunay三角形外接圆循环搜索的自然单元法(NEM)。采用175个节点对悬臂梁问题域进行离散，如图4(a)所示。对离散节点进行Delaunay三角化作为背景积分网格，同时该三角网格也是FEM的三角形单元，如图4(b)所示；

图5是采用本发明提出的Half-Edge NEM和三节点单元的FEM在y＝0处y方向位移分量数值解与精确解的对比分析图，从图中可以看出，Half-Edge NEM得到的位移数值解与精确解基本吻合，而三节点单元的有限元法得到位移数值解与精确解相差明显。

表1比较了Half-Edge NEM、MNEM和FEM在悬臂梁175个节点离散下各个求解步骤的计算时间，其中MNEM的r值按照文献蔡永昌等建议的节点排布的最大特征长度计算得到，r＝29。

表1悬臂梁数值模拟中基于各搜索方法的计算时间比较(单位：秒/s)

从表1可以看到，由于有限元法不涉及自然邻近点搜索的过程，且参与有限单元法形函数计算的仅为单元的顶点(该算例为三节点三角形单元)，形函数计算式为多项式计算，因此综合起来有限元法计算时间远低于自然单元法。自然邻近点搜索过程计算开销在总计算时间中占比最大，其次是形函数的计算，形函数形式为有理分式且参与计算的自然邻近点数目比有限元的单元节点多，其他计算步骤则相差不大。

在四种不同的自然邻近点搜索方法中，全局搜索的NEM的搜索时间达到了17.59s，效率最为低下；本发明提出的Half-Edge NEM的自然邻近点搜索时间仅为0.946s，总时间为1.644s，时间最短，效率最高，在一定程度上提高了自然邻近点的搜索效率，更加接近有限元的计算效率。从表中还可以看到MNEM的自然邻近点搜索时间为5.958s，总时间为6.653s；而NEM的全域自然邻近点搜索时间为16.850s，总时间为17.591s，该方法时间成本最大，效率最低。

上述仅为本发明的具体实施方式，但本发明的设计构思并不局限于此，凡利用此构思对本发明进行非实质性的改动，均应属于侵犯本发明保护范围的行为。

Claims (5)

1.一种基于半边数据结构的自然邻近点搜索的弹性仿真方法，其特征在于，包括如下步骤：

步骤一：建立计算对象模型，读入计算对象的材料和位移约束信息；

步骤二：根据离散节点划分Voronoi图与Delaunay三角网格，具体步骤如下：

1)采用Delaunay三角剖分库对离散节点进行剖分；

2)将Delaunay三角网格的每一条边分解为两条反向的有向矢量边，构建Delaunay三角网格的半边数据结构；

步骤三：搜索每个积分点的自然邻近点，具体步骤如下：

1)在每个Delaunay三角面片中分布高斯积分点；

2)循环所有Delaunay三角面片；

3)搜索某个Delaunay三角面片内高斯积分点的自然邻近点；具体为：根据外接圆准则，该三角面片的三个顶点为该面片内的高斯点的自然邻近点，将其存入高斯点的自然邻近点数组中；根据半边数据结构的邻接关系，用外接圆准则检验相邻的第二个三角面片，如果第二个三角面片的外接圆包含积分点则将其三个顶点无重复地存入自然邻近点的数组中；将搜索移动到与第二个三角面片相邻的第三个三角面片，重复上述步骤，若第三个三角面片的外接圆不包含高斯点，则将当前搜索面片移回到上一个三角面片；当第二个三角面片的所有相邻三角面片都搜索完毕，将当前搜索移回到该三角面片并从与该三角面片另一条边相邻的三角面片出发搜索；重复上述步骤，直到与该面片的三条边均相邻的面片搜索完毕；

4)根据Sibson插值或Non-Sibson插值算法计算形函数；

5)组装总体刚度方程；

6)施加载荷与边界条件；

步骤四：求解离散***的平衡方程，计算节点位移、应变、应力物理量，具体为：通过高斯消元法求解平衡方程，得到节点位移，通过物理方程与几何方程求得计算模型的应力、应变。

2.根据权利要求1所述的一种基于半边数据结构的自然邻近点搜索的弹性仿真方法，其特征在于，所述根据Sibson插值或Non-Sibson插值算法计算形函数，具体步骤为：

根据离散节点生成Voronoi图与二次Voronoi图；

采用Sibson插值构造形函数Φ_I(x)，即

进行形函数的构造，A_I(x)为计算点x的二次Voronoi图与其第I个自然邻近点的Voronoi图的相交面的面积，n为计算点x的自然邻近点数目；

采用Non-Sibson插值算法进行形函数的构造，即

式中s_J(x)为与节点J所对应的计算点x的二次Voronoi结构的边长，h_J(x)为计算点x与其自然邻近点J之间的距离。

3.根据权利要求2所述的一种基于半边数据结构的自然邻近点搜索的弹性仿真方法，其特征在于，组装单元刚度方程和总体刚度方程，具体为：

根据问题域Ω^h内的控制方程：

式中：σ为应力张量，b表示体力向量，Γ_u和Γ_t分别为问题域Ω^h内位移和力边界，u为位移向量，

为Γ_u上的位移向量，/>

为Γ_t上的面力向量，n为法向向量，/>

为微分算子；

令

和/>

分别为试探函数与检验函数的有限维子空间，则控制方程的Galerkin法变分弱形式为求/>

从而得：

式中，σ(u^h)和ε(v^h)分别为计算点x的应力与应变，计算点x处的位移试探函数u^h和检验函数v^h由相同的形函数线性插值而成，其插值公式为：

式中n为计算点x的自然邻近点个数；

得到矩阵形式的离散格式：

Ku＝f；

K＝∑K_e；

K_e＝∑K_IJ

K_IJ＝∫_ΩB_I ^TDB_JdΩ

式中，f为外力向量；K为总体刚度矩阵，K_e为自然相邻单元的单元刚度矩阵，I、J为自然相邻单元的节点；B_I为形函数偏导数矩阵，D为弹性矩阵；

其中，Φ_I,x为节点I处形函数对x的偏导数，Φ_I,y为节点I处形函数对y的偏导数，μ为泊松比。

4.根据权利要求3所述的一种基于半边数据结构的自然邻近点搜索的弹性仿真方法，其特征在于，所述施加载荷与边界条件，所述载荷的矩阵的离散格式为：

其中，f_I为载荷相对应的矩阵。

5.根据权利要求4所述的一种基于半边数据结构的自然邻近点搜索的弹性仿真方法，其特征在于，所述通过物理方程与几何方程求得应力应变，其具体过程为：

根据几何方程，求得节点应变：

ε＝Ba_e式中，B为应变矩阵，B＝[B₀B₁…B_I…B_m]，

其中m＝0，1，2，…为积分点的自然邻近点数，u_m和v_m为自然邻近点m在x，y方向的速度；

根据物理方程，求得节点应力：

σ＝Sa_e

式中S＝[S₀S₁…S_I…S_m]，

/>

Images