岩石在加载过程中的力学行为和破裂时空演化规律研究是地下岩体稳定性分级、采矿方法选择及采场结构参数设计的基础.目前用来模拟岩石损伤与变形行为的数值模拟方法主要包括连续介质模型和离散单元法(DEM).相对于连续介质模型而言,一方面,DEM借助模拟微观单元体间的力学行为来叠加生成宏观响应,不需要对模型预设复杂的数学本构关系;另一方面,DEM由于不涉及网格的划分,不必考虑划分尺寸的敏感性.基于离散单元法的颗粒流程序PFC2D[1]能够较好地模拟岩石在加载过程中的局部各向异性行为,重现微裂隙时空演化过程和最终断裂区,已被广泛应用于研究岩石的各项力学行为[2, 3, 4, 5].
在应用PFC2D解决特定问题之前,需要确定颗粒集合微观参数,包括颗粒尺寸分布、颗粒间接触刚度和强度等,计算模型微观参数的正确取值是获取精确模拟结果的关键.然而,岩石的力学行为通常借助如弹性变形常数、抗压强度和残余强度等宏观参数来描述.由于数值模拟微观参数和物理试验宏观响应的关系是非线性的,因而无法借助实验室所得宏观响应直接计算模拟所需的微观参数值.
目前国内外学者普遍应用“试错法”[3]来确定PFC2D数值计算模型微观参数,存在费时费力及盲目性大等问题,而且多参数确定过程中,纯粹的试错极难同时获得精确的微观参数值[3, 4, 5, 6].在应用过程中,文献仅提及“试错法”,在参数确定环节直接给出具体值,未对具体确定方法和判别准则作详细论述[7, 8, 9].
针对“试错法”存在的盲目性大和精确度低等问题,本文提出基于物理试验响应的PFC2D微观参数确定方法,以焦家金矿-450中段下盘糜棱岩为例,验证具体方法的可靠性和适用性.
1 方法详述颗粒流程序PFC2D数值模型微观参数包括物理几何参数和接触本构参数:物理几何参数通过具体模拟对象的基本特性确定,接触本构参数需要借助物理试验响应进行反演计算.
1.1 物理几何参数物理几何参数包括模拟试样高度H、宽度D、颗粒单元体最大半径Rmax、最小半径Rmin、粒径分布规律、颗粒集合孔隙率n、颗粒密度ρ,确定流程如图 1所示,具体过程如下.
1) H和D:按照模拟试样与物理试验岩样尺寸一致原则,通过测量方式确定H和D.
2) Rmax,Rmin、粒径分布规律和n:借助电镜扫描仪或显微仪器来观察试验岩样,通过绘制手段,获取统计层面每个晶粒的周长、面积、最大半径、最小半径、筛网尺寸、实际尺寸和误差等信息,绘制晶粒粒级分布曲线,确定Rmax,Rmin、粒径分布规律和n.
3) ρ:称取试验岩样质量、体积,去除孔隙所占比例,计算参数ρ.
1.2 接触本构参数PFC2D数值计算模型接触本构参数主要包括颗粒间接触模量Ec、正向切向刚度比Kn/Ks、正向强度Tn、切向强度Ts、正向强度偏差Tdn、切向强度偏差Tds、颗粒单元体之间的摩擦系数u.颗粒间的接触本构模型如图 2所示,A,B代表颗粒单元体,可以看出,颗粒集合刚度和强度分别由Kn,Ks和Tn,Ts,u决定[1].
借助单轴或常规三轴物理试验,获取岩石弹性模量E,泊松比μ,峰值强度бc,裂隙初始强度бci及峰后残余强度бca.在保证边界、加载及截止条件等与物理试验一致的基础上,开展数值模拟试验.通过按顺序不断更改各接触本构参数的方式,建立数值模拟微观参数与宏观响应之间的对应关系,按照数值模拟与物理试验宏响应相等原则,计算各类接触本构参数具体值.确定流程如图 3所示,具体过程可如下.
1) Ec和Kn/Ks:将模拟材料强度参数Tn,Ts设置为较大值(大约10倍的单轴抗压强度),改变Ec,建立Ec和E之间的关系,借助已知的E计算Ec.同理,改变Kn/Ks的取值,确定μ与Kn/Ks的关系,而Kn/Ks的改变会反过来影响E值,因此,需要借助迭代手段来确定微观参数Ec和Kn/Ks,其中迭代方法选用二分法.
2) Tn,Ts,Tdn和Tds:设置Tdn,Tds为0,通过改变Tn,Ts来建立бc与Tn,Ts的关系.值得注意,Tn/Ts的取值会影响矫正结果,必须保证Tn/Ts不变.数值模拟过程中,Tn/Ts值的大小控制了剪切和张破坏微裂隙生成数量,进而影响最终的破裂模式.一般情况下,小比率值可体现如钢铁等延性材料的力学行为,大比率值与岩石等脆性材料性能吻合[2].Fakhimi等[9]选择Tn/Ts=1:5来模拟孔洞岩石的力学行为,模拟结果与实际岩石物理试验一致,本文确定Tn/Ts=1:5.保证Tds/Ts=Tdn/Tn,同时增加Tdn,Tds,确定бci与Tdn,Tds的关系,裂隙初始生成强度可借助如声发射等其他手段获得.而此过程会影响бc,同样应用二分法迭代2)、以获参数Ts,Tn,Tdn和Tds.
3) u:离散元软件PFC2D在数值模拟过程中,颗粒间的破裂服从库伦准则,黏结破裂之前,材料强度由Tn,Ts决定,黏结破裂之后,材料强度主要由u确定.通过改变u的取值,建立岩石峰后行为与u的关系,通过数值模拟与物理试验的峰后残余强度值бca相等原则确定u.
2 实例分析基于焦家金矿-450中段下盘糜棱岩单轴抗压物理试验,确定对应的PFC2D数值计算模型微观参数,通过对比单轴、双轴模拟试验结果和物理试验响应,分别验证确定方法的可靠性和适用性.
2.1 物理试验简述将现场钻取的岩芯按照工程岩体试验标准进行单轴抗压试验,在加载过程中,利用声发射技术实时监测破裂事件[10, 11],加载系统如图 4所示,试样尺寸及物理参数、加载条件和试验结果见表 1.
按照PFC2D微观参数确定方法,获取数值计算模型微观参数,最终的物理几何和接触本构参数见表 2.糜棱岩的电镜扫描结果见图 5.可以看出,岩样晶粒辨析度较差,仅能从局部区域观察晶粒和孔隙大小及分布.为了尽量精确确定粒径分布规律,对同批试样的上下表面均进行扫描,将统计的整体结果作为Rmax,Rmin、粒径分布规律和n的最终取值.
借助确定的微观参数建立数值模型并进行单轴抗压模拟试验,将数值试验的应力应变曲线和最终破裂模式与物理试验响应进行对比,结果如图 6所示.图 6a可知,两者应力应变曲线在峰值强度之前基本重合,而峰值强度之后的差异较大,原因在于PFC2D是借助模拟微观单元体之间的力学行为来叠加生成宏观响应,不存在轴向应力的陡然下降.同时,由数值模拟裂隙数演化曲线可知,模拟试验的бci为51.4 MPa,与实验室声发射首个信号接收时的平均轴向应力值52.1 MPa几乎相等.图 6b~图 6c显示了岩石数值模拟和物理试验最终加载结束后的失稳破裂模式,图 6d展示了岩石在破裂过程中的微裂隙时空演化过程,圆圈代表微裂隙,白色至黑色代表微裂隙生成的时间前后顺序.可以看出,糜棱岩在单轴受载条件下,剪切断裂带主要位于岩石中央,方向大约与轴向应力成20°角.
借助矫正参数建立数值模型,实施不同围压下的双轴模拟试验,将模拟岩石不同围压下的强度值和最终破裂模式与常规三轴抗压物理试验响应进行对比.由强度曲线对比结果可知,实验室不同围压下的岩石强度值虽然存在一定的离散性,但其拟合曲线基本与数值模拟曲线重合,见图 7a.同时,图 7b展示了10 MPa围压条件下,岩石物理试验与数值模拟试验的最终破裂模式(其他围压下的结果一致),结果吻合性较好.双轴模拟试验与物理试验对比结果表明,基于物理试验的PFC2D微观参数确定方法具有较强的适用性,能够用来计算其他边界条件下的岩石力学行为.
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) |