基于伯格模型的非飽和土蠕變模型構(gòu)建

鄭 俊，王世梅3，周 輝，王 力3，鄒良超3，許匯源4，侯會明

(1.中國科學(xué)院武漢巖土力學(xué)研究所 巖土力學(xué)與工程國家重點試驗室，武漢 430071； 2.中國科學(xué)院大學(xué)，北京 100049； 3.三峽大學(xué) 三峽庫區(qū)地質(zhì)災(zāi)害教育部重點實驗室，湖北 宜昌 443002； 4.湖北省地質(zhì)局第七地質(zhì)大隊，湖北 宜昌 443100)

1 研究背景

滑坡是一種常見的地質(zhì)災(zāi)害，降雨、庫水位變動以及其耦合作用下會導(dǎo)致滑坡變形增加，甚至誘發(fā)滑坡失穩(wěn)破壞，滑坡變形發(fā)展一般具有一定的時間過程[1-4]。如美國Grand Coulee水庫在1941年蓄水后的12 a內(nèi)，先后發(fā)生滑坡500起，其中49%的滑坡發(fā)生在蓄水后2 a內(nèi)，51%發(fā)生在2～12 a之間[5]。據(jù)三峽庫區(qū)滑坡的監(jiān)測資料分析，大多數(shù)滑坡的位移都為一時間序列，均表現(xiàn)出明顯的蠕變特征[2-4]。由此可知，大多數(shù)滑坡變形具有蠕滑性質(zhì)，即其變形破壞具有時間效應(yīng)，是關(guān)于時間的函數(shù)。降雨以及庫水位變動對滑坡土體的影響較大，導(dǎo)致其處于飽和-非飽和狀態(tài)，因此滑坡土體變形和穩(wěn)定性與土體的飽和-非飽和狀態(tài)有關(guān)聯(lián)[6]。因此，可以認(rèn)為絕大多數(shù)滑坡的失穩(wěn)過程實質(zhì)為非飽和蠕變過程，考慮基質(zhì)吸力的非飽和土蠕變特性對于滑坡長期穩(wěn)定性研究極其重要。

目前已有不少國內(nèi)外專家學(xué)者進行了大量關(guān)于飽和土蠕變特性研究，并在試驗基礎(chǔ)上結(jié)合一定的理論建立了相關(guān)的蠕變模型[7-11]。相對于飽和土，非飽和土蠕變特性研究的成果較少。祝艷波等[12]為預(yù)測基質(zhì)吸力對巴東組碎屑土夾層蠕變行為的影響，開展了非飽和蠕變經(jīng)驗?zāi)Ｐ偷难芯俊riol等[13]在三軸蠕變試驗的基礎(chǔ)上結(jié)合Alonso等[14]提出的非飽和土彈塑性模型建立了考慮基質(zhì)吸力的非飽和土黏彈塑性模型。賴小玲等[15]在滑動帶土非飽和蠕變試驗基礎(chǔ)上，建立了各級吸力水平下滑動帶土的Mesri蠕變模型。

蠕變模型大致可以分為3類[16]，即經(jīng)驗-半經(jīng)驗?zāi)Ｐ?、基于一般流變理論的元件模型和黏彈塑性模型。?jīng)驗-半經(jīng)驗?zāi)Ｐ鸵话悴捎冒雽?shù)或雙對數(shù)坐標(biāo)系下應(yīng)變或應(yīng)變率與時間的線性關(guān)系來表示[17-20]?；谝话懔髯兝碚摰哪Ｐ椭饕ㄟ^元件組合，即從基本力學(xué)概念與理論出發(fā)，運用數(shù)學(xué)力學(xué)分析方法，用模型元件(牛頓黏性體N、虎克彈性體H、圣維南塑性體S)的組合來模擬土體的蠕變特性[21-23]。通過這類模型或模型組合形成的復(fù)合流變模型可以不同程度地模擬巖土流變特性。黏彈塑性模型則是在經(jīng)典塑性理論上發(fā)展起來的，認(rèn)為塑性勢和蠕變勢一致，從而建立黏彈塑性模型[24-26]。本文主要結(jié)合蠕變試驗數(shù)據(jù)特點，采用基于一般流變理論的元件模型進行組合，以構(gòu)建合適的蠕變模型。

