材料非線性有限元法

1、第四章 材料非線性有限元法 以上三章分別研究了線性彈性有限元法，材料非線性本構(gòu)方程和非線性方程組解法，本章就可以研究材料非線性有限元法了。 在材料非線性基本方程中，除第二章所述的本構(gòu)方程外，與線性彈性一樣，而非線性有限元法又歸結(jié)為一系列線性彈性問題。因此，只要在第一章中改用第二章的本構(gòu)方程，就可建立材料非線性有限元法的基本內(nèi)容。§4-1 非線性彈性有限元法 第二章提到，非線性彈性本構(gòu)方程與形變理論彈塑性本構(gòu)方程在形式上相同，所以與第二章一樣，這里也按塑性力學(xué)形變理論，研究非線性彈性有限元法，以便把二者統(tǒng)一起來。 1 非線性彈性基本方程 為了便于以后直接引用，這里列出全量形式的非線性彈

2、性（或形變理論彈塑性）基本方程，并用矩陣表示。 幾何方程： (1.14) 本構(gòu)方程：=D (2.13) 平衡方程：(在內(nèi)) (1.20) 邊界條件： （在A上） （1.22） （在A上） (1.23) 虛功方程： (1.28) 位能變分方程：=0 (1.31)其中 (1.32) (4.1) 2 非線性方程組的建立 由于虛功方程本身不涉及材料性質(zhì)，所以第一章由虛功方程得到的單元平衡方程（1.48）式和總體平衡方程（1.109）式完全適用于非線性彈性（或形變理論彈塑性）問題?？梢?，只要把非線性彈性（或形變理論彈塑性）本構(gòu)方程代入單元或總體的平衡方程，就可以建立非線性方程組。 （1）割線剛度方程 仿

3、照線性彈性有限元法，把（1.36）式代入（2.13）式后，再把（2.13）代入（1.48）式便得單元割線剛度方程，即 (4.2)其中單元割線剛度矩陣 (4.3)而割線本構(gòu)矩陣，如（2.14）式所示。 仿照（1.113）式的推導(dǎo)，同樣可得總體割線剛度方程即 (4.4)其中總體割線剛度矩陣 (4.5)而總體節(jié)點(diǎn)載荷P仍如（1.110）式所示。 由（4.5）式可知，總體割線剛度矩陣K取決于各單元的等效應(yīng)變 ；又由（2.5）式可知，等效應(yīng)變 是由應(yīng)變計(jì)算出來的；再由（1.36）和（1.106）式可知，應(yīng)變與總體節(jié)點(diǎn)位移U有關(guān)?？梢?，總體割線剛度矩陣K是總體節(jié)點(diǎn)位移U的函數(shù)，所以總體割線剛度方程（4.4

4、）式是一個(gè)非線性方程組。 必須指出，建立非線性方程組（4.4）式，只是為了說明非線性彈性（或形變理論彈塑性）有限元方程的非線性性質(zhì)。實(shí)際求解時(shí)并不用（4.4）式。因?yàn)榍蠼猓?.4）式要用直接迭代法，而正如 3-2指出，直接迭代法不但計(jì)算量太大，而且常常不收斂。 （2）切線剛度矩陣 由 3-2-3-6可知，在求解非線性方程組時(shí)，除上述直接迭代法外，都要用到切線剛度矩陣（至少要用到初始切線剛度矩陣k和K）。為此，這里討論一下建立非線性彈性（或形變理論彈塑性）有限元方程中的切線剛度矩陣問題。 由（1.109）和（3.11）式可知 (4.6)于是由（3.10）和（4.6）式可得 (4.7)由于由（2.

