基于Weibull分布的脆性巖石峰后應力-應變曲線特征及其統計損傷模擬方法

張 超，白 允，楊楚卿，雷 勇

（湖南科技大學巖土工程穩定控制與健康監測湖南省重點實驗室，湖南湘潭 411201）

引言

巖石工程已廣泛存在于公路、鐵路、礦山和土木等領域，各類巖石工程（如地下結構與隧道、邊坡、基礎工程，等等）的變形分析與計算已成為其結構設計不可或缺的重要內容，而合理的巖石本構模型是保證巖石工程穩定性控制和安全預測的關鍵。一般而言，應力水平是決定巖石表現為脆性變形還是延性變形的關鍵因素之一［1－2］，當圍壓低于脆延轉換臨界值時，巖石內部微裂紋迅速擴展、相互貫通，應力達到峰值點（極限強度）以后不斷發生跌落即應變軟化，其破壞模式表現為脆性；當圍壓高于脆延轉換臨界值時，巖石內部微裂紋仍發生擴展并相互貫通，但應力達到極限強度以后不發生跌落即應變硬化。由于巖石脆性變形破壞無明顯預兆，產生的危險性和破壞性極大，是巖石工程災害頻發的主要原因，因此，開展脆性巖石峰后應力-應變曲線特征及其統計損傷模擬方法研究對巖石工程結構災害防治工作具有重要的理論價值。

巖石剛性伺服力學試驗機研發成功以后，脆性巖石峰后應力-應變曲線得以實時測得，眾多學者對脆性巖石峰后變形力學特性的研究有了快速進展。根據巖石三軸壓縮試驗規程，可獲得不同圍壓作用下巖石全應力應變曲線，在此基礎上對巖石峰后破壞特征進行分析，揭示了巖石受荷變形破裂全過程的基本規律，尤其是巖石峰后應力跌落過程是連續的，完全破壞以后仍具有殘余承載能力［3］。三軸壓縮條件下巖石變形破壞全過程模擬方法研究也在不斷探索，由于巖石變形破壞伴隨著內部微裂紋不斷孕育、發展和相互貫通，當巖石發生損傷以后，巖石變形力學特性已超出彈性理論研究范疇，從連續介質力學角度分析巖石變形破裂的物理基礎已發生改變［4－5］，而巖石材料內部存在的大量微裂隙或微孔隙正是損傷力學理論所研究的范疇［6－7］，于是，基于有效應力概念和應變等效假設，將統計強度理論和連續損傷理論相結合，提出了統計損傷的概念，從而建立了脆性巖石統計損傷本構模型，其研究內容主要集中在巖石幾何損傷模型的構建方法［8］、巖石破壞準則和微元體強度隨機分布類型［9］、巖石發生損傷存在閾值的影響［10］以及本構模型參數的確定方法［11］等方面，從而不斷地完善巖石統計損傷軟化本構模型理論，使該類模型盡可能地能夠模擬脆性巖石變形破壞全過程。然而，目前該類模型理論曲線在描述巖石峰后應變軟化行為具有強烈的隨機性，即無法定量地模擬脆性巖石峰后應力-應變試驗曲線，這就導致難以利用該類模型預測脆性巖石達到極限強度以后的連續變形過程，但為文中開展脆性巖石峰后應力-應變曲線特征及其統計損傷模擬方法研究提供了一條可借鑒的途徑。

為此，本研究針對已有巖石統計損傷軟化本構模型在反映脆性巖石峰后應變軟化行為方面存在的缺陷，以微元體強度服從Weibull分布為前提，探討其分布參數對脆性巖石損傷演化規律產生的影響，提出一種新型幾何損傷模型，然后結合統計損傷的概念，建立能夠定量地模擬脆性巖石峰后應力-應變試驗曲線的統計損傷軟化本構模型，給出相應模型參數的確定方法，以期解決已有脆性巖石本構模型峰后理論曲線存在的強烈隨機性，進一步完善巖石統計損傷本構模型的研究內容與方法，為各類巖石工程結構變形與分析計算提供有力的本構模型理論依據，從而更好地進行巖石工程穩定性控制和安全預測。

