南溝門水庫滲流計算分析

龐耀立，董政良，張曉飛

（1.陜西省水務清潔能源集團有限公司，陜西 西安 710082；2.西安理工大學水利水電學院，陜西 西安 710048）

1 工程概況

南溝門水庫樞紐工程位于陜西省延安市黃陵縣境內，水庫壩址位于葫蘆河河口上游約3 km 處的寨頭河村南溝門附近。南溝門水利樞紐工程建設任務為工業和城鄉生活供水，兼顧灌溉和發電等綜合利用。南溝門水庫壩址斷面多年天然徑流量為1.31 億m3，100 年一遇洪峰流量為1130 m3／s，三日洪量為3410 萬m3；5000 年一遇洪峰流量為2910 m3／s，三日洪量為6430 萬m3。南溝門水庫屬Ⅱ等大（二）型工程，壩設計防洪標準為100 年一遇，校核洪水標準為5000 年一遇。工程區地震基本烈度為Ⅵ度。大壩是均質土壩，壩頂高程為852.0 m，大壩壩頂寬度為10 m，在其壩頂上游側設置了1.2 m 高的防浪墻，墻頂高程為853.2 m，大壩最底點開挖高程789.0 m，最大壩高為63.0 m，壩頂總長為504.4 m；壩體上游邊坡1∶2.75 和1∶3，在高程830 m 處設寬3.0 m 馬道；下游壩坡比均為1∶2.5，在高程835 m 和818 m 處設寬2.0 m 的馬道，在800.5 m 處設棱體排水平臺，最大壩底寬為351 m。壩體上游壩坡在高程812.5 m 以上采用厚0.4 m 的干砌石護坡，下游壩坡采用漿砌石網格內植草皮護坡。庫區地處陜北黃土高原強烈侵蝕區，河流切割基巖深達50 余米，兩岸溝谷縱橫，塬、梁、峁黃土地貌景觀明顯，相對高差200 m～300 m。河流階地發育不對稱，凹岸巖石裸露，地形陡竣；凸岸階地發育，第四系松散堆積覆蓋廣泛。庫盆周邊巖體完整，無大的構造現象，庫岸基巖頂板及地下水位高于水庫正常水位，無滲漏條件，不存在永久性滲漏問題[1］。

2 計算模型及參數

2.1 模型基本假定

土石壩滲流場與溫度場之間是相互影響、相互作用的，一方面，溫度影響水的密度和動力粘度，進而影響土的滲透系數，熱量的運移又會導致土中水分的遷移；另一方面，滲流作用引起土體中的熱量發生轉移，進而影響溫度場的空間分布。兩場耦合作用非常復雜，為了簡化求解同時又能體現關鍵要素，且能保證計算結果可靠的情況下，做以下基本假設：

1）土石壩屬于多孔連續介質，土石壩內的滲流是連續介質滲流；

2）水流運動過程中只受重力作用；

3）溫度不會導致水體發生相變；

4）不考慮因介質摩擦產生的熱量；

5）水流符合達西定律。

2.2 滲流場與溫度場耦合數學模型

不考慮流體與多孔介質之間熱平衡所需時間，即認為兩者之間的熱傳遞是瞬時完成的，忽略熱輻射、化學反應等引起的熱消耗，且流體密度僅和溫度及壓強有關，得到三維多孔介質地下水的水-熱運移耦合方程[2］：

溫度場：

式中：QG為溫度場的源匯項；Γ1和Γ2分別表示第一類邊界和第二類邊界；Φ0為第一類邊界上定溫度值；n 為第二類邊界上的外法線方向；Q0為第二類邊界上的單寬熱流量；T0為初始溫度分布；t 為時間；Ω 表示研究區域。

滲流場：

式中：Qs為滲流場的源匯項；P1為第一類邊界上水頭高度；P0為初始水頭分布；q1為第二類邊界上的通量；其他變量如前所述。

2.3 邊界條件及參數

2.3.1 邊界條件

滲流場方面，壩體上游定水頭為 852 m，邊界為上游坡面，壩體下游定水頭789 m，邊界為下游坡面，其余邊界均為零通量邊界。

溫度場方面，當地溫度隨時間周期性變化的傅里葉級數函數如下：

式中：T0為平均溫度；k 代表傅里葉級數的階數；Ak和Bk代表傅里葉級數的系數；L 為周期[3］。

這里采用陜西省延安市的氣溫數據（數據來源于中國氣象數據網站）并根據傅里葉函數來取5 階傅里葉級數就可以滿足精度要求。擬合得氣溫回歸模型系數為：T0=16.37、A1=-10.28、A2=-0.5417、A3=0.2667、A4=-0.2917、A5=-0.08822、B1=-7.133、B2=-0.2165、B3=0.3、B4=-0.0433、B5=0.08344。將系數帶入上式即可求得溫度隨時間周期性規律，溫度為16.37 ℃。

2.3.2 模型參數

由于本文采用COMSOL 軟件進行數值模擬，而COMSOL中自帶計算傳熱模塊是用熱量運移方程進行計算的，默認為導熱系數模型為體積平均模型，因此用PDE 模塊來代替COMSOL 自帶的傳熱模型進行計算，可以靈活修改導熱系數經驗模型，并通過修改相關的傳熱模塊來實現不同導熱系數條件下的土石壩流熱耦合模型求解。模型滲流場及溫度場的計算參數[4］，見表1，模型構圖見圖1[5］。

表1 土石壩流-熱耦合模型計算參數

續表2

圖1 模型結構示意圖

3 計算成果分析

生成的結果見圖2，可以發現上游壩體基本處于飽和狀態，帷幕灌漿區域軸線附近地帶的壩體中間部位處于飽和與非飽和的交界狀態，下游壩體則是大部分處于非飽和的狀態。

由圖2 可知由于帷幕灌漿區域滲透性比較低，等壓線相對密集，因此地下水滲流的流速極其緩慢，壩體浸潤線（圖中零壓力線）穿過帷幕灌漿區時急劇降低。該區域前壩體在庫水的作用下，其壓力水頭分布較為平緩。

圖2 2019 年6 月典型時刻孔隙水壓力水頭分布圖

在下游區域，壩體浸潤線逐漸降低變緩，與排水褥墊層的上表面高程基本一致，說明穿過帷幕灌漿區域的庫水主要通過排水褥墊層排出壩體，結合下游坡腳處的排水棱體起到了很好地防滲作用。因此，帷幕灌漿能有效降低壩體浸潤線，穿越帷幕灌漿這部分滲漏庫水能通過排水褥墊層和排水棱體排出，帷幕灌漿和排水棱體能有效保證土石壩的運行安全，計算結果合理。

4 結論

堤壩流熱耦合理論仍然是一個年輕的研究領域，模型和方法的發展至關重要。本文通過理論分析堤壩溫度場與滲流場相互影響，建立了土石壩飽和-非飽和流熱耦合的數學模型，并選取南溝門水庫攔河壩進行土石堤壩的滲漏研究，采用COMSOL Multiphysics 建立有限元模型來分析堤壩溫度場和滲流場的變化規律，并通過修改相關的傳熱模塊來實現不同導熱系數條件下的土石壩流熱耦合模型求解。發現滲流未達到穩定時，壩體浸潤線在經過帷幕灌漿區后依然急劇降低，下游浸潤線的高程與排水褥墊層基本相同，非飽和帶飽和度由壩體浸潤線至壩體外表面逐漸減小。結果可為今后陜西省土石壩工程的防滲提供參考依據。