0.   引 言

崩塌落石灾害是我国西部地区三大地质灾害之一^{[1 − 3]}。由于落石具有多发性、突发性、随机性、难预测、能量大等特点^{[4 − 6]}，常常会对山区公路、铁路等设施造成巨大的威胁^{[7 − 11]}。棚洞是防护公路、铁路等免受落石冲击破坏最有效的措施之一^[12]。典型的棚洞主要由混凝土梁柱、混凝土板以及覆盖在混凝土板上的土垫层构成。土垫层的主要作用是避免落石直接冲击混凝土板，缓冲落石的冲击能，减小冲击力，并将落石冲击力扩散到更大的范围（图1）。虽然现在已经有较多的棚洞设计规范，但是依然存在落石冲击穿透土垫层，导致棚洞主体结构破坏的现象^{[13 − 14]}。因此，开展落石冲击土垫层动力响应研究，有利于优化土垫层的设计，提升棚洞防护措施的有效性，增强崩塌落石灾害的防灾减灾能力。

图  1  崩塌落石与棚洞示意图

Figure  1.  Schematic diagram of rockfall and rock shed

下载: 全尺寸图片 幻灯片

落石冲击土垫层是一个非常复杂的过程，以快速加卸载、大变形、快速的能量转换和应力波传播为特征。国内外学者一直都在探索描述落石冲击力的理论计算方法。现在已形成的理论包括：赫兹弹性接触理论^[15]、赫兹弹塑性接触理论^{[16 − 19]}、能量守恒原理^{[20 − 21]}、冲量定理^[22]、地基承载力理论^[23]、BIMPAM流变学理论^[24]和基于Logistic函数的落石碰撞分析方法^[25]。

除理论研究外，国内外学者一直以来都在开展落石冲击土垫层的试验研究，包括小尺度、中等尺度试验和少量的原型试验。罗杰等^[26]采用试验研究了四种土壤（砂土、黏土、壤土和黄棕壤）的缓冲性能。研究表明，砂土的缓冲性能最佳。王林峰等^[27]基于小型的落石棚洞模型，研究了落石重量、下落高度和棚洞顶板倾角对落石冲击力的影响，研究表明落石重量是影响落石冲击力的首要因素，其次是下落高度，最后为棚洞顶板倾角。Calvetti等^[28]采用大尺度的试验研究了落石冲击土垫层的现象，研究表明垫层表层土的密度是影响落石冲击力的重要因素，而土垫层的倾斜角度影响不大。Kawahara和Muro^[29]研究了土垫层密度和厚度对落石冲击力的影响，研究表明落石的冲击力随垫层密度的增大而增大，土垫层与棚洞顶板的作用力随垫层厚度的增大而减小。

数值模拟是研究落石冲击土垫层的一种有效方法。数值模拟的优势在于：费用低，可重复性强，可以分析得到试验中无法获取的信息，可以节省大量的人力物力，避免仪器设备等被损坏。因为土垫层自身具有离散特性，因此离散单元法被广泛用于落石冲击土垫层的数值研究。王林峰等^[30]采用离散元软件（PFC^2D）研究了落石半径和密度对冲击力的影响。江巍等^[31]运用离散元软件研究了素填土、粉质黏土和砂质粉土的缓冲能力。Calvetti等^[28]应用离散元方法研究了落石冲击能量对落石冲击力的影响。Zhang等^[32]采用三维离散元法研究了落石冲击土垫层的反弹现象，分析了落石反弹与落石尺寸和垫层厚度的关系。上述研究表明，离散元法是研究落石冲击土垫层的一种有效方法。

综上所述，目前已经有较多关于落石冲击土垫层的成果，大多数的研究结果表明，土垫层厚度越大，落石冲击力越小，土垫层缓冲效果越好。但是，缓冲效果是否会一直增加，以及垫层厚度如何选择，现有研究还未回答。因此，本项研究拟采用离散元法探讨垫层厚度和下落高度对土垫层缓冲落石冲击力特性的影响，以期为土垫层的设计提供理论指导。

1.   落石冲击土垫层离散元数值模型

1.1   离散元理论简介

本项研究数值模拟采用开源离散元软件：ESyS-Particle^[33]。基于分子动力学的思想，离散单元法将土模拟为球形颗粒的集合体。在荷载作用下，颗粒可以发生平动和转动。在计算过程中，颗粒间被赋予一定的接触模型，两个相互接触的颗粒通过接触模型产生接触力。通过计算每个颗粒所受的合力（$ {{\boldsymbol{F}}_i} $）和合力矩（$ {T_i} $），并根据牛顿第二定律，采用显示积分的方法更新颗粒的速度和位置，如式（1）和式（2）所示。

式中：$ {m_i} $、$ {{\boldsymbol{r}}_i} $——第i个颗粒的质量和位置；

