离散元及其耦合方法在球磨机研究中的应用

朱洪民 王泽红 孙维波 齐晓楠

（东北大学资源与土木工程学院，辽宁沈阳110819）

球磨机作为选矿厂生产运行的关键设备之一，其能耗（电耗、材料成本）一般占全厂能耗的30%～60%［1］。球磨机具有结构简单、性能稳定、应用广泛的优点，同时具有能耗、钢耗较大且能源利用率低的缺点。国内外专家学者针对球磨机进行了大量实验研究工作，并取得了很多成果。近些年来，随着计算机性能的不断提高和一些新的数值分析方法的提出，数值分析研究在球磨机机理与参数优化中的应用得到了快速发展，并取得了一定的研究成果。

在球磨机机理与参数优化等研究中常用的数值模拟方法有离散元法（DEM）和离散元耦合方法。常见的离散元耦合方法有离散元与有限元耦合法（DEM-FEM）、离散元与光滑粒子流体力学耦合法（DEM-SPH）、离散元与计算流体力学耦合法（DEMCFD）以及离散元、光滑粒子流体力学和有限元耦合（DEM-SPH-FEM）等。数值方法的耦合又分为单向耦合方法和双向耦合方法。单向耦合方法只考虑一种（或几种）方法研究的对象对另一种（或几种）方法研究对象的影响，忽略反向影响；而双向耦合考虑多种方法研究的对象在设定步时内的交互影响。

1 常用数值分析方法简介

1.1 DEM

离散元法（Distinct Element Method，DEM）是由Peter Cundall于1971年提出的一种针对复杂非连续系统动力学问题的新型数值分析方法，其基本思想是把不连续体分离为刚性元素的集合，根据牛顿第二定律从刚性元素间接触力入手找出其接触的本构关系，对非连续、离散的刚性元素建立运动方程［2］。东北大学王泳嘉于1986年在第一届全国岩石力学数值计算及模型试验讨论会上首次引进介绍了该分析方法，随后在我国迅速发展起来［2-3］。Mishra于1992年首次在球磨机中应用二维DEM方法研究了衬板对钢球运动的影响并预测磨机功耗，预测结果与实验验证相吻合［4］。此后随着计算机计算能力的提高，三维DEM方法具有了可行性，推动了离散元法在球磨机介质运动规律及状态研究中的应用。

1.2 DEM-FEM

有限元法（Finite Element Method，FEM）于1943首次出现在Courant.R的论文中，在随后兴起的应用数学、现代力学及计算机科学相互渗透下，于1960年正式提出，其理论现今已逐渐完善成熟［5-6］。FEM的基本思想是将待研究连续介质对象的连续求解区域离散为若干有限个以一定方式相互连接的单元（小区域）。在节点上引进等效力以代替实际作用于单元上的外力，然后对单元进行力学分析，再整体分析，以化整为零、集零为整的思路使问题简化为适合数值解法的结构型问题。在研究球磨机时，该方法利用DEM分析球磨机内钢球和物料的各种运动学特性等信息，利用FEM分析球磨机衬板的应力和应变、轴承扭矩等信息。DEM-FEM的耦合类型分为2种［7］：单向耦合与双向耦合。球磨机研究中，单向耦合是指把DEM分析结果作为一种载荷条件输入到FEM分析中，通过FEM对设备进行力学分析，只考虑颗粒运动特性对设备的影响，不考虑设备的形变、应变等参数变化对颗粒运动的影响；双向耦合可以实现颗粒的运动特性与设备的变形等参数在一定计算步时内以载荷条件互相迭代，考虑了球磨机对颗粒运动影响的同时，也考虑了颗粒运动对球磨机的影响。在DEM-FEM耦合中，目前单向耦合较为成熟，双向耦合还需进一步深入研究［7］。在球磨机衬板优化设计与磨蚀预测研究中，衬板的变形量较小，对颗粒的影响不大，单向耦合同样取得了理想的效果。

1.3 DEM-SPH

