基于相場法的巖石三點(diǎn)彎曲細(xì)觀破壞數(shù)值模擬

李明耀, 李紹金, 左建平*, 王智敏, 彭磊, 薛喜仁

(1.中國礦業(yè)大學(xué)(北京)煤炭資源與安全開采國家重點(diǎn)實(shí)驗(yàn)室， 北京 100083; 2.中國礦業(yè)大學(xué)(北京)力學(xué)與建筑工程學(xué)院， 北京 100083)

巖石材料在復(fù)雜環(huán)境下的變形、破壞及失穩(wěn)等力學(xué)行為對巖體工程結(jié)構(gòu)的強(qiáng)度和安全至關(guān)重要。巖石內(nèi)部裂紋的萌生、擴(kuò)展和貫通等力學(xué)行為是其中的關(guān)鍵力學(xué)問題，是預(yù)測巖體工程的強(qiáng)度和穩(wěn)定性的關(guān)鍵因素之一，也是一直以來中外學(xué)者研究的熱點(diǎn)和難點(diǎn)。目前，采用數(shù)值模擬方法是研究巖石內(nèi)部裂紋擴(kuò)展行為的重要手段。

目前，針對裂紋擴(kuò)展問題發(fā)展提出了不同的數(shù)值模擬方法，如有限元法(finite element method)[1]、擴(kuò)展有限元法(extended finite element method)[2]、離散元法(discrete element method)[3]、位移不連續(xù)法(displacement discontinuity method)[4]、相場法(phase field method)[5]和彌散裂紋法(smeared crack method)[6]等。這些數(shù)值方法主要可分為離散方法和連續(xù)方法兩大類。離散方法將裂紋視為真實(shí)的或等效于物理情況的強(qiáng)不連續(xù)性裂紋。如余華中等[7]基于顆粒離散元顆粒流程序研究節(jié)理的宏細(xì)觀剪切力學(xué)行為以及不同法向應(yīng)力作用下微裂紋的發(fā)育及演化規(guī)律；Alneasan等[8]基于DDM模擬和分析混合加載模式下斷裂參數(shù)對巖石裂紋擴(kuò)展的影響；楊善統(tǒng)等[9]利用隨機(jī)分布Voronoi圖生成不同裂縫產(chǎn)狀巖石模型，采用剛性彈簧方法模擬了非貫通裂隙巖石細(xì)觀裂紋起裂和擴(kuò)展過程。離散方法能夠較為真實(shí)地反映裂紋形態(tài)，但需要先進(jìn)的計(jì)算能力才能處理較為復(fù)雜的裂紋擴(kuò)展問題，特別是對于流體流動(dòng)具有強(qiáng)非線性耦合問題的水力壓裂等[10-11]。連續(xù)方法則將裂紋視為連續(xù)分布形式的連續(xù)體，如Li等[12]基于FEM傳統(tǒng)重新網(wǎng)格劃分方法模擬和預(yù)測巖石在混合模式斷裂試驗(yàn)下裂紋擴(kuò)展的規(guī)律。XFEM模擬裂紋可以避免重畫網(wǎng)格從而得到了廣泛的應(yīng)用，如Yan等[13]基于擴(kuò)展有限元法研究了含預(yù)制裂紋的巖體中裂紋萌生、擴(kuò)展和貫通的規(guī)律；牛駿等[14]將擴(kuò)展有限元與嵌入式離散裂紋模型耦合的方法模擬頁巖水力壓裂，得到了頁巖水力裂縫內(nèi)部開度的動(dòng)態(tài)變化以及流體壓力對裂紋的影響。連續(xù)方法在連續(xù)體公式中規(guī)范了裂紋的不連續(xù)性，以避免捕獲與裂紋萌生、轉(zhuǎn)向和分叉相關(guān)的路徑演變的困難，通常以直接的方式應(yīng)用于計(jì)算。數(shù)值方法處理上的優(yōu)勢使得連續(xù)方法近年來得到了廣泛的應(yīng)用，但連續(xù)方法的缺點(diǎn)在于不可避免地難以獲得裂縫寬度和表面等信息。

