基于DEM-CFD耦合的辣椒清选仿真研究

张德俊， 张小明， 吴荻， 林蜀云， 张太华， 徐卫平*

（1.贵州师范大学机械与电气工程学院，贵阳 550001； 2.贵州省山地农业机械研究所，贵阳 550007）

辣椒产业是贵州最具代表性的传统经济优势产业之一，是贵州现代农业的重要经济支柱，同时，也是农民脱贫、创业增收的重要途径[1]。2020 年，贵州省辣椒种植面积36.3 hm2，产量724 万t，产值242亿元,生产、加工、销售规模居全国首位，贵州正从“辣椒大省”向“辣椒强省”迈进[2]。由于贵州山区辣椒种植面积大、田间地块小及人工收获成本高等问题，设计并制造了山地自走式履带辣椒收获机，通过机收降低了人工收获成本，但机收的辣椒含杂率较高，果实清选问题亟待解决。目前，谷物、花生等农产品已有成熟的风选设备并投入应用，辣椒收获机上虽然设有风选装置，但参数的选择主要是凭经验，辣椒风选效果不好。

清选方式一般有两种：一种是气流清选，主要是利用混合物的不同悬浮性质；另一种是利用气流和振动筛的联合作用来实现清选。气流清选是以空气为主要作用介质，在气流的作用下，根据颗粒密度和颗粒大小对各种颗粒材料进行空气分选的过程，基本原理为通过气流将较轻的颗粒水平的吹向远处，而较重的颗粒吹到近处[3]。

近年来，随着计算机技术的飞速发展，离散元法（discrete element method，DEM）和计算流体力学（computational fluid dynamics，CFD）在风选设备的分析中应用广泛。敬志臣等[4]通过Fluent 软件对小麦进行研究，采用旋风分离筒对分离过程进行仿真分析；李兴凯等[5]利用EDEM（extended DEM）软件对小区小麦联合收割机的脱粒机构进行仿真分析。使用CFD 技术可以模拟气流场，使用DEM 可以模拟粒子碰撞和接触过程。但是，它们只能以单一方式模拟清选过程，不能准确、全面地反映风选室中作物在气流、容器壁和颗粒碰撞接触作用下的复杂运动，DEM-CFD 耦合模拟技术将气流场、风选室和颗粒作为整体系统进行研究，可以准确分析出风选过程。张凤奎等[6]利用CFDDEM 耦合方法研究了红枣的负压气吸拾捡仿真并进行了分析；郭柄江等[7]通过EDEM 与Fluent 耦合方法研究了谷物风选，采用气吹与风筛组合的风选方式，不同风速参数下稻米清选装置工作性能的指标通过风选后稻米损失率和含杂率来进行衡量；丁贺贺等[8]通过DEM-CFD 耦合方法研究了谷物在垂直气流吹气下的仿真分离。因此，针对辣椒收获机收获的辣椒含杂率比较高的问题，设计了一种辣椒清选分离装置，采用DEM-CFD耦合方法，对朝天椒风选（吹气式）分离过程进行模拟仿真，旨在为设计及优化辣椒收获机的辣椒风选分离装置提供参考。

1 材料与方法

1.1 风选装置、原理及阻风效应

1.1.1风选装置 辣椒风选分离装置由进风口、进料口、清选室、出口Ⅰ、出口Ⅱ、出口Ⅲ组成，模型结构如图1 所示。使用SolidWorks 三维设计软件建模，将风选装置模型输出为XT 格式并导入ICEM CFD 网格划分软件中，添加进风口、物料入口、辣椒出口Ⅰ、辣椒茎秆及辣椒叶出口Ⅲ，设置完成后划分网格如图2所示。

图1 风选装置Fig. 1 Wind-sorting device

图2 风选装置网格Fig. 2 Grid of wind sorting device

1.1.2风选原理 在水平清选气流中，4 种力作用于物料颗粒，即水平风力FD、重力FG、阻力FR和浮力F′（图3）。空气密度很低，物料颗粒的体积也很小，因此相应的浮力F′很小，对气流中其他物料颗粒的径向运动速度影响不大，所以浮力F′可以忽略，因此，垂直方向上合力FGr与水平风力FD的合力等于F，合力F与风力FD的夹角为β。此时，物料颗粒沿合力F方向运动，轨迹大致呈抛物线。物料颗粒的最大运动距离与夹角β的值有关：夹角β越小，物料颗粒运动得越远，夹角β越大，颗粒运动距离越近。其中，Va为平均风速，V0为颗粒沉降速度，V为颗粒在合力作用下整体的速度。

