尾砂膠結充填體蠕變模型及在FLAC3D二次開發(fā)中的實驗研究

趙 奎，何 文，熊良宵，楊 欣，王曉軍

（1. 江西理工大學 工程研究院，江西 贛州 341000；2. 江西理工大學 鎢資源高效開發(fā)及應用技術教育部工程研究中心，江西 贛州 341000；3. 寧波大學 建筑工程與環(huán)境學院，浙江 寧波 315211；4. 重慶大學 土木工程學院，重慶 400030；

5. 江西理工大學 資源與環(huán)境工程學院，江西 贛州 341000）

1 前 言

由于在控制地壓、提高資源回收率以及促進礦山環(huán)保等諸多方面具有顯著的優(yōu)越性，充填采礦方法的應用日益廣泛[1]。近年來，研究人員在膠結充填體強度與充填材料配比、充填體力學性質、充填體穩(wěn)定性分析等方面進行了大量的研究，如通過正交試驗研究影響尾砂充填體強度的因素[2]；進行各種不同應力狀態(tài)下充填體破壞能耗試驗研究[3]；通過試驗建立充填體損傷模型[4]；研究布筋尾砂膠結充填體頂板力學性狀[5]，但在膠結充填體蠕變的研究方面鮮見報道。因此，通過試驗建立充填體蠕變模型，不僅具有一定的理論價值，對工程實際中充填體穩(wěn)定性計算分析也是十分必要的。

關于巖石和土的蠕變模型方面已有較多研究成果[6－7]，就巖石蠕變模型而言，材料的力學參數(shù)通常是隨時間等因素而變化的[6]。對于非定常參數(shù)蠕變模型，文獻[8]以西原模型為基礎，從試驗數(shù)據(jù)出發(fā)，運用非線性最小二乘法求得了不同應力水平下和不同時刻的蠕變參數(shù)；文獻[9]通過引入?yún)?shù)初始黏滯系數(shù)、參考時間、參考應力，建立了非定常的修正Burgers模型；文獻[10]通過對現(xiàn)有的非線性理論流變力學模型系統(tǒng)的總結，歸納出處理參數(shù)非線性的3種方式。

本文在巖石蠕變模型研究成果的基礎上，根據(jù)充填體蠕變試驗的結果，運用粒子群優(yōu)化算法，對蠕變模型的蠕變參數(shù)進行了辨識，研究了各蠕變參數(shù)對應力水平的敏感程度，提出并建立了充填體蠕變模型，應用 FLAC3D軟件二次開發(fā)了充填體蠕變模型，對開發(fā)的模型進行了蠕變試驗的數(shù)值模擬，模擬結果與試驗結果對比表明建立的蠕變模型是可行的，可為工程實際中充填體礦柱等蠕變計算提供參考依據(jù)。

2 非定常參數(shù)蠕變方程

2.1 蠕變試驗

試驗制備了配比為 1∶20（水泥與分級尾砂重量比）的充填體試件，水泥為32.5級普通硅酸鹽水泥，分級尾砂取自武山銅礦現(xiàn)場充填用尾砂，試件為圓柱形，尺寸（直徑×高）為φ 50 mm×100 mm，如圖1所示。在試件高度方向的中央粘貼垂直和水平方向的型號為 BX120－6AA應變片，用以測量充填體試件應變，應變儀為BZ2205型程控靜態(tài)電阻應變儀。

圖1 試件及應變片F(xiàn)ig.1 Specimen and strain gauges

蠕變試驗的加載方式通常有單級加載、分級增量加載，由于地下開采中充填體如礦柱等，常為單軸壓縮狀態(tài)，這里采用單軸分級增量循環(huán)加卸載方式。這種方式吸取了分級增量加載方式的優(yōu)點，并克服其缺點，能全面反映材料蠕變的加卸載過程[11]。試驗前測得試件的單軸抗壓強度約1.5 MPa，取強度的80%即1.2 MPa作為最后一級加載值，當位移增量小于0.001 mm/h時，施加下一級荷載。由于充填體材料強度較低，試驗采用自制的簡易分級增量循環(huán)加卸載裝置如圖2所示，試件所受實際壓力由壓力計讀出。

圖2 分級增量加卸載蠕變試驗現(xiàn)場Fig.2 Creep test under step loading and unloading

2.2 蠕變模型

試驗結果如圖3所示，充填體試件有顯著的瞬彈變形，而且變形隨時間的變化基本符合指數(shù)形式規(guī)律，可以采用Hoek-Kelvin（H-K）體的三元件模型（見圖4），分析不同應力水平下的蠕變試驗結果[12－13]。

