基于离散元的三七种子仿真参数标定与试验

于庆旭 刘 燕 陈小兵 孙 凯 赖庆辉

(1.农业农村部南京农业机械化研究所， 南京 210014； 2.昆明理工大学农业与食品学院， 昆明 650500)

0 引言

离散元法(Discrete element method，DEM)在农业装备中应用越来越广泛，应用离散元法研究散粒体动力学已成为一种发展趋势[1-3]。排种器在工作过程中，种子与种子间、种子与排种器之间的作用力非常复杂，可基于离散元法对排种器进行研究。三七是中国名贵药材，研究三七种子离散元模型和仿真参数，有助于离散元法在三七排种器研究中的应用。

为了提高离散元仿真试验的精度，在离散元软件中，需精确建立物料的离散元模型和准确定义仿真模型的物性参数[4-6]。物性参数主要包括本征参数和接触参数，本征参数是物体自身的特性参数，与外界因素无关，通常为固定值，如泊松比、剪切模量和密度等，可通过台架试验直接测得；接触参数包括碰撞恢复系数、静摩擦因数和滚动摩擦因数，是两个物体发生接触时的物性参数，与接触的两个物体有关系，可通过虚拟的仿真试验标定获得。

目前，国内外相关学者基于离散元对农业物料进行了大量的研究，已经完成土壤、稻谷、饲料、马铃薯、玉米等离散元模型的建立和仿真参数的标定。LANDRY等[7]建立有机肥料的离散元模型，通过剪切试验标定接触参数。COETZEE等[8]通过剪切和压缩试验分别标定玉米离散元模型的摩擦因数和刚度，通过筒仓卸料和斗装试验进行验证。UCGUL等[9]通过休止角和贯入度试验标定模拟土壤耕作所需的参数，得到的接触模型和接触参数可有效模拟无粘性土壤，预测土的受力和运动。MIYAMOTO等[10]通过休止角试验标定水稻离散元参数，通过分选机进行试验验证。刘彩玲等[11]基于三维激光扫描建立水稻的离散元模型，比常规离散元模型的仿真精度更高。刘文政等[12]在EDEM软件中建立微型薯的离散元模型，通过试验测定及仿真模拟相结合的方法，对微型薯颗粒离散元参数进行标定和校准。刘凡一等[13]基于响应面优化，以真实试验及不同参数组合下仿真得到的小麦颗粒堆休止角为响应值，标定了小麦离散元仿真参数。向伟等[14]构建土壤离散元仿真模型，基于土壤堆积试验，结合试验测定和EDEM 软件推荐的参数构建土壤仿真模型，以休止角为响应值，采用Design-Expert软件依次设计Plackett-Burman试验、最陡爬坡试验和Box-Behnken试验，完成土壤仿真物理参数标定及优化，通过成穴装置成穴的仿真试验与土槽试验的对比分析，验证黏壤土仿真模型的精准性。彭飞等[15]建立与颗粒饲料形态相近的离散元模型，基于注入截面法的休止角测定装置与方法，进行休止角的模拟与测定，标定了颗粒饲料离散元模型参数。鹿芳媛等[16]基于水稻芽种摩擦角试验与仿真测定，标定出不同含水率下水稻芽种离散元的主要接触参数。王云霞等[17]结合仿真试验和实际试验，对玉米种子颗粒模型的种间静摩擦因数和滚动摩擦因数两个关键参数进行标定，为了简化标定过程，建立数学回归模型主动寻找目标参数，使其在EDEM 中建立的玉米种子颗粒重新获得与真实颗粒相近的物理特性。

本文采用逆向工程技术，基于粘结颗粒模型，在EDEM软件中建立三七种子离散元模型；结合台架试验和仿真试验，标定三七种子与ABS塑料之间接触参数；通过堆积试验，基于二次回归正交旋转组合试验的响应面优化方法，确定EDEM仿真试验中最佳的三七种子之间接触参数；通过三七精密排种器进行验证试验，建立三七种子离散元模型，以期得到最佳的三七种子仿真参数。

