硬脆性巖石多尺度損傷蠕變模型及長期強度研究

趙倫洋，賴遠明，牛富俊，李鵬飛，朱其志

(1.華南理工大學亞熱帶建筑科學國家重點實驗室，廣東廣州，510641；2.華南理工大學華南巖土工程研究院，廣東廣州，511442；3.中國科學院西北生態環境資源研究院凍土工程國家重點實驗室，甘肅蘭州，730000；4.巴黎東部大學多尺度模擬與仿真實驗室，法國巴黎，77454；5.河海大學巖土力學與堤壩工程教育部重點實驗室，江蘇南京，210098)

隨著我國對地下空間開發利用和對地下資源需求的日益增大，大量的深部重大地下工程已陸續啟動實施，如地熱開發[1]、高放廢物地質處置庫[2]、川藏鐵路[3]等。大量研究表明[4?7]，深部高地應力硬脆性工程圍巖在承受小于其抗壓強度的荷載作用下，有可能產生較大的蠕變變形或誘發圍巖遲滯性巖爆災害，對工程的穩定性造成顯著的不利影響。由于硬脆性巖石材料的復雜性(非均質性、多相性、工程時效性、體積膨脹現象等)及工程環境的不確定性，建立硬脆性巖石蠕變本構模型可以為深部巖體工程的長期穩定性分析提供重要的理論依據。現有的巖石蠕變本構模型主要有元件模型和經驗模型兩大類。元件模型采用彈簧、滑塊和黏壺分別描述巖石的彈性、塑性和黏性行為，并通過它們之間的串并聯來模擬巖石的應力、應變與時間的關系，其中，西原模型、Bingham模型和Burgers 模型得到了較為廣泛的應用。這類模型雖然概念明確，表達直觀，但由于其參數確定和模型識別的復雜性以及對于加速蠕變階段描述的局限性，限制了其在工程上的應用。經驗模型是基于不同試驗條件下的試驗結果和現場觀察來確定應力、應變及時間函數關系的一種唯象模型，包括老化理論、時空理論、流動理論和績效理論等。由于對試驗的依賴性較大，經驗模型無法適應環境因素的復雜變化。此外，元件模型和經驗模型均是基于巖石受力變形的宏觀表征建立，屬于宏觀流變學范疇。目前，從硬脆性巖石蠕變細觀機理出發，建立能夠描述硬脆性巖石蠕變過程中非線性特征的蠕變模型成為巖石流變力學研究的重點和難點。

長期強度作為巖石流變力學范疇中又一重要的科學問題一直備受關注。巖石長期強度的定義為其強度隨時間而持續降低，并逐漸趨近于一個穩定收斂的低限定值[8]。目前確定巖石長期強度的方法主要有3 種：1)通過繪制應力?應變等時曲線簇，由曲線簇上的拐點確定長期強度[9]；2)根據穩態蠕變速率和應力水平的關系曲線，將曲線過渡到直線段的拐點確定為長期強度[10]；3)將體積應變開始反向時對應的裂紋損傷應力作為巖石的長期強度[11?12]。對于硬脆性巖石而言，裂紋損傷應力確定長期強度方法受到國內外學者的認同，MARTIN[11?12]通過對加拿大Lac du Bonnet花崗巖短期和長期強度進行系統研究，提出了以試樣開始擴容時裂紋損傷應力作為該試樣的長期強度；SZCZEPANIK 等[13]根據該理論得出花崗巖的單軸長期強度為常規單軸壓縮試驗峰值強度的70%～80%；SCHMIDTKE 等[14]的研究也驗證了這一方法的可行性；李良權等[15]基于向家壩砂巖的常規三軸壓縮試驗和三軸壓縮流變試驗成果，指出巖石的裂紋損傷應力可反映長期強度所在應力水平。盡管國內外學者對于長期強度進行了大量的試驗研究和理論分析，但是通過力學模型直接獲得硬脆性巖石長期強度解析表達的研究甚少。

