粉体混合过程及搅拌功率的DEM数值模拟和实验

宣颖，刘雪东，2，周成奇，刘文明，2，刘红梅，2，顾宇彤，彭涛

（1 常州大学机械工程学院，江苏常州 213164；2 江苏省绿色过程装备重点实验室，江苏常州 213164；3 安徽瑞旭搅拌设备有限公司，安徽广德 242200）

粉体混合机是一种对多种粉体物料进行混合作业的机械设备，广泛应用于石油、化工、饲料、制药、食品、材料等领域。在粉体混合机的研究中，搅拌功率是重要的参数之一，不仅直接反映粉体混合机的能量消耗，还为搅拌轴的强度计算提供理论依据[1-2]。目前，粉体混合时的搅拌功率没有准确可靠的计算公式，在实际工程应用中，搅拌功率的确定一般通过经验估算，具有较大的误差[3-4]。当选用的电动机功率不合适时，不仅会提高生产成本，还会降低电动机运行的功率因数[5]。因此，研究粉体混合机搅拌过程及搅拌功率，可以为粉体混合设备选择合适的电动机，以最小的能量输入来达到最好的混合效果，这对降低能耗、提高工业生产效率具有重要意义[6-7]。

目前，国内外学者针对粉体搅拌功率进行了大量实验和数值模拟。实验方面，马青山等[8]使用内外螺带和叶片式两种搅拌桨探究了搅拌型式、转速、粉体物性及加料量对搅拌功率的影响规律，针对卧式单轴搅拌反应器给出了功率准数关联式；杨敏官等[9]利用固体激光发生器和扭矩仪对偏心搅拌槽内颗粒的悬浮特性进行了实验研究，探究了功率消耗和混合时间的影响。数值模拟方面，Sato 等[10]对高剪切混合机内颗粒运动进行了数值模拟，对不同转速下搅拌桨的扭矩、颗粒速度分布及颗粒受力进行分析，发现颗粒动能决定搅拌扭矩和颗粒的运动行为；Siraj等[11]通过数值模拟探究桨叶各截面上的力场分布，发现颗粒间的作用力主要集中在桨叶前端和床层底部；Remy 等[12]通过数值模拟得到颗粒速度、颗粒温度分布及颗粒的混合动力学，在微观尺度上分析了叶片式混合机中颗粒的运动；Bao等[13]采用离散元法研究了叶轮结构、转速和填充水平对装有平桨叶和斜桨叶的圆筒混合器中颗粒流动和混合的影响，结果表明，三叶混合器的混合性能和效率优于两叶和四叶混合器。

国内外学者的实验或数值模拟重点多集中在搅拌功率与某种单一因素的定性关系研究上，没有给出明确的计算公式。鉴于此，本文采用数值模拟与功率测试相结合的方法，以直叶桨式粉体混合机为研究对象，研究粉体混合机中搅拌桨的转速、直径、桨叶数目等特性参数对功率消耗和扭矩的影响，并拟合出功率计算公式。最后，搭建粉体搅拌试验台进行实验验证，对比分析模拟与实验过程中粉体混合时搅拌功率和扭矩的变化特性，验证基于数值模拟所推导公式的准确性。

1 DEM数值模拟及其接触模型

1.1 DEM方法介绍

采用离散元方法（DEM）对粉体混合机中颗粒的运动和搅拌桨的扭矩进行分析。DEM 根据牛顿第二定律计算单个粒子的位置、速度和加速度以及作用于粒子的力，计算公式如式(1)、式(2)所示。

式中，X、m、t分别表示位置矢量、质点的质量和时间；Fc、g分别表示作用于质点的接触力和重力加速度；ωp、τc和I分别表示角速度、旋转运动的扭矩和粒子的惯性矩。混合机中的每一个颗粒受到的接触力来自三个方面：其他颗粒、搅拌桨叶和容器壁。

