有限元法讲义

《有限元法讲义》由会员分享，可在线阅读，更多相关《有限元法讲义（34页珍藏版）》请在装配图网上搜索。

1、有限单元法LAE图 3-1 简单的梁和桁架结构在这种情况下，对于梁单元、弹簧单元和桁架单元，我现考虑对图 3-1中结构的分析。在位移分析法中我们是把该结构看成是两个 梁单元、一个桁架单元和一个弹簧单元的分割体，第一步是计算对应于结构总体 自由度的单元刚度矩阵 们分别有K = El1L-6_z2 6 - L412_一 6 - Lu 一 Q6一 L 412一L2咖； U1,U2 ,U3 ,U4EzT一一e212K121212L 812一L1626 一 Lu-L2L67U,5U整个分割体的刚度矩阵可以由各个单元刚度矩阵通过直接刚度法有效地求 得。在这个过程中，结构刚度矩阵K是通过各单元刚度矩阵直接相

2、加而算得， 即K =工 Keii 而系统的平衡方程为KU = R式中U是系统的总体位移向量，R是作用在结构总体位移方向上的外力向量：UT =U U , Rt =Ir R R 1 7 1 2 7在求解结构的位移之前，我们需要利用边界条件U = 0和U = 0。这意味17着我们可以只考虑含有五个未知位移的五个方程，即KU = R，式中K是从K中删去第一和第七行以及第一和第七列后得到的，而Ut =U U U U U , Rt =10 一P 0 0 023456上面的讨论说明了有限元法的一些重要特点基本的处理过程是先把整个结 构看成为各个结构单元的分割体，计算对应于结构分割体总体自由度的各单元刚 度矩

3、阵，然后通过将各单元刚度矩阵叠加的方法形成结构刚度矩阵，求解单元分 割体的平衡方程组就得到单元的位移，然后利用它们来计算单元的应力。虽然原 来并不认为用位移法分析梁和衍架单元的分割体就是有限元分析，但以后我们将 会看到，实际上我们可以把桁架和梁单元看作两种特殊的有限元。在这个定义上， 上述的分析就是梁和衍架结构的有限元分析。而用另一种方法，我们可不通过求 解平衡微分方程，而是用虚功原理计算其刚度系数。导出表示弹性体平衡的相应方程的一个等效方法是利用虚位移原理，这个原理表示物体处于平衡的要求是：对于强加在该物体上的任意相容的微小的虚位移，总的内虚功应等于总的外虚功，即Jg T Q. dV =J

4、UTfBdV + J USTfsdS +V式中为 UiTFii3.1)fBfsFiXXXfB，fS =fs， F i =FiYYYfBfsFiZZZSVfB3.2)是作用在弹性体上的体力fB、表面力fs和集中力Fi。从未受荷载时的位置开 始的弹性体的位移以U表示，其中3.3)Ut =U V W相应于 U 的应变为8 T =18XX相应的应力为CT T = tXXYYYY8ZZZZY XYTXYY YZTYZY ZXT ZX3.4)3.5)亍和D表示虚应变和虚位移。分析的基本步骤和上述桁架和梁结构的分析步骤一样。图 3-2有限元平面分析254该问题的求解可按下列的步骤进行：假设每一单元节点i的两

5、个未知位移为u和八而u(X,Y)和v(X,Y)用 ii简单多项式函数来表示，其中的末知参数是单元节点位移。对于图 3-3中的单元未知位移是 u ,v , ,u ,v 。1 1 4 4(2) 利用虚位移原理计算每个单元对应于节点自由度的刚度矩阵。(3) 将各单元刚度矩阵叠加得到结构的刚度矩阵，利用边界条件解平衡方程 组求出节点位移，然后算出单元的应力。图 3-3 局部坐标系统中典型的二维四节点有限元3.1 利用虚位移原理建立有限元法的公式(一)平面应力分析的位移和应变位移的变换矩阵 为了便于说明，再考虑一个平面内荷载作用下悬臂板分析的例子(图 3-2)。该结构是处于平面应力的状态，因此可以将平桓

