柔性隔离层下多漏斗散体矿岩力链演化特征的离散元模拟

陈庆发，王少平，秦世康

1) 广西大学资源环境与材料学院，南宁 530004 2) 广西大学土木建筑工程学院，南宁 530004✉通信作者，E-mail: chqf98121@163.com

随着可持续发展理念的贯彻实施，当前采矿领域正朝着绿色开采、无废开采、协同开采等[1−3]方向发展. 在此背景下，笔者提出了一种大量放矿同步充填无顶柱留矿采矿方法[4]. 新采矿方法与传统留矿法的区别在于大量放矿前设置柔性隔离层，使得放矿过程中矿岩流动规律突破了传统放矿理论的描述范围. 为构建一套完整的“同步充填柔性隔离层下散体介质流理论”，笔者对大量放矿同步充填放矿过程做了一系列研究[5−12]. 研究发现矿岩颗粒由大量的非均匀块体组成，这些块体大小、形状均不相同，具有非连续和非均质的基本特征，在一定程度上表现出颗粒物质材料的本质属性. 颗粒物质材料从存在形态上可划分为颗粒气体、颗粒流体和颗粒固体，它们有着本质不同的动量传递机制，其中，后两者属密集颗粒物质体系. 而矿岩颗粒是一种具有固体和液体特性的特殊两相流体，从形态划分上属于密集颗粒物质体系，呈现出对外界微小作用的敏感性、非线性响应、自组织行为等复杂力学特性. 在矿岩颗粒体系内部，由于自身重力和外部荷载作用，颗粒物质在颗粒界面接触处产生接触力，接触力沿着颗粒接触网络传递，并最终形成完整的力链网络. 而力链作为连接颗粒物质宏观、细观特性之间的桥梁，通过力链研究能够进一步加强对柔性隔离层下散体介质流动规律的宏观认识，因此开展柔性隔离层下散体矿岩力链演化特征研究具有重要的理论意义.

近年来，力链作为颗粒物质力学中的研究热点，学者们分别从力学分析、室内试验和数值模拟三个方面展开了一系列研究. 在力学分析方面，Tordesilla等[13]从力链演化和相关运动学的角度研究了应变局部化试样非共轴性的微观力学根源，讨论了单剪试验和双轴试验中应变局部化试样宏观力学特性的关系. Hunt等[14]以结构力学的观点为基础，推测了颗粒介质中的力链屈曲模型.Socolar等[15]提出了基于有向力链网络的二维颗粒材料应力场理论，建立了不同方向力链密度的一般 Boltzmann方程. 在物理实验方面，Estep与Dufek[16]和Chen等[17]采用光弹性技术和数字图像处理方法分别研究了密集颗粒流中颗粒的位移和运动方向，力链的动态力学行为对颗粒流系统基底的影响，并对力链中的损伤点进行了预测分析.Hou等[18]建立了一个平行板剪切单元模拟了无限平行板的剪切运动，分析了力链的演化过程及影响因素，并利用力链的载荷分布率曲线描述了力链的受力行为和形态变化. 在数值模拟方面，Tordesilla[19]基于离散单元法从结构稳定性的角度定量研究了密集颗粒组合中剪切带内力链的屈曲程度和屈曲模式对应力、应变及能量耗散规律的影响. Sun等[20]采用离散单元法对11000个等粒径颗粒进行了二维单轴压缩试验，得到了力链长度分布的幂律形式. Tordesillas等[21]基于离散单元法开展了刚性平板冲头压入致密颗粒材料的数值试验，从压痕阻力、破坏、雷诺剪胀、未变形“死区”四个方面探讨了颗粒介质自组织线性、循环结构块的力链和接触网络循环. Zhang和Xu等[22−24]分别采用离散单元法研究了准静态双轴加载和动态冲击加载条件下，力链屈曲对试件整体力学响应及颗粒材料微观结构的影响，试样在不同时间段的剪切应变、偏应力和体积应变，以及力链强度、长度和方向系数与应力传递关系的定量表征.

综合分析上述文献，以上研究均未涉及放矿领域，并且学者们对力链方向、长度、强度和屈曲等已经展开了初步研究，但是多数研究处于定性描述层面，量化描述研究较少，且量化描述尚无统一的标准.由于力链演化特征的量化研究有助于揭示颗粒物质内部的作用机理，进而加强对其宏观力学行为的认识，因此，为了充分掌握柔性隔离层下散体介质的流动规律，本文基于离散元软件PFC开展了柔性隔离层下多漏斗散体矿岩力链演化特征的数值试验研究.

