隔离层下单漏斗放矿力链演变的参数敏感性分析

在协同开采理念的指导下，笔者于2010年提出一种新的采矿方法[1]，即大量放矿同步充填无顶柱留矿采矿法(简称同步充填留矿法)

该方法在放矿前预先于矿石堆上铺设柔性隔离层，然后从回风巷道充填废石，最后利用隔离层实现充填料与矿石同步下沉，进而从宏观和细观角度研究同步充填留矿法放矿过程散体介质体系的流动规律和接触力特性

研究发现，颗粒物质由大量离散的固体颗粒相互接触作用并形成复杂的体系[2]，在形态上包含固相、液相和气相三种

其中，固相和液相属于密集颗粒体系，对外界呈现出敏感性、非线性响应以及自组织行为等复杂的力学特性

采场内的矿石同样是由大量形态各异的固体颗粒组成，颗粒间密集分布，并且颗粒间的接触力通常大于间隙气体或液体，在传递荷载中起主要作用

颗粒相互接触、挤压，接触力沿着颗粒进行传递，从而形成准直线分布的链状路径(即力链)[3-5]

不同的力链相互交织，进而构成一个完整的力链网络并非均匀地贯穿于颗粒体系

矿石颗粒在流动时，力链不断进行着断裂重组的过程，其对外界条件的变化极其敏感，散体介质体系内任何一个参数发生变化，都可能影响力链网络的完整性和稳定性

针对力链复杂多变的演化特性与作用机理，国内外学者开展了一系列的研究，并取得了丰硕的研究成果

在力链演变规律方面，周勇等[6]基于散体介质双轴加载光弹物理试验装置和相似材料模拟试验，研究综放开采过程中覆岩断裂的演化特征与非连续覆岩关键层力链结构变化规律的关联性，并提取了力链作为采场覆岩矿压变化的判据

孙其诚等[7]以单轴压缩数值模拟试验为例，深入分析了土体颗粒骨架结构与力链形态之间的关系，指出强力链是决定颗粒体系宏观力学行为的关键因素

徐正红[8]构建了不同粒径的力链压曲变形分析模型，由模型提出物理和理论意义明确的塑性演化准则，并对力链模型进行了优化及参数影响效应分析

FANG等[9]通过光弹试验直观地观察了模型中力链的分布状态，讨论了应力拱的成拱机理，还在试验装置下方安装两根钢柱引起局部变形，并指出模型的几何参数和应力状态都会影响应力拱的形成

在离散元法应用方面，YANG等[10]研究了破碎颗粒材料在单向压缩过程中接触力和配位数的分布规律，发现颗粒内部接触力的分形分布是导致颗粒分形尺寸分布逐渐变化的主要因素

张炜等[11]利用离散元法模拟了金属粉末的高速压制过程，并对比了不同压制状态下的力链的变化规律，发现高速压制与普通压制过程力链的初始状态存在一定差异，且进行高速压制时力链特征参数的变化更为频繁

IMSEEH等[12]基于离散元法研究了三维砂土形态对砂土内部力链产生和演化的影响，对颗粒间接触力和力链的发育形态进行了准确的数值预测

PATINO-RAMIREZ等[13]将力链网络分解为单条力链，采用离散元法模拟了在颗粒状试样上移动的肋状界面，结果表明材料密度越大，力链集中度相对越大，对应的界面剪切强度就越大

此外，在正交试验分析方法应用方面，党洁等[14]选取隔离矿柱矿房侧壁的回采高度、侧壁厚度、采场长度、采场宽度等6种影响因素，建立了L8(27)正交试验设计表，构建了隔离矿柱的回采模型，模拟分析了不同试验方案中隔离矿柱回采过程的应力分布及位移变化情况，结果表明在6种影响因素中，矿房侧壁的回采高度对隔离矿柱的稳定性影响最为显著

黄德镛等[15]等根据岩体力学参数、充填体力学参数以及采场的布置特征建立了采场结构模型，利用正交设计方法对16种采场方案进行稳定性分析，并对不同方案进行模型优选，得出了最优采场结构参数尺寸和最佳充填配比

上述研究对力链的演化特性和作用机理进行了初步探索，同时也涉及了正交试验分析方法在矿业领域的应用，但关于影响力链动态变化因素及其敏感性的研究仍较少

