混合阶有限元的实现过程与详细分析:二维电磁场衰减
简述
关于二维问题的有限元数值模拟分析,之前分别介绍了一阶、二阶插值基函数、二阶叠层基函数的有限元实现过程,并测试验证了二阶有限元结果明显高于一阶,但是同时二阶求解的未知量个数也远大于一阶。在有些物理问题中,我们可能不要求研究区域中每个位置的精度都非常高,而只需要某些位置或者后处理结果精度高就可以了,因此整个区域使用高阶基函数可能会浪费计算和内存空间,因此混合阶有限元在优化计算效率方面就势在必行。例如本次测试的例子:二维电磁传播规律。在实际问题中,需要知道的是电磁波在传入研究区域后在入射表面位置的电磁场的变化情况(专业领域通过计算表面位置的视电阻率作为研究对象),下面是求解视电阻率的公式。
如图中,通过计算测试点(蓝色星号表示)的视电阻率的变化规律,来获得研究区域在测试面的视电阻率变化情况。
下面就详细介绍混合阶有限元的关键实现技术与结果分析对比。
1.叠层基函数与单元系数矩阵
红色框表示一阶所含有的内容,蓝色框表示二阶所含有的内容。可见高阶叠层基函数的特征一目了然:高阶系数项包含低阶系数项。
2.系数矩阵组装
混合阶有限元的两个关键点:重点1:对高阶部分的单元的棱边进行重新编号;重点2:对混合阶边界面的处理。
其他的单元系数矩阵组装和边界条件加载可以参考:二阶叠层基函数:二维电磁衰减数值模拟。这里用极简网格详细介绍两个重点内容的实现过程,网格如下图:
蓝色 区域为一阶、黄色 区域为二阶。给出了对应的节点编号和两个二阶三角形棱边编号。总共未知数为11个,四个单元。局部到全局编号为:
可以发现,由于存在两个三角形不是高阶,所以不能直接用所有的棱边编号,得重新对仅是高阶的三角形棱边进行编号。
第二个重点是在阶数分界面的位置需要特殊处理。
例如上述网格,阶数分界面含有两个节点编号3、4和一个棱边编号11,从一阶、二阶的单元系数矩阵不难看出,节点编号体现在一阶,而棱边编号表示的是二阶内容,阶数边界上要保持阶数一致,就得把高阶降至低阶统一。因为该棱边两侧最低阶为一阶,所以需要把二阶降至一阶,即把棱边表示的二阶内容去掉,与一阶统一。
具体在矩阵中实现就是对该棱边代表的行列去掉或者乘以大数置零处理。下面给出示意矩阵,可以看到系数矩阵的第11列的对角线的乘以大数以及右端项置零,就是对阶数边界的处理。
3.结果展示
a.以上述四个单元的网格为例,分别对一阶、二阶、混合阶(2小节的混合方式)计算结果,然后求解六个节点的误差结果,对比结果如下:
可见,1、2号节点由于是第一类强制边界,因此精度一致;3、4号节点在阶数分界面上,纯二阶精度最高,混合阶的精度整体高于一阶,而5、6号节点也是二阶精度最高、一阶精度最低。
b.设置5*6的网格单元的均匀介质,计算y=0位置处的视电阻率,测试混合阶有限元的优势。已知均匀电阻率为1的介质,在y=0处计算的理论视电阻率结果为1.其中混合阶y<10为2阶,其他为1阶,具体分布如下:
下面是66个网格节点的精度分布规律,可以看出,纯二阶精度明显要高于混合阶和一阶。而混合阶的精度整体精度略高于一阶,只有越接近y=0的时候,精度要更加高一些。
但是计算y=0位置的测点视电阻率的精度,发现混合阶几乎可以二阶媲美,远高于一阶计算精度。
可见,如果仅需要y=0处的视电阻率达到高精度,则不需要考虑每个节点的精度,这时候使用混合阶的优势就非常明显,就本例子而言,混合阶的未知数为151个,而二阶的未知数精度为231个,具体计算系数矩阵比例,二阶是混合阶多了2.3倍。c.设置5*10的网格,介质中存在不均匀电导率,分析混合阶不同数量下,精度与一阶、二阶的对比。如图所示,网格与测试的混合阶数量:
黑色框表示介质电阻率为100,蓝色、红色框等框表示不同混合阶,框内表示二阶部分。测试视电阻率如下:
一阶计算结果和二阶结果完全不一致,而混合阶在y<6区域使用二阶基函数就已经达到了纯二阶的计算精度。而此时混合阶y<6的未知数相比于二阶的231个未知数而言,只有119个。4.总结
1.介绍了二维三角形网格的混合阶有限元的实现,并详细介绍了其中关键步骤的实现过程;2.混合阶的整体精度是高于低阶,其精度提升的地方主要是远离一阶部分的二阶,因此仅需要对感兴趣的区域使用高阶即可。3.混合阶有限元在达到同等精度的时候,计算效率与内存资源利用往往远优于纯粹使用高阶基函数,这一点是混合阶有限元的优势所在。