1 三七种子离散元模型

1.1 轮廓模型

三七主要产自于云南省文山市，选用3年生成熟的文山三七种子，简称三七种子，三七种子含水率20%～60%，密度为929～1 132 kg/m3，长度5.2～7.2 mm，宽度为4.8～6.8 mm，高度为4.0～6.0 mm，平均直径为5.62 mm[18]，选取长宽高尺寸与平均值相近的三七种子建立轮廓模型，三七种子实物如图1a所示。

三七种子为非规则形状颗粒，为了精确建立轮廓模型，应用逆向工程技术，通过EinScan-S 3D扫描仪三维扫描三七种子外轮廓，采用Geomagic Wrap 3D软件处理三维扫描数据，得到种子点云数据，如图1b所示。

将点云数据导入UG软件中转换为多边形，进行裁剪多余曲面、删除钉状物、合并和光顺等操作，利用曲面逼近方法，对三七种子多边形模型进行优化处理，最终得到三七种子轮廓模型，如图1c所示。

图1 三七种子轮廓模型

1.2 离散元模型

1.2.1颗粒模型

在离散元仿真软件中，通常选用球形颗粒建立物料的离散元模型，球形颗粒的形状简单、规则，只有一个半径尺寸参数，而且球形颗粒之间只有一种接触状态，所以球形颗粒间的接触状态易于检测，可大大缩短检测计算的时间，有利于程序稳定运行。规则的类似球形种子可以用单个球形颗粒简化，图2a为单个球形离散元模型。

图2 颗粒模型

实际种子大多为不规则体，不规则种子之间的咬合和翻滚等特性，无法直接用单个球形离散元模型进行仿真模拟，为了真实仿真种子之间的特性，可以采用多球聚合模型[19-22](Multi-sphere method，MSM)和粘结颗粒模型[23-25](Bonded particle method，BPM)建立种子离散元模型。

多球聚合颗粒模型如图2b所示，通过若干个直径不等的球形颗粒重叠堆积完成，通常球形颗粒半径较大。该方法得到的种子离散元模型与种子实际外轮廓较为吻合，由于球形颗粒重叠堆积，在 EDEM 软件中被判定为一个独立体。

多球粘结颗粒模型如图2c所示，由多个直径相等的球形颗粒通过“粘结键”粘结而成，当球形颗粒半径越小时，越近似种子实际轮廓。球形颗粒彼此间是独立的个体，通常球形颗粒半径较小，通过多球粘结颗粒模型建立的种子离散元模型，适用于EDEM-CFD耦合仿真[26]。

1.2.2颗粒填充

选用多球粘结颗粒模型建立三七种子离散元模型，采用自动填充的方式，当球形颗粒直径越小时，粘结的球形颗粒数就越多。将三七种子轮廓模型导入到EDEM软件中，选择三七种子轮廓模型为球形颗粒生成工厂，球形颗粒接触模型为Hertz-Mindin无滑移模型。为了使球形颗粒填充得更加致密，设置较小的球形颗粒之间碰撞恢复系数、静摩擦因数和动摩擦因数，可加大球形颗粒填充速度。球形颗粒填充如图3所示。

图3 球形颗粒填充

如图3a所示，填充的球形颗粒半径为0.1 mm，球形颗粒数量为16 383；如图3b所示，填充的球形颗粒半径为0.3 mm，球形颗粒数量为1 582；如图3c所示，填充的球形颗粒半径为0.5 mm，球形颗粒数量为98。当球形颗粒直径越小时，越近似种子实际轮廓，但球形颗粒的数量越多，仿真计算时间越长，综合考虑球形颗粒半径选取0.5 mm。

