CN106295066A - 坝基帷幕防渗性能衰减的数值模拟方法 - Google Patents

Abstract

本发明涉及一种坝基帷幕防渗性能衰减的数值模拟方法，该方法以易发生溶解且作为帷幕主要组成的Ca(OH)₂作为目标组分，首先根据待模拟坝基的实际物理参数设置边界条件，然后将待模拟时段划分为多个子时段，计算每个子时段的渗流速度、孔隙度、目标组分浓度、渗透系数、弥散系数以及溶解速率，从而研究其溶解引起的帷幕防渗性能变化，以及由此引起的幕后渗流、离子浓度等的时空演变。本发明可以较好地反映帷幕孔隙度等细观结构的变化，并可以预测由帷幕引起的坝基渗流、离子浓度等的变化。

Description

坝基帷幕防渗性能衰减的数值模拟方法

技术领域

本发明涉及一种坝基帷幕防渗性能衰减的数值模拟方法，属于大坝安全监测技术领域。

背景技术

帷幕作为坝基防渗体系的重要组成部分，对大坝安全运行及防渗起着极其重要的作用。水利工程中有众多事故是由渗流所引发的，统计结果表明我国1981年241座大型水库的100宗事故中，由渗流问题而导致水库失事的占到30％～40％左右。因此，坝基帷幕的防渗问题得到越来越广泛的重视。

目前对帷幕防渗的研究主要集中在以下两个方面：(1)帷幕不同性态参数条件下对坝基渗流的影响：如通过数值模拟手段研究坝基帷幕不同深度、厚度、位置等因素对防渗效果的影响，以及通过建立帷幕随机缺损的数学模型模拟防渗体缺损大小、分布随机性对渗流场的影响；(2)对帷幕防渗效果及其耐久性进行评价：如通过单一指标，如水质等微观效应量，基于水化学图示法、饱和指数模型等方法对不同时期的帷幕防渗性能进行评价，以及通过建立幕后地下水位、排水量、扬压力等宏观指标和坝基地下水温度、pH值及其他水化学特征等微观指标相结合的多指标评价体系进行评价。但申请人检索发现，目前对水库运行过程中帷幕防渗性能衰减的研究还较为少见，尤其是考虑帷幕中水泥结石等组分的溶失对帷幕细观结构的改变等还未涉及。

发明内容

本发明要解决技术问题是：提供一种模拟大坝基础帷幕防渗性能衰减的方法，该方法可以为大坝的安全运行和维护提供依据从而减少因渗流引起的水坝失事事故。

为了解决上述技术问题，本发明提出的技术方案是：一种坝基帷幕防渗性能衰减的数值模拟方法，包括以下步骤：

1)根据待模拟的坝基的实际物理参数设置边界条件，并得到模拟前所述坝基的渗流速度U₀、孔隙度θ₀、目标组分浓度C₀、渗透系数K₀、弥散系数D₀以及溶解速率R₀；所述坝基包括基岩和帷幕两部分,所述目标组分浓度包括钙离子浓度和氢氧根离子浓度。如何设置边界条件为现有技术，本领域技术人员在看到坝基的实际物理参数后可以很容易地得到边界条件，不再赘述。

2)根据步骤1)设置的边界条件，建立渗流方程、溶质迁移方程和固相介质时变方程；

所述渗流方程为：

$\frac{\∂}{\∂x}\left(K\frac{\∂H}{\∂x}\right)+\frac{\∂}{\∂y}\left(K\frac{\∂H}{\∂y}\right)=\mathrm{\μ}\frac{\∂H}{\∂t}$

式中，H为水头，K为渗透系数，t为时间参数，μ为给水度，x、y分别为坝基的横坐标、纵坐标；

所述溶质迁移方程为：

$\frac{\∂\left(\mathrm{\θ}C\right)}{\∂t}=\frac{\∂}{\∂x}\[D_L\frac{\∂\left(\mathrm{\θ}C\right)}{\∂x}\]+\frac{\∂}{\∂y}\[D_T\frac{\∂\left(\mathrm{\θ}C\right)}{\∂y}\]-\frac{\∂\left(u_x\mathrm{\θ}C\right)}{\∂x}-\frac{\∂\left(u_y\mathrm{\θ}C\right)}{\∂y}+\mathrm{\θ}\·R$

式中，R为可溶组分的溶解速率，C为目标组分浓度；渗流速度U分解为水平方向的渗流速度u_x和垂直方向的渗流速度u_y；θ为孔隙度；弥散系数D分解为纵向弥散系数D_T和横向弥散系数D_L；x、y分别为坝基的横坐标、纵坐标；t为时间参数；

