深入淺出講清楚有限元法

1、(完好word版)深入淺出的講清楚有限元法(完好word版)深入淺出的講清楚有限元法23/23(完好word版)深入淺出的講清楚有限元法“有限元法基礎(chǔ)及應(yīng)用”增補(bǔ)講義（一）顧克秋2005年3月）一、引子彈簧單元與彈簧系統(tǒng)目標(biāo)：掌握失散構(gòu)造直接剛度法解析的原理和形式。認(rèn)識(shí)有限元位移法列式的形式和基本觀點(diǎn)。1、典型彈簧單元解析彈簧單元描述：彈簧的物理特點(diǎn)：圖1-12個(gè)節(jié)點(diǎn)：i,j圖1-2節(jié)點(diǎn)位移：ui,uj已知彈簧力位移關(guān)系：Fkk彈簧剛度節(jié)點(diǎn)力：fi,fjujui彈簧伸長(zhǎng)量單元自由度：2F彈簧力，拉伸為正考慮彈簧變形平衡時(shí)的條件和彈簧物理特點(diǎn)，獲得以下方程：fiFk(ujui)kuikujfjF

2、k(ujui)kui（1-1）kuj寫(xiě)成矩陣形式：fikkuifjkk（1-2）uj寫(xiě)成矩陣符號(hào)形式：fkd（1-3）1式（1-2）、（1-3）為彈簧單元的剛度方程，反響了單元特點(diǎn)：節(jié)點(diǎn)力與節(jié)點(diǎn)位移之間的關(guān)系。式中：k彈簧單元的剛度矩陣d單元節(jié)點(diǎn)位移列陣f單元節(jié)點(diǎn)力列陣（注意：?jiǎn)卧?jié)點(diǎn)力是節(jié)點(diǎn)對(duì)單元的作使勁）彈簧單元?jiǎng)偠确匠套h論：1）k有何特色？對(duì)稱(chēng)、奇怪、主對(duì)角元素恒正2）k中元素代表什么含義？剛度系數(shù)大小等于彈簧剛度；每列元素代表一端固定、另一端產(chǎn)生單位位移時(shí)加在彈簧單元上的節(jié)點(diǎn)力。3）上邊單元?jiǎng)偠确匠棠軌蚯蠼鈫?？為何？不能夠夠。剛度方程只是表征一個(gè)典型單元的彈性特點(diǎn)，單元水平上沒(méi)法確立單

3、元節(jié)點(diǎn)位移。只有把系統(tǒng)中所有單元特點(diǎn)集成后，在系統(tǒng)水平上才可能求出所有未知位移和反力。單元水平上，若已知單元的節(jié)點(diǎn)位移，可由剛度方程求出所有單元節(jié)點(diǎn)力重量。若節(jié)點(diǎn)力已知，單元節(jié)點(diǎn)位移不能夠確立，單元可作剛體運(yùn)動(dòng)（小位移）。這也是單元?jiǎng)偠染仃嚻婀中缘奈锢碇v解。2、彈簧系統(tǒng)整體解析原理以右圖的一個(gè)彈簧系統(tǒng)為例，研究怎樣由單元特點(diǎn)集成系統(tǒng)特點(diǎn)并建立對(duì)系統(tǒng)進(jìn)行求解的控制方程。由前面獲得的彈簧單元的剛度方程公式（1-2），分別寫(xiě)出2個(gè)彈簧單元的特點(diǎn)方程以下：圖1-3單元1（1-4）單元2（1-5）（注：右端節(jié)點(diǎn)力重量的下標(biāo)1，2為單元節(jié)點(diǎn)的局部編號(hào)，上標(biāo)是單元號(hào)）下邊按兩個(gè)方法完成系統(tǒng)特點(diǎn)的裝置和控制方

4、程的建立。并在特定條件下求解。1）由節(jié)點(diǎn)平衡方程導(dǎo)出：系統(tǒng)處于平衡時(shí)，考慮各節(jié)點(diǎn)（1，2，3節(jié)點(diǎn)）的平衡條件：節(jié)點(diǎn)碰到的外載荷與節(jié)點(diǎn)碰到與其連接的所有單元對(duì)其作使勁（單元節(jié)點(diǎn)力的反作使勁）之和等于零。所以有以下（節(jié)點(diǎn)）平衡方程（組）：Ff111Ff1f2221（1-6）Ff232把單元特點(diǎn)（1-4），（1-5）代入（1-6）獲得：F1k1u1k1u2F2k1u1（1-7）(k1k2)u2k2u3F3k2u2k2u3寫(xiě)成矩陣形式：（1-8）3或矩陣符號(hào)形式：KDF（1-9）式（1-8），（1-9）就是系統(tǒng)平衡方程，該方程建立了失散系統(tǒng)的外載荷與節(jié)點(diǎn)位移之間的關(guān)系，是求解節(jié)點(diǎn)位移的控制方程。KDF

