非飽和土滲流特性對庫岸邊坡穩(wěn)定性的影響

劉博，李江騰，王澤偉，廖峻，郝瑞卿

(中南大學 資源與安全工程學院，湖南 長沙，410083)

自從意大利瓦利昂水庫滑坡事件之后，各國學者開始意識到水庫庫岸邊坡的穩(wěn)定性對于工程建設(shè)和正常運營有著重大意義。水庫庫岸滑坡往往發(fā)生在水庫的庫水位發(fā)生急劇變化的階段，在日本約60%的水庫滑坡發(fā)生在水位驟降時，40%發(fā)生在水庫上升時。中村浩之等[1]對日本水庫滑坡進行分析與調(diào)研后認為引起庫岸邊坡滑坡的主要是由邊坡巖土體浸水、庫水位的急劇下降以及降雨等因素引起。鐘立勛[2]調(diào)查了Rasevelt 湖附近地區(qū)1941—1953年發(fā)生的一些滑坡，發(fā)現(xiàn)49%發(fā)生在1941—1942年的蓄水初期，30%發(fā)生在水位驟降10～20 m 的情況下，其余為發(fā)生在其他時間的小型滑坡，由此可見庫水位的漲落對于庫岸邊坡的穩(wěn)定性有著重大的影響，很多學者對此進行了大量研究，如鄭穎人等[3]基于Boussinesq 方程研究了庫水下降過程中地下水浸潤面的位置，并在此基礎(chǔ)上對比分析了根據(jù)經(jīng)驗概化浸潤面位置進行的穩(wěn)定性所得出的誤差及各種因素對庫水邊坡穩(wěn)定性的影響。廖紅建等[4]在基于Geo-slope軟件的基礎(chǔ)上，對不同滲透率和處于不同庫水位下降條件下的邊坡穩(wěn)定性進行了對比分析，并得出了滲透率和庫水下降速度對邊坡穩(wěn)定性的影響規(guī)律。Berilgen[5]在考慮土體模型為 hardening soil模型的基礎(chǔ)上，分析了邊坡穩(wěn)定性和變形性能受土體滲透系數(shù)、水位下降速度、水位下降高度比等因素的影響。黃茂松等[6?7]在考慮土的應(yīng)變軟化效應(yīng)和非飽和土強度準則的基礎(chǔ)上分析了非穩(wěn)定非飽和滲流對土壩穩(wěn)定性的影響。年廷凱等[8]基于abaqus軟件分析了基質(zhì)吸力對于邊坡穩(wěn)定性的影響。影響水庫庫岸邊坡穩(wěn)定性的因素有很多，如庫水位下降速度、下降高度、降雨條件、土體的飽和滲透系數(shù)、浸潤線位置、基質(zhì)吸力等。一些研究者[9]對這些因素對庫岸邊坡穩(wěn)定性的影響進行了大量分析，但是，對于土體的非飽和滲流特性對庫岸邊坡穩(wěn)定性的影響分析很少，為此，本文作者通過采用Fredlund和Xing提出的非飽和土水土特征曲線估算方程[10]對庫岸邊坡穩(wěn)定性受非飽和滲流特性的影響進行分析。

1 計算原理

1.1 土水特性曲線估算公式

土水特性曲線是描述非飽和土體中土體含水量和基質(zhì)吸力關(guān)系的曲線，這一曲線通常由實驗得到的數(shù)據(jù)與某一估算公式進行耦合確定公式中的待定常數(shù)而獲得?，F(xiàn)在常用的估算公式主要有Brooks等提出的公式、Van Genuchten公式和Fredlund和Xing提出的估算公式[10]，其中前兩者缺乏理論基礎(chǔ)，屬于經(jīng)驗公式，也不能準確地模擬整個基質(zhì)吸力范圍內(nèi)的非飽和土含水量；而Fredlund ＆ Xing估算公式則是以毛細定律為理論基礎(chǔ)，認為土體中孔隙的半徑與負孔壓成反比，結(jié)合土體孔隙的概率分布公式推導(dǎo)得來，所以，該估算公式更加科學，并且適用更加廣泛的土體情況和更廣的基質(zhì)吸力范圍，其待定參數(shù)的物理意義相對于前兩者也更加明確。Fredlund ＆ Xing估算公式如下：

