基于對(duì)流-擴(kuò)散理論的防滲帷幕溶質(zhì)溶蝕運(yùn)移分析

蔣明東，陶 駿，李永豐，任 鑫，李 驥，張貴金

(1.湖南平江抽水蓄能有限公司，湖南 平江 410400；2.長沙理工大學(xué)，湖南 長沙 410014；3.中國電建集團(tuán)中南勘測(cè)設(shè)計(jì)研究院有限公司，湖南 長沙 410014)

1 概論

帷幕灌漿是壩基防滲的常用方法，但由于施工難度，形成的帷幕體或多或少存在蜂窩、孔洞、孔隙等缺陷或薄弱段，在水頭作用或長期運(yùn)行后很容易產(chǎn)生滲漏通道，滲漏水或溶蝕帶出帷幕體中的氫氧化鈣會(huì)加劇缺陷，直至進(jìn)一步破壞防滲帷幕的連續(xù)性和完整性[1]。探索鈣離子在水工建筑物內(nèi)部的運(yùn)移規(guī)律，有助于評(píng)價(jià)防滲帷幕的耐久性。

以湖南省某抽水蓄能電站庫盆防滲工程為例，基于對(duì)流-擴(kuò)散理論，采用COMSOL多物理場(chǎng)耦合軟件，計(jì)算分析帷幕體的防滲能力，模擬研究各工況下帷幕體的溶蝕耐久性，為工程評(píng)價(jià)提供依據(jù)。

2 理論模型

2.1 滲流方程

研究壩體巖體內(nèi)部的壓力分布采用達(dá)西定律，公式如下：

(1)

(2)

2.2 溶質(zhì)運(yùn)移對(duì)流-擴(kuò)散方程

溶質(zhì)運(yùn)移的定義是地下水和土壤水中的溶質(zhì)在對(duì)流和擴(kuò)散等的共同作用下形成的物質(zhì)運(yùn)動(dòng)現(xiàn)象，故其公式使用對(duì)流-擴(kuò)散方程，方程如下：

(3)

(4)

式中，u—流體流動(dòng)速度矢量矩陣，m/s；cj—流體中j污染物的濃度，mol/L；εp—孔隙介質(zhì)孔隙率，取1；Sj—j物質(zhì)的源匯項(xiàng)。DD，j—溶質(zhì)分子擴(kuò)散系數(shù)，m2/s；De，j—分子彌散系數(shù)，m2/s。

3 數(shù)值模擬

3.1 工程簡(jiǎn)介和模型建立

湖南省某抽水蓄能電站庫壩為瀝青混凝土心墻堆石壩，壩高為48.44m，壩頂寬度為10.00m，瀝青混凝土心墻寬0.6m，正常蓄水位為43.44m。本堆石壩分為壩體和壩基，壩體由壓重區(qū)、堆石區(qū)、過渡區(qū)、瀝青混凝土心墻組成，最靠近上游的為壓重區(qū)，其次是堆石區(qū)，堆石區(qū)內(nèi)部設(shè)有2個(gè)瀝青混凝土心墻，寬度都為0.6m，下游側(cè)的瀝青混凝土心墻與堆石區(qū)之間還存在2個(gè)過渡區(qū)；壩基由全風(fēng)化、強(qiáng)風(fēng)化、弱風(fēng)化、微風(fēng)化花崗巖組成。風(fēng)化程度自上而下依次減弱。在壩基部分進(jìn)行帷幕灌漿，灌漿材料采用黏土和水泥，形成帷幕結(jié)石體厚0.8m，上部連接瀝青混凝土心墻，下部深入微風(fēng)化花崗巖層。

采用COMSOL多物理場(chǎng)耦合軟件進(jìn)行建模，使用軟件內(nèi)置的達(dá)西定律和多孔介質(zhì)中的稀物質(zhì)傳遞物理場(chǎng)。依據(jù)壩體典型斷面構(gòu)建二維模型，模型如圖1所示。網(wǎng)格采用COMSOL用戶控制網(wǎng)格進(jìn)行較細(xì)化劃分，形狀為自由三角形，最大單元大小為13.1，最小單元大小為0.0441，網(wǎng)格劃分結(jié)果如圖2所示，由圖2可知網(wǎng)格質(zhì)量良好。

