加载中…
个人资料
多物理场仿真
多物理场仿真
  • 博客等级:
  • 博客积分:0
  • 博客访问:204,426
  • 关注人气:267
  • 获赠金笔:0支
  • 赠出金笔:0支
  • 荣誉徽章:
正文 字体大小:

FEM之单元(1)---三角单元介绍

(2014-10-16 11:06:49)
标签:

科学数值计算

多物理场

偏微分

有限元

有限元单元

分类: FEM有限元方法
本文可以当做有限元的小品文看,目的是对有限元理论进行小结。

1.一阶三角形场函数
三角形单元因生成容易,计算简单,容易加密,成为有限元分析中最常用的单元。
针对一阶单元,即一个三角形有三个点。
对于任意一个三角形,假设三顶点坐标点为(x1, y1) (x2, y2) (x3, y3), 三点的场函数分别定义为点Fx1, Fy1, Fx2,Fy2,  Fx3,Fy3
三角形内任意一点的场函数可以表示为
Fx=a1+a2*x+a3*y         
Fy=b1+b2*x+b3*y
其中Fx为x方向的场函数,Fy为y方向的场函数,a1, a2, a3, b1, ,b2, b3为常数项。

由于已知三点坐标,即在三点上也满足上式。将三点的坐标带入上式,可得6个方程。例:
Fx1=a1+a2*x1+a3*y1
a1, a2, a3, b1, ,b2, b3共6个变量,6个方程,可以求出:
a1=((x2*y3-x3*y2)*Fx1+(x3*y1-x1*y3)*Fx2+(x1*y2-x2*y1)*Fx3)/(2*A)
b1=((x2*y3-x3*y2)*Fy1+(x3*y1-x1*y3)*Fy2+(x1*y2-x2*y1)*Fy3)/(2*A)

其中A=((x2*y3-x3*y2)+(x3*y1-x1*y3)+(x1*y2-x2*y1))/2
公式看起来比较繁琐,这些转换的目的只有一个:为了让场函数能用 三角形的三个顶点的坐标来表示。
将求解出来的 a1,a2,a3,b1,b2,b3 带入到
Fx=a1+a2*x+a3*y         
Fy=b1+b2*x+b3*y
中:
可以得到新的表达式:
Fx=N1*Fx1+N2*Fx2+N3*Fx3
Fy=N1*Fy1+N2*Fy2+N3*Fy3
其中
N1=(x2*y3-x3*y2)+(y2-y3)*x+(x3-x2)*y
以此类推,详细推导公式可以参考任意一本有限元书籍
总之最后的结果是:
三角形内任意一点的场函数可以用 三个顶点的坐标场函数,以及该任意点坐标表示,这样一来,只要我们求出了顶点的场函数的值,就可以通过插值计算出三角形内任意一点的场函数值。如果是矢量,需要两个表达式,标量只要一个表达式。

2.偏微分方程概念
有限元方法的目的就是求解偏微分方程,利用有限元方法求解偏微分方程主要有两种:
1. 变分原理的Ritz方法
2. 利用加权余值中的 伽辽金法(Galerkin weighted residual method)

一些基本概念
1.边界:
第一类边界,直接描述边界上结果,用于已知边界确切结果 比如边界上外力 F = 10,温度T=40
又叫 Dirichlet (狄利克雷)边界

第二类边界 不直接给出确切结果,用导数方式给出
又叫Neumann (诺伊曼/诺曼)边界

第三类边界可以看成是 第一类和第二类边界的叠加
又叫Robin  边界

2.常用标记:
FEM之单元(1)---三角单元介绍
上图截取自:
http://blog.sina.com.cn/s/blog_5701b67c0100x7fv.html


3. 不同物理场应用
平面力学问题:
场函数Fx,Fy取为位移
由平面应变为位移对坐标的导数
即x方向应变 = Fx'/x
y方向应变 = Fy'/y
可知 Fx=a2, Fy=b3 都为常数。所以一阶三角形为常应变单元,即一个单元内应变不发生变化,为了保证精度,所以在物体应变变化大的地方要加密网格。
弹性力学偏微分方程中的应力,应变通过变换也都可以用 场函数来表示,推导可参考有限元书(目前市面上关于有限元大都是力学方面的)。
由物体平衡时,物体整体势能最小,将势能函数取变分可得到2D静力平面问题的矩阵表达式,表达式中包含了刚度矩阵,位移向量和节点荷载向量,根据整理的公式,就可以直接用代码实现了。

平面热传递:
因为温度是标量,所以推导求解比力学问题要简单,场函数定义为:
F=a1+a2*x+a3*y    
仍然可以利用前面的推导,得出任意一点场函数的表达式。
平面温度场方程为:
FEM之单元(1)---三角单元介绍

以第一类边界为例,带入泛函计算,最后可得泛函表达式:
FEM之单元(1)---三角单元介绍
公式的推导参考 《有限单元法在传热学中的应用》
针对稳态,第二项可以去掉。

平面电磁
场函数F取静电势/磁势,电场/磁场
考虑如下二阶微分方程:
FEM之单元(1)---三角单元介绍
常用的二维拉普拉斯方程,泊松方程和赫姆霍兹方程都是上述方程的特殊形式
利用变分建立该函数的泛函,将场函数带入泛函,最后可推导出单个单元方程
[K]{F}-B=0
下一步和静力学一样,组装成总体刚度矩阵,求解
在电磁计算中,使用三角面片网格会出现伪解的现象,而且由于是多种材料,在处理边界时会比较困难,于是引入了矢量单元,将自由度赋在边而不是节点上(edge element),具体参考《电磁场有限元方法》
平面声学
场函数F 取声压
声学中需要求解波动方程
FEM之单元(1)---三角单元介绍
右边的P为声压。对于简谐振动,可以改写成 赫姆霍兹方程(Helmholtz)
得出该方程的泛函,也可以通过能量平衡原理导出泛函。取第一边界条件,场函数带入可求得
[K]*F=a[B]
之后可求得声压和频率

平面流体
流体主要求解Navier-Stokes方程
理论上三角形单元是可以用来求解流体N-S方程的,实际上由流体的特点,有限元通常使用四边形和六面体。考虑到计算效率,目前CFD软件用的最多的还是有限体积法。

对于三角形二阶单元(每条边上有一个点,一个单元共有6个点)以及高阶,推导过程类似。工程上二阶单元已经能很好满足要求。有些商业软件提供了高阶单元(p单元),高阶单元的好处是只需少量网格,针对某些特殊case,可以在不改变网格情况下,通过升阶提升精度。
有限元是求解偏微分方程一种方法,如果想发开有限元程序,求解偏微分方程是绕不过去的槛。不过好在科学家,数学家们已经做了很多研究,不用我们自己去推导头疼的公式,只要记住公式和结论都行了,试想 针对力学中二阶的四面体单元刚度矩阵有30*30=900个数据,实在没办法自己推导。但是作为有限元开发,了解推理过程对开发是非常有好处的。比如一般商业有限元程序中都会使用等参单元和高斯积分,积分点在单元内部而不在顶点,所以计算出的结果为高斯点的值,还有死锁,沙漏等诸如此类,这些都是与有限元算法本身特性相关的。

0

阅读 收藏 喜欢 打印举报/Report
  

新浪BLOG意见反馈留言板 欢迎批评指正

新浪简介 | About Sina | 广告服务 | 联系我们 | 招聘信息 | 网站律师 | SINA English | 产品答疑

新浪公司 版权所有