楊 戒，沈朝輝，吳禮舟

(地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護(hù)國家重點(diǎn)實(shí)驗(yàn)室，四川成都610059)

0 引 言

臨水邊坡在坡前水位降落時(shí)，由于水位驟降產(chǎn)生孔隙水滯留，坡內(nèi)水位高于坡外水位，坡內(nèi)外產(chǎn)生水壓差，常發(fā)生失穩(wěn)現(xiàn)象，如1963年的Vajont水庫失事事件[1]。如何真實(shí)反映水位下降過程對(duì)臨水邊坡的影響是一個(gè)值得深入研究的課題。

關(guān)于水位降落對(duì)邊坡穩(wěn)定性影響的研究[2]，目前除試驗(yàn)和單一滲流場(chǎng)理論外，較常用的研究方法是流-固耦合理論，通過考慮滲流場(chǎng)與應(yīng)力場(chǎng)的相互作用來分析坡體穩(wěn)定性，較符合實(shí)際。該理論需建立流-固耦合控制方程組，一般是先將有效應(yīng)力原理引入平衡方程，再將體應(yīng)變引入滲流連續(xù)性方程，最后將兩式聯(lián)立并結(jié)合幾何、本構(gòu)方程即得到流-固耦合控制方程組。但是，目前水位下降邊坡的流-固耦合分析理論大多將土體考慮為彈性材料進(jìn)行應(yīng)力場(chǎng)分析[3-7]，而土體實(shí)際上是一種彈塑性材料，也有不少學(xué)者進(jìn)行了研究[8-11]，但這些研究大都集中在利用強(qiáng)度折減法計(jì)算坡體的穩(wěn)定性系數(shù)上，鮮有對(duì)坡體內(nèi)部滲流、應(yīng)力、位移場(chǎng)耦合變化的研究。

本文基于非飽和流-固耦合理論，對(duì)水位下降邊坡進(jìn)行研究，建立了考慮塑性變形的水位下降邊坡模型，研究土體的彈塑性對(duì)滲流、應(yīng)力和位移多場(chǎng)耦合的影響，并與彈性情況的計(jì)算結(jié)果進(jìn)行比較，可為水位下降邊坡的穩(wěn)定性分析提供參考。

1 理論模型

1.1 基本假設(shè)

介質(zhì)為均質(zhì)各向同性材料；土骨架只發(fā)生微小應(yīng)變；不考慮土顆粒和孔隙水的壓縮，即?n/?t=?εv/?t，?ρf/?t=0。式中，n為孔隙率；t為時(shí)間；ρf為水的密度；εv為土骨架體積應(yīng)變；飽和度Sr與有效飽和度Se相等。忽略孔隙氣壓力的影響，不考慮土體的硬化。

1.2 滲流場(chǎng)控制方程

非飽和土的連續(xù)性方程為

·ρf[-kskr

(1)

式中，Hp為壓力水頭；g為重力加速度；ks為飽和滲透系數(shù)；kr為相對(duì)滲透系數(shù)；y為高程；θ為體積含水率。

結(jié)合假設(shè)和VG模型[12]中的儲(chǔ)水系數(shù)Cm對(duì)式(1)進(jìn)行變形，得

·ρf[-kskr

(2)

常用的VG模型[12]關(guān)系式為

(3)

式中，θr為殘余體積含水率；θs為飽和體積含水率；α、m、N為VG模型參數(shù)，其中m=1-1/N。

1.3 變形場(chǎng)控制方程

應(yīng)用畢肖普有效應(yīng)力原理[13]后的平衡方程為

-·(σ′-SepI)=[ρs(1-n)+ρfn]g

(4)

式中，ρs為土密度；σ′為有效應(yīng)力張量；g為重力加速度向量；I為單位矩陣。幾何方程為

ε=[(U)T+U]/2

(5)

式中，U=(u，v)，u、v分別為水平位移和豎向位移；ε為土骨架應(yīng)變張量。

c和φ表示粘聚力和內(nèi)摩擦角，θR(0≤θR≤π/3)表示羅德角[14]，摩爾-庫倫屈服函數(shù)為F，采用非關(guān)聯(lián)流動(dòng)法則，塑性勢(shì)函數(shù)Q取羅德角θR=π/3時(shí)所對(duì)應(yīng)的屈服函數(shù)，即Q=F(π/3)，這時(shí)的Q即為π平面上不規(guī)則六邊形的外接圓，收斂性良好。其關(guān)系式如下

(6)

綜上，式(2)、(4)為控制方程組，再結(jié)合幾何方程、本構(gòu)方程、初邊值條件就能實(shí)現(xiàn)考慮塑性的非飽和流-固耦合分析。當(dāng)應(yīng)變?cè)隽坎豢紤]塑性部分時(shí)，則土體簡(jiǎn)化為彈性材料。

2 數(shù)值模型及結(jié)果分析

2.1 實(shí)現(xiàn)方法

