一種模擬準(zhǔn)脆性斷裂問題的相場法能量退化法則

張 巍，莊夢如，朱瀟瀟，沈 蔚，陳 敏

（淮安市水利勘測設(shè)計(jì)研究院有限公司，江蘇 淮安 223005）

0 引 言

工程材料的裂縫形成與擴(kuò)展是材料失效和結(jié)構(gòu)破壞較為普遍的形式［1］，尤其是對水利工程中的防滲材料，如混凝土防滲體、土質(zhì)防滲體、防滲襯砌等，一旦出現(xiàn)裂縫，可能會威脅工程結(jié)構(gòu)的安全。如美國Teton大壩［2］由于高水壓引起的裂縫擴(kuò)展，導(dǎo)致潰壩事故；澳大利亞K?lnbrein 拱壩［3］由于壩踵裂縫出現(xiàn)了嚴(yán)重的滲漏，嚴(yán)重影響了工程效益?；炷?、巖石、非飽和黏性土及水泥砂漿等水利工程中常用防滲材料的裂縫擴(kuò)展多為準(zhǔn)脆性斷裂或脆性斷裂，量化這些裂縫的斷裂過程，對預(yù)防裂縫引起的材料失效具有重要的意義。

在材料準(zhǔn)脆性斷裂模擬方面，相場法基于斷裂能量變分原理，裂縫沿著使固體系統(tǒng)總能量最小的路徑擴(kuò)展，使得其可以較好地捕捉裂縫的開裂擴(kuò)展。Miehe［4］等人給出了材料脆性斷裂的相場法熱力學(xué)框架，并成功應(yīng)用于裂紋的靜態(tài)、動(dòng)態(tài)及多場耦合問題，此后相場法模擬準(zhǔn)脆性斷裂問題得到了廣泛重視與蓬勃發(fā)展［1］。在相場法基本理論中，相場變量是區(qū)分材料的斷裂區(qū)和未損傷區(qū)域的重要指標(biāo)，為了模擬斷裂區(qū)材料的力學(xué)性能變化，需要采用能量退化法則對相場變量與材料的力學(xué)特性進(jìn)行耦合，材料的力學(xué)響應(yīng)主要由能量退化法則控制［5］。因此，能量退化法則對相場法模擬材料的斷裂行為十分重要。

Bourdin［6］提出了一個(gè)二次冪函數(shù)型式的能量退化法則，由于其函數(shù)形式簡單、無需參數(shù)、易于編程實(shí)現(xiàn)，而被廣泛使用，后發(fā)展為了一種單參數(shù)冪函數(shù)型式的能量退化法則［7］。為了進(jìn)一步提高相場法在預(yù)測與裂紋形核以及現(xiàn)有裂紋擴(kuò)展相關(guān)的臨界載荷時(shí)的精度。Sargado［8］等運(yùn)用加權(quán)思想，采用冪函數(shù)型式的能量退化法則對指數(shù)型式的能量退化法則進(jìn)行修正。Wang［9］等基于廣義線性軟化法則，提出了一種適用于黏結(jié)模型的自動(dòng)校準(zhǔn)的能量退化法則。Lu［10］等從細(xì)觀損傷角度出發(fā)，運(yùn)用加權(quán)平均思想，建立了一種指數(shù)型式的材料宏觀能量退化法則。Borden［11］提出了一種單參數(shù)多次冪函數(shù)型式的能量退化法則，其可較好地捕捉材料在斷裂破壞前的線彈性行為。吳建營［1］在研究適用于脆性斷裂和準(zhǔn)脆性破壞的相場正則化內(nèi)聚裂縫模型時(shí)，提出了基于多項(xiàng)式和冪函數(shù)型式的能量退化法則。Zivkovic［12］采用等效塑性應(yīng)變，修正了二次冪函數(shù)能量退化法則的次數(shù)，形成了考慮塑性變形的能量退化法則。

