材科非線性有限元法

1、虛功方程是應(yīng)力邊界條件：n （二）s=s._ = p(3-9)TT _T _i godv 二fT pdv 亠 i f pdsvvs-(3-11)可見，除表示為物理方程的本構(gòu)關(guān)系外，其它基本方程與邊界條件均與線彈性問(wèn)題相同,且虛功方程（3-11）不涉及材料性質(zhì)，因而線彈性有限元的幾何關(guān)系和由虛功方程得到的單元和整體平衡方程完全適用于非線性彈住和彈塑性問(wèn)題，它們是f = N 3e(4-1)(4-2)BdV = Rev(4-3)寸e T T(c) B odV = Rve(4-4)將式（4-1）代入（3-2），式（3-2）代入（3-19），再將（3-19）代入（4-3），得ks（）八 Re(4-5)其

2、中單元割線勁度矩陣 k s為第四章材科非線性有限元法本章所研究的非線性有限元方法是材料非線性有限元。它的非線性是由本構(gòu)關(guān)系的非線性引起的。但它和線彈性有限元一樣，都屬于小變形問(wèn)題，因而關(guān)于形函數(shù)的選取、應(yīng)變矩 陣、應(yīng)力矩陣及勁度矩陣的形式都是相同的，不同的僅在勁度矩陣是按非線性彈性或彈塑性矩陣計(jì)算的，這是材料非線性有限元的基本內(nèi)容。4.1非線性彈性有限元法4.1.1非線性彈性有限元法的基本公式對(duì)于小變形的非線性彈性問(wèn)題，其基本方程為：平衡微分方程LT b p = 0(3-8)幾何方程= Lf(3-2)物理方程2 Ds (3-19)或d b = D T d e(3-20)式中，Ds和Dt的表達(dá)式

3、見（3-37）（ 3-40），它們是應(yīng)變強(qiáng)度 I的函數(shù)。邊界條件為(3-4)位移邊界條件：fks（）二TB DsBdV同樣可由（4-4 ）得整體平衡方程整體割線勁度矩陣Ks( S)注 Re TeKs( S 八(c ) ksC)ce(4-6)(4-7)(4-8)由于Ks與位移：有關(guān)，式（4-7）是一個(gè)非線性方程組。但實(shí)際求解時(shí)不用（4-7）式,因?yàn)榍蠼猓?-7）要用直接迭代法，這一方法不但計(jì)算量大，且常常不收斂。在求解非線性方程組時(shí)， 除直接迭代法要用割線勁度矩陣外，其它方法都要計(jì)算切線勁度矩陣。為此，須對(duì)非線性彈性有限元的切線勁度矩陣進(jìn)行研究。由式（4-4）得T T=送（c ） B cdV R

4、（4-9）e由上式及Newt on法迭代公式（2-7）可得切線勁度矩陣Kt( S)4.1.2求解的選代過(guò)程 不同的非線性方程組求解方法, 處理一般都分級(jí)按增量方法計(jì)算。T八(ce).vBeTT（ce） （ B DTBdV）ce（4-10）e迭代過(guò)程中的一些具體處理方法也不相同，但對(duì)荷載的1.荷載分級(jí)對(duì)于實(shí)際的工程問(wèn)題來(lái)說(shuō)，將荷載分級(jí)施加是相當(dāng)重要的，首先，這樣做可以模擬實(shí)際施工加載過(guò)程，進(jìn)行所謂”仿真的數(shù)值分析，這時(shí)彈塑性問(wèn)題顯得更為重要，因?yàn)椴煌募?載過(guò)程將得到不同的位移和應(yīng)力計(jì)算成果。其次，對(duì)荷載分級(jí)，容易使求解過(guò)程收斂。荷載的分級(jí)首先要考慮荷載的性質(zhì)，不同類型的荷載要根據(jù)實(shí)際施工加載情

