[拼音]：youxianyuanfa

[外文]：finite element method

把物体或结构整体所具有的域 V划分为有限多个被称为单元的子域V_n，以求得近似解的一种数值计算方法。划分整体的离散化思想，早在20世纪40年代就已经有人提出。50年代M.J.特纳等人首先提出对连续体进行离散化的分析法。近年来，在计算机的配合下，它已不同程度地推广应用到许多工程技术领域，成为处理工程计算问题的有效方法之一。

若单元是一根直的细长杆，可以为只沿杆的轴线方向有一个度量尺寸时，称为一维单元。图1所示的等厚度平面板被划分为28个小三角形平面板，它们很薄，其“中面”是二个度量尺寸的平面，称为二维单元。除三角形外，二维单元还可以是矩形、梯形、四边形或其他形状。若单元是长方形体、四面体、六面体、或其他形状的立体，具有三个度量尺寸，则称为三维单元。

离散后的物体或结构中相邻单元之间，通常只在若干点上相互连接，这些点称为节点；相邻单元在图形上相吻合的边不相接触，而只通过边上的节点传力。此外，由于计算工作的需要，单元内部也可设置内节点（图2）。所有荷载包括体力和面力都要移置到节点上；支座也要设置在节点处。节点可视为铰结或固结或其他状况的连接。如上所述，有限元法所分析的是与原物体或原结构的材料相同、形状接近而由许多个单元以一定的方式连接起来的离散物体，计算所得应力、变形等结果是近似的，但当划分的单元合理且数量很多时，即可得到非常逼近于真实的解。

在有限元法的应用中，最常用的计算法是以节点位移为基本未知量的位移法，另外还有力法和其他一些特殊的方法。采用位移法计算时，首要的问题是采用什么样的插值函数来建立单元的位移模式。

位移模式和插值函数

物体或结构离散化后，在单元中各点位移的变化可采用逼近于其原函数的近似函数来描述，这种近似函数称为位移模式，通常用插值多项式表示。如图1所示最简单的三节点三角形单元的位移模式，可表示为如下的线性插值多项式：

　 (1)

式中α_m(m＝1，2，…，6)是待定的参数，称为广义坐标。

取1、2、3三个节点为插值点，使在每个插值点上有：

　　(2)

式中x_i、 y_i为节点的坐标；u_i、v_i为节点线位移。由(2)式6个公式可定出6个广义坐标，于是可将(1)式所示的插值多项式改写为：

　　(3)

而位移模式列阵为：

式中δ^θ＝[u₁ v₁ u₂ v₂ u₃ v₃]^T，称为单元的节点位移矢量矩阵；/是二阶的单位矩阵；而

　　　　 (4)

式中A为三角形单元的面积；

(5)

在单元图中，三节点1、2、3的次序必须是逆时针转向。

由(4)式可见N_i(i＝1，2，3)是坐标(x，y)的函数，在数学上称为插值的基函数，在力学上由于它反映单元的位移形态，所以称为位移的形态函数，简称形函数。

在多节点三角形单元中，最常用的是在三角形123中三个边的中点增加三个节点4、5、6（图3）。在三角形内一点 P的位置，可用面积坐标L₁＝A₁/A、L₂＝A₂/A、L₃＝A₃/A表示。取 6个节点为插值点，位移模式可表示为完全的二次插值多项式，并类似于(3)式将插值多项式写为：

　　　　 (6)

式中的形函数可用面积坐标表示：

　　　(7)

(2)、(3)和(6)式中所示的节点线位移u_i、v_i称为节点参数。各节点的线位移的导数也可取作节点参数。单元中所有节点的节点参数的数目总和称为该单元的自由度。位移模式中广义坐标α_m(m＝1，2，…)的数目须等于或大于单元的自由度的数目。如图3所示6个节点的三角形单元，若1、2、3三个主节点各取如下6个节点参数：两个线位移u_i、v_i及其导数 ；而4、5、6三个副节点各取如下4个节点参数:两个线位移u_j、v_j和沿法线方向的两个导数，则该单元自由度的总数为3×6+3×4=30。位移模式可取完全四次插值多项式，此时广义坐标α_m(m＝1，2，…，30)的数目恰等于30。

更如图1所示三节点三角形单元，设若每个节点的节点参数都是u_i、v_i、、、、　(i=1，2，3)，则单元的自由度数等于3×6＝18，位移模式可取完全三次插值多项式：

(8)

式中广义坐标 α_m(m＝1，2，…，20)的数目为20，大于单元自由度数18。为了得解，必须增加两个包括广义坐标α_m的方程式称为约束方程式。

设有一多节点的矩形单元(图4)，称为一般拉族元，四边节点2(m+n)个，内节点(m-1)(n-1)个。取各节点为插值点，表示位移模式的插值多项式为：

