张友志 吴爱祥 王洪江 李 涛 孙 伟 闫其盼

(北京科技大学土木与环境工程学院，北京100083)

采矿塌陷区易诱发泥石流等地质灾害并造成生态破坏，严重威胁地表建筑稳定性和人员安全［1］。据不完全统计，全国因采矿引起的塌陷区有180多处，塌陷坑1 600多个，塌陷面积达1 150 km2之多。矿产资源开发带来的塌陷及其环境地质灾害，已经给矿区及矿业城市的可持续发展带来严重的影响［2］。尾砂作为一种廉价的矿山废料，大量用于矿山充填，既可以解决矿山充填材料的来源问题，又解决了矿山工业废料造成环境污染和占用耕地的问题［3］。目前处理矿山地表塌陷区的通常做法为采用矿山废弃物(煤矸石、粉煤灰、废石等)直接对塌陷区进行回填处置，但易造成回填体大规模塌陷、井下泥石流、地下水污染等次生灾害，引起矿石损失贫化。针对以上塌陷区处置方式的弊端，提出将矿山废弃物(尾砂、废石等)制备成混合处置体，采用柱塞泵通过管道泵送至塌陷区进行回填，以期达到有效控制活动塌陷区、实现井下安全、高效开采的目的。

颗粒流法是一种新兴的数值分析方法，在离散元理论基础上由 P.A.Cundall和 O.Strack［4］提出，广泛应用于岩石类材料基本特性、颗粒物质动力响应、岩石类介质破裂及其发展等基础性问题的研究［5］，适用于矿区采空区中的岩体断裂、坍塌、破碎和流动问题。该方法从微观角度出发，通过采用离散单元法来模拟圆形颗粒介质的运动及其相互作用，由平面内的平动和转动运动方程来确定每一时刻颗粒的位置和速度，能够自动模拟介质基本特性随应力环境的变化，实现岩土体对历史应力－应变记忆特性的模拟(屈服面变化、Kaiser效应等)，自动反映介质的连续非现行应力－应变关系、屈服强度和此后的应变软化或硬化过程等。目前，颗粒流法已广泛应用于岩石力学行为研究［6-10］，但在地表塌陷区回填处置方面的应用并不多见。

本研究以某矿地表塌陷区处置为工程背景，并基于颗粒流理论和PFC2D程序，采用Fish语言编制命令流，对塌陷区拟采用的回填工艺进行了数值模拟，得出不同回填工艺下混合处置体在塌陷坑底部采矿活动过程中的运移变化规律，研究结果可作为地表塌陷区回填处置的参考依据。

1 工程概况

某矿细脉带矿体自开采以来，井下贯通地表形成的塌陷坑已有10余个，塌陷坑相互嵌套，形状各异。如果不加处理，地表塌陷区仍有不断扩大的趋势，最大回填量将达到236.9万m3。同时该矿年产矿石量达220万t，尾矿产率96%。大量的尾矿以10%的浓度被排往尾矿库，使尾矿库的安全与环保的压力增大，并消耗大量的生产用水。另外，该矿地表碎石坑尚有35万m3的碎石，还有－3 mm的重介质尾矿砂，均可以用于回填塌陷区。为了充分利用废石、尾矿等工业固体废弃物，拟通过向尾矿中添加碎石的方式制备胶结与非胶结2种回填处置体，对地表塌陷区进行回填处理，达到消除塌陷区的目的，实现地表安全保障与环境保护的双重目标。

2 PFC2D数值模拟

在PFC2D的数值模拟中，颗粒和黏结体的相关参数与通常意义上的宏观参数存在较大的区别，岩石试件的力学参数不能直接输入到软件中，只有通过模拟力学实验，进行大量的数值模拟试验，根据数值试验的结果才能得到颗粒或者胶结材料的微观参数。当数值模拟的试件力学参数与实验室的宏观参数基本相符时，即可得到颗粒的微观参数［11-14］。只有将微观参数输入程序进行计算，才能使数值模拟与实际情况相符。通过室内实验对数值模型中表土层、稳定层、扰动带及胶结和非胶结处置体的宏观参数进行了测定，各类岩层的宏观、微观参数见表1。

表1 岩层宏观和微观参数对照Table 1 Contrast of rock macro and microscopic parameters