圖3 不同應力水平下試件蠕變曲線Fig.3 Creep curves of shale under different stress levels

圖4 H-K蠕變模型示意圖Fig.4 Sketch of H-K creep model

H-K體的三元件模型的蠕變方程為

擬合的數(shù)值A = 11.408 GPa，B = 0.00359 kPa-1，C = 1.025，相關系數(shù)R = 0.9998。EH和η變化較小，可直接求平均確定它們的取值，EH= 3.871 GPa，η= 41.102 GPa·h。

式中：σ為分級加載的應力，為一常量；ε為與時間有關的應變；EH、EK分別為H體和K體的彈性參數(shù)；η為材料的黏滯系數(shù)；t為蠕變時間。

2.3 蠕變方程的建立

粒子群算法是一種智能優(yōu)化方法，具有簡單易實現(xiàn)且收斂快、易與其他算法結合等優(yōu)點，可用于求解大量非線性、不可微和多峰值的復雜優(yōu)化問題，廣泛應用于巖土工程[14－17]，尤其是在巖石蠕變參數(shù)識別中廣泛使用[18－19]。

應用粒子群算法進行蠕變方程參數(shù)反演的關鍵是建立適應度函數(shù)，根據(jù)蠕變公式（1）和分級加卸載蠕變試驗數(shù)據(jù)，建立適應度計算公式：

式中：N為σ應力水平下應變監(jiān)測次數(shù)；εj為σ應力水平下第 j次應變監(jiān)測值；EH、EK、η分別為待反演參數(shù)，其余符號意義同前。

采用粒子群優(yōu)化算法參數(shù)反演的結果見表 1。從表中可以看出，運用H-K模型擬合的充填體蠕變方程相關系數(shù)均在0.97以上，表明該模型擬合充填體蠕變方程式是合適的；在不同分級應力水平加載下，參數(shù)EH、η變化較小，表明對應力水平變化不敏感；參數(shù)EK隨應力水平變化較大，表明該參數(shù)對應力水平變化較敏感。

表1 粒子群優(yōu)化算法反演結果Table1 Particle swarm optimization algorithm inversion results

為了更好地表征蠕變規(guī)律和應用于工程實際，參考文獻[10，15]采用式（3）對表1中反演的參數(shù)EK進行回歸擬合：

式中：A、B、C分別為待擬合常數(shù)。

3 蠕變模型在FLAC3D軟件中的實現(xiàn)

3.1 二次開發(fā)環(huán)境

FLAC3D軟件自定義本構模型開發(fā)必須使用C++語言，且編譯成DLL文件（動態(tài)鏈接庫），動態(tài)鏈接庫文件通過model load命令加載。

3.2 廣義Kelvin模型差分形式

二次開發(fā)過程中，需將式（1）一維蠕變本構方程擴展成為三維差分形式，文獻[18]給出了具體推導過程，本文略去推導，直接列出廣義Kelvin模型最終的三維差分增量形式方程：

式中：σij為應力張量；σm為球應力張量；Sij為偏應力張量；δij為kronecker符號；為新的球應力張量；為舊的球應力張量；K為體積模量；Δεkk為應變增量的第一不變量；為新的偏應變張量；為舊的偏應變張量；Δt為每一迭代步的蠕變時間增量；為新的偏應力張量；為舊的偏應力張量；為偏應變張量的增量；Gk為Kelvin體剪切模量；Ge為Hoek體剪切模量；η為黏滯系數(shù)。

3.3 充填體模型

通過對室內(nèi)充填體單軸蠕變試驗的研究，得到充填體的蠕變方程，在FLAC3D軟件中實現(xiàn)式（5）。

本文通過引入應力強度概念來實現(xiàn)軟件自動識別加載應力水平，定義應力強度為

式中：σ1、σ2、σ3為第一、第二、第三主應力。

由經(jīng)典彈塑性力學可得式（7）計算應力強度。

FLAC3D軟件中通過 ps->stnS指針讀取應力張量的各分量，根據(jù)式（7）則可求出應力強度σi。由應力強度的定義可知，單軸試驗時 σi=σ1。得到應力強度σi后，則可根據(jù)式（5）計算EK。因此，本文在H-K蠕變模型基礎上，通過引入應力強度概念，實現(xiàn)了軟件自動的讀取應力加載水平，即在FLAC3D軟件中實現(xiàn)了蠕變參數(shù)隨應力水平變化而變化的充填體蠕變模型。

3.4 開發(fā)模型的驗證

需要注意的是蠕變參數(shù)的換算問題。式（8）分別是H-K模型的一維本構和三維張量形式本構：