其中：wφ為水的體積分數(shù)；CΨ為函數(shù)的修正項，其使得水的體積分數(shù)在基質(zhì)吸力趨向于1 GPa時趨向于0；sφ為飽和水的體積分數(shù)；e為自然數(shù)；Ψ為負孔隙水壓力；iΨ為曲線拐點體積分數(shù)對應(yīng)的基質(zhì)吸力；φi為曲線拐點水的體積分數(shù)；Cr為對應(yīng)于殘余水的體積分數(shù)的基質(zhì)吸力有關(guān)常數(shù)；ks為經(jīng)過拐點的函數(shù)切線的斜率；a，m和n為擬合參數(shù)，通過與實驗數(shù)據(jù)對比確定，但這些參數(shù)具有一定的物理含義。由Fredlund ＆ Xing估算公式推導(dǎo)出的非飽和土滲透系數(shù)函數(shù)為[11]

其中：kw為非飽和土的滲透系數(shù)；ks為飽和土的滲透系數(shù)；kr為相對滲透系數(shù)；y為負孔隙水壓力算法的虛擬變量；j為函數(shù)描述的最小負孔隙水壓力；N為函數(shù)描述的最大負孔隙水壓力；φ′為方程的起始值。

將由式(1)～(7)計算得到非飽和土滲透系數(shù)代入土體的滲流公式中便可以得到非飽和土的滲流場變化。由于非飽和土的滲透系數(shù)與土體水的體積分數(shù)有關(guān)，在計算過程中，非飽和土體的滲透系數(shù)是不斷調(diào)整變化的。

1.2 流固耦合計算原理

在流固耦合計算過程中，認為土坡內(nèi)水的流動服從Darcy定律，其表達式為

式中：vw為流體的速度；s為飽和度；np為孔隙率；ρw為水的密度；uw為孔隙水壓力；k為滲透系數(shù)，該滲透系數(shù)既可以為飽和土的滲透系數(shù)，也可以為非飽和土滲透系數(shù)。

土體內(nèi)流體的連續(xù)性方程如下：

式中：n為邊界外法線方向向量。

滲流場與應(yīng)力場的耦合采用的是基于非飽和土有效應(yīng)力的BIOT固結(jié)理論。

1.3 邊坡穩(wěn)定性計算

通過流固耦合計算過程可以得到邊坡在某一時刻的滲流場與應(yīng)力場，在此基礎(chǔ)上對邊坡穩(wěn)定性的分析計算采用的是由 Zienkiewicz等[12]提出的有限元抗剪強度折減法，其相對于求解安全系數(shù)的極限平衡法具有不需要假定滑動面，適用于復(fù)雜地質(zhì)條件和表面形態(tài)邊坡的優(yōu)勢，現(xiàn)已得到廣泛推廣和應(yīng)用。該方法的強度折減方程如下：

式中：cr′和rφ′為折減后的強度參數(shù)；c′和φ′為實際強度參數(shù)；F為折減系數(shù)。對邊坡穩(wěn)定性的計算應(yīng)該是在單純的靜力條件下進行的，在折減計算過程中，將滲流場設(shè)定為不變值，這樣計算得到的才是滲流發(fā)生至某一時刻的安全系數(shù)。

2 模型計算及分析

2.1 計算模型

為研究非飽和土滲流參數(shù)對庫水下降條件下庫岸邊坡穩(wěn)定性的影響，建立如圖1所示的均質(zhì)土邊坡的平面模型。邊坡的高度為7 m，坡度為1:3，土體的強度參數(shù)參照文獻[13]中的取值：c′=10 kPa，φ′=20°，土體密度ρ=2 t/m3，彈性模量和泊松比分別為E=0.1 GPa，μ=0.3。邊坡土體在飽和狀態(tài)下的滲透率取為ks=1×10?4cm/s。非飽和土的滲流特性主要反映其土水飽和曲線的形式，在Fredlund和Xing提出的估算公式中存在3個待定參數(shù)，并且這3個參數(shù)對曲線的影響個不相同：a可以視為土體的進氣值，m與土體的剩余含水量相關(guān)；n決定了隨著基質(zhì)吸力增大土體含水量變化的速度。為此，參照文獻[10]，選取a為1，10，20，40和70 kPa，m取0.5，1.0和1.5，n取1.2，2.0，4.0和10.0，所得土體的土水特征曲線如圖2所示。其中圖2(a)所示為當m=1和n=4及不同a時的土水特征曲線；圖2(b)所示為當a=5 kPa和n=4及不同m時的土水特征曲線；圖2(c)所示為當a=5 kPa和m=1及不同n時的土水特征曲線。相對滲透率kr隨基質(zhì)吸力的變化曲線見圖3。

