考慮時(shí)效損傷的巖石分?jǐn)?shù)階蠕變本構(gòu)模型

蘇艷軍

(中國建筑東北設(shè)計(jì)研究院有限公司，沈陽 110166)

1 研究背景

在與巖體相關(guān)的工程結(jié)構(gòu)領(lǐng)域中，往往會(huì)研究預(yù)測(cè)巖石的長期變形和強(qiáng)度，為工程結(jié)構(gòu)長期穩(wěn)定性的把控提供一定參考[1-2]。蠕變本構(gòu)模型作為一種不斷進(jìn)步和發(fā)展的材料變形破壞模擬手段，是巖石蠕變特性研究中的核心內(nèi)容[3-4]。

國內(nèi)外對(duì)于巖石蠕變本構(gòu)關(guān)系的研究已有較多成果，Moghadam等[5]以鹽巖為研究對(duì)象，建立黏彈塑性蠕變本構(gòu)模型，提出一個(gè)非線性黏滯系數(shù)的牛頓體，構(gòu)建一個(gè)新的非線性黏彈塑性蠕變模型；Günther等[6]探究鹽巖穩(wěn)態(tài)蠕變速率與溫度及應(yīng)力之間的聯(lián)系，構(gòu)建了不同溫度、應(yīng)力條件下的鹽巖蠕變模型；王艷春等[7]進(jìn)行頁巖在不同溫度和pH值下的蠕變?cè)囼?yàn)，建立同時(shí)考慮溫度和pH值的蠕變本構(gòu)模型；黃方勇[8]在西原模型的基礎(chǔ)上，引入分?jǐn)?shù)階微積分理論，建立改進(jìn)后的分?jǐn)?shù)階蠕變模型；張玉等[9]在傳統(tǒng)Burgers模型的基礎(chǔ)上，引入損傷力學(xué)理論，建立改進(jìn)Burgers蠕變損傷模型；金俊超等[10]通過探索巖石峰后應(yīng)變軟化指標(biāo)與塑性變形之間的關(guān)系，構(gòu)建一個(gè)基于應(yīng)變軟化指標(biāo)的黏塑性體，與Burgers模型串聯(lián)得到新的非線性蠕變模型；張亮亮等[11]在廣義Burgers模型的基礎(chǔ)上，串聯(lián)一個(gè)具有非定常特征的黏塑性元件，得到一個(gè)能較好描述砂巖蠕變特征的本構(gòu)模型。

本文將Maxwell模型和黏塑性體串聯(lián)，得到結(jié)構(gòu)形式簡(jiǎn)練的基礎(chǔ)模型。采用連續(xù)損傷理論，根據(jù)能量損傷方式定義損傷變量，構(gòu)建考慮時(shí)效損傷的彈性體來反映巖石彈性應(yīng)變。引入分?jǐn)?shù)階微積分理論，得到分?jǐn)?shù)階軟體元件來描述巖石黏彈性應(yīng)變。為表征黏塑性應(yīng)變，在分?jǐn)?shù)階軟體元件基礎(chǔ)上進(jìn)行損傷演化，得到考慮時(shí)效損傷的分?jǐn)?shù)階黏塑性體。將改進(jìn)后的元件代入基礎(chǔ)模型，得到考慮時(shí)效損傷的巖石分?jǐn)?shù)階蠕變本構(gòu)模型，給出參數(shù)解析方法，利用泥質(zhì)板巖蠕變數(shù)據(jù)驗(yàn)證模型合理性和優(yōu)越性。探究損傷發(fā)展過程，分析參數(shù)敏感度，并通過紅砂巖、千枚巖單軸壓縮各向異性蠕變特性試驗(yàn)數(shù)據(jù)來驗(yàn)證模型適用性。

2 基礎(chǔ)模型及考慮時(shí)效損傷的彈性體

2.1 巖石蠕變特征及基礎(chǔ)模型擇取

