高壓輸水隧洞帷幕體防滲性能的數值模擬*

余 濤，章 光，高培培，李墨瀟，胡少華,3

(1.武漢理工大學 安全科學與應急管理學院，湖北 武漢 430070；2.湖北省水利水電規劃勘測設計院，湖北 武漢 430070；3.國家大壩安全工程技術研究中心，湖北 武漢 430010)

0 引言

為合理開發利用水能資源，我國建設1批大型、特大型水利水電工程，規模巨大，水文地質復雜，其中各類高壩普遍達200～300 m[1-2]。在高壩運行過程中，壩基帷幕體承受滲透壓力高達3 MPa，防滲排水系統之間巖體滲透坡降超過30，滲漏問題極為突出[3]。由于高壓輸水隧洞帷幕體長期處于高水頭、高滲壓環境，水與巖石、帷幕體相互作用發生溶蝕，使帷幕體防滲性能減弱甚至失效，導致壩基揚壓力升高、滲流量增加，對大壩穩定性和安全性造成嚴重影響[4-5]。針對高水頭下帷幕體性能演化帶來的具體問題提出合理優化治理方案，對水電站安全穩定運行意義重大。

近年來，帷幕體防滲性能研究主要集中在3個方面：1)通過帷幕體不同性態參數模擬對壩基滲流場的影響：高明軍等[6]通過建立2維滲流模型，分析帷幕體體型對壩基滲流的影響；陳書燕等[7]通過軟件模擬不同帷幕體厚度對滲流場的影響。2)利用水文地質參數對帷幕體時效性進行評價：霍吉祥等[8]通過幕后地下水質特征，采用水化學圖示法判定壩基帷幕體防滲性能及可能存在的防滲缺陷；劉勝等[9]通過壩基帷幕體滲透系數，定量分析帷幕體防滲可靠性。3)從水文地球化學基本理論出發，建立多物理場耦合模型分析帷幕體防滲效果：彭鵬等[10]依據連續介質動力學及化學動力學理論建立水-巖-帷幕體相互作用的多場耦合模型，定性及定量評價帷幕體防滲性能；張開來等[11]基于化學動力學原理建立滲透-耦合模型，確定帷幕體耐久性控制指標，評價帷幕體防滲性能。但目前針對復雜環境下帷幕體防滲性能研究相對較少，尤其對處于高水頭、高滲壓條件下帷幕體時效性的研究十分薄弱。

為研究復雜環境下帷幕體防滲性能衰減動態過程，本文將滲流模型、溶質運移模型與帷幕體溶蝕模型耦合，建立數值分析模型，采用有限元法得到可視化帷幕體微觀結構隨時間變化過程，并通過工程數據進行驗證，研究結果可為后續帷幕體時效性研究提供依據。

1 高壓輸水隧洞帷幕體多物理場耦合模型

高壓輸水隧洞帷幕體長期處于高水力梯度、高滲壓環境，水與帷幕體發生溶蝕作用，使內部溶質離子不斷析出，其中Ca2+是決定水泥基材料性能的重要部分，Ca2+析出會導致帷幕體微觀結構受損，防滲作用減弱甚至失效。本文基于連續介質動力學、溶質遷移及化學動力學原理建立多物理場耦合模型，表征高壓輸水隧洞帷幕體性能演變動態過程，并作如下3種假設：

1)帷幕體與周邊巖體均為等效連續介質，帷幕體與周圍巖體孔隙均勻分布在整個滲流區域內，并且服從達西定律。

2)帷幕體與周邊巖體為開放體系，與周圍環境存在Ca2+交換與能量交換。

3)地下水中Ca2+濃度改變，全部由帷幕體溶蝕擴散引起，不考慮周邊巖體溶解和沉淀引起的變化。

1.1 滲流模型

帷幕體與周邊巖體為等效連續介質，地下水在多孔介質內流動用達西定律描述，如式(1)所示：

(1)

式中：εp為孔隙率；P為壓力，Pa；t為時間，s；?為梯度；ρ為密度，kg/m3；u為滲流流速，m/s；Qm為源匯項，kg/(m3·s-1)；k為滲透系數，m/s。

帷幕體滲透系數本質受微觀結構影響，采用Kozeny-Carman方程表達帷幕體滲透系數與孔隙率間關系[12]，如式(2)所示：

k/k0=φ3(1-φ02)/φ03(1-φ2)

(2)

式中：k0和k分別為初始狀態與不同反應時刻下帷幕體滲透系數，m/s；φ0和φ分別為初始狀態與不同反應時刻帷幕體孔隙率。

1.2 溶質遷移模型

