基于DEM-CFD 耦合的谷物清选模拟研究

赵 磊,马学东,郭柄江,于海川

辽宁科技大学 机械工程与自动化学院,辽宁 鞍山 114051

谷物清选是农业生产中提高收获质量，降低谷物掺杂率的重要工序，是利用谷物、茎秆等成分之间的物理特性的差异，将混有短、碎茎秆、颖壳和尘土等细小杂物从谷物中分离出来的操作过程[1-3]。

农业中清选装置一般分为风选装置和风筛装置，其中风选装置是依靠气流将具有不同物理特性的物料分离开来，而风筛装置则是结合气流与振动筛的联合作用将物料分离。汉代史游《急就篇》有“碓石岂扇颓舂簸扬”说，此处之“扇”便是古代的风选装置。此外，元代王祯的《农书》、明末宋应星的《天工开物》中均有对清选装置构造及使用方法的详细介绍[4]。

由于谷物在清选过程中存在大量动能交换，且流体自身变化、颗粒间的碰撞、流体与颗粒之间相互影响所形成的耦合作用使得整个系统的物理特性极其复杂，因此在对谷物清选的模拟研究中，如果单纯采用DEM 或CFD 进行模拟[5]，则无法描述气流与谷物之间的相互作用，或只能将谷物视为多孔介质模型，不能准确计算谷物颗粒模型对气流的影响，宁新杰通过分析现有谷物清选装置的研究现状，发现单一DEM 模拟与CFD 模拟存在局限性，不能全面体现流场与物料之间的相互影响，并指出采用DEM-CFD 气固两相流耦合方法模拟将是未来清选理论研究的发展方向[6]。

Yuan 等利用DEM-CFD 气固耦合方法，模拟分析了稻谷脱粒混合料在圆筒筛清选过程中的运动行为及筛分特性，结果表明，入口气流速度对物料轴向平均速度、筛分质量影响显著[7]；江涛等基于DEM-CFD 耦合方法，对三种不同结构的风筛装置进行了仿真模拟，通过对比得出了清选效率最佳的结构方案[8]；王立军等采用DEM-CFD 耦合方法模拟了玉米籽粒在贯流阶梯式振动筛的筛分过程，确定了最优筛分参数，并通过对比试验，验证了模拟仿真的可行性[9]。上述研究均以风筛装置为研究对象，对于风选装置的研究还比较少。

由于风筛装置相比风选装置结构复杂、成本高昂、能耗较大，因此研究风选装置对于节能、降低成本有积极作用。本研究利用DEM-CFD 气固两相流耦合方法，以风选装置为研究对象，对谷粒、短茎秆、碎茎秆三元颗粒的清选过程进行仿真模拟，结合空气动力学，分析了三种物料在流场中的运动状态及分离机理，讨论了气流速度、气流倾角对谷物含杂率及夹带损失的影响。

1 模型描述

1.1 数学模型

CFD-DEM 进行耦合模拟时，主要耦合模型包括Lagrangian 模型和Eulerian 模型[10]，其中Lagrangian 模型使用的是单向流框架，不考虑颗粒相体积分数，Eulerian 模型采用多相流框架求解，有体积分数方程，并且考虑颗粒对流场的影响，结合气流穿过谷层时，谷物颗粒对气流的阻碍作用明显，因此采用Eulerian 耦合模型对谷物及其杂质的清选过程进行数值模拟。

Eulerian 模型中，流体的体积分数项和运动微分方程分别为[11]：

式中，ρ为气体密度；t为时间；u为流体流速；ε为气体的体积分数项；P为气体微元上的压强；g为重力；μ为粘滞系数；∇为哈密顿微分算子；S为动量源项。

动量源项S为作用在网格单元内气流阻力的总和，其表达式为[12]：

式中，Fi为第i个颗粒对气流的阻力，V为网格单元的体积。

1.2 颗粒接触碰撞模型

计算颗粒力学模型可以描述颗粒间的相互作用和接触力学行为。考虑到颗粒间的接触且颗粒速度基于接触力改变，本文采用软球干接触模型和Hertz-Mindlin(no slip)接触理论[13]。根据牛顿第二定律，第i个颗粒的运动方程为[14]：

式中：Vi和ωi分别为颗粒i的速度和角速度；Ii和mi分别为颗粒i的转动惯量和质量；g为重力加速度，P 为颗粒与气流相对运动时受到的作用力。

Fn,ij为法向分力、Ft,ij为切向分力，Tt,ij为切向力矩、Tr,ij为滚动摩擦力矩，根据三方程线性弹性-阻尼模型，将每种作用力简化为一个弹簧、一个阻尼以及一个滑动器，其中各种力和力矩的数学描述为[15]：

式中，k为刚度系数，V为颗粒的速度（矢量），δ为颗粒之间的位移变形（矢量），a为扭转变形（矢量），η为阻尼系数，f为摩擦系数，L为重叠量。下标ij表示颗粒i与颗粒j之间，i和j分别表示颗粒i和颗粒j，n表示法向，t表示切向，r表示周向或滚动，s表示滑动。