本文基于三峽庫區(qū)某滑坡滑帶土不同偏應(yīng)力水平、不同基質(zhì)吸力條件下的三軸蠕變試驗結(jié)果[27-28]，根據(jù)元件模型理論和非飽和土力學(xué)原理以及相關(guān)假設(shè)條件，類比飽和蠕變模型的構(gòu)建思想，構(gòu)建非飽和狀態(tài)下的擴展伯格蠕變模型；并根據(jù)蠕變數(shù)據(jù)中應(yīng)力-應(yīng)變關(guān)系呈現(xiàn)的非線性特征，借鑒雙曲線函數(shù)對構(gòu)建的非飽和狀態(tài)下的擴展伯格蠕變模型進行改進；并利用試驗數(shù)據(jù)對模型進行驗證，最后對模型進一步修正，并對修正后模型進行驗證。

2 非飽和蠕變模型構(gòu)建

本文采用擴展的伯格蠕變模型是在伯格模型的基礎(chǔ)上再串聯(lián)一個Kelvin體，也稱為M-2K模型，模型元件圖如圖1所示。

圖1 擴展的伯格體模型(M-2K)Fig.1 Extended Burger’s model(M-2K)

M-2K模型本構(gòu)方程為[29]

其中：

p1=El1El2，q1=El1El2β1，

q2=El1β1β3+El2β1β2，

對式(1)按應(yīng)變求解，可得到該模型的蠕變方程為[29]

(2)

元件理論中的元件都是一維的，為了實際工程應(yīng)用，可采用類似于彈性理論中將一維Hooke定理推廣到三維的方法，將元件模型推廣到三維，因此作如下假定[30]：

(1)材料的體積變形在受力瞬時完成，且體積變形為彈性，不隨時間發(fā)生變化。

(2)只有偏應(yīng)力張量引起蠕變，球應(yīng)力張量作用下的材料不發(fā)生蠕變。

(3)蠕變過程中，泊松比不隨時間發(fā)生變化。

在常規(guī)三軸蠕變試驗中，根據(jù)其受力狀態(tài)，軸向應(yīng)力為σ1，圍壓σ2=σ3，擴展的伯格蠕變模型的三維形式為[30]

(3)

其中：

σI=σ1-σ3， 2σII=σ1+2σ3。

一點的應(yīng)力狀態(tài)Sij可分解為偏應(yīng)力張量S″ij和球應(yīng)力張量S′ij之和[31]，即

(4)

式中：σm，δij分別為平均應(yīng)力和Kronecker符號；σkk為主應(yīng)力。

同時，在剪切應(yīng)力狀態(tài)下，其應(yīng)力狀態(tài)可采用3個主應(yīng)力σ1，σ2，σ3來表示，即

(5)

Fredlund等[32]認(rèn)為孔隙氣壓張量ua應(yīng)該被添加用于描述非飽和土中一點的應(yīng)力狀態(tài)，即雙應(yīng)力變量為基質(zhì)吸力張量和凈應(yīng)力張量。其中基質(zhì)吸力張量可表示為

s=ua-uw。

(6)

式中：s為基質(zhì)吸力張量；ua為孔隙氣壓力張量；uw為孔隙水壓力張量。

凈應(yīng)力張量可表示為

(7)

根據(jù)上述公式，式(3)可改寫為

(8)

式(8)即為非飽和狀態(tài)下的擴展伯格蠕變模型。

3 模型的改進

由于伯格蠕變模型為線彈黏性模型，其應(yīng)力-應(yīng)變關(guān)系為線性，但是從蠕變試驗數(shù)據(jù)可以看出，其應(yīng)力-應(yīng)變關(guān)系表現(xiàn)出明顯的非線性特征，因此為了將上述模型進一步拓展為非線性彈黏性模型，可借鑒經(jīng)驗蠕變模型中用來描述應(yīng)力-應(yīng)變關(guān)系的雙曲線函數(shù)對上述模型進行改進。