图3 水平气流中的颗粒受力分析Fig. 3 Force analysis of particles in horizontal airflow

颗粒沉降速度V0可表示如下[9]。

式中，d为颗粒的粒子直径，mm；ρs为颗粒密度，kg·m-3；ρ为风选的气流密度，kg·m-3；φ为阻力系数；g为颗粒重力加速度，m·s-2。

用沉降速度V0和风速Va的比值表示β的正切值。

从式(2)可以看出，在颗粒水平风速保持不变、颗粒直径大小相近的情况下，低密度颗粒的tanβ值较低，颗粒可以运动得更远，而高密度颗粒的tanβ值较高，颗粒运动得就会更近，不同的物料颗粒可通过气流按其密度、形状和大小进行分离。

1.1.3阻风效应 将多颗粒连续运动模型简化为颗粒叠加的粒子组模型（图4）。从图4 可明显看出，当气流中前端的流动力作用于前端粒子组时，气流受到前端粒子组的阻碍而形成纵向流，沿横向流过前端粒子组上、下表面。前端粒子组对后端粒子组的阻力屏障形成侧流，气流中的流动力将绕过后端粒子组。前端粒子群抑制气流方向的向前运动，风力不作用于后粒子群，产生阻碍效应（即阻力效应）影响后端粒子组的整体运动、受力和运动轨迹的改变。

图4 风选颗粒阻风效应简化原理Fig. 4 Simplified principle of wind-selected particle windblocking effect

1.2 DEM-CFD耦合模型构建

1.2.1Fluent 流体模型 采用欧拉耦合模型对辣椒及其杂质的风选分离过程进行模拟，流体体积分数项和流体运动微分方程表示如下[10]。

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

动量源项S为气流作用在单元网格内部气流阻力的总和[11]。

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

1.2.2EDEM 固体颗粒接触碰撞模型 在EDEM软件中对颗粒力学模型间的接触力学行为和相互作用力进行仿真分析。由于颗粒间的速度随接触力而发生变化，本文采用Hertz-Minderlin 接触理论和软球干接触模型[12]，根据牛顿第二定律获得第i个颗粒的运动方程[13]。

式中，g为重力加速度；P为颗粒在气流中运动时受的相互作用力；Ii和mi分别为颗粒i的转动惯量和质量；Vi和ωi分别为颗粒i的速度和角速度；Ft,ij为切向分力；Fn,ij为法向分力；Tr,ij为滚动摩擦力矩；Tt,ij为切向力矩。

根据弹性阻尼的线性三方程函数模型，每个力都以弹簧、阻尼和滑块的组合形式将其简化出来表示[14]。

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

颗粒在气流中受的相互作用力P表示如下。

式中，ρ为气流密度；kg为阻力系数，与物料结构形状、表面特性和雷诺数等有关；A为物料迎风面积；vw为物料速度；v为物料在气流中的相对速度；vq为气流速度。

1.3 DEM-CFD 耦合流程介绍及辣椒风选效果评价指标

本文用CFD-DEM 耦合方法，即Fluent 与EDEM的耦合，其流程如图5所示。

图5 DEM-CFD耦合流程图Fig. 5 Flow chart of DEM-CFD coupling

为深入探讨风速大小对辣椒风选效果的影响，采用控制变量法保证其他的各项条件参数相同，模拟辣椒风选的风速数值大小分别为15、20和25 m·s-1时的风选过程。

辣椒清选装置工作性能的主要衡量指标是风选后辣椒损失率（El）和含杂率（Ed），因此以辣椒损失率El、含杂率Ed为辣椒风选装置的性能评价指标[15-16]。

式中，W1为仿真过程中通过出口Ⅰ的辣椒质量；W2为仿真过程中通过出口Ⅱ的辣椒质量；M为仿真过程中通过出口Ⅰ的物料总质量。

1.4 模拟参数设置

通过WDW3100 微机控制电子万能试验机对朝天椒、辣椒叶及辣椒茎秆进行压破试验和拉伸试验求取泊松比与剪切模量[17-18]，密度用天平和量筒法测量；辣椒、辣椒茎秆和辣椒叶与钢板的摩擦接触参数采用斜板试验（分别为滑动法和滚动法）测量；辣椒、辣椒茎秆、辣椒叶与钢板之间的恢复系数采用碰撞法测量[19-20]。EDEM软件中材料的力学性能参数和材料颗粒的接触系数见表1和表2，辣椒清选分离装置材质为钢材。模拟中辣椒、辣椒叶、辣椒茎秆的生成比例设置为6∶1∶0.25，辣椒的生成速度设置为1 200个·s-1，辣椒叶的生成速度设置为200 个/s，辣椒茎秆的生成速度设置为50 个·s-1。EDEM中的仿真时间步长设置为29.6784%，即0.000 007 s，总仿真时间为6 s，Fluent 软件中的仿真模型采用标准的k-ε湍流模型，时间步长设置为EDEM 的100 倍，即时间步长设置为0.000 7 s，总步数为8 571步。

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

