彈性問題有限元方法

1、北京航空航天大學(xué)第3講 彈性問題的有限元方法北京航空航天大學(xué)3.1 彈性問題變分原理彈性問題變分原理3.2 連續(xù)體彈性問題的有限元求解過程連續(xù)體彈性問題的有限元求解過程平面問題平面問題 (三角形單元、四邊形單元三角形單元、四邊形單元)軸對稱問題軸對稱問題三維問題（四面體單元、六面體單元）三維問題（四面體單元、六面體單元）3.3 有限元分析中若干問題的探討有限元分析中若干問題的探討位移函數(shù)的構(gòu)造要求位移函數(shù)的構(gòu)造要求單元形狀函數(shù)的性質(zhì)單元形狀函數(shù)的性質(zhì)剛度矩陣（單元、整體）的特點剛度矩陣（單元、整體）的特點位移邊界條件的處理位移邊界條件的處理位移單元解的下限性位移單元解的下限性第3講 彈性問題的

2、有限元方法北京航空航天大學(xué)。 3.1 彈性問題變分原理彈性問題變分原理北京航空航天大學(xué)彈性勢能1()2xxyyzzxyxyyzyzxzxzUdV 012ijijdVUdVijijxzxzyzyzxyxyzzyyxxU21)(210單位體積的應(yīng)變能（應(yīng)變比能）：單位體積的應(yīng)變能（應(yīng)變比能）：北京航空航天大學(xué)微元變形能（正應(yīng)變）111() ()222xxxxdUF udydzdxdV uF北京航空航天大學(xué)外力做功ib()()pxyzxyzSWb u b vb w dVp up vp w dApiiiiSWbu dVpu dAip單位體積力面力北京航空航天大學(xué)彈性體系統(tǒng)的總勢能對于保守力場作用下的彈

3、性體系統(tǒng)總勢能：北京航空航天大學(xué)彈性體變分原理最小勢能原理1II212pTTpijijiiiiSTSdVbu dVpu dAdVdVdA u bu p北京航空航天大學(xué)最小勢能原理的等價性1II2pijijiiiiSdVbu dVpu dA 1II2pijklijkliiiiSDdVbu dVpu dA ppijiijiijklijkliiiiSijijiiiiSuuDdVb u dVp u dAdVb u dVp u dA 北京航空航天大學(xué)分部積分分部積分,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北京航空航天大學(xué)平衡方程力邊界條件等價于,ppppijijiiiiSiji jij jiiiiiSSij jiiij jiiSdVb u dVp u dAu l dAu dVb u dVp u dAbu dVlpu dA ,0pij jiiij jiiSbu dVlpu dA北京航空航天大學(xué)l以上證明過程說明：總可以找到滿足位移邊界條件的試函數(shù)(許可位移場)，在滿足幾何方程和物理方程的前提下，其結(jié)果可以精確滿足所剩下

5、的平衡方程和力邊界條件。l實際上，選擇的許可位移場有相當(dāng)?shù)木窒扌院兔つ啃浴Ｒ话愫茈y將真正精確的位移場包含在許可位移場中。最小勢能原理在許可位移場中求出最好的一組解。北京航空航天大學(xué)虛功原理關(guān)于虛位移原理外力作用下處于平衡狀態(tài)的彈性體，產(chǎn)外力作用下處于平衡狀態(tài)的彈性體，產(chǎn)生約束許可的微小虛位移（并同時在彈生約束許可的微小虛位移（并同時在彈性體內(nèi)產(chǎn)生虛應(yīng)變），外力在虛位移上性體內(nèi)產(chǎn)生虛應(yīng)變），外力在虛位移上所作的虛功等于彈性體內(nèi)各點的應(yīng)力在所作的虛功等于彈性體內(nèi)各點的應(yīng)力在相應(yīng)的虛應(yīng)變上所作的虛功。相應(yīng)的虛應(yīng)變上所作的虛功。北京航空航天大學(xué)dVxzxzyzyzxyxyzzyyxx)(dAwpvpu

6、pdVwbvbubzySxzyxp)()(dAupdVubdViSiiiijijp彈性問題中等價于最小勢能原理！北京航空航天大學(xué)比較：虛功原理和能量變分原理l虛功原理是理論力學(xué)上的一個根本性原理，可以用于一切非線性力學(xué)問題。l最小勢能原理只是虛功原理對彈性體導(dǎo)出的一種表述形式，但是對于線彈性問題，最小勢能原理的應(yīng)用非常方便。l能量變分原理方法可以很方便的擴(kuò)展到結(jié)構(gòu)位移場以外的不含非線性的領(lǐng)域，如求解熱傳導(dǎo)、電磁場、流體動力學(xué)等問題。北京航空航天大學(xué)關(guān)于集中力的說明1II2pijijiiiiiiSUWdVbu dVpu dAPu 1II2pijijiiiiSUWdVbu dVpu dA 體積力體

7、積力分布面力分布面力集中力集中力單獨考慮集中力單獨考慮集中力外力載荷北京航空航天大學(xué)3.2 連續(xù)體彈性問題的有限元求解過程首先看一個簡單的平面問題：p材料：低碳鋼體力：重力（密度 ）面力：p=1 N/mm2厚度：t等腰直角三角形腰長: l=20mm求：頂點處的位移？北京航空航天大學(xué)平面問題的有限元求解過程l幾何離散幾何離散：三角形單元或四邊形單元：三角形單元或四邊形單元 三角形單元三角形單元平面問題中最簡單的單元平面問題中最簡單的單元l單元特征分析單元特征分析構(gòu)造位移函數(shù)構(gòu)造位移函數(shù)單元應(yīng)變能單元應(yīng)變能單元外力功（單元等效節(jié)點力）單元外力功（單元等效節(jié)點力）l單元集成單元集成：系統(tǒng)的總勢能：系

8、統(tǒng)的總勢能l變分處理變分處理：系統(tǒng)的平衡方程（組）：系統(tǒng)的平衡方程（組）l應(yīng)用位移邊界條件應(yīng)用位移邊界條件求出節(jié)點位移求出節(jié)點位移l由節(jié)點位移由節(jié)點位移求出單元的應(yīng)變、應(yīng)力求出單元的應(yīng)變、應(yīng)力北京航空航天大學(xué)Step 1. 幾何離散采用3節(jié)點三角形單元662211vuvuvuq整體節(jié)點整體節(jié)點位移列陣位移列陣p體力：重力（密度 ）厚度：t0p p表面力單位體積力0gb112266xyxyxyPPPPPPP整體等效節(jié)整體等效節(jié)點力列陣點力列陣北京航空航天大學(xué)Step 2. 單元分析構(gòu)造單元位移函數(shù)l位移函數(shù)（模式）是指單元內(nèi)位移分布狀態(tài)，事先并不知道，合理選擇一種函數(shù)來逼近這種分布是有限元分析計

9、算過程中關(guān)鍵性的一環(huán)。l在實際應(yīng)用中普遍采用的是多項式函數(shù)，這是因為多項式函數(shù)的數(shù)學(xué)運算（微分和積分）比較方便，而且所有光滑函數(shù)的局部都可以用多項式來逼近。l關(guān)于多項式的項數(shù)和階次，要根據(jù)單元的節(jié)點自由度數(shù)和有關(guān)解的收斂性要求來確定。對于平面問題，位移函數(shù)如下：2212345622123456uxyxyxyvxyxyxy 北京航空航天大學(xué)123123uxyvxy112131212232312333uxyuxyuxy構(gòu)造位移函數(shù)：對u利用節(jié)點條件： i j m編號對應(yīng)關(guān)系：(局部 整體)先采用局部編號，最后換成整體北京航空航天大學(xué)111122223333111uxyuxyuxy1112233uu

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

11、c y uab xc y uAvab xc y vab xc y vab xc y vA北京航空航天大學(xué))(21ycxbaANiiii111232123233000000uvNNNuuNNNvvuv ( , )( , )ex yx yuNqN單元形狀函數(shù)矩陣qe 單元節(jié)點位移矩陣北京航空航天大學(xué) 0( , )( , )0 ( , )( , )( , )( , ) xeeyxyxu x yx yx yx yx yv x yyyx uNqBq應(yīng)變矩陣123123 0 0 0 0( , ) 0 0 0 0 xNNNx yNNNyyx BNStep 2. 單元分析應(yīng)變北京航空航天大學(xué)123123123

12、112233 0 0 01( , )0 0 0 2 bbbx ycccAcbcbcbBBBB312112233112233 0 0 01110 0 0 222 bbbcccAAAcbcbcbBBB北京航空航天大學(xué)21 0( , ) 1 0( , )11-0 0 2xxyyxyxyEx yx yD平面應(yīng)力：Step 2. 單元分析應(yīng)力( , )( , )( , )( , )eex yx yx yx yDDBqSq 123123( , ) x y SDBD BBBSSS2 2 (1)1-1- 22iiiiiiiibcEbcAcbSDB應(yīng)力矩陣平面應(yīng)變：用平面應(yīng)變彈性矩陣代入得到類似結(jié)果。北京航空航

13、天大學(xué)單元應(yīng)變能：111222TTTeeeeSUdVdVtdxdy D D111 222TTTeeeeTeeTeeeeSSUtdxdytdxdyqB DBqqB DBqqK qStep 3. 單元分析單元勢能111213212223313233 eeeeeTTeeeSeeetdxdytAKKKKB DBB DBKKKKKK單元剛度矩陣北京航空航天大學(xué)211 22114 (1) 22rsrsrsrseTrsrsrsrsrsrsb bc cb cc bEtKtAAb cc bc cb bB DB( ,1,2,3)r s 北京航空航天大學(xué)TTTTTTeeeeppeeeSSWdVdAdVdAu bu

14、pq N bq N pTTTTTTeeeeppeeeSSlWdVdAtdxdytdlqN bN pqN bN p TTTeepeeeeSlWtdxdytdlq PPN bN p單元等效節(jié)點力列陣單元外力功：北京航空航天大學(xué)KqqTeeUU2141eeTeeUqKq21擴(kuò)充疊加12TTUW q Kqq P41eTeWWq PTeeeW q P擴(kuò)充疊加Step 4. 單元集成系統(tǒng)勢能北京航空航天大學(xué)關(guān)于單元剛度矩陣的擴(kuò)充疊加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北京航空航天大學(xué)關(guān)于單元等效節(jié)點載荷列陣的擴(kuò)充疊加TeeeW q P12310 0emeTenniejWPq PuuuuuPPmij單元編號 mij北京航空航天大學(xué)Step 5. 變分處理變分處理12TTUW q Kqq P00 PKqPKq Step 6: 處理位移邊界條件并求解Step 7: 計算每個單元的應(yīng)變及應(yīng)力北京航空航天大學(xué)關(guān)于三角形單元l3節(jié)點三角形單元是常應(yīng)變（常應(yīng)力）單元。在應(yīng)變梯度較大的部位（亦即應(yīng)力梯度較大的部位），單元劃分應(yīng)適當(dāng)密集，否則不能反映真實的應(yīng)變變化而導(dǎo)

16、致較大的誤差。l提高計算精度的其它措施采用高精度三角形單元（2次單元、3次單元）采用四邊形單元（1次單元、2次單元）北京航空航天大學(xué)4節(jié)點四邊形單元北京航空航天大學(xué)12341234uxyxyvxyxy112131411212232422312333433412434444uxyx yuxyx yuxyx yuxyx y構(gòu)造位移函數(shù)：112131411212232422312333433412434444vxyx yvxyx yvxyx yvxyx y對u,v分別利用節(jié)點條件：北京航空航天大學(xué)1111112222223333334444441 1 1 1 uxyx yuxyx yuxyx yux

17、yx y112213344uuuuT11112222333344441 1 1 1 xyx yxyx yxyx yxyx yTTTT*1*TT的伴隨矩陣對于一般四邊形，逆矩陣的表達(dá)式比較復(fù)雜。北京航空航天大學(xué)1 12233441 1223 344uN uN uN uN uvN vN vN vN v12341234uxyxyvxyxy112213344uuuuT112213344vvvvT 1,2,3,4iiiiiNab xc yd xyi北京航空航天大學(xué)1122123431234344000 0000 0 uvuvNNNNuuNNNNvvuv ( , )( , )ex yx yuNqN單元形狀

18、函數(shù)矩陣qe 單元節(jié)點位移矩陣北京航空航天大學(xué)特例：4節(jié)點矩形單元0010020030041(1)(1)41(1)(1)4 1(1)(1)41(1)(1)4xxyyNabxxyyNabxxyyNabxxyyNab00,xy 矩形單元的重心坐標(biāo)12341234000 0( , )000 0 NNNNx yNNNNN北京航空航天大學(xué)12341234 0 000 0( , ) 0 000 0 xNNNNx yNNNNyyx BN( , )TeeeSStdxdyx y dxdyKB DBFTTeepeSltdxdytdlPN bN p北京航空航天大學(xué)l對于一般的四邊形單元，在總體坐標(biāo)系下構(gòu)造位移插值函

19、數(shù)，則計算形狀函數(shù)矩陣、單元剛度矩陣及等效節(jié)點載荷列陣時十分冗繁；而對于矩形單元，相應(yīng)的計算要簡單的多。l矩形單元明顯的缺點是不能很好的符合曲線邊界，因此可以采用矩形單元和三角形單元混合使用。更為一般的方法是通過等參變換將局部自然坐標(biāo)系內(nèi)的規(guī)格化矩形單元變換為總體坐標(biāo)系內(nèi)的任意四邊形單元（包括高次曲邊四邊形單元）。北京航空航天大學(xué)受均布內(nèi)壓作用的長圓筒軸對稱問題的有限元求解過程北京航空航天大學(xué)l研究軸對稱問題時通常采用圓柱坐標(biāo)系(r,z)，以z軸為對稱軸l由于對稱性: 4個應(yīng)力分量，4個應(yīng)變分量，2個位移分量rzzrrzzrwur北京航空航天大學(xué)3節(jié)點三角形軸對稱單元360度環(huán)形單元橫截面為3

20、節(jié)點三角形r z平面內(nèi)，單元節(jié)點位移和節(jié)點力為：北京航空航天大學(xué)構(gòu)造位移函數(shù)：利用節(jié)點條件：北京航空航天大學(xué)北京航空航天大學(xué)北京航空航天大學(xué)三維問題的有限元求解過程l離散時采用體單元：四面體或六面體l求解步驟和平面問題完全一樣l單元分析的時候?qū)⒍S擴(kuò)充到三維北京航空航天大學(xué)mmmjjjiiimjimjimjiwvuwvuwvuNNNNNNNNNwvu 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ( , , )( , , )ex y zx y zuNq形狀函數(shù)矩陣北京航空航天大學(xué) ( , , )( , , )( , , )( , , )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應(yīng)變矩陣北京航空航天大學(xué)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北京航空航天大學(xué) TTeepTeSeeedVdAWPN bN pq P單元等效節(jié)點力列陣北京航空航天大學(xué)3.3 有限元分析中的若干問題探討l單元位移函數(shù)的

