地應力環境對深埋硐室巖爆特性影響的數值模擬研究

趙紅亮，仇 巖，梁海安，楊澤平，張 娟

(東華理工大學 土木與建筑工程學院，江西 南昌 344000)

隨著礦山建設、水利水電和油氣開采等行業的蓬勃發展，我國地層的淺部資源正在日益枯竭，深部開采逐漸成為資源開發的常態。當巖體開挖深度超過700 m，工程開挖作業即進入深部地質環境，在深部地質環境下開挖的硐室被稱為“深埋硐室”。與淺埋硐室相比，深埋硐室具有埋置深度大、原巖應力和地應力高等顯著特征[1]。地應力賦存于巖體之中，是地下工程環境中最主要的指標之一，作為高地應力環境下特有的工程地質災害現象，巖爆是硐室邊緣部位應力集中超過了巖體的強度極限，從而導致巖體發生脆性破壞，并伴有圍巖片幫剝落、彈射拋擲等動力失穩的現象[2]，在煤炭行業被稱為沖擊地壓。輕微巖爆會使硐室圍巖產生片幫現象，而強烈的巖爆往往會使巖石爆裂，崩出的破碎巖石會造成設備損壞、工期延誤、人員傷亡等嚴重后果，巖爆問題已成為巖體地下工程和巖石力學領域的世界性難題。巖爆的產生與地應力聚集特性和地應力量級有著密切關系，一般來說，較高地應力區的巖體相較于較低地應力區的巖體，彈性模量更大，強度更高，儲存的彈性變形能更多，更易產生巖爆。

王學濱[3]等利用FLAC模擬了不同圍壓條件下圓形巷道的巖爆過程，模擬結果表明，隨著圍壓的增大，圓形巷道“V型”巖爆坑的破壞深度增加、圓形硐室附近剪切帶花樣的對稱性變差，在高壓狀態下，剪切帶花樣與滑移線網有相似之處；肖繁[4]等基于最大主應力的巖爆判別準則，預測了巷道周邊各部位的巖爆烈度等級，巷道不同部位的巖爆烈度等級由高到低依次為：底板兩側＞巷道兩幫＞拱頂邊緣；劉劍[5]等采用主成分分析法和改進貝葉斯(Bayes)判別法，選取巖石單軸抗壓強度、脆性系數、巖石彈性能量等指標，建立了巖爆綜合預測模型；陳衛忠[6]等按照地下工程開挖卸荷特點，開展了脆性花崗巖常規三軸、不同卸載速率條件下峰前、峰后三軸卸圍壓試驗，探討了巖爆巖石的變形破壞特征和巖爆形成力學機制；陳結[7]等對長期預測和短期預警方法及存在的問題、機器學習方法、沖擊地壓短期預警及長期預測的應用，以及傳統的經驗驅動和機理驅動與機器學習的數據驅動在沖擊地壓預測預警方面的聯系與差異做了總結；羅天敏[8]等對國家能源集團新疆能源有限責任公司下屬礦井沖擊地壓監測預警現狀進行了分析，提出了多參量集成監測預警平臺本地化建設方案；楊凡杰[9]等為了對巖爆進行更準確的預測，結合試驗和理論分析認為，巖爆機制的研究在本質上應是對巖體中裂紋動態擴展問題的研究，進而提出一個新的巖爆能量判別指標——單位時間相對能量釋放率指標；張修峰[10]等通過沖擊地壓多參量聯合監測預警方法研究、監測監控平臺系統開發和現場實踐，探索了新形勢下沖擊地壓監控預警新模式，實現了礦井與防沖研究中心、井上與井下、靜態數據與動態數據、預警結果與現場處置措施等智能聯動及自動上報；呂鵬飛[11]等提出一種基于粒子群算法(PSO)優化最小二乘支持向量機(LSSVM)預測方法，即沖擊地壓分級預測的PSO-LSSVM方法，該方法綜合考慮煤礦開采深度、地質構造、煤的堅固性系數、最大主應力、煤層傾角變化、煤厚變化、頂板巖層厚度、開采工藝、頂板和底板巖石強度共10項指標因素，構建沖擊地壓預測指標體系；陳衛忠[12]等開展脆性花崗巖常規三軸、不同控制方式、不同卸載速率條件下，峰前、峰后卸圍壓試驗，從能量的原理探討巖石破壞過程能量積聚-釋放的全過程，研究巖石的變形破壞特征、能量集聚-耗散-釋放特征和基于能量原理的巖爆判據；姚精明[13]等采用物理實驗和多重分形理論相結合的方法，研究了煤樣變形破壞過程中多重分形譜寬度與其沖擊傾向、釋放能量的內在關系；呂森鵬[14]對實際工程中巖爆的發生機理、預測預報模型開展試驗與理論研究、對巖爆防治對策進行實際探討，從而提出了可靠的預測方法和防治對策；陳玉濤[15]提出一種鄰近地質構造的不規則孤島工作面沖擊危險性評價方法；夏永學[16]等針對綜合指數法和多因素耦合法在進行沖擊危險性評價時存在臨界突變和權重量化困難等問題，通過因素分類、指數疊加和歸一化處理，建立了改進的綜合指數方法；宮鳳強[17]等基于現有的巖爆分級判據分析高地應力環境下圍巖的破壞現象，提出了高地應力“強度&應力”耦合判據及其定量分級標準，結果表明，基于巖石強度和地應力兩種要素提出的“強度&應力”耦合判據能夠在眾多工程實例檢驗結果中證實該判據的合理性；榮海[18]等對堅硬巖層結構特征、結構形成條件、結構失穩尺度、結構失穩釋放能量等進行了分析，并分別給出了相應的計算公式，確定堅硬巖層失穩對沖擊地壓的控制作用。

