CN102306396B - 一种三维实体模型表面有限元网格自动生成方法 - Google Patents

Abstract

本发明涉及一种三维实体模型表面有限元网格自动生成方法，具有网格质量、密度和尺寸控制性能好，划分过程收敛性好，计算效率高、稳健可靠的特点；实现步骤如下：(1.1)建立三维几何模型，读入模型表面STL中性文件，重建三角形面片的拓扑关系；(1.2)根据三角形面片与相邻面片的关系，将表面三角形面片进行分组，每组三角形面片为一个特征面；(1.3)确定每个特征面的边界，确定每一个特征面边界上的特征边；(1.4)根据网格密度分布，在每条特征边上生成网格节点，将特征面的边界用边界上的网格节点环表示；(1.5)选择切割面和最佳剖分面，并在切割线和和剖分线上生成网格节点；(1.6)在每一个特征面上生成表面网格，进行光滑处理。

Description

一种三维实体模型表面有限元网格自动生成方法

技术领域

本发明属于工程数值分析领域，涉及一种工程问题的三维实体模型表面有限元网格自动生成方法，尤其适用于板料冲压成形、体积成形、碰撞、汽车结构等工程问题有限元数值分析的有限元网格自动生成。

背景技术

有限元等数值方法已广泛应用于各种科学与工程问题的分析，在应用有限元方法分析之前，需要对分析对象(即实体几何模型)进行网格划分，即将所分析区域划分成有限数量的单元，这些单元只通过节点相互连接。网格划分是一个费时且容易出错的过程，网格质量对后续的数值分析精度具有很大影响。随着所分析工程问题的复杂程度及其所分析区域几何形状复杂性的不断提高，对网格生成技术的要求也越来越高。在过去的三十年里，网格划分技术取得了很大进步，且出现了很多网格生成软件。到目前为止，任意平面区域的网格划分技术已较成熟，任意曲面上的三角形网格和任意三维实体模型的四面体网格生成也比较成熟，而任意曲面上的四边形网格生成和任意三维实体模型六面体网格的生成技术还远远没有达到成熟状态，需要解决的问题还非常之多。

任意曲面和任意三维实体模型三维表面网格的生成包括三角形或四边形网格的生成方法具有广泛的应用领域。对于三维实体模型，在采用某些方法在其上生成三维实体单元之前，通常需要在实体模型的表面上首先生成表面网格，然后以此表面网格为基础再生成整个三维实体模型的体网格。例如，在采用Delaunay或波前法生成三维实体模型四面体网格之前，需要在实体模型的表面上先生成三角形网格；采用铺层法或扫描法生成三维实体模型的六面体网格之前，需要在实体模型表面或部分表面上先生成四边形网格。在很多结构问题和金属成形问题的数值分析中，需要分析的对象为薄板或者厚度方向尺寸相对较小的结构，例如建筑结构、汽车结构、汽车碰撞、冲压成形等问题的注释分析，对于这类工程问题，采用有限元方法对其进行数值分析时，只需要在分析对象(实体模型)的表面上生成相应的网格即可，不需要生成三维的体网格。

表面网格生成方法可以分为直接法和间接法。直接法可以直接在实体模型表面上生成三角形或四边形网格，例如表面三角形网格可以采用波前法生成，表面四边形网格可以采用Paving方法生成。直接法生成的网格质量较好，但生成算法比较复杂，由于直接法生成表面网格的收敛性难以保证，当网格密度过渡较大或者表面几何形状过于复杂时，经常会出现网格生成过程失败的情况。映射法是间接法中生成表面网格较早的一种方法，该方法首先将要划分的表面区域采用一定的映射关系映射到平面上，然后对平面上的映射区域进行二维网格划分，最后再将生成的二维网格反向映射到原三维表面上。映射法生成网格的过程较为简单，但是三维实体模型的表面与映射的二维区域差别很大，即使生成的二维网格质量很好，但是反向映射后的表面网格可能出现畸变或退化的单元，网格质量和网格密度难以控制。中国专利200810011852.5“基于三维实体模型的四边形有限元网格生成方法”中采用栅格法，栅格法是间接法中生成表面网格的另一种方法，该方法首先将表面围成的空间划分成栅格单元，然后将表面的栅格单元投影到表面上生成表面四边形网格。栅格法生成网格效率较高，但是生成的表面网格与栅格的取向及其大小有关，网格中某些单元的质量可能很差。因此，三维实体模型表面四边形和三角形网格自动生成方法一直是工程数值分析领域普遍关心的问题，对于可靠的、能够直接自动生成三维表面网格的方法在板料冲压成形、体积成形、碰撞、汽车结构、建筑结构等工程问题数值分析中具有迫切需求。

