适用于电磁场有限元计算的网格剖分算法

章春锋 汪 伟 吴天纬 安斯光

(中国计量大学机电工程学院 浙江 杭州 310016)

0 引 言

优异的适应场域边界几何形状以及媒介物理性质变异的能力，使有限元法成为各类电磁场、电磁波工程问题定量分析和优化设计的主导数值计算方法之一[1]。一般而言，有限元分析划分为建立模型 (前处理)、计算求解、结果处理和评定(后处理)。其中各阶段所用的时间占比分别为：40%～60%、5%～10%、30%～50%。因此，如何将实际工程中复杂场域离散为“有限单元”，这一剖分问题在数值分析过程是极其关键的一步，直接决定了最终结果的准确性、有效性[2]。

网格生成器诸如大型商业软件CAE/Maxwell，以映射法为理论基础[3]。其主要思想是，通过合适的映射函数将待剖分的几何域映射到参数空间形成规则的参数域；对规则参数域进行网格剖分；将参数域的网格反向映射回几何域，从而得到几何域的有限元网格。映射法是非全自动方法，必须通过人工交互方式，将剖分对象先剖分成具有简单拓扑关系的子域[4]。以映射法为核心算法的这些商业软件价格昂贵、操作复杂，且在生成大规模三角网格时速度较慢，往往是几乎不可访问的复杂代码，它们通常被用作“黑盒”，很难与其他代码或者软件进行交互，因此用户放弃了控制[5]。理解和使用网格生成算法的能力，好比掌握数据可视化、计算机图形中的几何建模的方法一样，是非常有价值的选择，值得深入研究[6-7]。

Persson-Strang算法是基于 MATLAB 的有限元网格自动生成算法[8]。与目前常用的Gmsh和Triangle开源网格剖分技术相比，Persson-Strang算法对一般的二维和三维的建模都能进行良好的处理，具有程序简短明晰、网格质量高、可移植性好的突出优点[9]。但是Persson-Strang算法并不完善，它可能不会终止并返回一个形状良好的网格，剖分速度相对较慢，且不易与实际应用相结合。目前，网格剖分存在剖分过程耗时长和网格质量不高的问题[10]。

本文主要从四个方面构造新的网格剖分算法：① 通过引入新的平滑函数减少迭代次数；② 利用坐标映射技术，避免由于微小形变而需重新剖分引起的计算资源浪费；③ 加入质量评估来作为终止条件，避免过度迭代；④ 对可能出现的重复节点进行删除处理。对于有限元应用解决了一些问题：确保边界节点位于边界上，提取有限元边界信息，将三角形节点逆时针编号处理。

1 网格剖分原理

1.1 Persson-Strang算法原理

Persson-Strang算法的简单性是由于使用有符号距离函数来指定和描述要网格化的几何区域。距离函数指定从空间中的任何点到域的边界的最短距离，其中函数的符号在区域外是正的，在内部是负的，在边界上是零。此定义用于判断节点是位于几何区域R2的内部还是位于几何区域的外部[8]。此外，距离函数d(x)的梯度指向边界的方向，允许外部的点有效地移回到域。几何域Ω可以定义为：

Ω={p=(x,y)∈R2|d(p)<0}

(1)

在该算法中，所需的网格大小由网格密度函数h(p)提供。网格密度函数h(p)是由用户指定的，其大致是区域的边界之间的距离。在几何域比较复杂的位置，可以使用自适应方法来实现局部特征变化[11]。对于均匀的网格，h(p)=1。

Persson-Strang算法采用基于力平衡的网格平滑函数来优化节点位置，该力平衡模型是三角形网格和二维桁架结构之间的类比。三角形的边对应于杆，顶点对应于桁架的节点。平衡问题对应一个非线性方程：

F(p)=0

(2)

找到一个定常解，满足式(2)。系统使用正向欧拉方法进行近似。在离散时间tk=kΔt时，近似解pk≈p(tk)更新为：

pk+1=pk+ΔtF(pk)

(3)

式中：Δt为方便离散所设置的时间分段值。

(4)

