陳 曦，張訓維，陳佳林，金 鋒，于玉貞

(1．北京交通大學 土木建筑工程學院，北京 100044；2. 清華大學 水沙科學與水利水電工程國家重點實驗室，北京 100084)

1 引 言

我國水庫大壩建設歷史悠久，全國已建成大小水庫總數(shù)量接近9萬座。 部分水庫的壩體由于長期服役，存在不同程度的隱患，環(huán)境的變化是誘發(fā)水庫風險的主要因素，例如，降雨入滲、水庫蓄水后庫水位驟降都可能引起水庫壩體安全系數(shù)下降，導致壩體危險產(chǎn)生。 許多學者對庫水位波動所引起的壩體滲流和壩體穩(wěn)定性變化進行了研究。 劉新喜等[1]基于 Richards飽和-非飽和滲流控制方程建立了有限元格式，分析了庫水位驟降時的瞬態(tài)滲流場，評價了滲流力作用下滑坡穩(wěn)定性分析的不平衡推力法，并認為滑坡的滲透系數(shù)和庫水位下降速度是影響壩體穩(wěn)定性的主要因素。張文杰等[2]基于GeoStudio中的Seep/W和Slope/W模塊對庫岸邊坡的滲流場和穩(wěn)定性進行了分析，其非飽和土的抗剪強度基于 Fredlund公式，穩(wěn)定性分析基于 Bishop極限平衡法，最后給出了坡體安全系數(shù)隨水位波動的變化曲線。陳曦等[3]認為非飽和土坡穩(wěn)定性分析方法的發(fā)展可分為3個階段，即時間無關的定性分析階段、時間相關的滲流和變形不完全耦合分析階段以及時間相關的滲流和變形完全耦合的分析階段，目前滲流穩(wěn)定性分析方法以第二或第三階段的分析方法為主。徐炎兵等[4]則基于多孔介質(zhì)理論，建立了非飽和土兩相流動與變形的耦合模型，開發(fā)了有限元分析程序，并采用經(jīng)典算例和離心機試驗對數(shù)值結果進行了驗證。

本文采用滲流與變形的不完全耦合分析方法，即先采用Richards方程進行非飽和滲流場的有限元求解，再基于非飽和滲流場和非飽和土的抗剪強度理論對壩體的穩(wěn)定性進行評價，通過數(shù)值算例分析了水庫庫水位波動對壩體穩(wěn)定性的影響。

2 非飽和土滲流Richards方程

2.1 非飽和土土-水特征曲線模型

常用的土-水特征曲線模型有Brooks-Corey 模型，Van Genuchten模型（簡稱VG模型）和Fredlund-Xing 模型。由于VG模型具有連續(xù)平滑的曲線，在進氣壓力值和接近殘余含水率狀態(tài)時具有較好的平滑過渡，本文采用VG土-水特征曲線模型，表達式為

式中：h為壓力水頭；θ、θr、θs分別為當前體積含水率、殘余體積含水率和飽和體積含水率；av、nv、mv為經(jīng)驗系數(shù)，mv通常簡化為mv=1- 1 /nv。VG模型通常是av、nv、mv的3參數(shù)的模型，當θs、θr不能直接給出時，VG模型可作為 5參數(shù)模型。式（1）也可表示為標準化的體積含水率Θ，即

非飽和滲流的水力傳導系數(shù)可采用Mualem公式進行描述：

式中：經(jīng)驗參數(shù)l一般取為1.0；ks為飽和水力傳導系數(shù)，為常數(shù)；kr(Θ)為相對滲透張量，是有效飽和度Θ的函數(shù)。

2.2 非飽和滲流有限元求解

非飽和土滲流建模的Richards方程可采用3種基本格式，即壓力水頭格式（h-form）、體積含水率格式（θ-form）和混合格式（mixed form），下面的壓力水頭格式的Richards方程應用更為廣泛，

式中：s為源匯項；η為比存儲量或土-水特征曲線的導數(shù)，即

式中：γw為水的重度；mw為土-水特征曲線的斜率。

對式（4）進行空間離散和時間差分，可得Richards有限元方程的修正Picard迭代公式：

