基于多項式型退化函數的脆性斷裂相場模型

余遠鋒, 鄭曉亞, 李鵬, 張中洲, 校金友

(1.西北工業大學 航空學院, 陜西 西安 710072; 2.西北工業大學 航天學院, 陜西 西安 710072;3.西安現代控制技術研究所, 陜西 西安 710065)

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

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

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

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

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

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

1 相場斷裂模型

1.1 斷裂能量泛函

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

圖1 二維裂紋拓撲

這里,可采用泛函[21]

(1)

表征尖銳的裂紋拓撲。裂紋面密度函數γ及其導數?γ/?d分別為

(2)

對于圖1中的相場斷裂問題,根據(1)式中引入的斷裂面函數Γ(d),可定義產生新的裂紋拓撲所需的功

(3)

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

根據線彈性理論,全局儲能泛函可表示為

(4)

儲能函數ψ描述了單位體積中物體存儲的體積能,對于由斷裂導致的能量降低,可表示為

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

(5)

式中,g(d)是能量退化函數,應滿足性質

(6)

參數k為殘余剛度系數,以保證數值的穩定性。ψ0為未損傷材料的應變能函數,對于各向同性材料,可表示為

(7)

(8)

(9)

1.2 變分相場模型

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

(10)

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

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

(11)

式中,能量驅動力f為

(12)

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

(13)

2 能量退化函數

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

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

(14)

和四次函數

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

(15)

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

(16)

這表明材料在達到極限應力時,即發生脆性斷裂時,已經發生了一定程度的損傷,且隨著多項式次數的改變,損傷程度發生了相應變化。

2.1 通用多項式型退化函數

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

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

(17)

根據退化函數性質((6)式),可得到

m=n+1

(18)

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

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

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

(19)

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

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

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

(20)

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

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

圖隨p的變化過程

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

2.2 均勻化解分析

(21)

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

忽略相場梯度Δd的影響,(21)式將變為

(22)

將(20)式的導數代入(22)式可得到均勻化解

(23)

在彈性極限處,存在相場值

de=0

(24)

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

(25)

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

(26)

(27)

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

圖5 臨界相場值dc隨著n的變化過程

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

將(27)式代入(23)式,或者對σ=g(d)E0ε關于ε求導(見附錄B),都可得到臨界應變

(28)

因此,可計算出臨界應力

(29)

為了表征彈性極限應力與臨界應力特性的關系,定義應力相對差值函數

(30)

可得到應力相對差值隨n變化的關系

(31)

如圖6所示。

圖6 應力特性隨n的變化過程

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

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

3 相場模型數值實現

3.1 單元離散

因為(11)式是關于位移場和相場的耦合方程,通常采用多場有限元方法求解。首先對位移場u及其相應的應變場ε進行離散化

(32)

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

同理,相場d及其相應的梯度場d可離散化為

(33)

相應的試探函數及其導數為

(34)

其相應的容許函數類和試探函數類如(35)式所示

?x∈Ω}

(35)

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

(36)

3.2 控制方程求解

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

對于弱耦合非線性問題,通常可以忽略非對稱項,即

(37)

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

式中

4 數值算例

由于不同次數的退化函數計算得到的彈性極限應力和臨界應力不同,而實際材料發生破壞時強度不會隨退化函數發生變化,因此,可根據(29)式計算出不同退化函數計算時的長度尺度參數l。

(40)

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

圖7 幾何模型及模擬結果

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

圖8 不同對應的載荷-位移曲線

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

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

5 結 論

本文采用一種通用的多項式型退化函數分析了AT2相場斷裂模型的力學特性,通過研究得到了以下結論:

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

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

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

附錄A

由(22)式可得到應變

(44)

對ε關于d求導,得到

(45)

對σ=g(d)E0ε關于d求導

(46)

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

(47)

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

(48)

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

(49)

整理化簡可得到一個關于d的三次多項式,但由于多項式較為復雜,此處不再展示,利用Matlab或者Mathematica求解(49)式,可得到一個關于n,p的多項式

d=d(n,p)

(50)

附錄B

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

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

(51)

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

(52)

對(52)式關于ε求導得

(53)

(54)

如(28)式所示。