22、構(gòu)造單元位移函數(shù)的構(gòu)造l單元形狀函數(shù)的性質(zhì)單元形狀函數(shù)的性質(zhì)l單元剛度矩陣的特點單元剛度矩陣的特點l整體剛度矩陣的特點整體剛度矩陣的特點l位移邊界條件的處理位移邊界條件的處理l位移單元解的下限性位移單元解的下限性北京航空航天大學(xué)單元位移函數(shù)的構(gòu)造滿足收斂性要求收斂收斂單元尺寸趨于零時，有限元解趨于真解單元尺寸趨于零時，有限元解趨于真解北京航空航天大學(xué)l準(zhǔn)則1：完備性包含常應(yīng)變項和剛體位移項如果在勢能泛函中所出現(xiàn)的位移函數(shù)的最高階導(dǎo)數(shù)是m階，則選取的位移函數(shù)至少是m階完全多項式。l準(zhǔn)則2：協(xié)調(diào)性相鄰單元公共邊界保持位移連續(xù)如果在勢能泛函中所出現(xiàn)的位移函數(shù)的最高階導(dǎo)數(shù)是m階，則位移函數(shù)在單元交界面

23、上必須具有直至(m-1)階的連續(xù)導(dǎo)數(shù)，即Cm-1連續(xù)性。北京航空航天大學(xué)l如果在單元交界面上位移不連續(xù)，表現(xiàn)為當(dāng)結(jié)構(gòu)變形時將在相鄰單元間產(chǎn)生縫隙或重疊，這意味著將引起無限大的應(yīng)變，這時必然會發(fā)生交界面上的附加應(yīng)變能補充到系統(tǒng)的應(yīng)變能中去，有限元解就不可能收斂于真正解。 北京航空航天大學(xué)多項式的Pascal模式l構(gòu)造一個單元的位移模式時，應(yīng)參考由多項式函數(shù)構(gòu)成的Pascal三角形或四面體從低階到高階多項式的項數(shù)由節(jié)點位移條件確定221 xyxxyy平面問題Pascal三角形三維問題Pascal四面體北京航空航天大學(xué)lC0型單元勢能泛函中所出現(xiàn)的位移函數(shù)的最高階導(dǎo)數(shù)是1階，在單元交界面上具有0階的

24、連續(xù)導(dǎo)數(shù)(平面問題單元、空間問題單元)。lC1型單元勢能泛函中所出現(xiàn)的位移函數(shù)的最高階導(dǎo)數(shù)是2階，在單元交界面上具有1階的連續(xù)導(dǎo)數(shù)（梁單元、板殼單元等）。北京航空航天大學(xué)關(guān)于非協(xié)調(diào)單元l當(dāng)單元的位移函數(shù)滿足完備性要求時，稱單元是完備的（通常較容易滿足）。當(dāng)單元的位移函數(shù)滿足協(xié)調(diào)性要求時，稱單元是協(xié)調(diào)的。l當(dāng)勢能泛函中位移函數(shù)的導(dǎo)數(shù)是2階時，要求位移函數(shù)在單元的交界面上具有C1或更高的連續(xù)性，這時構(gòu)造單元的插值函數(shù)往往比較困難。在某些情況下，可以放松對協(xié)調(diào)性的要求，只要單元能夠通過分塊試驗 (Patch test)，有限元分析的解答仍然可以收斂于正確的解。這種單元稱為非協(xié)調(diào)單元。分塊檢驗由首先提出，已經(jīng)證明它給出了收斂性的充分條件。