近年來，細觀試驗研究表明：微裂紋是硬脆性巖石非線性力學特性、蠕變變形和損傷的主要載體；微裂紋的擴展與當前的應力水平以及時間有著密切的關系[16?18]。在理論模型方面，SHAO等[19]指出硬脆性巖石內部微裂紋擴展方式包含：應力誘導的瞬時擴展和應力侵蝕導致的亞臨界(subcritical)擴展，并指出后一種方式引起了材料的時效變形行為，進而建立了與裂紋擴展有關的時效損傷本構模型，從細觀損傷角度為研究硬脆性巖石的時效特性提供了新的思路。在硬脆性巖石多尺度本構模型研究方面，朱其志等[20]基于Mori-Tanaka 均勻化方法和不可逆熱動力學理論，構建了多尺度力學損傷摩擦耦合模型；ZHAO等[21?22]從巖石細觀損傷累積和物理力學性質隨時間劣化的角度，建立了硬脆性巖石統一多尺度損傷模型，但是該模型對加速蠕變階段描述尚不充分，此外，模型參數的標定以及硬脆性巖石的長期強度研究均有待深入。本文作者在前期研究工作[19?23]的基礎上，從細觀角度出發，分析微裂紋擴展引起的宏觀變形響應，建立可充分描述硬脆性巖石蠕變三階段的多尺度損傷蠕變模型，并通過對摩擦損傷耦合和細觀時效損傷演化規律的研究，構建模型參數的跨尺度標定方法，推導硬脆性巖石長期強度的解析表達式，為深部巖體工程問題的長期穩定性和耐久性分析提供幫助。

1 多尺度損傷蠕變模型

微裂紋成核和擴展引起的損傷是導致硬脆性巖石變形和破壞的主要力學機制。研究表明[24?25]這類損傷包括在應力變化下的微裂紋成核和擴展產生的瞬時損傷以及常應力下微裂紋發生亞臨界擴展引起的時效損傷。據此，考慮到時間效應，對損傷d進行如下分解：

式中：di為瞬時損傷，dt為時效損傷。繼而基于均勻化方法、熱動力學理論和亞臨界擴展理論建立硬脆性巖石多尺度損傷蠕變模型。

1.1 Helmholtz自由能與狀態方程

硬脆性巖石宏細觀結構[18]示意圖如圖1所示。從硬脆性巖石的細觀結構出發(圖1(b))，選取如圖1(c)所示的代表性體積單元體(RVE)為研究對象，RVE的邊界為?Ω，整個RVE由線彈性固體基質和大量幣形微裂紋組成。根據ZHU等[26]的研究成果，含微裂紋RVE的Helmholtz自由能Π為基質彈性自由能和與微裂紋相關的塑性功之和：

圖1 硬脆性巖石宏細觀結構示意圖Fig.1 Schematic diagram of macroscopic and mesoscopic structure of hard brittle rock

式中：E為總應變；Ec為非彈性應變；Cm為各項同性基質彈性張量，Cm=2μmK+3kmJ，km和μm分別為基質體積壓縮模量和剪切模量，通過引入二階單位張量δ，四階張量Kijkl=(δikδjl+δilδjk)/2-δijδkl/3，Jijkl=δijδkl/3；Cb為四階塑性剛度張量，基于Mori-Tanaka均勻化方法可得其表達式[20]為

根據熱動力學基本理論，應變E、非彈性應變Ec、損傷d的共軛力分別為

式中：Σ為宏觀應力；Σc為作用在微裂紋表面的局部應力，可分解為偏應力Sc=K:Σc和球應力Pc=trΣc3；Dd為損傷驅動力。

1.2 摩擦準則和和瞬時損傷演化準則

在壓應力作用下，硬脆性巖石材料的破壞以壓剪類型為主，故可選取庫侖摩擦準則來描述非彈性變形的演化。假設裂紋面內局部應力分布是均勻的，庫侖摩擦滑動準則F可表示為[26]

式中：α為粗糙微裂紋面的摩擦因數。

出于簡化考慮，選取相關聯的流動法則，塑性勢函數G=F，根據正交法則計算非彈性應變增量

式中：λc為塑性乘子；V=Sc/‖Sc‖。

在不可逆熱動力學理論框架下，對于損傷行為，基于應變能釋放率的瞬時損傷演化準則可表述為[26]

式中：R(d)為損傷發展抗力，根據損傷抗力函數的基本特征[22]，其表達式如下：

式中：ξ=d/dc；dc為損傷臨界值，對應于材料的峰值應力處的損傷，當且僅當d=dc時，R(d)=R(dc)。

運用正交化準則計算瞬時損傷增量：

式中：λd為損傷乘子。結合式(4)和式(7)可以看出，損傷準則只與摩擦產生的非彈性應變Ec和當前損傷d有關，即損傷由微裂紋面上的摩擦滑移驅動。反之，損傷演化通過Cb對控制摩擦滑移的局部應力Σc產生影響。因此，摩擦和損傷的演化是強耦合的，可以通過聯立準則(5)和(7)的一致性條件方程求解塑性乘子λc和損傷乘子λd。

