基于DEM-FEM耦合的球磨机梯形衬板参数研究

李昆塬， 童昕， 李占福， 叶一青

(福建工程学院 机械与汽车工程学院， 福建 福州 350118)

球磨机是一种广泛应用于建材、化工、机械、选矿、机械等部门的磨料设备，衬板是球磨机中极其重要的一部分，又以梯形衬板使用较为广泛[1]。球磨机在工作过程中，衬板不仅受到磨料和钢球等的冲击和摩擦作用，还受到矿浆原料等化学物质的腐蚀，导致衬板磨损。选择梯形衬板的合理尺寸以改善衬板的抗磨损能力，提高衬板的强度一直是急需解决的工程问题。

M.Rezaeizadeh[2]利用离散元法研究了衬板数量、宽度、高度等参数对衬板磨损量的影响；M.Yahyaei[3-4]对衬板的轴向方向的磨损进行研究，并对其进行改进。N.Djordjevic[5-6]研究了衬板高

度、提升条数目、衬板与介质间的摩擦系数对磨矿效率的影响；P.Jonsén[7-8]利用DEM 与FEM 相互耦合为衬板的结构设计提供了新途径。本文以大型球磨机为例，利用DEM 与FEM相互耦合的方法进行仿真，对球磨机衬板的磨损进行研究，得出衬板的应力分布、变形程度，通过正交分析耦合结果来获得因素的显著性及衬板结构的优方案，进而提高球磨机的寿命和经济效益。

1 试验方法介绍

离散单元法(discrete element method,DEM)是在岩土工程领域里发展起来的一种模拟岩土块体，颗粒群力学过程的数值方法。试验使用的接触模型为Hertz-Mindlin(no slip) built-in[9]，颗粒在相互接触时会受到两种力，分别为法向作用力和切向作用力，法向力的依据是Hertz理论，切向力的依据是Mindlin-Deresiewicz[10]。有限元法(finite element)是将分析的机构或实体划分为有限个微小单元体的离散化方法，相邻两单元之间只通过节点相连接。

DEM-FEM的耦合分2种[11]：双向耦合和单向耦合。本文采用单向耦合，即把离散元中的分析结果作为一种载荷条件输入到有限元中，然后对设备进行静力学或动力学分析，不考虑设备的变形对颗粒运动的影响。具体步骤如图1所示。

图1 DEM-FEM耦合流程Fig.1 DEM-FEM coupling process

2 试验过程

2.1 建立正交试验

试验对象是有效直径为4 300 mm的某大型球磨机，其衬板厚度为100 mm，根据经验与实践[12]选择衬板的倾角、顶边长度、高度作为正交试验的因素，根据其结构尺寸范围设定3个水平，建立正交试验，正交试验水平因素表如表1所示，正交表为L9(34)的表格，如表2所示。

2.2 EDEM模拟仿真试验

建立球磨机的三维模型，为节省仿真时间，简化三维模型(采用周期性边界模型，即在有限空间内，颗粒从任何方向离开仿真区域将立即从相反方向重新进入)，只选取球磨机衬板中的一段(500 mm)，三维模型建好后保存为stp格式，导入EDEM软件中。球磨机内梯形衬板的材质为铸钢，研磨介质为钢球，磨料为矿石，颗粒与颗粒，颗粒与衬板，颗粒与钢球间的物理属性见表3。

表1 衬板正交试验因素水平表

表2 衬板L9(34)正交试验表

为了更好地模拟现实工况，试验选用粒径分别为80、100、200、300 mm，数量分别为15 800、2 700、255、25的矿石和一种粒径为120 mm、数量为6 220的钢球进行仿真。颗粒总数为25 000个，颗粒总体积约为球磨机内体积的35%，即球磨机的填充率约为35%，其中钢球约占15%，矿石约占20%，总质量为66 027 kg。该球磨机临界转速为14.45 r/min,转速率为74.46 r/min，则其工作转速为10.76 r/min，即1.126 8 rad/s,启动时间为1 s，落料时间为0.7 s，运行时间为9 s。至此，球磨机建模完毕，进行试验，最后导出Pressure文件。

2.3 ANSYS Workbench求解

在Workbench中建立DEM Solutions模块，将Pressure.axdt文件导入到Results中。建立静力系统模块，将Results与Setup连接起来，形成数据共享。导入球磨机衬板几何模型，进入Mechanical界面，定义球磨机衬板的材料、密度、泊松比；划分网格，考虑到网格数量和质量，将网格尺寸定为50 mm，网格节点数量为118 389个，网格单元数为22 710个。添加约束，载入Pressure，建立衬板的有限元模型，如图2。

2.4 有限元求解处理

在对衬板的磨损正交试验中，选择变形、应力、应变3个指标来判定衬板的磨损情况，共9组试验，每组重复3次，求解平均值。将试验结果记录在正交表中，记录结果如表4所示。

图2 衬板有限元模型Fig.2 Finite element model of the liner

颗粒模型密度/(kg·m-3)泊松比剪切模量/GPa滚动摩擦系数恢复系数静摩擦系数衬板78000．25700．50．20．01钢球78000．25700．50．20．01矿石31400．300．10．50．20．01

表4 正交试验结果Tab.4 Orthogonal test results

3 分析结果

3.1 正交试验-直观分析

