二维矢量有限元（结构化网格）实现流程

2月前浏览961

简述

在之前的文章，详细介绍了二维矢量有限元的实现过程，使用的是非结构化的三角形网格，在前一篇文章展示了一个几乎是结构化网格的实例后，发现结构化网格的矢量有限元也非常重要，因此本文采用相同的物理问题，具体详细介绍结构化（四边形）矢量有限元的实现过程。

物理规律的边值问题与对应的有限元问题不再介绍，可参考：二维矢量有限元的详细实现过程，这里给出最终需要离散的有限元问题：

1.四面体矢量基函数

不同于三角形矢量基函数的构造方法，结构化网格的四面体四条棱边本身就能完美体现电场的X和Y方向，因此很容易通过拉格朗日插值公式构造矢量基函数。

假设电场X、Y方向均以坐标轴的方向为正方向，如图所示，对于棱边上的Ex，Ey，通过拉格朗日插值公式，不难得出四面体网格内任意一点的Ex,Ey的插值公式：

对上述公式进行推导简化：

由此，得到四面体单元的插值基函数与插值公式。可以发现，基函数N的变化量与其方向是相互垂直的关系，是由简单的线性插值公式和与其垂直的方向组成。因此，根据梯度原理和旋度原理，很容易理解这四个基函数的梯度肯定为零，而旋度不为零。

上述的基函数特征正是矢量有限元的特征：弱解仅要求棱边的切向方向连续，而对法向方向不进行约束。这恰好满足电场的连续性规律：介质层分界面梯度方向不连续而切向方向连续。

2.单元系数矩阵推导

确定棱边基函数后，将该基函数带入有限元方程，可以得到旋度项：

四面体的面积等于长乘以高ab，写成方阵得到：

对于有限元方程的第二项积分，需要注意，矢量基函数具有方向性，有：

N1,N2表示X方向，N3,N4表示Y方向，因此这二者之间的乘以是等于零的，由此可以得到积分矩阵：

对上述矩阵的积分求解，与一维有限元基函数的推导方式一致的，具体可以参考最简单的一维有限元问题：求解cos函数分布，这里直接给出推导结果：

3.系数矩阵组装

为了详细的说明整个组装过程，这里以最简单的两个四边形网格为例，展示其中的网格单元的各个关系。

上述图中，黑色编号表示节点编号，红色表示棱边编号，圈红表示单元，因此可以得到每个单元与棱边与节点的关系。

其中单元与节点编号顺序是逆时针方向；单元与棱边编号是先长后宽，并以XY轴正方向作为棱边正方向，得到棱边编号：

在明确了单元-棱边-节点的关系后，就开始对每个单元系数矩阵进行组装，这里的物理模型参数：研究区域X=10米;Y=20米；频率10000Hz；电阻率为1欧姆米。

得到的单元1、2的系数矩阵是相同的，如下：

组装以后：

添加第一类边界条件后，最终有：

可以发现，系数矩阵matrix的对角线只有第二行没乘以大数处理过，其他均处理过，观察模型的网格剖分，的确仅第二行是内边界，其他七个边全部是外边界，而这些边界中，只有棱边1、5是沿着X方向的，所以体现在rhs中，只有1、5排有数值，其他边界为零。体现了在边界上电场仅沿着X方向传播。

4 求解结果

a.仅有两个四边形网格的棱边未知数结果：

不难发现，误差仅存在于2号棱边的解，这是由于其他棱边均是棱边，强加了第一类边界条件的原因。就数值上2号棱边的解基本上大体是一致的，误差还是非常大。

b.细化网格20*40，均匀介质，对比理论解

该结果的精度要明显高于三角形的矢量有限元，对比参考：二维矢量有限元的详细实现过程

c.研究区域内存在不同介质材料，中间区域的电导率为10欧姆米

对比三角形的结果，结果基本上可以吻合一致。对比参考：二维矢量有限元的详细实现过程

D 实际模型，具体参数结果参考：使用二维矢量有限元展示一个实际例子

中间区域为低阻10欧姆米，视电阻率结果：

中间区域为高阻1000欧姆米，视电阻率结果：

结果与三角形基函数矢量有限元的结果一致。但是计算时间要比三角形矢量基函数慢很多。

总结

1.从整个流程来看，四边形棱边基函数相对于三角形而言，更加容易理解与实现；

2.使用相同的网格，结构化矢量有限元计算精度要高于三角形棱边基函数结果，但是三角形基函数的耗时要远远低于四边形棱边基函数有限元。

来源：实践有限元

介绍几款二维网格剖分工具：三角形网格

简述本文列出作者用过的二维三角形网格剖分工具，以供参考。comsol（商业软件）Comsol作为多物理仿真的通用软件，二维网格剖分的功能也是多姿多样，提供给用户各种可选择的接口实现二维网格划分。网格结果可以输出各种格式，如果要自行解析文件还是比较麻烦，一般可以通过其他软件辅助，二次处理后生成容易解析的网格格式。Ansys-q2d（商业软件）q2d是集成于Ansys中的一款专业二维仿真软件，主要针对电磁相关的二维界面参数求解。专业性很强，网格划分质量不需要用户关注，自行使用自适应算法获得合乎仿真模型精度的网格。其生成的网格在业界是几乎是公认的最优网格。网格结果同样可以输出，解析格式不是太困难，但是其中仅包含网格，不含有边界、材料属性等信息。matlab-pdetool（商业软件）matlab是一款强大的数学商业软件，几乎搞科研算法的都有使用到，即使近些年被M国禁用。它里面也具有二维网格剖分功能，在命令行敲pdetool则可以弹出简单的二维网格剖分界面，可以进行简单的bool运算等建模功能。生成的网格可以保存为matlab脚本格式，很好读取其中的网格数据；Gmesh（开源代码）Gmesh的网格剖分可以说非常的完善，自成体系，拥有自己的界面，也提供了非常友好的接口，可以实现各种网格剖分、局部加密等功能。网格输出文件为msh格式，包含了建模的所有几何信息、材料信息、边界信息等。因此解析该格式相对不太容易。Triangle（开源代码）很完善的一套开源二维三角形网格剖分工具，用户可以通过命令的形式实现各种网格剖分，例如下图中区域划分，分区域加密等等。唯一的缺点是代码只支持X86,不支持X64,这一点需要特别注意。网格输出可供用户选择，支持单元、节点、棱边、邻单元等多种信息输出，也支持一些常用的格式输出。CDT（开源代码）相对于较为完善的tiangle,CDT则显得不是那么友善，似乎是专门给程序员使用，没有外部直接调用命令的方式，只能自己写接口实现输入参数和输出网格。不过适用于对网格的二次开发使用，例如其中不具备对单元属性的标记，不具备局部细化功能，似乎全局细化功能也不具备。网格信息的输出只用用户自己输出，并且只含有单元、节点信息。来源：实践有限元