真三軸壓縮下大理巖強度、變形與損傷特征數值

王浩然，王志亮，王昊辰，汪書敏

(合肥工業大學 土木與水利工程學院，合肥 230009)

自然界中巖體處于復雜應力狀態，在進行深部洞室開挖過程中，臨空面會釋放巖體中積累的應變能，造成應力重分布，導致巖體內部產生新裂紋，并且裂紋與裂紋之間擴展貫通，最終會造成巖體的失穩破壞[1-3]。以往的常規三軸壓縮試驗忽略了中間主應力(σ2)的影響，傳統的莫爾-庫侖和霍克-布朗準則對真三軸應力狀態下巖石力學表征難以適用，而中間主應力對巖石在真三軸應力狀態下的力學特性影響顯著[4]。近年來，許多學者對深埋大理巖進行真三軸壓縮試驗，以探究其在高地應力下的破壞機制。Zheng等[4]探討了錦屏大理巖在真三軸壓縮下的殘余強度特征，發現對于較高σ3(最小主應力)，隨著σ2的增大殘余強度先減小后增大，對于較低σ3，則隨著σ2的增大而逐漸減小；Feng等[5]對花崗巖、大理巖和砂巖等進行了真三軸循環加卸載試驗，指出巖石黏聚力對σ3和σ2不敏感，內摩擦角最小值不受σ3和σ2的影響，但其最大值受兩者的影響；Gao等[6]對含節理的錦屏大理巖進行真三軸壓縮試驗，發現節理大理巖強度對σ3的敏感性要大于對σ2的敏感性。

然而，在實驗室進行真三軸試驗需要先進且昂貴的實驗設備，試件的制備也費時耗力，并且試驗中很難觀察到細觀的破壞過程。近年來，數值模擬方法已成為傳統研究方法的有力補充。其中，離散單元法作為一個有效的數值分析工具，廣泛應用于巖石力學研究[7-8]。Yang等[9]研究了射孔角度和射孔布置對同步水力壓裂裂縫擴展機制和巖石破壞模式的影響，指出大主應力差下較小的射孔間距更有利于中心巖體的破裂。Lee等[10]對非平行雙裂隙花崗巖試樣進行單軸壓縮試驗，并且開展了顆粒流程序模擬，發現巖橋裂紋的貫通主要通過拉伸裂紋的擴展聚集。Zhang等[11]對含兩個平行裂隙的矩形類巖石試件進行單軸壓縮數值模擬，總結出巖橋貫通兩種新模式。由于實際工程問題復雜，簡化為二維問題后，其計算結果存在一定差距，故開展三維數值研究十分必要。Duan等[12]對砂巖展開真三軸壓縮離散元模擬，重點分析σ3、σ2對微裂紋取向、黏結力取向、組構和配位數演化的影響，指出隨著σ3水平的增加，σ2對顆粒尺度響應的影響逐漸減弱。Zhang等[13]采用平行黏結模型先標定大理巖的力學性能，然后建立真三軸壓縮數值模型，在數值巖樣上預制裂隙，研究裂隙大理巖在真三軸應力狀態下的破壞特征。Zhang等[14]采用黏結粒子模型研究花崗巖在真三軸應力狀態下的變形和強度特征，消除了端部效應，隔離了中間主應力的影響。

綜上，中間主應力對巖石力學特性影響較大，且目前對真三軸應力環境下的巖石內部裂隙擴展演化過程難以做到實時觀察。為此，擬采用三維離散元顆粒流程序(PFC3D)對大理巖試件開展真三軸應力狀態下的微觀裂紋擴展機制數值研究，并基于全應力-應變曲線特征，深入探析裂紋擴展過程與受壓變形過程當中的內在聯系，同時與實驗觀測結果進行比較，力求得出具有參考價值的結論。

1 數值模型的建立

1.1 數值模型

PFC5.0內嵌有十種接觸模型[8]，其中，平行黏結模型在模擬巖石類材料受壓破壞時，能夠沿著法向或者切向破壞，可較好地反映此類材料的破壞形式，已廣泛應用于巖石類材料相關問題分析中[15]，故本文選擇該模型進行后繼研究。大理巖試樣數值模型尺寸為50 mm×50 mm×100 mm，假定無內部缺陷。采用與文獻[16]常規三軸壓縮數值試驗相同的平行黏結模型細觀參數，大理巖顆粒的半徑在0.9～1.356 mm，避免顆粒尺寸對力學性能的影響。其中，真三軸數值模型包含26 018個顆粒，112 093個平行鍵黏結，4 504個線性黏結，如圖1所示。