近年來，相場法因其在模擬裂紋擴(kuò)展方面的優(yōu)勢被越來越多的學(xué)者關(guān)注。相場法從最初Francfort等[15]提出的斷裂變分理論以及Bourdin等[16-17]提出裂紋正則化理論，到Miehe等[18-19]提出了應(yīng)變譜分解方法和交替迭代法后，被廣泛應(yīng)用于模擬裂紋擴(kuò)展問題。相場法主要通過一系列的微分方程，引入序參量，利用彌散的方式描述裂隙邊界，避免了顯式追蹤裂紋面，用序參量的自動(dòng)演化獲取裂紋的路徑及位置，在模擬裂紋起裂、擴(kuò)展、分叉和多裂紋相互貫通等現(xiàn)象具有一定優(yōu)勢，為研究復(fù)雜裂隙擴(kuò)展問題提供了一種重要手段。

Molnar等[20]提出基于擴(kuò)散裂縫的速率無關(guān)變分原理的相場法應(yīng)用于研究具有復(fù)雜初始細(xì)觀結(jié)構(gòu)的巖石材料內(nèi)部的裂紋擴(kuò)展問題。在此基礎(chǔ)上，現(xiàn)通過模擬含偏置缺口玄武巖的三點(diǎn)彎曲實(shí)驗(yàn)，研究不同偏置距離缺口的玄武巖內(nèi)部裂紋擴(kuò)展過程及其影響因素，通過與室內(nèi)試驗(yàn)的宏觀和細(xì)觀結(jié)果對比分析表明，相場法能很好地模擬和預(yù)測具有復(fù)雜初始細(xì)觀結(jié)構(gòu)的巖石材料的內(nèi)部的裂紋擴(kuò)展問題，為研究巖石內(nèi)部裂紋擴(kuò)展問題提供一種重要手段。

1 相場法的基本理論及數(shù)值實(shí)現(xiàn)

相場法是在經(jīng)典Griffith斷裂力學(xué)理論的基礎(chǔ)上，利用最小能量變分原理來克服傳統(tǒng)斷裂力學(xué)理論框架需要復(fù)雜的裂紋擴(kuò)展判據(jù)的不足，是對傳統(tǒng)斷裂力學(xué)理論的繼承和發(fā)展。相場法理論認(rèn)為一個(gè)彈性體V，其總勢能Π由彈性應(yīng)變能和裂紋表面能兩部分組成：

(1)

式(1)中：V?Rn,其中n為維度，n=1，2，3；ψe為彈性應(yīng)變能密度；Γ為裂紋表面；ε(u)為應(yīng)變張量；Gc為材料的臨界能量釋放率。變分原理認(rèn)為，裂紋的起裂、擴(kuò)展和分叉等行為都應(yīng)使總勢能取最小值且遵循不可逆條件。

為了描述裂紋的分布情況，引入一個(gè)相場變量d(x)∈[0,1]，表示材料的損傷值，當(dāng)d=0時(shí)，表示材料完好，當(dāng)d=1時(shí)，表示材料完全破壞，即斷裂。裂紋表面函數(shù)Γl(d)可以表示為

(2)

式(2)中：γ(d,?d)為裂紋表面密度函數(shù)，形式為

(3)

式(3)中：lc為特征長度參數(shù)，控制裂紋的寬度。根據(jù)Miehe等[18-19]的研究成果，將式 (2) 改寫成歐拉變分形式

(4)

式(4)中：W={d|d(x)=1,?x∈Γ}。

彈性應(yīng)變能密度ψe可以表示為

ψe(ε,d)=g(d)ψ0(ε)

(5)

式(5)中：ψ0(ε)為彈性應(yīng)變能密度；g(d)為材料剛度削減函數(shù)，g(d)=(1-d)2+κ，其中參數(shù)κ為了保證當(dāng)材料局部發(fā)生完全破壞(d=1)時(shí)數(shù)值方法的穩(wěn)定性。對于線彈性完好材料，彈性應(yīng)變能密度可以表示為

(6)

式(6)中：C0為材料的彈性張量；應(yīng)變張量ε=[(?u)T+?u]/2。由于損傷的產(chǎn)生，彈性應(yīng)變能會(huì)隨著損傷的演化而衰減，通過對應(yīng)變能密度求解應(yīng)變張量的一階導(dǎo)數(shù)可得應(yīng)力張量為