颗粒与气流相对运动时受到的作用力P大小为：P=kgρAv2=kgρA(vq-vw)2(10)

式中，kg为阻力系数，与物体形状、表面特性和雷诺数有关；ρ为空气密度；A为物体的受风面积，即物料在气流方向的投影面积；v为物体与气流的行对速度；vq为气流速度；vw为物料速度。

1.3 几何模型构建

运用SolidWorks 软件创建风选装置模型，风选装置模型如图1 所示，装置采用单进风口结构，整体壁厚2 mm，箱体长度为200 mm，高度为160 mm，厚度为80 mm，进料斗宽口处长度为100 mm，宽度为70 mm，窄口处长度为50 mm，宽度为15 mm；进风口长度60mm，宽度60 mm。将模型导入ANSYS Workbench 中划分网格，网格如图2 所示。

图1 风选装置Fig.1 Winnowing device

图2 网格Fig.2 Grid

考虑到EDEM 软件自身建模的缺陷，所以选择谷物清选除杂中成分含量较高的谷粒、短茎秆、碎茎秆为研究对象。由于EDEM 中的颗粒均采用球形，故采用“多球丛聚法”对三种颗粒进行简化、重叠组合、填充[16]，如图3 所示，其中谷粒由13 个不同粒径的球体填充而成，短茎秆由2 个半径为2 mm 和19 个半径为1 mm 的球体填充而成，碎茎秆由36 个半径为0.5 mm 的球体填充而成。

图3 物料颗粒三维模型图Fig.3 3D models of material particles

1.4 模拟参数设置

EDEM 中物料颗粒的力学特性参数及接触系数如表1 和表2 所示[17,18]，风选装置材质选用钢，由于短茎秆和碎茎秆为同种物质，表1、表2 中均用茎秆代表。谷物、短茎秆与碎茎秆比例为4：1：0.25，设定谷粒生成速率为1200 个/s，短茎秆生成速率为300 个/s，碎茎秆的生成速率为75 个/s，仿真时间步长设置为瑞利（Rayleigh）时间步长的33.9182%，即4e-6 s，仿真时间总时长为10 s。Fluent 中模拟仿真采用标准的k-ε湍流模型。时间步长设定为EDEM 的100 倍，即4e-4 s。

表1 材料的物理参数Table 1 Physical parameters of the material

表2 材料的接触系数Table 2 Contact coefficient of the material

2 仿真模拟及分析

2.1 气流速度分析

在风选装置的入料口处设置颗粒工厂，谷物、短茎秆及碎茎秆在入料口处自由下落，进风口气流速度设置为5 m/s。图4为2 s 时风选装置内部的物料位置瞬态图，在出口2 和出口3处分别设置含杂率统计区域和损失率统计区域。从图4 中可以看出，物料落入气流作用区域后，在水平气流作用下3 种物料颗粒呈现出不同的运动轨迹，谷粒全部落入出口2 中，其中掺杂部分短茎秆，而出口3 收集到的全部为短茎秆和碎茎秆。

由于谷物、短茎秆、碎茎秆的空气动力特性不同，竖直下落的物料受到水平气流作用后会呈现不同的运动轨迹，物料在风选装置内受到自身重力G，空气浮力P´以及水平气流作用力P，三个力的合力为F，如图5 所示，物料将沿着F的方向运动，其运动轨迹为抛物线，物料的运动方向角为α。

图4 物料颗粒位置瞬态图(t=2 s)Fig.4 Material particle positions at 2 s

图5 物料颗粒受力图Fig.5 Diagram of material particle force

如果空气浮力P´忽略不计，则[19]：

当水平气流作用力P不变时，重力越大，α就越大，即物料颗粒的运动方向角α越大。空气动力学中，tanα为物料在流场中的飞行系数，物料的粒度、密度等物理性质不同，在同一气流中的飞行系数也不相同，当气流速度一定时，飞行系数越大的颗粒在气流驱使下所做的水平位移越大。由公式(11)可知，当α∈(0,π/2)时，物料的飞行系数与自身质量成反比，因此质量较大的谷粒下沉趋势明显，率先落入出口2 中，而质量相对较轻的短茎秆和碎茎秆在水平气流作用力的驱使下做平抛运动，落入出口3 中。

为了定量描述出口2 处的含杂率，引入短茎秆和碎茎秆的体积浓度作为含杂率的衡量标准，其中，短茎秆和碎茎秆的体积分数表达式为：

式中，Vs(s,t)，Vf(s,t)分别为区域s内t时刻短茎秆和碎茎秆的体积，Vt(s,t)为区域s内t时刻所有物料的体积。

图6(a)是通过统计分析后所得到谷物的含杂率，从图中可以看出，当水平气流速度为5 m/s 时，尽管出口2 处各时刻含杂率存在零点，但整体不稳定，波动较大，含杂率峰值出现在3 s 到4 s 之间，达到36.178%，平均含杂率为10.575%。利用水平气流清选谷物及其杂质的过程中，出口3 收集到的杂质中会夹杂部分谷粒，即存在夹带损失。由于出口1 为出风口，且仅有少量碎茎秆进入出口1，故定义出口3 收集到谷粒数量与出口2、出口3 收集到谷粒总数的比值为夹带损失率。计算结果显示，当气流速度设置为5 m/s，夹带损失率为0.066%，虽然夹带损失率比较理想，但出口2 处谷粒的含杂率峰值大，波动明显，且10.575%的平均含杂率表明水平气流速度为5 m/s 时的清选效果不佳。

