离散元法在球磨机仿真中的研究进展及应用

裴英杰　肖庆飞　石贵明　周强

摘要：球磨机作为矿石磨碎阶段的重要设备，對选矿厂提效降耗起着决定性作用。球磨机功耗与能耗分布、技术参数与衬板结构参数设计优化是影响球磨机磨矿效果的重要内容和技术难点，离散元法可以准确模拟颗粒运动的形变过程并记录颗粒的功率及能量分布，在球磨机仿真领域具有良好的应用前景。介绍了离散元法在球磨机仿真过程中应用比较广泛的颗粒模型、接触模型及联合仿真方法各自的特点及研究进展，以及离散元法在球磨机功率及能量分布研究、技术参数优化、衬板结构参数优化等方面的应用情况，归纳总结并指出了球磨机仿真模型精细化、联合仿真广泛应用、球磨机技术参数智能调节是未来球磨机仿真研究的发展趋势。

关键词：离散元法;球磨机;仿真;颗粒模型;接触模型;耦合方法

中图分类号：TD453文献标志码：A开放科学（资源服务）标识码（OSID）：

文章编号：1001-1277（2022）01-0069-09doi：10.11792/hj20220112

基金项目：国家自然科学基金地区科学基金项目（51964044）;矿冶过程自动控制技术国家（北京市）重点实验室开放基金项目（201905）;安徽省重点研究和开发计划面上攻关项目（201904a07020044）;云南省教育厅课题（2019J0738）

引言

球磨机作为矿石磨碎阶段的重要设备，具有结构简单、适应性强等优点，现已广泛应用于国内外各类矿山[1]。近年来，随着矿石“贫细杂”程度不断加深，选矿厂生产对球磨机设计提出了更高的要求。若球磨机参数设计不合理，会对磨矿效果产生消极作用，严重影响后续选别指标，并加剧钢材及电力的损耗[2-3]。为此，学者们进行了一系列的优化研究。在对球磨机研究的过程中，采用的研究手段主要分为2类：球磨机试验与模拟仿真[4-5]。在最初的研究中，基本研究过程是在理论计算后进行球磨机磨矿试验。近年来，随着计算机技术快速发展，模拟仿真方法越来越多地被应用到球磨机的研究中。目前，球磨机试验和模拟仿真这2种主要研究手段互相验证、相辅相成，共同推动球磨机研究领域的快速发展。

离散元法（Discrete Element Modelling，DEM）可以跟踪球磨机中单一颗粒的运动并准确模拟颗粒之间的碰撞，已成为球磨机仿真过程最常用的数值模拟方法[6]。为此，众多学者对球磨机功率及能量分布、技术参数及衬板结构参数优化等方面进行了广泛的模拟仿真研究，获得了大量研究成果，对后续研究人员了解球磨机内部的颗粒运动及优化各类参数具有积极的促进作用[7-9]。然而，受计算机能力限制，大多数仿真研究均经过模型简化处理，导致模拟过程存在一定误差[10]。适宜的建模仿真可在一定程度上降低仿真误差，提高球磨机仿真精确性。基于此，本文针对离散元法在球磨机仿真中的研究进展及应用进行了归纳总结，并提出了今后离散元法在球磨机仿真研究中的发展方向，以期为球磨机仿真研究提供一定的研究思路和理论指导。

1离散元法在球磨机仿真中的研究

离散元法于1971年由CUNDALL[11]基于分子动力学理论提出，最初被用来研究单颗粒与离散颗粒群之间的交互作用及其稳定性，由于该方法可以模拟岩体非均质、不连续和大变形等特点，逐渐开始应用于岩土力学领域。1992年，离散元法被MISHRA等[12]首次应用到球磨机衬板结构与功耗之间关系的研究中，但受到当时技术水平及计算机能力的限制，仅能进行二维建模仿真。1997年，AGRAWALA等[13]采用离散元法对比分析试验测量与预测所得球磨机功率的差异性，验证了三维离散元法用于球磨机仿真研究的可行性，为离散元法在球磨机仿真中的应用奠定了良好基础。

1.1颗粒模型

颗粒模型是指根据介质和矿石颗粒的形状及特性建立的仿真模型。自19世纪80年代以来，学者们提出包括球体、椭球体[14]、超二次曲面[15]、圆柱体[16]、多面体[17]、多聚球体[18]及多聚椭球体[19]等不同形状的离散元颗粒模型。为了达到不同的仿真要求，学者们选择的颗粒模型也各有差异。离散元法模拟仿真应用的部分颗粒模型见表1。

