王 力，王 琳，王世梅，郭 飛，鄒良超

(1.三峽大學(xué) 三峽庫區(qū)地質(zhì)災(zāi)害教育部重點(diǎn)實(shí)驗(yàn)室，湖北 宜昌 443002；2.三峽地區(qū)地質(zhì)災(zāi)害與生態(tài)環(huán)境湖北省協(xié)同創(chuàng)新中心，湖北 宜昌 443002)

1 研究背景

滑坡是一種常見的地質(zhì)災(zāi)害，而降雨和水庫蓄水是滑坡失穩(wěn)的主要誘因。大量實(shí)例表明，水庫蓄水導(dǎo)致滑坡變形演化失穩(wěn)并不是立即發(fā)生的，大都經(jīng)歷了一個(gè)時(shí)間發(fā)展過程[1-2]。如美國Grand Coulee水庫在1941年蓄水后的12 a內(nèi)，先后發(fā)生滑坡500起，其中49%起發(fā)生在蓄水后2 a內(nèi)，51%發(fā)生在2～12 a之間[1]；我國三峽水庫自2003年開始蓄水以來誘發(fā)變形的滑坡有151起，累計(jì)最大位移超過1 000 mm的大型滑坡達(dá)14個(gè)，絕大多數(shù)滑坡變形隨時(shí)間在不斷增加。如樹坪滑坡累計(jì)最大位移量達(dá)到4 300 mm，隨著庫水位周期性漲落，滑坡體地下水呈現(xiàn)周期性升降，地表位移變形速率也呈周期性降低和增大，但滑坡累計(jì)位移總體上呈階梯狀持續(xù)增大趨勢，即使庫水位和地下水位穩(wěn)定期間其變形仍在持續(xù)增大。這說明滑坡變形演化過程不僅與庫水變動(dòng)產(chǎn)生的滲流作用相關(guān)，還與蠕變效應(yīng)相關(guān)。實(shí)際滑坡變形演化預(yù)測常采用現(xiàn)場監(jiān)測和物理模型試驗(yàn)進(jìn)行，但這2種方法主要用于重大滑坡研究，實(shí)際應(yīng)用并不廣泛。具有蠕變特點(diǎn)的滑坡從啟動(dòng)到整體滑動(dòng)這一過程存在階段性的變形，基于非飽和流固耦合理論的數(shù)值方法對于此類滑坡的分析尚有不足，因此需將滲流與蠕變耦合來分析庫水作用下滑坡體內(nèi)實(shí)際力學(xué)過程。 目前國內(nèi)外對土體蠕變研究以試驗(yàn)為主[3]，針對蠕變的數(shù)值計(jì)算研究較少。謝寧和孫鈞[4]提出了一種能方便有效地求解土非線性蠕變問題的二次初應(yīng)變法，該方法在時(shí)步上仍采用與線性蠕變問題相同的黏性初應(yīng)變法，而在每一時(shí)步，采用非線性彈性初應(yīng)變法，通過反復(fù)迭代，求解非線性平衡方程，最終得到每一時(shí)步非線性蠕變的解。主要的商業(yè)數(shù)值計(jì)算軟件也提供了部分蠕變計(jì)算模型，如FLAC提供了經(jīng)典黏彈性模型(馬克斯韋爾體)、伯格蠕變模型、二分量冪定律等蠕變模型[5]；ABAQUS提供了3種可供用戶輸入相關(guān)參數(shù)來定義的蠕變模型：時(shí)間硬化蠕變模型、應(yīng)變硬化蠕變模型和Singh-Mitchell蠕變模型[6]。但是一般商業(yè)數(shù)值軟件中提供的蠕變模型相對較少，參數(shù)難以確定，且極少有考慮非飽和條件的蠕變模型或二次開發(fā)接口，因此對于許多實(shí)際工程問題往往需要自行編程求解。