水溶液中溶質遷移及因水與帷幕體發生化學反應而引起的溶質濃度變化,可通過Fick第二定律描述,如式(3)所示：

(3)

式中：C為地下水中溶質濃度,mol/m3；φ為孔隙率；u為達西流速,m/s；D為水動力彌散張量，m2/s；Rc為溶液中Ca2+增加速率，mol/(m3·s)。

彌散系數張量D主要受達西流速u影響，函數關系滿足式(4)～(5)[13]：

(4)

τ=-1.5tan{8.0(φ-0.25)}+2.5

(5)

式中：δ為Kronecker函數；αL、αT分別為橫向彌散系數和縱向彌散系數，m2/s；Dm為Ca2+在水中的擴散系數，mol/m3;τ為多孔介質曲折度。

1.3 帷幕體溶蝕模型

Ca(OH)2、水化硅酸鈣(CSH)等水化產物的溶解流失會導致帷幕體孔隙率增加。由于CSH凝膠溶解過程復雜，為簡化模型，將帷幕體中CSH凝膠溶蝕過程視為Ca(OH)2的溶蝕[14-15]。隨Ca(OH)2溶解，帷幕體孔隙率持續增大，其增長速率與帷幕體中Ca(OH)2溶解速率和摩爾體積有關[16]，如式(6)所示:

(6)

式中：Ms為Ca(OH)2的摩爾體積，cm3/mol；Rs為固相Ca(OH)2溶解速率，mol/(m3·s)。

固相Ca(OH)2溶解速率Rs是溶液中Ca2+濃度C的函數，如式(7)～(8)所示[17]：

(7)

ksp=0.012 5×109e-0.019T

(8)

式中：kd為溶解速率系數；ksp為固相Ca(OH)2的溶解積，與環境溫度T有關。溶液中Ca2+增加速率與帷幕體中固相Ca(OH)2溶解速率數值相等，即|Rc|=|Rs|。

1.4 耦合作用

通過構建數值分析模型并采用有限元法進行求解,以實現高壓輸水隧洞帷幕體防滲性能衰減數值模擬，模型由反映滲流場的滲流模型、反映化學場溶質遷移模型及反映帷幕體微觀結構的溶蝕模型組成。在高水頭的長期侵蝕和溶解作用下，帷幕體微觀結構改變，孔隙率增加，滲透系數增大，從而導致流速增大，加速Ca2+析出，溶蝕過程加快，最終導致帷幕體抗滲性和耐久性持續衰減。模型組成及耦合關系如圖1所示。

圖1 模型組成及耦合關系Fig.1 Model composition and coupling relationship

2 工程實例分析

2.1 工程概括

福建仙游抽水蓄能電站裝機容量1 200 MW，輸水系統如圖2所示。輸水系統由高壓岔管、排水廊道、斷層、節理密集帶及帷幕體組成。斷層F22與節理密集帶為軟結構面，兩者傾角陡且在高壓岔管上方相交，形成滲流通道。高壓岔管分岔前為鋼筋混凝土襯砌型式，分岔后采用壓力鋼管，模型整體長230 m，高150 m，帷幕體垂直布置于高壓水道和節理密集帶附近，長4.6 m，高70 m；壩基巖石主要為花崗斑巖，斑狀結構，塊狀構造，屬于堅硬巖類，其滲透系數均勻。

圖2 輸水系統計算模型Fig.2 Calculation model of water conveyance system

2.2 有限元模型

對模型網格進行精細劃分，并對邊界區進行加密處理，共生成24 153個單元，1 116個邊界元，62個頂點單元，平均單元質量0.931 5，網格質量全局逼近于最佳質量1，具體有限元網格如圖3所示。以1 a為時間步長，30 a為總時間，模擬高壓輸水隧洞帷幕體防滲性能衰減動態過程。

圖3 模型有限元網格Fig.3 Finite element mesh of model

2.3 計算參數及邊界條件

以仙游抽水蓄能電站輸水系統作為研究對象，模型水文地質參數見表1。模型上下游邊界條件及反應參數見表2。

表1 水文地質參數Table 1 Hydrogeological parameters

表2 模型上下游邊界條件及反應參數Table 2 Boundary conditions and reaction parameters of model

高壓岔管兩側取正常運行條件下水頭值703 m，右側垂直邊界為下游邊界，節理密集帶與斷層因其滲透系數大，地質構造復雜等特點，與排水廊道共同定義為潛在溢出型邊界，其余為隔水邊界；在溶質遷移場過程分析中，高壓岔管內部定義為化學場入口邊界，高壓水道內部水經檢測Ca2+濃度為0.001 mol/L，排水廊道、裂隙F22與節理密集帶定義為化學場出口邊界。根據滲流監測資料以及現場安全檢查資料，利用滲流監測數據進行反演分析以獲取材料的滲透系數[18]，初始孔隙率根據經驗數據給出。