光滑粒子流体力学法（Smooth Particle Hydrodynamics Method，SPH）由Gingold和Lucy等先后于1977年提出。该方法属拉格朗日型的一种无网格法，是最简单、最高效的无网格法之一，在与流体力学相关的领域中得到了很好的发展应用［6］。SPH方法［8，9］的基本思想是将整个连续介质流体离散化，用一系列任意分布的光滑粒子，在待求点的支持域内对物理量进行核近似，然后对所求解的微分方程进行粒子近似，进而通过求解解决问题。在球磨机研究中，该方法指的是用DEM模拟物料和钢球，用SPH研究矿浆；其单向耦合指的是只考虑矿浆对料球的影响，不考虑料球对矿浆的影响，或者只考虑料球对矿浆的影响而不考虑矿浆对料球的影响，而双向耦合考虑矿浆、料球在一定步时内的相互影响。运用DEMSPH耦合法研究球磨机时考虑到了湿磨工艺复杂矿浆条件对磨矿过程的影响，更加接近球磨机湿磨的工作状态，受到了许多研究工作者的关注。Cleary P.W.等［10］首次运用DEM-SPH单向耦合法模拟预测了磨矿过程矿浆体积、矿浆黏度、操作参数对矿浆的运动及分布的影响。近些年，计算机性能的迅速提升降低了数值模拟的计算成本，使球磨机应用三维DEM-SPH模型模拟计算具有了可行性，相关研究得到了一定发展但仍需不断深入。

1.4 DEM-CFD

现在国际上一般认为，英国气象学家Richardson L.F.于1910～1917年尝试用有限差分法迭代求解Laplace方程的方法来计算圆柱绕流和大气流动从而预报天气的工作标志着计算流体力学（CFD）的诞生［11］。CFD数值求解的基本思想是把原来在空间与时间坐标中连续的物理量的场（如浓度场、速度场等），用一系列有限个离散点（称为节点，node）上的值的集合来代替，通过一定的原则建立起这些离散点上变量值之间关系的代数方程（称为离散方程），求解所建立起来的代数方程以获得所求解变量的近似值［12］。近40多年来CFD计算方法发展特别迅速，以数值求解Euler方程和 RANS（Reynolds Averaged Navier-Stokes）方程为代表的CFD技术已经广泛应用于工程领域［13］。在球磨机研究中，DEM-CFD是指用DEM研究介质物料，用CFD研究矿浆，并将这2种方法在一定步时内建立力学关系从而实现矿浆与物料的交互影响。近几年来，随着CFD技术在选矿工程中的广泛应用，应对复杂流固体环境的DEM-CFD耦合模型进一步成熟，一些专家学者将DEM-CFD耦合技术应用到了湿式磨矿领域，研究矿浆对磨矿效果的影响，并取得了一些成就，相关研究正在不断深入。

1.5 DEM-SPH-FEM

在球磨机研究中，DEM-SPH-FEM指的是用DEM模拟球磨机中物料介质，用SPH模拟球磨机中矿浆，用FEM模拟球磨机腔体衬板，并将这3种方法通过力学关系耦合在一起，在一定步时内实现各个方法研究对象的交互影响或者单向影响。Jonsén P.等［14-15］在球磨机研究中使用了DEM-SPH-FEM的方法，并对模型、模拟结果进行了验证，取得了较理想的结果，相关研究刚刚开始仍需进一步深入。

2 单一离散元法在球磨机中的应用

单一离散元法相对于离散元耦合方法计算成本低、模型稳定、容易实现，在球磨机研究中应用相对更广泛。球磨机中单一离散元法的应用主要在介质运动状态及机理研究、磨机参数优化研究、球磨机功耗及能量分布等方面。

2.1 介质运动状态及机理研究

