收获期豌豆籽粒离散元仿真参数标定

张克平，侯传凯，孙步功，苏占科

(甘肃农业大学机电工程学院，甘肃 兰州 730070)

豌豆是一年生或越年生草本植物，是集粮、菜、肥、饲于一体的多用途作物，豌豆生长周期短、营养价值高，可作为豆粕和玉米等原料的替代物，从而有效降低饲料成本[1-2]。豌豆作物具有后熟特性，收获期豌豆籽粒含水率高且跨度大，机械化收获过程中易发生籽粒损伤现象，尤其在玉米/豌豆间作种植模式下，先熟的豌豆作物收获作业空间狭窄，通用谷物收获机械适用性较低，仍然以人工收获为主，劳动强度大且作业效率低。

豌豆籽粒和茎秆的物理和生物力学特性研究可为相关收获机械的研发提供理论依据。代治国等[3]对豌豆籽粒进行了压缩和剪切试验，探讨了豌豆籽粒的含水率与破碎负载、弹性模量、剪切力、剪切强度、硬度之间的函数关系；
王君等[4]通过对豌豆茎秆剪切和拉伸试验，发现豌豆茎秆的抗剪强度随含水率的增大呈先上升后下降的趋势，有利于确定豌豆收割的最佳时期；
Dange等[5]利用钟摆式动态试验仪对不同处理条件下的豌豆茎秆所需的力和能量进行了研究，得到切削速度、叶片斜角、含水量、茎径是影响豌豆茎秆切削能量和力的主要制约因素。近年来，离散元法在农作物收获机械设计和仿真分析中的应用日益广泛[6]，为了使仿真更贴近实际，分析过程中需要设置精确的农业物料参数，Balevicius等[7]通过豌豆和玻璃板的滑动试验测得两者之间的静摩擦系数；
Azadbakht等[8]通过冲击和摩擦试验探究了水分和能量对豌豆脱粒的影响，结果显示豆荚脱粒量随着冲击能和摩擦能的增大而增加，籽粒摩擦系数随其含水量的增加而增大，而对于收获期豌豆离散元仿真参数的标定鲜有报道。

本研究拟采用物理试验和仿真试验相结合的方法对豌豆籽粒离散元仿真参数进行标定，首先通过台架试验确定参数区间，依次通过 Plackett-Burman 试验、最陡爬坡试验以及 Box-Behnken 试验筛选对豌豆籽粒堆积角影响显著的参数并对其进行最优化处理，非显著参数则由离散元仿真来标定。研究结果可为豌豆机械化收获过程的离散元仿真分析提供参考，进而为豌豆收获机械的设计和优化提供依据。

1.1 豌豆籽粒物性参数

试验所用收获期豌豆于2021年7月30日采自甘肃省定西市渭源县试验田，利用电热鼓风干燥箱和分析天平(精度0.1 mg)测定其含水率为20.94%～25.69%。豌豆籽粒的物性参数主要包括百粒质量、三轴尺寸、密度、泊松比和弹性模量等，随机选取5组豌豆籽粒，每组100粒，利用电子天平(精度0.01 g)测其百粒质量并取均值；
利用数显游标卡尺(精度0.01 mm)测量三轴尺寸并取均值，利用量筒(精度1 ml)测定密度；
利用质构仪(TA.XT plus，距离精度0.001 mm，力量精度0.1 g)对豌豆籽粒开展压缩试验获取其泊松比和弹性模量[9-10]。豌豆籽粒的基本物性参数列于表1。

1.2 豌豆籽粒自然堆积角

豌豆籽粒的堆积角通过圆筒提升试验测量。将钢制圆筒(内径65 mm，高180 mm)垂直放置在水平钢板上，上部与自动伸缩杆相连，圆筒内装填满籽粒后以12 mm·s-1的速度向上提升[11]，豌豆籽粒自由下落在水平钢板上形成一个锥形颗粒堆，待其静止后用相机拍摄颗粒堆的正视照片(图1)，通过Origin 2018软件处理得到颗粒堆积角度，试验重复10次并取其平均值，得到豌豆籽粒堆积角为19.841°。

表1 豌豆籽粒基本物性参数Table 1 Basic physical parameters of pea grains

图1 豌豆籽粒堆积试验Fig.1 Pea grain stacking experiment

