Bandis经验公式修正及其应用

第36卷第12期 北京科技大学学报 Vol.36 No.12 2014年12月 Journal of University of Science and Technology Beijing Dec.2014 Bandis经验公式修正及其应用 祁建东”，金爱兵，王贺”，高艳华”，于咏妍” 1)北京科技大学土木与环境工程学院，北京1000832)北京矿治研究总院，北京100070 ☒通信作者，E-mail:jinaibing@vip.sina.com 摘要针对三维离散元程序中刚度参数合理选取方法及其应用开展研究.通过理论分析，给出了结构面剪切刚度及法向刚 度的理论公式：结合室内剪切试验研究和数值分析，得到了符合工程实际的剪切刚度修正Bnds经验公式及法向刚度参数经 验公式，从而明确提出了一种结构面刚度参数的选取方法：通过实际工程边坡稳定性分析，探究了基于上述方法选取的参数 对实际工程模拟的适用性.研究表明，岩体结构面刚度参数是应力的函数，修正Bdis经验公式能够较为完善地表征剪切破 坏前结构面的剪切刚度变化规律：利用编制FSH程序使刚度参数随应力变化而改变的方法，使3DEC软件对工程岩体变形特 征的模拟更加符合实际，从而验证了该方法的合理性. 关键词岩体：刚度；经验公式：剪切试验；数值分析：边坡稳定性 分类号TU452 Correction and application of the Bandis empirical formula QI Jian-dong",JIN Ai-bing,WANG He2,GAO Yan-hua,YU Yong-yan 1)School of Civil and Environmental Engineering,University of Science and Technology Beijing,Beijing 100083,China 2)Beijing General Research Institute of Mining and Metallurgy,Beijing 100070,China Corresponding author,E-mail:jinaibing@vip.sina.com ABSTRACT A reasonable method was discussed for choosing and applying stiffness parameters in a 3D distinct element code.The formulas of pre-peak shear stiffness and normal stiffness of the fracture were derived through theoretical analysis.In combination with data from indoor shearing experiment,a corrected Bandis empirical formula of shear stiffness and a normal stiffness empirical formula were developed,which were proved to be logical by numerical analytical study.At last,a reality engineering of slope slide was discussed as a further evidence to describe the applicability of stiffness parameters chosen by the method mentioned above.It is obvious that the stiffness parameters of the fracture are a function of stress whether shear stiffness or normal stiffness.The corrected Bandis empirical formula can reflect the behavior of shear stiffness in the pre-peak range well.It is practicable to make the stiffness of the frac- ture changed with stress by developing a FISH program (in 3DEC),and therefore,the simulation can more in line with the practical engineering.In conclusion,the method mentioned above is proved to be reasonable KEY WORDS rock mass:stiffness:empirical formulas;shear test:numerical analysis:slope stability 岩体结构面为岩体内具有一定方向、延展较大键因素之一.Nassir等同通过理论分析给出了节理 且厚度较小的面状地质界面口.在高陡岩质边坡 岩体刚度参数的表征公式：肖卫国等0利用经典弹 中，岩体的变形破坏主要受岩体内结构面尤其是软塑性理论提出了充填节理岩体增量型应力与位移的 弱结构面的空间组合关系和结构面的力学性质控 本构模型：陈卫忠等固和王芝银等分别提出了基 制回.合理的工程岩体参数是工程稳定性评价的关 于力学实验的结构面刚度参数确定方法 收稿日期：2013-10-10 基金项目：中央直属高校基本科研业务费资助项目(FRF-SD-12O02A) DOI:10.13374/j.issn1001-053x.2014.12.002:http://journals.ustb.edu.cn

