基于多項(xiàng)式型退化函數(shù)的脆性斷裂相場(chǎng)模型

余遠(yuǎn)鋒, 鄭曉亞, 李鵬, 張中洲, 校金友

(1.西北工業(yè)大學(xué) 航空學(xué)院, 陜西 西安 710072; 2.西北工業(yè)大學(xué) 航天學(xué)院, 陜西 西安 710072;3.西安現(xiàn)代控制技術(shù)研究所, 陜西 西安 710065)

材料失效是一個(gè)復(fù)雜的漸進(jìn)變化過程,在載荷的作用下,材料中的應(yīng)力和應(yīng)變不斷積累,在到達(dá)材料的破壞點(diǎn)之后,損傷開始發(fā)生,材料內(nèi)部形成微觀裂紋,隨著載荷繼續(xù)增加,裂紋不斷擴(kuò)展,逐步形成宏觀裂紋,最終導(dǎo)致材料破壞[1-2]。研究材料失效過程的分析方法可大致分為連續(xù)損傷力學(xué)(CDM)方法和斷裂力學(xué)(FM)方法[3-4]。連續(xù)損傷力學(xué)主要研究材料內(nèi)部微觀結(jié)構(gòu)缺陷的萌生、匯合、擴(kuò)展所造成的材料宏觀力學(xué)性能的局部劣化,通過對(duì)材料的本構(gòu)關(guān)系衰減研究來反映損傷的影響[5-6]。根據(jù)損傷機(jī)理的不同,可以建立相應(yīng)的損傷判別準(zhǔn)則和損傷演化方程,將損傷判別準(zhǔn)則引入到材料的本構(gòu)關(guān)系中,則可對(duì)整個(gè)分析過程的損傷程度進(jìn)行描述。斷裂力學(xué)主要研究宏觀裂紋對(duì)材料強(qiáng)度的影響,該方法對(duì)裂紋擴(kuò)展進(jìn)行模擬,以描述缺陷和破壞的累積[7]。

材料失效過程涉及非線性應(yīng)力、應(yīng)變變化,裂紋起始和擴(kuò)展過程等問題,一般難以進(jìn)行解析分析。因此,在結(jié)構(gòu)和力學(xué)研究領(lǐng)域,數(shù)值仿真模擬已經(jīng)成為分析材料失效過程的重要手段。數(shù)值分析方法通?？梢苑譃殡x散方法和彌散方法。離散方法如虛擬裂紋閉合技術(shù)(VCCT)[8]、擴(kuò)展有限元法(XFEM)[9]。VCCT方法在模擬裂紋擴(kuò)展時(shí),需要預(yù)先布置裂紋,且在裂紋附近細(xì)分或不斷重新劃分網(wǎng)格。XFEM方法通過增加位移的富集項(xiàng)對(duì)間斷項(xiàng)進(jìn)行描述,但其本身的收斂性較差。并且需要預(yù)先知道材料可能的破壞模式,定義損傷后相應(yīng)的開裂方向[9]。彌散的數(shù)值方法如近場(chǎng)動(dòng)力學(xué)方法[10]、相場(chǎng)方法[11],這些方法將裂紋彌散化,描述裂紋的擴(kuò)展過程。近場(chǎng)動(dòng)力學(xué)方法將固體力學(xué)中的變形協(xié)調(diào)方程和力平衡方程改寫為積分方程,從而描述裂紋的萌生、擴(kuò)展、分叉等問題[12-13]。然而,該方法難以與材料的本構(gòu)模型直接對(duì)應(yīng)。

相場(chǎng)方法研究起源于1990年下半年,從那時(shí)起,就成為了一個(gè)重要的研究方向[14]。該方法通過引入一個(gè)場(chǎng)變量(序參數(shù))——相場(chǎng),來描述系統(tǒng)中不同物理場(chǎng)之間的轉(zhuǎn)變,該變量在0和1之間變化,分別表示材料的完好狀態(tài)和完全破壞狀態(tài),從而對(duì)裂紋的不連續(xù)性質(zhì)進(jìn)行描述。利用該方法可以模擬復(fù)雜的斷裂過程,如裂紋的萌生、擴(kuò)展、分叉和結(jié)合等問題[15]。