　　　　　　(9)

(i=0，1，2，…，m；j=0，1，2，…，n)

式中插值的基函数（即形函数），一般用拉格朗日多项式：

N_ij＝L怢(ξ)L怹(η)　　　　　　　(10)

式中L怢(ξ)为拉格朗日乘子函数:

　　　　(11)

将上式中的_i换成η，i和m分别换成j和n，即得 (10)式中的L怹(η)。

对于图5所示的正方形单元，是最简单的4个节点的方形单元，当用线性拉族元的公式，得与四个节点相应的形函数为：

　 (12)

式中ξ_i、η_i是节点i的局部坐标，而该单元的位移模式为：

　　　 (13)

四节点正方形或矩形单元不能适应曲线边界和非正交的直线边界，也不便随意改变单元的大小，为此将正方形或矩形单元改变为图6a所示的任意四边形单元则可避免上述的缺点。

通过把整体坐标(x，y)变换为局部坐标（_ξ，η），图6b所示正方形单元的位移模式和形函数式分别如(13)和(12)式所示，同样适用于图6a所示的四边形单元，前者是为了后者而引出的单元称为基本单元或母单元，后者是计算时所实际采用的单元。

图6中所示的局部坐标（ξ，η）可依下列变换式转换为整体坐标：

　　　　 (14)

式中N_i(i=1，2，3，4)仍为如(12)式所列包含局部坐标ξ、η的形函数；x_i、y_i(i＝1，2，3，4)为四个节点的整体坐标值。

由于图6a、 b所示的两种单元的位移模式(13)和坐标变换式(14)中所用的形函数是等同的，所以实际采用的四边形单元称为等参数单元或同参数单元。

空间问题所采用的最简单的单元为图7所示的四面体单元，其位移模式类似于平面问题中三节点三角形单元的(3)至(5)式。

对于空间单元，可类似于图6b，取图8所示的八节点立方体作为基本单元，仿照式(13)和式(14)得位移模式和坐标变换式。

单元分析

如图7所示的空间四面体单元中，各节点上所承受的力U_i、V_i、W_i(i＝1，2，3，4)称为节点力，其矢量矩阵为

^θ＝[U₁ V₁ W₁…U₄ V₄ W₄]^T单元内的应力矢量矩阵为

单元内的应力和节点力与节点位移的关系式分别为

┿＝δ^θ　　　　　 　　　　　(15)

^θ＝δ^θ　　　　　 　　　　 (16)

以上二式中的和 分别称为单元的应力矩阵和刚度矩阵。所谓单元分析就是计算这两个矩阵以及将单元的体力和面力等效地移置到节点上。

根据弹性力学适用于小位移的几何方程和物理方程得如下的两个关系式：

ε＝δ^θ　　 　　　　　 　　　(17)

┿＝ε 　 　　　　　　(18)

式中ε为单元内的应变矢量矩阵；和分别为应变矩阵和弹性矩阵。单元的应力矩阵为

＝　　　　　　　　　 (19)

刚度矩阵和体力、面力移置到节点上的计量公式，可用如下最小势能原理导出。

在满足小位移几何方程和位移边界条件的所有允许的应变和位移中，实际的应变和位移必使弹性体的总势能为最小：

(20)

式中第一项为单元的应变能；=[uvw]^T，=[XYZ]^T分别为单元内各点的位移和体力的矢量矩阵；_b和分别为单元的某一边界面上的位移和面力的矢量矩阵。运用变分原理，得计算单元的刚度矩阵和体力、面力移置到节点上的节点荷载公式为

　　　　　　　　　(21)

　　　　　　　　　 (22)

　　　　　　　　　 (23)

对于平面问题，只须将以上三式中体积分转化为面积分，面积分转化为线积分。

对于等参单元，刚度矩阵公式和体力的移置公式只须将式(21)和式(22)中的dV变换为｜J｜dξdηdζ，式中

　　　　　　　　　(24)

此时积分符号应改为。对边界面上承受面力（如在ζ＝±1的ηζ面上承受匀布压力q）的荷载移置公式为

总体合成和边界条件的处理

把结构划分为许多单元，求出各个单元的刚度矩阵（简称单刚），再利用它们求出整个结构的刚度矩阵，称为总体刚度矩阵（简称总刚），才能达到求出结构中各节点的位移的目的。这种由部分求总合的工作（包括对边界条件的处理），称为总体合成，其方法有自由度编码法和变换矩阵法，前者的优点是在计算机上容易实现，缺点是当节点参数包括有线位移的导数时，则难适用。此时，当采用后者。

自由度编码法