2.4 計算結果驗證

幕后孔隙水壓監測值與模擬值對比如圖4所示。選取2處幕后孔隙水壓監測點P1(14，213.5)，P3(28，213.65)，計算從2014年至2019年間上游水位最接近模型水位的2個監測點的揚壓力折減系數，如式(9)所示：

圖4 揚壓力折減系數演變曲線Fig.4 Evolution curves of uplift pressure reduction coefficient

(9)

式中：α為揚壓力折減系數；hi為第i測點實測水位，m；hu為上游水位，m；hx為測點處基巖高程，m。

由圖4可知，監測數據揚壓力折減系數為負值，是因為仙游抽水蓄能電站運行期短，壩基巖石滲透性較小，帷幕體防滲性能較好；監測數據揚壓力折減系數呈上升趨勢，與仿真結果規律一致。監測數據在仿真結果曲線上下波動，仿真結果曲線與監測數據擬合曲線吻合較好，表明該模型能真實反應帷幕體防滲性能衰減過程。

3 結果分析與討論

3.1 滲流場運動規律

在已建立多物理場耦合模型基礎上，對工程計算模型進行離散化處理，通過軟件模擬分析揭示帷幕體滲流場隨時間變化規律。帷幕體滲透系數隨運行時間變化如圖5所示。抽水蓄能電站輸水系統30 a后帷幕體滲透系數變化如圖6所示。

圖5 帷幕體滲透系數隨運行時間變化Fig.5 Variation of permeability coefficient of curtain body with running time

圖6 30 a后帷幕體不同高程滲透系數變化Fig.6 Variation Law of permeability coefficient at different elevations of curtain body after 30 a

由圖5～6可知，帷幕體滲透系數隨運行時間逐漸增大，運行30 a后帷幕體滲透系數增大至1.95×10-8m/s，是初始滲透系數(1.0×10-8m/s)的1.95倍;帷幕體上游滲透系數高于下游，這與孔隙溶液運移作用有關。30 a后高程168 m處滲透系數演變曲線呈“U”型，說明帷幕體底端溶蝕呈環繞型溶蝕特征。

帷幕體滲透系數變化必然引起滲流形態變化，為對比30 a后帷幕體內滲流速度變化，在帷幕體中心線位置取中部上下段2點(7，204)、(7，199.5)，底部和頂部節理密集帶下側2點(7，168)、(7，228)進行分析，如圖7所示。

圖7 30 a后帷幕體內不同位置地下水流速變化Fig.7 Variation of groundwater velocity at different positions in the curtain body after 30 a

由圖7可知，隨滲透系數增大帷幕體內中心線位置各點流速逐漸增大，中部2點流速增長較快，約為初始流速2.5倍，說明帷幕體中部兩側微觀結構受損嚴重，防滲性能減弱；底部流速整體呈上升趨勢，但增長較為平緩；頂部節理密集帶下側流速增大幅度不大，趨勢亦不明顯。

3.2 溶質遷移場運移規律

帷幕體內水泥結石在高水頭，高滲壓環境下易被溶解流失，導致帷幕體強度和防滲性能降低。帷幕體抗滲性及耐久性劣化，取決于帷幕體中水化產物溶解過程。通過分析Ca2+濃度變化，能從空間及時間上定量研究帷幕體溶蝕程度。

運轉時間內計算區域Ca2+濃度分布，如圖8所示。由圖8可知，Ca2+帷幕體下游側濃度大于上游側濃度，由于帷幕體上、下部水泥結石逐漸溶解生成Ca2+并向兩側彌散，在滲流場作用下，Ca2+逐漸向下游側聚集擴散，最后進入節理密集帶及排水廊道。由于地下水滲流方向為右上方45°，影響帷幕體上部滲流侵蝕產生的橫向彌散與縱向彌散，導致Ca2+在帷幕體上部下游側匯集，最終使帷幕體兩側Ca2+濃度分布差異性增大，帷幕體頂端靠近排水廊道側位置Ca2+均處于高濃度狀態，然后經排水廊道排出。結果表明，上游側Ca2+濃度變化是由彌散作用引起，下游側Ca2+濃度變化是由對流作用引起，且對流作用大于彌散作用。

圖8 運轉時間內計算區域Ca2+濃度分布Fig.8 Distribution of calcium ion concentration in calculation area during operation time

為進一步探究幕后Ca2+擴散趨勢，在幕后同一水平方向調取3個坐標點d1(40，234.65)、d2(50，234.65)和d3(60，234.65)，各位置點Ca2+濃度變化如圖9所示。

