多级刚网结构拦阻碎屑流过程的数值模拟研究

苏嗣琦，樊赟赟

(东北大学 资源与土木工程学院， 辽宁 沈阳 110004)

碎屑流多发生在山区，是一种常见的地质灾害，其爆发时常具有高速、远程、冲击力强等特点。碎屑流爆发过程突然、历经时间短暂、来势凶猛、破坏力惊人[1]。碎屑流灾害常影响人类的正常生产生活并给经济财产带来巨大的损失，严重威胁着公路、铁路、输电线、石油管道等工程结构物的安全。因此，为应对碎屑流的冲击灾害，需要在受到威胁的基础设施附近采取防护措施，出于安全考虑防治工程应偏于保守[2]。由于布置不合理的防治工程结构损毁失效现象时有发生，为提高防护效率减少损失需要加强对碎屑流冲击作用、防护结构的拦阻作用、以及防护结构合理设计布局等进行深入研究。

碎屑流防治工程包括拦挡工程、排导工程及生物工程等方式[3]。其中拦挡工程是碎屑流灾害治理过程中一种最为常用的结构类型。试验表明拦挡结构能够大幅度的降低碎屑流流动的速度[4]，而固体颗粒间表现出的摩擦特性是碎屑流与防治结构相互作用中能量耗散的关键。近年来，研究人员Jiang等[5]开展了一系列的碎屑流冲击刚性挡墙的研究，通过改变坡体角度、坡面摩擦、坡面长度、颗粒属性等因素，揭示了碎屑流冲击过程中的冲击与堆积机制，发现了冲击力跟碎屑流的规模、冲击角度、厚度和速度有关，并初步探讨了冲击作用的量化估算方法。Liu等[6]利用压力计传感器进行了一系列水槽实验，探索颗粒对碎屑流冲击过程的影响。Zhou 等[7]提出了岩屑崩落运动会受地形条件、材料组成、滑床表面条件等因素的影响。王禄洲等[8]运用离散单元法给出了失稳块体不同运动阶段的运动形式，可为类似工程的防灾减灾提供理论依据。孙兴伟等[9]分析了西藏大古与街需水电站工程区的泥石流基本特征，通过室内物理模拟试验获取泥石流冲击力的试验结果，确定泥石流重度及坡度对冲击力的影响。眭静等[10]通过模型试验的设计，对冲击力的时间历程、分布范围、速度等试验结果进行分析，建立了碎屑流的冲击力学模型，推导了冲击力公式并验证了公式的合理性，基于试验对相关的影响参数进行了分析。模型实验对碎屑流运动研究已有众多科研成果，为工程建设及安全生产等方面做出了巨大贡献。但由于碎屑流运动过程中其内部机理十分复杂，环境复杂多变，使得目前对防护结构拦阻机理的研究仍不够透彻和全面。离散单元法能够描述颗粒材料的力学行为，适合于模拟防护结构拦阻碎屑流过程中的复杂相互作用。在室内试验研究的基础上，运用离散元法对拦阻过程进行数值模拟，可以实现对防护结构拦阻碎屑流过程深入的运动学与动力学特性分析。张睿骁等[11]以滑坡碎屑流为研究对象，通过对比模型试验和数值模拟结果，校正了三维离散元的计算参数，研究了挡板角度和高度对冲击力、最大水平运动距离的影响。李天话等[12]运用三维离散元法，对比试验探讨碎屑流不同级配对滑体堆积形态、运动速度及冲击力等动力学特征的影响。肖思友等[13]利用Hertz-Mindlin接触模型模拟，得到了碎屑流冲击柔性网时承力锚索的受力与变形。李坤等[14]通过物理模型实验及数值模拟的方法，以高速远程滑坡的分形粒径分布为依据，对不同粒径组成碎屑流的反粒序堆积特征及颗粒分选特征进行了研究。刘思杰等[15]采用离散元法对由不同粒径颗粒组成的岩屑滑坡碎屑流运动进行全程数值模拟，研究了碎屑粒径组成对岩石碎屑运动特性的影响。赵川等[16]采用离散元建立了三维碎屑流数值模型，对碎屑流运动特性进行深入探讨，研究碎屑流运移速度和颗粒接触数量关系，进一步定量分析高速远程碎屑流的动力学机制。近年来，数值方法因其成本低、应用方便、能够模拟碎屑流与拦挡结构相互作用等复杂问题而得到广泛应用。但由于离散元法对计算的要求较高，目前仅适用于小规模的碎屑流与防护结构相互作用的计算及力学机制的研究[17]。

