摘要：為研究撫順西露天礦邊坡中泥巖蠕變特征，利用ITZ-1型三軸剪切儀，采用陳氏分級加載法對撫順西露天礦邊坡中的巖樣進行蠕變試驗，獲得巖樣蠕變曲線。對試驗數(shù)據(jù)進行數(shù)值擬合得到改進的Burgers模型基本參數(shù)，在ABAQUS有限元分析軟件中對改進的Burgers模型二次開發(fā)編寫UMAT子程序。通過Python腳本結合Voronoi多邊形法建立巖石介觀結構模型，模擬得出巖樣蠕變曲線，分析試驗蠕變曲線及數(shù)值模擬得出的介觀尺度損傷。研究結果表明：撫順西露天礦邊坡中泥巖具有減速和等速蠕變階段特性，巖石介觀模型模擬所得蠕變結果與試驗結果誤差較小，泥巖中沿軸向生長的介觀結構更易發(fā)生損傷、破壞。所建巖石介觀模型可以更精確地模擬巖石介觀結構特征，為工程實際中的巖石力學仿真模擬提供參考。

關鍵詞：蠕變試驗；巖石介觀模型；改進Burgers模型；UMAT子程序

中圖分類號：TU458 文獻標志碼：A 文章編號：1008-0562（2024）02-0156-07

0 引言

蠕變特性是巖體重要的工程力學特性之一。大量工程實踐經(jīng)驗與理論研究表明，在大部分工程地質條件下，軟巖邊坡的破壞與失穩(wěn)不是在開挖時立刻發(fā)生的，而是在長期固定應力作用下產(chǎn)生隨時間不斷增長的蠕變應變，最終導致應變超出巖體的承受能力而破壞。巖體蠕變的本質原因是在載荷作用下的微裂紋及結構面發(fā)展。許多學者通過CT掃描、聲發(fā)射等技術觀測巖石介觀損傷現(xiàn)象，由此構建介觀損傷模型理論，并使用離散元（DEM）等軟件進行數(shù)值模擬仿真分析。

目前對巖石蠕變本構方程及其力學模型的研究主要從兩方面進行，一是蠕變試驗，通過試驗測得巖石各部分蠕變曲線特征，基于數(shù)值擬合等方法得出蠕變本構方程中的各項參數(shù)并以此體現(xiàn)巖石的蠕變特性。二是蠕變模型理論研究，通過對虎克彈簧、牛頓黏壺及圣唯南塑性體等基本元件串聯(lián)、并聯(lián)或對元件參數(shù)進行改進來描述各類巖石不同的蠕變特性。蠕變模型理論研究中，一部分從經(jīng)典蠕變模型（Maxwell模型、Burgers模型、Bingham模型等）的改進出發(fā)，對基本元件進行非線性參數(shù)改進或直接添加非線性元件以符合不同材料的蠕變性質；另一部分引入損傷理論、摩爾一庫倫理論等描述損傷等因素對巖石蠕變的影響。

在有限元計算中對蠕變介觀尺度損傷的探討還很少，大部分研究沒有考慮蠕變過程中巖石介觀結構的變化。本文運用三維Voronoi多邊形，使用Python腳本編程建立能夠描述巖石微裂紋、結構面、礦物質組成等介觀結構特征的模型，并通過編寫UMAT子程序，運用改進的Burgers模型對巖樣蠕變特性進行數(shù)值模擬，分析巖石介觀結構在蠕變條件下的損傷特性。

1 巖石分級加載蠕變試驗

1.1 試件制備及試驗方案