本文基于前期研究建立的適用于非飽和土的非線性蠕變模型及其三維形式[7]，采用試驗(yàn)數(shù)據(jù)對模型進(jìn)行驗(yàn)證，研究了基于該非飽和土蠕變模型的數(shù)值計(jì)算方法，并編制了滲流及蠕變耦合計(jì)算的有限元計(jì)算程序；利用該程序?qū)θ龒{庫區(qū)樹坪滑坡在庫水位變動(dòng)下的變形過程進(jìn)行了數(shù)值模擬，并與實(shí)際監(jiān)測資料進(jìn)行了對比分析。該項(xiàng)研究可望為分析非飽和土的長期強(qiáng)度和變形、預(yù)測在庫水和降雨作用下滑坡的長期穩(wěn)定性提供一種新理論和新方法。

2 基于伯格模型的非飽和蠕變模型

2.1 蠕變模型介紹

目前，宏觀蠕變模型大致可以分為3類：經(jīng)驗(yàn)-半經(jīng)驗(yàn)?zāi)Ｐ汀⒒谝话懔髯兝碚摰脑Ｐ?、黏彈塑性模型。本文根?jù)團(tuán)隊(duì)成員試驗(yàn)結(jié)果[7-9](不同基質(zhì)吸力的流變曲線)探究含有基質(zhì)吸力的變量——凈圍壓對非飽和土流變的影響規(guī)律;根據(jù)非飽和蠕變曲線特征，類比飽和蠕變模型的建模思想，將凈圍壓作為一個(gè)新的應(yīng)力變量加以考慮，研究凈圍壓與應(yīng)變之間的關(guān)系。從而將基質(zhì)吸力反映在蠕變模型中，得出能夠同時(shí)反映應(yīng)力-應(yīng)變-時(shí)間-基質(zhì)吸力的非飽和土蠕變模型[7]，即

2.2 模型參數(shù)求解及檢驗(yàn)

擴(kuò)展的伯格蠕變模型中含有負(fù)指數(shù)項(xiàng)，其參數(shù)一般需要采用非線性最小二乘法進(jìn)行回歸求得，但是利用非線性最小二乘法進(jìn)行回歸常常會(huì)得到不同的結(jié)果，即參數(shù)存在多組不同的匹配。在本文中，所建模型的參數(shù)都具有較明確的物理含義，因此可根據(jù)其物理意義及蠕變曲線特征先求一組接近實(shí)際值的初值，然后利用非線性回歸算法進(jìn)行回歸求得模型參數(shù)，以避免回歸過程中參數(shù)的跳躍，擴(kuò)展的伯格蠕變模型參數(shù)初值的求解公式為[10]

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

(2)

其中，

對式(2)取對數(shù)，有

ln(εij-A-2aC)=ln(-2aC)-Dt。

(3)

可得

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

(4)

式中A，B，C，D，E，F(xiàn)均為待求系數(shù)。

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

按上述參數(shù)求解方法求解模型的所有參數(shù)。限于文章篇幅,具體選取的數(shù)據(jù)分別為圍壓σ3=250 kPa、ua=200 kPa和圍壓σ3=300 kPa、ua=200 kPa共2組。求得應(yīng)力-應(yīng)變雙曲線函數(shù)的參數(shù)如表1所示。

表1 應(yīng)力-應(yīng)變雙曲線函數(shù)參數(shù)

注：Su為不排水剪切強(qiáng)度;(σ1-σ3)f為實(shí)際破壞剪應(yīng)力。

表2 各組圍壓下模型參數(shù)及相關(guān)系數(shù)

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

從表2可知，上述蠕變模型的模擬效果比較好，與試驗(yàn)數(shù)據(jù)的決定系數(shù)基本都在0.96以上，各級凈圍壓下，模擬蠕變曲線與實(shí)際試驗(yàn)值比較吻合，說明上述蠕變模型可以較好地模擬滑帶土的蠕變特征。

3 非飽和滲流與非飽和土蠕變耦合計(jì)算有限元程序編寫

3.1 求解策略