林阳辉等［16］采用DEM研究了转速及衬板磨损对球磨机介质运动轨迹的影响，发现衬板磨损后导致介质运动轨迹发生变化，降低了球磨机的磨矿效率，适当地提高转速可使球磨机维持衬板磨损前的磨矿效率。崔泽群［17］利用DEM简化球磨机为圆柱形水平滚筒并研究了其中二元颗粒轴向分层现象，发现了一种颗粒轴向分层的新模式，并解释了轴向分层的原因。董铭峰［18］、卢建坤［19］、李鸿程［20］分别利用DEM对球磨机介质运动状态进行了数值分析，并对球磨机转速率、充填率、料球比、衬板形状等参数进行了优化。夏恩品［21］、王卓［22］分别用DEM研究球磨机介质直径对磨机其他参数和物料运动状态的影响，并给出了介质直径及配比的选择方法。Cleary P.W.等［23］利用基于DEM的新计算模型研究了磨机的粉碎机制、颗粒形状演化和碰撞能量划分。Bian X.等［24］采用DEM对球磨机进行了数值模拟，并对球磨机介质的运动行为进行了定性和定量的验证，结果表明，衬板和挡板的扭矩受提升条数量、高度和磨机转速率的影响，解释了球磨机扭矩和功耗随提升条数量和高度及介质配比的变化规律。Cleary P.W.［25］基于DEM将微尺度模型嵌入到宏观尺度模型得到宽尺度模型，从而使计算中包含细颗粒变得可行。通过微尺度模型中粉体运动状态预测并证明了磨削作用来自相邻介质层之间的薄层粉体之间的相互挤压，描述了计算粉体效应有用功及碰撞能的方法，并说明了高粉体负载对介质流动行为缓冲作用的重要性。Weerasekara N.S.等［26］利用DEM模拟探究了球磨机大小和介质粒度分布对球磨机磨碎环境的影响，从而揭示磨碎过程及该过程能量消耗的本质，结果证明，脚趾冲击区域具有较高的冲击能量，大部分剪切区域有较大的切向能量，大部分的能量都是由中等大小的颗粒消耗的。Jiang S.等［27］利用DEM提出了颗粒重力加速度、磨机直径和旋转速度缩放比例关系，并在不同直径磨机上对物料运动进行了一系列缩放模拟，发现不同直径和转速的球磨机中颗粒的运动状态大体一致，从而定性验证了提出的放缩关系，结果证明，在有庞大数量颗粒的情况下利用相似原理模拟降低计算成本是一种有效的方法。Francisco P.等［28］利用DEM对小型球磨机进行了分析模拟，并与实验结果进行对比，得出可以通过球磨机负载扭矩信息准确表征球磨机的充填率，利用扭矩的频率域实现对球磨机充填率的监视。

2.2 磨机参数优化研究

陈友川等［29］用DEM对ϕ1 m×3 m球磨机内介质的运动状态进行了数值模拟，分别以短棒形和立方体等非球形物料为介质，分析了介质类型、介质配比、衬板高度、衬板数量、充填率等因素对介质运动状态的影响，为球磨机的参数优化设计提供了参考。Cleary P.W.［30］、Hlungwani O.等［31］利用二维 DEM分析研究了球磨机的充填率、转速、衬板形状和数量以及介质和物料性质（球配比、介质和物料的形状、钢球类型、物料粒度分布）对磨机内物料颗粒的运动情况和磨机功率产生的影响。Rezaeizadeh M.等［32］根据影响衬板磨损的主要参数建立了一种新的DEM预测模型，并利用该模型根据矿石类型、矿石与衬板之间的相对速度、衬板摩擦系数、衬板硬度来预测衬板磨损率，预测结果与伊朗Sarcheshmeh铜矿公司的磨机实测数据基本吻合。Powell M.S.等［33］利用DEM分析预测了衬板磨蚀过程，并阐述了衬板磨蚀后对磨矿效果的影响，提出衬板轮廓预测和磨削速率相结合是一种平衡衬板寿命和磨矿性能的方法，该方法将会成为衬板设计的重要工具。李臣等［34］利用DEM对设计的椭圆腔形结构球磨机与传统的圆筒腔形球磨机从转速差、粉磨效率等几个角度进行了数值模拟对比分析，得出研究的椭圆腔形结构球磨机较圆筒腔形球磨机磨矿效果有较大的提升。李笑同［35］在筒式球磨机中增加旋转导向板结构，并用DEM对结构改变后的球磨机进行仿真模拟，将振动特性和力学特性的变化同原球磨机进行了对比分析。边小雷［36］利用DEM研究了衬板提升条的高度和数量2个参数对球磨机转矩和球磨机功率的影响规律。王继生等［37］以直径6 m的大型球磨机作为研究对象，用DEM分别模拟分析了矩形衬板、梯形衬板和波形衬板对球磨机磨矿效率的影响。赵魏等［38］利用DEM模拟计算球磨机介质和矿料的相互运动状况，通过分析介质运动轨迹来判断磨机衬板及提升条设计的合理性。赵选恒等［39］运用DEM对典型的几种锥面分级衬板的分级作用进行数值模拟，表明双斜度分级衬板效果最佳，可以为衬板的选择提供理论依据。Xu L.等［40］利用DEM建立了剪切冲击能量模型，模拟分析与实验相结合研究了磨机转速对筒壁衬板的磨损速度和磨损分布的影响，最后根据能量利用和衬板磨损情况对转速和提升条形状进行了综合评估优化。

2.3 球磨机功耗及能量分布