6、理想化为二维平面应力有限元的 分割体，如图 3-2 所示。所需要的基本矩阵是对于分割体的每个单元的位移变换矩阵和应变-位移变 换矩阵。为了推导这些矩阵，我们着重研究如图 3-3所示的典型单元，并假设局 部单元位移u和v是以局部坐标变量x和y的多项式的形式给出：u( x,y ) = a+a xy + a xy(3.6)v(x,y) = P +P x + P y + P xy(3.7)1 2 3 4未知系数a ,P也称为广义坐标，它们可以用未知的单元节点位移u，u和1 4 1 4匕来表示。定义u T = u u12我们可以把式(3.6)和uuvvv3 4 1 2 3( 3.7)写成矩阵形式vl4u

7、( x,y) = Ga其中G=e0eu(x,y)=u(x,y)v(x,y)aT = Exa a a1234x y xyP1 P2 P3Pl43.8)( 3.9)( 3.10)( 3.11)3.12)对于单元的各个节点，方程(3.9)都一定是成立的。因此，对于所有四个节点，利用式(3.9)，我们有u = Aa其中A =人 00 A3.13)1xyx y111 11xyx yA =222 21xyx y3333_1x4y4x y443.14)因而广义坐标 a 是a = A -i 五3.15)其中A-i0 _A-i =_ 0A -i _将式(3.9)代入，可以得到u(x,y) = N(x,y)u其中

8、N( x,y ) = A-i(3.16)(3.17)3.18)我们用虚功原理来得到单元刚度矩阵。若是平面应力状态，单元应变为8 T = 8XX8YYY XY其中du8 = XX dx8YYdvdydvdu= + -丫 XXdx利用式3.17），我们得8 = Bu式中B= EA-1010y0000E=0000001x_001x010y单元应力是CT T =IXXCTYYTXY3.19)3.20)(3.21)3.22)(3.23)(3.24)假定是各向同性的线弹性材料，对于平面应力状态，有CT = D8其中g 0100巳2D = E1- g23.25)3.26)而E和卩是材料的弹性常数，E为杨氏模

9、量，p为泊松比。要进行实际的有限元分析，就需要求出式(3.18)和式(3.22)中分别给出的分 割体中每个单元的位移和应变-位移的变换矩阵，而算出每个单元的矩阵A，才 能求出这些矩阵。至今，我们所考虑的是通过单元局部节点位移来决定的每个有限元的位移。 然而，在推导整个单元分割体的平衡方程的时候，通过整个单元分割体节点位移 来表示单元的位移是方便的，即对单元m可写出u(m)( x,y ) = N(m)( x,y )U(3.27)而对于单元应变和应力类似地可写成(m)(x,y) = B(m)(x,y)U(3.28)(m)( x,y ) = D(m) (m)( x,y )(3.29)式(3.27)至

10、式(3.29)中的量与向量U有关，而U存储整个有限元分割体的总体 坐标系下的全部节点位移。作为一个例子，考虑N(2),即确定图3.2中悬臂板 理想化后的单元2 上位移矩阵。利用图 3.2 给出的整个单元分割体和图 3.3 给出 的单元的节点位移的定义，我们有uvuv3322UVUVUVUV112233440NNNN00N (2) =1317121600NNNN0023272226uvuv4411UVUVUV V5566779NNNN00 0 14181115( 3.30)NNNN00 024282125式中 N 是式( 317)中考虑单元2 时对应的单元。 ij3.2 建立有限元步骤归纳有限元

11、法与结构力学中的位移法相似。首先将连续体转化为离散化结构，即 将连续体代之以仅在节点互相连结的许多单元组成的结构。这种离散化结构类似 于结构力学中的桁架或刚架，但其中的单元不一定是杆件，可以是平面块体或空 间块体等。然后对此离散化结构按类似结构力学的位移法进行分析，步骤如下：1. 取每个单元的节点位移u作为基本未知数。2. 在单元内建立位移模式，u = Nu。3. 根据几何关系建立应变矩阵，= BU。4. 根据物理方程建立应力矩阵，a = De = DBU = SU。其中S = DB称为 应力矩阵。5. 根据虚位移原理建立单元刚度矩阵，K(m)=! BTDBdV(m)oV(m)6. 建立单元节

