首页/文章/ 详情

泊松方程三维二阶叠层有限元实现-四面体

2月前浏览729

简述

在高阶有限元中,有两种高阶基函数形式:一种是插值基函数,一种是叠层基函数。高阶插值基函数即使在上篇文章中介绍,它是依据线性基函数组装而成,形成的每个基函数在单元上都有具体的点一一对应,并且每个基函数均满足插值基本原则,细节参考:泊松方程三维二阶插值有限元实现-四面体

本文继续介绍二阶叠层基函数,它也是由线性基函数组装而成,特点是高阶基函数具有叠层的规律,即高阶基函数包含低阶基函数,并且其高阶部分的基函数在单元上并没有与之对应的点。这是与高阶插值基函数本质的差别。也正是由于叠层基函数的高阶包含低阶的规律,这使得高阶叠层基函数能够实现混合阶数有限元,而插值基函数则不能。

本文通过泊松方程的边值问题为例,实现三维的高阶叠层基函数有限元。泊松方程的有限元理论推导不再赘述,参考:泊松方程三维结构化有限元实现初探

1.四面体二阶叠层基函数

已知线性基函数表达式(推导可参考:泊松方程三维结构化有限元实现初探):

根据四面体示意图可以发现,5~10号位置也是使用棱边表示,但仅仅如此,并不意味着如插值基函数一样,在棱边上有具体的点。同样,棱边编号的顺序如下;    

由此给出二阶叠层基函数的表达式:

可以发现,二阶叠层基函数的N1~N4与一阶基函数是一样的,5~10阶的则是由对应棱边的两个节点基函数乘积而得。这可以理解为高阶项仅仅表示线性基函数组合的高阶,就好比泰勒级数展开的高阶项。

此时,四面体内任意一点p的插值公式同样为:

对应的基函数梯度项可以表示为:    

2.单元系数矩阵推导

观察二阶叠层基函数,仔细对比二阶插值基函数,我们发现有很多相似的地方:

除了低阶项完全不一样之外,N5~N10的形函数方式组装一样,仅仅相差了四倍而已,因此在系数推导的过程中,我们可以直接利用二阶插值基函数相同的结果。此同样把10*10的单元系数矩阵分成9个亚矩阵:

K11的推导非常简单,同时可以发现,K11与一阶基函数的系数矩阵是一样的,这正是由于二阶叠层基函数包含了一阶基函数的原因。    

K22,K23,K33系数矩阵与二阶插值基函数中的一致,仅仅相差16倍,因此可以直接获得(详细参考:泊松方程三维二阶插值有限元实现-四面体):

剩下的仅有K12,K13两个块矩阵需要推导,首先推导K12:

根据上述结果,很容易发现,由于N1~N4的偏导数等于常系数a,因此积分就非常容易获得,根据规律很容易得到K12以及K13项的结果:    

将所有亚矩阵整理,得到对x方向偏导数的10*10系数矩阵,对y、z方向的系数矩阵只需要将系数a改成b、c即可获得,如下:


右端项的推导相对也是很容易:

3.系数矩阵组装

使用的网格同样是tetgen剖分得到二阶网格,其中的注意细节参考:泊松方程三维二阶插值有限元实现-四面体。此外还需要注意,虽然生成的节点文件node包含棱边中点信息,实际这是不需要的,需要的仅是elem文件中对棱边顺序的排序。

边界条件加载与插值基函数也存在差异,因为叠层基函数的高阶项不表示网格上具体 位置的物理数值,因此在加载第一类边界条件的时候,低阶项会强加第一类边界,高阶项强加为零。    

例如,假设第一类边界条件如下,则高阶叠层基函数的第一类边界的加载方法:

完成组装系数矩阵和添加完边界条件后,最终组合成完整的系数矩阵方程可以表示如下,求解线性方程组即可获得数值结果。

4.结果展示

网格量2234,四面体顶点数目:593,顶点数+棱边中点数=3847。使用相同的网格求解的二阶叠层有限元泊松方程与上篇文章的二阶插值有限元对比如下:

上图中对比结果显示,一阶部分(红色 区域)二者的结果是一致的,而高阶部分完全不相同。再次说明了二者虽然均是二阶能提高计算精度,但是由于基函数不同,求解的数值在高阶部分也不同。

          

从图中可以发现,使用同等网格,二阶叠层有限元的结果与二阶插值基函数的结果可以比拟,从三维显示图的光滑程度上看也同样优于一阶基函数。

5.总结  

1.使用高阶叠层基函数实现的泊松方程的结果与高阶插值有限元结果具有同等的精度。

2.实现高阶叠层基函数需要注意的点:i 高阶部分仅使用棱边编号表示,并不代表某个点的物理数值;ii 第一类边界条件的加载也需要注意高阶项强加为零。


来源:实践有限元
理论
著作权归作者所有,欢迎分享,未经许可,不得转载
首次发布时间:2024-12-06
最近编辑:2月前
实践有限元
硕士 签名征集中
获赞 0粉丝 2文章 58课程 0
点赞
收藏
作者推荐

使用二维矢量有限元展示一个实际例子

简述二维电场从空气中入射地下空间。电场从高空中垂直入射到地面,电场在传播过程中遇到电阻率不同的物体会反射电磁波,通过在地面监测电阻率可以获得地下空间的实际电阻率变化情况。假设电场的振动方向是x方向,传播方向为y方向,因此在OXY平面,要模拟电场的传播规律,选择矢量有限元实现。具体实现参考:二维矢量有限元的详细实现过程1.基本信息频率:0.1Hz~1000Hz选择13个频点;视电阻率公式针对二维三角形网格矢量有限元,其插值任意一点公式为:2.网格整个仿真区域分为研究区域与外扩区域:研究区域是感兴趣的区域,实际做研究的区域;外扩区域是为了保证电场衰减下在边界面最小反射而添加的区域,外扩区域的范围大小与所研究频段的趋肤深度有关。趋肤深度,是定义电场在良导体内衰减会衰减1/e,即36.8%时,所深入的距离。本例子中的最大频率、最小频率的区域深度:因此,区域深度最好在0.1Hz的趋肤深度左右,如此在底面造成的反射才有可能尽量规避;同时在地表面为了保证1000Hz的精度,研究区域的网格尺度小于1000Hz的趋肤深度。由此获得符合实际物理情景的网格:3.仿真案例a.地下空间的电阻率为均匀电阻率100欧姆米。可以发现,在地表获得的视电阻率结果基本上保持在100欧姆米,与理论预期一致。在高频下的结果102欧姆米也基本上保证在精度范围内。b.在100欧姆米的地下空间存在一个10欧姆米的低阻物体,在实际情况中可能是矿等异常的资源。c.在100欧姆米的地下空间存在一个1000欧姆米的高阻物体。d.在100欧姆米的地下空间存在一个1000欧姆的高阻物体,一个10欧姆低阻的复杂情况。几个模型的结果基本上都能通过视电阻率反应地下异常体区域的电阻率情况。总结a.有限元在实际运用中的建模必须得考虑具体物理规律,保证仿真求解的结果能反应实际物理规律,这需要专业的知识与有限元算法特征的配合。b.可以发现本案例中虽然使用的三角形网格,但是是通过结构化网格实现的,造成了很多冗余的质量差的网格,实际上可以通过三角形网格生成软件生成更加优质的网格。来源:实践有限元

未登录
还没有评论
课程
培训
服务
行家
VIP会员 学习计划 福利任务
下载APP
联系我们
帮助与反馈