有限元基本理论及工程应用：第五章等参数单元（第一部分）

1、第五章 等参数单元(Isoparametric Elements)在前面的章节中我们已经认识了三角形单元和矩形单元。这两种单元的边均为直边，用直边单元离散曲边的求解域势必要用更多的单元数才能较准确地描述实际边界。本章将要介绍的等参数单元是目前应用最广的一类单元，可用这类单元更精确的描述不规则的边界。这类单元的出现不仅系统的解决了构造协调位移单元的问题，而且自然坐标系的描述方法也广泛为其他类型的单元所采用。等参数单元在构造形函数时首先定义一个规则的母体单元（参考单元），在母体单元上构造形函数，再通过等参数变换将实际单元与母体单元联系起来。变换涉及两个方面：几何图形的变换（坐标变换）和位移场函数的

2、变换，由于两种变换采用了相同的函数关系（形函数）和同一组结点参数，故称其为等参数变换。5-1四结点四边形等参数单元、 母体单元自然坐标和形函数(1,1)(-1,-1)342图5-1母体单元：边长为的正方形，自然坐标系，示于图5-1。取四个角点为结点，在单元内的排序为、。仿照矩形单元，可定义出四个形函数(5-1-1)显然有如下特点：（i）是，的双线性函数（ii）(iii)(5-1-2)、 实际单元与母体单元之间的坐标变换（） 坐标变换设xy平面上的实际单元e由母体单元经过变换F得到，即 且规定结点（i，i）与结点（xi, yi）对应（i=14）。这样的变换不只一个，利用（5-1-1）定义的形函数

3、即可写出这种变换中的一个(5-1-3)x,uy,v3241=-1=1=-1/2=0=1/2=1=1/2=0=-1/2=-1图5-2（5-1-3）所定义的变换有如下特点：x, y是，的双线性函数。沿母体单元中常数的直线（坐标线），x, y是的线性函数，对应于单元e中的一组直线，特别，单元e的一组对边1-2、3-4为直线。类似，中常数的另一组坐标线对应于单元e中的另一组直线。特别，e的另一组对边2-3、4-1也是直线，单元e为直边四边形。单元的其他直线（例如对角线1-3），变换到单元e中将是一条曲线（图5-2）（）Jacobi矩阵Jacobi行列式矩阵(5-1-4)称为变换的Jacobi矩阵。de

4、tJ称为变换的Jacobi行列式。一般情况下，J的元素和detJ都是，的函数。若detJ恒不为零（一般使它恒正），则J-1存在，变换F存在逆变换F-1。使单元e内的任一点（x, y）对应于单元内的一确定点（，）。此时称变换F为非奇异的。detJ称为变换特征量。detJ还具有明显的几何意义，如图5-3所示。设在（，）处detJ0在（，）附近取一边长为d，d的长方形。设此长方形与单元e内的一个小子区域d对应，可以证明，此小子域的面积d在略去高阶微量后有ed(,)dd图5-3(x,y)例图5-4所示的实际单元e为边长分别为2a、2b的矩形。结点坐标为：2a 2b 1 (c,d)4 2 3 xy 图5

5、-4则由（5-1-3），可得出坐标变换为同样得到：表明：当实际单元e为矩形时，经坐标变换得到的x, y是，的线性函数。Jacobi矩阵Jacobi行列式在单元内是常数。当结点序号按图5-4的转向排列时，detJ恒正。、 单元内假设的位移场对于平面问题，设沿总体坐标系的位移为u、v，结点（xi, yi）的位移为ui，vi实际单元e内的假设位移场（Trial function）取为(5-1-5)注意，这里u、v虽然是用点的自然坐标，表述的，但位移u、v (以及后面的单元刚度矩阵)却是对总体坐标系的。这与第二章中在单元局部坐标系下定义位移场的作法有区别。在坐标变换（5-1-3）和假定的位移场（5-1

6、-5）中使用的是同一套变换关系(形函数)，同一套变换参数（与（xi, yi）对应的结点位移(ui,vi)）满足这一特征的单元称为等参数单元。这样定义单元有不少优点，但也对我们提出了一些新问题。假定的位移场是，的双线性函数，当实际单元为矩形时，可表示成x, y的线性函数，假定的位移场u、v是x，y的多项式。但对一般单元而言，不能表示成x，y的多项式，因而位移场u、v不再是x，y的多项式，不能直接利用第四章的结果进行收敛性分析。、 收敛性分析（） 单元内位移场连续x、y、u、v都是，的双线性函数（连续函数）。只要Jacobi行列式detJ0，u、v就是x，y的连续函数。即在实际单元内u、v连续。（