表2 相关颗粒的材料接触系数Table 2 Material contact coefficient of related particles

2 结果与分析

2.1 辣椒颗粒运动轨迹分析

以风速大小为20 m·s-1的耦合仿真最终结果为例，辣椒颗粒的运动轨迹如图6 所示。可以看出，辣椒、辣椒叶和辣椒茎秆由颗粒工厂出料生成，然后自由运动下落并进入风选装置，当生成的颗粒通过气流通道流体区域时，颗粒质量较轻的辣椒叶颗粒和辣椒茎秆颗粒能随气流运动，被吹到出料口Ⅱ、Ⅲ从而下落，较重的辣椒在重力作用下继续自由下落，最终从出料口Ⅰ中掉出，从而实现辣椒、辣椒叶和辣椒茎秆的分离。从图6 中颗粒轨迹图可以直观地观察风选过程中颗粒轨迹，辣椒的红色轨迹线是一条略带斜度的直线，表明辣椒受风影响较小，水平位移较小，辣椒可以从出口Ⅰ自由落下；辣椒叶的绿色轨迹线是一条转弯的斜线，表明辣椒叶受风影响较大，辣椒叶随风扬动，从出口Ⅲ出来；代表辣椒茎秆的蓝线轨迹线是斜线，表示辣椒茎秆受风影响较明显，辣椒茎秆随风运动，从出口Ⅱ被吹出。

图6 辣椒颗粒的运动轨迹Fig. 6 Schematic diagram of the winnowing of pepper particles

2.2 风速20 m·s-1时辣椒风选效果分析

如图7 所示，在EDEM 软件中设置区域Ⅰ为辣椒含杂率的统计区域，统计0～6 s 内通过该区域的物料总质量M和辣椒质量W1；区域Ⅱ设置为辣椒损失率的统计区域，统计辣椒在0～6 s 内通过该区域的辣椒质量W2。

图7 物料颗粒位置瞬态图（t=1 s）Fig. 7 Transient diagram of material particle position （t=1 s）

从辣椒风选含杂率曲线（图8）可以看出，0～0.5 s 间的含杂率曲线在X轴上，这是因为该时间段为风选颗粒产生及开始下降的风选过程，经过0.5 s 的缓冲，颗粒经过风力风选位置并进行风选，风选后的辣椒颗粒经过辣椒含杂率统计区域并从出口Ⅰ掉落。通过风选模拟仿真分析可知，经过风力风选之后，辣椒颗粒的含杂率降明显，通过风选装置辣椒出口Ⅰ的辣椒平均含杂率为4.624%，通过风选装置其他出口的辣椒平均损失率（图9）为1.798%。

图8 风速20 m·s-1的辣椒含杂率Fig. 8 Impurity rate in pepper with wind speed of 20 m·s-1

图9 风速20 m·s-1的辣椒损失率Fig. 9 Pepper loss rate at wind speed 20 m·s-1

2.3 不同风速对风选效果的影响

从图10 可以看出，在较低风速（15 m·s-1）下，在辣椒出口处发现了较多杂质。辣椒颗粒分离度低；在20 m·s-1的较高风速下，混入的辣椒杂质很少，辣椒颗粒的分离度较高；而当风速为25 m·s-1时，几乎没有杂质混入辣椒出口处，然而，在茎秆出口处大量的辣椒被吹出，辣椒损失率过高。

图10 不同风速下的辣椒颗粒分布Fig. 10 Distribution of pepper particles at different wind speeds

不同风速下辣椒含杂率曲线及辣椒损失率曲线如图11、12所示，通过对辣椒在不同风速下的仿真结果分析，可知当风速大小为15 m·s-1时，辣椒含杂率峰值出现在4~5 s 之间，达到22.234%，风选后的辣椒平均含杂率为11.774%，辣椒平均损失率为0.959%；当风速大小为20 m·s-1时，辣椒含杂率峰值出现在5~6 s 之间，达到8.713%，辣椒平均含杂率为4.624%，辣椒平均损失率为1.798%；当风速大小为25 m·s-1时，辣椒含杂率峰值出现在2~3 s之间，达到3.554%，辣椒平均含杂率为1.676%，辣椒平均损失率为6.476%。根据模拟数据分析可知，在不考虑损失的情况下增加风速强度可以增加风选的效果。

