二维有限元混合网格的详细实现过程

实践有限元

22天前浏览117

简述

在二维有限元问题上，该公众号已经介绍了大量案例，从网格的角度上介绍了四边形网格的结构化有限元和三角形网格的非结构化有限元。

非结构化有限元的优势在于三角形网格形函数简单，可以很好的模拟任意形状物理模型的边界；缺点是低阶精度不足，必须耗费大量网格或者高阶基函数来提高精度。

结构化有限元的优势在相同阶数上精度要优于三角形网格，同时所需要的网格也远少于三角形；缺点则是无法精确模拟复杂物理模型，对于复杂边界只能通过不断细化网格来契合模型。

通过简单分析，二者各有优缺点，因此如果能将两种网格混合使用，或许就能够实现各自优点，避免各自缺点。本次文章介绍二维四边形网格和三角形网格的混合有限元实现过程，详细分析实现过程与精度对比分析。

1.边值问题

依旧使用二维电场在介质中衰减的物理模拟，具体边值问题参考：二阶叠层基函数：二维电磁衰减数值模拟，这里给出最终的有限元推导公式：

具体研究区域与物理参数如下：

2.基函数与单元系数矩阵

既然想要实现混合网格有限元，则必须给出三角形单元和四面体单元的基函数的单元刚度矩阵。只考虑线性插值的情况下有：

三角形基函数，三角形基函数使用面积坐标推导得出，具体推导参考：二阶叠层基函数：二维电磁衰减数值模拟

单元系数矩阵:

四边形基函数由一维线单元的基函数组合而成，具体推导参考：介绍一二维结构化有限元的常见边值问题

单元系数矩阵：

3.系数矩阵组装

以最简单的网格剖分为示例，分析求解过程中每个步骤的实现过程。网格示意图如下，两个三角形网格，一个四边形网格。

因此，可以得到对应的全局与局部单元编号表格。

按照上述列表，对有限元离散方程进行组装得到全局系数矩阵如下：

4.边界条件

添加边界条件，在y=0的边界上依旧强加第一类边界，使用乘以大数的方法实现。在y=20的边界上为第三类边界，分析发现，不管是三角形网格还是四边形网格，最后落在y=20的边上，都是1维线单元，因此处理方式不存在差异。

因此，参考一维二阶叠层有限元的基函数，可以得到一维线性基函数的单元系数矩阵维：

将上述边界条件，加载到对应的全局系数矩阵的位置，得到最终的系数矩阵与右端项，最终的全局系数矩阵：

右端项：

求解上述线性方程组，即可得到求解结果。

5.结果精度对比

对于上述的最简单模型进行对比，分别使用理论解对比纯三角形网格、纯四边形网格、混合网格的精度结果。

a.纯三角形网格计算结果

b.纯四边形计算结果

c.混合网格计算结果

可以看出，纯三角形网格的最大误差0.8高于纯四边形0.45;而混合网格的精度最大误差在0.8左右，混合网格在网格分界面上的精度要略高于纯三角形的最大误差。将网格尺寸继续加密，得到相对精确的结果，混合网格计算结果如下：

电场实部虚部可视化结果如下：

6.复杂模型测试

对于非规则模型而言，纯粹使用四边形网格就很难模拟模型的边界，例如下列模型：

在斜边上，只有通过三角形网格才可以做到完全的拟合模型边界。因此，粗略对下半区域使用三角形网格，而对上半区域使用四边形网格，求得结果：

总结

1.本文简单介绍了混合网格的有限元实现过程，对于这种方法的性能优势的具体体现还有待在更加复杂的模型中测试。

2.混合网格的有限元实现是可行的，初步总结规律即为复杂模型附近使用能完全拟合边界的三角形网格，在其他区域采用四面体网格。

来源：实践有限元

二阶叠层基函数：二维电磁衰减数值模拟