目前,分析材料失效的相場(chǎng)斷裂模型主要分為2類:第一類是以物理學(xué)中的Ginzburg-Landau理論為基礎(chǔ)[16-17],應(yīng)用非約束裂紋相場(chǎng)的Ginzburg-Landau型演化方程和非凸衰減函數(shù)描述材料的損傷演化過程,求解Ginzburg-Landau方程就可得到系統(tǒng)的位移場(chǎng)和損傷場(chǎng)。但該模型本質(zhì)上是完全黏性的,且主要適用于動(dòng)態(tài)斷裂[11]。第二類起源于斷裂力學(xué)的變分方法,是對(duì)經(jīng)典Griffith斷裂理論的延伸[18]。該理論由Francfort和Marigo[19]以及Bourdin等[20]提出,采用變分方法和Γ收斂準(zhǔn)則研究裂紋的演化過程,使研究者更易于理解其力學(xué)概念和推向工程應(yīng)用。該模型中用彈性勢(shì)能和裂紋面能描述系統(tǒng)的總勢(shì)能,對(duì)其總勢(shì)能求變分,得到力平衡方程和損傷演化方程,求解這2個(gè)方程,得到位移場(chǎng)和損傷場(chǎng),該理論進(jìn)一步由Miehe等[11]發(fā)展,并建立了兩類相場(chǎng)模型的統(tǒng)一。

Miehe等[21]又引入了歷史場(chǎng)概念,使裂紋相場(chǎng)演化滿足不可逆約束,為研究循環(huán)載荷下的裂紋擴(kuò)展過程奠定了重要的理論基礎(chǔ),并為相場(chǎng)理論的推廣和應(yīng)用做了重要鋪墊,使其能求解復(fù)雜的裂紋擴(kuò)展問題。在Miehe等[21]研究基礎(chǔ)上,相場(chǎng)方法被不斷改進(jìn)并被應(yīng)用到不同的研究領(lǐng)域。如脆性斷裂問題[22-23]、塑性斷裂問題[24-25]、動(dòng)態(tài)斷裂問題[26-27]、內(nèi)聚力斷裂問題[28-29]、材料疲勞問題[30-31]。同時(shí),研究材料范圍也不斷擴(kuò)大,如巖石類材料[32-33]、聚合物[34-35]、混凝土材料[36-37]、復(fù)合材料[38-39]。相比于其他的方法,如XFEM法、界面單元法,相場(chǎng)方法在研究材料失效過程中具有獨(dú)特的優(yōu)勢(shì),使其能對(duì)不連續(xù)的裂紋進(jìn)行描述,以模擬復(fù)雜形式的裂紋擴(kuò)展。

目前常用的相場(chǎng)模型主要有3種不同的選擇類型,即AT1和AT2模型,以及PFM-CZM模型。對(duì)于AT1模型和PFM-CZM模型,應(yīng)力-應(yīng)變關(guān)系存在線彈性階段,具有彈性極限應(yīng)力,即在載荷加載一段時(shí)間后,才會(huì)出現(xiàn)損傷。而AT2模型[40],采用常用的退化函數(shù)g2(d)=(1-d)2時(shí),應(yīng)力-應(yīng)變關(guān)系一直呈非線性形式,不存在彈性極限,即在載荷加載時(shí)損傷就立即發(fā)生,并且在材料到達(dá)極限應(yīng)力時(shí),損傷場(chǎng)已經(jīng)增大到d=0.25,在一定程度上影響了計(jì)算結(jié)果的準(zhǔn)確性。為了使AT2模型更好地刻畫材料的脆性斷裂過程,本文給出了一種通用的多項(xiàng)式型退化函數(shù),并分析退化函數(shù)對(duì)相場(chǎng)模型力學(xué)特性的影響。