目前，在球磨机仿真中最常用的颗粒模型主要是球体、多聚球体这2类。球体颗粒模型由于具有计算成本低的优势，常作为球磨机离散元仿真中的首选颗粒模型，而多聚球体颗粒模型具有可模拟真实颗粒形状的优点，常被用来模拟实际矿物，二者最大的差异表现在仿真时间及仿真精度上。BIBAK等[21]通过粒子填充算法对颗粒建模（步骤见图1），探究多聚球体颗粒模型和球体颗粒模型动力学之间的差异，研究发现多聚球体颗粒在高能冲击作用下的参与度较球体颗粒高5 %，多聚球体颗粒的模拟时间成本较球体颗粒增加35倍。BARRIOS等[27]为铁矿石球团建立了球体颗粒模型和多聚球体颗粒模型（4个球体），对比分析2种预测模型的休止角及磨矿功率。结果表明：采用多聚球体颗粒模型仿真测得的休止角角度略低于试验测得的休止角角度，预测磨矿功率与实测磨矿功率差值小于6 %;当仿真模型为球体颗粒模型时，休止角预测值误差过高，预测磨矿功率与实测磨矿功率偏差高达30 %。因此，在仿真条件允许的情况下，采用多聚球体颗粒模型仿真得到的试验数据及结论准确性更高。

在实际磨矿过程中，可能存在数十亿个不同粒度的矿粒在球磨机中同时运动，根据现阶段计算机性能，在离散元仿真中运行数千万个粒子颇具难度，模拟实际尺寸及数量的矿石颗粒几乎没有可能。IWASAKI等[28]研究发现，当球磨机内冲击能量的频率分布及其在研磨时间内的累积量在不同尺寸球磨机中一致时，可以成功实现球磨机及颗粒放大，节约仿真成本。POWELL等[29]将球磨机筒体长度按比例缩小，球磨机中颗粒数量由450万降至11万，模拟时间大幅降低，有效节约了仿真成本。因此，在确保仿真精度的前提下，对球磨机及颗粒进行适当缩放，将极大程度地减少离散元仿真所消耗的时间并有效提高计算效率。

1.2接触模型

离散元仿真過程即颗粒模型运动过程，颗粒运动时必然会引发颗粒之间的相互碰撞，对颗粒之间接触作用研究成为离散元仿真过程的重中之重。目前，在离散元仿真中主要采用的2种接触模型是硬球模型和软球模型。硬球模型的思想是ALDER等[30]在分子系统的相变研究中提出的，该模型认为颗粒间的碰撞是两体瞬时碰撞，且颗粒具有绝对刚度，但其常用于气体输送等稀疏颗粒系统，不适用于球磨机这类颗粒连续接触的致密颗粒系统。目前，球磨机离散元仿真中最常用的是由CUNDALL等[31]提出的软球模型，该模型通过对受力后发生的形变及断裂进行分析来处理颗粒间的接触力，颗粒间的接触为柔性接触，允许存在少量重叠，适用于准静态颗粒材料，对颗粒数量巨大的球磨机离散元仿真系统具有计算效率上的优势。

颗粒接触时接触力的计算方法，即颗粒接触本构模型，是离散元仿真颗粒接触模型的理论核心。目前，应用最广泛的离散元颗粒接触本构模型是CUNDALL等[31]提出的线性弹簧-阻尼器接触模型（见图2），该模型在接触点假设一组弹簧、阻尼器和滑块组件，用来表示2个颗粒之间的碰撞，同时考虑了离散颗粒之间接触力和非接触力的相互作用。在后续的离散元仿真研究中，学者们提出了多种不同的本构模型。JOHNSON等[32]基于Hertz-Mindlin模型，为颗粒添加一个额外的法向内聚力，并考虑接触物体的表面黏附能对接触半径的影响，提出了JKR模型。WALTON等[33]采用部分锁定弹性形变模型计算法向力，并依据MINDLIN等[34]的理论在切向加入了摩擦力对弹性体的作用，提出了半锁弹簧模型。IWASHITA等[35-36]分析了滚动摩擦机制，提出了抗滚动的离散元接触模型，其便于观察颗粒滚动的影响。

虽然颗粒接触本构模型众多，但并非全部都可以契合球磨机仿真研究。为了开发更适合于球磨机仿真的颗粒接触本构模型，有学者开展了相关研究。 IWASAKI等[9]忽略球碰撞时滚动摩擦的影响，提出了一种考虑球磨机转速率与球摩擦系数相关性的三维湿式球磨机接触力计算方法，并验证了该模型的有效性。FENG[37-38]基于接触体积模型开发出可以表示任意颗粒形状的节能法向接触模型，该模型可以自动确定任意颗粒形状的完整法向接触模型，包括接触法线方向、接触点（线）及作用力大小，无需任何其他假设，但使用该模型进行DEM仿真的整体计算效率并没有实际数据可以验证，需要后续学者进行相关工作。

1.3离散元耦合方法

球磨机仿真研究中将离散元法作为最常用的数值分析手段，但仅使用单一离散元法局限性较大。因此，学者们开发了不同的离散元耦合方法，其中最常用的是将离散元法与有限元（FEM）、计算流体力学（CFD）及光滑粒子流体力学（SPH）等方法进行耦合。由于以上方法均有各自的特点，它们与离散元法耦合后进行的研究也有所不同。

