碾壓混凝土重力壩廊道滲流測值分析與有限元數(shù)值模擬

周永紅，徐蘭玉，靳向波，徐淼

（1.山西汾河二庫管理有限公司，山西 太原 030012；2.水利部南京水利水文自動化研究所，江蘇 南京 210012；3.水利部水文水資源監(jiān)控工程技術研究中心，江蘇 南京 210012；4.河海大學，江蘇 南京 210098）

0 引言

碾壓混凝土重力壩廊道滲水在多個工程中被發(fā)現(xiàn)并得到一定程度的研究[1-3]。祁立友等[4]通過對桃林口水庫滲水觀測資料的分析，總結出了大壩滲水變化的一般規(guī)律。袁自立等[5]綜合大壩滲流觀測資料和有限元模型正反分析，得出了石漫灘碾壓混凝土重力壩滲流異常的成因。李榮等[6]通過ANSYS 軟件，采用等效連續(xù)介質模型模擬碾壓混凝土壩的層間滲流特性，探究思林水電站碾壓混凝土壩的滲流規(guī)律。孫亮[7]基于碾壓混凝土的滲流特性，對龍華河碾壓混凝土壩進行了三維滲流場分析。方衛(wèi)華等[8，9]借助COMSOL Multiphysics 軟件進行三維有限元多場耦合分析，評價汾河二庫攔河大壩除險加固效果。本文在已有的研究基礎上，基于汾河二庫廊道的滲流實測數(shù)據(jù)，在考慮帷幕和排水管等重要復雜滲控結構的基礎上，充分考慮壩體實際裂縫影響，研究了滲流場、應力場和溫度場三場耦合在COMSOL的實現(xiàn)方法。

1 工程概況

汾河二庫是一座以防洪為主，兼有供水、發(fā)電、旅游、養(yǎng)殖等綜合效益的大（2）型水利樞紐工程。攔河大壩為碾壓混凝土重力壩，壩頂長227.7 m，最大壩高88.0 m，河床中部設3 個溢流表孔，每孔凈寬12.0 m，堰頂高程為902.0 m。

由于歷史原因，在汾河二庫大壩主體工程竣工時，大壩壩體、壩肩固結和帷幕灌漿并未完成，一到冬季廊道內漏水嚴重。雖然2014年對壩肩、壩基及壩后連續(xù)墻等采取應急除險加固措施，但壩體內存在較長裂縫，廊道滲水現(xiàn)象依然存在，阻礙水庫效益的發(fā)揮。

2 廊道實測滲流數(shù)據(jù)分析

2.1 滲流觀測布置情況

2016年，為了解廊道內出現(xiàn)的異常滲流情況，在廊道內布置了23 個量水堰，用于觀測不同壩段滲流情況。堰上水頭測量采用振弦式4675LV-1傳感器，量程300 mm，精度為0.01%FS～0.02%FS。

2.2 滲流觀測數(shù)據(jù)與分析

2.2.1 滲流量過程及特征分析

壩體870 廊道的4 個測點WE-6～WE-9 在加固后滲流明顯，且呈現(xiàn)季節(jié)規(guī)律性，測點滲流量與氣溫和上游水位的過程線見圖1，測點總體特征值分析見表1。

由圖1和表1可知，測點的滲流量與溫度和上游水位具有較好的相關性，大致與溫度呈負相關，與上游水位呈正相關。滲流量的最大值普遍出現(xiàn)在冬季末至春季初，最小值普遍出現(xiàn)在夏季至秋季。

表1 測點滲流量特征值分析 ml/s

圖1 上游水位氣溫測點滲流量過程線

2.2.2 監(jiān)控模型

1）裂縫滲流的確定性模型

Witherspoon PA 根據(jù)裂隙巖體滲流的立方定律，導出粗糙裂縫中穩(wěn)定流流量Q的公式：

式中：G為常數(shù)，由裂縫的長度、寬度及水流粘滯系數(shù)確定；f為縫面特征系數(shù)，一般為1.04～1.65，當f=1.0 時，即為光滑情況；△H為裂縫上下游水頭差；b為裂隙的等效張開度。

2）統(tǒng)計模型

由于大壩裂縫的幾何和物理性參數(shù)均很難確定，一般不直接采用確定性模型分析大壩裂縫的滲流情況，而是根據(jù)確定性模型的形式，利用水庫水位、氣溫和裂縫滲流量的觀測數(shù)據(jù)建立裂縫滲流的統(tǒng)計模型。假定裂縫的下游水頭為常數(shù)，按下式計算水頭差△H：

式中：H為上游水位；α，H0為待定系數(shù)。

根據(jù)確定性模型的形式，基于水庫2014—2019年水位、氣溫和典型壩段廊道裂縫滲流量的實測數(shù)據(jù)建立裂縫滲流統(tǒng)計模型：

式中：△TF為橫縫寬度；△HF為水平縫寬度；T為觀測日氣溫。

綜合考慮式（1）—（4），按下式計算滲流量Q：

式中：c1，c2，c3為待定系數(shù)。

考慮到水溫和氣溫的關系實際上要比線性關系復雜，增加氣溫T4，多項式如下：

式中：a1～a15為待定系數(shù)。

870 廊道的4 個測點WE-6～WE-9 都具有良好的相關性，將4 個測點的滲流量合計為870 廊道的總滲流量進行回歸分析，擬合多項式的結果殘差為154.045，擬合度R2為0.777，擬合效果較好。