12、点荷载向量，R(m) = R(m) + R(m) - R(m) + R(m)。BSIC7. 建立平衡方程，KU = R，其中K = ZK(m)，R = ZR(m)，U是全部mm节点位移向量。8. 代入约束条件求解平衡方程，K U = R。aa a a位移模式中的N称为形态函数。在有限单元法中，各种计算公式都依赖于位 移模式，位移模式选择的恰当与否，与有限单元法的计算精度和收敛性有关。为 了使位移模式尽可能地反映物体中的真实位移形态，它应满足下列条件：(1) 位移模式必须能反映单元的刚体位移。(2) 位移模式必须能反映单元的常量应变。(3) 位移模式应尽可能地反映位移的连续性。 在单元之间，除了

13、节点处有共同的节点位移值外，还应尽可能反映在单元之间边界上位移的连续性。第一章、 杆件结构的有限单元法4.1 平面刚架结构平面刚架结构的有限单元法分析，通常采用两端固接的平面固结单元 (见图4-1)，现针对这种单元讨论有限元分析计算公式的建立。(a)坐标系，单元节点位移与节点力向量杆件结构的有限元分析需要建立单元局部坐标系与结构整体坐标系。对平面 固结单元(如图 4-1 所示)，结构矩阵位移法取结构的节点位移为基本未知量。单 元节点位移与基本未知量之间满足变形协调条件，在局部坐标和整体坐标下单元 节点位移、节点力分别记为u fT = Lv e u ve】(4.1)i ri i j jjUT =

14、 Uveuve】(4.2)i ri i j jjR fT = N /V9 M NV9 M】(4.3)i iijjjRt = IniVMNiijVMjj】(4.4)图4-1中各量值所示方向均为正方向。设单元的弹性模量为E,惯性矩为I 截面积为A,杆长为L。(b)位移模式，应力与应变矩阵 由材料力学知识可知,位移模式可以选择u(x) = a + a x(4.5)01v(x) = b + b x + b x2 + b x3(4.6)0123注意e=-dVdx可得单元的位移模式u(x) u =v(x)4.7)其中，形函数0X102x2x33 x 22x3一 X +0l12121303x 22x31 一

15、 + l2l30X 2X 3T -12单兀的线应变分为拉压应变和弯曲应变E两部分(略去剪切变形)，于是有其中=B=0b1 l0对应的应力为dudXd2v-ydx 24.8)6 y 12Xyl 2l 3一 41+6xyT IT一 6 y +12 xy一 12l3一2y+6xyT l2CT0CTb(c)单元刚度矩阵根据虚位移原理可知，在局部坐标系中CT =E E:卜 DB4.9)K(m) = f BTDBdV (m)4.10)V(m)具体表达式为EAT012EI136EI12lEA00l12EI6EI012EI1312136EI2EI06EI4EI12l12lsym4EI(d)坐标转换矩阵，整体坐

16、标下单元刚度矩阵EAT0从图4-1可知，局部坐标系相对于整体坐标系逆时针旋转了 8角。根据单元的节点位移或节点力在两个坐标系中的投影关系，可得坐标转换关系cosO sinB 0一 sinO cosO 04.11)_ 0 0 1_于是，单元节点位移或节点力的坐标转换矩阵(从整体坐标转到局部坐标)为O=4.12)根据正交矩阵的特性和矩阵的相关知识，可知在整体坐标中的单元刚度矩阵K(m) =O tK (m)O4.13)(e)节点荷载向量单元上的荷载按虚位移原理移置到节点上、得到等效节点荷载向量4.14)其中G是作用在杆上的均布荷载，Fi是集中荷载。整体坐标中的节点荷载同局部坐标中的关系为4.15)(

17、f)结构的平衡方程根据虚位移原理，结构的平衡方程表示为KU = R4.16)其中K = ZK(m)，U = Xu，R = ZR分别代表整体总刚度矩阵，总节点位移向量和总节点力向量。4.2 平面桁架结构平面桁架结构的有限元分析通常采用两端铰接的平面铰结单元，其坐标系如 图 4-2 所示。(a)坐标系，单元节点位移与节点力向量 杆件结构的有限元分析需要建立单元局部坐标系与结构整体坐标系。对平面 固结单元(如图 4-2 所示)，在局部坐标和整体坐标下单元节点位移、节点力分别 记为u T = U viiUT = U Viiv 1jV1(4.17)4.18)R fT = v QiiRt = P Pxi

