论文部分内容阅读
摘要：目前我国落石冲击力的计算主要有《铁路工程设计技术手册》推荐方法和《公路路基设计规范( JTJ13-95 ) 》推荐方法两种，但由于其计算的是落石的平均冲击力，其计算值比实际值小几倍甚至上百倍，从而导致了实际工程中被动防护结构的失效。本文通过对比与实际冲击力较符合的日本道路公团推荐算法的计算值、现场实测值、颗粒流模拟值，得出结论：通过PFC3D颗粒流计算落石冲击力，比公式法更精确，同时适用性也更广。
关键词：落石；冲击力；颗粒流
中图分类号：O319.56文献标识码：A
引言
珠海主要为剥蚀残丘地貌，再加上雨季台风频繁，因此，边坡的崩塌落石灾害时有发生。落石冲击力是设置被动防护结构的主要依据，国内求解落石冲击力的方法主要有路基规范法、隧道手册法、杨其新方法等，国外计算方法主要有日本道路协会方法、瑞士方法等。前两种方法是利用功能原理和冲量定理获得,杨其新法是根据室内试验进行拟合得到的计算公式，日本法和瑞士法是根据现场试验拟合得到的计算公式。
文献[1]对比了各类落石冲击力算法，发现国内几种算法的计算的仅仅是落石中的平均冲击力，从而导致计算出的落石冲击力严重偏小，而日本法和瑞士是基于现场试验所进行拟合得到的，相对来说更精确。其中特别值得注意的是，规范法的计算结果往往仅是日本法计算结果的几分之一甚至百分之一，从而导致按相关规范所建成防护结构不够安全，常常被落石损坏、从而丧失了保护作用。
从Pichler[2]的现场试验结果来看,与真实落石冲击力符合的最好的算法是日本道路公团法。但对比发现，日本法比实测值仍小30%—60% ,也就是说，它的计算结果也仅仅是可以参考而已，我们需要寻找更精确的算法。
1PFC3D算法、日本算法和Pichler现场试验结果的对比
1.1PFC3D力学原理
作为离散元的一种，三维颗粒流程序(Particle Follow Code PFC3D)的理论基础是Cundall [1979]提出的离散单元法，用微观颗粒及颗粒之间的接触关系，模拟宏观物体的物理性质，PFC3D有ball（颗粒）单元和wall（墙）单元，通过ball之间的组合，可以模拟岩土体的弹性模量，内聚力，内摩擦角，泊松比，单轴抗压强度等等；通过wall之间的组合，可以模拟三维地形及相关地形的力学参数，因此，对落石运动可以进行模拟。
1.1.1颗粒单元基本假设
1)颗粒单元为刚性体;
2)接触发生在很小的范围内，即点接触;
3)接触特性为柔性接触，接触处允许有一定的“重叠”量;
4) “重叠”量的大小与接触力有关，与颗粒大小相比，“重叠”量很小;
5)接触处有特殊的连接强度;
6)颗粒单元为球形。
1.1.2运动方程式
作用在刚性颗粒上的合力与合力矩决定颗粒的运动方式，其表现为质点移动与颗粒转动。质点的移动以位置、速度、及加速度描述; 而颗粒的转动以角速度及角加速度描述。其颗粒的运动方程式如下：
其中：为颗粒所受外力之总和；m为颗粒总质量；为重力加速度; 为颗粒所受的合力矩；为颗粒的转动惯量。
1.1.3接触本构关系
PFC通过每一个接触点的接触模型模拟材料的基本力学行为，其提供三种接触模型：
1）接触刚度模型（Contact-Stiffness Models）：定义了接触力和相对位移的弹性关系，包括线性接触模型（The Linear Contact Model）和简化的Hertz法则模型（The Hertz-Mindlin Contact Model）。
2）滑动模型（The Slip Model）：滑动模型是相互接触球体的一种固有特性，它没有法向抗拉强度，允许颗粒在抗剪强度范围内发生滑动。
3）连結模型（The Bond Model）： 包括两种连接：接触连接模型（The Contact-Bond Model）与平行连接模型（The Parallel-Bond Model）。接触连接认为连接只发生在接触点很小范围内，而平行连接发生在接触颗粒间圆形或方形有限范围内。接触连接只能传递力，而平行连接同时能传递力和力矩。
在落石过程模拟中，两接触物体之间的刚度为我们讨论的重点，我们主要采用接触刚度模型中的线性接触模型，其切向和法向的接触刚度计算公式为：
其中: [ A ] 与 [ B ]表示接触之两物体，
分别代表切线和法向刚度。
1.2日本道路公团推荐算法
日本道路公团基于落石冲击力现场试验的数据和 Hertz弹性碰撞理论[3]导出的落石自由落体情形下的最大冲击力计算公式为
式中
p为落石最大冲击力，单位是KN
M为落石质量，
为拉梅系数，在实验中其范围为1000-10000KN/m2，但论文[3]建议取1000KN/m2
H为落石的自由落高，单位是m。
1.3PFC3D算法、日本算法和Pichler现场试验结果的对比
Pichler在2004年做了相关现场实验，选取近似正方体的花岗岩，岩块的质量从10160kg—18260kg,自由落高从2m—18.85m。
在PFC3D的模拟中，质量为10160kg的岩块采用的参数为：
颗粒质量为10160kg（等效颗粒半径1m，等效密度为2425.5kg/m3），下落高度分别为2m和8.55m。其结果与日本法，数值模拟法对比见表1和图1。
在PFC3D的模拟中，质量为18260kg的岩块采用的参数为：
颗粒质量为18260kg（等效颗粒半径1.2158m，等效密度为2425.5kg/m3），下落高度分别为8.62m、18.85m和18.67m。其结果与日本法，数值模拟法对比见表2和图2。
其结果与日本法，数值模拟法对比如下：
表1 质量为10160kg的岩块冲击力计算结果
Table 1The contrast of rock impact force with 10160kg
图1质量为10160kg的岩块冲击力计算结果对比
Fig. 1The contrast of rock impact force with 10160kg
表2 质量为10820kg的岩块冲击力计算结果
Table 2The contrast of rock impact force with 10820kg
图2质量为18260kg的岩块冲击力计算结果对比
Fig. 2The contrast of rock impact force with 18260kg
小结：对比现场试验实测值，日本法计算值和PFC模拟法，我们发现，PFC模拟计算出的冲击力与实测冲击力最为符合，而PFC模拟计算冲击力时，不仅可以计算自由落体，也可以计算斜坡上落石的抛落，滑动，崩塌等，因此，适用性更广。
2 工程实例
珠海市华发山庄项目[4]位于珠海市香洲区迎宾北路的西侧及翠微东路的北侧，拟开发别墅和7层、11+1层、18+1层洋房。拟建场地北面、西面、南面三面环山，均为剥蚀残丘地貌，由于差异风化，拟建项目周围边坡存在多处孤石抛落，崩塌等地灾隐患，因此，设置主，被动防护装置，计算落石或颗粒流的冲击力就非常有必要。
2.1工况
华发山庄项目中危险边坡最高将近30m，一旦发生崩落，其影响范围和冲击力相当大。此地灾隐患为流动型崩塌，所处坡体主要由燕山三期花岗岩风化土体组成，边坡由人工削坡而成，坡高约25.6m，坡度63～65°，崩后坡度55～60°，崩落方向353°。
表3 边坡形状和危险岩土体体积
Table 3 slope shape and dangerous rock`s volume
表4 全风化花岗岩力学参数
Table 4mechanical parameters of fully weathered granite
2.2 PFC3D建模与计算
2.2.1三维模型的建立
根据上述实测数据，PFC三维建模如下：
边坡：坡高25.6m，坡度60°，坡脚两米处设有被动防护装置
崩滑体：质量17280kg，体积为9m3，由1125个半径为0.1m的小颗粒组成3×3×1的岩土块体。
图3PFC3D建模示意图
Figure 3 PFC3D modeling diagram
2.2.2力学参数的标定
参考华发山庄地灾评估报告和文献[5],可知该边坡全风化花岗岩的力学参数如下
表5 全风化花岗岩实测力学性质与PFC3D模拟值
table 5 The measured values and PFC3D simulation values
上述表格对比列出了全风化花岗岩的实测和数值模拟的宏观性质，对于PFC3D来说，岩土体的宏观性质跟颗粒之间的微观性质相关，现将模拟全风化花岗岩所需的微观参数罗列如下：
表6PFC3D对应全风化花岗岩微观参数表
Table 6 PFC3D microscopic parameters of weathered granite
2.2.3计算结果
1）崩塌体对平地冲击力分析
崩塌体刚撞击地面时，地面所受的冲击力最大，为3449KN,而横向冲击在竖向碰撞中并不占主导地位，同时我们也可以从图中得出，崩塌体从开始到第一次冲击地面总共用时3s，如下图所示
图4 崩塌体对水平地面的冲击
（蓝色为纵向冲击力，红色和黑色为横向冲击力）
Figure 4 the impact of collapse body on level ground
2）崩塌体对竖向被动防护结构的冲击力分析
对于崩塌体来说，撞击地面并不是灾害的结束，它将继续前行，并对建筑物或被动防护装置产生向前的冲力，如下图所示，被动防护墙所受最大冲击力为1495KN,约为第一次碰撞冲击力的一半。同时，对比图4和图5，我们可以看出，当崩塌体开始撞击时，由于碰撞后岩土体的破碎，冲击力比完整岩体冲击力小，但其随时间急剧变化，形成冲击波，在考虑建筑物的抗震性能时，冲击波式的瞬间多次冲击也是我们需要考虑的，而PFC的模拟计算，提供了此种可能。
图5 崩塌体对竖向建筑物的的冲击
（黑色为左右冲击力，蓝色为纵向，红色为横向冲击力）
Figure 5 the impact of collapse body on the vertical structures
3 结语
通过对比落石冲击力的几种计算方法，得出用PFC3D模拟法计算的落石沖击力与现场试验值符合的最好，并且有以下几个优点：
①PFC3D模拟法计算冲击力适用于任何地形：
目前流行的冲击力计算方法多基于自由落体运动，因此普适性不广，而PFC3D模拟法，可以计算任意地形下，落石或崩滑体对坡脚建筑的冲击力，适用性更广。
②PFC3D模拟法可以直观的得出落石或崩塌体的致灾范围和发生时间。
③PFC3D模拟法可记录冲击力随时间变化的曲线，为建筑物的动态防护提供一定的依据。
总之，引入数值模拟技术，对工程计算提供了新的可行性思路。
参考文献(References）
[1]叶四桥,陈洪凯,唐红梅. 落石冲击力计算方法的比较研究[J]. 水文地质工程地质,2010,02:59-64.叶四桥等 落石冲击力计算方法的比较研究
[2]B. Pichler;Ch. Hellmich;H.A. Mang.Impact of rocks onto gravel Design and evaluation of experiments[J].International Journal of Impact Engineering.2005(NO.5).
[3]Kawahara, S.;Muro, T..Effects of dry density and thickness of sandy soil on impact response due to rockfall.[J].Journal of Terramechanics.2006(No.3).
[4]曾新雄,龙明滔,陈开明等 广东省珠海市华发山庄边坡地质灾害危险性评估报告[R]2010
[5]于洪丹,陈卫忠,郭小红,贾善坡,曹俊杰,梁巍. 厦门海底隧道强风化花岗岩力学特性研究[J]. 岩石力学与工程学报,2010,02:381-387.
关键词：落石；冲击力；颗粒流
中图分类号：O319.56文献标识码：A
引言
珠海主要为剥蚀残丘地貌，再加上雨季台风频繁，因此，边坡的崩塌落石灾害时有发生。落石冲击力是设置被动防护结构的主要依据，国内求解落石冲击力的方法主要有路基规范法、隧道手册法、杨其新方法等，国外计算方法主要有日本道路协会方法、瑞士方法等。前两种方法是利用功能原理和冲量定理获得,杨其新法是根据室内试验进行拟合得到的计算公式，日本法和瑞士法是根据现场试验拟合得到的计算公式。
文献[1]对比了各类落石冲击力算法，发现国内几种算法的计算的仅仅是落石中的平均冲击力，从而导致计算出的落石冲击力严重偏小，而日本法和瑞士是基于现场试验所进行拟合得到的，相对来说更精确。其中特别值得注意的是，规范法的计算结果往往仅是日本法计算结果的几分之一甚至百分之一，从而导致按相关规范所建成防护结构不够安全，常常被落石损坏、从而丧失了保护作用。
从Pichler[2]的现场试验结果来看,与真实落石冲击力符合的最好的算法是日本道路公团法。但对比发现，日本法比实测值仍小30%—60% ,也就是说，它的计算结果也仅仅是可以参考而已，我们需要寻找更精确的算法。
1PFC3D算法、日本算法和Pichler现场试验结果的对比
1.1PFC3D力学原理
作为离散元的一种，三维颗粒流程序(Particle Follow Code PFC3D)的理论基础是Cundall [1979]提出的离散单元法，用微观颗粒及颗粒之间的接触关系，模拟宏观物体的物理性质，PFC3D有ball（颗粒）单元和wall（墙）单元，通过ball之间的组合，可以模拟岩土体的弹性模量，内聚力，内摩擦角，泊松比，单轴抗压强度等等；通过wall之间的组合，可以模拟三维地形及相关地形的力学参数，因此，对落石运动可以进行模拟。
1.1.1颗粒单元基本假设
1)颗粒单元为刚性体;
2)接触发生在很小的范围内，即点接触;
3)接触特性为柔性接触，接触处允许有一定的“重叠”量;
4) “重叠”量的大小与接触力有关，与颗粒大小相比，“重叠”量很小;
5)接触处有特殊的连接强度;
6)颗粒单元为球形。
1.1.2运动方程式
作用在刚性颗粒上的合力与合力矩决定颗粒的运动方式，其表现为质点移动与颗粒转动。质点的移动以位置、速度、及加速度描述; 而颗粒的转动以角速度及角加速度描述。其颗粒的运动方程式如下：
其中：为颗粒所受外力之总和；m为颗粒总质量；为重力加速度; 为颗粒所受的合力矩；为颗粒的转动惯量。
1.1.3接触本构关系
PFC通过每一个接触点的接触模型模拟材料的基本力学行为，其提供三种接触模型：
1）接触刚度模型（Contact-Stiffness Models）：定义了接触力和相对位移的弹性关系，包括线性接触模型（The Linear Contact Model）和简化的Hertz法则模型（The Hertz-Mindlin Contact Model）。
2）滑动模型（The Slip Model）：滑动模型是相互接触球体的一种固有特性，它没有法向抗拉强度，允许颗粒在抗剪强度范围内发生滑动。
3）连結模型（The Bond Model）： 包括两种连接：接触连接模型（The Contact-Bond Model）与平行连接模型（The Parallel-Bond Model）。接触连接认为连接只发生在接触点很小范围内，而平行连接发生在接触颗粒间圆形或方形有限范围内。接触连接只能传递力，而平行连接同时能传递力和力矩。
在落石过程模拟中，两接触物体之间的刚度为我们讨论的重点，我们主要采用接触刚度模型中的线性接触模型，其切向和法向的接触刚度计算公式为：
其中: [ A ] 与 [ B ]表示接触之两物体，
分别代表切线和法向刚度。
1.2日本道路公团推荐算法
日本道路公团基于落石冲击力现场试验的数据和 Hertz弹性碰撞理论[3]导出的落石自由落体情形下的最大冲击力计算公式为
式中
p为落石最大冲击力，单位是KN
M为落石质量，
为拉梅系数，在实验中其范围为1000-10000KN/m2，但论文[3]建议取1000KN/m2
H为落石的自由落高，单位是m。
1.3PFC3D算法、日本算法和Pichler现场试验结果的对比
Pichler在2004年做了相关现场实验，选取近似正方体的花岗岩，岩块的质量从10160kg—18260kg,自由落高从2m—18.85m。
在PFC3D的模拟中，质量为10160kg的岩块采用的参数为：
颗粒质量为10160kg（等效颗粒半径1m，等效密度为2425.5kg/m3），下落高度分别为2m和8.55m。其结果与日本法，数值模拟法对比见表1和图1。
在PFC3D的模拟中，质量为18260kg的岩块采用的参数为：
颗粒质量为18260kg（等效颗粒半径1.2158m，等效密度为2425.5kg/m3），下落高度分别为8.62m、18.85m和18.67m。其结果与日本法，数值模拟法对比见表2和图2。
其结果与日本法，数值模拟法对比如下：
表1 质量为10160kg的岩块冲击力计算结果
Table 1The contrast of rock impact force with 10160kg
图1质量为10160kg的岩块冲击力计算结果对比
Fig. 1The contrast of rock impact force with 10160kg
表2 质量为10820kg的岩块冲击力计算结果
Table 2The contrast of rock impact force with 10820kg
图2质量为18260kg的岩块冲击力计算结果对比
Fig. 2The contrast of rock impact force with 18260kg
小结：对比现场试验实测值，日本法计算值和PFC模拟法，我们发现，PFC模拟计算出的冲击力与实测冲击力最为符合，而PFC模拟计算冲击力时，不仅可以计算自由落体，也可以计算斜坡上落石的抛落，滑动，崩塌等，因此，适用性更广。
2 工程实例
珠海市华发山庄项目[4]位于珠海市香洲区迎宾北路的西侧及翠微东路的北侧，拟开发别墅和7层、11+1层、18+1层洋房。拟建场地北面、西面、南面三面环山，均为剥蚀残丘地貌，由于差异风化，拟建项目周围边坡存在多处孤石抛落，崩塌等地灾隐患，因此，设置主，被动防护装置，计算落石或颗粒流的冲击力就非常有必要。
2.1工况
华发山庄项目中危险边坡最高将近30m，一旦发生崩落，其影响范围和冲击力相当大。此地灾隐患为流动型崩塌，所处坡体主要由燕山三期花岗岩风化土体组成，边坡由人工削坡而成，坡高约25.6m，坡度63～65°，崩后坡度55～60°，崩落方向353°。
表3 边坡形状和危险岩土体体积
Table 3 slope shape and dangerous rock`s volume
表4 全风化花岗岩力学参数
Table 4mechanical parameters of fully weathered granite
2.2 PFC3D建模与计算
2.2.1三维模型的建立
根据上述实测数据，PFC三维建模如下：
边坡：坡高25.6m，坡度60°，坡脚两米处设有被动防护装置
崩滑体：质量17280kg，体积为9m3，由1125个半径为0.1m的小颗粒组成3×3×1的岩土块体。
图3PFC3D建模示意图
Figure 3 PFC3D modeling diagram
2.2.2力学参数的标定
参考华发山庄地灾评估报告和文献[5],可知该边坡全风化花岗岩的力学参数如下
表5 全风化花岗岩实测力学性质与PFC3D模拟值
table 5 The measured values and PFC3D simulation values
上述表格对比列出了全风化花岗岩的实测和数值模拟的宏观性质，对于PFC3D来说，岩土体的宏观性质跟颗粒之间的微观性质相关，现将模拟全风化花岗岩所需的微观参数罗列如下：
表6PFC3D对应全风化花岗岩微观参数表
Table 6 PFC3D microscopic parameters of weathered granite
2.2.3计算结果
1）崩塌体对平地冲击力分析
崩塌体刚撞击地面时，地面所受的冲击力最大，为3449KN,而横向冲击在竖向碰撞中并不占主导地位，同时我们也可以从图中得出，崩塌体从开始到第一次冲击地面总共用时3s，如下图所示
图4 崩塌体对水平地面的冲击
（蓝色为纵向冲击力，红色和黑色为横向冲击力）
Figure 4 the impact of collapse body on level ground
2）崩塌体对竖向被动防护结构的冲击力分析
对于崩塌体来说，撞击地面并不是灾害的结束，它将继续前行，并对建筑物或被动防护装置产生向前的冲力，如下图所示，被动防护墙所受最大冲击力为1495KN,约为第一次碰撞冲击力的一半。同时，对比图4和图5，我们可以看出，当崩塌体开始撞击时，由于碰撞后岩土体的破碎，冲击力比完整岩体冲击力小，但其随时间急剧变化，形成冲击波，在考虑建筑物的抗震性能时，冲击波式的瞬间多次冲击也是我们需要考虑的，而PFC的模拟计算，提供了此种可能。
图5 崩塌体对竖向建筑物的的冲击
（黑色为左右冲击力，蓝色为纵向，红色为横向冲击力）
Figure 5 the impact of collapse body on the vertical structures
3 结语
通过对比落石冲击力的几种计算方法，得出用PFC3D模拟法计算的落石沖击力与现场试验值符合的最好，并且有以下几个优点：
①PFC3D模拟法计算冲击力适用于任何地形：
目前流行的冲击力计算方法多基于自由落体运动，因此普适性不广，而PFC3D模拟法，可以计算任意地形下，落石或崩滑体对坡脚建筑的冲击力，适用性更广。
②PFC3D模拟法可以直观的得出落石或崩塌体的致灾范围和发生时间。
③PFC3D模拟法可记录冲击力随时间变化的曲线，为建筑物的动态防护提供一定的依据。
总之，引入数值模拟技术，对工程计算提供了新的可行性思路。
参考文献(References）
[1]叶四桥,陈洪凯,唐红梅. 落石冲击力计算方法的比较研究[J]. 水文地质工程地质,2010,02:59-64.叶四桥等 落石冲击力计算方法的比较研究
[2]B. Pichler;Ch. Hellmich;H.A. Mang.Impact of rocks onto gravel Design and evaluation of experiments[J].International Journal of Impact Engineering.2005(NO.5).
[3]Kawahara, S.;Muro, T..Effects of dry density and thickness of sandy soil on impact response due to rockfall.[J].Journal of Terramechanics.2006(No.3).
[4]曾新雄,龙明滔,陈开明等 广东省珠海市华发山庄边坡地质灾害危险性评估报告[R]2010
[5]于洪丹,陈卫忠,郭小红,贾善坡,曹俊杰,梁巍. 厦门海底隧道强风化花岗岩力学特性研究[J]. 岩石力学与工程学报,2010,02:381-387.