式中：eij=[pi,pj]为可变弹簧常数；f(pi,pj)取决于其当前长度‖eij‖和所需的长度rij。每一条边的力F取决于它当前的长度‖eij‖和理想长度rij之间的差值。对于均匀网格，rij是常数。但在许多情况下，在不同的区域有不同的尺寸是有利的[12]。轻微的非线性力函数可能会产生更好的网格，但分段线性力给出良好的结果，设置k=1，当‖eij‖=rij时，F=0，这是合理的。建议的边界处理方法，意味着没有任何点被迫停留在边界上，只是被阻止越过边界。

对于某些几何形状(如圆盘、椭球等)，创建带符号距离函数很容易。对于其他像凸多边形一样简单的几何域，距离函数很难得到表示。Persson-Strang算法通过引入固定点，可以使用简化的距离函数。除了比其他网格生成器更短、更简单之外，Persson-Strang算法还可以生成高质量的网格。

1.2 Persson-Strang可改进点

Persson-Strang算法是一种用于生成非结构化三角形网格的快速而灵活的工具，但是对于有限元应用而言存在一些主要缺点：

(1) 基于Laplacian函数的平滑算法仅描述了排斥力，并未考虑吸引力，造成迭代次数过多，剖分速度缓慢。

(2) 对于有限元方法(FEM)的数值模拟，生成合格的网格既复杂又耗时。在进行实际应用时，即使只有很小的变化，也需要为每组参数重复生成相应的有限元网格。重新网格化过程非常耗时。

(3) 算法执行缓慢和存在不终止的可能性。

(4) 有限元应用方面的不足。不能保证边界节点位于域的边界处，固定点和边界节点在网格剖分排序时可能出现重复。网格剖分核心函数Delaunay剖分生成的网格未逆时针排列编号。

2 改进的网格剖分算法

对于Persson-Strang网格剖分算法的不足之处，本文算法主要给出以下改进。

2.1 平滑过程

为确保有限元分析结果令用户满意，生成的有限元网格通常应具备以下条件：所有单元接近理想形状；变化梯度较大的地方网格密度应相应增大；粗细网格之间过渡均匀[10]。但通常情况下，在有限元网格自动生成器所生成的网格中总存在一些发生畸变的单元，剖分的几何域越复杂，畸形网格所占比例越大[13]。

本文研究通过平滑函数进行网格优化处理。平滑处理常用的Laplacian函数是应用最早同时也是最成熟的一种优化方法。其核心内容为：保持网格拓扑关系不变，将整个内部节点的位置摄动到由其相邻节点组成的多边形的质心处，使每个单元更接近于理想形状。将这个摄动过程遍历所有内部节点若干次，可较大幅度地提高网格质量[14]。经过Laplacian光滑处理的网格，不再具有Delaunay三角剖分的基本性质。

此外，化学中的Lennard-Jones函数考虑了吸引与排斥两种力两个方面的因素，给出分子间作用能表达式，被称为兰纳-琼斯势。兰纳-琼斯势函数是表示分子间相互作用势能的一种近似模型，如式(5)所示。

f(d)=d-13-d-14

(5)

式中：d是网格元素的边的长度。

Laplacian平滑函数只描述吸引行为，以帮助将节点分布到整个Ω域，但是由于只具有排斥力，不可避免出现迭代频繁；兰纳-琼斯势函数平滑函数效果受到数值不稳定的影响，有概率导致迭代不终止，陷入死循环[11]。综合两个函数的平滑特性，通过拟合方法，从数学上得到新的平滑函数NSF(New Smoothig Function)，NSF同时具有排斥力和吸引力，且在数值上具有稳定性：

f(d)=(1-d2)×ed3-d4

(6)

平滑函数如图1所示，其中：实点线为NSF平滑算法作用效果；虚线为Laplacian算法作用效果。

图1 平滑函数

(7)

2.2 自适应微小形变技术

提出一种稳定且简单的网格变形技术，用于快速准确地划分FEM网格。首先将边界网格定义为可伸缩变形的骨架，进而控制网格剖分生成的节点，不包括固定节点。更新所选边界网格的新坐标，就可以快速映射回细网格中节点的新位置。

