李坤蒙,李元辉,徐帅,安龙
东北大学 深部金属矿山安全开采教育部重点实验室, 辽宁 沈阳 110819
收稿日期: 2015-02-10
基金项目: 国家自然科学基金资助项目(51204031, 51274055); 国家科技支撑计划项目(2013BAB02B03); 中央高校基本科研业务费专项资金资助项目(N130401007).
作者简介: 李坤蒙(1990-), 男, 陕西西安人, 东北大学博士研究生; 李元辉(1968-), 男, 辽宁大石桥人, 东北大学教授, 博士生导师.

摘要: 针对目前应用“试错法”来确定PFC^2D数值计算模型微观参数所存在的盲目性问题，提出一种基于物理试验响应的参数确定方法.通过模拟具体对象的基本特性计算模型物理几何参数，借助物理试验响应反演计算颗粒间的接触本构参数.以焦家金矿-450中段下盘糜棱岩为例，在单轴抗压物理试验的基础上，确定相应的数值计算模型微观参数并进行单轴、双轴模拟试验，模拟结果与物理试验具有良好的一致性，验证了微观参数确定方法的可靠性和适用性.PFC^2D微观参数确定方法能够在保证精确模拟的基础上，大大降低前期建模工作量，为类似研究提供借鉴.
关键词: PFC^2D数值模拟微观参数确定试错法反演计算
Method to Determine Microscopic Parameters of PFC^2D Numerical Model
LI Kun-meng,LI Yun-hui,XU Shuai,AN Long
Key Laboratory of Ministry of Education on Safe Mining of Deep Metal Mines, Northeastern University, Shenyang 110819, China.
Corresponding author: LI Kun-meng, professor, E-mail: 13840061753@163.com
Abstract: Due to the blindness of the current “trial and error method” used in determining the microscopic parameters in the PFC^2D, a microscopic parameters determination method was put forward based on physical experiments. The physical and geometrical parameters of the simulation model were obtained by simulating the basic characteristics of real samples, and the contact constitutive parameters between particles were back calculated base on the physical test responses. Taking the mylonite of -450 level footwall in Jiaojia gold mine as an example, the micro parameters of the corresponding numerical model were determined, and uniaxial and biaxial simulation were carried out based on the uniaxial compressive physical test. The simulation results and physical responses have good consistency, which verifies the reliability and applicability of the proposed method. The PFC^2D microcosmic parameters determination method can ensure simulation accuracy, reduce the early modeling workload and offer references for similar researches.
Key words: PFC^2Dnumerical simulationmicroscopic parameter determinationtrial and error methodback calculation 岩石在加载过程中的力学行为和破裂时空演化规律研究是地下岩体稳定性分级、采矿方法选择及采场结构参数设计的基础.目前用来模拟岩石损伤与变形行为的数值模拟方法主要包括连续介质模型和离散单元法(DEM).相对于连续介质模型而言，一方面，DEM借助模拟微观单元体间的力学行为来叠加生成宏观响应，不需要对模型预设复杂的数学本构关系；另一方面，DEM由于不涉及网格的划分，不必考虑划分尺寸的敏感性.基于离散单元法的颗粒流程序PFC^2D[1]能够较好地模拟岩石在加载过程中的局部各向异性行为，重现微裂隙时空演化过程和最终断裂区，已被广泛应用于研究岩石的各项力学行为^{[2, 3, 4, 5]}.
在应用PFC^2D解决特定问题之前，需要确定颗粒集合微观参数，包括颗粒尺寸分布、颗粒间接触刚度和强度等，计算模型微观参数的正确取值是获取精确模拟结果的关键.然而，岩石的力学行为通常借助如弹性变形常数、抗压强度和残余强度等宏观参数来描述.由于数值模拟微观参数和物理试验宏观响应的关系是非线性的，因而无法借助实验室所得宏观响应直接计算模拟所需的微观参数值.
目前国内外学者普遍应用“试错法”^[3]来确定PFC^2D数值计算模型微观参数，存在费时费力及盲目性大等问题，而且多参数确定过程中，纯粹的试错极难同时获得精确的微观参数值^{[3, 4, 5, 6]}.在应用过程中，文献仅提及“试错法”，在参数确定环节直接给出具体值，未对具体确定方法和判别准则作详细论述^{[7, 8, 9]}.
针对“试错法”存在的盲目性大和精确度低等问题，本文提出基于物理试验响应的PFC^2D微观参数确定方法，以焦家金矿-450中段下盘糜棱岩为例，验证具体方法的可靠性和适用性.
1 方法详述颗粒流程序PFC^2D数值模型微观参数包括物理几何参数和接触本构参数：物理几何参数通过具体模拟对象的基本特性确定，接触本构参数需要借助物理试验响应进行反演计算.
1.1 物理几何参数物理几何参数包括模拟试样高度H、宽度D、颗粒单元体最大半径R_max、最小半径R_min、粒径分布规律、颗粒集合孔隙率n、颗粒密度ρ，确定流程如图 1所示，具体过程如下.
图 1(Figure 1)