因此，本文基于前人的研究成果，拟采用正交试验法和矩阵分析法，以大量放矿同步充填无顶柱留矿采矿法为背景，通过建立数值模拟放矿模型，研究在不同参数组合的影响下散体介质力链演变的规律及其敏感性，以期为力链研究提供一种新的思路并促进颗粒物质力学的发展

1放矿模型构建本文基于颗粒离散元法，利用PFC2D软件构建数值模拟放矿模型

在模拟放矿时，颗粒会产生位移与旋转，进而在颗粒间产生新的接触

因此，为准确反映颗粒间的物理力学性质，需选择合理的接触模型并确定相应的细观力学参数[16]

由于矿石颗粒体系是由大量非均匀块体所组成，其形态难以使用一个固定状态进行描述，利用PFC的Clump命令可以将数个粒径较小的颗粒组成不同的岩块，但在模拟过程中将所有矿块进行精确表征的难度较大，同时也会影响到后续的数据分析结果

在PFC内置的接触模型中，与其他模型相比，抗滚动线性接触模型增加了抗滚动因数[17]，通过调整因数的大小来描述颗粒形态对散体介质流动特性的影响

当颗粒相互接触时，该模型使接触位置增加了颗粒相对转动时呈线性增长的内力矩，这将降低颗粒的转动能力，可用来模拟非均匀的岩块体，与采场中矿岩的接触状态存在相似之处

因此，本文选取抗滚动线性接触模型来构建同步充填留矿法的数值模拟模型

由于选取随机形态的矿石颗粒进行模拟时会频繁出现漏斗堵塞的情况，而抗滚动线性接触模型可以抵消颗粒形态对矿石颗粒流动状态的影响，因此，本次试验采用固定粒径生成矿石颗粒

矿石颗粒的力学参数主要包括法向和切向刚度、摩擦因数以及抗滚动摩擦因数等

其中，根据Hertz接触理论和Mindlin-Deresiewicz接触理论分别确定法向刚度和切向刚度[18]；颗粒摩擦因数和抗滚动因数逐渐增加会使安息角表现出先增大后趋于稳定的趋势，且两种参数的不同组合均可获得相同的自然安息角，因此可通过矿石颗粒的自然安息角间接标定颗粒摩擦因数和抗滚动因数[19]；结合先期开展的散体介质流动规律物理试验[20]，经过多次调试并综合对比试验结果进行取值，矿石颗粒及墙体的初始力学参数设置如表1所示

表1墙体及初始矿石颗粒力学参数Table 1Mechanical parameters of wall and initial ore particlesWallInitial ore particleNormalstiffness/(N·m-1)Tangentialstiffness/(N·m-1)FrictioncoefficientNormalstiffness/(N·m-1)Tangentialstiffness/(N·m-1)FrictioncoefficientOre particledensity/(kg·m-3)Ore particleradius/m1×1071×1070.55×1075×1070.328000.008根据以上确定的接触模型、矿石颗粒力学参数、接触模型参数等，基于几何相似准则，结合散体介质流物理试验的模型尺寸[21]，将数值模拟模型的尺寸设为168 cm×128 cm，每个放矿口的间距为24 cm

利用PFC软件构建同步充填留矿法数值模拟模型，如图1所示

图1同步充填留矿法模型



Fig. 1Model of synchronous filling shrinkage method颗粒采用自重堆积法生成，数值模拟放矿试验的步骤如下：1) 颗粒间的初始接触模型设为线性接触模型，摩擦因数设为0.3；为便于观察放矿过程中矿石颗粒的流动规律，待初始模型平衡后，以10 cm为间隔将颗粒进行分层，并赋予不同的颜色

2) 利用Cubic命令在矿石颗粒上方生成一排细小颗粒，采用平行黏结模型使颗粒相互黏结进而生成柔性隔离层

隔离层的尺寸为长250 cm，颗粒粒径为0.15 cm，其初始细观力学参数设置如表2所示

表2隔离层初始力学参数Table 2Initial mechanical parameters of the isolation layerNormalStiffness/(N·m-1)TangentialStiffness/(N·m-1)Normal stiffness ofparallel bond/(N·m-1)Shear stiffness ofparallel bond/(N·m-1)FrictioncoefficientOre particledensity/(kg·m-3)Elasticmodulus ofparallel bond/PaRadius/m1×1071×1071×1061×1060.420005×1070.00153) 为抵消颗粒形态对矿石颗粒流动状态的影响，在放矿开始前，将颗粒接触模型由线性接触模型转变为抗滚动线性接触模型