在调用Delaunay三角形网格剖分后，进行网格优化之前，通过距离函数nfd(new function distance)定义新的边界网格框架，通过计算距离函数的数值梯度grad，将点投射到新的边界。

di(xi,yi)=nfd(xi,yi)

(8)

(9)

式中：deps为预设误差值，一般为10-3h0。

(10)

(11)

进行坐标映射之后，开始网格平滑优化过程，计算并存储每个节点的坐标,可通过重心公式判断是否需要细分，经验判断方式为：若重心在三角形内部，则不需要细分，若在外部，则需要细分。该算法既不需要网格重新剖分，也不需要添加新的未知参数。与现有方法相比，计算复杂度及其计算成本较低。该算法使用时，需关注网格质量q值，在形变较大时，仍须重新剖分以确保网格剖分的质量。

2.3 终止条件

该算法存在执行缓慢和不终止的可能性。实验表明通过附加的控制逻辑，可以使它更加健壮。终止标准基于当前迭代中节点的最大移动小于收敛精度dptol：

(12)

式中：Fr(pi)是每个节点pi的位移量；h0是初始网格尺寸大小。

由于一个高质量的网格有可能在达到该终止条件之前早已形成终止，标准应该包括一个质量评估，以避免迭代频繁。

在有限元应用中，误差上限取决于网格中的最小角度。要量化网格质量，通常使用的质量度量是最大内切圆半径rin与最小外接圆半径rout之间的关系，即：

(13)

式中：a、b和c是边长。根据经验法则，如果所有三角形的qmin>0.5，则网格是优的。在算法中加入该质量评估，以避免过度迭代。

2.4 问题节点处理

指定的固定点与应用概率拒绝方法后幸存下来的点可能重复，这将使得节点编号排序不正确，进而影响有限元求解。

对该问题进行如下处理：移除区域外的点，应用概率拒绝方法，再添加固定节点。在最后网格生成完成时，基于坐标判断是否重复，一般固定点位于边界处，在提取边界节点为有限元做准备后，在有限元解决计算过程中再次检查是否有重复节点。上述步骤在进行Delaunay三角剖分和平滑优化之前完成。

将三角形的任意两顶点交换，生成具有负面积的三角形排序，生成的三角形按标准逆时针形式排列；边界边在剖分域中只会出现一次，通过判断边的出现次数，实现边界节点筛选。

2.5 算法流程

本文算法流程如下：

步骤1初始化和参数设置：根据初始网格的大小h0，先把能涵盖欲划分区域的最大矩形划分为结构网格。根据距离函数定义，移除边界外的点。应用概率拒绝方法筛选点，当指定某些点要保留时，将保留的点加入，并删除重复的点。

步骤2Delaunay检索：执行Delaunay三角剖分并移除质心在外部的三角形。逆时针重新排列三角形顶点，形成不重复的边界边并提取。

步骤3微小形变处理：给出新的距离函数来定义新的边界框架，计算新距离函数的数值梯度，通过坐标映射技术将点投影到新的边界上。计算并存储坐标节点，防止出现节点重复。计算重心，判断是否需要细分。

步骤4网格平滑：计算和组装边缘力并移动节点。

步骤5边界节点：在边界上移动网格边界节点。网格的边界节点为边界边缘的端点。

步骤6终止标准：计算内部点的相对位移，并通过点来确定边长以计算网格质量q。如果在当前三角剖分中检测到大位移，则转到步骤2。如果在当前迭代中检测到大位移，则转到步骤3，如果q值均大于0.5，则终止。

步骤7有限元准备：清理并检查最终网格，输出包含边界节点序号、迭代次数等统计信息。

3 验证分析与应用

3.1 验证分析

实验条件：计算机配置为Intel(R)Core(TM)i5-7500 CPU @ 3.40 GHz。

取经典的二维复合几何网格生成实验提出的新的平滑函数对于网格剖分迭代次数以及网格生成质量的影响。采用新的平滑函数的网格剖分图如图2所示。

图2 NACA0012剖分

