有限单元法科普：到底什么是有限单元法(Finite Element Analysis)？

作者简介：

雪城第一98k， 美国雪城大学PhD在读；高级摸鱼运动员；有限元科研工作者；PUBG端游爱好者。 

由于个人的原因最近一直在回顾有限元的一些内容，再加上之前一直觉得自己写作能力有很大的问题，也算是自己给自己一个总结，就想非常认真的写一篇有限元的科普文章，为广大结构工作者能在使用有限元程序的时候能更加好的理解软件到底在算些什么，顺便看看自己的写作能力还有没有救， 也希望广泛 地 接受各位大神给我提点建议或者说一起讨论讨论。

在我学习有限元的过程中，也是参考了各行各业或者不同网站的资源，首先我觉得何晓明教授上传到B站上的《有限元基础编程》是非常不错的一个系列课程，但是需要你有一点有限元和Matlab以及编程的基础。

知乎上的Dr. Stein大神的文章也是非常不错的一个有限元资料，但是他讲的非常数学也非常专业，需要多读读。

再就是如果能翻墙的话Allan Bower 的网站也是非常好的资料并且给了不同版本的代码，Brown University毕业的学生跟ABAQUS软件有千丝万缕的联系（当然也是听说），去年上课的时候教授是Brown University毕业的大神，他说比较早期的ABAQUS是Brown University的学生参与了早期的设计，因此编程思路上跟Allan Bower教授网站上给的代码非常相似，包括形成总刚度矩阵的方式等。

如果真的想好好理解有限单元法，最好再去看看数值分析或者微分方程相关的基础知识。后面如果有什么想起来比较好的资源在补充。

首先，从本质上来说，有限元是一种求近似解的思路，最简单的明了的了解就是对于定义域在[0,π]上，f(x)=sin(x)是一条曲线，如果想求f(x)在定义域上的面积，最简单的办法就是积分。

但是如果我不知道具体的函数表达式，仅仅知道若干个点在定义域上所对应的函数值（f(0)=0, …,f(π/2) =1, …, f(π)=0）, 只要知道的点的个数足够多，那么我也可以把不同的点连成直线然后求面积，也可以得到一个近似的积分值。

那么，有限元到底是怎么应用到结构分析上的，在我看来首先要理解一个应力（stress），应变（strain），位移（displacement用u来表示）以及位移的导数 之间的关系。

用最简单的例子来解释，假设一个轴向受力单元（1D-bar element），材料的弹性模量是E，截面面积是A，受作用力P之前的位置是x1，x2，受力以后的位置是x1’，x2’，如图1：

△   图1

待解的未知量u1=x1’-x1，u2=x2’-x2，由于作用了一个力在单元上，单元会有一个拉长Δu = (x2’-x2)- (x1’-x1)，再带入材料力学（我隐约记得是第二章）给出的应变的表达式：

假设这个单元非常非常小，也就是对等式的右边取极限，就可以将应变转换为导数的形式：

在一维线弹性单元，应力应变之间遵循胡克定律：

再将上式稍微调整一下位置，再把P换成一般常用的字母小f来表示，Govern equation（实在不知道应该怎么翻译）可以表达为：

当然这是最最最简单的一种情况，复杂一点的单元没准会在后续中介绍（如果有的话，手动狗头）。

上面这个方程确实很好解，只需要把f移到右边，两边同时积分，就可以得到线弹性情况下的精确解，但是比如涉及到稳定性问题（stability），那么Govern equation就会变成：

我隐约记得之前某个结构群里讨论杆端弯矩对杆件失稳（buckling）有没有影响，如果是按照我的理解那肯定是有影响，但是似乎最后有人说了一句没有影响，然后问问题的大哥就说好的，那就是不影响。

现在回想一下，似乎生活节奏变快以后大家会自动过滤跟自己想法不同的声音，然后选择性相信自己想相信的内容，行吧，扯远了，如果对稳定性问题有疑虑，我强烈（夹带私货）的推荐你们看一下W.F.Chen和E.M.Lui的《Structural stability: theory and implementation》。