1 脆性巖石統計損傷本構模型理論及分析

目前，Lemaitre提出的巖石損傷本構基本關系式［6］在巖石損傷力學理論中具有代表性，也得到了廣泛應用，即

式中：σ和ε分別為名義應力張量和名義應變張量；C和D為巖石材料的彈性矩陣和損傷變量。

巖石材料作為一種天然地質材料，其內部構造不均質從而使巖石各微元體強度數值不盡相同，因此，基于統計強度理論假定各微元體強度F服從Weibull分布［7］，即

式中：m和F0為微元體強度的分布參數。根據巖石三軸壓縮試驗，當巖石軸向承受較低荷載時，巖石內部少量的、低強度的微元體發生破壞，陰影面積較小，如圖1（a）所示，隨著軸向荷載不斷增大，微元體發生破壞的數目不斷增多，陰影面積不斷增大，如圖1（b）所示，當巖石所受應力超過極限強度以后，大量微元體開始迅速發生破壞，巖石峰后應力快速跌落，表現為應變軟化，但仍有少量的、強度極大的微元體不能發生破壞，由此巖石進入殘余強度變形階段，但仍具有殘余承載能力。由此可見，陰影面積的變化可用于描述巖石變形破壞過程中微元體破壞累積過程，則陰影面積可用于定義損傷變量D，即

圖1 微元體強度F概率密度函數Fig.1 Probability density function of microelement strength F

將式（3）代入式（1）可得巖石統計損傷軟化本構模型，即

由此可知，要使該模型能夠模擬巖石變形破裂全過程，需對模型參數進行確定。該模型參數包括巖石材料彈性矩陣C，度量微元體強度F的相關參數以及分布參數m和F0，前兩類參數可由巖石基本力學試驗獲得，參數物理意義明確，屬于試驗類定常數，因此關鍵在于分布參數m和F0的確定方法。

目前，分布參數m和F0的確定方法主要有線性擬合法［12］、優化反分析法［13］和峰值點法［14］等。線性擬合法首先將本構模型方程變換為關于待擬合參數的線性方程，然后利用巖石三軸壓縮試驗曲線進行線性擬合從而獲得特定圍壓作用下分布參數m和F0的數值，該方法參數獲取便捷，但不具有物理意義，而且不同圍壓對巖石損傷本構模型產生不同的影響，導致特定圍壓下分布參數值對其余圍壓下巖石應力-應變曲線的適用性較差。優化反分析法是將不同圍壓作為不同工況并建立目標函數以及根據經驗確定的參數變化范圍，構成一個有約束的非線性規劃問題，利用單純形法求解分布參數，但需要基于大量的試驗數據且獲得的參數解并不唯一。峰值點法要求本構模型方程需同時滿足如下條件：

式中：σ1c和ε1c分別為脆性巖石所受應力達到極限強度時的最大主應力及相應的主應變。由此可見，只需本構模型理論曲線經過試驗峰值點且具有相應的極值特性即聯立上述方程從而求解分布參數m和F0，該方法處理簡捷，物理意義明確，得到了廣泛應用。于是，根據上述參數確定方法即可獲得巖石統計損傷軟化本構模型理論曲線［11］，如圖2所示，該類模型能夠較好地反映巖石峰前變形特征和應變軟化特征，但無法反映殘余強度變形特征。于是，曹文貴等［15］認為巖石發生損傷破壞后仍具有殘余強度，從而建立了新型巖石損傷模型，在基礎上引入統計損傷的概念，建立了考慮巖石統計損傷軟化本構模型，使該類模型還能夠反映巖石殘余強度變形特征。然而，上述研究都不能定量地描述脆性巖石峰后應力跌落過程，模型理論峰后曲線與試驗峰后曲線存在較大偏離現象，究其原因在于2方面：一方面是統計損傷本構模型理論本身存在問題，該問題非常復雜，這里暫不討論；另一方面是認為本構模型理論正確的前提下，其模型參數確定方法存在的缺陷所造成的，這是建立本構模型的重要組成部分，現對此進行分析。

