有限元网格划分

本文讨论了有限元网格的重要概念，包括单元的分类、有限元误差的分类与影响因素;并讨论分析结果的收敛性控制方法，并由实例说明了网格质量及收敛性对取得准确分析结果的重要性。同时讨论了一些重要网格控制的建议及其他网格设定的说明。

一、基本有限元网格概念

1.单元概述

几何体划分网格之前需要确定单元类型。单元类型的选择应该根据分析类型、形状特征、计算数据特点、精度要求和计算的硬件条件等因素综合考虑。为适应特殊的分析对象和边界条件，一些问题需要采用多种单元进行组合建模。

2.单元分类

选择单元首先需要明确单元的类型，在结构有限元分析中主要有以下一些单元类型：平面应力单元、平面应变单元、轴对称实体单元、空间实体单元、板单元、壳单元、轴对称壳单元、杆单元、梁单元、弹簧单元、间隙单元、质量单元、摩擦单元、刚体单元和约束单元等。根据不同的分类方法，上述单元可以分成以下不同的形式。

3.按照维度进行单元分类

根据单元的维数特征，单元可以分为一维单元、二维单元和三维单元。

一维单元的网格为一条直线或者曲线。直线表示由两个节点确定的线性单元。曲线代表由两个以上的节点确定的高次单元，或者由具有确定形状的线性单元。杆单元、梁单元和轴对称壳单元属于一维单元，如图1～图3所示。

二维单元的网格是一个平面或者曲面，它没有厚度方向的尺寸。这类单元包括平面单元、轴对称实体单元、板单元、壳单元和复合材料壳单元等，如图4所示。二维单元的形状通常具有三角形和四边形两种，在使用自动网格剖分时，这类单元要求的几何形状是表面模型或者实体模型的边界面。采用薄壳单元通常具有相当好的计算效率。

三维单元的网格具有空间三个方向的尺寸，其形状具有四面体、五面体和六面体，这类单元包括空间实体单元和厚壳单元，如图5所示。在自动网格划分时，它要求的是几何模型是实体模型(厚壳单元是曲面也可以)。

4.按照插值函数进行单元分类

根据单元插值函数多项式的最高阶数多少，单元可以分为线性单元、二次单元、三次单元和更高次的单元。

线性单元具有线性形式的插值函数，其网格通常只具有角节点而无边节点，网格边界为直线或者平面。这类单元的优点是节点数量少，在精度要求不高或

者结果数据梯度不太大的情况下，采用线性单元可以得到较小的模型规模。但是由于单元位移函数是线性的，单元内的位移呈线性变化，而应力是常数，因此会造成单元间的应力不连续，单元边界上存在着应力突变，如图6所示。

二次单元的插值函数是二次多项式，其网格不仅在每个顶点处有角节点，而且在棱边上还存在一个边节点，因此网格边界可以是二次曲线或者曲面。这类单元的优点是几何和物理离散精度较高，单元内的位移呈二次变化，应力呈线性变化，因此单元边界上的应力是连续的。但是在单元数量相同的条件下二次单元的节点数比线性单元的节点数多，模型的规模较大，如图7和图8所示。

三次单元的插值函数是三次多项式，其网格的每条边上存在两个节点，有些三次单元还具有内部节点。这类单元的离散精度更高，但是由于单元节点数较多，网格划分较为困难，模型规模很大，一般用于具有特殊精度要求的场合，如图9所示。

对于一阶和二阶单元，我们通常也称其为H单元。三阶及以上的单元，我们也称其为P单元，高阶次的P单元可以更好地拟合变形形状，特别对于曲率或者应力梯度变化较大的区域会较为真实的模拟，但会比H-单元有较多的运算量，如图10所示。

5.结构单元与非结构单元

根据单元能否离散成实际结构，可以将单元分为结构单元和非结构单元。能离散成实际结构的称为结构单元，如轴对称单元离散轴对称结构，杆、梁单元用于离散杆件结构，实体单元用于离散空间结构等，这些单元都属于结构单元。除此之外，还有一类单元并不用于实际结构的离散，而是在模型中模拟一些特殊的结构和边界条件，如质量单元用于实际的物体质量效应，弹簧和阻尼单元用于模拟结构的弹性支承和减振吸能部件，间隙和接触单元用于结构之间的相互接触作用，螺栓预紧力单元用于模拟螺栓的预紧力，刚体单元用于模拟节点之间的刚性连接等，这些单元称为非结构单元。

