土石坝坝基内部侵蚀规律研究现状述评

刘阳，侍克斌，毛海涛,2，胡洪浩

(1.新疆农业大学水利与土木工程学院，新疆乌鲁木齐830052；2.重庆三峡学院土木工程学院，重庆404100)

渗流潜蚀(Suffosion)又称内管涌，是指土体中的细颗粒在渗流作用下沿粗颗粒之间的孔隙通道移动流失的现象。管涌(Piping)被称为“翻花”或“翻沙”，是指在渗流作用下，土体细颗粒沿骨架颗粒形成的孔隙移动流失，形成土体贯穿的渗流通道的现象。从定义上来看，二者的成因相同，均是由于细颗粒流失，土体形成通道，最终破坏。实际上，两者发生部位完全不同。管涌又名溯源侵蚀或后向性侵蚀，其发展方向与渗透水流反向。管涌的发生开始于坝基下游某处薄弱点，当渗流水流流速较大时，该点的细颗粒首先被带出；随着水流的增强，较大的粒径从该口处不断被带出，出口处的空间沿渗流相反的方向扩张并形成通道；最终通道持续向上游发展，并且贯穿堤基，这就是管涌发展的全过程。然而，内管涌是由于土体中的细颗粒在渗流作用下首先达到启动条件，沿着大颗粒间的孔隙移动，并且沿着渗流方向形成通道；随着渗流流速加大，此通道在土体内部逐渐向外扩张，在一定条件下，此通道扩张至坝基外部，形成贯穿坝基的通道。可见，二者的发展方向正好相反，形成的条件也不尽相同。

内管涌常会引起水工建筑物或地基发生破坏，Foster[1]在2000年对世界范围的11 192座大型土石坝的破坏形式进行了统计分析，其中45%是由内外管涌引发的。如美国蒙大拿州的豪瑟(Hauser)湖大坝以及华盛顿州Elwha河大坝就是典型的由内管涌导致的大坝失事。在中国的大坝工程中，土石坝因其良好的适应变形能力、施工技术简单等优点，往往作为建坝的首选。据统计，2005年中国坝高30 m以上的大坝，土石坝所占比例高达67%[2]。由此可见，土石坝的安全问题不容忽视，而内管涌作为土石坝安全隐患的重要形式更是十分值得关注。

近年来，内管涌对大坝安全造成的重要影响逐渐成为讨论的热门话题。李喜安[3]、刘忠玉[4]等从理论方面推导内管涌发生的启动条件，Wan[5]、Fujisawa[6]、Jentsch[7]等从颗粒级配曲线推导内管涌发生的几何条件；Liang Y[8]、吴梦喜[9]等通过有限元模拟对大坝的安全进行计算；El Shamy[10]、Yuan Guo[11]等采用颗粒流软件对内管涌侵蚀过程进行模拟；Lin Ke[12]、倪小东[13]、A.Masoodi[14]等试图通过试验还原内管涌的侵蚀过程，将侵蚀过程可视化。可见，内管涌对大坝的安全影响问题已经引起广泛的关注，大量的学者试图通过各个方面对内管涌的侵蚀机制进行研究。

为了厘清内管涌发展现状及研究进展，便于进一步探究该领域前沿成果，本文将对内管涌的研究现状从理论、试验、数值模拟三大主流研究方向分别进行系统地梳理和分析。

1 理论研究

内管涌定义揭示了其发生必须满足2个条件：①水力条件，细颗粒能够在渗流作用下打破原有的力学平衡进行移动；②几何条件，细颗粒能够通过粗颗粒间的孔隙形成的通道移动。水力条件和几何条件一直贯穿于内管涌的全过程，当某一条件不满足时，内管涌即刻终止。图1所示，a状态到b状态即是渗透力作用下，部分细颗粒变成可流动细颗粒，逐渐流失的过程；随着时间推移，满足几何条件的土体细颗粒越来越少。由b状态到c状态可见，随着可动细颗粒沿着孔隙通道不断流走，最终只留下大颗粒骨架和极少量没有达到启动条件的细颗粒。因此，在理论研究方面，国内外绝大多数研究成果主要集中在水力条件和几何条件方面。