图6(b)给出了水平气流速度为7 m/s 时出口2 处清选所得到的谷物含杂率，由图可知，水平气流速度为7 m/s 时，各时刻含杂率零点较多，含杂率峰值存在于9 s 到10 s 之间，达到23.79%，计算后得到平均含杂率为2.162%，出口3 处的谷粒夹带损失率为0.351%。图6(c)为水平气流速度为9 m/s时，含杂率统计区域内收集到的谷物的含杂率，从图中可以看出，当水平气流速度为9 m/s 时，谷物含杂率波动较小，整体趋向于理想状态，峰值为15.914%，计算结果显示平均含杂率为0.307%，谷粒夹带损失率为1.275%。

图6 谷物含杂率Fig.6 Impurity rate of grain

3 组仿真的数据结果如表3 所示，通过对比表3 中的数据发现，气流速度调整为7 m/s 时，出口2 处谷物的平均含杂率相比气流速度为5 m/s 时下降了8.413%；夹带损失率较气流速度为5 m/s 时上升了0.291%。气流速度调整为9 m/s 时，清选过程中谷物含杂率标准差为7 m/s 的1/2，表明气流速度为9 m/s 时谷物的含杂率波动更小，平均含杂率降低了1.855%，夹带损失率上升了0.924%。由于气流速度的增加，即水平气流作用力P 变大，故物料运动轨迹与重力之间的夹角α与物料水平方向上的位移均变大，短茎秆和碎茎秆中夹杂的谷粒数量增加，因此增大水平气流速度可以降低谷物含杂率的同时，夹带损失率增大。

表3 仿真结果Table 3 Simulation results

实际生产中，为了降低粮食的含杂率和损失率，往往会对含有夹带损失谷粒的杂质进行二次清选，综合考虑三种气流速度下的含杂率和夹带损失率，结果表明气流速度为9 m/s 时的清选效果更优。

2.2 气流倾角讨论

为了探讨气流倾角对物料运动行为及清选质量的影响，将风选装置进风口由水平改为倾斜10°，即气流倾角β变为10°，图7(a)、(b)分别为气流速度为9 m/s 时，进风口倾角为10°条件下和进风口水平条件下的气流速度矢量分布图，从图中可以看到进风口倾斜时，仅入口处的气流呈现一定角度，相比水平条件下，进风口倾斜10°时，风选装置内部气流整体高度有所抬升，顶部内壁气流速度较大，流场整体差异不大。

图7 气流速度矢量分布Fig.7 Airflow velocity vector distribution

图8 v=9 m/s,β=10°时的谷物含杂率Fig.8 Impurity rate of grain at v=9 m/s,β=10°

图8 给出了气流倾斜角度为10°条件下，各时刻出口2 处的谷物含杂率。当气流倾角设置为10°时谷物的含杂率峰值为15.571%，平均含杂率为0.282%，计算后得到谷物的夹带损失率为1.583%。

通过对比进风口为水平条件下与进风口倾斜10°时的谷物平均含杂率及夹带损失率可知，平均含杂率降低了0.025%，夹带损失率上升了0.308%。根据空气动力学原理，相同条件下，物料在水平气流中的飞行系数要小于在倾斜气流中，即物料颗粒受气流作用所做的水平方向位移较大，因此相同气流速度下，气流倾斜10°时出口2 处谷物的含杂率下降，而出口3 处风选夹带损失率上升。

3 结论

运用DEM-CFD 耦合方法，对不同气流速度及气流倾斜角度下的谷物清选过程进行了仿真模拟，结合空气动力学分析了谷粒、短茎秆、及碎茎秆在流场中的受力状态及运动趋势，通过统计计算得到的不同参数下谷物的含杂率和夹带损失率，得到以下结论：

（1）水平气流速度不变时，物料自身重力越大，其运动方向角越大，且物料质量与其飞行系数成反比，即质量越小，飞行系数越大，物料水平方向位移越大；

（2）水平气流速度设置为5 m/s 时，谷物的平均含杂率为10.575%，夹带损失率为0.066%。通过对比水平气流速度分别为5 m/s、7 m/s、9 m/s 时的含杂率和夹带损失率发现，增大气流速度会降低谷物的含杂率，而夹带损失率随之增大；

（3）水平气流速度设置为9 m/s，气流方向角为10°时，谷物的平均含杂率为0.282%，夹带损失率为1.583%。相比水平气流条件下有所不同，由于物料在倾斜气流中的飞行系数要大于在水平气流中，即物料水平方向位移较大，因此相同条件下，气流方向角为10°时平均含杂率降低，而夹带损失率有所升高；

（4）DEM-CFD 耦合方法可以有效用于谷物清选过程的模拟分析，所得到的参数为风选装置的设计提供了参考依据。