彭汝軒

（廣東省水利水電科學研究院，廣東 廣州 510630）

1 有限元法計算原理

先鋒嶺水庫自建成后從不放空庫容，這使得水庫在蓄水過程中，其庫水位與浸潤線以上滑體一直處于非飽和狀態(tài)、以下滑體處于飽和狀態(tài)。當庫水位下降時，庫水位驟降會使壩體內(nèi)的浸潤線和孔隙水壓在短時間內(nèi)發(fā)生變化，從而對壩坡的穩(wěn)定性產(chǎn)生極其不利的影響。為了保證壩體與壩坡穩(wěn)定安全，根據(jù)基本的滲流理論——達西定律等分析，認為壩體內(nèi)的滲流變化符合達西定律[1]。為此，下文采用有限元方法分析多種邊界條件、多種材料的堤壩或土體的滲流情況。有限元法基本計算公式如下。

（1）滲流有限元分析基本方程為：

式中：K為滲透系數(shù)矩陣（m/d）；H為總水頭向量（cm）；M為單元儲水量矩陣；Q為流量向量（cm3/s）；t為時間（s）。

（2）滲透系數(shù)矩陣[K]為：

式中：K11，K12，K21，K22為滲透系數(shù)張量。

其中，當滲透主軸方向與所選x、y方向一致時，各項異性介質(zhì)表達式為：

式中：[λ]Γ為計算滲流邊界系數(shù)；Kx，Ky為2個坐標軸轉(zhuǎn)向滲透主軸；當設(shè)λ為旋轉(zhuǎn)矩陣時，其表達式為：

式中：a為最大透水系數(shù)的方向角（°），用于各項異性的材料中。

（3）由節(jié)點水頭求單元內(nèi)任意點水頭。其計算公式為：

式中：H為某節(jié)點（x、y）水頭（m）；N為插值函數(shù)的基函數(shù)的數(shù)量，是坐標x、y的函數(shù)，稱為形函數(shù)，其表達式為：

式中：A為單元面積（cm2）；系數(shù)a，b，c和三角形面積A均為節(jié)點坐標x、y的函數(shù)，其表達式為：

（4）由節(jié)點水頭求單元內(nèi)任意點流速。其計算公式為：

式中：V為單元內(nèi)任意點滲透流速（m/d）；當滲透主軸方向與所選的x、y方向一致時，其表達式為：

式中：Vx，Vy分別是x、y滲透主軸方向上的滲流流速（m/d）；Kx，Ky是x、y滲透主軸方向上的滲透系數(shù)。這里，需要注意的是單元內(nèi)任意點的流速均相等。

（5）單元剛度矩陣。根據(jù)有限元法得知，單元剛度矩陣可以表達為：

式中：[K]e為單元剛度矩陣；當單元內(nèi)應(yīng)變元素出現(xiàn)變化，則其表達式為：

2 實例分析

2.1 工程概況

先鋒嶺水庫位于平沙鎮(zhèn)孖髻山西側(cè)，壩址南距平沙鎮(zhèn)2 km，西距崖門水道約6 km，北距斗門區(qū)五山鎮(zhèn)5 km，東距斗門區(qū)五山鎮(zhèn)6 km。水庫三面環(huán)山，屬于低山丘陵地帶。水庫集雨面積3.82 km2，總庫容469.39 萬m3，正常蓄水位17.50 m，相應(yīng)庫容281.09 萬m3；設(shè)計洪水位18.08 m，相應(yīng)庫容386.59萬m3；核核水位18.75 m，相應(yīng)庫容420.93 萬m3。水庫功能以防洪、居民生活供水與農(nóng)田灌溉為主，是西區(qū)蓄淡防咸的樞紐工程，為小（1）型水庫，工程等別為Ⅳ等，主要建筑物為4級。大壩為均質(zhì)土壩，壩頂高程19.55 m、寬5.0 m、長470 m，最大壩高13.88 m，迎水坡比約1∶3.0，背水坡比1∶2.5及1∶3.0；迎水面壩肩設(shè)有高0.9 m、寬0.6 m 的漿砌石防浪墻，防浪墻頂高程20.35 m，壩腳設(shè)排水棱體；壩體填筑土料為黏性土、風化殘積土，密實度不均勻；壩基砂礫含量大，且石層的分布寬度在100.0 m 內(nèi)，厚度在10.0 m 左右，具有較強的透水性；壩體滲透系數(shù)最小值4.30×10-6cm/s，最大值1.18×10-5cm/s。大壩土層物理力學參數(shù)詳見表1，大壩斷面計算模型如圖1 所示。

