郭瑞凱，丁建華，趙 奎，谷睿智

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

充填體的分?jǐn)?shù)階微積分蠕變本構(gòu)模型及其在FLAC3D中的開(kāi)發(fā)應(yīng)用

郭瑞凱，丁建華，趙 奎，谷睿智

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

為了更加準(zhǔn)確地評(píng)估分析充填采礦法中實(shí)際礦山充填體礦柱的穩(wěn)定性，通過(guò)在廣義開(kāi)爾文體中引入分?jǐn)?shù)階黏性元件來(lái)模擬膠結(jié)充填體蠕變過(guò)程，根據(jù)試驗(yàn)數(shù)據(jù)對(duì)所得蠕變本構(gòu)方程利用爬山算法和全局優(yōu)化法相結(jié)合的方式進(jìn)行參數(shù)辨識(shí)，并將得到的蠕變本構(gòu)方程開(kāi)發(fā)進(jìn)巖土數(shù)值分析軟件FLAC3D中進(jìn)行數(shù)值計(jì)算。計(jì)算結(jié)果表明，改進(jìn)后的分?jǐn)?shù)階蠕變本構(gòu)方程很好地?cái)M合了試驗(yàn)數(shù)據(jù)，證明了改進(jìn)模型及其計(jì)算程序的正確性。該成果對(duì)礦山礦柱充填及穩(wěn)定分析具有實(shí)際指導(dǎo)意義。

膠結(jié)充填體；分?jǐn)?shù)階微積分；蠕變本構(gòu)模型；爬山算法；FLAC3D二次開(kāi)發(fā)

0 引言

隨著對(duì)能源需求量的增加以及礦山開(kāi)采強(qiáng)度的不斷加大，地表和地下淺部礦產(chǎn)資源逐漸減少，國(guó)內(nèi)外礦山開(kāi)始逐漸進(jìn)入資源的深部開(kāi)采[1]。充填體采礦法可以很好地控制地表的沉陷，解決因地下開(kāi)采引起的一系列難題[2]，同時(shí)在降低礦石損失貧化等方面也具有顯著的優(yōu)勢(shì)。有學(xué)者預(yù)測(cè)充填采礦法將成為21世紀(jì)最具潛力的采礦方法[3]，因此充填采礦法將被金屬礦山廣泛使用。近年來(lái)，地下礦山的開(kāi)采深度和強(qiáng)度在不斷增大，充填體的蠕變特性越來(lái)越受到重視，有關(guān)充填體蠕變特性的研究產(chǎn)生了大量的成果，趙奎等[4]將充填體制成直徑50 mm×100 mm的試樣，并在實(shí)驗(yàn)室進(jìn)行單軸蠕變?cè)囼?yàn)，同時(shí)利用Hoek-Kelvin模型來(lái)描述充填體的蠕變特性，將模型嵌入進(jìn)FLAC3D軟件中，很好地模擬了膠結(jié)充填體的蠕變特性。馬乾天等[5]通過(guò)單軸壓縮試驗(yàn)研究塊石膠結(jié)充填體的失穩(wěn)方式，并利用基于顆粒流理論的PFC程序?qū)Τ涮铙w蠕變進(jìn)行模擬，模擬結(jié)果與試驗(yàn)結(jié)果吻合良好。趙樹(shù)果等[6]對(duì)充填體進(jìn)行單軸分級(jí)蠕變?cè)囼?yàn)，并通過(guò)引入Burgers模型來(lái)描述其蠕變特性，結(jié)果表明蠕變理論曲線很好地吻合了試驗(yàn)曲線。