发明内容

针对现有三维实体模型表面网格生成方法存在的问题，本发明的目的是提供一种可靠的直接自动生成三维实体模型表面有限元网格的方法，该方法不仅能够生成三角形单元，也可以生成全四边形单元，能够有效地控制网格的质量和网格密度，提高有限元分析的效率和精度，同时对任意三维实体模型六面体网格的自动生成也具有重要作用。

本发明是通过下面的技术方案来实现的：

一种三维实体模型表面三角形和四边形网格自动生成方法，包括以下步骤：

(1.1)首先采用现有CAD设计软件，对所分析工程问题进行三维实体造型或者三维表面造型，获得实体模型或其三维表面几何模型，输出其三维表面几何模型的STL中性文件，根据此中性文件，重建STL中性文件中的三角形面片的拓扑关系；

(1.2)根据三角形面片与相邻面片的关系，将表面的三角形面片进行分组，每一组三角形面片称为一个特征面；

(1.3)确定每一个特征面的边界，相邻特征面的公共边界称为特征边，确定每一个特征面的边界上的特征边；

(1.4)根据网格密度分布，在每条特征边上生成网格节点；

(1.5)在每一个特征面上生成表面网格。

进一步地，本发明步骤(1.1)中还包括指定网格密度分布，本发明步骤(1.1)中的STL文件类型可以是文本或二进制格式，步骤(1.1)还包括以下步骤：

(2.1)对三角形面片的顶点重新编号，相同位置的顶点采用同一编号，确定每一个三角形面片的三个顶点编号，并按逆时针顺序给出；

(2.2)确定每一个三角形面片周围相邻面片的编号；

(2.3)计算每一个三角形面片与周围相邻面片之间的夹角。

(2.4)将拓扑关系重建后的三角形面片作为背景网格存放网格密度信息，网格密度存放到三角形面片顶点处；

(2.5)三角形面片上任意位置的网格密度值由三个顶点的密度值通过线性插值确定。

进一步地，本发明步骤(1.2)还包括以下步骤：

(3.1)如果三角形面片与相邻面片之间的夹角在100°～220°之间，则三角形面片与相邻面片属于同一组，确定每一组面片中各面片的编号；

(3.2)对每一组中的面片作进一步的分组，细分后每一组中的面片法矢之间的夹角小于45°，并且每一面片的法矢与该组所有面片的平均法矢夹角小于60°，每当有新加入小组的面片时，均需要重新计算小组的平均法矢。

(3.3)将同一组的面片沿该组平均法矢的反向投影到一个平面上，检查同一组的面片是否有重叠，将重叠的面片归为一组，并从所属的原组中删除。至此，属于同一组的面片称为特征面。

进一步地，本发明步骤(3.2)和(3.3)中的平均法矢是采用各面片的面积作为加权因子对各面片的法矢进行加权平均得到的。小组面片的平均法矢按下面公式计算：

$\stackrel{\‾}{n}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{A}_i{n}_i/\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{A}_i---\left(1\right)$

其中，

为小组面片的平均法矢，N为小组三角形面片的数目，A_i和n_i分别是小组面片i的面积和法矢。

进一步地，本发明步骤(1.3)还包括以下步骤：

(4.1)对每一特征面中的所有面片，统计面片中没有被别的面片所共享的边，这些边构成了特征面的边界，将这些边首尾相连形成封闭的环。

(4.2)统计两相邻特征面的公共边界，这些公共边界构成特征边，确定每一条特征边上顶点的组成，确定每一个特征面由哪些特征边组成。

进一步地，本发明步骤(1.4)还包括以下步骤：

(5.1)根据网格密度分布计算在每一条特征边上生成的网格节点数目；

(5.2)确定每一条特征边上网格节点的位置。

所述步骤(5.1)中计算在每一条特征边(即空间曲线)上生成的网格节点数目又包括以下步骤：

在计算生成节点数目之前，需要先沿着空间曲线对网格密度积分，其积分公式为

${\∫}_L\mathrm{\ρ}\left(s\right)\mathrm{ds}=R\≅N-1---\left(2\right)$

其中，R为积分结果，N为包括两个端节点在内生成的节点数，L为空间曲线，ρ(s)为空间曲线上的网格密度，ds为沿曲线的微分。