圖2 統計損傷軟化本構模型理論曲線Fig.2 Curves of statistical damage softening constitutive model

如上所述，脆性巖石統計損傷軟化本構模型參數［16－17］可分兩類：一類為巖石材料彈性矩陣C（如彈性模量E、泊松比μ、內摩擦角φ和黏聚力c）和度量微元體強度F的相關參數（如M-C強度準則α和κ），其數值主要取決于巖石物理性質，屬于基本力學參數，它們不會因本構模型或模型理論的改變而發生變化，因此，該類模型參數并非導致模型理論峰后曲線與試驗峰后曲線存在較大偏離的原因；另一類為微元體強度服從Weibull分布的參數m和F0，廣泛采用峰值點法進行確定，但峰值點法僅僅要求模型理論曲線需經過試驗曲線峰值點且滿足極值特性，由此獲得的2個方程無法對脆性巖石峰后應力-應變曲線特征產生任何約束或控制，致使模型理論峰后曲線存在強烈的隨機性，無法定量地模擬巖石峰后應力跌落過程，同時由該方程組獲得的分布參數m和F0數值具有唯一性，僅適用于模擬特定圍壓作用下脆性巖石應力-應變試驗曲線。

因此，將在探討分布參數m和F0對脆性巖石損傷演化規律產生的影響基礎上，提出一種新型脆性巖石損傷模型，然后結合統計損傷的概念，建立能夠定量地模擬脆性巖石峰后應力-應變試驗曲線的統計損傷軟化本構模型。

2 新型脆性巖石損傷模型

由前述分析可知，脆性巖石統計損傷軟化本構模型理論曲線的變化規律除受巖石材料基本力學參數的影響以外，還受損傷變量D變化的影響，但是巖石材料基本力學參數屬于試驗類定常數，而損傷變量D的變化規律主要與分布參數m和F0密切相關，因此，分布參數m和F0對模型理論曲線的變化規律產生重要影響，故基于式（3）探討分布參數m和F0對損傷變量D產生的影響，如圖3所示，可得：

圖3 分布參數m和F0對損傷變量D的影響Fig.3 Effects of distribution parameters m and F0 on damage variable D

（1）隨著微元體強度F的增大，損傷變量D不斷增大且曲線形狀表現為“S型”，損傷變量D的數值變化區間為[0，1]。

（2）分布參數F0保持不變，隨著分布參數m的增大，損傷變化率（?D/?F）不斷增大，D-F曲線形狀保持為“S型”并繞著定點A作逆時針轉動。根據巖石損傷本構基本關系式即式（1）可知，損傷變化率增大使（1-D）的數值快速降低，巖石峰后應力跌落加快，其破壞模式表現為脆性程度增大。

（3）分布參數m保持不變，隨著分布參數F0的增大，損傷變化率（?D/?F）不斷減小，D-F曲線形狀保持為“S型”且趨于向微元體強度F增大方向移動。根據巖石損傷本構基本關系式即式（1）可知，損傷變化率減小使（1-D）的數值慢速降低，巖石峰后應力跌落減慢，其破壞模式表現為脆性程度減小。

由此可見，分布參數m和F0對損傷變量D的變化規律產生重要影響，進而影響巖石峰后應力跌落即峰后應力-應變曲線特征，然而根據特定圍壓作用下巖石三軸壓縮試驗數據，由峰值點法確定的分布參數m和F0數值具有唯一性，使現有巖石統計損傷軟化本構模型針對特定圍壓情況僅能夠產生一條理論曲線，而且此條理論曲線由于峰值點法的原因而峰后表現出強烈的隨機性，無法定量地模擬巖石峰后應力跌落過程。然而，峰值點法處理簡捷，參數物理意義明確即分布參數m為巖石材料的均勻系數，分布參數F0為微元體強度統計均值，因此，本研究在繼承該方法優越性的前提下，經過不斷地嘗試和驗證，擬引入參數λ和η，提出新型脆性巖石損傷模型，即

