弹性问题有限元方法

1、北京航空航天大学第3讲 弹性问题的有限元方法北京航空航天大学3.1 弹性问题变分原理弹性问题变分原理3.2 连续体弹性问题的有限元求解过程连续体弹性问题的有限元求解过程平面问题平面问题 (三角形单元、四边形单元三角形单元、四边形单元)轴对称问题轴对称问题三维问题（四面体单元、六面体单元）三维问题（四面体单元、六面体单元）3.3 有限元分析中若干问题的探讨有限元分析中若干问题的探讨位移函数的构造要求位移函数的构造要求单元形状函数的性质单元形状函数的性质刚度矩阵（单元、整体）的特点刚度矩阵（单元、整体）的特点位移边界条件的处理位移边界条件的处理位移单元解的下限性位移单元解的下限性第3讲 弹性问题的

2、有限元方法北京航空航天大学。 3.1 弹性问题变分原理弹性问题变分原理北京航空航天大学弹性势能1()2xxyyzzxyxyyzyzxzxzUdV 012ijijdVUdVijijxzxzyzyzxyxyzzyyxxU21)(210单位体积的应变能（应变比能）：单位体积的应变能（应变比能）：北京航空航天大学微元变形能（正应变）111() ()222xxxxdUF udydzdxdV uF北京航空航天大学外力做功ib()()pxyzxyzSWb u b vb w dVp up vp w dApiiiiSWbu dVpu dAip单位体积力面力北京航空航天大学弹性体系统的总势能对于保守力场作用下的弹

3、性体系统总势能：北京航空航天大学弹性体变分原理最小势能原理1II212pTTpijijiiiiSTSdVbu dVpu dAdVdVdA u bu p北京航空航天大学最小势能原理的等价性1II2pijijiiiiSdVbu dVpu dA 1II2pijklijkliiiiSDdVbu dVpu dA ppijiijiijklijkliiiiSijijiiiiSuuDdVb u dVp u dAdVb u dVp u dA 北京航空航天大学分部积分分部积分,12() ijijiji jj iiji jijijij jiiji jij jiSdVuudVu dVudVu dVul dAu dV

4、puSSS说明：满足位移边界条件：0iu,pijijiji jij jiSdVul dAu dV 高斯定理高斯定理piji jiji jSSul dAul dA 0 on iuuS北京航空航天大学平衡方程力边界条件等价于,ppppijijiiiiSiji jij jiiiiiSSij jiiij jiiSdVb u dVp u dAu l dAu dVb u dVp u dAbu dVlpu dA ,0pij jiiij jiiSbu dVlpu dA北京航空航天大学l以上证明过程说明：总可以找到满足位移边界条件的试函数(许可位移场)，在满足几何方程和物理方程的前提下，其结果可以精确满足所剩下

5、的平衡方程和力边界条件。l实际上，选择的许可位移场有相当的局限性和盲目性。一般很难将真正精确的位移场包含在许可位移场中。最小势能原理在许可位移场中求出最好的一组解。北京航空航天大学虚功原理关于虚位移原理外力作用下处于平衡状态的弹性体，产外力作用下处于平衡状态的弹性体，产生约束许可的微小虚位移（并同时在弹生约束许可的微小虚位移（并同时在弹性体内产生虚应变），外力在虚位移上性体内产生虚应变），外力在虚位移上所作的虚功等于弹性体内各点的应力在所作的虚功等于弹性体内各点的应力在相应的虚应变上所作的虚功。相应的虚应变上所作的虚功。北京航空航天大学dVxzxzyzyzxyxyzzyyxx)(dAwpvpu

6、pdVwbvbubzySxzyxp)()(dAupdVubdViSiiiijijp弹性问题中等价于最小势能原理！北京航空航天大学比较：虚功原理和能量变分原理l虚功原理是理论力学上的一个根本性原理，可以用于一切非线性力学问题。l最小势能原理只是虚功原理对弹性体导出的一种表述形式，但是对于线弹性问题，最小势能原理的应用非常方便。l能量变分原理方法可以很方便的扩展到结构位移场以外的不含非线性的领域，如求解热传导、电磁场、流体动力学等问题。北京航空航天大学关于集中力的说明1II2pijijiiiiiiSUWdVbu dVpu dAPu 1II2pijijiiiiSUWdVbu dVpu dA 体积力体

7、积力分布面力分布面力集中力集中力单独考虑集中力单独考虑集中力外力载荷北京航空航天大学3.2 连续体弹性问题的有限元求解过程首先看一个简单的平面问题：p材料：低碳钢体力：重力（密度 ）面力：p=1 N/mm2厚度：t等腰直角三角形腰长: l=20mm求：顶点处的位移？北京航空航天大学平面问题的有限元求解过程l几何离散几何离散：三角形单元或四边形单元：三角形单元或四边形单元 三角形单元三角形单元平面问题中最简单的单元平面问题中最简单的单元l单元特征分析单元特征分析构造位移函数构造位移函数单元应变能单元应变能单元外力功（单元等效节点力）单元外力功（单元等效节点力）l单元集成单元集成：系统的总势能：系

8、统的总势能l变分处理变分处理：系统的平衡方程（组）：系统的平衡方程（组）l应用位移边界条件应用位移边界条件求出节点位移求出节点位移l由节点位移由节点位移求出单元的应变、应力求出单元的应变、应力北京航空航天大学Step 1. 几何离散采用3节点三角形单元662211vuvuvuq整体节点整体节点位移列阵位移列阵p体力：重力（密度 ）厚度：t0p p表面力单位体积力0gb112266xyxyxyPPPPPPP整体等效节整体等效节点力列阵点力列阵北京航空航天大学Step 2. 单元分析构造单元位移函数l位移函数（模式）是指单元内位移分布状态，事先并不知道，合理选择一种函数来逼近这种分布是有限元分析计

9、算过程中关键性的一环。l在实际应用中普遍采用的是多项式函数，这是因为多项式函数的数学运算（微分和积分）比较方便，而且所有光滑函数的局部都可以用多项式来逼近。l关于多项式的项数和阶次，要根据单元的节点自由度数和有关解的收敛性要求来确定。对于平面问题，位移函数如下：2212345622123456uxyxyxyvxyxyxy 北京航空航天大学123123uxyvxy112131212232312333uxyuxyuxy构造位移函数：对u利用节点条件： i j m编号对应关系：(局部 整体)先采用局部编号，最后换成整体北京航空航天大学111122223333111uxyuxyuxy1112233uu

10、uT112233111xyxyxyTTTT*1A2TA: 三角形面积*TT的伴随矩阵北京航空航天大学T23322332*3113311312211221x yx yyyxxx yx yyyxxx yx yyyxxTT111123*222123333123abcaaaabcbbbabccccTijmmjijmimjax yx ybyycxx北京航空航天大学11231212323123312aaaubbbuAcccu11231212323123312aaavbbbvAcccv同理可得：因此：1111222233331111222233331()()()21()()()2uab xc y uab x

11、c y uab xc y uAvab xc y vab xc y vab xc y vA北京航空航天大学)(21ycxbaANiiii111232123233000000uvNNNuuNNNvvuv ( , )( , )ex yx yuNqN单元形状函数矩阵qe 单元节点位移矩阵北京航空航天大学 0( , )( , )0 ( , )( , )( , )( , ) xeeyxyxu x yx yx yx yx yv x yyyx uNqBq应变矩阵123123 0 0 0 0( , ) 0 0 0 0 xNNNx yNNNyyx BNStep 2. 单元分析应变北京航空航天大学123123123

12、112233 0 0 01( , )0 0 0 2 bbbx ycccAcbcbcbBBBB312112233112233 0 0 01110 0 0 222 bbbcccAAAcbcbcbBBB北京航空航天大学21 0( , ) 1 0( , )11-0 0 2xxyyxyxyEx yx yD平面应力：Step 2. 单元分析应力( , )( , )( , )( , )eex yx yx yx yDDBqSq 123123( , ) x y SDBD BBBSSS2 2 (1)1-1- 22iiiiiiiibcEbcAcbSDB应力矩阵平面应变：用平面应变弹性矩阵代入得到类似结果。北京航空航

13、天大学单元应变能：111222TTTeeeeSUdVdVtdxdy D D111 222TTTeeeeTeeTeeeeSSUtdxdytdxdyqB DBqqB DBqqK qStep 3. 单元分析单元势能111213212223313233 eeeeeTTeeeSeeetdxdytAKKKKB DBB DBKKKKKK单元刚度矩阵北京航空航天大学211 22114 (1) 22rsrsrsrseTrsrsrsrsrsrsb bc cb cc bEtKtAAb cc bc cb bB DB( ,1,2,3)r s 北京航空航天大学TTTTTTeeeeppeeeSSWdVdAdVdAu bu

14、pq N bq N pTTTTTTeeeeppeeeSSlWdVdAtdxdytdlqN bN pqN bN p TTTeepeeeeSlWtdxdytdlq PPN bN p单元等效节点力列阵单元外力功：北京航空航天大学KqqTeeUU2141eeTeeUqKq21扩充叠加12TTUW q Kqq P41eTeWWq PTeeeW q P扩充叠加Step 4. 单元集成系统势能北京航空航天大学关于单元刚度矩阵的扩充叠加12310 0 0 0 0 00 0 0 0 0 00 0 01 0 0 020 0 eeemmmimjeTnneeeimiiijeejmjiUKKKq KquuuuuKKKK

15、KK1231 00 0 0 0 0 0nenjjuuuuumijmijeeTeeUqKq21单元编号 mij北京航空航天大学关于单元等效节点载荷列阵的扩充叠加TeeeW q P12310 0emeTenniejWPq PuuuuuPPmij单元编号 mij北京航空航天大学Step 5. 变分处理变分处理12TTUW q Kqq P00 PKqPKq Step 6: 处理位移边界条件并求解Step 7: 计算每个单元的应变及应力北京航空航天大学关于三角形单元l3节点三角形单元是常应变（常应力）单元。在应变梯度较大的部位（亦即应力梯度较大的部位），单元划分应适当密集，否则不能反映真实的应变变化而导

16、致较大的误差。l提高计算精度的其它措施采用高精度三角形单元（2次单元、3次单元）采用四边形单元（1次单元、2次单元）北京航空航天大学4节点四边形单元北京航空航天大学12341234uxyxyvxyxy112131411212232422312333433412434444uxyx yuxyx yuxyx yuxyx y构造位移函数：112131411212232422312333433412434444vxyx yvxyx yvxyx yvxyx y对u,v分别利用节点条件：北京航空航天大学1111112222223333334444441 1 1 1 uxyx yuxyx yuxyx yux

17、yx y112213344uuuuT11112222333344441 1 1 1 xyx yxyx yxyx yxyx yTTTT*1*TT的伴随矩阵对于一般四边形，逆矩阵的表达式比较复杂。北京航空航天大学1 12233441 1223 344uN uN uN uN uvN vN vN vN v12341234uxyxyvxyxy112213344uuuuT112213344vvvvT 1,2,3,4iiiiiNab xc yd xyi北京航空航天大学1122123431234344000 0000 0 uvuvNNNNuuNNNNvvuv ( , )( , )ex yx yuNqN单元形状

18、函数矩阵qe 单元节点位移矩阵北京航空航天大学特例：4节点矩形单元0010020030041(1)(1)41(1)(1)4 1(1)(1)41(1)(1)4xxyyNabxxyyNabxxyyNabxxyyNab00,xy 矩形单元的重心坐标12341234000 0( , )000 0 NNNNx yNNNNN北京航空航天大学12341234 0 000 0( , ) 0 000 0 xNNNNx yNNNNyyx BN( , )TeeeSStdxdyx y dxdyKB DBFTTeepeSltdxdytdlPN bN p北京航空航天大学l对于一般的四边形单元，在总体坐标系下构造位移插值函

19、数，则计算形状函数矩阵、单元刚度矩阵及等效节点载荷列阵时十分冗繁；而对于矩形单元，相应的计算要简单的多。l矩形单元明显的缺点是不能很好的符合曲线边界，因此可以采用矩形单元和三角形单元混合使用。更为一般的方法是通过等参变换将局部自然坐标系内的规格化矩形单元变换为总体坐标系内的任意四边形单元（包括高次曲边四边形单元）。北京航空航天大学受均布内压作用的长圆筒轴对称问题的有限元求解过程北京航空航天大学l研究轴对称问题时通常采用圆柱坐标系(r,z)，以z轴为对称轴l由于对称性: 4个应力分量，4个应变分量，2个位移分量rzzrrzzrwur北京航空航天大学3节点三角形轴对称单元360度环形单元横截面为3

20、节点三角形r z平面内，单元节点位移和节点力为：北京航空航天大学构造位移函数：利用节点条件：北京航空航天大学北京航空航天大学北京航空航天大学三维问题的有限元求解过程l离散时采用体单元：四面体或六面体l求解步骤和平面问题完全一样l单元分析的时候将二维扩充到三维北京航空航天大学mmmjjjiiimjimjimjiwvuwvuwvuNNNNNNNNNwvu 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ( , , )( , , )ex y zx y zuNq形状函数矩阵北京航空航天大学 ( , , )( , , )( , , )( , , )eex y zx y zx y

21、zx y z uNqBq 0 00 0 0 0 0 0 0 0 ( , , ) 00 0 ijxyNNzx y zyxzyzxB 0 0 0 0 0 0 0 0 0 0 0 0 0 0 mijmijmNNNNNNN应变矩阵北京航空航天大学eeTdVKB DB)1 (22-1 0 )1 (22-1 0 0 )1 (22-1 0 0 0 1 0 0 0 1 1 0 0 0 1 1 1)21)(1 (1ED对称单元刚度矩阵12TeeeeUq K q北京航空航天大学 TTeepTeSeeedVdAWPN bN pq P单元等效节点力列阵北京航空航天大学3.3 有限元分析中的若干问题探讨l单元位移函数的

22、构造单元位移函数的构造l单元形状函数的性质单元形状函数的性质l单元刚度矩阵的特点单元刚度矩阵的特点l整体刚度矩阵的特点整体刚度矩阵的特点l位移边界条件的处理位移边界条件的处理l位移单元解的下限性位移单元解的下限性北京航空航天大学单元位移函数的构造满足收敛性要求收敛收敛单元尺寸趋于零时，有限元解趋于真解单元尺寸趋于零时，有限元解趋于真解北京航空航天大学l准则1：完备性包含常应变项和刚体位移项如果在势能泛函中所出现的位移函数的最高阶导数是m阶，则选取的位移函数至少是m阶完全多项式。l准则2：协调性相邻单元公共边界保持位移连续如果在势能泛函中所出现的位移函数的最高阶导数是m阶，则位移函数在单元交界面

23、上必须具有直至(m-1)阶的连续导数，即Cm-1连续性。北京航空航天大学l如果在单元交界面上位移不连续，表现为当结构变形时将在相邻单元间产生缝隙或重叠，这意味着将引起无限大的应变，这时必然会发生交界面上的附加应变能补充到系统的应变能中去，有限元解就不可能收敛于真正解。 北京航空航天大学多项式的Pascal模式l构造一个单元的位移模式时，应参考由多项式函数构成的Pascal三角形或四面体从低阶到高阶多项式的项数由节点位移条件确定221 xyxxyy平面问题Pascal三角形三维问题Pascal四面体北京航空航天大学lC0型单元势能泛函中所出现的位移函数的最高阶导数是1阶，在单元交界面上具有0阶的

24、连续导数(平面问题单元、空间问题单元)。lC1型单元势能泛函中所出现的位移函数的最高阶导数是2阶，在单元交界面上具有1阶的连续导数（梁单元、板壳单元等）。北京航空航天大学关于非协调单元l当单元的位移函数满足完备性要求时，称单元是完备的（通常较容易满足）。当单元的位移函数满足协调性要求时，称单元是协调的。l当势能泛函中位移函数的导数是2阶时，要求位移函数在单元的交界面上具有C1或更高的连续性，这时构造单元的插值函数往往比较困难。在某些情况下，可以放松对协调性的要求，只要单元能够通过分块试验 (Patch test)，有限元分析的解答仍然可以收敛于正确的解。这种单元称为非协调单元。分块检验由首先提出，已经证明它给出了收敛性的充分条件。