所述固相介质时变方程为：

$\frac{\∂\mathrm{\θ}}{\∂t}=R\·M\·\mathrm{\θ}$

式中，M为Ca(OH)₂的摩尔体积，R为可溶组分的溶解速率，θ为孔隙度，t为时间参数。

3)将待模拟时段划分为n个子时段，其中第i个子时段的初始时刻为终了时刻为1≤i≤n，n≥2；

根据步骤1)得到的所述坝基帷幕的渗流速度U₀、孔隙度θ₀、目标组分浓度C₀、渗透系数K₀、弥散系数D₀以及溶解速率R₀，对步骤2)建立的渗流方程、溶质迁移方程和固相介质时变方程进行求解，依次计算得到第j个子时段的终了时刻的渗流速度、孔隙度以及目标组分浓度，1≤j≤n；同时利用渗透系数K与孔隙度θ之间的函数、弥散系数D与渗流速度U、孔隙度θ之间的函数以及可溶组分的溶解速率R与目标组分浓度C之间的函数，得到第j个子时段的终了时刻的渗透系数、弥散系数以及溶解速率；并将第j个子时段的终了时刻的渗流速度、孔隙度、目标组分浓度、渗透系数、弥散系数以及溶解速率作为第j+1个子时段的初始时刻的相应值；

其中，所述渗透系数K与孔隙度θ之间的函数为：

K＝K₀·(θ/θ₀)³[(1-θ₀)/(1-θ)]²

式中，K₀和θ₀分别为模拟前所述坝基帷幕的渗透系数和孔隙度；

所述弥散系数D与渗流速度U、孔隙度θ之间的函数为：

D＝α·|U|+D_m·θ

式中，α为弥散度，D_m为钙离子在水中的扩散系数；

所述可溶组分的溶解速率R与目标组分浓度C之间的函数为：

$R=A\·{(1-\frac{C_{{\mathrm{Ca}}^{2+}}\·C_{{\mathrm{OH}}^{-}}^2}{K_{sp}})}^n$

式中，A为反应动力学系数，n为动力学指数，K_sp为Ca(OH)₂的溶度积常数，为钙离子浓度，为氢氧根离子浓度；A，n，K_sp均为预设系数。

模拟前所述坝基的孔隙度θ₀可由坝基的几何参数直接得到。步骤1)中以上游库水位作为初始水头，在所述边界条件的约束下对坝基进行渗流和溶质运移的稳定态模拟，获得渗流稳定态时的渗流场和化学场分布情况，从而得到模拟前所述坝基的渗流速度U₀、目标组分浓度C₀、渗透系数K₀、弥散系数D₀以及溶解速率R₀。在所述坝基的孔隙度θ₀已知的前提下，对坝基进行渗流和溶质运移稳定态模拟为现有技术，不再赘述。

本发明的原理及带来的有益效果如下：

大坝在长期运行条件下，帷幕受到坝基水的侵蚀作用，其矿物组成的溶失使得帷幕细观结构发生改变，引起防渗性能衰减，可能造成渗流安全隐患。为反映帷幕防渗性能衰减机理及其造成的影响，本发明依据地下水动力学、溶质运移以及化学动力学的相关理论，建立了多物理场耦合的数学模型，以易发生溶解且作为帷幕主要组成的Ca(OH)2作为目标组分，研究其溶解引起的帷幕防渗性能变化，以及由此引起的幕后渗流、离子浓度等的时空演变。本发明的模拟结果表明，模型能较好地反映帷幕孔隙度等细观结构的变化，溶蚀部位主要发生在帷幕顶部和底部，其中顶部以贯穿型侵蚀为主，而底部则主要为表层侵蚀。同时，模型还可预测由帷幕引起的坝基渗流、离子浓度等的变化。

通过使用本发明的方法模拟大坝基础帷幕在某一时间段内的防渗性能衰减值，包括帷幕的渗流速度和孔隙度等，这样可以反映大坝基础帷幕防渗性能在该时段内的衰减过程及衰减程度，从而可以刻画该时段内坝基帷幕细观结构和防渗性能等随时间演变及空间分布情况。

本发明提出的坝基帷幕防渗性能衰减的数值模拟方法具有良好的精度，为直观、方便地了解坝基帷幕防渗性能的衰减发生的部位及程度提供可能，解决了因溶失引起帷幕防渗性能衰减的量化解析问题。本发明通过模拟帷幕中可溶组分溶解引起的帷幕孔隙度增大以及由此引起的防渗性能变化，可以为大坝的安全运行和维护提供依据从而减少因渗流引起的水坝失事事故。

