弹性力学论文：石灰岩细观力学特性的颗粒流模拟

第 31 卷增刊 2 岩 土 力 学 Vol.31 Supp.2 2010 年 11 月 Rock and Soil Mechanics Nov. 2010 收稿日期：2010-07-30 基金项目：国家自然科学基金资助项目(No. 40972191)；上海市教育委员会科研创新项目(No. 09YZ39)。 第一作者简介：徐金明，男，1963 年生，博士、教授、博士生导师，主要从事岩土和工程地质计算技术的教学和科研工作。Email: xjming@163.com 文章编号：1000－7598 (2010)增刊 2－0390－06 石灰岩细观力学特性的颗粒流模拟 徐金明 1，谢芝蕾 1，贾海涛 2 （1. 上海大学 土木工程系，上海 200072；2. 上海自然博物馆工程建设指挥部，上海 200041） 摘 要：岩体地区地质灾害的发生和发展取决于岩石细观组分的运动学行为。研究岩石运动学行为时通常将岩石作为整体研 究对象较多，而直接以细观组分为对象的研究较少。以石灰岩为例，根据室内试验获得的岩石力学性质指标，使用基于非连 续介质理论的颗粒流方法，将材料离散成刚性颗粒组成的模型，把颗粒细观变化与宏观力学特性联系起来，建立了石灰岩的 细观结构模型，获得了颗粒接触力、颗粒接触模量、接触连接强度和连接刚度比等细观力学参数。由于文中直接以细观成分 为研究对象、反映了岩石和岩体组成的本质特点，所得结论不仅对含裂隙岩石本构关系研究具有广阔的应用前景，而且对岩 体工程性质和地质灾害机制研究也具有重要的理论意义。 关 键 词：石灰岩；细观力学特性；颗粒流；模拟 中图分类号：TU 452 文献标识码：A Simulation of mesomechanical properties of limestone using particle flow code XU Jin-ming1，XIE Zhi-lei1，JIA Hai-tao2 ( 1. Department of Civil Engineering，Shanghai University, Shanghai 200072，China; 2. Shanghai Science and Technology Museum, Shanghai 200041，China) Abstract: The formation and development of geological disasters in rock area are dependant on the kinematic behaviors of rocks, especially of grains, fissures, and fillings in the rocks. In the conventional studies, rocks are generally treated as entireties and few concerns are concentrated on the individual meso-compositions in these rocks. Taking a limestone for example, macromechanical properties were obtained for the rock specimens of laboratory tests; and particle flow code in two-dimensions (PFC2D) was used for simulating the macromechanical properties of the rock material. In the simulation, the material was discretized as an assembly of rigid particles. The mesomechanical parameters, such as contact forces, contact modulus, normal contact strengths, and stiffness ratio, were obtained; and the mesostructural model was established for the limestone; connecting meso-level changes in particles with macromechanical properties. Because the individual compositions were taken as the direct objectives, reflecting the intrinsic features of rock materials or rock masses, the techniques presented herein may be of great significance in studying the constitutive law of fissured rocks, engineering properties of rock masses, and mechanism of geological disasters. Key words: limestone；mesomechanical property；particle flow；simulation 1 引 言 岩体地区地质灾害的发生和发展取决于岩石 的运动学行为、尤其是岩石中颗粒、裂隙、充填物 等细观组分的变化情况，常规宏观分析方法以岩石 整体为研究对象较多，直接以细观组分为对象进行 研究较少。基于非连续介质理论的颗粒流方法，将 材料离散成刚性颗粒组成的模型，把颗粒细观力学 参数与宏观力学特性联系起来，可以用于模拟颗粒 之间的相互作用和破裂面的形成扩展过程。 使用颗粒流方法对土的细观力学行为进行细 观模拟，多使用 PFC2D (particle flow code in 2-dimensions) 或 PFC3D (particle flow code in 3-dimensions)。周健[1]研究团队在这方面做了大量工 作（比如，模拟了不同水压下渗流引起砂土特性变 化的全过程）。 使用颗粒流方法进行岩石力学特性的细观模 拟也有一些报道。Potyondy[2]、Potyondy 和 Cundall[3] 增刊 2 徐金明等：石灰岩细观力学特性的颗粒流模拟 391 获得了 Lac du Bonnet(LDB)花岗岩的细观参数； Backstrom 等[4]研究了花岗岩细观模拟结果与宏观 应力-应变曲线的关系；Hsieh 等[5]研究了砂岩变形 机制，分析了细观组分含量与宏观力学特性之间的 关系；Fakhimi 等[6]用弹性模量、泊松比和单轴抗压 强度来确定砂岩的细观力学参数；Doležalová 等[7] 对比研究了碎裂状泥岩的细观模拟、室内试验、现 场试验结果；Stavropoulou[8]对大理岩芯样进行了细 观试验研究；徐文杰等[9]结合数字图像处理技术， 建立了土石混合体细观结构的细观数值计算模型。 此外，李凡[10]采用胶结在一起的小圆盘代替可 破碎颗粒并施加适当胶结，对岩土材料细观破损机 制进行了研究；Cho 等[11]认为，使用簇粒可以明显 改善岩石力学特性的预测能力；Cho 等[12]模拟了脆 性岩石材料的直剪试验过程；Kulatilake 等[13]研究 了节理几何参数(组数、密度、方向、大小、分布) 对含裂隙岩石单轴抗压强度的影响；Schopfer 等[14] 模拟了低围压高强度脆韧性层状岩层中正断层的发 展过程；王涛等[15]对矿岩自然崩落现象进行了细观 数值模拟。 本文拟以石灰岩为例，根据室内试验获得的岩 石力学性质指标，建立岩石的细观结构模型，使用 颗粒流数值模拟程序 PFC2D 获得细观力学参数。 本文方法对含裂隙岩石本构关系和岩体工程性质的 细观研究具有重要的参考价值。 2 细观力学参数的确定 2.1 原始数据 本文使用杨仕教等[16]的资料确定石灰岩的细 观力学参数，所用的主要试验数据见表 1，模拟时 岩石尺寸取为 42 mm×100 mm。 表 1 石灰岩单轴抗压强度试验数据(引自文献[16]) Table 1 Mechanical parameters of limestone under axial loading tests (from ref. [16]) 组号 弹性模量/GPa 峰值应力/MPa 泊松比 1 39.54 86.963 0.248 2 45.10 97.030 0.235 3 45.82 103.330 0.319 4 46.49 116.880 0.265 平均值 44.24 101.050 0.267 2.2 设置几何模型与物理模型 研究时，采用单位厚度的圆形颗粒来表示各组 分分布，由光滑无摩擦的“墙”来实现模型边界条件， 以石灰岩单轴抗压强度试验中弹性模量、峰值强度、 泊松比的平均值(见表 1)作为模拟目标。 PFC2D 中，颗粒接触模型分为基本模型和混合 模型两种。基本模型有刚度接触、平行接触和非连 接接触 3 类。由于平行连接可以传递颗粒之间的力 和力矩、平行连接强度和岩石材料强度有较好的对 应关系，颗粒之间的胶结物用平行连接来表征。 2.3 确定细观力学参数 岩石宏观参数不仅受到岩石中颗粒大小、形状 和分布的影响，还与颗粒和胶结物的变形和强度特 性有关。岩石细观力学参数的确定过程，就是根据 宏观参数(弹性模量、峰值应力、泊松比)不断调整 细观参数使模拟结果与试验结果误差达到最小的过 程。这一过程可以表述为： (1) 将平行连接法向连接强度σ 和切向连接强 度τ 设为一个较大的值(比如峰值应力)，将颗粒-颗 粒接触模量 cE 和平行连接弹性模量 cE 的比值 cE / cE 设置为 1.0，通过改变 cE 大小来选择岩石的 弹性模量； (2) 对颗粒法向刚度 kn 和切向刚度 ks、平行连 接法向刚度 nk 和切向刚度 sk ，将比值(kn/ks)/( n s/k k ) 和 n n/k k 设置为 1.0，通过改变 kn/ks来选择泊松比； (3) 将侧限压力 'xσ 设置为 0、将平行连接法向 切向强度均值比值(σ /τ )、法向强度差 sσ 和切 向强度标准差 sτ 的比值( sσ / sτ )设置为 1，将比值 ( sσ /σ )设置为一个比较小的数值(比如 0.12)，通过 改变σ 来选择峰值应力； (4) 经过步骤(1)～(3)，模拟效果通常已经比较 好；为了达到预设的误差精度，可以通过适当改变 上述(1)～(3)中的相关数值来实现(比如改变 cE 、kn、 ks、 nk 、 sk 、τ 、 sσ 、 sτ )。 初始细观力学参数在数值模拟过程中具有极 为重要的作用，模拟开始前必须首先设置。由于 Itasca 公司提供了一套算法并已通过程序代码包 Augmented FishTank 得以实现，下面将根据 Potyondy[2]、Potyondy和Cundall[3]、Cho等[11]与 Itasca 公司[17]的研究结果，说明算法中主要初始细观参数 的设置方法(取拉应力为正)。 （1）试验条件 试验条件主要通过设置变量 et2_ucs 来实现。 对于无侧限压缩试验和侧限压缩试验，这一变量分 别设为 1 和 0。 （2）试样生成控制参数 y 试样宽度 w：单位为 m，即 et2_xlen 的设 置值，根据实际试样宽度进行设置；研究 时取 42×10-3 m。 y 试样高度 h：单位为 m，即变量 et2_ylen 392 岩 土 力 学 2010 年 的设置值，根据实际试样高度进行设置； 研究时取 100×10-3 m。 y 最小颗粒半径 Rmin：单位为 m，即 et2_rlo 的设置值；研究时设为 0.275×10-3m。 y 粒径比 Rmax/Rmin：即 et2_radius_ratio 的设 置值，研究时设为 1.66。 y 墙刚度乘子 β：即 md_wEcfac 的设置值， 研究时设为 1.1。 y 内锁等向应力墙刚度乘子 σ0：单位为 Pa， 即 tm_req_isostr 的设置值，可取小于岩石 单轴抗压强度的 1%，研究时取-1.0×10-6 Pa。 y 内锁等向应力墙刚度乘子误差容许值 σ0t： 即 tm_req_isostr_tol 的设置值，研究时取 0.50。 y 最小非悬浮接触数 Nfmin：即 flt_def 的设置 值，研究时取 3。 y 剩余悬浮接触比 nfN：即 flt_remain 的设置 值，研究时设为 0.0。 （3）颗粒连接设置参数 y 颗粒密度 ρ：单位为 kg/m3，即 md_dens 的设置值，研究时取 2 630.0 kg/m3。 y 颗粒-颗粒接触模量 Ec：单位为 Pa，即 md_Ec 的设置值，研究时取 37×109 Pa。 y 颗粒刚度比 kn/ks：即 md_knoverks 的设置 值，研究时设为 2.46。 y 颗粒摩擦系数 μ：即 md_fric 的设置值，研 究时设为 0.5。 y 平行连接法向强度均值σ ：单位为 Pa，即 pb_sn_mean 的设置值、可用于模拟胶结物 的法向强度，研究时设为 77×106 Pa。 y 平行连接法向强度标准差 sσ ：单位为 Pa， 即 pb_sn_sdev 的设置值，研究时设为 10×106 Pa。 y 平行连接切向强度均值τ ：单位为 Pa，即 pb_ss_mean 的设置值、可用于模拟胶结物 的切向强度，研究时设为 77×106 Pa。 y 平行连接切向强度标准差 sτ ：单位为 Pa， 即 pb_ss_sdev 的设置值，研究时设为 10×106 Pa。 y 平行连接半径乘子λ ：即 pb_radmult 的设 置值、降低该值可以降低胶结物的数量与 连接刚度，研究时取 1.0。 y 平行连接弹性模量 cE ：单位为 Pa，即 pb_Ec 的设置值，研究时取 37×109 Pa。 y 平行连接刚度比 n s/k k ：即 pb_ knoverks 的设置值，研究时设为 2.46。 y 平行连接设置入口：即是否平行连接，将 md_add_pbonds 设为 1(平行连接)。 y 保存数据设置入口：即是否保存全部数据 设置，将 et2_prep_saveall 设为 1(保存)。 （4）模拟控制参数 y 侧墙刚度折减因子 βx：即 et2_knxfac 的设 置值，对于无侧限压缩试验、可不做设置； 对于侧限压缩试验，研究时设为 1.0；对 于巴西试验，设为 0.001。 y 竖墙刚度折减因子 βy：即 et2_knyfac 的设 置值，研究时设为 1.0。 y 目标侧限应力 'xσ ：单位为 Pa ，即 et2_wsxx_req 的设置值，根据实际侧限压 力进行设置；比如，实际侧限压力为 1MPa 时设置为-1.0×106 Pa。 y 目标竖向应力 'yσ ：单位为 Pa ，即 et2_wsyy_req 的设置值，研究时设为 -1.0×105 Pa。 y 墙伺服误差 ε：即 et2_ws_tol 的设置值， 研究时设为 0.01。 y 平台速度终值 ε：单位为 m/s，即 p_vel 的 设置值，研究时设为 5.0×10-2m/s。 y 平台加速度总循环数 Np：p_cyc 的设置值， 研究时设为 400。 y 平台加速多段数 Sp：p_stages 的设置值， 研究时设为 10。 y 初始裂纹标准 ciσ ：pk_ci_fac 的设置值、 即出现初始裂纹的应力与峰值应力的百 分比，研究时设为 0.10。 y 试验终止标准α ：et2_peakfac 的设置值， 研究时设为 0.80。 3 石灰岩离散元模拟结果 3.1 石灰岩细观参数 使用上面方法得到的石灰岩模型的主要细观 参数见表 2，宏观参数模拟结果见表 3。表 2 和表 3 中，抗拉强度 tσ 通过巴西试验获得，抗剪强度指标 c 和ϕ 通过不同侧限压力下的双轴试验获得。在具 体计算强度指标 c 和ϕ时，首先根据最小二乘法获 得峰值应力包线、然后根据这一包线的截距和斜率 计算指标 c和ϕ。 3.2 结果分析 下面说明使用石灰岩细观参数所得的一些颗 增刊 2 徐金明等：石灰岩细观力学特性的颗粒流模拟 393 粒流模拟结果。图 1 是石灰岩细观模拟时的颗粒几 何分布和平行连接分布。图 2 是石灰岩位移分布和 颗粒接触力分布。 表 2 石灰岩模型细观参数 Table 2 Model microproperties to match limestone 材料 变量参数 数值 试样 宽度/m 42×10 -3 高度/m 100×10-3 颗粒 最小半径 /m 0.275×10-3 粒径比 1.66 密度/(g/cm3) 2.63 接触模量/GPa 37 刚度比 2.46 摩擦系数 0.50 胶结物(平行连接) 法向强度均值/MPa 77 法向强度标准差/MPa 10 刚度比 2.46 切向强度均值/MPa 77 切向强度标准差/MPa 10 弹性模量 37 半径乘子 1.0 表 3 石灰岩宏观参数模拟结果 Table 3 Macroproperties of the model of limestone 颗粒 总数 弹性模 量/GPa 泊松比 侧限应力 /MPa 峰值应力 / MPa 抗拉强度 / MPa 凝聚力 / MPa 内摩擦 角/(º) 8 394 43.88 0.235 9 0.0 102.2 22.5 31.7 28.0 1.0 108.2 2.0 110.4 3.0 114.0 4.0 116.8 5.0 119.2 6.0 122.8 7.0 125.8 8.0 127.7 9.0 130.6 10.0 132.0 (a) 颗粒分布 (b) 平行连接分布 图 1 石灰岩细观模拟时的颗粒几何分布和平行连接分布 Fig.1 Grain distribution and parallel contacts of limestone in micro simulation. 关于细观参数对宏观参数的影响，下面作一简 要分析。分析时，宏观参数取弹性模量、峰值应力 和泊松比；细观参数取颗粒-颗粒接触模量 cE 、平 行连接法向连接强度σ 、平行连接法向刚度 nk ；一 个细观参数变化时，其它细观参数取表 2 中的数值。 (a) 位移分布 (b) 颗粒接触力分布 图 2 石灰岩细观模拟时的位移分布和颗粒接触力分布 Fig.2 Diasplacement and contact force distributions of limestone in meso-simulation. 表 4 为不同样品尺寸下的宏观参数。由表 4 可 以看出，试样宽度不变时，试样高度增加，弹性模 量、峰值应力和泊松比变化不大（基本在 10%以内 变化），变化规律性也不强。 表 4 不同样品尺寸下的宏观参数 Table 4 Macro properties for different specimen sizes 试样宽度 /mm 试样高度 /mm 弹性模量 /GPa 峰值应力 / MPa 泊松比 42 50 43.55 103.8 0.244 2 60 43.69 109.8 0.228 1 70 43.87 108.0 0.246 7 80 44.16 99.7 0.236 1 90 44.27 110.9 0.237 4 100 43.62 102.2 0.237 3 表 5 为颗粒-颗粒接触模量不同时的宏观参数。 由表 5 可以看出，颗粒-颗粒接触模量增加 62%（从 37.0 增加到 60.0 GPa），弹性模量相应增加 22%， 泊松比增加 15.0%，但峰值应力变化的规律性不强。 表 5 不同颗粒-颗粒接触模量下的宏观参数 Table 5 Macro properties for different particle-particle contact moduli 颗粒-颗粒接触模量 /GPa 弹性模量 /GPa 峰值应力 / MPa 泊松比 37.0 43.62 102.22 0.232 8 40.0 45.00 105.24 0.236 2 45.0 46.86 99.39 0.248 5 50.0 51.60 110.36 0.297 2 55.0 51.39 108.94 0.254 7 60.0 53.31 93.77 0.267 6 394 岩 土 力 学 2010 年 表 6 为平行连接法向连接强度不同时的宏观参 数。由表 6 可以看出，平行连接法向连接强度增加 30%（从 77 增加到 100 MPa），弹性模量和泊松比 变化不大，但峰值应力增加了 28%。 表 6 平行连接法向连接强度不同时的宏观参数 Table 6 Macro properties for different parallel normal contact strengths 平行连接法向连接强度 /MPa 弹性模量 /GPa 峰值应力 / MPa 泊松比 77 43.62 102.22 0.232 8 80 43.56 105.74 0.232 9 85 43.63 112.69 0.233 2 90 43.64 120.19 0.233 4 95 43.65 123.26 0.233 5 100 43.66 130.71 0.233 7 表 7 为平行连接刚度比不同时的宏观参数。由 表 7 可以看出，平行连接刚度比增加 50%（从 2.0 增加到 3.0），弹性模量降低 7.4%，峰值应力增加 2.8%，泊松比则增加 17.3%。 表 7 平行连接刚度比不同时的宏观参数 Table 7 Macro properties for different parallel stiffness ratios 平行连接法向刚度强度 /GPa 弹性模量 /GPa 峰值应力 / MPa 泊松比 2.0 44.26 99.77 0.213 3 2.2 45.14 98.43 0.222 3 2.4 43.95 99.75 0.253 4 2.6 42.87 102.55 0.237 8 2.8 41.89 102.19 0.244 2 3.0 40.99 102.60 0.250 2 4 讨 论 必须指出的是，在上述处理过程中，调整颗粒 细观参数花费时间很多，细观力学参数会受到细观 物理参数的影响，宏观力学性质参数也受到多方面 因素的影响，下面作一简要讨论。 （1）颗粒采用单位厚度的圆形颗粒来表示， 与实际情况不符，会导致内摩擦角降低(见文献[11] 和[17]的讨论)。为此，文献[17]建议采用可破碎颗 粒 cluster 来表征颗粒，而文献[11]建议采用不可破 碎颗粒 clump 来表征颗粒。到底采用哪一类颗粒来 模拟岩石颗粒，目前还没有定论。 （2）本次研究时，颗粒之间胶结采用平行接 触来模拟。实际岩石的宏观力学性质参数与胶结物 成分和胶结类型有一定关系。如何在细观结构模型 中准确地反映这一关系，也还需要进一步研究。 （3）本文采用巴西试验估计抗拉强度 tσ 。但 是，正如 Cho 等[11]所指出，模拟所得强度值比实际 抗拉强度要高出很多、即使采用不可破碎颗粒 clump也还是比实际值要大(很多文献都碰到类似问 题，比如文献[11])。作者认为，要解决这一难题， 可以在考虑颗粒接触时同时使用平行连接和刚度接 触、在模拟时考虑颗粒之间的转动问题。 （4）使用 PFC2D 进行岩石细观模拟时，确定 岩石细观物理力学参数比较麻烦、但又是必须进行 的一项工作。作者认为，可以根据已有宏观试验数 据及岩石和胶结物的不同特点，分别建立细观参数 和宏观参数的多对多关系，这对尽快、准确进行不 同岩石的细观模拟必将具有很好的参考价值。作者 正准备进行这一方面的研究。 5 结 语 （1）以石灰岩为例，根据室内试验获得的岩 石力学性质指标，使用基于非连续介质理论的颗粒 流方法，建立了石灰岩的细观结构模型。 （2）本文将石灰岩颗粒之间的连接使用平行 连接来表征，获得了石灰岩颗粒之间的连接强度、 刚度大小、弹性常数及其与宏观力学性质参数之间 的关系，同时对石灰岩模拟时初始细观参数的选取 方法进行了讨论。 （3）本文以石灰岩细观成分为研究对象、反 映了岩石和岩体组成的本质特点，对含裂隙岩石本 构关系和地质灾害机制研究具有重要的参考价值。 参 考 文 献 [1] 周健，姚志雄，张刚. 砂土渗流过程的细观数值模拟[J]. 岩土工程学报，2007，29(7)：977－981. ZHOU Jian, YAO Zhi-xiong, ZHANG Gang. Mesomechanical simulation of seepage flow in sandy soil[J]. Chinese Journal of Geotechnical Engineering, 2007, 29(7): 977－981. [2] POTYONDY D. A bonded-disk model for rock: Relating microproperties and macroproperties[C]// Discrete Element Methods: Numerical Modeling of Discontinua. Proceedings of the Third International Conference, Santa Fe: [s. n.], 2002. [3] POTYONDY D O, CUNDALL P A. A bonded-particle model for rock[J]. International Journal of Rock Mechanics and Mining Sciences, 2004, (41): 1329－ 增刊 2 徐金明等：石灰岩细观力学特性的颗粒流模拟 395 1364. [4] BACKSTROM A, ANTIKAINEN J, BACKERS T, et al. Numerical modelling of uniaxial compressive failure of granite with and without saline porewater[J]. International Journal of Rock Mechanics and Mining Sciences, 2008, (45): 1126－1142. [5] HSIEH Y M, LI H H, HUANG T H, et al. Interpretations on how the macroscopic mechanical behavior of sandstone affected by microscopic properties—Revealed by bonded-particle model[J]. Engineering Geology, 2008, (99): 1－10. [6] FAKHIMI A, CARVALHO F, ISHIDA T, et al. Simulation of failure around a circular opening in rock[J]. International Journal of Rock Mechanics and Mining Sciences, 2002, (39): 507－515. [7] DOLEŽALOVÁ M, KORÁN P. Micromechanical modeling of a dump material[J]. International Journal of Geomechanics, 2002, 2(1): 47－74. [8] STAVROPOULOU M. Modeling of small-diameter rotary drilling tests on marbles[J]. International Journal of Rock Mechanics and Mining Sciences, 2006, (43): 1034 －1051. [9] 徐文杰，胡瑞林，王艳萍. 基于数字图像的非均质岩土 材料细观结构 PFC2D 模型[J]. 煤炭学报，2007，32(4): 358－362. XU Wen-jie, HU Rui-lin, WANG Yan-ping. PFC2D model for mesostructure of inhomogeneous geomaterial based on digital image processing[J]. Journal of China Coal Society, 2007, 32(4): 358－362. [10] 李凡. 岩土材料破损特性的颗粒流研究[J]. 土木工程 学报，2007，40(9)：78－81. LI Fan. A study on the breakage properties of geological materials using particle flow simulation[J]. China Civil Engineering Journal, 2007, 40(9): 78－81. [11] CHO N, MARTIN C D, SEGO D C. A clumped particle model for rock[J]. International Journal of Rock Mechanics and Mining Sciences, 2007, 44: 997－1010. [12] CHO N, MARTIN C D, SEGO D. Development of a shear zone in brittle rock subjected to direct shear[J]. International Journal of Rock Mechanics and Mining Sciences, 2008, (45): 1335－1346. [13] KULATILAKE P H S W, MALAMA B, WANG J. Physical and particle flow modeling of jointed rock block behavior under uniaxial loading[J]. International Journal of Rock Mechanics and Mining Sciences, 2001, (38): 641－657. [14] SCHOPFER M P J, CHILDS C, WALSH J J. Localisation of normal faults in multilayer sequences[J]. Journal of Structural Geology, 2006, (28): 816－833. [15] 王涛，盛谦，熊将. 基于颗粒流方法自然崩落法数值模 拟研究[J]. 岩石力学与工程学报，2007，26(增 2): 4202 －4207. WANG Tao, SHENG Qian, XIONG Jiang. Research on numerical simulation of natural caving method based on particle flow method[J]. Chinese Journal of Rock Mechanics and Engineering, 2007, 26(Supp.2): 4202－ 4207. [16] 杨仕教，曾晟，王和龙. 加载速率对石灰岩力学效应的 试验研究[J]. 岩土工程学报，27(7)：786－788. YANG Shi-jiao, ZENG Sheng, WANG He-long. Experimental analysis on mechanical effects of loading rates on limestone[J]. Chinese Journal of Geotechnical Engineering, 2005, 27(7): 786－788. [17] Itasca Consulting Group, Inc. PFC2D (partial flow code in 2 dimension) version 3.1[R]. Minneapolis, Minnesota: Itasca Consulting Group, Inc., 2001.