第 7 章平面问题高阶单元

第7章平面问题高阶单元 7.1 位移模式阶次的选择在前面两章中讨论了平面问题三结点三角形单元，其位移模式的最高阶是坐标x、y的一次项。这种位移模式导致单元常应变、常应力特性，单元应变矩阵、应力矩阵、刚度矩阵均为常数矩阵，因此计算非常简单。但这种单元难以反映应力梯度的迅速变化。要想提高计算精度，必须细分网格，增加单元数和点数，因而加大输入数据的工作量。提高计算精度的另一条有效途径是采用高阶单元。由于高阶单元的应变、应力不再是常数，因此采用少量单元就

可能达到较高的精度。图7-1悬臂梁分别采用高、低阶单元计算就是一个典型的例子。可能达到较高的精度。图7-1悬臂梁分别采用高、低阶单元计算就是一个典型的例子。 P A B h 4h • • • • • • 常应变单元：A=0.866 B=0.619 悬臂深梁解析解：A=1.0 B=1.0 高阶单元：A=0.99 B=0.99 图7-1

1 x y x2 xy y2 x3 x2y xy2 y3 x4 x3y x2y2 xy3 y4 图7-2 多项式选择的怕斯卡三角形选择位移模式时，第2章提到要考虑解的收敛性，即要考虑到位移模式的完备性和协调性。实际操作中，一般应考虑位移模式的对称性。这是因为，有限元位移模式的选择实际是以帕斯卡（Pascal）三解形基础上的（如图7-2所示），由低价至高阶，顺序选取，组成多项式。多项式中的项数等于单元节点自由度数。如三节点三角形单元，位移模式取完全一次式，共3项。六节点三角形单元，位移模式取完全二次式共6项。如果某一阶次不能全取，则应按对称性原则适当选取。

例如在下节将要讨论的四 结点矩形单元中，位移模式不能取1，x，y，x2四项，也不能取1，x，y， y2四项，而应取1，x，y，xy四项。 1 x y x2 xy y2 x3 x2y xy2 y3 x4 x3y x2y2 xy3 y4 图7-2 多项式选择的怕斯卡三角形 7.2 四节点矩形单元图7-3示出的矩形单元，边长分别为2a和2b。取4个角点为节点，编号为i，j，l，m。将x轴和y轴置于单元的对称轴上。 1、位移函数单元的位移函数可取为：

y m l 图7-3 b x b j a a i （7-1）在上式表示的位移模式中，a1, a2, a3, a5, a6, a7， a8反映了单元的刚体位移和常应变。在单元的边界(x=±a 或y =±a)上（或），位移是按线性分布的。因此，相邻单元在公共边上的位移是连续的。这样，位移模式满足了解答收敛性的充分条件。在式（7-1）中代入节点位移和节点坐标后，可解出

（7-3） （7-2）各待定系数(a1… a8)。将这些系数再代入式（7-1），可得：式中形函数为：

（7-4） （7-5）（7-6）令在节点上的值为：则式（7-3）可简写为将位移函数写成矩阵形式，即有与式（2-20）相同的形式式中

（7-8） （7-9）（7-7）（7-10）其中，I为二阶单位矩阵。 2、应变矩阵根据几何方程，可得与式（2-25）同样的形式把应变矩阵[B]写成子矩阵形式其中

由此可见，[B]是、 的函数，即是x、y的函数。因此单元中的应变不再是常数。（7-11）（7-12）（7-13） 3、应力矩阵根据应力-应变关系，可以计算单元中的应力，得到式（2-28）同样形式应力矩阵[S]具有与式（2-29）同样形式将[S]写成子矩阵形式

其中 （2-33a）（7-14）上式对应平面应力情形。对于平面应变情形，只需将其中的E，作相应的改变即可。 4、单元刚度矩阵单元刚度矩阵可采用式（2-33a）进行计算

（7-16） 在四节点矩形单元中，[k]是一个8×8的矩阵。将[k]写成分块形式：其中的子矩阵[krs]2×2可由下式计算

（7-17） 上式对应平面应力情形。对于平面应变情形，只须将上式中的E、作相应的改变。 5、等价节点力单元体积力和表面力引起的节点力仍可用式（2-45）和（2-46）进行计算。

