考慮溫度效應(yīng)的花崗巖分數(shù)階蠕變模型研究

(成都理工大學 地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護國家重點實驗室，四川 成都 610059)

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

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

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

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

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

圖1 不同溫度下花崗巖的蠕變試驗曲線[8]Fig.1 Creep curves of granite under different temperatures

由圖1可以看出：在同一時間段，隨著溫度的升高，瞬時應(yīng)變不斷增加。若以彈性模量為研究對象，就會發(fā)現(xiàn)隨著溫度的上升，瞬時的彈性模量值會隨之減小；在相同的軸向應(yīng)力作用下，應(yīng)變率也隨著溫度的增長而增加[11]。

2 分數(shù)階蠕變模型

2.1 分數(shù)階微積分

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

(1)

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

2.2 分數(shù)階蠕變元件

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

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

(2)

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

令σ為一常數(shù)σ0，此時將式(2)兩邊都進行分數(shù)階積分。可描述蠕變行為的表達式如下：

(3)

此時，該元件即為Newton黏壺，代表理想流體； 當β=0時，該元件變?yōu)閺椈稍?，代表理想彈性體?？梢姡謹?shù)階蠕變元件是一種可以用來模擬介于理想流體與理想彈性體之間的材料模型[13-15]。

2.3 蠕變模型的建立

分數(shù)階黏彈塑性蠕變模型如圖2所示。

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

由圖2可知，此時若令總應(yīng)變?yōu)棣?，塑性體的應(yīng)變?yōu)棣舦p，分數(shù)階蠕變元件應(yīng)變?yōu)棣舦e，胡克體的應(yīng)變?yōu)棣舉，繼而可得出總應(yīng)變的表達式如下：

ε=εe+εve+εvp

(4)

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

(1) 胡克體：

(5)

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

(2) 分數(shù)階蠕變元件：

σ=η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為分數(shù)階蠕變元件的應(yīng)力。

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

εvp=0

(10)

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

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

(11)

(12)

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

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

(13)

即

(14)

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

(15)

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

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

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

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

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

由圖3可知，無論是彈性模量E0、黏性系數(shù)η0，還是求導(dǎo)階數(shù)β，都是隨著溫度的上升而呈現(xiàn)減小的趨勢。對E0而言，根據(jù)其物理意義，在保持軸向應(yīng)力不變的基礎(chǔ)上，溫度持續(xù)增加，應(yīng)變會呈現(xiàn)線性增長的趨勢。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)的蠕變方程，可作為花崗巖蠕變研究的一個參考。

4 分數(shù)階蠕變模型參數(shù)分析

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

4.1 彈性模量

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

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

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

4.2 黏性系數(shù)

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

圖5 黏性系數(shù)值對應(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可知，無論β如何取值，初始蠕變值均為一定值，即β對初始瞬時蠕變值影響較小，但對最終蠕變值的影響較大。當β較小時，應(yīng)變速率較低，加速蠕變階段不明顯，且更容易達到穩(wěn)定蠕變階段。

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

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

5 結(jié) 論

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

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