第10章(非線性有限元)(1)

1、第10章 非線性動(dòng)力有限元法110.1 幾何非線性問題的有限元法210.1.1 幾何非線性問題的牛頓迭代法210.1.2 典型單元的切線剛度矩陣410.2 材料非線性問題的有限元法810.2.1 彈/粘塑性問題的基本表達(dá)式810.2.2 粘塑性應(yīng)變?cè)隽亢蛻?yīng)力增量910.2.3 彈/粘塑性平衡方程1010.3 材料非線性問題的動(dòng)力有限元法1110.4 應(yīng)用舉例1410.4.1 粘彈粘塑性動(dòng)力有限元分析舉例14習(xí)題15第10章 非線性動(dòng)力有限元法當(dāng)機(jī)械結(jié)構(gòu)受到較大的外載荷，或受到持續(xù)時(shí)間較短的沖擊載荷作用時(shí)，結(jié)構(gòu)會(huì)產(chǎn)生過大的變形, 以至于必須考慮結(jié)構(gòu)幾何大變形對(duì)結(jié)構(gòu)整體剛度及固有頻率的影響，即所謂

2、的幾何非線性影響。另外, 對(duì)于多數(shù)非線性動(dòng)力學(xué)問題，還需要考慮材料非線性、接觸非線性等方面的影響。 非線性動(dòng)力學(xué)分析求解的基本方程有如下形式 (4.141)式中，為粘性效應(yīng)項(xiàng)，考慮阻尼、粘塑、粘彈等效應(yīng)。P為外部激勵(lì)。對(duì)于考慮各種非線性效應(yīng)的動(dòng)力學(xué)問題求解，需要對(duì)動(dòng)力學(xué)方程進(jìn)行直接時(shí)間積分。即非線性動(dòng)力有限元分析具有如下特點(diǎn)：（1）問題分析過程需要考慮時(shí)間積分效應(yīng)，不必做模態(tài)分析，不必提取固有頻率；（2）采用直接積分方法求解非線性動(dòng)力學(xué)方程，需要對(duì)時(shí)間作積分計(jì)算，因此計(jì)算量遠(yuǎn)遠(yuǎn)大于線性模態(tài)動(dòng)力學(xué)方法；（3）非線性動(dòng)力學(xué)分析中可以施加不同類型的載荷，包括結(jié)點(diǎn)力、非零位移、單元載荷；（4）在每個(gè)時(shí)

3、間步上，進(jìn)行質(zhì)量、阻尼、及剛度的集成，采用完整矩陣，不涉及質(zhì)量矩陣的近似；（5）可以同時(shí)考慮幾何、材料和接觸等多種非線性效應(yīng)。非線性動(dòng)力有限元分析程序常采用隱式Hilber-Hughes-Taylor法進(jìn)行時(shí)間積分運(yùn)算。這種方法適于模擬非線性結(jié)構(gòu)的動(dòng)態(tài)問題，對(duì)于沖擊、地震等激發(fā)的結(jié)構(gòu)動(dòng)態(tài)響應(yīng)以及一些由于塑性或粘性阻尼造成的能量耗散，隱式算法特別有效。隱式積分方法需要對(duì)剛度矩陣求逆計(jì)算，并通過多次迭代求解增量步平衡方程。隱式Hilber-Hughes-Taylor時(shí)間積分算法為無(wú)條件穩(wěn)定，對(duì)時(shí)間步長(zhǎng)沒有特別的限制。采用子空間法也可以對(duì)動(dòng)力學(xué)平衡方程作時(shí)間積分運(yùn)算。子空間法是提取模態(tài)分析得到的各階

