花生荚果离散元仿真参数标定*

徐效伟，魏海，颜建春，鲍国丞，杜元杰，谢焕雄

(1. 农业农村部南京农业机械化研究所，南京市，210014; 2. 江苏省农业机械试验鉴定站，南京市，210017)

0 引言

花生(peanut)，原名落花生，又名“长生果”“泥豆”“番豆”。花生是一种重要的油料作物，同时也是优质蛋白质来源。2019年全国花生种植面积约为4 633 khm2，产量为1 752 kt，花生种植面积和总产量常年居于世界前列[1-8]。近年来花生生产机械化率逐年提高，但相关机械较复杂[9-12]，通过传统试验进行前期的试验验证，费时费工，基于离散元法的EDEM软件在工农业中应用越来越广泛，因此运用EDEM软件进行前期仿真试验，可以有效缩短研发周期、降低研发成本。

目前，离散元仿真分析参数标定的研究对象多为肥料、土壤和种子等。韩树杰等[13]采用物理试验与仿真试验相结合的方法，以堆积角为响应值对散体厩肥离散元参数进行标定。邢洁洁等[14]对海南地区砖红土壤离散元接触参数进行仿真试验标定，以实际堆积角为目标值，进行仿真试验，获得最佳参数组合。于庆旭等[15]采用逆向工程技术对三七种子进行三维建模，并利用实际试验与仿真试验相结合的方法对三七种子的相关物理参数进行了标定，在后续的台架试验验证中，试验效果与标定后仿真结果一致。

花生荚果相关物理参数研究主要集中在摩擦系数和恢复系数等物理量。目前对于某一品种花生基本物理参数研究较少，离散元仿真分析在花生机械研发中应用较少。在花生荚果仿真物理模型参数标定及离散元仿真分析方面应用鲜有报道[16]。

本文以宛花2号花生荚果为研究对象，通过物理试验和仿真试验分析相结合的方法确定花生荚果基本物理参数和接触力学参数，并建立花生荚果离散元仿真三维模型[17]。使用EDEM软件(离散元法)对花生荚果颗粒堆积过程进行仿真试验，运用MATLAB软件进行图像处理测定花生荚果的堆积角，利用筛选试验、最陡爬坡试验及响应面试验确定花生荚果仿真参数间的最佳组合，以期为花生荚果采摘、烘干和脱壳等作业仿真提供参数依据。

1 材料和方法

1.1 试验材料与实验仪器

本文选用河南省主要种植品种宛花2号花生荚果为试验物料，试验所用花生荚果经初步清理和筛选，无损伤，无霉变，该花生荚果采用两段式收获方式(先挖掘晾晒，再捡拾收获)。试验主要用实验仪器有DWD型电子万能试验机(额定负荷5 kN、准确度等级0.5、位移分辨率0.01 mm、加载速度0.01～500 mm/min)、数显游标卡尺(精度0.01 mm)、自制摩擦系数测量仪、单轴倾角盒(精度0.05°)、DGF30/7-IA型电热鼓风干燥箱(温度0 ℃～300 ℃，电压220 V)、上皿电子天平(量程100 g，精度0.000 1 g)、FASTEC IMAGING生产的HiSpec5型高速摄像机、50 mL量筒(精度1 mL)、自制铁盒(容积1.047 625×10-3m3)。

1.2 物理参数测定

1.2.1 基本物理参数测定

花生荚果属于散粒体物料，其基本物理特性包括花生荚果外形尺寸(长X×宽Y×厚Z)、密度、含水率、容重、泊松比、弹性模量和剪切模量[18]。随机选取200颗花生荚果平均分为5组。利用量筒、上皿电子天平和数显游标卡尺测定花生荚果的外形尺寸和密度。根据GB/T 20264—2006《粮食、油料水分两次烘干测定法》，利用烘干箱和电子天平进行含水率测量[19-20]。用自制铁盒测定花生荚果的容重。

通过20次重复试验，测得花生荚果的外形尺寸(长X×宽Y×厚Z)、密度、容重和含水率如表1所示。

表1 花生荚果基本物理参数Tab. 1 Basic physical parameters of peanut pods