离散元法与有限元（DEM-FEM）耦合法通常采用DEM分析球磨机内矿石与介质的运动特性，采用FEM分析球磨机的衬板应力及扭矩。DEM-FEM耦合法分为单向耦合法和双向耦合法，目前单向耦合法已趋于成熟，而双向耦合法仍在进一步研究，因此大部分研究常采用单向耦合法分析磨矿介质、矿石颗粒和衬板提升条之间的相互作用。李昆塬等[39]采用DEM-FEM耦合法对Φ4.3 m球磨机进行模拟仿真，分析了球磨机中不同面角、高度和顶边长度的梯形衬板总变形程度与应力的关系。JONSḖN等[40]使用FEM对球磨机中的柔性橡胶衬板及提升条进行建模，使用DEM对磨矿介质进行建模，对比分析了仿真预测的挠度值与实际挠度值的关系。研究发现，DEM-FEM耦合法仿真获得的挠度曲线与测量值具有一致性，为后续学者采用DEM-FEM耦合法模拟球磨机中可形变的橡胶衬板提供了理论基础。

湿磨作业仿真需要对矿浆性质进行模拟，且由于入料矿石的粒径较小，无法使用单一离散元法对大量浆料进行建模。离散元法与计算流体力学（DEM-CFD）耦合法则可以通过描述物质经过网格坐标点时的密度、能量等来表示矿浆运动，且仿真精度高。BEINERT等[41]使用DEM-CFD耦合法对行星式球磨机进行数值研究，分析了磨矿介质的运动、每个时间段的接触及接触时的动能。MAYANK等[42]采用CFD对连续矿浆流建模，采用DEM对离散矿粒进行建模，通过DEM-CFD耦合法对球磨机内的矿石颗粒及矿浆进行模拟仿真，发现仿真结果与基于PEPT试验结果吻合度高。但是，该试验未考虑矿石颗粒与矿浆之间的相互作用，只考虑了矿石颗粒对矿浆的影响，模型需要进一步完善。

在SPH中粒子之间的关系根据当前时刻粒子的位置确定，其随时间变化而改变，因此该方法亦可用于流体等形变较大的物体，且不会产生类似于拉格朗日网格法（如FEM）严重的网格变形导致数值误差大的问题。离散元法与光滑粒子流体力学（DEM-SPH）耦合法也常被用于研究矿石颗粒、磨矿介质与矿浆之间的交互作用。SINNOTT等[43]使用DEM对矿石颗粒和介质建模，使用SPH对矿浆建模，提出一种球磨机中磨矿介质和矿浆完全耦合流动的预测模型，为后续学者研究球磨机中的固相和液相提供了理论基础。

为了更好地对球磨机运动进行仿真研究，JONSN等[44]将SPH、DEM和FEM 3种方法进行耦合，用DEM和SPH粒子表示介质及矿浆，同时使用FEM对球磨机结构进行建模，通过观察矿浆黏度和密度的差异对扭矩及矿浆压力分布的影响，考察球磨机模型及矿浆模型的可靠性，结果表明SPH-DEM-FEM耦合法为预测球磨机中高能区的体积及优化球磨机衬板提升条和运行参数提供了可能。

2离散元法在球磨机研究中的应用

随着离散元法研究的不断深入，其理论体系日渐成熟，该方法在球磨机研究中也得到了更加广泛的应用，其凭借非连续性的特点，将理论公式或常规试验无法解释的内容更加直观准确地表现出来。离散元法在球磨机研究中主要应用在功率及能量分布研究、技术参数优化、衬板结构参数优化等方面，具体研究流程见图3。

2.1功率及能量分布研究

球磨机运行过程通常需要消耗大量能量，由于磨矿过程中矿石对磨矿介质的破碎行为是随机的，磨矿介质可能与矿石发生碰撞，也可能与其他介质或衬板发生碰撞。当磨矿介质与其他介质或衬板发生碰撞的概率较高时，整个磨矿过程中大部分能量都在做无用功，导致能量利用率偏低，造成能源浪费。球磨机的功率及能量分布是由磨矿过程中多种因素决定的，衬板结构及磨矿技术参数均会对其产生影响，研究不同因素对磨矿过程中功率及能量分布的影响，不仅可以优化产品指标，对选矿厂节能降耗也至关重要。

KANO等[45]探究了球磨机磨矿过程中的冲击能量与磨矿速率常数之间的相关性，发现能量利用率与磨矿速率成正比。CLEARY等[46]对实验室球磨机建模，对比分析了不同料球比对球磨机功率消耗及能量利用方式的影响。研究发现，随着料球比的增加，矿石磨矿消耗能量比例呈线性增加，且随着矿石填充量的增加，磨矿产品粒度变得更细。IWASAKI等[9，47]基于离散元法探究了球磨機不同转速率与介质充填率对球磨机消耗能量的影响。研究发现，球磨机消耗能量不仅与介质充填率成正比，而且随球磨机速度的增加呈线性增加。

