基于分數階微積分的巖石硬性結構面蠕變本構模型

李任杰, 吉 鋒, 張津銘, 潘勇杰, 熊朝正

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

節理巖體普遍存在于巖體工程中,其力學性質等均較復雜,直接影響巖體的安全性。流變為巖石變形破壞中較為重要的一種表現形式,節理流變往往決定巖體的流變。巖體蠕變本構模型的研究是目前了解巖體蠕變特性的重點和熱點[1]。目前應用較為廣泛的本構模型構建方法主要為元件組合模型,其通過基本元件組合、改進以往的線性元件為非線性元件組合模型以及加入損傷參數等變量對蠕變過程進行擬合[2]。

孫海忠等[3]采用含分數階導數的Kelvin模型對軟土的黏彈性蠕變過程進行較好的模擬,并認為形式上簡單統一,且參數少,具有較為廣泛的應用前景。王紅娟等[4]于Kelvin體黏性元件中引入分數階函數,且彈性元件考慮彈性參數隨時間變化,對黏塑性體進行改進,構建變參數非線性黏彈性蠕變模型,經過驗證,模型模擬效果較好。陳家瑞等[5]通過引入分數階Kelvin模型,解決了整數解Kelvin的擬合效果不好的問題,最終證明模型能夠較好地模擬破碎泥巖的蠕變過程。郭佳奇等[6]將Kelvin模型中牛頓黏性體替換為含分數階微積分的FC元件,最后通過元件組合模型,形成基于分數階微積分的流變模型,通過與整數階流變模型相比較,其模擬效果較好。黃海峰等[7]基于Riemann-Liouville理論,于黏性元件中引入分數階微積分,使得更加符合巖石材料的流體性質,最終模型通過對紅層泥巖及相關試驗數據進行擬合,充分說明了本構模型的適用性。劉峻松等[8]基于分數階微積分理論,對蠕變本構模型中的基本元件進行改進,最終構建一種新的三元件蠕變損傷本構模型,通過辨識砂巖的試驗數據,證明了所建本構模型的合理適用性。唐皓等[2,9]利用分數階微積分的軟體元件替換西原模型中所有的黏壺元件,并將黏滯系數非定常化,得到分數階微積分改進的黃土西原模型,對黃土蠕變試驗數據進行辨識,結果證明模型的辨識效果較好。徐國文等[10]根據R-L分數階微積分理論,提出分數階微積分黏壺元件替換西原模型中的黏壺元件,結果表明模型能夠較好地對蠕變三階段進行反映。陳亮等[11]通過引入基于分數階微積分的軟體元件,對巖石蠕變的三階段進行很好的模擬。

上述研究將分數階微積分應用于蠕變本構模型中,且均取得較好的效果,但是多為完整巖石的蠕變本構模型,少有涉及貫通型硬性結構面和非貫通結構面的剪切蠕變本構模型研究。因此將分數階微積分應用于非貫通硬性結構面巖體剪切蠕變本構模型中具有較為重要的意義。

針對非貫通硬性結構面蠕變試驗及本構模型方面研究的缺乏,通過制作非貫通硬性結構面試樣,進行剪切蠕變試驗,揭示非貫通硬性結構面剪切蠕變特性,建立基于分數階微積分的巖石硬性結構面蠕變本構模型。

1 剪切蠕變試驗設計及成果

1.1 試驗設計及試驗設備

參考李銀平等[12]所提利用云母片模擬結構面的方法,模擬硬性結構面,進一步制作非貫通硬性結構面的試驗模型,進行剪切試驗及剪切蠕變試驗,揭示非貫通硬性結構面的剪切蠕變特性。試樣制作采用類巖石材料及其比例為水∶水泥∶標準砂=3∶8∶8的比例進行水泥砂漿的澆筑,待其固結24 h后進行拆模,并于水中養護20 d后進行相關試驗。試驗所用試樣制作流程如圖1所示。

圖1 模型制作流程Fig.1 Model making process

剪切蠕變試驗在地質災害防治與地質環境保護國家重點實驗室自主研發的YZJL-300型巖石剪切流變儀上進行。