大量試驗(yàn)研究表明，不同應(yīng)力狀態(tài)下巖石蠕變特征可歸納為圖1中的3類曲線[12]，其中曲線①、②均只有衰減和穩(wěn)定蠕變2個(gè)階段，不同之處在于應(yīng)力σ與長期強(qiáng)度σp的相對(duì)大小。曲線①的應(yīng)力條件為σ≤σp，穩(wěn)定蠕變階段速率為0；曲線②的應(yīng)力條件為σ>σp，穩(wěn)定蠕變階段速率為某一常數(shù)。曲線③在曲線②的基礎(chǔ)上，還出現(xiàn)了加速蠕變階段。t1和t2分別為衰減、穩(wěn)定和加速蠕變階段的分界時(shí)間點(diǎn)，t3為加速蠕變階段止點(diǎn)。

圖1 蠕變特征示意圖Fig.1 Creep characteristic curves

從圖1可知：巖石內(nèi)部同時(shí)賦存著黏彈塑性變形特征，由此模型中須同時(shí)含有彈簧體、牛頓體和塑性體。目前多數(shù)改進(jìn)后的蠕變本構(gòu)模型盡管能取得較好的擬合效果，但是通常存在參數(shù)眾多、物理意義不明確的缺陷。本文將Maxwell模型和黏塑性元件串聯(lián)，得到一個(gè)結(jié)構(gòu)簡(jiǎn)單、形式簡(jiǎn)練的元件模型(圖2)，將該模型作為基礎(chǔ)模型進(jìn)行改進(jìn)，以期得到較為全面、簡(jiǎn)練的蠕變本構(gòu)模型。

圖2 元件模型示意圖Fig.2 Model components

該元件模型的蠕變方程為

(1)

式中：σ為應(yīng)力；E1為瞬時(shí)彈性模量；ηα1和ηα2分別為模型第2、第3部分的黏滯系數(shù)；α1、α2分別為模型第2、第3部分黏滯體的階數(shù)；ε為模型總應(yīng)變；t為時(shí)間；σp為長期強(qiáng)度。

2.2 損傷變量定義

巖石在外界應(yīng)力、溫度、風(fēng)化、水等因素的作用下，內(nèi)部微裂紋發(fā)育、擴(kuò)展和貫通，造成巖石性能劣化的微結(jié)構(gòu)變化即為損傷[9,13]。由此在圖2模型的基礎(chǔ)上引入連續(xù)損傷理論，目前主要有幾何損傷和能量損傷2種方式定義損傷變量[13]，前者需測(cè)定結(jié)構(gòu)有效承載面積，后者基于彈性模量變化規(guī)律?？紤]到可操作性，采取能量損傷方式定義損傷變量D(t)。

研究表明，巖石蠕變變形過程中的彈性模量和黏滯系數(shù)往往隨時(shí)間的增長而減小[13]，許宏發(fā)[14]進(jìn)行泥質(zhì)板巖單軸壓縮蠕變?cè)囼?yàn)，得到蠕變變形過程中不同時(shí)刻的彈性模量，如表1所示。

表1中彈性模量隨著時(shí)間增長而遞減，時(shí)間超過100 d以后，衰減幅度明顯放緩。依據(jù)表1，通過能量損傷的方法作如下定義：①當(dāng)t=0時(shí)，D(t)=

表1 不同時(shí)間的彈性模量Table 1 Elastic modulus with time

0；②當(dāng)t→∞時(shí)，D(t)=1；③當(dāng)t∈(0,∞)時(shí)，D(t)隨時(shí)間增長逐漸趨于1?？紤]到彈性模型的衰減規(guī)律，由此建立如下負(fù)指數(shù)型損傷演化方程，即

D(t)=1-e-bt。

(2)

式中b為與巖石損傷相關(guān)的參數(shù)。

2.3 考慮時(shí)效損傷的彈性體

假設(shè)巖石材料為各向同性損傷，則彈性模量隨時(shí)間的損傷劣化可表示為

E(t)=E1[1-D(t)] 。

(3)

式中E(t)為隨時(shí)間變化的彈性模量。

當(dāng)加載應(yīng)力大于長期強(qiáng)度時(shí)，巖石內(nèi)部損傷才開始發(fā)展累積[9,13]。無損狀態(tài)下的巖石對(duì)加載應(yīng)力的響應(yīng)是瞬時(shí)的，通過Hooke體描述彈性應(yīng)變，即

(4)

式中ε1為改進(jìn)后模型中第1部分的應(yīng)變。

當(dāng)加載應(yīng)力超過長期強(qiáng)度時(shí)，依據(jù)Lemaitre[15]等效應(yīng)變?cè)碛?/p>