传统的能量评估指标均假定球沿轴向均匀分布，然而实际情况并非如此，若假定球沿轴向均匀分布无法准确地反映球磨机的磨矿效果。HUANG等[48]将试验研究和离散元法相结合，通过在球磨机中加入不同尺寸的钢球，发现球磨机中大小不同的钢球往往会彼此分开，而相同类型的钢球会沿轴聚集（见图4）。因此，有必要通过区域总能量来研究球磨机的能量分布情况，以分析球磨机中不同区域的磨矿效果。

学者们常将碰撞能量作为缩放仿真模型的参照，以达到提高仿真效率的目的。依据能量对球磨机筒体按比例缩放的方法，最早是由BOND[49]提出的基于能量的经验性缩放球磨机尺寸的方法。但是，该方法缺乏理论基础，最终磨矿产品粒度分布预测值与实际值存在约20 %的误差。学者们基于碰撞能量对缩放模型的研究均是在该模型的基础上不断完善，随着研究的不断深入，球磨机仿真模型的缩放模型研究已基本成熟。JIANG等[10]建立了不同半径球磨机中颗粒运动的DEM模型，对单个粒子的运动轨迹和能量变化进行分析，发现采用减小模型尺寸并改变颗粒粒度的方法，可以减小模拟计算的规模并提高计算效率。SHI等[50]开发了一种基于能量的球磨机尺寸减小模型，该模型采用基于能量的特定方程确定磨矿产品粒度分布，且充分考虑球磨机的几何形状、装料量及球磨机的生产能力等运行条件的影响。IWASAKI等[28]运用DEM对钢球的冲击能量进行数值分析，探究球磨机直径与转速率对冲击能量的影响。研究发现，在不同尺寸的球磨机中，当冲击能量的频率分布与能量累积量一致时，可以实现球磨机及颗粒的放大。

2.2技术参数优化

在球磨机运行过程中，影响磨矿效果的因素众多，但技术参数最易调整，通过优化技术参数改善磨矿效果是成本最低的方式。分析技术参数对球磨机中介质运动状态及磨矿产品粒度组成的影响，对降低选矿厂能耗及提高产品质量具有重要意义。离散元法对球磨机的不同技术参数进行微观模拟，可以更直观地了解球磨机内部运动情况，有助于对球磨机技术参数进行优化，具体优化研究流程见图5。

OROZCO等[51]使用键合单元法开展离散元仿真研究，以矿石平均粒度和比表面积为评价指标，对比分析了钢球数量与球磨机转速率对磨矿过程的影响（见图6）。研究发现，磨矿速度随相同尺寸钢球数量的增加而增大，但钢球数量过多会引起钢球之间的非弹性碰撞，增加能耗、降低磨矿效率。李腾飞等[8]通过DEM对球磨机建模，研究不同转速率下球磨机内介质的运动状态及不同粒度矿石与介质之间的碰撞情况，发现球磨机内介质的拋落角、下落高度、落回瞬间的动能及颗粒之间的碰撞能均与球磨机转速率有关，且转速率越高，越有利于破碎大尺寸及较硬矿石。d/d0—产品与给矿平均粒径比值S/S0—产品与给矿平均比表面积比值n—球磨机转数Nb—钢球数量

球磨机磨矿技术参数之间大多有交互作用，仅对一种参数进行分析通常很难达到优化的目的，因此需要多方面考虑。CLEARY[52]采用DEM探究了充填率、介质及矿石粒度分布对球磨机功耗及矿石运动状态的影响。孙鹏[53]采用DEM研究了充填率、转速率、钢球尺寸及料球比对球磨机单位功耗和介质与物料间冲击能量的影响。刘刚等[54]针对转速率、充填率、球料比等参数，以磨矿工作效率、破碎能力为性能指标，采用均匀试验方法，对球磨机内的物料进行了仿真运动分析及技术参数优化。

球磨机技术参数优化的另一个研究热点是球磨机预测模型。造成磨矿效果差的部分原因是缺乏监控及预测球磨机磨矿效果的适当方法，使得磨矿过程很难被控制。而离散元法能够有效预测球磨机磨矿效果及能耗[18]。PEDRAVES等[55]为有效控制球磨机磨矿过程，充分考虑球磨机转速率、待磨颗粒性质及充填率等因素，对DEM仿真生成的扭矩数据频谱进行分析，研究发现，球磨机的负载扭矩信号足以体现球磨机的负载水平，并提出了一种基于扭矩数据频谱分析球磨机充填率的方法。DE CARVALHO等[56]以矿石粒度分布、转速率、充填率、球磨机直径及介质尺寸为自变量，建立预测矿石破碎率的函数模型，研究发现，使用该模型的预测结果与HERBST等[57]提出的按比例关系放大程序得到的结果一致，预测效果明显。

2.3衬板结构参数优化

