(成都理工大學(xué) 地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護(hù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，四川 成都 610059)

溫度是影響巖石物理力學(xué)性質(zhì)的重要外因之一，尤其是對(duì)高放射性核廢料處理等工程有較大影響。由溫度導(dǎo)致的巖石力學(xué)性質(zhì)的變化規(guī)律是當(dāng)前巖石力學(xué)界研究的一個(gè)焦點(diǎn)問題。探究并獲得宏觀可用的考慮溫度變化的巖石蠕變模型對(duì)于指導(dǎo)巖石工程建設(shè)具有重要意義。

目前，溫度導(dǎo)致巖石蠕變行為機(jī)理的研究已得到諸多研究者的關(guān)注。Chan K S等[1]提出了一種考慮損傷斷裂的鹽巖蠕變模型，并通過蠕變?cè)囼?yàn)數(shù)據(jù)驗(yàn)證了模型合理性；郤保平等[2]進(jìn)行了層狀鹽巖蠕變特性試驗(yàn)，建立層狀鹽巖的穩(wěn)態(tài)蠕變率本構(gòu)方程；高小平等[3]對(duì)不同溫度處理后的鹽巖蠕變特性進(jìn)行研究，得到了巖石穩(wěn)態(tài)蠕變率本構(gòu)方程；周廣磊[4]等建立了溫度-應(yīng)力耦合作用下脆性巖石時(shí)效蠕變損傷模型，并在COMSOL的基礎(chǔ)上進(jìn)行二次開發(fā)，給出了溫度-應(yīng)力耦合作用下巖石時(shí)效蠕變損傷模型的數(shù)值求解方法；陳亮[5]等利用三維聲發(fā)射實(shí)時(shí)監(jiān)測(cè)信息，開展北山花崗巖的蠕變變形特征以及加載條件對(duì)其蠕變破壞過程的影響研究；楊春和[6]等給出了一個(gè)反映鹽巖蠕變?nèi)^程的非線性蠕變本構(gòu)方程；朱昌星[7]等在非線性黏彈塑性流變模型的基礎(chǔ)上，根據(jù)時(shí)效損傷和損傷加速門檻值的特點(diǎn)，建立了非線性蠕變損傷模型，并通過對(duì)板巖進(jìn)行試驗(yàn)研究，驗(yàn)證了模型的合理性；丁靖洋等[8-9]基于分?jǐn)?shù)階微分理論，利用變黏性系數(shù)Abel黏壺替代西原正夫模型中Newton黏壺的方法，建立了鹽巖蠕變損傷本構(gòu)模型。

由以上論述可知，在建立反映溫度影響的蠕變模型時(shí)，常有兩種代表性的做法：① 通過熱力學(xué)理論中關(guān)于能量損傷的辦法解決；② 一種經(jīng)驗(yàn)方法，即通過擬合蠕變參數(shù)隨溫度變化的關(guān)系，而后再代入蠕變模型中。相比之下，第一種方法的理論性更高，第二種做法則更為簡(jiǎn)單。而在巖石力學(xué)界，后者仍然是一種通用的研究手段。相應(yīng)的，采用后者時(shí)，蠕變本構(gòu)的選取對(duì)最終結(jié)果影響性較大，因而采用合理的蠕變模型尤為重要。本文采用后者建立花崗巖的溫度-蠕變耦合模型。基于已被廣泛認(rèn)可的分?jǐn)?shù)階蠕變模型，分析蠕變參數(shù)隨溫度的變化，然后獲得巖石溫度蠕變耦合模型，并對(duì)其模型參數(shù)展開詳細(xì)的討論。

1 不同溫度下花崗巖的蠕變特性

劉泉聲等[10]進(jìn)行了不同溫度下花崗巖的蠕變?cè)囼?yàn)。采用巖芯鉆取法將試件加工成直徑為50 mm，高度為100 mm的圓柱體。試驗(yàn)過程中，將三軸室的溫度按2℃/min的速率升至所需的溫度，保持該溫度4 h后開始加載試驗(yàn)。在溫度T=20 ℃，60 ℃，80 ℃，100 ℃，200 ℃，300 ℃條件下，保持軸向應(yīng)力為120 MPa，對(duì)三峽花崗巖進(jìn)行了一系列的單軸抗壓蠕變?cè)囼?yàn)。試驗(yàn)結(jié)果如圖1所示。

圖1 不同溫度下花崗巖的蠕變?cè)囼?yàn)曲線[8]Fig.1 Creep curves of granite under different temperatures

由圖1可以看出：在同一時(shí)間段，隨著溫度的升高，瞬時(shí)應(yīng)變不斷增加。若以彈性模量為研究對(duì)象，就會(huì)發(fā)現(xiàn)隨著溫度的上升，瞬時(shí)的彈性模量值會(huì)隨之減?。辉谙嗤妮S向應(yīng)力作用下，應(yīng)變率也隨著溫度的增長(zhǎng)而增加[11]。