设有包括1、2、…、5等5个节点的平面结构被划分为Ⅰ、Ⅱ、Ⅲ三个单元，根据结构的左竖边和下横边分别没有横向和竖向位移的边界条件，布置横向和竖向支承链杆(图9a)。

每个单元都有i、j、m3个顶点，6个节点位移u_i、v_i(i，j，m)，因而有6个自由度，编码为1、2、…、6，写作节点力F和节点位移 δ的下角标（图9b），于是各单元的节点位移矢量矩阵和节点力矢量矩阵分别为

δ^θ=[u_i v_i u_j v_j u_m v_m]^T=[δ₁ δ₂ δ₃ δ₄ δ₅ δ₆]^T

^θ=[U_i V_i U_j V_j U_m V_m]^T=[F₁ F₂ F₃ F₄ F₅ F₆]^T

式中u、U为横向的节点位移和节点力；v、V为竖向的节点位移和节点力;δ和F为不分横、竖而是按自由度编码的顺序来排列的节点位移和节点力。

对于整个结构，除去沿支承链杆方向的线位移为零不计外，共有5个节点位移，因而有5个自由度，编码为1、2、…、5，写作节点位移墹和节点荷载P的下角码(图9c)。

墹＝[墹₁墹₂墹₃墹₄墹₅]^T

＝[P₁P₂P₃P₄P₅]^T

单元节点的节点力包括加在该节点上的荷载 (即集中加在节点上的荷载、体力和面力移置到节点上的荷载)和其他汇交于该节点的单元加在该节点上的力，但对于整个结构来说，后者互相抵消，因而节点上只承受各类节点荷载的总和。

式中和分别为单刚和总刚。单刚中第r行第s列，总刚中第R行第S列的元素分别写作k_rs和K_RS。

总刚中的各元素都是由一个或若干个单元的单刚中的元素来贡献。单刚元素的下角标rs和总刚元素的下角标RS所代表的自由度两两相当时，该单刚元素k_rs便将贡献给K_RS。对照图9b和c可得单刚和总刚中两两相当的自由度的数码，如第Ⅰ、Ⅱ、Ⅲ单元中自由度的数码6、4、2与总刚中的数码2相当。总刚中的元素通过如上找出相当的自由度数码的方法，从若干个单元中相当的元素求和来一一算得，从而确定总刚。

变换矩阵法

设有一平面结构共有n个单元，把各单元的节点位移矢量矩阵δ、δ、…、δ按顺序排列(这些单元中相同的元素并未归并，都分别依序列入)， 得：

＝[δ₁^e δ₂^e…δ₃^e]^T　　　　　　　(27)

把各单元的单刚₁、₂、…、_n作为主对角线的子矩阵，其余为零，得如下的矩阵：

　　　　　　　　 (28)

若整个结构共有m个节点，δ₁、δ₂、…、δ_m为各节点的节点位移矢量矩阵，按顺序排列起来得

δ＝[δ₁δ₂…δ_m]^T　　　　　　　　 (29)

和δ两列阵之间的关系式为

＝δ　　　　　 　　　　 　 (30)

式中称为变换矩阵。

求出_s和后便可依下式列出总刚

＝^T_s　　　　　　　　　　　 (31)

列出了总刚和节点荷载矢量矩阵以后，并处理了边界条件，就可通过解式(26)那样的线性代数方程组把基本未知量──各节点的位移求出。于是各单元内的应变、位移和应力就可分别依前边列出过的公式算出。

收敛准则

有限元法的解答要求在单元的尺寸逐步取小时能够收敛于正确解答。

选取的位移模式必须满足如下的收敛准则：

（1）每个单元的位移一般包括本单元的形变引起的位移和本单元以外其他单元发生的形变而连带引起的位移两个部分。选用的位移模式必须能反映后一部分的位移，即位移与本单元的应变无关时，将使本单元产生刚体位移；

（2）每个单元的应变一般包括与各点的位置坐标有关的所谓变量应变和与该坐标无关的所谓常量应变两个部分。选用的位移模式必须能反映常量应变。

以上两条是能够收敛于正确解答的必要条件。

另一个条件是相邻单元之间的位移是协调的，即不出现图10所示的缝隙或重合。相邻单元位移的协调性是能够收敛于正确解答的充分条件。

特种方法

从有限元的基本方法派生出来的方法很多，如有限条法、边界元法、杂交元法、非协调元法和拟协调元法等，用以解决特殊的问题。

参考书目

1. 冯康、石钟慈：《弹性结构的数学理论》，科学出版社，北京，1981。
2. 徐次达、华伯浩：《固体力学有限元理论、方法及程序》，水利电力出版社，北京，1983。
3. 徐芝伦：《弹性力学问题的有限单元法》，水利电力出版社，北京，1978。
4. O.C.Zienkiewicz，The Finite Element Method，McGraw-Hill，London，1977.