σ=g(d)σ0=[(1-d)2+κ]σ0

=[(1-d)2+κ]C0:ε

(7)

彈性張量也可以表示為

C=g(d)C0

(8)

由此可見，由損傷變量直接影響著材料的彈性張量和應(yīng)力。

由于在有限元框架下由相場法直接求解出節(jié)點(diǎn)位移和相場具有強(qiáng)非線性，求解需要較多的迭代步，計(jì)算效率低。為了提高求解效率，可以對位移場和裂縫場解耦求解。研究表明，這種求解方法的計(jì)算效率可以大幅提升，同時(shí)可靠性高。因此，本文研究中采用Molnar等[20]中的解耦數(shù)值方法，在ABAQUS有限元軟件中通過編寫用戶子程序UMAT和UEL來實(shí)現(xiàn)基于相場法的裂紋擴(kuò)展數(shù)值模擬。

2 數(shù)值驗(yàn)證

為了驗(yàn)證本文數(shù)值方法的準(zhǔn)確性，首先模擬對比了兩個(gè)經(jīng)典數(shù)值實(shí)驗(yàn)：單邊缺口拉伸試驗(yàn)和對稱三點(diǎn)彎試驗(yàn)。通過對比裂紋擴(kuò)展過程和加載位移曲線，驗(yàn)證了本文采用的數(shù)值方法能夠準(zhǔn)確、合理地模擬裂紋擴(kuò)展過程。同時(shí)，分析了網(wǎng)格數(shù)量和長度參數(shù)(lc)對實(shí)驗(yàn)結(jié)果的影響，為后續(xù)研究奠定了基礎(chǔ)。

2.1 單邊缺口拉伸試驗(yàn)

單邊缺口拉伸試驗(yàn)是經(jīng)典的數(shù)值實(shí)驗(yàn)，該模型的幾何尺寸和邊界條件如圖1(a)所示，尺寸為1 mm×1 mm，其左邊中點(diǎn)預(yù)置水平裂紋Γ0,裂紋長度為0.5 mm。模型底部固定，在頂部施加豎直向上的位移，左右兩側(cè)自由。材料參數(shù)為彈性模量E=210 kN/mm2，泊松比ν=0.3，長度參數(shù)lc=0.007 5 mm，臨界能量釋放率gc=0.002 7 kN/mm。前500步位移均勻增量Δu=1×10-4mm，隨后位移均勻增量Δu=1×10-5mm完成整個(gè)模擬過程，總位移u=0.007 mm。數(shù)值分析的網(wǎng)格如圖1(b)所示，采用四節(jié)點(diǎn)單元，數(shù)量約8 000個(gè)，在裂紋可能擴(kuò)展的區(qū)域加密網(wǎng)格，加密區(qū)網(wǎng)格尺寸h≈0.003 mm。

圖1 單邊缺口拉伸試驗(yàn)?zāi)Ｐ虵ig.1 Single edge notched test model

圖2 裂紋路徑Fig.2 Crack pattern

圖3 力與位移曲線對比Fig.3 Comparison of force-displacement curve with the reference results

圖2和圖3分別展現(xiàn)本文所模擬的裂紋路徑和力與位移曲線與Molnar等[20]的對比結(jié)果，可以看出，兩者的結(jié)果幾乎一致，驗(yàn)證了本文采用的數(shù)值方法的準(zhǔn)確性。

2.2 對稱三點(diǎn)彎試驗(yàn)

對稱三點(diǎn)彎經(jīng)典試驗(yàn)的幾何尺寸及邊界條件如圖4(a)所示。彈性模量E=20.8 kN/mm2，泊松比是ν=0.3，長度參數(shù)lc=0.06 mm，臨界能量釋放率gc=5×10-3kN/mm。采用控制位移方式進(jìn)行加載，前360步位移均勻增量是Δu=1×10-4mm，隨后位移均勻增量Δu=1×10-5mm完成整個(gè)模擬過程，總位移0.06 mm。同樣在裂紋可能擴(kuò)展的區(qū)域進(jìn)行加密網(wǎng)格，如圖4(b)所示，網(wǎng)格采用四節(jié)點(diǎn)單元，數(shù)量約15 000個(gè)，加密區(qū)網(wǎng)格尺寸h≈0.006 mm。通過對比發(fā)現(xiàn)，本文的數(shù)值模擬與Molnar等[20]結(jié)果吻合度很高，幾乎一致(圖5和圖6)，進(jìn)一步驗(yàn)證了本文采用的數(shù)值方法的準(zhǔn)確性。