2 分?jǐn)?shù)階蠕變模型

2.1 分?jǐn)?shù)階微積分

分?jǐn)?shù)階微積分有多種定義。現(xiàn)采用最常見的Riemann-Liouville分?jǐn)?shù)階微積分定義方法。先給出分?jǐn)?shù)階積分定義[12]：設(shè)f在(0，+∞)上連續(xù)且屬于(0，+∞)的任何有限子區(qū)間內(nèi)都可積，由t>0，有Re(γ)>0，此時(shí)可得出：

(1)

式(1)即為函數(shù)f(x)的β階分?jǐn)?shù)階積分，其中，Γ(β)為Gamma函數(shù)。

2.2 分?jǐn)?shù)階蠕變?cè)?/h3>

由上述分?jǐn)?shù)階微積分的定義，將其引入到花崗巖蠕變本構(gòu)模型構(gòu)建之中。其本構(gòu)關(guān)系如下：

σ=ηdβε(t)/dtβ0≤β≤1

(2)

式中，η為黏性系數(shù)。

令σ為一常數(shù)σ0，此時(shí)將式(2)兩邊都進(jìn)行分?jǐn)?shù)階積分?？擅枋鋈渥冃袨榈谋磉_(dá)式如下：

(3)

此時(shí)，該元件即為Newton黏壺，代表理想流體； 當(dāng)β=0時(shí)，該元件變?yōu)閺椈稍?，代表理想彈性體?？梢姡?jǐn)?shù)階蠕變?cè)且环N可以用來模擬介于理想流體與理想彈性體之間的材料模型[13-15]。

2.3 蠕變模型的建立

分?jǐn)?shù)階黏彈塑性蠕變模型如圖2所示。

圖2 蠕變本構(gòu)模型示意Fig.2 A schematic diagram of creep constitutive model

由圖2可知，此時(shí)若令總應(yīng)變?yōu)棣?，塑性體的應(yīng)變?yōu)棣舦p，分?jǐn)?shù)階蠕變?cè)?yīng)變?yōu)棣舦e，胡克體的應(yīng)變?yōu)棣舉，繼而可得出總應(yīng)變的表達(dá)式如下：

ε=εe+εve+εvp

(4)

現(xiàn)分別表述圖2中當(dāng)總應(yīng)力為σ時(shí)各元件的應(yīng)力應(yīng)變關(guān)系。

(1) 胡克體：

(5)

式中，E0為胡克體中彈簧的彈性模量。

(2) 分?jǐn)?shù)階蠕變?cè)?/p>

σ=η0dβεve(t)/dtβ

(6)

(7)

(3) 黏塑性體中，摩擦滑塊應(yīng)力σp的值則是根據(jù)總應(yīng)力與屈服應(yīng)力的大小來決定的，即

(8)

式中,σs為屈服應(yīng)力。

根據(jù)元件串、并聯(lián)理論，可得出總應(yīng)力與彈性體的應(yīng)力和摩擦滑塊應(yīng)力之間的關(guān)系如下：

σ=σd+σp

(9)

式中,σd為分?jǐn)?shù)階蠕變?cè)膽?yīng)力。

當(dāng)σ<σs時(shí)，由式(8)和(9)可得σd=0，即

εvp=0

(10)

當(dāng)σ≥σs時(shí)，結(jié)合式(9)，可得本構(gòu)關(guān)系：

η1dγεvp(t)/dtγ+σs=σ

(11)

(12)

式中,η1為黏塑性體中Abel黏壺的黏性系數(shù)，γ為其求導(dǎo)階數(shù)。

由初始條件t=0時(shí)εvp=0，求解式(12)可得：

(13)

即

(14)

綜上所述，分?jǐn)?shù)階黏彈塑性蠕變模型的本構(gòu)方程可表示為

(15)

3 本構(gòu)模型的驗(yàn)證及參數(shù)分析

現(xiàn)將圖1中不同溫度對(duì)應(yīng)不同時(shí)刻的應(yīng)變值提取出來，代入式(15)。用差分進(jìn)化法擬合得到不同溫度下參數(shù)E0、η0以及β的取值(見表1)。

表1 彈性模量、黏性系數(shù)和求導(dǎo)階數(shù)隨溫度的變化Tab.1 Variation of elastic modulus, viscosity coefficient and derivation order

擬合結(jié)果與試驗(yàn)結(jié)果之間的吻合度較好，說明擬合得到的參數(shù)值較為合理。

將表1中的E0、η、β分別與溫度T進(jìn)行擬合，擬合結(jié)果如圖3所示。最終可以得到花崗巖不同時(shí)間段、不同溫度下的ε值。

