包衣棉种物性参数测定与离散元仿真参数标定

2022-05-16贺小伟郭文松王旭峰邢剑飞侯书林

中国农业大学学报 2022年6期

王龙贺小伟胡灿郭文松王旭峰邢剑飞侯书林

(1.中国农业大学工学院，北京 100083； 2.塔里木大学机械电气化工程学院，新疆阿拉尔 843300； 3.自治区教育厅普通高等学校现代农业工程重点实验室，新疆阿拉尔 843300)

2020年，新疆棉花播种面积达到2 500 khm，占全国播种面积的80%左右，覆膜精量播种技术是实现棉花大规模种植的关键。目前普遍使用的棉花播种机具仍存在漏播与重播现象，致使后期需要人工补苗与间苗，增加了人力劳动和经济成本。随着棉花种植向精细化智能化方向发展，需进一步提高棉花穴播器排种性能。

农作物种子大部分为散状颗粒物料，为提高播种机具工作性能，已有研究采用离散元法对取排种机理进行分析进而优化机具结构参数。张涛等利用离散元仿真得到了田间作业振动条件下气吸式排种器种子室内玉米种群的运动规律。赵淑红等设计了一种V型凹槽拨轮式导种部件，并运用EDEM软件分析籽粒在导种管内的滑移状态。胡建平等采用离散元法分析了各因素对磁吸板式精密排种器充种性能的影响。罗伟文等结合Box-Benhnken中心组合试验方法和EDEM离散元技术对小麦播种机碎秸导流装置参数进行设计与优化。可见利用离散元法研究农业颗粒物料在机构中的运动机理可行且高效。

离散元法中仿真模型与参数的准确性是影响仿真结果的重要因素，主要包括物料三维模型以及物料的本征参数(如密度、泊松比和剪切模量等)和接触参数(颗粒之间、颗粒与材料之间的摩擦因数和碰撞恢复系数)。一般物料本征参数直接采用试验所得的真实值，但由于三维模型与真实物料模型存在一定的差异，致使仿真接触参数与真实值存在一定误差，需先标定后才能进行仿真研究。刘彩玲等采用三维扫描法建立水稻离散元模型，发现多球聚合模型仿真精度比常规椭球体模型高。国内外已有研究建立了不同农作物颗粒物料的离散元单元模型，并利用试验及仿真相结合的方法对模型离散元参数进行了标定。将试验与离散元虚拟仿真相结合对农业物料接触参数进行标定能获得准确的仿真参数设定值。

棉花播种机具直接作业对象为包衣棉种，对脱绒棉种进行包衣处理能提高棉种发芽率和抗寒抗病能力，而包衣剂会对棉种表面特征参数造成一定的影响。因此，本研究拟采用逆向工程技术建立棉种离散元颗粒模型；结合物理与仿真落种试验，标定仿真试验中包衣棉种种间接触参数；利用小型棉花精密排种器进行排种试验验证标定的参数，以期为包衣棉种离散元仿真研究提供参考。

1 物理特性参数测定

1.1 三维尺寸及密度

选用南疆地区普遍种植的‘新陆中67号’未包衣及包衣棉种作为试验对象，包衣棉种种衣剂有效成分为福美双和甲基立枯磷，含量分别为10和5 g/100 kg。测得棉种的密度为0.981 g/cm，千粒质量为89.03 g，含水率(湿基)平均值为7.75%。

随机选取无损伤、外形较规则的200粒未包衣棉种，采用数显式游标卡尺(标康BK-318，精度为0.01 mm)测量三轴尺寸(图1)，棉种体积计算公式为：

图1 棉种的三轴尺寸Fig.1 Triaxial dimension of cotton seed

(1)

式中：

为棉种体积，mm；

为棉种长度，mm；

为棉种宽度，mm；

为棉种厚度，mm；

)。统计测量结果可得，棉种长度

为(9.00±0.56) mm，宽度

为(4.84±0.30) mm，厚度

为(4.25±0.28) mm，棉种体积

为(64.90±8.64) mm。

1.2 泊松比及剪切模量

包衣剂对棉种本征参数影响较小，因此直接选用未包衣棉种测定棉种泊松比和剪切模量。采用定义法测定棉种泊松比，即棉种在挤压过程中的横向应变与纵向应变的比值，其计算公式为：

(2)

式中：

为泊松比；

为棉种横向应变；ε为棉种纵向应变；Δ

为棉种横向变形量，mm；

为棉种横向原长度，mm；Δ

为棉种纵向变形量，mm；