虽然目前对碎屑流运动的研究已有很大进展，但因防护结构形式多样，应用环境复杂多变，对防护结构拦阻机理的研究仍存在不足。本文将利用离散元方法对多级刚网结构拦阻碎屑流过程进行数值模拟，重点研究不同条件下拦挡结构的受力特征，从而为进一步优化拦挡结构布局设计提供参考与建议。

1 数值模型及计算条件

本文采用基于离散单元方法的PFC3D软件进行数值模拟计算。结合本项研究的特点采用线接触和滑移模型，并采用clump单元模拟具有石英砂颗粒形态的碎屑颗粒，wall单元模拟边界，用linewall单元模拟刚性网状拦阻结构，接触模型如图1所示。目前除密度和摩擦系数等少数具有物理意义的细观参数外其他诸如接触刚度等细观参数与宏观参数没有确定的关系，这些细观参数主要以数值与试验结果吻合为标准对进行参数校准。

图1 接触模型示意图

离散单元间的接触力可分解为法向力和剪切力。

法向力和切向力增量可以分别表示为：

(1)

线弹性接触刚度模型是最为常用的接触模型，离散单元A和B间的法向接触刚度和切向接触刚度分别表示为：

(2)

式中:kn和ks分别表示两个离散单元的法向和切向刚度。

(2m[p]R[p]R[p]/5)}

(3)

Clump像刚体一样运动，其形心的动量Mi可用式(4)进行计算：

(4)

(j≠i,l≠k)

(5)

通过联立式(4)和式(5)便可得到clump的旋转运动求解。运动方程可以被写成方程(6)的形式：

(6)

通过求解clump单元在t+Δt时刻的质心位置可用式(7)表示：

(7)

在求得clump的平移和旋转速度之后，组成clump的每个粒子的速度可以用方程(9)表示为：

(8)

对所有离散单元不断进行上述迭代计算，从而可以最终实现对所研究问题的动态过程计算模拟。

2 数值验证与可靠性分析

本文在Jiang等[5]所作试验的基础上进行数值模拟分析，虽然该物理模型试验并未按照相似准则进行设计，但由于其较为准确的测量了下游刚性挡板的受力过程而被国内外诸多学者用于对碎屑流模拟的数值验证。在对刚性网状防护结构拦阻碎屑流进行系统研究之前，利用Jiang等所作试验对数值模拟的正确性验证与可靠性分析。

试验主要在如图2所示的试验装置中进行。试验中所用流通槽的高度为2.93 m，长度为0.30 m，宽度为0.35 m，流通槽底部与水平面之间的角度通过围绕底部固定点的旋转调节。在试验槽下游设置了垂直于试验槽底部的刚性挡板，刚性挡板的法向力可以通过传感器来测量。在试验前的初始时刻，将粒径为20 mm的石灰石砾石置于上游流动槽的方箱中。L和H分别是方箱的长度和高度。起动板被瞬间打开后碎屑材料开始在流通槽内流动，起动板与刚性挡板间的距离为2.19 m，试验流通槽的底部、侧面和刚性挡板的摩擦角可以从试验中获得，分别为25°、15°和21°。根据这些摩擦角可进一步确定出相应的摩擦系数。

图2 试验装置示意图[5]

由于实际的碎屑颗粒形状是不规则的，仅采用球形颗粒单元模拟碎屑流的运动过程有时效果不佳。为得到更好的模拟效果，在数值计算中采用4个相同大小的ball组成模拟石灰石砾石的clump单元，其中任意两个球间的距离相等。利用试验数据进行参数标定和优化后即可确定当球心距离D为球半径R的1.5倍时(即D=1.5R)模拟的效果最佳[19]，如图3所示。为使生成的clump单元和直径约为0.02 m的石灰石砾石的体积和密度相等，对球的半径和密度进行了校准，最终确定组成clump的球半径为0.066 m，密度为2.34×103kg/m3。为与试验的物理力学参数相一致采用了如表1所示的细观计算参数。