1 模型构建

利用颗粒离散元软件PFC构建同步充填留矿法试验模型，如图1所示. 试验中的接触模型选用抗滚动线性接触模型，颗粒生成采用自重堆积法.为实现柔性隔离层下同步充填放矿过程，颗粒的生成按如下三个步骤进行：

图1 数值试验中的颗粒模型Fig.1 Model of particles in the test

（1）初始颗粒的生成：通过 ball generate命令在墙体模型y轴方向0.08～130 cm范围内生成若干矿石颗粒，这些矿石颗粒的重力加速度g=9.81 m·s−2，其细观力学参数如表1所示. 为使散体介质体系内矿石颗粒尽快充填密实，初始颗粒的接触模型设置为线性接触模型，颗粒之间的摩擦系数取为0.3；同时为方便观察放矿过程中矿石颗粒的流动现象，待模型平衡后，以10 cm为间隔将颗粒赋予不同的颜色，并删除y轴正向128 cm以上的矿石颗粒.

（2）柔性隔离层的生成：为模拟柔性隔离层，利用Cubic命令，在生成的矿石颗粒上方生成一排长 250 cm，半径为 0.0015 cm 的细小颗粒，这些颗粒采用平行黏结，并赋予隔离层如表2所示的细观力学参数，放矿过程中柔性隔离层长度为178 cm.

（3）计算颗粒的生成：为模拟矿石颗粒大小及形状不同所带来的影响，将颗粒接触模型由线性接触模型变为抗滚动线性接触模型，颗粒的抗滚动线性摩擦系数为0.3（多次调试综合取值），此时散体介质体系内颗粒的细观力学计算参数如表3所示. 删除编号为4的底墙后，放矿口被打开，矿石颗粒从放矿口放出，矿石流动随即开始. 放矿过程中，每计算若干时步，关闭4号放矿口，并在矿石颗粒面上生成适量的充填废石颗粒（其细观力学参数与矿石颗粒细观力学参数相同），为实现同步充填效果，待模型在自重作用下解算平衡后，删除多余的充填废石颗粒，并再次打开4号放矿口，进入下一循环计算过程，直至矿石颗粒全部放出（计算且试验过程中放矿口的开闭方式为全部打开或全部关闭）.

数值模型中各颗粒间接触力采用Hertz-Mindlin理论进行计算，其法向、切向接触力如下[25]：

表 1 墙体及初始矿石颗粒力学参数Table 1 Mechanical parameters of walls and initial ore particles

表 2 隔离层相关参数Table 2 Parameters related to the isolation layer

表 3 矿石颗粒参数Table 3 Parameters of ore particles

墙体是与颗粒相互作用的具有任意可定义接触属性的一个平面，在PFC程序中通过赋予墙体如表1所示的刚度、摩擦系数等性质，生成了可以相互连接在一起的多边形，这些多边形即是数值模型的边界条件，它们将矿石颗粒约束在由其围成的区域范围内.

2 力链的程序识别

力链形成需要满足 3 个条件[25−27]：（1）接触力大小必须大于体系内平均接触力；（2）相邻接触之间的夹角必须小于等于角度阈值；（3）组成力链的颗粒个数必须大于等于3.

根据以上3个成链条件，对于力链识别设置3个判据：

（1）接触力判据.

散体介质体系内平均接触力大小：

接触力判据为：

（2）接触角判据.

在离散元模型中，接触为颗粒中心点之间的连线，假如颗粒A、B、C之间存在两个接触，颗粒A、B之间的接触为AB，其法向量为（αAB，βAB），颗粒B、C之间的接触为BC，其法向量为（αBC，βBC），接触AB与接触BC之间的夹角：

颗粒成链的角度阈值:

颗粒成链的接触角判据为：

（3）力链长度判据.

一条力链上颗粒的个数必须大于等于3，接触的个数必须大于等于2，所以颗粒成链的长度判据为：

式中，L为力链的长度，即一条力链上所含颗粒的个数；Nc为接触的个数，即一条力链上所含接触的个数.

利用PFC2D中的Fish语言输出不同放矿次数下散体介质体系内接触的编号、大小、位置坐标等信息. 基于以上3个判据，利用Matlab软件编制力链识别程序，遍历散体介质体系内所有接触，实现力链的自动检索及识别.

3 散体介质体系内力链演化机制

3.1 力链形态演化规律