2.1 碰撞恢复系数

2.1.1 豌豆籽粒-钢板碰撞恢复系数 碰撞恢复系数按照其表达的物理性质可分为3种，即运动学恢复系数、动力学碰撞恢复系数和能量恢复系数。运动学恢复系数定义为两物体在碰撞后质心的法向分速度与碰撞前的质心法向分速度的比值，此方法主要针对球状近刚性物体的低速碰撞问题，测定出只与材料相关的恢复系数的数值[12]。

试验装置及其原理如图2所示，颗粒自距离地面高2h处释放，在重力作用下与距离地面h并倾斜45°的钢板发生碰撞，最终籽粒落在初始释放位置水平距离L的地面[13]。由此可得，碰撞恢复系数可以简化为：

(1)

根据上述原理测量豌豆籽粒与钢板的碰撞恢复系数，释放高度取300 mm，自由下落过程重复10次，测量弹出水平距离获得最大值和最小值，根据式(1)计算豌豆籽粒-钢板的碰撞恢复系数取值范围为0.460～0.613。

2.1.2 豌豆籽粒间碰撞恢复系数 采用碰撞弹出法测量豌豆籽粒间的碰撞恢复系数[14]。测量之前制作豌豆籽粒碰撞板，取一块70 mm×55 mm白纸，将预先挑选好的豌豆籽粒用双面胶粘在纸上，尽可能减小豌豆籽粒间的空隙，制作的豌豆籽粒碰撞板如图3所示。试验时将制作好的碰撞板粘在图2所示测试装置的倾斜45°钢板上，重复2.1.1节中相同试验步骤测得豌豆籽粒间碰撞恢复系数范围为0.274～0.488。

2.2 静摩擦系数

2.2.1 豌豆籽粒-钢板静摩擦系数 静摩擦系数是物体所受到的最大静摩擦力与法向压力的比值，常用测量方法为斜面滑动法，其测试原理如图4所示，豌豆籽粒在斜面上静止时，受到重力mg(N)、静摩擦力f(N)、斜面对豌豆的支撑力N(N)，将籽粒的重力分解为垂直于斜面的力T(N)和平行于斜面的力F(N)，当斜面角度θ(°)小于滑动的临界角度时，平行于斜面的力F小于物体受到的静摩擦力f，物体在斜面上保持静止不动，随着斜面角度θ增大，增加到大于临界角度，斜面上的力F也增加到大于颗粒与斜面之间的静摩擦力，颗粒开始沿斜面向下滑动。故静摩擦系数fs与斜面角度θ的关系可按式(2)表示：

图2 碰撞恢复系数原理图(a)及其测试装置(b)Fig.2 Collision recovery coefficient schematic diagram (a) and its test device (b)

图3 豌豆籽粒测试板Fig.3 Pea grain test board

(2)

使用斜面滑动法测定豌豆籽粒与钢板的静摩擦系数，豌豆籽粒近似圆形，稍微施加外力就会做滚动运动，在试验过程中为了防止豌豆籽粒发生滚动，降低试验误差使参数与实际值更接近，将几个豌豆籽粒粘结成矩形，如图5所示，试验开始时，将钢板水平放置，把粘结成矩形的豌豆籽粒放置在钢板不转动的一侧，缓慢抬动钢板，使钢板绕转动轴侧做匀速圆周运动，当粘结籽粒有滑动的趋势时停止转动，并记录钢板与水平面的角度，试验重复10次，测得豌豆籽粒-钢板静摩擦系数范围为0.412～0.573。

2.2.2 豌豆籽粒间静摩擦系数 采用斜面滑动法测量豌豆籽粒间静摩擦系数[15]。如图6所示，将图3中制作好的豌豆籽粒测试板粘在匀速运动的钢板上，将一粒豌豆籽粒放在测试板上，使钢板匀速运动直到放上的豌豆籽粒开始滑动时停止运动，并记录钢板转过的角度，试验重复10次，测得豌豆籽粒间静摩擦系数范围为0.287～0.657。

图4 静摩擦系数测试原理图Fig.4 Schematic diagram of static friction coefficient test

图5 豌豆籽粒矩形板Fig.5 Pea grain rectangular plate

2.3 滚动摩擦系数