以T110塌陷区为模拟对象进行建模:模型长70 m，高84 m，塌陷坑深度为40 m。针对T110塌陷区的回填情况，拟采取完全胶结回填、完全非胶结回填、10 m厚度分层回填、20 m厚度分层回填及5 m、15 m厚度交替回填这5个方案对回填情况进行模拟。原始塌陷坑如图1所示，以完全胶结回填为例，回填效果模型如图2所示。

图1 原始塌陷坑Fig.1 Original collapse pit

图2 完全胶结回填效果模型Fig.2 Effect model of fully cemented backfill

不同回填方式的初始模型建成后，开始对塌陷坑底部的采矿活动进行模拟。按照该矿的年采矿量对塌陷坑底部矿石进行采出模拟，同时在回填体表层设置一系列监测点对回填体的沉陷位移进行监测，在塌陷坑底部设置测量圆对回填体与塌陷坑的接触状况及应力分布状况进行测量。监测点及测量圆布置情况见图3。

图3 监测点及测量圆布置Fig.3 Layout schematic of monitoring point and measuring round

3 数值模拟结果分析

以完全胶结回填、完全非胶结回填和5 m、15 m厚度交替回填为例，塌陷坑底部采矿活动结束后，不同回填方式下的处置体及塌陷坑运移情况分别如图4～图6所示。

图4 完全胶结回填方案运移情况Fig.4 Migration situation of fully cemented backfill scheme

图5 完全非胶结回填方案运移情况Fig.5 Migration situation of completely non-cemented backfill scheme

图6 5 m、15 m交替回填方案运移情况Fig.6 Migration situation of 5 m and 15 m alternative backfill scheme

3.1 监测点沉陷位移测量

塌陷区内处置体在回填过程中会出现不同程度下陷的情况，通过对处置体表层颗粒在重力方向的位移变化情况进行监测，可在很大程度上对处置体整体下陷情况做出评判，从而对回填方案的优劣情况做出评判。各方案在图3中5个监测点的沉陷位移曲线见图7。

从图7中可以得出:

(1)塌陷区内回填处置体在底部开挖过程中出现不均匀下降的现象，其中处置体与塌陷区边帮接触处的沉陷位移量较大，这是由于回填处置体和塌陷区边帮接触处不易形成稳定的胶结层造成的。另外处置体中心有较大沉陷位移量，这与工程实际也是相吻合的。

(2)完全胶结回填方案中处置体沉陷位移量最小，中心沉陷位移量为1.1 m，完全非胶结回填方案的沉陷位移量最大，中心沉陷位移量为2 m。分层回填方案的中心沉陷位移量在1.3～1.35 m，三者相差不大。从经济角度出发，优先推荐5、15 m分层交替回填方式。

3.2 回填体在底部开挖后的应力分布

图8为回填模型在塌陷坑底部矿石开挖后的接触力分布情况。颗粒间接触力的大小与线条粗细成正比例关系。从图8中可以看出，模型中稳定层两侧的颗粒间接触力均以竖向压力作用为主，说明在模拟计算过程中，扰动带内颗粒间存在一定的空隙空间，颗粒间的接触力较小，而稳定层内颗粒间接触力则相对较大。此外，在图8中的浅色标记区域内，颗粒间的接触力以横向作用力为主，形成一个横向应力拱，并与稳定层颗粒共同构成了一个连续的力学拱结构，这对阻止塌陷区持续下陷、维持塌陷区的稳定性起到了重要的稳固作用。因而在接触力分布图中，塌陷区下方横向力越明显，说明颗粒间横向压力作用越大，形成的横向应力拱也越稳固，对维持塌陷区的稳定性也就越有利。

图8 回填模型在底部矿石开挖后的接触力分布情况Fig.8 Distribution of contact force after ore being excavated at bottom

3.3 回填体与塌陷坑离层分析

随着塌陷坑内处置体的不断下移，处置体与塌陷坑底部岩层不一定能够良好接触，从而在接触处可能会产生一定的离层空间。接触处离层现象形成的空洞会对处置体的均匀下移起到阻碍作用，从而致使塌陷坑内处置体不能连续均衡稳定下移，对塌陷坑的稳定性造成影响。离层空间的大小是通过离层因子来判定的，在PFC2D中离层因子可通过对测量圆内的孔隙率进行测量获得，从而对不同方案的离层程度进行分析与评判:离层因子值较小，则说明测量圆内部颗粒分布比较均匀，即处置体颗粒与塌陷坑底部接触较为良好;而离层因子值较大则说明塌陷坑底部与处置体接触不佳。5个回填方案离层因子的测量数据见表2。

