第五章 其他常用单元的刚度矩阵

1、第五章 其他常用单元的刚度矩阵除了前面讲的一维、二维杆单元及三角形单元之外，有限元法中还根据分析对象的不同采用许多其他单元，如三棱圆环单元、等参数单元、平面四边形单元、四面体单元、六面体单元等等。鉴于学时所限，只介绍三棱圆环单元和等参数单元的刚度矩阵的求法，对其他单元同学们可查阅有关书籍。第一节 三棱圆环单元的刚度矩阵机器中许多零件如飞轮、缸体等在几何形状上具有共同点，即它们都是某一平面图形绕平面内某一轴线旋转而形成的回转体，此平面称为子午面。当回转体承受的载荷和支撑条件相对于该轴线也对称时，分析求解这类零件的应力、应变问题，称为轴对称问题。轴对称问题中，回转体内各点只有轴向和径向两个方向的位

2、移，一个三维问题就简化为二维问题。对这类零件的离散化可以在子午面内进行，最常用的是三角形截面的轴对称单元，简称为三棱圆环单元。如图 4-1 所示。1. 位移模式及形状函数由于轴对称的特点，不再用直角坐标系（,而用柱面坐标系（ ,z）描述物体。物体内任意一点只有沿 r 和 z方向的位移 u 和，而无 方向的位移。当纵剖面上三角形单元(e)的三个节点总码分别为 I、j、k 时，如图 4-1 所示，1相应的节点位移向量为T ( ) euwuwuwiijjkk与弹性力学平面问题中的三角形单元一样，采用线性位移模式，则u(r, z) w(r, z) r z123 rz456与平面问题的推导步骤完全相同，

3、可以得到与平面问题相似的结果：(r, z) ( )u e(e) ( , ) (e) (e)( , )r zN r zw(r, z)其中形状函数为：1N (r, z) (r z r z ) z r r z12jkkjjkkj1N (r, z) (r z r z ) z r r z22kiikkiik1N (r, z) (r z r z ) z r r z32ijjiijji2.应变与位移的关系（几何矩阵）轴对称问题中表示应变与位移关系的几何方程与弹性力学平面问题相似，所不同的是：单元内一点在径向产生的位移u r的圆 2 r u 后，其圆周长度变成。因此，在圆周方向的应变为2 (r u)2 (r

4、u) 2r2rur 由于各点在圆周方向上无位移，因而剪应变 和 均为vvrr2urur零。将应变写成向量的形式，则根据上式，可推导出几何方程 zrw zu wrz r z B(e)z00z0z0jkkjijN (r, z)N (r, z)( , )N r z00 jikrr0r0其中几何矩阵1 B 0rrikrkjjiij2 rzrzrzkjjkikkiji3.弹性方程和弹性矩阵D依照广义虎克定律，同样可以写出在轴对称中应力和应变之间的弹性方程，其形式为1 ( )urErz1 ( )uErz1 ( )uzEzr2(1 )r rzErz所以弹性方程为 D式中应力矩阵 Trzrz1001 弹性矩阵

5、 ED 1 0 )(1 2)1 200024.单元刚度矩阵 k(e)3与平面问题相同，仍用虚功原理来建立单元刚度矩阵，其积分式为 TkBD B (e)V在柱面坐标系中， 2 将代入，则 T(e)kBD B 2 k(e) BD B 2TV即为轴对称问题求单元刚度矩阵的积分式。与弹性力学平面问题的三角形单元不同，在轴对称问题B内有的元素（如等）是坐标rz的函N (r,z)ir数，不是常量。因此，乘积不能简单地从式 TBD B的积分号中提出。如果对该乘积逐项求 2 B D B T(e)k积分，将是一个繁重的工作。一般采用近似的方法：用三角 形形心的坐标值代替几何矩阵B内的 r 和 z 表示B在形心

