南水北调高填方段渠道渗流破坏机理研究

王平易，薛满宇，张 阳，吴多贤，刘汉声

(1.中国电建集团贵阳勘测设计研究院有限公司,贵州 贵阳 550081; 2.北京京能国际控股有限公司华北分公司，北京 100020；3.华北水利水电大学 水利学院，河南 郑州 450003)

南水北调工程是优化我国水资源时空分配的重要措施，在我国的宏观经济和社会发展中起着十分关键的作用。但由于南水北调工程规模大、跨度长，一旦失事将会产生严重的后果。渗流问题是常见的工程问题，自世界各国大型水利工程建成以来，都陆续出现了引发重大工程事故的渗流问题，有些事故不仅会造成经济损失，还会威胁到人们的生命。因此，对水利工程进行渗流与稳定分析是必不可少的。

胡元进[1]模拟三种不同工况下高边坡在降雨入渗条件下的稳定性，得出不同工况下,随着降雨入渗时间的不同,边坡应力应变反映出的不同分布特征。李广权等[2]对含有软弱下卧层的边坡稳定性进行数值模拟计算﹐探究了上部与软弱下层土体的抗剪强度、弹性模量、泊松比、密度对边坡稳定性的影响。张家发[3]建立了三维饱和、非饱和稳定渗流数学模型，并编拟了有限元计算程序，初步说明了土坝饱和非饱和渗流分析的意义。陈丽刚[4]通过有限元分析，研究了考虑水渗流作用的有限元强度折减边坡稳定安全系数的计算方法。李海艳等[5]采用饱和、非饱和渗流模型对土石坝渗流场进行数值模拟,分析了渗透系数的各向异性、渗透系数大小以及水位升降速度对饱和非饱和渗流场的影响。彭新峰[6]针对路面早期水损害以及降雨条件下边坡失稳问题，研究了在降雨期间和降雨停止后的路基路面渗流域饱和—非饱和渗流场的变化以及渗流场的变化对边坡稳定的影响。

本文主要对南水北调淅川段高填方渠段渗流场进行仿真分析，对渠道运行的稳定性进行分析，提供安全运行监测的参考依据，及时发现不安全因素，指导运行监测的工程措施，具有重大的经济和社会效益。

1 渗流计算原理及方法

从渗流模型中取出一个微分单元，其体积V=dxdydz。假设该单元中x方向的流速为vx，y方向的流速为vy，z方向的流速为vz。同时假设流体不可压缩，且在渗流过程中土壤的孔隙比保持不变。根据达西定律，可知三个方向上土体的渗流速度，进而可以推理得到三维渗流的连续性方程，如式(1)：

(1)

式中：kx、ky、kz分别为x、y、z方向上单元体的渗透系数，cm/s；φ为渗透水头函数，m。

该渠道模型形状复杂，为使计算结果精确，本文采用ABAQUS进行仿真模拟计算，并采用六面单元体进行网格划分，其渗流计算微分单元体格式如图1。

图1 渗流计算微分单元体示意图

2 工程概况和数值模型

选取南水北调中线工程淅川段高填方段渠道工程为实例进行研究，各层土的力学性能如表1所示。建立有限元模型，见图2。模型总共分为三层建立，给出了模型的整体有限元网格和渠道的网格划分。渠道有限元模型中列出了渠道及各层地基土层分布，并对渠道填土、各层地基土进行了组元分类，可实现程序计算的自动赋予渗流参数。

表1 淅川段典型断面各层土的力学性能

图2 淅川段危险段的有限元模型

3 渗流稳定分析

本文选取两种工况进行分析，分别是设计水位工况和加大水位工况。且每种工况采用两个时刻，分别是初始时刻和最终时刻。通过对于每种水位的两个时刻对比，明确在南水北调渠道运行过程中，渗流量、浸润线等随着时间的变化情况。并判断其是否对渠道基础造成管涌、流土等破坏。此外，对于最终时刻的两种水位也要进行比较，可以得出渗流量、浸润线、等水头线、以及流速分布等随着水位的增加的变化情况。

3.1 渗流浸润线与斜坡逸出点

设计水位工况的浸润线如图3(a)所示，从图中可以看出，浸润线基本上向坡脚倾斜，浸润线的位置相对不高。渠段外侧无水，渠段中浸润线的临水侧水位为145.103 m，渠段背面逸出点位于堤岸左侧和右侧。渠道两侧的逸出点高程较低。典型断面在渠道外坡脚的出逸高度很小，其左右出逸段比降分别为0.26和0.26。逸出段的渗透比降小于土的允许比降，不会发生渗透破坏。

加大水位工况的浸润线如图3(b)所示，其特征基本与设计水位特征相同，渠段中浸润线的临水侧水位为145.773 m，由于渠道和下层土渗透系数相差较小，水位加大后对渗流场分布影响不大。典型断面在渠道外坡脚出逸，左右出逸段的比降仍为0.26和0.26，不会发生渗透破坏。

