陳 曦，張訓(xùn)維，陳佳林，金 鋒，于玉貞

(1．北京交通大學(xué) 土木建筑工程學(xué)院，北京 100044；2. 清華大學(xué) 水沙科學(xué)與水利水電工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100084)

1 引 言

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

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

2 非飽和土滲流Richards方程

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

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

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

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

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

2.2 非飽和滲流有限元求解

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

式中：s為源匯項(xiàng)；η為比存儲(chǔ)量或土-水特征曲線的導(dǎo)數(shù)，即

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖4 不同水位下降速度時(shí)無(wú)心墻壩體安全系數(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位置時(shí)的自由液面Fig.5 Phreatic surface when water level reaches 26 m and 0 m, respectively

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

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

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

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

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

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

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

5 結(jié) 論

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

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

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

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

[1] 劉新喜, 夏元友, 練操, 等. 庫(kù)水位驟降時(shí)的滑坡穩(wěn)定性評(píng)價(jià)方法研究[J]. 巖土力學(xué), 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] 張文杰, 陳云敏, 凌道盛. 庫(kù)岸邊坡滲流及穩(wěn)定性分析.水利學(xué)報(bào), 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]. 巖土工程學(xué)報(bào), 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]. 巖土力學(xué), 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]. 巖土力學(xué), 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]. 巖土力學(xué), 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.