6、处计算出的矩阵B。其中(r,z)(r r r )(z z z )r k , z ijijk33只要单元尺寸不太大，经过这样处理引起的误差也不大。被积函数又成为常数，可以提出到积分号外面： T式中 三角形的面积。 k 2 BD B 2 BD B rT(e) T由式T k 2 B D B 2 BD B r(e)称的三角形单元，当形状、大小及方位完全相同而位置不同时，其刚度矩阵也不相同。距离主轴线越远的单元，其刚度4越大。这与平面问题不一样。二、等参数的刚度矩阵对一些由曲线轮廓的复杂结构，如果采用直角边单元进行离散，由于用直线代替了曲线，除非网格划分得很细，否则不能获得较高的精度；对另一些应力随坐标

7、急剧变化的结构，采用简单的常应力单元离散时，也必须划分成大量的微小单元，以保证足够的精度。为此引入一种高精度的单元等参数单元。它既能简化复杂单元划分的工作，又能在满足同样精度的要求时，大大减少使用的单元数。目前流行的大程序中较常用，它成功地解决了许多二维和三维的弹性力学问题。为导出等参数单元的刚度矩阵，首先要建立根据每个单元的形状确定的自然坐标系，然后将位移模式和形状函数都写成自然坐标的函数。一个单元在自然坐标系内的点余元整体坐标系内的点成一一对应的关系。通过映射，可以将整体坐标系中的图形转化为自然坐标系中的相应徒刑。例如可以将整体坐标系中的一个任意四边形（实际单元）映射到自然坐标系中成为一个

8、正方形（基本单元）。同样也可以将任意四面体、六面体（包括直边和曲边的）分别映射成正四面体和正六面体。这里只介绍较简单的一种平面问题的情况，将整体坐标系中的一个任意四边形映射成自然坐标系中的一个正方体，5不再叙述。1. 位移模式和形状函数图4-2 取其交点为原点，这两条直线分别为 和 的 和 值分别为 1自然坐标系。原坐标系 XOY 称为整体坐标系。在整体坐标系中，自然坐标系非正交，它由任意四边形的形状所确定。图 4-19如果将自然坐标系改画成直角坐标系，那么图 4-19（）中的任意四边形单元就成为图 4-19（）所示的正方形。上述两个四边形的点（包括顶点）一一对应，即它们之 Y和自然坐标、之间

9、的坐标转换式，即*X Y 1234 5678四边形四 个顶点 的坐标值在 XOY 坐标系中 分别为： 在 坐 标 系 中 相 应 为X ,Y , X ,Y , X ,Y , X ,Yo11223344。将有关数据代入*中的第一式，则有 1 , 1 , , X ,XX 1123421234X , 3123441234求解上述方程组得：6X X X X X X X X , 1234123441424 X X X XX X X X , 12341233444坐标变换方程成为11 XXX X 1 1 1 X41234同理11 Y Y 1 Y 1 Y 1 Y41234当引入函数后，坐标变换方程成为 ,N

10、i4 ,N XX Y iii14 ,N Yiii1式中1 , 1 1 N4iii变量 的正负号由相应节点的坐标值决定。例如 、 、ii 当i=4 时，因此，。1 1 , 1N4444下面再来研究函数的特性。 ,Ni对节点，相应的自然坐标值为（ -1，-11 X ,Y11中 很 容 易 看 出 ， 除 N =1 外 ，1 , 1 1 1N4iiiN =N =N =0。对其余 各节点也一样 。总而言之， 对节点234i(i=1,2,3,4)，除N=1 外，其余三个 N 值均为零。同时，不难i看出，即四个节点的 N 函数 , , , N, 1NNNi1234之和等于 1。函数具备上章所介绍的形状函数

11、应满足的条件， ,Ni7可作为本单元的形状函数。采用做形状函数，其位移模式为 ,NNi4 4 ,u , u vNviiiii1i14 X Y N,X对比和可以看出：ii44 N , u ,v N , vi1u4 iiii,NYi1i1iii1在这种实际单元（任意四边形）中，坐标变换式和位移模式不仅采用了相同的形状函数 ,Ni模型。这种性质的实际单元称为等参数单元。对用节点位移值 u（或 v 等）求单元内某一点位移量 u(或 vii等)的插值公式，只要将 u(或 v 等)换成 X(或 Y4 ,uu Niii1等) X(或Y 等)求相应点坐标 X(或 Y等)ii的插值公式。相反也是这样。2.几何矩