附图说明

下面结合附图对本发明作进一步说明。

图1是本发明实施例中待模拟坝基的结构示意图。

图2为模拟结束后帷幕的孔隙度二维空间分布示意图。

图3为模拟结束后帷幕典型断面孔隙度的一维空间分布示意图。

图4为模拟时段内帷幕防渗性能衰减程度的断面示意图。

具体实施方式

实施例

本实施例的待模拟的坝基二维模型如图1所示，大坝高150m，坝底宽60m，帷幕深100m，帷幕宽5m。其中，由ABCDEFIHG组成的区域Ω1为坝基基岩区域，由GHIF组成的区域Ω2为坝基帷幕。以大坝上游侧底部与坝基相连接处作为坐标原点(x＝0m,y＝0m)，指向下游侧为x轴正向，指向坝顶方向为y轴正向，建立坐标系。

本实施例中坝基矿物组成可视为惰性，并仅考虑帷幕组成中Ca(OH)₂的溶解作用。

本实施例的坝基帷幕防渗性能衰减的数值模拟方法，包括以下步骤：

1)根据待模拟的坝基的实际物理参数设置边界条件，并得到模拟前所述坝基的渗流速度U₀、孔隙度θ₀、目标组分浓度C₀、渗透系数K₀、弥散系数D₀以及溶解速率R₀，边界条件的设定为现有技术，不再赘述。所述坝基包括基岩和帷幕两部分，所述目标组分浓度包括钙离子浓度和氢氧根离子浓度。

模拟前所述坝基的孔隙度θ₀可由坝基的几何参数直接得到。在实施时，限定区域Ω1和Ω2中的初始水头为上游库水位，在边界条件的约束下对坝基进行渗流和溶质运移的稳定态模拟，获得渗流稳定态时的渗流场和化学场分布情况，从而得到模拟前所述坝基的渗流速度U₀、目标组分浓度C₀、渗透系数K₀、弥散系数D₀以及溶解速率R₀，即第1个子时段的初始时刻的渗流速度、孔隙度、目标组分浓度、渗透系数、弥散系数以及溶解速率。

在所述坝基的孔隙度θ₀已知的前提下，对坝基进行渗流稳定态模拟为现有技术，不再赘述。

2)根据步骤1)设置的边界条件，建立渗流方程、溶质迁移方程和固相介质时变方程；

所述渗流方程为：

$\frac{\∂}{\∂x}\left(K\frac{\∂H}{\∂x}\right)+\frac{\∂}{\∂y}\left(K\frac{\∂H}{\∂y}\right)=\mathrm{\μ}\frac{\∂H}{\∂t}$

式中，H为水头，K为渗透系数，t为时间参数，μ为给水度，x、y分别为坝基的横坐标、纵坐标。

渗流方程中边界条件设定为：边界AG和边界DE为第一类边界条件，即给定水头边界，分别等于上游库水位和下游库水位；边界AB、BC、CD、EF、FG作为第二类边界条件中的隔水边界。

所述溶质迁移方程为：

$\frac{\∂\left(\mathrm{\θ}C\right)}{\∂t}=\frac{\∂}{\∂x}\[D_L\frac{\∂\left(\mathrm{\θ}C\right)}{\∂x}\]+\frac{\∂}{\∂y}\[D_T\frac{\∂\left(\mathrm{\θ}C\right)}{\∂y}\]-\frac{\∂\left(u_x\mathrm{\θ}C\right)}{\∂x}-\frac{\∂\left(u_y\mathrm{\θ}C\right)}{\∂y}+\mathrm{\θ}\·R$

式中，R为可溶组分的溶解速率，C为目标组分浓度；渗流速度U分解为水平方向的渗流速度u_x和垂直方向的渗流速度u_y；θ为孔隙度；弥散系数D中涉及到的纵向弥散度α_L＝40m和横向弥散度α_T＝5m；x、y分别为坝基的横坐标、纵坐标；t为时间参数。

溶质迁移方程中边界条件设定为：边界AG为第一类边界条件，即给定浓度边界，等于上游库水中Ca²⁺浓度；边界DE为第二类边界，AB、BC、CD、EF、FG作为零通量边界。

所述固相介质时变方程为：

$\frac{\∂\mathrm{\θ}}{\∂t}=R\·M\·\mathrm{\θ}$

式中，M为Ca(OH)₂的摩尔体积，M＝33cm³mol^-1，R为可溶组分的溶解速率，θ为孔隙度，t为时间参数；