$ {I_i} $、$ {{\boldsymbol{\omega }}_i} $——第i个颗粒的转动惯量和转动速度。

本项研究中，颗粒间的接触模型采用无黏结摩擦模型，如图2所示。无黏结摩擦模型包括颗粒间的法向线性接触模型，见图2（a），切向接触模型，见图2（b），和抗转动接触模型，见图2（c）。为了模拟真实土颗粒间的非弹性碰撞，在法向接触模型中引入阻尼，如图2（a）所示。同时为了考虑真实土颗粒形状的影响，引入抗转动接触模型。

图  2  无黏结摩擦模型

Figure  2.  Cohesionless particle contacting model

下载: 全尺寸图片 幻灯片

根据图2所示的接触模型，两个颗粒间的接触力包括：法向接触力（$ {F_{\rm{cn}}} $）、法向阻尼力（$ {F_{\rm{cd}}} $）、切向接触力（$ {F_{\rm{cs}}} $）和滚动力矩（$ {M_{\rm{cr}}} $）。法向接触力由式（3）计算。

式中：$ {k_{\rm{cn}}} $——法向接触刚度；

$ {u_{\rm{cn}}} $——两个颗粒接触处的重叠距离。

法向接触刚度由式（4）计算。

式中：$ {E_{\rm{p}}} $——颗粒的杨氏模量；

$ {R_{\rm{A}}} $、$ {R_{\rm{B}}} $——两个接触颗粒的半径。

接触处的法向阻尼力由下式计算：

式中：$ \beta $——阻尼系数；

$ {m_{\rm{A}}} $、$ {m_{\rm{B}}} $——两个接触颗粒的质量；

$ {v_{\rm{cn}}} $——两个接触颗粒的法向相对速度。

接触处的切向接触力（$ {F_{\rm{cs}}} $）采用理想弹塑性模型，其线性阶段按增量的形式来计算：

式中：$ F_{\rm{cs}}^t $、$ F_{\rm{cs}}^{t - \Delta t} $——当前和前一个计算时步的切向力；

$ {k_{\rm{cs}}} $——切向接触刚度；

$ {u_{\rm{cs}}} $——两个颗粒在接触处的切向相对位移。

切向接触刚度由式7计算。

式中：$ {\upsilon _{\rm{p}}} $——颗粒的泊松比。

切向接触力的最大值由摩尔库仑定律控制，如式（8）。

式中：$ {\mu _{\rm{p}}} $——颗粒的摩擦系数。

滚动力矩用于考虑颗粒形状的影响，采用理想弹塑性模型，其计算方法为式（9）。

式中：$ M_{\rm{cr}}^t $、$ M_{\rm{cr}}^{t - \Delta t} $——当前时步和前一个时步的滚动力矩；

$ {k_{\rm{cr}}} $——抗滚动刚度，$ {k_{\rm{cr}}} = {k_{\rm{cs}}}({R_{\rm{A}}} + {R_{\rm{B}}})/2 $；

$ \Delta {\theta _{{\mathrm{r}}}} $——在一个计算时步内，两个接触颗粒的相对转 动角度。

最大的滚动力矩（$ M_{\rm{cr}}^{\max } $）定义为：

式中：$ {\eta _{\rm{p}}} $——颗粒塑性力矩系数。

1.2   数值计算模型

落石冲击土垫层的离散元数值模型如图3所示。该模型与文献[34]中的室内物理模型一致。数值模型由落石、土垫层和混凝土底座三部分构成。落石模拟为一个直径（D）为20 cm，质量为11.5 kg的球形颗粒。土垫层模拟为长1.0 m，宽1.0 m，厚度为H的立方形颗粒集合体。土垫层颗粒的直径均匀分布在1.0到2.0 cm之间。混凝土底座模拟为一层直径为1.0 cm的固定颗粒，该层颗粒不会发生平动和转动，但是可以与土垫层颗粒接触产生接触力。土垫层的生成过程包括两个步骤：首先在四个刚性墙和底座围城的矩形盒子内，随机生成规定半径范围内的颗粒；然后通过重力沉积作用形成指定厚度的颗粒层。颗粒之间的接触模型均为无黏结摩擦模型。数值模型的输入参数如表1所示。

图  3  落石冲击土垫层离散元模型

Figure  3.  Discrete element model of rockfall impacting the soil cushion layer

下载: 全尺寸图片 幻灯片

表  1  数值模型输入参数

Table  1.  Input parameters of the numerical model

变量	数值
土垫层颗粒直径/cm	[1.0, 2.0]
土颗粒密度/（kg·m−3）	2698.2
颗粒杨氏模量/MPa	1×102
颗粒泊松比	0.25
颗粒阻尼系数	0.01
颗粒摩擦系数	0.6
颗粒塑性力矩系数	0.15
计算时步/s	10−6
重力加速/（m·s−2）	9.81