图 1 物理几何参数确定流程 Fig. 1 Determination flow of physical and geometrical parameters

1) H和D：按照模拟试样与物理试验岩样尺寸一致原则，通过测量方式确定H和D.
2) R_max，R_min、粒径分布规律和n：借助电镜扫描仪或显微仪器来观察试验岩样，通过绘制手段，获取统计层面每个晶粒的周长、面积、最大半径、最小半径、筛网尺寸、实际尺寸和误差等信息，绘制晶粒粒级分布曲线，确定R_max，R_min、粒径分布规律和n.
3) ρ：称取试验岩样质量、体积，去除孔隙所占比例，计算参数ρ.
1.2 接触本构参数PFC^2D数值计算模型接触本构参数主要包括颗粒间接触模量E_c、正向切向刚度比K_n/K_s、正向强度T_n、切向强度T_s、正向强度偏差T_dn、切向强度偏差T_ds、颗粒单元体之间的摩擦系数u.颗粒间的接触本构模型如图 2所示，A，B代表颗粒单元体，可以看出，颗粒集合刚度和强度分别由K_n，K_s和T_n，T_s，u决定^[1].
图 2(Figure 2)

图 2 PFC2D微观接触本构模型 Fig. 2 PFC2D micro contact constitutive model

借助单轴或常规三轴物理试验，获取岩石弹性模量E，泊松比μ，峰值强度б_c，裂隙初始强度б_ci及峰后残余强度б_ca.在保证边界、加载及截止条件等与物理试验一致的基础上，开展数值模拟试验.通过按顺序不断更改各接触本构参数的方式，建立数值模拟微观参数与宏观响应之间的对应关系，按照数值模拟与物理试验宏响应相等原则，计算各类接触本构参数具体值.确定流程如图 3所示，具体过程可如下.
图 3(Figure 3)

图 3 接触本构参数确定流程 Fig. 3 Determination flow of contact constitutive parameters

1) E_c和K_n/K_s：将模拟材料强度参数T_n，T_s设置为较大值(大约10倍的单轴抗压强度)，改变E_c，建立E_c和E之间的关系，借助已知的E计算E_c.同理，改变K_n/K_s的取值，确定μ与K_n/K_s的关系，而K_n/K_s的改变会反过来影响E值，因此，需要借助迭代手段来确定微观参数E_c和K_n/K_s，其中迭代方法选用二分法.
2) T_n，T_s，T_dn和T_ds：设置T_dn，T_ds为0，通过改变T_n，T_s来建立б_c与T_n，T_s的关系.值得注意，T_n/T_s的取值会影响矫正结果，必须保证T_n/T_s不变.数值模拟过程中，T_n/T_s值的大小控制了剪切和张破坏微裂隙生成数量，进而影响最终的破裂模式.一般情况下，小比率值可体现如钢铁等延性材料的力学行为，大比率值与岩石等脆性材料性能吻合^[2].Fakhimi等^[9]选择T_n/T_s=1:5来模拟孔洞岩石的力学行为，模拟结果与实际岩石物理试验一致，本文确定T_n/T_s=1:5.保证T_ds/T_s=T_dn/T_n，同时增加T_dn，T_ds，确定б_ci与T_dn，T_ds的关系，裂隙初始生成强度可借助如声发射等其他手段获得.而此过程会影响б_c，同样应用二分法迭代2)、以获参数T_s，T_n，T_dn和T_ds.
3) u：离散元软件PFC^2D在数值模拟过程中，颗粒间的破裂服从库伦准则，黏结破裂之前，材料强度由T_n，T_s决定，黏结破裂之后，材料强度主要由u确定.通过改变u的取值，建立岩石峰后行为与u的关系，通过数值模拟与物理试验的峰后残余强度值б_ca相等原则确定u.
2 实例分析基于焦家金矿-450中段下盘糜棱岩单轴抗压物理试验，确定对应的PFC^2D数值计算模型微观参数，通过对比单轴、双轴模拟试验结果和物理试验响应，分别验证确定方法的可靠性和适用性.
2.1 物理试验简述将现场钻取的岩芯按照工程岩体试验标准进行单轴抗压试验，在加载过程中，利用声发射技术实时监测破裂事件^{[10, 11]}，加载系统如图 4所示，试样尺寸及物理参数、加载条件和试验结果见表 1.
图 4(Figure 4)

图 4 物理试验加载系统 Fig. 4 Physical test loading system

表 1(Table 1)