為了使三維張量模型和一維模型相對應，將三維張量模型簡化為

將式（8）中一維本構和式（9）中三維本構比較，并結合彈性模量和體積模量與剪切模量定義可得

運用 FLAC3D軟件模擬單軸蠕變時可參考式（10）轉換參數(shù)，即可實現(xiàn)軟件中三維本構和一維理論公式的對應。

進行單軸分級加載蠕變試驗模擬的步驟如下：

①通過模型調用命令，導入生成的動態(tài)鏈接庫文件來導入蠕變模型。

②設置蠕變分析的參數(shù)并調整蠕變時間步長，時間步長的設定必須滿足差分收斂的步長要求。

③施加分級載荷和設置約束條件，每一分級加載需先進行彈性部分計算以得到瞬時變形，再進行蠕變計算。

數(shù)值模擬計算結果和試驗結果如圖5所示。圖中，計算結果與試驗結果吻合，其中誤差的主要是試驗時彈性模量的差異性，而數(shù)值模擬時彈性模量取的是均值。

圖5 數(shù)值模擬與試驗結果對比Fig.5 Comparison between the model’s calculation and test results

4 結 語

通過運用粒子群優(yōu)化算法，對室內(nèi)充填體蠕變試驗數(shù)據(jù)進行了蠕變參數(shù)的辨識，發(fā)現(xiàn)在不同應力狀態(tài)下蠕變參數(shù)EH、η變化較小，蠕變參數(shù)EK對蠕變應力水平變化十分敏感。提出了充填體蠕變模型，并確定了蠕變模型中的參數(shù)，利用FLAC3D軟件二次開發(fā)了所建立的蠕變模型，用開發(fā)的模型進行了分級加載蠕變的數(shù)值模擬，計算結果與試驗數(shù)據(jù)吻合。

[1]周愛民. 礦山廢料膠結充填[M]. 北京: 冶金工業(yè)出版社,2007.

[2]李一帆,張建明,鄧飛,等. 深部采空區(qū)尾砂膠結充填體強度特性試驗研究[J]. 巖土力學,2005,26(6): 865－868.LI Yi-fan,ZHANG Jian-ming,DENG Fei,et al.Experiment study on strength characteristics of tailing cement backfilling at deep-seated mined-out area[J].Rock and Soil Mechanics,2005,26(6): 865－868.

[3]鄧代強,高永濤,吳順川,等. 復雜應力下充填體破壞能耗試驗研究[J]. 巖土力學,2010,31(3): 737－742.DENG Dai-qiang,GAO Yong-tao,WU Shun-chuan,et al.Experimental study of destructive energy dissipation properties of backfill under complicated stress condition[J]. Rock and Soil Mechanics,2010,31(3): 737－742.

[4]劉志祥,李夕兵,戴塔根,等. 尾砂膠結充填體損傷模型及與巖體的匹配分析[J]. 巖石力學與工程學報,2006,27(9): 1442－1446.LIU Zhi-xiang,LI Xi-bing,DAI Ta-gen,et al. On damage model of cemented tailings backfill and its match with rock mass[J]. Chinese Journal of Rock Mechanics and Engineering,2006,27(9): 1442－1446.

[5]趙奎,王曉軍,劉洪興,等. 布筋尾砂膠結充填體頂板力學性狀試驗研究[J]. 巖土力學,2011,32(1): 9－15.ZHAO Kui,WANG Xiao-jun,LIU Hong-xing,et al.Experimental study of mechanical behaviors of cemented tailings backfill roof with reinforcement[J]. Rock and Soil Mechanics,2011,32(1): 9－15.

[6]孫鈞. 巖石流變力學及其工程應用研究的若干進展[J].巖石力學與工程學報,2007,26(6): 1081－1106.SUN Jun. Rock rheological mechanics and its advance in engineering applications[J]. Chinese Journal of Rock Mechanics and Engineering,2007,26(6): 1081－1106.

[7]盧萍珍,曾靜,盛謙. 軟黏土蠕變試驗及其經(jīng)驗模型研究[J]. 巖土力學,2008,29(4): 1041－1044.LU Ping-zhen,ZENG Jing,SHENG Qian. Creep tests on soft clay and its empirical models[J]. Rock and Soil Mechanics,2008,29(4): 1041－1044.

[8]閻巖,王思敬,王恩志. 基于西原模型的變參數(shù)蠕變方程[J]. 巖土力學,2010,31(10): 3025－3035.YAN Yan,WANG Si-jing,WANG En-zhi. Creep equation of variable parameters based on Nishihara model[J]. Rock and Soil Mechanics,2010,31(10): 3025－3035.

