李飞翔， 王鹏， 王云飞， 葛越锋， 唐凯怿， 李得志

（洛阳智能农业装备研究院有限公司，河南 洛阳 471000）

气力式排种器因具有排种精度高、对种子损伤少等优点而被广泛应用。在实际排种作业中，种子处于气固两相环境中，运动情况较为复杂，离散元（discrete element method，DEM）与计算流体动力学（computational fluid dynamics，CFD）耦合法能够有效分析排种器中流场的变化以及种子运动的受力情况，为排种器的结构设计及参数优化提供依据[1-3]。为提高玉米的产量及质量，常常对玉米种子进行包衣处理。种子包衣是将杀虫剂、微生物肥料、着色剂等非种子材料包裹在种子表面，以提高种子抗逆性、抗病性，对农业机械化、现代化及农业可持续发展有着重要作用。研究人员基于离散元法对包衣种子已做了大量研究：闫聪杰[4]设计一种大豆包衣装置，并基于离散元法对其包衣过程进行分析；曲芳[5]针对现有大豆种子包衣搅拌装置包衣不均、破碎率高的问题，基于离散元法设计一种搅拌叶片并优化设计；胡建平等[6]以番茄磁粉包衣种子为对象，基于离散元法研究了种箱振动频率、振幅对种群运动规律及种箱供种性能的影响。

DEM-CFD 常采用 EDEM 与 Fluent 软件的耦合，要求EDEM 中颗粒体积小于Fluent 中网格体积。因此，EDEM 中常采用黏结颗粒模型（bonded particle model，BPM），用较小的颗粒黏结成物料颗粒模型[7]。BPM 需设置单位面积法向与切向刚度、极限法向与切向应力、黏接半径等参数[8]，现阶段对于玉米种子的仿真参数设置研究较少。基于离散元的Hertz-Mindlin 模型以堆积角为响应，王云霞等[9]以玉米颗粒种间静摩擦与滚动摩擦系数为因素，设计试验对其进行标定；王美美等[10]对玉米籽粒离散元仿真参数进行标定，为玉米籽粒离散元仿真提供参考；李秀清[11]基于离散元法对玉米粉碎过程进行分析，采用测定的楔形玉米籽粒建模，借鉴玉米籽粒压缩试验数据，将极限法向应力、极限切向应力作为固定值，选取法向刚度、切向刚度、黏接半径3个因素进行标定，得到玉米籽粒的仿真接触参数。

为得到玉米包衣种子BPM 离散元仿真所需的精确的种间接触参数，本研究将利用激光扫描仪得到的玉米包衣种子轮廓作为颗粒模板，采用BPM 以小颗粒黏结方式生成玉米包衣种子元颗粒模型，以堆积角作为响应值，基于Isight 软件的近似模型与DOE 联合模块，通过因素显著性筛选、最优逼近以及响应面优化对玉米包衣种子的仿真接触参数进行标定，以期为玉米包衣种子气力式排种器的两相耦合仿真提供参考。

1 材料与方法

1.1 试验材料获取

供试玉米品种为郑单958，纯度＞96%，净度＞90%，发芽率＞90%，含水率＞13.5%，由河南金博士种业股份有限公司提供。由于种子形状多样，将郑单958 玉米种子分为大扁形、长扁形、类球形3类（图1），比例为2∶3∶5。

图1 玉米包衣种子形状Fig.1 Shape of corn coated seed

1.2 接触模型

为方便与Fluent 软件耦合，EDEM 颗粒间接触选用BPM。该模型接触情况下，黏结键可以抵抗切向和法向运动，直至达到最大的法向和切向极限应力为止。黏结作用产生后，元颗粒所受的力和力矩设置为0，并根据每个时间步长进行增量调整[12-13]。

式中，δFn为法向力增量（N）；δFτ为切向力增量（N）；δMn为法向力矩增量（N·m）；δMτ为切向力矩增量（N·m）；Y为基础区域面积（m2）；Rb为黏结半径（m）；J为截面极惯性矩；Sn为法向刚度（N·m-3）；Sτ为切向刚度（N·m-3）；vn为法向速度（m·s-1）；vτ为切向速度（m·s-1）；wn为法向角速度（r·s-1）；wτ为切向角速度（r·s-1）；δt为时间步长（s）。

当法向和切向剪应力超过预定值σmax和τmax时，黏结就会断裂。

黏结应力与刚度值的不同设置代表了不同的材料属性，高的刚度值将产生高的结合力和应力。小颗粒在黏结之前，根据Hertz-Mindlin 模型相互作用，当达到键形成时间时，被定义的小颗粒将会黏结在一起形成元颗粒。

1.3 堆积角测定方法