圖1 邊坡數(shù)值計算模型Fig. 1 Numerical calculation model of slope

為了實現(xiàn)庫水位下降的條件下庫岸邊坡的模擬，本文通過Fortran語言編寫用戶子程序?qū)baqus軟件進行二次開發(fā)，從而實現(xiàn)了邊坡迎水面的邊界條件隨時間而變化。本文模擬了庫水位以1 m/d的速度勻速下降的情況。

2.2 參數(shù)a對邊坡穩(wěn)定性的影響

Fredlund ＆ Xing推導(dǎo)的估算公式中參數(shù)a能在一定程度上反映非飽和土體的進氣值(得到曲線的實際進氣值略小于a)，這一參數(shù)顯著影響了非飽和土的滲流特性。由圖2(a)可以看出：當a較大時，曲線的起始水平段變長，并且整個曲線趨于平緩，即飽和度隨基質(zhì)吸力增加平緩降低，這一影響趨勢同樣也反映在了相對滲透系數(shù)隨基質(zhì)吸力變化的曲線上。圖 3(a)中的相對滲透系數(shù)也有相同的變化趨勢。參數(shù)a對非飽和土滲流特性的影響將反映在土體的滲流場變化中，在庫岸邊坡穩(wěn)定性分析中，人們最關(guān)心的影響因素之一就是隨著庫水位下降土體中浸潤線的位置。浸潤線反映了坡體中孔隙水壓力的分布與超孔隙水壓力的耗散情況，可以間接地反映邊坡的穩(wěn)定性，對于土體參數(shù)及其他外部邊界條件相同的邊坡，浸潤線的位置越高，邊坡的穩(wěn)定性越差。對于本例中的模型，當庫水位以1 m/d的速度下降7 m后，不同a對應(yīng)的浸潤線位置如圖4所示。從圖4可以看出：隨著a的增大，浸潤線位置會顯著下降，即土體中超孔隙水壓力的耗散程度越大，邊坡越趨于穩(wěn)定。對比文獻[8]可知：a對浸潤線的影響與飽與土滲透系數(shù)對浸潤線的影響相類似。

圖2 a，m和n對土水飽和度的影響Fig. 2 Influence of a, m and n on soil-water saturation

安全系數(shù)F是反映邊坡穩(wěn)定性的直接指標，圖5所示反映了不同a對庫水位下降條件下邊坡安全系數(shù)的影響。由圖5可以看出：當a大于20 kPa時，邊坡的安全系數(shù)隨著水位下降是先減小后增大；當a小于20 kPa時，安全系數(shù)隨著水位下降將一直下降；參數(shù)a對安全系數(shù)的影響隨著水位下降深度的增加而不斷加大，在水位下降7 m后，a=70 kPa時的安全系數(shù)比a=5 kPa時的大23.14%。

圖4 a對土坡內(nèi)浸潤線的影響Fig. 4 Influence of a on phreatic line

圖5 a對邊坡安全系數(shù)F的影響Fig. 5 Influence of a on factor of safety of slope

2.3 參數(shù)m對邊坡穩(wěn)定性的影響

參數(shù)m對非飽和土剩余含水量起控制作用。從圖2(b)可以看出：當m越大時，非飽和土的剩余飽和度越低。而從圖3(b)可見：當m=0.5～1.5時，m對非飽和土的滲透系數(shù)影響很小。m對浸潤面位置的影響如圖6所示。從圖6可以看出：m對于邊坡浸潤面的影響也非常小。圖7所示為m對該庫岸邊坡安全系數(shù)的影響。由圖7可知：當m較小時對應(yīng)的安全系數(shù)略大于m較大時的安全系數(shù)，即非飽和土的剩余滲透系數(shù)對庫岸邊坡的穩(wěn)定性影響不顯著；當庫水位下降7 m后，最大安全系數(shù)僅比最小安全系數(shù)大2.6%。

2.4 參數(shù)n對邊坡穩(wěn)定性的影響