表 1 单轴抗压物理试验条件及结果 Table 1 Conditions and results of uniaxial compressive physical test 试样尺寸及物理参数 加载条件 结果 (H×D)/(mm×mm) ρ/(kg·m-3) V/(mm·s-1) бc/MPa E/GPa μ бci/MPa бca/MPa 100×50 2 783 0.001 74.5 34.8 0.26 52.1 16.1

表 1 单轴抗压物理试验条件及结果 Table 1 Conditions and results of uniaxial compressive physical test

2.2 数值模型参数按照PFC^2D微观参数确定方法，获取数值计算模型微观参数，最终的物理几何和接触本构参数见表 2.糜棱岩的电镜扫描结果见图 5.可以看出，岩样晶粒辨析度较差，仅能从局部区域观察晶粒和孔隙大小及分布.为了尽量精确确定粒径分布规律，对同批试样的上下表面均进行扫描，将统计的整体结果作为R_max，R_min、粒径分布规律和n的最终取值.
表 2(Table 2)

表 2 数值计算模型的微观参数 Table 2 Microscopic parameters of numerical model Rmax/Rmin Rmin/mm 粒径分布 ρ/(kg·m-3) n/% Ec/GPa Kn/Ks Tn/MPa Tn/Ts Tdn/Tn(Tds/Ts) u 1.66 0.02 均匀 2 630 9 60 4.5 45 1/5 1:10 0.5

表 2 数值计算模型的微观参数 Table 2 Microscopic parameters of numerical model

图 5(Figure 5)

图 5 岩样细观电镜扫描 Fig. 5 Meso SEM of rock sample

2.3 可靠性验证借助确定的微观参数建立数值模型并进行单轴抗压模拟试验，将数值试验的应力应变曲线和最终破裂模式与物理试验响应进行对比，结果如图 6所示.图 6a可知，两者应力应变曲线在峰值强度之前基本重合，而峰值强度之后的差异较大，原因在于PFC^2D是借助模拟微观单元体之间的力学行为来叠加生成宏观响应，不存在轴向应力的陡然下降.同时，由数值模拟裂隙数演化曲线可知，模拟试验的б_ci为51.4 MPa，与实验室声发射首个信号接收时的平均轴向应力值52.1 MPa几乎相等.图 6b~图 6c显示了岩石数值模拟和物理试验最终加载结束后的失稳破裂模式，图 6d展示了岩石在破裂过程中的微裂隙时空演化过程，圆圈代表微裂隙，白色至黑色代表微裂隙生成的时间前后顺序.可以看出，糜棱岩在单轴受载条件下，剪切断裂带主要位于岩石中央，方向大约与轴向应力成20°角.
图 6(Figure 6)

图 6 单轴抗压物理试验与数值模拟结果对比 Fig. 6 Comparison between results for uniaxial compressive physical test and numerical test (a)—轴向应力应变曲线； (b)—物理试验; (c)—数值模拟； (d)—数值模拟微裂隙演化.

2.4 适用性评价借助矫正参数建立数值模型，实施不同围压下的双轴模拟试验，将模拟岩石不同围压下的强度值和最终破裂模式与常规三轴抗压物理试验响应进行对比.由强度曲线对比结果可知，实验室不同围压下的岩石强度值虽然存在一定的离散性，但其拟合曲线基本与数值模拟曲线重合，见图 7a.同时，图 7b展示了10 MPa围压条件下，岩石物理试验与数值模拟试验的最终破裂模式(其他围压下的结果一致)，结果吻合性较好.双轴模拟试验与物理试验对比结果表明，基于物理试验的PFC^2D微观参数确定方法具有较强的适用性，能够用来计算其他边界条件下的岩石力学行为.
图 7(Figure 7)

图 7 不同围压下物理试验与数值模拟结果对比 Fig. 7 Comparison of results between physical experiment and numerical simulation under different confined pressure (a)—双轴试验强度曲线； (b)—10 MPa围压破裂模式.

3 结论1) 以物理试验模型与数值计算模型的基本特性一致为原则，针对具体研究对象，借助电镜扫描、称重和测量等手段，获取数值计算模型物理几何参数.
2) 通过有顺序地不断更改各接触本构参数的方式，建立数值模拟微观参数与宏观响应之间的线性关系，按照数值模拟与物理试验的宏观响应相等原则，计算各类接触本构参数具体值.
3) 以焦家金矿-450中段下盘糜棱岩为例，基于单轴抗压物理试验，确定对应的数值计算模型微观参数并进行模拟试验，模拟结果与物理试验响应具有良好的一致性，证实了微观参数确定方法可靠并具有较好的适用性.

参考文献