图11 不同风速下的辣椒含杂率Fig. 11 Impurity rate of pepper at different wind speeds

图12 不同风速下的辣椒损失率Fig. 12 Loss rate of pepper at different wind speeds

不同风速下的辣椒仿真结果如表3 所示。分析发现，当风速从15 m·s-1增大到20 m·s-1时，风选装置出口Ⅰ处辣椒的平均含杂率下降了7.150%，而平均损失率增加了0.839%；当模拟的风速再次从20 m·s-1增大到25 m·s-1时，辣椒的平均含杂率下降了2.948%，而平均损失率上升到4.678%。这是由于气流速度增加，即气流水平力FD增大，故合力F和风力FD之间的角度β及物料颗粒在水平方向上的位移变大，从而导致辣椒损失率变大，因此增加气流速度可以降低辣椒的杂质含量，同时也会导致辣椒的损失率增加。

表3 不同风速的仿真结果Table 3 Simulation results at different wind speeds

由以上比分析可知，风速的增加可以增强辣椒的风选效果，但是随着风速的增加，会导致辣椒从辣椒叶出口吹出的风险。综合考虑辣椒损失率及含杂率，在风选时应选择分离率较高及损失率较低的风选速度，即在20 m·s-1风速时分离效果最佳，此时，辣椒的含杂率为4.624%、损失率为1.798%。

3 讨论

邹道钟等[21]设计了一种适用于板椒收获机的清选系统，系统总体分为上、下清选机构，复脱机构以及风机风选机构，由于复脱机构中定、动弹齿的高速冲击，更容易损伤板椒果皮，增大了采收时的损伤；赵飞龙等[22]使用液压动力代替传统的机械动力来为风扇提供动力，采用PID（proportion integration differentiation）反馈控制速度回路可以快速响应动态负载变化，减少液体对风机回路的影响，可提高风选装置的使用寿命及提高风选装置的输入动力，能有效提高风选效果。

前人对辣椒的风选研究还停留在凭经验的机械设计及动力分析上，对风选装置中物料风选过程的运动特性、流场分布特性之间的关系研究较少。近年来，EDEM 与Fluent 软件耦合法逐渐成为研究流固两相间流动性的重要研究手段和方法[23-24]，国内外学者使用CFD-DEM 耦合方法研究农作物的风选，并取得了一些成果：Yuan 等[25]通过DEM-CFD 耦合方法研究了谷物和茎秆的运动行为和筛分特性；蒋恩臣等[26]模拟了联合收割机分离室中谷物的运动并进行了试验，表明DEMCFD 耦合方法适用于模拟谷物清选过程；王立军等[27]使用DEM-CFD 耦合方法研究了玉米清选并使玉米粒的损失率得到了降低。为了进一步研究辣椒的清选，本研究使用DEM-CFD耦合方法对辣椒清选装置进行仿真。通过非结构化的ICEMCFD 网格模型将仿真模型划分网格并导入EDEM和Fluent 软件中，并在2 个软件之间建立耦合关系以进行模拟仿真。Fluent 软件中流体的耦合模型采用欧拉模型，欧拉耦合模型使用多颗粒相流耦合框架方程来进行求解，有考虑体积能量分数，该方程考虑不同颗粒相间对流场的相互影响[28]，结合通过进入辣椒的表层进行循环的单向气流，辣椒层的颗粒对循环气流产生阻碍力的作用明显。EDEM 软件中的辣椒混合物颗粒在建模时，考虑到EDEM 软件本身自动模拟的一些不足，选择辣椒清选除杂过程中成分含量高的辣椒、辣椒叶和辣椒茎秆作为主要研究对象进行物料粒子建模。由于EDEM 软件中的所有粒子模型都用球体表示，因此采用“多球体类聚法”对3 种不同的研究粒子进行简化、重叠和填充，以代替实际风选的粒子模型[29]。本文采用CFD-DEM 耦合法模拟辣椒在不同风速下的清选分离过程，结合空气动力学分析辣椒、辣椒叶、辣椒茎秆在流场中的受力情况和运动趋势。上述研究表明，当风选速度不变时，物料颗粒本身的重力越大，其运动的方向角越大，则物料颗粒的水平位移就越大；此外，也验证了DEM-CFD 耦合方法对于研究辣椒清选分离是有效的，得到的清选分离参数可以为辣椒清选分离装置的设计和优化提供指导。