四节点矩形单元有限元分析

1、平面问题有限元分析平面问题有限元分析四节点矩形单元四节点矩形单元天津大学天津大学 建筑工程学院建筑工程学院 Tianjin University本节内容提要本节内容提要1 1、分析提高有限元法求解精度的、分析提高有限元法求解精度的途径途径2 2、简要回顾三节点三角形单元有、简要回顾三节点三角形单元有限元分析过程限元分析过程3 3、全面介绍四节点矩形单元有限、全面介绍四节点矩形单元有限元分析过程元分析过程4 4、总结、总结分析提高有限元求解精度的途径分析提高有限元求解精度的途径 Tianjin University一、三节点三角形单元的缺点l 三节点三角形单元精度低，收敛慢，由于单元内应力三节点

2、三角形单元精度低，收敛慢，由于单元内应力和应变均为常量，故在单元内不能很好地反映应力和和应变均为常量，故在单元内不能很好地反映应力和应变的变化。应变的变化。l 该单元只有三个节点，单元自由度少，单元位移插值该单元只有三个节点，单元自由度少，单元位移插值函数（位移模式）只能是线性函数，描述单元内位移函数（位移模式）只能是线性函数，描述单元内位移变化的能力差。变化的能力差。分析提高有限元求解精度的途径分析提高有限元求解精度的途径 Tianjin University 二、提高有限元求解精度的途径l 第一个途径是对某一种特定类型的单元采用网格划分第一个途径是对某一种特定类型的单元采用网格划分加密，依

3、靠单元的收敛性提高求解精度。加密，依靠单元的收敛性提高求解精度。l 第二个途径是对一定的单元网格和单元尺寸，采用高第二个途径是对一定的单元网格和单元尺寸，采用高精度单元来提高求解精度。精度单元来提高求解精度。分析提高有限元求解精度的途径分析提高有限元求解精度的途径 Tianjin University 三、建立高精度单元的原理和途径l 原理：提高单元位移插值函数多项式的阶次，从而提原理：提高单元位移插值函数多项式的阶次，从而提单元拟合局部区域位移、应力变化的能力。单元拟合局部区域位移、应力变化的能力。l 途径：增加单元的节点数目。途径：增加单元的节点数目。l 对于平面有限元问题，除三节点三角形

4、单元外，还可对于平面有限元问题，除三节点三角形单元外，还可以考虑六节点三角形单元和四节点矩形单元。以考虑六节点三角形单元和四节点矩形单元。 Tianjin University三节点三角形单元有限元分析过程三节点三角形单元有限元分析过程设位移函数设位移函数求位移函数中的未知量求位移函数中的未知量代入函数中代入函数中整理可得形函数整理可得形函数()i j mN、 （性质？性质？）几何方程求解应变（几何矩阵几何方程求解应变（几何矩阵 ） B物理方程求解应力物理方程求解应力（弹性矩阵（弹性矩阵 应力矩阵应力矩阵 ） DS运用虚功原理求解运用虚功原理求解 eK 由由 合成合成 ( (方法？方法？) )

5、 eK K建立节点荷载列阵建立节点荷载列阵 （方法？方法？ F处理位移约束条件（处理位移约束条件（方法？方法？）组成？组成？) ) Tianjin University四节点矩形单元有限元分析过程四节点矩形单元有限元分析过程一、四节点矩形单元位移函数12345678uxyxyvxyxy单元节点编号为单元节点编号为 k，l，m，n（逆时针）（逆时针）单元节点位移列阵为：单元节点位移列阵为： Tnnmmllkkevuvuvuvu设位移函数为设位移函数为：( , ) ( , ) x yf x y或写为或写为： Tianjin University四节点矩形单元有限元分析过程四节点矩形单元有限元分析过

