龐耀立，董政良，張曉飛

（1.陜西省水務(wù)清潔能源集團(tuán)有限公司，陜西 西安 710082；2.西安理工大學(xué)水利水電學(xué)院，陜西 西安 710048）

1 工程概況

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

2 計(jì)算模型及參數(shù)

2.1 模型基本假定

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

1）土石壩屬于多孔連續(xù)介質(zhì)，土石壩內(nèi)的滲流是連續(xù)介質(zhì)滲流；

2）水流運(yùn)動(dòng)過程中只受重力作用；

3）溫度不會(huì)導(dǎo)致水體發(fā)生相變；

4）不考慮因介質(zhì)摩擦產(chǎn)生的熱量；

5）水流符合達(dá)西定律。

2.2 滲流場與溫度場耦合數(shù)學(xué)模型

不考慮流體與多孔介質(zhì)之間熱平衡所需時(shí)間，即認(rèn)為兩者之間的熱傳遞是瞬時(shí)完成的，忽略熱輻射、化學(xué)反應(yīng)等引起的熱消耗，且流體密度僅和溫度及壓強(qiáng)有關(guān)，得到三維多孔介質(zhì)地下水的水-熱運(yùn)移耦合方程[2］：

溫度場：

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

滲流場：

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

2.3 邊界條件及參數(shù)

2.3.1 邊界條件

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

溫度場方面，當(dāng)?shù)販囟入S時(shí)間周期性變化的傅里葉級數(shù)函數(shù)如下：

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

這里采用陜西省延安市的氣溫?cái)?shù)據(jù)（數(shù)據(jù)來源于中國氣象數(shù)據(jù)網(wǎng)站）并根據(jù)傅里葉函數(shù)來取5 階傅里葉級數(shù)就可以滿足精度要求。擬合得氣溫回歸模型系數(shù)為：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。將系數(shù)帶入上式即可求得溫度隨時(shí)間周期性規(guī)律，溫度為16.37 ℃。

2.3.2 模型參數(shù)

由于本文采用COMSOL 軟件進(jìn)行數(shù)值模擬，而COMSOL中自帶計(jì)算傳熱模塊是用熱量運(yùn)移方程進(jìn)行計(jì)算的，默認(rèn)為導(dǎo)熱系數(shù)模型為體積平均模型，因此用PDE 模塊來代替COMSOL 自帶的傳熱模型進(jìn)行計(jì)算，可以靈活修改導(dǎo)熱系數(shù)經(jīng)驗(yàn)?zāi)Ｐ?，并通過修改相關(guān)的傳熱模塊來實(shí)現(xiàn)不同導(dǎo)熱系數(shù)條件下的土石壩流熱耦合模型求解。模型滲流場及溫度場的計(jì)算參數(shù)[4］，見表1，模型構(gòu)圖見圖1[5］。

表1 土石壩流-熱耦合模型計(jì)算參數(shù)

續(xù)表2

圖1 模型結(jié)構(gòu)示意圖

3 計(jì)算成果分析

生成的結(jié)果見圖2，可以發(fā)現(xiàn)上游壩體基本處于飽和狀態(tài)，帷幕灌漿區(qū)域軸線附近地帶的壩體中間部位處于飽和與非飽和的交界狀態(tài)，下游壩體則是大部分處于非飽和的狀態(tài)。

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

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

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

4 結(jié)論

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