依据已有对堆积角的研究[14-18]，以漏斗法对玉米包衣种子的堆积角进行测定，如图2所示，测定5次取其平均值，测得郑单958玉米包衣种子的堆积角为29.56o±0.27o。由于标定参数为颗粒间接触，为消除颗粒与部件的影响，采用亚克力材质加工落料底盘，底盘直径（L）为150 mm，底盘凹面深度为30 mm。待底盘落满种子后，继续落种；待玉米包衣种子溢出底盘后，观察落料中颗粒堆积高度的变化；待堆积高度无显著性变化时，用钢尺测定包衣种子的堆积高度（h），由式（9）计算玉米种子堆积角。

图2 堆积试验Fig.2 Stacking test

1.4 仿真试验

玉米种子为非规则形状颗粒，为精确建立种子轮廓模型以提高标定精准性，选择不同种类形状的种子，采用海克斯康Tango-S手持式激光扫描仪对种子的外轮廓进行扫描。得到种子的点云数据后，通过CATIA 软件对种子点云数据进行处理，最终得到玉米包衣种子的轮廓模型。

仿真模型采用与实际试验模型相同尺寸，按照筛选的大扁形、长扁形、类球形种子2∶3∶5的比例生成颗粒。颗粒堆积仿真试验如图3 所示，为节省仿真时间，先以圆柱底盘为虚拟工厂生成颗粒，待底盘表面落满颗粒后，导出仿真时间为0时的文件；再以漏斗为虚拟工厂，进行堆积角仿真试验，待底盘溢出一定数量颗粒后，停止生成颗粒。采用后处理自带的量角器，测量堆积角。

图3 玉米种子颗粒堆积的模拟仿真Fig.3 Simulation test model for repose angle

参照玉米种子BPM 参数的设置[19-22]，本研究各参数取值范围如表1 所示。仿真本征参数设定为玉米种子泊松比0.4，密度1 197 kg·m-3，剪切模量1.37×108Pa[23-26]。

表1 离散元仿真参数Table 1 Parameter in DEM simulation

1.5 标定方法

以玉米包衣种子实际堆积角为响应值，基于Isight 软件的近似模型和DOE 联合模块，通过因素显著性筛选、最优逼近以及响应面优化对玉米包衣种子的仿真接触参数进行标定。本研究的Plackett-Burman 设计以玉米包衣种子堆积角为响应值对BPM 参数的显著性进行筛选。低水平设定为最初原始水平，高水平设为低水平的2倍，为方便对试验误差进行分析，同时设定虚拟参数，试验参数如表2所示。利用Isight软件的近似模型与DOE模块联合对Plackett-Burman结果进行方差分析，运行前，添加ANOVA Table、主效应图、Pareto图对各参数的显著性进行分析。Plackett-Burman试验后，以对玉米包衣种子堆积角影响显著的参数进行最陡爬坡，以便能快速进入到最优值的附近区域。最陡爬坡试验根据Plackett-Burman试验所得的回归系数来确定爬坡步长。

2 结果与分析

2.1 轮廓模型扫描分析

图4 为玉米包衣种子扫描的轮廓模型，可以看出，通过点云处理后得到的不同种类种子轮廓与实际种子轮廓一致，能够用于仿真分析。

2.2 种子模型填充分析

通过激光扫描得到种子轮廓后，导入EDEM软件中作为颗粒模板。结合玉米气力式排种器研究进展[27-28]，此处以0.5 mm 的小球颗粒对玉米种子颗粒模板进行填充，如图5 所示。对颗粒模板进行上半部透明化显示，填充的颗粒与模板贴合程度良好，通过小颗粒间形成黏结键建立玉米种子元颗粒黏结模型，建立的颗粒模型与玉米包衣种子外形相似。

图5 玉米包衣种子仿真模型Fig.5 Simulation model of corn coated seeds

2.3 参数标定结果分析

2.3.1 Plackett-Burman 试验结果分析 Plackett-Burman 设计及结果如表3 所示，堆积角范围在23.03°～30.31°之间，各仿真值没有特别跳跃性情况。

表3 Plackett-Burman 试验设计及结果Table 3 Design and results of Plackett-Burman test

2.3.2 Plackett-Burman 试验显著性分析 主效应、Pareto结果如图6所示，可清晰观察各个因素对响应的影响趋势；ANOVA结果如表4所示，可知对玉米包衣种子堆积角影响极其显著的参数为：玉米种子-玉米种子静摩擦系数、法向刚度、切向刚度。其余参数结合已有研究进展[29-30]取值（玉米种子-玉米种子恢复系数0.18、玉米种子-玉米种子滚动摩擦系数 0.05、法向应力 1.1×107N·m-3、切向应力4.1×106N·m-3、黏结半径0.7 mm）进行最陡爬坡以及响应面试验设计。

