李修磊 ，陳洪凱 ，張金浩

（1. 重慶交通大學(xué)山區(qū)公路水運交通地質(zhì)減災(zāi)重慶市高校重點實驗室, 重慶 400074；2. 重慶交通大學(xué)土木工程學(xué)院, 重慶 400074）

巖石由于內(nèi)部大量隨機(jī)分布的微孔隙、微裂紋等初始缺陷，使得其力學(xué)性質(zhì)和破壞機(jī)制變得非常復(fù)雜. 為了研究巖石的力學(xué)變形特性，國外學(xué)者對各類巖石（砂巖、花崗巖、大理巖、頁巖等）開展了大量單軸和三軸試驗研究[1-4]，結(jié)果表明，巖石的變形破壞特征對圍壓具有很強(qiáng)的依賴性，構(gòu)建合理描述巖石變形破壞過程的本構(gòu)模型是巖石力學(xué)研究的主要內(nèi)容之一.

基于試驗研究，許多學(xué)者對巖石的力學(xué)變形特性進(jìn)行了大量的理論研究，也建立了相應(yīng)的本構(gòu)關(guān)系. Krajcinovic 等[5]首次將連續(xù)損傷和統(tǒng)計強(qiáng)度理論引入到巖石力學(xué)與破壞機(jī)制的研究中，隨后國內(nèi)外學(xué)者[6-10]也相繼建立了一系列的巖石統(tǒng)計損傷模型，基于Lemaitre應(yīng)變的等價理論[11]認(rèn)為巖石損傷的本質(zhì)是其內(nèi)部形成缺陷（如空洞、裂隙等），而缺陷的承載能力幾乎為0. 但該類模型有兩個方面的不足：一是無法反映荷載應(yīng)力較小時巖石的非線性彈性變形；二是難以準(zhǔn)確描述巖石變形峰值強(qiáng)度后殘余強(qiáng)度階段的變形特征. 針對以上兩方面不足，曹文貴等[12-13]基于微元強(qiáng)度的概念，認(rèn)為微元強(qiáng)度服從Weibull隨機(jī)分布，且微元損傷后以殘余強(qiáng)度的形式承擔(dān)荷載，考慮巖石破壞后的殘余強(qiáng)度變形特征，提出了能夠模擬巖石峰值后破壞變形特征的損傷統(tǒng)計模型；劉冬橋等[14]依據(jù)三軸試驗結(jié)果提出了巖石損傷變量的演化方程，基于微元強(qiáng)度概念發(fā)展了相應(yīng)的損傷本構(gòu)模型，用于描述巖石的應(yīng)變軟化；Li等[15]考慮了圍壓對巖石脆性指標(biāo)的影響，認(rèn)為微元強(qiáng)度服從Weibull分布、正態(tài)分布和冪函數(shù)分布，分別建立了對應(yīng)的巖石損傷本構(gòu)關(guān)系. 上述幾種模型均認(rèn)為巖石在屈服前呈線性彈性變形，與實際情況差別較大. 另外，曹文貴等[16]通過考慮荷載作用下巖石內(nèi)部空隙體積的變化和損傷閾值，建立了巖石的統(tǒng)計損傷模型，一定程度上反映了低應(yīng)力水平下巖石的非線性變形特征，但不能很好反映巖石破壞過程的變形規(guī)律. Xu等[17]雖考慮了巖石內(nèi)部空隙的壓密，并基于斷裂損傷理論建立了相應(yīng)的本構(gòu)模型，但峰值后的破壞變形與試驗結(jié)果差別較大.

基于上述分析，大部分巖石本構(gòu)模型能夠合理描述巖石部分階段的變形特征，很少有模型能夠準(zhǔn)確模擬巖石的變形全過程，關(guān)鍵在于現(xiàn)有模型對巖石的力學(xué)變形機(jī)理難以準(zhǔn)確反映，本文將根據(jù)荷載作用下巖石的變形全過程特征，分析巖石的力學(xué)變形機(jī)理，基于統(tǒng)計損傷理論，以期建立巖石的損傷統(tǒng)計本構(gòu)模型用于合理模擬巖石變形破壞的全過程.

