《1 前言》

1 前言

接触冲刷是指渗透水流沿着两种不同粒径组成的土层层面流动而带走细颗粒的现象。上下两地层的颗粒直径悬殊越大越易发生接触冲刷。在自然界中沿两种介质界面，如坝体与坝基、土坝与涵管等接触面流动而造成的冲刷，都属于此类破坏，一般认为接触冲刷在砂砾石层与粘土层的交界面处最易发生。笔者将从分析接触冲刷的物理过程及其地层渗透性概化模型出发，用分时段的稳定井流模型模拟接触冲刷发展的非稳定过程，探讨接触冲刷的发展规律。

《2 接触冲刷发生后地层渗透性概化模型》

2 接触冲刷发生后地层渗透性概化模型

接触冲刷的实质是渗透水流往往在两种不同粒径地层的触面附近水流最集中，因而最易将其中的细颗粒随水流带走。接触面附近水流最易集中的原因：a.沉积成因，两地层之间若存在沉积间断，因古风化而使接触面粗糙，空隙度增大，渗透性提高；b.固结成因，上下两地层固结的速度和程度有所不同，因而易在接触面形成微小的缝隙，大大提高接触面附近的渗透性能；c.结构物与土的刚柔不一而在它们之间产生隙缝；d.上覆为粘性土而下覆为砂或砂砾石层的双层结构地基，按正常的渗流计算，接触面附近的水力坡降最大，因而最易将细颗粒带走。接触面附近细颗粒的不断流走，进一步增大接触面附近的空隙度或缝隙，甚至将细颗粒层掏空，发育至河岸时形成集中渗漏通道，危及堤坝安全。

砂砾石层与粘土层的交界面处接触冲刷最易发生，以下讨论接触面附近细砂从砂砾石层中因冲刷而流失的接触冲刷发展过程。假设：

1）砂砾石层呈均质各向同性，在接触冲刷的发展过程中，砂砾石层中的细砂呈层状向下逐渐剥离而随水流流失，从而使砂砾石层分为上部的砾石层和下部的砂砾石层；

2）接触面呈水平状，将接触面因各种原因而增大的空隙或存在的缝隙概化为裂隙，其渗透性由裂隙的渗透性代替；

3）细砂和砾石均呈理想的等粒球状，砾石形成的空隙全部被细砂充填满，细砂流失后的砾石层渗透性由球体呈紧密堆积时的理论公式确定。

若砂砾石层的渗透系数为Ko，砾石层的渗透系数K1和裂隙的渗透系数K2为[1]

\(K_{1}=\frac{\beta n^{3}}{\lambda(1-n)^{2}} d^{2} \frac{\gamma_{\mathrm{w}}}{\mu}, \quad K_{2}=\frac{g b^{2}}{12 \nu}\)。        （1）

式中β为球体系数（圆球取β＝π／6），n为孔隙度，λ为决定于邻近颗粒的影响系数，d为颗粒直径，Yw为水的容重，μ为水的粘滞性，b为裂隙宽度，g为重力加速度，v为水的动力粘滞系数。

地层的渗透性概化模型如图1，任意时刻地层的渗透系数K可表述为

\(K=\left[K_{0}(M-\Delta H-b)+K_{1} \Delta H+K_{2} b\right] / M\),           (2)

式中符号见图1，其中M为渗透层厚度，砾石层厚度ΔH是关于时间以及渗透流速的函数。

《图1》

图1 接触冲刷前后的地层渗透性变化模型

Fig.1 The model of permeability change after contact scouring between the two layers

设在Δt时间内在接触面附近流出的细砂体积ΔV，将体积ΔV的细砂简化为在体积ΔxΔyΔH的砂砾石层中流出的细砂，即nΔxΔyΔH＝ΔV，得ΔH＝ΔV／nΔxΔy。n为砾石层的孔隙度，立方体排列时n＝0.4764，四面体排列时n＝0.2595，则

\(K=\left[K_{0}(M-\Delta V / n \Delta x \Delta y-b)+\right. \\\left.K_{1} \Delta V / n \Delta x \Delta y+K_{2} b\right] / M\)。           （3）

《3 计算模型》

3 计算模型

《3.1 直线透水边界附近井流计算理论》

3.1 直线透水边界附近井流计算理论

