陈彬鑫，聂兴山，田 淳
(1.太原理工大学 水利科学与工程学院，山西 太原 030024；
2.山西省水土保持科学研究所，山西 太原 030013)
为了提高水资源利用率，山西省水利厅将“淤地坝蓄水利用”作为水利扶贫的一项重要举措，明确要求加强小流域地表水开发，重点开展淤地坝防渗改造技术的研究[1]。为保证大坝在蓄水时能够安全运行，需对大坝改造前后的渗流和坝坡稳定性进行分析。大坝的渗流分析和坝坡稳定分析一般是分开进行的，即在进行渗流分析时不考虑大坝应力场的影响，在完成渗流分析后，将渗流分析得到的结果作为坝坡稳定计算的初始条件[2]。实际上，大坝蓄水运行时的渗流场和应力场一直是动态和相互影响的。在大坝上下游存在一定水头差的情况下，坝体内部会形成稳定的渗流。渗透水流产生的动水压力和静水压力，使得土体的有效应力分布发生变化，同时，有效应力反过来影响多孔骨架土体的孔隙比，从而改变土体的渗透系数[3]。李宗坤等[4]根据渗流理论和弹性模型建立流固耦合模型，结果表明考虑耦合的计算结果相比于传统方法的计算结果更加保守；
王嘉贵等[5]通过实测数据与计算结果的对比，发现流固耦合计算结果与实际情况更加接近。因此，为了真实反映大坝的运行状态，提高大坝安全性评价的科学性，在进行渗流和坝坡稳定计算时，应该考虑渗流场和应力场的耦合影响[6]。
咀儿上淤地坝下游有耕地和截潜流水源，不宜进行施工。本文在保持淤地坝下游坝坡不变的情况下，针对大坝上游提出了原土培厚、黏土斜墙、复合土工膜3 种防渗改造措施，利用ABAQUS 软件，对改造前后大坝的渗流和坝坡稳定性进行了安全性核算。分析计算结果，对不同改造方案的安全性和可行性做了评估，以期为类似防渗改造工程提供参考。
1.1 流固耦合的基本方程
在多孔土体介质中，假设土颗粒和水是不可压缩的，则土体的渗透系数k与孔隙率n存在下述经验关系：
式中：k0为初始渗透系数；
α为试验常数，软黏土的取值为4.0～6.0。
若体积应变 εv全部是由孔隙体积变化引起的，则荷载作用下孔隙率函数表达式为：
式中：n0为初始孔隙率。
针对具体的土体本构模型，应力场决定着体积应变，故土体介质的渗透系数与应力场的变化有关。渗流作用以渗透体积力的形式作用于土体介质，土体的平衡方程可以表示为：
式中：σij,j为应力场对空间坐标的偏导；
fi(h)为渗透作用产生的体积力，与水头相关；
Xi为荷载。
将渗流场数值模型和应力场数值模型及渗流和应力的经验关系进行组合，并以矩阵的形式表示，得到两场的耦合模型[7-8]：
式中：k为与渗透系数有关的矩阵；
h为 水头列向量；
s为与贮水率有关的矩阵；
∂h/∂t为水头对时间的偏导得到的列向量；
f为已知水头节点得出的常数列向量；
dσ为与土体自重有关的应力增量列向量；
dε为应变增量列向量；
D为弹（弹塑）性矩阵；
B为应变矩阵；
dw为节点位移增量列向量；
K为刚度矩阵；
w为节点位移列向量；
F为与渗透力有关的列向量；
X为荷载列向量。
式（4）表示应力场影响下的渗流场；
式（5）为土体的应力-应变关系；
式（6）为描述渗流场影响下的力的平衡方程。
1.2 流固耦合在ABAQUS 中的实现
在耦合计算中，ABAQUS 将节点位移和孔隙水压力作为节点自由度进行空间离散，得到平衡方程和渗流方程的矩阵形式，直接求解耦合方程[4]。求解流程如下：
（1）将整个求解域离散为特定的位移/孔压耦合单元，确定单元和节点的编号。
（2）给定定解条件，包括初始应力、位移边界、初始水头和流量边界等。为实现应力场对渗流场的影响，在设置土体介质的材料参数时，将渗透系数设置为与孔隙率有关的函数，如式（1）。
（3）在得到某个节点的位移和孔压时，单元上其他节点的位移和孔压通过内插函数获得，从而得到整个求解域的渗透矩阵和应力矩阵。
（4）求解整个区域的有限元方程，得到各个节点的位移和孔压。提取计算结果，计算其他变量的值。
咀儿上淤地坝位于偏关县城东南40 km 处的楼沟乡咀儿上村，所在沟道属黄河一级支流——县川河二级沟道，坝址位置为东经110°31′45.88″，北纬39°16′52.29″。该坝所在沟道为“U”型沟道，大坝为均质碾压土坝，洪水标准按20年一遇设计，设计淤积年限20年。原设计坝高22 m，设计洪水位20 m，设计总库容为97.57 万m3，控制流域面积4.63 km2。大坝坝顶宽3.5 m，坝顶长94 m，上游边坡坡度为1∶2.25/1∶2.00，下游边坡坡度为1∶2.00/1∶1.75。下游坝坡设有贴坡式反滤层，层高3.0 m，顶宽1.5 m。坝体的主要材料为中轻粉质低液限黏土，土粒密度2.72 g/cm3，孔隙比0.7。坝基主要为第四纪黄土，颗粒级配均匀、黏性大、密实，干后坚硬，土粒密度2.81 g/cm3，孔隙比0.65。淤地坝坝体完好，放水建筑物运行正常，坝体下游20 m 处有截潜流水源。
2.1 原坝体渗流和稳定性计算
原坝体有限元模型如图1 所示。根据咀儿上淤地坝横断面图，建立包括坝体和坝基的二维有限元模型，坝基向上下游各延伸30 m，向下深度取20 m，坝基总长度为155 m。对该模型进行网格划分，采用扫略式划分技术，并选择CPE4P 四节点位移/孔压耦合单元，网格最小尺度为0.5 m，最大尺度为4.0 m，得到网格总数5 297 个，节点5 525 个。模型计算共分为两个计算步，第一步为流固耦合计算得到坝体和坝基内的渗流和应力应变情况；
第二步进行强度折减，坝坡失稳时的折减系数即为大坝的坝坡稳定安全系数。
图1 原坝体网格模型Fig.1 Mesh model of original dam
计算时，假设坝体、坝基和防渗材料均为理想弹塑性材料，服从Mohr-Coulomb 屈服准则[9]。各材料具体物理参数见表1。
表1 主要材料物理参数Tab.1 Physical parameters of main materials
以设计洪水位（20 m）作为上游水位条件进行数值计算，并绘制大坝等势线及流速矢量图和水平位移分布图，如图2、3 所示（在图2 中绘出了模型验证中水力学法计算所得的浸润线）。由于上下游的水位差，坝体内部形成了稳定的渗流。浸润线的位置较高，下游出逸点高度为10.67 m，且受到基质吸力[10]的影响，部分水流流入浸润线上部。坝体作为渗流的主要通道，其单宽渗流量为4.380×10−7m2/s，部分渗透水由坝体流入坝基。渗透水在上游渗入段的最快流速为7.760×10−8m/s，渗入后流速有所减小，在下游出逸段附近又重新增大到6.357×10−8m/s 左右。蓄水后，大坝在水流影响较大的地方均产生了一定的水平位移。坝体上、下游坡的渗流进出段附近流速变化较快，渗透坡降较大，使得这两部分的位移最为明显，位移的最大值为0.141 m。
图2 原坝体等势线及流速矢量（单位：渗流流速，m·s−1；
等势线，m）Fig.2 Equipotential line and velocity vector of original dam（unit: seepage velocity in m·s−1; equipotential line in m）
图3 原坝体水平位移Fig.3 Horizontal displacement of original dam
本文以有限元强度折减法计算大坝稳定安全系数，采用数值计算不收敛作为坝坡失稳的判据[11]。该大坝在折减系数为1.084 时数值计算不收敛，对应的坝坡稳定安全系数为1.084。
2.2 模型验证
考虑渗流场和应力场的耦合影响，对原坝体进行渗流和稳定性分析，所得结果分别与传统水力学二段法和条分法进行对比（表2）。
由表2 可见，流固耦合计算所得出逸点高度相对水力学二段法提高了2.39%，单宽渗流量增加了3.30%，大坝渗流浸润线大致吻合；
流固耦合计算出的安全系数比条分法减小了6.27%。两者的计算结果相差不大，证明本文用流固耦合进行坝体的渗流和稳定性分析是可行的。
表2 验证结果对比Tab.2 Comparison of verification results
出逸比降是校核下游坝坡稳定安全的重要数据，根据谢斯塔可夫经验公式[12]计算得到咀儿上淤地坝下游坝坡反滤体以上部位的出逸比降为0.57，大于流土的临界渗透比降0.50，下游坡面在渗透力的作用下将产生局部破坏，极易危害下游的整体安全。根据碾压式土石坝设计规范，小型土石坝正常工况下的坝坡稳定安全系数需大于1.25，该坝下游坝坡的安全系数为1.084，不满足规范要求。
考虑到咀儿上淤地坝下游存在耕地和截潜流，仅对上游坡进行防渗处理。由于下游边坡系数未发生改变，下游坡的出逸比降不变。为使下游坡面不发生流土破坏，各方案根据其不同的浸润线高度在坝体下游坡设置反滤排水。特征方案的改造设计图如图4 所示，各方案具体说明如下。
图4 方案B、C 设计断面Fig.4 Design drawings of schemes B and C
（1）原土培厚：采用与淤地坝坝体相同的土料进行原土培厚。30 m 以下低坝的坝顶宽度一般在5～10 m，考虑到施工条件和工程量，培厚坝顶宽度至7.5 m，上游边坡坡度分别设置为1∶2.50、1∶3.00、1∶3.50（方案A1～A3）。
（2）黏土斜墙：采用黏土斜墙防渗体防渗（物理性质参数见表1）。土质斜墙防渗体应自上而下逐渐加厚，防渗体在顶部的水平宽度不宜小于3 m，在底部的水平宽度应大于1/5 水头。改造后的淤地坝的水头为20 m，故黏土斜墙方案（B）的黏土防渗墙在坝顶的厚度设为3 m，上游边坡设为1∶3.00，坝体由内向外依次为30 cm 过渡层、黏土斜墙、50 cm 保护层、20 cm 干砌石护坡。
（3）土工膜：在上游铺设土工膜防渗体（物理性质参数见表1）。25 m 以下土石坝常用土工膜厚度在0.4 mm 左右[13]，本文选择0.5 mm 复合土工膜对原淤地坝进行防渗改造。土工膜方案（C）的坝顶水平宽度为7.5 m，上游边坡坡度设为1∶2.25，坝体由内向外分别为30 cm 垫层、0.5 mm 复合土工膜、50 cm 保护层、20 cm 干砌石护坡。各方案详细数据见表3。
表3 各改造方案详情Tab.3 Details of each transformation scheme
对5 种改造方案分别进行数值计算。由于土工膜非常薄，在建立模型时采用等效土体法[14]对土工膜进行处理，即在数值计算中按当量渗透系数把厚度很小的土工膜等效成具有一定渗透系数且厚度较大的土质材料。本文采用的土工膜渗透系数为1×10−14m/s，等效的土体垂直厚度为1.11 m，等效渗透系数为2.22×10−11m/s。各方案计算结果特征值见表4。
表4 各方案渗流稳定分析结果Tab.4 Seepage stability analysis results of each scheme
（1）原土培厚3 种方案计算所得渗流情况和水平位移类似，限于篇幅仅给出A3 方案计算结果，见图5 和6 所示。可见，3 种原土培厚方案在逐渐放缓上游坝坡的过程中，坝体渗径延长，渗流出逸点高度逐渐降低，单宽渗流量有所减小，下游部分水力坡降相对原坝体也略有下降，但总体的防渗效果不明显，各方案下游坝坡的安全系数均不满足规范要求。该淤地坝在筑坝时并未采取有效的防渗措施，且其筑坝材料的渗透系数较大，抗剪强度偏小，仅通过培厚上游边坡难以达到理想的防渗效果。
图5 方案A3 等势线及流速（单位：流速，m·s−1；
等势线，m）Fig.5 Equipotential line and velocity vector of scheme A3（unit: seepage velocity in m·s−1; equipotential line in m）
（2）黏土斜墙防渗方案的等势线、流速矢量和水平位移如图7 和8 所示。黏土斜墙良好的防渗作用，使下游部分水力坡降相对原坝体明显减小，出逸点高度下降为3.00 m。改造后坝体的单宽渗流量减小为1.715×10−7m2/s，较原坝体下降了73.17%，上游渗入段的最大流速减小为3.129×10−9m/s，下游出逸段的最大流速为5.22×10−8m/s，在黏土斜墙后有些许水流从坝基回流至坝体。坝体下游部分的水力坡降大大降低，渗透动水压力的影响随之减弱，大坝的整体位移量减小，位移分布发生变化，静水压力成为上游坡面发生位移的主要因素。上游坡面的最大水平位移为0.094 m，坝体下游部分的最大水平位移为0.070 m。铺设黏土斜墙后下游坝坡的安全系数增大到1.598，满足了土石坝安全蓄水的规范要求。
图6 方案A3 蓄水时水平位移Fig.6 Horizontal displacement of scheme A3 when impounding water
图7 方案B 等势线及流速（单位：流速，m·s−1；
等势线，m）Fig.7 Equipotential line and velocity vector of scheme B （unit:seepage velocity in m·s−1; equipotential line in m）
图8 方案B 蓄水时水平位移Fig.8 Horizontal displacement of scheme B when impounding water
（3）复合土工膜防渗方案的等势线、流速矢量和水平位移如图9 和10 所示。土工膜优越的防渗效果使得坝体在土工膜处的浸润线骤降，渗流出逸点高度下降为0.857 m，坝体内部的水力坡降很小，几乎没有受到渗透动水压力的影响。坝体的单宽渗流量减小为8.414×10−9m2/s，上游渗入段的流速均在2.573×10−10m/s 左右，下游出逸段的最大流速为4.998×10−8m/s，坝基中有渗透水从土工膜后流入坝体。蓄水后大坝上游坡面受到静水压力的影响发生了较大的位移，位移最大值为0.093 m；
坝体下游部分的位移量较小，均小于0.056 m。由于大坝整体受到渗透水的影响很小，大坝的安全性较高，折减计算后得到下游边坡稳定安全系数为1.851，满足土石坝安全蓄水的规范要求。
图9 方案C 等势线及流速矢量（单位：渗流流速，m·s−1；
等势线，m）Fig.9 Equipotential line and velocity vector of scheme C（unit: seepage velocity in m·s−1; equipotential line in m）
图10 方案C 蓄水时水平位移Fig.10 Horizontal displacement of scheme C when impounding water
（4）原坝体及各方案潜在滑裂面如图11 所示。由于下游坡并未进行改造，各方案坝体失稳的滑裂面基本相同。在坝体失稳破坏的过程中，下游坡脚与坝基的交界处最先发生破坏，随着破坏程度的加深，破坏区域向坝体内部延伸，最终形成贯通上游坝坡和下游坝脚的圆弧状滑裂面。由于地基的抗剪强度大于坝体，且坝体中的渗流水对坝体材料的强度有着不利影响，发生失稳破坏时各方案的滑裂面并未延伸至坝基中。
图11 原坝体及各方案潜在滑裂面Fig.11 Potential slip surface of each scheme and original dam
针对咀儿上淤地坝，在保持下游边坡不变的前提下，讨论通过对上游边坡进行防渗改造，使其蓄水后满足渗流和稳定性要求。提出原土培厚、黏土斜墙、复合土工膜3 种防渗改造措施，共5 种方案，并考虑渗流场和应力场的耦合作用进行渗流和稳定性计算，对比分析后得出如下结果：
（1）对原淤地坝进行原土培厚，培厚坝顶宽度为7.50 m，放缓上游边坡系数至3.50 时，渗流量减小了18.29%，安全系数增大到1.190，可见仅采用上游原土培厚的工程措施时，防渗效果不明显，难以满足土石坝安全蓄水的要求。
（2）在对原淤地坝添加了黏土斜墙防渗体后，大坝水头下降显著，坝体受渗透动水压力的影响减小，使其稳定性有了明显提高。黏土斜墙方案的单宽渗流量为1.715×10−7m2/s，较原坝体下降了73.17%，安全系数增大到1.598。
（3）复合土工膜改造方案的单宽渗流量为8.414×10−9m2/s，下游坝体几乎不受渗透动水压力的影响，蓄水后其安全系数从原来的1.084 增大到1.851。
（4）黏土斜墙方案和复合土工膜方案的安全系数均大于1.25，满足土石坝安全蓄水的规范标准，两种改造方案理论上都是可行的。但黏土斜墙方案工程量相对较大，且黏土防渗体的施工容易受气候影响，复合土工膜改造方案具有更好的安全性、可行性，可以在条件允许的前提下优先考虑。
猜你喜欢坝坡淤地坝土工膜大中型淤地坝安全度汛“四预”模型的重点主题及其算据中国水利(2022年9期)2022-05-31陕西榆林地区无定河流域淤地坝遥感解译干旱区地理(汉文版)(2022年3期)2022-05-30不同铺设及联接型式的复合土工膜围堰应力变形特性长江科学院院报(2022年4期)2022-04-16土工膜和GCL防渗结构缺陷渗漏量试验研究中国水利水电科学研究院学报(2022年1期)2022-03-29淤地坝坝系工程除险加固施工时序安排探讨中国水土保持(2022年3期)2022-03-23库水位骤降偶遇地震作用的土石坝稳定分析三峡大学学报(自然科学版)(2021年3期)2021-04-29黄土高原地区淤地坝存在问题分析农业科技与信息(2020年15期)2020-12-18软体沼气工程聚乙烯土工膜技术性能与成本分析中国沼气(2020年1期)2020-07-28非稳定渗流场对黏土心墙坝坝坡及心墙稳定的影响分析价值工程(2019年34期)2019-12-18水源工程土石坝坝体渗流及坝坡稳定性分析科技创新与应用(2019年25期)2019-09-24