多层堤基结构管涌动态发展的数值模拟

王 霜，陈建生，钟启明

(1.河海大学水利水电学院，江苏 南京 210098；2.南京铁道职业技术学院铁道工程系，江苏 南京 210031；3.河海大学土木与交通学院，江苏 南京 210098；4.南京水利科学研究院岩土工程研究所，江苏 南京 210029)

2020年汛期，我国南方地区发生多轮强降雨，造成多地发生严重洪涝灾害，江西赣江大堤、安徽河段堤防等多地出现管涌险情。管涌是江河大堤在汛期的常见险情之一，是造成大坝失事的一种很重要的渗透破坏形式[1-3]，工程意义上的堤基管涌通常将接触冲刷、流土等渗透变形现象也归结其中。

多年来，众多学者采用试验、数值模拟等方法对管涌问题进行了广泛深入的研究。刘杰等[4-5]在室内进行了由渗流引起的堤基管涌试验，论证了管涌的发生以及发展过程对堤坝安全稳定性的影响。李广信等[6]采用彩色砂作为示踪材料，模拟了二元堤基结构及不同防渗墙深度下堤基内部管涌的发生与发展过程。罗玉龙等[7]对管涌发展机理提出了新认识：管涌是涉及孔隙水渗流、可动颗粒运移、多孔介质变形等众多复杂行为的多相多场耦合过程。Rochim等[8]在水力梯度控制或流量控制条件下，对4种不同的无黏性土进行了渗流管涌试验。研究表明，水力加载历史对管涌临界水力梯度有显著影响。明攀等[9]利用声发射技术实时监测堤防的管涌过程，发现管涌过程的水力参数和声发射参数具有相同的分布规律，可通过声发射参数对管涌过程进行判别。

随着计算机应用的迅速发展，数值模拟渐渐成为管涌研究的重要方法。殷建华[10]在国内首次采用有限元方法来研究管涌问题，计算了管涌发生区域的长度和渗透系数对堤身渗流的影响。倪晓东等[11]运用颗粒流模拟揭示了非稳定流作用下管涌发展的细观机制。Hu等[12]采用计算流体力学和离散元耦合的方法，对缺乏中间粒径的管涌型土和级配良好的土体进行了管涌侵蚀的微观和宏观研究。贾恺等[13]考虑了双层堤基管涌通道断面扩展和上溯与通道尖端坡降之间的相互影响，通过编制有限元程序模拟研究了管涌通道发展趋势。

目前，关于管涌的研究对象大多集中在双层堤基以及堤基内部土体。事实上，堤防总是沿江河而建，江河两岸由于历史上的河流冲刷、淤积和改道等，堤基结构存在很大差异，大体上可以概括为单层、双层、三层和多层4种类型[14-15]。对于不同的堤基结构，管涌的发生和发展过程存在很大差异。本文针对多层堤基结构，利用有限元分析软件对堤基管涌的动态发展过程进行模拟。颗粒流失影响局部土体渗透性，进而改变渗流场的分布以及堤基内各区域的水力梯度，继而又会影响土体颗粒流失，模拟过程中考虑了土颗粒流失和土体渗透性及其相互影响、相互耦合的过程。

1 管涌临界水力梯度确定

通过有限元数值模拟得到不同上游水头下堤基内部的水力梯度分布情况，将堤基内部水力梯度的分布与该区域土体发生破坏的临界水力梯度进行比较，从而确定在该上游水头下堤基内部是否有颗粒侵蚀，若无侵蚀，继续逐级增加上游水头，最终得到堤基内部土体发生侵蚀的上游水头，也就是管涌发生的临界上游水头，同时可以得到该水头下堤基内部颗粒侵蚀的位置与区域。

本文以细砂夹层埋入砂砾石深度为1 m的多层堤基为例，各层土体从上到下依次为黏土层、砂砾石层、细砂夹层、砂砾石层。砂砾石层以及细砂夹层的土料级配如图1所示，土体的物理参数如表1所示，上覆黏土的渗透系数为5×10-6m/s。

表1 堤基土体物理参数