18、yi图 4-2 平面桁架单元N j QjjjPP_ EAl0EAl0K(m)=0EA00EA0-l0l0_ 0000_xj yj(b)局部坐标中的单元刚度矩阵(c)坐标转换矩阵cos0sin0000=-sin0cos00000cos0sin0_ 00- sin 0cos0(d)整体坐标中的单元刚度矩阵K(m) = 0 tK(m)Q =C 2cEA C C c sl - C 2c-C Ccs其中 C = cos0 , C = sin0。cs(4.18)4.19)4.20)symC2s-C C C2c sc-C2 C C C2sc s s4.21)(4.22)第二章、 平面问题的有限单元法连续体的

19、离散化与三角形常应变单元(a)三角形单元的节点编码| y4/A/d.9 7/2J51图 5-1 平面应力问题我们用一组网格线把平板划分成若干三角形区域，每个三角形就算作是一个 单元，网格线的交点就算是节点。把这些单元和节点按一定的顺序编号，并把各 个节点沿坐标轴的位移取为基本未知量：对于节点1对于节点2对于节点iU T = U11U T =U23UTi=UU4(2i-1)U2i设有n个节点，则整体的总节点位移向量为Ut = U U U U U 12( 2i-1)2i2n对每个单元来说，它与三个节点发生联系，设其整体序号是i,j,k，依次标记 为 1,2 ,3 ， 则 此 单 元 的 节 点 位

20、 移 向 量 可 以 表 达 为 ( 参 看 图 5-2 )图 5-2 三角形常应变单元ut = u v u v u vU U1 1 2 2 3 3( 2i-1) 2i ( 2 j-1)单元的节点位移向量U是整体的总节点位移向量的一部分。U U U 2 j( 2k-1 )2k(b)位移模式和桁架不同，这里没有现成的用节点位移向量来表示单元内任一点位移的现 成公式，因而必须采用近似假设。我们知道，任一函数常常可以用有限项的幂级 数来迫近。所以不妨假设单元上任何一点的位移为u(x,y)=u(x,y) a +a x + a y畑3+a 4x+a：y5.1)也就是把高阶的项全都略去，只保留线性项，这自

21、然会带来截尾误差。从几何意义上来解释：函数u(x,y)可认为是定义在单元上的某种复杂曲面，而现在用一平面来代表它，这样做当然是有误差的。可是我们从直觉上可以想像得出来，若 把单元划分得很小，那末在足够小的范围内平面和曲面是相差不大的。所以可以 用式(5.1)来代表真正的位移函数而不致造成很大的误差。现在要用节点位移来表示式(5.1)中的系数a ,a a，将节点坐标代入 015该式，位移应为节点位移，即u1xya 1110u=1xya2221u1xya313325.2)由上式可以得到用节点位移表示待定系数的式子，之后再代回到式(5.1)之第式，就有U(X，y) = N1U1+ N 2 U2 +

22、N 3 U3其中5.3)1N = (a + b x + c y)1 2 A 1 1 15.4)a = x y 一x y ， b = y 一y ， c = -x + x (1 t 2 t 3 t 1) (5.5)12332123123x1x2x3y1y2y3为单元的面积。其他 N ,N 的式子可用下述方法得到：将上边四个式子小的下标按“循环交换”的规则依次更改。所谓循环交换就是（1 T 2 T 3 T 1）的交换方式。5.6)同理可得V(Xy) = N 儿 + N 2v 2 + N 3v 3综合起来，就是u(x,y)=u(x,y) _v(x,y)_=Nu其中_ N0N0N0N =123_ 0N0

23、1N20N5.7)5.8)通过以上处理方法，连续体上任一点处的位移皆可用若干选定的节点位移来 表示。一个寻求未知函数的问题已经转化成一个寻求若干个未知数的问题。这样， 我们就完成了将连续体离散化的关健一步。由此可见形状函数在有限单元法中占 据的重要地位。式（5.4）和式（5.6）表明，这种单元的形状函数矩阵是由三个节点处的形 状函数组合而成的。每个形状函数都是线性函数，在几何上代表一个平面。在所代表的节点处，函数的值为单位1，在另外两个节点处，函数的值为零。图 5-3是这类函数的图形表示（以N 为例）。图 5-3 形状函数平面（C）应力应变矩阵利用前面的知识，可知应变8x8y1xy其中应变矩阵

