绿豆种子离散元仿真参数标定与排种试验

张胜伟 张瑞雨 陈天佑 付 君,2 袁洪方,2

(1.吉林大学生物与农业工程学院， 长春 130025； 2.吉林大学工程仿生教育部重点实验室， 长春 130025)

0 引言

长期以来，绿豆作为小宗农作物，由于其种植地区分散、集中种植面积较小、机械化种植水平较低等问题，严重制约了我国绿豆产业发展[1]。为提高绿豆机械化种植水平，近年来国内学者对绿豆精密排种器进行了相关研究[2-3]。但与主粮作物排种器的研究相比，绿豆精密排种器设计及理论研究仍较少，且实际生产中多采用勺轮或窝眼轮等机械式排种器进行播种，存在播种精度低、伤种率高等问题，而气力式排种器可有效克服上述问题[2]。

近年来EDEM-Fluent气固耦合仿真方法在气力式精密排种器研究中已被广泛运用，可直观分析种子在复杂环境中的受力及运动规律，从而对排种器结构或工作参数予以优化[4-10]。在利用气固耦合仿真法对精密排种器进行仿真研究时，需定义仿真模型的本征参数(密度、剪切模量、泊松比等)和接触参数(滚动摩擦因数、静摩擦因数、碰撞恢复系数等)[11-16]，其本征参数与实际试验测量值基本一致，而接触参数难以通过实际试验获得，且接触参数是影响种子的流动特性与仿真准确性的重要因素，因此有必要对其接触参数进行标定[17-18]。

国内外学者已对玉米[19]、水稻[20]、小麦[21]、马铃薯[22]、三七[23]等种子的离散元仿真参数进行了标定，结果表明不同种子之间的接触参数具有较大差异，而对绿豆种子离散元仿真参数的标定鲜有报道。因此为应用离散元法优化设计气力式绿豆精密排种器，克服现有绿豆机械式排种器播种精度低、伤种率高等问题，以绿豆种子为研究对象，基于其物理本征参数，采用Hertz Mindlin with bonding粘结模型建立种子仿真模型，标定得到气力式绿豆精密排种器气固耦合仿真所需的离散元仿真参数，并利用台架试验验证其仿真模型与仿真参数的可靠性，从而为气力式绿豆精密排种装置的设计与仿真研究提供参考。

1 模型建立

1.1 绿豆种子物理模型

本文绿豆种子选取白绿9号，该品种具有抗病性强、质量好等优点。选取1 000粒种子平均分成5组，通过试验测得其含水率为(11±0.32)%，千粒质量为(60.14±0.429)g；密度为(1 245.9±1.946)kg/m3；泊松比为0.2；弹性模量为126 MPa[24]；根据公式[25]

(1)

式中G——绿豆种子剪切模量，MPa

E——绿豆种子弹性模量，MPa

μ——绿豆种子泊松比

计算得绿豆的剪切模量为52.5 MPa。

为建立绿豆种子的仿真模型，提高仿真精确度，本文随机选用150粒种子，通过游标卡尺测量其三轴尺寸(长度l、宽度w、厚度t)，并根据公式计算其等效直径dp与球形率φ。

(2)

(3)

测量与计算结果如下：绿豆的平均长度为5.19 mm，平均宽度为3.89 mm，平均厚度为3.93 mm，球形率为82.8%，等效直径为4.29 mm。

依据尺寸分布均值在CATIA中进行三维建模，将建好的三维模型导入EDEM软件，并运用API 颗粒替换法建立仿真颗粒粘结模型(粘结颗粒半径为0.3 mm，数量为207个，单位法向刚度为1×1010N/m3，单位切向刚度为1×1010N/m3，临界法向应力为5×1010Pa，临界切向应力为2×108Pa，颗粒间粘结半径为0.6 mm)，绿豆实物图与模型图如图1所示。

图1 绿豆种子实物及仿真模型Fig.1 Physical model and simulation model of mung bean seeds

1.2 接触模型选取

试验过程中，除颗粒与颗粒间接触，还会有颗粒与其他材料之间的作用力。本文中与绿豆接触的材料选择有机玻璃与Somos8000树脂，其参数如表1所示[26]。

表1 接触材料仿真参数Tab.1 Simulation parameters of contact materials

由于实际试验中绿豆种子与材料表面光滑且几乎无粘附力，因此EDEM仿真时选用Hertz-Mindlin 无滑移接触模型。