a)a状态 b)b状态 c)c状态图1 内管涌侵蚀示意

1.1 水力条件

内管涌是否发生主要取决于土颗粒被侵蚀的力学条件。提到内管涌的力学条件，就不得不先说河流运动学中的Shields起动拖曳力曲线。Shields在1963年对各类泥沙进行临界启动试验，得出启动剪应力与雷诺数的关系曲线。在内管涌的侵蚀研究中，对泥沙的启动条件的研究正是基于Shields曲线。

目前，对内管涌力学条件的推导，主要采用毛细管模型。Khilar[15](1985)首次将土体概化为由众多不同直径的毛细管组成的毛细管模型。在毛细管中，当颗粒所受毛管水流的剪切力超过临界剪切力，水流将其从管壁剥离将其带走。刘忠玉[4](2004)、李喜安[3](2010)、傅旭东[16](2011)等分别基于毛细管模型以临界剪切力确定侵蚀发生的方法。其中，李喜安在《渗流潜蚀作用临界发生条件的推导》一文中，采用毛细管模型(图2)对不同渗流方向细颗粒启动条件进行了推导。

图2 毛细管模型

当渗流方向向上时，细颗粒受到向上的拖曳力为：

F=6πμr[ν0(1-b2/B2)]

(1)

式中r——细颗粒半径；b——细颗粒中心到微型圆管轴线的垂直距离；B——微型圆管的半径；μ——颗粒-流体的表观黏滞系数。

细颗粒的水下重力为：

(2)

式中γs——细颗粒的重度；γw——水的重度。

设细颗粒与孔隙壁之间的黏聚力为C，则处于临界启动状态细颗粒的静力平衡条件为：

G′=F-C

(3)

将式(1)、(2)代入式(3)可得细颗粒启动条件为：

(4)

同理，当渗流方向水平时，细颗粒与管壁胶结处的抗剪强度为：

τ=σn·tanφ+C

(5)

(6)

φ为细颗粒与管壁胶结物的内摩擦角。当颗粒处于临界启动状态时，极限平衡条件为F=τ,则：

(7)

对于垂直向下的渗流，同上可得，细颗粒启动条件为：

(8)

内管涌的侵蚀是细颗粒流失过程，其受力分析须针对单个土颗粒。采用毛细管模型将圆管作为细颗粒周围的土体骨架，对单个土体细颗粒进行受力分析，能将复杂的情况简单化。因此，目前对水力条件的推导主要采用毛细管模型。但土体内部侵蚀作用发生的环境往往极其复杂，复杂的环境反映出复杂的渗流侵蚀启动过程，这一事实说明在推导内管涌侵蚀发生的临界条件时，应该根据具体情况，对症下药，显然仅仅一个毛细管模型并不能作为一个普遍适用的渗流侵蚀模型。

1.2 几何条件

当水力条件满足之后，细颗粒想要通过粗颗粒间形成的孔隙移动，其粒径R必须小于孔隙间的距离D，即：

R

(9)

因此，几何条件主要从土的颗粒组成和孔隙大小角度进行论述。在此基础上不难发现，外管涌和内管涌发生机理在几何条件上并无差别。关于外管涌几何条件的研究成果仍然适用于内管涌。20世纪50年代初，伊斯拖美娜提出，当土的不均匀系数Cu>20时，容易发生管涌破坏，这个观点在应用中存在一定的偏差，部分学者在后期大量的研究中进行改进和完善。刘杰[17]认为在Cu<5的不均匀土中，当细料含量P<25%时容易发生管涌，并给出了管涌土的临界水力梯度经验公式：

(10)