圖1 壩體壩基剖面模型圖

圖2 模型網(wǎng)格較細(xì)化劃分圖

3.2 參數(shù)及邊界條件設(shè)置

本研究需要壩基各區(qū)和壩體各區(qū)的孔隙率和水力傳導(dǎo)系數(shù)，根據(jù)工程設(shè)計(jì)報(bào)告和勘測(cè)報(bào)告參數(shù)可以得到水力傳導(dǎo)系數(shù)和壩體各區(qū)的孔隙率，壩基各區(qū)的孔隙率來源于相似工程，參數(shù)見表1。

表1 壩體巖體參數(shù)明細(xì)表

對(duì)于灌漿形成的不完全連續(xù)的帷幕結(jié)石體，采用正態(tài)分布描述帷幕體中的孔隙率，稱為隨機(jī)孔隙率，選取平均值為0.23，標(biāo)準(zhǔn)差為0.01，隨機(jī)種子為1130[15-17]。模型孔隙率整體圖及帷幕體孔隙率細(xì)部圖如圖3—4所示，設(shè)置壓力水頭為43.44m。在壩基地層的左右邊界和下邊界設(shè)置無限元域，即邊界法向方向無流出。模型內(nèi)所有區(qū)域內(nèi)部濃度初始值為0mol/m3，溶質(zhì)運(yùn)移的邊界為上游帷幕面0m到-37.71m段，溶質(zhì)濃度大小為20mol/m3。

圖3 模型孔隙率整體圖

圖4 帷幕體孔隙率細(xì)部圖

4 結(jié)果分析

4.1 帷幕防滲效果分析

壩基壩體壓力水頭計(jì)算結(jié)果如圖5所示。在壩基部分，水平方向上，上游側(cè)的壓力水頭大，最大可達(dá)到38m；下游側(cè)壓力水頭小，最小為0m，壓力水頭從上游到下游呈單調(diào)遞減；豎直方向上，花崗巖壓力水頭的大小為全風(fēng)化>強(qiáng)風(fēng)化>弱風(fēng)化>微風(fēng)化，觀測(cè)壓力分布圖中的顏色分布，大部分都是平順過渡，僅帷幕結(jié)石體處不是平順過渡，表明壓力水頭在此處發(fā)生突降；在壩體部分，水平方向的壓力水頭自上游至下游單調(diào)遞減，在豎直方向較均勻分布，上游壩體正常蓄水位以上部分壓力水頭很小，下游壩體因?yàn)樾膲Ψ罎B作用，壓力水頭極小。

圖5 壩基壩體壓力圖

上述為定性分析，為了在數(shù)量上進(jìn)行分析，對(duì)上游心墻帷幕面和下游心墻帷幕面的滲流流速和壓力分布情況進(jìn)行研究，計(jì)算結(jié)果如圖6所示，圖6中橫坐標(biāo)均為弧長，圖6(a)、(c)的弧長方向?yàn)閺南轮辽?，圖6(b)、(d)的弧長方向?yàn)閺纳现料?；圖6(a)、(b)中縱坐標(biāo)代表滲流流速，圖6(c)、(d)中縱坐標(biāo)代表水壓力。

圖6 上、下游帷幕面流速和壓力變化圖

壩基地層部分，由圖6(a)可知，上游帷幕面滲流速度為0～1.0×10-4，因?yàn)閴位貙痈鲄^(qū)的水力傳導(dǎo)系數(shù)最大為5.00×10-6m/s，水力傳導(dǎo)系數(shù)過小，表明該區(qū)在有水流過時(shí)，對(duì)其有較大的阻礙作用；該段在弧長為16.8、19、37.7m處出現(xiàn)有波動(dòng)，因?yàn)榇颂帪椴煌◢弾r地層的交界處，二者水力傳導(dǎo)系數(shù)不同，導(dǎo)致流速出現(xiàn)波動(dòng)，由圖6(b)可知，在弧長48.5～86.2m處，流速呈下降趨勢(shì)，在不同地層交界處發(fā)生波動(dòng)，上下游心墻帷幕面滲流流速差距不大。由圖6(c)、(d)可知，上下游心墻帷幕面壓力差距同樣不大。