颗粒与颗粒间、颗粒与几何体之间均采用振动方程进行模拟，采用Hertz-Mindlin 接触模型，其中，模型法向部分和切向部分的接触力分别从Hertz 理论和Mindlin 接触理论中得到[14]。应用Cundall 和Strack[15]提出的图1 所示接触模型计算接触力，该模型利用弹簧、缓冲器和摩擦滑块等机械元件表示接触力。弹簧在粒子碰撞时提供排斥力。缓冲器和摩擦滑块分别代表碰撞的能量耗散和摩擦特性。法向和切向的接触力Fcn和Fct的计算公式见式(3)～式(5)。

式中，δ、k、η和μ分别表示接触颗粒间位移、弹簧刚度、阻尼系数和摩擦系数；Vr分别表示颗粒间的相对速度。阻尼系数η[16]由式(6)计算。

式中，e表示恢复系数。

搅拌桨的扭矩T由式(7)表示为基于颗粒与搅拌桨叶碰撞的排斥力作用在旋转轴方向上的扭矩。

式中，Fcx,i和Fcy,i分别是在x和y轴方向上粒子与叶片之间产生的接触力；rx,i和ry,i分别是与叶片碰撞的粒子到叶片中心在x和y轴方向上的位移；nc是与桨叶碰撞的粒子总数。

1.2 仿真模型的建立

由于粉体混合机主要由筒体和搅拌桨构成，因此，为了提高计算效率，该粉体混合机三维模型包括筒体和搅拌桨，筒体内径为180mm，搅拌桨为直叶平桨，桨叶直径为d，具体结构如图2所示。

图2 粉体混合机仿真计算模型图（单位：mm）

1.3 DEM参数设置

采用颗粒力学行为离散单元法数值模拟软件EDEM 2.7.1，其相关参数设置[17-18]如表1、表2所示。搅拌桨的形式为直叶平桨，搅拌筒内部填充20000个球形聚甲基丙烯酸甲酯颗粒，按照上下填充的方式分为两种颗粒，各10000个，为了清晰可视混合过程，采用蓝白两色加以区分。在初始阶段，利用颗粒工厂生成颗粒，并设置颗粒处于静止堆积状态。颗粒直径均为4mm，达到搅拌筒有效容积的35%。混合过程中，搅拌桨逆时针方向旋转，转速设置为100r/min、200r/min、300r/min、400r/min、500r/min。

表1 EDEM颗粒、壁面参数设置

表2 EDEM接触参数设置

2 结果与讨论

2.1 颗粒的混合过程

图3为颗粒在粉体混合机内不同时刻的混合状态。以搅拌桨直径d为160mm 的四直叶平桨为例，探究不同转速下颗粒的混合过程。如图3所示，两种球形颗粒的直径均为4mm，颗粒数量为20000个。在初始静止状态下，对球形颗粒按上下方式进行填充，保证两种颗粒能够处在完全分离的状态。在搅拌桨的旋转作用下，粉体颗粒沿着搅拌桨的切向方向高速运动。由于受到较大离心力作用，颗粒被推向筒壁，部分颗粒沿筒壁向上运动。当颗粒向上运动到一定高度后，克服重力做功增多，同时颗粒与器壁产生碰撞损失能量，这两个因素使得颗粒在粉体混合机中心区域向下沉降，在搅拌桨再次作用下重新被抛起，切向运动与上升运动及重力沉降的结合使得颗粒处于螺旋状的往复运动状态，进而实现粉体颗粒的混合。

混合过程中不同时刻颗粒与壁面的碰撞频次如图4 所示。在混合的初始阶段，由于搅拌桨的作用，颗粒从静到动，被抛起撞向壁面，运动紊乱，与壁面高频地碰撞、反弹。1.04s时颗粒-壁面的碰撞次数达到最大值，为24415 次。之后颗粒-壁面碰撞次数开始下降，于1.17s 时达到最小值，为14205次。这主要是由于混合初始阶段，在桨叶启动瞬间，颗粒由静止状态被突然加速推向壁面，颗粒触碰到器壁并产生反弹，之后的重力沉降使得颗粒-壁面碰撞次数也因此减少。1.5s后颗粒-壁面碰撞次数趋于稳定，颗粒-壁面的碰撞次数大约在17108处波动，粉体混合机内颗粒运动达到相对稳定状态，4s后颗粒基本混合均匀。

