柳小波, 王连成
东北大学 资源与土木工程学院, 辽宁 沈阳 110819
收稿日期：2016-09-04
基金项目：“十二五”国家科技支撑计划项目(2015BAB15B01)。
作者简介：柳小波(1980-)，男，辽宁丹东人，东北大学副教授，博士。

摘要：基于离散元思想, 结合随机理论, 以崩落矿岩散体为研究对象, 建立了非均匀散体碰撞运动模型.该模型主要研究内容包括:模型假设条件、散体颗粒的随机生成、颗粒单元的稳定性分析、颗粒单元的初始移动、颗粒单元的碰撞过程以及颗粒单元运动的随机性等.根据所建立的非均匀散体碰撞运动模型, 利用C#编程语言研发了二维放矿仿真系统.该系统包括参数输入、模型生成、仿真模拟和数据统计四个功能模块.设计了无底柱分段崩落法放矿模拟方案, 所得结论与类似文献对比, 基本相同, 说明利用非均匀散体碰撞运动模型模拟散体流动是可行的.
关键词：非均匀崩落矿岩散体流动放矿仿真随机理论
Simulation Model of Granular Flow of the Non-uniform Caved Ore and Rock
LIU Xiao-bo, WANG Lian-cheng
School of Resources & Civil Engineering, Northeastern University, Shenyang 110819, China
Corresponding author: WANG Lian-cheng, E-mail: 13898110153@163.com
Abstract: By combining the discrete element method with the stochastic medium theory, the collision movement model of non-uniform granular media is proposed, which includes random generation of particles, stability analysis of granular units, initial movement of granular unit, collision process of granular unit and randomness of movement of particles and so on. Based on the established collision movement model of non-uniform granular media, a 2D drawing simulation system was developed using C# language, which includes four modules of parameter input, model generation, emulational simulation, and data statistics. A pillarless sublevel drawing simulation scheme was designed and was found to be consistent with the study from relevant literatures, indicating it is feasible to simulate discrete granule flow using collision movement model of non-uniform granular media.
Key Words: non-uniformcaved ore and rockgranular flowore-drawing simulationrandom theory
崩落采矿法是一种通过崩落围岩来管理地压的采矿方法, 因其具有效率高、成本低、适用范围广等优点, 被广泛应用于国内外矿山.在崩落采矿法中, 矿石在废石覆盖下进行放矿, 很容易造成贫化.矿石的损失贫化既浪费了国家的宝贵资源, 降低了矿山经济效益, 又加大了生产成本；而研究崩落矿岩散体流动规律是进行放矿贫化研究的基础.因此, 研究矿岩散体的流动规律, 实现放矿的仿真模拟, 优化采矿方法的相关参数, 对提高矿石回收率、提高矿山经济效益具有重要意义.
目前研究散体流动的理论或模型主要有离散元法、随机介质理论、九块模型等.基于离散元法编制的颗粒流软件PFC^2d及PFC^3d是根据牛顿第二定律等一系列复杂的力学计算来模拟颗粒体流动的^[1-4]; 1962年王泳嘉教授提出了散体移动的球体递补模型, 是随机介质理论在散体流动研究方面的重大突破; 1968年, 加拿大学者Jolly提出的散体移动九块模型在放矿仿真方面取得了重大进展^[5], 但九块模型的均匀性假设使其无法模拟非均匀散体的真实流动过程; 除此之外, 文献[6-7]中都对非均匀散体流动仿真模型提出了独到的见解.
本文以圆形颗粒为单元, 建立了非均匀散体碰撞运动模型, 以此实现散体流动过程的模拟.该模型结合了离散元思想和随机理论, 可以简单实用地模拟散体流动.该模型可以模拟大小不同的散体颗粒的流动问题, 详细地展现散体颗粒的运动过程, 兼有随机现象, 具有很好的发展前景.
1 非均匀散体碰撞运动模型的建立1.1 假设条件① 崩落矿岩散体为半径不同的圆形颗粒单元; ② 崩落矿岩散体颗粒为不可变形的刚体; ③ 不考虑散体碰撞过程中的变形和滑动, 即采用硬球模型; ④ 当多个颗粒满足移动条件时, 颗粒的移动先后具有随机性; ⑤ 不考虑摩擦力的影响, 球体颗粒之间的运动为滑动, 且满足动能定理; ⑥ 由于颗粒受自重力的影响, 不能向上运动.
1.2 散体颗粒的随机生成散体颗粒生成步骤如下.
1) 在一定矩形范围内(起点(x₀, y₀), 边长为a, b), 利用蒙特卡洛随机投点法^[8]生成一些点坐标(x_i, y_i), i=1, 2, …, n, 并以这些坐标为圆心, 以输入值r为半径生成圆, 其生成规则如下:
规则1:代表颗粒的圆不能超越边界, 即