n是控制非飽和土土水特性曲線拐點斜率的參數(shù)。由圖 2(c)可得：n越小，土水特性曲線越平緩；反之，則曲線越陡。同時，n減小也使得進氣值減小。由圖3(c)可得：當n較小時，非飽和土的相對滲透系數(shù)會隨著基質(zhì)吸力更快地減小為一個較小值；n=4和n=10的相對滲透系數(shù)曲線很相近，可以認為當n>4時，其對相對滲透系數(shù)曲線的影響很小。n對浸潤面和安全系數(shù)的影響如圖8和圖9所示。由圖8可得：當n為4和10時，兩者的浸潤面在斜坡段基本重疊。以n=10 時的浸潤線作為參照可以發(fā)現(xiàn)：隨著n減小，其對浸潤面位置的影響越來越大，而不是均勻變化的，這一趨勢與相對滲透系數(shù)曲線相同，也同樣反映在邊坡的安全系數(shù)中，其中最大安全系數(shù)比最小安全系數(shù)大7.5%。

圖6 m對土坡內(nèi)浸潤線的影響Fig. 6 Influence of m on phreatic line

圖7 m對邊坡安全系數(shù)F的影響Fig. 7 Influence of m on factor of safety of slope

圖9 n對邊坡安全系數(shù)F的影響Fig. 9 Influence of n on factor of safety of slope

2.5 結(jié)果分析

由上述計算結(jié)果可看出庫岸邊坡的安全系數(shù)、浸潤面位置與非飽和土的相對滲透系數(shù)與基質(zhì)吸力關(guān)系曲線之間存在相關(guān)性。相對滲透系數(shù)與基質(zhì)吸力關(guān)系曲線可以看作由2個水平段和1個傾斜段組成，起始水平段反映了非飽和土的進氣值，當該段曲線越長時，土體的進氣值越大；而中間的傾斜曲線反映了隨基質(zhì)吸力的增大，土體由飽和滲透率變化至剩余飽和度條件下非飽和滲透率的過程，曲線越斜率越大，則這一變化過程越快。曲線最后的水平段反映了剩余飽和度條件下的非飽和土滲透率。當初始水平段越短，同時傾斜段變化率越大時，因為土體的進氣值較小，隨著庫水位下降，邊坡中的飽和土將很快失水變?yōu)榉秋柡屯?，隨著庫水位進一步下降，土體中的基質(zhì)吸力進一步增大。由于曲線中傾斜段變化率大，非飽和土區(qū)域的滲透系數(shù)將迅速減小，使得水在這一區(qū)域滲流速度大幅度下降，大部分的水必須從飽和區(qū)域出滲流出邊坡，從而影響了邊坡整體的滲流場，浸潤面維持在一個較高的水平；隨著浸潤面下降，難以滲流的非飽和區(qū)域?qū)⒃龃?，從而進一步影響滲流速度，使得浸潤面更加難以下降。而較高的浸潤面位置反映了邊坡中的超孔隙水壓力處于一個較高水平，這將減小邊坡土體中的有效應(yīng)力，從而使得邊坡的穩(wěn)定性下降。反之，當起始水平段越長即進氣值越大時，曲線傾斜段越平緩，非飽和滲透系數(shù)相對基質(zhì)吸力減小較慢，從而對整個邊坡內(nèi)部的滲流速度的影響較小，這有利于超孔隙水壓力的耗散和浸潤面下降，從而也更有利于邊坡的穩(wěn)定。

3 結(jié)論

(1) 非飽和土的相對滲透系數(shù)與基質(zhì)吸力關(guān)系曲線對庫水位下降條件下庫岸邊坡穩(wěn)定性有重要影響。總體來說，曲線的起始水平段越短，傾斜段越陡，土坡內(nèi)浸潤面位置將越高，越不利于邊坡穩(wěn)定性。該曲線受a，m和n這3個參數(shù)不同程度的影響。

(2) 參數(shù)a與非飽和土的進氣值相關(guān)，對非飽和土的相對滲透系數(shù)與基質(zhì)吸力關(guān)系曲線中的起始水平段長度影響顯著；a越小，浸潤面位置越高，邊坡穩(wěn)定性越差。

(3) 參數(shù)m與非飽和土的剩余飽和度有關(guān)，但是，其取值對相對滲透系數(shù)與基質(zhì)吸力關(guān)系曲線影響很小，其對浸潤面和邊坡穩(wěn)定性的影響也很小。