对于特征边，特征边上顶点处的网格密度是已知的；对于切割线或剖分线，其形状是通过平面与三角形面片相交求得的，切割线或剖分线是由直线段组成的，故线段端点处的网格密度也可确定。在实际数值计算时，式(2)可以按下面的方法进行计算：用曲线上每一直线段两端点密度的平均值乘以该段的长度，然后对这些乘积求和，即可得到近似的积分结果R。

积分结果R为一实数，需要四舍五入到最接近的整数N-1。另外，如果要生成四边形单元，边界节点总数必须是偶数，所以，为满足这个要求，必要时还需要对N值进行偶数调整。

所述步骤(5.2)中确定每一条特征边上网格节点的位置，又包括以下步骤：

第1个和第N个节点的位置已经固定，但其余节点位置都是未知的。为了确定其余节点的位置，将节点的位置(相对起点的长度)作为一维变量，该变量满足下面的控制方程

ρ(s)ds＝Adξ＝C                                                    (3)

其中，ρ(s)为沿曲线的密度函数，ds为实际空间上的微分，dξ为计算空间上的微分(在计算空间上，节点是均匀分布的)，A为比例常数，C为常数。

将式(3)写成下面的两阶微分形式

$\frac{d}{\mathrm{d\ξ}}\left(\mathrm{\ρ}\left(s\right)\frac{\mathrm{ds}}{\mathrm{d\ξ}}\right)=0---\left(4\right)$

对于式(4)的微分方程，采用有限元迭代方法求得节点位置的数值解。将计算空间平均分成N-1个单元，共N个节点。经单元分析后，得到如下单元刚度方程：

$\mathrm{\ρ}\left(\frac{S_i+S_{i+1}}{2}\right)\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]\left(\begin{array}{c}{S}_i\\ S_{i+1}\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right)---\left(5\right)$

其中，S_i为第i个节点的位置，S_i+1为第i+1个节点的位置。

对所有单元的刚度方程进行组装，可得到整体刚度方程，将第一个节点位置S1＝0和第N个节点位置SN＝L(L为曲线的长度)作为边界条件，并假设其余节点的初始位置是均匀分布的，然后对整体刚度方程进行迭代求解，当求解过程收敛时，即可得到各节点的实际位置。

进一步地，对每一个特征面，本发明步骤(1.5)还包括以下步骤：

(6.1)特征面的边界用边界网格节点环来表示，外部边界的节点以逆时针顺序给出，内部边界的节点以顺时针顺序给出；

(6.2)选择合适的切割面，依次将内部边界与外部边界进行合并，计算各切割线的形状，根据网格密度分布在各切割线上生成网格节点，重新确定特征面边界上的网格节点环；

(6.3)选择最佳的剖分面，将特征面边界节点环围成的单连通区域剖分成两个区域，计算剖分线的形状，根据网格密度分布在剖分线上生成网格节点，对每一区域重新确定其边界节点环；

(6.4)对每一区域，重复步骤(6.3)，直到每一子区域不可再分解为止，即对于生成三角形单元，区域边界只有三个节点，对于生成四边形单元，区域边界只有四个节点或六个节点；

(6.5)对于生成四边形单元，将区域边界节点数目为6的区域按模板生成四边形单元；

(6.6)对生成的网格进行光滑处理。

所述步骤(6.3)中确定最佳的剖分面又包括以下步骤：

(7.1)剖分面与表面的交线称为剖分线；

(7.2)将区域的边界节点和区域内的三角形面片沿平均法矢反向投影到一个平面内，在该投影平面内，连接边界节点环内的两个节点，即可形成一条剖分线的投影线；

(7.3)确定最佳剖分线的投影线；

(7.4)确定最佳剖分面和剖分线。

所述步骤(7.3)中，确定最佳剖分线的投影线又包括以下步骤：

(8.1)对于四边形单元，如果节点i处的角度α_i小于60°，该节点就不能与其它的节点连接形成投影线；对于三角形单元，如果节点i处的角度α_i小于40°，该节点也不能与其它的节点连接形成投影线。

(8.2)如果连接节点i，j形成的四个投影角γ_ij1，γ_ij2，γ_ji1，γ_ji2中有一个小于30°；同理，对于三角形单元，如果连接节点i，j形成的三个投影角为20°，则节点i，j之间也不能形成投影线。

(8.3)如果连接节点i，j的线段与边界相交，则节点i，j之间就不能形成投影线。

(8.4)除了前面几种情况外，连接节点i，j就可以确定一条投影线。节点i，j配对后就不再与其它节点配对形成投影线。计算投影线ij的长度和投影线两侧子区域的面积比(≤1)。