COMSOL是一款基于偏微分方程組(PDE)建模的多物理場(chǎng)耦合軟件，本文通過修改固體力學(xué)和達(dá)西定律模塊中內(nèi)置PDE方程，實(shí)現(xiàn)用戶自定義控制方程的輸入，并利用如下函數(shù)模擬水位的下降。

坡面變化的壓力水頭邊界：當(dāng)h0-tv0

為驗(yàn)證水位變化模擬的正確性，利用本文彈塑性流-固耦合的方法，對(duì)文獻(xiàn)[15]中的水位下降邊坡試驗(yàn)進(jìn)行了計(jì)算，對(duì)壓力水頭隨水位下降的變化情況進(jìn)行了對(duì)比。邊坡尺寸參考文獻(xiàn)[15]，材料為中砂(滲透系數(shù)ks=2.4×10-4m/s)，試驗(yàn)的初始水位高度h0=0.55 m，水位下降速度v0=0.043 cm/s，坡面水位在第21 min時(shí)下降為0。以邊坡剖面圖的左下角為原點(diǎn)，圖1給出了靠近坡腳和遠(yuǎn)離坡腳的2個(gè)監(jiān)測(cè)點(diǎn)A(0.2，0.2)、B(1，0.2)的壓力水頭隨著坡面水位下降的消散情況。從圖1可以看出，計(jì)算值與試驗(yàn)值較為接近，說明本文的計(jì)算方法是可行的。

圖1 壓力水頭消散對(duì)比

2.2 模型簡(jiǎn)介

本文選取典型的坡度為1∶1的均質(zhì)土坡算例進(jìn)行分析，模型見圖2。其中，AB、BC、AF、EF為不透水邊界(流量q=0)；CD、DE為變化的壓力頭和荷載邊界；AB、EF限制水平位移；AF限制水平和豎向位移；a、b、c為3個(gè)特征點(diǎn)，分別代表坡腳、坡中、坡頂。模型參數(shù)：水密度ρf=1 000 kg/m3、土密度ρs=2 000 kg/m3、粘聚力c=25 kPa、內(nèi)摩擦角φ=25°、孔隙率n=0.4、彈性模量E=10 MPa、泊松比υ=0.3、VG模型參數(shù)α=1.06 m-1、N=1.395、殘余體積含水率θr=0.106、飽和體積含水率θs=0.469。

圖2 邊坡模型(單位：m)

2.3 計(jì)算結(jié)果分析

2.3.1滲流場(chǎng)

圖3給出了當(dāng)ks=10-6m/s、v0=1 m/d時(shí)坡腳處點(diǎn)(28，0)和遠(yuǎn)離坡腳處點(diǎn)(5，0)壓力水頭h隨時(shí)間t的變化情況。從圖3可知，坡面水位在第10 d時(shí)下降為0；在同一高程，離坡表越近，所對(duì)應(yīng)的壓力水頭越接近坡面水位；當(dāng)水位下降至末期水位時(shí)刻，壓力水頭曲線有明顯的轉(zhuǎn)折，且觀測(cè)點(diǎn)位置越接近坡表轉(zhuǎn)折越明顯；考慮彈性與考慮彈塑性得出的壓力水頭曲線幾乎重合，故實(shí)際工程(如土石壩工程)中僅關(guān)注滲流場(chǎng)的變化時(shí)，可將土體簡(jiǎn)化為彈性材料計(jì)算而不會(huì)產(chǎn)生較大誤差。

圖3 壓力水頭消散

2.3.2應(yīng)力場(chǎng)

圖4是正應(yīng)力σx在初始水位和末期水位的分布(ks=10-6m/s，v0=1 m/d)。從圖4可知，初始水位時(shí)，坡腳有應(yīng)力集中現(xiàn)象，σx值較大，坡頂后緣出現(xiàn)受拉；末期水位時(shí)，受拉區(qū)明顯擴(kuò)大，坡腳σx值明顯變大，彈塑性解的坡腳應(yīng)力集中現(xiàn)象還有向深部發(fā)展的趨勢(shì)，而彈性材料無這一現(xiàn)象，但彈性解與彈塑性解的差異在坡腳以外的區(qū)域并不明顯。圖5是剪應(yīng)力τxy在初始水位和末期水位的分布(ks=10-6m/s，v0=1 m/d)。從圖5可知，初始水位時(shí)，坡腳同樣也有應(yīng)力集中現(xiàn)象，τxy值較大；末期水位時(shí)，坡腳τxy明顯變大，彈塑性解坡腳應(yīng)力集中現(xiàn)象也有向深部發(fā)展的趨勢(shì)，若土體強(qiáng)度不足，則可能向深部發(fā)展甚至貫通整個(gè)坡體形成滑面，而彈性解無法反映這種情況。

圖4 正應(yīng)力 σx分布(單位：kPa)

圖5 剪應(yīng)力τxy分布(單位：kPa)

2.3.3位移場(chǎng)