球磨机衬板是矿石及介质提升的动力来源，球磨机的衬板结构对矿石颗粒的运动状态及磨矿效果有很大的影响。衬板的结构种类众多，根据材质和形状不同主要分为8类（见图7）：波形衬板、梯形衬板、阶梯形衬板、楔形衬板、平衬板、波纹形衬板、橡胶衬板、磁性衬板。在实际生产中，球磨机衬板提升条形状及数量等参数设计不合理时，会导致球磨机出现磨矿效率低、钢球消耗高及磨矿产品质量差等问题，对产品指标产生恶劣影响。

ROSALESMARḮN等[58]对4种不同面角衬板提升条进行离散元仿真，分析球磨机衬板几何形状对矿石破碎率和衬板提升条磨损率的影响。CHIMWANI等[59]采用离散元法对球磨机衬板形状进行了优化，发现沿球磨机轴向逐渐改变衬板提升条形状可以实现不同的破碎及磨矿要求。WU等[60]分别对不同尺寸球磨机的衬板磨损情况开展仿真试验研究，发现大型球磨机中钢球对衬板的冲击作用与小型球磨机相差无几。张学东等[61]基于仿生耦合理论和磨粒磨损相关理论，依据生物耐磨特征的尺寸参数和比例关系，分别对不同的非光滑衬板提升条形状和现有衬板提升条形状开展仿真对比试验，发现建立的生物非光滑衬板提升条较现有衬板提升条耐磨效果更好。

当现场衬板形状设计不合理或者磨损严重导致球磨机磨矿效果差时，若无法对球磨机内的产品进行实时观测，易造成严重的能源浪费，影响磨矿产品指标。因此，针对球磨机内衬板对介质及矿石颗粒的影响，实时精确控制球磨机内的磨矿过程及对衬板磨损进行预测具有重要意义。POWELL等[29]将磨矿速度作為预测球磨机衬板提升条磨损的影响因素，提出一种基于磨矿速度预测球磨机衬板磨损演变的方法。BOEMER等[62]对球磨机的轴向切片进行三维DEM仿真，开发了一种基于DEM预测衬板表面磨损演变过程的程序，使球磨机衬板仿真研究变得更直观，但衬板形状变化的精度受到计算机能力的限制。REZAEIZADEH等[63]基于影响衬板磨损的主要参数，开发了一种预测衬板磨损率的模型，该模型能根据球磨机物料的密度、几何特性和球磨机转速率的变化，预测不同衬板提升条的使用寿命。

综上可见，国内外学者应用离散元法在功率及能量分布研究、技术参数优化、衬板结构参数优化等方面积累了大量研究经验。在中国选矿厂提效降耗的背景下，该方法在球磨机仿真领域具有良好推广应用前景。

3离散元法在球磨机仿真中的发展趋势

目前，全球对矿石的需求量仍居高不下，矿产资源的开发利用逐渐趋向于集约化、规模化和智能化，颠覆传统的生产技术与研究方法不断革新。如前文所述，离散元法在球磨机设计优化中具有不可替代的作用，其加快了球磨机的发展进程及新产品研制的速度。随着计算机行业的快速发展及各种先进技术的应用，该方法必将展示出更强大的设计指导能力及试验替代能力。但是，目前该方法在球磨机仿真领域仍处于探索阶段，数据质量、测量精度等技术方面的短板依然饱受业界质疑，本文参照众多学者的研究结果，结合当前离散元法在球磨机仿真中的研究难题，梳理出在球磨机研究领域离散元法的模型精度、仿真方式及技术参数控制3个方面的发展趋势：

1）球磨机仿真模型的精细化发展。目前，球磨机仿真受到计算机能力限制，大多仿真研究对矿石颗粒的运动及衬板结构建立了简化模型。然而，随着计算机处理能力的提高、软件算法的优化，球磨机及磨矿系统的仿真模型将更加符合实际，具体表现为矿石颗粒形状及颗粒间的接触更符合实际、矿石颗粒粒度分布更确切、球磨机结构部件更真实、仿真效果更适合指导现场。

2）多种仿真方法联合仿真的广泛应用。根据现阶段研究，离散元法与其他数值仿真方法耦合的目的在于提高仿真精度、对球磨机进行综合分析。多种仿真方法联合仿真可以有效解决单一方法的局限性，对矿浆及球磨机结构进行综合研究，获得更加准确的仿真结果，为球磨机发展开拓新的视角。

3）球磨机技术参数的智能调节。目前，离散元法在球磨机中的仿真研究主要集中于能耗、技术参数及结构性能的优化研究。仿真模拟局限于球磨机的性能预测、衬板磨损失效及磨矿机理的研究，但对球磨机运行时的参数自适应调节研究较少。工业大数据、人工智能、数字孪生等智能技术能够为球磨机仿真研究提供新思想、新思路。通过建立球磨机磨矿参数的智能调节系统，可实时监测球磨机的运行状态，并动态调整技术参数，提高球磨机的磨矿效果及衬板和介质的使用寿命。