1 巖石變形力學(xué)機(jī)理與分析

1.1 巖石變形過程的力學(xué)機(jī)理

外部荷載作用下巖石的變形可劃分為兩部分：一是由巖石骨架產(chǎn)生；二是由巖石內(nèi)部空隙產(chǎn)生. 較小荷載應(yīng)力下巖石骨架和空隙部分將同時產(chǎn)生變形；當(dāng)荷載應(yīng)力增加到一定程度時，巖石內(nèi)部空隙閉合完成，此時巖石僅發(fā)生線彈性的骨架變形；荷載應(yīng)力進(jìn)一步增加超過巖石屈服強(qiáng)度后，巖石發(fā)生非線性的骨架變形，出現(xiàn)屈服、應(yīng)變軟化和完全破壞的變形階段. 巖石破壞變形全過程示意如圖1所示. 圖中：Δ ε1為骨架部分應(yīng)變與總應(yīng)變之差；Δ εa為初始時刻骨架部分應(yīng)變與總應(yīng)變的差值（即為巖石總的空隙應(yīng)變）；εa和 σa為空隙完全閉合時的軸向應(yīng)變和應(yīng)力；σc和 εc分 別 為 峰 值 應(yīng)力和峰值應(yīng)變；σ1?σ3為偏應(yīng)力，σ1為大主應(yīng)力，σ3為小主應(yīng)力（圍壓）； ε1為軸向應(yīng)變. 由圖1可知：OA為初始空隙壓密段、AB為線彈性變形段、BC為屈服段、CD為應(yīng)變軟化段和D點后為完全破壞階段.

圖1 巖石破壞變形全過程示意Fig. 1 Whole failure and deformation process of rocks

1.2 巖石內(nèi)部空隙部分變形的分析方法

為了分析荷載應(yīng)力下由巖石內(nèi)部空隙部分產(chǎn)生的變形，在巖石內(nèi)部取一個有代表性的柱體單元，如圖2所示，空白處表示空隙部分，陰影處表示骨架部分. 設(shè)加載前柱體單元初始總高度為h0，巖石骨架部分的高度為h0s，空隙部分的高度為h0g；若某荷載應(yīng)力 σ下巖石柱體單元總變形量為 Δh，巖石骨架和內(nèi)部空隙分別產(chǎn)生的變形量為 Δhs和 Δhg，則

圖2 巖石變形分析模型Fig. 2 The deformation analysis model of rocks

荷載應(yīng)力 σ作用下巖石的總應(yīng)變 ε、骨架部分應(yīng)變 εs以及空隙部分應(yīng)變 εg可分別表示為

由式（4）可知：若要分析荷載應(yīng)力 σ作用下的 εg，須先獲得空隙部分的變形量 Δhg. 為此，將 σ劃分為由n個等級組成的應(yīng)力增量 Δ σt逐級施加，即

將應(yīng)力增量 Δ σt作用下巖石骨架部分和空隙部分 的 變形量分別記作 Δhts和 Δhtg，相應(yīng)的應(yīng)變增量分 別記 作 為 Δ εts和 Δ εtg，巖石 的 總應(yīng) 變 增量 Δεt=Δεts+Δεtg，骨架部分和空隙部分的總變形量 Δhs和Δhg可 分別視為 Δhts和 Δhtg的 累加. 巖石總 應(yīng)變、骨架以及空隙部分的應(yīng)變可表示為

令應(yīng)力增量 Δ σt作用下空隙部分的應(yīng)變增量Δεtg所占巖石總應(yīng)變增量 Δ εt的比例為kt（即 Δεtg=ktΔεt），則 Δ εts=(1?kt)Δεt. 若應(yīng)力增量 Δ σt和骨架部分的應(yīng)變增量 Δ εts服從廣義Hook定律，則