綜上所述，現有的研究大都從影響巖爆產生的單一因素作為出發點，引入巖爆判別指標，判斷地應力對巖爆特征的影響，但是影響巖爆產生的因素具有多樣性，僅憑單一的巖爆判別指標研究地應力對巖爆特征的影響過于片面，因此筆者在前人研究的基礎上，利用硬巖在某個特定圍壓狀態下的試驗數據，擬合出了硬巖的極限儲能值經驗公式，從導致巖爆發生的能量角度和應力角度出發，引入了彈性變形能和巖石應力強度比作為巖爆發生的判別指標，通過FLAC3D數值模擬技術輸出巖爆判別指標云圖和監測曲線，結合相關巖爆傾向判據分析地應力環境對巖爆烈度等級和破壞區域的影響。

1 巖爆的評價指標和判據

1.1 彈性變形能指標

基于彈性變形能理論，巖體產生破壞主要是能量驅動的結果。當巖體可釋放的彈性變形能Ue達到巖體破壞所需的極限儲能U0時，也即Ue=U0時，巖體內部彈性變形能完全釋放，巖體產生靜態破壞，此時不會發生巖爆災害；當Ue=U0時，巖體發生動態破壞，誘發巖爆的產生。

彈性變形能作為表征巖爆現象發生的物理量，能夠結合相關巖爆判別準則來對巖爆發生的烈度等級進行分析。

巖體發生巖爆破壞的條件為

針對巖石的極限儲能值U0，筆者結合文獻[6，12，14]中有關花崗巖在加、卸載速率為0.1 MPa/s，不同控制方式下的極限儲能結果(表1)，擬合得到了其在特定圍壓狀態下的極限儲能經驗公式，擬合曲線如圖1所示。

圖1 巖石極限儲能與圍壓擬合曲線Fig.1 Fitting curve of rock limit energy storage and confining pressure

表1 花崗巖在不同控制方式下的極限儲能Table 1 Limit energy storage of granite under different control methods

極限儲能經驗公式為

彈性變形能物理表達式為

式中，σ1,σ2,σ3為單元體的3個主應力；μ為泊松比；E為彈性模量。