巖石的壓縮蠕變是巖體在長時間應力作用下應變逐漸增加的力學性質，采用巖石分級加載蠕變試驗可以充分體現(xiàn)泥巖的蠕變特性。試驗中巖樣取自撫順西露天礦開采邊坡工程現(xiàn)場，為黑色泥巖，通過3組X射線衍射（XRD）試驗得出該泥巖礦物組成及其質量分數(shù)分別為：高嶺石占35.4%～38.8%，微斜長石占21.8%～22.8%，石英占21.7%～26.6%，伊利石占11.3%～20.1%，其中高嶺石和伊利石屬于黏土礦物，這種泥巖質地軟，強度低，且由多種礦物結合而成，易出現(xiàn)結構面和微裂紋等影響強度的介觀結構。對試樣進行切割、磨平，過程中盡量減少對巖樣的破壞，以保持巖石原本的力學性質。為符合儀器尺寸，將取回的巖塊加工成尺寸為Φ58 mm×116 mm的圓柱體，試樣制備按照《水利水電工程巖石試驗規(guī)程》（SL/T 264-2020） 中有關規(guī)定進行。

試驗儀器采用TSZ-1型全自動三軸剪切儀，見圖1。儀器由圍壓控制室、軸向壓力系統(tǒng)（包括量力環(huán)、壓力室、軸向加載平臺）及自動數(shù)據(jù)采集系統(tǒng)組成。剪切蠕變試驗根據(jù)陳氏分級加載法，在前期測得試樣極限載荷后，取極限載荷的75%為最高應力條件。分為1 MPa、2MPa、3MPa、4MPa共4級連續(xù)加載，每級加載持續(xù)10 h左右，獲得分級應力加載條件下的流變曲線。

1.2 試驗結果及分析

在分級應力加載條件下根據(jù)試驗結果繪制出典型撫順西露天礦泥巖的蠕變曲線，見圖2。

由圖2可見，該泥巖試樣隨時間產(chǎn)生了明顯的蠕變現(xiàn)象。由于試驗儀器量程有限，巖樣未出現(xiàn)加速蠕變階段。在每級加載后巖樣首先經(jīng)過應力應變成正比關系的彈性變形階段，之后出現(xiàn)典型的蠕變階段，其中減速蠕變階段應變變化速度迅速下降，當應變變化速度下降到一固定值時，到達等速蠕變階段。分級加載試驗中，巖樣在較低應力水平作用下減速蠕變階段持續(xù)時間較短，等速蠕變階段蠕變速度較小，蠕變變形量較小，在1 MPa應力水平下經(jīng)過10 h后，巖樣軸向應變增長了2.19×10-4。隨著應力水平逐漸升高，蠕變變形量有逐漸增大的趨勢，減速蠕變階段持續(xù)時間增長。當應力增長至4 MPa時，巖樣軸向應變增長速度明顯加快，經(jīng)過10 h后應變增長了1.82×10-3。

巖石試樣的宏觀蠕變特性是其內部微損傷在等載荷作用下不斷演化的結果，表現(xiàn)為微裂紋的萌生、擴展和合并隨時間變化。等速蠕變階段的微裂紋數(shù)量較少并且主要集中在邊緣部分，之后微裂紋在巖樣全部范圍內快速增加，導致巖樣強度持續(xù)下降。對巖石介觀結構建模能夠考慮這種損傷，使巖石蠕變模擬更加精確，并直觀表現(xiàn)蠕變過程中這種介觀尺度的變化。

2 改進的Burgers模型

2.1 改進的Burgers模型

采用付凱敏等改進的Burgers模型進行泥巖蠕變模擬。改進的Burgers模型是采用元件模型法，將經(jīng)典Burgers模型中的外置黏壺元件改進為非線性元件所得，其模型示意見圖3。

改進的Burgers模型包含5個參數(shù)，其中，E1、E2為2個線性虎克彈簧元件的彈性模量，MPa；改進后外置黏壺的黏度系數(shù)為η1=AeBt，A、B為改進后黏壺的獨立系數(shù)；t為時間，s；η2為內置黏壺的黏度系數(shù)，MPa·h。由基本元件方程及Laplace變換可以得出改進的Burgers蠕變本構方程為

式中：ε為應變；σ為應力，MPa。