2 种子与材料间接触参数标定

2.1 种子与材料间碰撞恢复系数

如图2a所示，采用自由落体碰撞法进行种子与材料间的碰撞恢复系数标定，试验时，在距碰撞材料板H=150 mm的高度将种子释放，种子碰到材料板进行反弹，通过高速摄像系统测定最高反弹高度h。碰撞恢复系数ex可表示为种子与材料碰撞前后在碰撞接触点处的法向瞬时分离速度v1与瞬时接触速度v0之比，即种子与材料碰撞反弹最大高度h与初始下落高度H之比，计算公式为

图2 种子-材料碰撞恢复系数标定试验Fig.2 Calibration test of seed-material coefficient of restitution

(4)

试验重复5次，并通过式(4)计算恢复系数，结果取平均值。实际试验得到种子与有机玻璃碰撞最高反弹高度均值为27.42 mm，恢复系数为0.428；与Somos8000树脂碰撞最高反弹高度均值为24.26 mm，恢复系数为0.402。

仿真试验时，把除碰撞恢复系数以外的其他接触参数设为0。以碰撞恢复系数ex为因素，以最高反弹高度h为指标，选取绿豆与接触材料的碰撞恢复系数为0.1～0.9，间隔均取0.1，每组试验重复5次取均值。将试验结果绘制成散点图并拟合曲线(图3)。

图3 碰撞恢复系数与最高反弹高度拟合曲线Fig.3 Fitting curves of restitution coefficient and maximum springing height

拟合方程分别为

(5)

(6)

式中h1、h2——绿豆与有机玻璃、Somos8000树脂碰撞最高反弹高度

ex1、ex2——绿豆与有机玻璃、Somos8000树脂碰撞恢复系数

将实际试验所得的种子与材料间反弹高度均值代入式(5)、(6)，得ex1为0.445、ex2为0.434，将其分别代入EDEM仿真验证，仿真时最高反弹高度与实测最高反弹高度相对误差较小，因此将ex1定为0.445、ex2定为0.434。

2.2 种子与材料间静摩擦因数

种子与材料间静摩擦因数标定采用斜面滑动法，即质量为m的物体在角度为α的斜面上静止时，其重力可分解为平行于斜面的力F1和垂直于斜面的力F2，F1随着α增大而增大，当α大于物体滑动临界角时，则F1大于物体与斜面间的静摩擦力f，物体开始向下滑动。其中静摩擦因数计算式为

μf=f/F2=mgsinα/(mgcosα)=tanα

(7)

实际试验时将4粒绿豆种子粘结在一起(为防止种子滚动)放置于图4a所示接触材料的斜面左端，缓慢匀速提升放置种子的斜面端，当种子开始滑动时，记录斜面与水平面之间的夹角α，试验重复5次，结果取平均值。得到绿豆与有机玻璃、Somos8000树脂的滑动摩擦角均值分别为23.9°、30.8°，静摩擦因数分别为0.443、0.596。

图4 种子-材料静摩擦因数标定试验Fig.4 Calibration test of seed-material coefficient of static friction

仿真时在EDEM中以种子与材料间静摩擦因数μf为因素，斜面倾角α为指标进行试验(采用已标定完的碰撞恢复系数，其他接触参数设为0，选取静摩擦因数为0.1～0.9，试验水平间隔为0.1)。为避免种子在斜面上滚动，将4个粘结的颗粒替换模型放置在斜面一端，设置与实际试验相同的转速进行仿真，当种子开始滑动时，记录斜面倾角α，试验重复5次，结果取平均值。将试验结果绘制成散点图并拟合曲线(图5)。

图5 静摩擦因数与斜面倾角拟合曲线Fig.5 Fitting curves of static friction coefficient and inclination angle

拟合方程分别为

(8)

(9)

式中μf1、μf2——绿豆与有机玻璃、Somos8000树脂的静摩擦因数

α1、α2——有机玻璃板、Somos8000树脂板与水平面的夹角

把实际试验所测倾角分别代入式(8)、(9)，得μf1为0.458、μf2为0.556，将其分别代入EDEM验证，得到仿真滑动摩擦角与实际滑动摩擦角基本一致，因此将μf1定为0.458、μf2定为0.556。

2.3 种子与材料间滚动摩擦因数