表2 各回填方案离层因子Table 2 Separation factor of each backfill scheme

从表2中可以得出:

(1)完全胶结回填方案的离层因子最大，离层因子值为0.18，完全非胶结回填方案的离层因子最小，值为0.12。

(2)方案3、4、5 的离层因子分别为0.15、0.16 和0.148，三者相差不大。

4 结论

(1)塌陷坑底部矿石的开采会对塌陷处回填体造成影响，使之产生位移，完全胶结回填方案中回填处置体受到的影响最小，沉陷位移也最小，分层回填方案效果次之，完全非胶结回填方案位移最大，控制处置体位移的能力也最差。

(2)胶结作用下，回填处置体与周围岩体形成新的应力分布体系，所形成的横向连续应力拱结构可以阻止塌陷区下陷，从而对地表塌陷区稳定性起到较好的控制作用。

(3)回填处置体与塌陷坑底部岩层之间形成的离层空间对处置体的整体稳定均衡下移起到阻碍作用，胶结回填层厚度不小于5 m、非胶结回填层厚度5～15 m时，塌陷坑的稳定性最佳。

［1］ Olivier D，Harlalka A.Numerical study of the soil-structure interaction within mining subsidence areas［J］.Computers and Geotechnics，2010，37:802-816.

［2］ Hu R L，Yue Z Q，Wang L C，et al.Review on current status and challenging issues of land subsidence in China［J］.Engineering Geology，2004，76:65-77.

［3］ 刘同有.充填采矿技术与应用［M］.北京:冶金工业出版社，2001.Liu Tongyou.Filling Mining Technology and Application［M］.Beijing:Metallurgical Industry Press，2001.

［4］ Cundall P A，Strack O.Particle Flow Code in 2D［M］.Minnesota:Itasca Consulting Group，Inc.，1999.

［5］ 吴顺川，周 喻，高 斌.卸载岩爆试验及PFC3D数值模拟研究［J］.岩石力学与工程学报，2010，29(S2):4082-4088.Wu Shunchuan，Zhou Yu，Gao Bin.Study of unloading tests of rock burst and PFC3Dnumerical simulation［J］.Chinese Journal of Rock Mechanics and Engineering，2010，29(S2):4082-4088.

［6］ 王培涛，杨天鸿，柳小波.边孔角对无底柱分段崩落法放矿影响的颗粒流数值模拟研究［J］.金属矿山，2010(3):12-16.Wang Peitao，Yang Tianhong，Liu Xiaobo.Particle flow numeral simulation investigation on influence of lateral opening angle on ore drawing with sublevel pillarless caving［J］.Metal Mine，2010(3):12-16.

［7］ 刘增辉，高 谦，李 欣，等.回填废石层下放矿损贫预测及颗粒流数值分析［J］.金属矿山，2011(2):5-8.Liu Zenghui，Gao Qian，Li Xin，et al.Ore dilution prediction of ore drawing below backfilling of waste rock layer and particle flow numerical analysis［J］.Metal Mine，2011(2):5-8.

［8］ Potyondy D O，Cundall P A.A bonded-particle model for rock［J］.International Journal of Rock Mechanics ＆ Mining Sciences，2004，41(8):1329-1364.

［9］ Potyondy D O，Cundall P A.Simulating stress corrosion with a bonded-particle model for rock［J］.International Journal of Rock Mechanics ＆ Mining Sciences，2007，44(5):677-691.

［10］ Wang 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 ＆ Mining Sciences，2003，40:415-424.

［11］ Jeoungseok Y.Application of experimental design and optimization to PFC model calibration in uniaxial compression simulation［J］.International Journal of Rock Mechanics ＆ Mining Sciences，2007，44:871-889.

［12］ Okan S，Nuri A A.Numerical simulation of rock cutting using the discrete element method［J］.International Journal of Rock Mechanics ＆ Mining Sciences，2011，48:434-442.

［13］ Tawadrous A S，DeGagn D.Prediction of uniaxial compression PFC3Dmodel micro-properties using artificial neural networks［J］.International Journal for Numerical and Analytical Methods in Geomechanics，2009，33:1953-1962.

［14］ Kazerani T，Zhao J.Micromechanical parameters in bonded particle method for modeling of brittle material failure［J］.International Journal for Numerical and Analytical Methods in Geomechanics，2010，34:1877-1895.