7、） 刚体位移和常应变条件对于二阶问题，这个条件归结为假定的位移场中包括总体坐标的完全一次多项式。或者换一个提法：假定的位移场可以精确地表述任何一种线性变化的真实位移场。当试探函数直接用总体坐标的多项式描述时（像第四章所做的那样）采用前面一种提法是方便的。现在试探函数是用自然坐标表述的，则用后一种提法更合适一些。我们定义的形函数满足：(5-1-6)设真实位移场为x，y的线性函数将x，y按（5-1-3）代入，并利用 有注意到 （结点位移的真实值）y,vx,ueeMs图5-54则有类似有上述论证表明：只要所定义的形函数满足（5-1-6）（不管形函数的具体表达式如何），且坐标变换和假定的位移场使用同一

8、组形函数（等参数单元总是如此），那么这样假设的位移场一定能够精确地表述任何一种线性位移场，即刚体位移和常应变条件总可以得到满足。（） 协调性对于二阶问题要求穿过单元边界时位移连续。如图5-5所示，考虑一个实际单元e，它的母体单元为。以1-2边为例。沿1-2边常数，x、y、u、v都是的线性函数。设e边界上的M点与边界上的点对应，则M到结点的距离将是的线性函数。反过来也是的线性函数，因而u，v也是的线性函数，完全由这个边界上两个结点1、2的位移值u1、u2、v1、v2所决定。从另一相邻单元 e 看来，沿边1-2, u、v也是的线性函数。完全被结点、处的位移值所决定。从单元e和e 看来沿共同边界1-

9、2上的位移处处相同，即在边界上位移是连续的。对其他边界可用类似的方法加以证明。四结点四边形等参元的形状有较大灵活性，巧妙地解决了单元形状的灵活性和收敛条件（主要是协调条件）之间的矛盾。但是一般的四边形单元只能精确地再现线性变化的位移场，有限元空间Sh的次数k1=1。虽然能保证有限元解的收敛性，但精度不够满意。当实际单元是矩形时，是x、y的线性函数，假定的位移场将是x、y的二次多项式，但只完全到一次多项式，二次项不完全。这不完全的二次项有时可能改善精度，有时则不能。例如，在分析图5-6的“纯弯曲”应力场时，图（a）中的单元将比图（b）中的单元效果好，尽管还不能说满意。提高单元精度的一个途径是增加

10、结点个数，提高插值函数阶次。图5-6(b)（a）5、四结点单元的应用实例及相关限制条件 某求解域如图5-7(a)所示，若将该区域用个四结点等参元进行离散，母体单元如图5-7(b)所示。xy0.01124342332.03.05.05.03.02.00.0142y1432x(c)(a)12341(d)(b)图5-7从图中可以看出：、2号单元与母体单元的结点编号顺序一致，均为逆钟向，而号单元的编号顺序为顺钟向；、号单元为凸形单元，即连接任意两点结点的线段均在单元内部，而单元是非凸形单元，如连接结点、的线段不在单元内。下面讨论这些差别在母体单元与实际单元进行映射时的影响。在母体单元中形函数如式（5-

11、1-1）,坐标变换关系如式（5-1-3）。首先，计算出Jacobi矩阵中的各元素如下下面计算出各单元具体的变换关系及Jacobi行列式的值单元：各结点的坐标为，则Jacobi行列式是的线性函数，对所有的值（）Jacobi行列式的值恒为正，因此，母体单元与单元的变换是可逆的。单元：各结点的坐标为，则 Jacobi行列式的值沿着直线为零，母体单元中的阴影部分将映射到实际单元的阴影部分，这部分显然在实际单元之外。例如，母体单元中的点落在阴影部分，该点映射到了实际单元的。因此，母体单元与单元的变换不是可逆的。所以内角大于网格在任何单元中都是不允许的。一般来说，有限元网格中内角过大或过小都是不合适的。单

12、元：各结点的坐标为。则Jacobi行列式的值小于零表示：右手坐标系映射到左手坐标系，这种变换关系在有限元方法中也是不允许的。若将单元3的结点编号顺序改为逆钟向，即：等参变换成为可逆变换。5-2八结点四边形（直边或曲边）等参数单元8342图5-765、 母体单元形函数母体单元仍为边长为的正方形（图5-8）。自然坐标系，如图所示。在单元中配置八个结点，其中仍位于角点上，则位于各边中点。构造出八个形函数N1N8。如下：(5-2-1)验证可知，具备以下性质(5-2-2)构造出（5-2-1）的八个形函数是不太容易的，下一节将介绍构造这类形函数的一般方法。y76543218x0图5-、 实际单元和坐标变换

