有限元基本理论及工程应用：第三章最小势能原理和分片插值

1、第三章 最小势能原理和分片插值最小势能原理和分片插值是有限单元方法的核心内容之一，在这一章中我们将介绍最小势能原理的具体应用，同一力学问题的几种不同的表达方式及它们之间的联系；Ritz方法在单元内的应用；两种最常用的插值形式（Lagrange型和 Hermite型）。协调的位移型单元的收敛条件。3-1最小势能原理一个平衡问题，可以至少用以下叁种不同的方式加以描述：（i）平衡方程（ii）虚位移原理F(x,y)=0（iii）总势能取驻值（函数的极值问题）AB、 有限自由度系统质点系图3-1(a)为两个重分别为PA, PB的小球，由不计重量，弹性系数为k的弹簧相连，放置在光滑的曲面F(x,y)=0上

2、。该系统的平衡问题可由以下三种方法来描述：（）平衡方程（在平衡位置描述平衡）（a）分别取、为研究对象，建立一直角坐标系，如图3-1 (b)所示。由、两球组成的系统的平衡条件为：ABTNANBPBPATxyo其中，弹簧力以压力为正，即：（）虚位移原理 (在平衡位置附近描述平衡)系统的平衡位置可描述为 rA，rB或 x、y，x、y如图3-1 (c) 所示。两质点的虚位移分别为：AByxPBPATTrBrAo（c）图（b）r, r 或x,y,x,y由虚位移原理可得上式也称为微分型变分原理。进一步改写为(3) 总势能取驻值 （在全范围内描述平衡（寻找平衡位置）引入总势能函数由总势能函数取驻值，即：寻找

3、平衡位置。以上三种描述系统平衡的方法，对于有限自由度系统而言是等价的。然而，三种方法的特殊性是：平衡方程是在平衡位置上描述平衡问题；虚位移原理是在平衡位置附近描述平衡问题；而最小势能原理是在全局范围上描述（寻找）平衡问题。最小势能原理可以表述为：在所有满足给定位移边界条件和协调条件的位移中，满足平衡条件的位移使总势能取驻值，若驻值是最小值，则平衡是稳定的。作为一条基本原理，它的正确性应由事实加以检验，不必从理论上给以证明。但为了便于理解“位移边界条件“和“协调条件“的含义，下面通过三个例子对最小势能原理和平衡方程的等价关系加以验证。、无限自由度系统弹性体 （） 轴向受拉的直杆f(x)x, u图

4、设杆长为，截面积为，轴向分布载荷f(x)。x=0端固定，x=L端受端点集中力P。 设位移u(x)满足：（i）u(0)=0 (位移边界条件)（ii）u(x) 在O, L上连续（协调条件）(3-1-1)(iii) 使取最小值。若u(x)+u(x)为不同于u(x) 的另外一种位移分布函数，也满足上述的位移边界条件和协调条件，由(3-1-2)可推出而不难求得：因P(u) 取最小值，即的充分必要条件是：对任意满足(3-1-2)的u(x) 有(3-1-3)式(3-1-3)即势能驻值条件。若补充假定u(x) 存在、连续，则对（3-1-3）分部积分一次，并利用（3-1-2）, 可得到(3-1-4)(3-1-4

5、) 式对任意u(x) 都成立的充分必要条件是： （平衡方程）（力边界条件）由上述过程不难看出：由势能取驻值可以推出平衡方程。反之也对，说明两种描述方法在力学上等价。两种描述方法对边界条件的要求不同。用微分方程描述时，u必须满足：（平衡方程）（位移边界条件）（力边界条件）用最小势能原理描述时，要求函数满足位移边界条件 而力边界条件将作为势能取驻值的自然结果。当u(x) 使P精确取驻值时，平衡方程和力边界条件将精确地得到满足；当u(x) 使P近似取驻值时，平衡方程和力边界条件只能近似的得到满足。在这种意义下，力边界条件又称为自然边界条件(Natural Boundary Condition )，而

6、位移边界条件由称为强制边界条件（Essential Boundary Condition）。两种描述方法对函数的光滑程度（即可微性）要求不同。用微分方程描述时要求u(x)有连续的二阶导数（记作uC2(0, L)）。而用最小势能原理描述时，为了保证变形能存在，要求u(x)平方可积（记作uH(0, L)） H1为一阶导数平方可积的函数组成的函数空间，C2为二阶导数连续的函数组成的函数空间。nnnx, uP(x)(fx, fy)BDCAOy, v图nq(x)（） 平面应力问题正方形区域边长为a，厚度为t，受到体积力 (fx, fy), 边界AB固定。边界BC、CD自由。边界AD的法向应力为q(x)，

