碾压混凝土重力坝廊道渗流测值分析与有限元数值模拟

周永红，徐兰玉，靳向波，徐淼

（1.山西汾河二库管理有限公司，山西 太原 030012；2.水利部南京水利水文自动化研究所，江苏 南京 210012；3.水利部水文水资源监控工程技术研究中心，江苏 南京 210012；4.河海大学，江苏 南京 210098）

0 引言

碾压混凝土重力坝廊道渗水在多个工程中被发现并得到一定程度的研究[1-3]。祁立友等[4]通过对桃林口水库渗水观测资料的分析，总结出了大坝渗水变化的一般规律。袁自立等[5]综合大坝渗流观测资料和有限元模型正反分析，得出了石漫滩碾压混凝土重力坝渗流异常的成因。李荣等[6]通过ANSYS 软件，采用等效连续介质模型模拟碾压混凝土坝的层间渗流特性，探究思林水电站碾压混凝土坝的渗流规律。孙亮[7]基于碾压混凝土的渗流特性，对龙华河碾压混凝土坝进行了三维渗流场分析。方卫华等[8，9]借助COMSOL Multiphysics 软件进行三维有限元多场耦合分析，评价汾河二库拦河大坝除险加固效果。本文在已有的研究基础上，基于汾河二库廊道的渗流实测数据，在考虑帷幕和排水管等重要复杂渗控结构的基础上，充分考虑坝体实际裂缝影响，研究了渗流场、应力场和温度场三场耦合在COMSOL的实现方法。

1 工程概况

汾河二库是一座以防洪为主，兼有供水、发电、旅游、养殖等综合效益的大（2）型水利枢纽工程。拦河大坝为碾压混凝土重力坝，坝顶长227.7 m，最大坝高88.0 m，河床中部设3 个溢流表孔，每孔净宽12.0 m，堰顶高程为902.0 m。

由于历史原因，在汾河二库大坝主体工程竣工时，大坝坝体、坝肩固结和帷幕灌浆并未完成，一到冬季廊道内漏水严重。虽然2014年对坝肩、坝基及坝后连续墙等采取应急除险加固措施，但坝体内存在较长裂缝，廊道渗水现象依然存在，阻碍水库效益的发挥。

2 廊道实测渗流数据分析

2.1 渗流观测布置情况

2016年，为了解廊道内出现的异常渗流情况，在廊道内布置了23 个量水堰，用于观测不同坝段渗流情况。堰上水头测量采用振弦式4675LV-1传感器，量程300 mm，精度为0.01%FS～0.02%FS。

2.2 渗流观测数据与分析

2.2.1 渗流量过程及特征分析

坝体870 廊道的4 个测点WE-6～WE-9 在加固后渗流明显，且呈现季节规律性，测点渗流量与气温和上游水位的过程线见图1，测点总体特征值分析见表1。

由图1和表1可知，测点的渗流量与温度和上游水位具有较好的相关性，大致与温度呈负相关，与上游水位呈正相关。渗流量的最大值普遍出现在冬季末至春季初，最小值普遍出现在夏季至秋季。

表1 测点渗流量特征值分析 ml/s

图1 上游水位气温测点渗流量过程线

2.2.2 监控模型

1）裂缝渗流的确定性模型

Witherspoon PA 根据裂隙岩体渗流的立方定律，导出粗糙裂缝中稳定流流量Q的公式：

式中：G为常数，由裂缝的长度、宽度及水流粘滞系数确定；f为缝面特征系数，一般为1.04～1.65，当f=1.0 时，即为光滑情况；△H为裂缝上下游水头差；b为裂隙的等效张开度。

2）统计模型

由于大坝裂缝的几何和物理性参数均很难确定，一般不直接采用确定性模型分析大坝裂缝的渗流情况，而是根据确定性模型的形式，利用水库水位、气温和裂缝渗流量的观测数据建立裂缝渗流的统计模型。假定裂缝的下游水头为常数，按下式计算水头差△H：

式中：H为上游水位；α，H0为待定系数。

根据确定性模型的形式，基于水库2014—2019年水位、气温和典型坝段廊道裂缝渗流量的实测数据建立裂缝渗流统计模型：

式中：△TF为横缝宽度；△HF为水平缝宽度；T为观测日气温。

综合考虑式（1）—（4），按下式计算渗流量Q：

式中：c1，c2，c3为待定系数。

考虑到水温和气温的关系实际上要比线性关系复杂，增加气温T4，多项式如下：

式中：a1～a15为待定系数。

870 廊道的4 个测点WE-6～WE-9 都具有良好的相关性，将4 个测点的渗流量合计为870 廊道的总渗流量进行回归分析，拟合多项式的结果残差为154.045，拟合度R2为0.777，拟合效果较好。

3 有限元模拟及对比分析