(5)

式(5)即為考慮時(shí)效損傷的彈性體的本構(gòu)方程。取瞬時(shí)彈性模量E1=4.85 GPa，b=1×10-4，可得表1彈性模量衰減規(guī)律的擬合曲線，如圖3所示。

圖3 彈性模量擬合曲線Fig.3 Fitted curve of elastic modulus

由圖3可知，式(5)能較好地描述彈性模量衰減規(guī)律，將考慮時(shí)效損傷的彈性體替換基礎(chǔ)元件模型第1部分的彈簧體是可行的。

3 基于分?jǐn)?shù)階微積分的黏滯體和黏塑性體

3.1 Riemann-Liouville型分?jǐn)?shù)階微積分

分?jǐn)?shù)階微積分區(qū)別于傳統(tǒng)整數(shù)階微積分之處在于階數(shù)可為有理分?jǐn)?shù)、無理數(shù)甚至復(fù)數(shù)。分?jǐn)?shù)階微積分物理意義清晰，形式簡(jiǎn)練，能較好地描述非線性動(dòng)力行為，故將其引入來表征巖石蠕變的非線性非定常特征[16]。分?jǐn)?shù)階微積分根據(jù)時(shí)域的不同有Caputo型、Grunwald-Letnikov型、Riemann-Liouville型[16]。其中，Riemann-Liouville型相對(duì)簡(jiǎn)練，應(yīng)用較廣。函數(shù)f(t)在可積區(qū)間[0，t]的α(α>0)階Riemann-Liouville積分定義為

(6)

式中：f(t)為在可積區(qū)間[0，t]的某一函數(shù)；α為階數(shù)；s為用于Laplace變換的某一自變量；Γ(α)為Gamma函數(shù)，定義為

(7)

與積分相對(duì)應(yīng)，函數(shù)f(t)的α階微分定義為

(8)

式中n=[α]。

(9)

3.2 基于分?jǐn)?shù)階微積分的軟體元件

巖石材料是一種復(fù)雜的非線性材料，其性質(zhì)被認(rèn)為是介于理想固體和理想流體之間的某種狀態(tài)，根據(jù)牛頓體本構(gòu)方程的結(jié)構(gòu)特點(diǎn)，參考Blair[17]引入一種描述該狀態(tài)的軟體元件，其本構(gòu)方程為

(10)

式中：ε(t)為軟體元件的應(yīng)變；ηα為該軟體元件中的黏滯系數(shù)。當(dāng)α=0時(shí)，該軟體元件可描述理想固體；當(dāng)α=1時(shí)，該軟體元件退化為牛頓體；由此認(rèn)為當(dāng)α∈(0，1)時(shí)，該軟體元件可描述介于理想固體和理想流體之間的某種狀態(tài)。

當(dāng)應(yīng)力σ為恒定不變時(shí)，依據(jù)Riemann-Liouville分?jǐn)?shù)階微積分算子理論，對(duì)式(10)進(jìn)行分?jǐn)?shù)階積分得

(11)

式(11)即為基于分?jǐn)?shù)階微積分的軟體元件。

取σ=4 MPa，ηα1=0.5 GPa·d，依據(jù)式(11)得到不同α取值對(duì)應(yīng)的蠕變曲線，如圖4所示。

圖4 分?jǐn)?shù)階軟體元件蠕變曲線Fig.4 Creep curves of fractional software elements

由圖4可看出，當(dāng)α在區(qū)間[0，1]內(nèi)取不同量值，蠕變曲線表現(xiàn)出不同程度的非線性特征，隨著α值的減小，蠕變曲線的非線性特征愈加明顯。

3.3 分?jǐn)?shù)階黏滯體

圖4中蠕變曲線隨著α值的改變而不斷在某個(gè)范圍內(nèi)變化，能模擬圖1中的曲線①、②，由此說明該軟體元件能較好地反映巖石衰減、穩(wěn)定蠕變階段。鑒于該軟體元件具備描述巖石衰減、穩(wěn)定蠕變行為的靈活性，將該軟體元件替換基礎(chǔ)元件模型(圖2)中第2部分的牛頓體，利用其描述巖石蠕變黏彈性變形，于是有

(12)