回到上面的Govern equation，上面的等式实质上是一个x, y(x), y''(x)的微分方程，如果用待定系数法也可以得到精确解，但是解起来很麻烦解一道题差不多就用了好几页A4白纸。

因此，有些时候就需要提高计算速度，但是又可以得到相对精确的解，而且一定程度上最好能概念化结构化的交给计算机去做，有限单元法（Finite element analysis）就出现了。

首先，假设我们对于微分方程（Partial differential equation）：

且有边界条件：

的近似解解由一组基函数(basis function)的叠加组成：

为 基函数， 为待解的未知系数，除非我们假设的基函数跟实际解的形式相同（或包含关系），那么我们根据假设的解与精确解之间一定会有误差。

若将近似解带入原方程，原微分方程就会一些误差，用R来表示：

因此，加权余量法(Method of Weighted Residual) 要求我们若不能使微分方程的余量R=0，最起码要保证在定义域上的余量R尽可能 地 小。

为此，引入一个试探函数（test function）v，对微分方程的余量R与试探函数v的积在定义域上求积分，再令积分为0。

具体为什么这么做可以参考一下知乎上的大神Dr.Stein的系列文章：

变分法简介Part 1.（Calculus of Variations） - 知乎 (zhihu.com)

https://zhuanlan.zhihu.com/p/20718489

对于上式，我的理解是经过上面的处理，相对减弱了对于微分方程解的限制。

引用经典变分法的例子来说就是：

△   图2

假设一个小球从1点滚到2点，我们想要求小球最快下降的时间，那么时间既是坐标（x,y）的函数，同时也是小球运动路径的函数，那对于这个问题，精确解就是图2左边，两点之间直线最短。

但是，假设这个路径很难直接计算出来，那么我们稍微放松一下限制条件，小球反正1点2点定了，你要是不愿意沿着直线走那也行吧，但是只要走的别太离谱最后到达2点的时间大差不差就行，这样就可以更轻松的得到解，但是精度上可能大于第一种直接求得的解。

如果放到结构分析里来可以这么理解，一根杆件在P作用下的第一失稳模态（buckling mode）的精确解是y=Asin(πx/L)，但是如果我们假设他的解的形式是y=a0+a1*x+a2*x^2，也能得到一个大差不差的解，但是可能会略大于精确解，这种情况叫做stiffened effect（加强效应？我也不知道怎么翻译了，但是知道是啥事就行了）。

我个人的理解就是杆件一定会沿着它本身最不利的情况破坏，我们预设的情况如果能满足边界条件，那么就会求得一个略大的近似解。

OK，啰嗦半天就是希望大家能看懂我想表达的意思，再回到之前的Govern equation上，表达式中u(x)（trial function）是精确解，假设近似解是由一系列的基函数（Basis function）组成，即：

戴帽子的u表示的就是近似解，那么之前的Govern equation就可以重新表达为：

上式中对我来讲c一般是一个常数，因此材料的模量和截面性质一般是常数，就算涉及到材料非线性或者屈服，我们依然可以别的处理方法完成计算。

因此，上式的左边就可以重新写为：

再对其进行分部积分处理：

同时带入边界条件：

如果你好奇为什么这两项等于零，可以去认真的读一读我发在上面链接的文章，我们在这里讨论的是定端变分的问题，因此在定义域的端点上的变分一定会等于零。这里如果实在不好理解可以先放一放，回头再看可能就看懂了。

上面处理过最终的形式就化简为：

这样处理的好处就是让Trial function(u)和Test function(v)的derivative order（导数阶？）趋于一致。

在之前提过Test function(v)是由我们定，那好，我直接令他等于近似解 ，这就是传说中的伽辽金法 (Garlekin method) 。那么弱形式的终形态就变为：

当然，除了弱形式，还有强形式(Strong form)。由于弱形式用的比较多，这里就主要介绍弱形式。

到这里，如果带入我们之前的最早例子里面的参数c = EA，f = P，那么方程的左边就是应变能，右边就是外力做功。可以假设每一个单元都是一个小弹簧，结合回顾初中学胡克定律时候的内容，是不是感觉一切都联系起来了。