通过正交试验，将得出的结果进行分析，求其极差。根据极差R的大小来分析因素主次，比较算术平均值k1、k2、k3，确定其优方案。表5为变形的直观分析。从表中可以看出，对于总变形的影响因素主次为衬板的倾角>衬板顶边长度>衬板高度。从他们的极差值中可以看出，优方案的衬板倾角为22°，衬板顶边长度为130 mm，衬板高度为225 mm。同样地，由表6、7中可以得出，等效应力的影响因素主次为衬板高度>衬板的倾角>衬板顶边长度，优方案的衬板高度为200 mm，衬板的倾角为22°，衬板顶边长度为150 mm；等效应变的影响因素主次为衬板高度>衬板的倾角>衬板顶边长度，优方案为衬板高度为200 mm，衬板的倾角为22°，衬板顶边长度为150 mm。

表5 各因素的变形极差表

注：极差R=最大k值-最小k值

表6 各因素的等效应力极差表Tab.6 Equivalent stress range table

注：极差R=最大k值-最小k值

当筒体的半径和转速一定时，倾角的大小是确定的，即仿真速度为1.126 8 rad/s时，其倾角应为22°。顶边长度决定了相邻物料间的最小距离，过小的距离虽然可以使物料产生的冲击力集中，但加剧了衬板磨损。本次试验顶边长度为130、150 mm都为优方案的尺寸，对衬板有明显的磨损作用，在设计衬板的顶边长度时，应选择更大的尺寸，从本次试验来看，170 mm是最合适的。在衬板的高度对物料有一定的提升能力情况下，提高衬板的高度可以提升衬板的寿命，本次试验中，200、225 mm都为优方案，即会产生较大的变形量或者较大的应变或应力。因此，衬板高度为250 mm会有效缓解衬板的磨损，提高衬板的寿命。

表7、 8可看出，空列的极差大于任何一列的极差，从而可以推断存在不可忽略的交互作用，或者漏掉了对试验结果有重要性的其他影响因素。本文仅对衬板的结构尺寸进行研究，而现实中，球磨机的工作环境及其复杂，存在众多因素会对球磨机的衬板磨损产生影响，因而分析出来的结果是可信的。

表7 各因素的等效应变极差表Tab.7 Equivalent strain range table

注：极差R=最大k值-最小k值

表8 变形方差分析Tab.8 Variance analysis of deformation

注：F0.05(2,2)=19，F0.01(2,2)=99[12]

3.2 正交试验-方差分析

正交试验中，需要对试验结果进行方差分析，来弥补直观分析法的不足。表8是对最大变形的方差分析，通过对F值大小的比较，来判断因素的显著性。当F0.05

表9 等效应变方差分析

注：F0.05(2,8)=4.46，F0.01(2,8)=8.65

表10 等效应力方差分析Tab.10 Variance analysis of equivalent stress

注：F0.05(2,8)=4.46，F0.01(2,8)=8.65

4 结论

本次试验以衬板的磨损状况为目的进行研究，衬板的磨损量越小越好。对正交分析得出的结论做进一步分析，所得最优方案为衬板倾角22°，顶边长度170 mm，衬板高度250 mm。

运用有限元与离散元的耦合方法对衬板的结构尺寸进行研究，通过正交分析将变形、应力、应变作为衬板磨损的判定指标，发现衬板的倾角对

衬板的磨损具有显著的影响。因而在改善衬板磨损的情况时，应对衬板的倾角做进一步优化。而衬板的顶边长度和衬板高度虽然分析结果表明影响不显著，但同样存在优方案，可加以考虑。

球磨机的实际工况极其复杂，存在众多影响因素，如球磨机转速，球磨机内的化学环境等。本文只对衬板的倾角、顶边长度、高度进行研究，旨在提供一种可供参考的研究方案。

[1] 朱祥，王国强，邵明亮.球磨机梯形衬板设计[J].矿山机械,1993(5)：5-7.

[2] 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):1174-1181.

[3] Yahyaei M,Banisi S, Hadizadeh M. Modification of SAG mill liner shape based on 3-D liner wear profile measurements[J]. International Journal of Mineral Processing,2009,91(3):111-115.

[4] Yahyaei M, Banisi S. Spreadsheet-based modeling of liner wear impact on charge motion in tumbling mills[J]. Minerals Engineering,2010,23(15):1213-1219.

[5] Djordjevic N. Discrete element modelling of the influence of lifters on power draw of tumbling mills[J]. Minerals engineering,2003,16(4):331-336.

[6] Djordjevic N. Influnce of charge size distribution on net-power draw of tumbling mill based on DEM modeling[J]. Minerals Engineering.2005,18(3):375-378.

[7] Jonsén P, Pålsson B I, Tano K, et al. Prediction of mill structure behaviour in a tumbling mill[J]. Minerals Engineering,2011,24(3):236-244.

[8] Jonsén P, Pålsson B I, Häggblad H Å. A novel method for full-body modelling of grinding charges in tumbling mills[J]. Minerals Engineering,2012,33(6): 2-12

[9] Mindlin R D. Compliance of Elastic bodies in Contact[J]. Journal of Applied Mechanics，1949，16：259-268.

[10] 孙其城,王光谦.颗粒物质力学导论[M].北京：科学出版社,2009.

[11] 侯亚娟.基于离散元与有限元耦合的大型球磨机衬板性能研究[D].长春:吉林大学,2015.

[12] 李云雁,胡传荣.试验设计与数据分析[M].北京：化学工业出版社,2008.