壩體部分，由圖6(a)可知，在弧長48.5～86.2m段，滲流流速先增加后減??；在弧長37.7～71.6m段，增速接近線性；把弧長為71.6m對(duì)應(yīng)的轉(zhuǎn)折點(diǎn)命名為A，A點(diǎn)滲流流速最大，為1.10×10-5m/s；由圖6(b)可知，A點(diǎn)對(duì)應(yīng)的下游帷幕面處的滲流流速為1.5×10-8m/s，為A點(diǎn)滲流流速的1‰。由圖6(c)、(d)可知，在弧長37.7～82.6m段，上游帷幕面的壓力范圍為0MPa至6.07×10-2MPa，隨著位置降低增大；在弧長37.7m處達(dá)到頂峰，觀察下游帷幕面；在弧長0～37.7m段，水壓力大小為1.5×10-3MPa；在弧長37.7～72.5m段，心墻前的水壓力大于心墻后的水壓力。

4.2 溶質(zhì)運(yùn)移分析

4.2.1溶質(zhì)擴(kuò)散范圍分析

設(shè)置計(jì)算時(shí)長為50年，步長采取1年計(jì)算1次，觀察隨時(shí)間變化的溶質(zhì)擴(kuò)散圖，因?yàn)楹笃谧兓?xì)微，前期變化快，前2年內(nèi)將步長縮小1倍，重新計(jì)算，溶質(zhì)擴(kuò)散結(jié)果如圖7所示。計(jì)算結(jié)果顯示，在壩基地層的上游帷幕面析出的溶質(zhì)能夠運(yùn)移到壩體部分，運(yùn)移主要發(fā)生在前5年，5年之后運(yùn)移程度微小。

在壩體部分，由圖7(a)、(b)可知，前半年內(nèi)，壩基地層中帷幕結(jié)石體溶蝕析出的鈣離子向著下游壩坡面的方向運(yùn)移，在豎直方向上的運(yùn)移程度比水平方向上的運(yùn)移程度大，第1年內(nèi)，溶質(zhì)已經(jīng)擴(kuò)散了下游壩體1/4的部分；由圖7(c)可知，第1年到第2年內(nèi)，溶質(zhì)在豎直方向上的運(yùn)移程度比水平方向上的運(yùn)移程度大，第2年末，溶質(zhì)擴(kuò)散了下游壩體80%的部位；由圖7(d)可知，第2年末，溶質(zhì)已擴(kuò)散到了下游坡面，查看顏色圖例，接近坡面處的濃度為14～16mol/m3；由圖7(e)可知，第30年末，溶質(zhì)已運(yùn)移到整個(gè)下游壩體，大部分濃度都在18～20mol/m3；由圖7(f)可知，第50年末，溶質(zhì)可以從壩基地層析出。