圖4 對稱三點(diǎn)彎試驗(yàn)?zāi)Ｐ虵ig.4 Symmetric three point bending test model

圖5 裂紋路徑Fig.5 Crack pattern

圖6 力與位移曲線對比Fig.6 Comparison of force-displacement curve with the reference results

2.3 相場法相關(guān)參數(shù)的影響分析

2.3.1 網(wǎng)格影響分析

在驗(yàn)證了數(shù)值方法準(zhǔn)確性的基礎(chǔ)上，需要進(jìn)一步明確不同的網(wǎng)格數(shù)對數(shù)值模擬結(jié)果的影響。以上述對稱三點(diǎn)彎試驗(yàn)為分析對象，在其他參數(shù)不變的情況下，分別采用3種不同的網(wǎng)格劃分方式(網(wǎng)格數(shù)分別為6 914、19 358、194 280)，最終模擬的位移載荷曲線如圖7所示。結(jié)果表明，采用不同網(wǎng)格時(shí)的荷載位移曲線基本一致，由此可見，本文采用的相場法對網(wǎng)格數(shù)量不敏感。

圖7 不同網(wǎng)格的力與位移曲線Fig.7 The force-displacement curves of different mesh

2.3.2 長度參數(shù)(lc)影響分析

相場法通過引入正則化參數(shù)(即長度參數(shù)lc)來彌散裂紋寬度，其選擇對數(shù)值模擬結(jié)果有一定的影響。以對稱三點(diǎn)彎實(shí)驗(yàn)為例，中心加密區(qū)網(wǎng)格最小尺寸為0.015 mm，分別采用長度參數(shù)lc為0.015、0.03、0.06 mm進(jìn)行數(shù)值模擬。圖8所示為不同長度參數(shù)下的裂紋形態(tài)。通過對比發(fā)現(xiàn)，3個(gè)不同參數(shù)下，擴(kuò)展路徑幾乎一致，但lc越小，所模擬得到的裂紋寬度越接近真實(shí)裂紋。參照Miehe等[19]已驗(yàn)證結(jié)論，模型參數(shù)lc與加密區(qū)網(wǎng)格尺寸h的關(guān)系可近似表示為lc≥2h。如圖9所示為對應(yīng)的力與位移曲線，可以看出隨著長度參數(shù)lc的減小，最大荷載與失效位移均在增大。因此，在實(shí)際中選取lc時(shí)，除了滿足lc的定義lc→0以外，還應(yīng)當(dāng)通過試驗(yàn)和最小網(wǎng)格尺寸等方式來選取合理的lc，從而保證精度較高、誤差較小的裂紋數(shù)值模擬結(jié)果。

圖8 不同長度參數(shù)的裂紋形態(tài)Fig.8 Crack morphology with different length parameters

圖9 不同長度參數(shù)的力與位移曲線Fig.9 The force-displacement curves of different length parameter

3 含偏置缺口玄武巖三點(diǎn)彎曲的數(shù)值模擬

3.1 試樣描述

含不同偏置距離缺口的玄武巖三點(diǎn)彎室內(nèi)試驗(yàn)[21]的巖樣取自北京市門頭溝煤礦，埋深約410 m，X衍射分析結(jié)果顯示該巖樣為玄武巖，其內(nèi)部顆粒成分及體積分?jǐn)?shù)如表1所示。細(xì)觀試驗(yàn)觀測結(jié)果顯示，玄武巖試樣由礦物顆粒隨機(jī)分布并膠結(jié)而成，各礦物成分的彈性參數(shù)如表2所示。

表1 玄武巖X射線衍射結(jié)果[20]Table 1 X-ray diffraction results of basalt[20]

表2 礦物組分彈性參數(shù)[21]Table 2 Elastic parameters of mineral particles[21]

3.2 模型建立及彈性參數(shù)

