基于COMSOL Multiphysics的土石壩邊坡滲流穩定數值研究

李 輝，陳大雷，賈 璐

(1.淮安市水利工程建設管理服務中心，江蘇 淮安 223001；2.江蘇淮源工程建設監理有限公司，江蘇 淮安 223009)

我國是水資源大國，怎樣高效利用水資源問題一直是我們研究的熱點。在水利大壩中土石壩是最常見的，大約占比82.7%，因此加強對土石壩的研究意義重大。

土石壩由于其結構適用性強、經濟性好等原因被廣泛應用，同時由于其結構主體因素導致滲流破壞一直是土石壩面臨的問題，引起相關領域專家的關注[1- 3]。數值研究在滲流問題方面一直是研究的方向，津克維茨[4]提出有限單元法應用滲流計算，隨后有限元計算又逐漸將裂隙滲流問題、應力場問題等進行耦合[5]。周慶科[6]、郭仁東[7]、吳良驥[8]等人對有限單元法在滲流計算領域又進一步深入優化。任理[9]、彭華[10]等人對有限元計算的精度和收斂速度進行了改進，使我們的數值研究更加方便準確。李棟[11]、紀偉[12]和馬海兵[13]針對土石壩問題進行了相關數值研究。賀玉珍[14]和占平[15]分別使用ANSYS和AutoBANK軟件對不同情況下土石壩的滲流穩定性進行了相關研究。

本文在前人的研究基礎上，結合實際工程中土石壩的實際特點，對其在不同蓄水情況下滲流穩定特性進行研究。

1 工程概況

某反調節水庫主要由發電樞紐和反調節抽水系統兩部分組成。主體樞紐是Ⅳ等小(Ⅰ)型，由大壩主體、泄洪建筑物等建筑物組成；反調節抽水系統的泵站建筑物是3級建筑物，主要借助抽水泵從水庫匯流池抽取到前池，然后再加壓抽取到對應的輸水隧洞。該反調節水庫的總庫容、調蓄庫容、死水位庫容分別為：987萬、699萬、246.6萬m3；對應的水位分別為：921、920.3、890m。

2 模型建立

COMSOL Multiphysics軟件是瑞典COMSOL公司基于有限單元法開發的一款可以應用于壩體滲流的多物理場耦合分析軟件，該軟件通過求解多物理場對應的偏微分方程來實現多物理場之間的耦合計算，對應計算控制方程見公式(1)，本文涉及的公式較多，在此只列出主要的公式。此軟件可以很好地實現滲流場和邊坡穩定之間的精確計算[16]。

(1)

式中，ρ—流體密度，kg/m3；Cm—單位容水度；Se—有效飽和度；S—儲水系數，Pa-1；Kr—相對滲透系數；μ—流體動力黏度；D—高程，m；P—壓強，Pa；Q—流體源項，m3/s；t—時間，s；Ks—飽和滲透系數，cm/s。

本文根據工程實際情況確定計算需要的滲流參數，見表1。利用COMSOL Multiphysics軟件中的解析函數功能實現非飽和度情況下壓力水頭與滲透曲線和相對滲透曲線的關系曲線，分別如圖1—2所示。壩體計算采用的物理指標參數見表2。

表1 數值計算滲流參數

表2 壩體物理參數

圖1 非飽和土滲透系數曲線

圖2 非飽和土相對滲透系數曲線

使用COMSOL Multiphysics軟件進行建模，為了更準確地進行模擬，采用映射方法劃分網格，在不規則區域局部進行加密，整個模型的網格單元數為2607個，節點數為2686個。不同單元的質量不同，最小單元質量為0.45，平均單元質量為0.92，可以有效的保證計算的精度。建立的有限元模型如圖3所示。

圖3 有限元模型

本文使用COMSOL Multiphysics軟件中的多孔介質滲流模塊進行相關計算。上游邊坡設置為第一類邊界條件，正常蓄水位以上的土石壩壩坡表面和兩側及壩基底部設置為第二類邊界條件。在下游低于正常蓄水位的邊界設置為混合邊界條件。

本文分別對正常蓄水和死水位情況下的滲流穩定特性進行分析。

3 結果分析

3.1 正常蓄水位滲流分析

分別對正常蓄水位作用下的浸潤線、滲透水壓力分布云圖、滲透壓力矢量分布圖、等勢線以及壩體滲透比降分布云圖進行分析，如圖4所示。

圖4 正常蓄水位滲流分析

根據圖4可知，在正常蓄水位時，土石壩的壩體浸潤線從上游接近壩頂的壩坡面開始逐漸降低到下游壩坡面，浸潤線以上為非飽和區域，滲透水壓力為負值；以下為飽和區域，滲透水壓力為正值。對應的滲透水壓力分布云圖如圖4(b)所示，可以知道在最靠近上游的地基滲透水壓力最大，在浸潤線以上滲透水壓力均為負值。根據圖4(c)可知，在滲透水壓力和重力的綜合作用下，水流均向下游流動，上游的水流流速較小，由上游到下游流速先逐漸增大后逐漸減小，最大值出現在壩頂垂線附近。根據圖4(d)可知，等勢線從上游到下游逐漸降低。根據圖4(e)可知，在整個壩體內滲透比降變化不明顯，在結合槽底部出現了明顯的巨變，最大比降就出現在該位置，大小為2.89；在壩體內部的最大滲透比降出現在上游壩坡面浸潤線出現的位置，大小為0.9。

3.2 死水位滲流分析

分別對死水位作用下的浸潤線、滲透水壓力分布云圖、滲透壓力矢量分布圖、等勢線以及壩體滲透比降分布云圖進行分析，如圖5所示。

圖5 死水位滲流分析

根據圖5可知，在從正常蓄水位下降到死水位時，土石壩的壩體浸潤線發生了明顯的變化，在上游壩坡面死水位以上浸潤線出現了明顯的下降，但是在原有正常蓄水位時浸潤線的變化不大，主要體現在上游壩坡面附近。對應的滲透壓力也發生變化，在上游壩坡面死水位以上位置出現了明顯的負的滲透水壓力，其他位置變化不明顯。這主要是因為由于水位下降，在死水位以上的部分上游壩坡面，水庫內沒有水壓力，但是壩體內部依然存在水壓力，導致出流從壩體流向上游，出現了流向上游的流速(如圖5(c)所示)，導致此處的浸潤線下降。

由于上游壩坡面浸潤線的變化，壩體內等勢線也發生變化，最大值出現在上游壩坡面靠近壩頂位置附近。根據圖5(e)可知，最大滲透比降同樣出現在上游壩坡面靠近壩頂位置附近，最大值為1.12。

4 結論

本文使用COMSOL Multiphysics有限元軟件對不同水位情況下的土石壩滲透穩定特性進行分析，得出結論如下。

(1)在正常蓄水位時壩體的浸潤線從上游壩坡面靠近壩頂開始；下降至死水位時上游壩坡面浸潤線出現在死水位附近。

(2)在正常蓄水位時壩體的滲透水壓力負值主要出現在下游壩坡面附近；下降至死水位時在上游壩坡面也出現了滲透水壓力負值。

(3)在正常蓄水位時壩體的最大比降出現在結合槽底部；下降至死水位時最大比降出現上游壩坡面靠近壩頂位置附近。此時容易發生破壞需要特別加強保護。