视出渗口为承压完整抽水井，涌水量Q即为抽水井流量；设渗透层厚度为M，原始渗透系数为K0，导水系数T＝K0M；视河岸为直线透水边界，井中心至直线透水边界的垂直距离为a。根据镜像法原理，可设想在直线透水边界的另一侧与边界的垂直距离为-a的位置上存在一流量为-Q的注水井，见图2。

《图2》

图2 直线补给半无界井流计算示意图

Fig.2 A calculation chart of linear recharge of half space by the well flow

堤坝附近渗流场内任一点P（x，y）的水头降深S（x，y）可由叠加原理得[2]

\(\begin{array}{l}S(x, y)=\frac{Q}{2 \pi K_{0} M} \ln \frac{R}{r_{1}}- \\\frac{Q}{2 \pi K_{0} M} \ln \frac{R}{r_{2}}=\frac{Q}{2 \pi T} \ln \frac{r_{2}}{r_{1}},\end{array}\)       (4)

式中\(  r_{1}=\left((a-x)^{2}+y^{2}\right)^{1 / 2}, r^{2}=\left((a+x)^{2}+\right.   \left.y^{2}\right)^{1 / 2} \) 。

任一点的水头  H  和降深  S  的关系为

\(H=H_{0}-S\) 。       （5）

势函数为

\(\varphi=K H\) 。               （6）

将式（4）代入式（5），再将式（5）代入式（6），得

\(\varphi=K\left(H_{0}-S\right)= K H_{0}-\frac{K Q}{2 \pi T} \ln \frac{\left.\left((a+x)^{2}+y^{2}\right)\right)^{1 / 2}}{\left((a-x)^{2}+y^{2}\right)^{1 / 2}} \)          (7)

令\(\frac{(a+x)^{2}+y^{2}}{(a-x)^{2}+y^{2}}=\exp \frac{4 \pi T\left(K H_{0}-\varphi\right)}{K Q}=c\)，

则直接由式（7）可得等势线方程为

\(\left(x+\frac{(1+c) a}{1-c}\right)^{2}+y^{2}=\frac{4 c a^{2}}{(1-c)^{2} }\)          (8)

令φ为不同的值，可得到一系列的等势线。

《3.2 单元剖分及计算方法步骤》

3.2 单元剖分及计算方法步骤

单元的剖分按流网进行，即流网中的每个网格就是计算和剖分单元。根据前面的假设，渗透地层在发生接触冲刷前及初期是均质各向同性体，流线和等势线应互相垂直（图3）。另外，流网的绘制按相邻两等势线之间的水头差相等和相邻两流线之间的流函数差值相等的原则进行，则通过每个网格的流量相等。以相邻两流线之间的所有单元为一个单元组（图4），设有m个单元组，每个单元组中沿径向有n个单元，则总共有nm个单元，每个单元中的初始流量为管涌口流量Q的1／m。 

《图3》

图3 发生初期地下水流场分布图

Fig.3 Groundwater flow distribution at groundwater gushing beginning

《图4》

图4 第i单元组的单元剖分图

Fig.4 The divided unit ith

接触冲刷的发展过程是一个动态变化过程，随着时间的推移，包括渗透地层各单元接触面附近的细砂含量和渗透性以及动水压力等均将改变。为准确模拟接触冲刷的发展过程以了解接触冲刷发展规律，计算中将接触冲刷发展的非稳定动态过程分解为一系列时间间隔Δt的稳态过程，只要Δt的取值足够小，并假设在Δt，时间内各单元的渗透性K，近似为常数，那么这样的模拟过程完全符合接触冲刷本身的实际发展过程，其模拟结果完全是可靠的。模拟计算的步骤如下：

1）根据江水位和出渗口的地面高程和直线透水边界的井流公式，估算初期的涌水量Q，并平均分配到各单元组：

\(\begin{aligned} Q=& \frac{2 \pi K\left(H_{0}-H\right) M}{\ln \left(\frac{2 a}{r_{w}}\right)}, q_{i}=Q / m \\\end{aligned}\),      (9)

式中rw为出渗口半径，M为含水层厚度，H0为江水位，H为出渗口水位，qi为第i单元组的流量。

管涌发生初期通过镜像法求出任一点P（x，y）的水力梯度为

