基于EDEM的青稞接触参数仿真标定*

丁文波，朱继平，陈伟，张晋，袁栋，丁艳

(农业农村部南京农业机械化研究所，南京市，210014)

0 引言

青稞作为西藏等地区人民的主粮，在高海拔、寒冷地区的粮食作物种植中，具有不可替代地位。青稞生产机械化程度低，已经影响到青稞产业的发展和当地人民的生活质量，提高青稞生产的机械化水平具有重要的现实意义[1-2]。

青稞机械化播种作为其机械化生产的重要环节，受限于排种器的研究发展，排种器的研究难点主要在于青稞种子在排种器中的运动及受力情况复杂，难以进行深入研究。而EDEM离散元仿真技术已经广泛应用于农业领域[3-4]，能够有效直观的模拟机构的工作原理，有利于高效便捷的对各种机构进行深入的研究[5-8]。进行EDEM离散元仿真研究前需要建立物料的三维模型，设置物料的本征参数和接触参数。建立精确的颗粒三维模型和准确的设置仿真参数，可以提高仿真的精度。获取这些参数有直接测量法[9-10]和虚拟标定[11]两种方法，通常一些难以直接测量的材料接触参数需要通过虚拟标定的方法进行标定。

目前对主要粮食作物如水稻、小麦、大豆、玉米等都已进行相关的研究。张荣芳等[12]利用EDEM研究了不同填充球半径的水稻种子模型对颗粒间的动力学响应特性的影响，寻找种子模型最佳的填充球颗粒数量，结果表明填充颗粒半径为0.21 mm，仿真时长与仿真精度最优；张涛等[13]对大豆排种过程的离散元仿真参数进行研究，标定了大豆种子的本征参数、大豆种子与排种器的接触参数；崔涛等[14]通过离散元分析软件EDEM仿真试验对玉米种子休止角试验所得试验数据进行了验证，为排种器的设计与仿真分析提供理论依据。

目前对于青稞这方面的研究很少，本文以青稞藏青2000为研究对象，根据种子的形状特征采用分段取截面的方法建立其三维模型，并通过台架试验测定其部分参数，为其在仿真试验中的取值提供依据，最后以青稞休止角为指标，通过响应面法标定青稞的主要接触参数，以期为青稞后续的播种研究提供参考。

1 青稞的基本参数与三维模型

1.1 青稞的基本参数

青稞藏青2000由西藏自治区农牧科学院农业研究所于1995年以(藏青320×拉萨白青稞)与(喜拉19×昆仑164)复合杂交选育而成的青稞新品系，穗长方形、四棱长芒，黄白粒，具有良好的增产效果[15-16]。将经过筛选的青稞种子随机抽取五组测得平均千粒重46.5 g，含水率9.53%，密度10.2 g/cm3。随机选取100粒种子使用游标卡尺测得的平均长度7.75 mm、宽度3.76 mm、厚度2.67 mm。

1.2 青稞种子三维模型的建立

种子的三维模型与真实形状越接近，仿真的结果就会越准确，目前对于不规则的种子主要采取激光扫描的方式获取精确的三维模型，如于庆旭等[17]利用逆向工程技术得到三七种子的精确三维轮廓模型；刘彩玲等[18]基于三维激光扫描获取水稻种子精确的三维模型。对于小麦、青稞等近似椭球形，具有对称性的种子主要采用近似法建模[19]。精确的三维建模可提高仿真的精度，但精确建模难度大成本高。青稞藏青2000种子如图1(a)所示，与小麦外形相似，故采用近似法建模。为了提高建模的精度，随机抽取50粒种子，在长度方向进行等分获得等分截面如图1(b)所示，测量其各个截面的宽度k、厚度h的尺寸参数，计算得到各截面宽、厚尺寸平均值如表1所示。再利用CATIA软件中多截面实体功能，连接截面生成种子的一侧实体，由种子的对称性对一半实体进行镜像得到种子的三维模型如图1(c)。

表1 青稞截面尺寸Tab. 1 Section size of highland barley

(a) 青稞种子 (b) 等分截面图

2 台架试验与仿真标定