4结论

自离散元法首次在球磨机仿真中应用至今已有近30年的时间，在此期间研究人员对球磨机的离散元仿真进行了各种研究，加快了球磨机的发展进程及新产品研制的速度，使球磨机的理论研究及应用均取得了长足进展，但目前离散元法在球磨机仿真中仍然存在着模型精度低及仿真时间长等问题。因此，本文重点论述了各离散元颗粒模型、接触模型和联合仿真方法的特点，并从功率及能量分布研究、技术参数优化、衬板结构参数优化3个方面介绍了离散元法在球磨机仿真中的应用情况，提出球磨机仿真模型的精细化发展、球磨机联合仿真方法的广泛应用、球磨机技术参数的智能调节是离散元法在球磨机仿真应用中未来的发展趋势。随着工业大数据、人工智能、数字孪生等智能技术交叉融合，离散元法在球磨机仿真领域的智能创新应用，势必成为智慧矿山体系建设中不可或缺的重要技术手段。

[参 考 文 献]

[1]JAYASUNDARA C T，ZHU H P.Impact energy of particles in ball mills based on DEM simulations and data-driven approach[EB/OL].（2021-09-24）[2021-10-05].https：∥www.sciencedirect.com/science/article/abs/pii/S0032591021008561.

[2]WANG T，ZOU W J，XU R J，et al.Assessing load in ball mill using instrumented grinding media[J].Minerals Engineering，2021，173：107198.

[3]YU J W，QIN Y H，GAO P，et al.An innovative approach for determining the grinding media system of ball mill based on grinding kinetics and linear superposition principle[J].Powder Technology，2021，378（Part A）：172-181.

[4]GAO P，ZHOU W T，HAN Y X，et al.Enhancing the capacity of large-scale ball mill through process and equipment optimization：an industrial test verification[J].Advanced Powder Technology，2020，31（5）：2 079-2 091.

[5]MARIJNISSEN M J，GRACZYKOWSKI C，ROJEK J.Simulation of the comminution process in a high-speed rotor mill based on the feeds macroscopic material data[J].Minerals Engineering，2021，163：106746.

[6]耿凡，袁竹林，孟德才，等.球磨機中颗粒混合运动的数值模拟[J].热能动力工程，2009，24（5）：623-629，682.

[7]BIAN X L，WANG G Q，WANG H D，et al.Effect of lifters and mill speed on particle behaviour，torque，and power consumption of a tumbling ball mill：experimental study and DEM simulation[J].Minerals Engineering，2017，105：22-35.

[8]李腾飞，林蜀勇，张博，等.不同转速率下球磨机内钢球的碰撞研究[J].中南大学学报（自然科学版），2019，50（2）：251-256.

[9]IWASAKI T，YAMANOUCHI H.Ball-impact energy analysis of wet tumbling mill using a modified discrete element method considering the velocity dependence of friction coefficient[J].Chemical Engineering Research and Design，2020，163：241-247.

[10]JIANG S Q，YE Y X，TAN Y Q，et al.Discrete element simulation of particle motion in ball mills based on similarity[J].Powder Technology，2018，335：91-102.

[11]CUNDALL P A.A computer model for simulating progressive large-scale movements in blocky rock systems[C]∥ISRM.Proceedings of the symposium of the international society for rock mechanics.Lisbon：ISRM，1971：128-132.

[12]MISHRA B K，RAJAMANI R K.The discrete element method for the simulation of ball mills[J].Applied Mathematical Modelling，1992，16（11）：598-604.

[13]AGRAWALA S，RAJAMANI R K，SONGFACK P，et al.Mechanics of media motion in tumbling mills with 3d discrete element method[J].Minerals Engineering，1997，10（2）：215-227.

[14]HE S Y，GAN J Q，PINSON D，et al.Particle shape-induced axial segregation of binary mixtures of spheres and ellipsoids in a rotating drum[J].Chemical Engineering Science，2021，235：116491.

[15]YANG S L，WANG H，WEI Y G，et al.Flow dynamics of binary mixtures of non-spherical particles in the rolling-regime rotating drum[J].Powder Technology，2019，361：930-942.

[16]FENG Y T，HAN K，OWEN D R J.A generic contact detection framework for cylindrical particles in discrete element modelling[J].Computer Methods in Applied Mechanics and Engineering，2017，315：632-651.

[17]XIE C H，SONG T，ZHAO Y Z.Discrete element modeling and simulation of non-spherical particles using polyhedrons and super-ellipsoids[J].Powder Technology，2020，368：253-267.

[18]WEERASEKARA N S，LIU L X，POWELL M S.Estimating energy in grinding using DEM modelling[J].Minerals Engineering，2016，85：23-33.

[19]ZHANG B，REGUEIRO R，DRUCKREY A，et al.Construction of poly-ellipsoidal grain shapes from SMT imaging on sand，and the development of a new DEM contact detection algorithm[J].Engineering Computations，2018，35（2）：733-771.