为棉种纵向原长度，mm。

采用质构仪(型号TA.XT PlusC，英国Stable Micro Systems公司，测力精度0.000 1 N)测定棉种泊松比，先测得棉种横向与纵向长度后，将其水平放置于承压板上，确保棉种处于承压板中心位置，使用刚性平板(长×宽×厚为26 mm×26 mm×10 mm)进行加载，加载速率为5 mm/min，触发传感器记录的起始力为0.3 N。由于棉种在加载1～2 mm时会发生破裂，并伴随有破裂声音，棉种破裂时停止加载，并测量棉种横向与纵向长度。每组试验重复20次取平均值，可得新陆中67号棉种泊松比为0.27。

棉种的剪切模量可以利用弹性模量和泊松比求得。棉种剪切模量的计算公式为：

(3)

式中：

为棉种剪切模量，MPa；

为棉种弹性模量，MPa；

为棉种泊松比。

棉种为农业颗粒物料，在受到挤压时，内部应力分布较为复杂。根据ASAE S368.4 DEC2000 (R2017)《Compression Test of Food Materials of Convex Shape》标准可知，当对棉种进行刚性平板挤压时，其弹性模型计算公式为：

(4)

式中：

为棉种所受压缩接触力，N；

为棉种形变量，mm；

为棉种接触处最小曲率半径，mm；

′为棉种接触处最大曲率半径，mm；

为系数，量纲为1。

棉种曲率半径采用图像处理技术测定，图像采集装置主要包括：工业相机(型号MV-GE131GC-T，深圳市迈德威视科技有限公司，122万像素，CMOS传感器)，高清镜头(型号MV-LD-4M-G，深圳市迈德威视科技有限公司，焦距8 mm)，环形光源，支架等。将棉种水平放置于支架平面上，用工业相机对图像进行采集，利用MATLAB软件对图像进行处理，包括二值化、膨胀、腐蚀、边缘检测、轮廓提取和拟合等，可得棉种轮廓曲线及拟合曲线(图2)。棉种外形轮廓类似椭圆，基于最小二乘法将棉种轮廓进行拟合，得到拟合椭圆一般方程式为：

Bxy

(5)

式(5)一阶、二阶均可导，因此

(

)在(

(

))处的曲率半径为：

(6)

图2 棉种轮廓图像处理过程Fig.2 Image processing of cotton seed contour

本研究采用质构仪XT PlusC测定棉种弹性模量。先处理棉种的图像信息，根据式(5)和式(6)得到最小曲率半径和最大曲率半径后，将棉种水平放置在承压板上，确保其处于承压板中心位置，使用刚性平板进行加载，加载速率为10 mm/min，触发传感器记录的起始力为0.3 N，加载位移为2 mm。每组试验重复20次取平均值，可得‘新陆中67号’棉种剪切模量为14 MPa。

1.3 自然休止角

用排出法测定棉种休止角。棉种休止角装置采用透明有机玻璃制成，主要由上种箱、下种箱和挡板组成(图3(a))，上种箱和下种箱的长×宽×厚均为400 mm×200 mm×50 mm，挡板宽度为40 mm。有机玻璃的泊松比为0.5，密度为1 180 kg/m，剪切模量为177 MPa。试验时，先清理种箱，确保种箱清洁后插入挡板；在上种箱中均匀地加入适量棉种，约占上种箱总容积的3/4，且棉种群上表面尽量保持平齐；将挡板快速抽出，棉种在重力作用下沿开口处下滑，落入下种箱，待种箱中棉种稳定后，棉种堆积斜面与水平面的夹角即为棉种休止角

(图3(b))。

1.上种箱；2.棉种；3.挡板；4.下种箱 1. Seed feeding box; 2. Coated cotton seed; 3. Baffle; 4. Seed box γ为棉种休止角，(°)。γ is angle repose of cotton seed，(°).图3 棉种休止角测定试验装置及其试验结果Fig.3 Test device and result for measuring repose angle of cotton seeds

为减少人为测量操作过程中产生的误差，采用MATLAB软件对采集到的休止角图像进行处理，包括中值去噪、灰度、二值化及删除孤岛区域等，得到二值化图像(图4(b))；提取种群边界轮廓得到种群边界曲线，利用最小二乘法对边界曲线进行拟合(图4(c))，拟合直线的表达式为：

(7)

式中：

为图像垂直像素点值；

为图像水平像素点

图4 棉种休止角测定图像处理过程Fig.4 Image processing of repose angle measurement of cotton seed