4) 打开4号放矿口，矿石颗粒随即流出，每计算一定时步，关闭放矿口，并在矿石颗粒上方充填适量的废石颗粒，借此实现同步充填效果，待模型在自重作用下解算平衡后，删除多余的废石颗粒，并再次打开放矿口，进入下一循环，直至矿石颗粒全部放出

2正交试验设计正交试验设计是分式析因设计的主要方法[22]，具备快速、高效的优点

在进行设计时可优化试验的工况条件，从而找出最优的试验方案

在正交试验中，一个因素的改变通常会引起其他因素的变动，相较于单因素分析，多因素敏感性分析能够反映不同因素发生变化时对整体试验结果产生的综合影响，在一定程度上能弥补单因素分析的局限性[23]

同步充填留矿法与传统采矿方法的区别在于放矿前预先铺设柔性隔离层在矿石堆上，散体介质在流动过程中，因隔离层存在牵扯与控制作用，使散体介质受到充填料的非自由表面纵向荷载，隔离层的次生横向荷载以及采场边界限制条件等复合作用

因此，在考虑隔离层下散体介质流动过程影响力链演变规律的因素时，应将隔离层的力学性能参数作为主要影响因素之一

结合数值模拟试验的特点，将矿石颗粒粒径和摩擦因数，隔离层厚度和隔离层界面摩擦因数作为主要影响因素

综上所述，根据正交试验方法设计原理和同步充填留矿法的工艺特点，本文将隔离层厚度A、隔离层界面摩擦因数B、矿石颗粒摩擦因数C及颗粒半径D作为正交试验的四种影响因素

针对四种主要影响因素，选取四因素三水平正交试验设计表，在构建放矿模型时，已对矿石颗粒和隔离层力学参数的选取进行说明，结合矿石流动规律影响因素的敏感性相关研究[24]，各因素及对应水平表如表3所示，L9(34)正交试验设计表如表4所示

表3正交试验因素表Table 3Table of orthogonal test factorsLevelThickness,A/mInterfacial frictioncoefficient, BParticle frictioncoefficient, CParticle radius,D/m10.0030.20.20.00620.0040.50.50.00730.0050.80.80.008表4正交试验水平设计表Table 4Design table for orthogonal test levelTest No.Thickness,A/mInterfacial frictioncoefficient, BParticle friction coefficient, CParticle radius,D/m10.0030.20.20.00620.0030.50.50.00730.0030.80.80.00840.0040.20.50.00850.0040.50.80.00660.0040.80.20.00770.0050.20.80.00780.0050.50.20.00890.0050.80.50.0063力链演变的一般规律分析3.1力链宏观分布规律分析由于参数设置的水平不同，为避免在漏斗口出现堵塞的情况，编译防止漏斗堵塞的Fish函数

每个放矿节点利用History命令记录并输出颗粒接触信息，通过Geometry命令使放矿模型散体介质体系的力链宏观分布可视化

为描述力链宏观的一般变化规律，以正交试验4为例，选取第1次、第5次、第10次以及第15次放矿的力链宏观分布图进行分析

图2所示为正交试验4的力链宏观分布图(图中蓝色线条表示力链)

图2力链的宏观分布图



Fig. 2Macroscopic distribution of force chains: (a) First draw; (b) Fifth draw; (c) Tenth draw; (d) Fifteenth draw由图2可知，在放矿初期，力链分布较为均匀，主要集中在放矿模型的下部，模型上部的力链较为松散；随着放矿进行，隔离层在充填料荷载的作用下逐渐下沉，模型中部的力链逐渐向两侧移动，并以放出口中心线为对称轴沿着模型边壁向上延伸，不同力链之间的空隙逐渐缩小，部分力链出现较明显的集中效应；直至放矿后期，随着充填颗粒的数量逐渐增加，模型中未放出的矿石颗粒与隔离层相互作用，使隔离层上部的力链数量逐渐增加，且力链垂直于隔离层界面并向上蔓延，呈现较明显的方向性，模型中力链分布状态也表现出一定的对称性

3.2力链数量和长度变化特征分析力链数量在一定程度上可以反映散体介质内力链网络的变化过程