1 相場(chǎng)斷裂模型

1.1 斷裂能量泛函

在相場(chǎng)模型中,對(duì)于連續(xù)固體中的裂紋Γ,可采用一個(gè)有限寬度l將裂紋彌散化,使離散裂紋可以用某些連續(xù)函數(shù)模擬,如圖1所示。這樣就可引入輔助變量d(x)∈[0,1]來描述這種尖銳裂紋,d=0表示桿的完好狀態(tài),d=1表示桿的完全斷裂狀態(tài)。

圖1 二維裂紋拓?fù)?/p>

這里,可采用泛函[21]

(1)

表征尖銳的裂紋拓?fù)?。裂紋面密度函數(shù)γ及其導(dǎo)數(shù)?γ/?d分別為

(2)

對(duì)于圖1中的相場(chǎng)斷裂問題,根據(jù)(1)式中引入的斷裂面函數(shù)Γ(d),可定義產(chǎn)生新的裂紋拓?fù)渌璧墓?/p>

(3)

式中,Gc是Griffith型臨界能量釋放率。

根據(jù)線彈性理論,全局儲(chǔ)能泛函可表示為

(4)

儲(chǔ)能函數(shù)ψ描述了單位體積中物體存儲(chǔ)的體積能,對(duì)于由斷裂導(dǎo)致的能量降低,可表示為

ψ(ε,d)=[g(d)+k]ψ0(ε)

(5)

式中,g(d)是能量退化函數(shù),應(yīng)滿足性質(zhì)

(6)

參數(shù)k為殘余剛度系數(shù),以保證數(shù)值的穩(wěn)定性。ψ0為未損傷材料的應(yīng)變能函數(shù),對(duì)于各向同性材料,可表示為

(7)

(8)

(9)

1.2 變分相場(chǎng)模型

根據(jù)裂紋耗散泛函(3)式和儲(chǔ)能泛函方程(4),引入總的能量泛函

(10)

式中:b是區(qū)域B中的體積力;t是在?B上的表面外力。

對(duì)總能量泛函(10)式求變分,可得到以下強(qiáng)形式的相場(chǎng)模型控制方程

(11)

式中,能量驅(qū)動(dòng)力f為

(12)

為了保證損傷演化的不可逆性,引入歷史參考能量[21]

(13)

2 能量退化函數(shù)

為了模擬能量退化,采用常用的多項(xiàng)式函數(shù)，如前文的g2(d),以及三次函數(shù)[41]

g3(d)=3(1-d)2-2(1-d)3

(14)

和四次函數(shù)

g4(d)=4(1-d)3-3(1-d)4

(15)

采用g3(d)和g4(d)退化函數(shù)時(shí),都可使AT2相場(chǎng)模型得到彈性極限,且其臨界相場(chǎng)值分別為

(16)

這表明材料在達(dá)到極限應(yīng)力時(shí),即發(fā)生脆性斷裂時(shí),已經(jīng)發(fā)生了一定程度的損傷,且隨著多項(xiàng)式次數(shù)的改變,損傷程度發(fā)生了相應(yīng)變化。

2.1 通用多項(xiàng)式型退化函數(shù)

為了探究隨著函數(shù)次數(shù)的改變,退化函數(shù)對(duì)脆性斷裂特性的影響,采用一種通式來表達(dá)退化函數(shù)

gm(d)=m(1-d)n-n(1-d)m,m>n≥2

(17)

根據(jù)退化函數(shù)性質(zhì)((6)式),可得到

m=n+1

(18)

對(duì)于不同的m,gm(d)的變化趨勢(shì)如圖2所示。

圖2 gm(d)隨m變化形式