在水利工程及巖土工程領(lǐng)域，現(xiàn)多采用二次冪函數(shù)型式的能量退化法則。劉國威［13，14］等研究了巖石的雙平行翼型裂縫相交和動(dòng)態(tài)裂縫三維曲面擴(kuò)展等問題以及動(dòng)力水力壓裂作用下裂縫動(dòng)態(tài)擴(kuò)展。李鵬飛［15］等分析了包含不同巖橋傾角的預(yù)制雙裂隙巖石類材料在單軸壓縮作用下的損傷和破壞過程。易良平［16］等研究了煤砂互層中水力裂縫縱向延伸影響因素。Zhou［17］等分析了類巖石材料的壓剪耦合脆性斷裂，研究了不同地應(yīng)力狀態(tài)下巖石的水力壓裂行為［18］。Bilgen［19］等模擬了巖石材料的巴西劈拉試驗(yàn)。Zhang［20］等研究了不同試樣半徑及不同初始預(yù)制裂縫寬度對直裂縫半圓彎曲巖石試樣（NSCB）的裂縫擴(kuò)展路徑及斷裂韌度的影響。Zhang［21］等在考慮拉壓不等損傷的基礎(chǔ)上，分析了單裂隙人造巖石的初始裂縫擴(kuò)展路徑與位移-荷載曲線。Reinoso［22］等結(jié)合黏結(jié)模型，模擬了巖石的巴西劈拉試驗(yàn)及含裂縫巖石的拉伸與壓縮斷裂行為。此外，侯越［23］等采用了三次多項(xiàng)式型式的剛度法則，分析了砂漿裂紋相互作用失效的行為。Bilgen［24］等采用Borden 能量退化法則，研究了不同參數(shù)對巖石張拉斷裂行為的影響。Fei［25］等提出了拉伸-剪切雙相場模型，采用統(tǒng)一相場黏結(jié)模型［26］，研究了巖石裂縫的張拉-剪切斷裂耦合行為。Wang［27］等基于統(tǒng)一相場黏結(jié)模型框架，建立了復(fù)合型斷裂的相場模型，分析了巖石類材料的斷裂行為。

本文介紹了相場法模擬準(zhǔn)脆性斷裂問題的基本原理，詳細(xì)描述了冪函數(shù)能量退化法則與Borden 能量退化法則，提出了一種模擬準(zhǔn)脆性斷裂問題的相場法耦合能量退化法則，編制了二維相場法計(jì)算程序，研究了耦合能量退化法則對單邊缺口方形板拉伸斷裂行為的影響，進(jìn)而模擬了直裂縫半圓彎曲巖石試樣斷裂行為。

1 相場法簡介

1.1 相場法基本理論

對于準(zhǔn)脆性材料斷裂的模擬而言，相場模型可以分為AT2、AT1和內(nèi)聚裂縫模型等［1］，本文主要介紹AT2相場模型，其總勢能可以采用式（1）表示。

式中：ψε(ε，φ) 為Helmholtz 自由能，ψε(ε，φ)=g(c)ψε(ε)；ψε(ε)為彈性應(yīng)變能密度；ε為彈性應(yīng)變張量；φ為相場變量，表征了尖銳裂縫彌散為有限裂縫帶的拓?fù)湫问剑淙≈捣秶鸀椋?，1］，φ=1 表示裂縫，φ=0 表示材料完好；g(c)為能量退化法則函數(shù)，c=1-φ；Gc為臨界能量釋放率；l0裂縫尺度或者特征長度；Γ代表了非連續(xù)缺陷；Ω為區(qū)域外邊界；u為位移場；b和t分別為體力與面力，如圖1所示。

圖1 裂縫幾何正則化Fig.1 Schematic graph of fracture geometric regularization

總勢能對位移場和相場變分，同時(shí)引入自由能歷史最大狀態(tài)變量H=max[ψε(ε)]作為裂縫擴(kuò)展的不可逆條件，可得相場模型的控制方程。

式中：g＇(c)=?g(c)/?φ；σ為Cauchy 應(yīng)力張量?？刂品匠痰倪吔鐥l件可按式（3）進(jìn)行表示。

1.2 能量退化法則

在相場法基本理論中，能量退化法則耦合了相場與位移場，對材料斷裂行為的模擬至關(guān)重要。能量退化法則g(c)是隨著自變量c單調(diào)遞增的函數(shù)，且滿足以下條件［1，4］。

現(xiàn)有的AT2 相場模型中，能量退化法則多為冪函數(shù)能量退化準(zhǔn)則以及Borden 能量退化法則。冪函數(shù)能量退化法則表達(dá)式如式（6）所示。