5、況分級(jí)施加。為了更精細(xì)地模擬施工過(guò)程并使迭代過(guò)程收斂，每級(jí)荷載還可按2.2介紹的荷載系數(shù)法分成更小的增量。 圖4-1為三峽永久船閘的閘室橫剖面，為了分析閘室結(jié)構(gòu)的應(yīng)力狀態(tài)和圖4-1兩側(cè)邊坡的穩(wěn)定性，所施加的荷載按加載順 序有開挖荷載、襯砌自重、邊坡內(nèi)的滲壓和 閘室內(nèi)的水壓力。由于開挖荷載的計(jì)算需要 事先了解地應(yīng)力，在沒有地應(yīng)力量測(cè)值的情 況下，還要首先計(jì)算自重應(yīng)力。很明顯，對(duì) 于這個(gè)問(wèn)題，開挖是主要荷載，根據(jù)設(shè)計(jì)的 開挖施工程序再將它分為若干級(jí)增量，其它荷載也可根據(jù)需要再分級(jí)，目的就是為了更 精確地模擬加載過(guò)程和使求解過(guò)程收斂。2增量迭代法增量迭代法實(shí)際上就是非線性方程組求解的Euler-Ne

6、wt on法，即將荷載分成若干級(jí)增量，對(duì)每一級(jí)荷載增量，進(jìn)行迭代運(yùn)算。對(duì)第m級(jí)荷載，迭代的基本公式是血=(K T,m)( R m 一 F ；)=( KRm - 此)* i i(2-46)兇十=審+山式中Rm 一一第m級(jí)荷載增量施加后的總荷載，R m 第m級(jí)荷載增量，F(xiàn)m第m級(jí)荷載第i次迭代結(jié)束時(shí)的結(jié)點(diǎn)力，根據(jù)當(dāng)時(shí)的單元應(yīng)力二m按下式計(jì)算T TFm =無(wú)(ce) B olrndV(4-11)vm第m級(jí)荷載第i次迭代的不平衡力m =F m R m A切線勁度矩陣KT ,m根據(jù)式(4-9)由下式求出T t iKT,Z (ce) (fB Dt(寫B(tài)dV)ce(4-12)Lve如果已知第m級(jí)荷載增量時(shí)第

7、i次迭代的近似解決，相應(yīng)的應(yīng)變?yōu)?Bc e-；m(4-佝于是切線彈性矩陣 D t( 4)可以確定，由(4-12 )和(4-11 )分別計(jì)算(K tX和F m，最后利用式(2-46)求出=臨，這個(gè)迭代過(guò)程由i=0開始直至=臨足夠小，達(dá)到一定計(jì)算精度內(nèi)為零時(shí)終止。假定最后一次迭代i-l，則該級(jí)荷載增貴的最終位移為I3m二3m亠一 詁i =0(4-14)按(2-46)求出位移增量3m后，應(yīng)力增量為m 二 Dt( 4) Bce 3m(4-15)第m級(jí)荷載增量的最終應(yīng)力為Im 怖-i =0(4-16)一維問(wèn)題的迭代過(guò)程見圖4-2?？梢钥闯?，上一級(jí)荷載增量作用下，迭代收斂后仍然存在的不平衡力將轉(zhuǎn)至下一級(jí)。

8、2. 3初應(yīng)力法如果在彈性材料內(nèi)確實(shí)存在初應(yīng)力衍，則材料的應(yīng)力應(yīng)變關(guān)系為D 十衍(4-17)由上式及虛功原理可導(dǎo)出單元的結(jié)點(diǎn)力為eeT(4-18)F = k 3 + B %dV集合單元得出以下的有限元支配方程(4-19)式中，R0為由初應(yīng)力引起的等效結(jié)點(diǎn)荷載圖4-2e T TR；十八(ce)vBOodVe(4-20)對(duì)于非線性問(wèn)題，在求解非線性方程組時(shí)，可利用式（4-19）。此時(shí)，非線性的應(yīng)力應(yīng)變關(guān)系仍由式（4-17）表示，但 并非真正的初應(yīng)力，它是實(shí)際應(yīng)力b和彈性應(yīng)力 /的差，是一種虛擬的初應(yīng)力。從圖4-3可以看出，二e=D ；為彈性應(yīng)力，殆是由于應(yīng) 力應(yīng)變非線性而降低的應(yīng)力（ 衍為負(fù)值）。

9、初應(yīng)力法就是將初應(yīng) 力b看作是變化的，以此來(lái)反映應(yīng)力和應(yīng)變之間的非線性關(guān)系。 通過(guò)不斷地調(diào)整初應(yīng)力，使線彈性解逼近非線性解。前面所談的增量迭代法，在迭代的每一步都要重新形成整體 的切線勁度矩陣 K；,m并完整地求解一次線性方程組。如果將式（2-46 ）中的K T,m改用初始的切線勁度矩陣K T,m或K T,1 ,就得到常勁度的迭代公式。這相當(dāng)于求解非線性方程的 Euler修正Newt on法，（K T ）0是加載一開始時(shí)的切線勁度短陣，也 就是彈性勁度矩陣。K ；m可以認(rèn)為是第 m級(jí)荷載增量開始施加時(shí)的“彈性”勁度矩陣K:m 二 K m 八（ce）T vDBdVce（4-21）e式中的“彈性”

10、矩陣 D實(shí)際上是DT（蠱）。式（2-46）可寫成根據(jù)式（4-21）可得出下面的恒等式k m( c e)T B( /)mdv(4-23)其中(4-24)簡(jiǎn)稱彈性應(yīng)力。將式（4-22）的第二式代入第一式，并利用等式（4-23）及（4-11）可得km 陽(yáng)二 Rm （Ro）m（7=-遲 0廬 t（ om-（e）m）dv(4-24)與（4-20）相比，非線性求解用的虛擬初應(yīng)力為(4-26)（=血 一（/）m式（4-25）是直接計(jì)算總位移的迭代公式，一維的迭代過(guò)程見圖4-3。圖4-3設(shè)經(jīng)過(guò)M次迭代后在容許的精度范圍內(nèi)達(dá)到陷上器，則迭代結(jié)束，并由- 3m，按真實(shí)的本構(gòu)關(guān)系求出最后的應(yīng)力結(jié)果。在迭代過(guò)程中，o

11、m和（膚）m分別是由注按真實(shí)本構(gòu) 關(guān)系和”線彈性計(jì)算出來(lái)的應(yīng)力。（rZo）m在迭代過(guò)程中并不趨于零，而是趨于一常數(shù)矢量（圖4-3 （ b）中的M ）,它是與收斂解的真實(shí)應(yīng)力=bm和彈性應(yīng)力（be）M之差相應(yīng)的等效結(jié)點(diǎn)力。（匚0）m = ；：m是迭代結(jié)束時(shí)的虛擬初應(yīng)力。 因而，迭代過(guò)程也相當(dāng)于尋求一個(gè)合適的”初應(yīng)力場(chǎng)。0 mRm由式（4-25）可得第（i-1 ）次迭代的公式(Ro)將它與（4-25）相減得(4-27)這是求解位移增量的迭代公式。從圖4-3 （b）可見，在進(jìn)行第二次迭代（即 i=1 ）時(shí)，若按式（4-25 ）求總位移，方程 右端等于M2 A，若按式（4-27）求位移增量，方程右端為

12、M, M,二M2M2，對(duì)應(yīng)于1 IP圖4-3（ a）中的（s）m （線段 NN ）。第三次迭代（i=2 ）時(shí)，式（4-25）右端為M3 A，F(xiàn)fFFFF而式（4-27）右端為（RJ； （Rb）m =M2 M2 -M, M, = M 2M2，與（a）中的 N2N2所示的應(yīng)力對(duì)應(yīng)，經(jīng) M次迭代后，式（4-27）右端接近于零，此時(shí)，也趨于零。本級(jí) 的最后位移為Mm 二鵡=也.盡（4-28）i=0實(shí)際上，式（4-27 ）右端就相當(dāng)于失衡力。4.2彈塑性有限單元法彈塑性問(wèn)題屬于小變形，因而上節(jié)關(guān)于非線性彈性的一套基本公式及迭代運(yùn)算方法仍 然適用于彈塑性有限元，但有兩點(diǎn)必須強(qiáng)調(diào)。第一是應(yīng)力和變形呈非線性關(guān)

13、系，勁度矩陣與應(yīng)變有關(guān)， 是變化的，本構(gòu)方程必須用增量形式表示，非線性方程組的求解必須用增量方法，這與非線性彈性問(wèn)題是類似的。第二是彈塑性問(wèn)題中的材料本構(gòu)關(guān)系與應(yīng)力、變形的歷史有關(guān)，需要利用加卸載準(zhǔn)則判斷單元的材料究竟處于彈性狀態(tài)還是塑性狀態(tài)，以決定由應(yīng)變計(jì)算應(yīng)力時(shí)采用何種本構(gòu)關(guān) 系。在按照上節(jié)介紹的增量迭代法或初應(yīng)力法求出位移增量.: S后，應(yīng)變?cè)隽坑上率接?jì)算二 & = B=序二 Be乜 3（4-29）由應(yīng)變?cè)隽扣檀_定應(yīng)力增量 = C時(shí)，需要知道厶&和厶C之間的關(guān)系。彈塑性材料的本構(gòu)方 程為dDepd 廠(D - Dp)d (4-30)這是以無(wú)限小增量 de和d&的形式給出的。而在有限元數(shù)值