再将近似解带入我们之前假设的基函数的形式：

这种写法看着还是有点迷惑，但是作为仅仅是一篇科普文，那我换一种更加直接了当的说法：

上式中 是跟 相关， 是常数，那么我们的近似解就可以分离为一组常数（节点位移）乘以 的形式，这里的 就是传说中的型函数（Shape function）。

如果把 提出来，左右两边同时约去，再把求和号拿出来，最后的形式就可以表达为：

这样就把问题转移到求 上，注意一下， 是跟材料截面相关的常数， 是我们的未知量。重写成一种更为熟知的形式就是：

这样是不是亲切多了（手动狗头）。

在此再说一下型函数，型函数是我们根据已知推未知的一种方式，假设已知我们节点值，但是节点中间的值我们并不知道，因为我们已经将我们的问题离散化了，也就是把一个整体划分成很多小的单元(mesh)。

你可以理解成，已知两个人身高，第三个人的身高介于两人之间，你就可以根据Shape function大差不差的猜一下第三个人的身高。这听上去很有道理，如果你已知某个人身高介于姚明易建联之间，求出来的误差可能在几厘米，但是如果告诉你某个人身高介于姚明和郭敬明之间，你是不是想弄死我，因为我说了一句废话，而且可能误差也非常大，手动狗头。

故事从另一个角度告诉我们网格密度的正确性：如果两个节点之间比较接近，那么节点值之间的差值较小，那么你使用插值的误差就较小。

关键问题来了，怎么猜？你可以假设姚明和郭敬明之间的第三个人遵循一次函数之间的关系，也就是说从姚明头顶连一条线到郭敬明的头顶，那就可以根据离郭敬明近还是离姚明近来预测一下身高。那如果你假设姚明和郭敬明之间的第三个人遵循抛物线的关系，那你就需要谨慎了，特定情况下计算出来的身高可能比我郭还要矮。

故事从另一个角度告诉我们两个道理：

1、当你的计算结果看着不是那么对，需要考虑是不是单元选的有问题。

2、所有的有限元计算结果需要结合实际的物理或者结构常识来对比。

让我们回到一维轴向受力单元问题上：

△   图3

假设一个轴向受力构件，受一个大小为P的外力，我将其离散为3个单元4个节点，图3所示。我们的未知量就是在每个节点处的位移(u1, u2, u3, u4) ，边界条件 (Boundary Condition) ：假设最右端是固定端 (u1 = 0) ，求解其余三个未知量。

在这里我假设每个相邻的点的位移之间相互有联系，而且是线性的关系，线性关系非常好理解，在单元1中的位移u(x)由其左右两边的节点位移u1和u2线性插值得到，你可以理解成一个在节点1处u = u1+0*u2，但是随着向右边移动，直到节点2 u = 0*u1+u2，如图4所示：

△   图4

我 们假 设每个单元的长度是h，在此h=L/3，在此处i是型函数的编号，j是节点的编号，对于第一个单元来说，非零项仅存在于i,j=1,2。

将型函数带入单元刚度矩阵的表达式，就可以化简为：

对 于单元2，i,j=2,3，带入单元刚度矩阵

同样，单元3，i,j =2,3

△   图5

关于矩阵组合(matrix assemble)，如图5所示，如果实在是记不清了，麻烦回去看一下矩阵位移法。

终于快到最后了，对于等式右边在节点123处，没有外力，所以：

此时将计算结果带入全局刚度矩阵 (Global Stiffness Matrix)和外力矩阵：

如果这个时候尝试去求解矩阵，那么一定会提示你刚度矩阵不可逆，因为没有施加边界条件。

大家可以想象一下，一根轴向受拉构件，没有约束，受P大小的作用拉力，那岂不是会一直做刚体运动(rigid body movement)。因此，我们需要处理一下边界条件。

矩阵的本质就是由一系列的方程组成，如果想给节点1赋值，可以令u1的系数等于1，u2, u3, u4的系数等于0，然后令结果等于1，那么最终的矩阵就会变为：

后续剩下的内容就是非常简单的线性代数运算了~