(完整word版)材科非線性有限元法

1、43第四章材科非線性有限元法本章所研究的非線性有限元方法是材料非線性有限元。它的非線性是由本構關系的非線性引起的。但它和線彈性有限元一樣，都屬于小變形問題，因而關于形函數(shù)的選取、應變矩 陣、應力矩陣及勁度矩陣的形式都是相同的，不同的僅在勁度矩陣是按非線性彈性或彈塑性矩陣計算的，這是材料非線性有限元的基本內(nèi)容。4.1 非線性彈性有限元法4.1.1 非線性彈性有限元法的基本公式對于小變形的非線性彈性問題，其基本方程為：平衡微分方程LTbp = 0(3-8)幾何方程=Lf(3-2)物理方程2Ds(3-19)或db =DTde(3-20)式中，Ds和DT的表達式見（3-37）（ 3-40），它們是應變

2、強度I的函數(shù)。邊界條件為位移邊界條件：f虛功方程是應力邊界條件：n（二）s=s._= p(3-9)TT_T _i godv二fTpdv亠i f pdsvvs-(3-11)可見，除表示為物理方程的本構關系外，其它基本方程與邊界條件均與線彈性問題相同,且虛功方程（3-11）不涉及材料性質(zhì)，因而線彈性有限元的幾何關系和由虛功方程得到的單元和整體平衡方程完全適用于非線性彈住和彈塑性問題，它們是f = N3e(4-1)(4-2)BdV = Rev(4-3)寸寸e TT(c) B odV = Rve(4-4)將式（4-1）代入（3-2），式（3-2）代入（3-19），再將（3-19）代入（4-3），得ks

3、（）八Re(4-5)其中單元割線(3-4)44由于Ks與位移：有關，式（4-7）是一個非線性方程組。但實際求解時不用（4-7）式,因為求解（4-7）要用直接迭代法，這一方法不但計算量大，且常常不收斂。在求解非線性方程組時， 除直接迭代法要用割線勁度矩陣外，其它方法都要計算切線勁度矩陣。為此，須對非線性彈性有限元的切線勁度矩陣進行研究。由式（4-4）得TT=送（c ） B cdV R（4-9）e由上式及 Newt on 法迭代公式（2-7）可得切線勁度矩陣T八(ce).vBeTT（ce） （ B DTBdV）ce（4-10）e迭代過程中的一些具體處理方法也不相同，但對荷載的1.荷載分級對于實際的

4、工程問題來說，將荷載分級施加是相當重要的，首先，這樣做可以模擬實際施工加載過程，進行所謂”仿真的數(shù)值分析，這時彈塑性問題顯得更為重要，因為不同的加 載過程將得到不同的位移和應力計算成果。其次，對荷載分級，容易使求解過程收斂。荷載的分級首先要考慮荷載的性質(zhì)，不同類型的荷載要根據(jù)實際施工加載情況分級施加。為了更精細地模擬施工過程并使迭代過程收斂，每級荷載還可按2.2 介紹的荷載系數(shù)法分成更小的增量。 圖 4-1 為三峽永久船閘的閘室橫剖面，為了分析閘室結(jié)構的應力狀態(tài)和兩側(cè)邊坡的穩(wěn)定性，所施加的荷載按加載順 序有開挖荷載、襯砌自重、邊坡內(nèi)的滲壓和 閘室內(nèi)的水壓力。由于開挖荷載的計算需要 事先了解地應

5、力，在沒有地應力量測值的情 況下，還要首先計算自重應力。很明顯，對 于這個問題，開挖是主要荷載，根據(jù)設計的 開挖施工程序再將它分為若干級增量，其它荷載也可根據(jù)需要再分級，目的就是為了更 精確地模ks（）二TB DsBdV同樣可由（4-4 ）得整體平衡方程整體割線勁度矩陣Ks(S)注Re TeKs(S八(c ) ksC)ce(4-6)(4-7)(4-8)KT(S)4.1.2 求解的選代過程 不同的非線性方程組求解方法, 處理一般都分級按增量方法計算。圖 4-145擬加載過程和使求解過程收斂。2增量迭代法增量迭代法實際上就是非線性方程組求解的Euler-Newt on 法，即將荷載分成若干級增量，

