FEM之单元(1)---三角单元介绍
本文可以当做有限元的小品文看，目的是对有限元理论进行小结。
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.常用标记：
上图截取自：
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
仍然可以利用前面的推导，得出任意一点场函数的表达式。
平面温度场方程为：
以第一类边界为例，带入泛函计算，最后可得泛函表达式：
公式的推导参考 《有限单元法在传热学中的应用》
针对稳态，第二项可以去掉。
平面电磁
场函数F取静电势/磁势,电场/磁场
考虑如下二阶微分方程：
常用的二维拉普拉斯方程，泊松方程和赫姆霍兹方程都是上述方程的特殊形式
利用变分建立该函数的泛函，将场函数带入泛函，最后可推导出单个单元方程
[K]{F}-B=0
下一步和静力学一样，组装成总体刚度矩阵，求解
在电磁计算中，使用三角面片网格会出现伪解的现象，而且由于是多种材料，在处理边界时会比较困难，于是引入了矢量单元，将自由度赋在边而不是节点上（edge element），具体参考《电磁场有限元方法》
平面声学
场函数F 取声压
声学中需要求解波动方程
右边的P为声压。对于简谐振动，可以改写成 赫姆霍兹方程（Helmholtz）
得出该方程的泛函，也可以通过能量平衡原理导出泛函。取第一边界条件，场函数带入可求得
[K]*F=a[B]
之后可求得声压和频率
平面流体
流体主要求解Navier-Stokes方程
理论上三角形单元是可以用来求解流体N-S方程的，实际上由流体的特点，有限元通常使用四边形和六面体。考虑到计算效率，目前CFD软件用的最多的还是有限体积法。
对于三角形二阶单元（每条边上有一个点，一个单元共有6个点）以及高阶，推导过程类似。工程上二阶单元已经能很好满足要求。有些商业软件提供了高阶单元（p单元），高阶单元的好处是只需少量网格，针对某些特殊case，可以在不改变网格情况下，通过升阶提升精度。
有限元是求解偏微分方程一种方法，如果想发开有限元程序，求解偏微分方程是绕不过去的槛。不过好在科学家，数学家们已经做了很多研究，不用我们自己去推导头疼的公式，只要记住公式和结论都行了，力学中二阶的四面体单元刚度矩阵有30*30=900个数据，实在没办法自己推导。但是作为有限元开发，了解推理过程对开发是非常有好处的。比如一般商业有限元程序中都会使用等参单元和高斯积分，积分点在单元内部而不在顶点，所以计算出的结果为高斯点的值，还有死锁，沙漏等诸如此类，这些都是与有限元算法本身特性相关。
声明：原创文章，欢迎留言与我讨论，如需转载留言