雙曲線型應(yīng)力-應(yīng)變方程最早由Kondner[33](1963)提出，用來模擬土體在常速率軸向變形條件下的應(yīng)力-應(yīng)變特性，該等軸雙曲線可以寫為

(9)

由式(9)可以得出初始切線模量Eu為

(10)

其最終主應(yīng)力差(σ1-σ3)ult為

(11)

可以看出，當(dāng)應(yīng)變?yōu)闊o窮大時，雙曲線型應(yīng)力-應(yīng)變曲線才能達到最大主應(yīng)力差(σ1-σ3)ult，然而，試驗過程中土體實際破壞剪應(yīng)力(σ1-σ3)f往往在有限應(yīng)變εf狀態(tài)下即可達到。為了使雙曲線通過試驗中的破壞點[εf,(σ1-σ3)f]，特引入破壞比Rf[33]，其表達式為

(12)

將式(10)—式(12)代入式(9)得到

(13)

式中：c為黏聚力；φ為內(nèi)摩擦角。

在非飽和土中，基于Fredlund等[32]提出的非飽和抗剪強度理論，其抗剪強度凈應(yīng)力和吸力ua-uw作為應(yīng)力變量可表示為

τf=c′+(σ-ua)tanφ′+(ua-uw)tanφb。 (15)

式中φb為隨吸力ua-uw變量變化的內(nèi)摩擦角。

將上述雙曲線型應(yīng)力-應(yīng)變關(guān)系(式(13))引入非飽和擴展伯格蠕變模型(式(8))中，則改進的模型形式為

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

4.1 參數(shù)求解

由于應(yīng)力-應(yīng)變函數(shù)關(guān)系式中的參數(shù)都有具體的物理意義，其參數(shù)在式(16)中單獨求取。具體步驟如下[34]。

對式(13)進行如下變換：

(17)

(18)

(19)

由于擴展的伯格蠕變模型中含有負指數(shù)項，其參數(shù)一般需要采用非線性最小二乘法進行回歸求得；而在利用非線性最小二乘法進行回歸時，由于回歸結(jié)果存在不確定性，且本文中所建模型的參數(shù)都具有較明確的物理含義，因此可根據(jù)其物理意義及蠕變曲線特征先求一組接近實際值的初值，然后利用非線性回歸算法進行回歸求得模型參數(shù)，以避免回歸過程中參數(shù)的跳躍。擴展的伯格蠕變模型參數(shù)初值的求解方法如下：

εij=A+a[Bt+C(1-e-Dt)+E(1-e-Ft)] 。(20)

令t=0，公式右端等于A，即可以根據(jù)瞬時變形來確定A，當(dāng)時間足夠長時，右端后兩式項a[C(1-e-Dt)+E(1-e-Ft)]趨于常數(shù)a(C+E)，則公式可看作直線方程，直線斜率即為aB。選定初值時可近似假定C=E，D=F，即可把公式右端后兩項a[C(1-e-Dt)+E(1-e-Ft)]看作2aC(1-e-Dt)，在減速蠕變階段，式(20)可簡化為

εij=A+2aC(1-e-Dt) 。

(21)

將式(21)中常數(shù)移至等式左邊，在等式兩邊取對數(shù)，即ln(εij-A-2aC)=ln(-2aC)-Dt。

上式在對數(shù)坐標(biāo)系中為一直線方程，A已確定，D即為直線的斜率，再將D值代入式(21)即可解出C。由此可確定一組初值，再利用MatLab中l(wèi)sqcurvefit非線性回歸工具回歸出A，B，C，D，E，F(xiàn)的穩(wěn)定解。由于μ值對其他參數(shù)的影響很小[29]，本文中根據(jù)經(jīng)驗，μ值取 0.4。

圖2 不同凈圍壓條件下蠕變關(guān)系曲線Fig.2 Curves of under different net confiningpressures

表2 各組圍壓下模型參數(shù)及相關(guān)系數(shù)Table 2 Model parameters and correlation coefficient under different confining pressures

表1 應(yīng)力-應(yīng)變雙曲線函數(shù)參數(shù)Table 1 Parameters of stress-strain hyperbolic function

