1. 引言
落石灾害在我国山区铁路建设、运行过程中经常发生，导致列车受损，轨道损毁，人员伤亡。业界人士针对落石对结构的冲击效应展开系统研究。在研究初期，落石冲击效应的研究多以试验模型为主，如隧道设计技术手册 [1] 、公路路基设计规范 [2] ，西南交通大学(杨其新等) [3] 、日本道路协会 [4] 和瑞士学者 [5] 都是基于试验提出了落石冲击力的计算方法，由于这些方法有局限性，计算结果差异较大；后来，很多学者引入了有限单元法来模拟落石冲击过程；近几年，越来越多的学者开始采用离散元(Discrete Element Method，DEM)模拟落石冲击过程。
近年来，很多学者对于棚洞结构的落石冲击效应方面的研究有了很大的进展。裴向军，刘洋，王东坡等 [6] 通过有限元数值模拟落石冲击砂土垫层，得出砂土垫层可以有效减小落石对棚洞顶板冲击效应的结论，并经过优化得到合理垫层厚度；杨建荣、白羽、杨晓东等 [7] 提出了一种由金属柔性网和弹簧撑杆组成的柔性棚洞并采用有限元方法进行数值模拟，计算动力响应；王静峰、赵鹏、袁松等 [8] 基于大尺度模型试验和有限元数值模拟相结合的方法，对EPS和砂土垫层组成的复合垫层钢棚洞抗落石冲击性能展开了研究，得出复合垫层可以有效提高钢棚洞对落石的抵抗能力的结论；吴佳楠 [9] 同样以钢棚洞模型冲击试验和有限元数值模拟相结合的方法，对钢棚洞在EPS-细砂复合垫层下对落石的冲击效应展开研究，得出峰值加速度、最大应变与下落高度成正比的结论；王永东，周天跃，柴伦磊等 [10] 基于能量分析方法对落石冲击棚洞结构进行研究，考虑垫层弹塑性、垫层与落石之间摩擦力、棚洞自身耗能三方面影响因素，得到了冲击深度以及冲击压力的计算公式，还得出垫层成坑变形耗能最大达到50%、垫层与落石摩擦力耗能次之，达到30%以及棚洞自身耗能占比20%的结论；闫帅星，何思明，李新坡等 [11] 基于冲击动力学、修正的混凝土塑性特性及永久变形分阶段弹性接触理论，构建了钢筋混凝土棚洞在落石冲击作用下动力响应的一般解析方法，并通过有限元进行验证。
本文采用ITASCA公司开发的软件PFC来建立落石、框架棚洞、回填土等模型，模拟落石垂直下落在有回填土和无回填土两种情况下，对棚洞结构的冲击过程，通过在棚洞不同部位(如顶板、立柱、纵梁等)布置监测点来获取相应部位的应力–应变时程曲线，从而对框架棚洞的冲击效应展开研究并总结。
2. 落石冲击棚洞结构的力学模型和本构模型
本次建模包括物理模型和本构模型。物理模型包括山体、棚洞结构、回填土层、落石等分模型，棚洞结构组成部分主要有顶板、盖梁、支柱、基础等；本构模型主要讲述不同物理模型之间的相互作用，比如落石与回填土、回填土与棚洞结构之间的相互作用。本次建模主要采用接触本构模型，包括颗粒与颗粒之间的粘结和滑移模型以及颗粒与墙体之间的接触刚度模型 [12] 。
2.1. 物理模型
1) 棚洞结构模型
PFC3D模型中颗粒间通过平行粘结方式粘结，平行粘结模式有五个参数，其中平行粘结法向刚度pb_kn和平行粘结切向刚度pb_ks，单位是Pa/m，代表颗粒间经过平行粘结后的法向和切向刚度，影响粘结模型的变形特性，而平行粘结法向强度pb_ns和平行粘结切向强度pb_ss单位是Pa，代表的是颗粒间法向和切向粘结强度，影响粘结模型的强度特性；平行粘结半径放大系数pb_rad，代表的是平行粘结半径放大系数，粘结半径等于最小颗粒半径乘以这个半径放大系数。
棚洞结构由顶板、盖梁、支柱、基础组成，全部用球体颗粒单元通过平行粘结来模拟。平行粘结主要有五个微观参数，分别是pb_kn、pb_ks、pb_ns、pb_ss和pb_rad，使得模拟材料的参数属性与钢筋混凝土接近。
2) 回填土、山体和落石模型
回填土以及落石都用球体颗粒来模拟，山体坡面采用墙体单元模拟。
2.2. 力学本构关系
本次离散元数值模拟采用接触刚度模型，回填土颗粒单元采用接触粘结本构模型，棚洞颗粒单元采用平行粘结结构模型。
3. 棚洞结构建模
采用离散元颗粒流数值模拟软件PFC建模，模拟工况是有回填土(t = 1 m)和无回填土(t = 0)两种情况,落石重量2.0 t，垂直下落冲击棚洞结构。落石冲击时，需要考虑落石重力、回填土的摩擦力、回填土对落石的粘滞阻尼力和落石与回填土以及棚洞结构的碰撞力，阻尼选用粘滞阻尼，粘滞阻尼系数取值0.2。
落石冲击棚洞结构的物理模型和数值建模
1) 模型概述
棚洞结构设计图如图1所示，颗粒流PFC3D建模如图2所示。
(a) 横断面 
(b) 纵断面
Figure 1. Design drawing of shed tunnel structure [13]
图1. 棚洞结构设计图 [13]
(a) 无回填土(t = 0)
(b) 回填土(t = 1 m)
Figure 2. PFC3D model of rockfall and shed tunnel
图2. 落石、棚洞PFC3D模型
2) 离散元模型单元微观参数
模型微观参数 [14] 设置如表1所示。
3) 回填土参数标定
根据土力学 [15] 、工程地质手册 [16] 、弹性力学 [17] 、以及相关学术论文 [18] [19] 的方法进行回填土参数标定，标定结果如表2所示。
Table 2. PFC3D numerical simulation macroscopic parameter results of Triaxial test [8]
表2. 三轴试验PFC3D数值模拟宏观参数结果 [8]
4) 工况组合
本次数值模拟工况组合表3。
Table 3. Combination of working conditions [8]
表3. 工况组合 [8]
4. 棚洞结构的落石冲击效应研究
监测点布置情况见图3(a)~(b)，单元三维空间应力符号说明见图4，图中所示应力均是应力正方向，其中三个正应力σxx、σyy、σzz正值表示拉应力，负值表示压应力 [11] 。通过测量球(measure sphere)获取每个监测点的六个应力分量，分别是三个正应力分量σxx、σyy、σzz和三个切应力分量τxy、τyz、τzx，绘制不同监测部位的应力时程曲线，从而得出该部位的控制应力(最大应力值)。
(a) 
(b)
Figure 3. (a) Layout of monitoring points for roof; (b) Layout of monitoring points for Stringer, pillar and base
图3. (a) 顶板监测点布置图；(b) 纵梁、支柱和基础监测点布置图
Figure 4. Symbol of three-dimensional space stress [11]
图4. 单元三维空间应力符号 [11]
用三种方法研究落石冲击效应，三种方法介绍如下：
1) 单个监测点位应力时程曲线
在离散元数值模拟软件PFC中用fish语言可以实现用测量小球(measure sphere)读取监测点部位(例如选取监测点18冲击部位顶板中心)的六个应力分量，在模拟计算的过程中自动记录应力分量，绘制应力时程曲线，进而可以对该监测点的最大应力进行分析。
2) 条带峰值应力曲线
在框架棚洞结构选取顶板中心横向条带(监测点3，8，13，18，23，28，33)和纵梁纵向条带(36，37，38，39，40)，在PFC中通过Fish语言记录、读取并导出条带上各监测点的应力值，从中找出各监测点各应力分量的最大值(正值/负值)。
3) 应力峰值云图
在PFC软件中用fish语言记录、读取、导出所有监测点的应力数值，提取出各应力分量的最大值，在origin软件中导入相应数据绘制框架棚洞的峰值应力云图，进而在应力云图上找到在落石冲击下最大应力发生的部位，即危险部位。
4.1. 监测点应力-时程曲线
选取框架棚洞顶板横向条带几个监测点(3，8，13，18，23，18，33)进行分析和研究。框架棚洞顶板结构横向条带各监测点应力–时程曲线如图5所示。
Figure 5. Stress-time history dynamic curve of each monitoring point in the horizontal plane of the roof of the frame shed tunnel
图5. 框架棚洞顶板横向条带各监测点应力–时程动态变化曲线
从图5中可以得出以下结论：
1) 棚洞结构在落石冲击下，同一部位的各个应力分量变化情况不一致，例如图5(a)中显示，监测点3在落石冲击作用下Y向正应力突变最大，其他应力分量变化较小，说明该部位落石冲击效应以Y向正应力σyy为控制应力，同时为负值，如发生破坏就是受压破坏；
2) 棚洞结构在落石冲击下，不同部位的冲击效应不一致，如图5(a)、图5(c)、图5(k)、图5(m)中显示，顶板靠近端部位置(监测点3、8、28、33)都是Y向正应力σyy突变最大，而图5(e)、图5(g)、图5(i)中显示，顶板靠近中间(监测点13、18、23)落石冲击部位附近，则是τxy表现增大比较明显，说明顶板靠近冲击部位范围内，切应力τxy是控制应力，如果发生破坏则是剪切破坏，即顶板被落石砸穿或者砸破。
3) 回填土对落石冲击具有很大的缓冲作用。分别对比回填土厚土t = 0和t = 1 m各监测点的应力曲线，如图5(a)和图5(b)，有回填土时的最大应力明显小于无回填土，可以说明回填土对落石冲击效应有很大的缓冲作用。
4.2. 棚洞结构峰值应力分布
选取框架棚洞结构顶板结构横向条带几个监测点(3，8，13，18，23，18，33)进行分析和研究。框架棚洞顶板结构横向条带各监测点最大应力曲线见图6所示。
Figure 6. Maximum stress distribution curve of each monitoring point in the horizontal strip of theroof of the shed tunnel
图6. 顶板横向条带各监测点最大应力分布曲线
从图6中可以得出以下结论：
1) 在落石冲击下，棚洞结构不同部位的动力响应不一致。如图6(a)、图6(c)、图6(e)中显示，正应力分量突变最大的是Y向正应力σyy突变最大，最大值达到σyy = 120 MPa，发生在顶板端部位置；而图6(g)、图6(i)、图6(k)中显示；切应力突变最大的则是XY方向切应力τxy，发生在顶板中间靠近落石冲击部位附近，最大值达到τxy = 25 MPa，破坏形式将是剪切破坏，即顶板被落石砸穿或者砸破。该结论与4.1章的结论吻合。
2) 回填土对落石冲击具有很大的缓冲作用。分别对比回填土厚土t = 0和t = 1 m各部位峰值应力分布情况，如图6(a)和图6(b)，有回填土时最大应力明显小于无回填土，说明回填土对落石冲击具有很好的缓冲作用。该结论与4.1章的结论吻合。
4.3. 棚洞结构峰值应力云图
选取框架棚洞顶板结构(监测点1-35)进行分析，顶板的应力峰值云图见图7(a)~(x)。
从PFC3D软件中导出监测点的应力，再将数据导入Origin软件中进行处理，绘制应力云图，目前这种方法还在探索阶段。如何能更高效、准确地绘制离散元颗粒流的应力云图，是今后应该努力探索一个方向。
5. 结论
本文使用PFC3D软件建立落石冲击棚洞模型，模拟落石冲击过程，用Fish语言记录冲击过程中不同部位的应力分量，分别绘制了单个监测点位的应力时程曲线、条带应力峰值分布曲线、顶板结构峰值应力云图，在落石冲击下，研究棚洞结构的动态应力响应。可以得出以下结论：
1) 棚洞结构在落石冲击下，同一部位的正应力和切应力变化趋势不一致；不同部位的动态响应不同；具体分析见本文4.1、4.2。
2) 框架棚洞顶板在落石冲击作用下，落石冲击部位控制应力是XY方向剪切应力τxy，破坏形式是剪切破坏；顶板端部控制应力是Y向正应力σyy，破坏形式是受压破坏；具体分析见本文4.1、4.2。
3) 回填土对落石冲击棚洞结构具有非常明显的缓冲作用，具体分析见本文4.1、4.2；填筑回填土会影响棚洞顶板在落石冲击作用下的应力分布规律，见本文图7。
参考文献