（7-18） （2-45）（2-46）对本问题给定的位移函数，若体积力是重力的情形（设重度为），单元等价节点载荷列阵为：有了对单元的上述结果，便可应用第5章的方法组集结构刚度矩阵和节点荷载向量；求解节点位移；计算内力和应力。四节点矩形单元采用较高阶的位移模式，具有比三节点三角形单元较高的计算精度。但矩形单元也有缺点，

y j i • m • m • j i x 图7-4 一是不能适应斜线及曲线边界，二是不便于采用大小不同的单元。 7.3 六节点三角形单元 1、位移模式在三角形单元i, j, m的各边中点增设一个节点，使每个单元具有6个节点，得到图7-4所示的六节点三角形单元。这种单元具有12个自由度，可以采用完全二次多项式的位移模式：

（7-20） 所取位移模式反映了单元的刚体位移和常应变；单元内部是连续的；在单元边界上位移分量按抛物线变化，而每条公共边界上有3个公共结点，可以保证相邻两单元位移的连续性。因此，上述位移模式满足收敛的必要和充分条件。上述位移模式确定之后，可以用分析三节点三角形单元和四节点矩形单元相同的方法进行分析。得到形函数、应变矩阵、应力矩阵、单元刚度矩阵、等价节点力向量。但其过程十分繁复，采用面积坐标可以大大简化计算。

2、面积坐标 y （7-21） j · P m i x （1）定义对于一个三角形ijm（图7-5），三角形内任一点P（x,y）的位置，可以用如下的三个比值来确定： Ai Am Aj 图7-5

其中A为三角形ijm的面积，Ai, Aj, Am分别为三角形的Pjm, Pmi, Pijd的面积。这三个比值Li, Lj, Lm称为P点的面积坐标。（7-22）由于则由此可见，P点的三个面积坐标不是独立的。同时，面积坐标只是用以确定三角形内部某点的位置，因而是一种局部坐标。下面进一步给出面积坐标的几个性质。（2）面积坐标与直角坐标的关系在图7-5中，三角形Pjm的面积为

（7-23） （7-23a）（7-24）由式（7-23），式（7-21）化为将式（7-24）、（7-23a）和式（2-18）、（2-17）对比，可知，面积坐标就是三节点三角形单元的形函数 Ni、Nj、Nm。

将式（7-24）的3个式子分别乘以xi, xj, xm，然后相加，并利用关系式（7-23a），有同理（7-25）（3）面积坐标的导数公式根据面积坐标与直角坐标的关系，由复合函数的求导公式，有

（7-26） （4）面积坐标的积分公式下面给出面积坐标的幂函数积分公式。它们在计算单元刚度矩阵和等效结点载荷时有用。

（7-27） （7-28）在三角形单元上进行积分时，有在三角形某一边（设ij边，边长为l）上进行积分时，有 3、用面积坐标表示六节点三角形单元计算公式（1）形函数和位移表达式对应如图7-4所示的六节点三角形单元，形函数可用面积坐标表示为

y j i • m • m • j i x 图7-6 （7-29）现利用形函数的性质检验式（7-29）的正确性。先考虑三角形的角点，例如图7-6中的i点，有由式（7-21）（P16），有代入式（7-29），有

y j i • m • m • j i x 图7-7 说明形函数Ni在i点等于1，在其它节点等于0，因此是正确的。再考虑三角形的边中点，例如i点，面积划分如图7-7所示。显然有：由式（7-21）(P16)，有代入式（7-29） (P16) ，进一步说明式（7-29）所表示的形函数的正确性。

形函数确定后，单元中任意一点的位移可以表示为：形函数确定后，单元中任意一点的位移可以表示为：（7-30）（7-32）（7-33）其中（7-31）其中I为二阶单位阵，形函数由式（7-29）确定。（2）应变矩阵单元中的应变仍可表示为：

式中应变矩阵[B]为： （7-34）（7-35）其中

（3）应力矩阵 （7-36）（7-37）单元中的应力仍可表示为：式中[D]是弹性矩阵，由式（2-9）确定；应变矩阵由式（7-34）、（7-35）确定。根据矩阵乘法，可以给出用面积坐标表示的应力矩阵[S] （4）单元刚度矩阵单元刚度矩阵仍可表示为：根据[B]、[D]的表达式以及面积坐标的积分公式（7-27），可以求出[k]中元素的显式表示。由于较为繁复，这里就不列出详细结果。

