MATLAB上的有限元网格自动生成

刘 瑶，谭建国

(国防科技大学，湖南 长沙 410073)

0 引 言

一般地，可以将有限元分析[1]划分为：建立模型(前处理)、计算求解、结果处理和评定(后处理)。其中，各阶段所用的时间占比分别为：40%～60%、5%～10%、30%～50%。在数值分析过程中，网格划分[2]是极其关键的一步，可以说直接决定了最终结果的准确性、有效性。MATLAB[3]具有简洁的语言、丰富的基础库函数、强大的数据处理和可视化能力等突出特点，利用MATLAB进行有限元前处理是值得研究的问题。

DistMesh[4]是基于MATLAB的有限元网格自动生成程序，由加州大学伯克利分校的Per-Olof Persson和麻省理工学院数学系的Gilbert Strang开发，相较于MATLAB自带的PDE工具箱仅能求解比较简单的二维模型的缺陷，它可以生成三角形或四面体的网格，对一般的二维和三维的建模都能进行良好的处理，具有程序简短明晰、网格质量高、可移植性好等突出优点，在MATLAB有限元前处理方面应用广泛。

但是由于DistMesh用来筛选和优化节点的距离函数是解析表达式，对简单的模型来说可以得到良好的应用，用户却无法实现对复杂模型的构造，这限制了DistMesh的工程应用和推广。目前比较通用的三维建模方式[5]有网格建模、多边形建模、面片建模、非均匀有理B样条(NURBS)[6]建模等，而因为良好的数学特性和强大的建模功能，NURBS能够比传统的建模方式更好地控制物体表面的曲线度，从而创建出更逼真、生动的造型，成为CAD中几何定义的标准技术。通用的三维CAD软件也都采用NURBS建模方法。

文中基于MATLAB的NURBS工具箱[7]，对DistMesh进行改进，提出一种利用NUBRS曲线和曲面的方向来判断节点与曲线、曲面位置关系的方法，以完成对节点的筛选；再基于力的平衡原理对网格边的长度进行调整，以完成对节点位置的迭代优化。

1 网格生成器DistMesh

1.1 DistMesh设计原理

DistMesh基本原理是使用距离函数d(x,y)来定义几何边界，通过判断节点位置来筛选节点，再基于力的平衡来完成对网格的迭代优化。使用d(x,y)计算每个节点到边界的距离，规定在区域内部时的距离为负，外部为正，从而将在边界外的点去除。后续优化是基于三角形网格的边与桁架结构或者是线性弹簧之间的物理类比[8]。平面中的任何一组点都可以由Delaunay算法[9]进行三角化。在物理模型中，三角形网格的边(顶点之间的连接)对应于桁条，顶点对应于桁架的节点。

每一条边的力F取决于它的当前长度l及其理想长度l0的差值：

(1)

网格的疏密分布与几何形状的复杂程度有关[10]。在该算法中，所需的网格大小由网格密度函数h(x,y)提供。h(x,y)由用户指定，可以使用自适应方法[11]来实现局部特征变化，其大致是区域的边界之间的距离。自适应求解器可估计每个三角形中的误差，从而给出约束条件来精细化网格以获得更好的分布。

1.2 二维网格生成器算法描述

function[p,t]=distmesh2d(fd，fh，h0，bbox，pfix，varargin)

程序的输入参数如下：

(1)几何形状由距离函数fd给出。d=fd(p)返回值为每个节点位置P到其最近边界的带符号距离。

(2)所需的边缘长度h(x,y)由函数fh给出，为所有输入点返回h。

(3)h0是网格的初始大小。对于均匀网格h(x,y)=C。

(4)边界(最大外围矩形范围)。

(5)固定节点(必须出现在网格中的点)。

输出有：

(1)节点位置P：该N×2数组包含N个节点的x,y坐标。

(2)三角指标t：指定三角形的节点编号。

基本步骤为：

(1)在涵盖待划分几何形状的边界框中划分初始分布，初始分布可由用户任意规定。

(2)利用距离函数除去不在边界内的所有节点。

(3)进入主循环，节点位置不断迭代优化：

①去除质心在边界外的三角网格；

②计算理想的三角形边长，基于与实际边长的差值调整边长、移动节点；

③将移动到边界外的点拉回到边界上，避免节点流失。

(4)终止循环，终止条件由当前迭代中的节点最大移动距离决定。

1.3 存在的问题

距离函数是解析表达的，对于复杂图形很难构造，这大大限制了DistMesh的适用性。概括来说，算法中需要改进的地方有：步骤2，利用距离函数去除边界外的节点；进入循环后的第1步，用Delaunay方法划分为网格后，去除质心在边界外的三角网格；进入循环后的第3步，将移出的节点移回边界。

2 基于NURBS的网格生成方法

2.1 改进方法

(2)

图1 点的位置判断

