什么是有限元

Chapter 1

什么是有限元 by caoer

1.1 PDE

有限元是一种求解问题的数值方法，求解什么问题呢？－－求解PDE(偏微分方程). 那么PDE 是做什么用的呢？－－描述客观物理世界。我想如果这两个问题搞清楚了也就明白了为什么要用fem,fem 可以做那些东西。 PDE 可以描述很多物理现象，电磁，流体，换热，diffusion, 力学，河床变迁，物种竞争，股票金融，等等等等。。。。乃至整个宇宙，当然也不是所有的物理现象都可以用PDE 描述，如微观世界分子原子的运动等等，所以我从来都不建议用有限元方法仿真微观物质现象的原因，但也有PDE 应用于微观位置如possion 方程来解析plasma 的物理现象，这在量子物理里面用统计的方法过于庞大，泼松方程反而使问题简单而且能吻合实验，这些都是题外话就不多说了。除了PDE 以外，ODE 同样也可以描述客观世界，但ODE 多用于控制系统，很有一些线性PDE 的解法也都是将PDE 转化为ODE 来做解析解的

1.2 求解PDE

有了PDE 以后，问题是如何求解并得到结果，首先要说明的是不是所有的PDE 都有解的，往往有解的PDE 才有实际工程意义。对于数值解法，常用的是有限差分，有限元和谱方法，还有蒙特卡罗法。有限差分出现的较早，计算精度相对较高，但是费时，且模型形状必须规则，边界条件处理困难，好处是可以比较方便的控制计算精度，适用于流体类的仿真。有限元方法效率高且满足精度要求，边界条件容易处理，得到了广大的应用，尤其是固体领域。谱方法由于可以采用FFT 方法的来求解，使得程序有着精度高，收敛快的特点，也克服了有限元条件下使用高阶插值方程计算费时的缺点，常常使用periodic boundary condition，但也有越来越多的算法使得一类二类边界成为可能，适合微观尺度的PDE 解，谱方法和有限元结合产生的谱元法取两者之优点，使得应用前景非常好。蒙特卡罗法不是基于弱解形式的，随机数的多维采样最终得到统计上的结果，多用于金融分析。咱这里还是着重有限元解PDE ，顾名思义，有限元将整个计算几何模型划分为很多小的单元（element ）, 每个单元的含有一定数量的节点(node)，具体单个单元有多少节点，有对应的不同算法与差值方程，拿一个简单的线性4节点平面单元来说，每个单元包含4个节点，每个节点有对应的variable 值，比如简单固体力学问题，每个节点就有对应的位移值，热力学问题每个节点就有对应的温度值，等等。然后单元内部的variables 就通过差值方法计算得出弱解(weak formation) 是建立在变分法基础上的，通过这个方法将strong form 的PDE 转换为weak form, 使得有限元的求解成为可能，具体如何推导weak form 可以参考一些有限元书籍，如果一本基础的有限元书籍没有介绍如何推导weak form, 那么可以考虑选择换一本了。推导所得的弱解形式仍需要通过计算机来计算，Galerkin 方法推导出含有求和符号的公式，在计算机中多半以loop 的形式来计算这个量，往往这个量就是stiffness matrix中的component ，这个loop 多半是围绕gauss 积分点进行的。公式中还会存在积分计算，有限元方法多用gauss quaradure 的方法来计算，精度一般可以满足。也就说一般简单的限元计算中存在两次approximation, 一次是Galerkin 一次是gauss, 这也是很多人在计算完以后需要进行validation 的原因，同时对于非线性问题newton 法本身就是通过计算误差来判断收敛的，固体有限元常常使用能量增量作为newton 求解器的判据。单个单元的stiffness matrix计算完成后，还需要将所有单元的矩阵装配为一个大型的矩阵，然后进行线性代数计算。这个装配是很有技巧的，因为一般情况下stiffness matrix 是一个很大的稀疏矩阵，０值往往可以省去以节省计算量。一个由10x10x10个8节点矩形单元组成的模型会有11x11x11个node, 如果是热力学问题变量是温度T, 每个节点上有一个自由度，组装后的的stiffness matrix会有(11x11x11)x(11x11x11)=1771561之大，随着单元数或变量的增加，计算是惊人的，这样

一个大的矩阵求解显然不能用常规的方法(gaussian elimination)，这就是各种这样求解稀疏矩阵算法存在的原因了，有效且快速求解线代方程组是一个好的solver 的标准。

