孫琪皓 馬鳳山 趙海軍 郭 捷 馮雪磊

( ①中國科學院地質與地球物理研究所，中國科學院頁巖氣與地質工程重點實驗室 北京 100029)

( ②中國科學院地球科學研究院 北京 100029)

( ③中國科學院大學 北京 100049)

0 引 言

隨著人類對資源需求量的日益增加，能源開采逐漸向更深和更復雜的地層中進軍，深部地質工程所面臨的“三高和強擾動強時間效應”問題更加凸顯，巷道圍巖冒頂、幫縮、底鼓，分區(qū)化破裂，巖爆，突水等非線性物理力學現(xiàn)象更加突出( 何滿朝，2014) ，開挖后的巷道處于更加復雜的地質力學環(huán)境當中，其變形和破壞過程往往是在多個物理場耦合作用下發(fā)生的結果( 周宏偉等，2005) 。謝和平等( 2015) 詳細闡述了國內外關于深部與淺部的量化界定，何滿朝等( 2005) 對深部和淺步情況下巷道開挖引起的力學差異進行了詳盡分析，賀永年等( 2006) 分析了深部巖石工程研究中存在的若干問題，并揭示了巷道圍巖斷裂間隔分布的現(xiàn)象。周小平等( 2007) 研究了深部開采產生分區(qū)破裂化現(xiàn)象的機理，姜耀東等( 2004) 分析了深部巷道底鼓的不同類型和影響因素，并總結了深井軟巖巷道底鼓的防治措施。劉高( 2000) 論述了高應力軟巖的特點，總結了高應力軟巖巷道的變形類型及特征，分析了開挖前后地應力變化對巷道穩(wěn)定性的影響作用機理。李術才等( 2008) 通過對分區(qū)破裂化現(xiàn)象的現(xiàn)場監(jiān)測，分析了深部圍巖的破裂模式。還有眾多學者通過現(xiàn)場測量、物理模擬實驗、數(shù)值模擬等不同手段對巷道圍巖穩(wěn)定性機理進行了系統(tǒng)研究。