式(10)在中国国内研究管涌渗透破坏方面具有重要的指导意义。不均匀系数小的土，颗粒大小均匀，骨料与填料区别不明显，土体不能形成骨架结构，如果粒径较小，在较小的水头下，出口处小颗粒被水流带走，并逐渐发展为流土破坏。如果不均匀系数较大，曲率系数也较大，土体缺乏中间粒径，在土体结构中，大颗粒形成的骨架结构空隙较大，小颗粒在大空隙中处于自由状态，在较小的水力坡降下，小颗粒土体被带走，形成管涌破坏，这种破坏形式由于粗颗粒骨架还存在，还能承担一定的水力坡降，土体结构并未完全破坏。

1985年Kenney[18]指出无黏性砂的内部稳定性取决于土体的粗粒部分所组成的土体骨架能否阻挡细颗粒被渗透水流带出。Wan[5](2008)改进了以级配曲线判断土体内部稳定性的方法。Fujisawa等[6](2010)根据土体的渗透系数及孔隙率确定了土体孔隙的代表直径，认为小于代表直径的颗粒只要满足力学条件就能被带走。Jentsch等[7](2014)提出通过土体级配曲线判断土体内部骨架结构，以此分析土体内部稳定性。

以上从土的不均匀系数、颗粒级配角度对土体内部稳定性进行研究，说明土体内部的侵蚀机制(内管涌的发生机制)与土体本身的性质密切相关，也即是内管涌发生的几何条件很大程度上取决于土体本身的性质。

2 数值模拟

2.1 有限元数值模拟

有限元可以用有限的、相互关联的单元模拟无限的复杂体，使复杂的工程问题简单化。由于国内外对内管涌的研究起步较晚，目前利用有限元模拟内管涌的发展尚在探索阶段。

胡亚元[19](2013)基于Galerkin编制有限元程序，建立三相耦合的侵蚀控制方程，在轴对称情况下对土体内部管涌侵蚀过程进行数值模拟。张兴胜[20](2015) 利用有限元法计算内管涌对心墙坝造成的影响，结果证明，在内部管涌侵蚀的作用下，随着侵蚀系数α的增大，坝基产生的位移也会随之增大。Liang Y[8](2016)基于自适应移动网格，提出了一种模拟土体内部侵蚀的数值方法，该方法可以使流网自动适应土体的离散域，通过流网的变化模拟内管涌侵蚀路径。吴梦喜[9](2017)依据内管涌过程中泥沙输运的连续性方程和含沙水渗流的连续性方程，采用有限元法，通过顺序非线性耦合迭代，求解泥沙侵蚀过程中的渗流场和浓度场，并记录各时间步有限元节点的孔隙水压力、高斯点的泥沙浓度和体积侵蚀(淤积)率，得到内管涌的动态过程。

其中，吴梦喜将泥沙运输方程和含沙水的渗流方程进行耦合，实现对内管涌过程的动态模拟，同时采用单向双向耦合对大坝的渗流和变形进行计算。不仅在理论上进行了创新，而且能准确计算出内管涌造成大坝坝体和坝基的位移，直观地分析出内管涌侵蚀对大坝造成的影响程度。

根据质量守恒定律，泥沙输运的连续性微分(质量守恒)方程为：

(11)

式中Sv——孔隙水中泥沙体积浓度；φ——土体孔隙率；vi——渗流速度向量；q——单位时间内单位体积中泥沙的沉积体积。

假设泥沙的侵蚀速率与作用在土体的渗透坡降之间满足线性关系，即：

q=kd(J-Jc)

(12)

式中kd——土体侵蚀系数；J——土体的渗透坡降；Jc——砂粒启动所需的临界坡降。

含沙水的连续性方程：

(13)

式中θ——单位体积土体的孔隙体积含水率。

不同孔隙水饱和度和含沙浓度时的渗流，可以用下式表示：

νi=-Kij(β,Sν)kr(s)(ψ+z)

(14)

式中ψ——压力水头；vi——渗流速度向量；Kij——饱和介质渗透张量；β——泥沙体积侵蚀率；Sv——含沙浓度；kr(s)——相对渗透系数；z——参考平面高程。