[20]METZGER M J，GLASSER B J.Simulation of the breakage of bonded agglomerates in a ball mill[J].Powder Technology，2013，237：286-302.

[21]BIBAK Z，BANISI S.A combined physical and DEM modelling approach to investigate particle shape effects on load movement in tumbling mills[J].Advanced Powder Technology，2021，32（2）：916-930.

[22]李同清.球磨機矿磨介质动力学行为研究[D].徐州：中国矿业大学，2018.

[23]王肖江.基于离散元法的武山铜矿半自磨介质优化研究[D].昆明：昆明理工大学，2017.

[24]李鸿程，董为民，赵德化，等.矿物颗粒形状对球磨机磨矿性能的影响[J].中国粉体技术，2013，19（2）：70-73.

[25]GOVENDER N，RAJAMANI R，WILKE D N，et al.Effect of particle shape in grinding mills using a GPU based DEM code[J].Mine-rals Engineering，2018，129：71-84.

[26]CLEARY P W，MORRISON R D.Comminution mechanisms，particle shape evolution and collision energy partitioning in tumbling mills[J].Minerals Engineering，2016，86：75-95.

[27]BARRIOS G K P，DE CARVALHO R M，KWADE A，et al.Contact parameter estimation for DEM simulation of iron ore pellet handling[J].Powder Technology，2013，248：84-93.

[28]IWASAKI T，YABUUCHI T，NAKAGAWA H，et al.Scale-up methodology for tumbling ball mill based on impact energy of grinding balls using discrete element analysis[J].Advanced Powder Technology，2010，21（6）：623-629.

[29]POWELL M S，WEERASEKARA N S，COLE S，et al.DEM modelling of liner evolution and its influence on grinding rate in ball mills[J].Minerals Engineering，2011，24（3/4）：341-351.

[30]ALDER B J，WAINWRIGHT T E.Phase transition for a hard sphere system[J].Journal of Chemical Physics，1957，27（5）：1 208-1 209.

[31]CUNDALL P A，STRACK O D L.A discrete numerical mode for granular assemblies[J].Géotechnique，1979，29（1）：47-65.

[32]JOHNSON K L，KENDALL K，ROBERTS A D.Surface energy and the contact of elastic solids[J].Proceedings of the Royal Society A：mathematical，Physical and Engineering Sciences，1971，324：301-313.

[33]WALTON O R，BRAUN R L.Viscosity，granular-temperature，and stress calculations for shearing assemblies of inelastic，frictional disks[J].Journal of Rheology，1998，30（5）：949-980.

[34]MINDLIN R D，DERESIEWICZ H.Elastic spheres in contact under varying oblique forces[J].Journal of Applied Mechanics—Transactions of the Asme，1953，20（3）：327-344.

[35]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.

[36]ODA M，KONISHI J，NEMAT-NASSER S.Experimental micromechanical evaluation of the strength of granular materials：effects of particle rolling[J].Mechanics of Materials，1982，1（4）：269-283.

[37]FENG Y T.An energy-conserving contact theory for discrete element modelling of arbitrarily shaped particles：basic framework and general contact model[J].Computer Methods in Applied Mecha-nics and Engineering，2021，373：113454.

[38]FENG Y T.An energy-conserving contact theory for discrete element modelling of arbitrarily shaped particles：contact bolume based model and computational issues[J].Computer Methods in Applied Mechanics and Engineering，2021，373：113493.

[39]李昆塬，童昕，李占福，等.基于DEM-FEM耦合的球磨机梯形衬板参数研究[J].福建工程学院学报，2017，15（6）：523-527.

[40]JONSN P，PLSSON B I，TANO K，et al.Prediction of mill structure behaviour in a tumbling mill[J].Minerals Engineering，2010，24（3/4）：236-244.

[41]BEINERT S，FRAGNIRE G，SCHILDE C，et al.Analysis and modelling of bead contacts in wet-operating stirred media and planetary ball mills with CFD-DEM simulations[J].Chemical Engineering Science，2015，134：648-662.

[42]MAYANK K，MALAHE M，GOVENDER I，et al.Coupled DEM-CFD model to predict the tumbling mill dynamics[J].Procedia IUTAM，2015，15：139-149.

[43]SINNOTT M D，CLEARY P W，MORRISON R D.Combined DEM and SPH simulation of overflow ball mill discharge and trommel Flow[J].Minerals Engineering，2017，108：93-108.

[44]JONSN P，STENER J F，PLSSON B I，et al.Validation of a model for physical interactions between pulp，charge and mill structure in tumbling mills[J].Minerals Engineering，2015，73：77-84.

[45]KANO J，MIO H，SAITO F.Correlation of grinding rate of gibbsite with impact energy of balls[J].AIChE Journal，2000，46（8）：1 694-1 697.

[46]CLEARY P W，MORRISON R D.Understanding fine ore breakage in a laboratory scale ball mill using DEM[J].Minerals Engineering，2011，24（3/4）：352-366.

