岳慶河，郭清華，張大雷

(煙臺市城市水源工程運行維護心，山東 煙臺 264000)

老嵐水庫樞紐工程壩型為混凝土重力壩與黏土心墻土石壩結(jié)合的混合壩，使用ANSYS有限元分析軟件APDL語言編寫程序?qū)浖M行功能擴充，得到土石壩滲流場有限元分析程序，可以方便快捷地進行滲流數(shù)值模擬分析，研究老嵐水庫土石壩滲流特性，為滲流控制提供參考。

1 工程概況

老嵐水庫壩址位于煙臺市福山區(qū)外夾河中游，控制流域面積624km2，總庫容15792m3，興利水位44.8m，興利庫容8491萬m3，死水位34m，死庫容513萬m3，屬大(2)型水庫，樞紐工程于2020年10月開工，計劃2023年具備蓄水條件。工程建設(shè)任務(wù)以城市供水、防洪為主，兼顧灌溉和改善下游生態(tài)環(huán)境等綜合利用。壩址區(qū)屬丘陵河谷地貌，呈不對稱的寬U形河谷，兩岸壩肩地形起伏不大，邊坡巖體主要為花崗斑巖和二長花崗巖。壩型為混凝土重力壩與黏土心墻壩結(jié)合的混合壩，重力壩段布置于左岸及主河槽，黏土心墻壩段布置于右岸。黏土心墻砂殼壩壩頂高程50.5m，壩頂長439m，壩頂寬度8m，最大壩高27.5m，黏土心墻頂高程48.5m，頂寬3m，如圖1所示。

2 滲流程序開發(fā)

2.1 開發(fā)基礎(chǔ)

ANSYS提供了APDL參數(shù)化設(shè)計語言，允許用戶開發(fā)專用有限元分析程序，擴展了ANSYS使用范圍，二次開發(fā)滲流分析程序可充分利用ANSYS可視化的前后處理能力，方便了建模、模型檢察、結(jié)果查看及對結(jié)果的計算分析，避免了重復(fù)編制有限元程序的繁鎖[1]。通過對比分析，滲流場與溫度場的熱傳導(dǎo)都遵循質(zhì)量(或能量)守恒定律，滲流微分方程與熱傳導(dǎo)微分方程形式完全相同，參數(shù)及邊界條件相互對應(yīng)，為應(yīng)用ANSYS熱傳導(dǎo)模塊二次開發(fā)滲流場分析程序奠定了基礎(chǔ)[2- 5]。

不考慮土骨架變形及水壓縮性時，飽和—非飽和滲流微分控制方程為：

(1)

以溫度T為控制變量的熱傳導(dǎo)微分方程為：

(2)

式中，λx、λy—x、y方向的熱傳導(dǎo)系數(shù)，W/(m·℃)；c′—比熱系數(shù)，J/(kg·℃)；ρ—介質(zhì)密度，kg/m3。

圖1 樁號0+325處土石壩剖面圖

2.2 程序結(jié)構(gòu)

滲流計算時非飽和土滲透系數(shù)、浸潤線位置、逸出面邊界條件等需在計算過程中迭代確定[6- 7]。程序設(shè)計為兩層迭代結(jié)構(gòu)，內(nèi)層迭代確定非飽和土滲透系數(shù)，計算時保持邊界條件不變，滲透系數(shù)根據(jù)單元孔隙水壓力不斷調(diào)整，直到前后2次計算各單元孔壓變化差小于某一定值時計算收斂；外層迭代確定滲流逸出點位置，調(diào)整逸出面邊界條件范圍，滲流向外流通時計算收斂。完成計算后，利用ANSYS后處理模塊可以直觀的查看浸潤線位置、滲透比降、逸出點位置等信息，通過對下游坡節(jié)點熱流通量求和得到單寬滲漏量值。二次開發(fā)的滲流程序框圖如圖2所示。

圖2 滲流計算迭代程序框圖

3 滲流計算

3.1 網(wǎng)格劃分及滲透系數(shù)確定

在ANSYS中采用PLANE55單元，壩體部位取0.2m、壩基取0.5m進行網(wǎng)格劃分，共劃分119433個單元，239940個節(jié)點，壩體筑壩材料和壩基飽和滲透系數(shù)見表1。

表1 各土層飽和滲透系數(shù)

非飽和土相對滲透系數(shù)采用Van Genuchten-Mualem公式確定[8]。