为了提高离散元仿真试验的可靠性，本研究利用台架试验测量青稞种子与试验材料的静摩擦因数以及青稞种子与试验材料的碰撞恢复系数，为这两种离散元仿真参数的设置提供依据。对于不便直接测量的参数，采用台架试验与仿真试验结合的方式进行仿真标定，以种子的休止角为检验指标进行仿真试验，休止角采用圆筒提升法获得[20]。试验材料选择3D打印材料pla(聚乳酸)塑料[21]，方便为后续3D打印排种器进行排种机理研究做准备。

2.1 台架试验

2.1.1 斜面法确定静摩擦系数取值

通过斜面法测量青稞种子与pla材料的静摩擦系数，试验材料有pla3D打印平板如图2(a)所示(长150 mm、宽50 mm、厚2 mm)、青稞藏青2000种子、万能角度尺、多功能铁架台等。为了防止青稞种子在试验中发生滚动影响试验的结果，将2粒青稞种子粘在一起为一组，共设置五组如图2(b)所示。试验的原理：将pla平板置于多功能铁架台的可转动铁夹上，转动铁夹使得平板处于水平位置，将粘好的种子放置上面，然后缓慢转动铁夹使得平板逐渐倾斜，直到种子发生滑动时停止转动铁夹。用万能角度尺测量此时的平板的倾斜角，每组重复5次，完成平板倾角的测定，计算平均值求得平板倾斜角α=27.04°。根据公式求得青稞与pla平板的最大静摩擦系数x8=0.51。由于试验本身存在误差，故x8取值应该在0.51附近，取值范围定为0.3～0.7具有可信度。

(a) pla 3D打印平板 (b) 青稞种子粘结处理

2.1.2 台架试验确定恢复系数取值

对青稞—pla平板之间的碰撞恢复系数的测量是通过使青稞种子从一定的高度垂直落下，与放置于正下方呈45°放置的pla平板发生碰撞后，作抛物线运动，最后落入收集装置[22]，恢复系数测定装置如图3(a)，通过测量种子初始位置N点与碰撞点O的垂直高度h1、碰撞点O与落点P的垂直高度h2、碰撞点与落点的水平距离s1等位移，代入运动学方程[23]求解得到恢复系数。重复试验10次记录数据，通过计算求得碰撞恢复系数平均值x7=0.49。试验装置原理图如图3(b)。

(a) 恢复系数测定装置

2.1.3 圆筒提升试验确定休止角

采用圆筒提升法测量青稞休止角θ[24]。试验的材料有内径高度分别为21 mm、200 mm的玻璃圆筒、多功能铁架台、直径厚度分别为80 mm和3 mm的pla圆板、万能角度尺、游标卡尺等。将pla圆板悬空置于铁架上，然后将充满青稞颗粒的圆筒口对准圆板的中心倒置，再以0.02 m/s的速度缓慢提升圆筒，最后形成稳定的颗粒堆。通过万能角度尺测量得到颗粒堆休止角θ。重复10次试验求得平均休止角θ=25.1°，试验装置如图4。

图4 圆筒提升试验装置图

2.2 仿真模型的建立及参数设置

2.2.1 青稞的离散元仿真模型

将种子的三维模型(图1(c))导入EDEM软件中，采用多球聚合模型[25]对种子模型进行填充，通过若干个直径不等的球形颗粒重叠堆积完成，得到的种子离散元模型如图5所示，从图中可以看出颗粒模型与种子实际外轮廓较为吻合。

图5 青稞的离散元仿真颗粒模型

2.2.2 仿真参数设置

离散元颗粒堆积仿真试验所需要的基本参数有青稞和pla材料的本征参数及接触参数。其中已知的本征参数如表2所示：青稞的密度由试验测得为1.02 g/cm3；通过查阅文献[26]得到pla材料的泊松比0.29、剪切模量220 MPa、密度1.11 g/cm3。

表2 仿真参数的取值Tab. 2 Value of simulation parameters