式中：σi、σ′i和σr分別為脆性巖石表觀應力、有效應力和殘余應力。在分布參數m和F0數值保持唯一性的情況下，參數λ和η能否使在此基礎上建立的統計損傷軟化本構模型定量地模擬出脆性巖石峰后應力跌落過程，將在第4節進行分析和驗證。

3 脆性巖石變形破裂過程統計損傷模擬方法

根據損傷力學理論可知，脆性巖石可抽象為由損傷部分和未損傷部分組成，而未損傷部分的變形力學特性服從廣義虎克定律［18］，即

式中：E和μ分別為脆性巖石彈性模量和泊松比；σ′i、σ′j和σ′k分別為未損傷部分在i、j和k方向的有效應力；εi為脆性巖石在i方向的應變。由于三軸壓縮試驗條件下巖石試件在軸向荷載壓縮方向發生破裂變形并逐漸喪失承載能力，圍壓約束巖石試件發生側向變形并起到提高軸向承載能力的作用，因此可近似認為巖石試件損傷主要發生在軸向，忽略側向圍壓方向發生的損傷。于是，式（8）可改寫為

于是，將式（9）代入式（7）可得

其中，

由于脆性巖石變形破裂行為的本質是微元體發生破壞的數目不斷增多的過程，因此，本文在考慮脆性巖石力學特性受應力狀態的影響基礎上，將M-C強度準則［17］作為判斷微元體是否發生破壞的依據，即

其中，

式中：φy和cy分別為脆性巖石所受應力達到損傷閾值時的內摩擦角與黏聚力。根據損傷變量D變化規律可知，當脆性巖石僅發生彈性變形時，微元體未發生破壞，損傷變量D等于0；當脆性巖石發生損傷變形時，微元體開始發生破壞且破壞數目不斷增多，損傷變量D不斷增大；當脆性巖石變形持續增大而應力不再增加時，微元體破壞數目已達上限，損傷變量D等于1/λ并保持不變。于是，新型脆性巖石損傷模型下損傷變量D可表述為

式中：Fs為損傷變量D恰好達到1/λ時的微元體強度值。由于脆性巖石損傷演化過程是連續的，由式（15）可確定出Fs，即

于是，將式（15）代入式（10）可得脆性巖石變形破裂過程統計損傷模擬方法，即

要使式（17）能夠模擬出脆性巖石變形破裂過程，還需給出分布參數m和F0以及參數λ和η的確定方法。于是，本研究考慮峰值點法的優越性并繼承其參數確定方法，可得分布參數m和F0的解析表達式，即

其中，

由于參數λ和η變化對脆性巖石損傷演化規律及其損傷本構模型理論曲線產生的影響未知，因此暫無法給出其確定方法，需參數分析之后給出相應表達。值得注意的是，現有脆性巖石統計損傷軟化本構模型是最大主應力σ1和最大主應變ε1之間的理論曲線，而三軸壓縮條件下脆性巖石實測曲線是軸向應力σ1t和軸向應變ε1t之間的試驗曲線，兩者曲線之間存在明顯差異，需給出它們之間的理論關系，即

式中，εc為靜水壓力下巖石產生的初始應變，其數值為（1-2μ）σ3/E。于是，將式（22）和式（23）代入式（17）可得脆性巖石修正統計損傷軟化本構模型。

4 參數（λ和η）分析及本構模型驗證

YUMLU.M等［19］基于三軸壓縮試驗機給出了圍壓分別為3，5和8 MPa作用下砂巖軸向應力與軸向應變試驗曲線，如圖4所示，砂巖彈性模量E=27.0 GPa，泊松比μ=0.25，內摩擦角φy=50°，黏聚力cy=13.17 MPa。為了清楚地表明已有脆性巖石統計損傷軟化本構模型理論曲線存在的缺陷以及采用本文模型和方法所帶來的優越性，首先將參數λ和η均賦數值為1，本文模型則轉化為已有脆性巖石損傷本構模型［15］，將其模型理論曲線與試驗曲線進行比較，然后再對參數λ和η進行分析，討論其影響并給出其確定方法，最后對本文模型和方法的可行性與合理性進行驗證。

