王 猛，宋子楓，鄭冬杰，神文龍，勾攀峰,韋四江,2

(1.河南理工大學(xué) 能源科學(xué)與工程學(xué)院，河南 焦作 454000; 2.中國(guó)礦業(yè)大學(xué) 深部煤炭資源開采教育部重點(diǎn)實(shí)驗(yàn)室，徐州 江蘇 221116; 3.煤炭安全生產(chǎn)河南省協(xié)同創(chuàng)新中心，河南 焦作 454000)

準(zhǔn)確描述巖石強(qiáng)度與變形是評(píng)價(jià)工程巖體安全和穩(wěn)定性的基礎(chǔ)，現(xiàn)有應(yīng)力強(qiáng)度理論與破壞準(zhǔn)則很難有效描述巖石復(fù)雜的強(qiáng)度變化與整體破壞行為[1]。巖石破壞是能量驅(qū)動(dòng)下的一種狀態(tài)失穩(wěn)現(xiàn)象[2]，受載巖石變形過程伴隨能量演化，峰前以能量積聚為主，峰后主要表現(xiàn)為能量耗散與釋放，巖石變形破壞可看成是不同形式能量相互轉(zhuǎn)換的結(jié)果。因此，從能量角度研究巖石變形和失穩(wěn)，更接近其破壞本質(zhì)[3-5]。

基于能量守恒理論，國(guó)內(nèi)外學(xué)者對(duì)巖石加卸載過程能量演化規(guī)律開展了大量研究，取得了豐富的研究成果。具有代表性的有:謝和平等[6-8]闡述了巖石變形破壞過程能量耗散與損傷、能量釋放與整體破壞的關(guān)系，定義了基于能量耗散與釋放原理的巖石強(qiáng)度失效準(zhǔn)則和巖體整體破壞準(zhǔn)則。張志鎮(zhèn)等[9-11]研究了不同巖石類型能量演化的異同點(diǎn)，建立了巖石能量演化隨加載條件的自我抑制模型，揭示了巖石能量耗散的圍壓效應(yīng)。宮鳳強(qiáng)、李淼等[12-13]分別研究了動(dòng)靜態(tài)劈裂過程巖石能量構(gòu)成與耗散特征，揭示了動(dòng)靜荷載作用下巖石能量耗散機(jī)制及各向異性特征。與此同時(shí)，許江等[14]、楊磊等[15]、尹光志和蘇國(guó)韶[16-17]、李楊楊等[18]分別研究了含孔隙水砂巖、煤巖組合體、不同加載速率以及循環(huán)加卸載等條件下巖石能量積聚與耗散特征，豐富了巖石能量理論的研究成果。

現(xiàn)有研究常聚焦于室內(nèi)試驗(yàn)，多是針對(duì)特定加卸載條件、巖樣類型以及其他條件下(含水率、溫度[19]等)巖石能量演化特征的描述，研究對(duì)象著重關(guān)注巖石屈服、峰值等特征點(diǎn)，對(duì)于貫穿巖石變形破壞全過程的能量實(shí)時(shí)演化及分配規(guī)律認(rèn)識(shí)不足[9]。巖石能量積聚、耗散及其狀態(tài)間的轉(zhuǎn)化決定了其內(nèi)部裂隙擴(kuò)展和變形破壞，雖然從能量角度反演巖石破壞路徑對(duì)實(shí)驗(yàn)結(jié)果意義不大，但對(duì)于評(píng)估工程巖體區(qū)域穩(wěn)定、精準(zhǔn)定位加固位置至關(guān)重要。目前，實(shí)驗(yàn)室常用巖石破壞監(jiān)測(cè)手段有聲發(fā)射定位[20],CT掃描[21],SEM電鏡掃描[22],CCD相機(jī)[23]等，其中，僅有聲發(fā)射可以實(shí)時(shí)定位受載巖樣變形破壞路徑，但該方法多用于單軸壓縮試驗(yàn)，且存在操作復(fù)雜、誤差大等缺點(diǎn)，很難推廣應(yīng)用于工程巖體的監(jiān)測(cè)。分析工程巖體穩(wěn)定性時(shí)多是基于圍巖應(yīng)力場(chǎng)、位移場(chǎng)及支護(hù)場(chǎng)的耦合作用關(guān)系進(jìn)行綜合評(píng)判。