地下工程的安全問題常常與地下水的作用密不可分。自Biot( 1941) 進行三維固結問題的研究，奠定了流固耦合理論的研究基礎后，學者對滲流－應力耦合問題做了大量理論和試驗的研究工作，提出了一系列經驗公式和理論模型，總體上可以分為等效連續(xù)介質模型、離散裂隙網絡模型和雙重介質模型。近年來，巖體內部結構的損傷效應引起了廣泛的關注，損傷帶來了更加復雜的耦合效應，據(jù)此學者提出了許多基于斷裂力學和損傷力學的耦合理論。Shao et al. ( 2004) 建立了基于巖石微裂紋演化的各向異性損傷模型，并根據(jù)細觀力學分析提出了巖石宏觀力學參數(shù)及滲透張量與細觀損傷之間的關系。Oda et al. ( 2002) 進行試驗研究了巖石材料的相關力學參數(shù)與損傷變量，并考慮裂隙幾何張量分析了巖石滲透性的變化規(guī)律。榮傳新等( 2004) 根據(jù)彈塑性損傷理論、極值點失穩(wěn)理論分別推導出了滲流作用下巷道圍巖應力分布規(guī)律和巷道圍巖的穩(wěn)定性條件，并分析了突水事故發(fā)生的原因。盧應發(fā)等( 2007) 基于能量原理，通過Gibbs 自由能建立了一種正交損傷模型，并應用到流固耦合的研究中。鄭少河等( 2002) 基于自洽理論，研究了考慮損傷斷裂效應的巖體滲透張量演化過程，并建立了多裂隙巖體滲流損傷的理論模型。馮夏庭等( 2005) 通過自行研制的實驗裝置進行了應力－滲流－化學多場耦合的細觀力學試驗，對全過程顯微和監(jiān)測并對試驗現(xiàn)象進行了分析。楊天鴻等( 2001) 引入了滲透突跳系數(shù)這一概念，提出了巖石破裂損傷過程中的滲流－應力耦合機制。朱萬成等( 2009) 基于質量守恒和能量守恒定量，建立了考慮損傷的熱－流－力耦合模型，并分析了水力損傷的作用機制。沈振中等( 2009) 基于無單元法建立了巖體水力劈裂過程中的應力－滲流－損傷耦合模型，并對具有初始裂紋的平面應力模型進行了裂紋擴展的耦合分析。賈善坡等( 2009) 提出了與塑性損傷演化和滲流相關聯(lián)的Mohr-Coulomb 破壞準則，并采用該準則建立了考慮圍巖損傷和滲透性動態(tài)演化過程的流－固全耦合迭代計算模型。趙延林等( 2010) 利用斷裂力學和幾何損傷力學，基于滲透力對巖體的損傷和裂紋起裂的作用以及裂紋擴展導致滲流場改變的原理，建立了滲流－損傷－斷裂耦合模型，并對蓄水過程中庫岸變形機制進行了分析。王軍祥等( 2014) 給出了彈塑性應力－滲流－損傷耦合模型的數(shù)值解法，并根據(jù)智能反分析方法對不易測定的損傷參數(shù)進行了反演。卞康等( 2010) 利用Mazars 各向同性彈性損傷理論和線彈性斷裂力學理論，根據(jù)水工隧洞襯砌水壓致裂的實際過程建立了滲流－損傷－應力耦合模型。李利平等( 2011) 基于壓剪和拉剪雙重破壞準則，建立了可以適用于煤礦充填體的應力－滲流－損傷耦合方程，并采用有限元對斷層活化導致的突水機制進行了研究。翁其能等( 2013) 以代表體積單元( REV) 為研究對象，分析建立了長期處于高水壓狀態(tài)下隧道襯砌的滲流－應力－損傷耦合作用方程。冉小豐等( 2015) 考慮圍巖損傷演化特征與滲透性的動態(tài)變化特征，建立了滲流－應力－耦合模型，并對某油田鉆井過程中泥頁巖井壁的穩(wěn)定性進行了分析。劉黎等( 2016) 建立了采動煤巖體瓦斯?jié)B流－應力－損傷耦合模型，并設計程序對某煤礦瓦斯?jié)B流規(guī)律進行了模擬研究。劉仲秋等( 2012) 將圍巖和襯砌材料均視為彈塑性損傷材料，對錦屏二級水電站引水隧洞進行了從開挖支護到運行充水全過程的滲流－應力耦合模擬。目前，很多學者試圖通過分析巖石細觀缺陷結構的演化來解釋宏觀破裂的發(fā)展與耦合過程中現(xiàn)象，把宏觀和細觀參數(shù)通過一定的方式聯(lián)系起來，成為了當前解決工程問題和取得理論研究突破的焦點。

對于滲流－損傷－應力耦合問題已有的諸多研究，已經取得了一定的成果，但應用彈塑性損傷耦合模型進行巷道變形破壞的分析尚不多見，且現(xiàn)有模型仍需進一步改善以更好地反映巖石物理力學參數(shù)的動態(tài)變化特征尤其是峰后特征，以及巖體力學特性與水力學特性的不均勻性對巖體破壞的影響。本文基于前人的研究，充分考慮巷道開挖后破壞形態(tài)的影響因素，在巖體力學參數(shù)及水力學參數(shù)符合weibull 分布的假定下，考慮圍巖塑性屈服，滲透系數(shù)隨圍巖損傷演化以及圍巖力學參數(shù)弱化的過程，建立滲流－損傷－應力耦合作用數(shù)值模型，并利用多物理場數(shù)值模擬軟件，模擬分析巷道圍巖變形破壞的漸進性過程和關鍵部位的破壞演化過程，對比分析不同因素作用下巷道圍巖變形破壞模式的特征和區(qū)別，為合理開展地下工程和礦山工程提供可靠依據(jù)和科學指導。

1 滲流－損傷－應力耦合模型的建立

巷道開挖后，圍巖將發(fā)生應力重分布和應力回彈作用，當圍巖無法承受其作用時即發(fā)生塑性變形或破壞。前人已對巷道變形破壞類型及影響因素進行了系統(tǒng)地歸納和總結( 谷德振，1979; 張倬元，1981; 于學馥，1983; 王思敬，1984; 何滿朝等，2002; 孫廣忠，2004; 謝和平等，2015) 。