下载: 导出CSV 

| 显示表格

数值模拟过程中，落石被置于土垫层的正上方，并根据落石下落高度（h_f）设定初始速度（v₀）。初始速度和下落高度的关系如式（11）所示。

本项研究中，落石下落高度有5种，包括3.0，5.0，10.0，20.0，30.0 m。土垫层的厚度有4种，包括10.0，20.0，30.0，40.0 cm。因此，总共进行20组数值试验。为了评估土垫层的缓冲特性，提取了落石的冲击力峰值（$ F_{{\rm{block}}}^{\max } $），以及土垫层与底座接触面中心位置的峰值力（$ F_{\text{c}}^{\max } $）。土垫层与底座中心位置的接触力可以看作是土垫层与棚洞顶板中心位置的接触力。因此，$ F_{\text{c}}^{\max } $与$ F_{{\rm{block}}}^{\max } $的比值越小，表明土垫层的缓冲效果越好。

2.   计算结果分析

2.1   数值模型验证

通过与文献[34]报道的试验结果对比，本项研究首先验证了上述数值模型的有效性。图4给出了落石以3 m下落高度冲击30 cm厚土垫层情况下，落石冲击力和顶板中心力随时间演化曲线。由图中可以看出，从定性的角度，数值模拟结果能基本再现落石冲击力和顶板中心随时间的演化趋势；从定量的角度，数值模拟结果能再现落石冲击力峰值和顶板中心力峰值。因此，上述数值模型以及所选参数是有效的。

图  4  落石冲击力和顶板中心力随时间演化曲线（h_f = 3.0 m，H = 30.0 cm）

Figure  4.  Time evolution curves of the impact force and the roof center force (h_f = 3.0 m, H = 30.0 cm)

下载: 全尺寸图片 幻灯片

2.2   下落高度的影响

图5给出了落石冲击不同厚度土垫层情况下，落石冲击力峰值（$ F_{{\rm{block}}}^{\max } $）与下落高度（$ {h_{\rm{f}}} $）的关系。从图中可以看出，落石峰值冲击力随下落高度的增大而增大。垫层厚度小于落石直径时的峰值冲击力明显大于其它垫层厚度情况，而且随着下落高度的增加，越来越明显。从图中还能看出，无论垫层厚度为多少，落石峰值冲击力与下落高度都可以用统一的式（12）来表示。

图  5  落石冲击力峰值与下落高度的关系

Figure  5.  Dependence of the maximum rockfall impact force on the falling height

下载: 全尺寸图片 幻灯片

式中：$ {F_0} $和$ {h_0} $——拟合参数。

研究表明，$ {h_0} $取30.0 m，$ {F_0} $取下落高度为30.0 m的峰值冲击力时，可以达到较好的拟合效果。

图6给出了落石以不同下落高度冲击不同厚度土垫层情况下，顶板中心力峰值（$ F_{\text{c}}^{\max } $）与下落高度（$ {h_{\rm{f}}} $）的关系。从图中可以看出，顶板中心力峰值与下落高度呈线性关系。当垫层厚度小于落石直径时（H/D = 0.5），$ F_{\text{c}}^{\max } $随$ {h_{\rm{f}}} $的增长率（拟合直线的斜率）明显高于其它情况。随着垫层厚度的增加，$ F_{\text{c}}^{\max } $随$ {h_{\rm{f}}} $的增长率减小。当垫层厚度由0.5倍落石直径增加到1.0倍落石直径时，$ F_{\text{c}}^{\max } $随$ {h_{\rm{f}}} $的增长率由504.7 N/m减小到372.1 N/m。当垫层厚度增加到1.5倍落石直径和2.0倍落石直径时，$ F_{\text{c}}^{\max } $随$ {h_{\rm{f}}} $的增长率分别变化为87.0 N/m和为48.2 N/m。这表明，随着垫层厚度的增加，下落高度对顶板中心力峰值的影响逐渐减小。

图  6  顶板中心力峰值与落石下落高度的关系

Figure  6.  Dependence of the maximum roof center force on the falling height

下载: 全尺寸图片 幻灯片

2.3   土垫层厚度的影响

图7给出了落石冲击力峰值（$ F_{{\rm{block}}}^{\max } $）与垫层厚度和落石直径比值（H/D）之间的关系。从图中可以看出，随着垫层厚度的增加，落石的峰值冲击力减小。当H/D从0.5增加到1.0时，即垫层厚度从落石直径的0.5倍增加到1倍时，落石的冲击力峰值减小将近50%。当土垫层的厚度继续增加时（H/D > 1.0），落石峰值冲击力变化不大。并且，从图中可以看出，无论落石的下落高度为多少，即无论落石的冲击速度为多少，落石的峰值冲击力与土垫层厚度的关系均出现上述现象，即在土垫层厚度增加到1倍直径后，土垫层厚度对冲击力影响较小。