1.2 能量判別準則

陳衛忠[12]等定義巖爆能量指數為巖體單元的真實能量U與巖石極限儲能U0的比值，具體判據為

1.3 巖石應力強度比指標

HOEK[20]根據南非礦井的巷道破壞情況，利用巖石單軸抗壓強度和硐壁切向應力建立了表征巖爆與硐室圍巖應力關系的巖石應力強度比指標，即

式中，σθ為硐壁切向應力；σc為巖石單軸抗壓強度。

E.Hoek判別準則為

2 數值模擬結果分析

2.1 數值模型及巖石力學參數

硐室的長度遠大于其截面尺寸時，可視為平面應變問題，在FLAC3D程序中建立準三維數值計算模型，如圖2所示。

圖2 硐室三維數值計算模型Fig.2 Three-dimensional numerical model of chambers

模擬工況為：開挖隧硐埋深為1 500 m，模型尺寸為50 m×50 m，中央位置處為直徑10 m的圓形待開挖隧硐，左右兩側邊界施加水平方向位移約束，前后邊界施加水平方向位移約束，底部邊界施加垂直和水平方向位移約束，上部邊界施加上覆巖體自重應力。先對模型施加初始地應力場，當計算達到應力平衡狀態后，再開挖圓形硐室。假設巖體為均質各向同性，并采用FLAC3D內置的Morh-Coulomb本構模型。數值計算采用的花崗巖物理力學參數為：彈性模量E為50 GPa，泊松比為0.2，容重為26.9 kN/m3，黏聚力c為2.6，內摩擦角為43°。

2.2 不同地應力環境下主應力差分布特征

巖爆是高應力集中的結果，因此應著重關注硐室周圍的應力集中區域。巖爆產生的破壞類型多為壓剪破壞，導致巖體破壞的剪切應力可由主應力差間接表示。由于FLAC3D數值分析軟件無法直接輸出模型的主應力差，筆者基于主應力差理論，利用FISH語言進行軟件的二次開發，從而實現主應力差云圖的輸出。

圖3為不同側向壓力系數時的圓形硐室主應力差云圖?？紤]到實際工程當中面臨的地應力環境不同，筆者通過側向壓力系數λ為0.50，0.75，1.50，2.00時的模擬結果，分析地應力最大主應力方向對硐室巖爆產生的影響。

圖3 圓形硐室主應力差分布云圖Fig.3 Cloud map of main stress difference distribution in circular chambers

由圖3可知，受開挖擾動的影響，圍巖應力在一定范圍內重新分布，越靠近圓形硐室周邊區域，圍巖主應力差的分布不均勻性越顯著，隨著開挖遠離圓形硐室邊緣，主應力差分布不均勻程度明顯降低。

當地應力場中的水平應力小于豎向應力(λ＜1.0)時，硐室兩側的應力集中情況較為顯著，主應力差最大值位于硐室兩側，對應硐室兩側的剪應力集中程度較為明顯。應力集中區距離硐室邊緣約1.2 m，最大主應力差為68.762～73.772 MPa，該狀態下可能發生巖爆災害的區域主要為硐室兩側。

當地應力場中的水平應力大于豎向應力(λ＞1.0)時，硐室頂部和底部的應力由應力釋放狀態轉變為應力集中狀態，對應硐室頂部和底部區域的剪應力集中程度較為顯著，應力集中區距離硐室邊緣約3.1 m，最大主應力差為102.79～141.74 MPa，該狀態下可能發生巖爆的區域主要為硐室的頂部和底部位置。

2.3 不同地應力環境對硐室巖爆的發生及烈度等級的影響

2.3.1 彈性變形能變化特征及分析

基于彈性變形能理論得到巖體單元彈性變形能量公式(3)，利用FLAC3D內嵌FISH語言編寫計算函數，輸出圓形硐室彈性變形能云圖，通過巖體開挖后硐室周圍能量的變化，分析側向壓力系數對巖爆的影響。