2.2 模型參數(shù)辨識

利用Matlab軟件，使用改進的Burgers模型本構方程對試驗中巖樣分別在1 MPa、2 MPa、3 MPa、4MPa應力條件下得出的蠕變曲線進行數(shù)值擬合，擬合曲線見圖4。所有應力情況下擬合可決系數(shù)均在0.97以上，可見數(shù)值擬合效果較好。

由蠕變各階段分別擬合出所需元件的參數(shù)，最終所得各模型參數(shù)見表1。擬合所得參數(shù)具有一定規(guī)律性。由表1可見，泥巖巖樣彈性模量E1為452 MPa左右，η2隨應力增大而增大，巖樣的黏彈性增強，體現(xiàn)出在分級加載過程中蠕變的疊加性質。

3 巖石細觀模型蠕變數(shù)值模擬

3.1 UMAT子程序

使用有限元軟件ABAQUS進行泥巖蠕變模擬，并通過UMAT子程序賦予其材料屬性。UMAT子程序工作時在每個增量步中與ABAQUS主程序都進行數(shù)據(jù)交互，其中主程序將模型的初始應力和應變傳入子程序，再由子程序對傳入的數(shù)據(jù)進行處理，得出主程序需要的下一增量步中的應力和應變，并在模型每個單元中的積分點上進行計算，子程序與主程序協(xié)同工作過程見圖5。

根據(jù)ABAQUS中對用戶材料子程序的要求，需要將改進Burgers模型蠕變方程由一維推廣至三維，且采用增量形式來實現(xiàn)材料積分單元的迭代求解。在三維應力張量中，球張量僅對應巖體受靜水壓力作用下的應力狀態(tài)，而偏張量代表巖體蠕變效應，因此蠕變本構可借鑒塑性理論中的分析方法。將塑性力學中的米塞思等效應力定律代入蠕變本構方程中，得到新的蠕變方程為

式中：a0、a1、a2為材料常數(shù)；ε為等效應變；σ為等效應力，MPa。

對于蠕變應力與應變之間的關系，假設塑性流動定律仍成立，即

式中：f為類似于塑性勢平面的函數(shù)；Sij為偏應力張量；λc為塑性流動因子。

由式（2）～式（4）推出由米塞思塑性力學公式演化而來的蠕變本構方程為

蠕變產(chǎn)生的應變隨時間變化，因此其應力應變本構矩陣也與時間相關。采用初應變法，由蠕變應變增量代替增量步開始時的初應變，則應力增量可表示為

式中，Dijkl為雅克比矩陣。

3.2 建立巖石細觀模型

Voronoi多邊形是一種在一定空間內隨機構建多個多邊形的方法，其原理是在空間內生成隨機分布的離散點，再由離散點自動組成有序的三角網(wǎng)絡并編號，對每個三角形生成外接圓圓心，并按編號連接圓心即可構成多邊形。

由平面的Voronoi多邊形可進行三維推廣初步形成隨機的多塊體堆積模型，其結構符合巖石體內各種礦物質組成及結構面的特性。由于這種結構較為復雜，劃分普通網(wǎng)格后采用ABAQUS有限元計算難以收斂，因此需要對模型進行三維網(wǎng)格化處理，模型建立過程見圖6。運用Python語言編寫腳本并在ABAQUS中運行，最終創(chuàng)建出多個隨機多邊形顆粒組成的巖石細觀模型，見圖6 （c）。

建立巖石細觀模型后，首先對模型進行網(wǎng)格劃分，再通過Python腳本識別出每個顆粒之間連接的單元面，在這些單元面之間插入零厚度的Cohesive單元即可描述每個顆粒之間的黏結、摩擦和分離等力學性質，生成的Cohesive單元見圖7。

三維零厚度的Cohesive單元本構關系符合牽引分離準則，見圖8，其損傷起始階段準則采用二階損傷準則，即