式中：Δ εts,i為骨架部分在i方向上的應(yīng)變增量；E和μ分別為巖石骨架的彈性模量和泊松比，其中，E即為巖石應(yīng)力-應(yīng)變線彈性階段的變形模量；i=1，2，3，j=2，3，1，k=3, 1, 2，分別表示三維空間上的主應(yīng)力和主應(yīng)變的方向；Δ εi,t為應(yīng)力增量 Δ σt作用下i方向上的總應(yīng)變增量，余同；Δ σi,t為主應(yīng)力i方向的應(yīng)力增量，余同.

利用式（6）和式（7），對式（8）等號兩側(cè)進(jìn)行求和可得

式中：εi為i方向上的總應(yīng)變；K為空隙應(yīng)變比，如式（10）；Δ εi,tg為應(yīng)力增量作用下i方向上的空隙應(yīng)變增量，εi,g為i方向上的空隙應(yīng)變.

1.3 空隙應(yīng)變比K的確定

由巖石的軸向應(yīng)力-應(yīng)變試驗曲線（如圖3）可知：初始空隙壓密階段的軸向應(yīng)力-應(yīng)變曲線表現(xiàn)出明顯的非線性，此階段軸向剛度隨著應(yīng)力的增加在逐漸增大；線彈性階段的軸向剛度隨應(yīng)力的增加近似為定值；若不考慮應(yīng)力-應(yīng)變曲線的微小波動，則空隙完全閉合的應(yīng)力點對應(yīng)于由非線性增長向線性增長過渡的轉(zhuǎn)折點.

將圖3中線彈性階段（AB）的試驗數(shù)據(jù)進(jìn)行回歸分析，即可得到圖3中藍(lán)色直線的線性方程為

圖3 應(yīng)力-應(yīng)變試驗曲線Fig. 3 Stress-strain test curve

式中：b為斜率，其值等于巖石的彈性模量E，MPa；c為縱軸（偏應(yīng)力軸）上的截距，MPa；b和c均為正值.

巖石的軸向總空隙應(yīng)變 Δ εa等于圖3中藍(lán)色直線在橫軸（應(yīng)變軸）上的截距，可表示為

偏應(yīng)力水平相同時，利用式（11）求解得到的軸向應(yīng)變減去對應(yīng)的試驗軸向應(yīng)變，可得到應(yīng)力相同時圖3中藍(lán)色直線與試驗曲線之間的 Δ ε1，如式（13）.

利用式（13）求解并繪制 Δ ε1與偏應(yīng)力之間的關(guān)系曲線，如圖4所示. 由圖4可確定隨著偏應(yīng)力增大軸向應(yīng)變差初始等于0的A點，該點對應(yīng)著空隙閉合完成時的軸向應(yīng)變 εa，A點之前則對應(yīng)著初始空隙壓密階段，B點對應(yīng)著屈服應(yīng)力點.

圖4 軸向應(yīng)變差-偏應(yīng)力曲線Fig. 4 Axial strain difference -deviatoric sress curve

巖石骨架部分的應(yīng)變可用式（11）在c值為0時進(jìn)行求得，也就相當(dāng)于圖3中的藍(lán)色直線向左平移c/b個單位長度，使其通過坐標(biāo)原點. 巖石空隙部分的軸向空隙應(yīng)變ε1g如式（14）. 當(dāng) ε1≥εa時，巖石內(nèi)部空隙完全閉合，空隙應(yīng)變不再發(fā)生變化.

利用多功能試驗機(jī)對砂巖試樣（φ 50 × 100 mm）開展三軸壓縮試驗，得到了不同圍壓下的軸向應(yīng)力-應(yīng)變關(guān)系曲線，如圖5所示. 由應(yīng)力-應(yīng)變試驗曲線，利用式（11） ～ （14）可以得到 ε1g與 ε1之間的變化關(guān)系，如圖6所示. 由圖6可以看出：砂巖的空隙應(yīng)變隨著軸向應(yīng)變的增加而增大，且增加幅度逐漸減小，應(yīng)變水平較大時逐漸趨近于定值；應(yīng)變水平較小時，不同圍壓下砂巖試樣的空隙應(yīng)變非常接近；隨著軸向應(yīng)變的增加，不同圍壓下砂巖空隙應(yīng)變的差異性逐漸顯現(xiàn)，應(yīng)變水平較大時近似為定值的空隙應(yīng)變隨著圍壓的增加呈現(xiàn)出先增大后減小的變化趨勢.