6、對每一級荷載增量，進行迭代運算。對第m 級荷載，迭代的基本公式是血=(KT,m)(Rm一F；)=(KRm- 此)*i i(2-46)兇十=審+山式中Rm一一第 m 級荷載增量施加后的總荷載，Rm第 m 級荷載增量，F(xiàn)m第 m 級荷載第 i 次迭代結(jié)束時的結(jié)點力，根據(jù)當時的單元應力二m按下式計算T TFm=無(ce) BolrndV(4-11)vm-第 m 級荷載第 i 次迭代的不平衡力m=FmRmA切線勁度矩陣KT ,m根據(jù)式(4-9)由下式求出TTiKT,Z (ce) (fB DT(寫B(tài)dV)ce(4-12)Lve如果已知第 m 級荷載增量時第 i 次迭代的近似解決，相應的應變?yōu)?Bc e-

7、；m(4-佝于是切線彈性矩陣DT(4)可以確定，由(4-12 )和(4-11 )分別計算(KTX和Fm，最后利用式(2-46)求出=臨，這個迭代過程由 i=0 開始直至=臨足夠小，達到一定計算精度內(nèi)為零時終止。假定最后一次迭代i-l，則該級荷載增貴的最終位移為I3m二3m亠一 詁i =0(4-14)按(2-46)求出位移增量3m后，應力增量為m二 DT(4)Bce3m(4-15)第 m 級荷載增量的最終應力為Im怖-i =0(4-16)一維問題的迭代過程見圖4-2。可以看出，上一級荷載增量作用下，迭代收斂后仍然存在的不平衡力將轉(zhuǎn)至下一級。2. 3 初應力法46如果在彈性材料內(nèi)確實存在初應力衍，

8、則材料的應力應變關系為D十衍(4-17)47由上式及虛功原理可導出單元的結(jié)點力為eeTF = k3+ B%dV集合單元得出以下的有限元支配方程式中，R0為由初應力引起的等效結(jié)點荷載(4-18(4-19)圖 4-2eT TR；十八(ce)vBOodV(4-20)對于非線性問題，在求解非線性方程組時，可利用式（4-19）。此時，非線性的應力應變關系仍由式（4-17）表示，但并非真正的初應力，它是實際應力b和彈性應力 /的差，是一種虛擬的初應力。從圖 4-3 可以看出，二e=D；為彈性應力，殆是由于應 力應變非線性而降低的應力（衍為負值）。初應力法就是將初應 力b看作是變化的，以此來反映應力和應變之

9、間的非線性關系。 通過不斷地調(diào)整初應力，使線彈性解逼近非線性解。前面所談的增量迭代法，在迭代的每一步都要重新形成整體 的切線勁度矩陣K；,m并完整地求解一次線性方程組。如果將式（2-46 ）中的KT,m改用初始的切線勁度矩陣KT,m或KT,1,就得到常勁度的迭代公式。這相當于求解非線性方程的 Euler 修正 Newt on 法，（KT）0是加載一開始時的切線勁度短陣，也 就是彈性勁度矩陣。K；m可以認為是第 m 級荷載增量開始施加時的“彈性”勁度矩陣式中的“彈性”矩陣D實際上是DT（蠱）。式（2-46）可寫成48其中(4-24)簡稱彈性應力。將式（4-22）的第二式代入第一式，并利用等式（4

10、-23）及（4-11）可得與（4-20）相比，非線性求解用的虛擬初應力為（=血 一（/）m式（4-25）是直接計算總位移的迭代公式，一維的迭代過程見圖4-3。圖 4-3設經(jīng)過 M 次迭代后在容許的精度范圍內(nèi)達到陷上器，則迭代結(jié)束，并由- 3m，按真實的本構關系求出最后的應力結(jié)果。在迭代過程中，om和（膚）m分別是由注按真實本構 關系和”線彈性計算出來的應力。（RZO）m在迭代過程中并不趨于零，而是趨于一常數(shù)矢量（圖 4-3 （ b）中的 M ）,它是與收斂解的真實應力=bm和彈性應力（be）M之差相應的等效結(jié)點力。（匚0）m =；：m是迭代結(jié)束時的虛擬初應力。 因而，迭代過程也相當于尋求一個合

11、適的”初應力場。由式（4-25）可得第（i-1 ）次迭代的公式(Ro)將它與（4-25）相減得(4-27)根據(jù)式（4-21）可得出下面的恒等式Km(ce)T B(/)mdv(4-22)(4-23)Km陽二 Rm（Ro）m（7=-遲0廬T（om-（e）m）dv(4-24)(4-26)0mRm49這是求解位移增量的迭代公式。從圖 4-3 （b）可見，在進行第二次迭代（即 i=1 ）時，若按式（4-25 ）求總位移，方程 右端等于M2A，若按式（4-27）求位移增量，方程右端為M,M,二M2M2，對應于1 IP圖 4-3（ a）中的（s）m（線段NN）。第三次迭代（i=2 ）時，式（4-25）右端為

12、M3A，F(xiàn)fFFFF而式（4-27）右端為（RJ；（Rb）m=M2M2-M,M,= M2M2，與（a）中的N2N2所示的應力對應，經(jīng) M 次迭代后，式（4-27）右端接近于零，此時，也趨于零。本級 的最后位移為Mm二鵡=也.盡（4-28）i=0實際上，式（4-27 ）右端就相當于失衡力。4.2 彈塑性有限單元法彈塑性問題屬于小變形，因而上節(jié)關于非線性彈性的一套基本公式及迭代運算方法仍 然適用于彈塑性有限元，但有兩點必須強調(diào)。第一是應力和變形呈非線性關系，勁度矩陣與應變有關， 是變化的，本構方程必須用增量形式表示，非線性方程組的求解必須用增量方法，這與非線性彈性問題是類似的。第二是彈塑性問題中的

13、材料本構關系與應力、變形的歷史有關，需要利用加卸載準則判斷單元的材料究竟處于彈性狀態(tài)還是塑性狀態(tài)，以決定由應變計算應力時采用何種本構關 系。在按照上節(jié)介紹的增量迭代法或初應力法求出位移增量.: S 后，應變增量由下式計算二& = B=序二Be乜3（4-29）由應變增量厶確定應力增量 =C時，需要知道厶&和厶C之間的關系。彈塑性材料的本構方 程為這是以無限小增量de和d&的形式給出的。而在有限元數(shù)值計算中，由于荷載增量是以有限大小的形式R給出，所以，應力增量和應變增量都以有限大小的形式給出，設分別是厶e和也,可以利用數(shù)值積分的方法從式（4-30）得到 &和也e之間

14、的關系,.m丄七#心e=JDepd（4-31）m丄在荷載增量R作用前后，所考慮單元的高斯積分點上的介質(zhì)究竟處于什么狀態(tài)，可以用第三章介紹的應變空間內(nèi)的加卸載準則來判斷。dDepd廠(D - Dp)d(4-30)50 0加載幾何上就是當彈性應力增量d膚指向屈服面外側(cè)，使f ( c + d ,K) 0時為加載，材料處于塑性狀態(tài)。當彈性應力增量d膚指向屈服面內(nèi)側(cè)或與屈服面相切，使ef (bd ,K) 0(4-34)時為卸載和中性變載， 介質(zhì)處于彈性狀態(tài)。 值得注意的是，這兒用于判斷的應力增量d /是彈性應力增量，真正的應力增量de還未求出，當d&已知時，很容易由彈性本構關系求出d/，因而實際

15、上是以應變增量d&來判斷的。在有限元方法中，根據(jù)式(4-33)和(4-34),可由應變增量 ：&來判斷材料所處的狀態(tài)。設某個單元的某個積分點在荷載Rm作用下的應力為OmJ，內(nèi)變量為-mJ，施加荷載增量Rm后，產(chǎn)生的應變增量為A,彈性應力增量A = DA&。如果fm =f(mA一e,m)-0表示應力emj ： L ee在屈服面內(nèi)(彈性或卸載)或在屈服面上(中性變載)，沒有新的塑性應變產(chǎn)生，因而該點處于彈性狀態(tài)，可采用彈性本構關系由&求厶O-,即匚e =D=(4-35)如果fm4 =f(怖,54)=0fm=f(em 4e mJ)0說明該點始終處于塑性狀態(tài)，產(chǎn)生了新的

16、塑性應變。當應變增量厶&較小時，有e = Dep (4-36)當 ： &較大時，須按(4-31)式計算應力增量.：e。如果fm 4-f(em 4 , m 4)：0(4-32)(4-33)51fm= f(Om厶eSmJ0即從荷載增量施加前的彈性狀態(tài)進入施加后的塑性狀態(tài)，此時，可由條件f(emr ee,mJ =0(4-37)52確定應力增量中的彈性部分與總應力增量A c之比 r （0r1，見圖 4-4 （a）。也可以通過對屈服函數(shù)f采用線性內(nèi)插得到cm= 口m_1 +rDA也c=D + （1 _r）Dep（Cm）氐（4-40）=DA_（1_r）Dp（Cm）A從式（4-38）可知，只

17、要取r =1，式（4-39）和（4-40）都成為（4-35）式，表明增量R施加前后，材料狀態(tài)的變化是彈性到彈性，或塑性卸載或中性變載（圖4-4 （ b）,該點介質(zhì)的反應是純彈性的。如果取r =0,式（4-39）和（4-40）成為式（4-36）,表明材料處于塑性加載狀態(tài)（圖 4-4（ c）。當二較大，須采用式（4-31 ）或（4-39）計算；_c時，由于被積分的Dep或Dp很難寫 成顯式，因而可以把與積分對應的塑性變形部分再分為若干子增量，用分段線性的計算結(jié)果去逼近積分值，例如將塑性變形部分（1 -r）： 分為 N 個相等的子增量.：.： 即昶汙（i N則式（4-39）成為Nc = rD也 十瓦

18、Dep(C)j=(rD f bep(Nj=t在求出r之后，根據(jù)式（當丄很小時（要求圖 4-44-31 )可求出應力增量務i-r.一Dd丄m 1也 rD-rDepd(4-38)ir ,Dpd(4-39)R足夠?。?，上式可改寫為下面的近似公式T-(4-41)53整個計算過程見圖4-5。0由7X%：-1 +rzW* * g-J 0型云牛-確定r幾=%+2也&= OD+U幾1幾） 込j +&返544.3 特殊單元在非線性有限元中的應用在非線性有限元分析中，為了模擬材料的某些特殊性質(zhì)或結(jié)構上的一些特征，經(jīng)常要使用一些特殊單元。4.3.1 不抗拉單元對于由不抗拉材料，如破碎巖體、土體，可采

19、用不抗拉單元。此時要求出單元的應力，與材料內(nèi)的初應力迭加，再計算主應力。并檢驗是否出現(xiàn)拉應力，對出現(xiàn)的主拉應力需要進行應力轉(zhuǎn)移”。例如當二1.0時，認為二i不可能大于零，超過部分應轉(zhuǎn)移至周圍其它單元， 并按應力與主應力的關系，確定實際存在的應力再按F 二 (ce)TBTodVve求出結(jié)點力，進行迭代計算。類似地，可以定義低抗拉單元。4.3.2 層狀單元若巖體中有一組沒有膠結(jié)的平行層面，那么這組層面就決定了巖體不能承受拉力的方 向，因而方向明確。設層面的產(chǎn)狀為傾向，傾角一：，現(xiàn)以北方向N為整體坐標系的x正向，E 為y正向，按右手系，z向下為正。并令局部坐標系的z正向與層面的向上法線一致，x,y位

20、于層面內(nèi)。則z的方向余弦為l二cos：sin :m = sin sin：(4-42)n = - cos :由于這類介質(zhì)有強烈的各向異性，破壞的形式只可能是沿層面的剪切滑移或垂直層面的拉開。由有限元方法求出 b 后，根據(jù)彈性力學公式可求出層面的法向應力吃，當(yz 0 時,取 CTz-=0。按前面一樣的方法計算 F ,進行迭代計算。如果 bz =IK 、VB：ON10Nu2I式中，1尹N1(1-)2N1J(1)2頂面任一點的位移為門j0N103：0 N20Nj|U4|IV4因而，節(jié)理的相對位移匚3=許一= N 3式中*UT3 = U1v1u2v2u3v3u4v4式中A =(f)T=(K 4G)f

21、2G- 蔦Cc fa2c f .:cf (z)_l劇p閔p門 目閔p(4-44)(4-45)(4-46)I56關于勁度矩陣的推導有兩種做法：（1）設節(jié)理單元的應力和相對位移之間有如下關系：a = k込込S式中-Jo k0k是按局部坐標建立的，在整體坐標系中（圖 4-7）需要進行轉(zhuǎn)換，整體坐標系下的節(jié)理單 元勁度矩陣為TTkT其中，轉(zhuǎn)換矩陣為009T =0衛(wèi)0cossin日10=| _ _ sin日cos日k =廣NTkNdxU2_2ks102kn對稱ks02ks0kn02kn-ks0- 2ks02ks0-kn0-2kr02kn-2ks0-ks0ks02ks02kn0一kn0kn02kn一其中，ks和kn分別為節(jié)理的切向和法向勁度系數(shù)。0-N1-N200 N20N1-N20 N2