式中：m為冪函數(shù)能量退化法則的參數(shù)，m≥0。冪函數(shù)能量退化法則函數(shù)曲線均在g(c)=c2的下側(cè)，且為下凹曲線。Borden能量退化法則可按式（7）進(jìn)行表述。

式中：s為能量退化法則的參數(shù)，0＜s≤3。Borden 能量退化法則函數(shù)曲線為介于與g(c)=3c2-2c3（不包括）與g(c)=c3之間的曲線。

根據(jù)式（6）和式（7），s=2 時(shí)的Borden 能量退化法則與m=0 時(shí)的冪函數(shù)能量退化法則相同，均為g(c)=c2。參數(shù)s=2 的函數(shù)曲線為Borden 能量退化法則函數(shù)曲線上下限之內(nèi)的一條曲線，這表明Borden 能量退化法則函數(shù)曲線包括了冪函數(shù)能量退化曲線不能涉及的一塊區(qū)域，該區(qū)域?yàn)?＜s＜2的Borden 能量退化函數(shù)曲線。同時(shí)，s=3 時(shí)的Borden 能量退化法則與m=1 時(shí)的冪函數(shù)能量退化法則相同，均為g(c)=c3。參數(shù)s=3 的曲線為Borden 能量退化法則函數(shù)曲線的下限，而m=1 的曲線為冪函數(shù)能量退化法則函數(shù)曲線上限之下的一條曲線，這表明冪函數(shù)能量退化法則函數(shù)曲線包括了Borden 能量退化法則函數(shù)曲線不能涉及的一塊區(qū)域，該區(qū)域?yàn)閙＞1 的Borden 能量退化函數(shù)曲線。雖然冪函數(shù)與Borden 能量退化法則函數(shù)覆蓋域可以相互彌補(bǔ)，但是能量退化曲線的線型并不能較好融合，這導(dǎo)致材料的損傷演化路徑較為單一，無法全面反映材料的力學(xué)性質(zhì)。

2 耦合能量退化法則

為更為全面地模擬材料的斷裂行為，通過權(quán)重參數(shù)對冪函數(shù)與Borden 能量退化法則進(jìn)行耦合，形成的耦合能量退化法則，如式（8）所示。

式中：w為權(quán)重參數(shù)，0 ≤w≤1?？梢宰C明式（8）在0＜s≤3、m≥0 的情況下，滿足能量退化法則的式（4）和式（5）的要求。當(dāng)w=0時(shí)，式（8）為冪函數(shù)能量退化法則；當(dāng)w=1時(shí)，式（8）為Borden 能量退化法則。根據(jù)式（8），結(jié)合式（4）和式（5），可知參數(shù)s和m的取值范圍并不限于0＜s≤3 和m≥0，為分析方便，本文僅研究0＜s≤3和m≥0下的耦合能量退化法則。

圖2顯示了耦合能量退化法則參數(shù)s=0.1、m=4 時(shí)，不同參數(shù)w下，耦合能量退化法則的函數(shù)曲線。由圖2可知，耦合能量退化法則函數(shù)曲線涵蓋了冪函數(shù)與Borden 能量退化法則的所有區(qū)域，且函數(shù)曲線線型更加豐富，可用于較為全面地反映不同材料的能量退化行為。

圖2 耦合能量退化法則函數(shù)曲線（s=0.1、m=4）Fig.2 coupled degradation function curves（s=0.1、m=4）

3 數(shù)值實(shí)現(xiàn)

3.1 有限元離散

位移場和相場問題的殘差和可由式（9）和式（10）表示。

式中：g(c，k)=(1-k)g(c) +k；k僅為了數(shù)值計(jì)算需要，0＜k?1；Nui和Nφi為形函數(shù)；i表示每個(gè)單元的第i個(gè)節(jié)點(diǎn)。采用Newton-Raphson迭代求解，剛度矩陣Kuu和Kφφ可按式（11）和式（12）計(jì)算。

根據(jù)式（9）～（12）可對位移場與相場進(jìn)行交錯(cuò)耦合求解。

3.2 程序驗(yàn)證

本文編制了二維相場法計(jì)算程序，為了便于與其他學(xué)者的研究進(jìn)行比較，耦合能量法則參數(shù)取w=1 和s=2 或w=0 和m=0，此時(shí)耦合能量退化法則退化為二次冪函數(shù)能量退化法則。