如图6所示，种子与材料间滚动摩擦因数测定采用斜面滚动法。即在倾角为β的斜面板上以零初速度释放种子，种子沿斜面向下滚动一段距离S，由于受到滚动摩擦，种子最终滚落至水平面板上并静止，设种子水平滚动距离为L，假设种子为理想球体，在纯滚动过程中只受到滚动摩擦力影响，则通过能量守恒定律可得

图6 种子-材料滚动摩擦因数标定试验Fig.6 Calibration test of seed-material coefficient of rolling friction

mgSsinβ=mg(Scosβ+L)μs

(10)

式中μs——滚动摩擦因数

由于绿豆种子并非理想球体，当倾角β与斜面滚动距离S较大时，会导致种子在滚动过程中发生弹跳，影响试验结果准确性，当β与S较小时，种子滚动距离较小，不利于试验测量，因此经大量预试验调整，设定倾角β为20°、斜面滚动距离S为30 mm，以零初速度释放种子，使种子沿斜面向下滚动，待种子静止，测量种子水平滚动距离L。为减小试验误差，试验重复30次，结果取平均值。试验测得绿豆在有机玻璃、Somos8000树脂上的滚动距离均值分别为193、141 mm，滚动摩擦因数分别为0.046、0.061。

仿真试验采用相同的方法，EDEM中设置已标定完毕的碰撞恢复系数与静摩擦因数，其它接触参数设为0，以滚动摩擦因数μs为因素，选取因数范围0.01～0.09，试验水平间隔为0.01，以水平滚动距离L为评价指标，将试验结果绘制成散点图并拟合曲线(图7)。

图7 滚动摩擦因数与水平滚动距离拟合曲线Fig.7 Fitting curves of rolling friction coefficient and horizontal rolling distance

拟合方程分别为

(11)

(12)

式中μs1、μs2——绿豆与有机玻璃、Somos8000树脂的滚动摩擦因数

L1、L2——种子在有机玻璃板、Somos8000树脂板上的水平滚动距离

将实际测得绿豆在有机玻璃、Somos8000树脂上的水平滚动距离均值分别代入式(11)、(12)，得μs1为0.036、μs2为0.049，将其分别代入EDEM验证，得到仿真水平滚动距离与实际水平滚动距离相对误差较小，因此μs1为0.036、μs2为0.049。

3 种间接触参数标定

采用上述试验中自由落体碰撞法、斜面滑动法、斜面滚动法对种子之间的接触参数标定时，需将种子粘结成种子板，由于粘结的底板材料不同、种子板表面凹凸不平，对试验结果会产生较大影响。因此本文采用堆积角试验方法，以种间接触参数为因素，以实测堆积角与仿真堆积角的相对误差为指标，进行最陡爬坡试验、三因素五水平旋转组合设计试验，对试验结果优化得到种子之间的仿真接触参数。

3.1 实际堆积角试验

采用传统抽板法进行试验时，由于底部材料的不同，会导致种子堆积角产生较大差异，因此底面设置一直径为140 mm、高度为30 mm的圆筒，圆筒中先填满种子并刮平，距离圆筒上端80 mm处设置一漏斗，用挡板挡住漏斗下端面，漏斗中装适量种子，迅速抽掉挡板，在圆筒上端形成堆积锥角，装置材料分别采用有机玻璃、Somos8000树脂，厚度均为2 mm，经预试验发现两材料对堆积角影响无显著差别，因此选用有机玻璃进行堆积角试验。实际试验与仿真试验如图8所示。

图8 堆积角试验Fig.8 Stacking angle test

试验结束后用Matlab软件读取颗粒堆边缘图像(图9a)，对图像进行灰度化、二值化处理并提取二值化图像边界轮廓(图9b)，扫描轮廓边缘点并进行线性拟合，拟合直线(图9c)与水平面夹角即为实测堆积角θ。

图9 Matlab软件图像处理Fig.9 Image processed by Matlab

每组试验重复5次，结果取平均值，实际测得绿豆种子堆积角θ均值为22.1°。

3.2 最陡爬坡试验

为确定三因素五水平二次正交旋转组合设计试验的零水平及因素水平的最优值区间，进行最陡爬坡试验。在EDEM中设置已标定的种子与有机玻璃间的接触参数，以颗粒间碰撞恢复系数e、静摩擦因数μn、滚动摩擦因数μr为因素，以实测堆积角θ与仿真堆积角θ′相对误差σ为指标进行试验。