圖5 砂巖的三軸試驗結(jié)果Fig. 5 Triaxial test results of sandstone

圖6 軸向空隙應(yīng)變隨軸向應(yīng)變的變化規(guī)律Fig. 6 Variation of void strain with axial strain

利用式（14）對圖6中的試驗曲線進(jìn)行擬合，可得到空隙應(yīng)變比K隨 ε1變化的表達(dá)式為

式中：a1和a2均為計算參數(shù)，其中，a1為圖3中直線方程式（11）與軸向應(yīng)變軸的交點，即為c/b，a1可通過對巖石應(yīng)力-應(yīng)變曲線的線彈性段回歸分析獲得，a2可利用式（14）對確定的 ε1g- ε1關(guān)系曲線（見圖6）擬合獲得.

2 巖石損傷本構(gòu)模型的建立

2.1 損傷模型的構(gòu)建

基于Lemaitre [11]提出的應(yīng)變等價性，將巖石承受的宏觀名義應(yīng)力和凈應(yīng)力與其有效承載面積的減小建立聯(lián)系. 假定巖石由眾多微元均勻組成，所有微元在荷載作用下劃分為損傷和未損傷兩部分，且荷載由這兩部分共同承擔(dān). 圖7給出了巖石損傷的轉(zhuǎn)化過程. 圖中：空白部分的面積為S1表示未損傷部分，陰影部分的面積為S2表示損傷部分；σi為巖石整體受到的宏觀名義應(yīng)力：σ?i為未損傷部分受到的凈應(yīng)力；Rs為損傷部分受到的凈應(yīng)力. 初始時刻總面積S=S1，S2=0；巖石完全破壞時S=S2，S1= 0；損傷變量D如式（16）所示.

圖7 損傷轉(zhuǎn)換過程示意Fig. 7 Sketch of the damage transition process

取巖石微元進(jìn)行分析，由靜力平衡條件可得外部荷載為

聯(lián)立式（16）和式（17），可得

考慮巖石內(nèi)部空隙壓密的變形特征，根據(jù)廣義Hook定律，由式（9）和式（10）可得

將式（18）代入式（19）中，考慮常規(guī)三軸試驗條件下有 σ2=σ3，通過整理可得到反映巖石全應(yīng)力-應(yīng)變特征的本構(gòu)方程，如式（20）.

2.2 損傷部分殘余強(qiáng)度的計算

隨著外荷載的增加，巖石內(nèi)部有效微元數(shù)目逐漸減少，未損傷區(qū)域S1逐漸轉(zhuǎn)化為損傷區(qū)域S2，直到完全損傷. Menendez等[18-19]認(rèn)為巖石中微裂紋貫通形成剪切帶后，其強(qiáng)度主要依賴于剪切帶上的摩擦作用，黏聚力幾乎完全消失. 因而，可采用殘余強(qiáng)度Re來替代Rs，利用Mohr-Coulomb強(qiáng)度準(zhǔn)則計算：

式中：cr和 φr分別為巖石殘余強(qiáng)度對應(yīng)的黏聚力和殘余內(nèi)摩擦角.

由巖石三軸試驗結(jié)果可得不同圍壓下巖石的殘余強(qiáng)度Re，進(jìn)而得到殘余強(qiáng)度參數(shù).

2.3 損傷變量的演化

