2018, Vol. 39
2. 安阳工学院 机械工程学院，河南 安阳 455000
2. School of Mechanical Engineering, Anyang Institute of Technology, Anyang 455000, China
农作物籽粒种群属于典型的由大量离散单元颗粒组成的离散系统，采用有限元或流体等数学分析手段都很难对其进行描述，而采用离散元法可以直接获得离散物质大量复杂行为信息及不易测量的颗粒尺度行为信息[1-2]。近年来离散元法逐渐应用在农业生产过程中农业物料的运动状态模拟中[3-8]。
采用离散元法对颗粒进行建模和相关应用分析时，需要确定合适的材料和接触参数，才能确保离散元模拟结果预测的准确性[9-10]。有些学者通过直接测量得到颗粒的参数，Coetzee等[11]采用实验室剪切试验和压缩试验确定了材料的内摩擦角和刚度；Wang等[12]设计了一种测定不规则玉米碰撞恢复系数的自动控制平台，研究不同形状玉米对碰撞恢复系统的影响。然而直接测量单个颗粒的参数值并不总是可行和具有代表性，尤其是农业物料其参数更为复杂[13]。因此有些学者倾向于使用参数标定法，其涉及到改变参数值，直到模拟结果与相应的目标响应值之间取得最佳匹配[14]。
张荔等[15]以真实再现磨料颗粒堆积角为研究目的，在研究滚磨加工中磨料颗粒的离散元模拟参数时提出一种在已有正交试验基础上考虑参数间交互作用的参数标定方法。冯俊小等[16]通过虚拟堆积实验和滚筒内颗粒混合仿真与实体试验现象对比，标定离散元法的颗粒物性参数。张锐等[17]使用EDEM软件对颗粒之间相互作用参数进行标定，提出一种系统标定沙土颗粒相互作用参数的方法。武涛等[18]借助Box Behnken试验方法进行堆积角仿真试验，实现黏性土壤离散元模型的接触参数与接触模型参数标定。Coetzee[19]和Pasha等[20]通过研究表明不同形状颗粒包括不同球组合的颗粒在离散元模拟试验时采用的参数也不同，需要通过试验标定。由于玉米籽粒的离散元模型由多球填充组成，与真实玉米籽粒的光滑表面有一定差别[21]，为了使仿真试验更接近实际，确保离散元模拟玉米籽粒运动的准确性和有效性，因此需要对玉米籽粒离散元参数进行标定。
本文以玉米籽粒为研究对象，以堆积角为响应值，提出一种基于响应面法的离散元参数标定方法。将堆积角离散元模拟试验与实体试验对比，对玉米籽粒的离散元参数及它们之间的交互作用进行分析，标定其具体数值，为播种、脱粒等生产过程中玉米籽粒的运动状态研究奠定基础。
1 材料与方法 1.1 标定方法堆积角是表征散粒体物料颗粒流动、摩擦等特性的宏观参数[22-23]，操作方便、简单，农作物物料颗粒离散元模型多采用该方法作为验证试验[24-25]。
本文以堆积角作为表征玉米籽粒离散元参数的目标响应值，通过EDEM软件进行玉米籽粒堆积角仿真试验，不断调整玉米籽粒离散元参数使仿真结果与实体试验结果一致，实现对玉米籽粒离散元参数的标定。
1.2 实体试验试验采用玉米品种‘桥玉八号’，籽粒含水量（w）11.02%，密度1 197 kg·m–3。采用无底圆筒法测量玉米籽粒堆积角，选用直径60 mm、高180 mm的无底圆筒，侧面和底面均采用不锈钢材料。试验时将玉米籽粒通过漏斗落入无底圆筒内，填充结束后采用电机带动滑轮组以0.05 m·s–1的速度缓慢提升圆筒，最终玉米籽粒稳定地堆积在不锈钢平板中间。采用常规方法计算玉米籽粒的堆积角，测量玉米籽粒堆积物不同方位的直径及高度，多次测量取平均值。试验重复5次，堆积角在25.55°与27.46°之间，平均值为26.43°，标准误为0.36。
1.3 仿真试验 1.3.1 玉米籽粒模型将实体试验采用的玉米籽粒选取150粒，统计长宽厚等数据，得到玉米籽粒平均长度12.84 mm，标准误0.049，平均宽度10.24 mm，标准误0.048，平均厚度4.64 mm，标准误0.033。颗粒建模时将玉米籽粒外形简化处理为近似棱台形，其大端、小端和厚度等分别取真实玉米籽粒三维尺寸平均值，在Solidworks软件中完成玉米籽粒三维建模，生成.igs格式文件并导入到EDEM软件中作为手动填充模型的模板。
手动填充玉米籽粒模型(图1)由41个相互重叠的球颗粒组成，需要指定每个球颗粒的半径及球心X、Y、Z坐标以确定球颗粒间相互位置关系。由于手动填充的玉米籽粒模型较平滑，填充时球颗粒水平排列放置，由大端到小端球颗粒数目逐渐减少，果柄部位采用一个较小半径的球颗粒代替，其余部位球颗粒半径采用真实玉米籽粒厚度均值。
为减少仿真时间，对玉米籽粒离散元模型进行了简化处理，因此在仿真时需标定该模型的离散元仿真参数，以使玉米籽粒在宏观运动上与真实情况保持一致。
|
图 1 玉米籽粒仿真模型 Figure 1 Maize kernel simulation model |
采用EDEM软件进行仿真试验时，需要设置颗粒及边界相关参数。边界为不锈钢材料，质地均匀参数唯一确定，而玉米籽粒由于生物特性的复杂多变性，离散元参数数值并不唯一。除玉米籽粒密度为试验测量值外，其他参数通过查阅国内外文献[26-30]得到，归纳总结其变化范围如表1所示。
|
表 1 仿真参数设置 Table 1 The setting of simulation parameters |
在EDEM软件中导入与实体试验尺寸相同的无底圆筒模型，在圆筒底部添加与圆筒垂直的水平平面。为使圆筒内刚好填满颗粒且不影响堆积角测量，在圆筒上部设置一个倒置的锥形漏斗，颗粒工厂设置在圆筒上方，动态随机生成玉米籽粒落入圆筒内。当玉米籽粒填满圆筒时，用圆筒上方的圆柱体将多余的玉米籽粒去除，如图2所示，圆筒外部漏斗可避免刮下来的颗粒落在平板上影响堆积角的测量。然后设置圆筒以0.05 m·s–1的速度缓慢向垂直于不锈钢平面的上方移动，玉米籽粒逐渐流出，最终形成一个稳定的锥形物料堆。仿真采用的Rayleigh时间步长为25%，从玉米籽粒生成到圆筒上升达到玉米籽粒堆积稳定状态，仿真总时长为6 s。
|
图 2 试验仿真模型 Figure 2 Experiment simulation model |
采用PB试验可以从多个仿真参数中筛选出有显著效应的参数。根据EDEM仿真时用到的与玉米籽粒相关的参数，选取了8个参数变量，预留3个虚拟变量作误差分析，均匀分布在8个参数变量中。本文采用了N=12的PB设计表，以玉米籽粒堆积角为响应值，每个参数变量选取高、低2个水平(+1和–1)，如表2所示。
|
|
表 2 Plackett-Burman试验因素与水平1) Table 2 Plackett-Burman experimental factors and levels |
通过PB试验可以得到若干个对响应值影响较为显著的因素，为建立有效的响应面拟合方程，需通过最陡爬坡试验来快速、经济地逼近各显著因素的最佳区域。根据PB试验结果中显著因素的正负效应确定最陡上升路径方向，以因素效应值确定变化步长，其余不显著因素在PB试验中取平均水平。
1.4.3 响应面试验由最陡爬坡试验逼近最佳响应区域后，采用中心组合设计试验进行响应面分析，选取显著因素的5个水平进行试验设计，不显著因素取平均水平，设计中采用6个中心点。通过试验结果拟合响应面模型，得到模型的多元函数，再依据回归方程绘制响应面的三维立体分析图，并预测最接近实际堆积角的玉米籽粒仿真参数最佳组合。
1.4.4 优化验证试验为考察优化参数组合后得到的堆积角与实际堆积角之间的误差，采用优化后玉米籽粒仿真参数参照“1.3”进行堆积角仿真试验，并将仿真试验结果与实体试验结果对比，验证优化参数后仿真试验结果是否更接近实际情况。
2 结果与分析 2.1 PB试验结果PB试验设计及结果如表3所示，应用Design Expert软件对试验结果进行方差分析，由于PB试验不能区分主效应与交互效应，为避免遗漏重要因素，选用t检验法(P<0.1)选取影响显著的因素。
方差分析结果如表4所示，各因素对堆积角的影响从大到小依次为玉米–玉米静摩擦系数(P=0.009 1)、玉米泊松比(P=0.023 1)、玉米–玉米滚动摩擦系数(P=0.098 3)，其他因素对堆积角影响不显著(P>0.1)。因此在最陡爬坡试验时选取这3个显著因素。
|
|
表 3 Plackett-Burman试验设计及结果1) Table 3 Design and result of Plackett-Burman experiment |
|
|
表 4 Plackett-Burman试验方差分析1) Table 4 Variance analysis of Plackett-Burman experiment |
同时通过帕累托图(Pareto chart)(图3)不仅可以直观看到各因素的显著性排序，还可观察到因素的正负效应。玉米–玉米静摩擦系数和玉米–玉米滚动摩擦系数对堆积角具有正效应，而玉米泊松比对堆积角有负效应，因此在最陡爬坡试验时，玉米–玉米静摩擦系数和玉米–玉米滚动摩擦系数为上升路径，玉米泊松比为下降路径。
|
图 3 Plackett-Burman试验帕累托图 Figure 3 Pareto chart of Plackett-Burman experiment |
根据PB试验结果筛选出的3个显著因素及其爬坡方向，最陡爬坡试验设计及结果如表5所示。随着玉米泊松比的降低，玉米–玉米静摩擦系数和玉米–玉米滚动摩擦系数的增加，仿真得到的玉米籽粒堆积角与实体试验堆积角的误差先减少后增加，处理2误差最小，说明各因素最佳区域在处理2附近。因此选取玉米泊松比0.44，玉米–玉米静摩擦系数0.22、玉米–玉米滚动摩擦系数0.05为中心点进行后续的响应面试验。
|
|
表 5 最陡爬坡试验设计及结果1) Table 5 Design and result of the steepest ascent experiment |
根据PB试验得到的3个显著因素及最陡爬坡试验得到的中心点进行3因素5水平的中心组合设计试验，试验设计及结果如表6所示。利用Design Expert软件对表6试验数据进行二次多项式回归拟合，构建玉米籽粒仿真参数与堆积角之间的模型，得到的模型方程为：
|
|
表 7 中心组合设计试验模型方差分析1) Table 7 Variance analysis of model of central composite designs experiment |
Y=108.02–593.26A+585.63E–481.15G– 878.52AE–193.36AG–958.01EG+ 866.58A2–295.45E2+8 868.81G2。
该模型的方差分析结果如表7所示，该拟合模型P=0.006 1<0.01，说明在该模型中自变量与因变量相关性极显著。该模型方程拟合度较好，玉米泊松比、玉米–玉米静摩擦系数、玉米–玉米滚动摩擦系数及G2在95%置信区间内对玉米籽粒堆积角影响显著，拟失项P=0.389 0>0.05，差异不显著。精密度为9.425 7，相关系数R2=0.835 7，校正决定系数R2Adj=0.687 9，表明模型可信度高，预测值与实测值具有相关性。
|
|
表 8 中心组合设计优化模型方差分析1) Table 8 Variance analysis of modified model of central composite designs |
为进一步优化模型方程，在保证模型显著性好，相关系数高的前提下，剔除不显著因素，得到新的方程为：
Y=122.81–602.93A+537.73E–776.99G– 878.52AE+866.58A2–295.45E2+8 868.81G2。
对优化后模型重新进行方差分析(表8)，优化后拟合模型P=0.001 1，显著性提高，精密度=11.547 8，相关系数R2=0.820 5，校正决定系数R2Adj=0.715 8，较优化前增大，使模型拟合度提高。
|
|
表 6 中心组合设计试验设计及结果 Table 6 Design and result of central composite designs experiment |
优化后发现玉米泊松比与玉米–玉米静摩擦系数的交互作用显著(P=0.040 0<0.05)，应用Design Expert软件绘制了玉米–玉米滚动摩擦系数=0.050时，玉米泊松比与玉米–玉米静摩擦系数的双因子效应分析图(图4)。与玉米泊松比方向比较，玉米–玉米静摩擦系数效应面曲线(图4a)较陡，其等高线密度高于沿玉米泊松比方向移动的密度，表明玉米–玉米静摩擦系数对堆积角的影响较玉米泊松比显著，这与方差分析结果一致。等高线(图4b)形状类似马鞍形，说明玉米泊松比和玉米–玉米静摩擦系数之间交互作用显著。在玉米泊松比取确定值时，随着玉米–玉米静摩擦系数的增大，玉米籽粒堆积角呈增大趋势；在玉米–玉米静摩擦系数取确定值时，随着玉米泊松比的减小，玉米籽粒堆积角呈增大趋势；在玉米–玉米静摩擦系数值较高，玉米泊松比减小对玉米籽粒堆积角增大趋势的影响越发明显。
|
图 4 玉米泊松比和玉米–玉米静摩擦系数对堆积角影响的效应面及等高线图 Figure 4 Response surface and contour plot of the effects of corn Poisson's ratio and corn-corn static friction coefficient on repose angle |
由于模型的理想响应值应接近实际玉米堆积角26.43°，因此选择Point prediction优化模块预测分析结果，获得1组接近目标值的最优参数组合：玉米泊松比0.438、玉米–玉米静摩擦系数0.182、玉米–玉米滚动摩擦系数0.051，其余参数为中间水平，玉米弹性模量273.9 MPa、玉米–玉米恢复系数0.233、玉米–不锈钢恢复系数0.672、玉米–不锈钢静摩擦系数0.334、玉米–不锈钢滚动摩擦系数0.062 4。
2.4 优化验证试验结果为验证响应面法得到的最优参数的准确性和有效性，采用标定后的玉米籽粒离散元参数进行EDEM仿真试验，试验重复3次，得到的玉米籽粒堆积角分别为27.18°、25.87°和27.63°，平均值为26.89°，仿真值与实测值吻合良好。
3 结论1)通过响应面优化试验设计方法筛选出了3个对玉米籽粒堆积角有显著影响的因素：玉米泊松比、玉米–玉米静摩擦系数、玉米–玉米滚动摩擦系数，且玉米泊松比与玉米–玉米静摩擦系数有交互作用。
2)通过标定得到的玉米籽粒离散元参数为：玉米泊松比0.438、玉米–玉米静摩擦系数0.182、玉米–玉米滚动摩擦系数0.051、玉米弹性模量273.9 MPa、玉米–玉米恢复系数0.233、玉米–不锈钢恢复系数0.672、玉米–不锈钢静摩擦系数0.334、玉米–不锈钢滚动摩擦系数0.062 4。验证试验表明基于响应面法标定玉米籽粒离散元参数是可行的，可提高离散元仿真试验的准确性。
| [1] |
王国强, 郝万军, 王继新. 离散单元法及其在EDEM上的实践[M]. 西安: 西北工业大学出版社, 2010.
(  0) |
| [2] |
胡国明. 颗粒系统的离散元素法分析仿真: 离散元素法的工业应用与EDEM软件简介[M]. 武汉: 武汉理工大学出版社, 2010.
( 0) |
| [3] |
MA Z, LI Y, XU L, et al. Dispersion and migration of agricultural particles in a variable-amplitude screen box based on the discrete element method[J]. Comput Electron Agr, 2017, 142, Part A: 173-180. ( 0) |
| [4] |
张涛, 刘飞, 赵满全, 等. 基于离散元的排种器排种室内玉米种群运动规律[J]. 农业工程学报, 2016, 32(22): 27-35. DOI:10.11975/j.issn.1002-6819.2016.22.004 ( 0) |
| [5] |
BOAC J M, AMBROSE R P K, CASADA M E, et al. Applications of discrete element method in modeling of grain postharvest operations[J]. Food Eng Rev, 2014, 6(4): 128-149. DOI:10.1007/s12393-014-9090-y ( 0) |
| [6] |
LENAERTS B, AERTSEN T, TIJSKENS E, et al. Simulation of grain-straw separation by Discrete Element Modeling with bendable straw particles[J]. Comput Electron Agr, 2014, 101: 24-33. DOI:10.1016/j.compag.2013.12.002 ( 0) |
| [7] |
纪玉杰, 薛宏达, 李成华. 玉米精密排种器工作过程的离散元仿真分析[J]. 计算机仿真, 2013, 30(3): 393-397. ( 0) |
| [8] |
SARNAVI H J, MOHAMMADI A N, MOTLAGH A M, et al. DEM model of wheat grains in storage considering the effect of moisture content in direct shear test[J]. Res J Appl Sci Eng Technol, 2013, 5(3): 829-841. ( 0) |
| [9] |
RACKL M, HANLEY K J. A methodical calibration procedure for discrete element models[J]. Powder Technol, 2017, 307: 73-83. DOI:10.1016/j.powtec.2016.11.048 ( 0) |
| [10] |
COETZEE C J. Review: Calibration of the discrete element method[J]. Powder Technol, 2017, 310: 104-142. DOI:10.1016/j.powtec.2017.01.015 ( 0) |
| [11] |
COETZEE C J, ELS D N J. Calibration of granular material parameters for DEM modelling and numerical verification by blade–granular material interaction[J]. J Terramechanics, 2009, 46(1): 15-26. DOI:10.1016/j.jterra.2008.12.004 ( 0) |
| [12] |
WANG L, ZHOU W, DING Z, et al. Experimental determination of parameter effects on the coefficient of restitution of differently shaped maize in three-dimensions[J]. Powder Technol, 2015, 284: 187-194. DOI:10.1016/j.powtec.2015.06.042 ( 0) |
| [13] |
GRIMA A P, WYPYCH P W. Investigation into calibration of discrete element model parameters for scale-up and validation of particle–structure interactions under impact conditions[J]. Powder Technol, 2011, 212(1): 198-209. DOI:10.1016/j.powtec.2011.05.017 ( 0) |
| [14] |
COETZEE C J, ELS D N J. Calibration of discrete element parameters and the modelling of silo discharge and bucket filling[J]. Comput Electron Agr, 2009, 65(2): 198-212. DOI:10.1016/j.compag.2008.10.002 ( 0) |
| [15] |
张荔, 李文辉, 杨胜强. 滚磨光整加工中磨料颗粒堆积角的离散元参数标定[J]. 中国科技论文, 2016, 11(16): 1821-1825. DOI:10.3969/j.issn.2095-2783.2016.16.006 ( 0) |
| [16] |
冯俊小, 林佳, 李十中, 等. 秸秆固态发酵回转筒内颗粒混合状态离散元参数标定[J]. 农业机械学报, 2015, 46(3): 208-213. DOI:10.6041/j.issn.1000-1298.2015.03.030 ( 0) |
| [17] |
张锐, 韩佃雷, 吉巧丽, 等. 离散元模拟中沙土参数标定方法研究[J]. 农业机械学报, 2017, 48(3): 49-56. DOI:10.6041/j.issn.1000-1298.2017.03.006 ( 0) |
| [18] |
武涛, 黄伟凤, 陈学深, 等. 考虑颗粒间黏结力的黏性土壤离散元模型参数标定[J]. 华南农业大学学报, 2017, 38(3): 93-98. DOI:10.7671/j.issn.1001-411X.2017.03.015 ( 0) |
| [19] |
COETZEE C J. Calibration of the discrete element method and the effect of particle shape[J]. Powder Technol, 2016, 297: 50-70. DOI:10.1016/j.powtec.2016.04.003 ( 0) |
| [20] |
PASHA M, HARE C, GHADIRI M, et al. Effect of particle shape on flow in discrete element method simulation of a rotary batch seed coater[J]. Powder Technol, 2016, 296: 29-36. DOI:10.1016/j.powtec.2015.10.055 ( 0) |
| [21] |
MARKAUSKAS D, RAMÍREZ-GÓMEZ Á, KAČIANAUSKAS R, et al. Maize grain shape approaches for DEM modelling[J]. Comput Electron Agr, 2015, 118: 247-258. DOI:10.1016/j.compag.2015.09.004 ( 0) |
| [22] |
GUO Z, CHEN X, LIU H, et al. Theoretical and experimental investigation on angle of repose of biomass–coal blends[J]. Fuel, 2014, 116(15): 131-139. ( 0) |
| [23] |
吴爱祥, 孙业志, 刘湘平. 散体动力学理论及其应用[M]. 北京: 冶金工业出版社, 2002.
( 0) |
| [24] |
贾富国, 韩燕龙, 刘扬, 等. 稻谷颗粒物料堆积角模拟预测方法[J]. 农业工程学报, 2014(11): 254-260. DOI:10.3969/j.issn.1002-6819.2014.11.031 ( 0) |
| [25] |
GRIMA A P, WYPYCH P W. Development and validation of calibration methods for discrete element modelling[J]. Granul Matter, 2011, 13(2): 127-132. DOI:10.1007/s10035-010-0197-4 ( 0) |
| [26] |
周文秀. 玉米籽粒的物理力学特性研究[D]. 哈尔滨: 东北农业大学, 2015.
( 0) |
| [27] |
崔涛, 刘佳, 杨丽, 等. 基于高速摄像的玉米种子滚动摩擦特性试验与仿真[J]. 农业工程学报, 2013, 29(15): 34-41. DOI:10.3969/j.issn.1002-6819.2013.15.005 ( 0) |
| [28] |
DAI F, ZHAO W Y, HAN Z S, et al. Experiment on Poisson's ratio determination about corn kernel[C]// Trans Tech Publications Ltd. Advanced Materials Research. Switzerland: Trans Tech Publications Ltd., 2013: 799-802.
( 0) |
| [29] |
杨明芳. 基于离散元法的玉米排种器的数字化设计方法研究[D]. 长春: 吉林大学, 2009.
( 0) |
| [30] |
COŞKUN M B, YALҪIN N B, ÖZARSLAN C. Physical properties of sweet corn seed (Zea mays saccharata Sturt.)
[J]. J Food Eng, 2006, 74(4): 523-528. DOI:10.1016/j.jfoodeng.2005.03.039 ( 0) |