為增加能量理論模型分析工程問題時(shí)的適用性，筆者基于能量守恒和有限差分理論，推導(dǎo)巖石能量耗散有限差分方程式，采用FISH語(yǔ)言將其應(yīng)用于FLAC3D應(yīng)變軟化本構(gòu)模型。以室內(nèi)試驗(yàn)結(jié)果驗(yàn)證采用耗散能演化描述受載巖石變形破壞的合理性，進(jìn)一步討論能量耗散模型在深部高應(yīng)力巷道中的應(yīng)用，揭示工程巖體峰后變形破壞全過程耗散能演化特征與分配規(guī)律，為從能量耗散角度評(píng)估工程巖體安全和穩(wěn)定性提供理論依據(jù)，彌補(bǔ)FLAC3D軟件采用塑性區(qū)定性判斷地下工程圍巖破壞時(shí)的不足。

1 巖石能量耗散的室內(nèi)實(shí)驗(yàn)

1.1 巖石單、三軸壓縮試驗(yàn)

試驗(yàn)巖樣取自徐礦集團(tuán)三河尖煤礦吳莊區(qū)運(yùn)輸大巷，選取巖性均勻、結(jié)構(gòu)完整的取芯巖柱加工巖樣，采用GCTS-RTX3000巖石力學(xué)試驗(yàn)系統(tǒng)進(jìn)行單軸和常規(guī)三軸壓縮試驗(yàn)，測(cè)試巖石基本力學(xué)參數(shù)，圍壓等級(jí)設(shè)置為0,5,15和25 MPa，測(cè)試結(jié)果分別如圖1和表1所示。

表1 常規(guī)單、三軸試驗(yàn)結(jié)果Table 1 Results of uniaxial and triaxial tests

圖1 巖石應(yīng)力-應(yīng)變曲線Fig.1 Stress-strain curves of the rock samples

1.2 巖石能量耗散特征

巖樣在試驗(yàn)機(jī)外力作用下產(chǎn)生變形，假設(shè)該物理過程與外界沒有熱交換，即封閉系統(tǒng)，外力功所產(chǎn)生的總輸入能量為Wz，根據(jù)熱力學(xué)第一定律[6]得到

Wz=We+Wd

(1)

式中，We為巖石可釋放彈性應(yīng)變能;Wd為耗散能。

由式(1)可知，巖石總輸入能量Wz一部分作為可釋放的彈性能儲(chǔ)存在其內(nèi)部，另一部分則伴隨巖石變形破壞耗散和釋放，此過程不可逆[8-9]。巖石損傷程度越高，耗散能越大，殘余彈性能越少。

對(duì)于單軸壓縮試驗(yàn)，試樣總輸入能量Wz等于軸向輸入能量Wa;而常規(guī)三軸壓縮試驗(yàn)，軸向輸入能量Wa應(yīng)包含試樣總輸入能量Wz以及環(huán)向膨脹對(duì)液壓油做功所釋放能量Wh兩部分，為此，常三軸壓縮條件下巖石總輸入能量[6]可表示為

(2)

式中，ε1，ε3分別為軸向應(yīng)變和環(huán)向應(yīng)變;σ1，σ3分別為軸向應(yīng)力和圍壓。

根據(jù)彈性力學(xué)相關(guān)理論，巖石可釋放彈性能We[6]可由下式計(jì)算:

(3)

則巖石耗散能Wd可由下式求得

Wd=Wz-We

(4)