巖石內(nèi)部的缺陷會削弱其承載能力，這些缺陷在巖石內(nèi)部可以看作是隨機(jī)分布. 因此，從統(tǒng)計損傷的角度出發(fā)，認(rèn)為巖石損傷是一個連續(xù)的應(yīng)力過程，采用微元強(qiáng)度分布對巖石進(jìn)行定量分析. 現(xiàn)有研究大多認(rèn)為巖石微元強(qiáng)度服從Weibull函數(shù)分布[7-10,12,17,20]，本文采用相同的微元強(qiáng)度分布規(guī)律，相應(yīng)的概率密度函數(shù)為

式中：m和 ε0分別為形狀參數(shù)和尺寸參數(shù).

外部荷載作用下，應(yīng)變達(dá)到一定值時巖石內(nèi)部損傷區(qū)域為

由式（16）和式（23）可得巖石損傷變量的演化方程為

將式（15）和式（24）代入式（20）中，即可得到基于Weibull分布的巖石統(tǒng)計損傷本構(gòu)關(guān)系如下：

1） 當(dāng) ε1<εa時

2） 當(dāng) ε1≥εa時

為了得到上述Weibull分布的參數(shù)m和 ε0，對式（25）進(jìn)行整理后求對數(shù)，得到如下形式：

1） 當(dāng) ε1<εa時

2） 當(dāng) ε1≥εa時

式（26）右側(cè)可以看作是因變量，左側(cè)第一項lnε1可看作是自變量，m可視為斜率，?mlnε0為截距. 通過對現(xiàn)有不同圍壓下巖石的實測試驗數(shù)據(jù)，利用式（26）進(jìn)行線性回歸，可得到Weibull分布的參數(shù)m和 ε0.

圖8給出了不同m值對應(yīng)的損傷變量D隨ε/ε0的變化關(guān)系. 由圖可看出：D隨 ε/ε0的增加而增大； ε/ε0<1 時，m值越大意味著同等應(yīng)力水平下巖石的損傷程度越小，且 ε/ε0>1 時，m值越大意味著同等應(yīng)力水平下巖石的損傷程度越大，巖石越快達(dá)到完全損傷狀態(tài)（D= 1）.

圖8 不同m值對應(yīng)的損傷變量D隨ε/ε0的變化Fig. 8 Variation of damage variablesD with ε/ε0 for differentm

3 試驗驗證與討論

3.1 試驗結(jié)果分析

利用多功能電液伺服試驗機(jī)對砂巖試樣開展單、三軸壓縮試驗測，試樣直徑為50 mm，高度為100 mm. 試驗過程中采用位移控制，加載速率為0.002 mm/s，圍壓分別取0、10、20、30、40 MPa. 由試驗所得不同圍壓下砂巖的軸向應(yīng)力-應(yīng)變關(guān)系曲線（圖5）可以看出：不同圍壓下砂巖的應(yīng)力-應(yīng)變曲線具有明顯的階段性特征. 表1給出不同圍壓（σ3）下砂巖的彈性模量E、空隙完全閉合時的 εa、σc、εc和R. 可以看出：砂巖的E、σc、 εc和R隨著圍壓的增加均呈增大趨勢，而 εa隨圍壓的增加而減小. 由Morh-Coulomb強(qiáng)度準(zhǔn)則，可得到該砂巖峰值強(qiáng)度對應(yīng)的黏聚力（cc）和內(nèi)摩擦角（φc）分別為12.08 MPa和44.6°；殘余強(qiáng)度對應(yīng)的黏聚力（cr）和內(nèi)摩擦角（φr）分別為0.12 MPa和37.4°，巖石的泊松比+ μ = 0.21.

根據(jù)前述確定模型參數(shù)的方法，基于試驗結(jié)果得到了與空隙應(yīng)變比K相關(guān)的參數(shù)a1和a2，以及與損傷變量D相關(guān)的參數(shù)和m，見表2.

根據(jù)表2中對試驗數(shù)據(jù)進(jìn)行反演得到的模型參數(shù)，從而得到相關(guān)模型參數(shù)與圍壓 σ3