5、彈簧系統(tǒng)的構(gòu)造總剛度矩陣系統(tǒng)節(jié)點(diǎn)位移列陣系統(tǒng)節(jié)點(diǎn)載荷列陣議論：（1）K有那些特色和性質(zhì)？2）上述方程能求解嗎？由單元?jiǎng)偠确匠摊B加導(dǎo)出將單元1，2的剛度方程（1-4），（1-5）進(jìn)行增廣（擴(kuò)大到系統(tǒng)規(guī)模）：（1-10）（1-11）上述兩個(gè)矩陣方程疊加，得：（1-12）4上式中代入節(jié)點(diǎn)力平衡關(guān)系（1-6），就獲得與（1-8）同樣的節(jié)點(diǎn)平衡方程。上述兩種方法都一定考慮1）單元特點(diǎn)集成；2）失散構(gòu)造的節(jié)點(diǎn)上外載荷（系統(tǒng)外力）與節(jié)點(diǎn)力（系統(tǒng)內(nèi)力）的平衡。所以方程（1-8）的實(shí)質(zhì)是節(jié)點(diǎn)的力平衡關(guān)系，左側(cè)是由節(jié)點(diǎn)位移表示的（總）節(jié)點(diǎn)力，右側(cè)是節(jié)點(diǎn)所受外載荷。3）給定載荷和拘束條件下的求解設(shè)界線(xiàn)條件為：u10

6、（1-13）F2F3P則節(jié)點(diǎn)平衡方程（1-8）變化為：（1-14）該方程組睜開(kāi)后分為2個(gè)部分：第2，3個(gè)方程變化為：（1-15）第1個(gè)方程變化為：（1-16）先后解方程（1-15）、（1-16）獲得：（1-17）（1-18）5進(jìn)而解出了系統(tǒng)的未知位移和未知反力，并能夠進(jìn)一步求彈簧力。3、例題圖1-4所示一個(gè)3個(gè)彈簧的系統(tǒng)。k1100N/mm,k2200N/mm,k3100N/mm,P500N,u1u40求：（a）系統(tǒng)總剛度矩陣b）節(jié)點(diǎn)2，3的位移c）節(jié)點(diǎn)1、4的反力d）彈簧2中的力解：a）:寫(xiě)出各單元?jiǎng)偠染仃嚕簣D1-4應(yīng)用疊加法直接獲得系統(tǒng)總剛度矩陣：或：該總剛度矩陣特色：對(duì)稱(chēng)性、奇怪性、稀有

7、、非零元素沿主對(duì)角線(xiàn)呈帶狀散布。b）:參照前面的做法（1-8）式）和求出的總剛度矩陣，寫(xiě)出系統(tǒng)節(jié)點(diǎn)平衡方程：（1-19）考慮到位移界線(xiàn)條件：u1u40則平衡方程組（1-19）第2，3方程化為：求解上式得：（c）：7由（1-19）的方程1，4得：（d）：彈簧2內(nèi)力為：F2k22k2(u3u2)20032200(N)4、練習(xí)題對(duì)圖示彈簧系統(tǒng)，試用疊加法求其總剛度矩陣。并依據(jù)節(jié)點(diǎn)平衡方程的含義，試一試由各單元?jiǎng)偠染仃嚨脑刂苯訉?xiě)出總剛度矩陣的非零元素。（拉力）圖1-5二、桿單元目標(biāo)：經(jīng)過(guò)桿單元特點(diǎn)方程的建立，初步掌握有限元法單元解析的過(guò)程和原理。認(rèn)識(shí)桿系構(gòu)造解析的原理。81、等截面桿單元及其剛度矩陣

