1、 等参元的基本思想是：首先导出关于局部坐标系（Local coordinate, 或Natural coordinate, 自然坐标系）的规整形状的单元（母单元）的高阶位移模式，然后利用形函数多项式进行坐标变换，得到关于整体坐标系(Global coordinate)的复杂形状的单元（子单元），其中子单元的位移函数插值节点数与其位置坐标变换的节点数相等，位移函数插值公式与位置坐标变换式都采用相同的形函数与节点参数，这样的单元称为等参元。 等参数单元（Isoparametric elements）简称等参元，是根据特定方法设定的一大类单元，不一定具有相同的几何形状。因为等参元具有规范的定义原理和
2、较强的适应复杂几何形状的能力。在有限元理论中占有重要的地位。采用等参元，一方面能够很好地适应曲线边界和曲面边界，准确地模拟结构形状；另一方面，等参元一般具有高阶位移模式，能够较好地反映结构的复杂应力分布情况，即使单元网格划分比较稀疏，也可以得到比较好的计算精度。 (1) (1) 局部坐标系下的位移模式局部坐标系下的位移模式根据形函数的定义，在局部坐标系中，建立起几何形状简单且规整称之为母单元。下面以四边形单元为例说明等参元的基本概念。111111母单元是平面中的22正方形， 坐标形函数定义原点在单位形心上。单元边界是四条直线：,函数阶次相适应。对于具有线性形函数的四边形单元，共有四个节点，如图
3、6-1所示。三个，如图6-2。的单元，如图6-1所示，为保证用的未知量在相邻单元之间的连续性， 单元节点数目应与形如果是二次函数的四边形单元，单元每边的节点数为（a）母单元 （b）子单元图6-1 线性矩形单元及其平面坐标变换 对于如图6-1所示的线性4节点四边形等参元（Bilinear element），它在局部坐标系形函数如下：44443333222211111111, 1,N代入节点坐标值41100iN（i=1, 2, 3, 4） (6.1)isoparametricii00 其中 ， . 该形函数是定义成自然坐标下的归一化变量 的函数。,对于8节点二次四边单元，角点上的形函数为11141
4、0000iN （i=1, 2, 3, 4） (6.2a) 边中点上的形函数为Ni121120（i = 5, 7） （6.2b） Ni121120（i = 6, 8） （6.2c） 可见，线性四边形单元和8节点二次四边单元用局部坐标形函数表达的位移模式如下11,4,8nniiiiiiuNu vNvn (6.3)（2）等参坐标变换等参元需要用坐标变换把形状规整的母单元转换成具有曲线（面）边界的、形状复杂的单元。转换后的单元称为子单元。子单元在几何上可以适应实际结构的各种复杂外形。即可以采用各种形状复杂的子单元在整体坐标系中对实际结构进行划分。,zyx,子单元通过坐标变换映射成一个局部坐标系下的规整
5、的母单元。坐标变换是指在局部坐标和整体坐标关系。在这里，坐标变换关系利用形函数建立起来。例如，对于上述4节点线性四边形单元，有4112244141122441,iiiiiixNxNxNxNxyNyNyNyNy (6.4a)之间建立一一对应 (6.4b)对于上述8节点二次四边形单元，有8112288181122881,iiiiiixNxNxNxNxyNyNyNyNy ,iNiiyx ,其中，是用局部坐标表示的形函数，是节点i的整体坐标上式即 为平面坐标变换公式。 如图6-1和图6-2所示的二维单元的平面坐标变换，其中母单元是正方形，子单元变换成曲边四边形，且相邻子单元在公共边上的整体坐标是连续的
6、，且在公共节点上具有相同的坐标，即相邻单元是连续的。 （3） 两种坐标系的关系雅可比矩阵（Jacobian matrix） 局部坐标系和整体坐标系之间具有如下偏导数的关系。根据复合函数的求导法则，有yzNyyNxxNNiiii (6.5)(6.6) 上式写成矩阵形式 yNxNJyNxNyxyxNNiiiiii 其中，J称为雅可比矩阵（Jacobian matrix） 将式(6.4)的表达式代入雅可比矩阵式(6.7)，求得雅可比矩阵。对于式(7.4a)表示的双线性四边形单元，有 4411, iiiiiiNNxyxy(6.7)所以 1xyJ(6.9)例6-1：对一个四边形单元，在整体坐标系下的4个
7、节点（1, 2, 3, 4）的坐标分 别是(a, b), (c, b), (c, d), (a, d)，求其雅可比矩阵及其逆矩阵。 解：根据式(6.1)、(6.4a) 和式(6.7)可以得到雅可比矩阵各元素)1 ()1 ()1 ()1 (41432111xxxxxJ11yyxxJJ(6.8)因为 的表达式为4 , 3 , 2 , 1,iyxii 为4个节点在整体坐标系下的坐标，将各坐标值代入得到也可比矩阵，具体为0202c adbJ进一步，由式(6.8)可求得其逆矩阵为12020c ad bJ)1 ()1 ()1 ()1 (41432112yyyyyJ)1 ()1 ()1 ()1 (41432
8、121xxxxxJ)1 ()1 ()1 ()1 (41432121yyyyyJ就可以 以二维二次8节点单元为例，说明等参元分析的一般原理。8节点四边形等参元的位移模式为8181, ,iiiiiivNvuNu(6.10)其中，ui和vi是节点i的位移。ii,x y, 采用坐标变换使母单元的8个节点与等参元的8个节点整体坐标和局部坐标的变换式为yi）一一对应。整体坐标8181, ,iiiiiiyNyxNx(6.11)将上述等参元的位移模式代入弹性力学平面问题的几何方程，将会得到如下形式的、用应变矩阵B表示的单元应变分量计算式值（xi,128xeeyxyuxvyuvyxBBBB (6.12)为了求得
9、应变矩阵B，进行如下推导。由于形函数,iN是局部坐标的函数， 需要进行偏导数的变换 1iiiiNNxNNyJ(6.13)其元素根据坐标变换式确定，即88118811;iiiiiiiiiiiiNNxyxyNNxyxy(6.14) 将单元任一点的应变列阵代入平面问题的物理方程，得到单元应力列阵。再利用虚功原理，进一步推导出单元刚度矩阵1111d dd deTTt x yt kB DBB DB J 式中，t为单元厚度。注意上式中的积分限，在整体坐标系中的积分相应地转化成了局部坐标系下的积分，是积分限为-1到+1的定积分。 (6.15) (6.16) 在等参元分析中，单元外载荷的计算分析如下。设外载荷
10、包括信集中载荷、体积力、表面力等，可以写成如下形式11NNeeeeee RRFQPF Q P 其中，对于集中载荷，设单元任意点c作用有集中载荷TxyGGG，移置到单元有关节点上的等效节点载荷为(6.17) 1,2,8eee TiixiyicFFNiFG (6.18) 对于体积力，设单元上作用的体力为Txyp pP，移置到单元各节点载荷为11111,2,8TxeeeiixiyiiypPPN tdxdyNt J d dpi Pp 对于表面力，设单元某边界上作用的表面力为Txyq qq，则这条边上的等效载荷为xeee TiixiyiyqQQNtdsq Q (6.19)式中, 是单元作用有面力的边界域
11、, ds是边界域内的微段弧长。在上述分析的基础上，利用结构中所有等参元的单元刚度矩阵集成结构整体刚度矩阵。列写结构有限元方程、引入约束条件，进而进行结构整体分析。上的等效三个节点将原来在整体坐标系统下的三角形平面单元，通过坐标转换，可以转换到局部坐标系中去,设在整体坐标系下的节点分别为( , )( , )( , )iijjkkx y x y x y、，对应到局部坐标系 中的坐标分别为（-1，1）、（-1，-1）、图6-3 三角形单元在整体坐标系和局部坐标系下的对应变换设坐标转换方程式为 （1，-1），如图6-3所示。(6.21)将各点坐标值带入上式得到：1231231 1111 11110 0
12、 1111xxxa b cyyyd ef 由上式可解出：11231231 1111 111111111xxxxyyyy (6.22)所以11231231 111 1 10 0 11 1 1 1 1 1a b cx x xd e fy y y J(6.23)10011xabcydef (6.20)(6.25)对于三角形单元，坐标转换方程式分别对与 求偏微分有如下xaxbydye代入上式可得：deta bdxdyd dd e (6.26)平面三角形单元的应变能公式可以写成：1111121det2TTTTUdAJ dd B DBB DB(6.27)三角形单元的面积的积分表达式可以写成如下形式 111
13、11111det|xxdAdxdyd dd dyy J(6.24)具体表达式dd11T11eJDBBK(6.29)VUW0V设位能等于应变能与外力功之和，即，由最小势能原理, 可知 0UW上式可写成：(6.30)其中，K即为单元刚度矩阵，如下(6.31)式中为节点位移。将应变能U对求微分得到：1111det TUJ d d B DB(6.28)0eeRK 如图6-4所示的矩形单元，其边长分别为2a和2b，两边分别平行于x、y轴。若取该矩形的四个角点为节点，每个节点位移有两个分量，所以矩形单元共有8个自由度。采用平面三角形单元的分析方法，同样可以完成对这种单元的等参数有限元分析。这里引入一个局部
14、坐标系、，可以推出比较简洁的结果。 u2 (U2)u3 (U3)v3 (V3)v2 (V2)v1 (V1)v4 (V4)u1 (U1)u4 (U4)xoyo12432a2b图6-4矩形单元(1) 坐标变换与形函数 在图中，取矩形单元的形心为局部坐标系的原点，和轴分别与整体坐标轴x和y平行，两坐标系存在有以下的坐标变换关系byyaxx00(6.32)式中2)(2)(2)(2)(2)(2)(2)(2)(/142334124132043210yyyybxxxxayyyyyxxxxx(6.33)在局部坐标系中，节点i的坐标是(i , i )，其值分别为1。取局部坐标系中的位移模式为uv 1234567
15、8(6.34)将4个节点的局部坐标值代入上式，可以列出4个节点处的位移分量，得到一组关于8个未知参数1，2，8的方程组，由此可求得这8个未知参数，最终可以得到用形函数矩阵和节点位移表示的位移模式方程4141 ,iiiiiivNvuNu4/)1)(1(00iN(6.35)(6.36)(2) 单元刚度矩阵由弹性力学几何方程求得单元内任一点的应变分量为11111xyxyuuubaxvvvaybabuvuvuvabbayx (6.37) (6.38)1234eBBBB 将式代入（6.35）得式中000001011001411iiiiiiiiiNbbNaaabababNNabB (i=1, 2, 3,
16、4) (6.39) 可以得出用节点位移表示的单元应力，即1234e= DSSSS (6.40)式中iiS=D B (i=1,2,3,4) (6.41) 对于平面应力问题0000200111141111122iiiiiiibaEbaabab S (6.42)将单元刚度矩阵写成分块形式11121314212223243132333441424344ekkkkkkkkkkkkkkkkk (6.43)其中的子矩阵按下式进行计算TijijtdxdykB DB (6.44)设单元厚度t是常量。上式可以积分出如下显式表达式1111211111132324 11111112323Tijijijijijijij
17、ijijijijijijijtabd dbaabEtabba kB S (i, j =1,2,3,4) (6.45)对于平面应变问题，将上式中的E、分别换成E / 1- * 和 / 1-。eReR(3) 载荷列阵和整体结构有限元方程矩形单元的节点位移与单元节点载荷矩阵 之间的关系为eeek R (6.46)由于矩形单元有四个节点，所以单元载荷矩阵 有8个元素，即11223344eTxyxyxyxyFFFFFFFFR (6.47)与平面三角形单元相同，将各单元的、 e 和和Re 都扩充到ek维数,再进行叠加,可得到整个结构的平衡方程 (6.48)即， Kd=R整体结构自由度的(4) 矩形单元的特
18、性 矩形单元的位移模式比平面三角形单元的线性位移模式增添了项（相当于xy项），这种位移模式称为双线性模式。在这种模式下，单元内的应变分量不是常量，这一点可以从应变矩阵B的表达式中看出。 矩形单元位移模式中的1、2、3、5、6、7与平面三角形单元相同，反映了刚体位移和常应变。位移模式在单元的边界上( =1或 =1)，位移是按线性变化的，显然，在两个相邻单元的公共边界上，其位移是连续的。 但是，矩形单元存在一些明显的缺点：其一是矩形单元不能适应斜交的边界和曲线边界，其二是不便于对不同部位采用不同大小的单元。 20节点三维等参元的母单元是边长为2的20节点正方体单元，对应的是边界为曲面和曲边的六面体
19、子单元，如图所示。每个节点具有三个平动自由度， 图6-5 20节点空间等参元 根据等参元的概念，位移函数和几何坐标变换式应采用相同的形函数。20节点三维等参元的位移函数可表示为uNuvNvwNwiiiiiiiii ,120120120 (6.49)对于单元的20个节点分别写出20个形函数，具体表达式如下Niiiiiiiiiiiii1118111141111411114000000222200222200222200222 / (6.50)式中, iii000,节点5的坐标是（-1，-1，1）等。. i，i及i是节点i在局部坐标系中的坐标。例如，节点1的局部坐标是（-1，-1，-1），根据弹性力
20、学几何方程，可以得到单元应变列阵为 12122020 xyezxyyzzxuxvywzBBBuvyxvwzywuxzB (6.52)xNxyNyzNziiiiiiiii ,120120120 (6.51) 坐标变换关系可表示为000000;1, 2 , 2 0000iiiiiiiiiiNxNyNziNNyxNNzyNNzx B (6.53)其中，B是单元的应变矩阵，其分块形式如下 上式中的形函数Ni是局部坐标的函数。在对整体坐标求导时，类似于平面问题，根据复合函数求导数的规则，有以下用雅可比矩阵表达的关系式iiiiiiNNxNNyNNzJ1iiiiiiNNxNNyNNzJ (6.54)其中，三
21、维雅可比矩阵J为xyzxyzxyzJ (6.55)式中每个元素可以分别按如下公式求得202020111, , ., iiiiiiiiiNNNxyzxyz (6.56) 将单元应变代入空间问题的物理方程式，得到单元的应力1220TxyzxyyzzxeeeSSSDDBS (6.57)式中， S为应力矩阵 iiSD B（i=1, 2, , 20） (6.58)进而利用虚功原理，推导得到单元刚度矩阵为1,11,21,202,12,22,2020,120,220,20eTkkkkkkdxdydzkkkKB DB (6.59) 111,1111,2,20;1,2,20eTTi jijijdxdydzJ d
22、 d dij KBDBBDB (6.60)其中 在等参元的单元刚度积分过程中，矩阵中的每个元素都很复杂，必须采用数值积分方法加以完成，如高斯积分法。在这里对高斯积分法的基本原理进行介绍。 高斯积分法是计算复杂的定积分时通常采用的一种数值方法。对于一维定积分问题fd11 (6.61)可以近似的化为加权为题。所谓数值积分就是把定积分问题近似地化为加权求和问题，在积分区间选定某些点(称为积分点)，求出积分点处的函数值，然后再乘上与这些积分点相对应的求积系数（又称加权系数），再求和，所得的结果认为是被积函数的近似积分值。这种求积方法表达如下 niiifHdf111 (6.62)高斯积分法采用以上这种格
23、式，其中积分点坐标i及其对应的加权可以构造出二维和系数Hi如表6-1所示。逐次利用一维高斯求积公式三维高斯求积公式 高斯积分的阶数通常根据等参元的维数和节点数来选取。例如，平面4节点等参元可取2阶，平面8节点等参元可取3阶，空间8节点等参元可取2阶，而空间20节点等参元可取3阶。表6-1高斯积分法中的积分点坐标和加权系数积分点坐标i积分点数n加权系数Hi20.57735031.000000031.00000000.77459671.88888890.555555640.86113630.33998101.34785480.65214525 0.00000000.90617980.5384693
24、1.56888890.23692690.4786287fd dH H fijijjmin ,111111fd d dHHH fijkijkkljmin , , , 111111111 (6.63) (6.64)0.323000KN/mw现有一个受均匀分布载荷作用的薄板结构，几何尺寸及受力情况如图6-6(a) 所示，本节将采用4节点矩形等参元对其求解，读者应掌握具体的求解过程。设该结构的弹性模量E=210Gpa，泊松比板厚度t=0.025m，均布载荷试对该结构进行静力学分析。(a) (b) 图6-6 薄板结构力学模型及有限元模型 为了方便说明，我们将平板仅仅分解为两个单元，6个节点，如图6-6(
25、b)所示。分布载荷的总作用力平均分给节点3和节点6。每个节点力为求解过程如下:(1)对结构进行离散化3000 0.025 0.259.3752FKN节点坐标见表6-2 表6-2 节点的坐标节点坐标值123456x00.250.500.250.5y0000.250.250.25（2）求解单元的刚度矩阵 四节点四边形单元的单元刚度矩阵的表达式为1111d deTt kB DB J其中11122122JJJJJ而 )1 ()1 ()1 ()1 (41432111xxxxxJ)1 ()1 ()1 ()1 (41432112yyyyyJ)1 ()1 ()1 ()1 (41432121xxxxxJ)1 (
26、)1 ()1 ()1 (41432122yyyyyJ12341BBBBBJ应变矩阵中的每一项为 2212112111212212001,2,3,4iiiiiiiiiNNJJNNJJiNNNNJJJJB弹性矩阵D D为2101011002ED由此可以求得两个单元的刚度矩阵，分别为 2.5962 0.9375 -1.5865 -0.0721 -1.2981 -0.9375 0.2885 0.0721 0.9375 2.5962 0.0721 0.2885 -0.9375 -1.2981 -0.0721 -1.5865 K1 = 1.0e+006 * -1.5865 0.0721 2.5962 -0
27、.9375 0.2885 -0.0721 -1.2981 0.9375-0.0721 0.2885 -0.9375 2.5962 0.0721 -1.5865 0.93750 -1.2981-1.2981 -0.9375 0.2885 0.0721 2.5962 0.9375 -1.5865 -0.07211-0.9375 -1.2981 -0.0721 -1.5865 0.9375 2.5962 0.0721 0.28850.2885 -0.0721 -1.2981 0.9375 -1.5865 0.0721 2.5962 -0.9375 0.0721 -1.5865 0.9375 -1.2
28、981 -0.07218 0.2885 -0.9375 2.5962K2= 1.0e+006 *2.5962 0.9375 -1.5865 -0.0721 -1.2981 -0.9375 0.2885 0.0721 0.9375 2.5962 0.0721 0.2885 -0.9375 -1.2981 -0.0721 -1.5865 -1.5865 0.0721 2.5962 -0.9375 0.2885 -0.0721 -1.2981 0.9375-0.0721 0.2885 -0.9375 2.5962 0.0721 -1.5865 0.93750 -1.2981-1.2981 -0.93
29、75 0.2885 0.0721 2.5962 0.9375 -1.5865 -0.07211-0.9375 -1.2981 -0.0721 -1.5865 0.9375 2.5962 0.0721 0.28850.2885 -0.0721 -1.2981 0.9375 -1.5865 0.0721 2.5962 -0.9375 0.0721 -1.5865 0.9375 -1.2981 -0.07218 0.2885 -0.9375 2.5962（3）组集整体刚度矩阵利用直接组集法，可以将上述单元组集成整体刚度矩阵，由于共有6个节点，整体刚度矩阵是一 个12X12的方阵，具体值为 K=1.0
30、e+006 * 2.5962 0.9375 -1.5865 -0.0721 0 0 0.2885 0.9375 2.5962 0.0721 0.2885 0 0 -0.0721 -1.5865 0.0721 5.1923 0 -1.5865 -0.0721 -1.2981 -0.0721 0.2885 0 5.1923 0.0721 0.2885 0.9375 0 0 -1.5865 0.0721 2.5962 -0.9375 0 0 0 -0.0721 0.2885 -0.9375 2.5962 0Columns 1 through 70.2885 -0.0721 -1.2981 0.937
31、5 0 0 2.5962 0.0721 -1.5865 0.9375 -1.2981 0 0 -0.9375 -1.2981 -0.9375 0.5769 0 -1.2981 0.9375 -1.5865 -0.9375 -1.2981 0 -3.1731 0.9375 -1.2981 0.0721 0 0 -1.2981 -0.9375 0.2885 0.0721 0 0 0 -0.9375 -1.2981 -0.0721 -1.5865 0Columns 8 through 12 0.0721 -1.2981 -0.9375 0 0 -1.5865 -0.9375 -1.2981 0 0
32、0.9375 0.5769 0 -1.2981 -0.9375 -1.2981 0 -3.1731 -0.9375 -1.2981 0 -1.2981 0.9375 0.2885 -0.0721 0 0.9375 -1.2981 0.0721 -1.5865-0.9375 -1.5865 0.0721 0 0 2.5962 -0.0721 0.2885 0 0 -0.0721 5.1923 0 -1.5865 0.0721 0.2885 0 5.1923 -0.0721 0.2885 0 -1.5865 -0.0721 2.5962 0.9375 0 0.0721 0.2885 0.9375
33、2.5962 （4）引入边界条件 节点位移列阵为 112233445566TxyxyxyxyxyxyU U U U U U U U U U U UU节点力列阵为112233445566TxyxyxyxyxyxyF F FFF FFFFFFFF 由图6-6（b）可知 114422553366009.375,0,9.375,0 xyxyxyxyxyxyUUUUFFFFFFFF静力学的求解方程为 12 12KU F 将边界条件代入上述方程，即引入了位移边界条件和力边界条件。因为1节点和4节点位移为0，可以将方程6.66进行化简，提取整体刚度矩阵是的第3-6行、第9-12行的第3-6列、第9-12列作为子矩阵，得到新的求解方程8 8KUF22335566TxyxyxyxyU U U U U U U UU