未确定的仿真参数通过查找文献[27-28]及台架试验确定其取值的范围如表2所示。通过台架试验测得x8=0.51、x7=0.49。由于试验本身存在一定的误差，故以试验测定的值为中心值扩展，得到这两个仿真参数的取值区间，具有可信度。仿真颗粒的分布按照标准正态分布设置，平均值Mean：1；标准差StdDev：0.05；由青稞种子颗粒的运动特点，接触模型选择“Hertz-Mindlin(noslip)”；仿真中FixedTimeStep采用20%的瑞利时步；仿真中网格尺寸CellSize取3倍最小球形单元尺寸。

2.2.3 圆筒提升仿真试验

在圆筒提升试验中影响试验指标的因素有：青稞的泊松比x1、青稞的剪切模量x2、青稞—青稞的碰撞恢复系数x3、青稞—青稞的静摩擦系数x4、青稞—青稞的滚动摩擦系数x5、青稞—pla平板的滚动摩擦系数x6、青稞—pla平板的碰撞恢复系数x7、青稞—pla平板的静摩擦因数x8。试验的指标：台架试验中青稞的休止角θ。

利用三维软件建立出圆筒提升试验中圆筒和pla圆板的三维模型，将其保存为Stp格式导入EDEM软件中。在圆筒内设置颗粒工厂生成颗粒，经过预仿真试验确定以1 500个/s的速率生成颗粒，生成时间设置为1 s，可产生足够的试验颗粒，颗粒生成后自由落下填充圆筒。在颗粒生成结束后，给圆筒添加一个向上的直线运动，速度为0.02 m/s，随着圆筒的提升，颗粒缓慢从桶内流出，在pla圆板上形成稳定的颗粒堆，多余的颗粒则落入下方的容器中。仿真完成以后，进入EDEM后处理界面，利用软件自带量角功能直接测量颗粒堆的休止角θ1，如图6所示。

图6 青稞颗粒休止角测量

2.3 仿真模型的响应面试验设计

2.3.1 Plackett-Burman试验

通过Plackett-Burman试验筛选出显著性试验因素，并非所有试验因素都对试验指标颗粒堆积角有显著的影响，对于影响不显著的试验因素，无法进行仿真分析。故应用Design-Expert软件进行Plackett-Burman试验，仿真试验因素及水平如表3所示，筛选出x1～x8中对试验指标有显著影响的因素。仿真试验设计及结果如表4所示。

表3 Plackett-Burman试验因素及水平Tab. 3 Factors and levels of Plackett-Burman test

表4 Plackett-Burman 试验设计及结果Tab. 4 Design and results of Plackett-Burman test

通过Design-Expert的分析功能得到表5，从表5可以看出对青稞休止角影响最显著的参数分别是x4、x5、x8是影响休止角大小的主要因素，其他的影响因素对休止角的影响不显著。

表5 Plackett-Burman试验参数显著性分析Tab. 5 Analysis of significance of parameters in Plackett-Burman test

2.3.2 最陡爬坡试验

采用最陡坡试验可以较为快速的确定显著性因素的最佳取值范围，为后续的响应面试验提供取值依据，从而确定参数的最优值。由前文中Plackett-Burman试验可知显著性影响因素有x4、x5、x8，取值范围分别为0.1～0.6、0～0.1、0.3～0.7，等分取六个梯度值较为合适。其他因素对试验影响不显著，则都取中间值，分别取x1(0.125)、x2(5.5 MPa)、x3(0.5)、x6(0.05)、x7(0.5)。

最陡爬坡试验设计及结果如表6所示，结果表明随着x4、x5、x8的增大仿真试验的颗粒休止角也不断增大，仿真试验的颗粒休止角与台架试验测得的休止角之间的相对误差先减小后增大。在第二梯度时相对误差最小，故最优值在2号附近，则最优的取值区间在1号和3号之间。故表7中1、2、3号取值可分别作为后面Box-Behnken试验中-1、0、1水平的取值。

表6 最陡爬坡试验及结果Tab. 6 Design and results of steepest ascent test

2.3.3 Box-Behnken试验及回归模型