4.1 不考慮參數λ和η影響的脆性巖石統計損傷軟化本構模型分析

為了不考慮參數λ和η的影響，將參數λ和η均賦數值為1，則本文模型轉化為已有脆性巖石損傷本構模型［15］，根據已有方法可得其模型理論曲線σ1－ε1、修正模型理論曲線σ1t－ε1t以及砂巖損傷演化曲線D－ε1t，并將已有理論曲線與砂巖試驗數據進行比較，如圖4所示，可以看出：

圖4 已有模型理論曲線［15］（λ=η=1）與砂巖試驗數據對比Fig.4 Comparison of experimental data of sandstone and theoretical curves of existing model［15］（λ=η=1）

（1）已有模型理論曲線σ1－ε1能夠定性地反映砂巖線彈性、屈服硬化、應變軟化和殘余強度等4個變形階段，但存在明顯缺陷，對于峰值強度之前，隨著圍壓的增大，已有模型理論曲線與砂巖試驗數據的偏離程度不斷增大，對于峰值強度之后，已有模型理論曲線明顯偏離砂巖試驗數據，且與圍壓水平無關。

（2）修正模型理論曲線σ1t－ε1t能夠定量地描述砂巖峰值強度之前的試驗曲線，盡管峰值強度之后的修正模型理論曲線與砂巖試驗曲線的偏離程度減小，但兩者之間的差距仍然不可忽略或近似。

（3）砂巖損傷演化曲線D－ε1t形狀基本表現為“S型”，當砂巖未發生損傷即處于彈性變形階段時，損傷變量D為0，當砂巖進入屈服硬化和應變軟化階段時，損傷變量D由0逐漸增大，直至砂巖進入殘余強度變形階段，損傷變量D等于1并保持不變，可見D－ε1t能夠對砂巖階段性變形特征進行合理闡釋。

由此可見，不考慮參數λ和η影響的脆性巖石統計損傷軟化本構模型及其修正模型存在明顯缺陷，它們均無法很好地模擬脆性巖石峰后應力跌落過程，因此為了解決該問題，文中提出了新型脆性巖石損傷模型即式（7），進而在此基礎上建立了脆性巖石修正統計損傷軟化本構模型。

4.2 參數λ和η分析

為了探討參數λ和η的影響并給出其確定方法，文中采用控制單參數變量方式即將參數λ和η分別賦數值為1，考慮單參數的變化對脆性巖石損傷演化規律及其損傷本構模型理論曲線產生何種影響，分別如圖5和圖6所示。

圖5 參數λ對模型理論曲線和損傷演化規律的影響（η=1）Fig.5 Effects of parameterλon theoretical curves and damage evolution law（η=1）

圖6 參數η對模型理論曲線和損傷演化規律的影響（λ=1）Fig.6 Effects of parameterηon theoretical curves and damage evolution law（λ=1）

（1）不同圍壓作用下峰值強度之前模型理論曲線不受參數λ和η變化的影響，繼承了已有修正理論模型（λ=η=1）的優點，能夠很好地反映砂巖線彈性變形階段和屈服硬化變形階段，不同圍壓作用下峰值強度之后模型理論曲線受參數λ和η變化影響顯著。

（2）控制參數η不變（賦參數η數值為1），隨著參數λ的增大，峰值強度之后的模型理論曲線應力跌落過程加快，殘余應力降低，損傷變化率增大，損傷演化理論曲線形狀基本保持為“S型”并繞固定點A作逆時針轉動，從而表明參數λ與Weibull分布參數m對損傷演化規律產生的影響一致。

（3）控制參數λ不變（賦參數λ數值為1），隨著參數η的增大，峰值強度之后的模型理論曲線應力跌落過程減慢，殘余應力增大，損傷變化率增大，損傷演化理論曲線形狀基本保持為“S型”并趨于向軸應變減小方向移動，從而表明參數η與Weibull分布參數F0對損傷演化規律產生的影響接近相同。

