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

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

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

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

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

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

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

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

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

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

10、式中 (4.164)式中C為矩陣。根據(jù)的定義，可以將表示成任意一點位移的函數(shù)，引入形函數(shù)N后，可以得到 (4.165)對于(4.163)式進行微分，得 (4.166)因此， (4.167)B矩陣為 (4.168)這樣得到 (4.169) 另外，有 (4.170)利用矩陣C和列陣的性質(zhì)，得到 (4.171)式中I為三階單位矩陣，M是的六個應(yīng)力分量組成的矩陣。因此幾何剛度矩陣為 (4.172)故此，非線性三維單元的切線剛度矩陣為 (4.173)作為特例，可以直接寫出三角形單元的上述有關(guān)表達式。由三角形單元的位移模式 (4.174)其中，式中的等由結(jié)點坐標確定，為三角形單元的面積。 根據(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é)點坐標和結(jié)點位移有關(guān)。因此，線性剛度矩陣為： (4.181)式中為單元厚度。初始剛度矩陣為： (4.182)幾何剛度矩陣為： (4.183)式中 (4.184)因此，對于幾何非線性問題，平面三角形單元的切線剛度矩陣可以由式(4.173)求出。10.2 材料非線性問題的有限元法材料非線性是指材料的本構(gòu)方程是非線性

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

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

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

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

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

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

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

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