根据放矿过程不同节点获取的颗粒信息，包括颗粒编号、粒径、坐标、接触力等，利用PFC的Fish语言将颗粒接触信息输出；经过预先整理后，基于力链的识别判据[25-26]，编写力链的识别程序，检索体系内满足力链成链条件的矿石颗粒，记录并输出力链特征参数，实现力链信息的自动检索和识别

利用力链识别程序对试验4放矿过程的力链数目进行统计，结果如图3所示

图3力链数目的变化规律



Fig. 3Changing law of force chain number由图3可知，在放矿初期，散体介质原有的平衡被打破，导致体系内的力链数量在波动中逐渐减少；至放矿中期，随着同步充填工序进行，力链数量仍然保持着继续减少的趋势；在第13次放矿结束后，力链数量有所增加，最终稳定在761条左右

从整体来看，当放矿次数逐渐增加，体系内矿石颗粒因受到来自隔离层、自身重力及放矿口等因素的综合作用，初始游离的矿石颗粒数量逐渐减少，体系荷载逐渐被力链承担，体系内的短力链数量逐渐减少，长力链数量逐渐增加，因此力链数量总体呈现出下降的趋势

当模型内的大部分颗粒被放出，隔离层逐渐下沉至模型底部，其两侧未放出的矿石颗粒受到的挤压作用有所增强，颗粒间接触的密集程度强于放矿中期，这使得体系内的荷载逐渐趋于稳定，力链数量有所增加

通过归一化处理，计算每次放矿后散体介质内的力链数目，并统计力链长度的概率分布情况，结果如图4所示(因长度大于15的力链占比较小，故本文中暂不考虑长度大于15的力链)

图4力链长度的概率分布



Fig. 4Probability distribution of force chain length由图4可知，力链网络主要由短力链和长力链组成，其中短力链为主要组成部分

随着力链长度增加，其比例逐渐下降，力链长度的概率分布曲面呈指数形式递减

同时，可以看出长度为3的短力链数量比例随着放矿次数增加而逐渐减少，这表明在放矿初期，体系内的长力链分解成短力链，而后颗粒在流动过程中不断进行着断裂重组的过程，新组成的力链相互交织，从而生成长力链来维持体系的稳定

因此，随着放矿进行，短力链数量占比随之下降，而长力链数量占比逐渐提高

3.3力链强度分析力链强度是单条力链所有接触力强度的平均值，能够反映该条力链的承载能力

对于力链网络整体而言，需要统计颗粒体系中所有力链的强度，并计算平均值以表示力链网络的强度

将不同放矿次数下颗粒体系内部所有力链强度的平均值进行统计，对力链网络强度的变化规律进行初步分析

图5所示为力链网络的强度

图5力链网络的强度



Fig. 5Strength of force chain network由图5可知，根据力链强度的曲线变化规律，将力链网络强度的变化趋势分为三个阶段

其中，第一阶段为第1次放矿至第5次放矿结束，力链网络强度有所下降，但幅度较小；第二阶段为第5次放矿至第10次放矿结束，力链网络强度先下降后提升，且下降的幅度高于第一阶段；第三阶段为第10次放矿至第15次放矿结束，在该阶段内力链网络强度有所增加

为便于描述不同放矿次数的力链强度分布规律，将力链强度数据进行归一化处理，并计算体系内力链强度的概率分布，结果如图6所示

图6力链强度的概率分布



Fig. 6Probability distribution of force chain strength由图6可知，在放矿过程中，力链强度的分布概率呈现出相似的变化规律，先呈指数式上升，在f=0.7处达到峰值，然后呈指数式递减

同时，可以看出不同放矿次数力链强度的分布概率变化较小，力链强度主要集中分布于f=0.7处

3.4力链准直系数分析力链准直系数表示力链保持直线链状，传递外部荷载的能力强弱，能够反映力链网络的稳定性大小，对于研究力链准直性具有重要意义

统计单漏斗放矿条件下散体介质体系内所有力链的准直系数，并进行累加求和，根据力链数目计算平均值，绘制不同放矿节点下的力链准直系数分布曲线，如图7所示

图7力链准直系数的变化规律



Fig. 7Changing law of collimation coefficient of force chain由图7可知，根据力链准直系数的变化规律，将力链准直系数的变化趋势分为三个阶段

其中，第一阶段为第1次放矿至第5次放矿，该阶段力链准直系数在波动中逐渐下降；第二阶段为第6次放矿至第10次放矿结束，此阶段力链准直系数逐渐增加，但幅度较小；第三阶段为第11次放矿至第15次放矿，此阶段力链准直系数逐渐增加，且增加的幅度较大