6、程二、求解位移函数中的未知系数将节点坐标将节点坐标 (,),( ,),( , ),(, )ababa ba b代入函数中，并写成矩阵形式：代入函数中，并写成矩阵形式： 123456781000000001100000000110000000011000000001kklelmmnnuababvababuababvababAuababababvababuababv 解上述方程组可得解上述方程组可得： 1eA Tianjin University四节点矩形单元有限元分析过程四节点矩形单元有限元分析过程三、将所求 值代入位移函数中 1( , ) ( , ) ex yf x yA Tianjin Un

7、iversity四节点矩形单元有限元分析过程四节点矩形单元有限元分析过程四、整理位移函数可得形函数( , )k l m nN 1( , ) ( , ) ( , ) eex yf x yAN x y展开上式可得：展开上式可得： 0000( ,)0000kkllklmnklmnmmnnuvuvNNNNux yNNNNvuvuv四节点矩形单元有限元分析过程四节点矩形单元有限元分析过程其中，形函数为：其中，形函数为：)1)(1 (41)1)(1 (41)1)(1 (41)1)(1 (41byaxNbyaxNbyaxNbyaxNnmlk Tianjin University四节点矩形单元有限元分析过程四

8、节点矩形单元有限元分析过程单元位移插值函数可以由单元位移插值函数可以由单元形状函数单元形状函数与与节点位移值的节点位移值的乘积表示：乘积表示：, , ( , )l m niii kx yN x y即可以表示为：即可以表示为：( , )( , )( , )( , )( , )( , )( , )( , )kkllmmnnkkllmmnnuN x y uN x y uNx y uN x y uvN x y vN x y vNx y vN x y v Tianjin University四节点矩形单元有限元分析过程四节点矩形单元有限元分析过程( , )( , )( , )( , )( , )( ,

9、)( , )( , )kkllmmnnkkllmmnnuN x y uN x y uNx y uN x y uvN x y vN x y vNx y vN x y v由此可见，位移插值函数完全由形函数决定；因此抛开节点位移，由此可见，位移插值函数完全由形函数决定；因此抛开节点位移，只讨论形函数的性质，就可以了解单元的变形性质只讨论形函数的性质，就可以了解单元的变形性质。例如：例如：四节点矩形单元，若四节点矩形单元，若则由则由 可得：可得：因此可以看出，单元变形完全由形函数决定。因此可以看出，单元变形完全由形函数决定。1,0,klmnklmnuuuuvvvv, ,( ,)l m niiikx y

10、Nx y, ,( , )l m neiikkki kuN x yN uNv Tianjin University四节点矩形单元有限元分析过程四节点矩形单元有限元分析过程 另外，可以验证形函数另外两个性质：另外，可以验证形函数另外两个性质：（1 1） 同理对于其余三个形函数同理对于其余三个形函数 ,()1,(,)0,iiiijjN x yN xyiy,1()(1)(1)141()(1)(1)041()(1)(1)041()(1)(1)04kkklkkmkknnnabNx yababNx yababNx yababNx yab,(),( , ),( , )lmnNx yNx yNx y 。 Tia

11、njin University四节点矩形单元有限元分析过程四节点矩形单元有限元分析过程 Tianjin University（2 2） 即在单元内任意一点处的形函数之和等于即在单元内任意一点处的形函数之和等于1 1。,( , )1l m nii kN x y,1( , )(1)(1)(1)(1)(1)(1)(1)(1)4l m nii kxyxyxyxyNx yabababab1(1)(11)(1)(11)4xyyxyyabbabb 12(1)2(1)4xxaa1四节点矩形单元有限元分析过程四节点矩形单元有限元分析过程 Tianjin University五、几何方程求解应变将位移插值函数代入

12、几何方程中：将位移插值函数代入几何方程中： eexyyxBNxyyx00形函数矩阵经过微分算子矩阵作用后得到形函数矩阵经过微分算子矩阵作用后得到3 38 8几何矩阵：几何矩阵： ybxaybxaybxaybxaxaxaxaxaybybybybabB0000000041四节点矩形单元有限元分析过程四节点矩形单元有限元分析过程 Tianjin University六、物理方程求解应力由平面问题物理方程可得：由平面问题物理方程可得：其中：其中： eeDDBS21010(1)1002ED )(21)(21)(21)(21)(21)(21)(21)(21)()()()()()()()()()()()()

13、()()()()1 (42ybxaybxaybxaybxaxaybxaybxaybxaybxaybxaybxaybxaybabES因此，应力矩阵因此，应力矩阵 为：为： S四节点矩形单元有限元分析过程四节点矩形单元有限元分析过程 Tianjin University结论：结论： 对于平面四节点矩形单元，其单元上的应力、应变不对于平面四节点矩形单元，其单元上的应力、应变不再是常数，而是在一定程度上呈线性变化，即：再是常数，而是在一定程度上呈线性变化，即： 方向的方向的正应力和正应变随正应力和正应变随 坐标线性变化；坐标线性变化； 方向的正应力和正方向的正应力和正应变随应变随 坐标线性变化；剪应力

14、沿坐标线性变化；剪应力沿 坐标和坐标和 坐标均成线坐标均成线性变化。性变化。 因此，若在弹性体中采用相同数目的节点时，矩形单因此，若在弹性体中采用相同数目的节点时，矩形单元的精度要比常应变三角形单元的精度高。元的精度要比常应变三角形单元的精度高。 xyyxxy四节点矩形单元有限元分析过程四节点矩形单元有限元分析过程 Tianjin University七、运用虚功原理求解 eK 由虚功原理，节点力在节点的虚位移上所做的虚功应等由虚功原理，节点力在节点的虚位移上所做的虚功应等于单元内部应力在虚应变上所做的虚功，即内力虚功于单元内部应力在虚应变上所做的虚功，即内力虚功= =外力虚外力虚功，也即：功

15、，也即： ;WQ eTeTAFtdxdy将将 TTeTB和和 eBD代入上式，可得：代入上式，可得： TeeFBDBt A eTKBDB t A由此，可得：由此，可得：四节点矩形单元有限元分析过程四节点矩形单元有限元分析过程 Tianjin University 2222222222222222211111111111()(1)()(1 3 )()(1)(1)(1 3 )328348628681111111111()(1 3 )(1)(1)()(1 3 )()32868628341(1)ebaabbaabbaabbaababababababababEhKab 22222222222222222

16、21111111()(1)(1)(1 3 )()(1)3286862811111111()(1 3 )()(1)()32834862111111()(1)()(1 3 )328348111 ()328baabbaabbaababababababbaabbaababa 对 称2222221(1 3 )(1)6111()(1)32811()32babbaabab 四节点矩形单元有限元分析过程四节点矩形单元有限元分析过程 Tianjin University)1)(1 (41)1)(1 (41)1)(1 (41)1)(1 (41byaxNbyaxNbyaxNbyaxNnmlk1(1)(1)4rrrN

17、 引入无量纲坐标：引入无量纲坐标：(, , )rk l m n,xyab,rrrrxyab四节点矩形单元有限元分析过程四节点矩形单元有限元分析过程 Tianjin University 1xyxyuubxvvayabuvuvabyx由几何方程可得单元应变场表达式：由几何方程可得单元应变场表达式：可记为：可记为： eB四节点矩形单元有限元分析过程四节点矩形单元有限元分析过程 Tianjin University几何矩阵可表示成分块形式：几何矩阵可表示成分块形式： klmnBBBBB其中：其中： 0(1)01100(1)4(1)(1)rrrrrrrrrrrrrNbbNBaaabababNNab (

18、, , )rk l m n四节点矩形单元有限元分析过程四节点矩形单元有限元分析过程 Tianjin University由应力与应变关系，可得单元应力场表达式：由应力与应变关系，可得单元应力场表达式： eeDDBS应力矩阵可表示成分块形式：应力矩阵可表示成分块形式： klmnSSSSS其中：其中： 2(1)(1)(1)(1)4(1)11(1)(1)22rrrrrrrrrrrrrrbaESD Bbaabab (, , )rk l m n对于平面应变问题：对于平面应变问题：2,11EE四节点矩形单元有限元分析过程四节点矩形单元有限元分析过程 Tianjin University kkklkmkne

19、lklllmlnmkmlmmmnnknlnmnnKKKKKKKKKKKKKKKKK1112221224(1)rskkEtKkkab( , , )r sk l m n其中：其中： TrsrsKBDBt A即：即：单元刚度矩阵：单元刚度矩阵：四节点矩形单元有限元分析过程四节点矩形单元有限元分析过程 Tianjin University1112221224(1)rskkEtKkkab( , , )r sk l m n221112212222111(1)(1)3231()21()2111(1)(1)323rsrsrsrsrsrsrsrsrsrsrsrskbakabkabkab 其中：其中：四节点矩形单

20、元有限元分析过程四节点矩形单元有限元分析过程 Tianjin University八、由 合成 eK K刚度集成法：刚度集成法：首先求出各单元的贡献矩阵，然后将它们叠加首先求出各单元的贡献矩阵，然后将它们叠加形成整体刚度矩阵。但是由于编程时需先将各单元的贡献矩形成整体刚度矩阵。但是由于编程时需先将各单元的贡献矩阵储存起来，而贡献矩阵的阶数与整体刚度矩阵阶数相同，阵储存起来，而贡献矩阵的阶数与整体刚度矩阵阶数相同，因此占用非常大空间，不利于节约空间资源。因此占用非常大空间，不利于节约空间资源。单元定位数组法：单元定位数组法：将单元的节点位移编码按照节点顺序排成将单元的节点位移编码按照节点顺序排成

21、一行形成一维数组，利用各单元的定位数组，采用一行形成一维数组，利用各单元的定位数组，采用“边定位，边定位，边累加边累加”的方法。的方法。 四节点矩形单元有限元分析过程四节点矩形单元有限元分析过程 Tianjin University九、建立节点荷载列阵 = F +dEFF节点荷载列阵的组成：节点荷载列阵的组成： 其中，其中， 为节点荷载，为节点荷载， 为等效节点荷载。为等效节点荷载。 可按照虚功等效原则求解，即将单元内的荷载移置到节可按照虚功等效原则求解，即将单元内的荷载移置到节点上后，应当与原荷载所作虚功等效。点上后，应当与原荷载所作虚功等效。集中力、分布体力（均质等厚单元自重）、分布面力（

22、均布集中力、分布体力（均质等厚单元自重）、分布面力（均布侧压、侧压、X X方向均布荷载、方向均布荷载、X X方向三角形荷载）方向三角形荷载） FdFEFE四节点矩形单元有限元分析过程四节点矩形单元有限元分析过程 Tianjin University十、处理位移约束条件（1 1）降阶法：若第）降阶法：若第r r个自由度方向位移分量为个自由度方向位移分量为0 0，则将整体，则将整体刚度矩阵第刚度矩阵第r r行，第行，第r r列划掉，后一行上移，右一列左移，这列划掉，后一行上移，右一列左移，这样总刚减少一阶，未知数减少一个。样总刚减少一阶，未知数减少一个。 112233447032044200132

23、14021203274420021413212000304427032324021201321042003274021200214130uvuqhtvEtuqhtvuv11440uvuv例：例： 22337442413212034270322120130uqhtvEtuqhtv四节点矩形单元有限元分析过程四节点矩形单元有限元分析过程 Tianjin University（2 2）对角元素置）对角元素置1 1法：法： 11111211222122221212ininiiiiiniinnninnnnRKKKKRKKKKKKKKRKKKKR 11111121221222221212000inniiiiiiniinnnnnnniRKKKKKKKRKKKKKRKKKRK 1111112122122222120000100inniinnnnnnniRKKKKKKKRKKKKRKi例：已知位移边界条件例：已知位移边界条件 （