1.3 時效損傷演化準則

時效損傷變量dt是關于時間t和損傷d的函數，即

當t→∞時，dt→。為t時刻巖石細觀結構達到平衡狀態時的時效損傷。當dt<時，系統尚未達到平衡狀態，材料細觀結構將不斷向自平衡狀態發展。文獻[16]指出：細觀結構演化在動力學方面可以解釋為系統偏離平衡狀態的距離，即(-dt)，因此，時效損傷演化準則可用線性形式來描述[16?17]：

式中：γ為與材料性質有關的常數；?[0,∞)；

對于硬脆性巖石，當偏應力水平低于巖石長期強度時，巖石僅會發生衰減蠕變和穩態蠕變，并且穩態蠕變的速率很小，變形最終會趨于穩定；當偏應力水平高于巖石長期強度時，將會發生加速蠕變，并導致巖石迅速破壞。在微裂紋各向同性分布假定下，非彈性應變Ec可分解為

式中：標量β為非彈性體積應變，β=trEc；二階張量Γ=K:Ec為非彈性偏應變。

在應力空間，式(5)可表示為

式中：S=K:Σ；P= trΣ3。

由式(13)可知：當巖石受到的應力水平小于巖石初始屈服應力時(F<0)，損傷未發生；當應力水平高于初始屈服應力時(F=0)，在蠕變試驗過程中應力保持恒定，則和均為常數，即宏觀非彈性應變Ec與細觀損傷d的比值保持恒定。忽略彈性后效作用，這一理論結果與MOGI[27]和OHNAKA[28]通過試驗觀測到的結果相一致。因此，為反映硬脆性巖石在各蠕變階段的變形規律，對加速蠕變階段進行更深入的研究，在文獻[16]的基礎上引入細觀參數的表達式：

式中：B和n為模型參數，B>0，控制瞬態蠕變和穩態蠕變階段時效損傷的演化速度，n>0，控制加速蠕變階段時效損傷的演化速度。顯然，當損傷達到損傷臨界值dc時，加速蠕變啟動。

1.4 增量蠕變本構關系

在建立蠕變本構模型過程中，假定彈性應變全部是由應力變化產生，非彈性應變包括與時間無關的非彈性變形和與時間相關的非彈性變形，其中，與時間相關的非彈性變形由時效損傷所引起。

考慮時間效應后，當應力保持不變時，庫侖型摩擦準則式(5)依然適用。與時間相關的非彈性應變增量也可根據流動法則求得

式中：λct為與時間相關的塑性乘子。

根據摩擦準則的一致性條件：

求解可得：

綜上，蠕變本構關系的增量形式可寫為

式中：Sm=(Cm)-1為基質柔度張量。

2 長期強度研究

MARTIN等[12,29]研究發現：硬脆性巖石材料在壓縮荷載作用下的應力?應變關系可分為裂紋閉合階段、線彈性階段、穩定裂紋擴展階段和不穩定裂紋擴展階段[12]，分別對應巖石裂紋閉合應力Σcc、裂紋初始應力Σci、裂紋損傷應力Σcd和峰值強度Σf。裂紋損傷應力Σcd對應于巖石加載過程中體積應變?軸向應變曲線上體積變形的拐點。MARTIN等[12,14]研究表明：裂紋損傷應力Σcd可作為硬脆性巖石的長期強度。

由彈塑性理論可知，體積應變Ev包含彈性體積應變Ev,e和非彈性體積應變Ev,p

對于常規三軸壓縮試驗，忽略圍壓所產生的體積應變和裂紋閉合階段產生的非彈性應變，在軸向加載過程中

式中：P0為圍壓；Λc為偏壓引起的累計塑性乘子，表達式為[20,30]

ZHU[30]通過摩擦損傷耦合分析推導了在常規三軸加載情況下，摩擦準則(式(5))在主應力空間的表達形式為

將式(21)和式(22)代入式(20)并化簡可得

對損傷變量d進行求導，將=0(體積壓縮/膨脹轉化點)化簡，并結合邊界條件得

通過求解式(24)，得到滿足條件的解ξ∞(見圖2)，從而求得長期強度Σ∞對應的瞬時損傷為

圖2 損傷抗力和軸向應力演化示意圖Fig.2 Diagram of damage resistance and axial stress evolution

將式(25)代入式(22)，即可得到不同圍壓下硬脆性巖石長期強度的解析表達式為

3 模型驗證