Hesch［28］、Liu［29］和Zhou［30］對單邊缺口正方形板進(jìn)行了拉伸斷裂模擬，如圖3所示，模型底部固定約束，兩側(cè)水平約束，在頂部施加拉伸位移。計(jì)算模型參數(shù)：彈性模量E為210 GPa、泊松比v為0.3、臨界能量釋放率Gc為27 kJ/m2、特征長度l0為0.015 mm、k為1×10-7。單元特征長度為0.005 mm。分兩步進(jìn)行加載計(jì)算，第一步加載步長為5.0×10-6mm，加載至0.005 mm，第二步加載步長為1.0×10-6mm，直至方形板完全發(fā)生斷裂。

圖3 單邊缺口方形板拉伸斷裂幾何模型與邊界條件Fig.3 Geometry and boundary conditions for tension fracture

圖4（a）與4（b）顯示了裂縫的擴(kuò)展過程，裂縫沿著與加載方向垂直的方向進(jìn)行擴(kuò)展，擴(kuò)展路徑與Liu 等人的模擬較為接近。圖4（c）示顯示了荷載-位移曲線，就荷載峰值而言，本文結(jié)果與Hesch 的結(jié)果較為接近，在達(dá)到荷載峰值前，本文得到的荷載-位移曲線與Hesch、Liu 和Zhou 的研究基本一致，荷載超過峰值后，本文與Hesch、Liu 和Zhou 所得到的荷載-位移曲線線型基本相同。因此，本文編制的程序?qū)α芽p擴(kuò)展的模擬以及荷載-位移的捕捉是合理且準(zhǔn)確的，可用于分析材料的靜態(tài)斷裂問題。

圖4 裂縫擴(kuò)展路徑及荷載-位移曲線Fig.4 Fracture propagation and load-displacement curves of tensile fracture

4 耦合能量退化法則對單邊缺口方形板拉伸斷裂行為的影響

本節(jié)以3.2 節(jié)中單邊缺口方形板拉伸斷裂為例，研究了耦合能量退化法則參數(shù)對荷載峰值和荷載-位移曲線的影響。

圖5（a）顯示了當(dāng)s=0.1與m=4時(shí)，不同w值對單邊缺口方形板拉伸斷裂行為的荷載-位移曲線的影響，其中w=0與w=1.0分別對應(yīng)了m=4 的冪函數(shù)能量退化法則及s=0.1 的Borden 能量退化法則?？梢钥闯?，隨著參數(shù)w的增大，荷載峰值及其對應(yīng)的位移值都逐漸增加，且w=1.0 時(shí)的荷載峰值是w=0 時(shí)荷載峰值的2.45倍，這表明參數(shù)w對荷載峰值影響較大。同時(shí)，在荷載峰值前，荷載-位移曲線的非線性程度隨著參數(shù)w的增加而減小，當(dāng)w=1.0 時(shí)，荷載-位移曲線表現(xiàn)為線性，這說明了隨著參數(shù)w的增加，材料的斷裂行為由準(zhǔn)脆性斷裂向脆性斷裂轉(zhuǎn)變，即參數(shù)w影響了材料的斷裂進(jìn)程。

圖5（b）顯示了當(dāng)w=0.5 與m=4 時(shí)，不同s值下的單邊缺口方形板拉伸斷裂行為的荷載-位移曲線，隨著參數(shù)s的增加，荷載峰值變小。在w=0.5與m=4情況下，位移-荷載曲線均表現(xiàn)出非線性，隨著參數(shù)s的增加，荷載峰值前后的荷載-位移曲線的非線性程度增大。這說明參數(shù)s影響了材料的斷裂進(jìn)程，可以反映材料的準(zhǔn)脆性斷裂行為。

圖5（c）顯示了當(dāng)w=0.5 與s=0.1 時(shí)，不同m值下的單邊缺口方形板拉伸斷裂行為的荷載-位移曲線。隨著參數(shù)m的增加，荷載峰值有所降低，但相對圖5（a）與圖5（b）而言，荷載峰值下降幅度較小。同時(shí)，隨著參數(shù)m的增加，峰值荷載前的荷載-位移曲線非線性程度稍有增加，但是相對于參數(shù)w和s對荷載-位移曲線非線性的影響程度而言，參數(shù)m對荷載-位移曲線非線性行為影響較小。這表明參數(shù)w和s對單邊缺口方形板拉伸斷裂行為的影響大于參數(shù)m對單邊缺口方形板拉伸斷裂行為的影響。