1.4 后处理

其实对于最基本的有限元方法，求解得到的仅仅是variable 的值，如力学就是节点位移值，thermal problem就是温度值，流体就是位移速度加压力值，如果我们想知道应力或者应变怎么办呢？后处理系统里面个都会增加相应子程序来计算stress, strain, flux 等等。这也就是为何我们能看到各种各样的contour 的原因了，当然读者也可以自己加入计算各种量的子程序，如应变能密度什么的。关于什么是有限元就介绍到这里，仅仅是一些随笔和想法，具体的理论和推导需要自身实践与探索，本文行文仓促只是阐述自己对有限元的粗浅理解。有不对的地方还请指正。

下一章会谈及一些我曾经用过得软件。

Chapter 2 有限元软件的介绍与比较

有限元软件很多，有商业的，开源的，免费的，并行的，多物理场的，专业于某领域

【这里仅仅介绍一些笔者曾经用过的，遗漏的还请高手补充。】

2.1ANSYS

第一个学习的有限元软件，也是上世纪末本世纪初国内最流行的，为什么流行呢？应为ansys 的教程铺天盖地，所以大家都学习ansys 因为有教程上手快，那个时abaqus 的教程可谓是凤毛麟角，自然学习的人也就不多，当时组里导师都是用ansys ，自然我也就用了。如果一个有限元软件推广的好，那么他的教程一定要推广的好，fem 是一个专业性很强的软件，教程推广不够，销售自然不行。题外话，接下来说说对用他的感受。

当时学习ansys 是直奔apdl 去的，加上理论背景很弱所以学习过程异常艰苦。总得说来ansys 是一个很中规中具的软件，功能很强大，计算也很可靠，速度快，基本上固体力学的问题都能解决。作为一个工程软件还是很不错的。帮助文件很丰富，我也挑不出他什么缺点。只是感觉国内学习ansys 已经陷入了一个怪圈，用户喜欢就某一个应用技巧钻的很深，对于其基础的理论却忽略不视，这无不于大量ansys 教材的编写误导有关。是指算例，却不关心原理，悲哀。

2.2COMSOL （femlab ）

第一次用comsol 的时候它叫femlab2.0, 无奈的是只是作为一种兴趣学习的，因为时间精力问题最后也没有很深入，偶然机会１年后又开始认真的学习，逐渐体会到他的强大，我想说的是comsol 对于各种物理场的research 的确是一把利器，他可以任意的由用户输入ＰＤＥ, 计算结果并出contour, 这对于懂得有限元理论却又不想花大量的时间来编程的人来说，就是种bliss, 然而缺点就是速度太慢，这也是可以理解的，毕竟comsol 是将单元与pde 信息进行重组后计算，而不是很多商业有限元已经将Ｐde 固化在单元信息里了，但是其出色的weak form PDE模块是不可多得的有限元分析模块。还有其建模与后处理模块相当方便，比ansys 方便很多。至于计算的精确性笔者还抱有怀疑态度，细节就不提了。国内现在这个软件推的很好，大有赶超ansys 之势。值得一提的是comsol 的并行计算能力，采用的umfpack+mpich的构架，总体并行效率还是非常令人满意的。

2.3FEAP

这个其实是老大哥，虽然很多人不知道它。开源，不免费，但是很便宜，也有免费的学生版本，大家可以下载学习用。据说abaqus 和ansys 都是源于这套程序，abaqus 技术部里面的人都是feap 搞这个组里出来的，虽然是用fortran 编的，但是编程思路清晰，注解丰富，参考文档相近，实在是不可多得学习fem 的神器。不提供gui 的输入方式，命令的输入和apdl 和umat 方式很像，都是文本macro 方式。计算速度超快，带有后处理功能，提供强

劲的二次开发接口，可以自己编写子程序输出其他后处理软件如tecplot. 具有openmpi 并行能力，提供丰富的elastic, plastic ，hyperelastic 模块。缺点是之用于固体力学与热计算，当然可以自己开发其他单元，如电磁的。还有就是feap 只运行于linux 下，使用者需要知道如何make 程序。

2.4libmesh 和 deal II