采用式(2)～(4)分別計(jì)算巖樣總輸入能量Wz、彈性能We和耗散能Wd，如圖2所示。由圖2可知，巖樣總輸入能量隨軸向應(yīng)變?cè)黾又饾u增大，應(yīng)力達(dá)到峰值前主要以彈性能形式儲(chǔ)存在巖樣內(nèi)部，耗散能變化不大。當(dāng)巖樣加載至峰后階段，巖樣總輸入能量依然持續(xù)增加，但主要伴隨峰后破壞而耗散，此外，巖樣內(nèi)部?jī)?chǔ)存的大部分彈性能也將轉(zhuǎn)化成耗散能釋放，能量轉(zhuǎn)化程度與巖樣變形破壞程度正相關(guān)。

圖2 室內(nèi)試驗(yàn)巖樣能量演化規(guī)律Fig.2 Evolution law of energy of rock samples in laboratory test

2 FLAC3D能量耗散模型的開發(fā)

σ1<σ2<σ3

(5)

相應(yīng)的主應(yīng)變?cè)隽喀う舏[25]可分解為

(6)

單元體彈性應(yīng)變?cè)隽亢蛻?yīng)力增量之間的關(guān)系[25]可表示為

(7)

(8)

式中，α1,α2為切變模量G和體積模量K控制的巖石材料常數(shù)，其中，α1=K+4G/3，α2=K-2G/3。

基于Mohr-Coulomb的應(yīng)變軟化模型在峰后軟化過程控制強(qiáng)度參數(shù)描述材料的峰后軟化行為，對(duì)于剪切破壞，非關(guān)聯(lián)流動(dòng)法則[25]改寫為

(9)

式中，λs為塑乘因子;gs為剪切勢(shì)函數(shù)，其表達(dá)式為

(10)

式中，ψ為剪脹角。

因此，聯(lián)立式(6)，(7)和(9)，由總應(yīng)變表示的應(yīng)力增量Δσi表達(dá)式[25]為

(12)

(13)

式中，φ為內(nèi)摩擦角。

根據(jù)FLAC3D應(yīng)變軟化模型屈服準(zhǔn)則，巖體材料剪切破壞面上應(yīng)力[25]滿足:

(14)

式中，c為黏聚力。

拉伸破壞面上應(yīng)力[25]滿足:

ft=σ3-σt=0

(15)

式中，σt為抗拉強(qiáng)度。

式(14)和(15)分別為FLAC3D應(yīng)變軟化模型判定材料壓剪破壞和張拉破壞的失穩(wěn)準(zhǔn)則。

圖3 單元能量計(jì)算簡(jiǎn)圖Fig.3 Schematic diagram for the energy calculation

(16)

假設(shè)t時(shí)刻對(duì)應(yīng)第n循環(huán)，t+Δt時(shí)刻對(duì)應(yīng)第n+1循環(huán)，則第n+1循環(huán)內(nèi)總能量增量可表示為

(17)

第n循環(huán)至第n+1循環(huán)運(yùn)算的物理意義表示巖體受外力作用時(shí)，某一運(yùn)算時(shí)步內(nèi)單元網(wǎng)格應(yīng)力-應(yīng)變行為的差分運(yùn)算。t+Δt時(shí)刻單元總能量等于n+1循環(huán)總能量增量累加，即

(18)

不考慮峰后彈性模量衰減，單元可釋放彈性應(yīng)變能可由下式計(jì)算:

(19)

則單元耗散能計(jì)算式如下:

(20)

圖4 單元能量計(jì)算流程Fig.4 Energy calculation flow for the model cells

3 模型校驗(yàn)與運(yùn)算

3.1 巖樣模型建立

建立標(biāo)準(zhǔn)巖樣數(shù)值模型，尺寸φ×H=50 mm×100 mm，共劃分17 280單元，如圖5所示。模型底部固定位移邊界，四周施加環(huán)向邊界應(yīng)力模擬圍壓，上邊界通過施加軸向位移對(duì)巖樣加壓，加載速率2.5×10-5mm/step，計(jì)算時(shí)忽略巖石自重影響。

模型采用基于Mohr-Coulomb的應(yīng)變軟化本構(gòu)模型，以表1給出的巖樣力學(xué)參數(shù)對(duì)模型進(jìn)行初始賦值。巖石峰后軟化階段黏聚力和內(nèi)摩擦角隨應(yīng)變的變化數(shù)據(jù)由試驗(yàn)結(jié)果計(jì)算得到，見表2。