表4 Plackett-Burman 试验参数显著性分析Table 4 Analysis of significance of parameters in Plackett-Burman test

图6 玉米种子对颗粒堆积响应分析Fig.6 Response analysis of corn seed to grain accumulation

2.3.3 最陡爬坡试验结果分析 最陡爬坡试验设计及结果如表5所示，2号水平仿真堆积角与试验值的误差最小。由此选取2 号水平为中心点，1号、3号水平为低、高水平进行后续优化设计。

表5 最陡爬坡试验及结果Table 5 Design and results of steepest ascent test

2.3.4 Box-Behnken 试验结果分析 根据最陡爬坡试验结果，选取显著性参数进行3个水平试验设计。选取3个中心点对误差进行评估，结果如表 6 所示，基于 Isight 软件 RSM 优化模块，建立3个参数与堆积角间的二阶回归方程。

表6 Box-Behnken 试验及结果Table 6 Design and results of Box-Behnken test

Box-Behnken 试验方差分析结果如表7 所示，玉米包衣种间静摩擦系数(H)、法向刚度(J)、切向刚度（K）、法向刚度-切向刚度（JK）、静摩擦系数-法向刚度（HJ）、静摩擦系数二次项（H2）以及切向刚度二次项（K2）的P值＜0.05，表明这些参数对堆积角的影响显著。失拟项P=0.093＞0.05，变异系数CV=0.801%，表明模型良好，试验有较高的可靠性。决定系数R2=0.997 4，校正决定系数RAdj2=0.992 8，预测决定系数RPre2=0.958 9，表明模型能真实地反映实际情况。试验精密度Adep precision=46.973，表明模型良好的精确度。

表7 Box-Behnken 试验方差分析Table 7 ANOVA of quadraticl of Box-Behnken test

2.3.5 回归模型交互效应分析 根据优化回归模型方差分析结果可知，种间静摩擦系数-法向刚度（HJ）以及法向刚度-切向刚度（JK）交互项对玉米包衣种子堆积角影响显著（P＜0.05）。在玉米种间静摩擦系数（H）为0.2 以及切向刚度（K）为0.5 情况下，应用Isight 软件对静摩擦系数-法向刚度（HJ）以及法向刚度-切向刚度（JK）交互效应的响应曲面进行绘制，如图7所示，可以直观反映交互项对堆积角的影响。由HJ曲面可知，相对玉米包衣种间静摩擦系数（H），法向刚度（J）的效应面曲线比较陡，表明其对堆积角影响较为显著。由JK曲面可知，相对于法向刚度（J），切向刚度（K）的效应面曲线比较陡，表明其对堆积角影响较为显著。

图7 HJ与JK的交互效应Fig.7 Interaction effect of HJ and JK

2.3.6 试验优化结果分析 应用Isight软件RSM优化模块，以玉米包衣种子实际测定的堆积角29.56o为目标，对回归方程寻优求解得到玉米包衣种间静摩擦系数为0.269，法向刚度为2.54×108N·m-3，切向刚度为5.93×107N·m-3。用所得最佳参数组合进行堆积仿真试验，结果如图8 所示。仿真试验所得堆积角为29.85o，与实际试验值29.56o误差为0.98%，表明最佳参数组合可行。

图8 仿真与实际试验对比Fig.8 Comparison between simulation and actual test

3 讨论

现有堆积试验多采取物料直接堆积于材料表面，实际仿真标定中，物料与部件间的参数设置对堆积结果影响较大，为减小试验误差，堆积试验采用凹槽底盘，使得物料落满底盘后，在物料表面形成堆积。本研究选用郑单958 玉米包衣种子，在对种子进行分类筛选后，通过激光扫描仪对轮廓较好的种子进行扫描得到点云数据，采用逆向工程技术，通过CATIA 软件对点云数据处理，最终得到玉米种子仿真模型。

根据主效应和Pareto 图分析，筛选出对堆积角影响显著的参数：玉米种子-玉米种子静摩擦系数、法向刚度、切向刚度。基于RSM 模块，建立3个显著性参数与堆积角间的二次回归模型，由方差分析结果可知，交互项静摩擦系数-法向刚度（HJ）、法向刚度-切向刚度（JK）以及3个显著性参数的二次项（H2）对玉米包衣种子堆积角影响极显著。

以玉米包衣种子实际测定堆积角为目标，对回归方程进行寻优求解，得到参数的最佳组合为：种间静摩擦系数0.269，法向刚度2.54×108N·m-3，切向刚度5.93×107N·m-3。用所得最佳参数组合进行仿真试验，与实际值误差为0.98%，仿真与实际试验值差异不显著，表明标定所得参数可用于玉米包衣种子离散元BPM 仿真，为玉米气力式排种器的耦合模拟提供参考。