在應力和滲流的共同作用下，巷道圍巖即會處于損傷狀態(tài)，且損傷隨時間發(fā)生動態(tài)演化，同時損傷反作用于圍巖體，使其力學性質和水力學性質發(fā)生變化，在應力－損傷－滲流的耦合作用下，巷道圍巖產生漸進性的變形破壞。假設巖體介質為多孔彈性介質，下面給出耦合計算模型。

1.1 巖體非均質性的假定

處于地質環(huán)境中的巖體，其力學性質在空間上的分布是不均勻的，為了反映巷道圍巖的非均質性，引入weibull 統(tǒng)計分布函數(shù)來描述巖體的彈性模量、黏聚力、內摩擦角、滲透系數(shù)等物理力學參數(shù)，即:

式中，α 為巖體材料參數(shù)( 彈性模量等) ; α0為材料參數(shù)的均值; m 為均勻性系數(shù); φ( α) 為巖體力學性質α 的統(tǒng)計分布密度。

1.2 巖體力學參數(shù)弱化的表述

在應力和滲流環(huán)境中的巖體損傷實質上是巖體的力學參數(shù)不斷降低的過程( 江權，2008) 。地下工程巖體開挖后，由原先的三向應力狀態(tài)轉向兩向應力狀態(tài)，相當于完成了一次復雜的加、卸載過程，在此過程中產生應力重分布作用，應力集中或卸荷作用使得巖體產生不同程度的損傷，并導致一定范圍的屈服破壞，隨著損傷范圍和屈服范圍的動態(tài)擴大，應力將不斷發(fā)生重分布，進而又造成了更大范圍的損傷和屈服，并形成了損傷和應力重分布相互促進的動態(tài)過程。細觀上，在滲流、應力作用下圍巖內部原有的微小裂隙從閉合狀態(tài)趨于張開，并產生一定的擴展延伸，同時在拉應力和劈裂作用下圍巖可能促使新裂隙產生并發(fā)生裂隙間的貫通，使得圍巖性質產生明顯的劣化，從宏觀上可以認為巖體的彈性模量、黏聚力、內摩擦角等力學參數(shù)發(fā)生了變化( 黃勇等，2018; 趙晶等，2018) 。

從理論分析和現(xiàn)場實際工程經驗中，均可得知巖體的彈性模量隨著其受損程度的增加，而發(fā)生弱化不斷減小。根據(jù)典型的巖石單軸加載應力－應變曲線( 圖1) ，可以得出在彈性階段巖石的彈性模量基本保持不變或變化很小，進入屈服階段后，巖石的彈性模量開始下降，且其下降的速度由緩慢逐漸增加，在巖石達到峰值強度時變化速度不斷加快，直至巖石達到殘余強度時，巖石的彈性模量跌落到殘余值并基本不再發(fā)生變化。這一過程是加載進行時巖體的彈性模量隨著塑性應變的變化而發(fā)生的，而等效塑性應變可以用來較好的描述巖體屈服后的塑性程度，即:

圖1 巖石單軸加載應力－應變曲線圖Fig. 1 Rock stress-strain curve in uniaxial loading

因此，基于典型的巖石應力－應變曲線特征，考慮工程巖體彈性模量在加載過程中動態(tài)變化的特點，充分反映巖體力學參數(shù)的峰后變化特點以及應力跌落的突然性，簡化其弱化的過程如圖2，建立彈性模量關于巖體等效塑性應變的二次型弱化函數(shù)關系:

式中，E 為巖體的彈性模量; E0和Ed分別為巖體的初始彈性模量和到達殘余階段時的殘余彈性模量;為等效塑性應變;為彈性模量達到殘余值時的等效塑性應變。

圖2 彈性模量弱化過程曲線Fig. 2 Rock modulus of elasticity curve of weakening process

1.3 損傷變量的定義

為了定量描述巷道圍巖的損傷程度，定義損傷變量D?？紤]到巖體在達到屈服極限時，開始出現(xiàn)明顯的損傷現(xiàn)象并產生一定的塑性應變，因此本文將等效塑性應變作為損傷程度的判據(jù)，并通過上文彈性模量E 的弱化來定義損傷變量D，建立如下?lián)p傷演化方程:

考慮到巖體加卸荷過程中，當發(fā)生屈服作用時，微裂隙趨于張開，產生新裂隙并使裂隙間趨于貫通，從細觀角度上來看，巖體內部結構不斷分離，不連續(xù)面增多，巖體的完整性受到了破壞，表現(xiàn)為細觀參數(shù)黏聚力和內摩擦角的減小。據(jù)此建立簡化的巖體損傷過程中黏聚力c，內摩擦角φ 的動態(tài)損傷變化過程:

其中，C0和Cd分別為巖體的初始黏聚力和到達殘余階段時的殘余黏聚力; φ0和φd分別為巖體的初始內摩擦角和到達殘余階段時的殘余內摩擦角。

1.4 應力場控制方程

假定研究對象為理想彈塑性體，考慮流體的孔隙水壓力作用和損傷影響的本構關系如下:

式中，σij為應力分量; G( D) 為巖體的剪切模量，是損傷變量D 的函數(shù); λ( D) 為拉梅常數(shù)，其值大小為E( D) ν/［( 1 +ν) ( 1－2ν) ］; E( D) 為巖體的彈性模量，為損傷變量D 的函數(shù); ν 為泊松比( 變化忽略不計) ; εij為應變分量;為塑性應變分量; εkk為體積應變分量;為塑性體積應變分量; δij為Kronecker 系數(shù); α 為Boit 系數(shù); p 為孔隙水壓力。

式中彈性部分可根據(jù)胡克定律計算求解，塑性部分采用基于D－P 準則的理想彈塑性模型，其屈服條件為:

其中，αf、Kf為與巖體內摩擦角φ 和黏聚力c 有關的實驗常數(shù); J2為應力偏量第二不變量; I1為應力第一不變量。

1.5 滲流場方程

假設流體不可壓縮，無源穩(wěn)定情況下的滲流連續(xù)性方程為:

代入Darcy 定律方程，得到滲流場控制方程為:

通過室內進行的巖石全應力－應變過程滲透性試驗結果可知( 李世平，1995; 韓寶平，2000; 姜振泉，2001;李樹剛，2001) ，巖石在所受荷載逐漸增大的過程中，會發(fā)生細觀結構的損傷，其滲透率在此過程中會發(fā)生突跳性增長，因此，為表述這一現(xiàn)象引入滲透率突跳系數(shù)( 楊天鴻等，2001) ，同時考慮應力和損傷過程對滲透系數(shù)的影響，假定應力與滲透系數(shù)呈負指數(shù)關系，建立損傷演化過程中的滲流－應力耦合方程如下:

其中，k 和k0分別為滲透系數(shù)和滲透系數(shù)初值; P 為孔隙水壓力; ξ、α、β 分別為滲透系數(shù)突跳倍率、孔隙水壓力系數(shù)、耦合系數(shù)。

2 數(shù)值模擬結果分析

2.1 數(shù)值模型的建立

為了驗證耦合模型的合理性，建立巷道平面應變計算模型( 圖3) ，計算區(qū)域為寬( x 方向) 40 m，高( y 方向) 40 m 的正方形，巷道斷面采用半圓拱形，圓半徑為3 m，矩形區(qū)域寬6 m，高3 m。采用三角形單元網格剖分，共劃為6636 個域單元和228 個邊界單元。模型底部施加固定約束，初始地應力等效為邊界荷載施加于兩側和頂部，孔隙水壓力施加于模型頂部，巷道內邊界為自由界面且初始水頭為零。巷道圍巖物理力學參數(shù)及耦合參數(shù)如表1 所示。

圖3 數(shù)值模擬的計算網格模型Fig. 3 Calculate grid model of numerical simulation

2.2 力學參數(shù)弱化對巷道變形破壞影響分析

為了反映力學參數(shù)弱化對巷道圍巖損傷破壞情況的影響程度，采用以上計算參數(shù)，在相同的初始地應力和邊界條件下，分別對不考慮參數(shù)弱化和考慮參數(shù)弱化的兩種情況進行計算。不考慮滲流、應力對巖體損傷的影響，忽略加卸載過程中巖石物理力學參數(shù)的動態(tài)變化，并認為參數(shù)為常值時，計算結果如圖4 所示，當考慮力學參數(shù)按模型給出的方式弱化時，計算結果如圖5 所示。

