代古寺面板堆石壩滲流場有限體積法數值分析

張明亮，宿曉輝，吳正橋，張建濤，李浩瑾

(1.大連理工大學 建設工程學部水利工程學院，遼寧 大連 116033； 2.中水北方勘測設計研究有限責任公司，天津 300222)

滲流分析和滲流控制是面板堆石壩工程中極為重要的研究內容，國內外有許多由于滲流控制不當而引起破壞或潰決的工程，1993 年青海溝后面板堆石壩因滲透破壞導致大壩潰決[1]，墨西哥的Aguamilpa[2]大壩由于面板上的水平裂縫，滲漏量達到了257.7 L/s。由此可見，面板堆石壩中由面板和防滲帷幕組成的防滲系統對壩體和壩基的安全至關重要。目前大多數滲流數值分析采用有限單元法，張鳳財等[3-4]采用飽和穩定滲流有限元法研究混凝土面板壩在面板產生裂縫后的滲流場分布及滲流量的變化；潘正陽等[5]通過模擬不同深度、類型、位置、大小及分布的面板裂縫，獲得了面板裂縫對大壩滲流量及浸潤線等滲流要素的影響規律；何玲麗等[6]以 Richards 方程、運動波方程和有限元法為基礎，提出了一種降雨時考慮徑流流量補給的滑坡滲流數值模擬方法，正確反映降雨入滲對基質吸力大幅降低的作用。然而，有限元方法對計算機存儲需求高，計算工作量大，應用到大型工程滲流計算時耗時嚴重。有限體積法與有限元法相比，計算工作量小，并且在控制體積內局部絕對守恒的特性對流動問題處理具有天然的優勢。陳國芳等[7]提出了一種修正的中心型有限體積法，避免了標準有限體積方法的 “數值障礙”現象； Manzini 等[8]在非結構網格上用二階精確單元中心有限體積法數值求解二維滲流方程，提高了數值精度；郭影等[9]提出了一種滲流吸水誘發巖體強度弱化的有限體積數值計算方法，實現了滲流作用下孔隙內自由水向基質吸附水轉化過程的模擬，并實現了巖體模量與強度隨滲流時間逐漸弱化過程的模擬。

本文通過自主開發的三維有限體積法滲流計算軟件DUT-Seepage[10]，采用無系數矩陣數值迭代求解方法，對代古寺混凝土面板堆石壩三維復雜滲流場進行滲流特性分析，通過避免數值計算過程中的矩陣操作，減少計算工作量，并采用GPU-OpenACC 并行技術加快求解速度[10]，實現大型工程滲流模擬的快速求解。本文重點研究面板和防滲帷幕聯合作用下的滲流控制效果，并開展基巖滲透敏感性分析，研究基巖滲透系數對滲流量的影響，最后探討壩基擠壓破碎帶對面板堆石壩滲流場特性的影響，為其他存在地質缺陷的工程提供參考。

1 工程概況

代古寺水庫河道走向呈大“S”形，河谷狹窄，呈深切“V”形，為橫向或斜向河谷。兩岸壩肩邊坡巖體均為志留系板巖，均傾向坡內，自然邊坡穩定，整體穩定性較好。壩址區滑坡、泥石流不發育。壩址區未發現區域性斷層通過，構造行跡以裂隙及小規模斷層為主，多傾向坡里，對壩體抗滑穩定影響不大。為了加強抗滲性，沿左右基巖向下進行帷幕灌漿，減少滲漏量。工程擬定正常蓄水位1 804 m，最大壩高151 m。

本工程地質條件復雜，壩址區巖體以物理風化為主，JI38 擠壓破碎帶自壩（0+115）斷面往右岸延伸。壩址區風化程度不均，左、右岸強風化帶厚9.4～26.3 m，弱風化帶厚29.5～66.5 m。右岸巖體強卸荷帶水平深度15.0～18.0 m，弱卸荷帶水平深度大于35.0 m，左岸巖體強卸荷帶水平深度5.0～8.0 m，弱卸荷帶水平深度15.0～20.0 m。大壩最大斷面如圖1 所示。根據壩體填料分區原則，該面板堆石壩分區自上游至下游依次為：上游蓋重區（1B）、上游鋪蓋區（1A）、混凝土面板（F）、墊層料區（2A）、特殊墊層區（2B）、過渡料區（3A）、主堆石區（3B）、次堆石區（3C）和下游壩基排水區（3F）。

圖1 代古寺面板堆石壩最大斷面（單位：m）Fig.1 The maximum cross section of Daigusi concrete face rockfill dam (unit: m)