1.2.2 泊松比

宛花2号花生荚果形状不一，且外形不规则。随机选取10颗花生荚果，记录其长X和宽Y原始尺寸。利用DWD型电子万能试验机对整个花生荚果进行压力变形试验，通过测量花生荚果加载前后开裂极限处宽度方向和长度方向的变形量计算泊松比。试验过程中采用平板压头以10 mm/min速度沿花生荚果长度方向加载，当花生荚果出现破裂时停止加载。由电子万能试验机记录其长度方向变形量，数显游标卡尺记录花生荚果在宽度方向载荷开裂极限处的变形量。通过式(1)计算泊松比，结果取平均值。

(1)

式中：μ——花生荚果泊松比；

εX——宽度方向应变；

εY——长度方向应变；

ΔX——长度方向变形量，mm；

ΔY——宽度方向变形量，mm。

试验重复10次取平均值，测得花生荚果的泊松比为0.4。

1.2.3 弹性模量和剪切模量计算

花生荚果不是规则的圆柱体，若按照常规压缩试验法压缩整个花生荚果测定弹性模量误差较大，试验的难点在于确定压缩过程中花生荚果横截面面积。花生荚果弹性模量主要与荚果外壳相关，因此本试验选取形状近圆柱的花生荚果，并从中截取近圆柱体部分作为试样并量取试样高度，分别量取试样上、中、下三处内外直径，计算面积后取平均值作为花生荚果压缩过程中的横截面积。本文采用管件压缩试验测定弹性模量，试验时将试样放置在电子万能试验机平板上，并设定加载速度和加载时间分别为12 mm/min和10 s，使用平板压头，沿试样轴向压缩10 s后停机，从软件中获取试样压缩过程中载荷、位移相关数据。试验重复10次取平均值，通过式(2)计算得花生荚果试样弹性模量为17.46±5.72 MPa。试样压缩如图1所示。

(2)

式中：E——花生荚果弹性模量，MPa；

FN——压缩过程中的最大力，N；

l——试样原始长度，mm；

Δl——压缩过程中最大力下的变形量，mm；

A——试样压缩受力横截面积，mm2。

图1 花生荚果试样压缩试验

根据试验测定的泊松比和弹性模量，通过式(3)计算可得花生荚果剪切模量为6.5 MPa。

(3)

式中：G——花生荚果剪切模量，MPa。

1.3 接触参数测定

利用EDEM软件进行离散元仿真试验时所需要的接触参数有花生荚果之间和花生荚果与钢板之间的碰撞恢复系数、静摩擦系数和滚动摩擦系数，通过实际物理试验对相关基础参数进行测定，为离散元仿真试验提供参考。

1.3.1 碰撞恢复系数测定

碰撞恢复系数是衡量物体碰撞后恢复到原来形状的能力[18]。碰撞恢复系数定义为碰撞前后两物体沿接触处法线方向上的分离速度与接近速度之比[21-22]。即

(4)

式中：e——物体的恢复系数；

v1、v2——碰撞后物体1和2的速度，m/s；

v10、v20——碰撞前物体1和2的速度，m/s。

(5)

式中：h——花生荚果弹起后最大高度，mm；

H——花生荚果自由落体的高度，mm。

图2 恢复系数测试原理图

根据上述原理，本文采用HiSpec5型高速摄像机采集花生荚果跌落碰撞视频与照片，摄像机设置为150帧/s，试验装置主要包括高速摄像机、跌落架和间距为3 mm的方格背景纸，摄像机镜头与方格背景纸之间的距离为50 cm，跌落碰撞试验装置如图3所示。本文分别以Q235钢板和花生荚果整齐排列的荚果板作为碰撞接触材料，进行花生荚果跌落碰撞试验，用于测取花生荚果与Q235钢板和花生荚果之间的碰撞恢复系数。花生荚果从距离碰撞接触材料40 cm处自由下落，与被测接触材料碰撞后弹起，利用高速摄像机捕捉花生荚果试验过程。

图3 跌落碰撞试验装置

利用高速摄像机控制软件，将跌落碰撞试验的视频以图片形式导出，截取跌落碰撞过程图片，找出花生荚果碰撞后到达最高点处的图片，并数出背景方格纸的方格数，由式(5)计算得出花生荚果的碰撞恢复系数[23]。