基于大量放矿同步充填无顶柱留矿采矿法数值试验模型，利用PFC编制Fish语言模拟柔性隔离层下多漏斗放矿过程，每放出一定量的矿石颗粒，记录并输出矿石颗粒信息，并利用PFC2D中的Geometry命令，将散体介质体系内的强接触进行可视化，如图2所示.

图2 力链宏观分布特征Fig.2 Macroscopic distribution characteristics of the force chain

图2为多漏斗放矿过程中散体介质体系内力链的分布情况，图中蓝色线条为接触力大于散体介质体系内平均接触力的接触. 由于接触力大于散体介质体系内平均接触力的接触仅仅满足颗粒成链的接触力判据，所以图中显示的力链较实际情况偏多；但是这对散体介质体系内力链的宏观分布特征影响较小，因此其可以代表大量放矿同步充填过程中散体介质体系内力链的整体分布特征.

由图2可知，放矿初期，力链主要集中在隔离层下方，力链分布较为均匀，力链网络空隙较大；随着矿石颗粒的持续放出，隔离层同步弯曲下沉，隔离层上方力链逐渐增多，隔离层下方力链逐渐减少，力链主要集中在模型中部，力链网络空隙较大，存在局部的应力集中现象；放矿后期，隔离层上部力链占据主导地位，且呈现出明显的方向性，

主要沿着垂直隔离层的方向分布，隔离层下方力链数目较少，主要沿着垂直漏斗边壁方向分布.

3.2 力链数量和长度变化

根据颗粒成链的3个判据可知，强接触不一定能够组成力链.

图3所示为多漏斗放矿过程中强接触数量与力链接触数量在总接触数量上的占比，由图可知多漏斗放矿过程中强接触及力链接触占比均比较稳定，其中强接触占比稳定在33%左右，力链接触占比稳定在16%左右，上下波动幅度均不超过2%.

图3 多漏斗放矿过程强接触与力链接触占比Fig.3 Proportion of strong contact and force chain contact in multifunnel ore drawing process

图4为多漏斗放矿过程中力链总数的变化规律，由图可知随着放矿次数的增加力链总数不断波动减小，并在放矿后期稳定在790条左右.

图4 多漏斗放矿过程散体介质体系内力链数目变化规律Fig.4 Variation law of the number of force chains in the granular medium system during multi-funnel ore drawing

结合图3和4可知，放矿前期散体介质体系在外载荷作用下发生剧烈波动，散体介质体系内部力链不断断裂重组，并在不断断裂重组的过程中达到一个新的平衡，因此，放矿前期强接触占比及力链总数均波动较大，放矿后期，随着体系内荷载逐渐稳定，强接触及力链接触逐渐趋于稳定.

将每次放矿结束后散体介质体系内所含力链数目进行归一化处理，计算体系内力链长度的概率分布，最终结果如图5和图6所示（不考虑长度大于20的力链，因为其在力链总数中所占概率较小）. 通过观察图5可知，不同放矿次数下力链长度分布均呈指数形式递减. 对第8次放矿结束后散体介质体系内力链长度分布进行拟合，拟合相关系数R2=0.999，拟合函数为：

式中，a=32.26，b=−3.64.

利用上述拟合公式对其它放矿次数下的力链长度进行拟合，拟合效果均较好，相关系数均大于0.99.

由图5和图6可知，不同放矿次数下力链长度的概率分布几乎一致，均随着力链长度的增加呈指数形式递减；其中长度为3的力链数目最多，约占力链总数的55%，其个数随着放矿次数的增加逐渐减少；长度处在4～6之间的力链数目约占力链总条数的37%，其个数随着放矿次数的增加不断波动，但基本保持不变；长度大于6的力链数目较少，约占力链总数的8%，其个数随着放矿次数的增加逐渐增加. 这主要是因为放矿初期，体系原始应力被打破，冲击荷载的作用使得原本稳定的力链网络遭到破坏，长力链断裂弯曲形成短力链，因此，放矿初期散体介质体系内短力链数目较多，长力链数目较少；随着放矿过程的持续推进，体系内荷载逐渐稳定，新的力链分布网络逐渐形成，荷载逐渐被长力链分担，因此，随着放矿次数的增加长力链数目逐渐增多，短力链数目逐渐减少. 而在整个放矿过程中，中等长度力链不断断裂重组，因此其比例基本保持不变.

3.3 力链强度与方向分布