图3 渗流稳定场分布

3.2 渗流稳定性分析

图4显示了设计水位工况渗流场流速的分布。渠底仅一层较厚的粉质黏土层，黏土层的渗透系数较渠道回填土小，渠底向土层的渗流量较大。渠底向黏土层的下渗速度更快，增加了渠底的水流外渗。从图中可以看出，渠道中临水面的渗透速度也很大，渠底土层渗流场的流速远大于渠道中的流速，除此之外，渠底第一层土内和两侧坡脚处最大。由于两侧渠道断面的对称性，左右两侧坡脚外渗流速一样。图5显示了加大水位工况渗流场流速的分布。水位加大后，与设计水位时相同，两侧坡脚流速相同。

图4 设计水位工况渗流流速场分布

图5 加大水位工况渗流流速场分布

3.3 渗流量

经过分析计算，在设计水位下，渠道底部和坡脚的最大单宽渗流量为0.023 m3/(d·m)。在加大水位下，渠道底部和坡脚的最大单宽渗流量为0.025 m3/(d·m)。加大水位后，虽然提高了浸润线，加大了渗透比降，但由于初始的渗流量较小，在坡脚的渗流量增加也较小。

4 安全系数分析计算

有限元强度折减法是在外载一定的情况下，对土体的内聚力和内摩擦角进行相应的折减，进而得到相应的折减系数，当边坡出现塑性连续破坏或产生较大的位移的时，其折减系数与该边坡的稳定安全系数相等[7-8]。在实际分析中，通过假定不同的强度折减系数Fr，然后用折减后的参数进行有限元分析，然后绘制出折减系数随位移的变化曲线，观察曲线是否收敛，当土体达到临界破坏时，对应的强度折减系数Fr与边坡的安全系数Fs相等。有限元强度折减法的计算式(2)～式(3)如下[7-8]：

(2)

(3)

式中：cf为土体的临界内聚力，kPa；c为土体的内聚力，kPa；Fr为强度的折减系数；φf为土体的临界内摩擦角，(°)；φ为土体的内摩擦角，(°)。

本次研究应进行设置地应力平衡分析步，对土体在正常情况下的受力情况进行分析，然后设置通用静力分析步，对土体的强度进行折减，以此来求得边坡的稳定安全系数。通过对边坡模型失稳计算得出，特征点坡顶处Fr随U1的变化关系如图6所示。从图中能够看出，如果判断边坡稳定的标准为数值计算不收敛，那么Fr=1.52，即边坡的稳定安全系数是1.52，远大于极限平衡法中边坡的稳定安全系数，说明该判据对边坡的安全性估算过高。如果判断标准是特征点水平位移的拐点，那么Fr=1.24，即边坡的稳定安全系数是1.24，较接近极限平衡法中边坡的安全系数。比较上面两种判断标准，得出本论文应使用有限元强度折减法的第二种判据确定边坡的安全系数较合理，即通过特征部位坡顶的位移拐点来确定边坡稳定安全系数，且由此方法所确定的安全系数比用数值计算收敛评判的结果低。

图6 折减系数Fr变化曲线

图7给出了该边坡不同折减系数下的塑性区分布情况，如果边坡失稳时塑性区会发生相互贯通的现象，那么通过对比可以认为边坡失稳时的强度折减系数在0.94～0.99，取Fs=0.98，而极限平衡法中边坡的稳定安全系数是1。因此，由塑性区的贯通现象所确定的安全系数也合理。

图7 不同折减系数下的塑性区分布

土体失稳下的位移等值云图，如图8所示。从图中我们可以看出，滑移面的位置并不算清楚，因此需要通过增量位移等值云图(如图9所示)才能确定出滑移面的具体位置。与极限平衡法中的滑移面位置和形状一样，增量位移等值云图确定的滑移面基本呈现圆弧状且经过坡脚位置。

图8 土体失稳下的位移等值线云图

图9 增量位移等值云图

5 结 论

(1)本文通过对运行期高填方渠道危险渠段渗流场的仿真分析，可以认为当典型渠段处于设计流量和加大流量的正常运行条件下时，浸润线的位置不高，整体渗流量较小。外坡脚的出逸高度通常为零或者更小，并且计算出的出逸段渗透比率小于土壤的允许比率。外坡脚出逸部分的渗透比率也小于土壤的允许比率，且渠道渗流是稳定的。即当在设计流量的水位和加大流量的水位下运行时，该高填方渠段是安全可靠。

(2)运用ABAQUS软件建立有限元模型，通过有限元强度折减法来分析边坡稳定性，得出边坡实际失稳过程和模拟塑性区的发展过程相同，滑动面基本呈现圆弧状且经过坡脚位置；并证明以发生位移突变作为主要判断依据确定的边坡的安全系数较合理。