之所以把这两个程序放在一起，是因为都出自于同一个大师。c++开源免费程序，非常的robust ，但是还有一些小bug 存在，可以应用于更广泛的领域如流体，强劲在于adaptive meshing, 和remesh, 从而可以解决很多常规有限元无法解决的问题，如singularity 问题。libmesh 虽然是后起之秀，但是发展异常迅猛，无奈只有那20多个例子作为教程，对于初学FEM 的朋友那就是灾难了，而且似乎固体不是libmesh 的主攻方向，如果要在libmesh 上开发一个壳单元，就必须对整个内部插值环境有很好的了解，只能说难度很大。需要有c++背景以及有限元理论知识，也是笔者遇到门槛最高的一个有限元软件，推荐给专业人事。同样是在Ｌinux 下运行的软件。

2.5 ALGOR_pipeline

这个软件在也挺有名，但是国内用的少，主要是固流耦合这快做得比较好，就我所知一些有名的飞机和流体机械厂都用他，我主要是用pipeline 这个插件，也许名字记得不对，就是用来计算整个管路系统的震动的，前处理中设定管路，就可以计算出各阶振型与频率，很方便，据说准确性挺好。其主模块用的不多，不做讨论了。

基本上笔者主要使用的就是这几种有限元软件，我想如果能精通上面任何一种软件，学习其他的软件都是一个很快的过程。如果你是博士生且有一个比较长得研究学习时间，推荐用feap 和libmesh, 或者是基于如PETSc/blas求解器这样的自编程序， 如果是硕士，推荐comsol ，本科毕业设计，推荐ansys. 当然每个人水平周围环境都不一样，不能一概而论，总之多学一点东西没有坏处，每个软件都有其深厚的背景。

Chapter 3 有限元的数值方法

这里先略过有限元的几何建模与网格划分部分，因为一来不是数值方法的主要部分，二来我对这块也不是很精通，所以直接从数值解法入手。这部分可以说是ch4与ch5的总体概览部分。也是fem 的核心。

3.1 单元离散化与jacobian

划分好网格后，就意味着单元编号以及节点都已确定下来了，其实这步在有限元分析的一开始就设定好了。拿平面单元为例，如果选用平面8节点矩形单元来计算热力学问题，那么每个单元有8个节点温度值，要知道每个节点的位置并不是规则均匀分布在实际有限单元上的，人们通常将节点的global 坐标转换到local 坐标上，以方便计算推导，之后再通过jacobian 转换到global 上，这点和连续介质力学的reference 转换是同一道理。但是不是jacobian 的建立是不是必须的呢？不是，如果所有的插值方程都是基于global 的，那么local 空间可以省去，jacobian 自然不需要了，这种做法不被推荐，应该编出来的程序可扩展性非常差，而且计算量较大，所以在大型通用有限元里面都是计算jacobian 的。

3.2 运动方程与各种矩阵

得到了单个单元的8节点test function 或 grediant of test function的积分值以后（积分的计算用gauss 法），就需要联系所有的单元，装配成为一个整体的矩阵，这个矩阵就是stiffness matrix, 如果是一个4单元简单正方形区域，那么装配好的stiffness matrix 就是一个21x21矩阵，以k 标记，除此以外，如果是瞬态问题，也会有质量矩阵存在，多半是一个对角矩阵，其值一般是test function的积分值，如果存在damping 那么还会存在damping matrix, 一个有限元系统一般只存在这3个矩阵，然后通过一个运动方程将其连立起来，即

M*u_tt+D*u_t+K*u=F,F是source 。这就是将pde 转化为弱解再转化为有限元方程的最后形式。这个矩阵绝对是一个稀疏矩阵，对于稀疏矩阵传统的gaussian elimination方法似乎就显得非常笨拙了，所以各种各样求解稀疏矩阵的求解器就出现了。

3.3 Ax=b求解.

有了 M*u_tt+D*u_t+K*u=F 方程，接下来要做得就是求解他，如果是一个2x2的矩阵，手工即可计算得到，实际应用中，往往都是上万阶的矩阵。也许有人会问了u_tt和u_t是怎么解得的？瞬态问题由于有了时间变量的加入，所有必须要有对应的时间积分求解器，二阶的有newmarks, HHT，energy conserved 方法，一阶的有crank nicolson, 向前向后euler 法等等，这点在求解器部分会详述，经过时间求解器的计算后，运动方程还可以化简为 K_tilde＊u=F的形式，也就是线代求解器需要解的最后方程，求解后，各节点的u 值得到。完毕。

3.4 多场耦合方法