3 有限元模擬及對比分析

為了從理論上論證該典型壩段滲流監(jiān)測結果的準確性，進一步深入研究壩體裂縫、帷幕及壩基壩體排水管等作用下壩段滲流量的變化規(guī)律，利用有限元方法對存在問題的典型溢流壩段進行數(shù)值模擬。

3.1 模型的建立

針對汾河二庫除險加固前后的結構特點和幾何拓撲性質，考慮該壩段結構型式、壩體材料分區(qū)、基巖材料分區(qū)，以及帷幕灌漿、排水、固結灌漿的布置等因素，建立溢流壩段滲流場-應力場-溫度場三場耦合的雙重介質模型[9，10]。模型單元總數(shù)為350 717 個，計算網格如圖2所示。溫度場和滲流場方程采用瞬態(tài)形式，計算時長為720 d，步長為5 d。

圖2 除險加固后溢流壩段三維有限元模型網格剖分示意圖

滲流場計算時，上、下游水位以下為等水頭邊界，水位以上部分為混合邊界；除險加固前，底高程為870.0 m 的廊道為可能溢出邊界，底部廊道為滿水有壓邊界；除險加固后，所有廊道表面均為可能溢出邊界；其他邊界為不透水邊界。在COMSOL中可用透水層邊界模擬混合邊界條件。

位移場、應力場計算時，壩基底部視為固定約束，壩基兩側x方向為連桿約束，壩基及并縫灌漿高程（855.0 m）以下壩體左右岸方向（y方向）兩側為連桿約束。

溫度場計算時，壩基底部和壩基兩側為絕熱邊界條件，上下游水位以下為第一類邊界條件，混凝土與空氣接觸為第三類邊界條件?；鶐r和混凝土都有初始溫度，該地區(qū)地熱溫度為8.0 ℃，混凝土常溫為9.5 ℃。

3.2 計算參數(shù)及邊界條件

溢流壩段模型主要由壩體和地基組成，其材料性質及力學參數(shù)見表2和表3。

表2 壩體及地基材料力學參數(shù)指標

表3 壩體及地基材料熱學參數(shù)指標

3.3 計算工況

基于2018—2019年2年的氣溫及上游水位變化，文中取不考慮溢流壩段縱橫縫的情況為工況1，考慮溢流壩段縱橫縫情況為工況2。

3.4 計算結果及分析

1）三場分析

壩體部分的第一主應力分布、最低溫度分布、最高水位的滲流場水頭分別見圖3～5。在考慮裂縫之后，壩體增加了裂縫滲流引起的壓力，使壩體的第一主應力增加；滲流量的增加引起了壩體內水頭的抬高；應力的抬高及滲流速度的增加，加快了壩體內的熱交換及傳導速率，從而使溫度有了略微的提升。

圖3 壩體工況1 和2 的第一主應力圖（單位：MPa）

圖4 壩體工況1 和2 的溫度等值線圖（單位：℃）

圖5 壩體工況1 和2 的水頭等值線圖（單位：m）

2）滲流量分析

選取870 廊道工況1，2 條件下的滲流量計算結果與上游水位、氣溫進行對比，如圖6～9所示。

由圖6和圖7所示，在水位和氣溫的雙重影響下，有裂縫的滲流量大，工況2 相對于工況1 的變幅較大；由圖8和圖9可知，溫度對于裂縫滲流的影響更大，裂縫引起的滲流量與溫度有著明顯的負相關性。

圖6 滲流量計算值與水位對比

圖7 滲流量與氣溫對比

圖8 裂縫引起的滲流量與水位對比

圖9 裂縫引起的滲流量與氣溫對比

3）計算值與實測值對比分析

采用雙重介質滲流-應力-溫度三場耦合模型計算工況2 條件下870 廊道的滲流量，并與實測滲流量、上游水位及氣溫進行比較，如圖10～11所示。由圖10可知，采用耦合模型計算得到的滲流量與實測滲流量擬合趨勢性較好，但在滲流量最低時吻合性稍低，主要原因是有限元計算采用的溫度值為大氣溫度的擬合值，造成模型計算過程中溫度值的偏差。同時，由圖11可知，氣溫在20 ℃以上時，溫度對實測滲流量影響很小，基本可以忽略，但這在有限元計算結果中不能較好地體現(xiàn)，從而引起了溫度變化帶來的滲流量偏大。分析還表明，在水位和氣溫的共同影響下，壩體滲流和裂縫滲流共同作用引起了廊道滲流量的季節(jié)性變化。

圖10 滲流量實測值與計算值和水位對比

圖11 滲流量實測值與計算值和氣溫對比

4 結語

本文根據(jù)汾河二庫實際環(huán)境，以求解滲流量為目的，建立了水庫溢流壩段雙重介質三場耦合模型，并用COMSOL Multiphyics 模擬有無裂縫的各工作性態(tài)及廊道的滲流情況。該法有效融合了動態(tài)監(jiān)測數(shù)值和有限元計算結果，建立了壩體廊道內滲流量的在線監(jiān)測與計算模型，可模擬長周期內壩體滲流變化情況，為大壩滲流安全性態(tài)評價提供了新的思路和方法。該法模擬結果與實測值吻合度較好，但在滲流量最低時吻合性稍低，后期尚需進一步研究分析。