固相介质时变方程中边界条件设定为：边界GH、HI、IF、FG为第二类边界条件。

3)将待模拟时段划分为n个子时段，其中第i个子时段的初始时刻为终了时刻为1≤i≤n，n≥2；本实施例模拟大坝在运行一百年期间帷幕防渗性能的衰减过程，令时间间隔为10天，则总的模拟时段可划分为3650个子时段，即n＝3650。

根据步骤1)得到的所述坝基帷幕的渗流速度U₀、孔隙度θ₀、目标组分浓度C₀、渗透系数K₀、弥散系数D₀以及溶解速率R₀，对步骤2)建立的渗流方程、溶质迁移方程和固相介质时变方程进行求解，依次计算得到第j个子时段的终了时刻的渗流速度、孔隙度以及目标组分浓度，1≤j≤n；同时利用渗透系数K与孔隙度θ之间的函数、弥散系数D与渗流速度U、孔隙度θ之间的函数以及可溶组分的溶解速率R与目标组分浓度C之间的函数，得到第j个子时段的终了时刻的渗透系数、弥散系数以及溶解速率；并将第j个子时段的终了时刻的渗流速度、孔隙度、目标组分浓度、渗透系数、弥散系数以及溶解速率作为第j+1个子时段的初始时刻的相应值；

其中，所述渗透系数K与孔隙度θ之间的函数为：

K＝K₀·(θ/θ₀)³[(1-θ₀)/(1-θ)]²

式中，K₀和θ₀分别为模拟前所述坝基帷幕的渗透率和孔隙度；

所述弥散系数D与渗流速度U、孔隙度θ之间的函数为：

D＝α·|U|+D_m·θ

式中，α为弥散度，D_m为钙离子在水中的扩散系数，D_m＝1e^-9m²s^-1；

所述可溶组分的溶解速率R与目标组分浓度C之间的函数为：

$R=A\·{(1-\frac{C_{{\mathrm{Ca}}^{2+}}\·C_{{\mathrm{OH}}^{-}}^2}{K_{sp}})}^n$

式中，A为反应动力学系数，n为动力学指数，K_sp为Ca(OH)₂的溶度积常数，为钙离子浓度，为氢氧根离子浓度；A，n，K_sp均为预设系数。本实施例中n取4.5，A取1×10^{- 8}mol·L^-1·s^-1，K_sp取1×10^-5.33。

模拟结束后，参照图2可以看出，帷幕溶解主要发生在两个以下两个部位：1)帷幕浅层，即帷幕上部靠近坝体部位，帷幕溶失现象较为明显，自上游侧至下游侧孔隙度皆相对较大，形成一个较为明显的“贯穿型”渗流通道，且上游侧溶失面积较大；2)帷幕底部溶失也较为明显，但较顶部相对较弱，未形成自上游至下游的贯穿型渗流通道，仅沿帷幕表面向内部发生浅层溶蚀。从图2可以看出，模拟的结果与帷幕溶解的实测结果非常接近。

另外，分别选取帷幕内x＝1m、x＝2.5m和x＝4m处作为帷幕靠近上游侧、中心线和靠近下游侧的典型断面，作模拟时段结束时此三条断面处帷幕孔隙度的分布，如图3所示可以看出：1)帷幕顶部和底部处孔隙度相对较大，而帷幕中部孔隙度相对较小；2)靠近帷幕顶部部位处，无论是上游侧、中心线还是下游侧，孔隙度皆接近0.2，表明此部位的水泥结石已全部溶失，且自上游至下游侧溶失面积逐渐减小；3)帷幕中部，孔隙度自上游侧值下游侧逐渐减小，表明帷幕中部所受溶蚀自上游至下游逐渐减小；4)帷幕底部部位，上、下游侧的孔隙度明显大于帷幕中心线位置。

本实施例还可以计算获得帷幕防渗性能的衰减程度，选取帷幕中心线x＝2.5m处，作y＝0～-250m的典型断面，通过该断面的水流量如图4所示，通过坝基典型断面处的流量随时间呈明显增大的趋势，初始时刻通过断面流量为3.261356×10^-4m²/s，而在模拟时段末期流量为3.818986×10^-4m²/s，较初始时刻增大了约17.1％。

本发明不局限于上述实施例所述的具体技术方案，除上述实施例外，本发明还可以有其他实施方式。凡采用等同替换形成的技术方案，均为本发明要求的保护范围。

Claims (2)

1.一种坝基帷幕防渗性能衰减的数值模拟方法，包括以下步骤：

