第7章(等参数)

改公式号、图号第七章矩形单元和等参元矩形单元也是一种常用的单元，它采用了比平面三角形单元更高阶次的位移模式，因而可以更好地反映结构中的位移状态和应力状态。等参数单元则是一类用于曲线或曲面结构、具有特殊特性的单元。在本节中，除了对上述单元进行分析之外，还对结构有限元方程的求解方法及有关数值积分方法进行了简单介绍。7.1矩形单元的基本概念如图3-21所示的矩形单元，其边长分别为2a和2b，两边分别平行于x、y轴。若取该矩形的四个角点为结点，每个结点位移有两个分量，所以矩形单元共有8个自由度。采用平面三角形单元的分析方法，同样可以完成对这种单元的力学特性分析。这里引入一个局部坐标系、，可以推出比较简洁的结果。(1)坐标变换与形函数在图3-21中，取矩形单元的形心为局部坐标系的原点，和轴分别与整体坐标轴x和y平行，两坐标系存在有以下的坐标变换关系byyaxx00(3.146)式中，2)(2)(2)(2)(2)(2)(2)(2)(////////142334124132043210yyyybxxxxayyyyyxxxxx(3.147)其中，(xi,yi)是矩形单元结点i的整体坐标，i=1,2,3,4。在局部坐标系中，结点i的坐标是(i,i)，其值分别为±1。取局部坐标系中的位移模式为uv12345678(3.148)将4个结点的局部坐标值代入上式，可以列出4个结点处的位移分量，得到一组关于8个未知参数1，2，…，8的方程组，由此可求得这8个未知参数，再把这些参数代回式(3.148)中，得到用形函数矩阵和结点位移表示的位移模式方程4141,iiiiiivNvuNu(3.149)其中4/)1)(1(00iN(3.150)式中,0=i，0=i，i=1,2,3,4。(2)单元刚度矩阵u2(U2)u3(U3)v3(V3)v2(V2)v1(V1)v4(V4)u1(U1)u4(U4)xoyo12432a2b图3-21矩阵单元由几何方程求得单元的应变vbuavaubabvaubvbuaxvyuyvxuxyyx11111(3.151)将式(3.146)代入，得BBBBe1234(3.152)式中000011100141001iiiiiiiiibaababNbNaNaNbabB(i=1,2,3,4)(3.153)可以得出用结点位移表示的单元应力，即DSSSSe1234(3.154)式中SDBii(i=1,2,3,4)(3.155)对于平面应力问题SEabbabaabiiiiiii4111111211212000000(3.156)将单元刚度矩阵写成分块形式44434241343332312423222114131211kkkkkkkkkkkkkkkkke(3.157)其中的子矩阵按下式进行计算kBDBtdxdyijiTj(3.158)设单元厚度t是常量，则jijijijijijijijijijijijijTiijabbabaabEtddSBtabk311213112121311213111421111(i,j=1,2,3,4)(3.159)对于平面应变问题，将上式中的E、分别换成E/1-2和/1-。(3)载荷列阵和结构有限元方程四边形单元的结点位移与单元结点载荷列阵{R}e之间的关系为eeeRk(3.160)由于矩形单元有四个结点，所以单元载荷列阵{R}e有8个元素，即TyxyxyxyxeFFFFFFFFR}{44332211(3.161)例如，对于单元自重为W，移置到每个结点的载荷都等于四分之一的自重，其载荷列阵为TeWR}410410410410{(3.162)如果单元在一个边界上受有三角形分布的表面力，且在该边界上的一个结点处为零，而另一个结点处为最大，可将总表面力的三分之一移置到前一个结点上，而将其三分之二移置到后一个结点上。与平面三角形单元相同，将各单元的ek][、{}e和{R}e都扩充到整个结构自由度的维数，再进行叠加，可得到整个结构的平衡方程。即，[K]{d}={R}。(4)矩阵单元的特性矩形单元的位移模式比平面三角形单元的线性位移模式增添了项（相当于xy项），这种位移模式称为双线性模式。在这种模式下，单元内的应变分量不是常量，这一点可以从应变矩阵[B]的表达式中看出。矩形单元位移模式中的1、2、3、5、6、7与平面三角形单元相同，反映了刚体位移和常应变。位移模式在单元的边界上(=±1或=±1)，位移是按线性变化的，显然，在两个相邻单元的公共边界上，其位移是连续的。由矩形单元的应力矩阵表达式可以看出，矩形单元中的应力分量也都不是常量。其中，正应力分量x的主要项（即不与相乘的项）沿y方向线性变化，而正应力分量y的主要项则是沿x方向线性变化、剪应力分量xy沿x及y两个方向都是线性变化。因此，采用相同数目的结点，矩形单元的精度要比平面三角形单元的精度高。但是，矩形单元存在一些明显的缺点：其一是矩形单元不能适应斜交的边界和曲线边界，其二是不便于对不同部位采用不同大小的单元。6.2平面等参数单元利用边界都是直线的单元分析结构复杂的曲边或曲面外形，只能通过减小单元尺寸、增加单元数量进行逐渐逼近。如果这些单元的位移模式是线性的，问题的求解精度也受到限制。为了克服以上缺点，人们发展了等参数单元（等参元）的思想，一方面，这类单元能很好地适应曲线边界或曲面边界，相对准确地逼近结构形状；另一方面，这类单元可以具有较高阶的位移模式，能够更好地反映结构的复杂应力分布，即使单元网格划分比较稀疏，也可以得到较好的计算精度。等参元的基本思想是：首先导出关于局部坐标系的规整形状的单元（母单元）的高阶位移模式，然后利用形函数多项式进行坐标变换，得到关于整体坐标系的复杂形状的单元（子单元），其中子单元的位移函数插值结点数与其位置坐标变换的结点数相等，位移函数插值公式与位置坐标变换式都采用相同的形函数与结点参数，这样的单元称为等参元。（1）母单元及其形函数首先，根据形函数的定义，在局部坐标中，建立起几何形状简单且规整的单元，称之为母单元。对于一维问题，母单元采用局部坐标ξ，11。形函数的具体形式如下：1）线性单元（2结点）NN121212(3.163)2）二次单元（3结点）232112121NNN(3.164)对于二维问题，二维母单元是,平面中的2×2正方形，11，11，如图3-23所示，坐标原点在单位形心上。单元边界是四条直线：1，1。为保证用形函数定义的未知量在相邻单元之间的连续性，单元结点数目应与形函数阶次相适应。因此，对于线性和二次形函数，单元每边的结点数分别为两个和三个。（a）线性单元（b）二次单元图3-23二维母单元1)线性单元（4结点）的形函数如下：Ni11400（i=1,2,3,4）(3.165)其中，00ii.2)二次单元（8结点）的形函数为角点：Ni141110000（i=1,2,3,4）(3.166a)边中点：Ni121120（i=5,7）（3.166b）Ni121120（i=6,8）（3.166c）（2）等参坐标变换等参元需要用坐标变换把形状规整的母单元转换成具有曲线（面）边界的、形状复杂的单元。转换后的单元称为子单元。子单元在几何上可以适应实际结构的各种复杂外形。即可以采用各种形状复杂的子单元在整体坐标系中对实际结构进行划分。子单元通过坐标变换映射成一个局部坐标系下的规整的母单元。坐标变换是指在局部坐标,,和整体坐标xyz,,之间建立一一对应关系。在这里，坐标变换关系利用形函数建立起来。在整体坐标系中，子单元内任一点的坐标用形函数表示如下η-101ξ21η-101ξ(a)线性单元(b)二次单元图3-22一维母单元xNxNxNxyNyNyNyiiii,,,,,,11221122(3.167)其中，Ni,是用局部坐标表示的形函数，xyii,是结点i的整体坐标，上式即为平面坐标变换公式。图3-24表示了二维单元的平面坐标变换。母单元是正方形，子单元变换成曲边四边形，且相邻子单元在公共边上的整体坐标是连续的。以二次单元为例，两个相邻单元在公共边界上都是二次曲线（抛物线），而在三个公共结点上具有相同的坐标。因此，整个公共边界都有相同的坐标，即相邻单元是连续的。（a）母单元（b）子单元图3-24二维单元的平面坐标变换（3）两种坐标系的关系——雅可比矩阵讨论局部坐标系和整体坐标系下的偏导数的关系。根据复合函数的求导法则，有xyxyxyxy(4.58)上式写成矩阵形式xyJ(4.59)其中，J为雅可比（Jaccobi）矩阵xyxyJ(4.60)式(4.59)表示由x和y推导，的变换式，其逆变换式为1xyJ(4.61)其中，J-1是雅可比矩阵J的逆矩阵11yyxxJJ(4.62)（4）平面等参元的有限元分析上述二维二次单元的位移模式为8181,,,iiiiiivNvuNu(4.63)其中，ui和vi是结点i的位移。形函数的具体形式见式(4.54,4.55,4.56)。采用坐标变换使母单元的8个结点ii,与等参元的8个结点整体坐标值（xi,yi）一一对应。整体坐标xy,和局部坐标,的变换式为8181,,,iiiiiiyNyxNx(4.64)其中，Ni,与母单元的形函数相同，见式(4.54)、(4.55)、(4.56)。这样就确定了平面八结点曲边四边形等参元的几何形状和位移模式。将上述等参元的位移模式代入平面问题几何方程，得到单元应变分量的计算式128eeuxvyuvyxεBδBBBδ(4.65)其中，Te}{821δδδδ—单元结点位移列阵，iiivuδ，(i=1,2,…,8)。为了求得应变矩阵B，需要进行如下推导。由于形函数Ni,是局部坐标的函数，需要进行偏导数的变换1iiiiNNxNNyJ(4.66)式中雅可比矩阵的逆矩阵1J由式（4.62）给出，其元素根据坐标变换式确定，即iiiiiiiiiiiixNxxNxxNx