基于EDEM的稻种离散元摩擦参数标定

刘大为，段佳鹏，陈星宇，熊蒸科，王修善，谢方平

（1.湖南农业大学机电工程学院，长沙 410128；2.智能农机装备湖南省重点实验室，长沙 410128）

杂交水稻是我国一项重要的原创性粮食增产技术，为保障中国乃至世界粮食安全发挥了重要作用[1-4]。分植制种是我国培育杂交水稻常用的制种方法，即父母本分开播种和收获。其过程对劳动力需求量大，不利于杂交稻工程化制种机械化和轻简化发展。杂交水稻父、母本混播混收是降低制种成本的重要技术，其核心是父、母本种子的分离[5-6]。因此，研究杂交水稻父、母本种子在分选机械中的运动规律，建立较为精确的杂交水稻父、母本颗粒模型，对设计和优化杂交水稻种子分选机械，改善分选机械的作业性能尤为重要。

当前，多数研究者借助离散元法（discrete element method，DEM）开展收获、清选等农业工程领域的研究[7-10]。科学地运用DEM 方法需要准确定义仿真模型的物性参数，包括物料与材料的本征参数（如密度、泊松比、弹性模量等参数值）和物料与物料、物料与材料的接触参数（如静摩擦因数、滚动摩擦因数、恢复系数等参数值）。

在进行DEM 仿真时，物料和材料本征参数可以从已有研究中获取[11-12]，但是，由于建立的颗粒模型在形状上与实际情况存在一定的差异，使得接触参数与实测值存在误差，因此，需要对接触参数进行标定[13]。国内外已经进行了基于DEM 的相关仿真参数标定研究。王云霞等[13]通过建立一种能够主动寻找目标参数的数学回归模型，标定了玉米种子颗粒模型静摩擦因数和滚动摩擦因数；王韦韦等[14]通过Plackett-Burman 试验和Box-Behnken试验标定了玉米秸秆粉料致密成型离散元模型的秸秆粉料-粉料滚动摩擦因数和秸秆粉料-不锈钢板静摩擦因数；侯占峰等[15]通过物理试验和仿真试验相结合的方法标定了冰草种子的离散元仿真参数；BRYAN等[16]通过环剪试验和直剪试验标定了作物残体的摩擦系数、颗粒法向刚度和剪切刚度；CHEN 等[17]采用直接测量和校准试验相结合的方法，确定了玉米和小麦籽粒的离散元接触参数；GHODKI等[18]对实测结果和仿真模拟结果进行定性和定量分析，标定了大豆颗粒间的恢复系数和静摩擦因数等。综合国内外离散元参数标定的发展现状，鲜有关于杂交水稻种子离散元接触参数标定的相关文献报道。

由于颗粒间的静摩擦因数、滚动摩擦因数对颗粒的运动影响较大[19]，本研究以杂交稻（玉针香）父本种子和母本种子为对象，通过EDEM 分别建立颗粒模型，进行堆积角仿真试验，对堆积角实测试验和仿真试验结果进行对比分析，分别标定父、母本种子的种间静摩擦因数和滚动摩擦因数。首先，通过实测试验分别得到杂交水稻父、母本堆积角并采集图像，利用MATLAB 图像处理技术对种子堆积群外轮廓进行线性拟合处理，得到堆积角的实测值；其次，通过EDEM 软件，建立杂交水稻父、母本种子颗粒模型，以种间静摩擦因数和滚动摩擦因数为试验因素进行堆积角的全因素仿真试验；再次，将仿真试验结果与堆积角实测值进行对比，取相对误差最小组的种间滚动摩擦因数为标定值，再通过拉格朗日插值法确定种间静摩擦因数；最后以标定后的种间参数进行堆积角试验，验证该参数的可靠性。

1 材料与方法

1.1 堆积角实测试验

堆积角试验装置由万能试验机、不锈钢圆筒和不锈钢板等组成（图1）。不锈钢圆筒的内径和填充高度分别为45 mm 和135 mm，测量时，将不锈钢圆筒底部贴合不锈钢板，填充种子至指定高度后，放置在万能试验机上，通过计算机控制提升装置带动圆筒匀速提升，种子自然滑落并堆积，待种子堆积群稳定后，采集堆积角的图像，其斜面与水平面的夹角即为堆积角。

图1 圆筒提升试验装置Figure 1 Cylinder lifting test device