(8.5)继续寻找其它可能的投影线，并计算各投影线的长度和两子区域的面积比。

(8.6)对所有可能的投影线按长度和子区域的面积比从小到大排序，选择子区域面积比大于0.3且长度最短的投影线作为最佳的投影线；确定完最佳剖分线的投影后，即可确定剖分线的空间形状和位置。

本发明的有益效果是：直接在三维模型的表面上生成网格三角形单元和全四边形单元，网格节点的位置是根据网格密度的分布确定的，因此可以有效地控制网格质量、网格密度和网格单元尺寸，可提高有限元分析结果的精度；将模型表面分解成若干个特征面，对每一个特征面采用递归分解的方式划分网格，从而保证了网格划分过程的收敛性和收敛速度；是一种计算效率高、稳健可靠的三维实体模型表面和三维曲面网格生成方法。

附图说明

图1为三维表面三角形和四边形网格自动生成流程图；

图2为两相邻面片之间的夹角示意图；

图3为实体的特征面和特征边；

图4(a)为圆柱面的特征面和特征边的识别；

图4(b)为球面的特征面和特征边的识别；

图5为切割线的确定示意图；

图6为剖分线的投影线的确定示意图；

图7(a)为六节点区域映射模式1；

图7(b)为六节点区域映射模式2；

图7(c)为六节点区域映射模式3；

图8(a)为一个连杆锻造过程数值模拟中的三维实体模型表面网格划分实施例；

图8(b)为一个十字形实体模型结构表面四边形网格划分实施例。

具体实施方式

下面结合附图与实施例，对本发明做进一步说明。

图1为三维实体模型表面三角形和四边形网格自动生成流程图。根据图1所示，在三维表面自动生成三角形和四边形网格的流程如下：

读入三维表面几何造型STL中性文件，在该文件中三维实体模型或曲面模型是用一系列表面三角形面片来表示，对三角形面片的顶点重新编号，建立表面三角形面片的拓扑关系，包括三角形面片的顶点编号以及三角形面片周围单元编号等，计算相邻面片之间的角度；根据面片之间的角度关系，将三角形面片进行分组，每一组构成一个特征面；确定每一个特征面的周围边界以及相邻特征面的公共边界，相邻特征面的公共边界构成一条特征边，确定每一条特征边由哪些顶点组成以及每一个特征面有哪些特征边构成；根据表面的网格密度分布在各条特征边上生成网格节点；对每一个特征面分别进行网格划分。在每一个特征面上生成网格的过程是：将特征面的边界用边界上的网格节点环来表示，特征面的外部边界上节点以逆时针顺序给出，内部边界上节点以顺时针顺序给出；选择合适的切割面依次将内部边界与外部边界合并，计算各切割线的形状，根据网格密度分布在各切割线上生成网格节点，重新确定特征面的外部边界节点环，最后特征面转换成只有一个外部边界的单连通区域；选择合适的剖分面将区域剖分成两个子区域，计算剖分线的形状，根据网格密度分布在剖分线上生成网格节点，确定每一子区域的边界节点环，对每一子区域以递归的方式剖分成两部分，直到所有的子区域不可再分为止；对于生成四边形单元来说，当子区域的边界节点数目为6时，采用模板映射的方法生成四边形单元；当在一个特征面上生成完网格后，需要对网格节点进行光滑处理，将不在表面上的网格节点调整到模型表面上。

下面详细阐述在三维实体模型表面上生成三角形和四边形网格的具体步骤：

1.三维表面几何形状的读入和三角形面片拓扑关系的建立

STL格式文件是大部分CAD系统能够输出的一种几何造型交换文件，它是对实体模型或曲面模型表面进行三角形离散化得到的，仅包含所有三角形面片的三个顶点坐标和面片的法矢数据。由于STL文件中的三角形面片是无序的，为了识别三维表面的特征信息，需要建立三角形面片的拓扑关系。首先对三角形面片的顶点重新编号，相同位置的顶点采用同一编号，这样这可以确定每一个三角形面片的三个顶点编号，编号是按逆时针顺序给出的；接着采用遍历的方式寻找每一个三角形面片周围的相邻面片，即寻找与之共边的三角形编号，对于内部的三角形面片有三个相邻面片，边界上三角形面片有一个或两个相邻面片；最后需要计算相邻两个面片之间的夹角，相邻两个面片之间的夹角定义为：将两个面片沿着它们公共边的方向投影到一个平面上，两面片投影线之间的夹角，如图2所示。

