2-三角形有限单元法

第2讲三角形有限单元法－弹性力学平面问题的求解方程组其物理和几何方程可表示如下：εDσ.210101000称对DD在平面问题里，弹性力学的求解未知量为：TxyyxσTxyyxεTvuu物理方程几何方程L.uεxyyx00L(2.1)(2.2)－瑞利李兹法的位移函数对于右图所示的任意形状的分析物体，要找到一个满足位移边界条件的、全域的位移函数，是很困难的。其中的待定系数Am、Bm没有特定的含义。mmBAfVU,FDKmmmmyxFBvyxFAu,,2100mmBA－有限单元法的思路弹性力学方程组直接求解的困难在于难于找到一个全域的精确函数.在有限单元法里，这个问题可以通过定义分片插值的位移或应力函数来解决；ijmxy分析域上的单元和结点离散任意单元ijm－用结点位移作为待定系数假设单元的位移函数任意单元ijmmmjjiiuyxNuyxNuyxNu),(),(),(对于任意单元(i,j,m)，以结点位移(ui,uj,um)为待定系数，可以给出该单元的插值函数ybxbbu321:假设mmmjjjiiiybxbbuybxbbuybxbbu321321321:在各结点处满足321,,:bbb可求出??;?;mjiNNN如何在三角形单元上构造出上述的单元位移函数？？？－用结点位移作为待定系数假设单元的位移函数假设分析区域被划分成了n个结点、m个单元，对二维分析，本问题的求解未知量为n个结点上的位移(ui,vi)对于任意单元(i,j,m)，以结点位移(ui,uj,um)为待定系数，可以很容易给出该单元的插值函数：mmjjiiuyxNuyxNuyxNu),(),(),(其中Ni(x,y)为各结点坐标(xi,yi)为表示的函数。－三角形单元位移函数的构造过程yxu321设该三角形单元的位移插值函数为：将三个结点的坐标和待求结点位移代入上式得：mmmjjjiiiyxuyxuyxu321321321(2.3)(2.4)－三角形单元位移函数的构造过程求解该线性方程组(2.4)即可得到三个待定系数(β1β2β3)：其中将(β1β2β3)代入单元插值位移插值函数(2.3)，经过整理后得到：其中Ni称为单元的形函数：ycxbaAyxNiiii21),(),,(mjixxcyybyxyxamjimjijmmji(2.5)－三角形单元位移函数的构造过程mmjjiiuyxNuyxNuyxNu),(),(),(同理得到v方向的位移函数：mmjjiivyxNvyxNvyxNv),(),(),(将单元插值位移函数写成矩阵的形式：emjimjivuNaaaaNNNuiiivua(2.6)－矩阵形式的三角形单元位移函数简记为mmjjiimjimjiuuvuvuNNNNNNyxvyxu000000),(),(u其中，－单元的应变矩阵将单元插值函数(2.6)代入几何方程(2.1)得单元应变：eemjiBaaBBBεiiiiiiiiibccbAxNyNyNxN002100B(2.7)(2.8)其中应变矩阵：mmjjiimmjjiimjimjiuuvuvuxNyNxNyNxNyNyNyNyNxNxNxN000000.uLε简记为xyyx00L－单元的应变矩阵的特点从(2.8)式可以发现三角形单元的应变矩阵是常量矩阵，故对于存在应力梯度等问题的求解，三角形单元的精度较低。它的优点是原理简单，能较好地离散曲线边界。－单元的应力矩阵将单元应变(2.7)代入物理方程(2.2)得单元应力：(2.9)(2.10)eeSaDBaDεσ其中单元应力矩阵DBS－弹性体的总势能将单元位移(2.6)、应变(2.7)和应力(2.9)式代入上式可以求出单元势能和弹性体的总势能：(2.11)eSTTTfdddTufuσε21弹性体的总势能可表示为eSeeeeeeeeetdstdxdytdxdy21TTTTTTTNafNaaDBBa(2.12)－最小势能原理(2.13)根据最小势能原理与弹性力学求解体系的等价性，对总势能取驻值：(2.14a)0eeeeeFaK得整体平衡方程tdxdyeeBDBKTeesetdstdxdyTNfNFTT(2.14b)其中－有限元求解方程组(2.15)按对号入座的原则对单元方程进行组集，得到整体求解的结构控制方程：FaKK为2n×2n的方阵，F为2n的列阵。求解该线性方程组可以得到每个结点的位移(ui,vi)值；由(2.7)和(2.9)式可以求出单元内任意点的应力和应变。－有限元求解过程总结有限元法通过单元离散和最小势能原理，把定解条件下的微分方程组的求解巧妙地转化为线性方程组的运算，它的核心是单元分析；假定单元的位移函数[u]=[N]{a}求出单元应变ε和应力σ求出单元势能πe求出总势能π=∑πe用最小势能原理δπ=0得到整体求解方程[K]{a}={F}不同的单元类型或数值计算方法在上述各个具体步骤上会有一些自己的特色或特点，但是它们的的基本计算过程脱不开上述的框架。－瑞利－李兹法vs有限单元法mmBAfVU,FDKmmmmyxFBvyxFAu,,2100mmBA瑞利－李兹法有限单元法mmjjiiuyxNuyxNuyxNu),(),(),(iivugVU,0iuFUKijmxy2.分析实例－模型及网格划分划分为4个三角形单元，6个结点，故本题有限元求解的方程组的总规模为：xy2m2m2m2m02N/m2N/m(a)1121121212FaKxy1N/m123456(1)(b)(2)(3)(4)图(a)所示的分析物体，根据对称型取1/4作为分析模型，如图(b)所示结点自由度编号：结点号自由度号11;223;435;647;859;10611;12－单元弹性矩阵该问题可以简化为平面应力问题求解，假设弹性模量为E，泊松比0，则各单元的弹性矩阵为5.000101称对ED210101000称对DDxy1N/m123456(1)(b)(2)(3)(4)－单元应变矩阵3(3)25ijm任取一单元③，设其局部编号为，则其应变矩阵为：110110101000010001BiiiiiiiiibccbAxNyNyNxN002100B),,(mjixxcyybyxyxamjimjijmmjixy1N/m123456(1)(b)(2)(3)(4)3(3)25ijm－单元刚度矩阵单元③的刚度矩阵为：2000200110110110114000202211031021213TeKBDBtA349105634Et91056③tdxdyeeBDBKTxy1N/m123456(1)(b)(2)(3)(4)3(3)25ijm单元③的刚度矩阵2000200110110110114000202211031021213TeKBDBtA349105634Et91056③－单元刚度矩阵的组集按对号入座的原则将单元刚度矩阵Ke③叠加到整体刚度矩阵K中:************************************123456789101112123456789101112K－所有单元组集后的刚度矩阵再按同样的方法，最后得到本题的有限元求解方程为000000000010101120261114101602121031213016121016140061216012014665544332211vuvuvuvuvuvuEt称对－引入位移边界条件该问题的位移边界条件为0421uuu在上面的方程中消去对应的结点位移为0的项，得到最后的求解方程组为0000012200006121061206106224653321uuvuvvEt称对0654vvvxy1N/m123456(1)(b)(2)(3)(4)－求解位移、应变和应力从上式中求解得到结点位移有了结点位移后，将其代入可以求出单元应变；176.0176.0374.0088.0253.1253.31653321EtuuvuvveaB代入可以求出单元应力。D1N/m123456(1)(2)(3)(4)3.作业(1)设弹性模量为E=1，泊松比为0.5，重新求解上述分析实例。123(2)如下图所示的三角形单元，设各结点的坐标为节点1（0，0）、节点2（-1，0）、节点3（0，-1）；已知各节点的位移为u1=2、u2=1、u3=3，利用三角形单元的插值函数求点（-0.2,-0.2）处的位移值