滲透壓作用下黃山石灰石礦區炭質泥頁巖剪切蠕變特性分析

嚴 珊，胡 斌，魏二劍，馬利遙，田 磊

（1.武漢科技大學資源與環境工程學院，湖北 武漢 430081；2.冶金礦產資源高效利用與造塊湖北省重點實驗室，湖北 武漢 430081）

我國四川省黃山石灰石礦區大范圍存在含緩傾軟弱夾層的二疊系石灰巖地層[1]，經過調查發現，這類軟弱夾層主要是含炭質的薄層泥頁巖，由于我國西南地區雨水豐沛，軟弱夾層在雨水浸泡、疏干的長期循環作用下，其蠕變特性將發生顯著變化，軟弱夾層極易沿著破壞面發生壓剪破壞，從而導致礦山邊坡滑坡。

目前，國內外學者在巖石蠕變力學相關領域已經進行了大量研究：秦哲等[2]對經歷不同飽水-失水循環次數的巖樣進行下三軸蠕變試驗，從細觀角度分析了水巖作用對巖石的損傷機理；王瑋瑋等[3]對炭質泥巖開展了三軸壓縮蠕變試驗，探討了不同干濕循環次數下炭質泥巖孔隙的演化特征和蠕變規律；FUJII 等[4]開展了砂巖在飽水狀態與干燥狀態下的蠕變試驗，發現飽水狀態下，砂巖的環向應變更大；李鵬等[5]為探究含水率對軟弱夾層剪切蠕變特性的影響，開展了軟弱夾層在多種含水率下的剪切蠕變試驗；于永江等[6]通過開展不同含水率軟巖的直剪蠕變試驗，建立了能夠描述含水軟巖蠕變特性的非線性剪切蠕變模型。上述研究表明，目前對于水巖作用下巖石蠕變特性方面的試驗研究主要集中在試驗前對巖樣進行飽水-干燥循環作用以及巖樣的含水率方面，這與實際情況中巖體主要處于壓縮狀態下持續產生變形的同時，仍不斷經歷雨水浸泡、疏干的反復循環作用過程存在一定差異。

因此，以黃山石灰石礦區炭質泥頁巖為研究對象，開展滲透壓作用下的剪切蠕變試驗，為了研究其蠕變特性，基于傳統的西原模型，提出了一種能同時描述巖石黏彈塑性特性的非線性蠕變模型，并利用該模型進行剪切蠕變參數辨識。以此為基礎，運用大型有限元軟件FLAC3D對炭質泥頁巖的室內剪切蠕變試驗進行數值模擬，通過數值模擬方法獲得了室內剪切蠕變試驗中炭質泥頁巖的變形特征，這些結果拓展了試驗成果和實測資料，使之可以更合理地評價炭質泥頁巖在滲透壓作用下的剪切蠕變特性，也為眾多含二疊系軟弱夾層的露天礦山邊坡長期穩定性分析提供一定的理論支撐。

1 礦區工程地質條件

黃山石灰石礦山位于四川盆地邊緣低山～低中山地段。礦區地勢南高北低，最高標高為1 229 m，最低標高為640 m，相對標高589 m，屬低中山區。礦區氣候以四川盆地亞熱帶濕潤性季風氣候為主，全年陰濕多雨，雨季多集中在6 月到9 月。在大量降雨和周邊地震的影響下，雨水入滲流入邊坡內部軟弱夾層，造成軟弱夾層力學參數不斷劣化，內部節理裂隙貫通，加之軟弱夾層明顯的流變效應，極易引起滑坡，如圖1 所示的東采區滑坡、老鷹嘴滑坡等，危及礦山生產安全及礦山山腳下村民的生命安全。

圖1 黃山石灰石礦區數字化全貌圖Fig.1 Digital panorama of Huangshan limestone mining area

2 滲透壓作用下炭質泥頁巖剪切蠕變試驗

2.1 試驗概況

室內試驗的炭質泥頁巖試樣采自黃山石灰石礦區，室內試驗使用課題組研發的降雨滲流-爆破振動耦合模擬軟巖剪切流變試驗系統，如圖2 所示。將所取巖塊送至加工廠加工成150 mm×75 mm×75 mm的長方體試樣，試樣包含三個等孔距的滲流孔，孔距為37.5 mm，孔深度為37.5 mm，制備完成的試樣如圖3 所示。

圖2 軟巖剪切流變試驗系統Fig.2 Shear rheological test system of soft rock

圖3 軟巖試樣Fig.3 Soft rock specimen

2.2 試驗方案及結果分析

本次室內試驗是對四組炭質泥頁巖（J-1、J-2、J-3、J-4）進行直剪試驗，試驗結果見表1。再對另外四組炭質泥頁巖（J-5、J-6、J-7、J-8）進行不同滲透水/氣壓循環次數下的剪切蠕變試驗，分別為循環0 次、4 次、8 次、12 次。對試樣施加1 MPa 法向壓力和五級剪切應力的過程中交替施加0.5 MPa 滲透水/氣壓。以施加0.5 MPa 滲透水壓4 h 后再施加0.5 MPa 滲透氣壓4 h 為一次循環。需要說明的是，本次室內試驗設計的滲透水/氣壓的循環在施加第四級剪切應力前結束，在施加法向壓力和剪切應力的同時進行滲透水/氣壓的循環是本次試驗的一個創新點，更好地還原了實際情況中軟弱夾層主要處于壓縮狀態下持續產生變形的同時，仍不斷經歷雨水浸泡、疏干的反復循環作用過程。剪切蠕變試驗方案設計見表2 和圖4。剪切蠕變試驗法向壓力和滲透水/氣壓都為恒定值。