重建后的三角形面片还有一个重要的应用就是可以作为背景网格存放网格密度信息，根据表面曲率计算的网格密度或者用户指定的网格密度存放到三角形面片顶点处，三角形面片上任意位置的网格密度值可以由三个顶点的密度值通过线性插值的方式确定。

2.表面三角形面片的分组及特征面的确定

为了简化三维表面网格生成的难度，可以将三维几何表面分块，即将表面的三角形面片分组，每一组三角形面片构成一个特征面，然后在每一个特征面上生成网格单元。图3为内部有一方形通孔的立方体模型，该模型表面可以分成10个特征面，其中上下两个特征面内部各有一个方形孔洞。图3中的几何表面特征非常明显，特征边和特征面也容易识别。但对于图4(a)中的圆柱面和图4(b)中的圆球面，如果作为一个特征面，就不能直接生成网格，必须采用一定的规则，将这样的表面区域分解成若干个小的区域，这样的小区域仍称为特征面，小区域之间的边界仍称为特征边。

特征面的识别可以按照以下步骤进行：第一步，根据相邻面片之间的夹角对三维表面的三角形面片进行分组。如果相邻面片之间的夹角在100°～220°之间，这样的面片归为一组。第一组面片的搜寻可以从任一面片开始，如果它周围的面片与它的夹角在100°～220°之间，则将周围的面片加入该组，然后对新加入的面片，依次根据上面的方法决定其周围是否有新的面片加入该组。第一组面片搜寻完毕后，在剩余的面片中搜寻第二组面片，依次类推，直到所有的面片都分了组。对于图3所示的三维表面，总共有10组面片，这10组面片也就是最后的特征面。对于图4(a)，共有3组面片，而图4(b)，只有1组面片。

第二步，对第一步得到的各组面片作进一步的分组。对每一组面片，第1小组的搜寻从面积最大的面片T1开始，如果T1周围的面片与T1法矢之间的夹角小于45°，并且与第1小组所有面片的平均法矢小于60°，则将T1周围的面片加入第1小组。小组面片的平均法矢按下面的公式计算：

$\stackrel{\‾}{n}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{A}_i{n}_i/\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{A}_i---\left(1\right)$

其中，为小组的平均法矢，N为小组三角形面片的数目，A_i和n_i分别是小组面片i的面积和法矢。

每有新加入的面片，需要重新计算小组的平均法矢。按照同样的方法，依次对新加入的面片周围的面片进行处理，决定这些面片是否加入第1小组。第1小组搜寻完毕后，在剩余的面片中搜寻第2小组，依次类推。对于图4(a)中的圆柱面共细分为4个小组，图4(b)中的圆球面共细分为10个小组。

第三步，对每一小组的面片，沿着平均法矢的反向投影到一个平面上，如果小组中的面片的投影相互不重叠，则每一小组的面片作为一个特征面。如果有面片重叠，则将有重叠的面片单独归为一组作为一个特征面。

这样，经过上述三步就可以将任意三维表面上的面片进行分组，识别出各个特征面。确定每一个特征面后，搜寻每一个特征面的边界，包括外部边界和内部边界，内外边界都用边界顶点来表示，外边界顶点以逆时针顺序给出，内部边界顶点以顺时针顺序给出。然后搜寻相邻特征面的公共边界，公共边界构成特征边，确定每一条特征边有哪些顶点组成。这样，每一个特征面的边界都可以用特征边来表示。

3.空间曲线上的网格节点生成

在特征边上以及后面介绍的切割线和剖分线上生成网格节点的过程都是一样的，即在空间曲线上根据网格密度分布生成网格节点。在空间曲线上生成节点需要两个步骤。第一步需要计算生成节点的数目，第二步根据网格密度分布决定各个节点的位置。

要计算生成节点的数目，需要沿着空间曲线对网格密度积分，即

${\∫}_L\mathrm{\ρ}\left(s\right)\mathrm{ds}=R\≅N-1---\left(2\right)$

其中，R为积分结果，N为包括两个端节点在内生成的节点数，L为空间曲线，ρ(s)为空间曲线上的网格密度，ds为沿曲线的微分。

对于特征边，特征边上顶点处的网格密度是已知的；对于切割线或剖分线，其形状是通过平面与三角形面片相交求得的，因此，切割线或剖分线是由直线段组成的，线段端点处的网格密度也很容易求出。在实际数值计算时，式(2)可以按下面的方法进行计算：用曲线上每一直线段两端点密度的平均值乘以该段的长度，然后对这些乘积求和就可以得到近似的积分结果R。