2.3.1 豌豆籽粒-钢板滚动摩擦系数 滚动摩擦是指当一个物体在另一物体表面作无滑动的滚动或有滚动的趋势时，由于物体在接触部分受压发生形变而产生阻碍滚动的作用，豌豆籽粒与钢板的滚动摩擦系数通过斜面滚动法和能量守恒定律相结合测定。将豌豆籽粒从一固定倾斜角度θ为16°钢板以初速度0沿着钢板滚动，由于摩擦力存在，豌豆籽粒停在水平钢板的某一位置，豌豆籽粒在倾斜钢板上滚过距离为S=200 mm，在水平钢板上滚过距离为L。滚动摩擦系数μ可以由能量守恒定律求得(式3)，试验重复10次，测得豌豆籽粒-钢板滚动摩擦系数范围为0.068～0.139。

mgSsinθ=μmg(Scosθ+L)

(3)

2.3.2 豌豆籽粒间滚动摩擦系数 采用初始滚动时刻内近似能量守恒获得的滚动摩擦系数来预估真实值[11]。测试板由平放缓慢转动，物料受重力而滚动，如果在倾斜面上开始滚动的一微小角度α(°)，其速度由零增加，速度变化极小，此时假设动能为0。所以测量物料的重力势能等于其在滚动中损失的能量，由以上分析并简化可得滚动摩擦系数计算公式(式4)，试验重复10次取平均值，测得豌豆籽粒间滚动摩擦系数范围为0.337～0.627。

μ=tanα

(4)

3.1 仿真模型

3.1.1 豌豆籽粒离散元模型 豌豆籽粒的基本物性参数已经通过试验测定，由于其外形近似椭球形，为了便于在离散元仿真软件中建立豌豆籽粒的颗粒模型，前文已对其长度L、宽度K、厚度H尺寸测定，取其长轴长度a=L，短轴长度b=(K+H)/2[16]，由此得到豌豆籽粒模型长轴和短轴长度分别为8.220 mm和6.825 mm。在EDEM软件中通过多球颗粒组合功能建立豌豆籽粒的离散元模型，如图7所示。

图6 籽粒间静摩擦系数测试原理图Fig.6 Test schematic diagram of coefficient of static friction between grains

3.1.2 仿真参数 在建立仿真模型时，颗粒之间、颗粒-钢板之间的相互作用影响较大，在仿真参数标定试验时，为了提高软件仿真效率，采用固定尺寸的豌豆籽粒离散元模型[17]。本研究中通过台架试验和相关文献确定各仿真参数如表2所示。

3.1.3 圆筒提升模型 在EDEM中等尺寸建立无底圆筒以及矩形底板，圆筒上侧建一颗粒工厂，颗粒生成方式为dynamic，颗粒的生成速率为3 000个·s-1，生成数量不限，时间步长1.55×10-5s，网格尺寸为最小颗粒半径的3倍，待颗粒填满圆筒后停止继续生成颗粒，静置3 s待生成颗粒沉积，尽可能减小空隙使其达到平衡。以与台架试验相同的速度提升圆筒，颗粒下落在底板上形成颗粒堆，如图8所示。

3.2 影响显著性参数筛选试验

应用Design-Expert 10软件进行影响显著性参数筛选试验，以豌豆籽粒的堆积角为响应值，通过Plackett-Burman试验筛选出对豌豆籽粒堆积角影响显著的参数，如表3所示，试验变量的最大值、中间值、最小值分别编码+1、0、-1。仿真试验方案与结果列于表4，待仿真试验结束后，采用台架物理试验测堆积角的方法，测量豌豆籽粒的仿真堆积角。

图7 豌豆籽粒离散元模型Fig.7 Discrete element model of pea grain

表2 离散元仿真参数Table 2 Discrete element simulation parameters

利用Design-Expert 10软件对仿真堆积角结果进行分析，如表5所示，豌豆籽粒与豌豆籽粒静摩擦系数(E)和豌豆籽粒与豌豆籽粒滚动摩擦系数(F)对豌豆籽粒堆积角影响极显著，豌豆籽粒与豌豆籽粒碰撞恢复系数(D)则影响显著。而其余参数影响相对较小，对于影响较小的仿真参数标定则采用前述台架试验与之后模拟仿真台架试验相结合来完成。

3.3 非显著性参数仿真标定

