第四章平面问题有限元分析及程序设计

1、第四章第四章 平面问题有限元分析及程序设计平面问题有限元分析及程序设计4.1 平面问题单元离散4.2 平面问题单元位移模式4.3 平面问题单元分析4.4 平面问题整体分析4.5 平面问题有限元程序设计有限单元法及程序设计有限单元法及程序设计有限元网格划分的基本原则有限元网格划分的基本原则网格数目网格疏密单元阶次网格质量网格分界面和分界点位移协调性网格布局结点和单元编号网格自动剖分网格数量20万最小网格尺度150m最大网格尺度3500m平面问题单元：平面问题单元：平面应力：平面应力：三角形板三角形板平面应变：平面应变：三棱柱三棱柱平面问题结点：平面问题结点：平面问题约束：平面问题约束：平面问题荷

2、载：平面问题荷载：三角形单元三角形单元绞结点绞结点绞支座、链杆绞支座、链杆结点荷载和非结点荷载结点荷载和非结点荷载几个重要概念：几个重要概念：4.1 平面问题单元离散平面问题单元离散基本量和方程的矩阵表示基本量和方程的矩阵表示xyfff体积力体积力xyfff面力面力应力应力xyxy应变应变xyxy基本量和方程的矩阵表示基本量和方程的矩阵表示位移位移udv 物理方程物理方程2101011002xxyyxyxyE简写为简写为D只要知道了单元的位移函数，就可由几何方程求出应变，再由物只要知道了单元的位移函数，就可由几何方程求出应变，再由物理方程就可求出应力。理方程就可求出应力。几何方程：几何方程：4

3、.2 单元位移模式单元位移模式 Tyuxvyvxu Tmmjjiievuvuvu有限单元法：未知量是有限单元法：未知量是结点结点的位移分量的位移分量那么单元内任意一点的位移跟结点位移有什么关系呢？那么单元内任意一点的位移跟结点位移有什么关系呢？因此说，只要知道了位移场的分布，即可解决上述问题。因此说，只要知道了位移场的分布，即可解决上述问题。i (xi, yi)位移模式：单元位移场分布形式位移模式：单元位移场分布形式4.2 单元位移模式单元位移模式iujujvivmumvxyj (xj, yj)m (xm, ym),(yxfuu),(yxfvv),(iiuiyxfu ),(jjujyxfu )