14、計(jì)算中，由于荷載增量是以有限大小的形式 R給出，所以，應(yīng)力增量和應(yīng)變?cè)隽慷家杂邢薮笮〉男问浇o出，設(shè)分別是厶e和也,可以利用數(shù)值積分的方法從式（4-30）得到 &和也e之間的關(guān)系,.m丄七#心 e= JDepd （4-31）m丄在荷載增量 R作用前后，所考慮單元的高斯積分點(diǎn)上的介質(zhì)究竟處于什么狀態(tài)，可以用第三章介紹的應(yīng)變空間內(nèi)的加卸載準(zhǔn)則來(lái)判斷。 0加載幾何上就是當(dāng)彈性應(yīng)力增量d膚指向屈服面外側(cè)，使f ( c + d , k) 0時(shí)為加載，材料處于塑性狀態(tài)。當(dāng)彈性應(yīng)力增量d膚指向屈服面內(nèi)側(cè)或與屈服面相切，使ef ( b d , k) 0(4-34)時(shí)為卸載和中性變載， 介質(zhì)處于彈性狀態(tài)。 值得注

15、意的是，這兒用于判斷的應(yīng)力增量 d /是彈性應(yīng)力增量，真正的應(yīng)力增量de還未求出，當(dāng)d &已知時(shí)，很容易由彈性本構(gòu)關(guān)系求出d /，因而實(shí)際上是以應(yīng)變?cè)隽?d &來(lái)判斷的。在有限元方法中，根據(jù)式(4-33)和(4-34),可由應(yīng)變?cè)隽?： &來(lái)判斷材料所處的狀態(tài)。設(shè)某個(gè)單元的某個(gè)積分點(diǎn)在荷載R m作用下的應(yīng)力為OmJ，內(nèi)變量為- mJ，施加荷載增量 Rm后，產(chǎn)生的應(yīng)變?cè)隽繛锳 ,彈性應(yīng)力增量 A = DA &。如果fm = f( m A一 e ,m)-0表示應(yīng)力emj ： L ee在屈服面內(nèi)(彈性或卸載)或在屈服面上(中性變載)，沒有新的塑性應(yīng)變產(chǎn)生，因而該點(diǎn)處于彈性狀態(tài)，可采用彈性本構(gòu)關(guān)系由&