图3 Clump离散单元[19]

表1 数值模拟中选用的计算细观参数

为使模拟中的clump数量与实验中的砾石粒数近似相等，采用np=Vtγt/Vsγs计算clump的数量np。其中，Vt和γt分别为试验所用砾石的总体积和密度，而Vs和γs分别为clump团块的体积和密度。

本文选取文献[5]中编号为L44-H15-α40°的实验进行数值模拟计算和验证，如图4所示。试验编号所表示的试验条件为L=44 cm，H=15 cm，α=40°。通过离散元法的计算得到了在该编号试验中分别作用于刚性挡板上单宽法向力和力矩的变化，计算结果与试验数据的比较如图4所示。

图4 单宽法向力和力矩时程曲线(L44-H15- α40°)

从计算结果来看，尽管在某些细节上存在一定差异，但总体计算结果与试验数据基本吻合一致，从而也验证了计算的正确性和有效性，说明利用数值模拟进一步研究多级刚网结构拦阻碎屑流过程是可行且可靠的。数值模拟可以反映碎屑流对刚性挡墙的冲击过程，在此基础上进一步对刚性网状结构拦阻碎屑流的过程进行更深入的研究和讨论。

3 数值结果与分析

由于数值模拟计算可以得到拦挡结构拦阻碎屑流过程中的拦截效率、以及碎屑与拦阻结构的受力特征，这也启示我们可以利用数值模拟对一些难以应用现代技术测试的防护结构(如刚性网状结构)拦阻碎屑流过程进行研究。

为了进一步研究在不同设置条件下刚网结构的拦截效率及受力特征，在算例验证的基础上，在流通槽中的不同位置设置与流通槽相同宽度与高度的刚网结构，模拟刚性拦挡结构对碎屑流的拦阻效果，如图5所示。

图5 多级刚网结构在流通槽中的设置示意图

计算采用与数值验证算例具有相同参数的clump单元进行模拟，刚网结构则采用与侧墙参数相同的基于线接触模型和滑移模型的线墙单元进行模拟。在距离下游挡板L1=0.55 m处设置不同网孔尺寸单级刚网结构条件下计算得到了设置不同网孔尺寸条件下的拦截效率，如表2所示。

表2 不同网孔尺寸设置下单级刚网结构的颗粒拦阻效率

从表2数据可以看出，当采用较小的网孔尺寸时，利用单层刚网结构对碎屑流进行拦阻基本就可以达到较好的拦阻效果。当网孔尺寸达到二倍粒径0.04 m时刚网结构的拦截效率依旧较高为78.70%，这与我们直观感觉略有不同。在网孔尺寸大于颗粒粒径的情况下，虽然到达刚网结构的部分颗粒仍可能会因粒径小于网孔尺寸而穿透刚网结构，但在刚网结构前颗粒间通过相互挤压支撑所形成的稳定的传力结构可实现对颗粒的有效拦阻，因而较大网孔尺寸的刚网结构仍可拦截小于网孔尺寸的颗粒而获得较高的拦截效率。而当刚网结构的网孔尺寸超过0.04 m时，在刚网结构前由颗粒相互挤压形成的空间网状传力结构不再能有效阻止颗粒的持续穿透，因而此时利用网状传力结构拦截颗粒的力学机制便不再有效。由于单级刚网结构的受力特征已有文献详加阐述[20]，本文不再做过多讨论。在网孔较大的条件下，虽然单层刚网结构的拦截效率并不理想，但其仍然能够起到一定的拦阻和耗能的作用，因此可以进一步利用多级刚网结构对碎屑流进行拦阻。