表2 巖石峰后軟化參數(shù)Table 2 Post-peak parameters of rock samples

依據(jù)圖5建立的數(shù)值模型，表1和2所列巖石力學(xué)參數(shù)進(jìn)行數(shù)值分析，圖6給出了模擬巖樣應(yīng)力應(yīng)變曲線與試驗(yàn)結(jié)果的對(duì)比，由圖6可知，不同圍壓下巖樣應(yīng)力應(yīng)變曲線模擬結(jié)果較為匹配試驗(yàn)結(jié)果，單軸壓縮和低圍壓條件下巖樣脆性特征得到較好描述，表明所選取的巖石參數(shù)和數(shù)值本構(gòu)模型合理，可用于后續(xù)能量模塊的模擬分析。

圖5 巖樣數(shù)值模型Fig.5 Numerical model of rock samples

圖6 應(yīng)力-應(yīng)變曲線校驗(yàn)Fig.6 Verifications of the stress-strain curves

3.2 巖樣能量模型校驗(yàn)

考慮到耗散能是由巖石總能量減去彈性應(yīng)變能計(jì)算得到，限于篇幅，校驗(yàn)?zāi)芰磕Ｐ蜁r(shí)只校驗(yàn)總能量和耗散能兩個(gè)指標(biāo)。如圖7(a)所示，不同圍壓下巖樣總能量數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果基本吻合，表明FISH語(yǔ)言編寫的能量算法可靠性較強(qiáng)。

但需要指出的是，模擬巖樣峰后耗散能演化曲線與試驗(yàn)結(jié)果吻合度較好，但峰前存在差異性，如圖7(b)所示。峰前階段，巖樣試驗(yàn)耗散能主要呈緩慢線性增加;加載至峰值80%時(shí)，耗散能近似指數(shù)增加;一旦到達(dá)峰后某破壞點(diǎn)，耗散能演化為斜率較大的直線增長(zhǎng)，最終進(jìn)入殘余階段。數(shù)值模擬可近似反演耗散能演化趨勢(shì)，耗散能在峰后增長(zhǎng)曲線與試驗(yàn)結(jié)果匹配較好，且總量大致相等;但峰前階段存在顯著差異，模擬峰前耗散能近似等于0，主要原因是能量模型峰前假設(shè)不發(fā)生塑性變形，無(wú)法準(zhǔn)確描述孔隙、松軟巖體初始?jí)好茈A段和峰前塑性階段，導(dǎo)致圍壓越大，峰前耗散能相差越大。但是，對(duì)于地下工程巖體的分析更關(guān)注于峰后變形破壞，由于該模型可較好模擬巖石峰后能量耗散特征，可用于巖石峰后變形破壞的模擬。

圖7 能量模型校驗(yàn)Fig.7 Verifications of the energy model

3.3 巖樣峰后破壞路徑反演

圖8給出了不同圍壓下巖石破壞模式對(duì)比，其中，數(shù)值模擬給出的是巖石殘余階段耗散能密度演化云圖，對(duì)應(yīng)試驗(yàn)巖樣殘余階段破壞照片。由圖8可知，巖樣單軸壓縮試驗(yàn)呈縱向劈裂破壞，模擬巖樣耗散能集中區(qū)域位于巖樣中線偏上位置，并向周邊演化出一“八”字分支集中線，巖樣耗散能演化沿中線呈軸對(duì)稱分布，反映了巖樣單軸壓縮時(shí)的縱向劈裂特征，與試驗(yàn)結(jié)果吻合度較高。隨著圍壓增加，試驗(yàn)巖樣逐漸向剪切破壞演化，破斷角逐漸減小;模擬巖樣在5 MPa時(shí)耗散能密度集中成“y”型分布，長(zhǎng)邊集中程度顯著大于短邊，為巖樣主破斷面;隨圍壓進(jìn)一步增加，其破壞模式演化成斜切面，且圍壓越大，耗散能越大，破斷角越小，與試驗(yàn)結(jié)果基本吻合，表明采用耗散能密度演化揭示巖石峰后破壞行為是可行的。