將所求得的應(yīng)力-應(yīng)變雙曲線函數(shù)參數(shù)代入式(21)，按上述的方法求取參數(shù)初值，再利用非線性最小二乘法回歸，即可以得到模型所有參數(shù)。求解各組凈圍壓下蠕變模型的參數(shù)，結(jié)果如表2所示。

4.2 模型檢驗

為了檢驗?zāi)Ｐ偷臄M合效果，將以上參數(shù)代入蠕變模型，并對試驗數(shù)據(jù)進行模擬，求取模擬結(jié)果與試驗數(shù)據(jù)的擬合系數(shù)R2(詳見表2)，并繪制蠕變曲線與實際曲線進行對比，如圖3所示。

圖3 不同凈圍壓條件下時蠕變試驗值與模擬值對比Fig.3 Comparison of creep between simulation andtest under different net confining pressures

從表2和圖3可以看出，上述蠕變模型的模擬效果比較好，與試驗數(shù)據(jù)的相關(guān)系數(shù)基本都在0.96以上，各級凈圍壓下，模擬蠕變曲線與實際試驗值比較吻合，說明上述蠕變模型可以較好地模擬滑帶土的蠕變特征。

4.3 模型的進一步修正與驗證

考慮到模型參數(shù)與凈圍壓之間的線性相關(guān)性，可將其線性函數(shù)關(guān)系(式(22))代入原模型(式(16))，對模型進行修正，修正后的模型如式(23)。

(22)

(23)

圖4 蠕變試驗值與模擬值對比Fig.4 Comparison of creep between test and simulation

為了進一步驗證改進后的非飽和蠕變模型，選取另外一組試驗數(shù)據(jù)與模型預(yù)測值進行比較。該組試驗數(shù)據(jù)為σ′3=200 kPa，σ3=400 kPa，ua=200 kPa，偏應(yīng)力水平分別為0.18，0.35，0.50，0.65，0.80。試驗過程中所選的土樣為同一批土樣，忽略制樣差異，模型參數(shù)則可采用上述參數(shù)。繪制模擬曲線與試驗曲線，如圖4所示。

從圖4可以看出，蠕變試驗曲線與模型模擬值雖然不完全吻合，但是趨勢基本一致，說明本文所建立的非飽和土的蠕變模型合理有效，基本反映了非飽和土的蠕變特征。

5 結(jié) 論

(1)根據(jù)試驗蠕變曲線特征，選取元件模型中的伯格蠕變模型為基礎(chǔ)模型，并基于元件模型相關(guān)理論和非飽和土力學(xué)原理以及相關(guān)假設(shè)，推導(dǎo)出了適用于非飽和土的擴展的伯格蠕變模型的三維形式。根據(jù)蠕變曲線中應(yīng)力-應(yīng)變關(guān)系表現(xiàn)出明顯非線性，采用描述應(yīng)力-應(yīng)變關(guān)系的雙曲線函數(shù)對推導(dǎo)出的擴展的伯格蠕變模型進行改進，將以上非飽和蠕變模型進一步拓展為非線性彈黏性模型。

(2)根據(jù)模型參數(shù)的物理意義及蠕變曲線特征先求一組接近實際值的初值，然后利用非線性回歸算法進行回歸求得模型參數(shù)，通過參數(shù)的回代，對試驗數(shù)據(jù)進行模擬，結(jié)果表明改進后的模型對蠕變試驗數(shù)據(jù)的模擬效果較好，兩者相關(guān)系數(shù)都在0.95以上，說明改進后的非飽和蠕變模型能較好地反映三峽庫區(qū)某滑坡滑帶土的非飽和蠕變特性。

(3)考慮到模型參數(shù)與凈圍壓之間的線性相關(guān)性，將其線性函數(shù)關(guān)系代入到改進后的非飽和蠕變模型，并采用另外一組試驗數(shù)據(jù)對修正后的模型進行了驗證，結(jié)果表明：修正后的模型合理有效，能較好地描述某滑坡滑帶土的非飽和蠕變趨勢和主要特性。