根據(jù)礦物顆粒的物理力學(xué)性質(zhì)，將上述礦物成分分為兩類：①黏土和綠泥石視為基質(zhì)；②石英、斜長石，綠簾石視為礦物顆粒。根據(jù)Mori-Tanaka方法計(jì)算出復(fù)合后的顆粒與基質(zhì)及整個(gè)巖樣均勻化的彈性參數(shù)如表3所示。

實(shí)驗(yàn)中模型尺寸及邊界條件如圖10所示，其中C為偏置距離，預(yù)制缺口寬度為0.4 mm。為了減少計(jì)算量，數(shù)值模擬采用2D平面多尺度模型[23]，首先建立偏置距離2 mm的均質(zhì)化模型[圖11(a)]，研究了裂紋擴(kuò)展的可能區(qū)域[圖11(b)]，裂紋可能的擴(kuò)展路徑主要集中在加載點(diǎn)和預(yù)制缺口之間的區(qū)域。基于此，將模型分為主要區(qū)域和次要區(qū)域。主要區(qū)域?yàn)榱鸭y可能擴(kuò)展的區(qū)域，為了更詳細(xì)地研究裂紋擴(kuò)展行為，在該區(qū)域建立更能反映巖石真實(shí)細(xì)觀結(jié)構(gòu)的隨機(jī)礦物顆粒分布于基質(zhì)的模型，而在次要區(qū)域采用均勻化后的均質(zhì)化模型。圖12所示為偏置6 mm多尺度模型的幾何尺寸及邊界條件，中間為主要區(qū)域，兩邊為次要區(qū)域，主要區(qū)域的模型按照礦物組分的比例及力學(xué)性質(zhì)生成，次要區(qū)域參數(shù)為均質(zhì)后的參數(shù)。

表3 均勻化彈性參數(shù)[22]Table 3 The homogenized elastic parameters[22]

圖10 結(jié)構(gòu)幾何尺寸及邊界條件Fig.10 Geometry and boundary conditions of tested structure

圖11 均質(zhì)模型及裂紋路徑Fig.11 Homogeneous model and crack path

圖12 結(jié)構(gòu)幾何尺寸及邊界條件Fig.12 Geometry and boundary conditions of tested structure

3.3 數(shù)值模擬

以缺口偏置距離6 mm為例，采用三節(jié)點(diǎn)單元，網(wǎng)格總數(shù)36 000個(gè)，加密區(qū)網(wǎng)格尺寸h≈0.1 mm，根據(jù)3.3.2節(jié)已述，長度參數(shù)lc=0.2 mm。采用控制位移方式進(jìn)行均勻加載，總位移加載分600步完成，均勻位移增量為Δu=0.000 1 mm，總位移u=0.06 mm完成整個(gè)裂紋擴(kuò)展過程。參數(shù)臨界能量釋放率gc是相場法關(guān)鍵參數(shù)[24]，很難直接從試驗(yàn)測量，本文研究中選取偏置缺口6 mm室內(nèi)試驗(yàn)測定的荷載位移曲線進(jìn)行標(biāo)定，該室內(nèi)試驗(yàn)采用控制位移進(jìn)行加載，加載速率為1×10-4mm/s,詳細(xì)試驗(yàn)數(shù)據(jù)參見左建平等[21]。此外，根據(jù)左建平等[21, 25]所描述，裂紋擴(kuò)展大多是沿顆粒擴(kuò)展，且總是基于最小能量方式擴(kuò)展。因此，假設(shè)顆粒能量釋放率為基體的2倍，對于兩邊的次要區(qū)域，由于在均質(zhì)模型中沒有發(fā)生裂紋，所以假設(shè)次要區(qū)域的臨界能量釋放率遠(yuǎn)遠(yuǎn)大于主要區(qū)域的臨界能量釋放率，這樣可以讓裂紋更好地在基質(zhì)中擴(kuò)展。通過室內(nèi)試驗(yàn)數(shù)據(jù)標(biāo)定(圖13)，最終確定基體臨界能量釋放率gc為48×10-5kN/mm，即復(fù)合顆粒臨界能量釋放率gc=96×10-5kN/mm。