4、,(mmumyxfu ),(iiviyxfv ),(jjvjyxfv ),(mmvmyxfv 建立一个坐标系，如下图所示：建立一个坐标系，如下图所示：假定位移模式如下所示：假定位移模式如下所示：三个结点的位移也必定满足三个结点的位移也必定满足位移场函数，因此有：位移场函数，因此有：位移函数的选取是任意的，所选取的位移函数越接近于真实情况，所位移函数的选取是任意的，所选取的位移函数越接近于真实情况，所求得的形变和内力结果就越准确。求得的形变和内力结果就越准确。yxu321 最简单的位移场函数是线性函数，即最简单的位移场函数是线性函数，即：yxv654 iiiuyx321 iiivyx654 jj

5、juyx321 jjjvyx654 mmmuyx321 mmmvyx654 位移模式的选取位移模式的选取边界条件：边界条件：在三个结点也应满足位移场函数；在三个结点也应满足位移场函数；i 结点结点j 结点结点m 结点结点其中，其中， 、 、 、 、 、 是系数，由边界条件求得。是系数，由边界条件求得。1231234.2 单元位移模式单元位移模式4.2 单元位移模式单元位移模式写成矩阵形式123456100000011000=000110000001iiiiiijjjjjjmmmmmmxyuxyvxyuxyvxyuxyvyxu321 2A iycxbaNiiimmiyaxyx jjmiyb1y1

6、- jmixc1x1 j位移模式的选取位移模式的选取因此，因此， 、 、 是是 、 、 的线性函数；的线性函数；123iujumu同样，同样， 、 、 是是 、 、 的线性函数；的线性函数；455ivjvmv代入位移场函数，则代入位移场函数，则 是是 、 、 的线性函数，即：的线性函数，即：uiujumummjjiiuNuNuN其中，其中， 、 、 是系数，是是系数，是 、 的线性函数；的线性函数；iNjNmNxy可以求得：可以求得：mmjjiiyxyxyxA11121 其中：其中：注意：注意：i, j, m 必须是逆时针必须是逆时针排列，否则面排列，否则面积为负。积为负。(i, j, m )

7、同理，可求得同理，可求得 、 ，且下标可轮换，且下标可轮换 ；jNmNmmjjiivNvNvNyxv654 同理可得：同理可得：4.2 单元位移模式单元位移模式j上式也可以写成：上式也可以写成：形函数的性质：形函数的性质：mmjjiimmjjiyxyxyxyxyxyxN111111 iiNjiNmiNj i oiNoiN jiidsNAidxdyN(i, j, m) 、 、 表明了单元的位移表明了单元的位移形态（位移在单元的变化规律）形态（位移在单元的变化规律）iNjNmN称为形态函数，简称称为形态函数，简称形函数形函数 、 、 是坐标是坐标 (x 、y ) 的的线性函数；线性函数；iNjNm

8、Nim11/21/3j i21A311002131(i, j, m)4.2 单元位移模式单元位移模式位移函数：位移函数：mmjjiiuNuNuNyxu321 mmjjiivNvNvNyxv654 由于由于 、 、 是坐标是坐标 (x 、y ) 的线性函数，的线性函数，iNjNmN),(mjiuuu因此，因此， u、v 也是也是 、 的线性函数。的线性函数。),(mjivvvjimiujumujimivjvmv4.2 单元位移模式单元位移模式因此，因此， u、v 在坐标空间应该为一平面。在坐标空间应该为一平面。位移写成向量形式：位移写成向量形式： mmjjiimmjjiivNvNvNuNuNuN

9、vud mmjjiimjimjivuvuvuNNNNNN000000 eNdmmmjjjiiimmmjjjiiivNuvNuvNuvuNvuNvuN000000 N称为形函数矩阵。称为形函数矩阵。4.2 单元位移模式单元位移模式有限元分析中，应力转换矩阵、刚度矩阵都是依赖于位移模式建立有限元分析中，应力转换矩阵、刚度矩阵都是依赖于位移模式建立起来的，因此，位移模式必须能够反映弹性体的真实位移形态，才起来的，因此，位移模式必须能够反映弹性体的真实位移形态，才能得到正确的解答。能得到正确的解答。位移模式需要满足的条件：位移模式需要满足的条件：（1）位移模式必须能够反映单元的刚体位移；）位移模式必须

10、能够反映单元的刚体位移；（2）位移模式必须能够反映单元的常应变；）位移模式必须能够反映单元的常应变；（3）位移模式尽可能反映位移的连续性；）位移模式尽可能反映位移的连续性；必要条件必要条件充分条件充分条件yxu321 yxv654 xvvyuu0010 u40 v2 35刚体平动刚体平动刚体转动刚体转动作业：作业：P141 6-14.2 单元位移模式单元位移模式53531222xyy53534622yxx单元应变单元应变4.3 平面问题的单元分析平面问题的单元分析4.3.1 单元的应变向量单元的应变向量 xvyuyvxummjjiimmjjiimjimjivuvuvuxNyNxNyNxNyNy

11、NyNyNxNxNxN000000mmjjiimmjjiimjimjivuvuvubcbcbccccbbbA00000021由几何方程求。由几何方程求。4.3 平面问题的单元分析平面问题的单元分析4.3.1 单元的应变向量单元的应变向量 eB mmjjiimjimjibcbcbccccbbbAB00000021 iiiiibccbAB0021可以简写为可以简写为:mjiBBB其中其中: B是单元的应变矩阵，且是单元的应变矩阵，且:所以所以:(i, j, m)常量常量因此，单元应变是常数。因此，单元应变是常数。所以，三角形单元又称为所以，三角形单元又称为常应变单元常应变单元。4.3 平面问题的单

12、元分析平面问题的单元分析4.3.2 单元的应力向量单元的应力向量 eBDD eeSBD mjiSSSS iiiiiiibccbcbAES2121)1 (2221E1E单元应力单元应力由物理方程求。由物理方程求。其中其中:其中其中: S是单元的应力矩阵，且是单元的应力矩阵，且:平面应力平面应力:平面应变平面应变:(i, j, m)常量常量因此，应力也与坐标无关，所以单元应力是常数。因此，应力也与坐标无关，所以单元应力是常数。4.3 平面问题的单元分析平面问题的单元分析单元性质分析单元性质分析位移位移是是x、y的线性函数；的线性函数；误差是误差是 x、 y的二阶小量；的二阶小量；应变应变应力应力常

13、量；常量；相邻单元连续相邻单元连续相邻单元不连续相邻单元不连续误差是误差是 x、 y的一阶小量；的一阶小量；提高精度方法：提高精度方法：1）减小单元尺寸；）减小单元尺寸；2）采用高次位移函数，提高位移、应变和应力的精度；）采用高次位移函数，提高位移、应变和应力的精度；收敛速度和精度估计 若单元的插值函数是若单元的插值函数是完备完备而而协调协调的，当单元尺寸不的，当单元尺寸不断缩小而趋于零时，有限元解将趋于真正解。断缩小而趋于零时，有限元解将趋于真正解。 在有些情况下，如果用于有限元场函数近似解的多在有些情况下，如果用于有限元场函数近似解的多项式能精确地拟合真正解，则在有限数目的单元划分（项式能

14、精确地拟合真正解，则在有限数目的单元划分（甚至仅仅是一个单元）的条件下，也能得到精确的解答甚至仅仅是一个单元）的条件下，也能得到精确的解答。例如真正解是二次函数，若有限元的插值函数也包含。例如真正解是二次函数，若有限元的插值函数也包含了二次的完全多项式，则有限元解就能得到精确的解答了二次的完全多项式，则有限元解就能得到精确的解答。由此我们可以得到精度与单元尺寸的关系。例如位移。由此我们可以得到精度与单元尺寸的关系。例如位移可以展开成可以展开成TaylorTaylor级数：级数：iiiuuuuxyxy 这只是形式上的精度估计，并不能对有限元解的误差做出具体的估计。而后者在这只是形式上的精度估计，

15、并不能对有限元解的误差做出具体的估计。而后者在实际分析工作中更有用。一般可以通过两种途径解决：实际分析工作中更有用。一般可以通过两种途径解决：单元结点力单元结点力4.3 平面问题的单元分析平面问题的单元分析4.3.3 单元的单元的 刚度矩阵刚度矩阵 TmymxjyjxiyixemjieFFFFFFFFFF eeSBD Tmmjjiiemjievuvuvu Tmmjjiiemjievuvuvu* eB*单元结点位移单元结点位移单元应力向量：单元应力向量：给定一个虚位移：给定一个虚位移：单元虚应变：单元虚应变：i (xi, yi)iujujvivmumvxyj (xj, yj)m (xm, ym)

16、虚功原理虚功原理: :内力虚功等于外力虚功内力虚功等于外力虚功4.3 平面问题的单元分析平面问题的单元分析4.3.3 单元的单元的 刚度矩阵刚度矩阵 eTeF* tdxdyFATeTe* eATetdxdyBDBF eek tdxdyAT* tdxdyBDBAeTe* tdxdyBDBAeTTe* eATTetdxdyBDB*t 为单元厚度为单元厚度由于虚位移是任意给定的可能位移，故：由于虚位移是任意给定的可能位移，故： ek其中，其中， 是单元的刚度矩阵是单元的刚度矩阵 (6 6)；4.3 平面问题的单元分析平面问题的单元分析 emmemjemiejmejjejieimeijeiiekkkk

17、kkkkkk AtBDBksTrerssrsrsrsrsrsrsrsrersbbcccbbcbccbccbbAEtk21212121)1 (42 tdxdyBDBkATe tdxdyBDBAT tABDBT单元的刚度矩阵为：单元的刚度矩阵为：写成分块矩阵：写成分块矩阵：其中：其中：平面应力：平面应力：平面应变：平面应变：;,mjir ;,mjis 21E1E 211111122222211111122222211224 1iiiiiiiiijijijijimimimimiiiiiiiijijiijijm im iimimijijjijiebbccbccbbbccbccbbbccbccbbccb

18、ccbbb cc bccbbb cc bccbbbbccb cc bEtkA111122221111112222221111122222jjjjjjjjjmjmjmjmijijijijjjjjjjjjmjmjjmjmimimm im ijmjmmjmjmmmb bc cb cc bb bc cb cc bbccbccbbb cc bc cb bb cc bc cb bbbccb cc bb bc cb cc bb bc12111111222222mmmmmimimimimjmjmjmjmmmmmmmmmcb cc bbccbccbbb cc bc cb bb cc bc cb b写成元素矩阵：

19、写成元素矩阵：4.3 平面问题的单元分析平面问题的单元分析ersk单元刚度矩阵的特点：单元刚度矩阵的特点：1) 对称性：对称性：esrk2) 与单元尺寸无关，放大或缩小尺寸，单元刚度矩阵不变；与单元尺寸无关，放大或缩小尺寸，单元刚度矩阵不变；3）奇异性：它不存在逆阵）奇异性：它不存在逆阵4）主元（对角线元素）恒正）主元（对角线元素）恒正4.3 平面问题的单元分析平面问题的单元分析单元的刚度矩阵为：单元的刚度矩阵为： 2222222222222222102121002102121212121212121021002121021)1 (4 hlhlhlhhlhhlhlhlhhllhllhlhhlh

20、lhlhllhllhllhlhllAEtk单元的刚度矩阵为：单元的刚度矩阵为： kk作用在单元上的荷载，既有结点荷载，也有非结点荷载，作用在单元上的荷载，既有结点荷载，也有非结点荷载，因此需要将非结点荷载转换成等效的结点荷载。因此需要将非结点荷载转换成等效的结点荷载。4.3.4 等效结点荷载等效结点荷载等效结点荷载和原荷载在任何虚位移产生的虚功相等；等效结点荷载和原荷载在任何虚位移产生的虚功相等；刚体：刚体：等效原则：等效原则：4.3 平面问题的单元分析平面问题的单元分析结点荷载：结点荷载：非结点荷载：非结点荷载：直接集成到荷载列向量；直接集成到荷载列向量；等效成结点荷载；等效成结点荷载；原荷

21、载与等效结点荷载在任一轴上的投影之和相原荷载与等效结点荷载在任一轴上的投影之和相等，对任一轴的力矩之和也相等。等，对任一轴的力矩之和也相等。等效结点荷载向量：等效结点荷载向量： TLmyLmxLjyLjxLiyLixeLFFFFFFF静力等效静力等效4.3 平面问题的单元分析平面问题的单元分析 eTNyxvyxud*),(),( Tyxggf Tmmjjiievuvuvu* tdxdyfdFATeLTe*1）体积力的等效结点荷载）体积力的等效结点荷载单元结点为单元结点为 i, j, m，密度为，密度为 ，任意一点的体积力向量为：，任意一点的体积力向量为：假设单元各结点发生虚位移：假设单元各结点

22、发生虚位移：则单元内任意一点的虚位移为：则单元内任意一点的虚位移为：根据虚功原理：根据虚功原理： tdxdyfNATTe* tdxdyfNATTe*4.3 平面问题的单元分析平面问题的单元分析 tdxdyfNFATeL1）体积力的等效结点荷载）体积力的等效结点荷载因此，则有：因此，则有：一般情况下，重力与一般情况下，重力与y轴轴方向相反：方向相反： Tgf 00LmxLjxLixFFFdxdyNgtFAiLiygtA31因此，则有：因此，则有：LiyLmyLjyFFFgtA31 0000=000000TeLAiijAjmmTijmAFNfdxdy tNNNdxdy tNgNNgNgNgNdxd

23、y t结论：三个结点各承担结论：三个结点各承担总荷载的三分之一。总荷载的三分之一。4.3 平面问题的单元分析平面问题的单元分析2）分布荷载的等效结点荷载）分布荷载的等效结点荷载ixyjmqxqyq Tyxqqqf TeLqijFNfds t如图所示均布荷载，集度如图所示均布荷载，集度为为q，则：，则：则等效结点荷载为：则等效结点荷载为： *TeTeLqijFdfds t *TeTqijNfds t *TeTqijNfds tj itqFFxLjxLix210LmxFj itqFFyLjyLiy210LmyF 000=000TeLqijiixjyjijmmTxiyixjyjxmymijFNfds

24、 tNNqNds tqNNNq Nq Nq Nq Nq Nq Nds t结论：沿该荷载作用的结论：沿该荷载作用的边上的两个结点各承担边上的两个结点各承担总荷载的一半。总荷载的一半。如图所示线性分布荷载：荷载如图所示线性分布荷载：荷载与与x轴夹角为轴夹角为xyijm TeLqijTxiyixjyjxmymijFNfds tq Nq Nq Nq Nq Nq Nds tqcosxsqql1isNl 1122cossincossin003 23 23 23 2TeLqtlqtlqtlqtlF xyijmiljliiFqlslss TeLqijTxiyixjyjxmymijFNfds tq Nq Nq

25、Nq Nq Nq Nds t1isNl Fcos1Fsin1FcosFsin00TeiiiiLllllFtllll 4.4.1 整体刚度矩阵的集成步骤1 1、定位、定位单元单元结点结点编号编号 2 2、累加、累加整体整体结点结点编号编号 ijrs单元刚度系数单元刚度系数 整体刚度系数整体刚度系数4.4 整体刚度矩阵的集成整体刚度矩阵的集成eiikrrKeijkrsKejiksrKejjkssKmteimkrtKejmkstKemiktrKemjktsKemmkttK4.4.1 整体刚度矩阵的集成步骤例：例：P121P1214.4 整体刚度矩阵的集成整体刚度矩阵的集成单元号IIIIIIIV局部编码整体自由度编码i i159311 i2610412j j11395 J224106m m13759 m248610以以整体自由度编码整体自由度编码总刚的集成总刚的集成4.4.2 边界条件的处理方法4.4 整体刚度矩阵的集成整体刚度矩阵的集成1 1）划行划列法）划行划列法2 2）0 0、1 1置换法（填