圖8 巖石破壞模式校驗(yàn)Fig.8 Verifications of the failure modes of rock samples

工程現(xiàn)場(chǎng)，準(zhǔn)確掌握圍巖變形破壞路徑，精準(zhǔn)定位加固位置是控制其穩(wěn)定的前提。為此，反演巖石峰后耗散能演化過程可為揭示其主破壞區(qū)演化、評(píng)判巖石穩(wěn)定性提供新途徑。圖9給出了圍壓25 MPa時(shí)不同峰后目標(biāo)點(diǎn)巖樣耗散能密度演化云圖。巖樣初始進(jìn)入峰后階段(B點(diǎn))，經(jīng)歷裂隙孕育形成初始主裂隙面，并隨著軸向應(yīng)變?cè)黾?，主裂隙面耗散能逐漸增加，一旦巖樣進(jìn)入殘余階段(F點(diǎn))，主裂隙面耗散能密度增幅減小并趨于穩(wěn)定，表明巖樣此時(shí)已發(fā)生宏觀破斷，繼續(xù)加載將導(dǎo)致其結(jié)構(gòu)失穩(wěn)。

圖9 圍壓25 MPa巖樣峰后破壞路徑Fig.9 Post-peak destruction path of rock sample with a confining pressure of 25 MPa

4 能量耗散模型的應(yīng)用

將能量耗散模型應(yīng)用于分析吳莊區(qū)運(yùn)輸大巷，嘗試從能量耗散角度揭示巷道變形破壞路徑，為巷道加固提供依據(jù)。建立試驗(yàn)巷道三維數(shù)值模型，尺寸X×Y×Z=60 m×40 m×60 m，水平、底邊界限制位移，上邊界施加20 MPa載荷模擬巷道800 m埋深，依據(jù)地應(yīng)力測(cè)試結(jié)果，取側(cè)壓系數(shù)0.8。巷道斷面寬×高=5 m×4 m，模擬巷道錨桿支護(hù)參數(shù)與現(xiàn)場(chǎng)一致。

模擬采用應(yīng)變軟化模型，由于巖石實(shí)驗(yàn)參數(shù)無(wú)法直接應(yīng)用于工程巖體，模擬前需校驗(yàn)巖體參數(shù)?；趯?shí)驗(yàn)獲取的巖石強(qiáng)度參數(shù)衰減規(guī)律，采用RocLab軟件進(jìn)行巖體參數(shù)轉(zhuǎn)換，以獲取的初始巖體參數(shù)進(jìn)行賦參，以實(shí)測(cè)巷道變形作為已知值，校驗(yàn)獲取的巖體參數(shù)見表3,4。同時(shí)，為防止單元畸變中斷運(yùn)算，模擬時(shí)未設(shè)置set large命令。由于巷道周邊圍巖受力多處于三向不等壓狀態(tài)，計(jì)算單元能量時(shí)考慮中間主應(yīng)力，各應(yīng)力分量增量計(jì)算方法見式(12)，再由式(17)～(20)分別計(jì)算單元總能量、彈性能和耗散能。

表3 校驗(yàn)巖體參數(shù)(峰前)Table 3 Calibrated parameters of the rock masses

圖10給出了試驗(yàn)巷道模擬結(jié)果與現(xiàn)場(chǎng)對(duì)比圖。試驗(yàn)巷道開挖后，采用錨網(wǎng)索支護(hù)，巷道頂幫變形得到控制，頂幫變形均在200 mm以內(nèi);但底臌嚴(yán)重，最大底臌量接近1 000 mm，模擬巷道變形與現(xiàn)場(chǎng)實(shí)測(cè)數(shù)據(jù)基本吻合。

表4 校驗(yàn)巖體參數(shù)(峰后)Table 4 Calibrated parameters of the rock masses