测试Laplacian函数和NSF平滑函数对于Persson-Strang算法的影响，在使用式(3)的前提下，取四组h0进行测试对比，h0为初始定义的网格尺寸大小，结果如表1所示。使用NSF平滑函数优化网格到达平稳状态所需要的迭代次数是使用Laplacian平滑函数进行网格优化所需要迭代次数的50%。

表1 平滑函数对于迭代次数以及时间影响

通过FEATool Multiphysics实现并集成了统一的MATLAB网格生成框架，基于该工具箱，可以直接比较Persson-Strang算法和本文算法网格生成时间和网格质量(h0=0.02)。

如表2所示，经过四次剖分测试(L1，L2，L3，L4)，新的网格算法相较于原先速度得到提升，其生成的网格质量q值经检测均在0.7以上。

表2 CPU时间(复合几何) s

通过二维示例验证所提出的微小形变处理方法。如图3所示。假定内部圆形在变形过程中被放大。关于所提出的网格更新方法，首先需要定义包括所有可移动节点的边界网格。重新定义距离函数，可以根据新的边界网格(如图4(a)中黑色虚线所标)和不变的区域坐标来映射精细网格的新位置，得到图4(b)。

图3 微小形变示意图

图4 确立新的距离函数之后的剖分图

得到CPU运行时间，以显示本文方法在计算效率上的优点。此二维精细网格中有399个节点和680个元素。再生网格的CPU时间为7.263 s，而针对相同问题的建议方法的计算时间为6.024 s,缩短了大约17%。两个网格中每个元素的平均最大与最小边缘比分别为1.21和1.33。

3.2 应用于电磁场领域

同轴线是由两根同轴的圆柱导体构成的导行系统，内外导体之间填充空气或高频介质的一种宽频带微波传输线。同轴线主要以TEM模的方式广泛应用于宽频带馈线和元器件的设计中。如图5所示，使用本文算法剖分一个外导体边长为2 cm，内导体边长为1 cm，内外导体电位差为1 V的方同轴线，进行电磁场有限元处理，计算方同轴线间的电位分布，处理步骤如下[15-16]：

图5 方同轴线图

(1) 网格剖分。

(2) 读入网格剖分准备的数据文件。

(3) 计算总体系数矩阵K。

(4) 处理强加边界条件。

(5) 求解有限元方程。

使用本文算法对该方同轴线进行剖分，在几何域比较复杂的地方，通过网格密度函数来进行有针对性的细化，网格剖分图如图6所示。与Persson-Strang算法的速度和网格质量进行对比，结果如表3所示。

图6 网格剖分示意图

表3 不同算法剖分得到的结果

可以观察到，引入NSF平滑函数后的本文算法，迭代次数是原先Laplacian平滑算法的60.9%左右，运行速度得到显著提升，得到的网格质量是平稳的，网格质量优，而且相较于Persson-Strang算法，不会出现网格质量显著下降的情况。

进行电磁场有限元处理后得到的等位线分布图和电位三维曲面图如图7和图8所示。坐标轴设置单位长度为1 cm，得到以下四点坐标：P1(-0.85,-0.85)，P2(-0.75,-0.75)，P3(-0.65,-0.65)，P4(-0.55,-0.55)。选取四点来验证电位，如表4所示。

图7 等位线分布图

图8 电位分布三维曲面图

表4 不同算法对电位求解的影响V

可以看出，本文算法在二维规则形状的剖分处理上，相较于Persson-Strang算法，速度得到显著提升，并且维持了优异的网格质量；与Maxwell专业有限元软件剖分结果相比，本文算法经过有限元后处理返回得到的电磁场参数与商业软件得到的结果精确度相似，这验证了新的网格剖分算法网格剖分速度以及质量方面改进的有效性。

4 结 语

本文提出一种新的网格剖分算法，基于Persson-Strang算法网格剖分原理，通过引入新的平滑函数来减少网格迭代优化过程中迭代次数；通过在Delaunay三角形剖分处理后新定义网格边界框，利用梯度坐标映射技术，避免由于微小形变而须重新进行网格剖分引起的计算机资源浪费；通过设置质量评估来避免网格过度迭代，以此提升网格剖分速度；对网格节点进行编号，并提取边界信息为有限元处理做准备；最后应用于电磁场领域，验证了该方法的可行性。