目前對(duì)充填體蠕變本構(gòu)模型的研究大都采用傳統(tǒng)整數(shù)階元件組合的方式，然而對(duì)用分?jǐn)?shù)階流變?cè)M合來(lái)模擬充填體蠕變特性的研究?jī)H有少量研究成果。分?jǐn)?shù)階微分方程是將傳統(tǒng)微分方程中整數(shù)階次用分?jǐn)?shù)來(lái)替代，相比整數(shù)階微積分，分?jǐn)?shù)階微積分能夠更好地?cái)M合自然界中某些物理系統(tǒng)的動(dòng)態(tài)過(guò)程和記憶過(guò)程，因此逐漸有學(xué)者開(kāi)始將分?jǐn)?shù)階微積分引入到材料流變特性的研究當(dāng)中。Y Kawada[7]等研究了巖石和礦物黏彈性行為，并基于分?jǐn)?shù)階微積分用分?jǐn)?shù)階導(dǎo)數(shù)代表時(shí)間的延遲程度推導(dǎo)出其本構(gòu)方程。陳靜[8]基于分?jǐn)?shù)階微積分并結(jié)合巖石蠕變過(guò)程中的損傷效應(yīng)推導(dǎo)出花崗巖的蠕變本構(gòu)方程，并將研究成果應(yīng)用于錦屏水電站地下洞室群圍巖長(zhǎng)期穩(wěn)定性研究。吳斐[9]通過(guò)引入分?jǐn)?shù)階非線性黏壺元件較好地描述了巖鹽的蠕變?nèi)^(guò)程。

綜上所述，利用改進(jìn)的分?jǐn)?shù)階流變?cè)?lái)模擬充填體的流變性質(zhì)具有一定的實(shí)際意義。本文通過(guò)研究礦山采場(chǎng)充填體的蠕變特性，基于分?jǐn)?shù)階微積分推導(dǎo)出充填體在單軸加載下的蠕變本構(gòu)方程，并將得到蠕變本構(gòu)方程開(kāi)發(fā)進(jìn)巖土數(shù)值分析軟件FLAC3D中，對(duì)得到的充填體蠕變模型進(jìn)行蠕變?cè)囼?yàn)數(shù)值模擬，研究結(jié)果可應(yīng)用于礦山充填體礦柱的穩(wěn)定性分析。

1 蠕變?cè)囼?yàn)和蠕變模型

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

試驗(yàn)所用水泥和尾砂皆來(lái)自武山某礦，其中水泥和尾砂的配比為1∶20，將其澆注成50 mm×100 mm的圓柱，并將其置于自制的簡(jiǎn)易分級(jí)增量循環(huán)加卸載裝置上進(jìn)行單軸加載，按照其事先測(cè)定的單軸抗壓強(qiáng)度 1.5MPa，將其分成 210kPa、436kPa、662kPa、939 kPa、1 191 kPa五個(gè)級(jí)別進(jìn)行分級(jí)加載。加載過(guò)程中當(dāng)位移增量小于0.001 mm/h時(shí)，切換到下一個(gè)加載級(jí)別。加載得到的時(shí)間-應(yīng)變曲線如圖1所示。

圖1 單軸分級(jí)加載下充填體時(shí)間-應(yīng)變曲線Fig.1 Time-strain curve of filling body under uniaxial loading

由圖1可知充填體不同于普通巖石的脆性破壞、塑性破壞以及弱面剪切破壞，充填體表現(xiàn)出明顯的黏塑性特性，該類型蠕變形式存在明顯的減速蠕變和等速蠕變階段，無(wú)明顯的加速蠕變階段[10]。

1.2 蠕變模型

傳統(tǒng)的廣義開(kāi)爾文體由一個(gè)彈簧體和一個(gè)開(kāi)爾文體串聯(lián)組成，其力學(xué)模型如圖2所示。

圖2 廣義開(kāi)爾文體模型Fig.2 The generalized Kelvin model

圖2 中σ為應(yīng)力，MPa；k1和k2表示彈簧的彈性元件彈性參數(shù)，GPa；η表示黏壺的牛頓黏性系數(shù)，GPa·h；其蠕變方程如式（1）所示。