（5）等价节点力向量 （7-38）由于位移模式是非线性的，因此体积力和表面力引起的节点力向量不能采用静力等效原理进行分配，而应采用相应公式进行计算。单元体积力引起的等价节点力计算公式仍为：将由式（7-29）、（7-32）表示的[N]代入，并应用积分式（7-27），可以计算FVe。例如对于重力引起的FVe，有它表示各边中点承担单元重力的1/3。

单元表面力引起的结点力计算公式仍为： （7-39） y y j j ps pslh/6 i i • • 2pslh/6 m m • • 4pslh/6 m m l • • l pslh/6 j j ps pslh/6 i i x x 图7-9 图7-8 设在ij边上受有x方向的均匀分布力ps，对应的等价节点力向量为（图7-8）

如在ij边上受到x方向的三角形分布面力，其集度在i点为ps，在j点为0。对应的节点力向量为（图3-9） 它表示边中点承担载荷的2/3，载荷集度大的角节点承担1/3。六结点三角形单元中的应变、应力不为常量，因此可以应用于应力梯度较大的地方，精度较高。显然，其计算也较复杂。 7.4 四节点四边形等参数单元 1、等参数单元的概念

现在，我们从任意四边形单元着手，介绍等参数单元的概念。现在，我们从任意四边形单元着手，介绍等参数单元的概念。 y 3 4 2 1 x 图7-10 任意四边形单元前面讲到的四节点矩形单元虽然比较简单，但难以应用于斜线边界。图7-10所示四节点任意四边形单元容易适应这种边界，但要在整体坐标系内，写出它的统一的形函数又是相当复杂和困难的。

但是若能找到它与一个规则正方形的关系，就能写出它的统一的位移模式，这可以通过坐标变换来解决。在图7-10所示四边形单元上，用等分四边的两族直线分割该四边形，以两族曲线的中心（ =0、=0）为原点，沿、 增大的方向作轴和轴，并令四边的=±1、 =±1，就得出一组新坐标系 （图7-11）。 y x 3 = 1 = 1 4   = -1 2 = -1 1 这里，、是一种局部（单元）坐标，它只应用于单元范围内。而x，y是整体（结构）坐标，它适用于所有的单元。图中的任意四边形单元是研究对象，称为实际单元。图11 实际单元

为了得出实际单元的位移模式和局部坐标与整体坐标之间的变换关系，引入一个四节点的正方形单元，称基本单元（图7-12） 。 3 4 2 1 图7-12 基本单元 = 1  = 1 = -1  = -1 （7-40）参照式（7-2）和（7-3）P6，此基本单元位移函数可写为：

其中，形函数应为： 引入新变量 i、i（i=1,2,3,4) 基本单元的形函数被写成：

（7-41） （7-40）现在，把基本单元的位移模式（7-40）和形函数式（7-41）移用于图（7-11）所示的实际单元，则实际单元的位移模式取为：在结点处：在其它结点处：

利用形函数的上述性质，可以将任意四边形的整体坐标写成：利用形函数的上述性质，可以将任意四边形的整体坐标写成：（7-42）且，式中的形函数Ni仍由式（7-41）确定。而把式（7-41）中的、理解为图7-11所示实际单元的局部坐标，i、i便是实际单元中节点i的局部坐标。任意四边形单元中结点的整体坐标，如果它已知，那么(7-42)表示了局部坐标与整体坐标的变换

实际单元是任意四边形四节点单元，基本单元是正方形单元，可以认为：实际单元是对基本单元通过变换得来的。由于实际单元的位移模式中采用了基本单元等同的形函数，这个实际单元就称为等参数单元。实际单元是任意四边形四节点单元，基本单元是正方形单元，可以认为：实际单元是对基本单元通过变换得来的。由于实际单元的位移模式中采用了基本单元等同的形函数，这个实际单元就称为等参数单元。另一方面，式（7-42）表明了实际单元中局部坐标（、）与整体坐标（x、y）的一一对应关系，是一个坐标变换式。类似于本章3.2节进行的四结点矩形单元的特性分析，可以建立等参单元的应变矩阵、应力矩阵、刚度矩阵、节点力向量等的计算公式。与前面不同之处在于，在等参数单元法中，要将对整体坐标x、y的导数计算和积分计算转换为对局部坐标、的导数计算和积分计算。