散体矿岩体系中，矿石颗粒密集排布，外载荷及自身重力使得毗邻颗粒相互挤压、接触，形成诸多传递外载荷的路径. 在这些路径中，一些颗粒变形较大，传递较大份额的外载荷，形成强力链；一些颗粒变形微弱，传递小部分的外载荷，形成弱力链. 因此，通常把力链承载能力的强弱称为力链强度.

各个阶段散体介质体系内力链强度演变规律如图7所示.

由图7可知放矿初期散体介质体系内部力链强度发生了较大变化，第4次放矿结束之后力链强度波动范围逐渐减小并趋于稳定.

根据每次放矿结束后散体介质体系内力链强度的统计结果，计算体系内力链强度的概率分布，最终结果如图8所示. 通过观察图8可知，多漏斗放矿过程中力链强度的概率分布随着放矿次数的增加先呈指数形式上升再呈指数形式下降，并在0.7处出现一峰值. 对第1次放矿结束后散体介质体系内力链长度分布进行拟合，拟合相关系数R2=0.999，拟合函数为：

图5 多漏斗放矿过程力链长度概率分布Fig.5 Probability distribution of the force chain length in the multi-funnel ore drawing process

图7 多漏斗放矿过程力链强度演化规律Fig.7 Evolution law of the force chain strength in the multi-funnel ore drawing process

式中，A=2.528，B=2.279，C=4.788，D=2.307.

利用以上拟合公式对其他放矿次数下的力链强度进行拟合，拟合效果均较好，相关系数均大于0.99.

3.4 力链方向分布特征

将360°等分为36个区间，对每个区间内的力链数目及强度进行统计，进而求出多漏斗放矿过程中每个区间内力链的平均强度，最终得到多漏斗放矿过程中力链方向的分布演化规律如图9所示.

为了定量描述大量放矿同步充填过程中散体介质体系内部力链细观组构特征的演化规律，采用Rothenburg与Bathurst[28]提出的颗粒间接触力角度分布描述函数对放矿过程散体介质体系内力链强度与力链方向统计结果进行拟合：

式中，fn(θ) 为散体介质体系内力链强度的分布函数，f0为散体介质体系内所有力链强度的平均值，an为三角函数振幅即傅里叶系数，其值表示散体介质体系内力链方向分布的各向异性程度，θn为散体介质体系内力链分布的主方向，w为三角函数频率控制系数，三角函数的最小正周期T=2π/(π/w)=2w,n为接触法线方向的单位向量.

通过对比图9中红色与蓝色线条轮廓的相似程度可以看出，前6次拟合效果较好，后6次拟合效果较差，最终获得的拟合参数如表4所示.

综合图9和表4，由不同放矿次数下力链方向分布对比分析可得，前6次放矿过程中散体介质体系内力链方向分布的变化规律具有一致性，后6次放矿过程中力链方向分布的变化规律具有一致性. 放矿初始阶段，散体介质体系处于初始固结状态，系统受自身重力影响较大，力链主要沿铅垂方向分布，故散体介质体系内力链分布主方向θn在90°左右，力链方向分布形态近似花生状；此后随着矿石颗粒的持续放出，受回填废石上覆荷载及隔离层形状的影响，散体介质体系内部力链方向分布形态发生较大变化，局部应力集中现象明显，在图中表现为红色线条与拟合得到的蓝色线条相比在某些部位特别突出，拟合效果不好，且从图中可以明显看出力链分布主方向由一个演变为4个（铅垂方向、水平方向及与水平方向呈±60°夹角方向）. 对于力链分布各向异性程度表征参数an而言，第6次放矿结束之前，随着放矿次数的增加，an逐渐增大；第6次放矿结束之后，受隔离层横向摩擦以及隔离层形状的影响，力链主方向分布由一个演变为6个，所以体系内力链各向异性程度an逐渐减小，并在第10放矿结束之后达到最小值；第10次放矿结束之后，由于局部应力集中现象的影响，各放矿口上方颗粒并没有同时放出，隔离层形状的变化造成了与水平方向夹角为+60°夹角的力链明显增多，相应的各向异性系数也有所增加. 对于体系内力链平均强度f0而言，经过归一化处理，体系内力链的平均强度为1.00，所以其一直保持不变.

图8 多漏斗放矿力链强度的概率分布Fig.8 Probability distribution of the force chain strength in the multi-funnel ore drawing process

3.5 力链准直性分析

力链准直系数反映了力链是否能较好地形成直线链状结构. 将体系内所有力链的准直系数进行累加求和，并根据力链数目进行几何平均，得到每次放矿结束后散体介质体系内力链的准直系数，如图10所示.