如果有多个场存在怎么办呢？道理很简单，比如一个固体场一个温度场，两个pde 出来的运动方程是这样的 M1*u1_tt+D1*u1_t+K1*u1=F1和 M2*u2_t+K2*u2=F2,将两个场的运动方程线性进行叠加（线性耦合），得到 M1*u3_tt+(D1+M2)*u3_t+(K1+K2)*u3=F1+F2 ,再化简，得到 K_tilde＊u3=F3，使用求解器求解既得结果。这里只是简单的线性耦合，对于复杂物理现象的非线性耦合都是在这个基础之上进行变化。而且要注意的是，最后所得的K_tilde不一定是对称矩阵，也就意味着求解器必须要有解Unsymmetric matrix的能力。

3.5 边界条件加载

边界条件的加载，都是在所有单元矩阵装配完毕以后进行的，往往是variable 或是gradiant of variable的值被限定住了，对于低阶边界条件，如Dirichlet boundary condition, 常用的方法是增加一个惩罚量，这个惩罚量相当大>1e10，这样便使得u 被固定在人为设定的u0位置。对于natrual boundary condition, 也就是q_n的值如何确定，通常于变量的梯度值相关，如果q_n的表达式已知，那么只需要计算出q_n的精确值带入source 相即可，需要注意的是natural boundary condition要求比较精细的网格，否则会导致局部节点的计算失真。对于固体问题，q_n代表的是压力，具体边界上的力（或弯矩在板壳杆单元中）可以通过修改源相{F}来施加载荷。

Chapter 4 常用固体单元

记得自己刚学ansys 的时候就被自带的100多个单元模型搞的晕头转向，到底应该选择什么样的单元？这些单元有什么不同？这些问题一直困扰着我，其实把每一种单元都搞透是一件不太可能的事情，但有件事情可以做得就是对于几个大类的单元有个底层数值方法上的了解，这对于单元的选择与开发是至关重要的。这里仅以3d 单元作为例子。

有限元方法中以单个单元为基础，所有的单元都有相同的属性，这些单元组成了整个domain, 可以说了解了单元也就了解了整个模型，现代有限元软件都将pde 整合在了单元信息里面，还句话说，Pde 所推导出的弱解形式中非边界条件term 都体现在了单元中。

4.1 固体单元, 变量：位移u1,u2,u3

工程力学中最常用的单元，常用于连续介质的力学计算，一般分为small deformation 和 finite deformation 两个区域，记得ansys 里面有个大变形开关，一打开后计算变得惨不忍睹，所以选择small 还是finite 要看具体的问题, 此外对于不同的材料，如超弹，粘弹，弹塑材料都有不同的数学模型，这里选择的时候一定要慎重，不要想当然，最好要看看手册和理论，计算方法上常用的有displacement,mixed,enhanced strain方法，displacement 是最原始的，收敛效果不好，可以用于解决一般基础问题，推荐后面两种来解非常规材料。此外，对于同样的pde 模型，单个单元节点越多计算越精确，单元数目相同情况下，6面体单元比

4面体单元精度高。商业软件一般都考虑的很周全，如果是自编软件的话，就需要考虑很多问题了，如locking, buckling, convergence 等等，这些东西很有趣也很深奥，既需要深厚的理论功底又要有实际的编程能力，如果要做实际可用的大规模问题，对计算机cpu 的构架，cache 的大小，data layout 和编译器都需要很好的理解，才能编写出高质量的代码。所以编写一个small deformation 的固体单元对于科技发展的今天，已经是本科生的课堂设计实验课了。但对于大变形，多场耦合，粘弹性材料，多尺度，还是需要比较深厚的张量推导内力的，所以搞力学理论背景的要多了解了解计算机相关原理及建立高效编程习惯，搞计算机出生的要建立张量运算推导和物理现象的概念，那样就是有限元开发的有用人才，这句话我和广大有限元爱好者共勉。

4.2 热单元，变量：温度T .

是相对简单的一种单元，pde 本身也很简单，就是换热方程，节点自由度就是T , 由于温度T 没有方向性，所以T 只有1维量，计算量小。另外heat flux = gradient of T, 重点参看一下Fourier heat conduction equation.

4.3 壳单元，变量：u1 , u2 , u3 , θ1 , θ2 , （θ3）