從圖2可以看出,隨著m的增大,退化函數(shù)開始時(shí)平緩變化,隨后下降趨勢(shì)逐漸增大,在到達(dá)某個(gè)拐點(diǎn)后,下降趨勢(shì)開始變得緩慢,表明在開始階段材料性能下降緩慢,然后下降趨勢(shì)比較劇烈,并在某個(gè)拐點(diǎn)后,材料性能下降開始減緩。同時(shí)George等[42]也表示在材料開始發(fā)生損傷時(shí),由于材料內(nèi)部微裂紋的相互作用和抑制,在加載開始時(shí),材料損傷趨勢(shì)較慢,而在一定階段后,微裂紋間的相互阻止效應(yīng)減弱,損傷開始加快。因此,采用這種多項(xiàng)式型退化函數(shù),在一定程度上與George等人的分析較為吻合。

(19)

這將使裂紋在開始階段無法成核,因此,采用一種修正模型的退化函數(shù)

gm(d)=p[(1-d)m-(1-d)n]+

m(1-d)n-n(1-d)m,m>n≥2,p≥0

(20)

當(dāng)p=2,m=3,n=2時(shí),g3(d)恢復(fù)為經(jīng)典的退化函數(shù)g2(d)=(1-d)2。為了使相場(chǎng)模型在開始階段存在裂紋驅(qū)動(dòng)力,裂紋能夠成核,需要將p設(shè)置為一個(gè)非零正值,以保證在d=0時(shí)存在驅(qū)動(dòng)力。

圖3 gm(d)隨p的變化過程

圖隨p的變化過程

從圖3～4可以看出,當(dāng)p較大時(shí),如p=0.1,退化函數(shù)(17)式和(20)式的計(jì)算結(jié)果存在一定差別,而當(dāng)0

2.2 均勻化解分析

(21)

式中,應(yīng)力為σ=g(d)E0ε。

忽略相場(chǎng)梯度Δd的影響,(21)式將變?yōu)?/p>

(22)

將(20)式的導(dǎo)數(shù)代入(22)式可得到均勻化解

(23)

在彈性極限處,存在相場(chǎng)值

de=0

(24)

將(24)式代入(23)式,可得到彈性應(yīng)變和應(yīng)力極限

(25)

將(23)式代入σ=g(d)E0ε,并對(duì)d求導(dǎo)

(26)

(27)

從(27)式可以看出,在多項(xiàng)式型退化函數(shù)中,在n≥2情況下,臨界相場(chǎng)值恒滿足dc>0,這意味著材料在彈性階段后,還會(huì)存在一段非線性變化過程,因此,在到達(dá)脆性破壞極限時(shí),材料已經(jīng)發(fā)生了一定程度的損傷。dc隨n的變化形式如圖5所示。

圖5 臨界相場(chǎng)值dc隨著n的變化過程

從圖5可以看出,隨著n值的增大dc先增大后減小,在n=3時(shí),dc取得最大值,且當(dāng)n=2,4時(shí),兩者的臨界相場(chǎng)值dc相等。這說明對(duì)于采用gm(d)的相場(chǎng)模型,隨著n的增大,材料發(fā)生脆性斷裂時(shí)的損傷程度先增大然后逐漸減小,刻畫脆性斷裂的效果在增強(qiáng)。

將(27)式代入(23)式,或者對(duì)σ=g(d)E0ε關(guān)于ε求導(dǎo)(見附錄B),都可得到臨界應(yīng)變

(28)

因此,可計(jì)算出臨界應(yīng)力

(29)

為了表征彈性極限應(yīng)力與臨界應(yīng)力特性的關(guān)系,定義應(yīng)力相對(duì)差值函數(shù)

(30)

可得到應(yīng)力相對(duì)差值隨n變化的關(guān)系

(31)

如圖6所示。

圖6 應(yīng)力特性隨n的變化過程

從圖6可以看出,隨著n值的增大,彈性極限應(yīng)力與臨界應(yīng)力值都在減小,但彈性極限應(yīng)力減小的趨勢(shì)小于與臨界應(yīng)力的減小趨勢(shì),并且兩者之間的相對(duì)差值在不斷增大,說明隨著n值的增大,相場(chǎng)模型的線彈性變化階段在減弱,非線性變化趨勢(shì)在逐漸增強(qiáng)。