4、特征模態(tài)，并采用與線性模態(tài)動(dòng)力學(xué)分析方法相近的分析方式進(jìn)行求解。對(duì)于帶有微小非線性效應(yīng)的問題，如材料小范圍進(jìn)行入屈服、結(jié)點(diǎn)轉(zhuǎn)角不大的情況，子空間法效率比進(jìn)接積分法要高。此外，非線性動(dòng)力有限元分析還可以采用顯式動(dòng)態(tài)算法，如中心差分法。顯式時(shí)間積分算法為有條件穩(wěn)定，其臨界穩(wěn)定時(shí)間步長(zhǎng)限制了時(shí)間步長(zhǎng)的大小，與有限元模型最小單元尺寸、材料應(yīng)力波速等有關(guān)。顯式時(shí)間積分法適于模擬高速?zèng)_擊、接觸等問題。上述方法的選擇需要綜合考慮計(jì)算量、分析問題的規(guī)模、單元限制等多方面因素，需要豐富的有限元模擬的理論、經(jīng)驗(yàn)和實(shí)踐知識(shí)。以下以幾何非線性問題和材料非線性問題為例介紹非線性有限元法，其中粘彈粘塑性非線性材料問題的分

5、析是典型的非線性動(dòng)力有限元的求解思想。10.1 幾何非線性問題的有限元法幾何非線性問題一般是指物體經(jīng)歷大的剛體位移和轉(zhuǎn)動(dòng)，但固連于物體坐標(biāo)系中的應(yīng)變分量仍假設(shè)為小量, 即大位移小應(yīng)變情況。10.1.1 幾何非線性問題的牛頓迭代法 由數(shù)值分析技術(shù)可知，求解非線性方程組的數(shù)值方法的常規(guī)方法是Newton-Raphson法，即牛頓迭代法，這是一種近似線性化迭代求解方法。對(duì)于非線性方程，具有一階導(dǎo)數(shù)，在點(diǎn)作一階泰勒級(jí)數(shù)展開，它在點(diǎn)的線性近似為 (4.142)因此，非線性方程在附近似為線性方程： (4.143)當(dāng)時(shí)，由上式求得步的修正項(xiàng) (4.144)Newton-Raphson方法的迭代公式為 (4.

6、145)在幾何非線性有限元法中，結(jié)構(gòu)的剛度矩陣與其幾何位置有關(guān)，平衡方程由變形后的位形描述，因此，結(jié)構(gòu)的剛度矩陣是幾何變形的函數(shù)。設(shè)變形為, 結(jié)構(gòu)的平衡方程式 (4.146)為一個(gè)非線性方程組。記非線性方程 (4.147)用Newton-Raphson方法求的根時(shí)，迭代公式分別為 (4.148)其中, 滿足下式 (4.149)式中, 稱為切線剛度矩陣，表達(dá)式為 (4.150) 在每一個(gè)迭代步中，通過求解切線剛度矩陣，進(jìn)而用進(jìn)行迭代求解，稱為Newton-Raphson方法，又稱切線剛度法。牛頓法的收斂性是好的。但是某些非線性問題中，使用牛頓法迭代時(shí)，若出現(xiàn)奇異或病態(tài)，則對(duì)的求逆出現(xiàn)困難。關(guān)于這

7、一點(diǎn)也可以采用其它修正辦法，如引入阻尼因子。對(duì)于已經(jīng)建立的有限元方程，設(shè)表示內(nèi)為和外力矢量的總和，有 (4.151)式中, R為載荷列陣；為虛位移；為虛應(yīng)變用應(yīng)變的增量形式代入上式，消去項(xiàng)，可以得到非線性問題的一般平衡方程式為 (4.152)該式不論位移或應(yīng)變的大小與否均成立。在有限變形中，應(yīng)變和位移之間的關(guān)系是非線性的，即B矩陣是的非線性函數(shù)。但是，近似地可將進(jìn)行如下分解： (4.153)式中, 為線性應(yīng)變分析的部分； 為由非線性變形引起的，與有關(guān)。假定應(yīng)力應(yīng)變關(guān)系為線彈性，于是有 (4.154)式中 為材料的彈性矩陣； 為初應(yīng)變列陣；為初應(yīng)力列陣對(duì)于式(4.152)的非線性平衡方程式，可用