为了减少人工测量所带来的误差，通过MTALAB软件先后对采集到的堆积角图像进行二值化、边界提取、线性拟合等处理（图2），取左右两侧拟合直线与水平面夹角的平均值作为堆积角的实测值，重复4 次试验，得到父、母本种子堆积角的实测数平均值分别为36.07°和29.06°，标准差分别为0.50 和1.76，极差分别为0.94和4.28。考虑锥形堆表面处的平衡条件，结合库仑屈服条件，内摩擦角的正切应等于斜面倾角的正切，求得父、母本种间静摩擦因数值取值范围分别为0.51～0.74和0.43～0.62。

图2 父、母本种子堆积角图像处理Figure 2 Image processing of angle of repose for male and female seed

1.2 堆积角仿真试验

1.2.1 种子颗粒模型的建立 从杂交水稻父、母本种子中随机选取60 粒，其中父本、母本种子各30 粒，采用精度为0.02 mm 的游标卡尺分别对父、母本种子的长度、宽度和厚度，进行测量（图3），得到父本种子和母本种子的长、宽、厚的平均尺寸分别为10.76，2.97，2.02 mm和11.03，2.26，2.05 mm，因此，创建杂交水稻种子颗粒模型时，取父本和母本种子颗粒长、宽、高3个方向的尺寸分别为10.8，3.0，2.0 mm和11.0，2.3，2.1 mm。

图3 父、母本种子颗粒Figure 3 Male and female seed granules

在EDEM 软件中采用基础颗粒进行多球填充，建立杂交水稻父、母本种子颗粒的多球模型。本研究采取双列排布的方式，在满足宽度尺寸的前提下，所采用球体的直径均不大于种子厚度尺寸[19]。首先在AutoCAD软件中根据种子外形和尺寸绘制杂交水稻父、母种子多球模型的二维图，并记录各个球体单元的球心坐标和半径；然后在EDEM软件中通过手动输入球体球心坐标和半径值的方式，得到杂交水稻父、母本种子颗粒模型（图4），父、母本种子多球模型左右对称，分别由16和12个球体单元堆叠形成。

图4 父、母本种子多球模型Figure 4 Multi-bulb model of male and female seed

1.2.2 几何体建模 为了减少仿真试验的误差，圆筒的规格需与实测试验一致，即圆筒内径为45 mm，填充高度为135 mm。在EDEM 软件Geometry 中直接添加Cylinder（圆筒）和Polygon（底板），Cylinder 半径、高度分别为22.5 mm 和200 mm，Polygon 尺寸为500 mm×500 mm，材料均为不锈钢，并调整圆筒和底板的位置，使底板平面垂直于圆筒轴线。

1.2.3 接触模型的选取 EDEM 常用的离散元接触模型主要包括Hertz-Mindlin (No Slip)、Hertz-Mindlin with Bonding、Hertz-Mindlin with JKR 等。由于收获后的稻谷颗粒间的粘附力可以忽略，本次研究选择Hertz-Mindlin(No Slip)模型[21]。

1.2.4 杂交稻种子仿真参数的设置 稻种间滚动摩擦因数值通过在PVC 板上分别粘上单层父、母本种子，并以此做斜面，抬升斜面，测量单粒稻种发生滚动时的斜面倾角，求取相应的滚动摩擦因数值；静摩擦因数值取值范围依据前述堆积角实测试验得到；材料的本征参数和接触恢复系数等相关试验参数值参考文献[19]，父、母本种子及圆筒的仿真参数如表1和表2。由于EDEM 软件仿真时间步长对仿真效果的影响较大，仿真时间步长过大，物料群整体会产生旋涡状的运动状态和爆破现象，仿真时间过小，仿真结果虽然越接近理论值，但是会产生很大的计算量，极大地降低了仿真效率[19-22]。因此，本次仿真试验的仿真时间步长设置为Rayleigh 时间步长的25%，网格尺寸设置为最小颗粒半径的3倍，保存时间步长设置为0.01s。

表1 材料的本征参数Table 1 Intrinsic parameters of materials

表2 材料的接触参数Table 2 Contact parameters of materials