圖5 單邊缺口方形板拉伸斷裂的荷載-位移曲線Fig.5 Load-displacement curves of tensile fracture

綜上可知，耦合能量退化法則可以反映受拉作用下單邊缺口方形板的準(zhǔn)脆性斷裂以及脆性斷裂行為，尤其是參數(shù)w和s對單邊缺口方形板拉伸斷裂的進(jìn)程影響較大，它們對荷載峰值影響也較大。

5 直裂縫半圓彎曲巖石試樣斷裂模擬

直裂縫半圓彎曲（NSCB）試樣為國際巖石力學(xué)和巖石工程學(xué)會（ISRM）建議的測量巖石斷裂性能的方法之一［31］，已廣泛應(yīng)用于巖石、黏性土、混凝土及砂漿材料的斷裂性能的測試。本節(jié)模擬了NSCB試樣的I型斷裂試驗(yàn)，主要分析耦合能量退化法則各參數(shù)對裂縫擴(kuò)展路徑及荷載-位移曲線的影響。由于NSCB 試驗(yàn)，需要對試樣施加壓力，則應(yīng)對應(yīng)變能進(jìn)行分解，本文采用Miehe［32］等人提出的方法對對應(yīng)變能進(jìn)行分解。

直裂縫半圓彎曲試樣幾何模型（見圖6）與計(jì)算參數(shù)參考Zhou 等人［33］對巖石材料的研究，選取彈性模量E為92 GPa、泊松比v為0.21、臨界能量釋放率Gc為9.6 J/m2、材料特征長度l0為0.40 mm、計(jì)算參數(shù)k為1.0×10-9。模型單元的特征長度為0.083 mm。

圖6 直裂縫半圓彎曲試樣I型斷裂幾何模型與邊界條件Fig.6 Geometry and boundary conditions for NSCB

圖7顯示了不同耦合能量退化法則參數(shù)下的NSCB 試樣I型裂縫擴(kuò)展的荷載-位移曲線。當(dāng)s=0.1 與m=4 時(shí)，隨著參數(shù)w的增加，荷載峰值及其對應(yīng)的位移均會增大，參數(shù)w=1.0時(shí)荷載峰值是參數(shù)w=0 時(shí)荷載峰值的2.43 倍。同時(shí)，當(dāng)參數(shù)w的增量相同時(shí)，隨著參數(shù)w的增加，荷載峰值增量先減小后增加。此外，隨著參數(shù)w的增加，荷載峰值前的荷載-位移曲線由非線性特征變?yōu)榫€性，這表明隨著參數(shù)w的增加，初始裂縫的開裂擴(kuò)展會由準(zhǔn)脆性變?yōu)榇嘈裕瘩詈夏芰客嘶▌t中參數(shù)w可以反映NSCB試樣的準(zhǔn)脆性斷裂及脆性斷裂。

圖7 NSCB試樣I型斷裂的荷載-位移曲線Fig.7 Load-displacement curves of NSCB

當(dāng)w=0.5 與m=4 時(shí)，隨著參數(shù)s的增加，NSCB 試樣I 型裂縫擴(kuò)展的荷載峰值減小，參數(shù)s=0.1 時(shí)荷載峰值時(shí)s=3.0 時(shí)荷載峰值的1.38 倍。同時(shí)，當(dāng)參數(shù)s=3.0 時(shí)，荷載峰值前的荷載-位移曲線表現(xiàn)為非線性，即NSCB 試樣發(fā)生準(zhǔn)脆性斷裂。當(dāng)參數(shù)s=0.1 時(shí)，荷載峰值前的荷載-位移曲線表現(xiàn)為線性，即NSCB 試樣發(fā)生脆性斷裂。隨著參數(shù)s的增加，荷載峰值前的荷載-位移曲線由線性逐漸轉(zhuǎn)變?yōu)榉蔷€性，這表明參數(shù)s影響了NSCB試樣的I型斷裂進(jìn)程。