圖1 大壩計算模型

表1 大壩土層物理力學參數(shù)

2.2 浸潤線情況

為了能夠進一步分析先鋒嶺水庫庫水位下降期上游壩坡的穩(wěn)定性，還需觀察大壩正常蓄水位降落期、校核水位降落期的浸潤線情況。如，正常蓄水位從17.5 m 降落到死水位6.47 m 時，浸潤線在防滲墻的影響下開始出現(xiàn)下降，其下降情況如圖2所示；校核洪水位從18.75 m下降至死水位6.47 m時，浸潤線在防滲墻的影響下出現(xiàn)下降，其下降情況如圖3所示。

圖2 大壩正常蓄水位下降期浸潤線

圖3 大壩校核洪水位下降期的浸潤線

2.3 計算結(jié)果分析

2.3.1 滲流穩(wěn)定計算

本次滲流穩(wěn)定計算采用有限元法，計算軟件采用理正邊坡穩(wěn)定、滲流分析軟件，計算斷面選取大壩典型斷面，計算參數(shù)選取大壩土層物理力學參數(shù)。

（1）計算工況。為了能夠進一步分析水庫庫水位下降期上游壩坡的穩(wěn)定性，還需掌握大壩在各工況下的滲流量、滲流位置、滲流比降等情況。因此，根據(jù)《碾壓式土石壩設(shè)計規(guī)范》（SL274-2020），結(jié)合大壩實際運行中可能出現(xiàn)的不利情況，擬定下列工況進行滲流分析[2]，具體工況如下：①工況1：正常運用條件，正常蓄水位17.50 m 滲流分析，設(shè)計洪水位18.25 m 滲流分析，正常蓄水位17.50 m 降至死水位8.50 m 滲流分析；②工況2：非常運用條件Ⅰ，校核洪水位18.75 m 滲流分析，校核洪水位18.75 m 降至正常蓄水位17.50 m滲流分析，死水位8.50 m滲流分析；③工況3：非常運用條件Ⅱ，正常運用遭遇地震。本區(qū)地震動峰值加速度為0.10 g，地震動反應(yīng)譜特征周期為0.35 s，對應(yīng)的地震基本裂度為Ⅶ度。

（2）滲流量分析。本次大壩滲流量分析主要根據(jù)大壩滲透系數(shù)來復(fù)核斷面在各工況下的滲流量情況。經(jīng)計算，大壩加固前后滲流量詳見表2。

表2 大壩加固前后滲流量計算結(jié)果 m3/d

由表2 可知，大壩加固后在各種工況下滲流量減小，穩(wěn)定性提高。

（3）滲流計算分析。本次大壩滲流計算主要按下式進行：

式中：H1為上游水位（m）；H2為下游水位（m）；L為滲徑長度（m）。

經(jīng)計算，大壩在各種水位組合工況下的滲流計算結(jié)果如圖4—8所示，滲透比降計算成果詳見表3。

圖4 大壩正常蓄水工況滲流計算結(jié)果

圖5 大壩設(shè)計洪水工況滲流計算結(jié)果

圖6 大壩校核洪水工況滲流計算結(jié)果

圖7 大壩校核洪水位驟降至正常蓄水位滲流計算結(jié)果

圖8 大壩正常蓄水位驟降至死水位滲流計算結(jié)果

由表3 可知，在擬定工況下，大壩在正常蓄水位及以上水位運行時，滲流比降最大值出現(xiàn)在壩體與反濾棱體結(jié)合部位，其中最大滲透比降是校核洪水工況下的0.74，浸潤線溢出點均出現(xiàn)在棱體處，且棱體設(shè)有反濾層，可以防止壩體細小顆粒流出，由此認為加固后的大壩滲流安全滿足要求[3]。

2.3.2 壩坡穩(wěn)定分析