由于非结构单元非常抽象，使用起来有一定的难度，在设计仿真一体化分析里面通常会将其工程化，帮助使用者淡化其力学概念。

6.节点和单元的重要力学概念

针对前述的单元分类，此处要澄清关于节点和单元的一些非常有用的总结性概念。

◎有限元分析首先计算节点的位移量，接着再推算其对应单元的应变值，再计算积分点的应力。因此位移的准确性高于应变、应变高于应力;

◎当结构静力平衡时计算变形的单元是求得准确有限元分析结果的关键，因此线性计算中单元不可以变形过大，否则会造成求解失败;

◎网格质量概括来说，初始网格必须可呈现初始模型的几何形状，而且要足够“弹性”以符合静力平衡后的变形几何形状;

◎在预计会有应力梯度变化剧烈的位置上，为预测其准确变形情况，细小特征几何必须要更精确符合，以利于准确计算这些位置上的应力值;

◎在理想曲率边线与网格曲率边线之间的差距称之为离散误差。

二、有限元误差分析

1.有限元误差

有限元的误差主要来自两个方面，一是模型误差，一是计算误差。

模型误差是指将实际工程问题抽象为适合计算机求解的有限元模型时所产生的误差，即有限元模型和实际问题之间的差异。它包括有限元离散处理所固有的原理性误差，也可能包括几何模型处理、实际工况转化为模型边界条件时所带来的偶然性误差。

计算误差是指采用数值方法对有限元模型进行计算所产生的误差，误差的性质是舍入误差和截断误差。

模型误差包含离散误差、边界条件误差和单元形状误差，离散误差包含物理离散误差和几何离散误差。

2.离散误差

物理离散误差是插值函数和真实函数之间的差异，其大小与单元尺寸和插值多项式的阶次有关，单元尺寸减小也就是网格划分越密，插值函数的阶次增加，将使有限元的解收敛于精确解。

几何离散误差是指离散后的几何体与原有几何形状上的差异。对于由直线或者平面边界构成的规则结构，这类误差较小。对于具有复杂曲线或者曲面边界构成的结构，离散后会产生较大的形状误差。

本文下面通过SolidWorks Simulation来详细讨论物理离散误差与几何离散误差的具体操作细节。

三、收敛性及自动收敛方法

一般而言，网格拥有较多的单元，可得到较准确的结果。会有更多的节点可供计算，所以结果会较准确。较多的单元也就表示单元大小较小，所以物理离散误差可减小。实际分析上也有极限，在收敛性分析过程中网格尺寸一再缩减也不一定会对精确结果有帮助。

对一给定的几何而言，要达到收敛性的网格会与外部负载条件及边界约束条件有关(见图11)。在线性静态分析中，载荷大小不是收敛性的系数。

下面以SolidWorks Simulation的收敛性为例简单介绍收敛性的处理方法和技巧。

SolidWorks Simulation提供三种收敛性的技术。包括有手动控制收敛性和软件自动控制收敛性技术。其中自动控制收敛性的方法我们也称之为自适应方法，如图12所示，包括自动H自适应方法(H-adaptive)、自动P自适应方法(P-adaptive)。

1.自动H自适应方法(H-adaptive)

H方法的本质就是根据应力梯度的变化情况自动在应力梯度大的地方，根据预先规定的收敛准则，重新自动剖分网格，进行自动加密(见图13，原始网格与H自适应网格结果)。

SolidWorks Simulation的H方法具有以下特性。

1)适用于实体零件及装配体(仅支持实体单元)的静态分析研究;

2)在应变能误差较高的区域使用较小网格尺寸;

3)可以在应变能误差较低区域网格粗化(加大网格尺寸)，便于在后面的优化计算中降低计算规模，大大提高优化效率;

4)目标精度定义应变能量密度范数的精度等级默认值是98%，此处可以调整

能量密度范数的精度等级，一般情况下默认的精度可以达到分析的要求;

5)精度偏差设置(见图13);

◎精度偏差设置有局部(本地)和全局(整体);

◎滑动杆朝局部移动，指示程序以较少的单元取得精确的峰值应力结果;