8、Newton-Raphson方法進(jìn)行迭代求解。對(duì)該式微分，有 (4.155)不考慮初應(yīng)變和初應(yīng)力的影響，得并且這樣可得 (4.156)這里 (4.157)式中為通常的小位移的線性剛度矩陣。矩陣則是由于大位移引起，它可以寫成 (4.158)式(4.156)又可記成： (4.159)式中 (4.160)式中，是關(guān)于應(yīng)力水平的對(duì)稱矩陣，稱之為初應(yīng)力矩陣或幾何剛度矩陣。 因此，用Newton-Raphson方法迭代求解幾何非線性問題的步驟為： (1) 用線彈性解作為，即一次近似； (2) 通過定義求出，求出； (3) 確定切線剛度矩陣；(4) ， ； (5) 重復(fù)上述迭代步驟，直至足夠小。在這里，沒有

9、考慮載荷R可能由于變形而發(fā)生的變化，即在這里假設(shè)了載荷不因變形而改變其大小和方向，否則是非保守力作用下的大變形問題，在此不做討論。10.1.2 典型單元的切線剛度矩陣求解具體的幾何非線性問題時(shí)，必須計(jì)算單元的切線剛度矩陣。對(duì)于一般空間問題，無(wú)論位移和應(yīng)變大小，都可以利用應(yīng)變的基本定義寫出位移和應(yīng)變的關(guān)系式。用變形前的坐標(biāo)做為自變量，可以用位移定義如下大變形問題的應(yīng)變分量表達(dá)式 (4.161) 對(duì)于微小位移情況，可以略去二次以上的偏導(dǎo)數(shù)項(xiàng)，得到小變形時(shí)的應(yīng)變公式。在有限變形中，假設(shè)應(yīng)變?nèi)詾樾×?。?yīng)變和位移之間的關(guān)系為： (4.162)式中為線性應(yīng)變部分。對(duì)于非線性部分，可以寫成： (4.163)

10、式中 (4.164)式中C為矩陣。根據(jù)的定義，可以將表示成任意一點(diǎn)位移的函數(shù)，引入形函數(shù)N后，可以得到 (4.165)對(duì)于(4.163)式進(jìn)行微分，得 (4.166)因此， (4.167)B矩陣為 (4.168)這樣得到 (4.169) 另外，有 (4.170)利用矩陣C和列陣的性質(zhì)，得到 (4.171)式中I為三階單位矩陣，M是的六個(gè)應(yīng)力分量組成的矩陣。因此幾何剛度矩陣為 (4.172)故此，非線性三維單元的切線剛度矩陣為 (4.173)作為特例，可以直接寫出三角形單元的上述有關(guān)表達(dá)式。由三角形單元的位移模式 (4.174)其中，式中的等由結(jié)點(diǎn)坐標(biāo)確定，為三角形單元的面積。 根據(jù)式(4.16

11、4)把式(4.174)代入上式得 (4.175) 由(4.165)式可以知道 (4.176)根據(jù)定義，由式(4.163)確定的平面問題的C矩陣為 (4.177)這樣可以得到C的顯式為(4.178)故 (4.179)而由線性問題給出，即 (4.180) 至此，、和G都是常數(shù)矩陣，只與單元結(jié)點(diǎn)坐標(biāo)和結(jié)點(diǎn)位移有關(guān)。因此，線性剛度矩陣為： (4.181)式中為單元厚度。初始剛度矩陣為： (4.182)幾何剛度矩陣為： (4.183)式中 (4.184)因此，對(duì)于幾何非線性問題，平面三角形單元的切線剛度矩陣可以由式(4.173)求出。10.2 材料非線性問題的有限元法材料非線性是指材料的本構(gòu)方程是非線性

12、的。一般主要分為兩類: 一類是非線性彈性問題，如橡膠、塑料、巖石等，在加載時(shí)應(yīng)力應(yīng)變關(guān)系等性質(zhì)呈現(xiàn)非線性的物理現(xiàn)象，卸載時(shí)可逆。另一類是指材料的彈塑性問題。材料超過屈服極限后呈現(xiàn)出非線性。在機(jī)械結(jié)構(gòu)分析中，常見的本構(gòu)方程主要有線彈性和非線性彈性模型，其特點(diǎn)是應(yīng)力僅應(yīng)變的函數(shù)，加卸載規(guī)律相同。公式如下： (4.185)其中，對(duì)于線彈性材料為常數(shù)，對(duì)于非線性彈性材料，是的函數(shù)。此外，還有超彈性模型、次彈性模型、彈塑性模型等。對(duì)于彈塑性模型，可以認(rèn)為彈塑材料發(fā)生塑性變形時(shí)，其總應(yīng)變可以分解為兩部分： (4.186)即總應(yīng)變?yōu)閺椥詰?yīng)變和塑性應(yīng)變之和。加載時(shí)遵循一定規(guī)律，如Prandtl-Reuss方程