通過以上分析可知,采用多項(xiàng)式型退化函數(shù)時(shí),在材料承受載荷后,應(yīng)力-應(yīng)變首先呈線彈性變化,然后是伴隨著微小損傷的非線性變化,在達(dá)到材料破壞極限時(shí),已經(jīng)積累了一定程度的損傷。但相比較g2(d)函數(shù),這種多項(xiàng)式型函數(shù)能更好地描述材料脆性斷裂過程。

3 相場(chǎng)模型數(shù)值實(shí)現(xiàn)

3.1 單元離散

因?yàn)?11)式是關(guān)于位移場(chǎng)和相場(chǎng)的耦合方程,通常采用多場(chǎng)有限元方法求解。首先對(duì)位移場(chǎng)u及其相應(yīng)的應(yīng)變場(chǎng)ε進(jìn)行離散化

(32)

式中,插值矩陣N=[N1,…,Ni,…,Nn],幾何矩陣B=[B1,…,Bi,…,Bn]。

同理,相場(chǎng)d及其相應(yīng)的梯度場(chǎng)d可離散化為

(33)

相應(yīng)的試探函數(shù)及其導(dǎo)數(shù)為

(34)

其相應(yīng)的容許函數(shù)類和試探函數(shù)類如(35)式所示

?x∈Ω}

(35)

然后,可以得到以下弱形式的控制方程

(36)

3.2 控制方程求解

在求解(36)式時(shí),常用的方法有:整體牛頓法、分步算法、BFGS法。而由于所得到的相場(chǎng)模型控制方程的非線性,采用整體牛頓法計(jì)算時(shí)不容易收斂。為了改善計(jì)算的收斂性,可采用分步算法求解以上的相場(chǎng)斷裂模型控制方程,但這種方法求解效率較慢,計(jì)算復(fù)雜模型時(shí)非常耗時(shí)。由于BFGS法無需每次迭代都更新剛度矩陣,且計(jì)算過程較為簡(jiǎn)單,計(jì)算效率較高,被Wu等[43]應(yīng)用于控制方程的求解,在保證收斂性的同時(shí),提高了計(jì)算效率。

對(duì)于弱耦合非線性問題,通?？梢院雎苑菍?duì)稱項(xiàng),即

(37)

求解方程(37)可得到位移場(chǎng)和相場(chǎng)。

式中

4 數(shù)值算例

由于不同次數(shù)的退化函數(shù)計(jì)算得到的彈性極限應(yīng)力和臨界應(yīng)力不同,而實(shí)際材料發(fā)生破壞時(shí)強(qiáng)度不會(huì)隨退化函數(shù)發(fā)生變化,因此,可根據(jù)(29)式計(jì)算出不同退化函數(shù)計(jì)算時(shí)的長(zhǎng)度尺度參數(shù)l。

(40)

本節(jié)算例中模型尺寸為1 mm×1 mm,在模型的中間區(qū)域0.5 mm的位置有一條初始裂紋Γ,初始長(zhǎng)度a0=0.5 mm,下邊界固定,上邊界施加拉伸載荷u=0.01 mm,模型的邊界條件及載荷形式如圖7所示,其材料參數(shù)為E=210 kN/mm2,泊松比ν=0.3,能量釋放率Gc=2.7×10-3kN/mm。

圖7 幾何模型及模擬結(jié)果

計(jì)算中模型最小單元尺度選擇為h=0.5l,假定n=2,即m=3時(shí),ln=2=0.01,由于不同次數(shù)n時(shí)計(jì)算得到的臨界應(yīng)力相同,因此根據(jù)(40)式可等效計(jì)算出n=3,4,5時(shí)的尺度參數(shù)l,ln=3=0.005 3,ln=4=0.003 3,ln=5=0.002 3。

圖8 不同對(duì)應(yīng)的載荷-位移曲線

