FLAC3D中巖石能量耗散模型的開發(fā)與應用

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

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

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

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

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

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

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

1.1 巖石單、三軸壓縮試驗

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

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

圖1 巖石應力-應變曲線Fig.1 Stress-strain curves of the rock samples

1.2 巖石能量耗散特征

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

Wz=We+Wd

(1)

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

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

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

(2)

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

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

(3)

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

Wd=Wz-We

(4)

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

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

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

σ1<σ2<σ3

(5)

相應的主應變增量Δεi[25]可分解為

(6)

單元體彈性應變增量和應力增量之間的關系[25]可表示為

(7)

(8)

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

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

(9)

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

(10)

式中，ψ為剪脹角。

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

(12)

(13)

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

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

(14)

式中，c為黏聚力。

拉伸破壞面上應力[25]滿足:

ft=σ3-σt=0

(15)

式中，σt為抗拉強度。

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

圖3 單元能量計算簡圖Fig.3 Schematic diagram for the energy calculation

(16)

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

(17)

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

(18)

不考慮峰后彈性模量衰減，單元可釋放彈性應變能可由下式計算:

(19)

則單元耗散能計算式如下:

(20)

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

3 模型校驗與運算

3.1 巖樣模型建立

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

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

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

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

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

圖6 應力-應變曲線校驗Fig.6 Verifications of the stress-strain curves

3.2 巖樣能量模型校驗

考慮到耗散能是由巖石總能量減去彈性應變能計算得到，限于篇幅，校驗能量模型時只校驗總能量和耗散能兩個指標。如圖7(a)所示，不同圍壓下巖樣總能量數(shù)值模擬結(jié)果與試驗結(jié)果基本吻合，表明FISH語言編寫的能量算法可靠性較強。

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

圖7 能量模型校驗Fig.7 Verifications of the energy model

3.3 巖樣峰后破壞路徑反演

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

圖8 巖石破壞模式校驗Fig.8 Verifications of the failure modes of rock samples

工程現(xiàn)場，準確掌握圍巖變形破壞路徑，精準定位加固位置是控制其穩(wěn)定的前提。為此，反演巖石峰后耗散能演化過程可為揭示其主破壞區(qū)演化、評判巖石穩(wěn)定性提供新途徑。圖9給出了圍壓25 MPa時不同峰后目標點巖樣耗散能密度演化云圖。巖樣初始進入峰后階段(B點)，經(jīng)歷裂隙孕育形成初始主裂隙面，并隨著軸向應變增加，主裂隙面耗散能逐漸增加，一旦巖樣進入殘余階段(F點)，主裂隙面耗散能密度增幅減小并趨于穩(wěn)定，表明巖樣此時已發(fā)生宏觀破斷，繼續(xù)加載將導致其結(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 能量耗散模型的應用

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

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

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

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

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

圖10 試驗巷道變形破壞對比Fig.10 Comparisons on the roadway deformations

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

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

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

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

5 結(jié) 論

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

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

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