8、研究2節(jié)點(diǎn)等截面桿單元：?jiǎn)卧系牧W(xué)量和基本關(guān)系以下：L桿長(zhǎng)A截面積E彈性模量dudxEuu(x)桿單元位移(x)桿單元應(yīng)變(x)桿單元應(yīng)力圖2-1應(yīng)變位移關(guān)系：（2-1）單元節(jié)點(diǎn)位移：應(yīng)力應(yīng)變關(guān)系：（2-2）單元節(jié)點(diǎn)力：ffuiidfjuj下邊研究桿單元的單元特點(diǎn)。1）直接法導(dǎo)出桿單元特點(diǎn)采納資料力學(xué)基本知識(shí)對(duì)單元進(jìn)行力學(xué)解析：桿單元伸長(zhǎng)量：ujui（2-3）9桿應(yīng)變：LE桿應(yīng)力：EL桿內(nèi)力：EAEAFAkLL桿的軸向剛度：EAkL2-4）2-5）2-6）2-7）因?yàn)檩S向變形模式下，桿單元的行為與彈簧單元同樣，所以可比較彈簧單元的剛度方程（1-2），考慮到（2-7），直接寫(xiě)出桿單元的剛度方程

9、：fiEA11ui（9）fjL11uj2-8寫(xiě)成符號(hào)形式：fkd桿單元?jiǎng)偠染仃嚍椋篍A11k（2-10）L112）公式法導(dǎo)出桿單元特點(diǎn)步驟以下：（1）在單元上定義近似位移場(chǎng)把一個(gè)單元上的位移散布假設(shè)為簡(jiǎn)單多項(xiàng)式函數(shù)。有限元法中用插值法經(jīng)過(guò)節(jié)點(diǎn)位移重量作為待定參數(shù)來(lái)構(gòu)造單元位移函數(shù)。對(duì)圖2-1的桿單元，方便起見(jiàn)引入局部坐標(biāo)10因?yàn)樵摋U單元只有2個(gè)未知位移重量，所以單元上假設(shè)的簡(jiǎn)單位移函數(shù)采納一次多項(xiàng)式。故對(duì)單元的節(jié)點(diǎn)位移進(jìn)行線(xiàn)性插值。則簡(jiǎn)單定義出節(jié)點(diǎn)的插值函數(shù)以下：（2-11）所以單元上近似位移函數(shù)的插值形式為：（2-12）該位移函數(shù)也稱(chēng)為單元的位移模式，這里是線(xiàn)性位移模式。式（2-11）中的插

10、值函數(shù)又稱(chēng)為形狀函數(shù)，簡(jiǎn)稱(chēng)形函數(shù)。式（2-12）寫(xiě)成矩陣形式為：上式中N稱(chēng)為單元的形函數(shù)矩陣。ui（2-13）uNiNjujNd式（2-13）是有限元法中最重要的關(guān)系式之一，經(jīng)過(guò)該式把單元上的近似位移散布函數(shù)用節(jié)點(diǎn)位移來(lái)表示，為進(jìn)行單元層次上的解析打下了基礎(chǔ)。（2）單元應(yīng)變和單元應(yīng)力由桿一維變形的應(yīng)變位移方程（2-1）和單元的位移函數(shù)（2-13）求出單元的應(yīng)變散布和節(jié)點(diǎn)位移的關(guān)系：dudNdBd（2-14）dxdx式中：BdNi()Nj()1/L1/L（2-15）dx稱(chēng)為單元的位移應(yīng)變變換矩陣，簡(jiǎn)稱(chēng)應(yīng)變矩陣。由一維桿的應(yīng)力應(yīng)變關(guān)系（2-2），得單元應(yīng)力和單元節(jié)點(diǎn)位移的關(guān)系：EEBd（2-16）

11、（3）用彈性體的虛位移原理導(dǎo)出桿單元?jiǎng)偠确匠?1變形體的虛位移：設(shè)想在彈性體上發(fā)生的，知足位移允許條件（內(nèi)部連續(xù)，界線(xiàn)協(xié)調(diào)）的細(xì)小、任意位移場(chǎng)。能夠理解為某個(gè)位移場(chǎng)的細(xì)小擾動(dòng)（變分）。虛位移的特色：1）設(shè)想的，與真實(shí)位移沒(méi)關(guān)；2）幾何上是允許的：連續(xù)、協(xié)調(diào)；3）細(xì)小、任意大小。虛位移原理：彈性體受力平衡時(shí)，若發(fā)生虛位移，則外力虛功等于彈性體內(nèi)的虛應(yīng)變能（應(yīng)力在虛應(yīng)變上做的虛功）。下邊把虛位移原理應(yīng)用在所研究的桿單元上。定義桿單元的虛位移：節(jié)點(diǎn)虛位移單元虛位移單元虛應(yīng)變節(jié)點(diǎn)虛位移：duiuj單元虛位移：uNdd單元虛應(yīng)變：(u)Bddx那么，節(jié)點(diǎn)力虛功：dTf單元虛應(yīng)變能：TdVdTBTEBdd