式中ε2為改進(jìn)后模型中第2部分的應(yīng)變。

3.4 考慮時(shí)效損傷的分?jǐn)?shù)階黏塑性體

根據(jù)圖1曲線③可看出，在加速蠕變階段中曲線斜率不斷變化，巖石的黏滯系數(shù)將不再恒為常數(shù)，而軟體元件的黏滯系數(shù)是恒定的，故其無法描述巖石加速蠕變行為。因此采用連續(xù)損傷理論，依據(jù)Lemaitre[15]應(yīng)變等價(jià)性假說，通過損傷變量D(t)來描述黏塑性體黏滯系數(shù)的劣化，將軟體元件的本構(gòu)方程寫為式(13)。需注意的是，當(dāng)應(yīng)力水平超過了長期強(qiáng)度，損傷分?jǐn)?shù)階黏塑性體才會(huì)表現(xiàn)有黏塑性變形，由此應(yīng)將σ替換成(σ-σp)，即

(13)

式中ε3為改進(jìn)后模型中第3部分的應(yīng)變。

將式(2)代入式(13)可得

(14)

式中ηα2e-bt即為黏塑性體中隨時(shí)間劣化的黏滯系數(shù)。

當(dāng)σ為恒定不變時(shí)，依據(jù)Riemann-Liouville分?jǐn)?shù)階微積分算子理論，對(duì)式(14)進(jìn)行分?jǐn)?shù)階積分得

(15)

式中k表示數(shù)學(xué)求和符號(hào)的下界。式(15)即為考慮時(shí)效損傷的分?jǐn)?shù)階黏塑性體的本構(gòu)方程。

3.5 考慮時(shí)效損傷的巖石分?jǐn)?shù)階蠕變本構(gòu)模型

將分?jǐn)?shù)階黏滯體、考慮時(shí)效損傷的彈性體和分?jǐn)?shù)階黏塑性體串聯(lián)，參考基礎(chǔ)元件模型的組合方式得到

式(16)即為本文考慮時(shí)效損傷的巖石分?jǐn)?shù)階蠕變本構(gòu)模型。

4 模型驗(yàn)證及參數(shù)求解

4.1 參數(shù)求取方法

瞬時(shí)彈性模量E1取0 d時(shí)刻的4.85 GPa，式(16)第一式中的參數(shù)ηα1和α1通過一般的非線性最小二乘法便可擬合求解。由于式(16)第二式包含了∑符號(hào)，一般的非線性求解方法不適用，故采用雙參數(shù)的Mittag-Leffler函數(shù)[18]進(jìn)行模型表述，其中該函數(shù)表達(dá)式為

(17)

式中ω、n和z均為Mittag-Leffler函數(shù)的反演計(jì)算參數(shù)，其中ω>0、n>0、z為常數(shù)。

當(dāng)ω=1時(shí)，Mittag-Leffler函數(shù)的一般形式可以總結(jié)為[18]

(18)

將式(18)代入式(16)中的第2式可得當(dāng)σ>σp時(shí)的蠕變本構(gòu)方程，即

σ>σp。

(19)

粒子群算法(PSO)[19]隨機(jī)初始化每個(gè)粒子，每個(gè)粒子均表示一個(gè)可能性解，不斷更新粒子速度和位置，通過迭代評(píng)估每個(gè)粒子并得到全局最優(yōu)，基于MatLab編程實(shí)現(xiàn)粒子群算法，反演求解式(19)中最優(yōu)模型參數(shù)。

4.2 模型擬合及參數(shù)求解

為了驗(yàn)證本文所建模型的辨識(shí)效果，借助許宏發(fā)[14]泥質(zhì)板巖單軸壓縮蠕變?cè)囼?yàn)數(shù)據(jù)，由等時(shí)應(yīng)力-應(yīng)變曲線法[14]確定該巖石長期強(qiáng)度為11.50 MPa。通過本文所建模型進(jìn)行擬合，采用文獻(xiàn)[8]中未考慮時(shí)效損傷的分?jǐn)?shù)階蠕變模型進(jìn)行對(duì)比，得到擬合值與試驗(yàn)值對(duì)比曲線(將加載應(yīng)力標(biāo)識(shí)在曲線上)，如圖5所示。模型參數(shù)如表2所列。