Cleary P.W.运用三维DEM模拟计算了磨机的有用功率［14］。许利学等［41］通过DEM仿真得到球磨机碰撞能分布结果，利用该结果对破碎率计算公式中的相关参数进行求解计算，得到的破碎率函数与粉磨实验值基本吻合。史国军［42］利用DEM分析了球磨机介质尺寸、磨机转速率以及介质充填率等介质工作参数对球磨机介质运动规律与有用功率的影响。田秋娟等［43］用三维DEM研究ϕ5.5 m×8.5 m球磨机的衬板倾角、衬板高度和衬板数量等设计参数对球磨机功率与冲击能量分布的影响。郝万军等［44］用三维DEM研究大型球磨机的球料比对其工作性能的影响，从功率、比功率、碰撞能量谱、介质对衬板的冲击磨损消耗的能量占总消耗能量的百分比等方面揭示了球料比对球磨机工作效率的影响，并将其作为球磨机优化设计的参考依据。Cleary P.W.［45］利用DEM对直径4 m的球磨机的介质运动、磨矿效率进行了数值分析，探讨了钢球配比对磨矿效率的影响，模拟预测的功率与球磨机的额定功率基本一致。Wang M.H.等［46］将碰撞能、耗散能量和最大冲击能量与物料破碎联系在一起，利用DEM将能量信息与群体平衡模型相结合，预测了产品粒度随介质磨削时间的演变过程，并与实验数据进行了比较，结果表明，碰撞能不需要调整就能直接与颗粒破碎后的粒度分布建立联系，可以很好地预测产品的粒度。耗散能量和最大冲击能量也可以用于粒度预测，但需要选择函数仔细校准。Cleary P.W.等［47］利用DEM与实验结合研究了直径0.5 m的球磨机的磨碎过程，模拟分析获得的能量谱表明，对于典型的矿石，大部分碰撞不会对颗粒产生有效破碎，最好的操作条件下能量利用效率也不到10%。

3 离散元及其耦合方法在球磨机中的应用

球磨机研究中常见的离散元耦合方法有DEMFEM、DEM-SPH、DEM-CFD、DEM-SPH-FEM。多种方法耦合能实现各种方法优势特长的结合，能更加真实全面地模拟球磨机实际工作情况，相关研究具有极大的实际意义，已有很多开拓性的研究成果。但是由于耦合模型复杂不稳定、技术经验要求极高以及计算成本的剧增，相对于DEM，磨机中耦合方法的模拟仿真进展相对缓慢，仍需进行深入研究。

3.1 DEM-FEM

Jonsén P.等［48］建立了二维DEM-FEM双向耦合模型，分析了球磨机运转中衬板的应力和挠度变化，为衬板的结构优化设计提供了新方向。侯亚娟［7］利用DEM-FEM单向耦合方法对球磨机运转过程进行了模拟仿真，获得了球磨机衬板的磨损分布规律与球磨机运行时衬板的受力变化之间的关系。孙姗姗等［49］联合EDEM软件与ANSYS软件进行DEM-FEM单向耦合，研究了波形衬板的等效应力和分布，并根据分析结果设计出了新型波形衬板。汪滋润等［50］利用DEM-FEM单向耦合研究分析了球磨机在不同转速工作条件下，衬板应力的变化规律，为后续结构设计提供了数据支持。李昆塬［51］利用DEM-FEM单向耦合对直径为8.6 m的大型球磨机的梯形衬板进行了数值模拟，分析了衬板倾角、顶边长度和高度3个因素对球磨机梯形衬板变形、应变和应力的影响规律，并获得了最优方案。

3.2 DEM-SPH