本文在廣義開(kāi)爾文體模型的基礎(chǔ)上，引入一個(gè)改進(jìn)的分?jǐn)?shù)階微積分黏性元件，該元件用于描述充填體蠕變過(guò)程中處于彈性和黏彈性之間的狀態(tài)，并反映蠕變過(guò)程中時(shí)效性。改進(jìn)之后的模型如圖3所示。

圖3 改進(jìn)后的廣義開(kāi)爾文體模型Fig.3 The modified generalized Kelvin model

有關(guān)分?jǐn)?shù)階微積分的定義主要有三種，分別是Grümwald-Letnikov 定義、Riemann-Liouville 定義和Caputo定義，試驗(yàn)研究中選擇最常用的Riemann-Liouville定義。Riemann-Liouville對(duì)分?jǐn)?shù)階導(dǎo)數(shù)的定義如式（2）所示。

其中Γ（n-α）為伽馬函數(shù)，伽馬函數(shù)的定義是：

吳斐等[11]根據(jù)Riemann-Liouville定義推導(dǎo)出分?jǐn)?shù)階黏性元件的應(yīng)變與時(shí)間的關(guān)系如式（4）所示。

因此改進(jìn)后的充填體蠕變本構(gòu)方程見(jiàn)式（5）。

其中，k1、k2為彈性元件的彈性參數(shù)，GPa；η1、η2分別為改進(jìn)黏性元件和傳統(tǒng)黏性元件的黏性系數(shù)，GPa·h；β表示分?jǐn)?shù)階次，且 0≤β≤1；t為蠕變時(shí)間，s。

2 參數(shù)識(shí)別

有關(guān)本構(gòu)模型參數(shù)辨識(shí)的研究有很多，常見(jiàn)的有最小二乘法、遺傳算法、神經(jīng)網(wǎng)絡(luò)法等，本試驗(yàn)選擇比較簡(jiǎn)單但效率較高的爬山算法。爬山算法是通過(guò)模擬現(xiàn)實(shí)中的爬山過(guò)程，開(kāi)始隨機(jī)選取一個(gè)值進(jìn)行逼近，每次都在選定值的臨近的空間中選擇最優(yōu)解作為當(dāng)前解，直到局部最優(yōu)解。其爬山算法存在的缺點(diǎn)是比較容易在辨識(shí)過(guò)程中陷入局部最優(yōu)解，而局部最優(yōu)解并不一定滿足全局最優(yōu)解。因此本研究借助七維高科有限公司（7D-Soft High Technology Inc.）開(kāi)發(fā)的優(yōu)化分析軟件1stOpt，通過(guò)應(yīng)用爬山算法并結(jié)合1stOpt軟件內(nèi)嵌的通用全局優(yōu)化法對(duì)蠕變本構(gòu)方程（5）進(jìn)行參數(shù)辨識(shí)，克服了爬山算法的不足，辨識(shí)結(jié)果見(jiàn)表1。

表1 改進(jìn)模型擬合參數(shù)Tab.1 Fitting parameters of the modified model

3 蠕變本構(gòu)模型在FLAC3D中的開(kāi)發(fā)與實(shí)現(xiàn)

將自定義的蠕變本構(gòu)模型開(kāi)發(fā)進(jìn)FLAC3D中需要將式（5）轉(zhuǎn)化成三維的有限差分形式，研究根據(jù)陳靜[2]推導(dǎo)的結(jié)果得出充填體蠕變本構(gòu)方程的三維有限差分形式如式（6）所示。

式中：

式中：

式中：GM表示有彈性元件和改進(jìn)黏性元件組成的類Maxwell體的剪切模量。