2.2 搅拌桨直径变化对扭矩和功率的影响

图3 不同时刻颗粒的混合状态

图4 颗粒与壁面的碰撞次数

本节以粉体混合领域常用四直叶平桨为研究对象，通过改变搅拌桨直径研究桨径变化对粉体混合过程搅拌扭矩和功率的影响。分别建立搅拌桨直径d为80mm、100mm、120mm、140mm、160mm 的5种计算模型。模型除桨径不同外，其余几何参数均一致。在4s 内模拟混合20000 个直径为4mm 的聚甲基丙烯酸甲酯球形颗粒。

图5为4s时不同直径搅拌桨分别在粉体混合机中逆时针旋转时颗粒物料的速度分布云图。从图中可以看出，颗粒速度最大处在桨叶正前方，在桨叶后方速度较小，桨叶之间形成了颗粒速度梯度。当搅拌桨直径为160mm时，搅拌桨的作用范围最大，在4s 的有限时间内，物料的混合程度远高于其他直径搅拌桨，能够更好地带动周围颗粒甚至靠近搅拌筒壁的颗粒进行混合。当搅拌桨直径较小时，靠近搅拌筒壁的颗粒仅靠颗粒间的摩擦力不足以带动周围颗粒运动。在模拟的有限时间内，无法获得颗粒混合均匀时所需功率值。因此，在后续的模拟和实验中，选用直径为160mm的搅拌桨。

图5 100r/min下不同直径搅拌桨的颗粒分布云图

图6 不同直径搅拌桨的扭矩、功率随转速变化图

如图6(a)所示，同一转速下，扭矩值随着搅拌桨直径的增大而增大。在搅拌速度和装料量相同的情况下，搅拌桨直径越大，桨叶离搅拌筒壁越近，更能带动搅拌筒壁面附近颗粒的运动。与短桨叶相比，大直径的搅拌桨与颗粒的接触面积更大，移动相同颗粒床时，不仅颗粒对桨叶的作用力更大，而且桨叶受力的力矩也更长，因此扭矩值也会更大。由图6(a)还可以看出，当桨叶直径为80mm 时，扭矩值随着转速的增加有明显降低的趋势。桨叶直径较小时，桨叶所能接触到的物料有限，搅拌桨所施加的力不足以推动整个颗粒床进行混合流动。当搅拌桨直径为80mm时，搅拌桨尖端处与筒壁之间距离最大，搅拌过程中，有较多数量的颗粒处于搅拌死区。并且，当转速较高时，颗粒由于机械搅拌作用被搅拌桨击飞，向四周抛洒，并在筒壁处堆积，到达桨叶触及不到的地方。因此，与搅拌桨接触的颗粒物料反而变少，从而导致扭矩值减小。而当搅拌桨直径较大时，在转速较高的情况下，颗粒即使被搅拌桨击飞，也依旧落回到搅拌桨的作用范围内与其发生碰撞。因此，扭矩值随着转速的增加而增大。

图6(b)所示为功率随转速变化规律图，与图6(a)中扭矩的变化规律相似，在各个转速下，功率值都随着搅拌桨直径的增大而增大。当搅拌桨直径固定时，功率值随着转速的增加急剧增大。

采用MATLAB，考虑影响混合功率的因素，对功率数值模拟结果进行拟合，如式(8)、式(9)所示。

整理成量纲为式(10)。

式中，M为质量量纲；L为长度量纲；T为时间量纲。

由此得到式(11)。

基于模拟数据，通过MATLAB 拟合可得式(12)。

式(12)即为粉体混合时的搅拌功率公式。拟合相关系数R2=0.993，偏差较小。由公式可知，功率P可由转速N、搅拌桨直径d、搅拌筒内径D及粉体密度ρ等因素计算得出，并且功率P与转速N、搅拌桨直径d、搅拌筒内径D及粉体密度ρ都呈正相关。仅转速N变化时，功率P与N3呈正比；仅搅拌桨直径d变化时，功率P与d3.81呈正比。