简述在之前文章中介绍了二维情况下二阶插值基函数的实现过程，本文章继续介绍与之对应的另一种高阶基函数：高阶叠层基函数，高阶叠层基函数的特征在一维有限元的实现过程中已经详细介绍：叠层基函数允许在不同区域使用不同的阶数，也就是可以实现混合阶有限元。这是与插值基函数的重要区别。本次通过实现二维电场在介质中衰减规律的数值模拟，详细阶数二阶叠层基函数的实现过程，并给出常见单元系数矩阵，为后续介绍混合阶有限元打下基础。1.二维电磁衰减的边值问题该问题的边值问题与研究区域：对于均匀介质，解析解公式为：这里要注意电场Ey的物理意义，表示传播方向为y方向，振动方向与传播方向垂直，即x或者z方向，因为该问题是二维问题，因此只能表示其中某个方向。具体方向暂不讨论。2.有限元离散方程推导类似：介绍一个二维结构化有限元的常见边值问题，这里直接给出有限元方程为：在end边界上，是三类边界条件，二维积分落在边上降为一维积分，这部分使用一维的基函数获取单元系数矩阵，然后根据全局坐标组装到对应位置。3.二阶叠层基函数首先给出二维三角形通过面积坐标得到的线性基函数：二阶三角形叠层基函数则可以表示为：对应的单元三角形中编号为：可见由叠层基函数的性质，二阶基函数包含了一阶基函数，得到这两个基函数的梯度形式：4.单元系数矩阵推导推导过程类似：二维有限元实现详细过程：三角形网格二阶插值基函数这里给出梯度项高阶部分的详细推导过程：整理得到系数矩阵：对y方向的推导，只需要把上述a修改成b即可。对离散方程的第二项系数矩阵推导类似，这里直接给出推导结果：离散方程的第三项边界项，为一维二阶叠层有限元的基函数：至此，有限元方程中三个部分的单元系数矩阵就全部推导完成。5.系数矩阵组装以2*3网格作为示意图，如下显示：对应二阶叠层基函数的局部编号到全局编号的对应关系如下：按照上述的点、边映射关系，可以得到全局系数矩阵。将该映射关系与二阶插值基函数的映射关系做比较，可以发现两者在高阶的顺序上存在差别，这是因为基函数定义的顺序存在差别。此外插值基函数与叠层基函数的最大区别，插值基函数的高阶部分在网格上也有具体点，而叠层基函数的高阶使用棱边表示，不具体体现某个点的数值，体现的是高阶部分的精度。详细说明请参考：高阶插值基函数与叠层基函数细节对比。6.计算结果与精度分析a.1阶有限元5*10尺寸的网格下，均匀介质，电导率为1,频率1000Hz情况下，电场实部与虚部,误差图：b.2阶有限元5*10尺寸的网格下，均匀介质，电导率为1,频率1000Hz情况下，电场实部与虚部,误差图：从a、b的对比结果，相同网格下，二阶的叠层有限元的结果基本上比一阶有限元结果的精度高一个量级；c.2阶叠层有限元，10*20网格，介质体背景电导率1，中间存在高电导率100情况下，频率10000Hz的计算结果。下图为电导率的模型示意图：电场的计算结果如图所示：d.2阶叠层有限元，40*60网格，介质体背景电导率1，中间存在高电导率100情况下，频率10000Hz的电场分布：从c的结果分析，二维的优势既可以表示两个维度上的不一致，在该例子中，表示二维下电导率介质的不同，更能够模拟实际情况的电磁场；从d的结果，网格越多模拟的结果更加的光滑，也正是因为d图中的光滑结果，说明该模拟的电场应该是垂直与研究区域，也就是电场极化方向为Z方向，而不是Y方向。因为电场的法线方向在介质突变的分边界上是不连续的，不会出现如图中这种渐变光滑的情况。当电场极化方程为Z方向，不会出现电场振动方向垂直穿过分界面的情况，只会出现电场的切向方向穿过介质分界面，而电场的切向方向刚好又是连续的。所以判断电场极化方向为Z方向。7.总结a.本次内容通过实现了电磁衰减规律的模拟，详细介绍了二阶叠层有限元的详细过程；b.通过电场传播规律以及有限元求解结果，反推电场的极化方向，再次强调了解物理规律再选择合适的有限元求解方式的重要性；c.标量（节点）有限元，即未知数落在节点上，这种有限元在各个方向上均是自然连续的，这是标量有限元本身的特征，因此对于电场法线方向在介质分界面不连续的情况是无法模拟的。后续介绍矢量有限元则可以解决这个问题。来源：实践有限元