(4) 參數(shù)n和相對滲透系數(shù)與基質(zhì)吸力關(guān)系曲線的初始水平段長度和傾斜段的變化率都相關(guān)，隨著n變化，其對邊坡穩(wěn)定性和浸潤面位置的影響表現(xiàn)出非線性，隨n減小，其影響程度將越來越大。

(5) 對庫岸邊坡穩(wěn)定性產(chǎn)生重要影響的因素不只是飽和土滲流特性及庫水位下降速度等因素，土體的非飽和滲流特性對邊坡的穩(wěn)定性也有重要影響。

[1]中村浩之, 王恭先. 論水庫滑坡[J]. 水土保持通報, 1990,10(1): 53?64.Nakamura K, WANG Gongxian. On reservoir landslide[J].Bulletin of Soil and Water Conservation, 1990, 10(1): 53?64.

[2]鐘立勛. 意大利瓦依昂水庫滑坡事件的啟示[J]. 中國地質(zhì)災(zāi)害與防治學報, 1994, 5(2): 77?84.ZHONG Lixun. Enlightenments from the accident of vaiont landslide in Italy[J]. The Chinese Journal of Geological Hazard and Control, 1994, 5(2): 77?84.

[3]鄭穎人, 唐曉松. 庫水作用下的邊(滑)坡穩(wěn)定性分析[J]. 巖土工程學報, 2007, 29(8): 1115?1121.ZHENG Yingren, TANG Xiaosong. Stability analysis of slopes under drawdown condition of reservoirs[J]. Chinese Journal of Geotechnical Engineering, 2007, 29(8): 1115?1121.

[4]廖紅建, 盛謙, 高石夯, 等. 庫水位下降對滑坡體穩(wěn)定性的影響[J]. 巖石力學與工程學報, 2011, 24(7): 2071?2077.LIAO Hongjian, SHENG Qian, GAO Shihang, et al. Influence of drawdown of reservoir water level on landslide stability[J].Chinese Journal of Rock Mechanics and Engineering, 2011,24(7): 2071?2077.

[5]Berilgen M M. Investigation of stability of slopes under drawdown conditions[J]. Computers and Geotechnics, 2007,34(2): 81?91.

[6]黃茂松, 賈蒼琴. 考慮非飽和非穩(wěn)定滲流的土坡穩(wěn)定分析[J].巖土工程學報, 2006, 28(2): 202?206.HUANG Maosong, JIA Cangqin. Stability analysis of soil slopes subjected to unsaturated transient seepage[J]. Chinese Journal of Rock Mechanics and Engineering, 2006, 28(2): 202?206.

[7]HUANG Maosong, JIA Cangqin. Strength reduction FEM in stability analysis of soil slopes subjected to transient unsaturated seepage[J]. Computers and Geotechnics, 2009, 36(1/2): 93?101.

[8]年廷凱, 萬少石, 蔣景彩, 等. 庫水位下降過程中土坡穩(wěn)定強度折減有限元分析[J]. 巖土力學, 2010, 31(7): 2264?2269.NIAN Tingkai, WAN Shaoshi, JIANG Jincai, et al. Finite element analysis of slope stability under drawdown conditions by strength reduction technique[J]. Rock and Soil Mechanics,2010, 31(7): 2264?2269.

[9]蔣中明, 曾鈴, 付宏淵, 等. 降雨條件下厚覆蓋層邊坡的滲流特性[J]. 中南大學學報(自然科學版), 2012, 43(7): 2782?2788.JIANG Zhongming, ZENG Ling, FU Hongyuan, et al. Seepage characteristics of thick coating slope in raining condition[J].Journal of Central South University (Science and Technology),2012, 43(7): 2782?2788.

[10]Fredlund D G, Xing A. Equations for the soil-water characteristic curve[J]. Canadian Geotechnical Journal, 1994, 31(3): 521?532.

[11]Krahn J. Seepage modeling with SEEP/V: An engineering methodology[R]. Alberta, Canada: GEO-SLOPE International Ltd, 2004: 45, 51?52.

[12]Zienkiewicz O C, Humpheson C, Lewis R W. Associated and non-associated visco-plasticity and plasticity in soil mechanics[J].Geotechnique, 1975, 25(4): 671?689.

[13]Lane P A, Griffiths D V. Assessment of stability of slopes under drawdown conditions[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2000, 126(5): 443?450.