廣義擴(kuò)展有限元法的研究

1廣義擴(kuò)展有限元法自出版以來，xferm法被廣泛用于制定傳統(tǒng)的有限元法框架（gfm）和xferm法的兩個新方法，自成立以來，它引起了學(xué)術(shù)界的廣泛興趣。廣義有限元法將傳統(tǒng)有限元結(jié)點(diǎn)自由度廣義化,通過增加插值函數(shù)的階次,有效提高了計(jì)算結(jié)果的精度,且在非線性分析方面較傳統(tǒng)有限元法更為可靠。廣義有限元法的概念由Babushka首次提出,隨后Babuska、Strouboulis、梁國平和欒茂田等對廣義有限元法進(jìn)行系統(tǒng)研究,取得了許多重要成果。擴(kuò)展有限元法的優(yōu)勢則在于分析不連續(xù)問題,1999年,Belytschko教授在研究常規(guī)有限元框架下富集裂尖漸進(jìn)位移場形函數(shù)的主要項(xiàng)時,首次提出了擴(kuò)展有限單元法的思想。該方法以單位分解法為基礎(chǔ),將非連續(xù)結(jié)構(gòu)直接嵌入單元內(nèi)部,通過增加描述非連續(xù)性的附加函數(shù),模擬裂紋的不連續(xù)性及裂紋尖端的奇異性,避免了傳統(tǒng)有限元法分析不連續(xù)問題時需要重新剖分網(wǎng)格的缺點(diǎn)。擴(kuò)展有限元法問世后在國際上引起了極大關(guān)注,得到了快速發(fā)展和廣泛的應(yīng)用。在裂紋擴(kuò)展分析中,一般是依據(jù)斷裂力學(xué)的相關(guān)準(zhǔn)則判斷裂紋是否穩(wěn)定,解的精度往往成為數(shù)值模擬的關(guān)鍵;與傳統(tǒng)有限元法相比,擴(kuò)展有限元法本身對解的精度并沒有提高。本文結(jié)合上述兩種方法各自的特點(diǎn),提出了一種新的數(shù)值方法———廣義擴(kuò)展有限元法(GXFEM)。文中對擴(kuò)展有限元法中的非連續(xù)附加位移函數(shù)進(jìn)行了改進(jìn),闡述了廣義擴(kuò)展有限元法的基本原理,推導(dǎo)了相關(guān)公式,編寫了廣義擴(kuò)展有限元法程序,應(yīng)用所提出的廣義擴(kuò)展有限元法對裂紋擴(kuò)展過程進(jìn)行了計(jì)算模擬,有效提高了數(shù)值解的精度。2廣義差分和擴(kuò)展差分2.1廣義有限元分析傳統(tǒng)有限元中位移插值基函數(shù)取為Φ={φ1,φ2,…,φn},單元結(jié)點(diǎn)位移向量δe=[u1,u2,…,un]T,則單元任一點(diǎn)的位移可表示為而廣義有限元將傳統(tǒng)有限元的結(jié)點(diǎn)位移廣義化,認(rèn)為各結(jié)點(diǎn)位移具有任意多個廣義位移的展開式,即將結(jié)點(diǎn)位移表示為廣義自由度di的函數(shù),即式中de=[d1,d2,…,dn]T,F是由結(jié)點(diǎn)位移插值函數(shù)組成的矩陣。將式(2)代入式(1)得到新的插值逼近空間,即式中di=[di,1,di,2,…,di,kmi]T為結(jié)點(diǎn)的廣義位移,mi為結(jié)點(diǎn)位移展開式的項(xiàng)數(shù),k為維數(shù),Ni稱為廣義形函數(shù)。其與傳統(tǒng)形函數(shù)之間的關(guān)系為式中F(i)為結(jié)點(diǎn)位移的插值函數(shù),可任意選取,一般取為坐標(biāo)的多項(xiàng)式函數(shù)。以平面四結(jié)點(diǎn)等參單元為例,若結(jié)點(diǎn)位移插值函數(shù)取為零階,有則每個結(jié)點(diǎn)具有兩個自由度,即退化為傳統(tǒng)有限元法。若結(jié)點(diǎn)位移插值函數(shù)取為一階,有則每個結(jié)點(diǎn)具有6個自由度。若結(jié)點(diǎn)位移插值函數(shù)取為二階,有則每個結(jié)點(diǎn)具有12個自由度。有了單元的位移模式,便可進(jìn)行單元分析,類似傳統(tǒng)有限元一樣,形成廣義單元勁度矩陣和廣義單元荷載列陣,再進(jìn)行整體分析,形成整體結(jié)構(gòu)的廣義支配方程:式中K為廣義整體勁度矩陣,R為廣義整體荷載列陣,D為所求各結(jié)點(diǎn)的廣義位移。因?yàn)閺V義形函數(shù)中包含傳統(tǒng)有限元中的形函數(shù),故該方法的完備性與協(xié)調(diào)性自然滿足。并且由于廣義有限元法不同結(jié)點(diǎn)位移的多項(xiàng)式之間的協(xié)調(diào)性沒有要求,因此可以按照實(shí)際問題的需要去選取多項(xiàng)式的階次,具有很強(qiáng)的靈活性和適應(yīng)性。2.2附加非連續(xù)函數(shù)式中ui為節(jié)點(diǎn)i的自由度,Ni為傳統(tǒng)有限元單元中節(jié)點(diǎn)i的形函數(shù),ai是節(jié)點(diǎn)i的附加自由度,φi(x)為附加非連續(xù)函數(shù),是考慮非連續(xù)性結(jié)構(gòu)的影響而給節(jié)點(diǎn)i的豐富函數(shù)。在處理裂紋擴(kuò)展問題時,附加非連續(xù)函數(shù)的構(gòu)造是其中的一個重要環(huán)節(jié)。裂紋在有限元網(wǎng)格中有兩種情況:一種是把單元截成兩塊,另一種則是裂紋尖端終止于單元內(nèi)部如圖1所示。對于橫跨單元的裂紋,構(gòu)造的附加非連續(xù)函數(shù)必須體現(xiàn)位移在裂紋處的不連續(xù)性。此外,單元中的點(diǎn)距離裂紋越遠(yuǎn),裂紋對該點(diǎn)的影響就越弱。Belytschko教授在提出擴(kuò)展有限元法時,采用的是單位階躍函數(shù),難以反映裂紋影響的衰減效應(yīng)。本文將指數(shù)間斷函數(shù)作為附加非連續(xù)函數(shù),先定義最短距離d(x)為式中f(x)=0為裂紋的曲線方程,然后構(gòu)造的附加非連續(xù)函數(shù)φi(x)如下。該指數(shù)間斷函數(shù)的曲線如圖2所示。對于尖端終止于單元內(nèi)部的裂紋,構(gòu)造的附加非連續(xù)函數(shù)除了要求能反映裂紋處的位移不連續(xù)性外,還必須能體現(xiàn)裂紋尖端漸進(jìn)位移場的奇異性。無論是哪一種裂紋形式(Ⅰ,Ⅱ,Ⅲ),裂紋尖端位移場均可由下面幾個基函數(shù)組成的函數(shù)形式來表達(dá),即根據(jù)上述位移模式,便可類似傳統(tǒng)有限元法,進(jìn)行單元分析和整體分析,建立擴(kuò)展有限元的支配方程,求出主變量u和附加變量a。3廣泛擴(kuò)展限制法3.1不連續(xù)變量構(gòu)造廣義擴(kuò)展有限元的位移模式為與擴(kuò)展有限元相比,此處僅對形函數(shù)作了變動,其中的Niu和Nia分別為主變量的廣義形函數(shù)及附加不連續(xù)變量的廣義形函數(shù)。式中Fu(i)和Fa(i)分別為主變量的結(jié)點(diǎn)位移插值函數(shù)和附加不連續(xù)變量的結(jié)點(diǎn)位移插值函數(shù)。由于不同結(jié)點(diǎn)位移的插值函數(shù)之間的協(xié)調(diào)性沒有強(qiáng)制要求,理論上Fu(i)和Fa(i)可以隨意選取。為了方便數(shù)值計(jì)算,兩者均取為一階的多項(xiàng)式插值函數(shù),即則位移模式可表示為與式(6)相比,只需用F1(i)φi去替換Ni,就可以得到廣義擴(kuò)展有限元法的相關(guān)公式。3.2單元控制律設(shè)計(jì)等價(jià)于平衡條件的虛功方程為式中uh∈Uh,δuh∈U0h分別是真實(shí)位移和虛位移,b為體力,t為邊界Γt上的面力。將物理方程σ=C:ε和幾何方程ε=▽su代入,并進(jìn)行區(qū)域離散化后,得引入位移模式,并將單元結(jié)點(diǎn)自由度用整體結(jié)點(diǎn)自由度表示為ue=guu,ae=gaa,其中g(shù)u和ga為選擇矩陣。根據(jù)δue和δae的任意性,得到支配方程:式中u和a分別代表主變量和附加非連續(xù)變量,式中Bu和Ba分別為主變量的單元應(yīng)變矩陣和附加非連續(xù)變量的單元應(yīng)變矩陣,f和q分別代表主變量對應(yīng)的荷載列陣和附加非連續(xù)變量對應(yīng)的荷載列陣。3.3單元域的gauss積分keuu的積分方法采用傳統(tǒng)有限元法中相同的做法,而對于keua和keaa,由于Ba的非連續(xù)性,可采用分塊積分的辦法來進(jìn)行。具體處理方法是,對于被裂紋截?cái)嗟膯卧?可分別在截?cái)嗟膬蓚€子域中進(jìn)行Gauss積分,最后把各子域上的積分值相加即得整個單元域的積分值。需要注意的是,若子域?yàn)槲暹呅螘r應(yīng)再分割成一個三角形和一個四邊形,三角形的數(shù)值積分按退化的四邊形處理;對于裂紋尖端終止的單元,則首先用過裂紋尖端的垂直裂紋方向的直線把單元分成兩類單元,然后再在各子單元上數(shù)值積分,如圖3所示。此外,對于廣義擴(kuò)展有限元,在選用一階的結(jié)點(diǎn)位移插值函數(shù)時,由于單元勁度矩陣中各元素被積函數(shù)的階次提高,積分點(diǎn)數(shù)目也要有所增加,如對于平面四結(jié)點(diǎn)單元,宜對每個單元采用9個(3×3個)Gauss積分點(diǎn)進(jìn)行數(shù)值積分。3.4裂紋演化的計(jì)算過程采用J積分法計(jì)算裂尖應(yīng)力強(qiáng)度因子。對于線彈性材料,在準(zhǔn)靜態(tài)裂紋擴(kuò)展條件下,J積分即為裂紋擴(kuò)展能量釋放力,且與積分路徑無關(guān)。式中W為應(yīng)變能密度函數(shù),Ti=σijnj為應(yīng)力矢量,Г為環(huán)繞縫尖由裂紋下表面走向上表面的圍線。彈性材料的J積分與混合模型裂紋擴(kuò)展能量釋放率及應(yīng)力強(qiáng)度因子的關(guān)系為:對平面應(yīng)變問題,α=(1-ν2)/E,對平面應(yīng)力問題,α=1/E。純Ⅰ型裂紋或純Ⅱ型裂紋可由J積分直接求得相應(yīng)的應(yīng)力強(qiáng)度因子。對混合模式裂紋,由上式只能得到K1和K2組合值,可通過引入輔助平衡狀態(tài),計(jì)算相互作用能量積分得到分離的K1和K2值。裂紋演化的計(jì)算過程如下。(1)求出裂紋尖端的應(yīng)力強(qiáng)度因子,由斷裂判據(jù)判斷裂紋是否失穩(wěn);若裂紋失穩(wěn),則裂紋按擴(kuò)展的方向伸長一個距離ΔL,擴(kuò)展方向取為裂紋尖端最大周向應(yīng)力的方向。(2)進(jìn)行應(yīng)力釋放,法向以抗拉強(qiáng)度的大小釋放應(yīng)力,切向以抗剪強(qiáng)度的折減程度釋放應(yīng)力,釋放的應(yīng)力作為裂紋表面荷載等效到擴(kuò)展單元的節(jié)點(diǎn)上,重新進(jìn)行計(jì)算。(3)根據(jù)新求出的裂尖應(yīng)力強(qiáng)度因子判斷是否擴(kuò)展。若擴(kuò)展,則增大裂紋的擴(kuò)展距離;若不擴(kuò)展,則減小裂紋的擴(kuò)展距離。(4)對于每一荷載增量步如此往復(fù)進(jìn)行,直到裂紋穩(wěn)定為止。輸出該載荷步的裂紋擴(kuò)展方向和擴(kuò)展距離,然后繼續(xù)下一增量步的計(jì)算,這樣便可完成各增量步下的裂紋擴(kuò)展方向和擴(kuò)展距離,達(dá)到裂紋追蹤的目的。4計(jì)算與分析4.1廣義擴(kuò)展有限元法如圖4所示的平面薄板,板寬B=2m,高H=3m,板左側(cè)中部有一水平裂紋,裂紋長為a,板上下兩端受均布拉力σ=1kPa的作用。板的彈性模量E=10GPa,泊松比為0.167。當(dāng)ξ=a/B≤0.6時,此問題的解析解為:式中C為尺寸相關(guān)系數(shù),由下式給出,即利用廣義擴(kuò)展有限元法對該問題進(jìn)行數(shù)值模擬,單元劃分采用9×13形式,單元總數(shù)為117個,結(jié)點(diǎn)總數(shù)為140個。分別計(jì)算裂紋長度a為0.2m,0.3m,0.4m,0.5m和0.6m時對應(yīng)的應(yīng)力強(qiáng)度因子,并與解析解及擴(kuò)展有限元法數(shù)值解進(jìn)行比較,結(jié)果見表1。從表1給出的計(jì)算結(jié)果可以看出,廣義擴(kuò)展有限元法的數(shù)值解與解析解相對誤差控制在1.2%范圍以內(nèi),并且較擴(kuò)展有限元法的數(shù)值解更為接近解析解。4.2不同計(jì)算網(wǎng)格的裂紋擴(kuò)展過程圖5為一帶水平裂紋的懸臂梁,梁長L=6m,高H=2m,右端為固定端。初始裂紋水平布置在梁的左側(cè)中部,長a=2m。梁左側(cè)端部受相向的集中力作用,P1作用于梁的頂部,方向豎直向上,大小為10kN,P2作用于梁的底部,方向豎直向下,大小為15kN,彈性模量E取為25.5GPa,泊松比ν取0.17,材料的斷裂韌度為應(yīng)用廣義擴(kuò)展有限元法對該問題進(jìn)行數(shù)值模擬,分別采用五種不同的計(jì)算網(wǎng)格模擬該混凝土懸臂梁的裂紋擴(kuò)展過程:(1)15×9形式,(2)22×9形式,(3)30×9形式,(4)30×15形式,(5)50×15形式。圖6給出了不同計(jì)算網(wǎng)格的裂紋擴(kuò)展路徑。從圖6可以看出,第四種與第五種計(jì)算網(wǎng)格對應(yīng)的裂紋擴(kuò)展過程幾乎相同,表明當(dāng)計(jì)算網(wǎng)格精細(xì)到一定程度,裂紋的擴(kuò)展路徑將不再改變,可以較好的模擬裂紋的擴(kuò)展過程。5關(guān)于廣義擴(kuò)展有限元法的研究結(jié)果(1)廣義有限元法將傳統(tǒng)有限元結(jié)點(diǎn)自由度廣義化,通過提高插值函數(shù)的階次,有效地增加計(jì)算結(jié)果的精度,且在非線性分析方面較傳統(tǒng)有限元更為可靠;而擴(kuò)展有限元法通過在裂紋面單元引入附加不連續(xù)函數(shù)來構(gòu)造位移函數(shù),克服了常規(guī)有限元法分析不連續(xù)問題,特別是斷裂問題時的諸多缺點(diǎn)。本文結(jié)合廣義有限元法和擴(kuò)展有限元法各自的優(yōu)勢,提出一種新的方法———廣義擴(kuò)展有限元法(GXFEM),并利用該方法進(jìn)行裂紋分析和數(shù)值模擬。研究結(jié)果表明:利用廣義擴(kuò)展有限元法計(jì)算裂紋擴(kuò)展問題,不需要進(jìn)行過密的網(wǎng)格劃分,且網(wǎng)格在裂紋擴(kuò)展后無需重新剖分,具有相當(dāng)高的計(jì)算精度,是進(jìn)行結(jié)構(gòu)破壞演化分析的一種有效方法。(2)需要指出的是,與擴(kuò)展有限元法相比,廣義擴(kuò)展有限元法提高了位移插值函數(shù)的階數(shù),在單元數(shù)目相同的條件下,后者具有更高的計(jì)算精度,但也同時增加了計(jì)算工作量。作者曾對本文4.1中單邊裂紋平板受拉問題,基于不同網(wǎng)格,用兩種方法進(jìn)行了計(jì)算比較。結(jié)果表明,就該例而言,在達(dá)到相同的精度下,廣義擴(kuò)展有限元法比擴(kuò)展有限元法單純增加單元數(shù)目計(jì)算效率有所提高,但并不顯著。因?yàn)樵撍憷饕鞘芾瓎栴},高階位移插值函數(shù)的效果有限;對于彎曲問題,由于高階位移插值函數(shù)能更好地描述結(jié)構(gòu)的變形特征,在相同的精度下,廣義有限元計(jì)算效率比傳統(tǒng)有限元更為有效