◎滑动杆朝全局移动，指示程序取得整体零件刚度精确的结果，而不是应力结果。

6)若不确定，保持默认值即可;

7)网格粗糙化的目的是对应力梯度变化不大的区域，加大此处的单元尺寸，可以使用较小的网格得到较好的结果，同时也便于后继的优化求解。

图14展示了某机械零件，采用一阶单元不同单元大小，并采用自适应方法进行分析。方案1采用平均单元大小13.6mm 进行网格划分，然后采用H自适应网格划分，此时得到的最大应力点的应力是44MPa。方案2采用平均单元大小为3.4mm进行网格划分，然后采用H自适应网格划分，此时得到的最大应力点的应力是75.6MPa。两者之间有42%的差异。说明采用线性单元，使用H方法得不到准确的结果。

同样的模型，不做任何修改采用二阶单元进行网格划分，如图15所示，然后采用H方法，得到的应力误差小于3%，说明采用较大的全局二阶单元，然后采用H方法可以得到相当准确的应力结果。

2.自动P自适应方法(P-adaptive)

P方法的本质就是根据约束条件(如应变能)的变化情况自动在约束条件大的地方，根据预先规定的收敛准则，调整该处的单元形函数的阶次，在单元大小不变的情况下提高单元内部应力的精确性(见图16)。

SolidWorks Simulation的P方法具有以下特性。

1)适用于实体零件及装配体的静态分析研究，但装配体仅支持结合方式，不可以有其他接触存在;

2)收敛准则有总应变能、均方根合位移、均方根von Mises应力;

3)默认收敛准则是总应变能，均方根合位移及von Mises应力准则并不常用;

4)默认的设定通常就足够，由于系统通常会提前满足设定精度，因此最大p-order及最大循环数很少用到;

5)开始p-阶序起始于2，设为1会报错;

6)必须使用二阶单元为初始网格;

◎一般而言初始的网格尺寸影响很小(见图17);

◎Jacobian(雅可比检查)对误差有较大的影响，在局部位置无法反馈结果。建议将Jacobian检查设定在节点(见图18)。

四、手动收敛性检查

1.相对收敛性检查

在大多数复杂情况下很难通过自适应方法得到好的结果，必须通过相对收敛性检查得到收敛的结果(见图19)，其操作步骤如下。

1)执行多个分析研究，逐步调整加密网格，检查应力值的变化情形;

2)每次以2:1比例调整加细网格尺寸;

3)如果局部网格尺寸远小于整体网格尺寸，要留意扭曲失真的情况。

2.等值线质量检查(见图20)

1)应力等值线应该和连续几何体一样连续，使用不连续选项可以更清楚看到不连续的结果，如果几何体光滑连续而结果呈锯齿状，表明此处结果不好，需要加密网格或提高网格质量;

2)没有一个收敛的绝对测试标准，但是可以显示区域检查;

3)显示带有结果的单元结果可以显示他们如何影响结果。

3.误差估算方法1：能量范数值(见图21)

1)ERR能量范数误差绘图可以显示出相邻元素之应力值差异，理论上要愈小愈好;

2)并非绝对是收敛性测试方法，但会显示出要特别注意的高误差区域;

3)最佳的应用方式是在比较不同的设计过程。如果在某个区域的误差估算已经调整至最佳状况并维持一致，其分析结果的比较应该较具有实际意义。

4.误差估算方法2：节点和单元应力值比较(见图22)

1)节点解是临近单元的节点应力的平均值;

2)单元解是每个单元所有节点应力的平均值;

3)评定标准，理论上节点和单元应力值应该有较小的差异;

4)一般情况下，节点应力和单元应力的误差不允许超过5%。

五、手动收敛性控制和自动收敛性总结及建议

1)针对单一零件的分析

◎使用h-adaptive收敛方法，二阶质量单元及默认单元大小;

2)针对结合的装配体分析

◎使用p-adaptive收敛方法，二阶质量单元及默认单元大小;

◎如果有Jacobian误差发生，使用不同的初始网格尺寸执行可避免不正确的高应力值产生;

3)针对有连接接头或接触条件的装配体分析

◎使用传统的手动H-单元收敛方法;

◎使用二阶质量单元及默认单元大小;

◎使用初始网格控制以确保符合未变形的几何体; ◎使用局部网格控制在需要位置以达到收敛。(未完待续)