利用上述方法分别测取花生荚果与筛板、花生荚果之间的碰撞恢复系数，每个试验重复10次，测取结果取平均值，花生荚果与Q235钢板之间恢复系数0.275±0.064，花生荚果之间恢复系数0.287±0.043。

1.3.2 静摩擦系数测定

本文采用静力学原理对花生荚果摩擦系数进行测定，花生荚果在自制摩擦系数测量仪上的静力学分析如图4所示。由静力学原理可知，在花生荚果受力满足式(6)[24]；当斜面倾角逐渐增大到使花生荚果有下滑趋势时，满足式(7)。

图4 摩擦系数测量原理示意

(6)

μ1=f/N=tanθ

(7)

式中：θ——花生荚果的静摩擦角，(°)；

μ1——花生荚果静摩擦系数；

m——花生荚果质量，kg；

g——重力加速度，9.8 N/s2。

此时斜面倾角为花生荚果最大静摩擦角，利用该方法可以测定花生荚果的静摩擦系数。本试验利用自制摩擦系数测量仪(图5)和单轴倾角盒测量所需花生荚果与各接触材料间的静摩擦系数。测量花生荚果与Q235钢板的静摩擦系数时，将处理好的花生荚果平放在桌面上，用大头针固定花生荚果，先将单轴倾角盒放置于摩擦系数测量仪的铰接轴附近，避免影响花生滑动，再将花生荚果放置于摩擦系数测量仪上，保证测试花生荚果与接触材料面充分接触，缓慢转动摇把增大斜面倾斜角度，当花生荚果沿斜面突然下滑时，摇把停止转动，记录此时单轴倾角盒的数值。该值即为待测花生荚果在测试板上的静摩擦角。测定花生荚果间静摩擦系数时，先将花生荚果粘贴在平面板上，制成荚果板，然后重复上述试验过程，记录静摩擦角，计算荚果间的静摩擦系数。每组试验重复20次计算后取平均值，得到花生荚果与Q235钢板间静摩擦系数平均值为0.583±0.015，花生荚果间静摩擦系数平均值为0.823±0.041。

图5 摩擦系数斜面仪

1.3.3 动摩擦系数测定

动摩擦系数由自制摩擦系数测量仪测定，花生荚果在斜面抬升过程中出现滚动时斜面倾角为花生荚果动摩擦角[25]，每组试验重复20次，计算后取平均值，花生荚果与Q235钢板间动摩擦系数平均值为0.104±0.024，花生荚果间动摩擦系数平均值为0.283±0.043。由于花生荚果表面不规则，且表面粗糙，试验过程中系数变化较大，因此在后续的仿真试验过程中利用堆积角试验对动摩擦系数进行标定。

1.3.4 堆积角测定

本试验选用Q235钢圆筒，根据花生荚果尺寸，确定其内径和高度分别为15 cm和45 cm。试验时圆筒垂直放置在Q235钢板上，花生荚果经过初步筛选后填充满整个圆筒。随后将圆筒以0.05 m/s的速度匀速提升[26]，底板上形成近似锥型的花生荚果堆，试验重复10次。每次试验后拍摄荚果堆正视图像，利用Matlab对花生荚果堆积图像进行灰度处理、二值化处理、图像边界提取；利用Origin对边界数据进行拟合，得到花生荚果堆积角平均值为31.63°，如图6所示。花生荚果单侧堆积角边界拟合如图7所示。

(a) 原图像 (b) 灰度图像 (c) 二值化图像

图7 单侧堆积角边界拟合

2 花生荚果离散元模型建立及仿真参数标定

2.1 离散元仿真模型建立

利用Inventor软件按照花生荚果实际测量所得尺寸建立三维模型(颗粒模型)，将花生荚果模型转成.stl格式后倒入离散元仿真软件EDEM中。对于仿真模型的建立，其尺寸和形状对仿真结果影响较小，颗粒—颗粒、颗粒—钢板之间的相互作用对仿真结果影响较大[17]。本文通过EDEM软件中单球形颗粒对花生荚果模型进行填充，建立花生荚果的仿真模型，如图8所示。