值；

为拟合直线斜率；

为拟合直线截距。利用拟合直线的斜率即可求得棉种休止角

，计算公式为：

=artan

(8)

每组试验重复10次取平均值，可得未包衣、包衣棉种的自然休止角分别为36.79°和37.34°。表明脱绒处理后的棉种表面比较光滑，其流动性也较好，而包衣处理后，棉种表面粘附着包衣剂，增大了棉种表面的摩擦。

1.4 碰撞恢复系数

根据牛顿力学定律可知，碰撞恢复系数

为棉种碰撞被测材料前后的速度比值，其计算公式为：

(9)

式中：

为棉种碰撞被测材料前的速度，m/s；

为棉种碰撞被测材料后的速度，m/s；

为重力加速度，9.81 m/s；

为棉种碰撞后的回弹高度，m；

为棉种碰撞前的下落高度，m。

棉种碰撞恢复系数试验装置主要由华为mate30、有机玻璃板与坐标纸组成。试验时，将棉种从一定高度自由下落，与有机玻璃板碰撞后弹起，利用华为mate30慢动作视频拍下棉种整个碰撞与运动过程，慢动作视频设置为240帧/s；读取每一帧图片并记录棉种下落点、碰撞点和回弹最高点对应的位置尺寸，得到棉种的下落高度和回弹高度；根据公式(9)可得棉种与有机玻璃板之间的碰撞恢复系数。每组试验重复10次取平均值，可得未包衣、包衣棉种与有机玻璃之间的碰撞恢复系数分别为0.26和0.25。

1.5 摩擦因数

利用斜面仪搭建棉种与有机玻璃之间的摩擦因数试验装置(图5)。试验时，先调平底座，在斜面仪测试平面装上有机玻璃板，调节测试平面使角度指针指向0刻度。测定静摩擦因数时，为防止包衣棉种滚动，降低试验误差，将三粒尺寸较规则棉种粘结在一起，静置于斜面仪测试平面，匀速转动手摇轮，待粘结棉种开始下滑时停止转动，此时指针所指的角度值为滑动摩擦角，滑动摩擦角的正切值即为静摩擦因数。每组试验重复10次取平均值，可得未包衣、包衣棉种与有机玻璃之间的静摩擦因数分别为0.48和0.49。

1.支架；2.滑轮；3.棉种；4.手摇轮；5.斜面仪 1.Bracket; 2.Pulley; 3.Cotton seed; 4.Hand wheel; 5.Inclinometer图5 摩擦因数测定装置Fig.5 Determination device of friction coefficient

测定滚动摩擦因数时，将单粒棉种置于测试平面，确保棉种长轴与斜面倾斜方向垂直；匀速转动手摇轮，棉种沿斜面完全滚落时，记录指针所指的角度值为滚动摩擦角，滚动摩擦角的正切值即为滚动摩擦因数。每组试验重复10次取平均值，可得未包衣、包衣棉种与有机玻璃间的滚动摩擦因数分别为0.20和0.21。包衣棉种与有机玻璃的静摩擦因数和滚动摩擦因数均大于未包衣棉种。表明棉种包衣会对棉种表面摩擦特性产生影响，包衣棉种表面粘附的包衣剂增大了其表面粗糙度。

2 离散元模型建立与接触参数标定

2.1 包衣棉种离散元模型建立

为建立精确的包衣棉种仿真模型，本研究基于逆向工程技术，采用Capture MINI扫描仪扫描包衣棉种外形轮廓，获得包衣棉种点云数据(图6(a))；使用Geomagic Wrap 3D软件对点云数据进行去噪、压缩、封装与曲率修复处理后得到包衣棉种轮廓模型(图6(b))。在离散元仿真软件中，需要用球形颗粒来建立物料的离散元颗粒模型，而实际物料大多为不规则体，一般采用粘结颗粒模型(Bonded Particle Method)或多球聚合模型(Multi-sphere Method)建立物料模型。其中多球聚合模型对计算机配置要求不太高，仿真计算时间适中，综合考虑本研究采用多球聚合模型进行仿真试验。将包衣棉种轮廓模型导入EDEM软件中，依据棉种外形轮廓用13粒直径不同的球形颗粒重叠堆积建立包衣棉种多球聚合模型(图6(c))。

图6 包衣棉种仿真模型建立过程Fig.6 Establishment process of discrete element model of coated cotton seed

2.2 接触参数标定