由图10可知，多漏斗放矿过程中力链准直系数随着放矿次数的增加呈指数形式增长并逐渐保持稳定，利用式（11）对其进行拟合，结果如图10中的红色曲线所示. 拟合系数为0.945，拟合效果较好. 式（11）如下所示：

式中：y0=0.874，A=0.127，xc=9.033，w=0.833.

如果力链强度较大，准直性较好，就会使更多的颗粒连接在一起，因此长力链一般准直性较好、强度较大，短力链一般长度较小、强度也较小. 结合多漏斗放矿过程中力链长度、数量、强度、方向等的变化规律可知，放矿初期，散体介质体系主要受自身重力及放矿口等效荷载的影响，荷载分布均匀，每个力链所承受的荷载较小，所以放矿前期力链长度较短，数目较多，力链准直系数较小；随着大量放矿同步充填试验的持续进行，隔离层上方回填废石越来越多，相应的上覆载荷越来越大，所以隔离层下方力链长度逐渐增强，准直系数也逐渐增加.

图9 多漏斗放矿过程力链方向演化规律Fig.9 Evolution law of the force chain direction in the multi-funnel ore drawing process

表 4 多漏斗放矿过程力链方向分布的拟合结果Table 4 Fitting results of the force chain direction distribution in the multi-funnel ore drawing process

图10 力链准直系数变化规律Fig.10 Change law of collimation coefficient of the force chain

3.6 讨论

多漏斗放矿条件下，因各放矿口间的相互影响而产生交错、缺失等不同程度的变异，放出体并不是一个规则的椭球体，并且柔性隔离层曲线整体形态为正弦曲线. 放矿初期，散体介质体系主要受自身重力及放矿口等效荷载的影响，荷载分布均匀，整个隔离层界面近似在同一水平，并保持平缓下移，每个力链所承受的荷载较小，所以放矿前期力链长度较短，数目较多，力链准直系数较小.随着大量放矿同步充填试验的持续进行，隔离层上方回填废石越来越多，当最高矿石面下降到某一深度时，最高矿石面开始出现凹凸不平. 随着矿石面的继续下降，凹凸现象越来越明显，相应的上覆载荷越来越大，此时力链承受的荷载较大，所以隔离层下方力链长度逐渐增强，准直系数也逐渐增加，而柔性隔离层则出现较大的弯曲起伏现象，并呈圆弧形下移，直至放矿终了以波浪形悬浮于底部结构.

4 结论

（1）多漏斗放矿过程中强接触及力链接触占比均比较稳定，其中强接触占比稳定在33%左右，力链接触占比稳定在16%左右，上下波动幅度均不超过2%.

（2）力链总数随着放矿次数的增加不断波动减少，并在放矿后期稳定在790条左右. 不同放矿次数下力链长度的概率分布几乎一致，均随着力链长度的增加呈指数形式递减；其中长度为3的力链数目最多，约占力链总数的55%，其个数随着放矿次数的增加逐渐减少；长度处在4～6之间的力链数目约占力链总条数的37%，其个数随着放矿次数的增加不断波动，但基本保持不变；长度大于6的力链数目较少，约占力链总数的8%，其个数随着放矿次数的增加逐渐增加.

（3）放矿初期散体介质体系内部力链强度发生较大变化，第4次放矿结束之后力链强度波动范围逐渐减小并趋于稳定. 放矿过程中力链强度的概率分布随着放矿次数的增加先呈指数形式上升再呈指数形式下降，并在0.7处出现一峰值.

（4）前6次放矿过程中散体介质体系内力链方向分布的变化规律具有一致性，后6次放矿过程中力链方向分布变化规律具有一致性. 放矿初始阶段，力链主要沿铅垂方向分布，力链方向分布形态近似花生状；此后随着矿石颗粒的持续放出，散体介质体系内部局部应力集中现象明显，力链分布主方向由一个演变为4个（铅垂方向、水平方向及与水平方向呈±60°夹角方向）. 对于力链分布各向异性程度表征参数an而言，第6次放矿结束之前，随着放矿次数的增加，an逐渐增大；第6次放矿结束之后，体系内力链各向异性程度an随着放矿次数的增加逐渐减小，并在第10放矿结束之后达到最小值；第10次放矿结束之后，各向异性系数又随着放矿次数的增加逐渐增加. 体系内力链平均强度f0在放矿过程中一直保持在1.00.

（5）多漏斗放矿过程中力链准直系数随着放矿次数的增加呈指数形式增长并逐渐保持稳定.