圖1 數值建模

1.2 真三軸壓縮模擬加載方式

數值模擬中的加載應力路徑是仿照巖石真三軸壓縮試驗的加載方式[17-18]，如圖2所示，模擬的加載路徑由以下3個階段組成：

圖2 真三軸試驗的加載路徑

1)在靜水條件下(σ1=σ2=σ3)，通過伺服控制函數控制6個墻體以一定位移速率單調加載試樣，直到達到預設的σ3水平。

2)通過伺服控制函數控制x方向上兩個平行墻體保持σ3恒定，同時控制y、z方向兩對平行墻體以恒定位移速率繼續加載，保證σ1和σ2以相同的速率單調上升，直到達到預定的σ2。

3)通過伺服控制函數控制x、y方向上兩對平行墻體保持σ2、σ3恒定，控制z方向平行墻體保證σ1單調升高，直至峰值(σ1,peak)，然后繼續加載至σ1后下降到σ1,peak的70%時終止。

1.3 真三軸模擬應力路徑及破壞角的測量

對大理巖數值試樣分別進行7組不同應力控制的真三軸壓縮破壞試驗，即將最小主應力σ3分為5、20、50、80、100、120和150 MPa 7個應力水平，每個應力水平下，σ2的大小根據模擬過程中峰值強度確定，涵蓋了從σ2=σ3到σ2=σ1的整個范圍(表1)。

表1 大理巖真三軸模擬的破壞應力狀態和破壞面角

圖3展示出破壞面角度θ的測量方法，左邊為顆粒破壞后的fragment顯示(同種顏色的顆粒表示同種破碎塊體)，中間為微裂紋的空間分布(深色為剪切裂紋，淺色為拉伸裂紋)，右邊為顆粒位移矢量空間分布。當σ2=σ3時沒有明顯的破壞角；當σ2>σ3時，平行黏結鍵斷裂，無法傳遞力，故無法形成力鏈，則在平行黏結鍵斷裂處產生微裂紋，微裂紋的積累導致破裂面的產生。而且，破壞面平行于中間主應力方向，各個應力路徑下破壞角的范圍為60°～65°(見圖3(b)、表1)，對于恒定σ3，破壞角隨著σ2的增大而增大。

圖3 PFC計算結果及破壞角的測量

2 真三軸壓縮數值結果分析

2.1 中間主應力對真三軸破壞應力的影響

表1列出了破壞時的所有真三軸壓縮應力條件，對于每個最小主應力σ3，σ1,peak是σ2的函數，如圖4所示。可以看出，σ1,peak隨σ2的增加而增大，且當σ2達到某一固定平臺值時，超過此值σ1,peak逐漸下降，這種變化趨勢不是左右對稱的，相關硬巖真三軸壓縮試驗結果一致[19-20]。

由上可知，σ1,peak既是σ3的函數，又是σ2的函數。σ1,peak這種先升后降的趨勢符合二階多項式方程特征。曲線圖4中的虛線組成常規三軸壓縮(σ2=σ3)和常規三軸拉伸(σ2=σ1)兩種極限情況，σ1,peak與σ3的變化可以用下式來線性擬合，這實為傳統的三軸破壞判據，即Mohr-Coulomb準則。

圖4 σ1,peak與σ2在恒定σ3水平上的變化關系

對于常規三軸壓縮(σ2=σ3)

σ1,peak=150.84+3.68σ3,R2=0.999

(1)

對于常規三軸拉伸(σ1=σ2)

σ1,peak=169.28+5.40σ3,R2=0.999

(2)

對于給定的σ3，σ1,peak與σ2的關系曲線表現出強度上的差異性，這揭示了傳統的三軸壓縮試驗以及經典莫爾-庫侖和霍克-布朗等破壞準則存在不足。

取峰值應力50%處的割線模量計算巖樣的彈性模量E，圖5展示E在σ3應力水平下隨σ2的變化關系，可以看出，在每個最小主應力水平下，彈性模量隨著中間主應力的增加先急速增加，再逐漸趨近于平穩，并且最小主應力越大，平穩增加段所占比例就越大，表明最小主應力對彈性模量的增加有約束作用。采用Logistic函數可直觀表示E和σ2的變化關系，擬合方程列于圖5中，擬合度較高。