5、16），（1.36）和（1.106）式，并考慮到符號d和d分別是d和d，有 (4.8) (4.9) (4.10)所以把（4.8）-（4.10）式代入式（4.7）式便得總體割線剛度矩陣，即 (4.11)其中單元切線剛度矩陣 (4.12) （3）具有初應(yīng)變理論或初應(yīng)力的剛度方程 仿照線性彈性有限元法，把形式上相同的（3.101）式代入（2.13）式，并令=0或=0，再把（2.13）式代入（1.48）式便得單元?jiǎng)偠确匠蹋?(4.13)或 (4.14)其中單元?jiǎng)偠染仃嚭统鯌?yīng)變，初應(yīng)力節(jié)點(diǎn)載荷，仍分別如（1.50）和（1.53）、（1.54）式所示。但要強(qiáng)調(diào)，這里k的含義是單元初始切線剛度矩陣；中的初

6、應(yīng)變或中的初應(yīng)力隨迭代過程而變。 仿照線性彈性有限元法，同樣可得總體剛度方程，即 (4.15)或 (4.16)其中總體剛度矩陣和總體初應(yīng)變、初應(yīng)力節(jié)點(diǎn)載荷、在形式上均與線性彈性有限元法相同。 3 等效應(yīng)力、等效應(yīng)變關(guān)系 由（4.11）-（4.16）式可知，要建立并求解非線性彈性（或形變理論彈塑性）有限元方程，關(guān)鍵是要具體知道材料的本構(gòu)矩陣。而由（2.14）和（2.18）式可知，只要（2.15）和（2.19）式中的函數(shù)關(guān)系是已知的，那么本構(gòu)矩陣就是顯式的。 根據(jù)單一曲線假設(shè)，和的關(guān)系與單向拉伸時(shí)相同，即 (4.17)再考慮體積不可壓縮條件（），則 (4.18)其中取決于所采用的簡化模型。 理想塑

7、性（見圖4-1）： (4.19)線性強(qiáng)化塑性（見圖4-2）： (4.20)冪次強(qiáng)化塑性（見圖4-3）：， （4.21）4 迭代公式的具體化 由于非線性彈性（或形變理論彈塑性）有限元方程一般都寫成全量形式，所以這里只相應(yīng)的列出幾種迭代類型解法的具體迭代公式。（1）Newton-Raphson法 由（1.36）、（1.106）和（2.5）式以及（3.17）和（3.18）式，有 （4.22） （4.23） （4.24） (4.25) （4.26） （4.27）（2）初應(yīng)變迭代法 由（2.10）、（2.13）和（3.99）、（3.101）式可知 （4.28） （4.29） （4.30）所以仿照Newt

8、on-Raphson法，并考慮到（3.109）和（3.110）式，有 （4.31） （4.32） （4.33） （4.34） （4.35）（3）初應(yīng)力迭代法 由（2.10）、（2.13）和（3.118）、（3.120）式可知 （4.36） （4.37） （4.38）所以仿照Newton-Raphson法，并考慮到（3.126）和（3.127）式，有 （4.39） (4.40) (4.41) (4.42) (4.43)§4-2 非線性彈性手算例題為了熟悉非線性彈性(或形變理論彈塑性)有限元法及其非線性方程組的求解過程,這里以圖4-4(a)所示的彈塑性拉壓超靜定問題為例，用Newton-

9、Raphson法、初應(yīng)變迭代法、初應(yīng)力迭代法進(jìn)行手算。其中用Newton-Raphson法的求解作較詳細(xì)的敘述，以便了解非線性彈性(或形變理論彈塑性)有限元分析的全過程，而用其他方法的求解只給出主要計(jì)算過程和計(jì)算結(jié)果。在該拉壓超靜定桿中劃分的節(jié)點(diǎn)和單元如圖4-4（a）中所示。單元和分別由線性強(qiáng)化材料和線性彈性材料制成，如圖4-4（b）和（c）所示。兩個(gè)單元的截面積均為，長度均為，彈性模量均為，單元的強(qiáng)化模量為。節(jié)點(diǎn)2所受集中載荷，其中為單元的屈服極限。 1 2 3 圖4-4（a）1 用Newton-Raphson法求解（1）非線性有限元方程的形成 首先去掉兩固定端約束，用其約束反力和代替，所以