3.3.1 豌豆籽粒-钢板碰撞恢复系数仿真试验标定 仿真试验如图9所示，取释放高度为300 mm，因豌豆籽粒和钢板的静摩擦系数、滚动摩擦系数对籽粒碰撞弹出距离无影响，为避免干扰，除了豌豆籽粒与钢板碰撞恢复系数外，其他参数值均为0。将前述通过台架试验测定的豌豆籽粒和钢板碰撞恢复系数取值范围0.460～0.613，分为7组仿真试验，每组试验重复5次并取平均值，试验设计方案及结果列于表6，其中Y1为每组进行多次仿真试验得到籽粒与板碰撞后弹出距离的均值，每组弹出距离由EDEM软件后处理工具Ruler测量得到。

以豌豆籽粒与钢板碰撞恢复系数A为试验因素，弹出距离Y1为评价指标，对表6中的数据进行二次多项式曲线拟合见式(5)，二者关系可靠度较高。将Y1=170.7 mm代入式(5)得到A=0.505，在仿真试验中将豌豆籽粒与钢板碰撞恢复系数设定为该值，在释放高度为300 mm的条件下进行5次试验并取平均值，得到弹出距离为172.4 mm，与真实台架试验得到的弹出距离相对误差为0.996%。表明标定后的仿真试验结果和真实台架试验基本相同，所以碰撞恢复系数A=0.505。

图8 豌豆籽粒堆积角仿真模型Fig.8 Simulation model of pea grain stacking angle

图9 碰撞恢复系数仿真试验Fig.9 Simulation test of collision recovery coefficient

表3 Plackett-Burman试验参数Table 3 Parameters of the Plackett-Burman test

表4 Plackett-Burman 试验方案及结果Table 4 Plackett Burman test scheme and results

表5 Plackett-Burman 试验参数显著性分析Table 5 Significance analysis of Plackett-Burman test parameters

Y1=-231.37A2+517.17A-31.301
R2=0.9864

(5)

3.3.2 豌豆籽粒-钢板静摩擦系数仿真试验标定 在EDEM仿真中，采用多球组合的方式生成4个豌豆籽粒且其排列成矩形，如图10所示。仿真参数设置如下：豌豆籽粒与钢板碰撞恢复系数为已标定为0.505，豌豆籽粒与钢板静摩擦系数将台架试验确定的范围分为6组，试验时每组测5次并取平均值，试验方案与结果列于表7，为了避免干扰则其余参数均为0。

表6 碰撞恢复系数仿真试验方案与结果Table 6 Collision recovery coefficient simulation test scheme and results

图10 静摩擦系数仿真试验Fig.10 Static friction coefficient simulation test

表7 静摩擦系数仿真试验方案与结果Table 7 Static friction coefficient simulation test scheme and results

以豌豆籽粒与钢板静摩擦系数B为试验因素、倾斜角度Y2为评价指标，倾斜角度由EDEM软件后处理工具Protractor测量得到。对表7中的数据进行二次多项式曲线拟合见式(6)，二者关系可靠度较高。将Y2=25.5°代入式(6)得到B=0.462，在仿真试验中将豌豆籽粒与钢板静摩擦系数设定为该值，进行5次试验并取平均值，仿真试验倾斜角度为25.814°，与台架试验相对误差1.231%，表明标定后的仿真试验结果和真实台架试验值接近，所以静摩擦系数B=0.462。

Y2=3.8554B2+49.28B+1.9226R2=0.9972

(6)

3.3.3 豌豆籽粒-钢板滚动摩擦系数仿真试验标定 在EDEM软件中构建一倾斜角度为16°的长方形钢板(200 mm×300 mm)，与之相接一个400 mm×800 mm的水平钢板，如图11所示。仿真参数设置如下：碰撞恢复系数A=0.505，静摩擦系数B=0.462，将台架试验得到的滚动摩擦系数区间分为6份，试验方案与结果列于表8，为了避免干扰，其余参数均为0。

以豌豆籽粒与钢板滚动摩擦系数C为试验因素、滚动距离Y3为评价指标，每组试验的水平滚动距离由EDEM软件后处理工具Ruler测量得到。对表8中的数据进行二次多项式曲线拟合见式(7)，二者关系可靠度较高。将Y3=394.21 mm代入式(7)得到C=0.09，在仿真试验中将豌豆籽粒与钢板滚动摩擦系数设定为该值，进行5次试验并取平均值，仿真试验水平滚动距离388.617 mm，与台架试验相对误差1.439%，表明标定后的仿真试验结果和真实台架试验接近，所以滚动摩擦系数C=0.09。