之間的函數(shù)變化關(guān)系，如式（27）. 可知：隨著圍壓 σ3的增大，a1先增加后減小，a2近似線性減小，ε0近似為線性增加，m大致呈指數(shù)形式的減小.

由表1和表2可知：該模型中Weibull分布參數(shù) ε0略 大 于 εc，不 同 圍 壓 下 ε0?εc的 平 均 值 為0.103. 圖9給出了 ε0與 εc之間的關(guān)系曲線. 可以看出：兩者之間具有很好地一致性. 因而為了方便計算可用應(yīng)力-應(yīng)變曲線對應(yīng)的峰值應(yīng)變 εc+ 0.1來代替Weibull分布參數(shù) ε0.

表1 巖石三軸試驗參數(shù)Tab. 1 Triaxial test parameters for rock

表2 本文巖石統(tǒng)計損傷模型參數(shù)Tab. 2 Parameters of statistical damage model for rocks

圖9 Weibull分布參數(shù)ε0與εc的關(guān)系Fig. 9 Relationship between peak strainεc and the parameterε0 of Weibull distribution

3.2 模型驗證

以下采用砂巖試驗數(shù)據(jù)對本文所建巖石的損傷統(tǒng)計本構(gòu)模型進(jìn)行驗證，本文模型計算結(jié)果與試驗數(shù)據(jù)對比情況如圖10所示. 可以看出：本文模型在不同圍壓條件下的計算結(jié)果與試驗曲線均有較高的吻合程度，很好地反映荷載作用下巖石破壞變形的全過程，尤其是初始空隙壓密及屈服后各階段的非線性變形特征，從而驗證了本文所建巖石本構(gòu)模型的合理性.

圖10 本文模型計算結(jié)果與試驗曲線之間的比較Fig. 10 Comparison between the proposed constitutive model calculation values and experimental curves

以圍壓 σ3=30 MPa條件下砂巖的試驗數(shù)據(jù)為例，模型參數(shù)的取值見表2，比較已有模型和本文模型對試驗數(shù)據(jù)的模擬效果，如圖11所示. 可以看出：由于文獻(xiàn)[12]和文獻(xiàn)[20]建立的巖石損傷本構(gòu)模型均未考慮初始空隙壓密對巖石變形的影響，導(dǎo)致模型計算結(jié)果在初始段呈線性彈性變化，與試驗曲線差別較大；由于文獻(xiàn)[12]模型沒有考慮巖石破壞后的殘余強(qiáng)度，導(dǎo)致其殘余強(qiáng)度變形階段與試驗曲線有很大差別. 可見，本文所建模型能夠較好地模擬巖石應(yīng)力-應(yīng)變?nèi)^程5個階段的變形特征.

圖11 不同模型計算結(jié)果與試驗值的比較Fig. 11 Different model calculations versus experimental curve

4 結(jié) 論

1） 將巖石抽象為由實體骨架和內(nèi)部空隙兩部分組成的材料，分析了這兩部分的變形機(jī)理以及與巖石整體變形之間的關(guān)系，提出了空隙應(yīng)變比K的概念，建立了巖石的變形分析模型，奠定了荷載作用下巖石破壞變形全過程模擬方法研究的基礎(chǔ).

2） 利用三軸試驗結(jié)果，推導(dǎo)了空隙應(yīng)變比K的演化方程，結(jié)合巖石變形的分析方法，基于損傷統(tǒng)計理論并考慮巖石破壞后的殘余強(qiáng)度變形，建立了能夠合理描述巖石變形破壞全過程的損傷統(tǒng)計本構(gòu)模型，初始空隙壓密階段及隨后各階段的變形特征均能得到較好的反映，不同圍壓下所建模型與試驗結(jié)果均有較高的吻合度，相比現(xiàn)有模型更為合理.

3） 通過本文模型、現(xiàn)有相關(guān)模型和試驗曲線的比較，驗證了本文所建巖石損傷模型的有效性和合理性.

致謝：中國博士后科學(xué)基金（2018M633627XB）.