含水率和孔隙水压力混合表达的含沙水渗流微分方程：

(15)

此模型能够得出侵蚀过程中坝基每一个时刻的渗流场和浓度场，借此求出该时刻的侵蚀体积率，最终计算出侵蚀后坝体、坝基的位移和应力增量,是目前最新计算内管涌侵蚀的有限元模型。

众多学者采用有限元进行模拟时，为了模型建立的简单，通常假定当水库蓄水达到正常蓄水位时内管涌侵蚀开始发生，此时开始计算。事实上，在达到正常蓄水位之前侵蚀早已发生，这样计算得出的结果和实际结果有着一定的偏差。在今后对内管涌的有限元计算中，应该从实际发生侵蚀的时刻开始计算以保证计算结果的准确性和可靠性。

渗流的有限元计算，即使从最初的纯渗流模型发展到土-水耦合、三相耦合模型，仍然无法完全模拟出土的实际状态。且岩土体的物理参数存在随机性、模糊性和渗流的混沌性，因此有必要建立多场耦合模型。值得一提的是，目前工程中使用的耦合理论和本构方程绝大多数是各向同性的模型，如何建立能够反映土体实际情况并简单、实用多场耦合模型是值得注意和研究的问题。除此之外，为了合理地描述土层中多孔介质中多相或多种物理、化学场与环境的复杂相互作用过程，需深入研究多场耦合机制和耦合效应定量模型，使内管涌的研究更加科学和系统。

2.2 颗粒流数值模拟

内管涌的发展过程中，土体颗粒和水流相互作用的具体力学机制非常复杂，具有很大的随机性，影响的因素也很多。为了从颗粒尺度系统阐述内管涌发展的微观机理，众多学者逐渐选择以颗粒流模型模拟内管涌的发展过程。

El Shamy[10](2008)考虑被侵蚀的土颗粒、水流、土骨架之间的相互作用，采用N-S方程求解连续流体，用颗粒流软件模拟发生内管涌侵蚀的细颗粒，其数值模型见图3。

图3 N-S方程离散元模型

夏尔博士[21](2015)，第十六届欧洲土力学及岩土工程研讨会上展示了利用离散单元模型检测内管涌破坏的研究成果，证明了利用离散元检测内管涌的破坏情况是可行的。Y.Guo，S. M.[21](2017)利用耦合的CFD-DEM方法，模拟不同形状的土颗粒在渗流作用下的侵蚀情况，结果证明，颗粒的形状和角度对侵蚀的阻挡起着重要的作用。Yuan Guo[11](2017)将计算流体动力学(CFD)和离散元素法(DEM)的计算模型进行耦合，利用CFD-DEM对土体内部侵蚀进行模拟，对冲刷深度和土体侵蚀率，进行预测。证明了颗粒流软件能有效计算内管涌对大坝的危害程度。蒋中明[23](2018)利用颗粒流软件PFC3D对粗粒土的渗透过程进行模拟。并发现随着水头的增加被侵蚀的细小土颗粒，随渗流从上游向下游汇聚。再现了内部管涌侵蚀的全过程。张凤德[24](2018)等，采用PFC3D，对上覆黏土，下部为级配不良细砾的土层的坝基进行侵蚀模拟，分别模拟了在砾石层和砂石层由内部管涌侵蚀发展为管涌和流土的过程。Kuang Cheng[25](2018)等建立半解析的CFD-DEM模型，模拟内管涌侵蚀作用下，细颗粒的迁移过程，其模型见图4。

图4 CFD-DEM侵蚀模型

利用颗粒流方法能够将从微观的角度将土体的侵蚀过程模拟出来，直观地反映出侵蚀的机理。在内管涌侵蚀的研究中有着很大的优越性。但是想要用颗粒流方法模拟出整个大坝的侵蚀状况，需要设置的颗粒数一般远超百万，对计算机的要求非常高，现有的技术手段无法按照实际工程中实际颗粒组成进行整体模拟。因此，现在利用颗粒流软件研究内管涌侵蚀多停留在局部侵蚀，其结果是否能够代表整个大坝的侵蚀情况还有待考证，急需建立起微观与宏观之间的联系。

