撫順西露天礦泥巖蠕變特性試驗(yàn)及介觀模型研究

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

關(guān)鍵詞：蠕變試驗(yàn)；巖石介觀模型；改進(jìn)Burgers模型；UMAT子程序

中圖分類號：TU458 文獻(xiàn)標(biāo)志碼：A 文章編號：1008-0562（2024）02-0156-07

0 引言

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

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

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

1 巖石分級加載蠕變試驗(yàn)

1.1 試件制備及試驗(yàn)方案

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

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

1.2 試驗(yàn)結(jié)果及分析

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

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

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

2 改進(jìn)的Burgers模型

2.1 改進(jìn)的Burgers模型

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

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

式中：ε為應(yīng)變；σ為應(yīng)力，MPa。

2.2 模型參數(shù)辨識

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

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

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

3.1 UMAT子程序

使用有限元軟件ABAQUS進(jìn)行泥巖蠕變模擬，并通過UMAT子程序賦予其材料屬性。UMAT子程序工作時在每個增量步中與ABAQUS主程序都進(jìn)行數(shù)據(jù)交互，其中主程序?qū)⒛Ｐ偷某跏紤?yīng)力和應(yīng)變傳入子程序，再由子程序?qū)魅氲臄?shù)據(jù)進(jìn)行處理，得出主程序需要的下一增量步中的應(yīng)力和應(yīng)變，并在模型每個單元中的積分點(diǎn)上進(jìn)行計(jì)算，子程序與主程序協(xié)同工作過程見圖5。

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

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

對于蠕變應(yīng)力與應(yīng)變之間的關(guān)系，假設(shè)塑性流動定律仍成立，即

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

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

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

式中，Dijkl為雅克比矩陣。

3.2 建立巖石細(xì)觀模型

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

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

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

三維零厚度的Cohesive單元本構(gòu)關(guān)系符合牽引分離準(zhǔn)則，見圖8，其損傷起始階段準(zhǔn)則采用二階損傷準(zhǔn)則，即

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

（2）通過編寫UMAT子程序?qū)⒏倪M(jìn)的Burgers模型運(yùn)用于對巖石蠕變的數(shù)值模擬中，得出泥巖的蠕變曲線，所得蠕變應(yīng)變與試驗(yàn)誤差不超過2.1%。

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

基金項(xiàng)目：國家自然科學(xué)基金項(xiàng)目（51808351）;遼寧省教育廳基礎(chǔ)研究項(xiàng)目（lnjc202019）