圖5 E與σ2在恒定σ3水平上的變化關系

2.2 基于八面體理論的真三軸破壞應力分析

八面體理論認為，當巖石內部八面體剪應力達到某一臨界值時，巖石將進入破壞階段。為了研究大理巖試樣的空間破壞特性與中間主應力之間的關系，并且表1和圖4真三軸破壞時的應力狀態都可以由兩個主應力不變量之間的單一關系表示，即破壞時的八面體剪應力(τoct,f)和破壞時的八面體正應力(σoct,f)：

τoct,f=[(σ1,peak-σ2)2+(σ2-σ3)2+

(σ3-σ1,peak)2]1/2/3

(3)

σoct,f=(σ1,peak+σ2+σ3)/3

(4)

圖6(a)是在τoct,f-σoct,f區域中重新繪制破壞時的應力條件，可以看出，在大理巖數值試樣當中，τoct,f隨著σoct,f的升高而持續上升，盡管上升速率有所下降。因為所建數值模型沒有明顯缺陷，破壞時的主要原因是剪切而不是壓實[21]，與2.4節的結論相同。同時，為了更好地確定在恒定σ3應力條件下σ2在σ1和σ3之間的相對位置，引入應力比b[17]，對于確定的破壞應力狀態，可以獲得良好的數據基礎。應力比b定義為

b=(σ2-σ3)/(σ1-σ3)

(5)

對圖6(a)中所示的數據點進行二階多項式方程擬合，得到式(6)的τoct,f-σoct,f變化方程。

(6)

由圖6(b)可以看出，在τoct,f-σoct,f曲線區域中數據點表現出一定的分散性，但是分散中包含潛在規律，即對于每個級別的τoct,f，最小的σoct,f一致地出現在σ2=σ3時，最大的σoct,f一致地匹配σ2=σ1應力狀態。

在圖6(c)中，重新繪制了與圖6(a)和(b)相同的應力狀態數據點，從式(5)可看出,σ2=σ3和σ2=σ1時，應力比b=0和b=1。與此同時，圖6(c)為每個數據點分配了與其表示的b值相對應的顏色(側面顯示的顏色條)，易見低σ3下個體b的離散度非常低，σ2=σ3和σ2=σ1的擬合曲線(圖中虛線)明確地描繪b=0和b=1時的應力狀態，即分別對應巖石處于廣義三軸壓縮狀態(σ1>σ2=σ3)和廣義三軸拉伸狀態(σ1=σ2>σ3)。對比圖6(a)和(c)可以看出，對于0到1的任何其他b值，都可以得到一條擬合良好的曲線，并且對于恒定的b值，τoct,f隨著σoct,f線性增加，形成線性的破壞包絡線。

圖6 破壞時τoct,f隨σoct,f的變化

2.3 變形特征與裂紋擴展的聯系

圖7給出了σ3=150 MPa下的偏應力與3個主應變分量的關系。限于篇幅，僅討論一個最小主應力水平下的關系圖。易見隨σ2/σ1的增大，應變εy由拉應變轉為壓應變，并且壓應變也隨之增大。增加σ2/σ1值也會導致在最小應力方向上位移膨脹量的提升。在第1個拐點之前，3個主應力、主應變大小相同，壓縮方向一致且為正；在第1個拐點與第2個拐點之間，最大主應力和中間主應力大小相同，εz和εy由于兩個方向的壓縮而一致且為正，εx由于在該方向上開始擴張而為負；第2個拐點之后，εy與εz由重合變為分離，兩個方向上的應變變化較為明顯。

圖7 不同應力路徑的偏應力-應變關系

圖8展示出σ3在80 MPa下微裂紋數量、每10時步所產生的微裂紋數和偏應力與軸向應變的關系。對比常規三軸壓縮偏應力-應變曲線，真三軸大理巖破壞微裂紋的演化過程也可分為4個階段：第1個階段為線彈性階段(OB段)，其中3個主應力伺服至σ3時到達A點，數值試樣不斷被壓縮，包含初始壓密區，幾乎沒有微裂紋和聲發射事件產生，B點所對應的偏應力在40%σ1,peak附近；第2個階段為裂紋穩定擴展階段(BC段)，試樣聲發射事件呈穩定增長態勢，巖石損傷呈穩定速率上升，C點所對應的偏應力在80%σ1,peak附近；第3個階段為裂紋非穩定擴展階段(CD段)，D點對應峰值強度，從C點新生裂紋張開度明顯增大，巖樣出現不可恢復的塑性變形，并且隨著σ2增加該階段占比增大；第4個階段為峰后破壞階段(DE段)，D點之后，不同應力條件下微裂紋的擴展速率不盡相同，但是隨著宏觀裂隙的形成，該階段后期的微裂紋擴展速率逐漸降低。