从整个放矿过程来看，放矿结束后的力链准直系数高于放矿初期，如准直系数越大，则力链可以表现出更强的准直性，促使更多的矿石颗粒组成长力链

结合力链数量、长度、强度等参数在放矿过程中的变化特征可知，随着放矿次数增加，力链准直系数在波动中逐渐增大，使更多游离的强接触颗粒组成长力链的概率增大

因此，随着放矿进行，散体介质中短力链数量所占比例逐渐减少，长力链的比例逐渐增多，力链数目整体呈下降的趋势

4力链特征参数敏感性分析4.1评价指标选取通过对单漏斗放矿条件下力链演变的一般规律进行深入阐述，可知力链网络的变化受到不同参数的影响

因此，为进一步了解各因素对力链演变的作用规律，需建立评价指标体系进行量化分析

结合正交试验设计的特点以及力链参数特性，本文选取力链网络强度和力链准直系数作为评价指标

4.1.1力链网络强度力链强度为该条力链中所有接触力的平均值，可表示力链承受外部荷载能力的强弱

对于一个完整的力链网络，则需要统计所有力链强度的平均值表示力链网络强度，其值越大，则整个力链网络的稳定性越好

4.1.2力链准直系数力链具有较强的准直性，可保证在承受外部荷载时不发生明显的屈曲，并且能够较好地维持直线链状结构，使更多矿石颗粒相互接触，组成力链的一部分，从而形成长力链

同时，如力链具有较强的准直性，可使外部荷载在体系中顺利传递，保证整个力链网络维持稳定

4.2指标变化规律分析将不同试验在放矿过程中的力链强度进行统计，获得正交试验力链强度的变化规律如图8所示

图8正交试验力链强度的变化规律



Fig. 8Change law of force chain strength in orthogonal test由图8可知，在不同放矿条件下，不同试验的力链网络强度表现出相似的变化规律

在放矿初期，力链强度的波动幅度较小，放矿中期力链强度略微下降，直至放矿后期力链强度在波动中逐渐上升

此外，还可从中看出试验3、4、8在放矿过程中力链强度相对较大，试验2、6、7次之，试验1、5、9的力链强度较小

究其原因，在试验3、4、8中，矿石颗粒半径为0.008 m，同时试验3的矿石颗粒摩擦因数为0.8；当颗粒间摩擦因数取值较大，矿石颗粒在流动过程中产生接触作用更加剧烈，从而使相应颗粒组成的力链强度较大

结合试验2、6、7和1、5、9的力链强度演化规律可知，矿石颗粒半径对力链强度的变化规律影响较大，其次矿石颗粒摩擦因数对力链强度的变化也会产生一定影响，而在本次试验中隔离层厚度和隔离层界面摩擦因数对力链强度的演化规律影响程度较小

统计不同正交试验散体介质体系内所有力链的准直系数，并进行累加求和，根据力链数目计算平均值，得到不同试验的力链准直系数，绘制各试验在不同放矿节点下的力链准直系数分布情况，结果如图9所示

图9正交试验力链准直系数的变化规律



Fig. 9Change law of collimation coefficient in orthogonal test由图9可知，不同试验的力链准直系数在放矿前期波动幅度较小，在放矿中后期的波动幅度较大，整体呈缓慢上升的趋势

力链准直系数可以从侧面反映其他力链特征参数的变化规律，如准直系数越大，则力链可表现出较强的准直性，使更多的矿石颗粒相互连接形成长力链

各试验散体介质体系的力链准直系数保持在0.8~0.95之间，其中试验3、5、7的力链准直系数较大，试验2、4、9次之，试验1、6、8相对较小

通过与各试验的参数条件进行对比，发现试验3、5、7的矿石颗粒摩擦因数均是0.8，试验3、5、7的矿石颗粒摩擦因数是0.5，试验1、6、8的矿石颗粒摩擦因数则是0.2

由此可知，在单漏斗放矿条件下，散体介质体系的力链准直系数大小受矿石颗粒摩擦因数影响较大；颗粒摩擦因数越大，准直系数相对也越大

4.3矩阵优化分析通过计算可得到不同放矿条件下力链网络强度和力链准直系数的平均值，并与数值模拟试验的四种影响因素及对应的水平组成正交试验直观分析表，如表5所示