7、切向应力为p(x)。势能P和强制边界条件为：uAB=vAB=0利用以及格林公式（其中LX，LY为区域之边界的外法线n的方向余弦）。势能P的驻值条件为：(3-1-5)注意到沿边界，外法线n的方向余弦为ABBCCDDALX-1010LY0-101以及沿线AB：u=v=0则(3-1-5)对任意u,v都成立的充分必要条件为沿BC：y=xy=0沿CD：x=xy=0沿AD：y=q,xy=p关于边界条件的讨论以及u, v可微性的讨论与上例相似，这里从略。（） 梁的平面弯曲q(x)x Q图4Mv对于图3-4所示的弯曲问题，总势能和强制边界条件为势能驻值条件(3-1-6)对上式分部积分两次，并注意到由于必须满足

8、强制边界条件 v(0)=v/(0)=0 则有(3-1-7)使(3-1-7)对任意v(x) 都成立的充分必要条件是：（平衡方程）（自然边界条件）微分方程的阶数为。关于v、v 的边界条件为自然边界条件，关于v、v 的边界条件为强制边界条件。当用微分方程描述时要求v(x)有四阶的连续导数vC4(0,L)。用最小势能原理描述时，为保证变形能存在，只要求v(x)平方可积vH2(0,L)。本节讨论的三个例题，可作为维数不同，阶数不同的问题的代表。现把一些重要结论归纳如下表。其中，“微分方程阶数”是以位移为基本未知量来计算，“可微性要求”是对最小势能原理而言的。问题微分方程阶数强制边界条件协调条件可微性要求

9、杆的拉伸关于u的边界条件u连续u平方可积平面问题关于u, v的边界条件u, v连续u, v平方可积梁的弯曲关于v, v的边界条件v, v连续v平方可积3-2 Ritz法 (有限元方法的基础之一)由于有限单元方法可以理解为在单元（子域）内应用的Ritz法。Ritz法是一种求近似解的常用方法，它的基本步骤是：（i）选一组满足强制边界条件、协调条件和可微性要求的函数、n这组函数被称为基函数。（ii）假定近似解（试探函数trial function）的形式为nn(iii) 将试探函数作为近似解代入描述问题的能量泛函中，由泛函取驻值，即vL/4PvLvABxC图定出系数1n。从而得到近似解。下面以图3-

10、5所示的简支梁为例，求解在集中力P作用下的变形，分别采用不同的基函数来介绍Ritz法。解法1：取多项式形式的基函数：显然，1、1、2、2 连续，1、2平方可积，且取 为试探函数，,为待定常数，则描述这一问题的能量泛函可写成由能量泛函的驻值条件可得到一代数方程组解代数方程组可得试探函数中的待定常数点位移的近似值为解法：取正弦函数为基函数这样选择的基函数同样满足协调性、可微性以及强制边界条件的要求。可计算出：令可求得点位移该问题的精确解为1x103x102x1004x1ux0u2u3u4u10ux图x3214两种方法求得的点位移绝对值小于精确值。正弦的基函数，使支座处弯矩为零的条件（不属于强制边界

11、条件）也得到满足。尽管Ritz法本身并不要求这一点，但是第二个近似解的精度显然比第一个要好得多。不难看到，当选取多项式或三角函数为基函数时，由于这些函数可以求导任意次，协调条件和可微性要求很容易得到满足，我们要做的事仅为设法满足强制性边界条件。这在一维情况下很容易做到，但在二维、三维情况下，除一些形状规则的区域（例如矩形、园、扇形、长方体等）外却很难实现。因而古典的Ritz法的应用受到了限制。3-3分片插值形式的基函数和试探函数（解的收敛性与插值函数的选取关系很大）、 一维Lagrange型插值图3-6为一轴向受拉的直杆，截面积和轴向分布载荷f可以是x的函数。因而轴向位移u(x) 可能是x的复

12、杂函数。将区间0, L分成若干（这里是三个）子区间（单元），编号为 。取端点和分点为结点，编号为14。坐标为x1x4。(1)基函数1 当j=i0 当ji定义基函数14。满足(i) (ii) 设基函数在单元内是x的一次函数。(2)试探函数的形式取为基函数的线性组合根据i的定义显然有：u(x)是x的分段线性函数，且系数ui恰好代表结点i的位移值，相互之间是独立的。这样分段（片）定义的试探函数的一个显著优点是：强制边界条件很容易得到满足。例如u(0)=0的条件只要简单地令结点的位移u1=0即可以实现。而且允许我们在任何方便的时候（例如组装总体刚度矩阵时）引入这些边界条件。由于强制边界条件问题已经有了