借助巖土數(shù)值分析軟件FLAC3D預(yù)留的端口，在Visual Studio 2010中進(jìn)行充填體蠕變本構(gòu)模型的編寫，編寫內(nèi)容主要包括頭文件Gamma函數(shù)、基類Constituive Model等，其中式（6）、（7）和（8）將被編寫進(jìn)基類ConstituiveModel的成員函數(shù)run（）中，生成可供FLAC3D計(jì)算調(diào)用的dll.文件，并將該文件嵌入到FLAC3D軟件包中，在計(jì)算過(guò)程中通過(guò)config cppudm命令即可對(duì)蠕變模型進(jìn)行調(diào)用。具體內(nèi)容許多文獻(xiàn)[12-16]已經(jīng)報(bào)道，本文不再贅述。

數(shù)值計(jì)算模型為與試驗(yàn)試件相同尺寸的圓柱，如圖2所示，模擬蠕變過(guò)程中利用FLAC3D內(nèi)嵌的fish語(yǔ)言編寫時(shí)步函數(shù)time_fun，計(jì)算過(guò)程中模型底部在y軸方向上施加約束，模型頂部施加垂直載荷，因模擬單軸加載，故模型側(cè)向不施加約束，蠕變分析時(shí)步設(shè)置為1×10-3，監(jiān)測(cè)點(diǎn)為圓柱頂部的中心原點(diǎn)。

圖4 數(shù)值計(jì)算模型Fig.4 Model of numerical calculation

4 數(shù)值模擬計(jì)算結(jié)果與試驗(yàn)數(shù)據(jù)的對(duì)比

FLAC3D模擬結(jié)果如圖5所示，將數(shù)值計(jì)算結(jié)果從FLAC3D中導(dǎo)出并與試驗(yàn)值進(jìn)行對(duì)比，由圖5和圖6可知，數(shù)值計(jì)算結(jié)果很好地?cái)M合了充填體的減速和等速蠕變階段，證實(shí)了本試驗(yàn)所開(kāi)發(fā)程序的正確性。

圖5 FLAC3D模擬充填體分級(jí)加載Fig.5 Hierarchical loading of FLAC3Dsimulated filling body

圖6 模擬數(shù)據(jù)與試驗(yàn)數(shù)據(jù)對(duì)比Fig.6 Comparison between simulated data and experimental data

5 結(jié)論

通過(guò)試驗(yàn)與數(shù)值模擬相結(jié)合，通過(guò)二次開(kāi)發(fā)，使理論得到的該配比下充填體分?jǐn)?shù)階微積分蠕變本構(gòu)方程更好地應(yīng)用于實(shí)際礦山中充填體礦柱的穩(wěn)定性分析，得出以下結(jié)論：

（1）該配比下的膠結(jié)充填體具有顯著的黏塑性特征，存在明顯的減速蠕變和等速蠕變階段。

（2）利用擬合軟件1stOpt中內(nèi)嵌的爬山算法和全局優(yōu)化法對(duì)改進(jìn)的蠕變本構(gòu)方程進(jìn)行參數(shù)辨識(shí)，辨識(shí)結(jié)果表明擬合參數(shù)合理有效，因此將爬山算法與全局優(yōu)化法相結(jié)合可以提高其精度和尋優(yōu)速度，對(duì)于參數(shù)辨識(shí)具有一定的意義。

（3）將數(shù)值計(jì)算結(jié)果與試驗(yàn)結(jié)果進(jìn)行對(duì)比，發(fā)現(xiàn)改進(jìn)的分?jǐn)?shù)階蠕變本構(gòu)方程可以較好地?cái)M合該配比下充填體的蠕變過(guò)程，并體現(xiàn)出了蠕變過(guò)程的時(shí)效性，從而證實(shí)了所開(kāi)發(fā)模型的適應(yīng)性和正確性。

[1] 何滿潮，謝和平，彭蘇萍，等.深部開(kāi)采巖體力學(xué)研究[J].巖石力學(xué)與工程學(xué)報(bào)，2005，24（16）：2803-2813.HE Manchao，XIE Heping，PENG Suping，et al.Study on rock mechanics in deep mining engineering[J].Chinese Journal of Rock Mechanics and Engineering，2005，24（16）：2803-2813.