第 36 卷 第 12 期 2014 年 12 月 北京科技大学学报 Journal of University of Science and Technology Beijing Vol． 36 No． 12 Dec． 2014 Bandis 经验公式修正及其应用 祁建东1) ，金爱兵1) ，王 贺2) ，高艳华1) ，于咏妍1) 1) 北京科技大学土木与环境工程学院，北京 100083 2) 北京矿冶研究总院，北京 100070  通信作者，E-mail: jinaibing@ vip． sina． com 摘 要 针对三维离散元程序中刚度参数合理选取方法及其应用开展研究． 通过理论分析，给出了结构面剪切刚度及法向刚 度的理论公式; 结合室内剪切试验研究和数值分析，得到了符合工程实际的剪切刚度修正 Bandis 经验公式及法向刚度参数经 验公式，从而明确提出了一种结构面刚度参数的选取方法; 通过实际工程边坡稳定性分析，探究了基于上述方法选取的参数 对实际工程模拟的适用性． 研究表明，岩体结构面刚度参数是应力的函数，修正 Bandis 经验公式能够较为完善地表征剪切破 坏前结构面的剪切刚度变化规律; 利用编制 FISH 程序使刚度参数随应力变化而改变的方法，使 3DEC 软件对工程岩体变形特 征的模拟更加符合实际，从而验证了该方法的合理性． 关键词 岩体; 刚度; 经验公式; 剪切试验; 数值分析; 边坡稳定性 分类号 TU452 Correction and application of the Bandis empirical formula QI Jian-dong1) ，JIN Ai-bing1)  ，WANG He2) ，GAO Yan-hua1) ，YU Yong-yan1) 1) School of Civil and Environmental Engineering，University of Science and Technology Beijing，Beijing 100083，China 2) Beijing General Ｒesearch Institute of Mining and Metallurgy，Beijing 100070，China  Corresponding author，E-mail: jinaibing@ vip． sina． com ABSTＲACT A reasonable method was discussed for choosing and applying stiffness parameters in a 3D distinct element code． The formulas of pre-peak shear stiffness and normal stiffness of the fracture were derived through theoretical analysis． In combination with data from indoor shearing experiment，a corrected Bandis empirical formula of shear stiffness and a normal stiffness empirical formula were developed，which were proved to be logical by numerical analytical study． At last，a reality engineering of slope slide was discussed as a further evidence to describe the applicability of stiffness parameters chosen by the method mentioned above． It is obvious that the stiffness parameters of the fracture are a function of stress whether shear stiffness or normal stiffness． The corrected Bandis empirical formula can reflect the behavior of shear stiffness in the pre-peak range well． It is practicable to make the stiffness of the frac￾ture changed with stress by developing a FISH program ( in 3DEC) ，and therefore，the simulation can more in line with the practical engineering． In conclusion，the method mentioned above is proved to be reasonable． KEY WOＲDS rock mass; stiffness; empirical formulas; shear test; numerical analysis; slope stability 收稿日期: 2013--10--10 基金项目: 中央直属高校基本科研业务费资助项目( FＲF--SD--12--002A) DOI: 10． 13374 /j． issn1001--053x． 2014． 12． 002; http: / /journals． ustb． edu． cn 岩体结构面为岩体内具有一定方向、延展较大 且厚度较小的面状地质界面［1］． 在高陡岩质边坡 中，岩体的变形破坏主要受岩体内结构面尤其是软 弱结构面的空间组合关系和结构面的力学性质控 制［2］． 合理的工程岩体参数是工程稳定性评价的关 键因素之一． Nassir 等［3］通过理论分析给出了节理 岩体刚度参数的表征公式; 肖卫国等［4］利用经典弹 塑性理论提出了充填节理岩体增量型应力与位移的 本构模型; 陈卫忠等［5］和王芝银等［6］分别提出了基 于力学实验的结构面刚度参数确定方法．

·1576· 北京科技大学学报 第36卷 近年来，离散单元法作为一种显式求解的数值 式中，f(r/r,)为修正函数，A、B和C为常数.具体 方法被广泛应用于节理岩体的计算中).蒋坤和夏 修正过程后文将详细讨论，在此不赘述。 才初图利用离散元法基于不同的刚度参数模型对 1.2法向刚度 岩体边坡的稳定性进行了评价.3DEC程序可以用 实验研究表明，影响结构面法向刚度的因素主 于模拟静态或动态荷载下节理岩体的力学响应回. 要包括：结构面之间的实际初始接触面积、自相关距 巨能攀等0利用3DEC探讨了岩质边坡块体失稳 离及有效高度分量：结构面粗糙程度：充填物的厚 动态过程并进行了稳定性评价. 度、类型及力学性质等回 目前，针对数值模拟中刚度参数选取的研究并 根据文献2]所述，对于不包含充填物的结 不多见，在3DEC模拟分析时大多数学者将其简化 构面而言，法向应力σ。与法向变形△V符合经验 为常数，这种假设显然与工程实际有较大出入，模拟 公式 结果的准确性亦有待商榷.本文结合工程实例，通 lgon=p+q'△V (5) 过理论分析结合室内试验给出了结构面刚度参数的 式中，p和g为拟合常数 合理选取方法：通过数值实验，对Bandis剪切刚度 法向刚度K。反映结构面产生单位法向变形的 经验公式进行了修正，并编制FISH程序实现了边 法向应力梯度，即 坡稳定性分析中刚度参数随法向应力改变而改变， K.=aAV;-lge (6) 从而使数值分析更加符合工程实际，为类似工程数 0.3434 值分析提供参考 2结构面刚度参数试验研究 1结构面刚度参数理论分析 2.1结构面直接剪切试验 1.1剪切刚度 为探究结构面剪切刚度及法向刚度的变化规 对剪切破坏前的结构面剪切变形特性而言， 律，参考JTGE41一2005《公路工程岩石试验规 Kulhawy给出了其r-d,方程： 程》进行了结构面直接剪切试验 1 -=m· -+n. (1) 试验设备采用北京科技大学岩石力学实验室研 d 制的岩石弱面直剪仪.试验岩芯通过蜡封保持了天 式中，r为剪应力，d为剪切位移，m=1/K,K.为初 然结构面的性状特征（本次试验岩样结构面无充填 始剪切刚度，n=1/斤T为r-d双曲线水平渐近 物、面壁较为平整且起伏较小，因此可忽略表面性状 值. 不同对变形特性的影响)，经过切割加工、混凝土浇 剪切刚度K,即单位变形内的应力梯度为 筑及适当养护后，进行如下试验： k.()( (2) (1)将达到养护龄期的混凝土包裹试样放入直 剪仪，固定法向及剪切加载设备，安装位移计并将其 式中，K为剪切刚度系数，σ。为法向应力，n为剪切 调零校准： 刚度指数，R为破坏系数，R=T/r山 (2)首先施加法向载荷至0.2MPa,待位移读数 Bandis等通过实验研究引入函数r/r,修正 稳定后，读取法向位移.然后施加切向载荷，等间距 式(2)得到经验公式： 读取位移(d)一剪应力(r)组合数据： 2 K=K(o.）y1- ,T≤T。 (3) (3)当剪切位移变化速率显著增大且总剪切位 移达到试件直径的10%时，完成一次剪切试验； 式中，r。为抗剪强度峰值. (4)重复步骤(1)~(3)，仅改变法向载荷为 式(3)即为应用较为广泛的Bandis经验公式， 0.4、0.6和0.8MPa,完成一件试样剪切试验； 它能在一定程度上表征剪切刚度的实验变化规律. (5)重复步骤(1)~(4)，直至完成所有试验 但是，经过笔者反复地数值实验表明该公式对剪切 2.2结构面刚度参数拟合 刚度规律的反映并不完善，而需要重新修正，得到修 2.2.1剪切刚度 正Bandis经验公式： 将上述剪切试验数据进行拟合，得到剪切刚度 K=K(a）"1-fr/r,)·R]2, 参数如表1所示. f(r/,)=A(x/r,)2+B(r/r)+C, (4) 将拟合结果代入式(3)，得到剪切刚度K的 T≤Tp Bandis经验公式：