根据Plackett-Burman试验与最陡爬坡试验，确定了显著性影响因素有x4、x5、x8及其在Box-Behnken试验中-1、0、1水平的取值。其他因素的取值与最陡爬坡试验中相同。采用Design-Expert软件进行Box-Behnken试验设计并建立回归模型，共进行17次仿真试验。Box-Behnken试验设计及结果如表7所示。根据模型得到其二次多项式方程如式(1)。

θ1=26.67+5.12x4+1.51x5-0.29x8-

0.63x4x5+0.15x4x8+0.15x5x8+

0.61x42-0.81x52+0.20x82

(1)

表7 Box-Behnken试验设计及结果Tab. 7 Design and results of Box-Behnken test

仿真试验的方差分析如表8所示，模型的p值小于0.01该回归模型的拟合性好；x4、x5等因素的p值均小于0.01，说明这两个因素都对颗粒休止角具有极显著性的影响；x8的p值小于0.05说明对休止角具有显著性影响；交互项x4x5、x42、x52的p值都小于0.01说明其对休止角也具有极显著性影响，交互项x4x5对休止角的影响如图7所示，x4、x5、都对休止角起正效应作用；交互项x4x8、x5x8、x82的p值都大于0.05，对颗粒的休止角没有显著影响；决定系数R2=0.998 2，校正决定系数R2adj=0.996 0，二者都接近1，表明拟合方程具有较高的可靠性。

在保证模型显著、失拟项不显著的基础上，剔除模型中的不显著回归项，对模型进行优化[29]得到新的回归方程如式(2)。

θ1=26.76+5.12x4+1.51x5-0.29x8-

0.63x4x5+0.62x42-0.80x52

(2)

最后，将台架试验测定的青稞休止角代入到Design-Expert中，求得3个显著性参数的最优值，x4(0.19)、x5(0.01)、x8(0.43)。

表8 Box-Behnken 试验设计二次多项式模型方差分析Tab. 8 ANOVA of quadratic polynomial model of Box-Behnken test

图7 x4与x5交互作用

2.4 试验验证

将Box-Behnken试验得到的三个显著性参数的最优值代入到EDEM仿真模型中，进行仿真验证试验。仿真参数设置：青稞—青稞的静摩擦系数x4=0.19、青稞—青稞的滚动摩擦系数x5=0.01、青稞—pla平板的静摩擦因数x8=0.43，其他非显著性因素与之前相同，重复5次重复试验得到颗粒的平均休止角θ1=25.26°，与实际青稞休止角25.1°相对误差0.64%，相对误差很小，说明标定的参数可信，可以用于青稞离散元仿真，为后续青稞播种等研究提供参考。

3 结论

1) 通过Plackett-Burman试验，从8个试验因素中找出了主要影响试验指标青稞休止角的三个显著性因素：青稞—青稞的静摩擦系数x4、青稞—青稞的滚动摩擦系数x5、青稞—pla 3D打印材料的静摩擦因数x8。其他因素为非显著性因素，对试验指标的影响不显著。

2) 根据最陡爬坡试验确定3个显著性参数的最优值取值范围，为后面进行Box-Behnken试验提供高低水平的取值，以此为基础建立显著性参数与休止角之间的二次回归模型并对其进行优化，试验表明：3个显著性因素对颗粒休止角具有显著的影响，此外青稞与青稞的静摩擦系数的二次项、青稞与青稞的滚动摩擦系数的二次项、二者的交互项对颗粒休止角也有显著影响。

3) 取响应面试验得到青稞—青稞的静摩擦系数x4为0.19、青稞—青稞的滚动摩擦系数x5为0.01、青稞—pla平板的静摩擦因数x8为0.43。其余非显著性参数取值：青稞的泊松比x1为0.18、青稞的剪切模量x2为5.5 MPa、青稞—青稞的碰撞恢复系数x3为0.5、青稞—pla平板的滚动摩擦系数x6为0.05、青稞—pla平板的撞恢复系数x7为0.5。将以上参数导入EDEM中进行仿真验证试验与实际台架试验的结果进行比较，平均休止角θ1=25.26°与实际青稞休止角25.1°相对误差0.64%，误差很小，表明本研究标定的青稞离散元仿真参数具有较高的可靠性，可以为后续的青稞仿真研究提供参考。