圖9 30 a后幕后Ca2+濃度分布Fig.9 Distribution of calcium ion concentration behind curtain after 30 a

點d1計算曲線接近防滲帷幕頂端位置Ca2+濃度增加較快，后呈緩慢下降趨勢；由點d2、d3計算曲線可知，此類較遠區域需經過一段時間后，Ca2+才能完成擴散和遷移，Ca2+濃度開始增大。最終d2、d3Ca2+濃度與d1類似，呈緩慢下降趨勢。

3.3 帷幕體微觀結構分布規律

模擬結束后，帷幕體孔隙率分布如圖10所示。隨溶蝕時間增加，帷幕體中水化產物不斷溶解流失，孔隙率逐漸增大，但這種變化具有時空變異特性。由圖10可知，在早期階段，整個帷幕體孔隙率基本相同，增幅較小；在20 a運轉時間后，帷幕體不同部位出現不同程度溶蝕；隨溶蝕時間時間增加，差異化越來越明顯。受滲流場影響，在帷幕體頂端部位下方約6.4 m處，滲流場水溶液一部分從防滲帷幕流出到周圍圍巖，另一部分直接匯入周圍節理密集帶，因此頂端部位未受侵蝕，帷幕體性能較好。帷幕體溶蝕主要發生在以下2個位置：

圖10 帷幕體模擬時段孔隙率分布Fig.10 Porosity distribution of curtain body in simulation period

1)帷幕體中部兩側孔隙率較大，溶蝕現象明顯，是帷幕體防滲性能衰退最快區域。帷幕體上游側溶蝕面積明顯大于下游側，說明帷幕體溶蝕可能和地下水壓力、滲流速度有關，帷幕體中部接近高壓水道，其水頭壓力高于帷幕體上部和底部水頭壓力；靠近上游側動水壓力大于下游側動水壓力。

2)帷幕體底部溶蝕較明顯，溶蝕面積相對于中部較小，是一種“環繞型”溶蝕，其特征是溶蝕由表及里，帷幕體與地下水接觸，從表面開始逐漸溶解，隨之被水侵蝕脫落，內部開始發生溶蝕反應。造成這種現象的主要原因跟高壓水繞過帷幕體的滲流運動有關。

為分析帷幕體內孔隙率分布，分別選取帷幕內x為6.1，8.1，10.1 m處作為帷幕上游側、中心線和下游側典型斷面，模擬時段結束3條斷面處孔隙率分布如圖11所示。

圖11 30 a后防滲帷幕典型斷面孔隙率分布Fig.11 Porosity distribution at typical section of anti-seepage curtain after 30 a

由圖11可知,帷幕體頂部孔隙率基本無變化，無論是上游側、中心線還是下游側，孔隙率皆接近0.08，表明此部位Ca(OH)2幾乎未發生溶解反應；帷幕中部，孔隙率相對其他部位最大，在x=6.1 m上游側斷面處，出現“陡坎”式激增,這是由于地下水滲流方向為右上方45°，該角度對應帷幕體上游側滲透坡降較大，溶蝕程度加重，孔隙率相對較大,帷幕體中部上游側至下游側孔隙率逐漸減小；帷幕底部，帷幕體上游側孔隙率相對其他位置最大，而下游側孔隙率比中心線位置大，并最終趨于平衡。

4 結論

1)基于連續介質動力學、溶質遷移及化學動力學原理，建立多物理場耦合模型，研究高壓輸水隧洞帷幕體在實際運轉過程在滲流場、化學場及微觀結構耦合作用下溶蝕破壞規律，能在時空真實反應高壓輸水隧洞帷幕體物理化學特性變化過程。

2)高壓輸水隧洞帷幕體防滲效果隨時間增加逐漸減弱，主要因為帷幕體中諸如Ca(OH)2，水化硅酸鈣(CSH)等可溶性物質溶失所致，可溶性物質溶解導致帷幕體滲透系數增大，孔隙率增加。

3)地下水滲流方向影響帷幕體滲流侵蝕產生Ca2+的橫向彌散與縱向彌散過程，是導致Ca2+遷徙并富集的主要影響因素；帷幕體侵蝕破壞主要發生在帷幕體中部與高壓岔管接觸部位，而底部溶蝕主要以由表面向內部溶蝕為主，其溶蝕面積小于中部兩側位置。

4)通過與水電站監測資料及滲流數據進行對比，模擬計算結果符合其自然規律，可將此模型應用于其他工況中，通過與實際各個監測量進行對比，評估帷幕體抗滲性及耐久性。