10、 （a）由于 （b）所以 （c） （d）設(shè)單元內(nèi)任意一點(diǎn)位移為而所以其中幾何矩陣 （e）單元節(jié)點(diǎn)位移向量這里的單元都是單向應(yīng)力狀態(tài)，即 （f） （g）所以由（4.24）、（4.25）和（4.26）式可得 （h） （i）把（c）（g）式代入（h）和（i）式，并考慮到，有 （j） （h）把（a）、（j）和（k）式代入（3.17）式，并考慮到，有 （）最后由邊界條件和，用第一章所述的消行降階法便得約束處理后的非線性有限元方程，即 （m）（2）切線剛度迭代公式的建立 由于所以無量綱迭代公式為 （n） （o） （p） （q） （r） （s）（3）節(jié)點(diǎn)位移和單元應(yīng)力的求解 按迭代公式（n）（s）式，求節(jié)點(diǎn)

11、位移和單元應(yīng)力的迭代過程和計(jì)算結(jié)果如表4-1所示： 表4-1n101.5001.5001.250-1.50022.7500.1661.6661.333-1.66633.00002 用初應(yīng)變迭代法求解 若令（j）式中，則 （A）仿照用Newton-Raphson法的求解，由（3.103）式有 （B）約束處理后可得 （C）因?yàn)橹挥袉卧哂兴苄?，所以初?yīng)變節(jié)點(diǎn)載荷只有由單元提供。由于所以 （D）由圖4-4（b）可知因此 （E）于是由（n）、（C）（E）和（3.109）式可得無量綱迭代公式，即 （F） （G） (H) （I） （J）按迭代公式（F）（J）式，求節(jié)點(diǎn)位移和單元應(yīng)力的迭代過程和計(jì)算結(jié)果如表

12、4-2所示： 表4-2n1001.5001.500-1.50020.5000.2501.7501.250-1.75030.2500.1251.6251.375-1.62540.3750.1881.6881.313-1.68750.3130.1561.6561.343-1.65660.3430.1721.6721.329-1.67270.3290.1641.6641.335-1.66480.3350.1681.6681.333-1.66890.3330.1661.6661.333-1.666100.3333 用初應(yīng)力迭代法求解 仿照用初應(yīng)變迭代法求解，有 （K）約束處理后可得 （L）由于只有單元

13、具有塑性，所以初應(yīng)力節(jié)點(diǎn)載荷只由單元提供。由于所以 （M）于是由（n）、（L）、（M）和（3.126）式可得無量綱迭代公式，即(N) （O） （P） （Q）按迭代公式（N）（Q）式，求節(jié)點(diǎn)位移和單元應(yīng)力的迭代過程和計(jì)算結(jié)果如表4-3所示： 表4-3n101.5001.250-1.50020.1251.6251.312-1.62530.1561.6561.328-1.65640.1641.6641.332-1.66450.1601.6671.333-1.66760.166§4-3 彈塑性有限元法在§4-1中討論了基于形變理論的彈塑性有限元法。這里將討論基于流動(dòng)理論的彈塑性有限

14、元法。為了跟蹤加載歷史求出位移、應(yīng)變和應(yīng)力的全量，基于流動(dòng)理論彈塑性有限元方程只能取增量形式。同非線性彈性（或形變理論彈塑性）問題一樣，第一章由虛功方程得到的平衡方程也完全適用于流動(dòng)理論的彈塑性問題。所以由（1.109）、（3.11）和（3.59）式便得總體增量剛度方程，即 （4.44）其中總體切線剛度矩陣 （參考（4.11）式） （4.45）單元切線剛度矩陣 （參考（4.12）式） （4.46）（4.46）式中的切線彈塑性矩陣，按（2.63）式（正則屈服面）或（2.72）和（2.88）式（非正則屈服面）計(jì)算。幾種常見材料模型的切線彈塑性矩陣已在§2-3 §2-5中給出顯示