在非飽和蠕變過程中，隨著土體中含水量的變化，基質(zhì)吸力對蠕變性質(zhì)影響較大，從本文所建立的模型也可以看出基質(zhì)吸力在蠕變模型中為一個(gè)重要參數(shù)。而土體中的基質(zhì)吸力與土體中含水量的變化息息相關(guān)，即與土體中的非飽和滲流過程密切相關(guān)。因此在求解非飽和蠕變模型時(shí)，需要同時(shí)進(jìn)行非飽和滲流計(jì)算。

圖1 滲流與蠕變耦合計(jì)算流程

考慮到蠕變模型與非飽和滲流方程(Richards方程)的相對獨(dú)立性及兩者都與時(shí)間相關(guān)，可以將非飽和蠕變模型和非飽和滲流方程分開求解，即：在同時(shí)間步中先計(jì)算非飽和滲流方程，并將計(jì)算得到的負(fù)孔隙壓力(基質(zhì)吸力)及滲透壓力荷載加入土體的蠕變計(jì)算中，計(jì)算該時(shí)刻的蠕變變形；然后進(jìn)入下一時(shí)步，依次循環(huán)，直到在指定時(shí)間域計(jì)算完畢。滲流與蠕變耦合計(jì)算流程如圖1所示。值得指出的是該求解策略考慮了非飽和滲流對蠕變變形的影響，而忽略了蠕變變形對非飽和滲流的影響。

3.2 滲流蠕變耦合主要子程序

3.2.1 INPUT子程序

INPUT子程序主要讀取單元、節(jié)點(diǎn)、參數(shù)、滲流初始條件、邊界條件及控制信息，為計(jì)算處理好所有的計(jì)算參數(shù)。

3.2.2 UNSATURATED_SEEPAGE子程序

UNSATURATED_SEEPAGE子程序?yàn)橛?jì)算土體的非飽和滲流過程的主程序，計(jì)算出各個(gè)時(shí)間步各節(jié)點(diǎn)的水頭、負(fù)孔隙水壓力及滲透力。

3.2.3 CREEP子程序

CREEP子程序?yàn)槿渥冇?jì)算的主程序，調(diào)用非飽和滲流的計(jì)算結(jié)果，計(jì)算出各個(gè)時(shí)間步的位移場及應(yīng)力應(yīng)變場。

3.2.4 SOLVE子程序

SOLVE子程序主要用于求解組裝后的線性方程組，本程序采用共軛梯度法。由于在本程序中非飽和滲流計(jì)算需要在各個(gè)時(shí)步進(jìn)行多次迭代，求解各時(shí)間步位移場時(shí)也需要依次求解，調(diào)用頻率較高，其計(jì)算效率及正確性直接決定本程序的總體效率和正確性。

3.2.5 OUTPUT子程序

OUTPUT子程序主要輸出程序的計(jì)算結(jié)果，具體為在每時(shí)間步計(jì)算完畢后將計(jì)算結(jié)果如總水頭、孔隙水壓力、位移場、應(yīng)力應(yīng)變場按照Tecplot軟件的文本格式輸出相關(guān)信息，以便用Tecplot軟件打開，進(jìn)行后處理。

為了驗(yàn)證本程序的正確性和可靠性，分別選用具有試驗(yàn)數(shù)據(jù)的非飽和滲流算例及非飽和三軸蠕變試驗(yàn)進(jìn)行計(jì)算，以檢驗(yàn)有限元程序的正確性。

3.3 非飽和滲流有限元計(jì)算程序驗(yàn)證

非飽和滲流的驗(yàn)證算例采用赤井浩一等[12]做的模型試驗(yàn)。試驗(yàn)?zāi)Ｐ蜑?15 cm×23 cm×33 cm(長×寬×高)的均質(zhì)沙槽。沙槽中介質(zhì)的飽和滲透系數(shù)為0.33 cm/s，體積含水量θ與負(fù)孔隙水壓力h及非飽和相對滲透系數(shù)kr的關(guān)系如表3所示。

