考慮軟巖參數(shù)隨時(shí)間變化的非定常西原蠕變模型參數(shù)辨識(shí)及數(shù)值模擬

邴卿德, 孫騰閣, 張保良, 田忠喜, 張明賢

(1. 山東建筑大學(xué) 資產(chǎn)管理處, 山東 濟(jì)南 250101; 2. 濰坊醫(yī)學(xué)院 后勤管理處, 山東 濰坊 261053;3. 聊城大學(xué) 建筑工程學(xué)院, 山東 聊城 252000)

軟巖蠕變威脅施工和運(yùn)營(yíng)期地下工程的安全。 在軟巖中開挖洞室時(shí), 圍巖因蠕變變形而擠壓支護(hù)等結(jié)構(gòu), 嚴(yán)重時(shí)會(huì)造成襯砌開裂等破壞, 蠕變變形和斷裂的問題在巖土工程竣工后仍時(shí)有發(fā)生[1], 因此探明軟巖蠕變規(guī)律對(duì)保證工程安全具有重要意義。

軟巖蠕變變形表示方法通常分為經(jīng)驗(yàn)公式和元件模型[2-4]。經(jīng)驗(yàn)公式的適用性較低,而軟巖蠕變變形的元件模型方法以理想化的彈性元件、黏性元件、 塑性元件等為基本元件,將這些基本元件進(jìn)行串、 并聯(lián)組合,并根據(jù)基本元件的組合方式推導(dǎo)軟巖蠕變相應(yīng)的應(yīng)力-應(yīng)變關(guān)系方程。根據(jù)應(yīng)力-應(yīng)變關(guān)系方程,可以推導(dǎo)蠕變以及松弛方程,概念清晰且思路明確,在軟巖蠕變數(shù)值分析中應(yīng)用較廣泛[5]。Sterpi等[6]在Kelvin-Voigt模型基礎(chǔ)上串聯(lián)一個(gè)黏塑性元件,得到了非線性黏彈-黏塑性蠕變模型。顏海春等[7]基于Shvedov-Bingham模型研究了隧道圍巖塑性區(qū)的蠕變,分析不同位置的剪應(yīng)力和體積變化特征,應(yīng)用元件模型推導(dǎo)了應(yīng)力-應(yīng)變曲線及蠕變區(qū)半徑的解析解。張冬梅等[8]為了預(yù)測(cè)軟土隧道地表的長(zhǎng)期沉降,采用了彈性元件與Kelvin元件串聯(lián)的三參量黏彈性模型,結(jié)果表明,軟土隧道的沉降具有持續(xù)時(shí)間長(zhǎng)、 沉降量大的特點(diǎn)。趙旭峰等[9]采用三參量廣義Kelvin模型對(duì)廈門東通道海底隧道三軸壓縮蠕變?cè)囼?yàn)數(shù)據(jù)進(jìn)行擬合,得到了較準(zhǔn)確的蠕變參數(shù)。閻巖等[10]基于西原蠕變模型研究各蠕變參數(shù)與應(yīng)力、 時(shí)間的關(guān)系,得到了變參數(shù)的蠕變方程,但是該方程未考慮彈性模量隨應(yīng)力的變化規(guī)律。潘曉明等[11]將西原模型引入ABAQUS軟件中進(jìn)行二次開發(fā),經(jīng)過與單軸壓縮蠕變?cè)囼?yàn)結(jié)果對(duì)比發(fā)現(xiàn),蠕變模型程序較準(zhǔn)確。

大多數(shù)研究假定軟巖相關(guān)參數(shù)為定值,而在實(shí)際中,隨著蠕變的產(chǎn)生,軟巖內(nèi)部結(jié)構(gòu)會(huì)改變,軟巖蠕變的參數(shù)也會(huì)改變。為了更準(zhǔn)確地描述軟巖蠕變規(guī)律,應(yīng)注意軟巖各項(xiàng)參數(shù)隨應(yīng)力及時(shí)間的變化,并且不同類型軟巖的蠕變也有區(qū)別。

