三维混合阶有限元实现的关键技术-霍姆霍兹方程
简述
众所周知，仿真的根本目的是高效准确的获得物理规律，从有限元算法的角度，高阶有限元可以大大提升计算精度，但是会导致计算效率降低;低阶虽然计算效率高，但是精度不够。而混合阶有限元则是结合二者优点。具体而言，在同一个模型中，某些重要的区域采用高阶基函数，而不重要的区域采取低阶基函数，如此即可避免计算资源的浪费，又可提高计算精度。
混合阶有限元的具体理论就是使用叠层基函数的高阶基函数包含低阶基函数的特性。本文将重点介绍实现过程中的一些关键技术点。
边值问题使用霍姆霍兹方程，具体有限元方程的推导这里不再展示，请参考文章：三维高阶叠层与插值有限元细节对比-霍姆霍兹方程，这里给出有限元方程：
1.基函数与单元系数矩阵
二阶叠层基函数的理论这里不再赘述，具体参考三维高阶叠层与插值有限元细节对比-霍姆霍兹方程与泊松方程三维二阶叠层有限元实现-四面体，文章中有详细推导，这里直接给出结果。
二阶叠层基函数表达式如下：
有限元方程中第一项积分结果的单元系数矩阵为：
第二项积分结果的单元系数矩阵为：
从单元系数矩阵也能看出，二阶的系数矩阵包含了一阶系数矩阵，这也是混合阶有限元能实现的根本原因。
2.关键技术
a.棱边编号重组
不同于纯粹的二阶网格可以直接使用软件生成棱边的编号有限元由于同时存在1阶，2阶，因此需要的棱边总数不再是所有四面体的棱边数，而是标记为二阶网格的棱边数，是小于等于四面体的总棱边数的，因此原本的二阶网格的棱边编号不再适用。
具体的处理方法：根据原有的棱边编号，映射新的棱边编号。首先挑选出二阶四面体；然后把二阶四面体原有的棱边编号全部集 合起来，得到的就是标记的四面体的原有棱边编号集 合；去掉其中重复编号，剩下的就是实际二阶网格的棱边数目，对这些编号重新编号，即得到新的棱边编号，这时候也得到旧编号与新编号的映射关系；最后把标记为二阶的四面体单元中棱边的编号更新，通过旧编号去查找新的编号，一一替换即可。
两个四面体接触位置为三角形，因此需要把三角形的阶数降低为一阶。通过对比一阶、二阶基函数在四面体中的位置不难发现，一阶落在四面体顶点上，二阶的高阶部分落在棱边，因此需要把棱边部分降阶，即置零，表示高阶的棱边处理为“不处理”的状态，即可实现分界面上阶数对齐。
3.结果展示
针对下述网格，将单元中心坐标Z小于5的四面体设置为二阶单元，其余网格设置为一阶单元。通过控制Z小于程度，即可演示二阶与一阶网格的在Z方向的比例。
a.当频率为1e4Hz,网格全为1阶的时候，数值解与理论阶对比：
b.当频率为1e4Hz,网格全为二阶的时候，数值解与理论阶对比：
c.当频率为1e4Hz,当网格Z<5为二阶，数值解与理论阶对比：
d.当频率为1e4Hz,网格Z<8为2阶的时候，数值解与理论阶对比：
编号 | 未知数 | 最大误差% |
全1阶 | 593 | 0.557 |
全2阶 | 3847 | 0.01136 |
Z < 5 | 2286 | 0.3282 |
Z < 8 | 3196 | 0.2708 |
为了说明混合阶更大的优点，考虑当频率升高到1e6Hz的时候，这时候电场在z=0的附近衰减的很快，梯度变化很大，导致这附近的误差更大，如下图：
e.当频率为1e6Hz,全部为一阶的时候，数值解与理论阶对比：
f.当频率为1e6Hz,全部为二阶的时候，数值解与理论阶对比：
g.当频率为1e6Hz,当网格Z<5为二阶，数值解与理论阶对比：
对比分析：1阶的计算结果误差已经非常的明显，达到了0.15%左右，对比而言，二阶的结果要好得多，最大误差在0.05%不到，但是使用的未知数达到3847。而使用混合阶后，对比完全2阶的误差结果，趋势基本上一致，最大误差也低于0.05%，但是未知数降低到了2286个。可见本例中，混合阶有限元在达到与全二阶几乎同等精度的情况下，未知数减少，计算效率得到了改善。
同时，从这个例子也可以得出一个混合阶使用规律，高阶使用梯度变化大的区域，而低阶使用在梯度变化小的区域。
总结
1.混合阶有限元实现关键技术点：对全部编号的重排；对混合阶数分解面的处理；
2.一个混合阶有限元阶数分布规律：仿真场梯度变化大的区域使用高阶，梯度变化小的区域使用低阶，分配得当，混合阶能达到与全高阶几乎一致的精度。