在正式试验前，通过预试验确定本次仿真试验的因素和水平（表3）。进行仿真试验时，待几何体生成后，创建颗粒工厂生成颗粒填满圆筒至135 mm 高度；设定圆筒提升速度为0.05 m·s-1，运动时间1.6 s[23]。颗粒将从圆筒底部逐渐流出，在底板上形成堆积群，整个过程如图5，待堆积群稳定后，采集堆积群在XOY 和XOZ 两个平面的投影图像，通过MATLAB 进行图像处理，得到堆积角的角度，重复3次试验，取3次试验结果的平均值作为堆积角的仿真值。将堆积角仿真值与实测值进行对比，计算出相对误差：

表3 仿真试验因素水平Table 3 Simulation test factor and level table

图5 堆积角形成过程示意图Figure 5 Schematic diagram of formation process of repose angle

2 结果与分析

2.1 堆积角仿真试验

根据表3 所安排的试验因素水平进行3×3 全因素试验，试验结果如表4 和表5。显然，杂交水稻父、母本种子堆积角的目标值已包含在试验结果中，表明本次试验因素水平安排合理。同时，父本种子堆积角仿真结果中P2 的试验结果相对误差最小，母本种子堆积角仿真结果中Q2 的试验结果相对误差最小，由此可以确定P2 和Q2 组试验的仿真参数值最接近目标值，取该两组的种间滚动摩擦因数水平作为标定值。通过仿真结果可以看出，种间滚动摩擦因数一定的情况下，堆积角随着种间静摩擦因数的增大呈现增大的趋势。在工程应用中，处理离散数据时，通常采用插值法进行预测[25]。对给定的n个插值点x1,x2,...,xn及对应的函数值y1,y2,...,yn，利用Lagrange 插值多项式公式，则对插值区间内任意x的函数值通过式（2）求得：

表5 堆积角全因素仿真试验结果（母本种子）Table 5 The results of all-factor simulation test of repose angle(female seeds)

因此，本研究通过Matlab编程，采用拉格朗日插值法分别对P2和Q2进行插值，得到父、母本种子的种间静摩擦因数分别为0.575和0.453。

综上，父本种子间的滚动摩擦因数和静摩擦因数分别为0.11 和0.575；母本种子间的滚动摩擦因数和静摩擦因数0.07和0.453。

2.2 验证试验

为了验证本次试验所标定的仿真参数的准确性，对已标定参数进行堆积角仿真验证分析。分别用标定后的杂交水稻父、母本种间滚动摩擦因数和静摩擦因数进行仿真试验，重复3 次，得到父本种子的堆积角为36.31°，相对误差为0.67%，母本种子的堆积角为29.19°，相对误差为0.45%（表6）。可见，仿真结果与实测值无显著差异，本研究可以为杂交水稻种子分选仿真过程提供较为准确的颗粒模型和种间参数。

表6 验证试验结果Table 6 Validation test results

3 讨论与结论

混播制种是当前公认的提高杂交水稻制种效率的有效途径之一，种间滚动摩擦因数和静摩擦因数是影响杂交水稻种子在分选机械内运动的主要因素，利用离散元法开展杂交水稻种子分选影响因素研究具有现实意义。稻种间接触参数研究的主要方法是通过圆筒提升和方盒抽板两种方式形成堆积角，并建立二元回归方程组求解滚动摩擦因数和静摩擦因数两个目标参数，由于离散元仿真试验加载时间往往很长，相对复杂的试验安排会增加试验次数，时间成本高，在一定程度上降低了研发效率。目前，已有不少学者利用离散元法研究物料在清选装置的运动情况，但鲜有见到杂交水稻种子分选机械的仿真研究。

本研究针对杂交水稻（玉针香）父、母本种子，结合堆积角物理实验和仿真试验，通过对比分析和拉格朗日插值，标定了杂交水稻父、母本种子颗粒模型的种间滚动摩擦因数和静摩擦因数，其结果分别为0.11 和0.575，0.07 和0.453。以标定后的种间参数进行堆积角验证试验，试验结果显示，与物理试验相比，父、母本堆积角相对误差分别为0.67%和0.45%，表明本次研究建立的杂交水稻父、母本种子颗粒模型及标定的参数较为精准。故而本研究能够为杂交水稻种子分选仿真过程提供较为准确的颗粒模型和种间摩擦参数。此外，考虑到杂交水稻品种众多，外形尺寸和力学性质差异较大，若实测结果与仿真结果相对误差较大，应当重新设计试验求解种间参数。