16、求厶O-,即匚 e = D= (4-35)如果fm4 = f(怖,54)=0fm = f( em 4e m J)0說(shuō)明該點(diǎn)始終處于塑性狀態(tài)，產(chǎn)生了新的塑性應(yīng)變。當(dāng)應(yīng)變?cè)隽扣?較小時(shí)，有e = D ep (4-36)當(dāng) ： &較大時(shí)，須按(4-31)式計(jì)算應(yīng)力增量.：e。如果fm 4 - f ( em 4 , m 4 ) ： 0fm = f( Om厶 eSmJ0即從荷載增量施加前的彈性狀態(tài)進(jìn)入施加后的塑性狀態(tài)，此時(shí)，可由條件f(emr ee, mJ =0(4-37)確定應(yīng)力增量中的彈性部分與總應(yīng)力增量A c之比r （0r=N十也込山邛/ArZkt*%nT1+3q 爲(wèi)i +的 *0由 7X%：-1

17、 +rzW* * g-J 0 型云牛-確定r幾=%+2也&= OD+U幾1幾）込 j +&返何4.3特殊單元在非線性有限元中的應(yīng)用在非線性有限元分析中，為了模擬材料的某些特殊性質(zhì)或結(jié)構(gòu)上的一些特征，經(jīng)常要使用一些特殊單元。4.3.1不抗拉單元對(duì)于由不抗拉材料，如破碎巖體、土體，可采用不抗拉單元。此時(shí)要求出單元的應(yīng)力，與材料內(nèi)的初應(yīng)力迭加，再計(jì)算主應(yīng)力。并檢驗(yàn)是否出現(xiàn)拉應(yīng)力，對(duì)出現(xiàn)的主拉應(yīng)力需要進(jìn)行應(yīng)力轉(zhuǎn)移”。例如當(dāng)二1.0時(shí)，認(rèn)為二i不可能大于零，超過(guò)部分應(yīng)轉(zhuǎn)移至周圍其它單元， 并按應(yīng)力與主應(yīng)力的關(guān)系，確定實(shí)際存在的應(yīng)力再按F 二 (ce)T BT odVve求出結(jié)點(diǎn)力，進(jìn)行迭代計(jì)算。類似地，

18、可以定義低抗拉單元。4.3.2層狀單元若巖體中有一組沒有膠結(jié)的平行層面，那么這組層面就決定了巖體不能承受拉力的方 向，因而方向明確。設(shè)層面的產(chǎn)狀為傾向，傾角一：，現(xiàn)以北方向N為整體坐標(biāo)系的x正向，E為y正向，按右手系，z向下為正。并令局部坐標(biāo)系的z正向與層面的向上法線一致，x , y位于層面內(nèi)。則 z的方向余弦為l 二 cos ： sin :m = sin sin ：(4-42)n = - cos :由于這類介質(zhì)有強(qiáng)烈的各向異性，破壞的形式只可能是沿層面的剪切滑移或垂直層面的拉開。由有限元方法求出 b后，根據(jù)彈性力學(xué)公式可求出層面的法向應(yīng)力吃，當(dāng)(yz 0時(shí),取CTz - =0。按前面一樣的方

19、法計(jì)算 F ,進(jìn)行迭代計(jì)算。如果bz =IK 、Vb： o N1 0Nu2 I式中，1 尹N1(1- )2N1 J(1)2頂面任一點(diǎn)的位移為門 j 0N1 03：0 N20Nj|U4|IV4因而，節(jié)理的相對(duì)位移匚 3=許一 = N 3式中(4-44)(4-45)(4-46)*Ut3 = U1v1u2v2u3v3u4v40-N1-N200 N20N1-N20 N2關(guān)于勁度矩陣的推導(dǎo)有兩種做法：（1）設(shè)節(jié)理單元的應(yīng)力和相對(duì)位移之間有如下關(guān)系：a = k 込 S式中由虛功原理可得節(jié)理單元的勁度矩陣為。-Jo k0k =廣 n T kNdxU2_2ks102kn對(duì)稱ks02ks0kn02kn-ks0- 2ks02ks0- kn0-2kr02kn- 2ks0-ks0ks02ks02kn0一 kn0kn02kn 一其中，ks和kn分別為節(jié)理的切向和法向勁度系數(shù)。k是按局部坐標(biāo)建立的，在整體坐標(biāo)系中（圖 4-7）需要進(jìn)行轉(zhuǎn)換，整體坐標(biāo)系下的節(jié)理單 元?jiǎng)哦染仃嚍門TkT其中，轉(zhuǎn)換矩陣為0 09T =0衛(wèi)0cossin 日 10=| _ _ sin 日cos 日關(guān)于kn和ks的取值，可以通過(guò)試驗(yàn)確定它們的取值。一般來(lái)說(shuō)，當(dāng)圖4-7-/ - 0取