积分结果R为一实数，需要四舍五入到最接近的整数N-1。另外，如果要生成四边形单元，边界节点总数必须是偶数个，所以为了满足这个要求，必要时还需要对N值进行调整。

确定了要生成节点的数目N后，下一步就要确定这些节点的位置。除了第1个和第N个节点的位置固定之外，其余的节点位置都是未知的。可以将节点的位置(相对起点的长度)作为一维变量，则该变量应满足下面的控制方程

ρ(s)ds＝Adξ＝C                                            (3)

其中，ρ(s)为沿曲线的密度函数，ds为实际空间上的微分，dξ为计算空间上的微分(在计算空间上，节点是均匀分布的)，A为比例常数，C为常数。

式(3)可写成下面的两阶微分形式

$\frac{d}{\mathrm{d\ξ}}\left(\mathrm{\ρ}\left(s\right)\frac{\mathrm{ds}}{\mathrm{d\ξ}}\right)=0---\left(4\right)$

对于上面的微分方程，可以采用有限元迭代技术求得节点位置的数值解。将计算空间平均分成N-1个单元，共N个节点。经过单元分析，可得到单元刚度方程如下

$\mathrm{\ρ}\left(\frac{S_i+S_{i+1}}{2}\right)\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]\left(\begin{array}{c}{S}_i\\ S_{i+1}\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right)---\left(5\right)$

其中，S_i为第i个节点的位置，S_i+1为第i+1个节点的位置。

对所有单元的刚度方程进行组装可得到整体刚度方程，将第一个节点位置S1＝0，第N个节点位置SN＝L(曲线的长度)作为边界条件，其余节点的初始位置假定是均匀分布的，对整体刚度方程进行迭代求解，当求解过程收敛时就可得到各节点的位置。

4.区域内外边界的合并

如果特征面有内部孔洞，需要将内部边界与外部边界进行合并，生成只有一个外边界的单连通区域。如图5所示的特征面就是一个包含内部孔洞的区域。可以选择一个合适的切割面，在区域上切割出一条厚度为零的缝隙，将内外边界连接起来。内外边界的合并可以按照下面的步骤进行：首先将区域的边界节点和三角形面片沿平均法矢的反向投影到一个平面上。在投影平面上，搜寻内外边界上距离最近的一对网格节点，连接这两个节点，连线即为切割线在投影面上的投影线，从投影面片中寻找与连线相交的面片，计算与面片边的交点位置，进而可计算在实际表面上交点坐标，这样切割线的形状和位置也就确定了；下一步就是根据网格密度在切割线上生成网格节点；最后将内部边界和切割线上的节点按顺序插入外边界节点环中，形成只有一个边界的节点环。如果特征面内有多个孔洞，可以采用上述方法依次合并内部孔洞。

5.最佳剖分面的确定

由于表面网格生成是采用区域递归分解方式进行的，必须选择合适的剖分面将区域一分为二。剖分面与表面的交线称为剖分线。剖分面的选择将对最后生成的网格质量有很大的影响。最佳剖分面可以按下面的方法进行确定：将区域的边界节点和区域内的三角形面片沿平均法矢反向投影到一个平面内，在投影平面内，通过确定最佳剖分线的投影线，进而确定最佳剖分面和剖分线。

在投影面内，边界节点的投影如图6所示，连接边界节点环内的两个节点就可以形成一条剖分线的投影线。最佳的投影线按下面的步骤进行确定：

(1)对于四边形单元，如果节点i处的角度α_i小于60°，该节点就不能与其它的节点连接形成投影线；对于三角形单元，如果节点i处的角度α_i小于40°，该节点也不能与其它的节点连接形成投影线。

(2)如果连接节点i，j形成的四个投影角γ_ij1，γ_ij2，γ_ji1，γ_ji2中有一个小于30°，节点i，j之间就不能形成投影线；同理，对于三角形单元，如果连接节点i，j形成的三个投影角为20°，节点i，j之间也不能形成投影线。

(3)如果连接节点i，j的线段与边界相交，则节点i，j之间就不能形成投影线。

(4)除了前面几种情况外，连接节点i，j就可以确定一条投影线。节点i，j配对后就不再与其它节点配对形成投影线。计算投影线ij的长度和投影线两侧子区域的面积比(≤1)。

(5)继续寻找其它可能的投影线，并计算各投影线的长度和两子区域的面积比。

