基于PFC2D-DFN的自然崩落法数值模拟研究

doi:10.11872/j.issn.1005-2518.2019.02.189

基于PFC2D-DFN的自然崩落法数值模拟研究

方传峰^,¹^,²^,³, 王晋淼¹^,⁴, 李剡兵⁵, 贾明涛¹^,⁴

1. 中南大学资源与安全工程学院，湖南长沙 410083

2. 长沙矿山研究院有限责任公司，湖南长沙 410012

3. 金属矿山安全技术国家重点实验室，湖南长沙 410012

4. 中南大学数字矿山研究中心，湖南长沙 410083

5. 云南迪庆有色金属有限责任公司，云南香格里拉 674400

Numerical Simulation Research of Natural Caving Method Based on PFC2D-DFN

FANG Chuanfeng^,¹^,²^,³, WANG Jinmiao¹^,⁴, LI Shanbing⁵, JIA Mingtao¹^,⁴

1. School of Resources and Safety Engineering，Central South University，Changsha 410083，Hunan，China

2. Changsha Institute of Mining Research Co. ，Ltd. ，Changsha 410012，Hunan，China

3. State Key Laboratory of Safety Technology of Metal Mines，Changsha 410012，Hunan，China

4. Center of Digital Mine Research，Central South University，Changsha 410083，Hunan，China

5. Yunnan Diqing Nonferrous Metal Co. ，Ltd. ，Shangri-La 674400，Yunnan，China

收稿日期: 2018-02-01 修回日期: 2018-05-30 网络出版日期: 2019-04-29

基金资助:

国家重点研发计划项目“深部集约化开采生产过程智能管控技术”. 编号：2017YFC0602905

Received: 2018-02-01 Revised: 2018-05-30 Online: 2019-04-29

作者简介 About authors

方传峰（1990-），男，山东泰安人，硕士研究生，从事矿山开采数值模拟分析研究工作15084850237@163.com , E-mail：15084850237@163.com

摘要

采用二维颗粒流软件PFC2D构建矿山模型，引入符合矿山实际节理产状的DFN网络，基于室内岩石试验及矿山实际崩落情况，标定模拟所用细观参数。以某工程应用为实例，模拟分析崩落破坏机理与崩落区演化趋势。结果表明：裂纹发起于推进面与顶板交界处，在拱形区域内密集扩展、贯通，矿岩破裂以节理拉伸破坏为主，矿堆以岩石剪切破坏为主；崩落区前期呈拱形演化，崩透地表后引起邻近高陡边坡滑动；部分拉底仅造成新开挖顶板上方局部崩落。该研究结果可对矿山崩落采矿起到指导作用。

关键词： 自然崩落法 ; 数值模拟 ; PFC2D软件 ; DFN ; 拉底 ; 岩石破裂 ; 力学分析 ; 崩落演化

Abstract

Surface subsidence and rock burst may occur during mining in mines using natural caving method，which affect the safety of mine production and the life safety of underground workers.Therefore，it is necessary to research the evolution law of caving area with the advance of undercutting.The research process of this paper is as follows.Firstly:using two-dimensional discrete element software PFC2D to construct a geometric model consistent with the actual topography of a mine.Secondly，according to the complete rock parameters obtained from laboratory mechanical experiments，by means of simulating physical experiments，the meso-parameters are continuously adjusted to approximate the actual physical parameters，the meso-mechanical parameters used in PFC2D are calibrated.In addition，according to the statistic law of joint sets distribution obtained on the spot，choose an appropriate joint generation method，make use of discrete fracture network DFN to reconstructe a joint model consistent with the statistical law.Then，based on the caving height of a certain period obtained by mine drilling television，the simulation results are approximated to the monitoring values by dichotomy method，and the corresponding joint strength parameters are back analyzed.Finally，the whole mine model is coupled with the joint model，and the particles and joints are assign property.After stress balance，the excavation is simulated.Statistic simulation results，the conclusions are as follows:Cracks first develop at the intersection of the undercutting advance working face and the roof，then extend，develop and penetrate in the arch region，resulting in crushing and collapse of the original rock.The main failure types of ore and rock are joint tensile failure，while shear failure often occurs in caved ore and rock.In the early stage，the caving area develops arch upward，which causes the adjacent steep slope to slide to the valley after penetrating the surface.In some undercutting steps，the caving area only occurs in the local area above the newly excavated undercutting region，which is consistent with the change trend of caving quality and roof height.