13、实际单元e的八个结点坐标为（xi, yi）（图5-9），则母体单元到实际单元e的坐标变换取为(5-2-3)Jacobi矩阵当Jacobi行列式 detJ0时，（5-2-3）规定的变换是非奇异的。（5-2-1）所定义的形函数对于变量或变量来说次数都不超过。沿母体单元中常数的直线（坐标线），根据（5-2-3），x、y将是的二次函数，因而对应于实际单元e中的一族曲线（图5-3）。母体单元中常数的直线将对应于实际单元e中的另一族曲线。在一般情况下单元e将是曲边四边形。当实际单元e为矩形，且结点位于各边中点时，变换（5-2-3）的右端退化为、的一次多项式，反过来、也可表示为x、y的线性函数。、 单元内假

14、设的位移场单元内的位移场（即试探函数）采用与坐标变换相同的一套形函数（5-2-1），对平面问题则认为单元e内有(5-2-4)在（5-2-4）中，u、v是、的多项式，共有项。如果注意到（5-2-1）定义的形函数中、的三次项只有2和2（没有3和3），展开后必归为以下形式的位移场。(5-2-4)完全到、的二次多项式，三次项不完全。一般情况下u、v不是x、y的多项式，但当、可以表示成x、y的线性函数（单元e为矩形，且结点位于各边中点）时u、v将是x、y的多项式，且完全到x、y的二次多项式。、 收敛性分析当detJ0时可保证单元内位移场连续，（5-2-2）则保证了刚体位移和常应变要求。只剩下一个协调性问

15、题。0图5-10y, v76543218x，usMe(a)(b)4现以结点1-5-2所在的边为例（图5-10）。沿此边常数，根据（5-2-5），u、v将是的二次函数。设边上一点与单元边上的点对应，到结点的弧长为，则的坐标将是的单值函数(S)，这说明，沿1-5-2边，u、v是(S) 的二次函数，完全由这条边上三个结点：、处的位移u、u、u、v、v、v决定。从而保证穿过这段边界时位移u、v的连续性。这样八结点等参数单元满足了收敛性所要求的各项条件。、 精度分析由于性质（5-2-2），八结点等参数单元（以及四结点等参元）能够精确地再现任何线性分布的位移场，因而它们的有限元空间Sh次数k1至少为。当实

16、际单元e为矩形，且边结点位于所在边中点时，假定的位移场将包含x、y的完全二次多项式，这种情况下有限元空间Sh的次数k1为。位移的误差为(h3) ，应力的误差为(h2)。图5-118122227654321八结点等参元的边界可以为曲边，因而可以更好地逼近区域的曲线边界。从逼近曲线边界来讲对提高精度有利，但单元形状的畸变（单元为正方形，单元e为曲边四边形！）却有可能损失插值精度。只有在单元畸变不大的情况下，才能保证位移和应力的误差分别为(h3)和(h2)。所谓畸变小是指：(i)连接直边四个角点所成直边四边形的内角远大于0o，远小于180o。（ii）上述直边四边形接近平行四边形或矩形。即图5-11中

17、的 1=(h2)(iii)实际单元（直边或曲边）的边结点与上述直边四边形中点的偏离较小。即图5-10中的2=(h2)对于工程问题中通常采用的、不是有意加密的网格而言，上述条件很难同时得满足。在这种情况下达不到位移(h3) 和应力(h2)的精度，但一般总比四结点单元精度好。图5-12中，显然 (a) 形状比 (b) 形状的单元畸变小，因而精度也比 (b) 好。图5-1268(a)(b)342178754326515-3四九结点等参数元通用分析程序中选用的往往不是固定的四结点或八结点单元，而结点数目可由用户在一定范围（例如四九）内任意选择。这样不仅为用户提供了不同精度的单元，而且为用户同时使用不同

18、精度的单元时提供了过渡方法（图5-13）图5-13四结点单元（低精度单元）五结点单元（过渡单元）八结点单元（高精度单元）仍取边长为的正方形为母体单元，自然坐标系、取法同前面一样。在母体单元中安排了九个位置，用户可以根据自己的需要在这些位置上安排结点（图5-14）。其中位于角点，这四个节是必需的。8位于各边中点，位于形心，后面这五个位置可以自由安排结点。从上一节介绍的两种等参数单元可知，在等参数单元中形函数起着关键的作用。下面重点讨论不同情况下构造形函数的方法。4图5-1494xyxy1,43413ee图5-15(b)(c)(a)、 四结点单元即第5-1节中的单元，形函数如（5-1-1）所示。母