图1 堤基内土体级配曲线

多层堤基中，由于渗透系数较小的细砂夹层的存在，将堤基的砂砾石层分成了细砂夹层上部与细砂夹层下部两部分。堤基渗透破坏发生的标志主要有两种形式：细砂夹层上部，管涌口下端砂砾石发生竖向侵蚀破坏，只有管涌口处颗粒首先沿竖向发生侵蚀形成空洞，才能为水平向的颗粒运移提供足够的空间，侵蚀才能沿水平向逐渐向上游发展；细砂夹层下部，细砂夹层冲破，发生竖向流土破坏，同样细砂夹层竖向流土破坏为水平向发展提供了空间。多层堤基管涌发生的标志主要取决于这两种渗透破坏形式哪一种的临界水力条件先得到满足。

对于管涌口下端砂砾石竖向侵蚀破坏，其水力梯度Jc[16]为

(1)

式中：d为管涌时被冲蚀细粒粒径，根据试验观测可取质量分数小于3%颗粒对应的粒径d3，cm；k为砂砾石的渗透系数，cm/s；n1为砂砾石的孔隙率。将砂砾石土料的各参数代入，得到砂砾石颗粒发生竖向侵蚀的临界水力梯度为0.36。

(2)

细砂夹层及其上部土颗粒发生竖向流土破坏的临界水力梯度J为

(3)

将堤基土体的物理参数代入式(2)计算可得，当细砂夹层埋深为1 m时，细砂夹层及上部土体被冲破所需的细砂层底部竖向水力梯度为2。

本文采用逐级提升水头的方式，模拟管涌集中渗漏通道的发生与发展过程。当堤基管涌口下端砂砾石的水力梯度小于0.36或者细砂夹层底部水力梯度小于2时，渗透破坏未发生，颗粒不发生流失，该种情况文中未做分析。

2 多层堤基管涌发展过程的模拟

2.1 上游水头H=14 m

采用有限元软件对细砂夹层埋深为1 m的多层堤基内部渗流场进行模拟，模型以向堤内方向为x轴正方向、向上为y轴正方向建立坐标系，堤基向两边无限延伸，模拟的是半无限边界，上游侧水头的边界条件作用于模型左侧，下游侧水头的边界条件作用于模型右侧的黏土层表面，假定为零势面。本文不考虑上部堤防渗流对堤基内部渗流场产生的影响，假定堤防为不透水体，同时假定距堤脚26 m处为黏土层薄弱处，薄弱区域直径为4 m，即在模型中将x轴坐标为84～88 m的黏土层定为管涌口，渗透系数设定为0.5 m/s,计算结果如图2和图3所示。

图2 多层堤基管涌口形成后水头等值线(单位：m)

由图2可以看出，细砂夹层底部水头分布较均匀，断面流量较小，管涌口下端的细砂夹层承担了较大的水头，细砂夹层底部y方向的水力梯度最大，最大水力梯度达到2.14(图3(a))，超过了计算所得细砂夹层冲破的临界水力梯度2，细砂夹层冲破，发生管涌破坏。而管涌口下端黏土层与砂砾石层界面砂砾石所承担的水头差非常小，竖向水力梯度几乎为零(图3(b)y=28 m处)，小于颗粒竖向运动的允许水力梯度，从而认定管涌口下端黏土层与砂砾石层界面竖向侵蚀破坏并未发生。因而，对于细砂夹层埋深为1 m的多层堤基，管涌破坏发生的标志为细砂夹层冲破，临界上游水头为14 m。

图3 不同区域y方向水力梯度

由于细砂夹层的颗粒粒径较小，管涌发生时，假设管涌口下端的细砂夹层颗粒被冲散，全部带出管涌口。细砂夹层下部砂砾石的冲蚀深度同样根据竖向水力梯度来确定，由图3(b)可以看出，细砂夹层下部离细砂夹层越近，竖向水力梯度越大，细颗粒侵蚀粒径也就越大，因而可以假设细砂夹层冲破时夹层下部1 m深度的砂砾石和上部的砂砾石被冲散，细颗粒被全部带出，粗颗粒沉积下来。