3 试验

在实际工程中，内管涌的发生及发生过程无法观测，相关数据更是无从获取，因此，多年来大量的学者开展试验研究，主要目的是将内管涌过程可视化，发现其物理现象。Stewart[26](2002)在对大坝的安全研究中，采用CT技术，进行试验，对内管涌侵蚀进行观测。B.Burns[27](2006)提出利用微型电阻率阵列探头利用电阻率的变化对大坝内管涌的侵蚀过程进行监测，试验模型见图5。试验证明电阻率的变化能与内管涌的扩散过程相对应，这种远程传感方法在内管涌的监测中是有效的。

图5 电阻率阵列探头及侵蚀试样

Bendahmane[28](2008)设计了图6所示的试验装置，该装置能够使土体处于三向受压，一定程度上还原了土体的原状受压状态，反映土体内部侵蚀规律。

图6 三轴侵蚀实验装置

Lin Ke[12](2012)通过开展一维向上渗透试验，研究内管涌侵蚀对土体强度的影响。结果证明土体内部细颗粒损失随着施加的水力梯度增加而增加，且内管涌侵蚀会导致土体的抗剪强度降低。YA Sharif[29](2015)通过试验研究土体组成对坝基侵蚀的影响，结果表明，黏土的含量直接影响侵蚀速率。倪小东[13](2018)等引入多色填充砂，对土体的侵蚀过程进行追踪，直观地揭示了内部管涌侵蚀的发展过程。A.Masoodi[14](2018)通过短程摄影测量法对内管涌侵蚀试验过程进行记录，依据影片的高分辨率，分析内管涌侵蚀程度。

通过试验研究内管涌的侵蚀机制，能够将较好的还原整个大坝的实际状况。直接观察到大坝内管涌侵蚀的发生、发展、溃坝的全过程。但是，即使是可视化试验仍然很难在试验中观测到土体内部侵蚀的全过程，且在试验模拟中的关键参数往往由于技术手段难以获取，很难将试验定量化分析。

4 结论与展望

内管涌是威胁土石坝和堤坝安全的一个主要潜在问题，本文从理论、数值模拟和试验3个方面对国内外内管涌的研究现状进行归纳总结。现将目前研究中存在的问题进行讨论如下：在理论研究方面，现在尚无对内管涌准确的定义；在数值模拟方面，有限元数值模拟，尚无法做到完整的计算全过程；颗粒流数值模拟，由于对计算机设备要求的严格，能够观察到侵蚀的微观现象，却难以模拟整个工程，工程实用性价值有限。在试验研究方面，能够通过其可视化特点观察到侵蚀规律，但难以获得孔隙率、应力应变等关键参数，对大坝整体安全性进行分析。利用CT、射线等方法监测整个大坝造价高昂，对试验设备和试验技术的要求较高。高清数码摄像技术难以把控全局。

鉴于上述观点，结合国内外有关内管涌侵蚀的相关会议，建议今后对大坝内管涌侵蚀的研究可从以下方面开展。①理论方面，从内管涌发生的几何条件和力学条件进行理论推导，得出明确的内管涌触发机制。数值模拟在当前内管涌的研究中占有重要地位，利用数值模拟可以解释侵蚀试验中难以解释的现象。完善有限元数值模型，准确、完整地对整个大坝侵蚀全过程进行计算；改进现有的颗粒流软件，降低对计算机性能的要求，模拟整个大坝侵蚀的全过程。依托试验开展数值模拟，二者相互验证，保证结果的准确性。②试验方面，结合各种先进技术，用CT、射线等方法观察内管涌侵蚀发生时细颗粒变化特征，同时采用高清摄像追踪细颗粒的发展过程，结合砂砾石侵蚀前后的三轴试验和侧限压缩试验参数结果，综合评定细颗粒流失对大坝的安全影响。