(1)

规则2:代表颗粒的圆之间不能相交, 即

(2)

此时得到的结果见图 1a, 需要进一步压实(图 1b).
图 1(Fig. 1)

图 1 未压实和压实的散体颗粒Fig.1 Uncompacted and compacted particles (a)—未压实颗粒；(b)—压实颗粒.

2) 由于散体颗粒受重力作用, 因此处于不稳定状态的散体颗粒会在重力作用下产生位移, 散体颗粒移动的过程中遵循规则1和规则2, 直到所有的散体颗粒全部达到稳定状态.具体的流动过程见下文.
3) 重复以上操作, 当每一颗粒都达到稳定状态, 即如图 1b所示, 生成完毕.
崩落法放矿是覆盖下放矿, 因此采用两种不同灰度的圆盘代表矿石和废石, 其中浅色圆盘代表废石, 深色圆盘代表矿石, 两种圆形单元分两次形成并在重力下压实，完成初始状态的模拟.
1.3 颗粒单元的稳定性分析在重力的作用下, 散体颗粒具有向下移动的趋势, 能否产生移动需看其下方是否有阻碍.当要移动的颗粒下方有阻碍时, 则需要利用力学判据对其可否移动、运动的方向和过程作进一步判断.其判断过程如下.
以颗粒圆心为坐标原点建立直角坐标系, 将整个区域划分成四个象限.根据力的合成方法, 若合力为零, 则认为该颗粒单元处于稳定状态.假设颗粒的第三象限有其他单元或边界阻碍其运动, 则认为该颗粒在阻碍单元未移动之前不会向该方向移动, 第四象限则同理.如果颗粒第三象限和第四象限都有阻碍, 则认为该颗粒是稳定的.因此, 首先通过探索判断颗粒是否有接触点.以建立的坐标系为基础, 当颗粒的周围存在一个接触点时, 若接触点在y轴上时, 该颗粒单元处于稳定状态(如图 2所示), 否则不稳定; 当颗粒第三、四象限存在接触点时, 该颗粒稳定(见图 3a), 当颗粒第三象限和第一象限存在接触点时, 若颗粒1圆心与颗粒2、颗粒3圆心连线在逆时针方向所成夹角小于180°, 该颗粒稳定, 否则不稳定(见图 3b); 当颗粒第三象限存在接触点且右侧与边界接触时, 该颗粒稳定(如图 3c); 当颗粒第二象限和第四象限存在接触点时, 判断方法同理.
图 2(Fig. 2)

图 2 颗粒有一个接触点Fig.2 Particle has one point in contact

图 3(Fig. 3)

图 3 颗粒的接触状态Fig.3 Contact status of particle (a)—颗粒1在第三、四象限存在接触点；(b)—颗粒1在第一、三象限存在接触点; (c)—颗粒1靠近边界.

1.4 颗粒单元的初始移动处于不稳定状态的颗粒单元, 必然会发生移动.对于下方没有阻碍的颗粒, 颗粒将做自由落体运动, 而对于有接触阻碍的颗粒, 情况则如图 4所示, 不考虑颗粒之间的摩擦, 颗粒1由静止状态开始绕颗粒2做圆周运动, 初速度为v₀(v₀=0), 运动过程中由重力mg的分力F_n提供向心力, 当某一刻F_n不足以提供做圆周运动的向心力时, 颗粒1脱离颗粒2做抛物运动.
图 4(Fig. 4)

图 4 颗粒的初始移动分析Fig.4 Analysis of initial movement of particles

颗粒做圆周运动时, 有以下公式成立:

(3)

(4)

(5)

式中:m为颗粒单元1的质量; F_n为重力的分力提供的向心力; h为颗粒单元1的垂直方向位移; r为颗粒单元1的半径.
通过分析可得, 当

(6)

成立时, 即当

(7)

颗粒1脱离颗粒2，作抛物运动，此时速度为

(8)

1.5 颗粒单元的碰撞过程崩落的矿岩散体在重力下放出, 当其具有一定的速度并与其他颗粒接触时必然会发生碰撞, 其碰撞类型主要包括:与边界发生碰撞, 与稳定颗粒发生碰撞, 不稳定颗粒之间发生碰撞.下面以不稳定颗粒之间发生碰撞(见图 5)为例进行推导, 其过程如下.
图 5(Fig. 5)