\(\begin{aligned}J(x, y)=\frac{Q}{\pi T} \frac{a^{2}}{r_{1}^{2} r_{2}^{2}}=\\& \frac{Q}{\pi K_{0} M} \frac{a^{2}}{\left[(a-x)^{2}+y^{2}\right]\left[(a+x)^{2}+y^{2}\right]}\end{aligned}\) 。            (10)
在主渗漏通道（y＝0）上的水力坡降为

\(J_{x}=\frac{Q}{\pi K_{0} M} \frac{a^{2}}{\left(a^{2}-x^{2}\right)^{2}}\)  ，          (11)

作出y＝0方向上水力梯度J，随x的变化曲线。

2）判定各单元内的渗透流速是否达到或超过颗粒的起动流速。由于在Δt，时段内单元的渗透系数为常数，故通过各单元中的渗透系数与水力梯度可以求出该单元的渗透流速，如果流速达到了细颗粒的起动速度，则将有细颗粒被带出该单元。

细砂的起动速度v可按下式计算：

\(B v^{2}+M v-N=0 \)      (12)

式中dsk为等效粒径，y，为土粒容重，n为水的容重，a为土体颗粒的形状系数，ds为土粒直径，

\(\begin{array}{c}M=\frac{6 A v(1-n)}{d_{\mathrm{s} k}}, N=\frac{1}{6}\left(\frac{\gamma_{\mathrm{s}}}{\gamma}-1\right) \alpha n^{3} g d_{\mathrm{s} } \\\end{array}\)

作为层流运动（一般流速很缓慢）时，可令B＝0，于是得到v＝N／M。

3）计算Δtk时间内各单元中细砂量的改变。由于单元组的划分是以流线为界的，因此细砂颗粒的移动只能在单元组的上下游相邻单元之间发生。设各单元向下游单元流出的砂量与该单元的流量成正比，即Q砂＝c1Q水，式中c1为常数。该单元同时又接受了上游单元的来砂量，因此，在Δtk时间内第ij单元中

实际砂量的改变量为

\(\begin{array}{c}\Delta V_{i j}=V_{i j}-V_{i, j-1}\end{array}\) ，     (13)

式中 Vij为第ij单元流出的细砂量体积，Vi,j-1为由第（i，j-1）单元流入第ij单元的细砂体积。

4）计算各单元在Δt末的渗透系数。据式

1）得

\(\begin{array}{c}K_{i j}=\left[K_{0}\left(M-\frac{\Delta V_{i j} \Delta t}{n \Delta S_{i j} \Delta L_{i j}}-b\right)+\right. \\\left.K_{1} \frac{\Delta V_{i j} \Delta t}{n \Delta S_{i j} \Delta L_{i j}}+K_{2} b\right] / M \\\end{array}\)，              （14）

5）根据涌水量Q估算水头损失Δh。公式[3]为Δh＝c2Q2，其中c2为井损常数，可通过实验确定。

6）计算新的出渗口水头。在新的时间单元Δtk＋1中，由于各单元组的平均渗透性和出渗口水位的变化，涌水量已经发生了改变，可通过各单元组的平均渗透性\(\bar{K}_{i}\)和新的管涌口水位计算新的流量。

考虑径向第i单元组，假设流过第ij单元的流量为

\(q_{i}=K_{i j} \frac{\Delta H_{i j} \Delta L_{i j}}{\Delta S_{i j}} M\),

式中ΔHij为第ij单元的水头差，ΔLij为第ij单元垂直于流线方向的平均长度，ΔSij为第ij单元的平均渗径长度，从而可求得流过该单元的水头变化

\(\Delta H_{i j}=\frac{q_{i} \Delta S_{i j}}{K_{i j} \Delta L_{i j} M}\)。

由于沿着流线方向，流过各单元的流量相等，而总水头变化值（ΔHj）等于各个单元水头变化值之和，即\( \Delta H_{i}=\sum \Delta H_{i j}\)。

考虑井损后的水头为\( \Delta H_{i}=\sum \Delta H_{i j}+\Delta h \)。

7）计算新的涌水流量。第i单元组的平均渗透系数为K；和新的流量为

\(\begin{array}{c}\bar{K}_{i}=\sum_{j=1}^{n} \Delta S_{i j} / \sum_{j=1}^{n} \frac{\Delta S_{i j}}{K_{i j}} \\q_{i}=\bar{K}_{i} M \Delta H \sum_{j=1}^{n} \Delta L_{i j} / n \sum_{j=1}^{n} \Delta S_{i j}\end{array}\)