同理，对于三维实体，如图2所示的球面，其法线方向是u方向和v方向的向量积。曲面法线朝外，从曲面外一点pout到其最近点pntm引向量，由于最近点pntm处法线和所引向量的夹角大于90°，二者数量积必为负；反之，对于曲面内的点数量积应为正。同理，当存在内表面时，若内表面法线方向朝内，则内表面内的点(即边界外的点)向最近点所引向量与最近点处法线的数量积为负，边界内的点为正。故可以定义外表面法线方向朝外，内表面法线方向朝内。

图2 通过点与曲面位置关系判断其内外

2.2 算法描述

二维平面指定调用函数：function[p,t]=nrbdistmesh2d(crvout,crvin,fh,h0,bbox,pfix)，函数的输入参数由fd变成了crvout和crvin；用户只需要使用NUBRS工具箱进行内外环曲线的建模。

三维实体指定调用函数：function[p,t]=nrbdistmesh3d(srfout,srfin,fh,h,box,fix),同理，用户只需要使用NUBRS工具箱进行内外曲面的建模。

基本步骤如下：

(1)先对涵盖用户待划分区域的外围矩形(长方体)创建初始分布：二维平面为三角形网格，三维实体为四面体网格。

(2)判断外环曲线(曲面)的方向：基于上述原理，先在边界外确定一点，通过该点在边界外这一既定结论反推曲线(曲面)的方向。

对于三维实体，从外围长方体表面上找一点pout，这一点必在边界外；找到外表面上最近点pntmout，求出pntmout所在处法线方向wout；再从pout向pntmout引向量，计算这个向量与法线方向的数量积zout；若zout小于0，说明外表面法线方向朝外，反之，改变外表面方向。

(3)判断用户是否输入内环(内表面)：如果存在，应当先定义其方向；否则直接执行第5步。

(4)定义内环曲线(内表面)的方向：内环曲线为顺时针方向，内表面法向向内。

对于平面，与外环原理相同，将外围矩形上的点替换为刚刚得到的外环上的A点，此点必在内环外；当这个点指向最近点的向量与最近点处导数的向量积z分量为正，说明内环是顺时针；否则，变换方向。对于三维实体，也从外表面上找点，符合要求的内表面，得到的数量积应该是正的，否则将内表面变换方向。

(5)除去边界外的点。与原方法计算距离的正负不同，这里仍然利用向量积(数量积)。

对于平面，因为满足条件的点都在曲线左侧或在曲线上，故而向量积的Z分量小于0的都是在边界外的，需要去除；对于三维实体，因为外表面法线朝外，内表面法线朝里，满足条件的点都在曲面法线方向的相反侧，则从节点向边界上最近点引向量，与最近点处的法线方向是相同的(小于90°)，故而数量积小于0的都是需要去除的。

(6)进入循环，进行节点的迭代优化。

①利用网格质心与内外边界的位置关系，去除质心在边界外的网格；

②计算理想的网格边长，基于与实际边长的差值调整边长、移动节点；

③将移动到边界外的点拉回到边界上，把边界外点相对应的边界上最近点赋给边界外点即可。

(7)终止迭代循环，与原方法相同。

3 验证与分析

3.1 网格质量评价指标

对于三角形网格，网格质量[12]常根据三角形最大内切圆与其最小外接圆的半径比率来判断：

(3)

其中，a,b,c是三角形的边长。对于正三角形来说，q=1；对于退化三角形，q=0。如果所有的三角形网格都有q>0.5，可认为是高质量的。

3.2 二维平面改进结果

(1)单位圆和带孔正五边形的优化过程如图3和图4所示，从左至右分别是迭代1次、10次、100次的效果。从图中可以看出，边界处的三角形网格变化尤为明显，边界也趋于清晰理想化，最终整体网格质量很高，q都在0.99以上。

图3 单位圆网格优化过程

图4 带孔正五边形网格优化过程

3.3 三维实体改进结果

图5为球体的四面体网格优化过程的剖面图，由形状千奇百怪最后趋于均匀的正四面体。

图5 球体四面体网格优化过程

3.4 结果分析

以上所有二维的例子q的平均值超过0.98，显著好于拉普拉斯平滑的Delaunay细化网格[14]。和原来的DistMesh进行对比，如表1所示，其网格质量仍然保持在高水平。

表1 改进算法与DistMesh网格质量对比

4 结束语

文中介绍了基于MATLAB的网格自动生成器DistMesh的主要原理，由于其核心的距离函数是解析表达的，对于复杂图形很难构造，故基于MATLAB中的NUBRS工具箱，将NURBS建模嵌入到算法中，解决了DistMesh不能实现复杂模型构造的缺陷。在二维平面的改进中，利用点到边界曲线上最近点处的向量，与最近点处切线方向向量的向量积，来定义曲线方向、判断点与曲线的位置关系，以筛选节点；在三维实体的改进中，利用点到边界曲面上最近点处的向量，与最近点处法线方向向量的数量积，来定义曲面的法线方向、判断点与曲线的位置关系，以筛选节点；基于力的平衡原理的迭代循环，通过对比理想的与实际的网格长度来不断调整网格边长，完成对网格的优化。验证结果表明，改进算法在保证网格生成质量的前提下，提高了DistMesh的适用性。