例：实际单元的结点整体坐标如图(a)中括号内数字所示，基本单元的结点局部坐标如图(b)中括号内数字所示。例：实际单元的结点整体坐标如图(a)中括号内数字所示，基本单元的结点局部坐标如图(b)中括号内数字所示。 y 4(－1,1) 3(1,1) 3(1,2) 4(0,1) 0 0(3/4,3/4) x 1(0,0) 2(2,0) 1(－1,－1) 2(1,－1) 图(a) 实际单元图(b) 基本单元 (1) 试验证基本单元上的结点局部坐标与实际单元上对应点的整体坐标的对应关系。 (2) 求基本单元的局部坐标原点( ),在实际单元上的整体坐标(x,y)是多少？

解（1）以3结点为例，根据形函数的性质：在3结点处：应用式（7-42）说明了由基本单元上结点的局部坐标可映射出实际单元上对应的结点整体坐标（2）由于，由（7-40）P33 得：

代入（7-42） 由上例可知：利用（7-42）在基本单元上任意一点，都可以在实际单元上找到一个对应点的坐标(x,y), 这样就把实际单元与基本单元紧密地联系起来。反之，则比较困难，这是因为形函数是一个二次函数。为了避开这个困难，一般都假定基本单元上已知点去求实际单元上的对应点。

2、应变矩阵 （7-8）（7-43）（7-9）（7-44）单元的几何方程与式（7-8）、（7-9）相同，即：式中 ?

（7-45） 这里采用记号由于形函数式（7-41）是用局部坐标、给出的，将、看作x、y的函数，根据复合函数的求导规则，有：上式可记为：

上式右边第一个矩阵称为雅可比（Jacobi）矩阵： （7-48） (7-46) 其逆矩阵为： (7-47) 式中|J|为雅可比行列式由式（7-42）p35，有

（7-49） （7-50）（7-51）由式（7-41）p34，有由式（7-45）p41，有

和（7-52）式中分别由式（7-47）和（7-50）确定。从而由式（3-43）、（3-44）确定出应变矩阵[B]。 3、应力矩阵应力矩阵仍由下式得到 4、单元刚度矩阵单元刚度矩阵是一个8×8的矩阵，仍为

由于[B]是用局部坐标系、给出的，坐标变换时有面积微元公式 （7-53）因此，[k]可由下式计算 5、等价节点力向量（1）体积力设单元的体积力是，则等价节点力公式为

（7-54） （7-55） (2) 表面力设单元的某边（如对应的=±1）上作用有表面力ps=[psx psy]T，则该边上节点的节点力为

上式中，用到坐标变换时线积分微元公式(当=常量时) p37 对于=±1的边界表面力，等价节点力的计算过程完全一样。 6、关于高斯（Gauss）积分在计算单元刚度矩阵式（7-53）和等价节点载荷向量式（7-54）、（7-55）时，由于被积函数比较复杂，通常可采用数值积分。即在单元上选择某些点，称为

积分点，求出被积函数在这些积分点上的数值，再用一些权函数乘这些函数值后求和，就可得到近似积分值。高斯积分法是数值积分法中具有较高精度的方法。现简要介绍如下：积分点，求出被积函数在这些积分点上的数值，再用一些权函数乘这些函数值后求和，就可得到近似积分值。高斯积分法是数值积分法中具有较高精度的方法。现简要介绍如下：（7-56）（1）一维高斯积分公式式中Hi对应于积分点i的权函数，对于n=2~4个积分点，其坐标i和权函数Hi列于表7-1中。

（7-57） 表3-1 高斯求积公式中的积分点坐标和权系数（2）二维高斯积分公式

积分点j、i和权函数Hi、Hj同样可按表3-1进行取值。积分点j、i和权函数Hi、Hj同样可按表3-1进行取值。（7-58）（7-59）（7-60）对式（7-53）进行高斯积分，有对式（7-54）进行高斯积分，有对式（7-55）进行高斯积分，有

第 7 章 平面问题高阶单元