有些具有曲面外壳的材料，其一个方向的变形要远远小于其他方向的。 也有一种大变形的壳单元，只含有5个变量, 少一个θ或是u, 但是有可能在计算的时候不收敛，所以在使用不同的shell 单元时，一定要多看看说明，遇到不收敛的情况也不用慌张，多从手册分析解决问题。壳单元很多，最著名的要属MITC 了，bathe 在84年开发的，影响至今还在，基本上目前所有的壳单元(除了degenerated solid method）都是基于MITC 的。壳单元（理论）可以说是近100年来固体力学工作者最令人神往的地方，因为起广泛的应用使得大家不得不对器感兴趣，小到人的眼球，大到屋顶，重要到导弹潜艇，到处都是壳啊。。。。壳太复杂了，以至于有很多人尽量避开他去研究板，要说目前最高深的板壳理论还是钱伟长大师的非线性板壳理论（个人理解），可惜在现在西方的文献里，大家都只提von karman nonlinear plate/shell theory, 钱老的研究可以说是提前的人类文明至少50年（还是个人理解），因为就目前来说nonlinear 板壳的应用还是凤毛麟角。将nonlinear shell/plate理论变成有限元单元的，目前笔者还没看到，等待某位大牛出现将其编出。

4.4 框架单元, 变量:u1,u2,u3,θ1,θ2,θ3

常用于结构计算，如弯曲，剪切变形。也就是材料力学里面各种梁的计算，这里pde 方程是关于u 的4阶grediant, 解的方法是将二阶gradiant of u 看作是Moment,3阶看作是shear stress ，这类单元通常是2节点单元，一头一尾相互连接。

索单元和板单元就不赘述了，可以参看相关资料，与框架单元大同小异。索单元甚至更简单，和热单元类似。还有一些薄膜单元点单元什么的应用很少，就不多说了。

看到论坛上很多人开口闭口二次开发，其实所谓二次开发很大程度上都集中于新单元的开发，如果你有新的pde ，要想计算实现它，就必须写出新的单元代码，从而实现计算。其他方面的二次开发无外乎一些后处理或是mesh 上面的新功能，而这些商用有限元基本成熟，所以再提二次开发之前，先想想现有的单元库里面有没有我需要的单元，如果没有再二次开发它，而这个过程是漫长而检苦的。

Chapter 6 有限元的未来

基本上有限元数值方法的最核心内容已经差不多说完了，对应不同的问题人们题出了很多很好的方法，但这些方法都是基于我之前所述的理论基础之上的。接下来谈谈我对有限元将来发展趋势的看法，或许能和打算或正在从事fem 的朋友有些共鸣。

6.1 并行计算

随着多核cpu 进入个人电脑市场，并行软件已经是大势所趋，那么并行fem 软件作为大型科学计算软件，并行趋势极为明朗，各大fem 公司分分研制推出并行版，但是并行版的计算准确性与效率还需要经受考验，由于并行机制的设定，每个cpu 返回的计算结果不同步，会导致latency ，矩阵的blocks 分解如何建立行至有效的ghost node/halos, 减少communication time, 对于各种并行构架的适应性，都是大型有限元并行化要解决的问题，不过当Nvidia 公司推出基于Fermi 构架的GPU 以后，科学计算界为之兴奋，但其是否能催生出有效的并行通用有限元程序还是一个挑战。而且目前个人电脑上并行计算的speedups 并不能提高很多，有2-3倍就已经很不错了。所以并行fem 仍需经历考验，但一定是大势所趋。

6.2 多尺度模型计算

高性能显微设备的诞生，人们对于微观尺度的物质越发感兴趣，fem 能否描述微观世界呢？纳米级材料，微流场，原子的diffusion,dislocation 等等，量子力学由于考虑的电子的作用，使得计算量过于庞大，分子动力学方法虽然已获得了一定成功，但是终归不如fem 来的有效和方便。描述不同尺度下的材料，俨然以成为众多学者们研究的主要方向，在坚凯围棋的第6版有限元书中增加了此方面内容，也有不少学者结合md 和mc 的方法，来建立多尺度仿真系统，都有一些进展。

动网格，边界元，无网格等等都是fem 的发展或联系方向，本人不是很了解，不赘述，欢迎补充指正。

推荐书籍：(建议按顺序阅读)

1. 《Introduction to the finite element method，3rd edition》，by Reddy

2. 《The finite element method for engineers， 4th edition》， by Huebner

3. 《Computational inelasticity》， by Simo

4. 《Introduction to the nonlinear finite element》， by Reddy

5. 《The Finite Element Method, 6th edition》， by O.C. Zienkiewicz