式中：下標n、m分別為時間步指標和修正Picard方法迭代步指標；M、C分別為質(zhì)量矩陣和水力傳導特征矩陣。

3 非飽和土抗剪強度與穩(wěn)定性分析

進行非飽和土體的穩(wěn)定性評價時需要采用合理的非飽和土抗剪強度理論，Sheng等[5]認為，盡管目前所提出的非飽和土抗剪強度理論和公式較多，但大多可以看作是 Bishop和 Blight強度理論和Fredlund強度理論的變體，各種非飽和土強度理論的主要差別在于參數(shù)或變量的選取不同。Bishop和Blight強度理論一般可表示為

式中：c′、φ′分別為土的有效黏聚力和有效內(nèi)摩擦角；χ為與飽和度(Sr)相關的參數(shù)；Fredlund強度理論考慮了一個描述基質(zhì)吸力對剪切強度貢獻的角度φb，則式（8）可表示為

式中：抗剪強度公式中第三項是基質(zhì)吸力的貢獻，可將其合并到真黏聚力，形成“假黏聚力”c′，即

對于 ta nφb，本文采用 Fredlund-Oberg公式，即

土體穩(wěn)定性分析方法較多，其中有限元強度折減法實際應用效果良好[6]。為了獲得應力和應變分布、位移場，塑性分布等更多信息，本文采用有限元強度折減法進行土體穩(wěn)定性分析，仍為基于單網(wǎng)格強度折減技術，為獲得更高的求解效率，可以采用雙網(wǎng)格強度折減技術[7]。

4 水位波動條件下壩體滲流和穩(wěn)定性分析

4.1 壩體臨水坡坡度大于背水坡坡度的算例

USSA為開發(fā)的非飽和土滲流與穩(wěn)定性分析軟件，用戶界面采用Python腳本語言，主計算程序采用FORTRAN語言，并通過算例驗證了計算程序的可靠性[8]。如圖1所示，算例1為均質(zhì)土壩，其左右坡坡度分別為1︰1.8和1︰1.4；土壩左側初始水位為0 m，右側初始水位為7.05 m。在此初始條件下達穩(wěn)態(tài)滲流后，右側水位從初始水位以1 m/d的較快速度下降，最后水位降為0 m。圖2為水位下降過程中壩體邊坡安全系數(shù)的變化過程，實線和實心標識為真正的邊坡安全系數(shù)。在水位下降不超過4.11 m時，左側邊坡為危險邊坡，并先于右側邊坡發(fā)生破壞；當水位下降超過4.11 m時，右側邊坡先于左側邊坡發(fā)生破壞。虛線和空心標識為猜測的左右側邊坡安全系數(shù)演化過程，即隨著水位下降，左側邊坡的安全系數(shù)緩慢增加，右側邊坡的安全系數(shù)由于水壓的降落而迅速減少，并在水位下降超過4.11 m時，開始小于左側邊坡的安全系數(shù)。水位下降過程中，左側邊坡安全系數(shù)緩慢增加歸因于非飽和區(qū)比例的逐漸增加；右側邊坡則相對復雜，由于水壓降落的影響要高于非飽和區(qū)比例逐漸增加的影響，導致右側邊坡安全系數(shù)的總體變化趨勢是減小的。左右側坡體安全系數(shù)變化曲線相交點約為FOS≈1.56。

圖1 均質(zhì)土壩有限元網(wǎng)格劃分Fig.1 Finite element mesh of homogeneous soil dam

圖2 水位降低過程中均質(zhì)土壩安全系數(shù)的變化曲線Fig.2 Variation curves of safely factor with the drawdown of water level for homogeneous dam

4.2 壩體臨水坡坡度等于背水坡坡度的算例

算例2的整個區(qū)域由3種土組成，即壩體土體、壩體心墻材料和壩基土體，如圖3所示為整體分析區(qū)域的有限元網(wǎng)格劃分。3種土體的彈塑性力學參數(shù)見表 1，其中φ、φ分別為土的內(nèi)摩擦角和剪脹角。非飽和土土性參數(shù)見表2。3種土的不飽和自然重度和飽和重度分別為17 kN/m3和20 kN/m3，其中壩基和心墻材料采用非飽和土數(shù)據(jù)庫 UNSODA中編號為3 091的砂質(zhì)壤土和編號為1 183的黏土[9]。