本文中對(duì)經(jīng)典西原模型參數(shù)進(jìn)行非定?；?建立一種考慮軟巖參數(shù)隨時(shí)間變化的非定常西原蠕變模型(簡(jiǎn)稱本文模型),并通過單軸壓縮蠕變?cè)囼?yàn)研究軟巖蠕變規(guī)律,基于單軸壓縮蠕變?cè)囼?yàn)結(jié)果,對(duì)本文模型進(jìn)行參數(shù)辨識(shí)。通過ABAQUS軟件的二次開發(fā)功能,將本文模型嵌入其中進(jìn)行軟巖蠕變過程的模擬驗(yàn)證。

1 模型建立

1.1 經(jīng)典西原模型

經(jīng)典西原模型由彈性元件、 Kelvin元件和黏塑性元件串聯(lián)而成,如圖1所示。該模型由E1、E2、η1、η2、σs這5個(gè)參數(shù)進(jìn)行描述,其中E1、E2分別為彈性元件、 Kelvin元件的彈性模量,η1、η2分別為黏塑性元件、 Kelvin元件的黏性系數(shù),σs為黏塑性元件的屈服強(qiáng)度。

E1、 E2—彈性元件、 Kelvin元件的彈性模量;η1、 η2—黏塑性元件、 Kelvin元件的黏性系數(shù);σs—黏塑性元件的屈服強(qiáng)度; σ—軸向應(yīng)力。

經(jīng)典西原模型根據(jù)受力大小分為2種模型。 當(dāng)軸向應(yīng)力σ<σs時(shí), 經(jīng)典西原模型中的黏塑性元件為剛體, 不會(huì)發(fā)生變形, 經(jīng)典西原模型變?yōu)閺V義Kelvin模型, 本構(gòu)方程[2]為

(1)

蠕變方程[2]為

(2)

當(dāng)σ≥σs時(shí),經(jīng)典西原模型發(fā)生不穩(wěn)定蠕變,經(jīng)典西原模型變?yōu)椴袼鼓Ｐ?本構(gòu)方程[2]為

(3)

蠕變方程[2]為

(4)

1.2 本文模型

經(jīng)典西原模型中彈性元件的彈性模量E與黏性元件的黏性系數(shù)η均為定值, 但是軟巖的性質(zhì)會(huì)隨著蠕變發(fā)生變化, 定常的經(jīng)典西原模型無法準(zhǔn)確描述軟巖的非線性加速蠕變規(guī)律。 為了使經(jīng)典西原蠕變模型更準(zhǔn)確地描述軟巖蠕變規(guī)律, 本文中采用隨時(shí)間變化的彈性元件的彈性模量E1(t)以及黏塑性元件、 Kelvin元件的黏性系數(shù)η1(t)、η2(t)替換定常元件的E1、η1、η2, 建立本文模型, 如圖2所示。

E1(t)—隨時(shí)間t變化的彈性元件的彈性模量;E2—Kelvin元件的彈性模量; η1(t)、 η2(t)—隨時(shí)間t變化的黏塑性元件、 Kelvin元件的黏性系數(shù);σs—黏塑性元件的屈服強(qiáng)度; σ—軸向應(yīng)力。

1.2.1 非定常彈性元件

考慮經(jīng)典西原模型中的彈性元件隨時(shí)間而變化,該元件變化代表軟巖的瞬時(shí)變形。當(dāng)σ<σs時(shí),軟巖不發(fā)生塑性損傷,軟巖的彈性模量為常數(shù),此時(shí)彈性元件的變形為定值。當(dāng)σ≥σs時(shí),軟巖進(jìn)入屈服狀態(tài),彈性模量減小,軟巖的彈性變形逐漸增大。為了得出軟巖發(fā)生蠕變后彈性模量的變化規(guī)律,對(duì)經(jīng)歷不同軸向應(yīng)力、蠕變時(shí)間條件下未破壞的軟巖進(jìn)行單軸壓縮蠕變?cè)囼?yàn),測(cè)量彈性模量。軟巖蠕變后軟巖沒有明顯的彈性階段,國(guó)際巖石力學(xué)學(xué)會(huì)給出了3種定義非線性巖石彈性模量的計(jì)算方法,本文中采用割線模量作為軟巖的彈性模量,即取軸向應(yīng)力σ=σ50/2(其中σ50為軟巖的抗壓強(qiáng)度)處的割線模量作為蠕變后軟巖的彈性模量,即