圖5 不同加載應(yīng)力下應(yīng)變擬合值與試驗(yàn)值對(duì)比曲線Fig.5 Comparison between fitting values and test values of strain under varied loading stress

表2 模型參數(shù)Table 2 Model parameters

由圖5和表2可看出，引用的未考慮時(shí)效損傷的分?jǐn)?shù)階蠕變模型的辨識(shí)能力尚可，平均R2為0.967 2，但對(duì)穩(wěn)定蠕變階段快結(jié)束部分以及衰減蠕變階段的擬合存在較大偏差。本文所建模型擬合精度較高，平均R2達(dá)到了0.991 1，能較精準(zhǔn)地識(shí)別巖石蠕變曲線，尤其是加速蠕變階段。比較兩者模型，本質(zhì)上的不同在于本文所建模型有機(jī)結(jié)合了宏觀的蠕變變形和微細(xì)觀的損傷發(fā)展，不僅能較好地實(shí)現(xiàn)對(duì)蠕變量值的擬合，還可從“內(nèi)外結(jié)合”的角度模擬巖石蠕變變形全過程。

4.3 巖石損傷發(fā)展

當(dāng)應(yīng)力超過長期強(qiáng)度時(shí)，巖石材料內(nèi)部損傷才開始累積。泥質(zhì)板巖蠕變?cè)囼?yàn)中，在第4—第6級(jí)加載等級(jí)下，加載應(yīng)力超過了長期強(qiáng)度，從第4級(jí)加載起，損傷開始累積，巖石內(nèi)部微裂紋不斷發(fā)育、延展，至第6級(jí)加速蠕變階段末，微裂紋貫通，巖石出現(xiàn)宏觀變形破壞。本文以第4級(jí)加載起，根據(jù)式(2)和表2中模型參數(shù)，繪制損傷累積曲線，如圖6所示，利用玻爾茲曼[7-8]疊加處理第4—第6級(jí)蠕變曲線，得到分級(jí)加載形式，同時(shí)繪入圖6中。

圖6 損傷累積曲線Fig.6 Accumulated damage curves

由圖6可知，損傷曲線隨著時(shí)間增長而逐漸累積，至最后一級(jí)應(yīng)力加載結(jié)束，損傷變量趨近于1。在第4、第5級(jí)以及第5、第6級(jí)應(yīng)力加載之間，損傷變量出現(xiàn)了小幅的降低，分析其原因可能是應(yīng)力提升的過程中，巖石內(nèi)部微裂紋在外力增強(qiáng)作用下出現(xiàn)瞬間閉合現(xiàn)象，巖石應(yīng)變此時(shí)隨之小幅降低。待加載應(yīng)力穩(wěn)定后，巖石應(yīng)變不斷增長，損傷變量繼續(xù)累積。

4.4 參數(shù)敏感度分析

加載應(yīng)力的大小關(guān)系到損傷機(jī)制的觸發(fā)，分?jǐn)?shù)階黏塑性體的階數(shù)α2與巖石加速蠕變行為緊密關(guān)聯(lián)，故分析加載應(yīng)力和α2對(duì)蠕變發(fā)展的影響，判斷其敏感度。以泥質(zhì)板巖最后一級(jí)加載的模型參數(shù)為例，分別以加載應(yīng)力和α2作為因變量，保持其它參數(shù)不變，繪制參數(shù)敏感度曲線，如圖7所示。

圖7 參數(shù)敏感度曲線Fig.7 Curves of parameter sensitivity

分析圖7(a)，巖石應(yīng)變隨著加載應(yīng)力的提升而遞增，加載應(yīng)力影響整個(gè)蠕變過程的應(yīng)變，對(duì)衰減蠕變階段的曲線斜率有一定影響，但是不影響穩(wěn)定、加速蠕變階段的曲線斜率，由此可認(rèn)為加載應(yīng)力主要影響巖石衰減蠕變階段速率和整體應(yīng)變，故以最大應(yīng)變作為敏感度判斷依據(jù)。當(dāng)加載應(yīng)力從15.5 MPa提升到18.5 MPa時(shí)，最大應(yīng)變從0.016 6增加到0.028 6，加載應(yīng)力平均每1 MPa的增加引起了最大應(yīng)變0.40%的變化。