13、妥善的解决办法，我们的注意力将转向协调条件和可微性问题。（）协调性和可微性上面定义的i(x)和u(x)都存在着“尖点“，光滑程度不高。但是：（i） i(x)和u(x)在单元内连续，在结点处也连续；（ii）i(x)和u(x)在单元内连续，在结点处可能不连续。但只有有限的跳跃量。在区间0, L上平方可积。i(x)和u(x)属于同一类型的函数。对于轴向受拉杆（二阶问题），u(x)满足最小势能原理对协调性和可微性的要求。由及uNjuiujx0x,uxixjijNixixj1xxixjui1xxixj图可求得u、u、u、u的值，从而得到一个近似解。（）Lagrange插值i(x)、u(x)都涉及这样一个

14、问题：由两个结点上的函数值在单元内确定一个线性变化的函数。图3-7为一个一般性的单元，两个结点i、j的坐标为xi、xj，假定单元内u(x)是x的线性函数由定出1、2，可得到其中Ni、Nj称为形函数，它们在单元内是x的线性函数，且满足uNluiujx0xlxixj1Nixixj1xxixj1xxixj图8x0Nj0xlxlul0每个形函数由分子和分母两部分组成，分子保证了一个结点的形函数在其他结点处为，而分母的选择则恰好使得这个形函数在自己的结点个取值为。掌握了这些特点就可用“凑”的方法得到形函数的表达式。如果在每个单元内在增设一个结点l就可以假定在每个单元内u(x)是x的二次函数。形函数也是x

15、的二次函数。若结点为i、j、l，则可以用“凑”的方法得出各形函数：试探函数图3-8是它们的图形。用插值点的函数值构造的插值函数通常称为Lagrange插值。q(x)xx3214v1x11x1rad21rad1rad4x1rad3x、 一维Hermite型插值2x13x114x图3-9为一根梁，横向载荷q和截面惯性矩I可以是x的函数。因而挠度v是x的复杂函数。梁的弯曲是四阶问题，试探函数v及v 应在0, L上连续。将0, L分为若干（这里仍是三个）子区间（单元），编号为。取端点和分点为结点，编号为14。图9（）基函数定义基函数1(x)4(x)、1（x）4(x) 满足(i) i(x)、 j（x）在

16、单元内是x的三次函数1 当 j=i0 当 ji（ii） 1 当 j=i0 当 ji(2)试探函数根据i、i的定义可知，v(x)是x的分段三次函数，且满足系数vi、vi恰为结点处v、v之值。这些值相互之间是独立的。这样定义的试探函数保证了v、v在单元内连续，在结点处也连续，满足四阶问题所要求的协调条件。由于试探函数是分段定义的强制边界条件v(0)=v(0)=0可以简单地令v1(0)=v1(0)=0得以实现。（）可微性i(x)、i(x)和v(x)是分段定义的三次函数，函数本身及一阶导数在0, L上连续。二阶导数是分段线性函数，在结点处不连续，但只有有限跳跃量，在0, L上平方可积。满足最小势能原理

17、对试探函数可微性的要求。如图3-9（续）所示。vv1xv1v2v3v4v2v3v4vx图3-9（续）（）Hermite插值定义i(x)、i(x)和v(x) 都涉及这样一个问题：由函数及一阶导数在结点处的值确定单元内的一个三次函数。图3-10为一个一般性单元。它的长度为L,结点为i、j。i为局部坐标原点,与x平行。设梁的挠度函数(试探函数)为hihjHiHjiiijixjjjxxx11rad1rad1图3-10则以结点坐标i=0、j=L代入则有定出4个常数代入挠度函数并整理成的形式，即hi、 Hi、 hj、 Hj的图形已示于图3-10。显然基函数的非零部分正是由它们拼成的。用插值点的函数值及导数

18、值构造的插值函数通常称为hermite插值。 3-4常应变三角元的理论依据qqt图2-8现在重新考虑图2-8所示的平面应力问题。仍把所研究的区域分为个单元。单元和节编号方案不变。、基函数1 当 j=i0 当 ji定义基函数1(x,y)9(x,y)满足（i）i(xj,yj)=ij(ii)i(x,y)在每个单元内是x、y的线性函数。、试探函数根据i(xj,yj)的定义，在结点上有系数ui、vi恰为结点上的位移值。、在每个单元内，u、v是x、y的完全一次多项式，可以由三个结点上的函数值唯一确定，可由第二章的（2-3-2）表示。显然在单元内u、v连续，且其导数为常数。、u,vxyij图3-11沿单元边

19、界（例如i、j边），u、v按线性变化，完全由这条边上两个结点上的函数值ui、vi、uj、vj所决定，故穿过单元边界时u、v连续（如图3-11所示）。但穿过单元边界时其导数却一般不连续，有有限的跳跃量，但在内它们平方可积。所选择的试探函数满足平面应力问题（二阶问题）对协调性和可微性的要求。可用Ritz法求近似解。、总势能将假定的位移场写成其中为一个218的长方矩阵，在一个单元内有列全零（一个单元只有6个结点自由度，而上式是18个自由度）。该问题的总势能函数可以表达为(3-4-1)4、势能取驻值即强制边界条件 可以通过对K、F划行划列或在组装K、F的过程中加以实现。5、单元刚度矩阵在一个单元内有列