E=σ/ε50,

(5)

式中ε50為軸向應(yīng)力為σ50/2時(shí)的應(yīng)變。

表1所示為不同軸向應(yīng)力和蠕變時(shí)間時(shí)軟巖的彈性模量。

表1 不同軸向應(yīng)力和蠕變時(shí)間時(shí)軟巖的彈性模量

當(dāng)軸向應(yīng)力大于軟巖屈服強(qiáng)度時(shí)才會(huì)發(fā)生蠕變,可以推斷彈性模量的減小與軸向應(yīng)力與屈服強(qiáng)度之差有關(guān),并且隨著蠕變時(shí)間的延長(zhǎng),軟巖內(nèi)部裂縫發(fā)展程度加重,彈性模量也會(huì)減小,即彈性模量與蠕變時(shí)間相關(guān),因此彈性模量受軸向應(yīng)力與屈服強(qiáng)度之差σ-σs和蠕變時(shí)間t共同影響,即

E(t)=f[(σ-σs)t]。

(6)

利用單軸壓縮蠕變?cè)囼?yàn)得到軟巖的應(yīng)力-應(yīng)變曲線,如圖3所示。從圖中可以看出,當(dāng)軸向應(yīng)力約為9 MPa時(shí),軟巖的徑向應(yīng)變突然增大,原因是此時(shí)軟巖發(fā)生屈服,軟巖內(nèi)部出現(xiàn)微裂縫,因此可以判斷軟巖的屈服強(qiáng)度為9 MPa。

根據(jù)式(6)及表1可得不同(σ-σs)t對(duì)應(yīng)的軟巖的彈性模量,如表2所示,其中當(dāng)(σ-σs)t為0時(shí),彈性模量為初始彈性模量。

表2 不同(σ-σs)t對(duì)應(yīng)的軟巖的彈性模量

圖4所示為彈性模量隨(σ-σs)t的變化。從圖中可以推導(dǎo)得出

E(t)=A+Bexp[-(σ-σs)t/C],

(7)

σ—軸向應(yīng)力; σs—屈服強(qiáng)度; t—蠕變時(shí)間。

式中A、B、C為待定系數(shù)。

對(duì)圖4中的數(shù)據(jù)點(diǎn)進(jìn)行擬合,得到

E(t)=2.68+0.861exp[-(σ-σs)t/185.4] 。

(8)

由式(8)可知,當(dāng)蠕變時(shí)間趨于無窮大時(shí),軟巖的彈性模量接近2.68 GPa,即該軟巖的極限彈性模量為2.68 GPa,因此非定常彈性元件的本構(gòu)方程為

σ=E1(t)ε1,

(9)

蠕變方程為

(10)

式中ε1為非定常彈性元件的應(yīng)變。

1.2.2 非定常Kelvin元件

黏性系數(shù)的衰減符合指數(shù)形式[12-13],即

η(t)=η0exp(-λt),

(11)

式中η0、λ為待定參數(shù)。非定常Kelvin元件的本構(gòu)

方程為

表2為改善因子a對(duì)最優(yōu)值的影響。結(jié)果顯示，當(dāng)其他參數(shù)不變時(shí)，隨著a的增加，N*顯著增加，m*基本保持不變，全周期期望維修費(fèi)率逐漸變小。這是因?yàn)閷?duì)MC實(shí)施的不完全維修效果越好，則只依靠不完全維修就能夠更好的維持系統(tǒng)運(yùn)行，替換維修的需求就越小，維修費(fèi)率也就會(huì)變小。

(12)

利用分離變量法求定積分

(13)

整理得到非定常Kelvin元件蠕變方程為

(14)

式中:ε2為非定常Kelvin元件的應(yīng)變;λ1為待定參數(shù)。

1.2.3 非定常黏塑性元件

非定常黏塑性元件的本構(gòu)方程[14-15]為

(15)

求解微分方程得到蠕變方程為

(16)

式中λ2、η20為待定參數(shù)。

ε=ε1+ε2+ε3=