表3 土體負(fù)孔隙水壓力h、體積含水量θ、非飽和相對滲透系數(shù)kr數(shù)據(jù)

現(xiàn)選取一種工況的試驗(yàn)結(jié)果與計(jì)算結(jié)果相比較。即左右兩端初始水位都在離槽底10 cm處，左端水位瞬時(shí)上升至離槽底30 cm處。由于本程序維度為二維，取模型長-高截面為計(jì)算截面，剖分四邊形網(wǎng)格如圖2所示，共計(jì)節(jié)點(diǎn)768個(gè)，單元693個(gè)。

圖2 沙槽試驗(yàn)計(jì)算模型網(wǎng)格

圖3 不同時(shí)刻水位線分布

圖4 試樣計(jì)算模型

3.4 非飽和蠕變有限元計(jì)算程序驗(yàn)證

非飽和蠕變有限元計(jì)算程序的驗(yàn)證采用一組非飽和三軸蠕變試驗(yàn)數(shù)據(jù)，為了便于比較，仍采用σ′3=200 kPa，σ3=400 kPa，ua=200 kPa，偏應(yīng)力水平為0.50，即偏應(yīng)力188 kPa。試驗(yàn)用土樣尺寸Φ6.18 cm×12 cm，經(jīng)過軸心選取矩形截面為計(jì)算截面，剖分四邊形網(wǎng)格如圖4，共計(jì)節(jié)點(diǎn)325個(gè)，單元288個(gè)。

計(jì)算參數(shù)采用表1求取的模型參數(shù)，計(jì)算時(shí)長為15 000 min，時(shí)間步長參考試驗(yàn)采樣周期，以模型頂部中間節(jié)點(diǎn)(節(jié)點(diǎn)編號為35)作為考察點(diǎn)，繪制該點(diǎn)應(yīng)變曲線并與試驗(yàn)值進(jìn)行比較，結(jié)果如圖5所示。

圖5 應(yīng)變程序計(jì)算值與試驗(yàn)值及模型模擬值比較

從圖5可以看出,程序計(jì)算值與三軸蠕變試驗(yàn)值及模型模擬值基本吻合，其中程序計(jì)算值略偏大，可能是計(jì)算時(shí)考慮了土樣自重應(yīng)力，說明編寫的非飽和蠕變有限元計(jì)算程序正確可靠。

4 算 例

4.1 滑坡簡介

選取三峽庫區(qū)樹坪滑坡為模擬對象，樹坪滑坡為一個(gè)古崩滑堆積體。樹坪滑坡分布高程約為65～500 m,縱向長度約為800 m，橫向?qū)挾燃s為700～900 m，滑坡前緣突入長江，剪出口高程約為65～68 m，滑體厚約40～70 m，體積約為2 600×104m3?；滦螒B(tài)總體呈下陡上緩斜坡，坡度約為22°～35°。

選取滑坡變形較大的1#滑體軸向剖面為計(jì)算截面，計(jì)算模型后邊界取至滑坡后緣460 m高程處，前邊界取至54 m高程處。對計(jì)算截面進(jìn)行四節(jié)點(diǎn)矩形單元網(wǎng)格剖分，共剖分3 855個(gè)單元，4 010個(gè)節(jié)點(diǎn)(如圖6所示)。模型中的左右邊界條件為水平約束，底部約束條件為雙向約束。地層自上而下分別為滑體(上部以粉質(zhì)黏土為主，下部碎石含量較多)、滑帶及基巖。

圖6 樹坪滑坡數(shù)值計(jì)算模型

4.2 計(jì)算參數(shù)

