尾砂膠結(jié)充填體蠕變模型及在FLAC3D二次開(kāi)發(fā)中的實(shí)驗(yàn)研究

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

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

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

1 前 言

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

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

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

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

2.1 蠕變?cè)囼?yàn)

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

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

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

圖2 分級(jí)增量加卸載蠕變?cè)囼?yàn)現(xiàn)場(chǎng)Fig.2 Creep test under step loading and unloading

2.2 蠕變模型

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

圖3 不同應(yīng)力水平下試件蠕變曲線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(chǔ) = 11.408 GPa，B = 0.00359 kPa-1，C = 1.025，相關(guān)系數(shù)R = 0.9998。EH和η變化較小，可直接求平均確定它們的取值，EH= 3.871 GPa，η= 41.102 GPa·h。

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

2.3 蠕變方程的建立

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

應(yīng)用粒子群算法進(jìn)行蠕變方程參數(shù)反演的關(guān)鍵是建立適應(yīng)度函數(shù)，根據(jù)蠕變公式（1）和分級(jí)加卸載蠕變?cè)囼?yàn)數(shù)據(jù)，建立適應(yīng)度計(jì)算公式：

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

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

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

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

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

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

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

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

3.2 廣義Kelvin模型差分形式

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

式中：σij為應(yīng)力張量；σm為球應(yīng)力張量；Sij為偏應(yīng)力張量；δij為kronecker符號(hào)；為新的球應(yīng)力張量；為舊的球應(yīng)力張量；K為體積模量；Δεkk為應(yīng)變?cè)隽康牡谝徊蛔兞?；為新的偏?yīng)變張量；為舊的偏應(yīng)變張量；Δt為每一迭代步的蠕變時(shí)間增量；為新的偏應(yīng)力張量；為舊的偏應(yīng)力張量；為偏應(yīng)變張量的增量；Gk為Kelvin體剪切模量；Ge為Hoek體剪切模量；η為黏滯系數(shù)。

3.3 充填體模型

通過(guò)對(duì)室內(nèi)充填體單軸蠕變?cè)囼?yàn)的研究，得到充填體的蠕變方程，在FLAC3D軟件中實(shí)現(xiàn)式（5）。

本文通過(guò)引入應(yīng)力強(qiáng)度概念來(lái)實(shí)現(xiàn)軟件自動(dòng)識(shí)別加載應(yīng)力水平，定義應(yīng)力強(qiáng)度為

式中：σ1、σ2、σ3為第一、第二、第三主應(yīng)力。

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

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

3.4 開(kāi)發(fā)模型的驗(yàn)證

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

為了使三維張量模型和一維模型相對(duì)應(yīng)，將三維張量模型簡(jiǎn)化為

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

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

進(jìn)行單軸分級(jí)加載蠕變?cè)囼?yàn)?zāi)M的步驟如下：

①通過(guò)模型調(diào)用命令，導(dǎo)入生成的動(dòng)態(tài)鏈接庫(kù)文件來(lái)導(dǎo)入蠕變模型。

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

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

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

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

4 結(jié) 語(yǔ)

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

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

[2]李一帆,張建明,鄧飛,等. 深部采空區(qū)尾砂膠結(jié)充填體強(qiáng)度特性試驗(yàn)研究[J]. 巖土力學(xué),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]鄧代強(qiáng),高永濤,吳順川,等. 復(fù)雜應(yīng)力下充填體破壞能耗試驗(yàn)研究[J]. 巖土力學(xué),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]劉志祥,李夕兵,戴塔根,等. 尾砂膠結(jié)充填體損傷模型及與巖體的匹配分析[J]. 巖石力學(xué)與工程學(xué)報(bào),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]趙奎,王曉軍,劉洪興,等. 布筋尾砂膠結(jié)充填體頂板力學(xué)性狀試驗(yàn)研究[J]. 巖土力學(xué),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]孫鈞. 巖石流變力學(xué)及其工程應(yīng)用研究的若干進(jìn)展[J].巖石力學(xué)與工程學(xué)報(bào),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]盧萍珍,曾靜,盛謙. 軟黏土蠕變?cè)囼?yàn)及其經(jīng)驗(yàn)?zāi)Ｐ脱芯縖J]. 巖土力學(xué),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]. 巖土力學(xué),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].中南大學(xué)學(xué)報(bào)(自然科學(xué)版),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ù)非線性理論流變力學(xué)模型研究進(jìn)展及存在的問(wèn)題[J]. 巖石力學(xué)與工程學(xué)報(bào),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]袁海平,曹平,萬(wàn)文,等. 分級(jí)加卸載條件下軟弱復(fù)雜礦巖蠕變規(guī)律研究[J]. 巖石力學(xué)與工程學(xué)報(bào),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]丁志坤,呂愛(ài)鐘. 巖石黏彈性非定常蠕變方程的參數(shù)辨識(shí)[J]. 巖土力學(xué),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]呂愛(ài)鐘,丁志坤,焦春茂,等. 巖石非定常蠕變模型辨識(shí)[J]. 巖石力學(xué)與工程學(xué)報(bào),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]陳云敏,魏新江,李育超. 邊坡非圓弧臨界滑動(dòng)面的粒子群優(yōu)化算法[J]. 巖石力學(xué)與工程學(xué)報(bào). 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]. 巖土力學(xué),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]. 巖石力學(xué)與工程程學(xué)報(bào),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]蘇國(guó)韶,馮夏庭. 基于粒子群優(yōu)化算法的高地應(yīng)力條件下硬巖本構(gòu)模型的參數(shù)辨識(shí)[J]. 巖石力學(xué)與工程學(xué)報(bào). 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黏彈塑性流變模型的二次開(kāi)發(fā)研究[J]. 巖土力學(xué),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]羅潤(rùn)林,阮懷寧,朱昌星. 基于粒子群一最小二乘法的巖石流變模型參數(shù)反演[J]. 遼寧工程技術(shù)大學(xué),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.