24、B = B1d udXd v荡d ud v+ - dydxckbk=Bu；k = 1,2,3(5.9)5.10)5.11)可见由于采用了线性的形状函数，单元上的应变(以及应力)为常数。所以这种单 元也可称之为常应变(或常应力)单元。这当然是近似的，这种应力实际上只是真实应力的某种平均值。对于平面应力问题，应力为b T = L Q T = D8(5.12)x y xyg 01 00 口2其中g为材料的泊松比，E为材料的弹性模量。对于平面应变问题，只要将弹性矩阵D中的E用需、g用佥代替即可。将应变代入，可得q = DBU = SU其中 S 称为应力矩阵(针对平面应力情形)S3S = DB = S1

25、Sk 2( 1 g 2 )Abk1gbgk1g c2kgckck巳b2k；k = 1,2,3(5.13)5.14)(5.15)(d)单元刚度矩阵根据虚位移原理，由式(3.39)可知，单元的刚度矩阵为K(m) = J BTDBdV(m)V(m)当平面问题的板厚度为h时，上式改写为K(m)= hJ BTDBdA(m)A( m )5.16)此处的标记(m)均表示第m个单元。将式(5.9)到式(5.15)代入上式，可得Ksym11K21KK22KKK ( m ) =313233式中每个子矩阵都是2x 2阶的，其下标与单元的节点编号相对应，具体的计算公式如下：Ehy 4( 1 - g 2 )Ab b +

26、 _ c cgb c + _ bei j 2 i ji j 2 j i1- 1-b c + b c c c + b bj i 2 i j i j 2 i j；i,j = 1,2,3(e)节点荷载向量对于连续体来说，在单元的边界上有着分布作用的边界力，它们是其他单元 或是结构的外部环境对该单元的作用力。而在节点上，一般却没有集中的节点力， 这和杆系结构不一样。不过出于习慨在物理意义上我们常常可以设想把这些边界 力向节点处简化，形成一组等价的虚拟的节点力。显然当单元尺寸很小时用这种 等价节点力代表真实的边界力不致引起很大的误差。由第三章的知识，我们知道三节点三角形单元的节点荷载向量表示为Rt =R

27、 R R =yY X Y1 2 3 1 1 2 2X Y33面给出几种情况下的节点荷载向量的表达式：(1)设在三角形单元的(兀,y)点受到集中荷载FtRt = Lf NF NF N F N F N F1 x 1 y 2 x 2 y 3 x 3 y=Fx则5.17a)(2)设单元受y方向的体力作用，例如受重力作用。其合力W =-yAh，y为容重。Rt = -310 1 0 1 0 15.17b)(3)设单元的1 -2边上受有沿x方向的集中力F，其作用点距离节点1,2x的距离分别为 l ，l ， l l +l ，贝1212Rt F_ l20 li 00 0(5.17c)xll(4) 设单元的1 -

28、 2边上受有沿x方向的三角形分布荷载,荷载在1点的集 度为q，在2点的集度为零，边长为l，贝URT = 2qlh3 0 1 0 0 0(5.17d)(5) 设单元的1 - 2边上受有均匀分布的侧压力q，则R = 2qhh(y! - y2 )(X 一 X! ) G 一 y ) (X2 - X! ) 0 01 U5.1 面积坐标在采用较精密的三角形单元时，利用面积坐标，可以大大简化荷载向量、应 力矩阵和刚度矩阵矩阵的推导工作。在如图 5-4 所示的三角形单元中，任意一点 P 的位置，可以用如下的三个比值来确定X图 5-4 面积坐标L i ,i AJmmA5.18)其中A为三角形ijm的面积，A、A

29、、A分别为三角形Pjm、Pmi、Pij的面i j m 积。这三个比值就称为 P 点的面积坐标。注意，三个面积坐标并不是互相独立 的，即L + L + L 1(5.19)i j m 这里所引用的面积坐标，只限于用在一个三角形单元之内，在该三角形之外并没有定义，因而是一种所谓局部坐标。与此相反，以前所用的直角坐标X和y ,则是一种所谓整体坐标，它通用于所有单元的总体也就是通用于整个结构物。根据面积坐标的定义，在图5.4中不难看出，在平行于jm边的一根直线上的所有各点，都具有相同的L坐标而且这个坐标就等于“该直线至jm边的距离” i与“节点i至jm边的距离”的比值。图中示出L的一些等值线。 i根据面