15、形式。由于流動(dòng)理論彈塑性有限元方程只能取增量形式，所以在第三章所述的解法中，只有增量類型的解法才能用來對它進(jìn)行求解，其中常用的是Euler-Cauchy法、Euler一次迭代法、初應(yīng)變增量法和初應(yīng)力增量法。1 逐步求解公式的具體化 這里列出一求得的第I增量步的具體公式。（1）Euler-Cauchy法 對比（3.59）、（3.61）式和（4.44）（4.46）式，并仿照非線性彈塑性有限元法，Euler-Cauchy法逐步求解公式為： （4.47） （4.48） （4.49） （4.50）Euler一次迭代法：同樣，仿照非線性彈塑性有限元法，Euler一次迭代法逐步求解公式為： （4.51） （

16、4.52） （4.53） （4.54） （4.55）（2）初應(yīng)變增量法 仿照§4-1中的初應(yīng)變迭代法，并考慮到（3.102）和（3.114）（3.116）式，有 （4.56） （4.57） （4.58） （4.59）對于Mises模型，若把塑性應(yīng)變增量直接看作是初應(yīng)變增量，并寫成有限增量形式，則由（2.126）、（2.135）和（2.136）式可得： （4.60）初應(yīng)變增量法：仿照§4-1中的初應(yīng)力迭代法，并考慮到（3.121）和（3.135）（3.137）式，有： （4.61） （4.62） （4.63） （4.64）2 Mises模型和Drucher-Prager模型

17、（1）單元彈、塑性狀態(tài)的確定 由于同一增量步的不同單元可能處于不同的狀態(tài)彈性或塑性，而不同狀態(tài)下單元的應(yīng)力增量應(yīng)按不同公式計(jì)算，所以要計(jì)算（4.47）（4.64）式中的單元應(yīng)力增量，必須參照§2-2中的加卸載準(zhǔn)則，確定第I增量步單元的彈、塑性狀態(tài)。為了確定第I增量步單元的狀態(tài)，可先假定該單元在此增量步內(nèi)不產(chǎn)生新的塑性變形（彈性、卸載或中性變載），于是有： （4.65）據(jù)此可以計(jì)算第i步終了時(shí)的加載函數(shù)值或屈服函數(shù)值。為了計(jì)算方便，可把加載函數(shù)寫成另一種形式。由（2.102）、（2.148）和（2.165）式可知，Mises等向強(qiáng)化加載函數(shù)值為： （4.66）Mises隨動(dòng)強(qiáng)化（線性）

18、加載函數(shù)值為： （4.67）Drucher-Prager理想塑性加載函數(shù)（即屈服函數(shù)）值為： （4.68）參照§2-2中的加卸載準(zhǔn)則，若則應(yīng)力在加載面內(nèi)（彈性或卸載）或在加載面上（中性變載），說明上述假定是正確的。我們把這種單元簡稱為彈性單元。若則應(yīng)力在加載面外，這當(dāng)然是不可能的，說明上述假定是錯(cuò)誤的。但這種情況下的單元只有兩種可能狀態(tài)：全塑性或彈塑性，而單元的確切狀態(tài)又取決于第i-1步終了時(shí)（或第i步開始時(shí)）的單元狀態(tài)。若第i-1步終了時(shí)單元是塑性的，即則第i步單元是全塑性的（有一種塑性狀態(tài)進(jìn)入另一種塑性狀態(tài)，即完全加載狀態(tài)），我們把這種單元簡稱為塑性單元。若第i-1步終了時(shí)單元是

19、彈性的，即則第i步單元是彈塑性的（有彈性狀態(tài)進(jìn)入塑性狀態(tài)，即部分加載狀態(tài)），我們把這種單元簡稱為彈塑性單元（或過度單元）。（2）單元應(yīng)力增量的計(jì)算 對于彈性單元，顯然有： （4.69）對于塑性單元，當(dāng)單元應(yīng)變增量較小時(shí)，有： （4.70）當(dāng)單元應(yīng)變增量較大時(shí)，單元應(yīng)力增量應(yīng)當(dāng)用切線模量按增量法求下列積分的近似值： （4.71）對于彈塑性單元，單元應(yīng)力增量由彈性和塑性兩部分組成。若用r表示彈性應(yīng)力增量與總應(yīng)力增量之比參看圖4-5（a）和相應(yīng)的一維問題圖4-5（b）,則可通過滿足加載條件確定這個(gè)比例因子。確定出r以后，當(dāng)單元應(yīng)變增量較小時(shí)，有 （4.72） （4.73）當(dāng)單元應(yīng)變增量較大時(shí)，單元應(yīng)