2.3 搅拌桨桨叶数对扭矩和功率的影响

由前节可得，搅拌桨直径d为160mm的结构可以在模拟时间内获得较好的混合结果。因此，保持搅拌桨直径参数不变，且维持相同的混合时间及相同的颗粒粒度，仅改变搅拌桨叶数目，来研究搅拌桨桨叶数对扭矩和功率的影响。通过比较不同桨叶数目来分析在粉体混合机中桨叶数目的改变对扭矩和搅拌功率的影响。桨叶数n分别为1、2、3、4、5，采用控制变量法，在4s 的模拟时间内混合20000 个直径为4mm 的球形颗粒，结果如图7所示。

由图7(a)可知，同一转速下，随着桨叶数目增加，扭矩值呈增大的趋势。桨叶数目越多，与颗粒的接触面积越大，相同力矩下受力越大，从而导致扭矩增大。当桨叶数目不变时，扭矩值随着转速增加而增大。然而，在转速达到300r/min附近时，随着转速的增加，同种桨型的扭矩值增速变缓。图8为四叶平桨在转速为500r/min时，粉体混合机内颗粒进行混合运动时的速度矢量图。图8(a)为搅拌桨刚开始转动时的速度矢量图。箭头方向表示粉体混合机内颗粒的运动方向。从图8(b)可知，混合初期，搅拌桨开始转动，转速较高，粉体混合机内的颗粒在短时间内被赋予比较大的动能，使得这些颗粒分成两部分运动。大部分颗粒被搅拌桨击飞，速度较快，向上运动；少部分颗粒速度较慢，依然停留在桨叶前方，随着搅拌桨的推动向前运动。因此，当转速较高时，更多的颗粒不断被搅拌桨击飞向上运动，与搅拌桨接触的颗粒数目反而减少，表现为搅拌桨高速转动下的扭矩值增速变缓。图8(c)可以看出，速度较高的颗粒碰到容器壁或者附近颗粒时会产生碰撞反弹，使得颗粒在运动过程中损失动能，最终因为动能的损失及重力作用而沉降，进入图8(d)所示的稳定混合状态。

图7 不同桨叶数目搅拌桨的扭矩、功率随转速变化图

功率随转速的变化的规律图如图7(b)所示，可以看出功率的变化与图7(a)中扭矩的变化规律相似。在转速等其余因素都相同的情况下，粉体混合机的功率随着桨叶数目的增加而增大。当搅拌桨桨叶数目固定时，随着转速的增加，功率值急剧增加。并且通过MATLAB 拟合综合模拟数据可得到基于上述情况的公式，推导过程与2.2 节中类似，最终得到满足上述情况的功率公式，见式(13)。

式(13)即为粉体混合时的搅拌功率公式。拟合相关系数R2=0.99，偏差较小。由公式可知，功率P与转速N、粉体密度ρ及桨叶数目n都呈正相关。仅转速N变化时，功率P与N3呈正比；仅桨叶数目n变化时，功率P与n2.27呈正比。

2.4 综合影响因素

为了更加深入探究直叶桨式粉体混合机中特性参数对搅拌功率的影响，对搅拌桨转速N、搅拌桨直径d、搅拌筒内径D、粉体密度ρ及桨叶数目n与功率值的特性关系进行综合考虑。基于模拟数据对搅拌功率进行多维度探究，并推出较为综合的搅拌功率公式。

图8 固定转速下的y轴截面上的速度矢量变化图

在粉体搅拌过程中，综合考虑影响搅拌功率的主要因素包括N、D、d、ρ、n等，则有式(14)。

推导过程与2.2 节中类似，最终得到满足上述情况的功率计算公式，见式(15)。

