考慮初始空隙壓密的巖石變形全過程本構(gòu)模型

李修磊 ，陳洪凱 ，張金浩

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

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

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

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

1 巖石變形力學機理與分析

1.1 巖石變形過程的力學機理

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

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

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

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

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

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

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

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

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

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

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

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

1.3 空隙應變比K的確定

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

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

圖3 應力-應變試驗曲線Fig. 3 Stress-strain test curve

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

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

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

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

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

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

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

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

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

利用式（14）對圖6中的試驗曲線進行擬合，可得到空隙應變比K隨 ε1變化的表達式為

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

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

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

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

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

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

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

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

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

2.2 損傷部分殘余強度的計算

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

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

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

2.3 損傷變量的演化

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

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

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

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

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

1） 當 ε1<εa時

2） 當 ε1≥εa時

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

1） 當 ε1<εa時

2） 當 ε1≥εa時

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

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

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

3 試驗驗證與討論

3.1 試驗結(jié)果分析

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

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

根據(jù)表2中對試驗數(shù)據(jù)進行反演得到的模型參數(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)系曲線. 可以看出：兩者之間具有很好地一致性. 因而為了方便計算可用應力-應變曲線對應的峰值應變 ε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)模型進行驗證，本文模型計算結(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所示. 可以看出：由于文獻[12]和文獻[20]建立的巖石損傷本構(gòu)模型均未考慮初始空隙壓密對巖石變形的影響，導致模型計算結(jié)果在初始段呈線性彈性變化，與試驗曲線差別較大；由于文獻[12]模型沒有考慮巖石破壞后的殘余強度，導致其殘余強度變形階段與試驗曲線有很大差別. 可見，本文所建模型能夠較好地模擬巖石應力-應變?nèi)^程5個階段的變形特征.

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

4 結(jié) 論

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

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

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

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