圖10 試驗(yàn)巷道變形破壞對(duì)比Fig.10 Comparisons on the roadway deformations

模擬巷道塑性區(qū)主要集中于兩幫和底板，基于軟件拉剪屈服破壞準(zhǔn)則，僅可以區(qū)分底幫圍巖變形破壞方式，但很難量化底幫圍巖的破壞程度，也就無(wú)法準(zhǔn)確定位巷道后期加固區(qū)域。而利用開發(fā)的能量耗散模型，從巷道耗散能集中區(qū)分布可實(shí)現(xiàn)對(duì)巷道圍巖主失穩(wěn)區(qū)域的實(shí)時(shí)定位，同時(shí)，依據(jù)圍巖耗散能集中程度可以量化巖體變形破壞程度。如圖10所示，巷道底板破壞后，底角位置分別存在兩個(gè)耗散能集中區(qū)，均呈條狀分布，一條貫通與底角，誘發(fā)底角剪切破壞;另一條則向底板中部演化，加劇底臌。

基于3.3節(jié)，研究開發(fā)的能量耗散模型除可定位圍巖主破壞位置，還可以通過控制運(yùn)算時(shí)步反演巷道變形破壞路徑。圖11給出了不同時(shí)步對(duì)應(yīng)的巷道圍巖耗散能密度演化云圖，運(yùn)算平衡時(shí)對(duì)應(yīng)的耗散能云圖如圖10所示。由圖11可知，巷道開挖后耗散能集中區(qū)位于巷道兩底角處，表明直角相比于拱形更易產(chǎn)生剪切破壞，這與前人研究結(jié)論相吻合。隨著時(shí)間增加，巷道底角破壞后，耗散能集中區(qū)域?qū)⑾蛏畈繃鷰r轉(zhuǎn)移，并逐漸演化成帶狀分布，一端向巷道兩幫轉(zhuǎn)移，一端朝向巷道底板方向，誘發(fā)巷道底臌和兩幫收斂。對(duì)于破壞后的巷道，如對(duì)底幫尤其是底角圍巖不加以控制，劇烈底臌將進(jìn)一步引發(fā)巷幫松動(dòng)破壞，甚至危及頂板穩(wěn)定，最終誘發(fā)巷道整體災(zāi)變失穩(wěn)。

圖11 巷道破壞路徑演化模擬Fig.11 Simulated results of the failure paths of roadway

對(duì)于試驗(yàn)巷道，后期采用了巷幫短錨索和底板錨桿(索)加固，巷道大變形得到有效控制。限于篇幅，具體參數(shù)不再贅述。綜上，研究開發(fā)的能量耗散模型對(duì)于反演巷道變形破壞路徑和定位巷道主破壞位置提供了一條新途徑，可為巷道失穩(wěn)預(yù)測(cè)預(yù)報(bào)和加固設(shè)計(jì)提供了參考，具有較強(qiáng)的實(shí)用價(jià)值。

5 結(jié) 論

(1)基于能量平衡和有限差分理論，推導(dǎo)了巖石耗散能有限差分方程式，采用FISH語(yǔ)言將其寫入FLAC3D應(yīng)變軟化模型，補(bǔ)充了軟件能量計(jì)算模塊。通過與室內(nèi)試驗(yàn)對(duì)比，該模型可有效描述巖石峰后變形特征及破壞路徑。

(2)利用開發(fā)的能量耗散模型分析了深部巷道能量演化規(guī)律，采用耗散能演化不僅可定位巷道主破壞位置，同時(shí)可有效反演圍巖全過程變形破壞路徑，為巷道變形預(yù)測(cè)及加固設(shè)計(jì)提供了一條新途徑。

(3)巖石能量耗散模型在FLAC3D中的實(shí)現(xiàn)，擴(kuò)展了FLAC3D軟件的適用范圍，一定程度上彌補(bǔ)了常規(guī)采用塑性區(qū)定性判斷地下工程圍巖破壞及其穩(wěn)定的不足，具有重要的理論意義和現(xiàn)場(chǎng)實(shí)用價(jià)值。