管涌口在涌砂过程中，由于细砂夹层颗粒被冲散，带出管涌口，上部砂砾石层颗粒下沉，管涌口下部形成与细砂夹层相等厚度的1 m深空洞，空洞的渗透系数与管涌口保持一致，设为0.5 m/s。细砂夹层上部被冲散的砂砾石下沉填充了原细砂夹层的位置，与细砂夹层下部1 m深度范围内的冲蚀砂砾石形成了一个侵蚀区域，侵蚀区域如图4所示，在管涌口下方y=25～27 m粗实线区域内，侵蚀区域内细颗粒流失，涌出的细颗粒引起相应区域孔隙率发生改变，则

(4)

式中：nt为管涌侵蚀区的孔隙率；Vvt为涌砂体积；Vv为管涌侵蚀区初始孔隙体积；V为管涌侵蚀区总体积；n0为初始孔隙率。

根据Kozeny-Carman方程[17]，土体渗透系数与孔隙率有如下关系：

(5)

土体未发生侵蚀前，土体的初始孔隙率为n0，初始渗透系数为k0，则

(6)

(7)

根据上述假设，管涌口下端砂砾石中的细颗粒全部流失，同时结合土体的初始颗粒级配，计算可得颗粒发生流失后的土体孔隙率nt以及土体等效粒径Dht，代入式(7)得管涌口下端侵蚀区域砂砾石的渗透系数为2.8×10-2m/s。

图4为细砂夹层冲破后x方向的水力梯度等值线图，可以看出，水力梯度最大的区域发生在管涌口前端细砂夹层底部附近。为了进一步研究多层堤基管涌发生后管涌通道的位置以及管涌通道是如何发展的，需要对不同区域土体颗粒发生运移的临界水力梯度进行分析和确定。

图4 细砂夹层冲破后x方向水力梯度等值线

对于堤基内部砂砾石颗粒在水平渗流下所引起的水平管涌，其临界水力梯度JH[16]为

(8)

式中：α为形状系数，根据明兹[18]的试验资料，对于各种砂粒α=1.16～1.17，对于锐角颗粒α=1.50～1.67；d为管涌时被冲动的细粒粒径，根据试验观测可取d3；tanφ为砂砾石的内摩擦系数；γs1为砂砾石的容重;d3K为等效粒径，可采用下式计算[18]：

(9)

式中：A为层流阻力系数，A=5.1[18]。将式(9)代入式(8)并采用一些常用数据γw=9.81 kN/m3、γs1=26 kN/m3、g= 981 cm/s2、ν= 0.011 6 cm2/s、α=1.165，可得管涌临界水力梯度的近似公式：

(10)

将砂砾石土料的各参数代入式(10)，可得砂砾石颗粒发生水平向管涌的临界水力梯度为0.24。

图5(a)为黏土层与砂砾石层界面x方向水力梯度分布图，可见水力梯度在管涌口边界x=84 m处有所减小，该处水力梯度的减小主要受管涌口的影响，管涌口处的渗透系数很大，因而水头差较小，水力梯度很小；在管涌口前端，由于管涌口处的水头为0，水头在管涌口前端发生突变，水力梯度最大，达到了0.29，大于颗粒发生水平向管涌的临界水力梯度0.24，因而在黏土层与砂砾石层界面颗粒发生侵蚀，管涌通道在此界面最终发展到x=79 m处。将管涌口前端的水力梯度0.29代入式(10)反算，可得到颗粒粒径小于0.15 mm的颗粒被侵蚀。为了简化计算，假设管涌通道上不同位置处的颗粒侵蚀均相同，侵蚀颗粒均为粒径小于0.15 mm的颗粒。将颗粒侵蚀后的土体孔隙率nt以及等效粒径Dht代入式(7)，同时结合土体未被侵蚀时的初始参数，最终可计算出管涌通道的渗透系数为1.68×10-3m/s。

图5 细砂夹层冲破后不同区域x方向水力梯度