图8 花生荚果仿真三维模型

在EDEM仿真模拟试验中采用Hertz-Mindin模型来模拟花生荚果在干燥仓内的运动特性及物料之间的相互作用。花生荚果仿真模拟堆积角试验模型按照实际堆积角测定装置进行绘制，如图9所示。在导管上方口径处建立颗粒工厂，用于生成花生荚果颗粒，颗粒采用动态生成方式，生成速率为4 800粒/s，一共生成1 200粒，颗粒尺寸采用固定形式。仿真试验总时间为5 s，时间步长为2.88×10-7s，网格尺寸大小为最小颗粒半径的3倍[27]。

图9 堆积角测定装置

仿真开始时，花生荚果颗粒从导管上方口径处的虚拟颗粒平面生成并下落，设置颗粒初始速度为2 m/s，0.25 s后全部花生荚果生成完毕，导管以0.05 m/s的速度向上运动，经过5 s仿真时间后全部花生荚果颗粒在底部底板上静止，形成花生荚果堆积角。

2.2 离散元仿真参数标定

2.2.1 确定显著性影响参数

利用Design-Expert软件进行Plackett-Burman试验设计，试验参数依照实际物理试验测定结果，以实际花生荚果堆积角为目标值，通过Plackett-Burman试验筛选出对目标值具有显著性影响的参数[28]。分别将表2中的8个试验参数的最大、最小值编码为水平+1、-1。Plackett-Burman试验方案及结果如表3所示。

表2 Plackett-Burman 试验参数范围表Tab. 2 Plackett-Burman test parameter range table

表3 Plackett-Burman试验方案及结果Tab. 3 Plackett-Burman test protocol and results

在Design-Expert软件中对试验结果进行方差分析，得到仿真试验参数显著性结果，如表4所示。花生荚果—花生荚果滚动摩擦系数的P值<0.01，对仿真模拟试验有极其显著影响；花生荚果—花生荚果静摩擦系数、花生荚果—钢板静摩擦系数的P值<0.05，对仿真模拟试验有显著影响；其他仿真试验参数的P值>0.05，对仿真模拟试验影响较小。

表4 Plackett-Burman试验参数显著性分析Tab. 4 Significance analysis of Plackett-Burmantest parameters

2.2.2 最陡爬坡试验设计

在Plackett-Burman试验的基础上，对筛选出的3个显著性因素(花生荚果—花生荚果静摩擦系数、滚动摩擦系数，花生荚果—钢板静摩擦系数)进行最陡爬坡试验，以离散元仿真分析堆积角与实际堆积角的相对误差作为评价标准，确定仿真参数最优范围。仿真试验中，花生荚果泊松比取0.4，剪切模量取6.5 MPa；花生荚果—花生荚果碰撞恢复系数取0.275，花生荚果—钢板碰撞恢复系数取0.275，花生荚果—钢板滚动摩擦系数0.125。最陡爬坡试验设计方案及结果如表5所示。仿真试验结果显示，3号试验时，相对误差最小，因此可以确定最优参数范围在3号附近，将3号作为中心点，2、4号分别当作低、高水平进行后续的Box-Behnken响应面试验。

表5 最陡爬坡试验设计方案及结果Tab. 5 Steepest climbing test design scheme and results

2.2.3 Box-Behnken试验设计

以最陡爬坡试验中3号当作中心点(0)，2号、4号分别当作低(-1)、高水平(+1)在Design-Expert软件中进行显著性参数的Box-Behnken试验，设计方案及结果如表6所示。仿真试验中其他非显著性参数均按照最陡爬坡试验中所用参数。

表6 Box-Behnken 试验设计方案及结果Tab. 6 Box-Behnken experiment design scheme and results

在Design-Expert软件中对Box-Behnken试验结果进行二次响应面回归分析，得到离散元仿真试验堆积角二次回归方程

θ=151.578-3 367.919Dx+214.719Ex+

4 273.733Gx+525DxEx-1 444.444DxGx-

344.444ExGx+2 523.056Dx2-771.389Ex2-

2 701.944 4Gx2

