老嵐水庫土石壩滲流數(shù)值模擬研究

岳慶河，郭清華，張大雷

(煙臺(tái)市城市水源工程運(yùn)行維護(hù)心，山東 煙臺(tái) 264000)

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

1 工程概況

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

2 滲流程序開發(fā)

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

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

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

(1)

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

(2)

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

圖1 樁號(hào)0+325處土石壩剖面圖

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

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

圖2 滲流計(jì)算迭代程序框圖

3 滲流計(jì)算

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

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

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

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

(3)

3.2 不同上下游水位工況計(jì)算

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

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

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

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

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

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

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

滲漏問題是導(dǎo)致土石壩潰壩事件的重要原因，假設(shè)土石壩可能引發(fā)滲漏問題的不利工況，研究其對(duì)滲流場(chǎng)變化影響[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)得到工況⑧。列出了部分計(jì)算圖形如圖7—8所示，求得關(guān)鍵部位滲流水力要素見表3。

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

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

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

4 成果分析

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

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

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

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

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

5 結(jié)語

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