圖4為不同側向壓力系數下，圓形硐室的彈性變形能分布云圖，當λ＜1.00時，彈性變形能積聚區位于硐室兩側，彈性變形能的釋放區域位于硐室的頂部和底部；當λ=0.50 時，彈性變形能最大值為1.939 5 MJ；當λ=0.75 時，彈 性 變形 能最 大 值為1.838 4 MJ。當λ＞1.00時，彈性變形能的積聚區域由硐室的兩側轉變為硐室的頂部和底部；當λ=1.5時，彈性變形能最大值為4.766 9 MJ；當λ=2.00時，彈性變形能最大值為5.326 MJ。

圖4 圓形硐室彈性變形能分布云圖Fig.4 Cloud map of elastic deformation energy distribution in circular chambers

由彈性變形能分布特征可知，地應力環境的變化會影響硐室周邊彈性變形能積聚區的分布及數值的變化，相應地發生巖爆破壞的部位和烈度等級也會發生變化。當λ＜1.00時，彈性變形能的積聚區域位于硐室兩側，對應硐室兩側發生巖爆的可能性較硐室其他部位大；當λ＞1.00時，彈性變形能積聚區域由硐室兩側位置轉移至硐室的頂部和底部，對應硐室頂部和底部兩處位置發生巖爆的可能較硐室其他部位大。

為了進一步探究側向壓力系數λ對硐室巖爆特征的影響，在圓形硐室周圍布置監測點，監測巖體開挖后彈性變形能的變化，得到硐室周圍彈性變形能變化曲線，結合巖爆能量判據對圓形硐室巖爆發生的破壞區域和烈度等級進行分析。由于模擬巖體為均質各向同性，開挖圓形硐室上下左右均對稱，因此在圓形硐室頂部、腰部分別布置一排5個監測點，如圖5所示。

圖5 圓形硐室監測點布置示意Fig.5 Layout diagram of monitoring points in circular chamber s

由圖6(a)和圖4(a)可知，當λ=0.50時，硐室頂部和底部的能量整體處于彈性變形能釋放狀態，巖體內部不存在彈性變形能積聚區，且各監測點的能量均未超過巖石的儲能極限值，因此硐室頂部和底部不會發生巖爆災害。

圖6 λ=0.50 時硐室周圍監測點彈性變形能變化曲線Fig.6 Variation trend of the elastic deformation energy curves of the monitoring point around the chambers when λ=0.50

由圖6(b)可知，6號監測點巖體單元的能量U6=0.68 MJ，巖體單元能量超過巖石的極限儲能值U06=0.42 MJ，巖爆能量指數U6/U06=1.6，因此硐室腰部6號監測點巖體單元存在巖爆風險；7號監測點巖體單元的能量U7=0.72 MJ，巖體單元能量超過巖石的極限儲能值U07=0.50 MJ，巖爆能量指數U7/U07=1.4，因此硐室腰部7號監測點巖體單元存在巖爆風險；8號監測點巖體單元的能量U8=0.63 MJ，巖體單元能量超過巖石的極限儲能值U08=0.62 MJ，巖爆能量指數U8/U08=1.0，因此硐室腰部8號監測點巖體單元存在巖爆風險。結合巖爆能量判據可知，硐室兩側表現為IV級強烈巖爆傾向。應力的調整使得能量積聚區向巖體內深部轉移，由巖體內部9號和10號監測點彈性變形能增量的變化，結合圖4(a)可知，圓形硐室腰部兩側的彈性變形能積聚區距離 硐室邊緣約3.8 m，影響范圍約2.5 m。

圖7為λ=0.75時硐室周圍監測點的彈性變形能變化曲線。

圖7 λ=0.75 時硐室周圍監測點的彈性變形能變化曲線Fig.7 Variation trend of the elastic deformation energy curves of the monitoring point around the chambers when λ=0.75