[2] 史俊偉，陳章良，孫玉峰，等.充填開(kāi)采控制地表沉陷可靠性模糊綜合評(píng)判[J].中國(guó)鎢業(yè)，2015，30（4）：21-26.SHI Junwei，CHEN Zhangliang，SUN Yufeng，etal.Fuzzy comprehensive evaluation of reliability of surface subsidence control by filling mining[J].China Tungsten Industry，2015，30（4）：21-26.

[3] 胡 華，崔明義.高水材料硬化體特性及其充填體力學(xué)作用機(jī)理分析[J].礦業(yè)研究與開(kāi)發(fā)，2001，21（5）：23-25.HU Hua，CUI Mingyi.Characteristics of high water content hardening body and analysis on the mechanical mechanism of fill body[J].Mining Research and Development，2001，21（5）：23-25.

[4] 趙 奎，何 文，熊良宵，等.尾砂膠結(jié)充填體蠕變模型及在FLAC3D二次開(kāi)發(fā)中的實(shí)驗(yàn)研究 [J].巖土力學(xué)，2012（增刊1）：112-116.ZHAO Kui，HE Wen，XIONG Liangxiao，et al.Testing study of creep model for tailing cemented backfill and its secondary development based on FLAC3D[J].Rock and Soil Mechanics，2012（supply1）：112-116.

[5] 馬乾天，張東煒.基于顆粒流的塊石膠結(jié)充填體短時(shí)蠕變特性研究[J].礦業(yè)研究與開(kāi)發(fā)，2015（7）：68-71.MA Qiantian，ZHANG Dongwei.Study on short-time creep characteristics of rocky cemented filling body based on PFC[J].Mining Research and Development，2015（7）：68-71.

[6] 趙樹(shù)果，蘇東良，鄒 威.充填體分級(jí)加載蠕變?cè)囼?yàn)及模型參數(shù)智能辨識(shí)[J].礦業(yè)研究與開(kāi)發(fā)，2016，36（6）：54-57.ZHAO Shuguo，SU Dongliang，ZOU Wei.Study on creep properties of backfill under step loading and intelligent identification of constitutivemodelparameters[J].MiningResearchandDevelopment，2016，36（6）：54-57.

[7] KAWADA Y，YAJIMA T，NAGAHAMA H.Fractional-order derivative and time-dependent viscoelastic behaviour of rocks and minerals[J].Acta Geophysica，2013，61（6）：1690-1702.

[8] 陳 靜.硬巖時(shí)效變形和破壞機(jī)制及其工程行為分析 [D].北京：中國(guó)科學(xué)院大學(xué)，2013.CHEN Jing.Time-dependent deformation and fracture mechanisms of hard rock and engineering behavioural analysis[D].Beijing：The University of Chinese Academy of Sciences，2013.

[9] 吳 斐，謝和平，劉建鋒，等.分?jǐn)?shù)階黏彈塑性蠕變模型試驗(yàn)研究[J].巖石力學(xué)與工程學(xué)報(bào)，2014，33（5）：964-970.WU Fei，XIE Heping，LIU Jiangfeng，et al.Experimental study of fractional viscoelastic-plastic creep model[J].Chinese Journal of Rock Mechanics and Engineering，2014，33（5）：964-970.

[10] 楊 欣.充填體蠕變本構(gòu)模型及其工程應(yīng)用[D].贛州：江西理工大學(xué)，2011.YANG Xin.Creep constitutive model of filling body and its engineering application[D].Ganzhou：Jiangxi University of Science and Technology，2011.

[11] 吳 斐，劉建鋒，武志德，等.鹽巖的分?jǐn)?shù)階非線性蠕變本構(gòu)模型[J].巖土力學(xué)，2014，35（增刊 2）：162-167.WU Fei，LIU Jianfeng，WU Zhide，et al.Fractional nonlinear creep constitutive model of salt rock [J].Rock and Soil Mechanics，2014，35（supply2）：162-167.