为研究多级刚网结构布置对碎屑流动力过程的影响，在流通槽内设置了多级不同网孔尺寸的刚网结构，其中为了保证最后较高的拦截效率，末级刚网结构采用了能够实现较高拦截效率的网孔尺寸。经过计算得到了在不同多级刚网结构设置下的拦截效率，如表3所示。

表3 不同网孔尺寸设置下多级刚网结构的颗粒拦阻效率

从表3所示结果可以看出为保证最终有效的拦截效率，Case_B系列的工况中所设置的最后一级网孔尺寸都比较小，所以在多级刚网结构拦阻碎屑流的计算中拦阻效率都比较高，对下游刚网结构网孔尺寸的控制是提高碎屑流拦截效率最有效的措施之一。在末级孔径相同的条件下设置三级刚网结构的工况相比设置二级刚网结构的工况拦截效率更高，其主要的原因是碎屑流的能量在多级刚网结构拦阻条件下逐级减小，从而使颗粒快速在网前形成堆积，并且也使通过相互碰撞从刚网结构上方通过的颗粒也有所减少，在本文的算例中由于末级刚网结构尺寸均较小从而使这种差别并不明显。在最下游刚网结构位置相同的条件下，增加间距使处于中游和上游的刚网结构位置更高，这意味着碎屑流释放的势能要更小，或者说碎屑颗粒达到刚网结构时的速度更小，因而由高速运动穿透刚网结构到下游的颗粒也得以减少，体现在整体的拦截效率上则是随着刚网结构间距的增加拦截效率增大。

为了便于考察多级刚网结构拦阻碎屑流过程的特点，选取能够突出拦阻过程多项特征的Case-B7的多级刚网结构拦阻碎屑流过程进行研究，得到在不同时刻的运动速度及形态过程如图6所示。

图6 Case-B7工况的碎屑流运动过程

从图6可见，当颗粒到达刚网结构后，碎屑流的速度显著减小，粒径小于网孔尺寸的颗粒穿过刚网结构向下游运动，穿透上一级刚网结构的颗粒在下一级刚网结构的拦阻作用下可发生反弹、减速和继续穿透运动。在刚网结构的拦阻作用下网前的颗粒通过传力网络形成稳定的支撑以拦阻颗粒，后续的颗粒在网前颗粒所形成静力区影响下进行爬升运动，部分颗粒通过碰撞作用从网上运动向下游。随着拦阻过程的进行，流在网前达到新的高度后仍会有部分颗粒穿过刚网结构向下游运动，从而可见颗粒按照自下而上的顺序陆续穿过刚网结构，从较高处穿透刚网结构的颗粒会以抛体运动的形式向下游运动。在间距有限的条件下刚网结构的受力通过网间颗粒建立关联，而当间距足够大时，各级结构间的颗粒不再将两者紧密关联，刚网结构的受力则相对独立。

设置多级刚网结构可使各级结构共同分担碎屑流的冲击力，以此增加工程结构的可靠性，因而需要更进一步的考察各级结构的受力特征。为比较在间距设置相同条件下网孔尺寸对冲击过程的影响，计算Case_B1-B4工况的刚网结构受力时程曲线，如图7所示。

图7 碎屑流冲击各级结构的时程曲线(Case_B1-B4)

从图7中可以看到，对于2级刚网结构而言，虽然Case_B1和Case_B2工况的拦截效率几乎相同，但两种网孔尺寸刚网结构设置条件下各级结构网的受力分配并不相同，Case_B2下游刚网结构相比其上游刚网结构受到碎屑流的峰值和稳定力都比较大，这说明2级刚网结构对拦挡碎屑流起到的耗能作用相对较小，从而使下游刚网结构的受力相对更大，Case_B2的设置没有很好的起到逐级耗能以均衡各级防护结构受力的作用，而起到最终防护作用的末级刚网结构所受力过高是较为危险的。Case_B1上游和下游刚网结构受到碎屑流作用的峰值和稳定力相差不大，相比Case_B2更加合理。Case_B3和Case_B4设置了3级刚网结构对碎屑流进行拦阻，对于Case_B3而言由碎屑流冲击造成的受力主要由前两级刚网结构承担，而末级的受力则非常小，由此可见虽然末级防护结构的保障作用十分重要，但Case_B3的设置没有发挥对碎屑颗粒的拦阻作用在实践中会造成浪费。从受力上来看Case_B3工况和Case_B1工况的冲击时程是类似的，起到拦阻颗粒作用的也主要是孔径为0.04 m和0.02 m的两层刚网结构。在Case_B4和Case_B1的比较中，由于Case_B4在上游孔径为0.08m的刚网结构首先对碎屑流进行拦阻并有效耗散其能量，所以Case_B4在中游和下游的刚网结构两工况中所受的峰值和稳定力相对Case_B1较小，显然在条件允许的情况下在上游位置多设置一层孔径较大的刚网结构进行耗能是十分必要的。在综合了用料用量、受力均衡和安全可靠的原则下，在Case_B1-B4系列中Case_B4显然是较为理想的。