20、力增量應(yīng)當(dāng)用切線模量按增量法求下列積分的近似值： （4.74）顯然，當(dāng)r=0時(shí)，此式即為（4.71）式。（3）比例因子r的確定 假定用或表示彈性應(yīng)力偏張量，于是把（4.66）（4.68）式中的或換 2 ri-1 i-1 10 (a) i-1 ri-1 ri-1 i-1(b)圖45為或，并令其值為零便得 （4.75）其中A、B和C，Mises模型等向強(qiáng)化材料為（4.76）（4.77）（4.78）Mises模型隨動(dòng)強(qiáng)化（線化）材料為 （4.79）（4.80）（4.81）Drucker-Prager模型為理想塑性材料為（4.82）（4.83）（4.84）（4）積分近似值的計(jì)算（4.71）和（4.74

21、）式中的積分很難積得顯式。為此，我們可以把積妥對應(yīng)的塑性變形部分再分為若干子增量，用分段線性的計(jì)算結(jié)果去逼近積分結(jié)果。由于（4.71）式是（4.74）的r=0的特殊情況，所以這里只討論（4.74）式中積分的近似計(jì)算問題。若把，再分為N個(gè)相等的子增量（例如圖46中N5），則有（4.85）123456(1-r)i-1 圖46由于初始子增量步開始時(shí)的單元應(yīng)力為（4.86）第j子增量步開始時(shí)的單元應(yīng)力為（4.87）所以（4.88）（5）第i增量步的單元計(jì)算過程在第I增量步，已知的是和，要計(jì)算的是和。若用表示判定單元應(yīng)變增量大小的標(biāo)準(zhǔn)，則計(jì)算過程如圖47所示。3Tresca模型由于Tresca模型的屈服

22、面比較復(fù)雜，所以這里只限于平面應(yīng)力問題的Euler-Cauchy解法。由§25中的和可知，Tresca模型彈塑性矩陣各元素均與主應(yīng)力方向有關(guān)。當(dāng)單元開始進(jìn)入塑性狀態(tài)以后，主方向隨載荷的增加而改變，因而每一增量步都要重新對它進(jìn)行計(jì)算。為了使應(yīng)力滿足加載條件，可以按下述插值方法對它進(jìn)行計(jì)算第i增量步的主方向值。（1）按計(jì)算主應(yīng)力和塑性功用（塑性單元）和（彈性單元）以及它們的組合（彈塑性單元）形成，求得位移增量以后，有 （4.89）（4.90）（4.91）（4.92）（4.93）（4.94）（4.95）（2）按計(jì)算主應(yīng)力和塑性功仿照（1）中的計(jì)算，有（4.96）（4.97）（4.98）（4

23、.99）（4.100）（4.101）（3）按插值方法計(jì)算主方向由、和按一種插值方法求，使得用它求得的應(yīng)力滿足加載條件。對于圖26中的邊6，有（4.102）§44彈塑性手算例題我們?nèi)砸?#167;42中的彈塑性拉壓問題為例，用Euler-Cauchy法、 Euler一次迭代法、初應(yīng)變增量法和初應(yīng)力增量法進(jìn)行手算。為了反映應(yīng)力與加載歷史的相關(guān)性，除Euler-Cauchy法外，均用敘述方式。（1）用Euler-Cauchy法求解仿照§42中用Newton-Raphson法的求解，有（R）約束處理后可得（S）若以總載荷為參考載荷，并按兩個(gè)增量步取，則有于是無量綱逐步求解公式為（T） >0返