式中：（·）為Macaulay括號，用以確保Cohesive單元不在純壓狀態(tài)下破壞；σxo、σyo、σzo為單元變形完成后3個維度的峰值名義應力，MPa，其中x、y、z為坐標軸的3個維度。

圖8中δo表示相對于損傷開始時的有效位移，m；δf為完全失效時的有效位移，m；G為能量，N/m。

當單元中名義應力滿足式（7）時，損傷演化開始，在三維應力狀態(tài)下混合模態(tài)損傷響應符合Benzeggagh-Kenane斷裂準則，即

式中：GD為因失效耗散的能量，N/m；Gy-Gy+Gz，GZ=Gx+Gz;η為材料參數(shù)。參考文獻[21]，GDx取40N/m，GDy取400 N/m，η取1.5。

損傷變量D通過線性損傷演化表示為

3.3 巖石蠕變數(shù)值模擬

為驗證改進的Burgers模型、UMAT子程序及巖石細觀模型的正確性，將使用ABAQUS計算得到的巖樣蠕變曲線與試驗所得曲線進行對比。

通過巖石細觀模型腳本建立Φ58mm×116 mm的圓柱體模型并劃分為3 650個六面體單元，單元類型為C3D8，由UMAT子程序確定巖石材料屬性，插入的Cohesive單元類型為COH3D8。模型底端固定，在頂面施加軸向的載荷。在彈性分析之后，使用VISCO黏性分析步得出泥巖蠕變10 h后軸向應變數(shù)據(jù)，繪制蠕變曲線，見圖9。

由圖9可以看出，在各級應力狀態(tài)下模擬所得的泥巖蠕變曲線與試驗所得的曲線吻合程度較高。數(shù)值模擬所得蠕變曲線更光滑，且蠕變應變較試驗結果稍低，誤差小于2.1%。試驗結果會受到不同程度的人為和外界環(huán)境影響，而數(shù)值模擬條件為理想條件，因此二者存在一定差距。

在4 MPa應力水平下介觀模型顆粒間Cohesive單元損傷情況見圖10，可見Cohesive單元受到不同程度損傷，且在單軸壓縮狀態(tài)下沿軸向的單元比沿橫向單元損傷程度更高，傾向于發(fā)生剪切破壞，即沿軸向生長的微裂紋較橫向生長的微裂紋更易發(fā)生損傷，巖石蠕變后出現(xiàn)軸向裂紋而發(fā)生破壞，這也符合大部分巖石蠕變試驗的規(guī)律。巖石介觀模型精細化地模擬了巖體中的微裂紋等介觀結構，反映出的巖石損傷特性及蠕變等與試驗結果相符，說明模型能夠較好地模擬工程實際中巖石的應力集中和介觀損傷。

4 結論

通過對撫順西露天礦泥巖進行單軸蠕變試驗得到巖樣蠕變特性，并通過Python腳本程序及UMAT子程序驅動ABAQUS有限元模擬軟件進行巖石介觀模型數(shù)值模擬，分析巖石蠕變介觀結構損傷特性，得出如下結論。

（1）撫順西露天礦邊坡中的泥巖具有典型的減速和等速蠕變階段，隨著應力狀態(tài)的改變，由于分級加載蠕變變形的疊加效應，其黏性系數(shù)增大，蠕變特征更加明顯。

（2）通過編寫UMAT子程序將改進的Burgers模型運用于對巖石蠕變的數(shù)值模擬中，得出泥巖的蠕變曲線，所得蠕變應變與試驗誤差不超過2.1%。

（3）通過構建的巖石介觀模型數(shù)值模擬發(fā)現(xiàn)，Cohesive單元損傷集中在沿軸向單元上，因此巖石在單軸蠕變條件下沿軸向生長的微裂紋及結構面更易發(fā)生破壞。

基金項目：國家自然科學基金項目（51808351）;遼寧省教育廳基礎研究項目（lnjc202019）