20、全零，故k1有行、列全零。ki就是第二章中“补零扩充“后的单元刚度矩阵。（3-4-1）的推导过程给出了由单元刚度矩阵ki组装总体刚度矩阵的另一种解释：总变形能等于单元变形能之和。去掉ki中全零的行和列，可得到一个66的方阵k。若单元三个结点的编号为i、j、k，当单元内位移按线性分布时，由（2-3-3）有其中B的表达式为（2-3-4）。在单元内为常数。单元变形能单元刚度矩阵与（2-3-5）完全相同载荷q的等效结点力归结为计算下列两个积分单元单元有限元方法可以看成采用分片插值形式的Ritz法。由于试探函数采用分片插值形式。即使在区域形状比较复杂的情况下，强制边界条件也很容易得到满足。但所选择的试探

21、函数必须满足协调性和可微性要求。这是最早出现的关于有限元方法的理论论证。3-5收敛条件一般说来，用分片插值形式定义的试探函数很难做到与问题本身的真实解（精确解）完全吻合，因而有限元解一般都是近似解。我们希望在网格逐步加密、单元尺度无限变小时有限元解能收敛到真实解。为了保证收敛性，各单元内假定的位移场（试探函数）应满足以下条件（）假定的位移场在单元内连续。这是最小势能原理所要求的，也是很容易做到的。（）能够描述任何一种常应变状态（常曲率）。q(x)归纳出这一条件的理由是：当单元无限细分时，真实解在每个单元内都将接近某种常应变状态。如果有限元解能够无限接近总势能的真实最小值，随之也就逼近了真实解。

22、满足这一条件的另一个好处是，如果真实解在整个区域上为一常应变状态或者应力变化比较平缓，那么大尺度的单元即给出精确解或较好的近似解。（）包括足够的刚体位移模式。在假定单元位移场时不考虑位移约束条件（强制边界条件），这些条件在组装总体矩阵时再加以考虑。如果单元位移场缺少某种刚体位移模式，实际意味着加上了额外的约束反力。实现条件（）和（）是不困难的。对于轴向受拉杆，只要包含完全一次多项式即可。其中1为刚体平移，2x为常应变项。对于平面应力问题，只需包含x、y的完全一次多项式即可，即其中第一个括号内为三个刚体型位移：二个平移一个旋转。第二个括号内为常应变项；对于梁的弯曲，只需包含x的完全二次多项式其中

23、1+2x为刚体型位移，一个平移、一个旋转。3x2为常应变（常曲率）项（）协调条件对二阶问题要求穿过单元边界时位移（试探函数）连续。对四阶问题则要求穿过单元边界时位移及其一阶导数都连续。这是最小势能原理所要求的。可以证明，上述条件（）（）是有限元解收敛到真实解的充分条件。进一步研究表明：有些不满足条件（）的单元仍然可以收敛到真实解，这种单元通常称为非协调单元，它们不能再以最小势能原理为理论基础。非协调元能否收敛到真解的硬指标是：考查其能否通过“分片检验”。“分片检验”的概念将在第六章中介绍。对于协调元来讲，以最小势能原理为理论依据，保证能够通过分片检验。因此，条件（）并不是保证收敛性的必要条件，

24、而条件（）（）则属于必要条件。3-6其他形式的二维分片插值（Lagrange型插值）0kjnmx, uy, vlis图3-12一个具体的单元由三个要素所决定(i) 单元形状。(ii) 结点的配置和结点参数的选取。对于Lagrange型插值总是取结点处的函数值（即不取导数值）作为结点参数。0kjnmy, vli图3-13hprt(iii)插值函数（一般为多项式）的具体形式。这些多项式的系数要能够被结点参唯一确定。同时还要考虑收敛条件的要求。、高阶三角元图3-12为六结点三角元。结点i、j、k位于三角形项点，结点l、m、n位于各边中点。单元内假定位移场的形式为x、y的完全二次多项式可以由六个结点上的u、v值唯一确定。显然u、v在单元内连续，且包含了x、y的完全一次多项式，收敛条件（）（）得到满足，沿单元边界（例如ilj边）x、y是s的一次函数，故u、v是s二次函数，完全被这条边上三个结点处的函数值所决定，故协调条件（）也得到满足。图3-13为一个10结点的三角元，单元内假定的位移场是x、y的完全三次多项式。同样可以满足收敛条件（）（）。0y,v2b2ax,ujilk0yx图、矩形单元