由圖3可知，無論是彈性模量E0、黏性系數(shù)η0，還是求導(dǎo)階數(shù)β，都是隨著溫度的上升而呈現(xiàn)減小的趨勢(shì)。對(duì)E0而言，根據(jù)其物理意義，在保持軸向應(yīng)力不變的基礎(chǔ)上，溫度持續(xù)增加，應(yīng)變會(huì)呈現(xiàn)線性增長(zhǎng)的趨勢(shì)。3種參數(shù)的擬合結(jié)果如下：

E0=-91.197T+167980

(16)

η0=1.23098×1011·T-3.99813

(17)

(18)

將式(16)、(17)、(18)聯(lián)立代入式(15)可得：

(19)

式(19)即為考慮溫度效應(yīng)的蠕變方程，可作為花崗巖蠕變研究的一個(gè)參考。

4 分?jǐn)?shù)階蠕變模型參數(shù)分析

為研究各參數(shù)對(duì)蠕變模型的影響規(guī)律，采用控制變量的方法，分別對(duì)彈性模量E0、黏性系數(shù)η0和求導(dǎo)階數(shù)β進(jìn)行分敏感性分析。模型參數(shù)取值：σ=120 MPa，E0=145 GPa，β=0.25，η0=65 GPa·h0.25。

4.1 彈性模量

分析彈性模量時(shí)，黏性系數(shù)η0和求導(dǎo)階數(shù)β保持不變，可得到不同E0狀態(tài)下的蠕變曲線。由圖4可知，E0決定了瞬時(shí)蠕變強(qiáng)度以及整體蠕變強(qiáng)度。隨著時(shí)間的增加，應(yīng)變率基本不變；隨著E0的降低，達(dá)到同一應(yīng)變值的時(shí)間逐漸增加。這表明在蠕變的前兩個(gè)階段，較小的E0值，會(huì)縮短從加速蠕變到穩(wěn)態(tài)蠕變的時(shí)間，但對(duì)蠕變的整體趨勢(shì)影響較小。

圖3 參數(shù)擬合結(jié)果比較Fig.3 Comparison of the fitting parameters

圖4 彈性模量變化值對(duì)應(yīng)變的影響Fig.4 The influence of change elastical of modulus on strain

4.2 黏性系數(shù)

分析黏性系數(shù)時(shí)，彈性模量E0和求導(dǎo)階數(shù)β保持不變。圖5可知，隨著η0的增加，花崗巖進(jìn)入穩(wěn)態(tài)蠕變的時(shí)間大幅降低，且當(dāng)η0值增加幅度較大時(shí)，這種現(xiàn)象尤其明顯。此外，從瞬態(tài)蠕變過渡到穩(wěn)態(tài)蠕變的時(shí)間也會(huì)大大縮短。

圖5 黏性系數(shù)值對(duì)應(yīng)變的影響Fig.5 The influence of the viscosity value on strain

4.3 求導(dǎo)階數(shù)

保持彈性模量、黏性系數(shù)不變，改變求導(dǎo)階數(shù)，得到不同求導(dǎo)階數(shù)下巖石蠕變曲線。由圖6可知，無論β如何取值，初始蠕變值均為一定值，即β對(duì)初始瞬時(shí)蠕變值影響較小，但對(duì)最終蠕變值的影響較大。當(dāng)β較小時(shí)，應(yīng)變速率較低，加速蠕變階段不明顯，且更容易達(dá)到穩(wěn)定蠕變階段。

綜上所述，彈性模量對(duì)花崗巖蠕變強(qiáng)度的變化起主要作用，卻并不影響總體蠕變趨勢(shì)。就黏性系數(shù)與求導(dǎo)階數(shù)而言，在溫度影響下，隨著參數(shù)值的降低，黏性系數(shù)會(huì)增加巖石進(jìn)入穩(wěn)態(tài)蠕變的時(shí)間，后者則恰恰相反。綜合來看黏性系數(shù)影響較大。

圖6 求導(dǎo)階數(shù)值對(duì)應(yīng)變的影響Fig.6 The influence of derivation order value on strain

5 結(jié) 論

(1) 根據(jù)分?jǐn)?shù)階黏彈塑性蠕變模型擬合不同溫度下巖石蠕變?cè)囼?yàn)結(jié)果，獲得了考慮溫度效應(yīng)的蠕變方程，對(duì)類似高放射性廢物的處置工程等提供了參考。

(2) 對(duì)分?jǐn)?shù)階黏彈塑性蠕變模型的蠕變參數(shù)進(jìn)行了敏感性分析，揭示了該模型參數(shù)對(duì)蠕變的影響規(guī)律。研究表明：溫度的上升導(dǎo)致蠕變參數(shù)值減小，其中彈性模量與黏性系數(shù)的值均與花崗巖到達(dá)穩(wěn)態(tài)蠕變階段的時(shí)間呈負(fù)相關(guān)關(guān)系，而求導(dǎo)階數(shù)值則相反。