[9]熊良宵,楊林德,張堯. 巖石的非定常Burgers模型[J].中南大學學報(自然科學版),2010,41(2): 679－684.XIONG Liang-xiao,YANG Lin-de,ZHANG Yao.Non-stationary Burgers model for rock[J]. Journal of Central South University(Science and Technology),2010,41(2): 679－684.

[10]夏才初,金磊,郭銳. 參數(shù)非線性理論流變力學模型研究進展及存在的問題[J]. 巖石力學與工程學報,2011,30(3): 454－463.XIA Cai-chu,JIN Lei,GUO Rui. Nonlinear theoretical rheological model for rock: A review and some problems[J]. Chinese Journal of Rock Mechanics and Engineering,2011,30(3): 454－463.

[11]袁海平,曹平,萬文,等. 分級加卸載條件下軟弱復雜礦巖蠕變規(guī)律研究[J]. 巖石力學與工程學報,2006,25(8): 1575－1581.YUAN Hai-ping,CAO Ping,WAN Wen,et al. Study on creep rules of soft and intricate ore-rock under step load and unload[J]. Chinese Journal of Rock Mechanics and Engineering,2006,25(8): 1575－1581.

[12]丁志坤,呂愛鐘. 巖石黏彈性非定常蠕變方程的參數(shù)辨識[J]. 巖土力學,2004,29(增刊1): 37－40.DING Zhi-kun,Lü Ai-zhong. Parameter identification ofnon-stationary creep equation of rocks[J]. Rock and Soil Mechanics,2004,29(Supp. 1): 37－40.

[13]呂愛鐘,丁志坤,焦春茂,等. 巖石非定常蠕變模型辨識[J]. 巖石力學與工程學報,2008,27(1): 16－21.Lü Ai-zhong,DING Zhi-kun,JIAO Chun-mao,et al.Identification of non-stationary creep constitutive models of rock[J]. Chinese Journal of Rock Mechanics and Engineering,2008,27(1): 16－21.

[14]陳云敏,魏新江,李育超. 邊坡非圓弧臨界滑動面的粒子群優(yōu)化算法[J]. 巖石力學與工程學報. 2006,25(7):1443－1449.CHEN Yun-ming,WEI Xin-jiang,LI Yu-chao. Locating non-circular critical slip surfaces by particle swarm optimization algorithm[J]. Chinese Journal of Rock Mechanics and Engineering,2006,25(7): 1443－1449.

[15]高瑋. 基于粒子群優(yōu)化的巖土工程反分析研究[J]. 巖土力學,2006,27(5): 795－798.GAO Wei. Back analysis algorithm in geotechnical engineering based on particle swarm optimization[J].Rock and Soil Mechanics,2006,27(5): 795－798.

[16]田明俊,周晶. 巖土工程參數(shù)反演的一種新方法[J]. 巖石力學與工程程學報,2005,24(9): 1492－1496.TIAN Ming-jun,ZHOU Jin. New algorithm for parameter inversion in geotechnical engineering[J]. Chinese Journal of Rock Mechanics and Engineering,2005,24(9): 1492－1496.

[17]蘇國韶,馮夏庭. 基于粒子群優(yōu)化算法的高地應力條件下硬巖本構模型的參數(shù)辨識[J]. 巖石力學與工程學報. 2005,24(17): 3029－3034.SU Guo-shao,FENG Xia-ting. Parameter identification of constitutive model for hard rock under high in-situ stress condition using particle swarm optimization algorithm[J].Chinese Journal of Rock Mechanics and Engineering,2005,14(17): 3029－3034.

[18]褚衛(wèi)江,徐衛(wèi)亞,楊圣奇,等. 基于 FLAC3D黏彈塑性流變模型的二次開發(fā)研究[J]. 巖土力學,2006,27(11):2005－2010.CHU Wei-jiang,XU Wei-ya,YANG Sheng-qi,et al.Secondary development of a viscoelasto-plastic rheological constitutive model of rock based on FLAC3D[J]. Rock and Soil Mechanics,2006,27(11):2005－2010.

[19]羅潤林,阮懷寧,朱昌星. 基于粒子群一最小二乘法的巖石流變模型參數(shù)反演[J]. 遼寧工程技術大學,2009,28(5): 750－753.RUO Lun-lin,RUAN Huai-ning,ZHU Chang-xing.Parameter inversion of rock creep model based on PSO-least square method[J]. Journal of Liaoning Technical University,2009,28(5): 750－753.