徐國文， 何 川， 代 聰， 王士民

(西南交通大學(xué)交通隧道工程教育部重點(diǎn)實(shí)驗(yàn)室，四川 成都610031)

巖石的流變特性是巖石內(nèi)在的時(shí)間特性，為了更準(zhǔn)確地描述巖石的三階段流變變形特性，許多學(xué)者考慮巖石流變過程中的損傷，建立流變本構(gòu)模型時(shí)引入了損傷力學(xué)理論.

高賽紅等在流變模型中考慮了硬化-損傷效應(yīng)［1］;陳衛(wèi)忠等考慮了黏滯系數(shù)損傷對(duì)巖石結(jié)構(gòu)長期特性的影響［2］;朱昌新等研究了巖石蠕變的時(shí)效損傷［3］;K S Chan 等將損傷理論引入鹽巖流變的研究中［4］;楊圣奇等將巖石蠕變損傷演化特性分為2 個(gè)階段，考慮了損傷對(duì)巖石流變特性的影響［5］.

由于實(shí)際工程的需要，一些學(xué)者借助商業(yè)軟件開發(fā)平臺(tái)進(jìn)行模型的二次開發(fā).以FLAC3D 軟件為例，楊文東等將Burgers 損傷模型用于邊坡長期安全的研究［6］;黃明等對(duì)隧道圍巖在含水條件下的劣化進(jìn)行了分析［7］.

本文在五參數(shù)廣義Kelvin 模型的基礎(chǔ)上，基于Lemaitre 應(yīng)變等效原理，提出了考慮巖石蠕變過程中材料劣化效應(yīng)的廣義Kelvin 蠕變損傷模型，并對(duì)該模型進(jìn)行了二次開發(fā);結(jié)合FLAC3D 數(shù)值計(jì)算程序及Matlab 平臺(tái)，采用粒子群算法、模擬退火算法與FLAC3D 相結(jié)合的方法對(duì)已有試驗(yàn)數(shù)據(jù)進(jìn)行反演.

1 廣義Kelvin 蠕變損傷模型

廣義Kelvin 模型(圖1(a))由1 個(gè)彈性體和2 個(gè)Kelvin 體構(gòu)成，不能反映巖石的塑性特性和加速蠕變過程.因此，將流變參數(shù)的損傷特性和M-C(摩爾-庫倫)元件［8］引入該模型，建立廣義Kelvin蠕變損傷模型(圖1(b)).

圖1 2 種流變本構(gòu)模型Fig.1 Rheological constitutive models

(1)σ≤σf時(shí)

本構(gòu)方程為:

式中:ε 為總應(yīng)變;σ 和σf分別為總應(yīng)力和屈服應(yīng)力;Ei和ηi(i =1，2)分別為2 個(gè)Kelvin 損傷體的彈性模量和黏滯系數(shù);Eh為彈性損傷體的彈性模量;Dt=1 - e-at為損傷變量［6］，其中a 為損傷參數(shù)，t 為時(shí)間.

由式(1)得到一維蠕變方程:

三維軸向蠕變方程為:

式中:ε1為軸向蠕變;K 為體積模量;G 為剪切模量;σ1和σ3分別為最大與最小主應(yīng)力;Gb和Hb(b=1，2)分別為2 個(gè)Kelvin 損傷體的三維剪切模量和三維黏滯系數(shù).

(2)σ ＞σf時(shí)

模型總應(yīng)變?yōu)?

式中:εij為總應(yīng)變;為彈性損傷體的應(yīng)變;和分別為2 個(gè)Kelvin 損傷體的應(yīng)變;為塑性體的應(yīng)變.

偏量行為由式(5)描述:

式中:Sij為偏應(yīng)力.

式中:δij為克羅內(nèi)克函數(shù)(當(dāng)i =j 時(shí)，δij=1;當(dāng)i≠j時(shí)，δij=0);λ 為塑性指示因子;為塑性體應(yīng)變.

式中，g 為與屈服準(zhǔn)則對(duì)應(yīng)的勢(shì)函數(shù).

塑性體采用帶拉伸截止限的摩爾-庫倫屈服準(zhǔn)則，屈服函數(shù)

采用非關(guān)聯(lián)流動(dòng)準(zhǔn)則，有勢(shì)函數(shù)

體應(yīng)力速率

式中:σV為體應(yīng)力;eV為總的體應(yīng)變.

模型三維本構(gòu)方程的差分形式為:

式中:Δt 為時(shí)間步;Ab=1 +GbΔt/(2Hb);Bb=1 -GbΔt/(2Hb);和分別為Δt 內(nèi)新、舊應(yīng)力偏量;和分別為Kelvin 體的新、舊應(yīng)變偏量.

式中:

Δeij為Δt 內(nèi)的總應(yīng)變?cè)隽繛棣 內(nèi)的塑性應(yīng)變?cè)隽?

體應(yīng)力的差分表達(dá)式為:

同理，Kelvin 體新的體應(yīng)變

2 模型的程序?qū)崿F(xiàn)與驗(yàn)證

基于FLAC3D，用C + +語言對(duì)模型進(jìn)行了二次開發(fā).廣義Kelvin 蠕變損傷模型考慮了巖石的塑性和損傷特性，現(xiàn)用算例驗(yàn)證模型的正確性.

(1)塑性力學(xué)特性

對(duì)模型的塑性特性進(jìn)行分析. 在廣義Kelvin蠕變損傷模型(圖1(b))的基礎(chǔ)上去掉1 個(gè)Kelvin損傷體，且不考慮模型其余元件力學(xué)參數(shù)的損傷，此時(shí)廣義Kelvin 蠕變損傷模型退化為考慮塑性的Kelvin 模型(圖2). 模型參數(shù):K =30 GPa，Gh=40 GPa(Gh為三維剪切模量)，G1= 60 GPa，H1=80 GPa·d，c = 12 MPa，φ= 30°，ψ = 10°，σte=1 MPa.計(jì)算模型取邊長為1 m 的正方體，頂面施加20 MPa 均布荷載，其余各邊施加相應(yīng)的位移約束.圖3 為蠕變30 d 后的塑性區(qū).可見，2 種模型的塑性區(qū)范圍一致，且主要集中在隧道洞周.

圖2 考慮塑性的Kelvin 模型Fig.2 Kelvin model considering plasticity

圖3 塑性區(qū)分布Fig.3 Distribution of plastic zone

(2)損傷力學(xué)特性

采用圖4 所示模型進(jìn)行損傷特性驗(yàn)證.該模型尺寸(水平×豎向)為5 cm×10 cm，底部固定豎直方向位移，頂部施加30 MPa 均布?jí)毫?，記錄點(diǎn)A 的豎向位移.

三參數(shù)廣義Kelvin 模型(圖2)的黏彈性參數(shù)和廣義Kelvin 蠕變損傷模型的黏彈塑性參數(shù)取值同上. 新增元件參數(shù)取值:a = 0. 05 d-1，G2=60 GPa，H2=80 GPa·d. 圖5 給出了2 種模型點(diǎn)A豎向位移隨時(shí)間的變化曲線.

圖4 數(shù)值模型Fig.4 A numerical model

圖5 點(diǎn)A 的豎向位移Fig.5 Vertical displacement of point A

可見，三參數(shù)廣義Kelvin模型僅能反映衰減蠕變階段(圖5(a)中OA)及穩(wěn)態(tài)蠕變階段(圖5(a)中AB);而廣義Kelvin 蠕變損傷模型由于考慮了塑性與損傷，可反映巖石的非線性蠕變階段(圖5(b)中AB)和加速蠕變階段(圖5(b)中BC)，更符合巖石的特性.

3 模型參數(shù)智能辨識(shí)

3.1 參數(shù)辨識(shí)方法

將模擬退火算法［10］引入粒子群優(yōu)化算法［11］，稱為PSO(particle swarm optimization )-SA(simulated annealing)算法.具體步驟:

步驟1 建立目標(biāo)函數(shù)

(1)獲得試驗(yàn)數(shù)據(jù){(t1，ε1)，(t2，ε2)，…，(tn，εn)}，確定設(shè)計(jì)變量R=(K，Gh，G1，G2，H1，H2，a);

(2)用FLAC3D 軟件建立數(shù)值模型，采用廣義Kelvin 蠕變損傷模型，將設(shè)計(jì)變量值代入計(jì)算，得到與(1)中試驗(yàn)數(shù)據(jù)對(duì)應(yīng)的計(jì)算數(shù)據(jù){(t1，～ε1)，(t2，)，…，(tn，)};

(3)建立目標(biāo)函數(shù)F(R)，

步驟2 參數(shù)辨識(shí)

(1)在Matlab 環(huán)境下，采用PSO-SA 優(yōu)化算法對(duì)參數(shù)進(jìn)行初始化. 調(diào)用FLAC3D，將初始參數(shù)作為模型輸入，求對(duì)應(yīng)的變形量，通過式(17)求得每個(gè)粒子的適應(yīng)值，然后返回Matlab 程序，取優(yōu)更新個(gè)體與群體極值.

