考慮溫度效應的花崗巖分數階蠕變模型研究

(成都理工大學 地質災害防治與地質環境保護國家重點實驗室，四川 成都 610059)

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

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

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

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

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

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

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

2 分數階蠕變模型

2.1 分數階微積分

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

(1)

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

2.2 分數階蠕變元件

由上述分數階微積分的定義，將其引入到花崗巖蠕變本構模型構建之中。其本構關系如下：

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

(2)

式中，η為黏性系數。

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

(3)

此時，該元件即為Newton黏壺，代表理想流體； 當β=0時，該元件變為彈簧元件，代表理想彈性體?？梢?，分數階蠕變元件是一種可以用來模擬介于理想流體與理想彈性體之間的材料模型[13-15]。

2.3 蠕變模型的建立

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

圖2 蠕變本構模型示意Fig.2 A schematic diagram of creep constitutive model

由圖2可知，此時若令總應變為ε，塑性體的應變為εvp，分數階蠕變元件應變為εve，胡克體的應變為εe，繼而可得出總應變的表達式如下：

ε=εe+εve+εvp

(4)

現分別表述圖2中當總應力為σ時各元件的應力應變關系。

(1) 胡克體：

(5)

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

(2) 分數階蠕變元件：

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

(6)

(7)

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

(8)

式中,σs為屈服應力。

根據元件串、并聯理論，可得出總應力與彈性體的應力和摩擦滑塊應力之間的關系如下：

σ=σd+σp

(9)

式中,σd為分數階蠕變元件的應力。

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

εvp=0

(10)

當σ≥σs時，結合式(9)，可得本構關系：

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

(11)

(12)

式中,η1為黏塑性體中Abel黏壺的黏性系數，γ為其求導階數。

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

(13)

即

(14)

綜上所述，分數階黏彈塑性蠕變模型的本構方程可表示為

(15)

3 本構模型的驗證及參數分析

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

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

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

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

由圖3可知，無論是彈性模量E0、黏性系數η0，還是求導階數β，都是隨著溫度的上升而呈現減小的趨勢。對E0而言，根據其物理意義，在保持軸向應力不變的基礎上，溫度持續增加，應變會呈現線性增長的趨勢。3種參數的擬合結果如下：

E0=-91.197T+167980

(16)

η0=1.23098×1011·T-3.99813

(17)

(18)

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

(19)

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

4 分數階蠕變模型參數分析

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

4.1 彈性模量

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

圖3 參數擬合結果比較Fig.3 Comparison of the fitting parameters

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

4.2 黏性系數

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

圖5 黏性系數值對應變的影響Fig.5 The influence of the viscosity value on strain

4.3 求導階數

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

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

圖6 求導階數值對應變的影響Fig.6 The influence of derivation order value on strain

5 結 論

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

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