(17)

2 參數(shù)辨識(shí)

基于軟巖全自動(dòng)單軸蠕變伺服設(shè)備進(jìn)行軟巖單軸壓縮蠕變?cè)囼?yàn),研究破碎帶軟巖蠕變力學(xué)特性,對(duì)不同條件下軸向蠕變規(guī)律進(jìn)行探討。為了保證巖樣性質(zhì)和天然含水量不受外界環(huán)境的影響,將制得的標(biāo)準(zhǔn)巖樣進(jìn)行蠟封保存,巖樣的底面直徑、 高度分別為50、 100 mm,手工磨制軟巖蠕變巖樣如圖5所示。單軸壓縮蠕變?cè)囼?yàn)是巖樣在不同軸向應(yīng)力作用下的單軸壓縮蠕變?cè)囼?yàn),需要在恒溫、 恒濕條件下進(jìn)行。 加載速率為0.375 MPa/min, 變形穩(wěn)定后保持恒定軸向應(yīng)力, 直至設(shè)定蠕變時(shí)間。 在軸向應(yīng)力加載階段首先進(jìn)行預(yù)壓,試驗(yàn)機(jī)的壓頭與巖樣上端面接觸;然后逐漸增大軸向應(yīng)力直至設(shè)定應(yīng)力并保持不變,每隔1 h記錄應(yīng)變隨時(shí)間的變化規(guī)律;最后,通過單軸壓縮蠕變?cè)囼?yàn)數(shù)據(jù)得出軟巖蠕變力學(xué)特性。

(a)巖樣1 (b)巖樣2 (c)巖樣3 (d)巖樣4

圖6所示為不同軸向應(yīng)力時(shí)軟巖的蠕變規(guī)律。從圖中可以看出: 隨著施加軸向應(yīng)力的增大,軟巖蠕變變形增加,初始衰減蠕變持續(xù)時(shí)間和變形進(jìn)入穩(wěn)態(tài)蠕變時(shí)間延長(zhǎng)。當(dāng)軸向應(yīng)力為8、 9 MPa時(shí),衰減蠕變階段之后,蠕變曲線基本為水平直線,穩(wěn)態(tài)蠕變階段速率趨于0,蠕變變形不明顯。當(dāng)軸向應(yīng)力大于10 MPa時(shí),軟巖出現(xiàn)加速蠕變階段,并且軸向應(yīng)力越大,加速蠕變出現(xiàn)越早。

圖6 不同軸向應(yīng)力時(shí)軟巖的蠕變規(guī)律

利用三軸壓縮試驗(yàn)測(cè)得軟巖的物理力學(xué)參數(shù),如表3所示。當(dāng)軸向應(yīng)力大于σs時(shí),黏塑性元件才會(huì)發(fā)生塑性蠕變,即軟巖進(jìn)入屈服狀態(tài),因此σs取值為軟巖的屈服強(qiáng)度,即σs=9 MPa。

采用本文模型對(duì)不同軸向應(yīng)力時(shí)的單軸壓縮蠕變?cè)囼?yàn)數(shù)據(jù)進(jìn)行擬合。 圖7所示為不同軸向應(yīng)力

表3 利用三軸壓縮試驗(yàn)測(cè)得軟巖的物理力學(xué)參數(shù)

(a)軸向應(yīng)力為8 MPa

(c)軸向應(yīng)力為10 MPa

(d)軸向應(yīng)力為12 MPa

時(shí)本文模型擬合曲線與單軸壓縮蠕變?cè)囼?yàn)數(shù)據(jù)。 從圖中可以看出, 擬合曲線與單軸壓縮蠕變?cè)囼?yàn)數(shù)據(jù)基本吻合。 不同軸向應(yīng)力時(shí)本文模型的擬合結(jié)果如表4所示。 由表可知, 本文模型擬合曲線的相關(guān)系數(shù)大于或等于0.997, 表明本文模型擬合曲線與單軸壓縮蠕變?cè)囼?yàn)數(shù)據(jù)擬合精度較高。 由于不同類型軟巖對(duì)應(yīng)的非定常參數(shù)不同, 因此對(duì)于其他類型的軟巖可以通過單軸壓縮蠕變?cè)囼?yàn)擬合得出。