(6)对所有可能的投影线按长度和子区域的面积比从小到大排序，选择子区域面积比大于0.3且长度最短的投影线作为最佳的投影线。

确定了最佳剖分线的投影后，就可以确定剖分线的空间形状和位置。然后根据网格密度分布在剖分线上生成网格节点，确定每一个子区域的边界节点环，以递归的方式对每一个子区域进行剖分，只到子区域不可再分为止。

如果生成四边形且子区域边界节点数目为6，可以采用模板映射的方式生成四边形单元。图7为六节点区域映射模板的三种基本模式，对于模式1连接相对边上的两个节点，将区域分解成两个四边形单元；对于模式2和3，分别在区域内部增加1个和2个新节点，将区域分解成三个和四个四边形单元。

6.网格光滑处理

在每一个特征面上生成网格之后，需要对内部网格节点进行光滑处理以提高网格的质量。每一个内部节点用其周围节点坐标的平均值来代替，此过程迭代2～5次就可以收敛。光滑之后的内部节点可能偏离原表面，最后需要将节点投影到表面来调整节点的位置。

7.三维实体模型表面网格自动生成实例

图8给出了两个表面网格自动生成的例子。其中图8(a)为一个连杆锻造过程数值模拟中的三维实体模型表面网格划分实施例，图中给出了模型表面三角形网格划分结果；图8(b)为一个十字形实体模型结构表面四边形网格划分实施例，图中给出了模型表面四边形网格划分结果。

Claims (7)

1.一种三维实体模型表面有限元网格自动生成方法，采用现有CAD设计软件，对所分析工程问题进行三维实体造型或者三维表面造型，获得三维实体模型或其三维表面几何模型，输出其三维表面几何模型的STL中性文件，其特征是，该方法还包括以下步骤： 

(1.1)读入STL中性文件，重建三角形面片的拓扑关系； 

(1.2)根据三角形面片与相邻面片的关系，将表面的三角形面片进行分组，每一组三角形面片构成一个特征面； 

(1.3)确定每一个特征面的边界及每一个特征面的边界上的特征边； 

(1.4)根据网格密度分布计算在每一条特征边上生成的网格节点数目，确定每一条特征边上网格节点的位置； 

(1.5)在每一个特征面上生成表面网格，并进行光滑处理，生成三维实体模型的网格； 

所述步骤(1.2)中将表面的三角形面片进行分组的步骤如下： 

(3.1)根据相邻面片之间的夹角对三维表面的三角形面片进行分组：如果三角形面片与相邻面片之间的夹角在100°～220°之间，则三角形面片与相邻面片属于同一组，确定每一组面片中各面片的编号； 

(3.2)对每一组中的面片作进一步的分组，细分后每一组中的面片法矢之间的夹角小于45°，并且每一面片的法矢与该组所有面片的平均法矢夹角小于60°，每当有新加入小组的面片时，均需要重新计算小组的平均法矢； 

(3.3)将同一组的面片沿该组平均法矢的反向投影到一个平面上，检查同一组的面片是否有重叠，将重叠的面片归为一组，并从所属的原组中删除，至此，属于同一组的面片称为特征面； 

所述步骤(3.2)和(3.3)中的平均法矢是采用各面片的面积作为加权因子对各面片的法矢进行加权平均得到的，小组面片的平均法矢按下面公式计算： 

其中，为小组面片的平均法矢，N为小组三角形面片的数目，A_i和n_i分别是小组面片i的面积和法矢，i=1,2,…,N。 

2.如权利要求1所述的一种三维实体模型表面有限元网格自动生成方法，其特征是，所述步骤(1.1)中，重建三角形面片的拓扑关系的步骤如下： 

(2.1)对三角形面片的顶点重新编号，相同位置的顶点采用同一编号，确定每一个三角形 面片的三个顶点编号，并按逆时针顺序给出； 

(2.2)确定每一个三角形面片周围相邻面片的编号； 

(2.3)计算每一个三角形面片与周围相邻面片之间的夹角； 

(2.4)得到三角形面片与相邻面片的关系。 

3.如权利要求1所述的一种三维实体模型表面有限元网格自动生成方法，其特征是，所述步骤(1.3)中，确定每一个特征面的边界及每一个特征面的边界上的特征边的步骤如下： 

(4.1)对每一特征面中的所有面片，统计面片中没有被别的面片所共享的边，这些边构成了特征面的边界，将这些边首尾相连形成封闭的环； 

