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

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

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

0 引言

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

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

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

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

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

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

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

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

1.1 滲流模型

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

(1)

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

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

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

(2)

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

1.2 溶質(zhì)遷移模型

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

(3)

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

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

(4)

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

(5)

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

1.3 帷幕體溶蝕模型

Ca(OH)2、水化硅酸鈣(CSH)等水化產(chǎn)物的溶解流失會導致帷幕體孔隙率增加。由于CSH凝膠溶解過程復雜，為簡化模型，將帷幕體中CSH凝膠溶蝕過程視為Ca(OH)2的溶蝕[14-15]。隨Ca(OH)2溶解，帷幕體孔隙率持續(xù)增大，其增長速率與帷幕體中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的函數(shù)，如式(7)～(8)所示[17]：

(7)

ksp=0.012 5×109e-0.019T

(8)

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

1.4 耦合作用

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

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

2 工程實例分析

2.1 工程概括

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

圖2 輸水系統(tǒng)計算模型Fig.2 Calculation model of water conveyance system

2.2 有限元模型

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

圖3 模型有限元網(wǎng)格Fig.3 Finite element mesh of model

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

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

表1 水文地質(zhì)參數(shù)Table 1 Hydrogeological parameters

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

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

2.4 計算結果驗證

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

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

(9)

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

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

3 結果分析與討論

3.1 滲流場運動規(guī)律

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

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

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

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

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

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

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

3.2 溶質(zhì)遷移場運移規(guī)律

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

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

圖8 運轉(zhuǎn)時間內(nèi)計算區(qū)域Ca2+濃度分布Fig.8 Distribution of calcium ion concentration in calculation area during operation time

為進一步探究幕后Ca2+擴散趨勢，在幕后同一水平方向調(diào)取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計算曲線可知，此類較遠區(qū)域需經(jīng)過一段時間后，Ca2+才能完成擴散和遷移，Ca2+濃度開始增大。最終d2、d3Ca2+濃度與d1類似，呈緩慢下降趨勢。

3.3 帷幕體微觀結構分布規(guī)律

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

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

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

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

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

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

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

4 結論

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

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

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

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