为了从理论上论证该典型坝段渗流监测结果的准确性，进一步深入研究坝体裂缝、帷幕及坝基坝体排水管等作用下坝段渗流量的变化规律，利用有限元方法对存在问题的典型溢流坝段进行数值模拟。

3.1 模型的建立

针对汾河二库除险加固前后的结构特点和几何拓扑性质，考虑该坝段结构型式、坝体材料分区、基岩材料分区，以及帷幕灌浆、排水、固结灌浆的布置等因素，建立溢流坝段渗流场-应力场-温度场三场耦合的双重介质模型[9，10]。模型单元总数为350 717 个，计算网格如图2所示。温度场和渗流场方程采用瞬态形式，计算时长为720 d，步长为5 d。

图2 除险加固后溢流坝段三维有限元模型网格剖分示意图

渗流场计算时，上、下游水位以下为等水头边界，水位以上部分为混合边界；除险加固前，底高程为870.0 m 的廊道为可能溢出边界，底部廊道为满水有压边界；除险加固后，所有廊道表面均为可能溢出边界；其他边界为不透水边界。在COMSOL中可用透水层边界模拟混合边界条件。

位移场、应力场计算时，坝基底部视为固定约束，坝基两侧x方向为连杆约束，坝基及并缝灌浆高程（855.0 m）以下坝体左右岸方向（y方向）两侧为连杆约束。

温度场计算时，坝基底部和坝基两侧为绝热边界条件，上下游水位以下为第一类边界条件，混凝土与空气接触为第三类边界条件。基岩和混凝土都有初始温度，该地区地热温度为8.0 ℃，混凝土常温为9.5 ℃。

3.2 计算参数及边界条件

溢流坝段模型主要由坝体和地基组成，其材料性质及力学参数见表2和表3。

表2 坝体及地基材料力学参数指标

表3 坝体及地基材料热学参数指标

3.3 计算工况

基于2018—2019年2年的气温及上游水位变化，文中取不考虑溢流坝段纵横缝的情况为工况1，考虑溢流坝段纵横缝情况为工况2。

3.4 计算结果及分析

1）三场分析

坝体部分的第一主应力分布、最低温度分布、最高水位的渗流场水头分别见图3～5。在考虑裂缝之后，坝体增加了裂缝渗流引起的压力，使坝体的第一主应力增加；渗流量的增加引起了坝体内水头的抬高；应力的抬高及渗流速度的增加，加快了坝体内的热交换及传导速率，从而使温度有了略微的提升。

图3 坝体工况1 和2 的第一主应力图（单位：MPa）

图4 坝体工况1 和2 的温度等值线图（单位：℃）

图5 坝体工况1 和2 的水头等值线图（单位：m）

2）渗流量分析

选取870 廊道工况1，2 条件下的渗流量计算结果与上游水位、气温进行对比，如图6～9所示。

由图6和图7所示，在水位和气温的双重影响下，有裂缝的渗流量大，工况2 相对于工况1 的变幅较大；由图8和图9可知，温度对于裂缝渗流的影响更大，裂缝引起的渗流量与温度有着明显的负相关性。

图6 渗流量计算值与水位对比

图7 渗流量与气温对比

图8 裂缝引起的渗流量与水位对比

图9 裂缝引起的渗流量与气温对比

3）计算值与实测值对比分析

采用双重介质渗流-应力-温度三场耦合模型计算工况2 条件下870 廊道的渗流量，并与实测渗流量、上游水位及气温进行比较，如图10～11所示。由图10可知，采用耦合模型计算得到的渗流量与实测渗流量拟合趋势性较好，但在渗流量最低时吻合性稍低，主要原因是有限元计算采用的温度值为大气温度的拟合值，造成模型计算过程中温度值的偏差。同时，由图11可知，气温在20 ℃以上时，温度对实测渗流量影响很小，基本可以忽略，但这在有限元计算结果中不能较好地体现，从而引起了温度变化带来的渗流量偏大。分析还表明，在水位和气温的共同影响下，坝体渗流和裂缝渗流共同作用引起了廊道渗流量的季节性变化。

图10 渗流量实测值与计算值和水位对比

图11 渗流量实测值与计算值和气温对比

4 结语

本文根据汾河二库实际环境，以求解渗流量为目的，建立了水库溢流坝段双重介质三场耦合模型，并用COMSOL Multiphyics 模拟有无裂缝的各工作性态及廊道的渗流情况。该法有效融合了动态监测数值和有限元计算结果，建立了坝体廊道内渗流量的在线监测与计算模型，可模拟长周期内坝体渗流变化情况，为大坝渗流安全性态评价提供了新的思路和方法。该法模拟结果与实测值吻合度较好，但在渗流量最低时吻合性稍低，后期尚需进一步研究分析。