為驗證所提出的多尺度損傷蠕變模型的準確性和有效性，本文采用該模型模擬向家壩砂巖在不同圍壓下的三軸壓縮蠕變試驗。李良權等[15]在圍壓分別為3，5 和7 MPa 下對向家壩砂巖進行了常規三軸壓縮試驗和三軸壓縮蠕變試驗，并對向家壩砂巖的長期強度進行了分析，其常規三軸壓縮試驗和三軸壓縮蠕變試驗的圍壓和偏壓加載速率均為7.5 MPa/min，其中，蠕變試驗采用多級加載方法，詳細的試驗方案和試驗曲線見文獻[15,31]。

3.1 模型參數跨尺度標定

本模型涉及的5個力學參數和3個時效參數均可通過室內試驗和理論推導確定。

1)基質的彈性模量Em和泊松比νm可通過常規三軸壓縮試驗應力?應變曲線線性段確定。

2)由式(22)得到巖石在受壓狀態下的強度為

顯然，式(27)可以與摩爾?庫侖強度準則直接關聯，繼而建立細觀參數(α和R(dc))與材料宏觀強度參數(內摩擦角?和黏聚力c)的跨尺度關系

3)損傷臨界值dc對應于應力峰值處的損傷，可通過聲發射試驗確定，即在常規三軸壓縮試驗時記錄不同圍壓下峰值應力處聲發射信號的累積計數。在本構模型中，dc在應力?應變關系上的宏觀表現為控制峰前和峰后非彈性變形，不影響巖石的峰值強度。LOCKNER[32]指出，dc與圍壓大致呈線性關系。針對向家壩砂巖，由于其圍壓差別不大，本文將dc取為常數。運用上述方法最終確定向家壩砂巖的5個力學參數如表1所示。

表1 基本力學參數Table 1 Basic mechanical parameters

4)當巖樣受到長期強度荷載時，瞬時損傷和時效損傷滿足如下關系

結合式(14)和式(28)，可以得到模型參數

將表1中的力學參數代入式(24)和式(25)，計算得到長期強度對應的瞬時損傷d∞=0.62。將dc和d∞代入式(30)，得到B=0.9。

5)蠕變參數γ可通過三軸壓縮蠕變試驗中瞬態蠕變的速率來確定[16]，其取值與材料的性質有關，一般，軟巖的γ取值在10-5～10-6量級，硬巖的γ取值在10-4量級，針對向家壩砂巖取γ=10-4。蠕變參數n可通過三軸壓縮蠕變試驗加速蠕變持續時長來確定，為簡便取n=1。

3.2 單級加載下長期強度預測

通過式(26)可以得到在單級加載情況下圍壓為3，5和7 MPa砂巖的長期強度。表2所示為多尺度模型確定的長期強度與李良權等[15]運用應力?應變等時曲線簇法和穩態蠕變速率?應力水平的關系曲線法確定的長期強度對比。從表2可以看出：除圍壓為3 MPa下多尺度模型對長期強度的預測結果較常用方法確定的長期強度偏小，圍壓為5 MPa 和7 MPa下多尺度模型預測的長期強度與上述2種方法得到的結果相一致。

表2 向家壩砂巖在不同圍壓下的長期強度Table 2 Long-term strengths of Xiangjiaba sandstone for different confining pressure MPa

圖3所示為不同圍壓、長期強度荷載下砂巖的蠕變演化模擬曲線。從圖3可以看到典型蠕變的三階段特征，即瞬態蠕變、穩態蠕變和加速蠕變，加速蠕變發生時間在200～250 h之間。

圖3 不同圍壓下長期強度荷載下砂巖的蠕變演化模擬曲線Fig.3 Simulation curves of creep evolution of sandstone under long-term strength with different confining pressures

圖4所示為圍壓5 MPa時，長期荷載作用下砂巖的軸向蠕變、軸向蠕變速率以及細觀損傷演化的模擬曲線。從圖4可見：細觀損傷演化曲線與蠕變曲線具有明顯的一一對應關系，尤其是根據臨界損傷可以很直觀地在三階段蠕變曲線中找到加速蠕變啟動點。此外，還可以發現蠕變應變的3個階段與應變速率緊密聯系。

圖4 圍壓5 MPa時長期強度荷載下砂巖的軸向蠕變、軸向蠕變速率以及細觀損傷演化的模擬曲線Fig.4 Simulation curves of creep,creep rate and microscopic damage of sandstone under long-term strength with confining pressure of 5 MPa

