基于对流-扩散理论的防渗帷幕溶质溶蚀运移分析
2022-05-07蒋明东李永丰张贵金
蒋明东,陶 骏,李永丰,任 鑫,李 骥,张贵金
(1.湖南平江抽水蓄能有限公司,湖南 平江 410400;2.长沙理工大学,湖南 长沙 410014;3.中国电建集团中南勘测设计研究院有限公司,湖南 长沙 410014)
1 概论
帷幕灌浆是坝基防渗的常用方法,但由于施工难度,形成的帷幕体或多或少存在蜂窝、孔洞、孔隙等缺陷或薄弱段,在水头作用或长期运行后很容易产生渗漏通道,渗漏水或溶蚀带出帷幕体中的氢氧化钙会加剧缺陷,直至进一步破坏防渗帷幕的连续性和完整性[1]。探索钙离子在水工建筑物内部的运移规律,有助于评价防渗帷幕的耐久性。
以湖南省某抽水蓄能电站库盆防渗工程为例,基于对流-扩散理论,采用COMSOL多物理场耦合软件,计算分析帷幕体的防渗能力,模拟研究各工况下帷幕体的溶蚀耐久性,为工程评价提供依据。
2 理论模型
2.1 渗流方程
研究坝体岩体内部的压力分布采用达西定律,公式如下:
(1)
(2)
2.2 溶质运移对流-扩散方程
溶质运移的定义是地下水和土壤水中的溶质在对流和扩散等的共同作用下形成的物质运动现象,故其公式使用对流-扩散方程,方程如下:
(3)
(4)
式中,u—流体流动速度矢量矩阵,m/s;cj—流体中j污染物的浓度,mol/L;εp—孔隙介质孔隙率,取1;Sj—j物质的源汇项。DD,j—溶质分子扩散系数,m2/s;De,j—分子弥散系数,m2/s。
3 数值模拟
3.1 工程简介和模型建立
湖南省某抽水蓄能电站库坝为沥青混凝土心墙堆石坝,坝高为48.44m,坝顶宽度为10.00m,沥青混凝土心墙宽0.6m,正常蓄水位为43.44m。本堆石坝分为坝体和坝基,坝体由压重区、堆石区、过渡区、沥青混凝土心墙组成,最靠近上游的为压重区,其次是堆石区,堆石区内部设有2个沥青混凝土心墙,宽度都为0.6m,下游侧的沥青混凝土心墙与堆石区之间还存在2个过渡区;坝基由全风化、强风化、弱风化、微风化花岗岩组成。风化程度自上而下依次减弱。在坝基部分进行帷幕灌浆,灌浆材料采用黏土和水泥,形成帷幕结石体厚0.8m,上部连接沥青混凝土心墙,下部深入微风化花岗岩层。
采用COMSOL多物理场耦合软件进行建模,使用软件内置的达西定律和多孔介质中的稀物质传递物理场。依据坝体典型断面构建二维模型,模型如图1所示。网格采用COMSOL用户控制网格进行较细化划分,形状为自由三角形,最大单元大小为13.1,最小单元大小为0.0441,网格划分结果如图2所示,由图2可知网格质量良好。
图1 坝体坝基剖面模型图
图2 模型网格较细化划分图
3.2 参数及边界条件设置
本研究需要坝基各区和坝体各区的孔隙率和水力传导系数,根据工程设计报告和勘测报告参数可以得到水力传导系数和坝体各区的孔隙率,坝基各区的孔隙率来源于相似工程,参数见表1。
表1 坝体岩体参数明细表
对于灌浆形成的不完全连续的帷幕结石体,采用正态分布描述帷幕体中的孔隙率,称为随机孔隙率,选取平均值为0.23,标准差为0.01,随机种子为1130[15-17]。模型孔隙率整体图及帷幕体孔隙率细部图如图3—4所示,设置压力水头为43.44m。在坝基地层的左右边界和下边界设置无限元域,即边界法向方向无流出。模型内所有区域内部浓度初始值为0mol/m3,溶质运移的边界为上游帷幕面0m到-37.71m段,溶质浓度大小为20mol/m3。
图3 模型孔隙率整体图
图4 帷幕体孔隙率细部图
4 结果分析
4.1 帷幕防渗效果分析
坝基坝体压力水头计算结果如图5所示。在坝基部分,水平方向上,上游侧的压力水头大,最大可达到38m;下游侧压力水头小,最小为0m,压力水头从上游到下游呈单调递减;竖直方向上,花岗岩压力水头的大小为全风化>强风化>弱风化>微风化,观测压力分布图中的颜色分布,大部分都是平顺过渡,仅帷幕结石体处不是平顺过渡,表明压力水头在此处发生突降;在坝体部分,水平方向的压力水头自上游至下游单调递减,在竖直方向较均匀分布,上游坝体正常蓄水位以上部分压力水头很小,下游坝体因为心墙防渗作用,压力水头极小。
图5 坝基坝体压力图
上述为定性分析,为了在数量上进行分析,对上游心墙帷幕面和下游心墙帷幕面的渗流流速和压力分布情况进行研究,计算结果如图6所示,图6中横坐标均为弧长,图6(a)、(c)的弧长方向为从下至上,图6(b)、(d)的弧长方向为从上至下;图6(a)、(b)中纵坐标代表渗流流速,图6(c)、(d)中纵坐标代表水压力。
图6 上、下游帷幕面流速和压力变化图
坝基地层部分,由图6(a)可知,上游帷幕面渗流速度为0~1.0×10-4,因为坝基地层各区的水力传导系数最大为5.00×10-6m/s,水力传导系数过小,表明该区在有水流过时,对其有较大的阻碍作用;该段在弧长为16.8、19、37.7m处出现有波动,因为此处为不同花岗岩地层的交界处,二者水力传导系数不同,导致流速出现波动,由图6(b)可知,在弧长48.5~86.2m处,流速呈下降趋势,在不同地层交界处发生波动,上下游心墙帷幕面渗流流速差距不大。由图6(c)、(d)可知,上下游心墙帷幕面压力差距同样不大。
坝体部分,由图6(a)可知,在弧长48.5~86.2m段,渗流流速先增加后减小;在弧长37.7~71.6m段,增速接近线性;把弧长为71.6m对应的转折点命名为A,A点渗流流速最大,为1.10×10-5m/s;由图6(b)可知,A点对应的下游帷幕面处的渗流流速为1.5×10-8m/s,为A点渗流流速的1‰。由图6(c)、(d)可知,在弧长37.7~82.6m段,上游帷幕面的压力范围为0MPa至6.07×10-2MPa,随着位置降低增大;在弧长37.7m处达到顶峰,观察下游帷幕面;在弧长0~37.7m段,水压力大小为1.5×10-3MPa;在弧长37.7~72.5m段,心墙前的水压力大于心墙后的水压力。
4.2 溶质运移分析
4.2.1溶质扩散范围分析
设置计算时长为50年,步长采取1年计算1次,观察随时间变化的溶质扩散图,因为后期变化细微,前期变化快,前2年内将步长缩小1倍,重新计算,溶质扩散结果如图7所示。计算结果显示,在坝基地层的上游帷幕面析出的溶质能够运移到坝体部分,运移主要发生在前5年,5年之后运移程度微小。
在坝体部分,由图7(a)、(b)可知,前半年内,坝基地层中帷幕结石体溶蚀析出的钙离子向着下游坝坡面的方向运移,在竖直方向上的运移程度比水平方向上的运移程度大,第1年内,溶质已经扩散了下游坝体1/4的部分;由图7(c)可知,第1年到第2年内,溶质在竖直方向上的运移程度比水平方向上的运移程度大,第2年末,溶质扩散了下游坝体80%的部位;由图7(d)可知,第2年末,溶质已扩散到了下游坡面,查看颜色图例,接近坡面处的浓度为14~16mol/m3;由图7(e)可知,第30年末,溶质已运移到整个下游坝体,大部分浓度都在18~20mol/m3;由图7(f)可知,第50年末,溶质可以从坝基地层析出。
坝基地层由4层不同风化的花岗岩层组成,每层的运移程度都不一样。观察可得,微风化、弱风化、强风化花岗岩区溶质扩散范围小,全风化花岗岩区溶质扩散范围大。产生这种现象是因为孔隙率不同,强风化、弱风化、微风化花岗岩的孔隙率小,相对密实,扩散范围较小;全风化花岗岩孔隙率大,扩散范围较大。溶质在坝基地层的运移没有在坝体中运移急剧。由图7(a)可知,前半年内,全风化花岗岩区的溶质运移很快,其他花岗岩区的溶质运移微小,对比坝体的溶质运移情况发现,前半年内,溶质在全风化花岗岩内部运移的程度比在下游坝体堆石区运移的程度大;由图7(b)可知,第1年末,与图7(a)不同,下游坝体堆石区的运移程度超过全风化花岗岩区的运移程度;由图7(c)—(f)可知,下游侧的全风化花岗岩地层剖面图大致为一个直角三角形,全风化花岗岩的溶质运移图同样趋向于直角三角形,前2年内,溶质在微风化,弱风化,强风化花岗岩层运移程度微小,第2年到第5年内,开始加速运移,在第5年后,依然产生相对较大的运移,在30年后趋向稳定。
4.2.2溶质扩散通量分析
扩散通量定义为单位时间内垂直通过扩散方向的单位面积的物质的流量。扩散通量可以表现溶质在模型中任意点扩散的速度,扩散通量计算结果如图8所示。
由图8可知,扩散通量图与溶质运移图中扩散边界完全吻合,由图8(a)可知,溶质即将运移时,在帷幕结实体的顶端和底端,呈现深红色,其余部分呈现土黄色,观察颜色图例表,得到在两端扩散通量最大,达到10-8mol/(m2·s),其余部分为7×10-9mol/(m2·s),表明溶质在两端运移最快;由图8(b)可知,在半年后,帷幕结石体大部分变成了青色,在坝体与坝基交界处是深红色,在花岗岩层交界处变成绿色,观察颜色图例表,可知在帷幕结石体顶端溶质运移最快,在花岗岩层的交界处次之,在帷幕结石体的大部分区域最慢,这些区域通量下降到4×10-9mol/(m2·s);由图8(c)可知,扩散通量图在下游坝体处消失,说明此部分不存在运移,表示此时溶质已扩散到下游坡面;对比图8(c)、(d)可知,扩散通量在十年后趋向稳定,不再变化。
图8 溶质扩散通量图
5 结论
(1)采用正态分布方法随机孔隙率定义防渗帷幕结石体,可以有效反映灌浆后结石体的真实情况,为同类型研究提供了新的思路。
(2)溶质运移速度第1、2年较快,此后开始变缓,第1年溶质扩散了下游坝体的1/4,第1年扩散范围达到80%,此后扩散范围缓慢增加,30年后趋于稳定。
(3)在坝基地层部分,溶质第1年运移速度较快,此后开始变缓;溶质的扩散范围在全风化花岗岩地层最大,在强、弱风化层运移次之,在微风化岩层运移范围最小;50年后坝体下游有溶质析出,导致防渗帷幕防渗性能降低。