2 滲流計算方法和原理

三維穩態滲流控制微分方程[10]為：

式中：H為總水頭（m）；直角坐標軸x,y,z為滲透主方向；Kxx,Kyy,Kzz為沿滲透主方向的滲透系數。在本次工程滲流計算分析中，假定壩體各部分材料均為各向同性材料。

本文提出采用基于非結構化網格有限體積法對滲流控制方程進行離散求解，采用預處理共軛梯度法（PCG）+ GPU 并行技術對求解過程進行加速。采用Cell-Vertex（CV）方法形成控制體，并在控制體上對控制方程進行積分求解?？刂企w如圖2 所示，控制體由圍繞節點P的邊中心點與相應的形心點連接而成，此控制體形成方法會在一定程度上降低網格奇異性的影響，保證計算結果的準確性。

圖2 圍繞節點P 的控制體示意Fig.2 Diagram of control volume around node P

采用有限體積法對滲流控制方程進行離散求解：

對方程（2）應用高斯散度定理：

式中：積分區域cv為控制節點i的控制體；積分區域S為控制體的閉合曲面；K=Kxx=Kyy=Kzz為滲透系數；H為總水頭；V為控制體體積；ncell為與控制節點P相關聯的單元數目；ΔSci為控制體在四面體單元i中的對應邊界面。對于給定的四面體單元，可以采用閉合曲面積分公式dS=0計算nΔSci。

根據對封閉控制體所有邊界面法向量積分為0，可以得到：

將式（4）代入式（3）可以得到：

式中：nΔSpi為含有控制節點P的四面體單元中P點所對應的面和法向量的乘積；(?H)i為含有控制節點P的四面體單元的中心梯度；Ki為含有節點P的四面體單元的中心滲透系數，其四面體計算單元如圖2 所示。

四面體單元的中心梯度采用格林公式計算：

四面體單元的中心滲透系數通過體積加權法計算：

式中：Hi為四面體單元頂點i處的變量值；Si為四面體單元頂點i所對應的面與法向量之積；V為四面體單元體積；Km為四面體單元頂點滲透系數；Vm為單元頂點m控制體的體積。

最終得到滲流控制方程的離散形式：

采用預處理共軛梯度法對式（8）進行迭代求解，其求解步驟為：

（1） 給定求解域初始值H0；

（2） 通過初始值H0，求得初始殘差向量r0，并利用預處理矩陣M?1得到p0；

式中：M=diag(a11,a22,···,ann)

（3） 循環迭代計算，當計算收斂誤差小于給定允許值時，循環迭代停止。

3 三維滲流場分析

3.1 滲流計算模型邊界及參數

滲流計算模型邊界類型主要包括三類：（1）水頭邊界，包括壩址區上游水位線以下的上游壩坡、水庫岸坡和庫底及壩址區下游水位線以下的下游壩坡、岸坡和河道；（2）逸出邊界，包括壩體下游坡面在下游水位線以上的區域、壩址區下游水位線以上的左、右岸邊坡面、模型下游側截取邊界面在下游水位線以上的區域；（3）不透水邊界，包括模型上游側截取邊界面、模型左右岸兩側截取邊界面和模型底部截取邊界面。該面板堆石壩的三維滲流計算模型如圖3 所示。

圖3 三維滲流計算模型Fig.3 Three-dimension seepage calculation model

根據滲透及反濾試驗，并結合地質、水文及鉆孔壓水實驗資料，對模型中各介質的滲透系數取值如表1 所示。正常蓄水工況下，壩體上游水位1 804.00 m，下游水位1 675.00 m。

表1 面板堆石壩各材料分區滲透系數Tab.1 Permeability coefficient of each material partition of concrete face rockfill dam單位： cm/s

3.2 三維滲流場計算結果分析

圖4 為滲流計算域的浸潤面，圖5 為沿壩軸線方向孔隙壓力分布，圖6 為壩防滲結構體的總水頭分布。由計算結果可以看出：

圖4 計算域滲流浸潤面Fig.4 Seepage infiltration surface of the calculation domain

圖5 沿壩軸線方向不同典型剖面壓力水頭分布Fig.5 Pressure head distribution of different typical cross sections along dam axis

圖6 面板堆石壩不同部位的總水頭分布Fig.6 Total water head distribution diagram of different partitions of concrete face rockfill dam

（1）面板壩滲流場水頭分布規律合理，總水頭等值線在面板、趾墻、帷幕等處較為密集，上游水頭由這些部位承擔，滲流控制系統具備很好的防滲效果。