圖6給出了v0=1 m/d時(shí)，3個(gè)特征點(diǎn)a、b、c坡面總位移隨時(shí)間t的變化情況。從圖6可知，當(dāng)ks=10-4m/s時(shí)，彈性解與彈塑性解在點(diǎn)a、b、c都沒有差異；當(dāng)ks=10-5m/s時(shí)，彈性解與彈塑性解在點(diǎn)a、b、c有微小差異；當(dāng)ks=10-6m/s時(shí)，彈性解與彈塑性解在點(diǎn)a、b、c出現(xiàn)差異，且差異從a點(diǎn)至c點(diǎn)逐漸增大。因此，當(dāng)水位下降速度一定時(shí)，總位移的彈性解與彈塑性解的差異與特征點(diǎn)位置和飽和滲透系數(shù)有關(guān)：ks越小，彈性解與彈塑性解的差異越大，尤其在坡腳。

圖6 總位移變化

圖7為水位下降末期的坡面水平位移分布。從圖7可知，ks=10-4m/s時(shí)，彈性解與彈塑性解無差異；ks=10-5m/s時(shí)，彈性解與彈塑性解出現(xiàn)差異，且在坡腳尤為顯著，彈塑性解在坡腳處出現(xiàn)了尖點(diǎn)，而彈性解沒有出現(xiàn)尖點(diǎn)；ks=10-6m/s時(shí)，彈性解與彈塑性解的差異又進(jìn)一步擴(kuò)大。因此，ks越小，彈性解與彈塑性解的差異越大，尤其在坡腳，結(jié)論與圖6相似。

圖7 末期水位坡面水平位移

圖8 點(diǎn)c的總位移

圖9 等效塑性應(yīng)變

由于坡腳應(yīng)力集中較為明顯，且上述分析中考慮彈塑性對(duì)坡腳的影響比其他區(qū)域大，因此繼續(xù)對(duì)坡腳進(jìn)行分析，圖8給出了坡腳c點(diǎn)(ks=10-6m/s)在不同水位下降速度v0的總位移變化。從圖8可知,較大的v0會(huì)產(chǎn)生較大的總位移；隨著水位的逐漸下降，c點(diǎn)彈性解與彈塑性解的差異也逐漸變大(v0越大差異越明顯)；當(dāng)v0=0.1 m/d時(shí)，彈性解與彈塑性解的差異基本消失；v0=0.5 m/d時(shí)，彈性解與彈塑性解在水位下降7 m時(shí)出現(xiàn)差異，且隨著水位的下降，差異逐漸變大；v0=1 m/d時(shí)，彈性解與彈塑性解在水位下降6 m時(shí)出現(xiàn)差異，且在相同水位高度v0=1 m/d時(shí)的差異比v0=0.5 m/d時(shí)大。因此，當(dāng)ks為定值時(shí)，v0越大，彈性解與彈塑性解的差異就越大。

邊坡破壞可看作是塑性區(qū)逐漸發(fā)展、擴(kuò)大直至貫通而進(jìn)入完全塑流狀態(tài)、無法繼續(xù)承受荷載的過程。圖9給出了邊坡塑性區(qū)的動(dòng)態(tài)發(fā)展過程(ks=10-6m/s、v0=1 m/d)。由圖9可知，t=0時(shí)，坡內(nèi)還未產(chǎn)生塑性區(qū)，整個(gè)坡體處于穩(wěn)定狀態(tài)；隨著水位的下降，塑性區(qū)首先在坡腳出現(xiàn)并逐漸向深部發(fā)展，由于巖土體具有一定強(qiáng)度，塑性區(qū)發(fā)展未貫通，邊坡仍處于穩(wěn)定狀態(tài)；若巖土體強(qiáng)度不足，則塑性區(qū)將由坡腳向坡體深部貫通，從而影響邊坡的穩(wěn)定性；塑性區(qū)主要集中在坡腳，這也可解釋為何將土體分別考慮為彈性和彈塑性材料時(shí)，在坡腳體現(xiàn)的差異最大；而隨著水位下降，塑性區(qū)逐漸擴(kuò)大，從而也將導(dǎo)致彈性和彈塑性材料之間的偏差逐漸增大。

3 結(jié) 語

本文基于非飽和土流-固耦合理論，對(duì)水位下降的邊坡進(jìn)行了研究，得出以下結(jié)論：

(1)實(shí)際工程(如土石壩工程)中僅關(guān)注滲流場(chǎng)的變化時(shí)，可將土體簡(jiǎn)化為彈性材料計(jì)算而不產(chǎn)生較大誤差。

(2)彈塑解對(duì)應(yīng)力場(chǎng)的影響主要在坡腳，彈塑性解能反映水位下降時(shí)坡腳應(yīng)力集中向深部發(fā)展的現(xiàn)象，而彈性解無法反映這種情況。

(3)彈塑解對(duì)臨水邊坡位移場(chǎng)的影響與土體飽和滲透系數(shù)和水位下降速度有關(guān)。水位下降速度為定值時(shí)，飽和滲透系數(shù)越小，彈性解與彈塑性解的差異越大，尤其在坡腳；飽和滲透系數(shù)為定值時(shí)，水位下降速度越大，彈性與彈塑性之間的差異越大。