Keywords： natural caving method ; numerical simulation ; PFC2D software ; DFN ; undercutting ; rock failure ; stress analysis ; caving evolution

PDF (9888KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

方传峰, 王晋淼, 李剡兵, 贾明涛. 基于PFC2D-DFN的自然崩落法数值模拟研究[J]. 黄金科学技术, 2019, 27(2): 189-198 doi:10.11872/j.issn.1005-2518.2019.02.189

FANG Chuanfeng, WANG Jinmiao, LI Shanbing, JIA Mingtao. Numerical Simulation Research of Natural Caving Method Based on PFC2D-DFN[J]. Gold Science and Technology, 2019, 27(2): 189-198 doi:10.11872/j.issn.1005-2518.2019.02.189

自然崩落法自1895年应用于Pewabic铁矿起，经过不断研究，在矿体可崩性方面取得了巨大进展，已突破最初使用条件，并不断发展成熟^[1,2]。该方法是在拉底空间上方，利用矿体自身软弱结构面，在其重力及次生构造应力作用下，诱导矿体失稳崩落，从而省却大量凿岩爆破工程及费用^[3]。自然崩落法是一种低成本、高效率且安全性好的采矿方法，得到广泛的推广应用^[4]。

然而，自然崩落法也存在地表破坏严重、井下气浪冲击等缺陷^[5]，有必要研究自然崩落的影响因素与发展规律，研究方法包括经验推导法、观测法、相似物理试验法和数值模拟法等^[6,7,8,9]。随着计算机技术的发展，数值模拟法已成为国内外学者研究自然崩落力学响应特征的有效手段。王涛等^[10]采用PFC2D分析矿岩在崩落过程中的裂纹扩展及接触力分布情况；王连庆等^[11]采用PFC2D数值模拟方法分析了自然崩落法的崩落规律；曾庆田等^[12]采用FLAC^3D与PFC2D分析了多斗均匀放矿与放矿步距的合理性；Vyazmensky等^[13]根据实例统计、经验图分析和有限—离散元耦合模拟，研究了影响地表沉降崩落角的多种因素；Woo等^[14]将FLAC^3D与InSAR技术相结合，评价、预测矿块崩落法引起的地表变形，并修正模拟参数。

综上所述，研究自然崩落力学响应特征的数值模拟主要借助于连续介质算法、离散介质算法或二者耦合^[15]。与连续介质算法相比，离散介质算法可直观、形象、动态地模拟、跟踪矿体开裂和崩落过程，分析崩落过程的力学响应特征，具有更广泛的应用前景^[16]。传统离散元模拟方法未将对崩落影响巨大的随机节理网络引入模型，仅将颗粒自身作为崩落单元或将节理简化为若干组平行面，难以真实体现矿山构造。本文以西南地区某铜矿为例，采用PFC2D软件建立模型，参考节理分布规律引入随机节理网络，依据岩石试验数据及矿山监测数据标定岩石及节理细观力学参数，并基于实际拉底进度模拟开挖，研究矿山自然崩落的矿岩破裂机理与崩落规律。

1 矿山模型构建

1.1 整体矿山模型

根据地形等高线生成某矿山地形几何集，利用Dimine建模软件的“线附着”功能生成研究范围的地形轮廓线，见图1，对轮廓线进行“光顺”和“抽稀”等处理操作。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 矿山地形及地形轮廓线

Fig.1 Mine topography and terrain contour

轮廓线保存为DXF格式，将其导入PFC2D转化为墙体Wall，在稍大于模型范围建立矩形墙体，墙体内部均匀填充颗粒，利用PFC2D的放射法判断颗粒位置，删除模型以外的颗粒及辅助墙体，生成包含270万颗粒的矿山模型。3720水平生成线状几何体，模拟该水平坑道，3736-3751水平从模型中心向两帮延展，每30 m分为一拉底步，模型见图2。

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 PFC2D矿山模型

Fig.2 PFC2D mine model

1.2 节理网络生成

3720水平坑道调查获取节理相关数据，筛选出优势节理组，并统计节理倾角、倾向、间距和连续度的分布规律，矿山节理分为2个主节理组（1、2组）和2个次节理组（3、4组），这4组节理共占节理总数的78.5%，各节理组的产状和连续度均服从正态分布，节理间距呈对数正态分布，具体结果见表1和表2。

表1 优势节理组产状统计结果

Table 1 Statistical results of dominant joint set occurrence

组号	占比/%	倾向/（°）						倾角/（°）
组号	占比/%	最大值	最小值	均值	方差	分布规律	最大值	最小值	均值	方差	分布规律
1	36.36	166	120	143.17	126.89	正态分布	41	79	71.78	21.75	正态分布
2	32.74	196	239	216.21	95.13	正态分布	44	79	71.86	19.83	正态分布
3	5.41	329	306	317.07	35.58	正态分布	51	75	69.22	21.58	正态分布
4	3.99	31	48	39.99	23.69	正态分布	56	77	70.51	28.00	正态分布

新窗口打开| 下载CSV

表2 节理间距与连续度统计结果

Table 2 Statistical results of joint spacing and continuity(m)

项目	最小值	最大值	均值	标准差	分布规律
节理间距	0.00	1.38	0.093	0.15	对数正态分布
节理连续度	0.05	10.61	2.350	0.37	正态分布

新窗口打开| 下载CSV

DFN是模拟离散裂隙网络的概率统计分布函数，其结合节理密度、随机种子可生成用于研究节理岩体的DFN实体^[17]。PFC2D中DFN生成步骤为：依据矿山产状及连续度数据生成若干组节理样本，依据节理调查所用方法选择恰当的节理密度指标生成节理模型。PFC2D/3D包含P32、P21、P10和P10geom等节理密度指标。其中P10geom的机理为：当随机生成的节理与某条指定直线单位长度的交点个数达到设定值时，停止节理生成。节理数据由坑道调查所得，因此节理生成密度指标以P10geom为宜。

考虑计算机运行内存，兼顾体现节理空间形态，将实测节理间距扩大至1.0 m，并假设节理在空间均匀分布，4组节理的三维倾向与二维倾角的对应关系按照式（1）换算：

\{\begin{array}{l} d i p = d i p_{0} & d i p d \leq 180 ° \\ d i p = 180 ° - d i p_{0} & d i p d > 180 ° \end{array}

(1)

式中：dip为PFC2D所用倾角；dip₀为实测倾角；dipd为实测倾向。

此外，另有21.5%的节理不属于上述4组节理，对该部分节理按倾角均匀分布处理。各组节理与坑道单位长度交点数（即P10geom指标）按照式（2）换算。DFN各组节理倾角参数及P10geom指标见表3。

表3 DFN产生所用参数

Table 3 Parameters of DFN generated

节理分组	节理与坑道单位长度交点数/个	DFN节理倾角/（°）
节理分组	节理与坑道单位长度交点数/个	最小值	最小值	均值	标准差
1	0.36	40	70	71.78	4.67
2	0.33	100	136	108.14	4.45
3	0.05	104	130	110.78	4.64
4	0.04	56	77	70.52	5.29

新窗口打开| 下载CSV

n_{i} = \frac{\frac{L}{s p} \times d_{i}}{L} = \frac{d_{i}}{s p}

(2)

式中：L为坑道长度（m）；n_i为第i组节理与坑道的单位长度交点个数；d_i为第i组节理的数量百分占比；sp为节理间距（m）。

共生成节理约8万条，节理空间分布如图3所示。将节理与线状几何体交点的X坐标数据导出，统计模型的节理间距分布规律，确认其与实际情况相符，见图4。

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 DFN生成离散节理局部放大图

Fig.3 Local enlarged diagram of discrete joints generated with DFN

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 DFN节理间距统计

Fig.4 Statistical chart of DFN spacing

2 模拟参数标定与设置

PFC2D数值模拟共选用3种本构模型：线性本构模型用于墙与球之间的接触，且将二者之间的摩擦系数设置为零，作为滚动支撑；平滑节理模型不考虑沿节理局部颗粒接触的方向，节理两侧颗粒间的接触调用平滑节理模型；线性平行接触模型可同时传递力与力矩，提供了在2个接触实体之间类似混凝土材料的力学行为，用于模拟完整岩石。各本构模型力学计算原理详见文献[18]。

2.1 完整岩体细观参数标定

颗粒细观参数与宏观参数无法直接对应，需大量模拟与变形、强度有关的数值试验，并将其结果与室内试验结果对比，最终反分析标定 ^[19]。根据室内试验所得岩石物理力学参数，选用线性平行接触模型不断调整细观参数，使模拟结果与试验数据吻合。

模拟单轴压缩试验，调整接触模量Ec、平行粘结接触模量pb_Ec，获取岩石弹性模量E，调整平行粘结抗拉强度pb_ten，获取单轴抗压强度，改变刚度比 $K_{n} / K_{s}$ 、 ${\bar{K}}_{n} / \bar{K_{s}}$ ，调整岩石泊松比 $υ$ 。获取应力—应变曲线[图5（a）]，标定岩石单轴抗压强度及弹性模量。处理岩石破坏前的应变部分[图5（b）]，标定岩石泊松比。

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 单轴压缩试验曲线图

Fig.5 Plots of uniaxial compression test

在颗粒赋予的微观参数与单轴压缩试验一致的条件下，微调平行粘结黏聚力pb_coh、平行粘结摩擦角pb_fa，模拟三轴伺服压缩试验，试验曲线图如图6所示。

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 三轴伺服压缩试验曲线图

Fig.6 Plots of triaxial servo compression test

根据轴向应力 $σ_{1}$ 及相应侧向应力 $σ_{3}$ ，以 $\frac{σ_{1} + σ_{3}}{2}$ 为横坐标， $\frac{σ_{1} - σ_{3}}{2}$ 为纵坐标，在直角坐标系上用最小二乘法绘制最佳关系曲线。在最佳关系曲线上选定若干组对应值，以 $\frac{σ_{1} + σ_{3}}{2}$ 为圆心， $\frac{σ_{1} - σ_{3}}{2}$ 为半径，绘制应力莫尔圆，最终确定岩石抗剪强度参数。最终模拟结果与实测值对比（表4）可知：二者误差较小，细观参数能较好地反映完整岩石宏观力学特性，具体细观参数见表5。

表4 岩石实际物理力学参数与岩石模拟物理力学参数对照表

Table 4 Comparison of actual physical mechanics parameters of rock with physical and mechanical parameters of rock simulation

对比参数	试验结果	模拟值	误差/%
弹性模量/GPa	54.58	54.29	0.53
泊松比	0.27	0.26	3.70
抗压强度/MPa	127.96	125.93	1.60
黏聚力/MPa	22.06	24.83	12.50
摩擦角/（°）	47.8	49.00	2.50

新窗口打开| 下载CSV

表5 模拟采用细观参数

Table 5 Microscopic parameters used in the sumulation

参数类型	参数名	数值
密度/（kg·m^-3）		2 700
线性参数	弹性模量/GPa	48.0
	刚度比	1.0
	摩擦系数	1.0
平行粘结参数	弹性模量/GPa	9.0
	刚度比	2.5
	抗拉强度/MPa	75.0
	黏聚力/MPa	25.0
	摩擦角/（°）	75.0
	半径系数	1.5

新窗口打开| 下载CSV

2.2 节理参数反演

节理直剪试验所用试样在钻孔提取时已将填充物冲走，且提取的岩样分布较为集中，代表性不足，加之此次模拟节理间距扩大10倍，因此，模拟直剪试验的方法标定节理参数并不适用。在选用平滑节理模型的条件下，根据矿山实际崩落情况反演确定节理参数。

当拉底跨度约为90 m时，矿山发生大规模崩落，即在模型拉底第三步时发生了大面积崩落，根据矿山溜井拍摄的监控视频估计崩落高度约为30 m。假定节理刚度按照一般取值，黏聚力为抗拉强度的10倍，摩擦角为25°，选择较大黏聚力Coh₁和较小黏聚力Coh₂，采用不断二分法逼近，在满足第一、二步不发生大规模崩落的情况下，寻找使第三步发生大面积崩落且崩落高度约为30 m的黏聚力值，反演流程见图7。

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7 节理参数标定流程图

Fig.7 Flow chart of joint parameters calibration

图7流程运行前，需先反演模型地应力。由于PFC2D仅能从两侧边界墙Wall或柔性墙Ball进行挤压，水平应力无法在含有复杂地形的模型内部均匀传递。因此，仅在重力作用下平衡模型。考虑模拟动态崩落，将局部阻尼设置为0；由于崩落下的矿岩碰撞造成能量损失，法向临界阻尼设置为0.3。最终确定的节理力学参数见表6。

表6 节理参数

Table 6 Parameters of joint

参数名称	数值	参数名称	数值
法向刚度/GPa	6.0	摩擦角/（°）	25
切向刚度/GPa	3.0	抗拉强度/MPa	0.7
黏聚力/MPa	7.0

新窗口打开| 下载CSV

3 结果计算及分析

模型自拉底中段向两侧扩展，每步拉底30 m，模拟过程中不出矿，直至崩落区崩透地表时模拟结束。共模拟11步，处理模拟结果，并于下文详细分析崩落破坏机理与崩落区演化特征。

3.1 崩落破坏机理分析

追踪第一步拉底时裂纹扩展情况，结果如图8所示（图中黑线为裂纹）。裂纹首先出现在推进面与顶板交界处[图8（a）、8（b）]，此处为应力集中区；之后裂纹向上发展并相交，呈拱形[图8（c）、8（d）]；随后裂纹在顶板上方拱形区间内密集产生[图8（e）]；最终基于岩石破碎程度决定是否崩落[图8（f）]。

图8

新窗口打开| 下载原图ZIP| 生成PPT

图8 第一步拉底裂纹扩展

Fig.8 Crack propagation in the first step

裂纹在前两步拉底过程中发育、延伸，相互连通[图9（a）]；拉底第三步，紧贴拉底层的裂纹贯通，破碎严重，发生大规模崩落[图9（b）]，裂纹扩展成拱形，并在未崩落区间向拱形两侧的外围延伸。

图9

新窗口打开| 下载原图ZIP| 生成PPT

图9 裂纹扩展与破裂类型统计

Fig.9 Crack propagation and statistics of fracture types

统计拉底过程中裂纹产生数量、类型以及第八步拉底时裂纹空间分布情况，见图9（c）~9（d），原岩崩落主要发生节理拉伸破坏，拱形面底端原岩主要发生节理剪切破坏，崩落下的矿岩在冲击作用下主要发生岩石剪切破坏。

3.2 崩落区演化特征

随着拉底的推进，分析崩落区演化形态、顶板崩落高度和崩落质量情况，其中崩落区演化形态见图10，顶板高度和崩落质量统计结果如图11所示。

图10

新窗口打开| 下载原图ZIP| 生成PPT

图10 崩落区形态演化特征

Fig.10 Evolution characteristics of caving area shape

图11

新窗口打开| 下载原图ZIP| 生成PPT

图11 崩落区统计图

Fig.11 Statistics of caving area

由图10可知，前期拉底并未立刻引起顶板大规模崩落，崩落始于第三步拉底。随拉底向两侧推进，采场跨度增加，崩落区呈拱形向上、向四周发展，在第十一步时崩透地表，引起地表塌陷和邻近高陡边坡向谷地方向滑动。

由图11可知，因第一步与第二步未发生大规模崩落，所以崩落质量与顶板高度几乎为0；拉底前期崩落质量与顶板高度随拉底推进均匀增加，后期出现若干跳动，原因是后期某些拉底步骤仅在新开挖拉底区上方局部破裂，随着应力释放，其他位置岩石再次趋于稳定，所以质量增加较少且顶板高度几乎不变，但裂纹在此期间充分发育又促使下一步拉底时产生大规模崩落，造成崩落质量与顶板崩落高度突然增大。由于最后一步崩透地表，所以崩落质量与顶板高度相较前一步显著增加。

3 结论

采用PFC2D内置DFN功能，结合某矿山节理分布数据，生成与实际矿山相符的节理裂缝网络；通过模拟室内试验，反演标定完整岩石细观力学参数，并依据矿山实际崩落情况反演分析模拟所用节理细观强度参数。分析模拟崩落结果，所得结论如下：

（1）裂纹首先发育于拉底推进面与顶板相交位置，之后在拱形范围内延伸、发展、贯通，造成原岩破碎并崩落；崩落破坏类型主要为节理拉伸破坏，且崩落下的矿岩多数会发生剪切破坏。

（2）崩落区前期呈拱形向上发展，崩透地表后引起邻近陡坡向谷地滑动；某些拉底步骤，崩落区仅发生在新开挖拉底层上方局部范围，与崩落质量及顶板高度的变化趋势一致。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

冯兴隆

自然崩落法矿岩工程质量数字化评价及模拟技术研究

[D].长沙：中南大学，2010.