表5单漏斗放矿正交试验直观分析表Table 5Visual analysis table of single funnel drawing orthogonal testTest No.Thickness,A/mInterfacial friction coefficient, BParticle friction coefficient, CParticle radius,D/mForce chain strength,PCollimation coefficient,δ10.0030.20.20.006163.5900.83920.0030.50.50.007232.7050.88730.0030.80.80.008296.6720.89040.0040.20.50.008293.6830.88550.0040.50.80.006188.4850.89860.0040.80.20.007216.1070.84670.0050.20.80.007240.6500.89480.0050.50.20.008272.5220.84590.0050.80.50.006176.7830.887K1230.989232.641217.406176.286Visual analysis of PK2232.758231.237234.390229.821K3229.985229.854241.936287.626R2.7732.78724.530111.340Optimal schemeA2B1C3D3K10.8720.8730.8430.875Visual analysis of δK20.8760.8760.8860.876K30.8750.8750.8940.873R0.0040.0030.0510.003Optimal schemeA2B2C3D2由表5可知，对于力链网络强度而言，最优方案是A2B1C3D3；对于准直系数而言，最优方案是A2B2C3D2

针对多指标正交试验问题，利用矩阵分析法[27-28]，分别计算两个指标的权矩阵，进一步优选最佳的方案

矩阵分析模型是一个包含指标层、因素层以及水平层的数据结构模型，具体的分层内容如下

4.3.1试验评价指标层该层数据对应的正交试验有l个因素，每个因素有m个水平，单个因素Ai在第j个水平对应指标的平均值为kij

如试验评价指标Kij值越大，对试验效果越好，则令Kij=kij；反之则令Kij=1/kij

以此为基础，建立矩阵如式(1)所示

(1) 4.3.2因素层因素层矩阵，可构建矩阵如式(2)所示

(2) 4.4.3水平层正交试验中可包含多个因素，若因素Ai的极差为si，令，构建矩阵如式(3)所示

(3) 通过定义构建以上三层矩阵，进而建立影响试验指标的权矩阵；令，其中，，后者是因素A1对应第一水平指标值占该因素所有水平指标之和的比；，该项为因素A1极差占所有因素极差之和的比；两项比值的乘积能够表示因素对指标值大小的影响程度，同时也能反映该因素的极差大小

通过综合计算，可得到各因素水平对试验指标的影响权重大小，进而可根据权重大小优选出最佳方案，筛选出各因素对指标影响大小的主次顺序

在正交试验中，第一个试验指标为力链网络强度P，其值越大越好，矩阵分析过程如下

(4) (5) (6) (7) 同理，第二个试验指标为力链准直系数δ，其值越大越好，权矩阵计算结果如下：(8) 经过计算正交试验的两个指标值权矩阵，进一步得出总的权矩阵，总权矩阵为两个指标权矩阵的平均值，计算结果如下：(9) 经过综合计算，绘制正交试验各因素所占总权重的分布曲线，如图10所示

图10因素权重分布图



Fig. 10Weight distribution of different factors由图10及权重计算结果可知，各因素对正交试验指标影响的主次顺序是CDAB

因素A中影响权重最大的是A2，因素B是B2水平，因素C是C3水平，因素D是D3水平

因此，单漏斗正交试验的最优方案为A2B2C3D3

对应因素的水平即厚度0.004 m，隔离层界面摩擦因数0.5，颗粒摩擦因数0.8，颗粒半径0.008 m，此时散体介质体系力链网络在放矿过程中最为稳定

5结论1) 在放矿初期，力链网络分布均较为均匀，主要集中在模型下部，模型上部的力链较为松散；随着放矿进行，模型中部的力链逐渐向两侧移动，并以放出口中心线为对称轴沿着模型边壁向上延伸，部分力链出现较明显的集中效应；直至放矿后期，隔离层上部的力链数量逐渐增加，并垂直于隔离层界面向上蔓延，呈现较明显的方向性

2) 随着放矿次数增加，力链数量先逐渐减少后逐渐增加，最终稳定在761条左右

力链长度的概率分布曲面呈指数形式递减，长度为3的短力链数量所占比例逐渐减少，但仍占据着力链网络的主要部分，长力链数量所占比例逐渐增加，力链数量相比放矿初期有所减少

3) 力链强度在第一阶段有所下降，但幅度较小，在第二阶段先下降后提升，且下降的幅度高于第一阶段，在第三阶段有所增加