1)根据待模拟的坝基的实际物理参数设置边界条件，并得到模拟前所述坝基的渗流速度U₀、孔隙度θ₀、目标组分浓度C₀、渗透系数K₀、弥散系数D₀以及溶解速率R₀；所述坝基包括基岩和帷幕两部分，目标组分浓度包括钙离子浓度和氢氧根离子浓度；

2)根据步骤1)设置的边界条件，建立渗流方程、溶质迁移方程和固相介质时变方程；

所述渗流方程为：

$\frac{\∂}{\∂x}\left(K\frac{\∂H}{\∂x}\right)+\frac{\∂}{\∂y}\left(K\frac{\∂H}{\∂y}\right)=\mathrm{\μ}\frac{\∂H}{\∂t}$

式中，H为水头，K为渗透系数，t为时间参数，μ为给水度，x、y分别为坝基的横坐标、纵坐标；

所述溶质迁移方程为：

$\frac{\∂\left(\mathrm{\θ}C\right)}{\∂t}=\frac{\∂}{\∂x}\[D_L\frac{\∂\left(\mathrm{\θ}C\right)}{\∂x}\]+\frac{\∂}{\∂y}\[D_T\frac{\∂\left(\mathrm{\θ}C\right)}{\∂y}\]-\frac{\∂\left(u_x\mathrm{\θ}C\right)}{\∂x}-\frac{\∂\left(u_y\mathrm{\θ}C\right)}{\∂y}+\mathrm{\θ}\·R$

式中，R为可溶组分的溶解速率，C为目标组分浓度；U为渗流速度，可分解为水平方向的渗流速度u_x和垂直方向的渗流速度u_y；θ为孔隙度；D为弥散系数，可分解为纵向弥散系数D_T和横向弥散系数D_L；x、y分别为坝基的横坐标、纵坐标；t为时间参数；

所述固相介质时变方程为：

$\frac{\∂\mathrm{\θ}}{\∂t}=R\·M\·\mathrm{\θ}$

式中，M为Ca(OH)₂的摩尔体积，R为可溶组分的溶解速率，θ为孔隙度，t为时间参数；

3)将待模拟时段划分为n个子时段，其中第i个子时段的初始时刻为终了时刻为1≤i≤n，n≥2；

根据步骤1)得到的所述模拟前坝基帷幕的渗流速度U₀、孔隙度θ₀、目标组分浓度C₀、渗透系数K₀、弥散系数D₀以及溶解速率R₀，对步骤2)建立的渗流方程、溶质迁移方程和固相介质时变方程进行求解，依次计算得到第j个子时段的终了时刻的渗流速度、孔隙度以及目标组分浓度，j∈{1,2,…,n}；同时利用渗透系数K与孔隙度θ之间的函数、弥散系数D与渗流速度U、孔隙度θ之间的函数以及可溶组分的溶解速率R与目标组分浓度C之间的函数，得到第j个子时段的终了时刻的渗透系数、弥散系数以及溶解速率；并将第j个子时段的终了时刻的渗流速度、孔隙度、目标组分浓度、渗透系数、弥散系数以及溶解速率作为第j+1个子时段的初始时刻的相应值；

其中，所述渗透系数K与孔隙度θ之间的函数为：

K＝K₀·(θ/θ₀)³[(1-θ₀)/(1-θ)]²

式中，K₀和θ₀分别为模拟前所述坝基的渗透系数和孔隙度；

所述弥散系数D与渗流速度U、孔隙度θ之间的函数为：

D＝α·|U|+D_m·θ

式中，α为弥散度，D_m为钙离子在水中的扩散系数；

所述可溶组分的溶解速率R与目标组分浓度C之间的函数为：

$R=A\·{(1-\frac{C_{{\mathrm{Ca}}^{2+}}\·C_{{\mathrm{OH}}^{-}}^2}{K_{sp}})}^n$

式中，A为反应动力学系数，n为动力学指数，K_sp为Ca(OH)₂的溶度积常数，为钙离子浓度，为氢氧根离子浓度；A，n，K_sp均为预设系数。

2.根据权利要求1所述的坝基帷幕防渗性能衰减的数值模拟方法，其特征在于：步骤1)中以上游库水位作为初始水头，在所述边界条件的约束下对坝基进行渗流和溶质运移的稳定态模拟，获得渗流稳定态时的渗流场和化学场分布情况，从而得到模拟前所述坝基的渗流速度U₀、目标组分浓度C₀、渗透系数K₀、弥散系数D₀以及溶解速率R₀。