表4 不同軸向應(yīng)力時(shí)考慮軟巖參數(shù)隨時(shí)間變化的非定常西原蠕變模型的擬合結(jié)果

3 數(shù)值模擬

本文模型直接應(yīng)用于實(shí)際工程中較困難,因此借助數(shù)值軟件實(shí)現(xiàn)本文模型的工程應(yīng)用。ABAQUS軟件不僅能提供標(biāo)準(zhǔn)的有限元分析程序,而且具有良好的開放性。使用ABAQUS軟件的用戶子程序接口編輯非標(biāo)準(zhǔn)的分析程序,在實(shí)際工程中進(jìn)行廣泛應(yīng)用。本文中利用CREEP子程序?qū)⒈疚哪Ｐ投x到ABAQUS軟件中,對(duì)單軸壓縮蠕變進(jìn)行數(shù)值模擬,驗(yàn)證本文模型在數(shù)值軟件中應(yīng)用的可行性。

依據(jù)本文模型建立ABAQUS軟件中的計(jì)算單元,如圖8所示。計(jì)算單元為高度是100 mm且底面直徑是50 mm的圓柱體,底面為固定約束,其余面為自由面,網(wǎng)格為六面體單元。計(jì)算單元的材料參數(shù)取值見表3、 4。通過在計(jì)算單元頂面施加軸向恒定荷載模擬軟巖單軸壓縮蠕變?cè)囼?yàn)。

圖8 計(jì)算單元尺寸及邊界

圖9所示為不同軸向應(yīng)力時(shí)計(jì)算單元應(yīng)變模擬結(jié)果與單軸壓縮蠕變?cè)囼?yàn)數(shù)據(jù)。由圖可知:當(dāng)軸向應(yīng)力為9 MPa時(shí),巖樣發(fā)生等速蠕變;當(dāng)軸向應(yīng)力為12 MPa時(shí),巖樣發(fā)生加速蠕變。當(dāng)軸向應(yīng)力為9、 12 MPa時(shí),蠕變各階段的模擬結(jié)果與試驗(yàn)結(jié)果基本吻合。說明本文模型的蠕變方程可以在數(shù)值軟件中準(zhǔn)確應(yīng)用,并進(jìn)行相關(guān)工程計(jì)算。

(a)軸向應(yīng)力為9 MPa

(b)軸向應(yīng)力為12 MPa

4 結(jié)論

基于經(jīng)典西原模型,通過對(duì)彈性元件、 Kelvin元件、黏塑性元件的相關(guān)參數(shù)進(jìn)行非定?；?建立考慮軟巖參數(shù)隨時(shí)間變化的本文模型,并采用單軸壓縮蠕變?cè)囼?yàn)對(duì)本文模型進(jìn)行驗(yàn)證與參數(shù)辨識(shí),基于ABAQUS軟件,實(shí)現(xiàn)了軟巖蠕變的模擬,得出以下主要結(jié)論:

1)采用單軸壓縮蠕變?cè)囼?yàn)測(cè)得了不同軸向應(yīng)力水平和蠕變時(shí)間時(shí)軟巖的彈性模量,并推導(dǎo)彈性模量隨(σ-σs)t的變化公式。

2)采用本文模型對(duì)不同軸向應(yīng)力時(shí)的軟巖單軸壓縮蠕變?cè)囼?yàn)進(jìn)行擬合及參數(shù)辨識(shí),擬合曲線的相關(guān)系數(shù)不小于0.997,表明本文模型擬合精度較高,并且對(duì)軟巖蠕變的初期蠕變、等速蠕變、加速蠕變3個(gè)階段描述較準(zhǔn)確,驗(yàn)證了本文模型的準(zhǔn)確性。

3)利用CREEP子程序?qū)⒈疚哪Ｐ颓度階BAQUS軟件中,針對(duì)單軸壓縮蠕變過程進(jìn)行數(shù)值模擬,蠕變各階段的模擬結(jié)果與試驗(yàn)結(jié)果基本吻合,驗(yàn)證了本文模型在數(shù)值軟件中應(yīng)用的可行性。