我对有限元法的理解（修订版1.0）

Posted 2020-08-12 caicailiu

            温故而知新。                                          

                                                    --《论语》 

 

这里只讨论基于位移（注1）的有限元方法，它最终建立的是关于位移为未知量的方程组。也就是说，在推导过程中涉及到的其他未知物理量都要转化为位移的表达。

推导过程中用到的关键物理定律，虚功原理或最小势能原理。由于这两者在某些情形下是等价的，这里只以虚功原理为例进行说明。所谓的虚功原理，是说对于静态平衡的系统，力（注2）在虚位移上所做的功等于应力在虚应变上所作的功[1]。

本段来具体解释虚功原理中涉及到的概念和方程。先来看功。在标量的情形下，功等于力和其作用距离的乘积，而在向量的情形下，功等于力与物体位移（这两者都是向量）的内积，标量的结果可以看作它的一种特殊情形。值得一提的是，由于体力，面力以及外力和物体的具体形状都给定，所以虚功原理中的力可以看作是已知的。再来看应力和应变。直观地理解，应力可以看作物体的内力（注3）（鲍老师语），应变是在应力作用下引起的变形的一种描述，它包括物体内一点邻域内变形前后长度和角度的变化量[2]。在弹性力学中给出了应力和应变之间的关系，称为本构关系，它是有限元推导过程中一个必不可少的关系（注4）；由于讨论的是基于位移的有限元，如果能够给出应变与位移之间的函数关系（称为几何关系，推导过程中另一个关系），加上刚提到的本构关系，那么应力在虚应变上所作的功就可以表述为位移的函数。最后来看“虚”这个概念，虚位移指的是符合约束条件的任意的无穷小位移。

根据虚功原理列出相应的方程组后，经过移项等处理后，得到关于位移的函数乘以虚位移等于零的结果。注意到虚位移的任意性，得到关于位移的函数为零（这一点可以给出证明[3]，它也是数学分析中一个经典的证明，有兴趣的同学可以自己证明下）。这就是有限元要建立的最终结果。

细心的童鞋可能会发现，上面的叙述只是一个理论上的框架，其中某个未知变量关于位移的函数只是一笔带过。那么，具体到编程实现上，该如何给出相应的函数呢？首先，给出整个物体上关于位移（它是一个连续量）的函数在绝大部分情况下是不可行的。退而居其次，在离散的空间中进行定义，如果离散的程度很小，就可以近似连续情形下的结果。这是数值求解实际问题的一个通用的步骤，称之为连续问题离散化。具体到有限元中，就是将一个连续区域划分成彼此不相交的区域（二维情形下是三角形或者四边形，三维情形下是四面体或者六面体），称之为单元或者网格，单元之间通过节点相连，有限元最终建立的就是关于单元节点处位移的方程组。而单元中除了节点外其他的点处的位移，可以基于节点处的位移值进行插值得到。其次，即便将问题域离散化，那么采用什么样的函数（称之为形函数）来拟合单元内某点处的位移呢？最容易想到的，自然是多项式函数[4]，主要是它容易理解，计算简单，逼近的精度也能达到要求的优点。当然，函数系的选取要结合具体的问题。如果边界条件是周期边界条件，经常使用三角函数来进行逼近。现在有了单元内任一点的位移表达式，那么相应的应变和应力也就可以表达出来，自然地，应力在虚位移上所作的功也可以表示出来。注意的是，这里要使用积分的表达式，而且最终的结果中是关于位移的二次项表达。

再来看等式的另一端，目的是将之表达为节点处位移的函数，再根据虚位移的任意性，就可以将问题转化为关于位移的方程组。仅以体力为例进行说明，给定单元内任一点处的位移表达式，根据等功原理，即体力在单元上所作的功将其等价为集中力力在单元节点处所作的功，就将体力的功转化到单元节点处的表达。