由圖7(a)和圖4(b)可知，當λ=0.75時，圓形硐室頂部和底部的能量整體處于彈性變形能釋放狀態，巖體內部不存在彈性變形能量積聚區，且各監測點的能量均不超過巖石的儲能極限值，因此硐室頂部和底部不會發生巖爆災害。

由圖7(b)可知，6號監測點巖體單元的能量U6=0.94 MJ，巖體單元能量超過巖石的極限儲能值U06=0.90 MJ，巖爆能量指數U6/U06=1.0，因此硐室腰部6號監測點巖體單元存在巖爆風險；7號監測點巖體單元的能量U7=1.03 MJ，巖體單元能量超過巖石的極限儲能值U07=0.95 MJ，巖爆能量指數U7/U07=1.1，因此硐室腰部7號監測點巖體單元存在巖爆風險；8號監測點巖體單元的能量U8=1.03 MJ，巖體單元能量超過巖石的極限儲能值U08=0.98 MJ，巖爆能量指數U8/U08=1.1，因此硐室腰部8號監測點巖體單元存在巖爆風險。結合巖爆能量判據可知，硐室兩側表現為IV級強烈巖爆傾向。應力的調整使得能量積聚區向巖體內部轉移，由巖體內部9號和10號監測點彈性變形能增量的變化，結合圖4(b)可知，圓形硐室腰部兩側的彈性變形能量積聚區域距離硐室邊緣約3.4 m，影響范圍約2.0 m。

圖8為λ=1.50時硐室監測點的彈性變形能曲線。

圖8 λ=1.50 時硐室周圍監測點彈性變形能變化曲線Fig.8 Variation trend of the elastic deformation energy curves of the monitoring point around the chambers when λ=1.50

由圖8(b)和圖4(c)可知，當λ=1.50時，圓形硐室腰部兩側的能量整體處于彈性變形能釋放狀態，巖體內部不存在彈性變形能量積聚區，因此腰部兩側不會發生巖爆災害。

由圖8(a)可知，1號監測點巖體單元的能量U1=2.20 MJ，巖體單元能量超過巖石的極限儲能值U01=1.24 MJ，巖爆能量指數U1/U01=1.8，因此硐室頂部1號監測點巖體單元存在巖爆風險；2號監測點巖體單元的能量U2=2.35 MJ，巖體單元能量超過巖石的極限儲能值U02=1.35 MJ，巖爆能量指數U2/U02=1.7，因此硐室頂部2號監測點巖體單元存在巖爆風險；3號監測點巖體單元的能量U3=2.40 MJ，巖體單元能量超過巖石的極限儲能值U03=1.40 MJ，巖爆能量指數U3/U03=1.7，因此硐室頂部3號監測點巖體單元存在巖爆風險；結合巖爆能量判據可知，硐室頂部和底部表現為IV級強烈巖爆傾向。應力的調整使得能量積聚區向巖體內部轉移，由巖體內部4號和5號監測點彈性變形能增量的變化，結合圖4(c)可知，圓形硐室頂部和底部的彈性變形能量積聚區域距離硐室邊緣約4.1 m，影響范圍約4.6 m。

圖9為λ=2.00時硐室監測點的彈性變形能曲線。

圖9 λ=2.00 時硐室周圍監測點彈性變形能變化曲線Fig.9 Variation trend of the elastic deformation energy curves of the monitoring point around the chambers when λ=2.00

由圖9(b)和圖4(d)可知，當λ=2.00時，圓形硐室腰部兩側的能量整體處于彈性變形能釋放狀態，巖體內部不存在彈性變形能量積聚區，因此腰部兩側不會發生巖爆災害。