图5(b)为细砂夹层底部与砂砾石层界面x方向水力梯度分布图，可见细砂夹层底部水力梯度大于临界水力梯度的区域为x=79～84 m，因而管涌通道在此界面形成，且通道在细砂夹层底部发展到x=79 m处。为了简化计算，假设细砂夹层底部和上部的砂砾石颗粒侵蚀相同，细砂夹层底部的通道渗透系数同样采用1.68×10-3m/s。在该水头下，细砂夹层上部和下部的砂砾石均发生侵蚀，也就是深层管涌破坏与浅层管涌破坏同时作用。

对于堤基内部细砂夹层在水平渗流下所引起的流土破坏，细砂颗粒基本上可认为是均一粒径的，故式(8)中的d3K=d，于是得到水平向流土破坏的临界水力梯度JB[16]计算公式：

(11)

式中：tanφ2为细砂的内摩擦系数，细砂的内摩擦角φ2可取28°。将细砂的各参数代入式(11)可得细砂夹层发生水平向流土破坏的临界水力梯度为0.54；根据细砂夹层水平向流土破坏的临界水力梯度与细砂内摩擦角之间的关系，同样可以计算出细砂夹层发生竖向流土破坏的临界水力梯度为1.0。

细砂夹层冲破后，细砂夹层底部的水力梯度最大(图5)，细砂夹层如果发生流土破坏，发生在细砂夹层底部的可能性最大，然而，由图5(b)细砂夹层底部x方向水力梯度分布可知，细砂夹层底部水力梯度小于发生水平向流土破坏的临界水力梯度，因此，在该上游水头下细砂夹层不会发生水平向的流土破坏。

2.2 上游水头H=18 m

在上游水头为14 m时，管涌通道形成，管涌通道在细砂夹层上部和下部均发展到x=79 m处，通道的渗透系数为1.68×10-3m/s，管涌口下端被侵蚀的砂砾石的渗透系数为2.8×10-2m/s，保持上游水头14 m不变，管涌通道不继续发展，此时提升上游水头至18 m，得到管涌通道形成后在该水头下x方向的水力梯度等值线如图6所示。由于通道内颗粒侵蚀的粒径较小，渗透系数增大的幅度不大，因而水力梯度最大的区域仍发生在管涌口前端细砂夹层底部附近，而不是管涌通道前端。

图6 管涌通道形成后x方向水力梯度等值线

土体破坏区域的分析与上述分析过程相同，此处不再赘述。图7(a)为黏土层与砂砾石层界面x方向水力梯度分布图，可见水力梯度在管涌通道前端x=79 m处减小，然后在管涌通道区域内又会有小幅增大，在管涌口边界水力梯度再次减小，由此可知，水力梯度的变化规律与该界面各区域土体渗透系数的不同有关，水力梯度在渗透系数增大处发生减小；在18 m上游水头下，管涌通道在此界面继续发展到x=70 m。

图7(b)为细砂夹层底部与砂砾石层界面x方向水力梯度分布图，同样地，水力梯度在管涌通道前端和管涌口边界减小；在18 m上游水头下，细砂夹层底部砂砾石中的管涌通道继续发展到x=73 m，细砂夹层未发生水平向的流土破坏。

图7 管涌通道形成后不同区域x方向水力梯度

由于管涌通道内水力梯度分布并不相同，因而颗粒侵蚀的粒径不同，最终使得通道内各处的渗透系数不同，为了简化计算，假设土体侵蚀后，管涌通道内的渗透系数保持1.68×10-3m/s不变；同样地，管涌口下端的砂砾石中有颗粒不断地被带出管涌口，也有水平方向来的颗粒不断补充，其渗透系数是复杂多变的，为了简化计算，假设管涌口下端砂砾石的渗透系数为2.8×10-2m/s不变，在下文的数值模拟中均保持管涌通道的渗透系数和管涌口下端砂砾石的渗透系数不变，不再重复说明。

2.3 上游水头H=22 m