12、VdTBTEBdVdVVV對(duì)桿單元應(yīng)用虛功原理，那么上述節(jié)點(diǎn)力（外力）虛功等于虛應(yīng)變能，所以有以下關(guān)系：dTfdTBTEBdVd（2-17）V（2-18）考慮到d的任意性，從上式能夠獲得：12fBTEBdVdkdV上式就是桿單元的剛度方程，桿單元的剛度矩陣為：kBTEBdV（2-19）V這就是單元?jiǎng)偠染仃嚨耐ㄊ?，其?dǎo)出原理和計(jì)算方法可實(shí)行到其余種類(lèi)的單元。詳細(xì)計(jì)算以下：（2-20）顯然，與前面直接法獲得的單元?jiǎng)偠染仃嚕?-8）式同樣。3）桿單元議論只有拉伸、壓縮變形的桿單元在局部坐標(biāo)系下是一維問(wèn)題，2節(jié)點(diǎn)單元只有2個(gè)節(jié)點(diǎn)位移重量單元有2個(gè)自由度，單元?jiǎng)偠确匠?、剛度矩陣?階。單元?jiǎng)偠染仃囋氐?/p>

13、物理意義：設(shè)單元?jiǎng)偠确匠虨椋篺ik11k12ui（2-21）fjk21k22uj令：ui1（2-22）uj0帶入（2-21）獲得：fik（2-23）11fjk21上式表示，單元?jiǎng)偠染仃嚨谝涣性鼐褪钱?dāng)單元節(jié)點(diǎn)位移知足式（2-22）時(shí)的單元節(jié)點(diǎn)力重量。若是能想法求出此時(shí)的節(jié)點(diǎn)力，就獲得第一列的剛度系數(shù)。13所以，一般地，單元?jiǎng)偠染仃嚨牡趇（i=1,2）列元素表示當(dāng)保持單元的第i個(gè)自由度位移為，其余自由度位移為0時(shí)，施加在單元上的節(jié)點(diǎn)力重量。（也能夠用此方法直接導(dǎo)出桿單元的剛度矩陣元素，試練習(xí)）單元?jiǎng)偠染仃噷?duì)稱(chēng)、奇怪、主元恒正。4）例題例1求圖示桿中的應(yīng)力。解：構(gòu)造劃分為個(gè)桿單元，單元之間在節(jié)點(diǎn)鉸

14、接。依據(jù)桿單元?jiǎng)偠染仃嚨墓椒謩e寫(xiě)出兩個(gè)單元的剛度矩陣為：圖2-2參照前面彈簧系統(tǒng)的解析方法，裝置桿系統(tǒng)的有限元方程（平衡方程）以下：220u1F1EA（2-24）231u2F2Lu3F3011考慮圖2-2中的拘束u1u30和載荷狀況后，方程（2-24）變化為：2200F1EA1u2P23（2-25）L110F30則上式的第2個(gè)方程為：（2-26）14求解該方程后獲得系統(tǒng)的位移解：u1PL0u21（2-27）3EAu30計(jì)算應(yīng)力：?jiǎn)卧?1E1E1Eu2u1EPL0PLLL3EA3A單元22E2E2Eu3u2E0PLPLLL3EA3A提示：1）本例中單元應(yīng)力的計(jì)算采納了資料力學(xué)中的方法,與采納有

15、限元單元應(yīng)力公式EEBd的結(jié)果同樣，請(qǐng)考據(jù)。2）對(duì)錐形桿，單元截面積能夠用平均值，進(jìn)而轉(zhuǎn)變?yōu)榻票绢}的問(wèn)題求解。3）求應(yīng)力從前需要先求出節(jié)點(diǎn)位移，所以本方法稱(chēng)為有限元位移法。4）若是桿上受連續(xù)散布的軸向載荷或節(jié)點(diǎn)之間受軸向集中載荷，解析時(shí)能夠依據(jù)虛功相等的原則先把單元上的載荷等效移置到節(jié)點(diǎn)上。例2已知：圖2-3求：桿兩端的支反力。單元、節(jié)點(diǎn)的定義如圖2-3。先檢查桿右端（節(jié)點(diǎn)3）與墻壁可否接觸。15計(jì)算右端的自由伸長(zhǎng)：所以，右端縫隙將閉合。系統(tǒng)平衡方程為：（2-28）引入以下載荷和位移界線(xiàn)條件：則有限元平衡方程（2-28）成為：（2-29）分別出第二個(gè)方程：（2-30）即：（2-31）解得：所