北 京 科 技 大 学 学 报 第 36 卷 近年来，离散单元法作为一种显式求解的数值 方法被广泛应用于节理岩体的计算中［7］． 蒋坤和夏 才初［8］利用离散元法基于不同的刚度参数模型对 岩体边坡的稳定性进行了评价． 3DEC 程序可以用 于模拟静态或动态荷载下节理岩体的力学响应［9］． 巨能攀等［10］利用 3DEC 探讨了岩质边坡块体失稳 动态过程并进行了稳定性评价． 目前，针对数值模拟中刚度参数选取的研究并 不多见，在 3DEC 模拟分析时大多数学者将其简化 为常数，这种假设显然与工程实际有较大出入，模拟 结果的准确性亦有待商榷． 本文结合工程实例，通 过理论分析结合室内试验给出了结构面刚度参数的 合理选取方法; 通过数值实验，对 Bandis 剪切刚度 经验公式进行了修正，并编制 FISH 程序实现了边 坡稳定性分析中刚度参数随法向应力改变而改变， 从而使数值分析更加符合工程实际，为类似工程数 值分析提供参考． 1 结构面刚度参数理论分析 1. 1 剪切刚度 对剪切破坏前的结构面剪切变形特性而言， Kulhawy［11］给出了其 τ － dh方程: 1 τ = m·1 dh + n． ( 1) 式中，τ 为剪应力，dh为剪切位移，m = 1 /Ksi，Ksi为初 始剪切刚度，n = 1 /τult，τult为 τ － dh双曲线水平渐近 值． 剪切刚度 Kst，即单位变形内的应力梯度为 Kst = τ dh = Kj ( σn ) nj ( 1 － Ｒf ) 2 ． ( 2) 式中，Kj为剪切刚度系数，σn为法向应力，nj为剪切 刚度指数，Ｒf为破坏系数，Ｒf = τ /τult ． Bandis 等［12］通过实验研究引入函数 τ /τp修正 式( 2) 得到经验公式: Kst = Kj ( σn ) n (j 1 － τ·Ｒf τ ) p 2 ，τ≤τp ． ( 3) 式中，τp为抗剪强度峰值． 式( 3) 即为应用较为广泛的 Bandis 经验公式， 它能在一定程度上表征剪切刚度的实验变化规律． 但是，经过笔者反复地数值实验表明该公式对剪切 刚度规律的反映并不完善，而需要重新修正，得到修 正 Bandis 经验公式: Kst = Kj ( σn ) nj ［1 － f( τ /τp )·Ｒf ］2 ， f( τ /τp ) = A( τ /τp ) 2 + B( τ /τp ) + C， τ≤τp { ． ( 4) 式中，f( τ /τp ) 为修正函数，A、B 和 C 为常数． 具体 修正过程后文将详细讨论，在此不赘述． 1. 2 法向刚度 实验研究表明，影响结构面法向刚度的因素主 要包括: 结构面之间的实际初始接触面积、自相关距 离及有效高度分量; 结构面粗糙程度; 充填物的厚 度、类型及力学性质等［12］． 根据文献［12］所述，对于不包含充填物的结 构面而言，法向应力 σn与法向变形 ΔVj符合经验 公式 lgσn = p + q·ΔVj ． ( 5) 式中，p 和 q 为拟合常数． 法向刚度 Kn反映结构面产生单位法向变形的 法向应力梯度，即 Kn = σn ΔVj = qσn lge = qσn 0. 3434． ( 6) 2 结构面刚度参数试验研究 2. 1 结构面直接剪切试验 为探究结构面剪切刚度及法向刚度的变化规 律，参 考 JTG E41—2005《公路工程岩石试验规 程》［13］进行了结构面直接剪切试验． 试验设备采用北京科技大学岩石力学实验室研 制的岩石弱面直剪仪． 试验岩芯通过蜡封保持了天 然结构面的性状特征( 本次试验岩样结构面无充填 物、面壁较为平整且起伏较小，因此可忽略表面性状 不同对变形特性的影响) ，经过切割加工、混凝土浇 筑及适当养护后，进行如下试验: ( 1) 将达到养护龄期的混凝土包裹试样放入直 剪仪，固定法向及剪切加载设备，安装位移计并将其 调零校准; ( 2) 首先施加法向载荷至 0. 2 MPa，待位移读数 稳定后，读取法向位移． 然后施加切向载荷，等间距 读取位移( dh ) --剪应力( τ) 组合数据; ( 3) 当剪切位移变化速率显著增大且总剪切位 移达到试件直径的 10% 时，完成一次剪切试验; ( 4) 重复步骤( 1) ～ ( 3) ，仅改变法向载荷为 0. 4、0. 6 和 0. 8 MPa，完成一件试样剪切试验; ( 5) 重复步骤( 1) ～ ( 4) ，直至完成所有试验． 2. 2 结构面刚度参数拟合 2. 2. 1 剪切刚度 将上述剪切试验数据进行拟合，得到剪切刚度 参数如表 1 所示． 将拟合结果代入式( 3) ，得到剪切刚度 Kst 的 Bandis 经验公式: · 6751 ·

第12期 祁建东等：Bandis经验公式修正及其应用 ·1577· K=1.45a88. 表2岩石物理力学参数 Table 2 Physico-mechanical characteristics of rock [1-0.3050.+0.029j(0.360.+0.03》 岩石密度， 体积模量， 剪切模量， 摩擦角， 黏聚力， (7 p/(kg'm-3) K/GPa G/GPa /() c/MPa 表1剪切刚度参数拟合 2220 8.45 7.71 10.19 Table 1 Shear stiffness parameters fitting 利用3DEC内置函数GEN划分网格，共生成 Go/ Tp/ K./ Tulk/ 40930个有限差分域. MPa MPa (GPa-m-1)MPa 边界条件：在z=60及x=-30处施加不同应 0.2 0.10 0.55 0.104 力边界，其余边界面在各自法向量方向施加速度边 0.4 0.14 0.73 0.169 1.45 0.63 界v=0. 0.60.21 1.08 0.215 3.1基于Bandis经验公式数值实验 0.8 0.28 1.30 0.315 合理选取结构面刚度是数值分析结构面变形特 2.2.2法向刚度 性的保障.首先采用Bandis经验公式选取剪切刚度 利用实测数据通过共轭梯度法拟合得到lgσ。- 及法向刚度参数如表3所示，并据此进行了三组不 △V关系式： 同应力组合的数值实验.分析中利用History函数 lg0m=7.34△V-1.77. (8) 监测记录结构面中心点处的剪切位移. 结合式(6)及式(8)得到结构面法向刚度K。: 表3基于Bandis经验公式法选取的刚度参数 7.340m Table 3 Stiffness parameters for numerical experiments according to the K.=0.3434 (9) Bandis formula 3 Bandis经验公式修正 K/ K/ o。/MPa T/MPa (MPa“mm-l) (MPa'mm-) 通过基于室内试验结果的数值实验，利用反演 0.05 0.317 4.275 法对Bandis经验公式进行修正.本文中数值实验采 0.20 0.10 0.005 4.275 用三维离散元程序3DEC.该程序是基于粗糙结构 0.07 0.489 8.336 面的假设而开发的，虽然不能探究结构面不同性状 0.40 0.14 0.008 8.336 特征对其变形特性的不同影响，但对结构面的变形 0.11 0.595 12.611 特性进行定性半定量的分析是可行的. 0.60 0.21 0.009 12.611 综合考量边界效应影响及简化计算后建立模型 0.07 1.123 16.886 为60mm×60mm×120mm立方体，如图1所示. 0.14 0.755 16.886 0.80 0.21 0.304 16.886 0.28 0.011 16.886 当法向应力σ。=0.8MPa保持不变时，等差取 结构面 剪应力r直至r=T。=0.28MPa共4种组合，数值 实验结果如图2所示. 当法向应力σ。依次取为0.2、0.4、0.6和 0.8MPa时，选取各自对应的剪应力r=0.5r,共四种 组合，数值实验结果如图3所示 当法向应力o。依次取为0.2、0.4、0.6和 0.8MPa时，选取各自对应的剪应力T=T,共四种组 图1数值实验模型 Fig.1 Numerical experiment model 合，数值实验结果如图4所示. 从图2~图4可以看出：基于Bandis经验公式 选用刚性块体模型模拟砂岩试件，结构面选用 的刚度参数在剪应力相对较小（即，≤0.5r,)时，能 Coulomb滑移模型.砂岩物理力学参数见表2. 够得到与室内试验结果相近的剪切变形数值模拟

第 12 期 祁建东等: Bandis 经验公式修正及其应用 Kst = 1. 45σ0. 63 n [ · 1 － τ 2 ( 0. 305σn + 0. 029) ( 0. 336σn + 0. 033 ] ) 2 ． ( 7) 表 1 剪切刚度参数拟合 Table 1 Shear stiffness parameters fitting σn / MPa τp / MPa Ksi / ( GPa·m － 1 ) τult / MPa Kj nj 0. 2 0. 10 0. 55 0. 104 0. 4 0. 14 0. 73 0. 169 1. 45 0. 63 0. 6 0. 21 1. 08 0. 215 0. 8 0. 28 1. 30 0. 315 2. 2. 2 法向刚度 利用实测数据通过共轭梯度法拟合得到 lgσn － ΔVj关系式: lgσn = 7. 34ΔVj － 1. 77． ( 8) 结合式( 6) 及式( 8) 得到结构面法向刚度 Kn : Kn = 7. 34σn 0. 3434． ( 9) 3 Bandis 经验公式修正 通过基于室内试验结果的数值实验，利用反演 法对 Bandis 经验公式进行修正． 本文中数值实验采 用三维离散元程序 3DEC． 该程序是基于粗糙结构 面的假设而开发的，虽然不能探究结构面不同性状 特征对其变形特性的不同影响，但对结构面的变形 特性进行定性半定量的分析是可行的． 综合考量边界效应影响及简化计算后建立模型 为 60 mm × 60 mm × 120 mm 立方体，如图 1 所示． 图 1 数值实验模型 Fig． 1 Numerical experiment model 选用刚性块体模型模拟砂岩试件，结构面选用 Coulomb 滑移模型． 砂岩物理力学参数见表 2． 表 2 岩石物理力学参数 Table 2 Physico-mechanical characteristics of rock 岩石密度， ρ /( kg·m － 3 ) 体积模量， K /GPa 剪切模量， G /GPa 摩擦角， φ/( °) 黏聚力， c/MPa 2220 8. 45 7. 71 43 10. 19 利用 3DEC 内置函数 GEN 划分网格，共生成 40930 个有限差分域． 边界条件: 在 z = 60 及 x = － 30 处施加不同应 力边界，其余边界面在各自法向量方向施加速度边 界 v = 0． 3. 1 基于 Bandis 经验公式数值实验 合理选取结构面刚度是数值分析结构面变形特 性的保障． 首先采用 Bandis 经验公式选取剪切刚度 及法向刚度参数如表 3 所示，并据此进行了三组不 同应力组合的数值实验． 分析中利用 History 函数 监测记录结构面中心点处的剪切位移． 表 3 基于 Bandis 经验公式法选取的刚度参数 Table 3 Stiffness parameters for numerical experiments according to the Bandis formula σn /MPa τ /MPa Kst / ( MPa·mm － 1 ) Kn / ( MPa·mm － 1 ) 0. 20 0. 05 0. 317 4. 275 0. 10 0. 005 4. 275 0. 40 0. 07 0. 489 8. 336 0. 14 0. 008 8. 336 0. 60 0. 11 0. 595 12. 611 0. 21 0. 009 12. 611 0. 07 1. 123 16. 886 0. 80 0. 14 0. 755 16. 886 0. 21 0. 304 16. 886 0. 28 0. 011 16. 886 当法向应力 σn = 0. 8 MPa 保持不变时，等差取 剪应力 τ 直至 τ = τp = 0. 28 MPa 共 4 种组合，数值 实验结果如图 2 所示． 当法 向 应 力 σn 依 次 取 为 0. 2、0. 4、0. 6 和 0. 8 MPa时，选取各自对应的剪应力 τ = 0. 5τp共四种 组合，数值实验结果如图 3 所示． 当法 向 应 力 σn 依 次 取 为 0. 2、0. 4、0. 6 和 0. 8 MPa时，选取各自对应的剪应力 τ = τp共四种组 合，数值实验结果如图 4 所示． 从图 2 ～ 图 4 可以看出: 基于 Bandis 经验公式 的刚度参数在剪应力相对较小( 即 τ≤0. 5τp ) 时，能 够得到与室内试验结果相近的剪切变形数值模拟 · 7751 ·

·1578 北京科技大学学报 第36卷 18r 值；当剪应力？增大并趋近于T,时，数值模拟则会 出现较大偏差.这说明Bandis经验公式得出的刚度 15 参数仅在一定程度上反映了结构面变形特性，并不 12 完善，欲将其应用于3DEC刚度参数选取中，须对其 数值实验结果 进行合理修正 3.2基于数值实验结果的Bandis经验公式修正 6 经过大量数值实验分析，就3DEC数值分析时 试验实测结果 的剪切刚度参数选取而言，式(2)应重新引入修正 3 函数f(x/:,)才能更好地符合室内试验得出的变化 - 规律. 0.05 0.10 0.15 0.200.250.30 经过反复试算及参数反演，给出了当法向应力 /MPa c.=0.8MPa恒定时应力组合对应的K值，详见 图2 。=0.8MPa、不同，下试验及模拟结果对比 表4 Fig.2 Comparison of tests and numerical results for different r with =0.8 MPa 表4基于Bandis经验公式法选取的刚度参数(an=0.8MPa) Table 4 Stiffness parameters for numerical experiments according to the Bandis formula =0.8 MPa) 0.21 剪应力， 剪切刚度， 法向刚度， 0.20 试验实测结果 /MPa K/(MPa'mm-1) K./(MPa'mm-1) 0.19 0.07 1.030 16.886 0.18 0.14 0.730 16.886 0.21 0.440 16.886 0.16 0.28 0.150 16.886 0.15 数值实验结果 0.14< >0.14 通过比较分析，修正函数f(r/r。)应为二次多 0.2 0.12 0.4 0.10 项式形式，即式(4)中所示. 0.6 0.08 /MPa 0.8 .06 本例通过最小二乘法拟合上述数值结果得到： 0.04 f(r/:,）=0.18(r/r,)2+0.09(r/r,)+0.45. 图3不同o.及r=0.5rp组合下试验及模拟结果对比 (10) Fig.3 Comparison of test and numerical results under different 3.3基于修正Bandis经验公式数值实验 groups of o and T=0.5Tp 针对该工程利用修正Bandis经验公式(4)求得 剪切刚度参数见表5.取法向应力σ.依次取为0.2、 18 0.4、0.6和0.8MPa时，选取各自对应的剪应力T= 0.5r,共四种组合：取法向应力σ.依次取为0.2、 12 0.4、0.6和0.8MPa时，选取各自对应的剪应力T= 数值实验结果 试验实测结果 T,共四种组合.通过数值实验，得到模拟结果如图5 所示 从图5易知，基于修正Banids经验公式进行的 数值实验与室内试验结果一致，从而说明经过修正 0.4 的经验公式是合理的 02 03 0.4 0.2 3.4修正的Bandis经验公式应用讨论 0.6 G/MPa 0.8 70.1 T/MPa 基于3DEC的结构面变形特性分析中，难以实 0 现结构面剪切刚度取值随切应力？的变化而改变. 图4不同σ。及r=T组合下实验及模拟结果对比 选取r=T时式(4)计算所得K作为该应力条件下 Fig.4 Comparison of test and numerical results under different group 的结构面剪切刚度参考值，此时剪切刚度值可近似 of oa and T=Tp 简化为法向应力的函数：

北 京 科 技 大 学 学 报 第 36 卷 图 2 σn = 0. 8 MPa、不同 τ 下试验及模拟结果对比 Fig． 2 Comparison of tests and numerical results for different τ with σn = 0. 8 MPa 图 3 不同 σn及 τ = 0. 5τp组合下试验及模拟结果对比 Fig． 3 Comparison of test and numerical results under different groups of σn and τ = 0. 5τp 图 4 不同 σn及 τ = τp组合下实验及模拟结果对比 Fig． 4 Comparison of test and numerical results under different group of σn and τ = τp 值; 当剪应力 τ 增大并趋近于 τp时，数值模拟则会 出现较大偏差． 这说明 Bandis 经验公式得出的刚度 参数仅在一定程度上反映了结构面变形特性，并不 完善，欲将其应用于 3DEC 刚度参数选取中，须对其 进行合理修正． 3. 2 基于数值实验结果的 Bandis 经验公式修正 经过大量数值实验分析，就 3DEC 数值分析时 的剪切刚度参数选取而言，式( 2) 应重新引入修正 函数 f( τ /τp ) 才能更好地符合室内试验得出的变化 规律． 经过反复试算及参数反演，给出了当法向应力 σn = 0. 8 MPa 恒定时应力组合对应的 Kst 值，详见 表 4． 表 4 基于 Bandis 经验公式法选取的刚度参数( σn = 0. 8 MPa) Table 4 Stiffness parameters for numerical experiments according to the Bandis formula ( σn = 0. 8 MPa) 剪应力， τ /MPa 剪切刚度， Kst /( MPa·mm － 1 ) 法向刚度， Kn /( MPa·mm － 1 ) 0. 07 1. 030 16. 886 0. 14 0. 730 16. 886 0. 21 0. 440 16. 886 0. 28 0. 150 16. 886 通过比较分析，修正函数 f( τ /τp ) 应为二次多 项式形式，即式( 4) 中所示． 本例通过最小二乘法拟合上述数值结果得到: f( τ /τp ) = 0. 18( τ /τp ) 2 + 0. 09( τ /τp ) + 0. 45． ( 10) 3. 3 基于修正 Bandis 经验公式数值实验 针对该工程利用修正 Bandis 经验公式( 4) 求得 剪切刚度参数见表 5． 取法向应力 σn依次取为 0. 2、 0. 4、0. 6 和 0. 8 MPa 时，选取各自对应的剪应力 τ = 0. 5τp共四种组合; 取法向应力 σn 依次取为 0. 2、 0. 4、0. 6 和 0. 8 MPa 时，选取各自对应的剪应力 τ = τp共四种组合． 通过数值实验，得到模拟结果如图 5 所示． 从图 5 易知，基于修正 Banids 经验公式进行的 数值实验与室内试验结果一致，从而说明经过修正 的经验公式是合理的． 3. 4 修正的 Bandis 经验公式应用讨论 基于 3DEC 的结构面变形特性分析中，难以实 现结构面剪切刚度取值随切应力 τ 的变化而改变． 选取 τ = τp时式( 4) 计算所得 Kst作为该应力条件下 的结构面剪切刚度参考值，此时剪切刚度值可近似 简化为法向应力的函数: · 8751 ·

第12期 祁建东等：Bandis经验公式修正及其应用 ·1579· 表5基于修正的Bandis经验公式法选取的刚度参数 于边帮开挖线，该范围内滑坡体平面形似扇形， Table 5 Stiffness parameters for numerical experiments according to the 横宽55m,纵长40m,总滑坡体体积约16500m3, modified Bandis formula 详见图6. O./MPa /MPa K./(MPa*mm-1) K.(MPa.mm-1) 0.05 0.268 4.275 0.20 0.10 0.061 4.275 0.07 0.432 8.336 0.40 0.14 0.089 8.336 0.11 0.552 12.611 0.60 0.21 0.108 12.611 0.14 0.738 16.886 图6滑坡现场资料图 0.80 0.28 0.152 16.886 Fig.6 Picture of the landslide site 4.23DEC模型及参数选取 3.0 根据工程实际建立边坡分析模型为200m× 2.4 100m×100m,本构模型采用刚性块体模型，结构面 试验实测结果 采用Coulomb滑移模型.初始应力场取为重力场. 1.8 岩体力学参数详见表6. 12 表6岩体物理力学参数 数值实验结果 Table 6 Physico-mechanical characteristics of rock mass 0.6 岩石密度， 体积模量，剪切模量，摩擦角，黏聚力， >0.30 岩性 K/GPa G/GPa p/(o）c/MPa 0.2 0.24 p/(kg'm-3) 0.4 .1 16 0.12 t/MPa 砂岩 2220 7.21 3.04 33 6.00 /MPa 0.8 0.06 1.00 泥岩 2520 5.88 2.26 25 7.16 图5不同σ。及：组合下实验及模拟结果对比 Fig.5 Comparison of tests and numerical results under different 根据上文讨论令？=T。,结合式(10)和式(11) groups of o and 并做适当折减7一W后给出针对该工程的节理剪切 刚度经验公式： Ka=D.Ki()"). (11) K=0.0826. (12) 式中，D为常数 通过3DEC内置的FISH语言编制程序，自定义 式(11)能够模拟剪切破坏前结构面的剪切变 结构面变形本构模型.利用自定义的函数自动修改 形特征，且不产生不可接受的误差，在数值分析中也 结构面的剪切及法向刚度参数，具体步骤如下： 易于实现。对于变形误差要求不甚严格的工程模拟 (1)建模过程中利用D对结构面进行标记，记 计算中，可以实现变形特征的半定量分析. 为ID=101,ID=102,D=103: 4工程实例 (2)利用FISH定义函数J_St,该函数根据结构 面的法向应力结合式(9)和式(12)计算出剪切及法 在矿产资源开发中，露天开采作为主要的生产 向刚度值并分别赋予各结构面，其流程见图7. 模式被广泛应用.随着开采强度的不断加大，高 (3)编制FISH函数C_S,利用该函数实现模拟 陡边坡日益增多，随之而产生的边坡安全问题日益 过程中每循环200步自动调用J_St函数，从而实现 突出s.岩体中的结构面无疑是边坡稳定的主 结构面剪切及法向刚度参数的自动调整. 要不利因素之一.本节利用编制FISH函数实现法 4.33DEC数值计算及结果分析 向及剪切刚度经验公式(9)和(11)的工程应用，以 4.3.1模拟工况 期为同类型的边坡稳定性分析提供借鉴 为分析产生滑坡的主要影响因素，本次模拟分 4.1工程概况 两部分： 某露天矿首采区东北帮砂岩地层范围内的 (1)模型包含两个近似垂直边坡面结构面D= 次典型结构面控制型滑坡.滑动方向基本垂直 101(倾向215°，倾角75)及ID=102(倾向325°，倾

第 12 期 祁建东等: Bandis 经验公式修正及其应用 表 5 基于修正的 Bandis 经验公式法选取的刚度参数 Table 5 Stiffness parameters for numerical experiments according to the modified Bandis formula σn /MPa τ /MPa Kst /( MPa·mm － 1 ) Kn /( MPa·mm － 1 ) 0. 20 0. 05 0. 268 4. 275 0. 10 0. 061 4. 275 0. 40 0. 07 0. 432 8. 336 0. 14 0. 089 8. 336 0. 60 0. 11 0. 552 12. 611 0. 21 0. 108 12. 611 0. 80 0. 14 0. 738 16. 886 0. 28 0. 152 16. 886 图 5 不同 σn及 τ 组合下实验及模拟结果对比 Fig． 5 Comparison of tests and numerical results under different groups of σn and τ Kst = D·Kj ( σn ) nj ． ( 11) 式中，D 为常数． 式( 11) 能够模拟剪切破坏前结构面的剪切变 形特征，且不产生不可接受的误差，在数值分析中也 易于实现． 对于变形误差要求不甚严格的工程模拟 计算中，可以实现变形特征的半定量分析． 4 工程实例 在矿产资源开发中，露天开采作为主要的生产 模式被广泛应用［14］． 随着开采强度的不断加大，高 陡边坡日益增多，随之而产生的边坡安全问题日益 突出［15--16］． 岩体中的结构面无疑是边坡稳定的主 要不利因素之一． 本节利用编制 FISH 函数实现法 向及剪切刚度经验公式( 9) 和( 11) 的工程应用，以 期为同类型的边坡稳定性分析提供借鉴． 4. 1 工程概况 某露天矿首采区东北帮砂岩地层范围内的 一次典型结构面控制型滑坡． 滑动方向基本垂直 于边帮开挖线，该范围内滑坡体平面形似扇形， 横宽 55 m，纵长 40 m，总滑坡体体积约 16500 m3 ， 详见图 6． 图 6 滑坡现场资料图 Fig． 6 Picture of the landslide site 4. 2 3DEC 模型及参数选取 根据工程实际建立边坡分析模型为 200 m × 100 m × 100 m，本构模型采用刚性块体模型，结构面 采用 Coulomb 滑移模型． 初始应力场取为重力场． 岩体力学参数详见表 6． 表 6 岩体物理力学参数 Table 6 Physico-mechanical characteristics of rock mass 岩性 岩石密度， ρ /( kg·m － 3 ) 体积模量， K /GPa 剪切模量， G /GPa 摩擦角， φ/( °) 黏聚力， c/MPa 砂岩 2220 7. 21 3. 04 33 6. 00 泥岩 2520 5. 88 2. 26 25 7. 16 根据上文讨论令 τ = τp，结合式( 10) 和式( 11) 并做适当折减［17--18］后给出针对该工程的节理剪切 刚度经验公式: Kst = 0. 082σ0. 63 n ． ( 12) 通过 3DEC 内置的 FISH 语言编制程序，自定义 结构面变形本构模型． 利用自定义的函数自动修改 结构面的剪切及法向刚度参数，具体步骤如下: ( 1) 建模过程中利用 ID 对结构面进行标记，记 为 ID = 101，ID = 102，ID = 103; ( 2) 利用 FISH 定义函数 J_St，该函数根据结构 面的法向应力结合式( 9) 和式( 12) 计算出剪切及法 向刚度值并分别赋予各结构面，其流程见图 7． ( 3) 编制 FISH 函数 C_S，利用该函数实现模拟 过程中每循环 200 步自动调用 J_St 函数，从而实现 结构面剪切及法向刚度参数的自动调整． 4. 3 3DEC 数值计算及结果分析 4. 3. 1 模拟工况 为分析产生滑坡的主要影响因素，本次模拟分 两部分: ( 1) 模型包含两个近似垂直边坡面结构面ID = 101( 倾向 215°，倾角 75°) 及 ID = 102( 倾向 325°，倾 · 9751 ·