剪切蠕變試驗基于單軸壓縮試驗和直剪試驗的成果上,采用陳氏分級加載法,進行非貫通硬性結構面的剪切蠕變試驗。對連通率k=(L4L5)/(L1L3)=0.64試樣加載不變的正應力σ,由低到高加載剪應力τ,每級剪應力如表1所示,每級剪應力持續時間為48 h,直至試樣破壞。

表1 剪切蠕變試驗參數

1.2 試驗成果

通過剪切蠕變試驗,在應力加載瞬間或增大至下一級應力時,非貫通硬性結構面試樣有明顯的瞬時變形現象。同時,剪切蠕變試驗過程中出現了較為完整的蠕變三階段,即衰減蠕變階段、穩速蠕變階段、加速蠕變階段,如圖2所示。

圖2 剪切流變試驗曲線Fig.2 Shear rheological test curve

2 分數階微積分軟體元件

2.1 分數階微積分簡介

分數階微積分的定義方式多樣,而最常用且應用廣泛的定義方式為Riemann-Liouville。分數階微積分核心在于將整數階微積分推廣至分數、復數層面。

函數f(t)的β階積分為[9]

(1)

式(1)中:t、t0、τ為自變量；β為階數。

分數階微分定義為

(2)

式(2)中:β∈(n-1,n],n為正整數,且β>0；Γ(β)為Gamma函數,即

(3)

式(3)中：Re代表取復數的實部。

2.2 基于分數階微積分建立的軟體元件

實際上,巖石為一種介于理想彈性體和理想流體之間的一種流體材料。

理想彈性體(Hooke體)的應力σ與應變ε之間呈正比例關系，即

(4)

式(4)中:τ為應力,MPa；ε為應變；E為彈性模量,GPa。經過改寫可以寫為[2,9,13]

(5)

理想流體(Newton體)應力-應變-時間關系為

(6)

式(6)中:η為黏性系數；t為時間。經過改寫可以寫為[2,9,13]

(7)

由于巖石材料在黏彈性階段既非理想彈性體也非理想流體的性質,恰好與式(5)、式(6)中的微分階數0和1相對應,因此其巖石的本構關系可能為

(8)

式(8)中：β指元件階數,β=0或1時,正好對應于理想狀態下的彈性體和流體。因此,該軟體元件(圖3)可應用于巖石材料變形的描述。

圖3 軟體元件Fig.3 Software element

當應力τ恒定時,該元件即可用來描述巖土材料的蠕變特性。對式(8)進行分數階積分,依據R-L分數階微積分算子理論,可得:

(9)

3 加速蠕變元件的提出

目前針對于加速蠕變過程的本構元件的提出以改進現有的基本元件、基本元件組合、引入M-C(摩爾-庫侖)塑性體元件來對加速蠕變過程進行模擬。

張清照等[14]、Zhang等[15]依據結構面剪切蠕變曲線的特點,提出一個非線性黏性元件,如圖4所示,以更好地對結構面剪切蠕變全過程進行模擬。

圖4 非線性黏性元件Fig.4 Nonlinear viscous element

非線性黏性元件的應力、應變關系為

(10)

即非線性黏性元件的本構模型為

(11)

當巖石應力達到屈服應力條件下,巖石蠕變開始進入加速蠕變階段。此時,蠕變應變速率開始快速增長,黏壺元件黏滯系數開始逐漸降低,依據Jean Lemaitre[8,16]應變等價性假說,對式(11)中的η進行非定常化,加速蠕變元件經非定常化之后的本構方程為

(12)

式(12)中：b為與流變有關的參數。

4 蠕變本構模型構建

通過對剪切蠕變試驗進行分析,剪切蠕變試驗體現了瞬時變形、衰減蠕變過程、穩速蠕變過程、加速蠕變過程。根據此試驗現象,建立了一個由彈簧元件、黏彈性元件、一個能夠反映加速蠕變過程的非線性黏性元件通過串并聯的方式組成基于分數階微積分的巖石硬性結構面蠕變本構模型,如圖5所示。

圖5 基于分數階微積分蠕變本構模型Fig.5 Creep constitutive model based on fractional calculus

由本構模型的并聯型關系可知,圖5中本構模型的總應變等于圖中1、2、3部分的總和,即

ε=ε1+ε2+ε3

(13)

(1)彈性體(胡克體)的應力-應變關系滿足Hooke定律,其應力-應變關系應滿足:

(14)

(2)分數階微積分Kelvin體中含基于分數階微積分所建立的軟體元件,能夠更好地對巖石的黏彈性變形過程進行模擬,其應力-應變關系應滿足:

(15)

式(15)中：εve、εe、ε分別為Kelvin體、Kelvin體中彈性元件、Kelvin體中黏性元件所產生的應變。

將式(8)和式(14)代入式(15),可得黏彈性體的應力-應變-時間關系為

(16)

經過轉換,即可得：

(17)

對式(17)進行拉普拉斯變換后可得

(18)

式(18)中：s表示拉普拉斯變換空間復變量。

假設,τ/η1=a,E2/η1=b,代入式(18)并整理得：

(19)

進一步通過對式(19)進行拉普拉斯的逆變換,即可得出：

(20)

(3)模型加速蠕變元件的本構模型為

(21)

式(21)中：λ的意義與式(12)中n的意義相同。

根據本構模型總應變與各元件之間的應變關系,即式(14),得所用的本構模型為

(22)

式(22)中：第1式適用范圍為在巖石蠕變過程未出現加速蠕變階段的過程；第2式適用于在巖石蠕變過程中出現加速蠕變階段的巖石蠕變。

5 模型參數辨識及驗證

對模型參數辨識及驗證主要采用1stOpt數學優化軟件,基于麥夸特法M-L算法和通用全局算法,利用不同條件下蠕變試驗數據對上述本構模型中的參數進行擬合,擬合參數如表2所示。通過分析本構模型中待定參數,若巖石蠕變過程未出現加速蠕變階段時,其所需擬合的參數量為5個,分別為E1、E2、η1、β、n；若巖石蠕變過程出現加速蠕變階段,其所需要擬合的參數量為8個,分別為E1、E2、η1、η2、β、n、b、λ。總體上,模型中所需擬合的參數相比較于現有的多數模型參數減少,同時根據R2,試驗數據與本構模型之間的相關性較高。

表2 硬性結構面蠕變本構模型參數

注:擬合參數小于1×10-10時取為0。

圖6 蠕變試驗值與預測值對比Fig.6 Comparisons between creep test values and predicted values

圖6(a)、圖6(b)所示為結構面連通率k=0.64、σ=8.40 MPa下剪切蠕變試驗與擬合值對比。對比非貫通結構面巖樣的剪切蠕變試驗值與預測值,本文中的本構模型能夠較好地對非貫通結構面剪切蠕變特性的3個階段進行模擬,之間的平均相關性系數R2為0.94。圖6(c)、圖6(d)分別為利用本文所建本構模型對文獻[1,17]中貫通型硬性結構面巖石的蠕變試驗數據進一步驗證,最終通過比較試驗值與擬合值,表明本構模型能夠較好地對蠕變的3個階段進行模擬,之間的平均相關性系數R2為0.97,說明本文所構建基于分數階微積分的巖石硬性結構面蠕變本構模型也可適用于貫通型硬性結構面。經過對所建本構模型的驗證,本構模型能夠較好地反映出含硬性結構面巖樣(巖體)的蠕變特性。

6 結論

通過元件組合方式構建適用于巖石含硬性結構面蠕變本構模型,在傳統Kelvin模型中引入基于分數階微積分的軟體元件,同時引入一種能夠描述結構面蠕變的加速元件,將加速蠕變元件的黏滯系數進行非定常化,構建基于分數階微積分的巖石硬性結構面蠕變本構模型。

通過連通率k=0.64、σ=8.40 MPa、不同剪應力下的剪切蠕變試驗數據對所構建的本構模型進行辨識,結果表明,構建的本構模型能夠較好地對非貫通硬性結構面巖體剪切蠕變特性的3個階段進行模擬,其平均R2為0.94左右。

利用含節理泥板巖及三峽船閘區巖石的貫通型硬性結構面試樣的蠕變試驗數據,進一步對本文所建的本構模型進行驗證,結果表明,本構模型能夠較好地對貫通性結構面的試驗數據進行辨識,其平均R2為0.97。

綜上,本文所構建基于分數階微積分的巖石硬性結構面蠕變本構模型能夠較好地反映含硬性結構面巖體的蠕變特性。