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

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

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

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

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

1 工程概況

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

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

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

2 滲流計(jì)算方法和原理

三維穩(wěn)態(tài)滲流控制微分方程[10]為：

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

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

圖2 圍繞節(jié)點(diǎn)P 的控制體示意Fig.2 Diagram of control volume around node P

采用有限體積法對(duì)滲流控制方程進(jìn)行離散求解：

對(duì)方程（2）應(yīng)用高斯散度定理：

式中：積分區(qū)域cv為控制節(jié)點(diǎn)i的控制體；積分區(qū)域S為控制體的閉合曲面；K=Kxx=Kyy=Kzz為滲透系數(shù)；H為總水頭；V為控制體體積；ncell為與控制節(jié)點(diǎn)P相關(guān)聯(lián)的單元數(shù)目；ΔSci為控制體在四面體單元i中的對(duì)應(yīng)邊界面。對(duì)于給定的四面體單元，可以采用閉合曲面積分公式dS=0計(jì)算nΔSci。

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

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

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

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

四面體單元的中心滲透系數(shù)通過體積加權(quán)法計(jì)算：

式中：Hi為四面體單元頂點(diǎn)i處的變量值；Si為四面體單元頂點(diǎn)i所對(duì)應(yīng)的面與法向量之積；V為四面體單元體積；Km為四面體單元頂點(diǎn)滲透系數(shù)；Vm為單元頂點(diǎn)m控制體的體積。

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

采用預(yù)處理共軛梯度法對(duì)式（8）進(jìn)行迭代求解，其求解步驟為：

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

（2） 通過初始值H0，求得初始?xì)埐钕蛄縭0，并利用預(yù)處理矩陣M?1得到p0；

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

（3） 循環(huán)迭代計(jì)算，當(dāng)計(jì)算收斂誤差小于給定允許值時(shí)，循環(huán)迭代停止。

3 三維滲流場分析

3.1 滲流計(jì)算模型邊界及參數(shù)

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

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

根據(jù)滲透及反濾試驗(yàn)，并結(jié)合地質(zhì)、水文及鉆孔壓水實(shí)驗(yàn)資料，對(duì)模型中各介質(zhì)的滲透系數(shù)取值如表1 所示。正常蓄水工況下，壩體上游水位1 804.00 m，下游水位1 675.00 m。

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

3.2 三維滲流場計(jì)算結(jié)果分析

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

圖4 計(jì)算域滲流浸潤面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）面板壩滲流場水頭分布規(guī)律合理，總水頭等值線在面板、趾墻、帷幕等處較為密集，上游水頭由這些部位承擔(dān)，滲流控制系統(tǒng)具備很好的防滲效果。

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

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

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

4 基巖滲透系數(shù)對(duì)滲流量的影響

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

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

從表2 可知，隨著基巖滲透系數(shù)的增大，左、右岸及壩基滲透量也呈現(xiàn)逐漸增大的趨勢(shì)，尤其是當(dāng)微風(fēng)化基巖滲透系數(shù)增大時(shí)，壩基滲透量出現(xiàn)急劇增大現(xiàn)象。

通過與相近壩高的國內(nèi)外面板堆石壩工程進(jìn)行滲流量[11-17]比對(duì)發(fā)現(xiàn)，該面板堆石壩的滲流量處于中等偏低范圍，總體偏于安全，對(duì)比結(jié)果如圖7 所示。

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

5 JI38 擠壓破碎帶影響分析

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

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

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

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

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

6 結(jié) 語

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