式中n为单元组中的单元剖分个数。

8）重复以上1-7的步骤，进行下一个时间段的计算，直到管涌口的水位或流量基本达到稳定止。

《4 算例分析》

4 算例分析

假设地基土层上覆为粘土层，粘土层水平分布，下伏为砂砾石层，该砂砾石层裂隙宽度2cm，其渗透系数140.94 cm/s，砂砾石层的渗透系数0.078  cm/s，大颗粒直径25mm，渗透系数36 cm/s，水头差5m，地层厚度1m，接触冲刷出渗口距离河岸10m。接触冲刷出渗口的半径假设为0.5m。

《4.1 接触冲刷发展过程中水头和渗透性变化》

4.1 接触冲刷发展过程中水头和渗透性变化

通过模型计算，得到接触冲刷发生前后地层中水头的变化，接触冲刷发生初期地层中的渗透系数相等。地层水头变化参见图5a，近似为偏心的圆形分布，靠近堤坝一侧的水力梯度较大。随着接触冲刷的发展，接触面附近的细砂被水流带出，地层中的渗透系数发生变化，出渗口附近的渗透系数逐渐变大并逐步向江堤方向发展，出渗口附近的流场发生畸变，水头分布如图5b所示，出渗口附近的水位升高，靠近堤坝一侧的渗透系数明显增大。当接触冲刷进一步发生直至水量基本稳定后，水头分布参见图5c，出渗口的水位进一步升高，靠近堤坝一侧的接触面附近细砂都基本上被带出，最终形成了集中渗漏通道。接触冲刷发生前后的地层渗透性见图6。

《图5》

图5 接触冲刷发生后地层中的水头分布

Fig.5 The groundwater head distribution after the contact scoring

模型研究发现，接触冲刷从发生到形成集中渗

《图6》

图6 接触冲刷发生后的地层渗透性

Fig.6 The change of permeability distribution when contact scouring emergence

漏通道经历的时间很短，计算得到的经历时间一般仅为几分钟到十几分钟，由于计算中没有考虑颗粒间的摩阻力、水流与颗粒间的摩擦力、孔隙通道形状的复杂性等原因，实际情况可能比计算的经历时间要长一些，但一般形成集中渗漏通道的时间也只会增加几倍，也就是说，一旦接触面附近的水力梯度超过临界值，就可能产生接触冲刷破坏，这个破坏过程经历的时间非常短，在洪水期应注意加强堤坝附近的巡视，特别是靠近堤坝段的巡视，发现隐患，应立即进行处理，否则一旦形成接触冲刷破坏，在几分钟到几十分钟内就可能在地层中形成通道，最终造成溃堤事故。

《4.2 接触冲刷发展过程中沿流线方向上水力坡降的变化》

4.2 接触冲刷发展过程中沿流线方向上水力坡降的变化

计算发现，沿各流线方向上水力坡降的变化以出渗口到河堤之间的垂直方向上为最大、最典型，如图7。初始时刻（t＝0）是一条双曲线，水力梯度J从河堤到出渗口呈增大趋势，接近口处骤然增大，并在出渗口附近的水力梯度最大（图7中t＝0线），接触冲刷发生后，一旦接触面附近的砂涌出，水力梯度随即降低（如图7中t＝Δt线）。随着时间的推移和细砂的进一步涌出，水力梯度的峰值向河堤方向移动并呈增大趋势。这说明接触冲刷的发展方向是溯流而上并加速发展，进一步说明了发生接触冲刷后及时堵漏抢险的重要性。

《图7 》

图7 接触冲刷发生后主渗漏通道上 J-x变化曲线

Fig.7 The J-x curve in the main pathway of the layer after contact scouring emergence

4.3 接触冲刷形成贯通性集中渗漏通道与出渗口位置及地层渗透性的关系

如果其他条件不变，只是在接触冲刷发生初期，砂砾石地层的渗透系数较小，设为0.0096cm／s，渗透系数如图8所示，虽然在破坏处渗透系数也发生改变，但不会形成渗漏通道。如果砂砾

《图8》

图8 接触冲刷渗透系数分布（砂砾石地层渗透系数0.0096 cm/s）