模型計算結果顯示，在考慮力學參數(shù)弱化過程的滲流－應力耦合模型計算中，塑性屈服破壞區(qū)要明顯大于未考慮力學參數(shù)弱化時的情況，這主要體現(xiàn)了在加卸載過程中圍巖體處于損傷演化的過程，不連續(xù)面的延伸、貫通和產生對于巖體整體質量的損傷作用無法忽略，反映出巖體損傷動態(tài)變化的過程對于巷道圍巖的變形破壞的程度可產生較大的影響，表明了在滲流－應力耦合計算過程中引入力學參數(shù)弱化的必要性。

2.3 不同深度下巷道變形破壞分析

地下工程圍巖體的變形破壞常常是漸進性發(fā)展的，由于巖體的非均質性、幾何結構特點、應力分布不均勻等原因，應力集中程度最高的部位往往成為巷道圍巖變形破壞的突破口，在局部產生破壞后應力重分布作用將應力轉移至其他部位，再次引起其他部位的逐一破壞，最終導致大范圍的失穩(wěn)破壞。長期以來，人們對自重應力引起的變形破壞進行了大量的研究，但隨著地下工程的深度不斷增加，尤其是受到褶皺和擠壓性斷層等構造因素影響的區(qū)域，水平地應力可能會成為影響巷道穩(wěn)定性的主導因素?？紤]到不同深度地應力下巷道圍巖變形破壞過程的差異性，根據(jù)側壓力系數(shù)的取值，分別設置3 種情形進行計算分析。情形一: 水平地應力等于垂直地應力，側壓力系數(shù)λ =1; 情形二: 水平地應力大于垂直地應力，側壓力系數(shù)λ = 2; 情形三: 水平地應力小于垂直地應力，側壓力系數(shù)λ =1/2。計算組合詳細的地應力取值見表2。

表1 數(shù)值模擬的計算參數(shù)Table 1 Physical and mechanical parameters for numerical simulation

圖4 不考慮參數(shù)弱化的塑性屈服云圖Fig. 4 Plastic yield contours without considering parameter weakening

圖5 考慮參數(shù)弱化塑性屈服云圖Fig. 5 Plastic yield contours considering parameter weakening

表2 不同側壓力系數(shù)下的地應力Table 2 Initial stress with different lateral pressure coefficients

當側壓力系數(shù)等于1時，巷道圍巖塑性屈服區(qū)的變化過程如圖6 所示。從圖中可以看出，巷道底角處最先產生屈服破壞，隨著時間的推移，頂拱附近開始產生明顯的環(huán)形損傷屈服區(qū)，并逐漸向邊墻延伸最后塑性區(qū)發(fā)展到整個巷道周邊，呈現(xiàn)出四周近乎均勻的屈服破壞狀態(tài)。這表明在水平地應力與垂直地應力趨于相等的情況下，底角處因其棱角狀容易形成應力集中區(qū)，切向應力過大導致該處首先產生屈服，繼而頂板由于在重力作用下所受拉應力相對較大，從頂板開始不斷向兩側邊墻發(fā)生屈服，因為水平和垂直地應力均發(fā)揮主導作用，所以巷道圍巖最后的狀態(tài)大致為沿巷道周邊形成了近于均勻的屈服狀態(tài)。

當側壓力系數(shù)為2 時，巷道圍巖塑性屈服區(qū)的變化過程如圖7 所示。從圖中可以看出，巷道頂拱處最先產生塑性屈服，隨著時間的推移，水平應力向深部傳遞的趨勢十分明顯，頂拱和底板處產生均出現(xiàn)較大程度的屈服并不斷發(fā)展擴大，最后在頂拱和地板處形成了不規(guī)則梯形狀的屈服區(qū)。這表明了在這種地應力水平下，水平應力過高導致頂板處產生較大的剪切應力，剪切應力集中使得頂拱處成為應力集中系數(shù)最高的部位，因此首先產生屈服破壞，繼而底部由于同樣承受較大的剪切應力，故屈服破壞主要沿頂拱、底板處延伸發(fā)展，而兩側邊墻的變化程度很小，可見在水平地應力主導的情況下，頂、底部是最容易發(fā)生破壞的部位，頂部容易發(fā)生剪切破壞、楔形冒落等現(xiàn)象，底部容易產生底鼓、褶皺體破壞等現(xiàn)象。

