深厚覆盖层是指堆积于河谷中，厚度大于30 m的第四纪松散沉积物[1]，其内部常存在块石架空层，建坝后易形成强渗透作用，为坝基发生内部侵蚀提供条件．深厚覆盖层、坝基或堤身内部的渗透变形问题实质是侵蚀型内管涌问题，在侵蚀型内管涌发生发展过程中细颗粒始终在孔隙中移动，因此从颗粒细观层面研究内部侵蚀过程是难点和重点．针对侵蚀型内管涌带来的严重后果，学者对土体侵蚀型内管涌或内部失稳的影响因素进行了大量研究．文献[2]探究了土壤抗侵蚀能力与临界反滤层设计之间的相互关系．此外，研究人员又发现其他影响因素：内部因素，如GSD曲线和裂缝通道[3]；外部因素，如围压[4]和水力梯度．在物理测量方面，文献[5]采用CT扫描进行侵蚀型内管涌试验观测，但CT扫描费用昂贵，于是文献[6]通过采用显微摄像可视化跟踪与数字图像处理技术，成功地对砂土管涌的整个过程进行跟踪记录．得益于计算机技术的发展，学者开始采用数值模拟研究内部侵蚀机理．内部侵蚀宏观特性常用有限元法(FEM)[7]或有限差分法[8]模拟；而颗粒运移等细观特性常用离散单元法开展相关研究，该方法可以单独考虑每个颗粒的运动[9]，例如文献[10]应用PFC3D对不同层间系数下砂土管涌的基料-滤层系统进行离散元数值模拟，从细观角度揭示了管涌发展过程中颗粒的运动特性和反滤层防治机理．最近，其他数值方法也被引入来解决内部侵蚀特殊问题，如与DEM相结合的计算流体动力学[11](CFD)、孔隙尺度有限体积法[12]等．侵蚀型内管涌建模实质是流固耦合问题，当前应用DEM-CFD法进行流固耦合的案例较多．文献[13]基于泥石流启动过程中砂土的非饱和特性，应用PFC3D模拟流滑型泥石流的破坏形态和启动过程，同时与原有未考虑非饱和特性的数值模拟方法比对，证明了该方法的可靠性．文献[14]基于DEM-CFD方法模拟隧道上方土体内侵蚀的发生发展过程，并从细观层面呈现土颗粒团力链网络随侵蚀发展的动态演化过程，其变化特征符合实际规律．本课题组前期曾提出结合PLIF与MRI技术重构室内模型的方法，并验证了采用这一方法的可行性[15]．本研究基于颗粒-孔隙尺度流固耦合(DEM-CFD)分析模型开展不同水力梯度无黏性土的渗流仿真试验，从细观角度探讨不同水力梯度下粒间应力传递路径的演变趋势，综合细料流失率、喉径分布变化情况等，揭示细料受力与侵蚀型内管涌发生的内在关联，探究地层中侵蚀型内管涌时空演化过程．1 水力梯度对侵蚀型内管涌发生发展影响性试验本研究开展侵蚀型内管涌数值模拟，主要聚焦颗粒间接触力分布及喉管等细观尺度影响因素，因此采用DEM-CFD法开展研究，通过PFC软件自带的CFD模块实现流固耦合．1.1　模型生成及工况介绍考虑避免尺寸效应影响，并允许出现细料沉降、堵塞等效应，建立了200 mm×100 mm×100 mm的数值模型槽，上游边界与Y=0平面重合，流体网格单元尺寸为10 mm×10 mm×10 mm，共划分为2 000个单元，渗流方向沿y轴正方向，如图1所示，图中：P为水头；P0为初始水头．10.13245/j.hust.230082.F001图1生成模型及测量球布置示意图采用填充扩展法生成土颗粒，粒径放大次数为5次，达到平衡的标准是土样系统最大不平衡力小于1×10-5 N，此时可认为颗粒处于自然堆积状态，以模拟覆盖层土．模型共计生成7.073 5×104个颗粒，其中：粗颗粒2 829个，粒径为8~10 mm，构成骨架颗粒；细颗粒6.790 6×104个，粒径为1~2 mm，填充在粗颗粒孔隙中，细料质量分数为10%，孔隙率为0.35．建立5个测量球监测和记录模型的孔隙率变化情况，测量球的球心位于模型槽的中轴上，间距为5 mm，从左到右依次编号1~5，如图1(c)所示．颗粒间的接触模型选用线性接触，不考虑颗粒间黏结强度，颗粒与墙体的法向、切向刚度均设置为1×108 N/m，粗、细骨料内摩擦系数分别为1.25和0.70，墙体摩擦系数为0.5，流体密度与动力黏滞系数分别为1×103 kg/m3，1.01×10-3 Pa∙s，重力场默认9.8 m/s2，DEM和CFD计算时步均设置为1×10-8 s．为探究水力梯度对于内部侵蚀发展的影响，本研究设置了5组不同工况，如表1所示．10.13245/j.hust.230082.T001表1A组工况基本参数参数A1A2A3A4A5水力梯度0.500.751.001.251.50上覆荷载/MPa0.50.50.50.50.51.2　试验结果及分析各工况颗粒迁移过程相似，仅对A1，A3，A5工况颗粒最终迁移结果进行展示，红色框选区域用于统计出流口附近局部范围的颗粒迁移量，用于表征颗粒的逃逸能力，如图2所示．10.13245/j.hust.230082.F002图2不同水力梯度细颗粒迁移最终结果1.2.1　颗粒流失量变化分析图3统计了细颗粒流失数占颗粒总数比例的变化规律，图中：e为细颗粒流失数占总颗粒数比例；b为时步．可以看到：在系统达到稳定之前，细颗粒流失保持较高速率，表明侵蚀型内管涌的发生时间极短，且无明显破坏现象，发展极为隐蔽；随着水力梯度的增大，细颗粒流失量达到稳定时间逐渐提前，侵蚀型内管涌的破坏速度加快．结合图2红框区域可知：水力梯度越小，颗粒流失总量越大，其原因在于当水力梯度较小时，颗粒流速也较小，颗粒能够有序通过较为狭窄的通道，而当水力梯度较大时，细颗粒流速也较大，有少量颗粒同时堆积于一个狭窄通道中，造成部分流失通道堵塞，导致流失总量减少．10.13245/j.hust.230082.F003图3细颗粒流失数占比变化规律1.2.2　颗粒运移区域分析模拟结束后，将细颗粒运移距离分段，统计不同水力梯度下各运移距离范围的颗粒数，如表2所示，表中l为运移距离．水力梯度增大，颗粒运移速度加快，故而从表2中可以发现：水力梯度越大，运移距离大于100 mm的颗粒数量越多，小于100 mm的颗粒数量则越少．10.13245/j.hust.230082.T002表2不同水力梯度下各组运移距离的颗粒数量水力梯度l/m(0.01，0.05](0.05，0.10](0.10，0.15]0.150.127.5102.6050.4810.3130.327.4723.1530.9170.4140.527.0643.3211.1520.6480.726.1723.0191.5430.8391.025.7252.8631.5881.006103统计1号、3号和5号测量球中细颗粒数随时步的变化情况，如图4所示．从图中可以看出：1号和3号测量球中的细颗粒含量呈下降趋势，5号测量球则呈先上升再下降的趋势，试样细颗粒的流失主要体现在1号测量球中，占比约为65.5%，3号球的颗粒损失量相对较少，占比约为28.5%，5号球中的细颗粒几乎没有损失．10.13245/j.hust.230082.F004图4测量球中细颗粒数量随时步变化规律随时间推移，下游出流口附近应出现细颗粒堆积现象，而5号测量球中细颗粒数几乎不变，与实际不符．本研究将5号测量球沿渗流方向分为1~5层，统计5种工况下各层细颗粒数均值在试验前后的变化量，如表3所示．可以发现：试验结束后1-3层中细颗粒数较试验前均明显增加，第4层细颗粒数变化较小，而第5层细颗粒数明显减少，表明前4层中出现了细颗粒堆积情况，第5层的颗粒发生流失．因此前4层颗粒变化情况是符合实际的，而第5层下游由于没有设置墙体阻碍颗粒运移，出流口大小为整个系统的xoz面面积，无法滞留大量细颗粒，因此没有出现颗粒堆积现象．10.13245/j.hust.230082.T003表3试验前后5号测量球各渗流层细颗粒数变化渗流层12345细颗粒数变化量415420384202-133 32 水力梯度对力链网络影响性分析2.1　对力链数量影响性分析试验开始，一固定上覆荷载施加到试样上，导致颗粒间接触力陡增，尤其是粗颗粒间接触由于承担主要应力传递功能，接触力增量较大，接触总数变化规律如图5(a)所示．在渗流力作用下，粗细颗粒间产生大量接触力从而导致接触总数不断增加，水力梯度越大，新增接触力越多，直到系统稳定，粒间应力调整结束，接触总数不再大量增加，与细颗粒流失量变化规律相似．图5(b)中强力链数变化规律与接触总数类似，强力链数先快速增加后逐步趋于稳定，但增幅却远小于接触力总数．10.13245/j.hust.230082.F005图5强力链相关参数随时步变化规律强接触占接触总数比例的变化规律如图5(c)所示．强接触占接触总数的比例随时步呈下降趋势，表明游离态细颗粒形成的接触大多是弱接触，对主要应力的承载作用很小，只有少部分颗粒形成强力链参与应力传递，这些颗粒构成骨架颗粒支撑整个试样的结构稳定性，其他颗粒围绕在强力链周围支撑强力链不发生断裂．当细颗粒流失量逐渐增加时，细颗粒在力链中的参与度降低，部分参与强力链构成的细颗粒流失或转化为弱接触力，对强力链的支持作用减弱，强力链发生断裂的可能性增加，导致结构破坏的风险提高，这与侵蚀型内管涌的发展机理符合．在系统稳定之后，强接触占接触总数的比例不再发生大的波动，此时水力梯度越大，渗流力就越大，为抵抗渗流力而参与到应力传递的颗粒数越多，从而导致强力链数量增加．由于细颗粒间产生大量弱接触，使得平均接触力降低，而系统认为大于平均接触力即为强接触力，因此成为强接触的门槛降低，导致部分弱接触被判定为强接触．通过比对A1组与A5组在试验结束后强接触的变化情况发现：当以A1组的平均接触力标准判定A5组时，A5组强接触数量小于A1组，但接触力整体高于A1组，表明水力梯度越大，接触力链在断裂重组后，力链网络产生更为集中的应力传递路径以抵抗更大的渗流力．弱力链的支持作用是有限的，强力链承担的应力份额增大，系统失稳破坏的风险提高，因此水力梯度的增大会加剧侵蚀型内管涌的发生发展．2.2　对接触力分布方向影响性分析随着细颗粒持续流失，试样系统内部颗粒位置重新调整，影响力链网络结构．选择接触力与xoy，xoz及yoz平面的夹角展示接触力方向，限于篇幅，仅对部分工况初始、末尾时刻的接触力各方向分布情况进行介绍，如图6所示，图中j为接触力数量．10.13245/j.hust.230082.F006图6部分工况不同时刻各方向接触力数量统计图对比不同工况接触力方向分布情况，发现：在水力梯度较小处，接触力与xoy平面的夹角大于45°的部分更多，且集中于70°以上，此时上覆荷载对接触力方向的影响更明显，接触力主要传递上覆荷载引起的应力；与xoz平面的夹角大于45°的部分更多，与yoz平面的夹角则是小于45°的部分更多，说明接触力的分布更接近y轴，这是因为水力梯度的作用使得接触力向y轴聚集以抵抗由水力梯度产生的渗流力．随着水力梯度的增大，虽然接触力与xoy平面的夹角大于45°的部分仍相对更多，但从A3组接触方向分布图中能够看出，接触力方向主要集中在50°~60°，与之前的主要集中在70°以上相比，接触力向45°方向聚集．这表明水力梯度对接触力方向的影响逐渐明显，力链网络正在进行结构调整以平衡水力梯度与上覆荷载的影响，同时，接触力与这两个平面的夹角也正朝向y轴收拢，这也印证了水力梯度对接触力分布方向的影响逐渐明显．在水力梯度增大到一定程度后，接触力与xoy平面的夹角分布逐渐变扁，表明夹角小于45°的接触逐渐增多，最终数量反超大于45°的部分．同时，与xoz平面的夹角大于45°的接触变多，与yoz平面的夹角小于45°的接触也变多，表明接触力向y轴集中的程度增加，这是因为随着水力梯度的增大，渗流力也同步增大，试样原有的力链不足以承担渗流力，因此为了传递增大部分的渗流力，更多细颗粒参与进来生成新的接触力和力链，从而导致力链网络发生调整．3 细颗粒流失概率理论适用性分析3.1　水力梯度对颗粒流失概率影响分析Dallo[16]对许多渗流模型进行了测试，最终提出1个导致向前运动喉管F和0.33个侧向运动喉管0.33S的模型(表示为1F+0.33S)，因此细颗粒通过喉管，发生流失的概率为P=p+0.33(1-p)p2，(1)式中p为CSD曲线中直径大于d1的部分的喉管累积体积分数．采用Dallo的理论，根据文献[17]提出的颗粒-空隙弦方法进行分层，确定试样共分为28层，每层厚度为7.2 mm．采用可视化软件Avizo获取试样喉径分布，如图7所示．细颗粒粒径在1~2 mm范围，取间隔为0.1 mm将细颗粒粒径划分为10个组，分别计算这10个粒径对应的P值，如表4所示，表中d为粒径．10.13245/j.hust.230082.F007图7CSD曲线图10.13245/j.hust.230082.T004表4细颗粒粒径及对应的P值d/mm1.11.21.31.41.51.61.71.81.92.0P0.7800.7570.7140.6500.6130.5790.4910.4290.3970.373通过常规方法获取A组工况流失细颗粒的GSD曲线，与Dallo理论计算得出的GSD曲线对比情况如图8所示，图中u为小于某粒径的细颗粒质量分数．从图中可以看出：随着水力梯度的增大，GSD曲线逐渐上凸，流失细颗粒的GSD曲线无法用Dallo的概率理论很好地表示出来．10.13245/j.hust.230082.F008图8A组试验流失细颗粒GSD曲线在Dallo的理论研究中，细颗粒的流失是颗粒重力作用的结果，而水力梯度引起的渗流力与颗粒重力同属于体积力，本研究预测水力梯度对于细颗粒的流失概率的影响与颗粒重力相似．本研究尝试从颗粒细观承受的体积力这一角度，通过比较颗粒所受重力与渗流力之间的关系，探究水力梯度对细颗粒流失概率的影响．1F+0.33S模型决定了P的计算式，因此当采用渗流力对该计算式进行修正时，同样可以从水平向流失概率与竖向流失概率进行．为确定重力与水力梯度共同作用下不同流失方向在流失概率中的权重，根据颗粒流失总量统计了5组工况下各粒径组的流失概率P，如表5所示．10.13245/j.hust.230082.T005表5修正理论与模拟试验各粒径组P对比统计表d/mmA1A2A3A4A5理论试验理论试验理论试验理论试验理论试验1.20.6750.6790.6820.6870.7230.7270.7720.7660.7870.7821.40.6440.6480.6680.6640.6810.6880.6580.6660.6600.6651.60.6580.6590.6310.6340.5830.5870.5410.5410.5210.5271.80.6530.6500.6080.6130.5310.5360.5040.5120.3900.3942.00.6030.6060.5950.5910.4960.49504110.4130.3010.309选取两个权重系数α和β，其中：α为水平向体积力对细颗粒流失影响所占权重；β为竖向体积力对细颗粒流失影响所占权重．根据A组工况不同粒径组的细颗粒流失概率理论值，确定α和β．然而，由于α和β随着水力梯度的增大发生变化，直接分析α和β与水力梯度之间的关系难度较大，因此补充设置5组(B1~B5)不同工况，通过开展不同渗流方向内侵蚀数值模拟确定α和β在特定条件下的值，并预测其与水力梯度间的关系．数值模型由单一粒径的粗细颗粒组成，粗颗粒直径为6 mm，细颗粒直径为1 mm．各工况下的细颗粒流失概率指标如表6所示．10.13245/j.hust.230082.T006表6B组试验颗粒流失概率指标统计表分组渗流力/(kN∙m-3)土的容重/(kN∙m-3)理论P值试验P值αβB118.0(竖直向下)-18.00.630.840.021.36B20.0-18.00.630.650.000.97B318.0(竖直向上)-18.00.630.020.340.00B418.0(水平向右)-18.00.630.720.331.03B536.0(水平向右)-18.00.630.871.150.16从表6中可以看出：在重力条件固定、渗流力大小相同的情况下，当渗流力与土的重力方向一致时，细颗粒流失概率最大，当渗流力与土的重力方向相反时，细颗粒流失概率最小．因此，本研究认为细颗粒流失是由水平向体积力与竖向体积力共同作用的结果，而每一个单独方向的体积力遵循Dallo的理论，于是将Dallo的理论计算式修正为P=β[p+0.33(1-p)p2]+0.33αp．(2)细颗粒的流失方向取决于α和β的大小关系．根据式(2)可以确定B组各工况下α和β的特定值，见表6．当渗流力与重力方向一致时，β值最大，α值最小；当渗流力与重力方向相反时，β值最小；当渗流力与重力方向垂直时，α值最大．由此可以得出：β值受竖向力的影响较为明显，而α值则受水平力的影响更明显．同时可以看出：当细颗粒流失方向确定之后，与流失方向相同的体积力作为主导因素占有更多权重，而另一个方向的同等大小的体积力大约占到了主导体积力权重的1/3．因此，本文推测权重的计算公式为α=(GH/GV)/3；(3)β=1+(GH/GV)/3，(4)式中：GH为水平向体积力；GV为竖向体积力．根据B5组可推得：当水平向体积力大于竖向体积力时，该计算公式不再适用，因为颗粒的流失方向发生改变．在将α和β的值对换之后获得了较为一致的结果，所以对于式(2)进行修正，有P=β1[p+0.33(1-p)p2]+0.33α1p；(5)α1=(GH/GV)/3；(6)β1=1+(GH/GV)/3，(7)式中：α1为主导体积力对细粒流失影响权重；β1为非主导体积力影响权重．式(5)~(7)完成了对Dallo理论公式的修正．通过式(5)对A组模拟的细颗粒流失概率进行预测，预测结果如表5理论组所示，可以看出预测值与实际值较为一致．3.2　改进理论适用性验证试验为验证修正后的Dallo细颗粒流失概率理论，对前人所做的部分试验进行整合．文献[18]利用PFC构建了一个具有图9(a)所示GSD曲线的数值土壤样本，采用Delaunay镶嵌法获得CSD曲线，通过模拟渗流区内基土颗粒的渗流过程，给出了可在12，36，60 mm厚度下通过渗流区的基土颗粒GSD曲线，图中D为渗流区厚度．文献[19]进行了物理模型试验测试，获得可以透过3.3，5.0，14.5 cm渗流厚度的基土颗粒GSD曲线，如图9(b)所示，这些结果可用于验证新的概率理论．10.13245/j.hust.230082.F009图9改进理论与既有研究流失细颗粒GSD曲线比较土层间距由前面提到的数学方法得出，当D=12，36，60 mm时，文献[18]土层间距分别为2.02，6.05，10.08 mm，当D=3.3，5.0，14.5 cm时，文献[19]土层间距均为4.19 mm．统计试样每个粒径组的细颗粒数量，获得可通过任意层的每个粒径组的粒子数量．通过计算，可以得到流失细颗粒的GSD曲线，并与理论计算结果进行比较．由图9可知，通过新的概率理论预测流失细颗粒的GSD曲线与既有研究的结果具有较好的一致性，表明该概率理论是较为准确的．4 结论a．水力梯度的增加会对土体内部抗渗稳定性产生不利影响．随着水力梯度的增加，细颗粒流失速度加快，侵蚀型内管涌发展速度也随之增加；同时，试样内部产生大量新接触，但基本以弱接触为主，试样平均接触力降低，导致成为强接触力的门槛降低，部分弱接触力被判定为强接触力．按照初始阈值判定试验结束后的接触力分布发现：能够承担应力传递的强接触数量减少且分布方向更加集中，表明力链网络产生更为集中的应力传递路径以抵抗增大的渗流力，在细颗粒大量流失对强接触的支撑作用减弱后，强力链断裂的概率增大，从而导致土样发生内部侵蚀的风险增大．b．Dallo提出的细颗粒流失概率理论具有局限性，忽略了其他体积力对细颗粒流失概率的影响，不同方向体积力对细颗粒流失概率具有一定的影响权重，权重较大者决定细颗粒流失的方向；本研究沿用Dallo的研究思路确定了不同方向体积力影响权重计算公式，对原有概率理论进行修正，并根据前人的研究数据对修正的细颗粒流失概率理论进行验证，结果表明修正的理论能够较好地预测流失细颗粒的粒度分布情况．