13、，而卸載時(shí)為彈性。對(duì)于應(yīng)力足夠大時(shí)的金屬、土壤、巖石等材料，都有此類特征。在非線性動(dòng)力有限元分析時(shí)，具有黏性特性的模型十分重要。對(duì)于黏彈性模型，一般包括松馳（指突加應(yīng)變作用下應(yīng)力逐漸減少）和蠕變（指突加應(yīng)力作用下應(yīng)變會(huì)逐漸增加）。典型的松馳模型是如下Maxwell模型： (4.187)蠕變模型如Voigt-Kelvin模型為： (4.188)而彈/粘塑性模型是指材料的塑性變形與時(shí)間有關(guān)，本構(gòu)方程中出現(xiàn)非齊次的時(shí)間微分。下面以彈/粘塑性材料的變形分析為例，說明非線性動(dòng)力有限元的基本步驟。即彈/粘塑性材料的變形過程中要考慮時(shí)間效應(yīng)，材料開始屈服后，塑性流動(dòng)、應(yīng)力和應(yīng)變均與時(shí)間有關(guān)。10.2.1 彈

14、/粘塑性問題的基本表達(dá)式 假設(shè)總應(yīng)變分離成彈性應(yīng)變和粘塑性應(yīng)變，即 (4.189)其中表示對(duì)時(shí)間的求導(dǎo)?？倯?yīng)力率取決于彈性應(yīng)變率，有 (4.190)式中，D是彈性矩陣。粘塑性性質(zhì)為 (4.191)式中，是單向屈服應(yīng)力，是硬化參數(shù)的函數(shù)。設(shè)粘塑性應(yīng)變率僅取決于當(dāng)前的應(yīng)力，如下式所示： (4.192)式中，是塑性勢(shì)，是流動(dòng)參數(shù)，項(xiàng)對(duì)于是正單調(diào)增量函數(shù)，即 (4.193)在這里只討論的情況。函數(shù)的兩種常用形式為 (4.194)或 (4.195)式中，為常數(shù)。10.2.2 粘塑性應(yīng)變?cè)隽亢蛻?yīng)力增量 對(duì)于式(4.192)所表示的應(yīng)變定律，定義時(shí)間間隔出現(xiàn)的應(yīng)變?cè)隽繛?，?(4.196)其中為常數(shù)， 。上

15、式中的可用如下近似公式表達(dá) (4.197)式中, . 取決于應(yīng)力水平，可導(dǎo)出顯式公式，具體內(nèi)容參見有關(guān)文獻(xiàn)歐文，塑性有限元。對(duì)于應(yīng)力增量，有 (4.198)用位移增量表示時(shí)，有 (4.199)10.2.3 彈/粘塑性平衡方程 在任何瞬時(shí)的平衡方程應(yīng)滿足 (4.200)在時(shí)間增量中，平衡方程由增量形式給出，即 (4.201)在時(shí)步中出現(xiàn)的位移增量能為 (4.202) (4.203)這里，是切線剛度矩陣，使用如下形式，即 (4.204)把位移增量代回式(4.199)，可以得到應(yīng)力增量. 再有 (4.205) (4.206)且有 (4.207) (4.208)應(yīng)力增量的計(jì)算是基于增量平衡方程(4.2

16、01)的線性化形式，累積所有這樣的應(yīng)力增量得到總應(yīng)力是不正確的，并不會(huì)真實(shí)地滿足平衡方程(4.200)。為此，可以采用計(jì)算殘余平衡力的方法進(jìn)行迭代求解，即 (4.209) 對(duì)于幾何非線性問題，由位移計(jì)算得出。然后，殘余平衡力迭加到下一個(gè)時(shí)間步的載荷增量上。在彈/粘塑性分析中，時(shí)間步長(zhǎng)的選擇十分重要，限于篇幅，這里不再加以討論。10.3 材料非線性問題的動(dòng)力有限元法在材料非線性問題的動(dòng)力有限元分析中，首先考慮材料的非線性本構(gòu)關(guān)系，此外還要考慮是否包含大位移、大應(yīng)變等因素，在實(shí)際工程分析中，需要對(duì)具體不同問題具體處理，目前沒有統(tǒng)一的方法。在這里只限于考慮小應(yīng)變條件下的材料非線性問題，且對(duì)有限元?jiǎng)悠?/p>

