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

（淮安市水利勘測設計研究院有限公司，江蘇 淮安 223005）

0 引 言

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

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

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

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

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

1 相場法簡介

1.1 相場法基本理論

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

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

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

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

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

1.2 能量退化法則

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

現有的AT2 相場模型中，能量退化法則多為冪函數能量退化準則以及Borden 能量退化法則。冪函數能量退化法則表達式如式（6）所示。

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

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

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

2 耦合能量退化法則

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

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

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

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

3 數值實現

3.1 有限元離散

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

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

根據式（9）～（12）可對位移場與相場進行交錯耦合求解。

3.2 程序驗證

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

Hesch［28］、Liu［29］和Zhou［30］對單邊缺口正方形板進行了拉伸斷裂模擬，如圖3所示，模型底部固定約束，兩側水平約束，在頂部施加拉伸位移。計算模型參數：彈性模量E為210 GPa、泊松比v為0.3、臨界能量釋放率Gc為27 kJ/m2、特征長度l0為0.015 mm、k為1×10-7。單元特征長度為0.005 mm。分兩步進行加載計算，第一步加載步長為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）顯示了裂縫的擴展過程，裂縫沿著與加載方向垂直的方向進行擴展，擴展路徑與Liu 等人的模擬較為接近。圖4（c）示顯示了荷載-位移曲線，就荷載峰值而言，本文結果與Hesch 的結果較為接近，在達到荷載峰值前，本文得到的荷載-位移曲線與Hesch、Liu 和Zhou 的研究基本一致，荷載超過峰值后，本文與Hesch、Liu 和Zhou 所得到的荷載-位移曲線線型基本相同。因此，本文編制的程序對裂縫擴展的模擬以及荷載-位移的捕捉是合理且準確的，可用于分析材料的靜態(tài)斷裂問題。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

6 結 論

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

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

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

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