常用单元的刚度矩阵

rur r u r =-+=πππεθ22)(2由于各点在圆周方向上无位移，因而剪应变θr v 和r v θ均为零。

将应变写成向量的形式，则{}⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧∂∂+∂∂∂∂∂∂=⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧=r w z u z w r u r u rz z r γεεεεθ根据上式，可推导出几何方程{}[]{})(e B ϕε=其中几何矩阵[]⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡∆=ij jikiikjkkj ji ik kj k j i ijkjjkz r z r z rr r r r z r N r z r N r z r N z z z B 0000),(0),(0),(00021 3.弹性方程和弹性矩阵[D]依照广义虎克定律，同样可以写出在轴对称中应力和应变之间的弹性方程，其形式为[])(1θσσσε+-=z r r u E [])(1z r u E σσσεθθ+-=[])(1θσσσε+-=r z z u Erz rz Er τμ)1(2+=所以弹性方程为{}[]{}εσD = 式中应力矩阵{}{}T rz z r τσσσσθ=弹性矩阵[]⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡-----+=221000010101)21)(1(μμμμμμμμμμμμED 4.单元刚度矩阵[])(e k与平面问题相同，仍用虚功原理来建立单元刚度矩阵，其积分式为[][][][]dV B D B k VT e ⎰=)(在柱面坐标系中，drdz dV π2=将drdz dV π2=代入[][][][]dV B D B k VT e ⎰=)(，则[][][][]rdrdz B D B k T e ⎰⎰=π2)(即为轴对称问题求单元刚度矩阵的积分式。

与弹性力学平面问题的三角形单元不同，在轴对称问题中，几何矩阵[B]内有的元素（如rz r N i ),(等）是坐标r 、z 的函数，不是常量。

因此，乘积[][][]B D B T 不能简单地从式[][][][]rdrdz B D B k T e ⎰⎰=π2)(的积分号中提出。

如果对该乘积逐项求积分，将是一个繁重的工作。

一般采用近似的方法：用三角形形心的坐标值代替几何矩阵[B]内的r 和z 的值。

用[]B 表示在形心),(z r 处计算出的矩阵[B]。

其中3)(,3)(k j i k j i z z z z r r r r ++=++=只要单元尺寸不太大，经过这样处理引起的误差也不大。

被积函数又成为常数，可以提出到积分号外面：[][][][][][][]∆==⎰⎰r B D B rdrdz B D B k TTe ππ22)(式中∆——三角形的面积。

由式[][][][][][][]∆==⎰⎰r B D B rdrdz B D Bk TTe ππ22)(可以看出，两轴对称的三角形单元，当形状、大小及方位完全相同而位置不同时，其刚度矩阵也不相同。

距离主轴线越远的单元，其刚度越大。

这与平面问题不一样。

二、等参数的刚度矩阵对一些由曲线轮廓的复杂结构，如果采用直角边单元进行离散，由于用直线代替了曲线，除非网格划分得很细，否则不能获得较高的精度；对另一些应力随坐标急剧变化的结构，采用简单的常应力单元离散时，也必须划分成大量的微小单元，以保证足够的精度。

为此引入一种高精度的单元——等参数单元。

它既能简化复杂单元划分的工作，又能在满足同样精度的要求时，大大减少使用的单元数。

目前流行的大程序中较常用，它成功地解决了许多二维和三维的弹性力学问题。

为导出等参数单元的刚度矩阵，首先要建立根据每个单元的形状确定的自然坐标系，然后将位移模式和形状函数都写成自然坐标的函数。

一个单元在自然坐标系内的点余元整体坐标系内的点成一一对应的关系。

通过映射，可以将整体坐标系中的图形转化为自然坐标系中的相应徒刑。

例如可以将整体坐标系中的一个任意四边形（实际单元）映射到自然坐标系中成为一个正方形（基本单元）。

同样也可以将任意四面体、六面体（包括直边和曲边的）分别映射成正四面体和正六面体。

这里只介绍较简单的一种平面问题的情况，将整体坐标系中的一个任意四边形映射成自然坐标系中的一个正方体，并导出单元刚度矩阵。

其它种单元的映射，可依次原理进行。

不再叙述。

1. 位移模式和形状函数图4-2中的任意四边形单元上，作连接对边中点的直线，取其交点为原点，这两条直线分别为ξ和η轴，并令四条边上的ξ和η值分别为1±，建立一个新的坐标系，称之为该单元的自然坐标系。

原坐标系XOY 称为整体坐标系。

在整体坐标系中，自然坐标系非正交，它由任意四边形的形状所确定。

图4-19如果将自然坐标系改画成直角坐标系，那么图4-19（a ）中的任意四边形单元就成为图4-19（b ）所示的正方形。

上述两个四边形的点（包括顶点）一一对应，即它们之间相互映射。

因此，需要写出整体坐标X 、Y 和自然坐标ηξ、之间的坐标转换式，即ξηαηαξααξηαηαξαα87654321+++=+++=Y X *四边形四个顶点的坐标值在XOY 坐标系中分别为()()()()44332211,,,,,,,Y X Y X Y X Y X ：在ηξo 坐标系中相应为()()()()1,1,1,1,1,1,1,1----。

将有关数据代入*中的第一式，则有43214432134321243211,,αααααααααααααααα-+-=+++=--+=+--=X X X X求解上述方程组得：4,44,443214432134321243211X X X X X X X X X X X X X X X X -+-=++--=-++-=+++=αααα坐标变换方程*成为()()()()[]4321111141X X X X X ξηηξξηηξξηηξξηηξ---+++++--+++--=同理()()()()[]4321111141Y Y Y Y Y ξηηξξηηξξηηξξηηξ---+++++--+++--=当引入函数()ηξ,i N 后，坐标变换方程成为()()ii i ii i Y N Y X N X ηξηξ,,4141∑∑====式中()()()ηηξξηξi i i N ++=1141,变量ηξ、的正负号由相应节点的坐标值i i ηξ、决定。

例如当i=4时，1,144=-=ηξ，因此，()()()411,4ηξηξ+-=N 。

下面再来研究函数()ηξ,i N 的特性。

对节点()11,1Y X ，相应的自然坐标值为（-1，-1）。

从式()()()ηηξξηξi i i N ++=1141,中很容易看出，除N 1=1外，N 2=N 3=N 4=0。

对其余各节点也一样。

总而言之，对节点i(i=1,2,3,4)，除N i =1外，其余三个N 值均为零。

同时，不难看出()()()()1,,,,4321=+++ηξηξηξηξN N N N ，即四个节点的N i 函数之和等于1。

函数()ηξ,i N 具备上章所介绍的形状函数应满足的条件，可作为本单元的形状函数。

采用()ηξ,i N 做形状函数，其位移模式为()()i i i i i i v N v u N u ηξηξ,,,4141∑∑====对比()()ii i ii i Y N Y X N X ηξηξ,,4141∑∑====和()()i i i i i i v N v u N u ηξηξ,,,4141∑∑====可以看出：在这种实际单元（任意四边形）中，坐标变换式和位移模式不仅采用了相同的形状函数()ηξ,i N ，而且具有相同的数学模型。

这种性质的实际单元称为等参数单元。

对用节点位移值u i （或v i 等）求单元内某一点位移量u(或v 等)的插值公式()i i i u N u ηξ,41∑==，只要将u(或v 等)换成X(或Y等)，便成为利用节点值X i (或Y i 等)求相应点坐标X(或Y 等)的插值公式。

相反也是这样。

2.几何矩阵[B]由于几何矩阵[B]通过对位移求偏导数而得出，所以首先必须利用复合函数求导的规则得出下述公式[]⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧∂∂∂∂=⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧∂∂∂∂⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧∂∂∂∂+∂∂∂∂∂∂∂∂+∂∂∂∂=⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧∂∂∂∂y u X u J u u Y Y u X Xu Y Y u X X u u u ηξηηξξηξ或写成式中[]⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡∂∂∂∂∂∂∂∂=ηηξξY XY X J ，此式称为雅可比矩阵。

为了将几何矩阵[B]写成变量ηξ、的函数，必须将式[]⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧∂∂∂∂=⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧∂∂∂∂y u X u J u u ηξ改写成 []⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧∂∂∂∂=⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧∂∂∂∂-ηξu u J y u X u 1，同理，[]⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧∂∂∂∂=⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧∂∂∂∂-ηξv v J y v X v 1 从表示单元内各点位移与其应变关系的几何方程可知：{}TY v Xv Yuu v u X YY X ⎭⎬⎫⎩⎨⎧∂∂∂∂∂∂∂∂⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎭⎬⎫⎩⎨⎧⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡∂∂∂∂∂∂∂∂=ξε0110100000010将式[]⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧∂∂∂∂=⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧∂∂∂∂-ηξu u J y u X u 1和[]⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧∂∂∂∂=⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧∂∂∂∂-ηξv v J y v X v 1合并，则[][]TTv v uuJ J Y v Xv Yu Xu ⎭⎬⎫⎩⎨⎧∂∂∂∂∂∂∂∂⎥⎦⎤⎢⎣⎡=⎭⎬⎫⎩⎨⎧∂∂∂∂∂∂∂∂--ηξηξ11对单元（e ），任意一点的位移u ，v 对自然坐标ηξ，的偏导数可利用上式求出，写成矩阵形式为：[]{}()e p TN v v u u ψηξηξ=⎭⎬⎫⎩⎨⎧∂∂∂∂∂∂∂∂ 式中{}(){}Te v u v u v u v u 11111111=ψ[][][][][][][][][][][][][][][][][]⎥⎦⎤⎢⎣⎡=pppppppppN N N N N N N N N 4321432100000000对于i=1,2,3,4 []Ti i ip N N N ⎭⎬⎫⎩⎨⎧∂∂∂∂=ηξ将[][]TTv v u uJ J Y v Xv Yu Xu ⎭⎬⎫⎩⎨⎧∂∂∂∂∂∂∂∂⎥⎦⎤⎢⎣⎡=⎭⎬⎫⎩⎨⎧∂∂∂∂∂∂∂∂--ηξηξ110和[]{}()e p TN v vu u ψηξηξ=⎭⎬⎫⎩⎨⎧∂∂∂∂∂∂∂∂代入{}TY v Xv Y uu v u X Y Y X ⎭⎬⎫⎩⎨⎧∂∂∂∂∂∂∂∂⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎭⎬⎫⎩⎨⎧⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡∂∂∂∂∂∂∂∂=ξε01101000000100，则可得出表示在整体坐标系中位移和应变关系的几何方程：{}()[](){}()e e e B ψε= 式中的几何矩阵[B]是自然坐标ηξ，的函数：[][][][]p N J J B ⎥⎦⎤⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=--1100011010000001也可利用[]Ti iip N N N ⎭⎬⎫⎩⎨⎧∂∂∂∂=ηξ求得的[]ip N 以及()()ii i ii i Y N Y X N X ηξηξ,,4141∑∑====和[]⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡∂∂∂∂∂∂∂∂=ηηξξY XY X J 求出[]J ，[][][][][]{}Tp p p p Y Y Y Y X X X X N N N N J ⎥⎦⎤⎢⎣⎡=432143214321。