采用上述模型參數(shù)模擬了偏置距離為1、3、4 mm的三點(diǎn)彎試驗(yàn)，并與上述室內(nèi)試驗(yàn)中相對應(yīng)模型實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對比，如圖14所示，本文數(shù)值模擬結(jié)果與巖石三點(diǎn)彎實(shí)驗(yàn)的力與位移曲線擬合較好，對巖石的強(qiáng)度和峰后脆性破壞有較好的預(yù)測。因?qū)嶒?yàn)中巖石巖樣存在初始裂紋，導(dǎo)致實(shí)驗(yàn)曲線彈性階段具有一段壓密階段，需進(jìn)一步考慮巖石壓密階段，但關(guān)鍵點(diǎn)彈性模量、峰值強(qiáng)度基本吻合。未達(dá)到破壞荷載前，力與位移曲線呈線性關(guān)系，達(dá)到破壞荷載后，曲線開始跌落，直至試樣失效，曲線揭示巖石破壞過程為彈脆性破壞過程。從圖15看出數(shù)值模擬中的峰值強(qiáng)度變化趨勢與實(shí)驗(yàn)數(shù)據(jù)一致，隨著偏置缺口距離的增加，巖石的破壞荷載也逐漸增加。

圖13 能量釋放率的標(biāo)定Fig.13 Calibration of the parameter

圖14 力與位移曲線與其試驗(yàn)結(jié)果對比Fig.14 Contrast force and displacement curves with experiments

圖15 峰值荷載與偏置距離關(guān)系Fig.15 Relationship between peak load and bias distance

圖16 偏置6 mm時(shí)裂紋擴(kuò)展Fig.16 Crack propagation for the 6 mm offset case

3.4 裂紋擴(kuò)展的影響因素分析

以缺口偏置距離6 mm巖樣的裂紋擴(kuò)展數(shù)值模擬為例，裂紋擴(kuò)展的整個(gè)過程可分為4個(gè)階段，如圖16所示：①當(dāng)位移u=0.037 mm時(shí)，缺口開始出現(xiàn)損傷；②位移增加到位移u=0.040 mm時(shí)，缺口處的損傷加深，裂紋萌生；③當(dāng)位移u=0.047 mm時(shí)，裂紋開始向頂部快速擴(kuò)展；④直至位移u=0.06 mm時(shí)裂紋擴(kuò)展至整個(gè)巖樣，導(dǎo)致巖樣失效。裂紋擴(kuò)展過程與圖13的力-位移曲線各個(gè)階段基本一致，但是，從細(xì)觀尺度而言，巖石中的顆粒大小、位置分布以及物理力學(xué)性質(zhì)都會(huì)影響到巖石中裂紋的擴(kuò)展行為。為此，本節(jié)將從細(xì)觀尺度詳細(xì)分析裂紋擴(kuò)展的影響因素。

3.4.1 礦物顆粒分布的影響分析

巖樣中的礦物顆粒是隨機(jī)分布的，為了進(jìn)一步研究礦物顆粒不同分布對裂紋擴(kuò)展的影響規(guī)律，在礦物組分體積分?jǐn)?shù)一定，排除礦物顆粒尺寸過大過小而引起的離散性，建立了礦物顆粒尺寸一致的3種隨機(jī)分布模型，裂紋擴(kuò)展路徑如圖17所示。模擬結(jié)果可以明顯看出，受礦物顆粒分布隨機(jī)性的影響，3個(gè)巖樣模型的裂紋擴(kuò)展路徑各有不同，但大方向仍是從缺口處萌生，向頂部擴(kuò)展。力與位移曲線如圖18所示，3條曲線在彈性階段完全吻合，從裂縫萌生開始出現(xiàn)不同，裂縫萌生的起始位置和后期擴(kuò)展過程均表現(xiàn)出一定的差異，但相差不大。因此，在礦物組分比例一定，顆粒尺寸一致的情況下，顆粒的分布會(huì)改變裂紋擴(kuò)展路徑而對力-位移曲線僅有較小影響。

圖17 3個(gè)不同模型下裂紋路徑：相同的顆粒體積分?jǐn)?shù)及大小，不同顆粒位置Fig.17 Three models of different crack paths: same particle volume fraction and size, different particle location

圖18 不同模型的力與位移曲線Fig.18 Force and displacement curves of different models