保持上游水头18 m不变，管涌通道在细砂夹层上部和下部分别发展到x=70 m和x=73 m后不再继续发展，此时提升上游水头至22 m，图8为该上游水头下x方向的水力梯度等值线图，水力梯度最大的区域仍然发生在管涌口前端细砂夹层底部附近。由图9(a)可知，水力梯度同样在管涌通道前端和管涌口边界减小，在管涌通道内增加；通过水力梯度分布可以看出管涌通道继续在黏土层与砂砾石层界面发展到x=42 m处。

图8 管涌通道发展过程中x方向水力梯度等值线

图9 管涌通道发展过程中不同区域x方向水力梯度

图9(b)为细砂夹层底部与砂砾石层界面x方向水力梯度分布图，同样地，水力梯度变化规律与细砂夹层上部类似；然而，在22 m上游水头下，细砂夹层底部的水力梯度最大值达到了0.54，正好等于由式(9)计算出的细砂发生水平向流土破坏的临界水力梯度，因此细砂夹层在x=82 m处发生水平向的流土破坏，假设发生流土破坏的厚度为细砂夹层厚度的一半，同时，细砂夹层底部砂砾石中的管涌通道继续发展到x=48 m。

当细砂夹层底部发生流土破坏后，上部的细砂和砂砾石发生沉降，使得黏土层与砂砾石层之间形成间隙，设间隙处的渗透系数与管涌口一致，取0.5 m/s，间隙的产生使得黏土层与砂砾石层界面变得相对薄弱，更易于侵蚀的发生。细砂夹层上部和下部管涌通道分别发展到x=42 m和x=48 m，且细砂夹层发生流土破坏后，保持上游水头22 m不变，得到间隙产生后x方向水力梯度等值线如图10所示，可见最大水力梯度发生在流土破坏后的细砂夹层前端。

图10 间隙产生后x方向水力梯度等值线

由图11可以看出，细砂夹层上部和下部砂砾石层中的管涌通道继续发展最终贯通。在本文的数值模拟中，假定相同粒径的颗粒被侵蚀，且管涌通道形成后，其渗透系数不变，由于侵蚀的颗粒粒径较小，通道的渗透系数相对较小，管涌通道贯通后还存在较大的水力梯度，在水流的作用下，通道中的土体继续侵蚀，管涌通道逐渐扩宽。

图11 间隙产生后不同区域x方向水力梯度

图12为间隙前端D截面y方向水力梯度，可以看出，间隙下部细砂夹层处y方向的水力梯度最大，达到了1.15，大于细砂夹层发生竖向流土破坏的临界水力梯度1.0，因此，细砂夹层在竖向发生流土破坏，黏土层与砂砾石层之间的间隙逐渐向深部和上游发展。

图12 间隙前端D截面y方向水力梯度

因而，在多层堤基管涌发展过程中，细砂夹层的存在将土体分为了多个区域，每个区域破坏发生的条件不同，在细砂夹层埋深为1 m的多层堤基结构中，通道发展主要经历了细砂夹层上部和下部砂砾石层的水平向管涌破坏、细砂夹层水平向流土破坏和竖向流土破坏等破坏形式，是深层管涌破坏与浅层管涌破坏同时作用，产生的管涌通道深度很大，对大坝安全产生重大影响。

3 结 论

a.多层堤基结构存在渗透系数较小的细砂夹层，对于细砂夹层埋深为1 m的结构，管涌发生的标志为管涌口下端细砂夹层冲破，发生竖向流土破坏，颗粒大量流失，堤基内部渗流场重新分布。

b.管涌发展过程中，堤基内部土颗粒不断流失，局部土体渗透性随之发生变化，渗透性的变化影响渗流场的分布以及堤基内各区域的水力梯度，继而又会影响土颗粒的流失，因而土颗粒流失和土体渗流是一个相互影响、相互耦合的过程。

c.管涌发展过程中存在多种类型的破坏形式，主要经历了细砂夹层上部和下部砂砾石层的水平向管涌破坏、细砂夹层水平向流土破坏和竖向流土破坏等破坏形式，是深层管涌破坏与浅层管涌破坏同时作用，产生的管涌通道深度很大，对大坝安全产生重大影响。