(2)計(jì)算粒子更新前、后的適應(yīng)值，更新前后適應(yīng)值之差為適應(yīng)值變化量J. 若J ＜0，則粒子位置得到更新;若exp(-J/θ)＜γ(其中θ 為溫度，γ 為取值［0，1］的隨機(jī)數(shù))，粒子位置同樣得到更新，否則不更新.若接受新值，降溫θ 變成kθ(k 為退火參數(shù)，k∈(0，1))，否則不降溫，并返回步驟2，直至滿足精度要求.

3.2 參數(shù)辨識(shí)結(jié)果

用粉砂巖軸向蠕變數(shù)據(jù)［12］和砂板巖蠕變數(shù)據(jù)［13］進(jìn)行參數(shù)辨識(shí).

M-C 體的材料參數(shù)由常規(guī)試驗(yàn)確定，具體可參考文獻(xiàn)［12，14］.以K、Gh、G1、G2、H1、H2和a 為待求參量. 粒子群搜索空間:K = Gh= G1= G2∈(0，500 GPa)，H1=H2∈(0，500 GPa·d )，a∈(0，1);群體規(guī)模為100，最大搜索次數(shù)為100.模擬退火算法參數(shù):初始溫度θmax=5 000 ℃，溫度最低取值θmin=0.01 ℃，溫度退火參數(shù)k =0.9. 反演得到的參數(shù)見表1.

圖6 為粉砂巖、砂板巖反演曲線與室內(nèi)試驗(yàn)結(jié)果的比較，圖7 為誤差函數(shù)曲線.從圖6 可見，反演曲線與試驗(yàn)結(jié)果吻合，說明模型能正確反映該巖石的各蠕變階段. 從圖7 可見，對(duì)于粉砂巖，PSO-SA優(yōu)化算法在18 個(gè)迭代步內(nèi)即達(dá)到了較高的收斂精度;對(duì)于砂板巖，由于其流變曲線較復(fù)雜，PSO-SA優(yōu)化算法在40 步左右才達(dá)到全局最優(yōu). 誤差曲線具有平臺(tái)段，原因在于，為了解決PSO 算法易進(jìn)入局部最優(yōu)的缺陷，PSO-SA 優(yōu)化算法有一定接受適應(yīng)值較大粒子個(gè)體的概率，從而使得其本身也有可能進(jìn)入局部最優(yōu)［15］.但從結(jié)果看，只要達(dá)到一定迭代步，算法可以獲得較好結(jié)果.

表1 流變參數(shù)反演結(jié)果Tab.1 Inversion results of rheological parameters

圖6 反演結(jié)果與試驗(yàn)結(jié)果比較Fig.6 Comparison of inversion and test results

圖7 誤差曲線Fig.7 Error curves

反演算法的參數(shù)選擇參考了已有研究成果，是基于經(jīng)驗(yàn)進(jìn)行的選擇，因此產(chǎn)生了一定誤差. 在智能算法中，如何優(yōu)化算法自身參數(shù)以及算法理論的完善，還需要進(jìn)行深入研究.

4 結(jié) 論

本文在五參數(shù)廣義Kelvin 模型基礎(chǔ)上，基于Lemaitre 應(yīng)變等效原理，提出了廣義Kelvin 蠕變損傷模型，并基于FLAC3D 軟件提供的接口對(duì)模型進(jìn)行了二次開發(fā)，獲得以下結(jié)論:

(1)廣義Kelvin 蠕變損傷模型考慮了材料塑性和損傷特性，可以表征巖石的非線性和加速蠕變階段.

(2)采用粒子群算法、模擬退火算法與FLAC3D 相結(jié)合的智能方法，對(duì)已有試驗(yàn)數(shù)據(jù)進(jìn)行反演，該模型能較好模擬低應(yīng)力下巖石的兩階段蠕變和高應(yīng)力下巖石的三階段蠕變效應(yīng).

(3)流變曲線越復(fù)雜，反演算法收斂速度越慢.但只要選定合適的參數(shù)值，反演算法在40 代左右就可以收斂到全局最優(yōu)解.

［1］ 高賽紅，曹平，汪勝蓮，等. 改進(jìn)的巖石非線性黏彈塑性蠕變模型及其硬化黏滯系數(shù)的修正［J］. 煤炭學(xué)報(bào)，2012，37(6):936-943.GAO Saihong，CAO Ping，WANG Shenglian，et al.Improved nonlinear viscoelasto-plastic rheological model of rock and its correction of hardening coefficient of viscosity［J］. Journal of China Coal Society，2012，37(6):936-943.