從圖8可以看出,相比于二次退化函數(shù),采用的多項(xiàng)式型退化函數(shù)使材料在斷裂發(fā)生前更好地保持線彈性變化。在初始加載時(shí),由于材料內(nèi)部微裂紋的相互抑制作用,材料損傷增加地非常緩慢,模型保持了線彈性變化,很好地描述材料的脆性斷裂性能,如g3(d),g4(d)函數(shù)得到的載荷-位移曲線。

同時(shí),從圖8可以發(fā)現(xiàn),隨著多項(xiàng)式函數(shù)次數(shù)的增大,由于其函數(shù)下降趨勢(shì)逐漸增加,材料內(nèi)部微裂紋間的相互抑制效應(yīng)減弱,使材料損傷加快,引起材料的斷裂。所以雖然在計(jì)算中假定不同退化函數(shù)最終破壞的載荷相同,即得到位移-載荷曲線應(yīng)該相同,但隨著函數(shù)次數(shù)的增加,其下降趨勢(shì)增加,即材料損傷趨勢(shì)增大,會(huì)誘發(fā)材料產(chǎn)生一定程度的微觀損傷,引起最終的斷裂破壞,如g5(d),g6(d)函數(shù)得到的載荷-位移曲線。通過該算例也說明了描述材料性能下降的退化函數(shù),其變化趨勢(shì)在一定程度上也影響了材料損傷和斷裂變化過程。

5 結(jié) 論

本文采用一種通用的多項(xiàng)式型退化函數(shù)分析了AT2相場(chǎng)斷裂模型的力學(xué)特性,通過研究得到了以下結(jié)論:

1) 通過一維簡(jiǎn)單模型分析得到了一般形式的彈性應(yīng)變和應(yīng)力、臨界相場(chǎng)值等公式,發(fā)現(xiàn)隨著函數(shù)次數(shù)的增大,臨界相場(chǎng)值先增大,然后減小,表明當(dāng)函數(shù)次數(shù)增大時(shí),理論上材料發(fā)生斷裂時(shí)的損傷程度更小。

2) 相比于常用的二次退化函數(shù),多項(xiàng)式退化函數(shù)使相場(chǎng)模型具有彈性應(yīng)力極限,特別是當(dāng)退化函數(shù)的次數(shù)較小時(shí),能很好地描述脆性斷裂行為,使材料在斷裂前保持線彈性變化,直至發(fā)生脆性斷裂。

3) 退化函數(shù)的變化趨勢(shì)會(huì)影響材料發(fā)生斷裂的時(shí)間。在函數(shù)次數(shù)較小時(shí),如m=3,4時(shí),模型斷裂前很好地保持線彈性變化。隨著函數(shù)次數(shù)的增大,如m=5,6時(shí),其函數(shù)下降趨勢(shì)逐漸增加,導(dǎo)致材料內(nèi)部微裂紋間的相互抑制效應(yīng)減弱,誘發(fā)材料更早地發(fā)生損傷,引起最終的斷裂。但這種影響相比于二次退化函數(shù),還是比較微弱。

附錄A

由(22)式可得到應(yīng)變

(44)

對(duì)ε關(guān)于d求導(dǎo),得到

(45)

對(duì)σ=g(d)E0ε關(guān)于d求導(dǎo)

(46)

將(45)式和(44)式代入(46)式

(47)

-2d[g′(d)]2+g(d)[-g′(d)+dg″(d)]=0

(48)

將(41)～(43)式代入(48)式,得到

(49)

整理化簡(jiǎn)可得到一個(gè)關(guān)于d的三次多項(xiàng)式,但由于多項(xiàng)式較為復(fù)雜,此處不再展示,利用Matlab或者M(jìn)athematica求解(49)式,可得到一個(gè)關(guān)于n,p的多項(xiàng)式

d=d(n,p)

(50)

附錄B

將退化函數(shù)帶入應(yīng)力表達(dá)式σ=g(d)E0ε中,得到

σ=[m(1-d)n-n(1-d)m]E0ε

(51)

將(23)式代入(51)式,可得到

(52)

對(duì)(52)式關(guān)于ε求導(dǎo)得

(53)

(54)

如(28)式所示。