水位下降邊坡的流-固耦合彈塑性分析

楊 戒，沈朝輝，吳禮舟

(地質災害防治與地質環境保護國家重點實驗室，四川成都610059)

0 引 言

臨水邊坡在坡前水位降落時，由于水位驟降產生孔隙水滯留，坡內水位高于坡外水位，坡內外產生水壓差，常發生失穩現象，如1963年的Vajont水庫失事事件[1]。如何真實反映水位下降過程對臨水邊坡的影響是一個值得深入研究的課題。

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

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

1 理論模型

1.1 基本假設

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

1.2 滲流場控制方程

非飽和土的連續性方程為

·ρf[-kskr

(1)

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

結合假設和VG模型[12]中的儲水系數Cm對式(1)進行變形，得

·ρf[-kskr

(2)

常用的VG模型[12]關系式為

(3)

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

1.3 變形場控制方程

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

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

(4)

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

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

(5)

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

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

(6)

綜上，式(2)、(4)為控制方程組，再結合幾何方程、本構方程、初邊值條件就能實現考慮塑性的非飽和流-固耦合分析。當應變增量不考慮塑性部分時，則土體簡化為彈性材料。

2 數值模型及結果分析

2.1 實現方法

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

坡面變化的壓力水頭邊界：當h0-tv0

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

圖1 壓力水頭消散對比

2.2 模型簡介

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

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

2.3 計算結果分析

2.3.1滲流場

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

圖3 壓力水頭消散

2.3.2應力場

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

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

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

2.3.3位移場

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

圖6 總位移變化

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

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

圖8 點c的總位移

圖9 等效塑性應變

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

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

3 結 語

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

(1)實際工程(如土石壩工程)中僅關注滲流場的變化時，可將土體簡化為彈性材料計算而不產生較大誤差。

(2)彈塑解對應力場的影響主要在坡腳，彈塑性解能反映水位下降時坡腳應力集中向深部發展的現象，而彈性解無法反映這種情況。

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