二阶矢量有限元实现：二维电场在介质中的衰减

简述

二维一阶的矢量有限元在之前的文章中已有实现，本文在高斯积分求解的基础上，实现三角形网格的叠层二维二阶矢量有限元。并分析其计算精度与计算效率。

二维电磁在介质中的传播问题的边值问题这里不再详细介绍，具体可以参考：

1.二阶矢量基函数

对于矢量基函数而言，首先依旧是得到三角形网格的线性基底面积坐标：

由此，可以得到一阶矢量旋度基函数表达式为：

对于二阶矢量基函数，则是保留完全的一阶矢量基函数，加上二阶矢量基函数的旋度部分，进而得到：

对于第三项面矢量基函数而言，虽然节点组合可以得到三组面基函数，但是由于其中任意一项与其他两项是呈线性相关的，因此仅需要考虑任意两项。具体展开后，每个三角形由八个矢量基函数组成，如下：    

可以发现，前三项基函数与一阶的矢量基函数是一致的，这也是叠层矢量基函数的特征，高阶项保留低阶项，以便于使用混合阶矢量有限元。

2.系数矩阵的推导

同样，由于引入了高斯积分来求解具体的系数矩阵，因此不需要完全推导N*N的系数矩阵，只需要推导得到其表达式即可。

因此，根据上述二阶矢量基函数，需要进一步推导的内容是矢量基函数的旋度。第一部分一阶矢量基函数的旋度直接给出结果：

第二部分的矢量基函数本身由梯度组成，因此也很容易得出其旋度等于零。

第三部分面矢量基函数进行推导，难度也不大：

最终得到三部分的旋度表达式如下：    

接下来的求解过程就交给高斯积分完成，最终求解得到每个三角形的系数矩阵，具体高斯离散求解如下所示：

3.系数矩阵的组装与边界条件

对于二阶矢量有限元，组装过程中需要注意对应棱边编号、正负以及对应面的方向确定。这一点是非常重要的。基本原理就是保证同一条边在全局上是具有唯一的方向的，同一个面在全局方向也是唯一确定的。

同时也需要把这些边的全局方向映射到局部单元上，以确保根据母单元顺序可以唯一的组装好系数矩阵。对于母单元而言，三角形的排列顺序为：

对于第一类边界条件，其加载原理与叠层节点的有限元一致，在一阶部分使用强加第一类边界，对于二阶部分强加为零。同时矢量有限元需要将节点位置的场值映射到棱边切向分量，具体实现可参考：    

4.结果展示

测试模型为10m*20m的矩形,频率为1000Hz的电场垂直入射Y方向，介质体电导率为1欧姆米，模拟电场在介质中的衰减过程。

Eg1.网格横纵剖分5*10*2共计100个三角形网格

其中棱边165个，未知数共计165*2+100*2=530个。下面是具体的三角形网格剖分。

系数矩阵的分布如下，可以看出两个棱边部分与两个面部分在系数矩阵中体现的很分明，总体上该系数矩阵还是大型稀疏矩阵。    

求解该矩阵方程后，得到结果，将二阶结果、一阶结果的电场与理论解进行对比，如下：

对比发现，相比于一阶矢量有限元，二阶有限元计算精度提高了30多倍，当然二阶的计算时长是一阶的2.79倍。二阶精度提升在与理论解对比上也是非常的明显。

Eg2.剖分网格为20*20*2=800个网格，共计4080个未知数。

在模型中间放置一个低电导率，得到的电场衰减分布规律如下：    

可以明显看出，电场Ex的实部虚部在介质分界面有明显的突变情况。这也体现了矢量基函数的特点：仅保证切向连续，不保证法向连续。

总结

1.成功实现了二维二阶矢量有限元，相对于一阶而言计算精度有明显提升；    

2.二阶矢量有限元实现的关键点需要注意单元于自由度之间的映射关系以及确保每个自由度全局方向唯一性。

3.矢量基函数的构成方面，可以理解为一阶矢量基函数为二阶节点基函数的旋度部分；二阶矢量基函数为二阶节点基函数的旋度、梯度部分加上三阶节点基函数的旋度部分。