等参数单元.doc

第四章 等参数单元为了方便应用和提高计算精度，目前许多实用程序采用了等参数单元，取得了较好的效果。本章从平面问题的任意四边形单元入手，介绍等参数单元的一些基本概念，并根据工程实际应用需要，重点介绍空间六面体等参数单元分析。4-1 等参数单元的概念图4-1 任意四边形单元一、平面等参数单元1. 四节点四边形等参数单元在平面问题的有限单元法中，最简单和最常用的是三节点三角形单元，其次是四节点矩形单元。由于三节点三角形单元采用的是线性位移模式，是实际位移分布的最简单逼近形式，求解精度受到限制。而四节点矩形单元，由于它的位移函数是坐标的二次函数，单元内的应力不是常量而是线性变化的，所以能比简单三角形单元较好地反映实际应力变化情况。但是矩形单元不能适应曲线边界和非正交的直线边界，也不便随意改变大小，应用范围受到很大限制。如果采用任意四边形单元(如图4-1所示)，而仍采用矩形单元的位移模式，则基本上能够克服矩形单元的上述不足之处。但是此时在不平行于坐标轴的边界上，由于y=kx(k0)，位移函数为是坐标的二次函数，这样就不能由边界上二个节点的位移来唯一地确定该边界上的位移，故位移的连续性将得不到保证，其变形协调条件就得不到满足。采用坐标变换可解决这一矛盾，现说明如下。图4-2 不恰当的单元划分图4-3 矩形单元在图4-1所示的任意 四边形单元上，用等分四边的两族直线分割该四边形。以两族直线的中心为原点(=0)，并令四边上的值、值分别为1,这样就得到一个新的坐标系，单元上的任一点都取一个新的坐标(，)。这里的，是一种局部坐标，只适用于一个单元的范围内。与此相反，原坐标x，y则是一种整体坐标，和以前一样地通用于所有单元的整体。为确保局部坐标，和原坐标x,y有一一对应关系，即存在确定的坐标变换关系，应使任意四边形不能大歪斜，它的任意 一条边的延伸线不能再分割单元(如图4-2)。简单地说，就是要避免不必要的不规则，对四边形和六面体等参单元来讲，应使内角接近90，且总是在0180。为了较简单地得出任意四边形符合收敛准则的位移函数及局部坐标和整体坐标之间的变换式，我们引入一个四节点正方形单元(图4-3)。由第二章知，此单元的位移函数为 (4-1)并且有 (4-2)其中形函数 (i=1,2,3,4) (4-3)而(i，i)为节点i的坐标，ui，vi为节点i的位移分量。式(4-2)说明该矩形单元在采用此种位移模式满足常应变准则。并且，由于位移函数是，的双线性函数，使单元位移在每条边上是或的线性函数，可以被该边上二节点的位移唯一确定，从而保证了相邻单元的位移连续性所以，式(4-1)及(4-2)所确定的位移模式在图4-3坐标下是满足收敛条件的。现在，把式(4-1)所示位移模式和式(4-3)所示的形函数移用于图4-1所示的四边形中的局部坐标，ui、vi理解为该四边形节点处的实际位移分量。则显然可见，式(4-1)在四个节点处给出节点位移，也就是(4-1)式给出了单元位移的某种模式；并且，在单元的四条边上，位移是线性变化的，从而保证了位移的连续性。因此，式(4-1)就是所需要的正确的位移模式。式(4-1)得给出了用局部坐标表示的单元位移模式，而实际运算时，需要的是用整体坐标表达的位移函数。因此还必须给求出整体坐标和局部坐标间的坐标变换式。仿照位移模式(4-1)，把坐标变换式取为 (4-4)显然可见，该式在四个节点处给出了节点的整体坐标，并且局部坐标系上的一平行于坐标轴的直线，映射到整体坐标上，正是一条过相应二条边上相应等分点的一条直线。例如，局部坐标上的一条直线，由于形函数的双线性特性，映射到整体坐标上仍为一条直线，而且这条直线两端点分别是局部坐标点(,1)与(，-1)的变换对应点。同样由于形函数的双线性特性，在=1变换得到的直线上，整体坐标x,y均是的线性函数，因此上述两端点一定就是这两条边的相应的的等分点。局部坐标系中的=1(或1)的直线，映射到整体坐标系中为任意四边形的相应边线。因此，式(4-4)正是需要的坐标变换关系，还能保证变的得到的单元间即不相互重叠，又不裂开。这样，我们看到，单元的位移函数式和坐标函数式(坐标换式)具有完全相同的形式，它们用同样数目的节点值为参数，并以完全相同形函数Ni(，)为系数。当参数取节点位移时，就得到单元的位移函数式，当参数取为节点坐标值的，就得到单元坐标的表达式。这样的单元称为参数单元，简称等参元。由前几章分析中知道，平面问题的三角形单元，矩形单元及空间问题的四面体单元、六面体单元的位移函数式和坐标式均有相同的形式，只不过作为系数的形函数是用整体坐标表示而已。所以它们均是等参数单元，但一般不这样讲。平时说的等参元是指经过坐标变换得到的单元。图4-4 八节点曲边四边形单元图4-5 八节点矩形单元在上面等参元的位移模式及坐标变换式的建立过程中，直接借用了规则的正方形单元的位移模式表达式。这种方法不仅适用于简单的线性形函数，还适应于高次形函数。在这里，局部坐标系中被 借用的规则单元称为基本单元或母单元，而由其变换得到的单元为实际单元或子单元。此外可以证明：如果位移模式在基本单元中能够充分反映刚体位移和常量应变，则它在等参单元中也能充分反映刚体位移和常量应变。2. 八节点曲边四边形等参数单元。为提高计算精度，还可以采用更多节点的等参单元。通常用得较多的是八节点曲边四边形等参单元(图4-4)。利用上面介绍的坐标变换方法，可以导出八节点等参单元的位移函数。 现取一边长为2的正方形八节点单元为母单元，如图4-5所示。其位移函数为 (4-5)其中形函数为 (4-6)式中的i，i为节点i的坐标。由于形函数是坐标，的双二次函数，位移在和为常数的单元边线上，只是或的二次函数，可以被边界上三个节点的位移唯一确定，这就保证了相邻单元的位移的连续性。同时，因为有即单元满足“常应变准则”，所以以上位移模式对母单元是完全完备协调条件的。若将其以局部坐标的形式应用于实际单元，也是能够满足收敛准则的。和四节点等参数单元一样，根据单元位移函数公式(4-5)就可以写出，由局部坐标(，)变换到整体坐标(x,y)的坐标变换公式： (4-7)由于式(4-5)能保证单元位移的连续性，(4-7)式也就可以保证变换得到的实际单元的连续性。为看清局部坐标系中八节点正方形单元映射到整体坐标系中的形状，只需看任意一条边线就可以了。如图4-4中的284边，其在局部坐标系中是=1，将其代入(4-7)式，可以得到284边在整体坐标系下的参数方程 (4-8)其中是参数，消去，可知这曲线为一条抛物线。此外，由形函数的性质知，2，4，8三点变换到整体坐标系中即为各自实际坐标决定的点。所以，(4-8)式所决定的抛物线必然通过实际坐标系中的284三点，且为此三点唯一决定。同理，其他三条边也是被相应三个节点唯一决定的抛物线。因此，在整体坐标系中的实际单元形状是由过八个节点的四条抛物线围成的曲边四边形。在实际计算中用到的只是每个单元在整体坐标下八节点的坐标(xi，yi)(i=1，2，8)。因此，在整体坐标下进行单元划分时，只须给出每个单元八个节点的实际坐标值，而不必画出抛物线曲边的精确形状，示意画出即可。和任意四节点四边形等参单元一样，为保证坐标变换能顺利进行、减小计算误差等，一般要求单元不要太歪斜，各边线上的中间点应尽量处于角点连线的中央部位。以上分析方法，可以推广到具有更高精度的更多节点的等参单元。但在实际应用中，平面八节点等参单元的精度已经足够了。图4-6 八节点六面体单元二、空间等参单元对于空间问题的等参数单元，可以由平面问题等参数单元分析中叙述的方法加以推广。1. 八节点六面体等参数单元为分析任一空间八节点六面体等参数单元，取如图4-6所示的局部坐标下边长为2的正六面体八节点单元为基本单元，位移函数为 (4-9)其中形函数表达式为 (i=1,2,，8) (4-10)并且有 (4-11)和 (4-12)坐标变换式可写为： (4-13)由(4-9)和(4-10)可以看到，局部坐标下，当一个自变量(坐标)固定时，单元的位移函数是另外二个坐标的双线性函数。所以在六面体的每一个侧面上，单元位移函数完全由其上的四个节点的位移值所唯一确定。因此，在相邻单元的公共面上位移是连续的。母单元上位移函数的连续性，同样保证了坐标换式(4-13)的连续性及整体坐标单元的位移函数的连续性。整体坐标系中实际单元的具体形状，同样可以通过分析单元上某些特征线、面的映射情况来说明。如母单元上任意平行某一坐标轴的直线，由于在其上两个局部坐标值是固定的，实际单元中其映射线上的实际坐标x,y,z均为变化着的局部坐标的线性函数，因此该线仍为一直线，原来线上的某一等分点，仍为相应的等分点。所以，实际单元每条棱边为两相应节点的直线连线。在局部坐标下，单元某一侧面上二条相对直线棱边的某一等分点的连线(平行于某一坐标轴)，也必对应着整体坐标系下相应棱边上对应等分点的连线，如图4-7。不过，由于实际单元的对应棱边一般是不平行的，局部坐标下的每一侧面经过坐标变换在整体坐标下的映射面一般也不再是一平面，而是由二组直线所组成的双曲抛物面。 图4-7 空间八节点六面体单元 图4-8 二十节点正六面体单元2. 二十节点曲六面体等参数单元基本单元取如图4-8所示的边长为2的二十节点正六面体单元。仿照前述方法可以得出单元的位移函数和坐标变换式为： (4-14) (4-15)其中形函数的统一表达式为 (4-16)二十节点的曲六面体等参单元在整体坐标下的形状是相当复杂的。例如棱边3 19 7，它在局部坐标下的方程=1，=1；代入(4-15)式，就得到其在整体坐标下的方程为消去，可以看出其表示两个抛物线柱面的交线，它完全由此交线上的三个节点坐所决定。因而，在整体坐标下，单元的棱边一般不再是直线，而为一二次曲线；而其侧面也是由两族这样的二次曲 线所交织构成的曲面。和平面八节点曲边四边形等参单元一样，对空间六面体等参数单元，也只需要给出各节点的精确坐标值，而不用了解单元的精确形状。为使坐标变换顺利进行 ，保证计算精度，也要求空间等参单元的形状不要太歪斜，节点的布置应尽量接近基本单元上节点的分布。按照上述方法，可以推广到具有n个节点的空间等参单元，其位移函数和坐标变换式可以写成一般形式： (4-17) (4-18)其中Ni(，)是局部坐标，的显函数，在不同的单元中取不同的形式。4-2 空间等参数单元的数学分析在进行空间等参数单元的力学分析时，需要用到各个形函数对整体坐标的导数，局部坐标系中微分体的体积及微分面的面积，以及局部坐标面的法向余弦。现在来分别导出这些参数的表达式。一、形函数对整体坐标的导数由上节(4-17)和(4-18)式知，等参数单元的形函数Ni及x,y,z均只是局部坐标，的显函数，所以，利用下式 (x,y,z)无法求出形函数Ni对整体坐标x,y,z的导数。但是，我们注意到等等，所以有 (4-19)从而有 (4-20)这里 (4-21)称为雅可比矩阵，为了求得这个矩阵，必须将(4-18)式代入上式。于是得 (4-22)图4-9 局部坐标微分体能够保证式(4-20)顺利运算的充要条件就是雅可比矩阵行列式J的值在整个单元上任意点均不为零，这也是等参单元可以进行坐标变换的必要条件。具体地讲，这一条件就是前面所说的，在进行单元划分时，每个单元顶角的配置应限制在0180，一般应尽量控制在90左右；节点的配置也应尽量和基本单元相似。这样，则既保证坐标变换的条件，又保证得到足够的计算精度。二、局部坐标系中微分体的体积试在单元内的任一点P，沿局部坐标、的方向分别作微分矢量、 (分别对应局部坐标的d，d，d)，在实际空间形成一平行六面体(图4-9)。由于在方向，只是坐标在变化，而及支持不变，所以微分矢量在整体坐标轴上的投影分别为同样可得矢量及在整体坐标轴上的投影：由矢量运算，可知此三矢量构成之平行六面体的体积为此即为局部坐标系中微分体ddd映射到整体坐标系中相应微分体的体积。将上式化简，最后得。 (4-23)即 (4-24)式中J为前述雅可比矩阵的行列式，即 (4-25)三、局部坐标系中微分面的面积局部坐标系中平行于坐标 面的微分面dd，映射到整体坐标系中即为图4-9中矢量、所构成的平行四边形。由矢量运算知，该四边形的面积为和矢量积的模，即 (4-26)由前面导出的微分矢量在各坐标轴上的投影，极易求得： (a,b,c；,) (4-27)其中 (，)也可求得： (a,b,c;,) (4-28)其中 (，) 将(4-27)式和(4-28)式代入(4-26)式，可得 (a,b,c;,)此即为局部坐标系中微分面dd映射到实际坐标系中微分面的面积，即 (4-29)四、局部坐标面的法向余弦微分矢量与所构成平面的法线方向，即为矢量的方向，也就是局部坐标面的法线方向(但一般并不是坐标轴的方向，因为整体坐标下，一般并不是正交坐标)。利用矢量运算公式，有 (4-30)其模即为与所成平面的面积，也就是上面所求dA。因此，可由上列各式求得局部坐标面的法向余弦表达式： (，) (4-31)4-3 空间等参数单元的力学分析现在来对n个节点的空间等参数单元进行力学分析，从而建立这种单元的载荷列阵，应力矩阵和刚度矩阵。一、等效节点载荷计算空间等参数单元的载荷移置，仍遵循静力等效原则，因而前面几章介绍的普遍公式仍然适用。1. 集中力当单元在任意点受有集中力时，等效节点载何为 (4-32)但在这里 (4-33)而矩阵N可用三阶单位矩阵I及各个形函数表示成为 (4-34)由于其中Ni是由受力点的局部坐标、表示的，为计算方便，单元划分时应将集中力作用点取为节点。否则 ，由整体坐标求局部坐标很不方便。2. 体力单元受有分布体力，其强度为，此时有上式中N是局部坐标的显函数，且单元的真实形状及边界描述起来非常困难，一般将上式化成局部坐标积分比较方便。利用(4-24)式，有 (4-35)对于变量的分布体力，须将体力分量X、Y、Z表示成局部坐标，的函数，再进行积分。3. 面力当单元在某一边界上(例如在=1的面上)，受有分布的面力时，可利用(4-34)式的积分求得等效节点载荷应用(4-29)式，可以将上式表示成为 (4-36)对于变量的分布面力，须将面力分量、表示成局部坐标的函数，再进行积分。如果上述分布面力为静水压力q，则将(4-31)式代入，然后再代入(4-36)式，得 (4-37)当面力作用在的面上或1的面上时，只须将(4-37)中的，轮换，即可得出各种状况下的载荷列阵表达式。二、应力矩阵将位移函数表达式(4-17)代入空间问题的几何方程，可得单元内的应变表达式 (4-38)其中是单元的节点位移列阵，即 (4-39)而 (i=1,2,，n) (4-40)和以前一样，单元内的应力可以表示成 (4-41)将其中应力矩阵S写成分块形式 (4-42)则有 (i=1,2,，n) (4-43)三、刚度矩阵将单元上的节点力表示成为则由虚功方程可以得出和以前一样的公式其中单元刚度矩阵为 (4-44)单元刚度矩阵是一个3n3n的矩阵，可以写成分块形式： (4-45)其中每一个子矩阵是33方阵，为 (4-46)(i=1,2,n;j=1,2,n)有了单元的载荷列阵和刚度矩阵之后，即可利用以前多次用过的方法，建立那些求解节点位移所需的节点平衡方程，在求出节点位移之后，再按(4-41)式求单元的应力。应当说明是，对于等参单元中的任一点，由它的局部坐标求出它的整体坐标时，只须直接应用(4-18)式；但是，要由整体坐标算出局部坐标，则须利用(4-18)式，求解三元非线性联立方程组，这是非常繁琐的。因此，在计算单元中的应力时，一般是先设定一组局部坐标（一般是节点或积分点），代入应力矩阵以求出应力；并根据这组局部坐标算出相应的整体坐标，从而可知所求应力是单元(构件)那一点处的应力。对于单元中以整体坐标标定的一点（除非是节点或积分点等特殊点），要计算应力则是非常复杂和不现实的。计算节点处的应力，则比较方便。因为其局部坐标是已知的，且非常简单，不外乎0，1以及简单的分数，等等。因此，在整理应力计算结果时，一般宜采用绕点平均法。但应当指出，如果单元的形态比较差，则计算求得的节点处的应力的表征性也就较差，通过平均后其表征性可能仍然不够好。对于这样的单元，就应整理单元内点处的应力，或由内点处的应力来外插求节点的应力。对于等参数单元，单元内点的应力精度优于边界点的精度，角节点的应力精度最差。4-4 高斯积分法及其应用由4-3知，在计算空间等参数单元的载荷列阵及刚度矩阵时，需用到如下形式的定积分：； 其中被积分函数f(，)一般是很复杂的，即使能够得出它的显式，其积分也是很繁的。因此，一般用数值积分来代替函数的定积分，即，在积分区域内按一定规则选出一些点，称为积分点，算出被积函数f(，)在这些积分点处的值，然后再乘以相应的加权系数并求和，作为近似的积分值。数值积分的方法有多种，其中高斯积分法可以用相同的积分点数达到较高的精度，或者说用较少的积分数达到同样的精度。一、高斯积分法首先介绍一维积分的高斯公式，应有 (4-47)其中f(i)是被积函数在积分点i处的数值，Hi为加数系数，n为积分点数目。可以证明，对于n个积分点，只要选取适当的加数系数及积分点位置，能够使(4-47)式在被积分函数为不超过(2n-1)次多项式时精确成立。例如，n=1时 (a)不论f()的次数是0还是1，只需取H=2，上式均是精确成立的。因为 (b) (c)当n=2时，能保证(4-47)式精确成立所允许的多项式的最高次数是3，此时，f()的通式为 (d)其精确积分为 (e)数值积分为 (f)为了在C0C3取任意值(包括取零值在内)时公式(f)是精确的，显然应有， ， 所以，应取同样，对于不超过五次的多项式，只要取n=3即可保证得到精确的积分值。对于n=1到n=6时，高斯求积公式中积分点坐标i及加权系 Hi的数值列于表4-1。以一维高斯积分公式(4-47)为基础，极易导出二维及三维公式。求二维重积分的数值时，可以先对进行积分，这时将当作常量，于是由(4-47)式得到 (g)再对进行积分，得出 (h)将式(g)代入，即得或改写为 (4-48)这就是二维的高斯积分公式。表3-1 高斯求积公式积分点坐标与加权系数点数n坐标i加权系数Hi1234560.0.57735 02691 896260.77459 66692 41483 1.00000 00000 000000.86113 63115 940530.33998 10435 848560.90617 98459 386640.53846 93101 056830.00000 00000 000000.93246 95142 031520.66120 93864 662560.23861 91860 831972.1.00000 00000 000000.55555 55555 555550.88888 88888 888880.34785 48451 374540.65214 51548 625460.23692 68850 561890.47862 86704 993660.56888 88888 888890.17132 44923 791790.36076 15730 481390.46791 39345 72691同样，可以求得三维高斯积分公式： (4-49)式(4-48)及(4-49)中的n，m，l是分别关于变量，的积分点数目，由各个自变量在被积函数中可能出现的最高次数分别决定，一般并不要求相同。但为应用方便，常常在各个方向取相同的积分数，即 (4-48a) (4-49a)由前面的推导可见，当在每个方向取n个积分点时，只要多项式被积函数中自变量的次数m2n-1，则用高斯求积公式求得的积分值是完全精确的。反过来，对于m次多项式的被积函数，为了积分值完全精确，积分点的数目必须取。二、合适积分点数的确定在空间等参数单元分析中，计算节点载荷及刚度矩阵时要用到高斯积分公式，下面以20节点的等参数单元为例，来探讨一下相应公式中积分点的数目。分布体力的等效节点载荷列阵计算见公式(4-35)，即 (4-35)由式(4-16)可见，形函数Ni对每个局部坐标来说，一般都是二次式，因而在元素中，每个局部坐标都可能以二次幂出现，同时，由(4-15)式又知，在整体坐标的表达式中，每个局部坐标也都可能以二次幂出现，根据(4-21)式，在中每个局部坐标可能出现的次数是5。据此，当体力分量为常量时(例如体力为均质单元自重)，(4-35)式中的被积函数对于每个局部坐标来说，可能出现的最高幂次是m=5+2=7。为了使积分值完全精确，每个坐标方向积分点的数目应为。对于三维积分，其