由此可見，參數λ和η對損傷演化理論曲線的影響與Weibull分布參數m和F0對損傷演化理論曲線的影響是等效的，這表明基于峰值點法確定的Weibull分布參數m和F0數值在保持不變的情況下，通過改變參數λ和η的數值來影響損傷演化理論曲線變化規律，進而使統計損傷軟化本構模型能夠模擬脆性巖石峰值強度之后的應力跌落過程，這一研究思路是可行的。

根據參數λ和η分析可知，參數λ和η的變化對不同圍壓作用下峰值強度之前模型理論曲線不產生影響，但對峰值強度之后模型理論曲線變形特征起著關鍵作用。由于巖石峰值后區變形特征可由脆性指標和殘余應力進行表征，因此可提出參數λ和η的確定方法如下：

張超等［20］綜合考慮全應力應變曲線峰前峰后對巖石脆性程度的共同影響，提出巖石的新脆性指標B L，可表示為

其中，

式中：ε1r為巖石所受應力為殘余應力時的應變。B L數值越大，代表巖石脆性程度越高，由式（24）～式（25）可得ε1r的表達式，即

將式（26）代入式（12）可得巖石所受應力恰好達到殘余應力σr時的Fs，即

于是，將式（27）代入式（15）可得微元體強度達到Fs時的巖石損傷變量Ds，即

根據連續性損傷條件，由式（15）可得參數λ的確定方法，即

根據巖石殘余強度變形階段特征可得

于是，將式（17）第3式代入式（30）可得參數η和參數λ的數值相等。

4.3 考慮參數λ和η影響的脆性巖石修正統計損傷軟化本構模型驗證

根據上述參數的確定方法可得出文中理論模型相關參數，如表1所示，將其代入考慮參數λ和η影響的脆性巖石修正統計損傷軟化本構模型，可得不同圍壓下砂巖變形破壞全過程理論曲線，并將其與砂巖試驗數據進行比較，如圖7所示。

圖7 本文修正理論模型曲線與砂巖試驗數據對比Fig.7 Comparison of experimental data of sandstone and theoretical curves of proposed model

表1 理論模型相關參數Table 1 Related parameters of theoretical model

由此可見，本研究修正理論模型能夠準確地描述砂巖峰值強度之前線彈性和屈服硬化等變形階段，繼承了已有修正理論模型（λ=η=1）的優點，也能夠較好地模擬出砂巖峰值強度之后試驗曲線即應變軟化和殘余強度變形階段，盡管理論曲線應變軟化階段與砂巖試驗數據尚存些許偏離，究其原因，本文模型和方法不適用于脆性指標趨于極大的極端情況，也未考慮砂巖初始損傷所帶來的影響，但它明顯縮小砂巖峰值強度之后應力跌落過程與試驗數據的偏離程度，表明本文模型和方法具有一定的合理性與可行性。

5 結論

針對已有脆性巖石統計損傷軟化本構模型在模擬脆性巖石峰值強度之后應力跌落過程存在強烈的隨機性，提出新型脆性巖石損傷模型，在此基礎上建立脆性巖石變形破裂過程統計損傷模擬方法，可得如下結論：

（1）在模擬三軸壓縮試驗條件下軸向應力-軸向應變曲線時，需對脆性巖石統計損傷軟化本構模型進行修正，即對最大主應力-最大主應變與軸向應力-軸向應變兩者之間的關系進行轉換。

（2）Weibull分布參數m和F0能夠使理論模型很好地模擬出脆性巖石峰值強度之前變形過程，但無法對脆性巖石峰值強度之后變形過程起到約束作用，從而使峰后理論模型曲線存在強烈的隨機性。

（3）參數λ和η變化對峰值強度之前理論模型曲線不產生影響，但對峰值強度之后理論模型曲線產生強烈影響，它們與分布參數m和F0對損傷演化規律產生的影響等效，即參數λ與分布參數m對損傷演化規律產生的影響一致，參數η與分布參數F0對損傷演化規律產生的影響接近相同。

（4）本文修正理論模型不僅能夠準確地脆性巖石峰值強度之前變形過程，也能夠明顯縮小脆性巖石峰值強度之后應力跌落過程與試驗曲線的偏離程度，表明本文模型和方法具有一定的合理性與可行性。