圖6 側壓力系數(shù)為1 時的破壞過程Fig. 6 Damage process under the lateral pressure coefficients of 1

圖7 側壓力系數(shù)為2 時的破壞過程Fig. 7 Damage process under the lateral pressure coefficients of 2

圖8 側壓力系數(shù)為1/2 時的破壞過程Fig. 8 Damage process under the lateral pressure coefficients of 1/2

當側壓力系數(shù)為1/2 時，巷道圍巖塑性屈服區(qū)的變化過程如圖8 所示。從圖中可以看出，在巷道底角處因應力集中首先產生屈服之后，巷道兩側的邊墻拐角部位也開始產生損傷破壞，并沿兩側邊墻不斷擴展，最后覆蓋整個邊墻區(qū)域，在兩側邊墻處形成了不規(guī)則條帶狀的屈服區(qū)，而在頂、底板處變形相對較小。這表明了在垂直地應力占主導作用的情況下，巷道的破壞變形過程有著明顯的區(qū)別，在圍巖兩側更容易產生剪切應力集中，因此兩側邊墻是最可能產生破壞的部位，有可能發(fā)生表層剝落、脆性錯斷甚至在頂角處產生大塊楔形體坍落、巖爆的現(xiàn)象。

由損傷屈服面積的變化過程曲線( 圖9) 可以看出，在地應力水平不同情況下的加載過程，損傷區(qū)面積均呈現(xiàn)不斷增加的趨勢，且變化速度均為不斷加快。這表明隨著加載過程的進行，巖體損傷破裂的程度不斷加劇，耦合作用效果逐漸加強，在多物理場耦合的作用下，巖體的力學參數(shù)不斷弱化，巖體質量逐漸降低，使巖體產生屈服破壞形成損傷區(qū)的速度變得更快。

3 結 論

本文首先總結了前人關于巖體損傷與滲流－應力耦合的研究。在充分考慮了巖體的非均質性、塑性階段的屈服過程以及巖體力學參數(shù)隨加卸載狀態(tài)發(fā)生動態(tài)弱化的基礎上，建立了滲流－損傷－應力耦合的數(shù)值計算模型。將建立的滲流－損傷－應力耦合模型應用于巷道圍巖變形破壞的模擬，主要結論如下:

(1) 考慮力學參數(shù)弱化的耦合模型與其他模型相比，得出了更加嚴重的圍巖變形破壞結果，充分反映了巷道開挖后應力重分布對于圍巖造成的不可逆轉的損傷作用，損傷必然造成力學參數(shù)的弱化，這種損傷是必然存在且不可忽視的，每一時步動態(tài)損傷所造成的結果必然會對下一時步的加卸載過程造成影響，產生疊加損傷效應，這將對巷道圍巖的穩(wěn)定性產生更加不利的影響，而這在以往的模型中并沒有充分體現(xiàn)。由于材料損傷對于工程問題的重要性，應從損傷狀態(tài)變化的角度上考慮滲流－損傷－應力耦合問題，才能提出更能保證安全和符合實際的解答。

圖9 側壓力系數(shù)為1、2、1/2 時損傷面積的變化過程Fig. 9 Change of damage area under the lateral pressure coefficients of 1、2、1/2

(2) 應用本文提出的模型對處于不同深度和地應力水平下的半圓拱形巷道穩(wěn)定性進行分析，通過其屈服破壞的漸進過程，得出在滲流－損傷－應力耦合作用下，產生破壞失穩(wěn)的過程是由緩慢到逐漸加快的，多物理場耦合作用的效果是逐漸加強的，當側壓力系數(shù)趨向于1 時，巷道周邊最終均會發(fā)生相近程度的屈服破壞，當側壓力系數(shù)大于1 時，巷道頂拱和地板最終會出現(xiàn)似梯形狀的大面積破壞，當側壓力系數(shù)小于1 時，屈服破壞區(qū)則主要分布在巷道兩側邊墻處，呈不規(guī)則條帶狀。結論對于地下工程進行開挖后，能正確掌握并控制不同地應力下巷道發(fā)生破壞的關鍵部位，并有針對性地采取經濟合理的支護措施來加固圍巖有著重要意義。