为研究在相同网孔设置条件下多级刚网结构间距设置对冲击过程的影响，计算Case_B4-B7的刚网结构受力时程曲线，如图8所示。

图8 碎屑流冲击各级结构的时程曲线(Case_B4-B7)

由图8可见在设置相同网孔尺寸的三层刚网结构条件下， 不同的间距将会使各级结构的受力有很大的不同。 在间距较小时， 穿透上级刚网结构的颗粒快速被下一级刚网结构拦阻并减速， 而减速的颗粒将对上游陆续穿透的颗粒进行拦阻， 甚至可能造成颗粒对上游刚网结构的反向作用力， 从而使上一级刚网结构所受的力大幅减小， 这在Case_B4和Case_B5工况中第1级荷载的受力中显得尤为明显。随着刚网结构间距不断增加，在Case_B5-B7工况中受到下游刚网结构的影响减弱，被第2级拦阻的颗粒对第1级刚网结构的影响不再那么显著，因而峰值受力相差不大，但稳定受力则随下游堆积空间的增加而减小。在第2级刚网结构中Case_B4的单宽受力相比其他工况要大一些，其主要的原因是由于间距过小所限，处于前两级刚网结构间的颗粒将荷载传递给第2级，这与该工况第1级受力偏小是对应的。Case_B7的第2级受力峰值并不明显是由于随间距增加，第2级所处位置更靠上游，颗粒达到刚网结构时的速度相对较小，并且在经过第1级拦截后减速颗粒在第2级网前的爬升运动不显著。同时由于网间距足够大，在网前更小的堆积高度及下游颗粒的被动作用力使其稳定受力也较小。在第3级刚网结构中几个工况的受力峰值均有所减弱，这是由于经过前两级刚网结构的拦阻作用后，碎屑流的速度大为减小，因而导致碎屑流在第3级刚网结构上的爬升运动并不显著，因而所受冲击力的峰值并不明显。而且随着间距的增加第3级网受力依次减小，这是由于稳定的网前堆积高度随下游空间的增加而减小。从结果可见，在相同多级刚网结构孔径设置条件下，改变刚网结构的间距将会对各级结构的受力产生较大的影响。在Case_B4-B7系列中从各级受力均衡的角度考虑，Case_B4和Case_B5是设置较为合理的工况。

4 结 论

本文采用离散单元法对多级刚网结构拦阻碎屑流过程进行了数值模拟，在进行数值验证的基础上，结合多级刚网结构拦阻下碎屑流的运动特征研究了分级、网间间距和各级网孔尺寸布置对拦截效率及各级结构受力的影响。研究表明在多级刚网结构的拦阻作用下，碎屑流将产生反弹、穿透和爬升等运动形式，这些运动是各级结构受力产生波动和峰值等特征的主要原因。在刚网结构间距有限的条件下，网间颗粒将作为传力介质将各级结构的受力进行关联，进而使结构的受力过程变得复杂。基于离散元法的数值计算可以得到多级刚网结构的拦截效率、各级结构的受力过程等防护结构设计中的重要指标，而这些指标常不能通过简单的计算而得到。在测试条件受限时利用数值模拟所得到的结果可为多级刚网结构的网孔尺寸和位置布局等优化设计提供参考。