圖3 土壩有限元網(wǎng)格劃分Fig.3 Finite element mesh of earth dam with core wall

表1 壩體土體彈性和強度力學參數(shù)Table 1 Elastic and strength parameters for dam soils

表2 壩體土體非飽和土性參數(shù)Table 2 Unsaturated hydraulic parameters for soils of dam

首先模擬沒有心墻的情況，壩體在右側水位為38 m，左側水位為0 m條件下達穩(wěn)定滲流。分別采用0.01 m/d和 1 m/d兩種不同的水位下降速度，分析了水位下降過程中壩體安全系數(shù)的變化，如圖 4所示。對于無心墻壩體，當水位下降速度較低時（如0.01 m/d），壩體的安全系數(shù)呈緩慢穩(wěn)定增長趨勢；但是當水位下降速度較快（1 m/d時），安全系數(shù)基本維持在某一值（2.27～2.28）附近，直到水位下降了34 m或下降到4 m左右，安全系數(shù)開始明顯減小?？捎蓤D5來理解上述變化，對于左側坡體，水位下降到26 m時和水位下降到0 m時的自由液面相差較小，而在右側坡體內(nèi)則相差較大。當右側坡體缺少足夠的水壓力保護時，右側坡體的安全系數(shù)開始小于左側坡體，這與算例1的情況類似，只不過由于左右側坡體坡度不同，左右側坡體安全系數(shù)變化曲線的相交點的位置也有所不同。

圖4 不同水位下降速度時無心墻壩體安全系數(shù)變化曲線Fig.4 Variation curves of safely factor for different drawdown speed of water level for dam without core wall

圖5 水位以1 m/d速度下降到26 m和0 m位置時的自由液面Fig.5 Phreatic surface when water level reaches 26 m and 0 m, respectively

考慮心墻時，心墻與壩體材料見表1、2。壩體在右側水位為38 m，左側水位為0 m條件下達穩(wěn)定滲流，此時壩體、心墻和壩基的體積含水率的分布陰影圖和壓力水頭分布如圖6所示。從圖中可以看出，自由液面以下土體都達到各自的飽和體積含水率，心墻的飽和體積含水率（0.521）明顯高于壩體和壩基土體的飽和體積含水率。

同樣，壩體右側水位從38 m分別以0.01 m/d和1 m/d的速度下降。圖7為水位下降過程中安全系數(shù)的變化曲線。從圖中可以看出，由于心墻的存在，壩體安全系數(shù)略高于無心墻的情況。對于不同的水位下降速度，在水位下降了34 m或下降到4 m左右，安全系數(shù)都有明顯的減小，且在水位下降為0 m時都有回升。圖8為以0.01 m/d的水位下降速度，水位分別下降至26 m和0 m時壩體內(nèi)體積含水率陰影圖。由于滲透性較低，心墻會對左側坡體的自由液面下降起到一定的阻尼作用。

圖6 水位38 m時的穩(wěn)態(tài)滲流場Fig.6 Finite element mesh of earth dam with core wall

圖7 不同水位下降速度時心墻壩體安全系數(shù)變化曲線Fig.7 Safely factor curves of dam with core wall for different drawdown speed of water level

圖8 水位下降速度為0.01 m/d時壩體內(nèi)體積含水率陰影Fig.8 Contour of volumetric water content in dam for water level drawdown speed of 0.01 m/d

以左右兩側為0 m水位條件下達穩(wěn)定滲流，在此條件下分別以 0.01 m/d和1 m/d的速度蓄水。圖 9為水位上升過程中安全系數(shù)的變化曲線。可見，在水位上升初期，即約為0～8 m時，安全系數(shù)會有所提高，表明這個期間由于水壓的施加，右側邊坡的安全度一直增加，而左側邊坡的安全度隨著滲流的發(fā)生，非飽和區(qū)逐漸減小，自重增加，導致其安全系數(shù)逐漸減小。在達到約8 m水位后，安全系數(shù)開始小于右側邊坡。因此，水位上升過程中安全系數(shù)變化曲線亦可看成左右側坡體安全系數(shù)交叉獲得的曲線。需要說明的是，水位上升速度較快（即1 m/d）時，由于心強的隔水和防滲作用，左側坡體滲流自由液面較低，安全系數(shù)下降較為平緩，而水位上升速度較慢（即0.01 m/d）時，安全系數(shù)下降更為明顯。