17、衡方程進(jìn)行時(shí)間上的數(shù)值積分，以實(shí)現(xiàn)動(dòng)力分析。根據(jù)虛功原理，可以導(dǎo)出結(jié)構(gòu)材料在時(shí)刻的動(dòng)平衡方程。平衡方程與材料本身的性質(zhì)無(wú)關(guān)。對(duì)于每一個(gè)結(jié)點(diǎn)應(yīng)滿足相應(yīng)的平衡方程： (4.211)其中為內(nèi)阻力： (4.212)體力一致力為(為作用的體力矢量) (4.213)式中，為時(shí)刻的形函數(shù)。慣性力為(為密度) (4.214)阻尼力為(為阻尼系數(shù)陣)： (4.215)邊界面力的一致力為(為面力矢量)： (4.216)式中, 為受面力作用的邊界與單元邊界重合的部分。利用Gauss-Legendre乘積方法，依據(jù)單元形函數(shù)，可以對(duì)上述各式進(jìn)行數(shù)值計(jì)算。設(shè)材料的粘彈性本構(gòu)方程中各部分如下： 1) 彈性應(yīng)變 (4.21

18、7)其中為彈性矩陣，設(shè)為常數(shù)矩陣。 2) 粘彈性應(yīng)變?cè)隽空硰椥詰?yīng)變率可以表示為 (4.218)即 (4.219)在時(shí)刻為 (4.220)且有 (4.221)采用如下插分格式 (4.222)其中，為插分系數(shù)，如0，0.5，1等。則粘彈性應(yīng)變?cè)隽颗c應(yīng)力增量之間的關(guān)系為 (4.223) 3)粘塑性應(yīng)變?cè)隽慨?dāng)滿足一定的屈服準(zhǔn)則時(shí)，材料產(chǎn)生粘塑性變形。采用類似的插分方程 (4.224)其中用Taylor級(jí)數(shù)來(lái)近似 (4.225)其中由屈服準(zhǔn)則和粘塑性流動(dòng)準(zhǔn)則確定 (4.226) 采用Ducker-Prager準(zhǔn)則： (4.227)式中由材料粘聚力和內(nèi)摩擦角定；為應(yīng)力第一不變量；為應(yīng)力偏量的第二不變量，當(dāng)

19、取關(guān)聯(lián)的流動(dòng)法則時(shí)有 (4.228)這樣，由式(4.224)和式(4.225)得到 (4.229) 4)總應(yīng)變?cè)隽坑蓱?yīng)變分解可知 (4.230)其中總應(yīng)變?cè)隽靠梢杂山Y(jié)點(diǎn)位移得到 (4.231)式中, 為位移應(yīng)變矩陣；為單元結(jié)點(diǎn)位移增量。由式(4.223)、式(4.229)和式 (4.230)可得應(yīng)力增量為： (4.232)是粘彈粘塑性模量 (4.233)5) 動(dòng)態(tài)有限元?jiǎng)偠染仃囋谌我粫r(shí)刻，可以證明，該時(shí)刻的切線剛度矩陣為 (4.234) 3. 顯式時(shí)間積分法任一時(shí)刻的動(dòng)平衡方程可以寫成如下矩陣形式： (4.235)式中為總質(zhì)量矩陣，為總阻尼矩陣，為內(nèi)力的總矢量，為作用的體力、面力的一致結(jié)點(diǎn)力矢量。為簡(jiǎn)便起見，用中心差分公式對(duì)上述動(dòng)平衡方程進(jìn)行離散化。加速度、速度分別為： (4.236)