主体框架搭建完毕，简言之，最终建立了一个关于单元节点处位移的方程组。剩下的一些关键问题，比如从单元刚度矩阵到总体刚度矩阵的构建（涉及到坐标系的旋转变换，计算机图形学的基础内容），线性方程组的求解，积分的数值运算，力的叠加等等，其中第一个问题是编程时必须考虑的问题，总刚矩阵最后的结果使用稀疏矩阵来存储，后两个问题是数值计算研究的主要内容。至于说力的叠加，对每个单元使用虚功原理，求出单元上每个节点相应的内力，添加到节点力的矩阵中。也就是说，具体的实现就成了一个数学问题。          上述问题如果不加任何限制的话，在数学上是有无数多组解的（如果是线性方程组，对应的系数矩阵奇异）。为了得到唯一的解，需要增加约束条件，通常在有限元法中对边界施加，称为边界条件。以传热问题为例，边界条件又分为边界节点的温度固定（第一边界条件），边界节点处温度的梯度固定（第二边界条件），或者两者的组合（第三边界条件）。不较真地可以这样认为，上述问题在数学上都是可解的。     尽管从理论上问题是可解的，但是从程序实现上仍然有很多地方需要注意。由于方程组的系数矩阵是一个大规模的对称稀疏矩阵，在存储时采用三个向量来存储，一个用来存储某一行的非零元素的个数，一个用来存储该行非零元素所在的列，一个用来存储对应的某行非零元素的具体数值[5]；在将边界条件考虑进去以后，对应的方程组也要相应地进行修改，对于第一边界条件，直观的想法自然是直接对其赋值，但是在程序的具体实现过程中，使用矩阵乘大数法；至于大规模稀疏线性方程组的求解，以前常常以迭代法为主，随着并行和多核的发展，目前直接法也得到了越来越多的关注；而针对大规模的非线性问题，则使用Newton-Rapshon方法将其转化为一系列线性方程组的求解。     感谢耐心地看到这里，最近一个关于有限元的边界条件问题仍旧困扰着我，一般的教材讲到有限元边界条件时，都是在问题求解区域的边界上施加限制， 那么能否在一部分边界上施加条件，一部分内部节点上施加条件呢？从程序实现的角度来分析是可行的，另一个棘手的问题是这样施加的条件，有唯一解吗？更进一步，施加什么样的条件，数学上能够保证问题有唯一解？希望大家不吝赐教。       总结：有限元法用到的知识点，虚功原理（刚度矩阵的构造和整体方程的构造），转换矩阵（局部坐标系和整体坐标系的转换），叠加原理（力的叠加，单元刚度矩阵构造整体刚度矩阵），边界条件（实际问题及唯一解的需求）。
    附：线性问题中三角形单元涉及到各个量的转换关系图[6]：

                  

 

 

 

注1：假设通过求解线性方程组得到位移的结果，那么根据几何关系，就可以得到应变的分布，再根据本构关系，能够得到应力的分布。

注2：这里的力包括体力和面力。体力，指的是分布在物体体积上的力，如重力等，它在有限元中对应于体积载荷，推导过程中要将其转化为对应的单元节点力，文中有述；面力，即作用在物体表面上的力，在有限元中对应于边界分布的情形，和体力一样，推导时也要将其转化为对应的单元节点力。         

注3：关于应力和内力的转换，依旧使用虚功原理来实现，只不过这个时候将其应用于某一个单元内，而不是整个求解区域。这样就可以求解出该单元每个节点的内力向量，它表达为一个矩阵（线性变形下为常数）乘以应力。节点内力和位移之间的关系通过单元刚度矩阵来表达[6]。

注4：在一维情形下，弹簧的受力（应力）拉伸或者压缩（应变，变化前后的长度差和原始长度的比值）满足本构关系，即f=kx。在二维情况下，应力和应变之间的转换，Omega=M*Delta，只不过这里应力和应变是向量形式，对应的M是矩阵形式。如果M是常数矩阵的话，那么相应的力学问题成为线性问题，如果M中的元素是关于位移的函数时，则该问题是非线性问题。  

 

参考文献：

[1] http://zh.wikipedia.org/wiki/%E8%99%9B%E5%8A%9F#.E8.99.9B.E5.8A.9F.E5.8E.9F.E7.90.86

[2] 《物理学与偏微分方程第二版（上册）》，第五章：弹性力学，李大潜等著，北京：高等教育出版社，2006年

[3] 《微分方程数值解法（第四版）》，李荣华等著，北京：高等教育出版社，2009年

[4] http://zh.wikipedia.org/wiki/%E9%80%BC%E8%BF%91%E7%90%86%E8%AE%BA    

[5] Helfenstein R, Koko J. Parallel preconditioned conjugate gradient algorithm on GPU[J]. Fuel & Energy Abstracts, 2012, 236(15):3584-3590.  4.1 Matrix storage                 

[6] http://wenku.baidu.com/link?url=wNatQMjgOjH-viTSBwPu-xNT1a46mPl9q0bf6logQQL0iXVl6urOlzEoQ0PwlwnB4LujELtOcR1PixcZY__TB2lMBFwVDVNagX7kKfHvcgu  (同济大学孙远韬)第7讲_三角形单元      P24 以及 Step3 单元刚度矩阵