在离散元仿真过程中，材料的本征参数和接触参数十分重要，其中本征参数包括材料的密度、泊松比、弹性模量等，接触参数包括材料之间的碰撞恢复系数、静摩擦因数和滚动摩擦因数。棉种形状较不规则，包衣棉种种间的接触参数无法直接采用试验测定，因此采用物理试验与仿真试验相结合的方法进行标定。根据棉种物理特性参数、包衣棉种的碰撞恢复系数和摩擦因数测定结果，确定仿真试验中棉种本征参数和包衣棉种与有机玻璃之间的接触参数(表1)。

根据棉种颗粒之间的粘结性质，仿真过程中颗粒接触模型选取Herz-Mindlin (no slip)接触模型。为了更真实地模拟实际试验过程，在仿真试验中棉种颗粒模型体积按试验统计的体积分布值生成，棉种颗粒数量、试验过程尽量与物理试验一致。仿真试验结束后采集棉种休止角图像，同样采用MATLAB软件处理图像，可得仿真休止角

′。休止角相对误差

计算公式为：

(10)

式中：

为物理试验所得包衣棉种自然休止角，由1.3节可知，为37.34°；

′为仿真试验所得的休止角，(°)。

表1 棉种本征参数和包衣棉种与有机玻璃接触参数
Table 1 Intrinsic parameters of cotton seed and contact parameters between coated cotton seed and plexiglass

参数Parameter数值Value棉种密度/(g/cm3)Density of cotton seed0.981棉种泊松比Poisson’s ratio of cotton seed0.27棉种剪切模量/MPaShear modulus of cotton seed14包衣棉种与有机玻璃碰撞恢复系数Restitution coefficient between coated cotton seed and plexiglass0.25包衣棉种与有机玻璃静摩擦因数Static friction coefficient between coated cotton seed and plexiglass0.49包衣棉种与有机玻璃滚动摩擦因数Rolling friction coefficient between coated cotton seed and plexiglass0.21

在仿真试验过程中，以包衣棉种种间碰撞恢复系数

、静摩擦因数

和滚动摩擦因数

为待标定参数，休止角相对误差为试验指标，进行最陡爬坡和中心组合试验。根据大量的仿真预试验结果及相关文献，确定包衣棉种种间接触参数的最陡爬坡试验设计方案，进行仿真试验后得到休止角和休止角相对误差(表2)。

休止角相对误差值随试验因素值的增大呈先减小后增大的趋势，在第4组时达到最小，为2.17%，表明棉种种间接触参数的最优组合在第3组和第5组之间。因此选取第3组、第4组和第5组试验因素确定仿真试验因素编码(表3)，根据中心组合仿真试验方案进行仿真试验得到休止角相对误差见表4。

表2 最陡爬坡试验设计方案及其仿真试验结果
Table 2 Scheme and results of steepest ascent experiment

序号No.试验因素 Test factorX1X2X3休止角/(°)Repose angle休止角相对误差/%Relative error of repose angle10.070.150.0629.0122.3120.100.180.0832.3113.4730.130.210.1035.385.2640.160.240.1238.152.1750.190.270.1439.365.4260.220.300.1643.7517.1770.250.330.1845.3221.37

注：、、分别为包衣棉种种间碰撞恢复系数、静摩擦因数和滚动摩擦因数，表3、表4、表5同。

Note: ， and respectively for the restitution coefficient, static friction coefficient and rolling friction coefficient between coated cotton seeds. The same as in the
Table 3,
Table 4 and
Table 5.

表3 仿真试验因素编码
Table 3 Factors and codes of simulation experiment

编码Code试验因素 Test factorX1X2X3-1.680.110.190.09-10.130.210.10 00.160.240.12 10.190.270.14 1.680.210.290.15

根据仿真试验结果，运用Design-Expert 12.0软件进行多项式回归分析，可得包衣棉种种间接触参数对休止角相对误差影响的回归模型为：

55-0

934 5

487 9

+ 0

302 5

-1

+ 0

617 7

(11)

对回归模型进行显著性检验及方差分析(表5)，可得回归模型的拟合度极显著(

<0.01)。包衣棉种种间碰撞恢复系数和静摩擦因数的交互项(

)对休止角相对误差的影响不显著(

>0.05)，其他各项对模型的影响均显著，表明相关接触参数对响应值的影响是非线性关系，存在二次关系。模型的失拟项

=0.229 1>0.05，不显著，表明模型无其他影响指标的主要因素存在。回归方程的拟合优度

和Adjusted

分别为0.974 2和0.951 0，表明回归模型方程的预测值与实际值拟合度较好，自变量对因变量的解释程度较高。