(4.2)统计两相邻特征面的公共边界，这些公共边界构成特征边，确定每一条特征边上顶点的组成，确定每一个特征面由哪些特征边组成。 

4.如权利要求1所述的一种三维实体模型表面有限元网格自动生成方法，其特征是，所述步骤(1.4)中确定每一条特征边上网格节点的位置的步骤如下： 

第1个和第N个节点的位置已经固定，但其余节点位置都是未知的；为了确定其余节点的位置，将节点的位置即相对起点的长度作为一维变量，该变量满足下面的控制方程 

ρ(s)ds＝Adξ＝C            (1)

其中，ρ(s)为沿曲线的密度函数，ds为实际空间上的微分，dξ为计算空间上的微分，A为比例常数，C为常数； 

将式(1)写成下面的两阶微分形式 

对于式（2）的微分方程，采用有限元迭代方法求得节点位置的数值解；将计算空间平均分成N-1个单元，共N个节点；经单元分析后，得到如下单元刚度方程： 

其中，Si为第i个节点的位置，Si+1为第i+1个节点的位置； 

对所有单元的刚度方程进行组装，可得到整体刚度方程，将第一个节点位置S₁ =0和第N个节点位置S_N =L，L为曲线的长度作为边界条件，并假设其余节点的初始位置是均匀分布的，然后对整体刚度方程进行迭代求解，当求解过程收敛时，即可得到各节点的实际位置。 

5.如权利要求1所述的一种三维实体模型表面有限元网格自动生成方法，其特征是，所 述步骤(1.5)中，在每一个特征面上生成表面网格的实现步骤如下： 

(6.1)特征面的边界用边界网格节点环来表示，外部边界的节点以逆时针顺序给出，内部边界的节点以顺时针顺序给出； 

(6.2)选择合适的切割面，依次将内部边界与外部边界进行合并，计算各切割线的形状，根据网格密度分布在各切割线上生成网格节点，重新确定特征面边界上的网格节点环； 

(6.3)选择最佳的剖分面，将特征面边界节点环围成的单连通区域剖分成两个区域，计算剖分线的形状，根据网格密度分布在剖分线上生成网格节点，对每一区域重新确定其边界节点环； 

(6.4)对每一区域，重复步骤(6.3)，直到每一子区域不可再分解为止，即对于生成三角形单元，区域边界只有三个节点，对于生成四边形单元，区域边界只有四个节点或六个节点； 

(6.5)对于生成四边形单元，将区域边界节点数目为6的区域按模板生成四边形单元； 

(6.6)对生成的网格进行光滑处理。 

6.如权利要求5所述的一种三维实体模型表面有限元网格自动生成方法，其特征是，所述步骤（6.3）中确定最佳的剖分面又包括以下步骤： 

(7.1)剖分面与表面的交线称为剖分线； 

(7.2)将区域的边界节点和区域内的三角形面片沿平均法矢反向投影到一个平面内，在该投影平面内，连接边界节点环内的两个节点，即可形成一条剖分线的投影线； 

(7.3)确定最佳剖分线的投影线； 

(7.4)确定最佳剖分面和剖分线。 

7.如权利要求6所述的一种三维实体模型表面有限元网格自动生成方法，其特征是，所述步骤(7.3)中，确定最佳剖分线的投影线又包括以下步骤： 

(8.1)对于四边形单元，如果节点i处的角度α_i小于60°，该节点就不能与其它的节点连接形成投影线；对于三角形单元，如果节点i处的角度α_i小于40°，该节点也不能与其它的节点连接形成投影线； 

(8.2)对于四边形单元，如果连接节点i，j形成的四个投影角γ_ij1、γ_ij2、γ_ji1、γ_ji2中有一个小于30°，节点i，j之间就不能形成投影线；同理，对于三角形单元，如果连接节点i，j形成的三个投影角为20°，节点i，j之间也不能形成投影线； 

(8.3)如果连接节点i，j的线段与边界相交，则节点i，j之间就不能形成投影线； 

(8.4)除了步骤（8.1-8.3）的情况，连接节点i，j就可确定一条投影线；节点i，j配对后 就不再与其它节点配对形成投影线，计算投影线ij的长度和投影线两侧子区域的面积比； 

(8.5)继续寻找其它可能的投影线，并计算各投影线的长度和两子区域的面积比； 

(8.6)对所有可能的投影线按长度和子区域的面积比从小到大排序，选择子区域面积比大于0.3且长度最短的投影线作为最佳的投影线；确定完最佳剖分线的投影后，即可确定剖分线的空间形状和位置。 

Images