Fig.8 The permeability distribution of the layer (the gravel permeability of the layer is 0.009 6 cm/s)

石层的渗透系数仍为0.078cm／s，而把出渗口与河堤的距离改为15m，其他条件不变，其结果如图9。由渗透系数分布图可以看出，接触冲刷能否形成通道不仅与接触冲刷发生前夹层的渗透系数有关，而且还与堤坝后薄弱处至堤坝的距离有关。如果薄弱处离堤坝越远，则发生接触冲刷的可能性越小。

《图9》

图9 接触冲刷地层渗透系数分布（出渗口距河岸15m时）

Fig.9 The permeability coefficient distribution between layers after contact scouring

《5 工程实例分析》

5 工程实例分析

1998年，长江全流域发生大洪水，长江中游某段高水位（22.82m）持续时间长达2个月，长江干堤险象环生。溃口地段地处长江南岸I级阶地前缘，地势平坦，地面高程一般在17.66～18.79m。堤内距堤脚约20～30m沿堤方向展布有水塘，深度大于5m。1998年溃口抢险时，水塘大部分已用块石、炉渣等填平，目前仅在西南角残存一小水塘。堤外滩宽约30～90m，地面高程18.8～19.55m。溃口堤段内地表分布有人工填土层，其下为第四系全新统河流冲积层和冲湖积层，具二元结构，厚30～35m。上部为棕黄色、灰褐色粉质壤土、粉质粘土，厚度10～11m，其间多夹有粉细砂夹层；下部为灰色、灰黄色粉细砂、中粗砂砾、卵石层，厚20～25m。下伏基岩为下第三系新余群泥钙质砾岩及泥质粉砂岩。大堤地层示意图见图10。

《图10》

图10 大堤地层示意图

Fig.10 The layers structure of the dike

发生渗透破坏时江水位22.82m，渗透破坏处的高程18.8m，渗透破坏处距离大堤60m，砂砾石层厚度35m。计算发现，该段大堤溃堤的主要原因是：堤基土上部的10～11m粉质粘土层透水性较弱，为相对隔水层。下部砂层、卵砾石层为强透水层，是主要含水层。洪水期间，江水位高于上部粘性土底面高程时，地下水具有承压性质。假设不计水头损失，考虑最不利的情况，粘土层底部的水头差有14m，全部作用于粘土层底部水压力p＝1.0x14kPa，水塘处粘土层较薄弱，渗透破坏有可能在此发生，此处粘土层厚7m，水塘中有3m深的河水，粘土层底部的有效应力为σ'＝1.67x7＋1.0x3＝14.69kPa，粘土层底部的水压力小于有效应力。经过计算发现地下水的坡降不大可能超过临界水力坡降，因此该段的溃堤不是由于下部的砂卵石层引起的。从资料中可以看出，在上覆粘土层中有多层粉细砂夹层，夹层的渗透系数较大（k＝0.01 cm/s），经过计算发现，夹层中的水力坡降超过临界水力坡降，造成渗透破坏，由于细砂的起动速度较小，因此细砂在水流速度较小的情况下就被带走，在交界面处形成一小裂隙，随着小裂隙的进一步发展，在水流的冲刷作用下，粘土层的结构将发生破坏，从而造成溃堤。因此，该段大堤渗透破坏主要是由于粘土层中的夹层遭到接触冲刷而引起的破坏。

《6 结语》

6 结语

通过堤坝基薄层砂砾石层渗透性的概化和计算模型的建立和对接触冲刷发展过程的模拟，可以得到如下重要结论：

1）堤坝基接触冲刷的发生、发展是一个十分复杂的过程，笔者提出了针对均质各向同性地层渗透性的概化和计算模型，通过算例模拟了接触冲刷发展全过程，得到的接触冲刷发展规律具有一定理论意义。

2）接触面附近细砂颗粒的移动取决于渗透流速的大小是否达到颗粒的起动速度，而接触冲刷的发展方向是溯流而上并加速发展。

3）随着时间的推移和接触冲刷的发展，细砂从周围向出渗口至垂直河堤的主渗透方向上汇集并由出渗口向河堤方向扩展的同时，地层的渗透性增
大，水流集中可以形成集中渗漏通道，地层的渗透性和出渗口的位置对集中渗漏通道的形成有一定影响。