图 5 不稳定颗粒之间发生碰撞Fig.5 Collision between unstable particles

如图 5所示, 颗粒单元a以速度v_a¹与速度为v_b¹的颗粒单元b发生碰撞.根据模型假设, 在不考虑能量损失的情况下, 认为颗粒之间发生完全弹性碰撞.两颗粒圆心连线即为法线方向.将速度v_a¹和v_b¹分解为法线方向和切线方向, 不考虑滑动的情况下, 切向速度保持不变.由于不考虑能量损失、摩擦力和转动等复杂情况, 则根据动量守恒定律和能量守恒定律, 有以下式子成立:

(9)

(10)

联立以上两式可求得碰撞后的速度为

(11)

(12)

根据矢量合成定律, 即可求得碰撞后速度的大小和方向.碰撞后的颗粒在重力作用下不断改变速度方向和大小, 经过一系列碰撞后最终达到稳定或放出漏口.
1.6 颗粒单元运动的随机性对于同时满足运动条件的散体颗粒来说, 它们的运动必然具有先后顺序.本文认为, 颗粒单元的运动顺序具有随机性, 先运动的颗粒速度大于后运动的颗粒, 而相同的时间间隔也必然运动得更远.如图 6所示, 颗粒1和颗粒2同时满足运动条件, 则二者的运动顺序具有随机性, 若颗粒1先运动到虚线位置, 必然阻碍颗粒2的运动, 由此颗粒1较颗粒2先行向下运动.在概率赋值的过程中, 半径较小的颗粒先运动的概率也应大于半径较大的颗粒.通过对不稳定颗粒进行有规律的随机概率赋值, 便可以得到随机的散体颗粒运动顺序, 以此达到随机理论与离散运动相结合的目的.
图 6(Fig. 6)

图 6 颗粒单元运动的随机性Fig.6 Randomness of movement of granular

2 放矿仿真系统的开发与实现2.1 功能模块划分根据功能需求, 将放矿仿真系统分成四个功能模块, 分别为:输入参数模块、模型生成模块、模拟流动模块和数据统计模块，如图 7所示.
图 7(Fig. 7)

图 7 仿真系统功能模块划分Fig.7 Function module division of simulation system

2.2 界面设计根据功能模块的划分, 设计二维放矿仿真系统包括四个部分, 分别是参数输入、程序控制部分、数据统计部分和流动过程实时展示区域；界面效果如图 8所示, 黑色矩形代表放矿进路.
图 8(Fig. 8)

图 8 程序界面设计图Fig.8 Design of program interface

3 放矿仿真系统的实验对比3.1 无底柱分段崩落法结构参数优化概述无底柱分段崩落法具有生产能力大、结构简单、机械化程度高，以及生产作业安全等优点, 应用范围十分广泛.长期以来, 由于采矿设备的限制, 我国的无底柱分段崩落法采用的结构参数偏小, 采切比过高致使生产成本居高不下, 而生产能力也无法得到有效提高, 严重制约了我国无底柱分段崩落法的发展.通常情况下, 对矿石损失贫化影响较大的结构参数包括分段高度、进路间距和崩矿步距, 三者是互相联系和互相制约的.因此, 设计合理的结构参数具有十分重要的意义.
3.2 模拟放矿方案选择本次实验选取了分段高度和进路间距两个因素, 共设计9组实验, 对放矿过程进行仿真模拟, 以达到获得最优结构参数的目的.本文参照文献[9]的实验参数设计了放矿方案, 具体的实验方案及实验结果见表 1.
表 1(Table 1)

表 1 实验方案及实验结果Table 1 Experimental scheme and results 方案 进路间距 /m 分段高度 /m 矿石回收率/% 废石混入率/% 1 13 10 40.1 3.0 2 13 14 43.8 3.4 3 13 18 41.6 3.5 4 16 10 44.6 2.9 5 16 14 49.1 1.8 6 16 18 46.0 1.6 7 19 10 42.4 2.1 8 19 14 43.4 0.7 9 19 18 40.0 1.7

表 1 实验方案及实验结果 Table 1 Experimental scheme and results

本次模拟采用参数如下:分段个数2个, 进路个数2个, 岩石层厚10 m, 单位体积岩石混入率33%, 进路尺寸3 m×3 m, 粒径分布符合高斯分布, 半径0~30 cm占30%, 30~40 cm占40%, 40~50 cm占30%.由于第一个分段为不正常分段, 因此, 统计数据以第二分段为准.
3.3 模拟结果分析根据模拟得到的实验结果, 利用Matlab建立分段高度H和进路间距B两个因素与矿石回收率φ的回归方程^[10]:

(13)

利用该回归方程预估各种参数条件下的矿石回收率, 说明矿石回收率与分段高度和进路间距的变化关系.当分段高度一定时, 即可分析进路间距对矿石回收率的影响, 如图 9c所示.文献[9]室内放矿实验和PFC数值模拟结果如图 9a和9b所示.当进路间距一定时, 实验结果如图 10所示.
图 9(Fig. 9)

图 9 进路间距对矿石回收率的影响Fig.9 Effect of drift interval on ore recovery rate (a)—物理实验; (b)—PFC数值模拟; (c)—设计方案模拟.

图 10(Fig. 10)

图 10 分段高度对矿石回收率的影响Fig.10 Effect of sublevel height on ore recovery rate (a)—物理实验; (b)—PFC数值模拟; (c)—设计方案模拟.

根据实验数据分析, 可以得到以下结论:
1) 当分段高度一定时, 进路间距对矿石回收率的影响基本呈现二次曲线的变化关系, 且呈先增大再减小的变化趋势.
2) 当进路间距一定时, 分段高度对矿石回收率的影响基本呈现二次曲线的变化关系, 且呈先增大再减小的变化趋势.
3) 在本实验设计的参数下, 当分段高度为14 m, 进路间距为16 m时, 矿石的回收率最高.
4) 本文得到的结论与文献[9]研究所得结论基本一致, 这表明本研究是合理的.
4 结论1) 本文将随机理论和离散元思想相结合, 建立了一种研究散体流动的新方法——非均匀散体碰撞运动模型.该方法既能模拟不同粒径的颗粒在重力作用下相互碰撞的流动情况, 又能模拟非均匀散体的流动, 同时还具有随机的特点, 是一种研究非均匀散体流动的全新尝试.
2) 基于非均匀散体碰撞运动模型, 本文通过对关键方法的设计和系统模块的划分, 设计研发了二维放矿仿真系统，并参照相关文献设计了无底柱分段崩落法放矿实验；得到的结果与该文献研究结论基本一致, 说明利用非均匀散体碰撞运动模型进行放矿实验模拟是可行的.
参考文献

[1]	Han K, Feng Y T, Owen D R J. Sphere packing with a geometric based compression algorithm[J].Powder Technology, 2005, 155(1): 33–41.DOI:10.1016/j.powtec.2005.04.055
[2]	Liu Z, Zhou N, Zhang J X. Random gravel model and particle flow based numerical biaxial test of solid backfill materials[J].International Journal of Mining Science and Technology, 2013, 23(4): 463–467.DOI:10.1016/j.ijmst.2013.07.001
[3]	Brady B H G, Brown E T. Rock mechanics for underground mining[M]. London: Allen & Unwin, 1985: 334-339.
[4]	Itasca Consulting Group Inc. Particle flow code (version 3.0)[M]. Minneapolis: ICG, 2004: 22-26.
[5]	Jolly D. Computer simulation of the movement of ore and waste in an underground mining pillar[J].The Canadian Mining and Metallurgical, 1968, 67(2): 854–859.
[6]	李昌宁. 非均匀矿岩散体放矿的计算机模拟[J].有色金属, 2002, 54(2): 98–103. ( Li Chang-ning. Computer simulation of heterogeneous loose body drawing[J].Nonferrous Metals, 2002, 54(2): 98–103.)
[7]	柳小波. 散体流动时空演化仿真模型的研究与应用[D]. 沈阳: 东北大学, 2009. ( Liu Xiao-bo.Study on simulation model of time-space evolution of granular flow and its applications[D]. Shenyang:Northeastern University, 2009.http://cdmd.cnki.com.cn/Article/CDMD-10145-2010256950.htm)
[8]	Metropolis N, Ulam S. The Monte Carlo method[J].Journal of the American Statistical Association, 1949, 247(44): 115–129.
[9]	吴爱祥, 武力聪, 刘晓辉. 无底柱分段崩落法结构参数研究[J].中南大学学报(自然科学版), 2012, 15(5): 1845–1850. ( Wu Ai-xiang, Wu Li-cong, Liu Xiao-hui. Study on structural parameters of sublevel caving[J].Journal of Central South University (Science and Technology), 2012, 15(5): 1845–1850.)
[10]	Sen S K, Shaykhian G A. Matlab tutorial for scientific and engineering computations[J].Nonlinear Analysis Theory Methods and Applications, 2009, 71(12): e1005–e1020.DOI:10.1016/j.na.2009.01.069