（1）計算工況。根據(jù)《碾壓式土石壩設(shè)計規(guī)范》（SL274-2020）有關(guān)規(guī)定，壩坡穩(wěn)定計算斷面選擇最大壩高斷面，計算參數(shù)采用大壩土層物理力學參數(shù)，計算工況擬定下列幾種工況[4]：①工況1：正常運用條件，正常蓄水位17.50 m 迎水坡穩(wěn)定分析，正常蓄水位17.50 m 背水坡穩(wěn)定分析，設(shè)計洪水位18.25 m迎水坡穩(wěn)定分析，設(shè)計洪水位18.25 m背水坡穩(wěn)定分析，正常蓄水位17.50 m降至死水位8.50 m迎水坡穩(wěn)定分析；②工況2：非常運用條件Ⅰ，校核洪水位18.75 m 迎水坡穩(wěn)定分析，校核洪水位18.75 m 背水坡穩(wěn)定分析，校核洪水位18.75 m 降至正常蓄水位17.50 m 迎水坡穩(wěn)定分析，死水位8.50 m 迎水坡穩(wěn)定分析，死水位8.50 m 背水坡穩(wěn)定分析；③工況3：非常運用條件Ⅱ，正常運用遭遇地震。本區(qū)地震動峰值加速度為0.10 g，地震動反應(yīng)譜特征周期為0.35 s，對應(yīng)的地震基本裂度為Ⅶ度。

（2）計算方法。水庫壩坡安全復(fù)核抗滑穩(wěn)定計算采用剛體極限平衡法?？紤]到水庫大壩是均質(zhì)土壩，本工程采用計及條塊間作用力的簡化畢肖法計算圓弧滑動壩坡穩(wěn)定（如圖9 所示），采用有效應(yīng)力法計算壩坡土體抗剪強度，并對庫水位下降期不同工況下壩坡抗滑安全穩(wěn)定系數(shù)進行計算。

圖9 圓弧滑動條分法計算示意

瑞典圓弧法計算公式為[5]：

式中：K為滑動安全系數(shù)；W為土條重力（kN/m3）；b為土條寬度（m）；Q，V分別為地震水平、垂直慣性力（向上為負，向下為正）（kN/m3）；u為作用于土條底面的孔隙壓力（kN）；a為條塊重力線與通過此條塊底面中點的半徑之間的夾角（°）；c'為土條底面的有效應(yīng)力抗剪強度（kPa）；φ'為土的有效內(nèi)摩擦角（°）；Mc為水平地震慣性力對圓心的力矩（kN·m）；R為圓弧半徑（m）。

根據(jù)《碾壓式土石壩設(shè)計規(guī)范》（SL274—2020）相關(guān)規(guī)定，壩坡穩(wěn)定計算采用剛體極限平衡法進行，計算軟件采用GeoStudio2018，并利用滲流計算得到的孔隙水壓力和浸潤線等成果，附加到壩坡穩(wěn)定計算模型中，滑弧計算方法為簡化畢肖普法。根據(jù)不同水位工況，自動搜索最危險滑弧，如圖10—15所示。

圖10 正常蓄水工況迎水坡最危險滑弧

圖11 正常蓄水位驟降至死水位迎水坡最危險滑弧

圖12 正常蓄水位驟降至死水位遭遇地震迎水坡最危險滑弧

圖13 設(shè)計水位工況遭遇地震迎水坡最危險滑弧

圖14 死水位工況迎水坡最危險滑弧

圖15 校核洪水位驟降至正常蓄水位迎水坡最危險滑弧

從圖10—15 可以看出，水庫大壩上（迎水）游壩坡在正常運用和非常運用條件下的抗滑穩(wěn)定安全系數(shù)隨著庫水位的下降而趨于平穩(wěn)，并不存在不利水位。大壩各工況下的壩坡穩(wěn)定計算結(jié)果，詳見表4。

表4 水庫大壩穩(wěn)定計算結(jié)果

由表4 可知，水庫大壩壩坡在各種工況下的抗滑穩(wěn)定最小安全系數(shù)值均大于相關(guān)規(guī)范要求，因此水庫大壩抗滑穩(wěn)定是安全的。

3 結(jié)論

本文基于有限元法基本原理，對水庫大壩在庫水位下降期的滲流變化進行模擬分析，經(jīng)計算發(fā)現(xiàn)，不同工況下壩體滲流比降最大值出現(xiàn)在壩體與反濾棱體結(jié)合部位，其中最大滲透比降是校核洪水工況下的0.74，浸潤線溢出點均出現(xiàn)在棱體處。根據(jù)壩體計算出的滲流結(jié)果，采用極限平衡原理對壩坡不同水位工況下的穩(wěn)定性進行分析，結(jié)果表明，在庫水位下降期，壩體滲透系數(shù)越小浸潤線越低，迎水坡的抗滑穩(wěn)定安全系數(shù)越大，因此水庫大壩抗滑穩(wěn)定安全系數(shù)滿足相關(guān)規(guī)范要求。