表1 軟巖瞬時強度參數Table 1 Instantaneous intensity parameters of soft rock單位：MPa

表2 剪切蠕變試驗方案設計Table 2 Shear creep test protocol design

圖4 剪切蠕變試驗方案設計Fig.4 Protocol design of shear creep test

圖5 為不同循環次數下的剪切蠕變曲線。由圖5 可知，在前四級剪切應力作用下，炭質泥頁巖經歷了兩個階段，即減速蠕變階段和穩定蠕變階段，在第五級剪切應力作用下，炭質泥頁巖則經歷了完整的蠕變三階段。在減速蠕變階段中，炭質泥頁巖的剪切應變剛開始隨時間增大較快，但其應變率隨時間逐漸減小。在穩定蠕變階段中，曲線近似成直線，炭質泥頁巖的剪切應變隨時間勻速增加。在加速蠕變階段中，炭質泥頁巖的剪切應變隨時間的增長而加速增加，其應變率也隨時間迅速增加，直至炭質泥頁巖發生破壞。此外，炭質泥頁巖巖樣在每一級剪切應力加載的瞬間都會產生瞬時應變，且隨著滲透水/氣壓循環次數的增加，瞬時應變也隨之增加。在前四級剪切應力作用下，其剪切蠕變變形都非常小。以滲透水/氣壓循環0 次為例，在施加第一級剪切應力1.78 MPa 時，炭質泥頁巖巖樣達到穩定蠕變階段時發生的剪切位移為0.021 mm，在施加第二級剪切應力2.16 MPa 時，炭質泥頁巖巖樣達到穩定蠕變階段時發生的剪切位移為0.025 mm，在施加第三級剪切應力2.54 MPa 時，炭質泥頁巖巖樣達到穩定蠕變階段時發生的剪切位移為0.029 mm，在施加第四級剪切應力2.91 MPa 時，炭質泥頁巖巖樣達到穩定蠕變階段時發生的剪切位移為0.032 mm。

圖5 不同循環次數下的剪切蠕變曲線Fig.5 Shear creep curves at different cycle times

2.3 炭質泥頁巖的長期強度

巖體的長期強度是評價工程長期穩定性的重要參量[7]，目前應用最廣泛的是穩態流變速率法，該方法具有計算簡單、直觀方便的優點。根據剪切蠕變試驗曲線計算不同循環次數下、各級剪切應力下炭質泥頁巖巖樣的穩態蠕變速率，并采用指數函數擬合得到巖樣的穩態蠕變速率-剪切應力關系曲線，如圖6 所示。根據該方法確定了在滲透壓作用下炭質泥頁巖剪切蠕變的長期強度，炭質泥頁巖的長期強度即為由穩定蠕變過渡到加速蠕變的極限應力。炭質泥頁巖的蠕變破壞強度即炭質泥頁巖發生剪切破壞時的強度，由試驗即可確定，結果見表3。

表3 軟巖長期強度與蠕變破壞強度參數Table 3 Long-term strength and creep failure strength parameters of soft rock

圖6 蠕變速率與剪切應力關系曲線Fig.6 Relationship curve between creep rate and shear stress

3 改進的西原模型及參數辨識

3.1 改進的西原模型

傳統的西原模型由虎克體、黏彈性體以及黏塑性體串并聯而成[8]，如圖7 所示。由于傳統的西原模型難以描述炭質泥頁巖巖樣非線性加速蠕變階段的蠕變特性，在傳統西原模型的基礎上，提出了一種能同時描述巖石黏彈塑性的非線性蠕變模型，該模型在西原模型的基礎上串聯了一個非線性黏塑性體，并引入了兩個改進的非線性牛頓體黏壺，如圖8 所示。E1為瞬時彈性模量，E2、E3均為黏彈性模量，η1、η2、η3均為黏滯系數，ε1、ε2、ε3、ε4為圖8 中各部分對應的應變，當剪切應力小于長期強度τs時，該非線性黏塑性體不發揮作用，當剪切應力大于長期強度τs時，巖樣進入加速蠕變階段，該非線性黏塑性體發揮作用[9]。

圖7 西原模型Fig.7 Nishihara model

圖8 改進的西原模型Fig.8 Improved Nishihara model

夏才初等[10]研究指出直接將蠕變參數中的黏滯系數η替換為η（t）得到的結果將是錯誤的蠕變方程。因此，本文參照相關文獻[11-15]對黏滯系數η2、η3進行等效的非線性化處理，針對圖8 中ε3部分和ε4部分非線性黏塑性體構造相應的本構方程，見式（1）和式（2）。