(3)

3.2 不同上下游水位工況計算

計算3種設(shè)計工況，分別為：①上游興利水位44.8m、下游相應(yīng)最低水位32.4m；②上游設(shè)計洪水位46.27m、下游相應(yīng)水位34.82m；③上游校核洪水位48.35m、下游相應(yīng)水位36.88m。為了對比分析，增加2種假設(shè)工況，分別為：④上游校核洪水位48.35m，下游水位34.82；⑤上游校核洪水位48.35m，下游水位32.4。列出了部分計算圖形如圖3—6所示，求得關(guān)鍵部位滲流要素見表2。

圖3 工況③總水頭等值線圖

圖4 工況③孔隙水壓力等值線圖

圖5 工況③下游坡滲透比降圖

圖6 工況③壩基中粗砂層滲透比降圖

表2 不同水位組合下關(guān)鍵部位滲流要素表

3.3 假設(shè)不利工況計算

滲漏問題是導(dǎo)致土石壩潰壩事件的重要原因，假設(shè)土石壩可能引發(fā)滲漏問題的不利工況，研究其對滲流場變化影響[9- 10]。在工況②的基礎(chǔ)上，分別假設(shè)防滲帷幕防滲效果差(調(diào)整帷幕滲透系數(shù)為1×10-7m/s)得到工況⑥、心墻及截水槽防滲效果差(調(diào)整心墻滲透系數(shù)為2.9×10-7m/s)得到工況⑦、反濾排水體失效(調(diào)整反濾排水體滲透系數(shù)為3.5×10-5m/s)得到工況⑧。列出了部分計算圖形如圖7—8所示，求得關(guān)鍵部位滲流水力要素見表3。

圖7 工況⑥總水頭等值線圖

表3 假設(shè)不利工況下關(guān)鍵部位滲流水力要素表

圖8 工況⑦總水頭等值線圖

4 成果分析

(1)各種工況下游壩坡浸潤線均與下游水位平齊，黏土心墻降低浸潤線高度效果顯著，出逸點位置均在下游水位線處，無自由出逸面；根據(jù)地質(zhì)勘查成果，壩殼砂礫料允許滲透比降0.45、壩基粗砂層允許滲透比降0.22，出逸點及壩基砂滲透比降均未超過允許值。

(2)滲透比降大值集中在心墻及防滲帷幕中，壩殼砂料及壩基其它區(qū)域滲透比降較??；下游壩坡滲透比降最大點發(fā)生在出逸點位置，壩基砂層滲透比降最大點發(fā)生在壩坡與壩基接觸點處。

(3)下游反濾排水體失效對總水頭分布影響不大，但壩基砂層滲透比降增大，說明此時部分滲流繞經(jīng)壩基排出，下游壩殼及壩基材料良好的排水性有效防止了浸潤線抬升。

(4)隨著上下游水位差增大，滲漏量也逐步增大；截滲墻失效后，總水頭線向兩側(cè)擴散，導(dǎo)致兩側(cè)地基中滲透比降增大；黏土心墻及截滲墻失效后單寬滲漏量增加明顯，說明兩個部位是控制滲漏的關(guān)鍵部位，施工過程中重點加強質(zhì)量控制，建成運行后進行重點監(jiān)測。

(5)心墻土石壩斷面設(shè)計合理，正常工況滲流水力要素計算結(jié)果符合要求，部分壩體區(qū)域失效對壩體整體的滲透安全影響不大，說明設(shè)計斷面適應(yīng)性好，安全保障程度高。

5 結(jié)語

基于ANSYS熱分析模塊二次開發(fā)得到的滲流程序，能夠很好實現(xiàn)土石壩滲流場數(shù)值模擬，可廣泛應(yīng)用于土石壩斷面設(shè)計、優(yōu)化及相關(guān)研究；對老嵐水庫土石壩滲流特性的研究成果，為工程建設(shè)期質(zhì)量控制及建成后滲流安全監(jiān)測提供指導(dǎo)。但由于地質(zhì)條件的復(fù)雜性，所建立土石壩滲流模型及相關(guān)參數(shù)特別是非飽和土相對滲透系數(shù)確定是否合理，應(yīng)在工程建成后與實際滲流觀測資料對比、分析，根據(jù)實測資料完善模型、優(yōu)化參數(shù)，還可進一步研究建立三維模型對老嵐水庫土石壩滲流特性做更進一步研究。