经查阅文献[1,14-16,21]发现，大部分粮谷类种子之间的碰撞恢复系数在0.15～0.85之间，静摩擦因数在0.1～0.6之间，滚动摩擦因数在0.01～0.08之间，因此方案设计与试验结果如表2所示。

表2 最陡爬坡试验方案设计与结果Tab.2 Scheme and results of the steepest ascent experiment

由表2可知，采用第2组组合试验仿真得到堆积角与实测堆积角相对误差最小，为1.81%，因此分别选用第1组、第2组、第3组作为三因素五水平旋转组合设计试验编码值。

3.3 二次正交旋转组合设计试验

以颗粒间碰撞恢复系数、静摩擦因数、滚动摩擦因数为试验因素，以实测堆积角θ与仿真堆积角θ′相对误差Y为指标，进行三因素五水平二次正交旋转组合设计试验，并利用Design-Expert软件对试验数据进行回归分析，利用响应面法寻找最佳参数组合。其因素编码如表3所示。试验方案与结果如表4所示，A、B、C为因素编码值。

表3 试验因素编码Tab.3 Experiment factors and codes

表4 试验方案与结果Tab.4 Experiment scheme and results

各参数对堆积角相对误差的方差分析如表5所示。由表5可知，3个试验因素中，颗粒间静摩擦因数、滚动摩擦因数对试验指标影响均极显著，而颗粒间碰撞恢复系数对试验指标无显著影响，分析其原因可能是：①试验前为排除底部碰撞材料的影响，底部圆筒中先充满颗粒，而颗粒形成的表面凹凸不平，试验时颗粒碰撞在颗粒群表面上更易受其形状影响，限制了其碰撞后反弹运动。②由于颗粒群内部孔隙率较大，且呈现松散状态，相比密度均匀无孔隙的材料板，颗粒群的吸能性更好，颗粒碰撞在颗粒群上能量耗散会更快，颗粒碰撞后运动速度急剧下降。

表5 堆积角相对误差的方差分析Tab.5 Analysis of variance for the relative error in stacking angle

去除对指标影响的不显著项后，通过Design-Expert软件得拟合回归方程为

Y=6.37-5.98B-2.94C+1.95AB+
4.12A2+4.85B2+2.26C2

(13)

由表5可知，该回归模型极显著，且模型的失拟项不显著，说明拟合的回归方程具有高度可靠性，能够准确反映各显著因素与堆积角相对误差之间的关系。

根据上述试验结果及回归方程，将非显著因素——颗粒间碰撞恢复系数设为中间水平，以相对误差最小为目标，寻找颗粒间静摩擦因数与滚动摩擦因数的最优组合，其设定目标函数及约束条件为

(14)

最终得到当颗粒间碰撞恢复系数为0.3、静摩擦因数为0.23、滚动摩擦因数为0.03时，堆积角实测值与仿真值相对误差为3.91%。将上述参数代入EDEM仿真验证后，结果表明与实际试验误差较小，因此将上述3个参数值作为其初步标定结果。

4 验证试验

4.1 试验装置

为验证上述标定的仿真参数能够应用于排种器的仿真模拟中，进行验证试验，实际试验采用文献[2]所设计的自吸式绿豆精密排种器(种腔为有机玻璃，种盘为Somos8000树脂)，选取含水率为11%的白绿9号种子。试验台选择吉林大学农机实验室的JPS-12型排种器试验台(图10)。

图10 排种器试验台Fig.10 Test bed of seed-metering device1.真空计 2.台架 3.自吸式排种器 4.驱动电机 5.高速摄像机 6.摄像机控制系统 7.试验台数据采集显示器 8.试验控制台

由于该排种器利用机械结构产生负气压进行吸种，因此仿真试验采用EDEM-Fluent流固耦合方式进行仿真，为提高流体域网格划分质量与耦合仿真速度，将排种器结构简化，省略种腔内部不规则结构，与实际种腔(图11a)相对比，简化模型只保留种盘、种腔、负压吸气口3部分，所建立的排种器简化模型如图11b所示。仿真时当种子被吸附在型孔中转过吸气口掉落时，仍可回到种腔继续应用于仿真，目的是为了生成少量的种子而提高仿真速度。

图11 排种器模型Fig.11 Model of seed-metering device

利用ICEM软件对图11b对应的流体域进行网格划分，其中种腔与负压吸气口对应流体域采用结构化网格划分，型孔流体域采用非结构四面体网格，将种腔流体一侧壁面设为inlet，负压吸气口流体一侧壁面设为outlet，种腔流体与型孔流体接触面、型孔流体与负压吸气口流体接触面均设置interface面，结果如图12所示。