分析圖7(b)，巖石應(yīng)變隨著α2的增大而遞增，α2影響穩(wěn)定、加速蠕變階段的應(yīng)變和曲線斜率，幾乎不影響衰減蠕變階段的應(yīng)變和曲線斜率，由此認(rèn)為α2主要影響穩(wěn)定、加速蠕變階段的蠕變應(yīng)變及速率，故以蠕變速率(從46.96 d—202.03 d取近直線段斜率)和最大應(yīng)變和同時(shí)作為敏感度判斷依據(jù)。當(dāng)α2值從0.34增大到0.37時(shí)，最大應(yīng)變從0.016 6增加到0.023 5，蠕變速率從3.78×10-5d-1增加到7.49×10-5d-1，α2平均每0.01量綱的增加引起了最大應(yīng)變23.00%的變化，α2平均每0.01量綱的增加引起了蠕變速率123.67%的變化。

總體上，蠕變曲線隨著加載應(yīng)力和α2的改變而表現(xiàn)出一定規(guī)律性的變化，α2的敏感度高于加載應(yīng)力。

4.5 模型適用性

為驗(yàn)證本文所建模型模擬不同巖石蠕變行為(尤其是加速蠕變階段)的適用性，引用文獻(xiàn)[20-21]中紅砂巖和千枚巖單軸壓縮各向異性蠕變特性試驗(yàn)研究得到的相關(guān)蠕變?cè)囼?yàn)數(shù)據(jù)，分別通過本文所建模型和文獻(xiàn)[8]引用模型進(jìn)行擬合，得到擬合值與試驗(yàn)值對(duì)比曲線，如圖8所示。

圖8 不同巖石在不同加載應(yīng)力下的應(yīng)變擬合值與試驗(yàn)值對(duì)比曲線Fig.8 Comparison between fitting values and test values of strain under varied loading stress

由圖8可看出，本文所建模型識(shí)別效果較好，擬合精度較高，R2分別為0.992 4和0.989 3。引用模型的R2分別為0.975 5和0.953 6，對(duì)于穩(wěn)定、加速蠕變階段的曲線識(shí)別存在較大偏差。

本文所建考慮時(shí)效損傷的巖石分?jǐn)?shù)階蠕變本構(gòu)模型中彈性模量E1可通過試驗(yàn)確定，僅有ηα1、α1、ηα2、α2和b這5個(gè)參數(shù)需計(jì)算求解，模型結(jié)構(gòu)簡(jiǎn)單、形式簡(jiǎn)練、參數(shù)較少、物理意義較明確，具有較強(qiáng)的蠕變曲線識(shí)別能力，對(duì)不同類型巖石蠕變行為的描述具有較好的適用性。

5 結(jié) 論

(1)本文引入連續(xù)損傷理論，基于彈性模量隨時(shí)間衰減規(guī)律，根據(jù)能量損傷的方式以負(fù)指數(shù)型函數(shù)定義損傷變量，構(gòu)建考慮時(shí)效損傷的彈性體，進(jìn)行數(shù)據(jù)擬合驗(yàn)證該損傷演化方式的可行性。

(2)引入Riemann-Liouville型分?jǐn)?shù)階微積分算子理論，構(gòu)建具有非線性特征的分?jǐn)?shù)階軟體元件。利用該軟體元件作為分?jǐn)?shù)階黏滯體描述巖石黏彈性應(yīng)變，在該軟體元件的基礎(chǔ)上，結(jié)合損傷演化規(guī)律，建立考慮時(shí)效損傷的分?jǐn)?shù)階黏塑性體。

(3)聯(lián)合分?jǐn)?shù)階黏滯體、考慮時(shí)效損傷的彈性體和分?jǐn)?shù)階黏塑性體，建立一個(gè)新的考慮時(shí)效損傷的分?jǐn)?shù)階蠕變本構(gòu)模型。給出參數(shù)解析方法，利用泥質(zhì)板巖蠕變數(shù)據(jù)驗(yàn)證模型合理性和優(yōu)越性。損傷曲線隨著時(shí)間增長而逐漸累積，參數(shù)α2的敏感度高于加載應(yīng)力。通過紅砂巖、千枚巖單軸壓縮各向異性蠕變特性試驗(yàn)數(shù)據(jù)驗(yàn)證了模型的適用性。