30、积的计算公式，可知1 XA = 1 Xij1 Xmyyjjym(a + b X + c y),2 i ii(i,j,m )5.20)同式（5.4）相比，可知三角形单元的形状函数N，就是对应的L。用矩阵表示 ii为Li1aibici1LabcX( 5.21)j2 AjjjLmambmcm_ y _将上三式分别乘以 X ,X ,X 然后相加，可得ijmx L + x L + x L = X i i j j m m 同样，用 y 坐标，可得yL +yL +y L =yi i j j m m5.22)二式表明了直角坐标与面积坐标之间的关系，写成矩阵形式_ 1_ 111X=XXXijm_ y _y1-

31、IyjjymLi Lj L m对面积坐标函数求导数时，可以利用下面的公式ddxd1 b I2A cL ibmcmdLddL5.23)dLm求面积坐标的幂函数在三角形单元上的积分值时，可应用积分公式D LflLbLc dxdy =i j mAa!b!c!(a + b + c + 2)!2A5.24)求面积坐标的幂函数在三角形莱一边上的积分值时，可应用积分公式f LaLbds =ijl5.25)其中l为ij边的长度。(a)位移模式假定位移模式按照二阶多项式变化，即u(x,y)=u(x,y) I a +a x + a y +a x2 +a xy + a y2I = I 012345v(x,y) I

32、I a +a x + a y + a x2 +a x +a y2 6 7 8 9 10 115.26)图 5-5 六节点三角形单元上式中包含着十二个待定系数，为了把它们用节点位移表示出来，每个单元应当 有六个节点。通常的作法是在三个顶点之外再取三个边的个点，如图 5-5所示。对应的形状函数矩阵是N = BNNNNN234565.27)其中N = Nf ，并且iI0N IIi1n. =(g.1 + g 2x + g 3y + g 4x2 + g 5xy + g 6y2)i2A2i1i 2i3i4i5i6式中六个系数与三角形的几何尺寸有关i1=a (a 一 A);g = b (2a - A)i i

33、i 2i ig = c (2a 一 A);g = b2i 3 i ii 4ig = 2b c ; g = c 2i 5i i i 6i(i = 1,2,3 )、 g= 2a a;g =2(a b一 ab)i1i - 3 i - 2 i 2i - 3 i - 2i - 2 i - 3g= 2( a c一 ac );g= 2b bi3i-3i-2i-2i-3i4i-3i-2g= 2( b c一 bc );g= 2c ci5i-3i-2i -2i-3i 6i-3i-2(i =4,5,6)式中a,b,c的下标仍然按照(1T 2 t 3 T1)循环。(b)应变矩阵B = BBBBBB1234561式中B

34、 - 1g0 2 g0 g0 i20gx+i4LL 0g+ yi5LL 02 gi 2A 2i 32A2i 52A2i 6L g 3g 2L g 52 gi42 gi6g 5对应的应变矩阵为(5.29)可见应变是线性变化的，并可以表示为B = B + B x + B y。ii 1ix iy(c)单元刚度矩阵同理，单元刚度矩阵表达为K (m ) =SKY22KK3233KKK424344KKKK52535455KKKK62636465112131415161式中任意一个子矩阵可以写成K =1ijV(M )i1ix=刀工Krq jr=1,x,y q=1,x,yBT DB dV (m )V(m )

35、i j(B + B x+B y)T D(B V(M ) i1ix iyj166+B x+ B y)dV(m) jx jyijrqdV (m )V(m)5.30)其中Krq =BTDB(r =1,x,y;q=1,x,y)ij ir jq此外，计算时应注意到，J rqdV (m)分别为三角形单元的体积、静矩、惯性矩V(m)和惯性积等等。5.2 矩形单元三角形单元的优点之一是它 的“适应性”， 任何复杂边界的 弹性体，总是可以划分成三角 形。可是在规则边界的情况下， 显然划分成矩形更加方便。另 外，许多计算实例已经证明，矩 形单元的计算精度也比三角形 单元好。这是因为在每个小矩形 范围内，矩形单元有

36、连续变化的 应力场，而对应的两个三角形单元却只有阶梯变化的应力场。所以，在复杂边界 的情况下，同时使用三角形和矩形单元将是可取的计算方案。与三角形单元不问，不可能采用完全的多项式作为位移函数。因为一次的完 全多项式只有3个待定系数而矩形单元却有4个角点。为使二者的数目一致应在 二次的多项式中选取补充项。设我们补充兀2的项，则在矩形的上、下边界上位 移曲线是二次的，它不可能由两角点处的位移来唯一地确定，因而是不相容的。同样的道理，也不可以取y 2的项。所以，应取u( x,y ) = c + c x + c y + c xy1 2 3 4( 5.41 )v(x, y)=c +c x+c y+c x