19、体单元如图5-15（a），一般实际单元如图(b)。若四边形单元的两点结点（例如和）重合为一点，将得到图 (c) 的三角形单元在重合的结点上detJ0，但是进一步对位移场分析可知，这时得到的正是常应变三角元，只要在积分时避开重合结点，就不会遇到什么困难。4图5-16、 五结点单元设在位置处再增加一个结点，则成为五结点单元(图5-16)。(1) 构造形函数为了保证在（0, -1）处为，在其他三条边上为零，可以构造出(5-3-1)(2) 修改结点、当前的形函数四结点单元形函数N1、N2在(5、5)的值为为此必须对它们进行修改(5-3-)而N3、N4则维持原状(5-3-)（5-3-1）（5-3-3）就

20、是图(5-16)五结点单元的形函数。当第五个结点取在位置时可用类似的步骤得出形函数。4图5-17、 六结点单元设第六个结点取在位置(图5-17)。（） 构造形函数(5-3-4)（） 修改结点、当前的形函数由（5-3-2）和（5-3-3）可知在（-1,0）之值为均将它们修改为(5-3-5)N2、N3不必修改(5-3-6)注意，这时N1、N5、N4已同（5-2-1）相应的形函数相同。再加一个边结点仍然要做类似的两个事：补充定义一个边点形函数，修改两个角点当前的形函数。当位置均设置结点时，就成为第5-2节的八结点单元，形函数即（5-2-1）、 九结点单元4图5-18设第九个结点取在形心处(图5-18

21、)。（） 构造形函数 N9(,)(5-3-7)要求N9(,)在四边上为零，在形心(0,0)处为。可取（）修改结点当前的形函数，由（5-2-1）有将N1(,) 修改为(5-3-8)类似(5-3-9)修改N5(,),(5-3-10)类似(5-3-10)上述以四结点单元为基础，随着结点增加，逐步扩充、修改形函数的方法很容易由计算机去执行。如果有必要，将结点个数扩充到16也是容易实现的。一旦形成了形函数，单元分析的工作就可以按统一的步骤进行。5-4等参数单元的单元分析公式本节以八结点单元为例，讨论等参数单元的单元分析公式。将单元内位移场表示为矩阵形式其中坐标变换为、 Jacobi矩阵J及其逆阵J-1(

22、5-4-1)(5-4-2)若detJ0则J-1存在上式中的各元素的分子、分母一般情况下均为,的多项式，当detJ为常数时（实际单元为矩形时），上式中的各元素退化为多项式。若实际单元近似为矩形，则detJ近似为常数时，各元素为近似的多项式（与插值多项式阶数一致的多项式）、 求导公式合并成得到（下式必须以解析的形式给出各元素的表达，以便后续的积分能顺利进行）、 应变几何矩阵其中B称为几何矩阵，由八个子块组成(5-4-4)几何矩阵是从单元刚度矩阵中分离出的一个任务，在具体计算中有重要的用途。、 单元刚度矩阵单元变形能单元刚度矩阵(5-4-5)、 单元载荷的等效结点力（） 体积力（fx, fy）等效结

23、点力(5-4-6)例当实际单元也为正方形时（与母体单元一样），若fx为均布力，则将其等效到个结点上的情况为图5-19所示。, y44，xfx图5-19由虚功可知具体计算时，对本例中的正方形单元而言，只需计算两类结点上的等效结点力（角点和边中点）第、结点上的等效结点力与第结点上的等效结点力相等。第、结点上的等效结点力与第结点上的等效结点力相等。qlmnyx0(1,-1)(1,1)(-1, 1)(-1,-1)图5-20（）边界力的等效结点力(i) 当实际单元e的形状十分规则时，分析可以直接在总体坐标下进行。例设实际单元e为边长为的正方形，在结点l、m、n所在边上作用均布压力q, x, y轴的方向如图2-20所示，单元厚度为t。求q的等效结点力解沿lmn边位移v是x的二次函数，三个形函数为三个量的比值为1：，既不是：也不是：。(ii)当单元形状比较复杂时，必须借助于母体单元。例在单元e的一边lmn上作用均布压应力p，与母体单元上的边对应，单元厚度为t。求p的等效结点力。图5-21nx,upmley,v25(a)(b)解沿边lmn=1弧微分边界内法线的方向余弦为，等效结点力（5-4-7）（）温升载荷的等效结点力初应力虚位移（虚应变）由虚功等效的原则有其中（5-4-8）注意上述的积分