［2］ 陳衛(wèi)忠，王者超，伍國軍，等. 鹽巖非線性蠕變損傷本構(gòu)模型及其工程應(yīng)用［J］. 巖石力學(xué)與工程學(xué)報(bào)，2007，26(3):467-472.CHEN Weizhong，WANG Zhechao，WU Guojun，et al.Nonlinear creep damage constitutive model of rock salt and its application to engineering［J］. Chinese Journal of Rock Mechanics and Engineering，2007，26(3):467-472.

［3］ 朱昌星，阮懷寧，朱珍德，等. 巖石非線性蠕變損傷模型的研究［J］. 巖土工程學(xué)報(bào)，2008，30(10):1510-1513.ZHU Changxing，RUAN Huaining，ZHU Zhende，et al.Non-linear rheological damage model of rock［J］.Chinese Journal of Geotechnical Engineering，2008，30(10):1510-1513.

［4］ CHAN K S，BRODSKY N S，F(xiàn)OSSUM A F，et al.Damage-induced non-associated inelastic flow in rock salt［J］. International Journal of Plasticity，1994，10(6):623-642.

［5］ 楊圣奇，徐鵬. 一種新的巖石非線性流變損傷模型研究［J］. 巖土工程學(xué)報(bào)，2014，36(10):1846-1854.YANG Shenqi，XU Peng. A new nonlinear rheological damage model for rock［J］. Chinese Journal of Geotechnical Engineering，2014，36(10):1846-1854.

［6］ 楊文東. 壩基軟弱巖體的非線性蠕變損傷本構(gòu)模型及其工程應(yīng)用［D］. 濟(jì)南:山東大學(xué)巖土與結(jié)構(gòu)工程研究中心，2008.

［7］ 黃明，劉新榮，鄧濤. 基于含水劣化特性的隧道圍巖時(shí)效變形數(shù)值計(jì)算［J］. 巖土力學(xué)，2012，33(6):1876-1882.HUANG Ming，LIU Xinrong，DENG Tao. Numerical calculation of time effect deformations of tunnel surrounding rock in terms of water degradation［J］.Rock and Soil Mechanics，2012，33(6):1876-1882.

［8］ 袁海平，曹平，許萬忠，等. 巖體黏彈塑性本構(gòu)關(guān)系及改進(jìn)Burgers 蠕變模型［J］. 巖土工程學(xué)報(bào)，2006，28(6):796-799.YUAN Haiping，CAO Ping，XU Wanzhong，et al.Visco-elastic-plastic constitutive relationship of rock and modified Burgers creep model［J］. Chinese Journal of Geotechnical Engineering，2006，28(6):796-799.

［9］ Itasca Consulting Group，Inc. Fast Lagrangian analysis of continua in three dimensions (version 3.0):user's manual［R］. Minnesota: Itasca Consulting Group，Inc.，2003.

［10］ 史峰，王輝，胡斐，等. Matlab 智能算法30 個(gè)案例分析［M］. 北京:北京航空航天大學(xué)出版社，2011:101-134.

［11］ 李麗，牛奔. 粒子群優(yōu)化算法［M］. 北京:冶金工業(yè)出版社，2009:62-83.

［12］ 黃明. 含水泥質(zhì)粉砂巖蠕變特性及其在軟巖隧道穩(wěn)定性分析中的應(yīng)用研究［D］. 重慶:重慶大學(xué)土木工程學(xué)院，2010.

［13］ 蔣昱州，張明鳴，李良權(quán). 巖石非線性黏彈塑性蠕變模型研究及其參數(shù)識(shí)別［J］. 巖石力學(xué)與工程學(xué)報(bào)，2008，27(4):832-839.JIANG Yuzhou，ZHANG Mingming，LI Liangquan.Study on nonlinear viscoelasto-plastic creep model of rock and its parameter identification［J］. Chinese Journal of Rock Mechanics and Engineering，2008，27(4):832-839.

［14］ 陳裕民. 錦屏一級(jí)水電站建基巖體力學(xué)參數(shù)選取研究［D］. 成都:成都理工大學(xué)環(huán)境與土木工程學(xué)院，2010.

［15］ 劉衍民，牛奔. 新型粒子群算法理論與實(shí)踐［M］.北京:科學(xué)出版社，2013:45-68.