37、y这个函数虽然含有二次项，但在单元的任何一条边界上都是线性的，所以只 要两点就可以唯一的确定位移，因而满足相容条件。根据上式可以写出：N=(1 - 2a)( 1 -缶；N=X( 1 - 2b)； n=x2b ； n=(i -壮若改用自然坐标表示(见图5-8), g= Xi a，n = y； b可得11N = 4 (1 弋)(1 -n)； N2 = 4(1Y)( 1 -n )11n = -（i+g）（i+n）； n = -（i）（i+n）3 4 4 4 或者概括成1N=4（1+软）（1+灯）；i=，45.42)单元的应变为& = Bu 式中UT =U1B = B1UU34B B 34UUU567

38、U 为单元节点位移向量。8g (1+nnii0n( 1+gg )bii对于平面应力问题，n (1+gg)b ；lIg( 1+nn )aii单元的刚度矩阵为iiiii = 1,2,3,45.43)K(m)=K11K21K31K41sym42K33K K4344式中K =!ijBT DB dV (m ) V(m ) i jEh12( 1 -m 2)(3+nn )卩 +茫料n (3+gg )ji j2p i j j33就n +(1 -m尤nj i 2i j33n4n +t/1-m 兀 n(i,j =1,2,3,4)(5.44)j 2jnn (3+gg )(3+nn )pj i j2 i j1 j其中

39、p = ba , h为单元厚度。对于平面应变问题按照前述同样处理。从式（5.43）可以看出，在单元内部，应力（应变）不是常数。它虽然不是完全的一次多项式表达的，但却是一种线性变化的规律。6.1 等参单元等参单元是在1966年公开引入的（由Irons）,它们使得非矩形四边形单元 的存在成为可能。等参单元最明显的特点是单元可以有曲边，并且它们有特殊的 坐标系（称为自然坐标系,。等参单元在模拟具有曲线边界的结构和从粗网格到 细网格作的单元过渡时是很有用的。等参单元具有多方面的适用性；在二维和三 维弹性分析、壳体分析以及非结构的应用中已被证明是有效的。等参单元的构造 过程初看起来似乎生疏，但并不复杂。

40、一旦掌握了一种类型单元的构造过程，就 很容易应用到大多数其它类型的等参单元上。单元节点确定两件事情1节点位移向量（亦称节点自由度）表示单元内部的位移u，具有u = NU ；2节点坐标向量定义了单元内部任意一点的坐标（兀,（用向量x表 示），具有x = N。矩阵N和NN都是自然坐标的函数，如果上述两项中的节点集相同，且N和庐也相同，则这个单元是等参的。6.2 平面线性等参单元首先，需要再次说明，不同的等参单元终究是类似的。本节的列式可直接推 广到更复杂的单元：主要的变化是增加节点和利用不同的形函致。矩形单元，最大的不足在于其对边界的要求高，即只能是规则的、平行于整 体坐标轴的。如果能有任意四边形

41、的单元，这样的问题就解决了，但任意的四边 形单元的形函数的构造却不容易，等参单元能够很好的解决这样的困难。图 6-2 四边形等参单元映射关系图6-2中的坐标轴和n通过四边形对边中点，它们不需要正交，也不必平行于x或y轴。3坐标的方位由节点号确定。就是说，对于如下定义的形函数n,节点1在=n = -i处,节点2在 = -n = 1处等等。整体坐标和位移为;卜用；-其中敷=1X- y1y2y3x y 44u v 440 N0N34N0 N023有时将式（6.9）写成如下形式更为方便x =N x u = S Nu 亍 * i和 S y=SN y v=SN vi ii i其中，各个形函数为1n = -

42、（i+gg）（i+nn）； i = 1,2,3,4 i 4ii（6.9）6.10a）（6.10b）（6.10c）（6.11）（6.12）同矩形单元的完全一样。通常，u和v分别是平行于x和y方向的位移，在一般情形下它们不平行于g和n。一个特殊情形是其边平行于整体坐标xy的矩形，此时E和n变成无量 纲的形心坐标。下述推导类似上节的推导。为了求得单元刚度矩阵K,我们需要B矩阵，我 们不能用x和y来表示B矩阵，而用g和n来写出它，这就要求引用下列导数的 坐标变换。设e是某一x和y的函数。那么由复合求导法则得到wede=dx 鬲 dxdy _空dx de或者fe,w Ie,y6.13)6.14)其中J为