圖9 不同水位上升速度時心墻壩體安全系數(shù)演化曲線Fig.9 FOS curves for different rise speed of water level

5 結 論

（1）左右壩體邊坡各存在一條安全系數(shù)變化曲線，真實的安全系數(shù)變化曲線是這兩條安全系數(shù)變化曲線交叉獲得的較小的部分。

（2）水壓對壩體邊坡的穩(wěn)定性具有重要的影響，通常高水位有益于臨水邊坡的穩(wěn)定性，而臨水邊坡的穩(wěn)定性通常由水壓和非飽和區(qū)比例變化共同支配。

（3）背水面邊坡的分析相對簡單，在水位變化過程中背水面邊坡的非飽和區(qū)的比例變化所支配。

（4）心墻具有對水位變化所產(chǎn)生的背水面坡體滲流場變化的阻尼作用。

[1] 劉新喜, 夏元友, 練操, 等. 庫水位驟降時的滑坡穩(wěn)定性評價方法研究[J]. 巖土力學, 2005, 26(9): 1427－1431.LIU Xin-xi, XIA Yuan-you, LIAN Cao, et al. Research on method of landslide stability valuation during sudden drawdown of reservoir level[J]. Rock and Soil Mechanics, 2005, 26(9): 1427－1431.

[2] 張文杰, 陳云敏, 凌道盛. 庫岸邊坡滲流及穩(wěn)定性分析.水利學報, 2005, 36(12): 1510－1516.ZHANG Wen-jie, CHEN Yun-min, LING Dao-sheng.Seepage and stability analysis of bank slopes[J]. Journal of Hydraulic Engineering, 2005, 36(12): 1510－1516.

[3] 陳曦, 劉春杰. 逐步完善的非飽和土邊坡穩(wěn)定性有限元分析方法[J]. 巖土工程學報, 2011, 33(增刊1): 380－384.CHEN Xi, LIU Chun-jie. Staged development of finite element methods for stability of unsaturated soil slopes[J].Chinese Journal of Geotechnical Engineering, 2011,33(Supp.1): 380－384.

[4] 徐炎兵, 韋昌富, 李幻, 等. 非飽和土滲流與變形耦合問題的有限元分析[J]. 巖土力學, 2009, 30(5): 1490－1496.XU Yan-bing, WEI Chang-fu, LI Huan, et al. Finite element analysis of coupling seepage and deformation in unsaturated soils[J]. Rock and Soil Mechanics, 2009,30(5): 1490－1496.

[5] SHENG D, ZHOU A N, FREDLUND D G. Shear strength criteria for unsaturated soils[J]. Geotechnical and Geological Engineering, 2011, 29(2): 145－159.

[6] 于玉貞, 林鴻州, 李榮建, 等. 非穩(wěn)定滲流條件下非飽和土邊坡穩(wěn)定分析[J]. 巖土力學, 2008, 29(11): 2892－2898.YU Yu-zhen, LIN Hung-chou, LI Rong-jian, et al.Stability analysis of unsaturated soil slope under transient seepage flow state[J]. Rock and Soil Mechanics, 2008,29(11): 2892－2898.

[7] CHEN X, WU Y, YU Y, et al. A two-grid search scheme for large-scale 3-D finite element analyses of slope stability[J]. Computers and Geotechnics, 2014, 62: 203－215.

[8] 陳曦, 于玉貞, 程勇剛. 非飽和滲流 Richards 方程數(shù)值求解的欠松弛方法[J]. 巖土力學, 2012, 33(增刊 1):237－243.CHEN Xi, YU Yu-zhen, CHENG Yong-gang. Underrelaxation methods for numerical solution of Richards equation of variably saturated flow[J]. Rock and Soil Mechanics, 2012, 33(Supp.1): 237－243.

[9] NEMES A, SCHAAP M G, LEIJ F J. The UNSODA unsaturated soil hydraulic database(version 2.0)[M].Riverside, CA: US Salinity Laboratory, 1999.