式(15)即为粉体混合时的搅拌功率公式。拟合相关系数R2=0.95，偏差较小。由公式可知，功率P可由转速N、搅拌桨直径d、搅拌筒内径D、粉体密度ρ及桨叶数目n等因素计算得出，且功率P与转速N、搅拌桨直径d、搅拌筒内径D、粉体密度ρ及桨叶数目n都成呈相关。

2.5 实验验证

为了验证以上数值模拟的准确性，根据上文的分析结果，采用搅拌桨直径d为160mm 的四叶平桨，搭建粉体混合机的实验装置。按照实验装置尺寸进行1∶1 数值模拟建模和混合过程模拟，以便更好地分析模拟结果。如图9所示，试验台主要包括OS40-Pro LCD 数控顶置式电子搅拌器、筒体、桨叶、实验室用容器夹等。其中OS40-Pro LCD 数控顶置式电子搅拌器是集电动机、转速仪、扭矩传感器为一体的装置，转速范围可在50～2200r/min之间调节，控制精度为±3r/min，LCD 显示屏通过mode 按钮进行切换可直接显示转速和扭矩值等相关信息，扭矩测量精度为1N·cm。为了便于观察粉体混合机内颗粒运动状态，筒体材料选用亚克力管材，无色透明，并且具有一定的机械强度。如图10 所示，采用等量的红、棕两种直径4mm 的聚甲基丙烯酸甲酯颗粒放入试验设备中，数量共20000个，采用上下方式进行填充。

控制桨型及填料高度一致，分别进行转速为100r/min、200r/min、300r/min、400r/min 和500r/min时的混合实验。为了定量比较实验和模拟的顶面颗粒分布，将图像分成26 个块，并进行图像分析以识别每个块中的红色粒子。根据红色粒子的面积，得到每个块中红色粒子的比例。模拟与实验的相对偏差可用式(16)表示。

图9 粉体混合机实验设备

图10 实验与模拟初始填充状态

式中，Eri为相对偏差；psi为模拟中红色颗粒的比例；pei为实验中红色颗粒的比例。

图11 中最右侧一栏显示了每个块中的Eri即红色颗粒的模拟比例与实验值的偏差，模拟值和实验值之间的差异随着色块变轻而减小。模拟的颗粒分布与实验结果吻合较好。

图11 实验和模拟表面颗粒分布对比

从图12 中可知，实验结果与模拟数据变化趋势类似，模拟与实验的扭矩和功率值较为接近，且实验所得结果略高于模拟数据。进一步，对实验结果和模拟数据进行误差分析，结果表明，实验值与模拟值偏差范围在7%～10%。这是因为在实际的混合实验中，粉体混合机的消耗功率包括电动机空转时所消耗的功率和移动颗粒床所需的功率，并且在实验过程中发现，颗粒和壁面间产生了一定程度的静电效应，而静电效应使得筒体壁面对颗粒产生了吸附作用，增加了搅拌阻力。而数值模拟则是一种理想状态，忽略了电动机空转功率及静电的存在。

图12 不同转速下实验与模拟结果对比

以搅拌桨直径d为160mm的四叶平桨为例，对实验测得的功率值与基于模拟数据拟合的公式计算得到的功率值进行对比。由图13 可以看出，实验值与计算值相比偏差小于10%，说明所推导的功率计算公式具有一定可靠性。

3 结论

采用DEM 模拟的方法研究直叶桨式粉体混合机中粉体混合时搅拌功率和扭矩的变化特性，比较了粉体混合机中搅拌桨的转速、直径、桨叶数目等结构特性对粉体混合时搅拌功率和扭矩的影响，并使用MATLAB 软件拟合出功率计算公式。最后，在模拟结果的基础上进行颗粒混合实验，对比分析两种过程下粉体混合机内搅拌功率的差异，验证拟合公式的准确性，得出以下结论。

图13 不同转速下实验与模拟结果对比

（4）后续研究将继续针对颗粒粒径、颗粒形状、颗粒密度及充填量等因素开展，期待拟合出更具工程普遍意义的粉体搅拌功率的计算公式。