图12 流体域网格划分图Fig.12 Diagram of grid of fluid area

4.2 试验方法

实际试验以排种盘转速为因素，以漏播率、重播率、合格率为指标，进行单因素试验，确定排种轴转速分别为90、100、110、120、130 r/min，其对应的排种盘转速分别为27.8、30.9、34、37、40.1 r/min。每组试验重复3次，结果取均值。

仿真试验将图11b排种器简化模型保存为stp格式导入EDEM，种腔材料设置为有机玻璃，种盘设置为树脂，其相关材料及接触参数参照上述标定结果进行设置，随后在种腔内生成120粒颗粒模型，并保存成0 s时刻的文件。颗粒替换成功，打开保存好的0 s时dem文件，并启动耦合按钮Show Coupling Server，将划分的网格文件导入Fluent 17.0软件中，根据不同的种盘转速及负压进行相关参数设置。仿真试验以种盘转速为因素(排种盘转速分别为27.8、30.9、34、37、40.1 r/min时，利用真空计测得相对应的实际负压分别为-1.84、-2.62、-3.36、-4.15、-4.93 kPa)，以漏吸率(未充填种子的型孔数除以120)、重吸率(充填两粒及以上的型孔数除以120)、单粒率(充填一粒种子的型孔数除以120)为指标(由于仿真速度与时长的限制，本次试验只统计120个型孔)，进行单因素试验，每组试验重复3次，结果取均值。

4.3 试验结果

如图13所示，截取仿真过程中的一些时刻，0.39 s时3个型孔已经分别成功充填3粒种子，无产生漏吸与重复吸种现象，仿真进行到0.64 s时开始出现漏吸，仿真0.98 s时，种子已经转过负压区产生掉落现象，1.25 s时，型孔中已经产生种子重复吸种的现象。

图13 仿真过程图Fig.13 Diagram of simulation process

仿真与台架试验对比结果如图14所示。

图14 仿真与台架试验对比结果Fig.14 Comparison result of simulation test and bench test

从图14中可以看出，随着种盘转速的提升，仿真试验与台架试验的3个指标变化趋势基本一致。漏吸率、漏播率均随着种盘转速加快而下降，两者最大相对误差为4.71%；重吸率、重播率均随着种盘转速加快而上升，最大相对误差为4.94%；单粒率、合格率均随着种盘转速加快而下降，最大相对误差为0.98%。台架试验的性能指标总体上均劣于仿真试验，由于仿真试验比较理想，台架试验存在振动、种子在种管中碰撞等因素，因此存在一定差距是合理的，但性能指标总体变化趋势一致，且误差不大于5%，说明标定试验结果可靠。

5 结论

(1)运用离散元仿真软件EDEM，分别采用自由落体碰撞法、斜面滑动法、斜面滚动法对绿豆种子与接触材料(有机玻璃、Somos8000树脂)间的碰撞恢复系数、静摩擦因数、滚动摩擦因数进行了标定。绿豆与有机玻璃碰撞恢复系数、静摩擦因数、滚动摩擦因数分别为0.445、0.458、0.036，绿豆与Somos8000树脂碰撞恢复系数、静摩擦因数、滚动摩擦因数分别为0.434、0.556、0.049。

(2)采用堆积角试验方法，以种间接触参数为因素，以实测堆积角与仿真堆积角的相对误差为指标，进行最陡爬坡试验、三因素五水平旋转组合设计试验，以实测堆积角与仿真堆积角相对误差最小为约束条件，取非显著因素——颗粒间碰撞恢复系数为0.3，对试验结果寻优得到绿豆种间静摩擦因数、滚动摩擦因数分别为0.23、0.03。

(3)选用自吸式绿豆精密排种器为试验装置，以种盘转速为试验因素，台架试验以漏播率、重播率、合格率为指标，仿真试验以漏吸率、重吸率、单粒率为指标，在EDEM中设置已标定的仿真参数，进行单因素对比验证试验，结果表明：仿真试验漏吸率与台架试验漏播率最大相对误差为4.71%，重吸率与重播率最大相对误差为4.94%，单粒率与合格率最大相对误差为0.98%，性能指标总体变化趋势一致，且误差均不大于5%，说明标定的仿真参数可用于离散元仿真试验。