[12] 藍(lán) 航，姚建國(guó)，張華興，等.基于FLAC3D的節(jié)理巖體采動(dòng)損傷本構(gòu)模型的開(kāi)發(fā)及應(yīng)用[J].巖石力學(xué)與工程學(xué)報(bào)，2008，27（3）：572-579.LAN Hang，YAO Jianguo，ZHANG Huaxing，et al.Development and application of constitutive model of jointed rock mass damage due to mining based on FLAC3D[J].Chinese Journal of Rock Mechanics and Engineering，2008，27（3）：572-579.

[13] 陳育民，徐鼎平.FLAC/FLAC3D基礎(chǔ)與工程實(shí)例[M].北京：中國(guó)水利水電出版社，2009.

[14] 周宏偉，王春萍，段志強(qiáng)，等.基于分?jǐn)?shù)階導(dǎo)數(shù)的鹽巖流變本構(gòu)模型[J].中國(guó)科學(xué)：物理學(xué) 力學(xué) 天文學(xué)，2012，42（3）：310-318.ZHOUHongwei，WANGChunping，DUANZhiqiang，etal.Rheological constitutive model of salt rock based on fractional derivative[J].Scientia Sinica Physica，Mechanicaamp;Astronomica，2012，42（3）：310-318.

[15]朱曉鵬.基于節(jié)理巖體損傷本構(gòu)模型的FLAC3D二次開(kāi)發(fā)及應(yīng)用[D].北京：中國(guó)地質(zhì)大學(xué)，2015.ZHU Xiaopeng.The secondary development of FLAC3DBased on Jointed rock mass damage model and its application[D].Beijing：China University of Geosciences，2015.

[16] 高文華，劉 正，張志敏.基于FLAC3D的粉砂巖壓縮蠕變?cè)囼?yàn)數(shù)值模擬研究[J].土木工程學(xué)報(bào)，2015（3）：96-102.GAO Wenhua，LIU Zheng，ZHANG Zhimin.Numerical simulation study on compression creep experiment of siltstone based on FLAC3D[J].China Civil Engineering Journal，2015（3）：96-102.

Fractional Calculus Creep Constitutive Model of Filling Body and its Development and Application in FLAC3D

GUO Ruikai,DING Jianhua,ZHAO Kui,GU Ruizhi
(Faculty of Resources and Environmental Engineering,Jiangxi University of Science and Technology,Ganzhou 341000,Jiangxi,China)

For an accurate analysis on the stability of filling pillars,this paper simulated the creeping process of cemented backfill by applying fractional viscous components in the generalized Kelvin system.Parameters were identified on the creep constitutive equation by combining Hill Climbing algorithm and global optimization method.The creep constitutive equation is developed into numerical analysis software FLAC3Dfor rock and soil numerical calculation.The calculated results show that the modified fractional creep constitutive equation fits the experimental data well,which testifies the correctness of the improved model and its calculation program.

cemented filling body;fractional calculus;creep constitutive model;hill climbing algorithm;FLAC3Dsecondary exploration

TD853；TU443

A

10.3969/j.issn.1009－0622.2017.05.005

2017-08-28

江西省科技計(jì)劃項(xiàng)目（20143ACG70010）；江西省主要學(xué)科學(xué)術(shù)和技術(shù)帶頭人培養(yǎng)計(jì)劃（20153BCB22004）；江西省“5511”優(yōu)勢(shì)科技創(chuàng)新團(tuán)隊(duì)（20165BCB19012）

郭瑞凱（1990-），男，河北邢臺(tái)人，碩士，主要從事巖石力學(xué)方面的研究。

趙 奎（1969-），男，安徽六安人，教授，主要從事巖石力學(xué)測(cè)試技術(shù)方面的研究。

（編輯：游航英）