填充完成后，在 EDEM 软件后处理界面中，导出98颗球形颗粒的半径尺寸、ID编号和中心坐标参数。将EDEM仿真时间归零，在后处理界面中仅选择颗粒数据，导出Simulation Deck文件，文件格式为XML，将98颗球形颗粒的半径尺寸、ID编号和中心坐标参数编辑保存到该文件中，得到新的XML格式文件。在新建的EDEM中导入新的XML格式文件，在EDEM软件的Particles项中，自动填充生成三七种子离散元模型，后续的EDEM仿真试验可直接调用，三七种子离散元模型如图4所示。

图4 三七种子离散元模型

2 接触参数标定

由于不规则种子的各向异性，通过试验直接测量接触参数的真实值，该数据变化差异较大，同时有些参数很难通过试验直接测量；而且不规则种子的离散元模型与实际尺寸存在差异，仿真接触参数和真实接触参数存在误差，所以通过试验测量的真实接触参数不直接应用到仿真试验[12]。为了提高离散元仿真试验的可靠性，本文结合台架试验和仿真试验，标定仿真试验中的接触参数。前期设计采用3D打印，3D打印材料为ABS塑料，本文对三七种子与ABS塑料之间和三七种子之间的接触参数进行标定。常用的接触参数标定试验有碰撞弹跳试验、斜面滑移试验、斜面滚动试验和堆积试验等[27-28]。标定种子和材料的接触参数时，种子含水率对接触参数有重要的影响，选取含水率为45%～55%的三七种子。在EDEM软件中进行仿真试验，参数如表1所示。

表1 仿真试验参数

2.1 三七种子与ABS塑料之间接触参数标定

2.1.1碰撞恢复系数

图5 碰撞恢复系数标定试验

碰撞恢复系数是物体在接触碰撞变形后恢复能力的参数，采用碰撞弹跳试验标定三七种子与ABS塑料之间的碰撞恢复系数，如图5所示，将ABS塑料板水平放置，将三七种子从高度H1=100 mm处自由下落，与ABS塑料板发生碰撞，通过高速摄像记录三七种子最高弹起高度，台架试验重复5次，计算平均值，得最高弹起高度实测值h1=19.62 mm。

三七种子与ABS塑料之间静摩擦因数x2和滚动摩擦因数x3，以及三七种子之间碰撞恢复系数X1、静摩擦因数X2和滚动摩擦因数X3，对弹跳高度没有影响，为了避免干扰，在EDEM仿真试验中，x2、x3、X1、X2和X3的值均为0，经过预仿真试验，取三七种子与ABS塑料之间碰撞恢复系数x1范围为0.50～0.80，取步长0.05，进行7组仿真试验，每组试验重复5次取平均值，试验设计方案与结果如表2所示，表中y1为最高弹起高度仿真值。

表2 碰撞恢复系数仿真试验设计方案与结果

为得到仿真试验中三七种子与ABS塑料之间碰撞恢复系数与最大弹起高度的关系，对表2中的试验数据进行曲线拟合，得到2次多项式拟合曲线如图6所示，曲线方程为

(1)

图6 仿真试验碰撞恢复系数与最大弹起高度拟合曲线

式(1)的决定系数R2=0.988 4，接近1，表明该方程拟合的可靠度高。将台架试验的最大弹起高度实测值19.62 mm代入式(1)，求得x1=0.611，输入EDEM软件中进行验证，试验重复5次，计算平均值，测得最大弹起高度仿真值为20.27 mm，与实测值的相对误差为3.29%，由此表明标定后的仿真结果与台架试验基本一致。确定EDEM仿真试验中，三七种子与ABS塑料之间碰撞恢复系数x1=0.611。

2.1.2静摩擦因数

静摩擦因数是物体所受的最大静摩擦力与法向压力的比值，可以采用斜面法和拖重法进行测试。使用斜面滑移试验标定三七种子与ABS塑料间的静摩擦因数，试验如图7所示。试验过程中为防止三七种子滚动，降低试验误差，将4颗三七种子粘结成正方形。初始时ABS塑料板上水平放置，将粘结种子放置在ABS塑料板上，使ABS塑料板绕着一边缓慢匀速转动，粘结种子开始滑移时停止转动，记录倾斜板倾斜角，台架试验重复5次，计算平均值，求得倾斜板倾斜角实测值γ1=26.60°。