按照本文所建立的非飽和蠕變模型對樹坪滑坡的長期蠕變變形進(jìn)行數(shù)值模擬，需要的非飽和滲流參數(shù)及非飽和蠕變模型參數(shù)數(shù)目較多，若一一通過試驗(yàn)測定或通過參數(shù)反演則工作量巨大。本文主要分析滑坡長期變形趨勢，參數(shù)主要通過以往資料進(jìn)行工程類比以及根據(jù)經(jīng)驗(yàn)值多次試算來確定。樹坪所在的三峽庫區(qū)存在大量涉水型滑坡，根據(jù)附近泄灘滑坡及千將坪滑坡等相關(guān)資料[13-15]以及本團(tuán)隊(duì)相關(guān)非飽和三軸蠕變試驗(yàn)參數(shù)，確定的滑體、滑帶及滑床基巖的非飽和滲流參數(shù)分別為3.0,0.003,0.000 216 m/d。

滑體、滑帶及滑床基巖的土-水特征曲線及滲透性函數(shù)分別如圖7—圖9所示。

圖7 滑體土-水特征曲線及滲透性函數(shù)

圖8 滑帶土-水特征曲線及滲透性函數(shù)

圖9 滑床土-水特征曲線及滲透性函數(shù)

滑體、滑帶及滑床基巖的非飽和蠕變模型參數(shù)如表4 所示。

表4 滑坡非飽和蠕變參數(shù)

4.3 計(jì)算工況

樹坪滑坡外荷載主要有自重荷載和因三峽庫區(qū)庫水位升降及降雨入滲而產(chǎn)生的動(dòng)水壓力。由于每年降雨情況具有較大的隨機(jī)性和離散性，本文暫不考慮降雨工況。因此計(jì)算工況主要考慮庫水位的變動(dòng)情況。由于每年庫水調(diào)度情況大致相近，因此計(jì)算工況選取1 a時(shí)間內(nèi)庫水位實(shí)際調(diào)度周期(2010年6月15日—2011年6月15日)，庫水位如圖10所示。為了計(jì)算方便，忽略庫區(qū)水位小范圍的波動(dòng)，對庫水位的變動(dòng)進(jìn)行線性簡化，即庫水位開始在145 m左右穩(wěn)定一段時(shí)間，然后從145 m上升至 160 m，之后從160 m降至145 m，再從145 m上升至 175 m，然后穩(wěn)定大約75 d，再從175 m降至145 m。

圖10 2010年6月15日—2011年6月15日三峽庫區(qū)水位變動(dòng)

4.4 初始條件

滑坡滲流邊界假設(shè)底部為不透水邊界，滑坡前緣與后緣給定水頭邊界，坡面為透水邊界。計(jì)算應(yīng)力應(yīng)變場時(shí)，假定滑床基巖底部邊界上節(jié)點(diǎn)水平及豎向位移約束為0，前后兩側(cè)節(jié)點(diǎn)水平位移為0?；聺B流場初始值采用庫水位穩(wěn)定在145 m時(shí)的滲流場，初始應(yīng)力應(yīng)變場則采用在自重應(yīng)力及初始滲流場作用時(shí)計(jì)算所得的應(yīng)力應(yīng)變。為了使得非飽和滲流計(jì)算與蠕變計(jì)算同步，在非飽和計(jì)算中，基質(zhì)吸力對孔隙水壓力有顯著的影響，所以設(shè)置不同時(shí)刻的孔隙水壓力，根據(jù)滲流場、位移場在初始時(shí)刻孔隙水壓力分布云圖來分析滲流場的變化。初始狀態(tài)滲流場及位移場孔隙水壓力分布云圖如圖11所示。

圖11 初始時(shí)刻滑坡孔隙水壓力分布云圖

初始應(yīng)力主要有自重應(yīng)力及初始滲透壓力，通過程序計(jì)算可得到x和y方向初始位移分布云圖如圖12所示。

圖12 x和y方向滑坡初始位移

4.5 計(jì)算結(jié)果

計(jì)算時(shí)為了使得非飽和滲流計(jì)算與蠕變計(jì)算同步，時(shí)間步長取為1 d，為了分析不同時(shí)刻滲流場，分別輸出時(shí)間t為1,30，60,120，240，360 d時(shí)的孔隙水壓力,如圖13所示。