當(dāng)w=0.5 與s=0.1 時(shí)，隨著參數(shù)m的增加，NSCB 試樣I 型裂縫擴(kuò)展的荷載峰值與荷載-位移曲線幾乎不變，這表明參數(shù)m對NSCB試樣I型裂縫擴(kuò)展的影響較小。

綜上可知，耦合能量退化法則可以反映NSCB試樣I型裂縫擴(kuò)展的準(zhǔn)脆性行為和脆性行為，通過改變耦合能量退化法則的參數(shù)可以控制NSCB 試樣I型裂縫的擴(kuò)展進(jìn)程，尤其是參數(shù)w和參數(shù)s對荷載峰值與荷載-位移曲線影響都較大。

圖8、9分別顯示了當(dāng)s=0.1、m=4、w=0.2和s=0.1、m=4、w=0.8時(shí)NSCB 試樣I 型裂縫擴(kuò)展在荷載峰值與荷載降為0 的相場變量，其可以反映裂縫擴(kuò)展情況。在荷載峰值時(shí)，兩種參數(shù)下NSCB 試樣的相場區(qū)域較小，s=0.1、m=4、w=0.2 時(shí)NSCB 試樣的相場區(qū)域小于s=0.1、m=4、w=0.2時(shí)NSCB 試樣的相場區(qū)域，且s=0.1、m=4、w=0.2 時(shí)NSCB 試樣的相場變量最大值0.79 大于s=0.1、m=4、w=0.2 時(shí)NSCB 試樣的相場變量最大值0.59。在荷載降為0時(shí)，兩種參數(shù)下NSCB試樣的相場變量最大值達(dá)到1.0，裂縫沿著初始裂縫直線擴(kuò)展，且s=0.1、m=4、w=0.2時(shí)NSCB 試樣的相場區(qū)域?qū)挾却笥趕=0.1、m=4、w=0.8 時(shí)NSCB 試樣的相場區(qū)域?qū)挾?，而兩種參數(shù)下NSCB 試樣的相場變量最大值區(qū)域基本一致。

圖8 s=0.1、m=4、w=0.2時(shí)的NSCB試樣I型裂縫擴(kuò)展Fig.8 Fracture propagation of NSCB（s=0.1、m=4、w=0.2）

圖9 s=0.1、m=4、w=0.8時(shí)的NSCB試樣I型裂縫擴(kuò)展Fig.9 Fracture propagation of NSCB（s=0.1、m=4、w=0.8）

耦合能量退化法則各參數(shù)較難通過試驗(yàn)進(jìn)行直接測量，可根據(jù)試驗(yàn)數(shù)據(jù)，通過反分析得到。通過上述分析可知，參數(shù)w和s對NSCB試樣I型裂縫擴(kuò)展影響較大，參數(shù)反演過程中，可著重考慮這個(gè)兩個(gè)參數(shù)。同時(shí)，為對NSCB試樣I型裂縫擴(kuò)展進(jìn)行合理地模擬，精確反映荷載-位移曲線的型式，耦合能量退化法則中引入?yún)?shù)w和s的十分必要。

6 結(jié) 論

（1）提出了耦合能量退化法則，用于較為全面地反映不同材料的斷裂行為，耦合能量退化法則函數(shù)曲線覆蓋區(qū)域較廣，其豐富了能量退化法則函數(shù)曲線線型。

（2）編寫了基于能量退化法則的二維相場法計(jì)算程序，以單邊缺口正方形板拉伸斷裂為例，驗(yàn)證了程序的準(zhǔn)確性。

（3）耦合能量退化法則可以反映單邊缺口方形板拉伸的準(zhǔn)脆性斷裂和脆性斷裂行為，耦合能量退化法則參數(shù)s與m越大，拉伸荷載峰值越低；增加參數(shù)w，會提高拉伸荷載峰值。

（4）模擬了直裂縫半圓彎曲巖石試樣斷裂過程，耦合能量退化法則能夠反映NSCB 試樣I 型裂縫擴(kuò)展的準(zhǔn)脆性行為和脆性行為。不同耦合能量退化法則參數(shù)下巖石試樣的裂縫擴(kuò)展過程不同，耦合能量退化法則參數(shù)m對荷載峰值影響較小，而參數(shù)s與w對荷載峰值影響較大，荷載峰值與參數(shù)s負(fù)相關(guān)，而與參數(shù)w正相關(guān)。