圖5所示為圍壓5 MPa時不同軸向荷載作用下砂巖的蠕變曲線。從圖5可以看到，當軸壓為155 MPa時，不會發生加速蠕變，即不會產生蠕變失效；當軸壓為157.1 MPa時，發生蠕變失效，且失效時間約為240 h；當軸壓為160 MPa 時，蠕變破壞時間縮短約40 h。說明當荷載低于長期強度時，巖石不會發生蠕變失效，而當荷載高于長期強度時，隨著軸壓增大，蠕變失效時間在逐漸減小，這與眾多試驗觀察結果相一致。

圖5 圍壓5 MPa時不同軸向荷載下砂巖的蠕變模擬曲線Fig.5 Simulation curves of creep evolution of sandstone under different axial load with confining pressure of 5 MPa

當圍壓為5 MPa，長期強度荷載下n=0.5，1.0，1.5 和2.0 時多尺度蠕變模型獲得的砂巖軸向蠕變和細觀損傷演化曲線如圖6所示。從圖6可以看出，參數n僅對加速蠕變段產生影響，隨著n增大，加速蠕變階段減少。換言之，參數n對蠕變失效時間有很大的影響。

圖6 不同參數n下砂巖的軸向蠕變和細觀損傷模擬曲線Fig.6 Simulation curves of creep and microscopic damage of sandstone with different parameter n

3.3 三軸壓縮蠕變試驗模擬

為進一步對多尺度損傷蠕變模型進行驗證，采用3.1 節的模型參數對圍壓為5 MPa 和7 MPa 時向家壩砂巖三軸壓縮蠕變試驗應變?時間全過程曲線進行數值模擬(圖中藍色曲線)，并與試驗數據進行對比，結果如圖7所示。從圖7可以看出，模型可以較為完整地描述不同圍壓下硬脆性巖石分級加載蠕變的力學行為，包括在低應力水平下衰減蠕變和穩態蠕變的特征、在高應力水平下加速蠕變破壞的特征。值得指出的是，用分級加載的方式，上一級加載的應力水平會對巖樣造成不同程度的損傷，而隨著加載級數的增多，損傷逐漸增加。在圖7(a)的模擬中，軸向應力為155 MPa等級下僅經歷約50 h后加速蠕變發生。而從圖5可以看出，在單級加載情況下，軸向應力為155 MPa等級下不會發生加速蠕變破壞。說明從巖石材料隨時間損傷累積和物理力學性質劣化的角度出發建立的多尺度蠕變損傷模型很好地考慮了多級加載下損傷的累積效應。此外，根據本文預測的長期強度可以對分級加載蠕變試驗的荷載分級提供一個可靠的參考依據。例如，在圍壓為5 MPa的分級加載蠕變試驗中，將軸向應力為155 MPa的荷載等級持續時間延長至約150 h，若巖樣未破壞，再加荷載至165 MPa；在圍壓為7 MPa的分級加載蠕變試驗中，在荷載等級為177 MPa與187 MPa之間建議增加一個持續時長為75 h的182 MPa荷載等級。此外，結合圖7的分析結果，通過調整參數n可以提高模型預測結果的準確性，如圖7中紅色曲線所示。

圖7 向家壩砂巖三軸壓縮蠕變試驗結果[15]與數值模擬對比Fig.7 Comparison between numerical simulation and experimental data[15]of triaxial compression creep test of Xiangjiaba Sandstone

4 結論

1)基于均勻化方法建立多尺度損傷蠕變模型是可行的，和通常的宏觀損傷力學蠕變模型相比，本模型所體現的力學機制更為明確，而且所含參數少，物理意義明確，并均可通過室內試驗和理論推導予以確定。

2)通過摩擦損傷耦合分析和細觀時效損傷演化規律研究，推導了硬脆性巖石長期強度的解析表達式。通過對向家壩砂巖長期強度預測和對比，驗證了運用多尺度模型確定長期強度的可行性，為多級蠕變試驗的荷載分級提供理論依據。

3)模型可以較為完整地描述硬脆性巖石在不同圍壓下分級加載和單級加載蠕變的力學行為，包括在低應力水平下衰減蠕變和穩態蠕變的特征、高應力水平下發生加速蠕變破壞的特征。通過多級加載和單級加載蠕變試驗的模擬對比表明模型可以較為準確地考慮到應變歷史對損傷演化的影響，說明基于損傷累積和物理性能隨時間劣化的角度出發建立蠕變模型的合理性。

4)需要指出的是本文所提出的是一種基于細觀力學的各向同性簡化模型，而在實際工程中，巖體的應力狀態往往較為復雜，各向異性損傷以及與微裂紋張開閉合有關的單邊效應都需要予以考慮。