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

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

0 引言

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

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

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

1 高壓輸水隧洞帷幕體多物理場(chǎng)耦合模型

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

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

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

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

1.1 滲流模型

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

(1)

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

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

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

(2)

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

1.2 溶質(zhì)遷移模型

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

(3)

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

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

(4)

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

(5)

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

1.3 帷幕體溶蝕模型

Ca(OH)2、水化硅酸鈣(CSH)等水化產(chǎn)物的溶解流失會(huì)導(dǎo)致帷幕體孔隙率增加。由于CSH凝膠溶解過程復(fù)雜，為簡化模型，將帷幕體中CSH凝膠溶蝕過程視為Ca(OH)2的溶蝕[14-15]。隨Ca(OH)2溶解，帷幕體孔隙率持續(xù)增大，其增長速率與帷幕體中Ca(OH)2溶解速率和摩爾體積有關(guān)[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有關(guān)。溶液中Ca2+增加速率與帷幕體中固相Ca(OH)2溶解速率數(shù)值相等，即|Rc|=|Rs|。

1.4 耦合作用

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

圖1 模型組成及耦合關(guān)系Fig.1 Model composition and coupling relationship

2 工程實(shí)例分析

2.1 工程概括

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

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

2.2 有限元模型

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

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

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

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

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

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

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

2.4 計(jì)算結(jié)果驗(yàn)證

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

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

(9)

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

由圖4可知，監(jiān)測(cè)數(shù)據(jù)揚(yáng)壓力折減系數(shù)為負(fù)值，是因?yàn)橄捎纬樗钅茈娬具\(yùn)行期短，壩基巖石滲透性較小，帷幕體防滲性能較好；監(jiān)測(cè)數(shù)據(jù)揚(yáng)壓力折減系數(shù)呈上升趨勢(shì)，與仿真結(jié)果規(guī)律一致。監(jiān)測(cè)數(shù)據(jù)在仿真結(jié)果曲線上下波動(dòng)，仿真結(jié)果曲線與監(jiān)測(cè)數(shù)據(jù)擬合曲線吻合較好，表明該模型能真實(shí)反應(yīng)帷幕體防滲性能衰減過程。

3 結(jié)果分析與討論

3.1 滲流場(chǎng)運(yùn)動(dòng)規(guī)律

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

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

帷幕體滲透系數(shù)變化必然引起滲流形態(tài)變化，為對(duì)比30 a后帷幕體內(nèi)滲流速度變化，在帷幕體中心線位置取中部上下段2點(diǎn)(7，204)、(7，199.5)，底部和頂部節(jié)理密集帶下側(cè)2點(diǎn)(7，168)、(7，228)進(jìn)行分析，如圖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)中心線位置各點(diǎn)流速逐漸增大，中部2點(diǎn)流速增長較快，約為初始流速2.5倍，說明帷幕體中部兩側(cè)微觀結(jié)構(gòu)受損嚴(yán)重，防滲性能減弱；底部流速整體呈上升趨勢(shì)，但增長較為平緩；頂部節(jié)理密集帶下側(cè)流速增大幅度不大，趨勢(shì)亦不明顯。

3.2 溶質(zhì)遷移場(chǎng)運(yùn)移規(guī)律

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

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

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

為進(jìn)一步探究幕后Ca2+擴(kuò)散趨勢(shì)，在幕后同一水平方向調(diào)取3個(gè)坐標(biāo)點(diǎn)d1(40，234.65)、d2(50，234.65)和d3(60，234.65)，各位置點(diǎn)Ca2+濃度變化如圖9所示。

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

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

3.3 帷幕體微觀結(jié)構(gòu)分布規(guī)律

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

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

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

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

為分析帷幕體內(nèi)孔隙率分布，分別選取帷幕內(nèi)x為6.1，8.1，10.1 m處作為帷幕上游側(cè)、中心線和下游側(cè)典型斷面，模擬時(shí)段結(jié)束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ā)生溶解反應(yīng)；帷幕中部，孔隙率相對(duì)其他部位最大，在x=6.1 m上游側(cè)斷面處，出現(xiàn)“陡坎”式激增,這是由于地下水滲流方向?yàn)橛疑戏?5°，該角度對(duì)應(yīng)帷幕體上游側(cè)滲透坡降較大，溶蝕程度加重，孔隙率相對(duì)較大,帷幕體中部上游側(cè)至下游側(cè)孔隙率逐漸減小；帷幕底部，帷幕體上游側(cè)孔隙率相對(duì)其他位置最大，而下游側(cè)孔隙率比中心線位置大，并最終趨于平衡。

4 結(jié)論

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

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

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

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