图7 静摩擦因数标定试验

三七种子与ABS塑料之间滚动摩擦因数x3，以及三七种子之间碰撞恢复系数X1、静摩擦因数X2和滚动摩擦因数X3，对倾斜板倾斜角没有影响，为了避免干扰，在EDEM仿真试验中，x3、X1、X2和X3的值均为0，采用已经标定的参数：三七种子与ABS塑料之间碰撞恢复系数x1=0.611，经过预仿真试验，取三七种子与ABS塑料之间静摩擦因数x2范围为0.35～0.65，取步长0.05，进行7组仿真试验，每组仿真试验重复5次取平均值，试验设计方案与结果如表3所示，表中y2为倾斜板倾斜角仿真值。

表3 静摩擦因数仿真试验设计方案与结果

为得到仿真试验中三七种子与ABS塑料的静摩擦因数与倾斜角的关系，对表3中的试验数据进行曲线拟合，得到2次多项式拟合曲线如图8所示，曲线方程为

图8 仿真试验静摩擦因数与倾斜角拟合曲线

(2)

式(2)的决定系数R2=0.995 6，接近1，表明该方程拟合的可靠度高。将台架试验的倾斜角实测值26.60°代入式(2)，求得x2=0.473，输入EDEM中进行验证，试验重复5次，计算平均值，仿真试验测得倾斜角仿真值为26.64°，与实测值的相对误差为0.14%，由此表明标定后的仿真结果与台架试验一致。确定EDEM仿真试验中，三七种子与ABS塑料之间静摩擦因数x2=0.473。

2.1.3滚动摩擦因数

滚动摩擦是指一个物体在另一个物体表面作无滑移的滚动或有滚动趋势时，两物体在接触部分产生的形变对滚动的阻碍作用。使用斜面滚动试验标定三七种子与ABS塑料的滚动摩擦因数，试验如图9所示。将三七种子放置于倾斜角γ2=35°的倾斜面上，在固定高度H2=30 mm，以初速度为0沿着倾斜面向下滚动，种子最终滚落至水平面上，待种子在水平面上完全静止，测量种子的水平滚动距离，台架试验重复5次，计算平均值，求得水平滚动距离实测值s=61.06 mm。

图9 滚动摩擦因数标定试验

三七种子间的碰撞恢复系数X1、静摩擦因数X2和滚动摩擦因数X3，对水平滚动距离没有影响，为了避免干扰，在EDEM仿真试验中，X1、X2和X3的值均为0，采用已经标定的参数：三七种子与ABS塑料之间碰撞恢复系数x1=0.611和静摩擦因数x2=0.473，经过预仿真试验，取三七种子与ABS塑料之间滚动摩擦因数x3范围为0.04～0.1，取步长0.01，进行7组仿真试验，每组试验重复5次取平均值，试验设计方案与结果如表4所示，表中y3为水平滚动距离仿真值。

表4 滚动摩擦因数仿真试验设计方案与结果

为得到仿真试验中三七种子与ABS塑料之间滚动摩擦因数与水平滚动距离的关系，对表4中的试验数据进行曲线拟合，得到2次多项式拟合曲线如图10所示，曲线方程为

(3)

图10 仿真试验滚动摩擦因数与水平滚动距离拟合曲线

式(3)的决定系数R2=0.985 9，接近1，表明该方程拟合的可靠度高。将台架试验的水平滚动距离实测值61.06 mm代入式(3)，求得x3=0.067，输入EDEM中进行验证，试验重复5次，计算平均值，仿真试验测得水平滚动距离仿真值为62.99 mm，与实测值的相对误差为3.16%，由此表明标定后的仿真结果与台架试验基本一致。EDEM仿真试验中，三七种子与ABS塑料之间滚动摩擦因数x3=0.067。