[1]	Itasca Consulting Group.Online manual table of contents of particle flow code in 2 dimensions[EB/OL].[2015-02-10].http://www.itascacg.com/.(2)
[2]	Huang H,Lecampion B,Detournay E.Discrete element modeling of tool-rock interaction I:rock cutting[J].International Journal for Numerical and Analytical Methods in Geomechanics,2013,37(13):1913-1929.(2)
[3]	张社荣,孙博,王超,等.双轴压缩试验下岩石裂纹扩展的离散元分析[J].岩石力学与工程学报,2013,32(2):3083-3091.(Zhang She-rong,Sun Bo,Wang Chao,et al.Discrete element analysis of crack propagation in rocks under biaxial compression[J].Chinese Journal of Rock Mechanics and Engineering,2013,32(2):3083-3091.)(3)
[4]	Wang N C,Tannant D D,Lilly P A.Numerical analysis of the stability of heavily jointed rock slopes using PFC2D[J].International Journal of Rock Mechanics and Mining Sciences,2003,40(3):415-424.(2)
[5]	Fakhimi A,Carvalho F,Ishida T,et al.Simulation of failure around a circular opening in rock[J].Rock Mechanics and Mining Sciences,2002,39(4):507-515.(2)
[6]	Zhang S R,Sun B,Wang C,et al.Influence of intermediate principal stress on failure mechanism of hard rock with a pre-existing circular opening[J].Journal of Central South University,2014,21(4):1571-1582.(1)
[7]	Cai M,Kaiser P K,Moriok H,et al.FLAC/PFC coupled numerical simulation of AE in large-scale underground excavations[J].International Journal of Rock Mechanics & Mining Sciences,2007,44(4):550-564.(1)
[8]	王培涛,杨天鸿,朱立凯,等.基于PFC2D岩质边坡稳定性分析的强度折减法[J].东北大学学报(自然科学版),2013,34(1):127-130.(Wang Pei-tao,Yang Tian-hong,Zhu Li-kai,et al.Strength reduction method for rock slope stability analysis based on PFC2D[J].Journal of Northeastern University(Natural Science),2013,34(1):127-130.)(1)
[9]	Fakhimi A.Application of slightly overlapped circular particles assembly in numerical simulation of rocks with high friction angles[J].Engineering Geology,2004,74(1/2):129-138.(2)
[10]	Kentaro O,Masayasu O.Crack classification in concrete based on acoustic emission[J].Construction and Building Materials,2010,24(12):2339-2346.(1)
[11]	Caroline C G,Sergei S,Ian G M,et al.Comparison of polarity and moment tensor inversion methods for source analysis of acoustic emission data[J].International Journal of Rock Mechanics and Mining Sciences,2010,47(1):161-169.(1)

2

本文献在全文中的定位：
... 不必考虑划分尺寸的敏感性.基于离散单元法的颗粒流程序PFC^2D[1]能够较好地模拟岩石在加载过程中的局部各向异性行为，重现微裂隙时 ...

... K_s和T_n，T_s，u决定^[1] ...



2

本文献在全文中的定位：
... 时空演化过程和最终断裂区，已被广泛应用于研究岩石的各项力学行为^{[2, 3, 4,} ...

... 现如钢铁等延性材料的力学行为，大比率值与岩石等脆性材料性能吻合^[2].Fakhimi等^[9]选择 ...



3

本文献在全文中的定位：
... 泛应用于研究岩石的各项力学行为^{[2, 3, 4, 5]} ...

... 目前国内外学者普遍应用“试错法”^[3]来确定PFC^2D数值计算模型微观参数，存在费时费力及盲目性大 ...

... 而且多参数确定过程中，纯粹的试错极难同时获得精确的微观参数值^{[3, 4, 5,} ...



2

本文献在全文中的定位：
... ^{[2, 3, 4, 5]}. ...

... 错极难同时获得精确的微观参数值^{[3, 4, 5, 6]} ...



2

本文献在全文中的定位：
... ^{3, 4, 5]}. ...

... ^{[3, 4, 5, 6]}.在应用过程中，文献仅提及“试错 ...



1

本文献在全文中的定位：
... ^{4, 5, 6]}.在应用过程中，文献仅提及“试错法”，在参数确定环节直接给出具 ...



1

本文献在全文中的定位：
... 数确定环节直接给出具体值，未对具体确定方法和判别准则作详细论述^{[7, 8, 9]} ...



1

本文献在全文中的定位：
... 体确定方法和判别准则作详细论述^{[7, 8, 9]}. ...



2

本文献在全文中的定位：
... ^{[7, 8, 9]}. ...

... 脆性材料性能吻合^[2].Fakhimi等^[9]选择T_n/T_s=1:5来模拟孔洞岩石的力学行为，模 ...



1

本文献在全文中的定位：
... 进行单轴抗压试验，在加载过程中，利用声发射技术实时监测破裂事件^{[10, 11]}，加载系统如 ...



1

本文献在全文中的定位：
... 利用声发射技术实时监测破裂事件^{[10, 11]}，加载系统如图 4所示，试样尺寸及物 ...