（2）在混凝土面板堆石壩內，面板上下游側水位下降最大，上游水流經過混凝土面板后，總水頭由上游正常蓄水位1 804.00 m 降到1 680.50 m 左右，經墊層區、砂礫堆石區至排水區時，總水頭降至 1 675.15 m，經過混凝土面板后，壩體最大斷面逸出點水力坡降為0.008，小于允許水力坡降0.1，滿足正常運行要求，其余各填筑分區水力坡降均小于允許水力坡降。

（3）通過計算壩軸線斷面的流量得到最終壩基滲流及繞壩滲流量，其中通過壩基和帷幕的滲流量為586 m3/d，兩岸繞壩滲流量為82 m3/d。計算結果表明左、右岸滲漏量很小，繞壩滲流并不明顯。

（4）通過GPU-OpenACC 并行技術[10]，對本案例千萬級計算單元網格的滲流計算僅需3～5 h，能夠滿足工程滲流時間要求。

4 基巖滲透系數對滲流量的影響

由于面板堆石壩基巖實際滲透系數難以確定，為了分析基巖滲透系數的變化對面板堆石壩滲流量的影響，本文另外選取了2 組不同的基巖滲透系數，計算得到不同基巖滲透系數對應的滲流量，計算結果見表2。

表2 正常蓄水位下不同基巖滲透系數壩體滲流量Tab.2 Seepage flow with different bedrock permeability coefficients under normal storage level

從表2 可知，隨著基巖滲透系數的增大，左、右岸及壩基滲透量也呈現逐漸增大的趨勢，尤其是當微風化基巖滲透系數增大時，壩基滲透量出現急劇增大現象。

通過與相近壩高的國內外面板堆石壩工程進行滲流量[11-17]比對發現，該面板堆石壩的滲流量處于中等偏低范圍，總體偏于安全，對比結果如圖7 所示。

圖7 相近工程滲流量對比Fig.7 Comparison of engineering seepage flow

5 JI38 擠壓破碎帶影響分析

本工程基巖中存在JI38 擠壓破碎帶，寬度自左岸向右岸逐漸加寬，壩（0+285）斷面的破碎帶寬度較大，約11 m。由于JI38 擠壓破碎帶的滲透系數未知，故選取3 組不同的滲透系數分析JI38 擠壓破碎帶對浸潤線位置的影響，并將計算結果與GEO-Studio 計算結果進行對比。

建立面板堆石壩（0+165）二維滲流分析模型，針對壩（0+165）斷面處存在的JI38 擠壓破碎地質缺陷，在正常蓄水位1 804.00 m、尾水位1 675.00 m 工況下，開展JI38 破碎帶滲透敏感性分析。JI38 擠壓破碎帶的滲透系數分別取1.00×10?1、1.00×10?4和1.00×10?5cm/s，壩體其余部分滲透系數參照表1。

分別采用DUT-Seepage 與GEO-Studio 程序，分析不同工況下JI38 破碎帶滲透系數對壓力水頭和總水頭分布的影響，計算結果見圖8～10。對比可發現DUT-Seepage 與GEO-Studio 的計算結果較一致；不同滲透系數下的大壩浸潤線位置幾乎不變，但當滲透系數為1.00×10?1cm/s 時，壩基的壓力水頭及總水頭分布產生較大變化，并且JI38 破碎帶周邊的水力坡降從0.40 增至0.75，增大了壩基滲透破壞風險，建議對 JI38 破碎帶進行加固處理，降低滲透系數，避免壩基發生滲透破壞。

圖8 不同工況下JI38 破碎帶滲透系數對孔隙水壓力分布的影響Fig.8 The effect of permeability coefficient of JI38 crush zone on pore water pressure distribution under different cases

圖10 不同工況下JI38 破碎帶滲透系數對水力梯度的影響Fig.10 The effect of permeability coefficient of JI38 crush zone on hydraulic gradient distribution under different cases

6 結 語

本文利用三維有限體積法滲流計算分析軟件DUT-Seepage，采用無矩陣數值算法及GPU-OpenACC 并行技術對代古寺混凝土面板堆石壩進行了三維滲流數值模擬分析。三維滲流結果顯示代古寺面板堆石壩工程的防滲系統，包括面板與防滲帷幕有效阻斷了上游水流的滲透，具有較好的防滲效果，堆石壩主體無滲透破壞的風險；當基巖滲透性較大時，大壩整體滲漏量約為2 100 m3/d，處于同類工程的偏低水平；若擠壓破碎帶滲透系數過大時，存在地基滲透破壞風險，建議對壩基破碎帶進行加固處理。