2.2 三七种子之间接触参数标定

三七等小粒径种子之间接触参数不易测量标定，采用传统的种子板进行测量标定时，测量结果的相对误差较大。种子在自由下落和堆积成型过程中，相互之间有碰撞力和摩擦力等，种子之间碰撞恢复系数、静摩擦因数和滚动摩擦因数均影响种子堆积的形状，所以本文通过堆积试验，对比台架试验和仿真试验标定三七种子之间参数，基于二次回归正交旋转组合试验的响应面优化方法，确定EDEM仿真试验中最佳的三七种子之间接触参数。

2.2.1堆积试验

堆积台架试验如图11a所示，在EDEM中创建仿真试验，如图11b所示。漏斗、挡板和底板材料均为ABS塑料，漏斗下端面距底端底板距离H3=75 mm，用挡板挡住漏斗下端面，将漏斗中装满三七种子，迅速抽掉漏斗下端挡板，三七种子自由下落，三七种子在底板上堆积，通过图像处理技术测量堆积角，图像处理如图12所示，对图像进行二值化处理，通过边缘检测提取轮廓曲线，对轮廓曲线的边缘点进行线性拟合，拟合直线与水平面的夹角即为堆积角。每组台架试验重复5次取平均值，求得堆积角实测值θ=30.9°。

图11 三七种子在ABS塑料板上堆积试验

在EDEM仿真试验中，采用已标定的三七种子与ABS塑料之间接触参数，选取三七种子之间碰撞恢复系数、静摩擦因数和滚动摩擦因数作为试验因素，仿真试验堆积角与台架试验堆积角的相对误差作为试验指标，堆积角相对误差δ计算式为

(4)

式中θ′——堆积角仿真值，(°)

θ——堆积角实测值，(°)

2.2.2最陡爬坡试验

采用最陡爬坡试验确定二次回归正交旋转组合试验因素的0水平及最优值区间，最陡爬坡试验方案与结果如表5所示，三七种子堆积图如图13所示，图13a～13g分别为最陡爬坡试验a～g组仿真试验的三七种子堆积图，图13h为台架试验的三七种子堆积图。

图12 图像处理

表5 最陡爬坡试验方案及结果

图13 三七种子堆积图

由表5可得，堆积角相对误差先减小后增大，d组仿真试验的相对误差最小，及最优值区间在d组试验附近，因此选取c、d和e组试验因素作为二次回归正交旋转组合试验的-1、0和1水平的因素。

2.2.3二次回归正交旋转组合试验

为寻求EDEM仿真试验中三七种子之间碰撞恢复系数、静摩擦因数和滚动摩擦因数的最佳参数组合，进行三因素二次回归正交旋转组合试验，仿真试验因素编码如表6所示，表7中试验因素编码A、B和C分别为三七种子间碰撞恢复系数、静摩擦因数和滚动摩擦因数编码值，仿真试验设计方案与结果如表7所示，试验结果为仿真堆积角与实测堆积角相对误差Y。

表6 仿真试验因素编码

表7 试验方案与结果

采用Design-Expert 8.0.6软件对试验数据(表7)进行多元回归拟合，得到堆积角相对误差Y的回归方程为

Y=1.73-1.28A+1.72B-0.51C+0.53AB+0.54AC-0.006 25BC+0.74A2+0.70B2+1.15C2

(5)

回归方程的显著性检验如表8所示。由表8可知，模型的拟合度极显著(P<0.01)。三七种子间静摩擦因数和滚动摩擦因数的交互项(BC)的P>0.1，对堆积角相对误差的影响不显著，其他各项的P检验均显著，说明相关试验因素对响应值的影响不是简单的线性关系，存在二次关系。失拟项P=0.253 9，不显著，说明无其他影响指标的主要因素。回归方程的决定系数R2=0.93，说明回归方程的预测值与实际值拟合良好。影响堆积角相对误差的因素由大到小为三七种子间静摩擦因数、碰撞恢复系数和滚动摩擦因数。