圖13 不同時(shí)刻滑坡孔隙水壓力分布云圖

從圖13可以看出，不同時(shí)刻的滲流場變化主要集中于庫水位在145～175 m之間，而其他部位變化不是特別大，主要是因?yàn)閹焖恢饕?45～175 m之間波動(dòng)，對滑體上部影響較為有限。

圖14 不同時(shí)間滑坡位移云圖

圖14給出時(shí)間t分別為30，60，120，240，360 d時(shí)x和y方向滑坡位移分布云圖。

從圖14可以看出，滑坡的變形主要集中在滑體部位，隨著時(shí)間的推移，x和y方向的位移都逐漸增大，x方向最大位移從初始時(shí)刻的0.7 m左右增大至1.2 m左右，y方向位移從初始時(shí)刻的0.15 m左右增大至0.35 m左右。

為了驗(yàn)證計(jì)算結(jié)果，選取主滑截面上具有監(jiān)測資料的點(diǎn)，將變形計(jì)算結(jié)果與實(shí)際監(jiān)測資料進(jìn)行趨勢對比(為便于對比，變形都選取凈變形，即實(shí)際變形與初始時(shí)刻變形的差)，如圖15所示。

圖15 變形計(jì)算結(jié)果與實(shí)際監(jiān)測結(jié)果趨勢對比

從圖15可以看出，監(jiān)測值出現(xiàn)一些波動(dòng)，而計(jì)算值比較平滑，是因?yàn)楸O(jiān)測采用的是多點(diǎn)位移伸縮計(jì)，所測得數(shù)據(jù)在局部存在一定波動(dòng)誤差。計(jì)算結(jié)果與實(shí)際變形的監(jiān)測值盡管在數(shù)值的絕對大小上存在一定差異，但是變形趨勢比較接近，即1 a時(shí)間里，變形逐漸增大，且隨著庫水位的升降劃分為3個(gè)階段：①隨著庫水位的上升，變形逐漸增大；②庫水位保持在175 m左右時(shí)，蠕變變形逐漸進(jìn)入衰減階段；③庫水位下降時(shí)，蠕變變形又逐漸增大，且蠕變增大幅度要比上升時(shí)大。這個(gè)過程相對于庫水位的變動(dòng)存在一定的滯后性，這可能是因?yàn)閹焖坏难杆偕颠^程中，滑坡體內(nèi)滲流場改變存在滯后性。

5 結(jié) 論

(1)本文研究了改進(jìn)后蠕變模型的參數(shù)求取方法，求解出了模型的所有參數(shù)，對試驗(yàn)數(shù)據(jù)進(jìn)行模擬，并采用另一組數(shù)據(jù)對模型進(jìn)行修正，結(jié)果表明：修正后的模型合理有效，能較好地描述該非飽和滑帶土的蠕變趨勢和主要特征。

(2)在對非飽和蠕變模型進(jìn)行有限元數(shù)值分析時(shí)，需要同時(shí)計(jì)算非飽和滲流過程。根據(jù)蠕變模型與非飽和滲流方程的相對獨(dú)立性及兩者都與時(shí)間相關(guān)性，確定了在各個(gè)時(shí)步將非飽和蠕變模型和非飽和滲流方程分開求解的求解策略，并用一個(gè)實(shí)際算例進(jìn)行計(jì)算分析，將計(jì)算結(jié)果與實(shí)際監(jiān)測結(jié)果進(jìn)行對比。結(jié)果表明：計(jì)算結(jié)果與實(shí)際變形的監(jiān)測值趨勢比較接近，本文所建立的蠕變模型合理，編制的程序可靠，可用于實(shí)際滑坡的長期變形計(jì)算。

滑坡土體變形是滲流與蠕變共同作用的結(jié)果，滲流產(chǎn)生的滲透力將會(huì)影響到土體蠕變特性，土體蠕變使得孔隙發(fā)生變化亦將影響到土體滲流，下一步將重點(diǎn)研究蠕變與滲流的耦合作用。