16、有位移解為：16（2-32）依據(jù)上式位移解，從系統(tǒng)平衡方程（2-28）的第1，3個(gè)方程分別求出支反力以下：解畢。2、2-D和3-D空間中的桿單元（平面和空間桁架單元）1）2-D空間中的桿單元2-D空間中建立桿單元的基本思路是依據(jù)前面在桿的一維局部坐標(biāo)系下建立的單元特點(diǎn)方程經(jīng)過(guò)坐標(biāo)變換，變換為2-D整體坐標(biāo)系下的方程，同時(shí)獲得坐標(biāo)變換后的單元?jiǎng)偠染仃?。而系統(tǒng)整體解析的原理和方法與一維狀況同樣。圖2-4（1）變換圖2-4為一個(gè)桿單元及其局部坐標(biāo)系與2-D整體坐標(biāo)的關(guān)系。節(jié)點(diǎn)的位移重量和節(jié)點(diǎn)力重量在2-D局部坐標(biāo)系x-y下描述，桿節(jié)點(diǎn)i具17有2個(gè)自由度：位移重量為ui，vi；節(jié)點(diǎn)力重量為fxi，f

17、yi此中只有x方向的位移重量和節(jié)點(diǎn)力重量用來(lái)描述單元特點(diǎn)。節(jié)點(diǎn)上的位移和節(jié)點(diǎn)力向量在2-D局部坐標(biāo)系與2-D整體坐標(biāo)系下的變換以下：稱(chēng)為方向余弦。上述變換的矩陣形式：（2-33）符號(hào)形式：（2-34）diTdi此中l(wèi)m（2-35）Tlm稱(chēng)為向量的坐標(biāo)變換矩陣，是單元特點(diǎn)坐標(biāo)變換的基本元素。顯然是正TT交矩陣，即：T1TT（2-36）所以，由（2-33）可得單元節(jié)點(diǎn)位移向量的坐標(biāo)變換式以下：（2-36）18或：dTd此中T0T0T比較（2-37）獲得單元節(jié)點(diǎn)力的坐標(biāo)變換式：Tf2）2-D空間剛度矩陣下邊能夠?qū)С鰡卧獎(jiǎng)偠确匠毯蛦卧獎(jiǎng)偠染仃嚨淖鴺?biāo)變換式。已經(jīng)知道桿的一維局部坐標(biāo)系下的剛度方程為：把該

18、方程擴(kuò)大到2-D局部坐標(biāo)系x-y下的4階形式：1010uifxiEA0000vifyiL1010ujfxj0000vjfyj符號(hào)形式：kdf引入單元列陣變換式（2-37），（2-39）得：kTdTf考慮到變換矩陣T的正交性，獲得：TTkTdf或：kdf2-37）2-38）2-39）2-40）2-41）2-42）2-43）2-44）2-45）19此中：式（2-45）就是2-D整體坐標(biāo)系下剛單元的剛度方程，k就是二維空間桿單T（2-46）kTkT元?jiǎng)偠染仃?，其?jì)算式以下：（2-47）此中方向余弦由桿節(jié)點(diǎn)整體坐標(biāo)求得。（3）單元應(yīng)力由單元應(yīng)力計(jì)算公式（2-16）和位移向量變換得：即：（2-48）（4

19、）思慮與議論怎樣從二維空間整體坐標(biāo)系下桿單元?jiǎng)偠确匠蹋?-45），依據(jù)剛度矩陣元素的物理意義，直接導(dǎo)出整體坐標(biāo)系下桿單元?jiǎng)偠染仃囉?jì)算式（2-47）？2）例題平面桁架由2根同樣的桿組成（E，A，L）。求：1）節(jié)點(diǎn)2位移2）每根桿應(yīng)力20求出每個(gè)單元在整體坐標(biāo)下的剛度矩陣：?jiǎn)卧?，（i,j）=(1,2)：45，lm22圖2-5k1T1Tk1T11100T01011001EA22110000001100L22001110100011001100000011u1v1u2v21111EA1111（2-49）2L11111111135，l2，m222k2T2Tk2T21100T01011001EA22110000001100L220