Box-Behnken试验方差分析结果如表7所示，由分析结果可知Dx、Ex、DxGx、Dx2、Ex2、Gx2对堆积角影响极其显著；Gx、DxEx对堆积角影响显著；ExGx对堆积角影响不显著。该回归模型P<0.000 1，失拟项P=0.735 2>0.05，表明该回归模型拟合度较好，无失拟现象发生。回归方程决定系数R2=0.987 5，校正决定系数R2=0.971 5，接近于1；回归方程变异系数CV=1.18%。综上分析数据所示，该回归模型极其显著，能够反映真实情况，可用于目标堆积角的预测分析。

表7 Box-Behnken试验回归模型方差分析Tab. 7 Variation analysis of Box-Behnkentest quadratic model

2.2.4 仿真参数标定与试验验证

通过Design-Expert软件中的优化模块，以实际试验堆积角(31.63°)为目标值，对二次回归方程式(8)进行优化求解，从若干最优参数组合中选择一组与实际物理试验测定数堆积相近的参数组合。花生荚果—花生荚果静摩擦系数、花生荚果—花生荚果滚动摩擦系数、花生荚果—钢板静摩擦系数分别为0.74、0.24、0.58，其他非显著性仿真参数取实际物理试验所测定的平均值。

为验证标定后花生荚果仿真参数的准确性，将上述标定后的物理参数作为离散元仿真参数，进行3次仿真模拟试验验证，得到花生荚果堆积角分别为32.29°、33.58°、31.75°，对仿真模拟试验数据与实际试验数据进行独立样本T检验，得到P=0.169>0.05，检验结果表明实际试验堆积角与仿真参数标定后的仿真模拟试验堆积角无显著性差异，且实际试验堆积角平均值31.63°与仿真模拟试验堆积角平均值32.54°的相对误差为2.877%，试验对比如图10所示，验证了仿真模拟试验的准确性。

(a) 实际物理试验

同时对比了其他物料参数标定时所用试验方法及标定结果，如表8所示。进一步表明标定后的参数可为相关研究提供参考。

表8 部分物料试验方法及结果比较Tab. 8 Partial material test methods and comparison of results

3 结论

1) 通过实际试验测得花生荚果基本物理参数，花生荚果外形尺寸为(27.575±2.049) mm×(14.39±0.955) mm×(13.646±0.914) mm，密度为465.08±26.51 kg/m，含水率为(19.28±0.335)%，容重为306.865±9.222 kg/m3，泊松比为0.4，弹性模量为17.46±5.72 MPa，剪切模量为6.5 MPa；采用HiSpec5型高速摄像机和自制摩擦系数测量仪测得花生荚果间的碰撞恢复系数、静摩擦系数、滚动摩擦系数的平均值分别为0.287±0.043、0.823±0.041、0.283±0.043；花生荚果与Q235钢板间的碰撞恢复系数、静摩擦系数、滚动摩擦系数的平均值分别为0.275±0.064、0.583±0.015、0.104±0.024。

2) 以实际试验测得的物理参数作为仿真模拟试验参数选择依据，开展Plackett-Burman试验，筛选出对堆积角存在显著性影响的因素为：花生荚果—花生荚果静摩擦系数、滚动摩擦系数，花生荚果—钢板静摩擦系数，并进一步通过最陡爬坡试验确定显著性因素的取值范围。

3) 开展Box-Behnken试验，建立了堆积角与显著性因素之间的二次回归方程，并以实际物理试验堆积角(31.63°)为目标值对方程进行求解，得到最佳仿真模拟参数：花生荚果—花生荚果静摩擦系数、滚动摩擦系数、花生荚果—钢板静摩擦系数分别为0.74、0.24和0.58。

4) 对试验分析后确定的最佳仿真参数进行仿真模拟试验，对取得的仿真模拟值与实际试验值进行独立样本T检验，得到P=0.169>0.05，表明实际试验堆积角与仿真模拟试验休止角无显著性差异，且相对误差为2.877%。验证了仿真模拟试验的准确性。通过对比其他物料参数标定时所用方法及试验结果，进一步表明标定后的参数可为相关研究提供参考。