圖8 軸向應變與偏應力、微裂紋數和聲發射的關系

從圖8可以看出，隨著σ2的增大，試樣峰后表現出由塑性破壞轉向脆性破壞的趨勢。當σ3保持不變時，隨著σ2增大CD階段所占比重增加，并且剪切裂紋數目慢慢大于拉伸裂紋數目，破壞方式從拉伸破壞轉變為拉剪混合破壞。對于恒定的σ3，不同應力路徑下σ2的大小影響微裂紋的產生速度，即損傷的變化速率，且σ2越大則峰后損傷速率下降得越快。

2.4 微裂紋形成對宏觀破裂面的影響

圖9為不同應力路徑下微觀裂紋傾向傾角與其相對數量關系的赤平極射投影。對比圖9(a)與(d)可看出，初始圍壓對微觀裂紋的傾向傾角有較大影響：初始圍壓較低時微觀裂紋傾向分布較均勻，裂紋主要平行于最大主應力方向，類似單軸壓縮時裂紋的赤平極射投影[22]；σ3=σ2=80 MPa屬于初始圍壓較高的狀態，裂紋的傾向開始往中間主應力方向(90°)和最小主應力方向(0°)上集中，并且裂紋的傾角慢慢變大。從圖9(b)和(c)可以看出，因σ2>σ3時出現明顯的破壞角，裂紋傾向以及數量分布從中間主應力方向逐漸向最小主應力方向擴展，并且隨著中間主應力的增大，裂紋傾角也慢慢偏離最大主應力方向，最終擴展的角度形成實際的宏觀破壞角(位于60°～65°)。

圖9 細觀裂紋赤平極射投影

從圖10(a)可以看出，對于低σ3應力水平，拉伸裂紋占比較大；當σ3>50 MPa之后，隨著σ2增大，剪切裂紋占比開始增大，破壞方式從拉伸破壞向拉剪混合破壞轉變。在圖10(b)和(c)中，對于恒定的σ3，微裂紋數目隨σ2先下降然后再上升，巖樣損傷發展呈現出“對勾”型變化趨勢，即開始下降快，而上升速度滯緩。當σ3=σ2時，隨著σ3的增大，微裂紋數目快速上升；而當σ1=σ2時，隨著σ3的增大，微裂紋數目上升速率較為緩慢，后期微裂紋數目幾乎沒有變化。從這兩圖中的虛線可以看出高中間主應力抑制損傷的上升趨勢。

圖10 峰后70%σ1,peak時微裂紋數目與σ2變化關系

圖11對比了不同真三軸應力路徑下大理巖數值模擬破壞形態與試驗后的宏觀破壞形態[4]。主裂紋為剪切破壞，次生裂紋為拉伸破壞，二者吻合較好，且與文獻[5-6]中大理巖的破壞形態基本相同。此外，圖9與10的結果也體現了大理巖真三軸宏觀破壞方式，故本文所建模型是合理的。

圖11 宏觀破裂與模擬破壞形態對比

3 結 論

1)基于平行黏結模型的三維PFC細觀數值模擬，能較好地從微觀角度觀察大理巖的破壞模式，中間主應力σ2對巖樣的力學響應影響顯著，且微裂紋的赤平極射投影可動態展示微裂紋的傾向傾角和數量分布規律。

2)八面體剪切應力與八面體正應力可準確地擬合出此大理巖真三軸壓縮條件下的破壞應力，隨著σ2的增加，大理巖試樣破壞模式由拉伸破壞向拉剪破壞轉變，并且應力-應變曲線的峰后段由塑性狀態逐漸向脆性狀態過渡。

3)隨著σ2的增加，偏應力與3個主應變曲線將出現明顯拐點，中間主應變由拉伸向壓縮轉變；根據應力-應變曲線形狀，可將微裂紋演化過程分為4個階段，巖石損傷演化受σ2影響顯著，呈現出“√”型的變化趨勢。