文恩杨１，李玉华１，牛子孺１，侯加林１，２
（１．山东农业大学机械与电子工程学院，山东泰安 ２７１０１８；２．山东省农业装备智能化工程实验室，山东泰安 ２７１０１８）
摘 要：为确定蒜种离散元模型仿真参数，采用多球聚合的方法建立了蒜种离散元模型，并对可实现测量的模型参数（蒜种－钢板恢复系数、蒜种－蒜种恢复系数、蒜种－钢板静摩擦因数、蒜种－蒜种静摩擦因数）进行了物理试验测定。对不易测量的模型参数（蒜种－钢板滚动摩擦系数、蒜种－蒜种滚动摩擦系数）进行了仿真试验标定，应用、建立了两种不易测量的接触参数与蒜种堆休止角、密相区域圆半径的二次回归模型，最终求解寻优得到参数最优值。通过上述试验方法，得到蒜种离散元模型参数：蒜种－钢板恢复系数为０．５１１，蒜种－蒜种恢复系数为０．４８７，蒜种－钢板静摩擦因数为０．４７３，蒜种－蒜种静摩擦因数为０．５０３，蒜种－钢板滚动摩擦因数为０．１１１，蒜种－蒜种滚动摩擦因数为０．１０８。最后，在得到的接触参数下进行仿真验证试验，结果表明：蒜种堆休止角、密相区域圆半径的仿真值与实际值的相对误差分别为１．３６％、１．５０％，无明显差异，排种功率、排种速度的仿真与试验值在整体变化趋势上具有一致性，验证了蒜种离散元模型与仿真试验的有效性，可为大蒜播种机械仿真设计与优化提供参考。
关键词：离散元；蒜种；参数标定；休止角；密相区域圆半径
中图分类号：Ｓ１８３ 文献标识码：Ａ文章编号：１００３－１８８Ｘ（２０２１）０５－０１６０－０８
０ 引言
目前，我国大蒜几乎全为人工种植，面临着劳动强度大、人工成本高、播种效率低等问题，迫切需要机械化作业［１－３］。由于传统试验方法的限制，蒜种在播种机械各作业流程中的运动并未得到系统的理论分析，蒜种与播种机械作用机理的模糊性为大蒜播种机械化的发展带来了局限性。近年来，计算机技术与离散单元法的进步发展为明晰蒜种在播种机械中的动力学行为带来了可能，为大蒜播种机械的设计及优化提供有力的理论依据；然而，不正确模型参数值的采用将导致仿真模型发生失真现象，则建立准确有效的蒜种模型对于离散单元法的应用尤为重要。
目前，国内外学者在农业物料颗粒参数标定方面做了很多研究：王云霞等基于休止角试验仿真标定了玉米种子间静摩擦因数和滚动摩擦因数［４］；韩燕龙采
收稿日期：２０１９－１０－１０
基金项目：山东省农机装备研发创新计划项目（２０１７ＹＦ００１－０９）；国家现代农业产业技
术体系特蔬菜创新团队项目（ＣＡＲＳ－
２４－Ｄ－０１）；山东省农业重大应用技术创新项目
（ＳＤ２０１９ＮＪ００４）
作者简介：文恩杨（１９９４－），男，山东聊城人，硕士研究生，（Ｅ－ｍａｉｌ）２４０５２１００５０＠ｑｑ．ｃｏｍ。
通讯作者：侯加林（１９６３－），男，山东高密人，教授，博士生导师，（Ｅ－ｍａｉｌ）ｊｌｈｏｕ＠ｓｄａｕ．ｅｄｕ．ｃｎ。用系数标定方法测量了米粒间的滚动摩擦因数［５］；刘文政等利用试验与模拟相结合的方法标定了微型马铃薯仿真参数［６］；袁全春等通过休止角物理试验标定了机施有机肥散体颗粒的模型参数［７］，而对于蒜种颗粒离散元仿真参数标定方面的研究未见报道。蒜种相对于上述类球体或类椭球体的物料，具有形状不规则且体积较大的特征，目前尚未有测定蒜种滚动摩擦因数的有效方法，因此蒜种接触参数的标定成为建模过程中的难题。
以山东苍山大蒜为研究对象，以多球聚合的方法建立蒜种离散元模型，选取“Ｈｅｒｚ－Ｍｉｎｄｌｉｎ（ｎｏｓｌｉｐ）”接触模型，建立物理试验测定模型与仿真试验标定模型，并应用ＣｅｎｔｒａｌＣｏｍｐｏｓｉｔｅＤｅｓｉｇｎ建立了两种不易测量的接触参数与蒜种堆休止角、密相区域圆半径
的二次回归模型，最终求解寻优得到参数最优值。同时，并结合蒜种堆体试验与排种试验对蒜种离散元模型参数进行验证，通过试验与仿真误差分析验证标定方法的有效性，为轮廓复杂物料的离散元参数标定提供参考。
１ 参数标定模型建立
１．１ 蒜种模型
选取山东苍山大蒜作为蒜种离散元模型的原型，
该品种产量高，蒜头为四六瓣，具有蒜种特征的代表性。为保证蒜种出芽的一致性及蒜苗质量，蒜农在播种前通常会对蒜种进行有条件地筛选。为此，按照“蒜形饱满、大小相近”的原则，挑选出５００粒蒜种进行外形特征的统计，蒜种外形特征如图１所示。针对蒜种的长、宽、厚及鳞芽长度等特征，通过数显游标卡尺（型号ＭＮＴ－１５０Ｔ，精度０．０１ｍｍ）对蒜种样本进行测量并统计，得出其平均值。蒜种外形特征尺寸如表１
un38.3
所示。
图１ 蒜种外形特征图
Ｆｉｇ．１ Ａｐｐｅａｒａｎｃｅｃｈａｒａｃｔｅｒｉｓｔｉｃｓｏｆｇａｒｌｉｃｓｅｅｄ
表１ 蒜种外形特征尺寸
Ｔａｂｌｅ１ Ｃｈａｒａｃｔｅｒｉｓｔｉｃｄｉｍｅｎｓｉｏｎｓｏｆｇａｒｌｉｃｓｅｅｄ
ｍｍ
蒜体长度
Ｌ
蒜体宽度
ｂ
蒜体厚度
ｄ
鳞芽长度
Ｎ
蒜种总长
３４．５０１７．９４１６．２９１７．０７５１．５７
根据蒜种外形特征统计数据，建立蒜种离散元模型，蒜种作为农业散体物料，其形状不规则，外形轮廓复杂。为提高蒜种离散元模型参数的准确性，在ＥＤＥＭ中采用多球聚合的方法建立模型，蒜种模型已接近真实物料外形，如图２
所示。
图２ 蒜种离散元模型及实物图
Ｆｉｇ．２ Ｄｉｓｃｒｅｔｅｅｌｅｍｅｎｔｍｏｄｅｌｏｆｇａｒｌｉｃｓｐｅｃｉｅｓａｎｄｔｈｅｒｅａｌｐｉｃｔｕｒｅ
仿真试验过程中，设备材料选用农业机械中常用的Ｑ２３５钢，通过试验和查阅相关文献［８］，得到蒜种与Ｑ２３５钢的本征参数，如表２所示。
表２ 蒜种、Ｑ２３５钢本征参数
Ｔａｂｌｅ２ ＩｎｔｒｉｎｓｉｃｐａｒａｍｅｔｅｒｏｆｇａｒｌｉｃｓｅｅｄａｎｄＱ２３５ｓｔｅｅｌ
材料泊松比
剪切模量／ＭＰａ
密度
／ｇ·ｃｍ
－３
蒜种０．２３２３．８２
１．０８Ｑ２３５钢
０．２８
８．２×１０
４
７．８５
１．２ 仿真接触模型
蒜种种皮表面较为干燥，表面黏着力较小，且假设颗粒在堆积过程中位移、力、速度等参数的变化是通过颗粒与接触体之间产生的微小交叠量值的不同确定。基于以上假设，为准确模拟蒜种间及蒜种与材料间的力学运动，采用Ｈｅｒｔｚ－Ｍｉｎｄｌｉｎ（ｎｏｓｌｉｐ）模型［９］模拟仿真蒜种间及蒜种与机器间的相关运动过程。其中，法向力分量基于Ｈｅｒｔｚ接触理论，切向力模型基
于Ｍ
ｉｎｄｌｉｎ－Ｄｅｒｅｓｉｅｗｉｃｚ的理论。２ 蒜种离散元模型参数标定
２．１ 蒜种离散元模型参数物理试验测定２．１．１ 蒜种－钢板弹性恢复系数测定
利用蒜种的弹性碰撞试验［１０］
对蒜种－钢板间的
恢复系数进行测定，为简化计算，将试验设计为钢板水平静置，蒜种在钢板正上方自由落体。假设钢板在碰撞前后速度均为零，则蒜种－钢板间的弹性恢复系数表达式为
ｅ
十七大报告＝ｕｖ＝２ｇＨ槡１２ｇＨ槡
０
＝
Ｈ１
Ｈ０槡
（１）
式中 Ｈ１
—蒜种颗粒碰撞后第一次弹跳最大上升高度；
Ｈ０
—蒜种初始下落高度。由于蒜种自身质量分布不均匀，在下落的过程中会有轻微的偏转，为减少蒜种外形对试验的干涉，尽量达到蒜种与钢板点接触的理想状态，试验前需对蒜种进行预处理。将蒜种的尾部及蒜尖处切除少许，剩余的蒜种部分呈类扇形柱体，在进行弹性碰撞试验时，使蒜种初始状态为背部弧面向下，使其自由落体至水平钢板上，利用高速摄像机（德国ＯｐｔｒｏｎｉｓＧｍｂＨ制造，ＣＲ３０００ ２）记录下蒜种颗粒碰撞后第一次弹跳最大上升高度Ｈ１
。
为最大限度保证蒜种在下落过程中姿态的一致性，得到准确的试验结果，选取蒜种初始下落高度为２００ｍｍ，取３０粒蒜种，每粒蒜种进行５次试验，分别记录颗粒碰撞后第一次弹跳最大上升高度Ｈ１，如图３所示。在５组试验中选取ｈ最大值作为该粒蒜种的试验
粉末冶金网
值，将３０粒蒜种的试验值求平均值，得试验结果Ｈ１为５２．１６ｍｍ，利用公式计算得到弹性碰撞恢复系数为０．５１１，选取此值作为蒜种－钢板间仿真恢复系数，以待
验证。
（ａ） 初始 （ｂ） 碰撞 （ｃ） 上升
图３ 蒜种－钢板间恢复系数测定试验
Ｆｉｇ．３ Ｍｅａｓｕｒｅｍｅｎｔｔｅｓｔｇａｒｌｉｃｓｅｅｄ－ｓｔｅｅｌｐｌａｔｅｃｏｅｆｆｉｃｉｅｎｔｏｆｒｅｓｔｉｔｕｔｉｏ
２．１．２ 蒜种－蒜种间弹性恢复系数测定
测定蒜种间的弹性恢复系数，需要将蒜种与蒜种进行碰撞，观察到蒜种侧面平坦，在蒜种相互碰撞过程中，可以最大程度实现蒜种沿法线方向反弹，故以此为试验改进依据，在蒜种－钢板弹性恢复系数
测定模型的基础上改进蒜种－蒜种弹性恢复系数测定模型，如图４所示。在钢板上固定高为１０ｍｍ的圆柱器皿，在器皿内固定一粒蒜种ａ，使蒜种ａ侧面与器皿上边缘在同一水平面内，按照上文所述方法，按两粒１
组取３０组蒜种，并对其中１粒蒜种ｂ进行预处理；将处理后的蒜种ｂ静置于蒜种ａ正上方２００ｍｍ处并使其自然垂落，每组蒜种进行５次试验。利用高速摄像机记录颗粒ｂ碰撞后第一次弹跳最大上升高度Ｈｂ，在５组试验中选取Ｈｂ最大值作为该粒蒜种的试验值，将３０组蒜种的试验值求平均值，得试验结果Ｈｂ为４７．４３ｍｍ。利用式（１）计算得弹性碰撞恢复系数为０．４８７，选取此值作为蒜种－钢板间仿真恢复系数，
以待验证。
（ａ） 初始 （ｂ） 碰撞 （ｃ） 上升
图４ 蒜种－蒜种恢复系数测定试验
Ｆｉｇ．４ Ｍｅａｓｕｒｅｍｅｎｔｔｅｓｔｏｆｇａｒｌｉｃｓｅｅｄ－ｇａｒｌｉｃｓｅｅｄｃｏｅｆｆｉｃｉｅｎｔｏｆｒｅｓｔｉｔｕｔｉｏｎ
２．１．３ 蒜种－钢板静摩擦因数测定
采用斜面法测量蒜种－钢板静摩擦因数，原理：质量为ｍ的蒜种的重力可分解为两个力，
一个为平行于斜面的力Ｆ，另一个为垂直于斜面的力Ｎ。当斜面的倾角α很小时，
Ｆ小于蒜种与斜面间的静摩擦力ｆ，蒜种在斜面上保持静止状态；随着斜面倾角α的增大，Ｆ也越来越大；当Ｆ大于滑块与斜面间的静摩擦力ｆ时，蒜种开始运动并顺着斜面下滑。
在蒜种开始滑动的瞬间，有
ｆ＝Ｆ＝ｍｇｓｉｎα
强度极限
（２）
Ｎ＝ｍｇｃｏｓα（３）ｆ＝μＮ
（４）
因此，蒜种与钢板间静摩擦因数μ＝
ｆ
Ｎ
＝ｔｇα。试验时，将钢板置于水平状态，钢板上静置有单粒蒜种，缓慢旋转手柄，使钢板缓慢上升，防止振动。当蒜种开始下滑时，利用数显角度尺（精度０．０５°，量程０°～３６０°）记下此时的平面倾角α，重复试验３０次，取平均值α
＝３０．３０°，进而计算得蒜种与钢板间静摩擦因数μ＝０．４７３；选取蒜种－钢板间仿真静摩擦因数为
０．４７３，以待验证。
２．１．４ 蒜种－蒜种静摩擦系数测定
由于蒜种表面轮廓复杂且不规则，通过蒜种体来测量蒜种与蒜种间静摩擦因数显得尤为不易。观察发现：苍山蒜种侧面表皮表面平坦且易提取，故将蒜种侧面表皮切割成小矩形块，用双面胶粘贴排列在玻璃板上，将玻璃板固定在斜面仪上，蒜种置于蒜种表皮面一端，缓慢旋转手柄，使钢板缓慢上升，防止振动；当蒜种开始下滑时，记下此时的平面倾角α，重复
试验３０次，取平均值α＝２６．７０，进而计算得蒜种间静摩擦因数μ＝０．５０３。选取蒜种之间仿真静摩擦因数为０．４７３，
以待验证。
图５ 蒜种－蒜种静摩擦因数测定试验
Ｆｉｇ．５ Ｍｅａｓｕｒｅｍｅｎｔｔｅｓｔｏｆｇａｒｌｉｃｓｅｅｄ－ｇａｒｌｉｃｓｅｅｄｃｏｅｆｆｉｃｉｅｎｔｏｆ
ｓｔａｔｉｃｆｒｉｃｔｉｏｎ
２．２ 蒜种离散元模型参数仿真实验标定２．２．１ 堆体试验
休止角是表征颗粒物料流动、摩擦等特性的宏观参数，其静摩擦因数、滚动摩擦因数等接触参数对结果影响显著。本文采用空心筒法进行休止角测量试验，试验装置由空心筒（直径１３０ｍｍ，筒高２６０ｍｍ，底板长５００ｍｍ，宽５００ｍｍ）组成，装置材料为Ｑ２３５钢。在空心筒中填充３００粒蒜种，静置后，以０．０５ｍ／ｓ的速度缓慢提升空心筒；蒜种从圆筒底部逐渐漏出形成休止角，待蒜种稳定后，蒜种堆斜面与底板形成夹角即为蒜种堆的休止角，使用高清照相机采集测量试验的蒜种堆体图像，利用ＭａｔＬａｂ对图像依次进行灰度处理、二值化处理、孔洞填充，提取轮廓曲线；最后，对轮廓曲线进行直线拟合，如图６（ｂ）所示。其中，横纵坐标分别为图像水平像素点和垂直像素点，不具有实际量纲。拟合直线斜率即为蒜种堆的休止角正切值，继而得到蒜种休止角。该试验重复进行１０次，取平均
值，最终得蒜种堆休止角为２
６．５１°
。（ａ）
实物图
（ｂ） 拟合直线图图６ 休止角图像处理
Ｆｉｇ．６ ＩｍａｇｅａｎａｌｙｓｉｓｂｙＭａｔｌａｂｆｏｒｒｅｐｏｓｅａｎｇｌｅ
为明晰滚动摩擦因数对蒜种堆积体的影响，引入密相区域圆半径作为另一个评价指标，在进行蒜种堆
体试验时发现，堆体俯视投影面不是连续的圆面，而是在堆体边缘蒜种比较分散。其中，靠近边缘区域是连续的圆面，且堆体边缘分散蒜种数目约占总数的１０％，将包含９０％蒜种总数的圆作为密相区域圆，测量其半径ｒ，计算１０次试验平均值，得密相区域圆半径为１４１．２７ｍｍ。蒜种堆俯视图如图７
所示。
图７ 蒜种堆俯视图Ｆｉｇ．７ Ｖｅｒｔｉｃａｌｖｉｅｗｏｆｇａｒｌｉｃｐｉｌｅ
２．２．２ 仿真试验标定
在ＳｏｌｉｄＷｏｒｋｓ中建立钢制圆筒和底板的三维模型，并导入ＥＤＥＭ２０１８中作为结构模型，将建立的蒜种离散元模型作为颗粒模型。仿真试验时，设置空心筒提升速度５０ｍｍ／ｓ，生成颗粒总量为３００个，生成速率１５００个／ｓ，固定时间步长取瑞利时步的２０％，数据保存时间间隔０．０１ｓ，网格尺寸取２倍最小球形单元尺寸。蒜种堆体试验通过改变蒜种间滚动摩擦因数及蒜种－钢板间滚动摩擦因数，观察休止角与密相区域圆半径，其余接触参数均采用上述测定试验获得的标定值。
经过大量预仿真试验，确定蒜种间滚动摩擦因数仿真试验水平范围为０．１００～０．１２０，蒜种－钢板间滚动摩擦因数仿真试验水平范围为０．１００～０．１２０。根据此范围进行试验［１１］，堆体仿真试验因素的水平值设置如表３所示。各取５个水平，试验方案及结果如表４所示。共进行１３次试验，包括５次中心点的重复试验。
表３ Ｃｅｎｔｒａｌｃｏｍｐｏｓｉｔｅｄｅｓｉｇｎ因素水平表
Ｔａｂｌｅ３ Ｆａｃｔｏｒｓａｎｄｌｅｖｅｌｓｔａｂｌｅｏｆｃｅｎｔｒａｌｃｏｍｐｏｓｉｔｅｄｅｓｉｇｎ编码
因素
Ａ（蒜种间滚动摩擦因数）Ｂ（蒜种－钢板滚动摩擦因数）１．４１４０．１０１０．１０１
１０．１０３０．１０３
００．１１００．１１０
－１０．１１７０．１１７
－１．４１４０．１１９０．１１９
表４ Ｃｅｎｔｒａｌｃｏｍｐｏｓｉｔｅｄｅｓｉｇｎ方案及结果
Ｔａｂｌｅ４ Ｓｃｈｅｍｅａｎｄｒｅｓｕｌｔｓｏｆｃｅｎｔｒａｌｃｏｍｐｏｓｉｔｅｄｅｓｉｇｎ
序号
试验因素
ＡＢ休止角
／（°）
半径
／ｍｍ
蔡星辰
１－１－１２３．３３１５１．９２２１－１２８．３１１３５．９４３－１１２４．０５１４９．３５４１１２８．９１１３４．２６５－１．４１４０２１．１８１６０．２８６１．４１４０３０．４２１３０．１９７０－１．４１４２２．５４１５４．８７
续表４
序号
试验因素
ＡＢ
休止角
／（°）
半径
／ｍｍ８０１．４１４２８．７６１３４．６８
９００２６．３１１４２．８９
１０００２６１４３．８６
１１００２６．８４１３９．２６
１２００２６．５２１４１．２４
１３００２６．７１４０．９６
应用软件ＤｅｓｉｇｎＥｘｐｏｒｔ８．０．６［１２］对试验结果进行分析，得到二次回归模型，对该二次回归模型进行方差分析和回归系数显著性检验，结果如表５、表６所示。由表５可知：发现两试验模型均为极显著，失拟项分别为０．１４３７、０．１６６０均大于０．０５，说明回归方程与实际拟合中非正常误差所占比例小，拟合性较好；两试验的变异系数分别为１．５３、１．４６，故试验可靠性较好；决定系数分别为０．９９０９、０．９７６７，校正决定系数分别为０．９７８３、０．９４４１，均接近于１，故所得回归方程可靠度较高。则有
θ＝２６．４７＋３．２７Ａ＋２．２０Ｂ－０．０３０ＡＢ－
０．２３Ａ２－０．３１Ｂ２－１．８７Ａ２Ｂ－０．８１ＡＢ２
ｒ＝１４１．６１－１０．６４Ａ－７．１４Ｂ＋０．２２ＡＢ＋
１．２６Ａ２＋１．０３Ｂ２＋６．０８Ａ２Ｂ＋２．８７ＡＢ２
（５）
表５ 蒜种堆休止角回归方程方差分析
Ｔａｂｌｅ５ Ｖａｒｉａｎｃｅａｎａｌｙｓｉｓｏｆｇａｒｌｉｃｐｉｌｅｒｅｐｏｓｅａｎｇｌｅｒｅｇｒｅｓｓｉｏｎｅｑｕａｔｉｏｎ方差来源均方和自由度均方Ｆ值Ｐ值
模型８７．５９７１２．５１７８．２１＜０．０００１
同济图书馆Ａ４２．６９１４２．６９２６６．８２＜０．０００１
Ｂ１９．３４１１９．３４１２０．９１０．０００１
ＡＢ０．００４１０．００４０．０２３０．８８６６
Ａ２０．３７１０．３７２．３２０．１８８６
Ｂ２０．６５１０．６５４．０６０．０９９９
Ａ２Ｂ６．９９１６．９９４３．６７０．００１２
ＡＢ２１．３１１．３０８．１４０．０３５７
残差０．８５０．１６
失拟项０．３６１０．３６３．２９０．１４３７
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议