不同放矿节点力链强度分布概率呈现出相似的变化规律，先呈指数式上升，在f=0.7处达到峰值，然后呈指数式递减

力链强度的分布概率变化较小，力链强度主要集中分布于f=0.7处

4) 正交试验表明，四种影响因素对试验指标影响的主次顺序是颗粒摩擦因数C、颗粒半径D、隔离层厚度A、隔离层界面摩擦因数B

因素A中影响权重最大的是A2水平，因素B是B2水平，因素C是C3水平，因素D是D3水平

对应的最优方案为A2B2C3D3，即厚度0.004 m、隔离层界面摩擦因数0.5、颗粒摩擦因数0.8、颗粒半径0.008 m，此时散体介质力链网络在放矿过程中最稳定

参考文献

陈庆发, 吴仲雄. 大量放矿同步充填无顶柱留矿采矿方法: CN101864957A[P

2010-10-20.CHEN Q F, WU Z X. Great ore drawing simultaneous filling non-top column shrinkage mining method: CN101864957A[P

2010-10-20.王光谦, 孙其诚. 颗粒物质及其多尺度结构统计规律[J

工程力学, 2009, 26(S2): 1-7.WANG G Q, SUN Q C. Granular matter and the scaling laws[J

Engineering Mechanics, 2009, 26(S2): 1-7.ZHANG W, ZHOU J, ZHANG X J, et al. Quantitative investigation into the relation between force chains and stress transmission during high-velocity compaction of powder[J

Journal of the Korean Physical Society, 2019, 74(7): 660-673.ZHANG W, ZHOU J, ZHANG X J, et al. Quantitative investigation on force chain lengths during high velocity compaction of ferrous powder[J

Modern Physics Letters B, 2019, 33(10): 1950113.NARUMI T, SEE H, SUZUKI A, et al. Response of concentrated suspensions under large amplitude oscillatory shear flow[J

Journal of Rheology, 2005, 49(1): 71-85.周 勇, 李 飞, 韩现刚. 综放开采覆岩破断与顶煤力链结构演化特性分析[J

采矿与安全工程学报, 2018, 35(4): 751-755, 825.ZHOU Y, LI F, HAN X G. Analysis on evolution characteristics for overburden and chain of top-coal using fully mechanized top coal caving mining[J

Journal of Mining & Safety Engineering, 2018, 35(4): 751-755, 825.孙其诚, 辛海丽, 刘建国, 等. 颗粒体系中的骨架及力链网络[J

岩土力学, 2009, 30(S1): 83-87.SUN Q C, XIN H L, LIU J G, et al. Skeleton and force chain network in static granular material[J

Rock and Soil Mechanics, 2009, 30(S1): 83-87.徐正红. 颗粒物质的力链压曲变形及本构模型研究[D

杭州: 浙江大学, 2011: 21-50.XU Z H. Research of the buckling deformation of force chain and constitutive model of granular materials[D

Hangzhou: Zhejiang University, 2011: 21-50.FANG Y G, GUO L F, HOU M X. Arching effect analysis of granular media based on force chain visualization[J

Powder Technology, 2020, 363: 621-628.YANG Y, CHENG Y M. A fractal model of contact force distribution and the unified coordination distribution for crushable granular materials under confined compression[J

Powder Technology, 2015, 279: 1-9.张 炜, 萧伟健, 袁传牛, 等. 离散元法铁粉末压制中粒径分布对力链演化机制的影响[J

力学学报, 2022, 54(9): 2489-2500.ZHANG W, XIAO W J, YUAN C N, et al. Effect of particle size distribution on force chain evolution mechanism in iron powder compaction by discrete element method[J

Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(9): 2489-2500.IMSEEH W H, JARRAR Z A, ALSHIBLI K A. Influence of sand morphology on interparticle force and stress transmission using three-dimensional discrete- and finite-element methods[J

Journal of Engineering Mechanics, 2021, 147(10): 04021081.党 洁, 张晨洁, 郭生茂, 等. 隔离矿柱回采稳定性影响因素的正交试验分析[J

金属矿山, 2014(5): 24-26.DAND J, ZHANG C J, GUO S M, et al. Orthogonal Tests Analysis on the Influencing Factors for the Stoping Stability of Isolated Pillars[J

Metal Mine, 2014(5): 24-26.黄德镛, 贾子月, 吕世玮, 等. 基于正交设计试验方案的云南某铜矿采场结构参数优化研究[J

有色金属工程, 2023, 13(4): 111-119.HUANG D Y, JIA Z Y, Lü S W, et al. Study on optimization of stope structure parameters of a copper mine in Yunnan based on orthogonal design test scheme[J

Nonferrous Metals Engineering, 2023, 13(4): 111-119.刘 欢, 任凤玉, 何荣兴, 等. 模拟矿岩散体的PFC细观参数标定方法[J

金属矿山, 2018(1): 37-41.LIU H, REN F Y, HE R X, et al. Calibration methods of the PFC microscopic parameters for simulating the loose ore rock[J

Metal Mine, 2018(1): 37-41.IWASHITA K, ODA M. Rolling resistance at contacts in simulation of shear band development by DEM[J

Journal of Engineering Mechanics, 1998, 124(3): 285-292.孙其诚, 王光谦. 颗粒流动力学及其离散模型评述[J

力学进展, 2008, 38(1): 87-100.SUN Q C, WANG G Q. Review on granular flow dynamics and its discrete element method[J

Advances in Mechanics, 2008, 38(1): 87-100.尹成薇, 梁冰, 姜利国. 基于颗粒流方法的砂土宏—细观参数关系分析[J

煤炭学报, 2011, 36(S2): 264-267.YIN C W, LIANG B, JIANG L G. Analysis of relationship between macro-micro-parameters of sandy soil based on particle flow theory[J

Journal of China Coal Society, 2011, 36(S2): 264-267.陈庆发, 陈青林, 仲建宇, 等. 柔性隔离层下单漏斗散体矿岩流动规律[J

工程科学学报, 2016, 38(7): 893-898.CHEN Q F, CHEN Q L, ZHONG J Y, et al. Flow pattern of granular ore rock in a single funnel under a flexible isolation layer[J

Chinese Journal of Engineering, 2016, 38(7): 893-898.陈庆发, 陈青林, 秦世康. 同步充填柔性隔离层下散体介质流理论[M

长沙: 中南大学出版社, 2019: 26-32.CHEN Q F, CHEN Q L, QIN S K. Loose medium flow theory under a flexible isolated layer in synchronous filling mode[M

Changsha: Central South University Press, 2019: 26-32.王学深. 正交试验设计法[J

山西化工, 1989(3): 53-58.WANG X S. Orthogonal experimental design method[J

Shanxi Chemical Industry, 1989(3): 53-58.甘健胜. 概率论和数理统计[M

北京:清华大学出版社, 2005: 219-222.GAN J S. Probability theory and mathematical statistics[M

Beijing: Tsinghua University Press, 2005: 219-222.CHEN Q F, ZHAO F Y, CHEN Q L, et al. Orthogonal simulation experiment for flow characteristics of ore in ore drawing and influencing factors in a single funnel under a flexible isolation layer[J

JOM, 2017, 69(12): 2485-2491.付龙龙, 周顺华, 田志尧, 等. 双轴压缩条件下颗粒材料中力链的演化[J

岩土力学, 2019, 40(6): 2427-2434.FU L L, ZHOU S H, TIAN Z R, et al. Force chain evolution in granular materials during biaxial compression[J

Rock and Soil Mechanics, 2019, 40(6): 2427-2434.张 炜, 周 剑, 于世伟, 等. 双轴压缩下颗粒物质接触力与力链特性研究[J

应用力学学报, 2018, 35(3): 530-537, 687.ZHANG W, ZHOU J, YU S W, et al. Investigation on contact force and force chain of granular matter in biaxial compression[J

Chinese Journal of Applied Mechanics, 2018, 35(3): 530-537, 687.周玉珠. 正交试验设计的矩阵分析方法[J

数学的实践与认识, 2009, 39(2): 202-207.ZHOU Y Z. A matrix analysis of orthogonal design[J

Mathematics in Practice and Theory, 2009, 39(2): 202-207.闫红杰, 夏 韬, 刘 柳, 等. 高铅渣还原炉内气液两相流的数值模拟与结构优化[J

中国有色金属学报, 2014, 24(10): 2642-2651.YAN H J, XIA T, LIU L, et al. Numerical simulation and structural optimization of gas-liquid two-phase flow in reduction furnace of lead-rich slag[J

The Chinese Journal of Nonferrous Metals, 2014, 24(10): 2642-2651.