澳大利亚学者 Cleary P.W.和 Morrison R.D.［52］采用DEM-SPH单向耦合对ϕ1.8 m×0.6 m半自磨机进行了模拟仿真，使用DEM对矿物颗粒进行建模，通过SPH对矿浆进行建模，模拟结果反映了磨机内矿浆流动的详细信息，包括了不同区域的预测、矿浆室的充填率、矿浆流的最高点、物料的运动负荷；对粉碎室内物料流量对磨矿效率的影响进行模拟，分析了侧壁在产生复杂介质三维再循环模式方面的重要性及其对轴流的影响，揭示了侧壁提升条在促进矿浆运输流动方面的积极影响。Cleary P.W.［53］基于DEMSPH双向耦合法提出并演示了一种可用于复杂工业过程的密集矿浆与粗颗粒物相互作用的模型，并将该模型运用在球磨机湿磨中，与单向耦合模型结果进行了对比，分析总结了改变耦合力大小对矿物颗粒和矿浆负载的影响以及改变充填率对介质运动状态和介质分布的影响。Sinnott M.D.等［54］利用DEMSPH耦合法研究了溢流型球磨机排料和滚筒筛的矿浆流流动对磨机性能的影响，建模中粗粒矿石和球磨机介质模型使用了DEM，矿浆模型使用了SPH，DEM-SPH双向耦合模拟分析了磨机内复杂固液环境中的物料运动信息，并揭示了滚筒筛运输分级性能和滚筒筛挡板几何形状之间的联系。研究发现，在磨粉室里由于矿浆的加入，水拖拽力可以降低肩部区域和脚趾区域的大小，并在脚趾区域上方形成一个大的“水洼”。滚筒筛中有挡板，可以提高筛选性能。使用倾斜度较大的螺旋挡板会减少料浆在筒体内的滞留时间，既增加了矿浆和球的垂直提升，又增加了滚筒筛筛选的有效面积，从而使物料能及时排除而提高磨矿效率。

3.3 DEM-CFD

Mayank K.等［55］利用DEM方法模拟固体颗粒，用CFD方法模拟连续流体，提出了一种基于离散元与CFD单向耦合的三相计算模型。用C++程序进行建模，通过CFD网格得到介质的体积和速度，并对矿浆对介质的拽力作用进行了分析，将与模拟规模大小相等的实验球磨机等效PEPT（正电子辐射颗粒跟踪技术）实验结果分别与DEM-CFD单向耦合、单一CFD模拟、单一DEM模拟结果比较，结果发现，DEMCFD耦合模型具有单一DEM和单一CFD均不能实现的功能，能够有效地模拟计算介质和矿浆的自由表面轮廓，同时还能计算球磨机质量剖面的轴向中心。

3.4 DEM-SPH-FEM

Jonsén P.等［14，15］提出了一种新的研究球磨机湿磨的耦合模型，利用DEM模拟物料颗粒、SPH模拟矿浆、FEM模拟球磨机衬板，并用3种方法建立耦合的计算模型预测矿浆流体的压力，描述了介质冲击波在球磨机中的传播规律，研究了矿浆的黏度和密度对矿浆阻尼作用和球磨机扭矩的影响，与装有精确转矩仪的小型批量球磨机的磨蚀试验结果进行了比较，验证了该模型和方法的有效性，将模拟介质的运动状况与真实介质运动状况的高速视频进行比较表明，数值计算结果与实验测量结果吻合较好。因此，利用此方法预测球磨机中介质的冲击力分布、预测磨机中高能区域大小、优化操作条件具有可行性。

4 展望

数值分析在球磨机中的广泛应用对球磨机降低能耗、降低钢耗、提高磨矿效率提供了新的思路与途径，并成为研究解决球磨机实际问题的一种便捷、低成本的方法。随着计算机性能以“摩尔定律”的速度不断迅猛发展，计算成本不断降低，计算精度不断提高，数值模拟将更加接近实际，更加广泛地被应用在球磨机领域。数值分析在球磨机中的应用虽然取得了许多进展，但由于球磨机磨矿过程影响因素多，数学模型只是实际问题的近似表达，模拟过程和结果依旧具有一定的局限性，下面几个方面仍需要继续深入研究。

（1）多种耦合计算方法的耦合研究。选矿厂中常用的湿磨是在固液环境下进行的，而衬板又是连续刚性材料，单一的离散元法不足以接近真实的环境，需要多种方法联合耦合，而多种耦合环境复杂，模型不稳定，耦合参数难以确定等问题依旧需要进一步深入研究。

（2）注重实验与模拟的有机结合，以解决大型球磨机设计选型、参数优化、装球制度等实际问题。

（3）数学模型理论和算法的研究比较少，目前的模型相比实际情况过于简单，忽略了一些较重要因素的影响作用。研究更接近球磨机实际情况、基于破碎机理的新数学模型或者基于目前算法的不断优化、二次开发等都具有极大的实际意义。

（4）集成通用软件的开发。模拟过程需要多学科的技术支持，对研究者的要求极高，集成通用模拟软件的开发，将降低模拟的门槛，使模拟仿真在球磨机设计生产中更加方便快捷地推广开来。

（5）数值仿真结果的实验验证不够充分，特别是球磨机湿磨过程中物料、矿浆、介质运动行为的实际数据难以检测，不能精确地验证实验结果，模拟中需要开发高新技术对实验结果进行验证。