式中：τ0、τs分別為剪切應力和長期強度；ts為巖樣由穩定蠕變階段轉為加速蠕變階段對應的時刻；n為蠕變指數；Δt為單位時間，主要用來保持量綱一致，取1 h。故改進的西原模型一維蠕變方程見式（3）和式（4）。

3.2 剪切蠕變模型參數辨識

根據炭質泥頁巖在滲透壓作用下的剪切蠕變試驗曲線，識別改進的西原剪切蠕變模型相應的參數，其參數識別結果見表4。由表4 可知，模型表達式與試驗數據擬合效果很好，相關性系數高，R2均在0.96以上，如圖9 所示。由此可知，改進的西原剪切蠕變模型能夠很好地模擬炭質泥頁巖完整的蠕變三階段。因此，所建立的改進的西原剪切蠕變模型是合理的。

表4 改進的西原剪切蠕變模型擬合參數Table 4 Fitting parameters of improved Nishihara shear creep model

圖9 試驗曲線與擬合曲線對比圖Fig.9 Comparison of test curve and fitting curve

采用指數函數擬合得到炭質泥頁巖的彈性模量-滲透水/氣壓循環次數關系曲線與黏滯系數-滲透水/氣壓循環次數關系曲線如圖10 所示。受滲透水/氣壓循環次數的影響，炭質泥頁巖的平均瞬時彈性模量、平均彈性模量、平均黏滯系數均不斷減小，進一步說明了降雨反復作用影響下炭質泥頁巖力學強度降低，水理性質變差，從而影響礦山邊坡穩定性。

圖10 模型參數變化規律Fig.10 Variation pattern of model parameters

4 數值模擬分析

4.1 數值計算模型

將自定義的蠕變本構模型進行二次開發，數值計算模型的形狀采用長方體，尺寸與炭質泥頁巖試樣的尺寸相同，長寬高分別為150 mm×75 mm×75 mm，共8 803 個單元，3 385 個節點。采用的邊界條件為試樣底部約束法向位移，試樣上半部分的左面約束橫向位移。采用的加載條件與室內試驗中滲透水/氣壓循環4 次下的加載條件相同，試樣上表面施加1 MPa 的法向應力，試樣下半部分的右面施加剪切應力，剪切應力為1.69～3.20 MPa，應力梯度為0.38 MPa，分五級加載，每級剪切應力加載時間為24 h，共120 h。采用自定義的蠕變本構模型計算，計算參數見表4。

4.2 數值模擬結果分析

圖11 為數值模擬計算結果圖。由圖11 可知，最大位移值位于試樣的右下角部分，由于試樣上半部分的左面施加了橫向約束，所以試樣的左上角部分位移值為0，且隨著剪切應力的不斷增加，X軸方向的位移也在不斷增大。在第一級剪切應力作用下，最大位移值為0.18 mm；在第二級剪切應力作用下，最大位移值為0.32 mm；在第三級剪切應力作用下，最大位移值為0.47 mm；在第四級剪切應力作用下，最大位移值為0.62 mm；在第五級剪切應力作用下，最大位移值為1.17 mm。X軸方向的位移云圖的變化趨勢與實際情況也比較相符，圖11（b）可以看到明顯的剪切變形。

圖11 數值模擬計算結果Fig.11 Numerical simulation calculation results

對試樣上（0.15，0，0.037 5）點進行位移監測，得到如圖12（a）所示的位移監測圖。由于在試驗過程中采用分級加載方式，剪切蠕變呈臺階型，需利用Boltzmann 迭加原理對試驗數據進行二次處理，處理后的剪切蠕變曲線如圖12（b）所示。由此可見，使用FLAC3D數值模擬軟件得到的剪切蠕變曲線與室內試驗剪切蠕變曲線總體趨勢相同，從而驗證了室內試驗和所建立的非線性黏彈塑性本構模型的正確性。

圖12 室內試驗結果與數值模擬結果對比Fig.12 Comparison of indoor test results and numerical simulation results

5 結論

1）在滲透水/氣壓的循環作用下，隨著剪切應力的增加，炭質泥頁巖試樣的剪切蠕變曲線呈現了完整的蠕變三階段。不同的滲透水/氣壓循環次數下，試樣的含水率不同，剪切蠕變曲線也會有所區別。

2）基于室內滲透壓作用下炭質泥頁巖剪切蠕變試驗得到的剪切蠕變曲線，在西原模型的基礎上，串聯了一個非線性黏塑性體，并引入了兩個改進的非線性牛頓體黏壺，建立了改進的西原剪切蠕變本構模型，能夠很好地擬合炭質泥頁巖試樣蠕變的三個階段，具有較好的實用性。

3）運用FLAC3D軟件對炭質泥頁巖試樣以相同的室內試驗條件進行數值模擬計算，結果表明：最大位移值位于炭質泥頁巖試樣的右下角部分，可以看到明顯的剪切變形。使用FLAC3D數值模擬軟件得到的X軸方向的位移監測曲線與室內試驗剪切蠕變曲線趨勢相同，都呈臺階狀且包含完整的蠕變三階段，與實際結果比較吻合。