[47]穆海芳，何康，韓君，等.基于离散元法的球磨机介质运动的仿真研究[J].济宁学院学报，2019，40（2）：12-16.

[48]HUANG P，DING Y F，WU L，et al.A novel approach of evaluating crushing energy in ball mills using regional total energy[J].Powder Technology，2019，355（6）：289-299.

[49]BOND F C.The third theory of comminution[J].Transactions of AIME Mining Engineering，1952，193：484-494.

[50]SHI F N，XIE W G.A specific energy-based size reduction model for batch grinding ball mill[J].Minerals Engineering，2015，70：130-140.

[51]OROZCO L F，NGUYEN D H，DELENNE J Y，et al.Discrete-element simulations of comminution in rotating drums：effects of grinding media[J].Powder Technology，2020，362：157-167.

[52]CLEARY P W.Charge behaviour and power consumption in ball mills：sensitivity to mill operating conditions，liner geometry and charge composition[J].International Journal Mineral Processing，2001，63（2）：79-114.

[53]孙鹏.基于离散元素法的岩土球磨机仿真与分析[J].软件，2019，40（8）：145-148，189.

[54]刘刚，张冠宇，李婧锡，等.基于离散元的球磨机工作参数优化[J].矿山机械，2014，42（7）：79-82.

[55]PEDRAVES F，NORNIELLA J G，MELERO M G，et al.Frequency domain characterization of torque in tumbling ball mills using DEM modelling：application to filling level monitoring[J].Powder Technology，2018，323：433-444.

[56]DE CARVALHO R M，TAVARES L M.Predicting the effect of operating and design variables on breakage rates using the mechanistic ball mill model[J].Minerals Engineering，2013，43/44：91-101.

[57]HERBST J A，FUERSTENAU D W.Scale-up procedure for continuous grinding mill design using population balance models[J].International Journal of Mineral Processing，1980，7（1）：1-31.

[58]ROSALES-MARḮN G，ANDRADE J，ALVARADO G，et al.Study of lifter wear and breakage rates for different lifter geometries in tumbling mill：experimental and simulation analysis using population balance model[J].Minerals Engineering，2019，141：105857.

[59]CHIMWANI N，BWALYA M M.Using DEM to investigate how shell liner can induce ball segregation in a ball mill[J].Minerals Engineering，2020，151：106311.

[60]WU M，WANG V.Modeling ball impact on the wet mill liners and its application in predicting mill magnetic liner performance[J].Minerals Engineering，2014，61：126-132.

[61]張学东，董为民，周海燕，等.具有仿生特征的球磨机提升条耐磨性数值模拟[J].有色金属（选矿部分），2017（6）：56-62.

[62]BOEMER D，PONTHOT J P.A generic wear prediction procedure based on the discrete element method for ball mill liners in the cement industry[J].Minerals Engineering，2017，109：55-79.

[63]REZAEIZADEH M，FOOLADI M，POWELL M S，et al.A new predictive model of lifter bar wear in mills[J].Minerals Engineering，2010，23（15）：1 174-1 181.

作者简介：裴英杰（1997—），男，山东平度人，硕士研究生，研究方向为碎磨理论与工艺;昆明市五华区一二一大街文昌路68号，昆明理工大学国土资源工程学院，650093;E-mail：z179180802@163.com

通信作者，E-mail：13515877@qq.com，13888532894

裴英杰，肖庆飞，石贵明，周强（1.昆明理工大学国土资源工程学院; 2.矿冶过程自动控制技术国家（北京市）重点实验室;3.玉溪师范学院化学生物与环境学院

Research progress and application of discrete element modeling in ball mill simulationPei Yingjie，Xiao Qingfei，Shi Guiming，Zhou Qiang

（1.School of Land Resource Engineering，Kunming University of Science and Technology;

2.State（Beijing） Key Laboratory of Process Automation in Mining & Metallurgy Process;

3.School of Chemical Biology and Environment，Yuxi Normal University）

Abstract：As important equipment in the ore grinding and crushing stage，the ball mill plays a decisive role in improving efficiency and reducing consumption of the concentrator.Ball mill power consumption and energy consumption distribution，the design optimization of technical parameters and liner structural parameter are important content and technical difficulties that affect the grinding effect of ball mills.Discrete element modeling can accurately simulate the deformation process of particle motion and record the power and energy distribution of particles，so that in the field of ball mill simulation，it has good application prospects.The paper introduces the characteristics of the widely used particle model，contact model，and joint simulation method of the discrete element modeling in the ball mill simulation process and their research progress，as well as  the application of ball mill power and energy distribution research，technical parameter optimization，liner structural parameter optimization，summarizes and points out that the refinement of the ball mill simulation model，the wide application of joint simulation，and the intelligent adjustment of ball mill technical parameters are the future ball mill simulation research development trend.

Keywords：discrete element modeling;ball mill;simulation;particle model;contact model;coupling method