壩基地層由4層不同風(fēng)化的花崗巖層組成，每層的運(yùn)移程度都不一樣。觀察可得，微風(fēng)化、弱風(fēng)化、強(qiáng)風(fēng)化花崗巖區(qū)溶質(zhì)擴(kuò)散范圍小，全風(fēng)化花崗巖區(qū)溶質(zhì)擴(kuò)散范圍大。產(chǎn)生這種現(xiàn)象是因?yàn)榭紫堵什煌瑥?qiáng)風(fēng)化、弱風(fēng)化、微風(fēng)化花崗巖的孔隙率小，相對(duì)密實(shí)，擴(kuò)散范圍較?。蝗L(fēng)化花崗巖孔隙率大，擴(kuò)散范圍較大。溶質(zhì)在壩基地層的運(yùn)移沒有在壩體中運(yùn)移急劇。由圖7(a)可知，前半年內(nèi)，全風(fēng)化花崗巖區(qū)的溶質(zhì)運(yùn)移很快，其他花崗巖區(qū)的溶質(zhì)運(yùn)移微小，對(duì)比壩體的溶質(zhì)運(yùn)移情況發(fā)現(xiàn)，前半年內(nèi)，溶質(zhì)在全風(fēng)化花崗巖內(nèi)部運(yùn)移的程度比在下游壩體堆石區(qū)運(yùn)移的程度大；由圖7(b)可知，第1年末，與圖7(a)不同，下游壩體堆石區(qū)的運(yùn)移程度超過全風(fēng)化花崗巖區(qū)的運(yùn)移程度；由圖7(c)—(f)可知，下游側(cè)的全風(fēng)化花崗巖地層剖面圖大致為一個(gè)直角三角形，全風(fēng)化花崗巖的溶質(zhì)運(yùn)移圖同樣趨向于直角三角形，前2年內(nèi)，溶質(zhì)在微風(fēng)化，弱風(fēng)化，強(qiáng)風(fēng)化花崗巖層運(yùn)移程度微小，第2年到第5年內(nèi)，開始加速運(yùn)移，在第5年后，依然產(chǎn)生相對(duì)較大的運(yùn)移，在30年后趨向穩(wěn)定。

4.2.2溶質(zhì)擴(kuò)散通量分析

擴(kuò)散通量定義為單位時(shí)間內(nèi)垂直通過擴(kuò)散方向的單位面積的物質(zhì)的流量。擴(kuò)散通量可以表現(xiàn)溶質(zhì)在模型中任意點(diǎn)擴(kuò)散的速度，擴(kuò)散通量計(jì)算結(jié)果如圖8所示。

由圖8可知，擴(kuò)散通量圖與溶質(zhì)運(yùn)移圖中擴(kuò)散邊界完全吻合，由圖8(a)可知，溶質(zhì)即將運(yùn)移時(shí)，在帷幕結(jié)實(shí)體的頂端和底端，呈現(xiàn)深紅色，其余部分呈現(xiàn)土黃色，觀察顏色圖例表，得到在兩端擴(kuò)散通量最大，達(dá)到10-8mol/(m2·s)，其余部分為7×10-9mol/(m2·s)，表明溶質(zhì)在兩端運(yùn)移最快；由圖8(b)可知，在半年后，帷幕結(jié)石體大部分變成了青色，在壩體與壩基交界處是深紅色，在花崗巖層交界處變成綠色，觀察顏色圖例表，可知在帷幕結(jié)石體頂端溶質(zhì)運(yùn)移最快，在花崗巖層的交界處次之，在帷幕結(jié)石體的大部分區(qū)域最慢，這些區(qū)域通量下降到4×10-9mol/(m2·s)；由圖8(c)可知，擴(kuò)散通量圖在下游壩體處消失，說明此部分不存在運(yùn)移，表示此時(shí)溶質(zhì)已擴(kuò)散到下游坡面；對(duì)比圖8(c)、(d)可知，擴(kuò)散通量在十年后趨向穩(wěn)定，不再變化。

圖8 溶質(zhì)擴(kuò)散通量圖

5 結(jié)論

(1)采用正態(tài)分布方法隨機(jī)孔隙率定義防滲帷幕結(jié)石體，可以有效反映灌漿后結(jié)石體的真實(shí)情況，為同類型研究提供了新的思路。

(2)溶質(zhì)運(yùn)移速度第1、2年較快，此后開始變緩，第1年溶質(zhì)擴(kuò)散了下游壩體的1/4，第1年擴(kuò)散范圍達(dá)到80%，此后擴(kuò)散范圍緩慢增加，30年后趨于穩(wěn)定。

(3)在壩基地層部分，溶質(zhì)第1年運(yùn)移速度較快，此后開始變緩；溶質(zhì)的擴(kuò)散范圍在全風(fēng)化花崗巖地層最大，在強(qiáng)、弱風(fēng)化層運(yùn)移次之，在微風(fēng)化巖層運(yùn)移范圍最?。?0年后壩體下游有溶質(zhì)析出，導(dǎo)致防滲帷幕防滲性能降低。