43、雅可比矩阵x，匸yJ11J 12xyJJJ1-刃2122仲,y从式(6.13)可以得到逆关系6.15)上述关系是一般的情况。对于本章的平面单元，e就是“或卩，雅可比阵j 的系数的数值与单元的尺寸、形状和单元所处的方位(节点坐标和节点位置编号)有关。对于目前讨论的 4 节点单元。从式(6.11)有J = x = N x + N x + N x + N x(6.16)11 上1,g 12,g 23,g 34,g 4其它的J ,J ,J具有类似的表达式。如果x = g和y，则J = J-1 = I。12 21 226.17)现在，可以计算B矩阵了。0J-1uu,nvvi理丿6.18)00NNi,nN

44、4N4 ,n0000N4 ,gN4 ,n6.19)u1000,xxu =0001VyvY0110,xxyv,yB 矩阵就是式(6.17)、(6.18)和(6.19)中三个矩阵的乘积。单元刚度矩阵为K(m) = J BTDBdV(m)V(m)=J1J1BTDBhJddn-1 -16.20)其中J是雅可比矩阵的行列式值，h为单元厚度。在平面问题中6.21)J = detJ = J J 一J J11 22 21 12雅可比行列式J是单元内位置的函数，并且等于从面积dd 变到dxdy的放大倍 数。矩阵 B 中的一个典型系数是与节点坐标有关的，并且在分子和分母上都是 变量gn的多项式。因此，式（6.20

45、）中的积分必须以数值方法求出。对于单元 荷载向量中的体力和初始应力部分（见式（3.36）的R ,R部分）都需要体积积bI分也必须采用数值积分（见下节）（NTf 一BtCt ）hJdgdn（6.22）一1 一ibi我们可以根据伴随矩阵求得式（6.18）中雅可比矩阵的逆矩阵J一 Jz2212（6.23）一 J J21 11式（6.20）和（6.22）中的积分微元的变换，是来源于如下关系。在直角坐 标中dV = dx x dy - dz = dxi x dyj - dzk = dxdydz图 6-3 坐标变换关系6.24)由于gn为非正交坐标系（见图6-3），所以dV = hdg x dn - k其

46、中dg = Xdg i + 彷 dg j “ 淹d i淹 dx労二 旳二亦dm 而dn* J代入到式（6.24），可得w=hdA=h I詈墙即g叽6.25)6.3 8 节点平面等参单元将单元的边节点加到线性单元上，可以得到等参族较高次的单元。图6-4 所 示8节点单元，边缘的形状和位移按照g 2或2变化。在二次单元中的形状和位移的插值多项式是以不完全三次多项式为基础的，但是，在这里我们不打算用广义位移替换节点位移的方法建立形函数N,而是用一种更直接、更有明显的物理意义的方法来构造形函数。我们来考察各种模式并把它们组合起来，在图6-5a中，节点位移按照g二次方、按照线性变化。除掉在节点5之处为1

47、外，在其它节点处皆为零，因而形函数N = u，类似的，5a形函数N = u。在图6-5c中，位移u按g和线性变化，而在节点5和8分别8b取12。这样在节点 1 将节点 5 和 8 的位移的一般减去，便可以得到节点 1 的形函数N = u1 (d )11u u u(c) 2 (a ) 2 (b)6.29)其中 u= (1g 2)(1 -n )2(a)u 二(i-g)( i-m)2(b)u = (i-g)(i-m) 4(c)表 6-2 平面二次等参单元的形函数4个线性边3, 2, 1或0个线性边，其余为二次边N包括14节点i加节点5加节点6加节点7加节点8N (i-g)(i-m)i4一 N /25一 N /28N (i+g)(i-m)24一 N /25一 N H6N (i+g)(i+m)34一 N /26-N /2N (i-g)(i-m)44-竹/2一 N /28N5(i-g 2)(i-m)2N6(i+g)(i-m 2)2N7(i-g 2)(i-m)2N8(i+g)(i-m 2)2