表8 回归方程方差分析

注： *表示影响显著(P<0.05)，** 表示影响极显著(P<0.01)。

2.2.4最优参数

利用Design-Expert 8.0.6软件的优化模块，以堆积角相对误差的最小值为目标，对回归方程进行求解，分析响应曲面，对回归模型进行寻优。目标及约束条件方程组为

(6)

得到多组最佳参数，最终选取三七种子间碰撞恢复系数为0.492、静摩擦因数为0.202和滚动摩擦因数为0.083。利用该最佳参数进行仿真堆积角试验，与实测试验结果在角度和堆形上均具有较高的相似性，堆积角相对误差仅为0.51%，说明所得最佳参数准确可靠，可用于后续的EDEM仿真试验。

3 验证试验

为了进一步验证三七种子离散元模型和仿真参数的可靠性，选取三七精密排种器进行试验验证，该排种轮两侧布置窝眼孔，可同时播种2行，窝眼孔充入1粒种子即为充种合格[29]，以排种器的充种合格率、漏播率和重播率作为试验指标，种层高度等试验条件均一致，在不同排种轮转速的试验条件下对比各试验指标的实测值和仿真值。

选取文山三七种子进行台架试验，选用直径为5.0～6.0 mm，排种器采用3D打印，在JPS-12型计算机视觉排种器性能检测试验台上搭建台架试验，验证试验的台架试验如图14a所示。将排种器的简化模型、三七种子离散元模型和标定的接触参数导入离散元软件EDEM中进行仿真试验，仿真试验如图14b所示。

图14 验证试验

根据GB/T 6973—2005 《单粒(精密)播种机试验方法》实施，每组试验统计2行，每行连续测量100个窝眼孔，每组试验重复4次取平均值。验证试验方案和结果如表9所示。

在不同排种轮转速的试验条件下，分别对比各试验指标的实测值和仿真值，由表9可得，排种器充种合格率、漏播率和重播率的相对误差均小于5%，排种器充种合格率、漏播率和重播率相对误差的平均值分别为0.64%、2.98%和2.74%，表明该三七种子离散元模型和接触参数可用于离散元仿真试验。

4 结论

(1)采用逆向工程技术，通过三维扫描仪器，得到三七种子轮廓模型；基于粘结颗粒模型，在EDEM软件中，通过自动填充方式，在三七种子轮廓模型中填充98颗半径为0.5 mm的球形颗粒，建立了三七种子离散元模型。

(2)结合台架试验和仿真试验，在EDEM软件中标定三七种子与ABS塑料之间接触参数，通过碰撞弹跳试验，得到三七种子颗粒模型与ABS塑料之间碰撞恢复系数为0.611；通过斜面滑移试验，得到三七种子颗粒模型与ABS塑料之间静摩擦因数为0.473；通过斜面滚动试验，得到三七种子颗粒模型与ABS塑料之间滚动摩擦因数为0.067。

表9 验证试验方案与结果

(3)通过堆积试验，利用Design-Expert 8.0.6软件处理试验数据，基于二次回归正交旋转组合试验的响应面优化方法，确定EDEM仿真试验中三七种子之间最佳的接触参数，得到三七种子之间碰撞恢复系数、静摩擦因数和滚动摩擦因数分别为0.492、0.202和0.083。

(4)采用三七精密排种器进行试验验证，在JPS-12型计算机视觉排种器性能检测试验台上进行台架试验，在离散元软件EDEM中进行仿真试验。在不同排种轮转速的试验条件下，分别对比排种器充种合格率、漏播率和重播率的实测值和仿真值，得到充种合格率、漏播率和重播率的相对误差均小于5%，排种器充种合格率、漏播率和重播率相对误差的平均值分别为0.64%、2.98%和2.74%，进一步验证了离散元仿真试验的可靠性。