以休止角相对误差的最小值为优化目标，对回归方程进行寻优求解，得到目标函数及非线性约束条件方程组为：

(12)

表4 仿真试验中心组合试验设计方案与结果
Table 4 Scheme and results of the central composite design experiment

试验序号Test No.试验因素水平值 Test factor level valueX1X2X3休止角相对误差/%Relative error of repose angle1-1(0.13)-1(0.21)-1(0.10)5.2621(0.19)-1-14.833-11(0.27)-14.51411-15.985-1-11(0.14)6.8361-110.937-11110.8081115.429-1.68(0.11)0(0.24)0(0.12)4.14101.68(0.21)002.64110(0.16)-1.68(0.19)06.361201.68(0.29)010.571300-1.68(0.09)3.6514001.68(0.15)5.59150000.94160001.46170002.17180001.69190000.97200002.05

注：括弧内为试验因素值。

Note: The values in brackets are the test factor values.

表5 休止角相对误差回归模型的方差分析
Table 5 Variance analysis of regression model of repose angle relative error

方差来源Soruce of variation平方和Sum of square自由度Degree of freedomF值F-valueP值P-value模型 Model156.44941.97<0.000 1** X111.93128.800.000 3** X218.61144.92<0.000 1** X33.2517.850.018 7* X1X20.7311.770.213 2 X1X318.97145.81<0.000 1** X2X38.12119.610.001 3** X125.50113.280.004 5** X2283.841202.44<0.000 1** X3215.97138.550.000 1**残差 Residual4.1410失拟项 Lack of fit2.7752.020.229 1误差 Error1.375总和 Sum160.5819

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

Note: ** and * indicate significance at <0.01 and <0.05, respectively.

根据方程组(12)可得最佳包衣棉种种间接触参数组合，包衣棉种种间碰撞恢复系数

、静摩擦因数

和滚动摩擦因数

分别为0.19、0.23和0.13。利用该最佳参数组合进行仿真休止角试验，得到休止角相对误差为0.72%，将包衣棉种物理试验与仿真试验结果进行对比(图7)，可知休止角堆形上相似度较高，表明该最佳参数组合下的仿真结果准确可靠，可为后续EDEM仿真试验研究过程中的接触参数设定提供参考。

图7 包衣棉种物理试验与仿真试验结果对比Fig.7 Result comparison between physical experiment and simulation experiment of coated cotton seed

3 试验验证

为进一步验证包衣棉种离散元模型和仿真参数的可靠性，搭建小型棉花精密排种器排种试验装置(图8)。选取‘新陆中67号’包衣棉种进行台架试验，取种轮采用3D打印，窝眼孔的直径为7 mm，深度为10 mm。将排种器简化模型、包衣棉种离散元模型和标定所得最佳接触参数组合导入离散元软件进行仿真试验(图9)。

1.电源；2.取种轮；3.种箱；4.CL57C驱动器； 5.57闭环步进电机；6.步进电机控制器 1.Power; 2.Seed picking wheel; 3.Seed box; 4.CL57C drive; 5.57 closed loop stepper motor; 6.Stepper motor controller图8 棉种排种试验装置Fig.8 Cotton seed metering experiment device

图9 包衣棉种排种仿真试验Fig.9 Seed metering simulation experiment of coated cotton seed

根据GB/T 6976—2005《单粒(精密)播种机试验方法》，窝眼孔中充入1粒种子即为合格，窝眼孔中未充入种子即为漏播。选取取种轮转速为10、15、20、25和30 r/min进行试验，连续测量100个窝眼孔，每组试验重复3次取平均值，可得不同取种轮转速下各试验指标的台架试验和仿真试验结果(表6)，取种轮充种合格率和漏播率相对误差平均值分别为1.34%和3.93%，均小于5%，表明该包衣棉种离散元颗粒模型和接触参数可用于离散元仿真试验。

表6 不同转速下排种试验的合格率和漏播率
Table 6 Qualified rate and missed seeding rate of seed metering experiment at different speeds

转速/(r/min)Speed合格率/% Qualified rate漏播率/% Missed seeding rate仿真试验Simulationexperiment台架试验Benchexperiment相对误差Relativeerror仿真试验Simulationexperiment台架试验Benchexperiment相对误差Relativeerror1084.6883.731.137.657.943.651582.8681.671.469.589.154.702079.5580.421.0812.4712.923.482576.4375.181.6615.3814.843.643070.8271.791.3520.5719.744.20