图11 滚动摩擦系数仿真试验Fig.11 Simulation test of rolling friction coefficient

表8 滚动摩擦系数仿真试验方案与结果Table 8 Rolling friction coefficient simulation test scheme and results

Y3=51621C2-16112C+1427
R2=0.9984

(7)

3.4 最陡爬坡试验

由表5可知豌豆籽粒间静摩擦系数E和碰撞恢复系数D对堆积角效应为正值，滚动摩擦系数F的效应为负值，对这3个显著性参数进行最陡爬坡试验，以堆积角相对误差为评价标准，确定仿真参数最优范围，试验方案与结果列于表9，可知4号的相对误差最小，故将4号作为后续 Box-Behnken响应面试验的中心点。

3.5 Box-Behnken 试验及回归模型

根据最陡爬坡试验结果，开展Box-Behnken堆积角仿真试验，试验设计方案及结果列于表10。利用Design-Expert 10软件对表10中的数据进行分析，得到豌豆籽粒堆积角θ与3个显著性参数的二阶回归方程为：

θ=42.81395-75.81167D+42.64176E-

85.93713F+114.62916DE-235.76584DF-

79.45014EF+166.07626D2-

41.62253E2+226.31540F2

(8)

表9 最陡爬坡试验方案与结果Table 9 Test scheme and results of steepest climb

表10 Box-Behnken 试验设计方案及结果Table 10 Box-Behnken test design scheme and results

表11 Box-Behnken 试验回归模型方差分析Table 11 Variance analysis of Box-Behnken test regression model

图12 DE和DF交互效应图Fig.12 Interaction effect diagram of DE and DF

由图12中DE交互效应图可知，随着籽粒间恢复系数D和静摩擦系数E的增加，籽粒堆积角增加，当两参数中任意一个取较大值时，籽粒堆积角随另一参数增加而增大的速率越明显。从DF交互效应图趋势可知，当籽粒间恢复系数D取定值时，籽粒堆积角随滚动摩擦系数F的增大而减小，当滚动摩擦系数F取定值时，籽粒堆积角随恢复系数D变化波动较小。

3.6 仿真参数最优组合与试验验证

在Design-Expert 10软件中，为了得到使仿真试验与物理试验最接近的最优参数组合，以物理试验堆积角19.841°为目标，对回归方程(8)优化求解，得到最佳参数组合；
豌豆籽粒间碰撞恢复系数、静摩擦系数、滚动摩擦系数分别是0.364、0.519、0.444。以此参数作为仿真参数进行多组仿真试验，得到堆积角仿真均值19.714°与实际堆积角相对误差为0.64%。

1)豌豆籽粒间、豌豆籽粒-钢间碰撞恢复系数范围分别为0.274～0.488、0.460～0.613；
静摩擦系数范围分别为0.287～0.657、0.412～0.573；
滚动摩擦系数范围分别为0.337～0.627、0.068～0.139。

2)对豌豆籽粒堆积角影响显著的参数有豌豆籽粒-豌豆籽粒碰撞恢复系数、豌豆籽粒-豌豆籽粒静摩擦系数、豌豆籽粒-豌豆籽粒滚动摩擦系数；
而豌豆籽粒-钢板碰撞恢复系数、豌豆籽粒-钢板静摩擦系数、豌豆籽粒-钢板滚动摩擦系数对堆积角影响不显著。

3)以实际堆积角为目标，通过最陡爬坡试验和Box-Behnken试验得到显著性参数的最优组合：豌豆籽粒-豌豆籽粒碰撞恢复系数为0.364，豌豆籽粒-豌豆籽粒静摩擦系数为0.519，豌豆籽粒-豌豆籽粒滚动摩擦系数为0.444；
非显著性参数通过仿真试验标定：豌豆籽粒-钢碰撞恢复系数0.505、静摩擦系数0.462、滚动摩擦系数0.090，与实测值的相对误差分别为0.966%、1.233%和1.439%。以标定的参数进行堆积角仿真试验，得到最优组合下堆积角均值为19.714°，与台架试验相对误差为0.64%。