图  7  落石冲击力峰值与标准化垫层厚度的关系

Figure  7.  Dependence of the maximum rockfall impact force on the normalized soil cushion thickness

下载: 全尺寸图片 幻灯片

图8给出了落石以不同下落高度冲击不同厚度土垫层情况下，顶板中心力峰值（$ F_{\text{c}}^{\max } $）与垫层厚度和落石直径比值（H/D）的关系。从图中可以看出，无论落石下落高度（冲击速度）为多少，随着垫层厚度的增加，顶板中心力峰值不断减小，$ F_{\text{c}}^{\max } $与H/D呈负指数幂函数关系，表明$ F_{\text{c}}^{\max } $随H/D减小的速度不断变小。相比于0.5D的情况，当垫层厚度增加到一倍落石直径时（H = 1.0D），顶板中心力峰值减小64%；当垫层厚度增加到1.5D时，顶板中心力峰值减小86%；当垫层厚度增加到2.0D时，顶板中心力峰值减小92%。因此，垫层厚度从1.5D增加到2.0D仅仅使顶板中心力峰值减小6%。由此可见，在垫层厚度增加到1.5D后，继续增加垫层的厚度，土垫层缓冲效果（顶板中心力的减小量）增加不明显。结合土垫层厚度对落石冲击力峰值的影响，可以得出垫层厚度取落石直径的1.5倍较为合适。

图  8  顶板中心力峰值与标准化垫层厚度的关系

Figure  8.  Dependence of the maximum roof center force on the normalized soil cushion thickness

下载: 全尺寸图片 幻灯片

图9给出了落石以不同下落高度冲击不同厚度的土垫层情况下，顶板中心力峰值与落石冲击力峰值的比值（$ {{F_{\mathrm{c}}^{\max }}/{F_{{\rm{block}}}^{\max }}} $）与下落高度的关系。从图中可以看出，在垫层厚度为0.5倍落石直径情况下，$ {{F_{\mathrm{c}}^{\max }} / {F_{{\rm{block}}}^{\max }}} $随着下落高度的增大而减小。对比图4和图5，可以发现，这是由于在垫层厚度小于落石直径的情况下，落石下落高度对冲击力峰值的影响高于对顶板中心力的影响。当土垫层的厚度增大到落石的直径的1.5倍时（H/D = 1.5），对于同一厚度垫层，$ {{F_{\mathrm{c}}^{\max }} /{F_{{\rm{block}}}^{\max }}} $基本上不随下落高度变化，表明此时，垫层的缓冲效果不受落石下落高度的影响。此外，对于H/D = 1.0、1.5和2.0情况下的$ {{F_{\mathrm{c}}^{\max }} / {F_{{\rm{block}}}^{\max }}} $平均值分别为0.097、0.034和0.02。顶板中心力峰值与落石冲击力峰值的比值随垫层厚度的增大而减小，表明垫层缓冲作用随垫层厚度的增大而增大。当H/D从1.0增加到1.5时，$ {{F_{\mathrm{c}}^{\max }}/{F_{{\rm{block}}}^{\max }}} $减小0.063；当H/D从1.5增加到2.0时，$ {{F_{\mathrm{c}}^{\max }}/ {F_{{\rm{block}}}^{\max }}} $仅减小0.014。表明，在土垫层厚度增加到1.5倍落石直径后，继续增加垫层的厚度，垫层的缓冲效果增加不明显。

图  9  顶板中心力峰值和落石冲击力峰值的比值与落石下落高度的关系

Figure  9.  Dependence of the ration between the maximum roof center force and the maximum rockfall impact force on the falling height

下载: 全尺寸图片 幻灯片

3.   结论

基于离散单元法，建立落石冲击土垫层的数值模型，研究不同厚度土垫层缓冲落石冲击力的特性，得到以下结论：

（1） 在土垫层厚度一定的情况下，落石冲击力峰值与落石下落高度呈幂函数关系；顶板中心力峰值与下落高度呈线性正相关关系。

（2） 在下落高度一定的情况下，顶板中心力峰值与垫层厚度呈负指数幂函数关系；随着垫层厚度的增加，落石冲击力峰值减小，当垫层厚度增加到落石直径的1.0倍之后，落石冲击力峰值与垫层厚度无关。

（3） 随垫层厚度的增大，顶板中心力峰值与落石冲击力峰值的比值减小，垫层缓冲效果增大；当垫层厚度增加到落石直径1.5倍之后，垫层缓冲效果增加不明显。垫层厚度建议取值为落石直径的1.5倍。