12、阵B由于几何矩阵B通过对位移求偏导数而得出，所以首先必须利用复合函数求导的规则得出下述公式 u u Xu Y u u XYXu或写成 Juu XX u Yu Y yXY 式中，此式称为雅可比矩阵。 J XY 为了将几何矩阵B写成变量的函数，必须将式 、8 u u 改写成 Xu Ju y u v u v ，同理， XuXv J 1 J 1uv y y从表示单元内各点位移与其应变关系的几何方程可知： 0X1 0 0 0 uu uvv T 0 0 0 0 1 Y Y XYv0 1 1 0YX u v u v 将式和 合并，则 XuXv J 1 J 1uv y y u uuvv T10uv TJv 1

13、XY XY0 J对单元（e u，v 对自然坐标 的偏 ，导数可利用上式求出，写成矩阵形式为： u uv Tv N e p式中 T uevuvuvuv11111111 0000NNNN N 1 2 p 3 p 4 pp0N0N0N0Np4 p1p2 p3 p对于i=1,2,3,4 NN TN ii 将u和 uuY Xvv T10uvvTJ 1XY0 J9 代入uuv Tv e N p 0X1 0 0 0T，则可得出表示 uuuYvv 0 0 0 0 1 YXYv0 1 1 0YX在整体坐标系中位移和应变关系的几何方程： eeBe式中的几何矩阵B是自然坐标 的函数：，1 0 0 0 N0 J 1B

14、 0 0 0 1 0J1p0 1 1 04 X Y N,X也可利用 求得的 以及和NNTiiN Ni1iiip4 ,NYiii1XY 求出 ， 。 XXXYX TJ JJ N1234NNNXY1 p2 p3 p4 pYYY1234 3.单元刚度矩阵 ek设单元板厚为 t，根据虚功方程有：， ke B T D B tdAA此式中几何矩阵B和弹性矩阵D都已求出。因为几何矩阵B中的变量是自然坐标面积 。，在实际单元中任取一点 ，其整体坐标位 ，其相应的自然坐标为 。过 p 点做 的等值线，同时做，的等值线，围成一小块微分面积 dA，如图 4-20d， d（ pqrs 4-20（）10所示。实际上，取

15、得很小，因此该四边形可视为平行， dd四边形。若相邻的两边用向量表示，则两者的乘积等于 a b该平行四边形的面积 dA。图 4-20dA a b sin a b若, a a i a j b b i b j，则xyxyaabbdA a i a j b i b j xyxxyxyy为了求出的值，要先写出 a和 b 两端节点 、sa ，a ,b ，bxyxy的坐标值。点p： X, ,YYY,Xpp点q：点：X X , ,d ,Y dqqX X, ,dY, Y dss利用泰勒技术展开并略去高阶项，可得X , , XdXXdX , , Xdd对，也可写出相应的展开式。利用式 , , , YdYdXX ,

16、 , XdXd可得： , , XdXdXYa X X ,d aYYdX ，xqpyqpYb X X ,d bYYdxspysp11将此式代入式 得到：aabbdA a i a j b i b j xyxyxyxyXY，简写为aabbdA x d ddA J ddxyXYy 单元刚度矩阵为：，这个积分可以采用“数值方法”，用 11 B T D B t J d d e k11高斯求积分公式很方便的求出，在此不作介绍。例：求如图所示四边形的雅可比矩阵。解：求雅可比矩阵可在整体坐标系中进行，也可以在实际单元的局部坐标系中进行。为便于计算，本例在局部坐标系中进行。对单元（1将四个节点的自然坐标值（ -1，-11，-11，1（-1，1）代入下式：，再将所得到的值及四个节点1 , 1 1 NN4iiii实际单元在局部坐标系的坐标值（-3-23-232（-3，2）代