3.4.2 能量釋放率的影響分析

巖石中礦物的力學(xué)性質(zhì)，尤其是臨界能量釋放率的不同也會(huì)影響到巖石中裂紋的擴(kuò)展規(guī)律。為深入分析能量釋放率的影響規(guī)律，研究以下兩種情況，A：基質(zhì)能量釋放率小于顆粒能量釋放率；B：基質(zhì)能量釋放率大于顆粒能量釋放率。

圖19 不同gc下u=6 mm裂紋路徑圖Fig.19 Path of crack u=6 mm under different gc

圖20 不同能量釋放率的力與位移曲線Fig.20 Force and displacement curves of different energy release rates

針對第一種情況，在主要區(qū)域中基質(zhì)選取4個(gè)不同的能量釋放率，A1:gc=48×10-5kN/mm；A2：gc=38×10-5kN/mm；A3：gc=28×10-5kN/mm，A4：gc=18×10-5kN/mm，顆粒能量釋放率保持gc=98×10-5kN/mm不變，以偏置距離2 mm模型為例，最終裂紋形態(tài)如圖19所示。為更好地研究裂紋擴(kuò)展路徑的規(guī)律，顯現(xiàn)出部分顆粒。可以看出，裂紋增長時(shí)均繞過礦物顆粒在基質(zhì)中擴(kuò)展，與實(shí)驗(yàn)觀察到的沿晶擴(kuò)展一致。此外，可以觀察到，由于4種基質(zhì)的能量釋放率不同，裂紋在基質(zhì)中的擴(kuò)展形態(tài)也各有差異，影響了裂紋的最終形態(tài)。Zhou等[26]研究表明，均質(zhì)模型中能量釋放率不會(huì)改變裂紋的最終形態(tài)，而本文在考慮了礦物顆粒后，裂紋路徑在不同能量釋放率時(shí)最終形態(tài)發(fā)生改變，意味著裂紋擴(kuò)展行為受巖石內(nèi)部細(xì)觀結(jié)構(gòu)影響。圖20所示的4種情況下的力與位移曲線趨勢幾乎一致，不同能量釋放率下曲線斜率相同，峰值隨gc的增加而增加。

針對第二種情況，保持基質(zhì)能量釋放率gc=48×10-5kN/mm不變，選取兩種不同的顆粒能量釋放率，B1：gc=48×10-5kN/mm，B2：gc=8×10-5kN/mm。同樣以偏置距2 mm模型為例，最終裂紋形態(tài)如圖21所示。結(jié)果表明，隨著顆粒能量釋放率的降低，裂紋擴(kuò)展路徑逐漸向顆粒內(nèi)擴(kuò)展，不再沿顆粒邊界在基質(zhì)中擴(kuò)展，而出現(xiàn)穿晶擴(kuò)展。可見，裂紋擴(kuò)展路徑總是趨于最小耗能形式。巖石中礦物顆粒能量釋放率的不同會(huì)影響裂紋擴(kuò)展的速度、最大荷載以及裂紋的形態(tài)。

圖21 不同gc下裂紋路徑圖Fig.21 Path of crack under different gc

4 結(jié)論

(1)考慮了巖石內(nèi)部初始細(xì)觀結(jié)構(gòu)，在裂紋擴(kuò)展區(qū)建立了具有隨機(jī)分布礦物顆粒的多尺度模型，在一定程度上更合理地描述了巖石內(nèi)部細(xì)觀結(jié)構(gòu)。

(2)利用相場法從細(xì)觀尺度預(yù)測和分析了巖石三點(diǎn)彎曲試驗(yàn)過程中的裂紋起裂、擴(kuò)展和貫通的破壞過程，分析了細(xì)觀結(jié)構(gòu)分布和力學(xué)特征對裂紋擴(kuò)展行為的影響規(guī)律。

(3)與室內(nèi)宏細(xì)觀試驗(yàn)對比結(jié)果表明，相場法能很好地模擬和預(yù)測巖石內(nèi)部裂紋萌生、擴(kuò)展和貫通的破壞過程，具有一定的應(yīng)用前景。在未來的研究中，可以進(jìn)一步向三維問題和多場耦合問題拓展。