由圖9(a)可知，1號監測點巖體單元的能量U1=2.50 MJ，巖體單元能量超過巖石的極限儲能值U01=0.98 MJ，巖爆能量指數U1/U01=2.6，因此硐室頂部1號監測點巖體單元存在巖爆風險；2號監測點巖體單元的能量U2=2.60 MJ，巖體單元能量超過巖石的極限儲能值U02=1.10 MJ，巖爆能量指數U2/U02=2.4，因此硐室頂部2號監測點巖體單元存在巖爆風險；3號監測點巖體單元的能量U3=2.70 MJ，巖體單元能量超過巖石的極限儲能值U03=1.30 MJ，巖爆能量指數U3/U03=2.1，因此硐室頂部3號監測點巖體單元存在巖爆風險；結合巖爆能量判據可知，硐室頂部和底部表現為IV級強烈巖爆傾向。應力的調整使得能量積聚區向巖體內部轉移，由巖體內部4號和5號監測點彈性變形能增量的變化，結合圖4(d)可知，圓形硐室頂部和底部的彈性變形能量積聚區域距離硐室邊大約4.5 m，影響范圍約5.8 m。

2.3.2 巖石應力強度比變化特征及分析

利用巖石應力強度比來表征巖爆與硐室圍巖應力的關系，基于FISH語言編寫巖石應力強度比公式(4)的計算函數，通過提取應力集中區監測點的最大切向應力，替代彈性分析中硐壁最大切向應力的解析解，輸出圓形硐室巖石應力強度比云圖，從力學角度研究不同側向壓力系數λ對硐室巖爆特征的影響。

圖10為不同側向壓力系數λ下圓形硐室巖石應力強度比云圖。

圖10 不同側向壓力系數λ 下圓形硐室巖石應力強度比云圖Fig.10 Rock strength stress ratio nephogram of circular chambers under different lateral pressure coefficients

由圖10可知，當λ＜1.00時，巖石應力強度比的最大值位于硐室兩側；當λ=0.50時，巖石應力強度比最大值為1.139 2；當λ=0.75時，巖石應力強度比最大值為0.997 46，結合E.HOEK巖爆判據可知硐室兩側部位表現為IV級強烈巖爆傾向。當λ＞1.00時，巖石應力強度比的最大值位于硐室頂部和底部；當λ=1.50時，巖石應力強度比最大值為1.505 8；當λ=2.00時，巖石應力強度比最大值為2.060 1，結合E.HOEK巖爆判據可知硐室頂部、底部表現為IV級強烈巖爆傾向。

3 結 論

(1) 當λ＜1.00時，應力集中區域和彈性變形能聚集區位于硐室兩側，隨著側向壓力系數λ的增加，應力集中區的主應力差最大值降低，能量聚集區的彈性變形能數值逐漸下降。應力重新分布后，圍巖內部彈性變形能聚集區向巖體內部轉移的距離由3.8 m縮小至3.4 m，對應巖爆產生的風險變小，破壞深度減小。

(2) 當λ＞1.00時，應力集中區域和彈性變形能聚集區域位于硐室頂部和底部，隨著側向壓力系數λ的增加，應力集中區的主應力差最大值增大，能量聚集區的彈性變形能數值逐漸增加。應力重新分布后，圍巖內部彈性變形能聚集區向巖體內部轉移的距離由4.1 m增加至4.5 m，對應巖爆產生的風險變大，破壞深度增大。

(3) 不側向壓力系數λ下，硐室存在巖爆災害的風險等級均表現為IV級強烈巖爆。當λ＜1.00時，產生巖爆的部位在硐室兩側，隨著側向壓力系數λ的增大，該部位的巖石應力強度比逐漸減小，對應巖爆發生的風險減小，破壞范圍縮??；當λ＞1.00時，產生巖爆的部位在硐室頂部和底部，隨著側向壓力系數λ的增大，巖石應力強度比值由1.5增大至2.1，對應巖爆產生的風險變大，破壞范圍增大。

筆者借鑒國內相關學者對巖石儲能機制的相關試驗研究，重新擬合了不同圍壓下花崗巖的極限儲能經驗公式，引用相應的巖爆指數能量判別準則用于表征巖石受載過程中的儲能特征和能量轉化，可以較好地從能量角度解釋巖石工程失穩破壞機制，但限于巖石力學試驗的復雜性，對加卸載速率的影響考慮還不夠充分，尚有不足之處。