南夏豆25收获期植株离散元建模方法与参数标定

谭云峰，陈 霖，胡 森，王 键，陈治帆，吕小荣

（四川农业大学机电学院，四川 雅安 625014）

大豆作为我国重要农业作物之一，是植物油和植物蛋白主要来源，对我国油脂安全供给、畜牧业发展及国家粮食安全具有重要意义[1-2]。大豆机械化收获是实现大豆生产全程机械化重要环节。作为收获机核心部件，脱粒装置对整机工作质量起决定性作用[3]，收获过程损失降低经济效益。茎秆是大豆收获过程中主要杂质，其物理特性对收获过程影响较大。脱粒过程复杂，试验中较难获得物料状态变化数据。因此，有必要建立还原度较高的大豆茎秆模型，为脱粒过程中茎秆与脱粒装置交互作用研究提供可靠依据。

目前已有针对植株仿真试验的相关研究。大多采用圆形颗粒堆砌成刚性长杆模型，但植株脱粒过程中出现弯曲变形、破碎等现象，传统刚性植株模型在仿真计算过程中被视作无任何位移变形的刚体。刘基、谢干等使用茎秆几何参数建立刚性大豆、小麦茎秆模型，利用离散元仿真研究脱粒过程中籽粒破碎情况、功耗影响因素[4-5]。王万章等在多球面填充法构建刚性茎秆模型基础上，基于Hertz－Mindlin with bonding接触模型构建小麦柔性茎秆模型，该接触模型可承受一定量法向和切向位移，直至到达临界法向和切向力，模型可弯曲和破碎[6]；赵吉坤等利用EDEM 软件中Hertz-Mindlin with Bonding 模型建立水稻秸秆模型，在秸秆弯曲力学特性试验基础上，对模型参数进行标定[7]；Wang 等利用空心圆柱弹性键建立水稻植株离散元模型并在颗粒水平上准确描述柔性水稻植株在多重相互作用条件下动态行为[8]；Mao 等使用空心圆柱键建立柔性管模型并验证可有效再现模型流动行为，用于建立新秸秆模型[9]；Leblicq 等建立DEM中茎秆弯曲模型，可用于单个茎秆弯曲行为验证[10]。Yu 等基于离散元方法，建立玉米穗数值模型并用于玉米脱粒过程仿真，仿真结果与试验数据吻合[11-12]。张锋伟等早期研究结果表明，分析玉米秸秆破碎揉丝过程中受力，可表征玉米秸秆力学特性，结合力学特性试验对玉米秸秆进行颗粒黏结模型黏结参数校核[13]。刘禹辰等通过对玉米秸秆进行压缩、玉米秸秆外表皮与内穰拉伸和剪切等力学试验，采用离散元法进行力学特性仿真，将黏结键断裂情况作为玉米秸秆受外力时破裂程度评定指标[14]。张开飞等针对大豆秸秆进行弯曲、剪切力学试验发现，载荷加载速度对大豆秸秆弯曲应力及剪切应力产生影响，同一株大豆秸秆不同部位其表现强度也有较大差异[15]。Öztürk等研究大豆茎秆切削性能，切削力随茎秆直径增加而线性增加，不同高度间茎秆切削力无显著差异[16]。基于离散元法建立的茎秆模型已广泛应用于现代农业装备研究，但大豆茎秆柔性离散元模型建立与黏结参数标定研究仍未见报道。

针对实际脱粒过程中大豆茎秆出现弯曲、破碎等情况，本文以南夏豆25（其在四川省平坝、丘陵及低山区广泛种植）为研究对象，通过弯曲物理试验测得大豆茎秆参数极限断裂载荷和极限断裂位移，利用EDEM 软件中Hertz-Mindlin with bond⁃ing 接触模型，构建可破碎柔性大豆植株离散元模型，以极限断裂载荷和极限断裂位移为响应值，将其弯曲仿真试验与物理试验对比，分析大豆Bonding 参数及其交互作用，并对其Bonding 参数进行标定。以籽粒损失率为指标，进行仿真与物理对比试验，验证模型准确性。

1 材料与方法

1.1 大豆植株基本物性参数确定

本研究大豆植株取样于四川省仁寿县玉米-大豆间套作种植基地，品种为南夏豆25，为获取大豆植株几何尺寸，利用游标卡尺（精度：0.01 mm）进行测量，样品数量为100株，结果如表1所示。

表1 大豆植株几何尺寸Table 1 Plant geometry of soybean

1.2 大豆茎秆物理弯曲试验

收获期大豆茎秆含水率为40%~50%，但大豆茎秆含水率对大豆茎秆黏聚力、流动函数值并无显著影响，仅对其自身力学物理性质有影响[17]。为测得大豆茎秆极限破碎位移和极限破碎载荷，选用直径为7.13~7.75 mm茎秆为试验材料，在三点弯曲夹具质构仪（精度为0.001 mm、0.001 N、刀刃直径为10 mm，购自上海保圣实业发展有限公司）上作弯曲试验，试验过程见图1。试验重复20次取平均值，得到大豆茎秆极限位移和载荷，位移载荷曲线见图2。由图2可见，随刀刃下压，茎秆开始弯曲发生弹性变形，大豆茎秆所承受载荷迅速增大，当达到最大载荷98.994 N 时，大豆茎秆加载位移为3.559 mm，茎秆与仪器接触点四周产生纵向裂横，伴随断面增大的同时，载荷力减小，直到茎秆下表面出现横向断裂，载荷下降至零。

图1 大豆植株弯曲过程Fig.1 Bending process of soybean plant

图2 大豆植株弯曲力-位移曲线变化Fig.2 Change of curve of bending force-displacement of soybean plant

2 收获期大豆植株离散元仿真模型创建与仿真试验

2.1 离散单元法基本原理及接触模型选取

离散单元法最初由美国学者Cundall P.A 在20世纪70 年代初提出。离散单元法是将不连续体分离为刚性单元体，各单元体满足运动方程，运用迭代方法进行循环迭代计算刚性单元体的运动方程，得到研究对象宏观运动规律[18-21]。

Hertz-Mindlin with bonding V2接触模型采用球体之间镶嵌重叠将颗粒组合体连接，球体单元以Bond 键相互粘结方式构造，黏结点可承受来自外界的切向和法向位移，黏结断裂需满足达到最大法向和切向剪切应力。Hertz-Mindlin with bonding V2 接触模型适用于颗粒破碎、分解的问题。根据大豆茎秆可弯曲、破碎特性，采用该接触模型研究大豆植株离散元参数，其模型原理图和Bond 键示意分别如图3和4所示。

图3 Bonding 模型原理图Fig.3 Schematic diagram of the Bonding model

图4 Bond键示意图Fig.4 Bond bond diagram

2.2 大豆植株黏结离散元模型建立

通过前期预试验，得出收获期不同含水率豆荚破荚力分布在0.76~6.64 N，脱粒作业中脱粒齿提供打击力为40~70 N，远高于豆荚破荚力，实际脱粒作业中豆荚破开籽粒离开豆荚即视为完成脱粒作业，为简化模型，可将豆荚简化为籽粒。考虑仿真模型计算体量，将大豆植株视为圆柱体，基于Hertz-Mindlin with bonding V2 黏结接触模型，利用半径3.5 mm 球形粒子相互黏结构建大豆茎秆离散元Bonding V2模型。大豆籽粒形状与椭球体相似，因此在建立颗粒模型时将大豆种子简化为椭球体，测得大豆籽粒尺寸数据后，利用五球填充方法构建大豆籽粒仿真模型，可反映大豆籽粒接触状态[22-23]。为分析脱粒装置脱粒分离能力和夹带损失率，籽粒与茎秆间同样可利用Bond 键黏结。籽粒与茎秆质量比例等于喂入物料草谷比3∶2，收获时期大豆植株离散元模型如图5所示。

图5 收获期大豆植株离散元模型Fig.5 Discrete element model of soybean plant at harvest

2.3 仿真试验参数设置

利用离散元仿真软件EDEM 2020 建立与物理试验一致的大豆茎秆弯曲数值模拟试验，如图6所示。仿真试验数据保存间隔为0.001 s，总时长预设为10 s，时间步长为8.7846×10-7s，网格最小颗粒半径划分为3倍。通过预试验及文献得到仿真中用到的材料及接触参数[24-25]，数值如表2所示。

图6 大豆茎秆弯曲数值模拟试验Fig.6 Numerical simulation experiment of soybean stem bending

表2 离散元仿真参数Table 2 Discrete element simulation parameters

2.4 试验方?法设计

Bonding 模型由单位面积法向刚度、切向刚度、临界法向应力、临界切向应力和黏结半径比例5个参数决定最终黏结键强度。为探究各参数对极限位移及载荷影响的显著性，以极限断裂位移Y1和极限断裂载荷Y2为目标，以单位面积法向刚度、单位面积切向刚度、临界法向应力、临界切向应力和黏结半径比例为标定参数。为获取对试验结果有显著效应的参数，采用Plackett-Burman试验方法，观察显著因素t值正负效应，确定最陡爬坡试验路径方向，对仿真过程中其余不显著因素取中间值；基于上述Plackett-burman 试验结果，进行最陡爬坡试验，确定显著仿真因素最佳范围；最后采用Box-Behnken试验进行响应面优化分析，获取试验仿真最佳参数。

3 仿真试验结果与分析

3.1 Plackett-burman试验

Plackett-burman 试验5 个因子为单位面积法向刚度、单位面积切向刚度、临界法向应力、临界切向应力和黏结半径（X1、X2、X3、X4、X5），选用N=11 的Plackett-burman 试验设计，预留6 个虚拟项作误差分析，试验因素水平如表3所示。

表3 Plackett-burman试验因素水平Table 3 Plackett-burman trial factor levels

Plackett-burman 试验结果如表4所示，为得到各因素显著效应情况，利用Design-expert 11.0 进行方差分析和t检验选取影响显著的因素，结果分别如表5、图7所示。

图7 帕累托图Fig.7 Pareto chart

表4 Plackett-burman试验结果Table 4 Plackett-burman trial results

表5 Plackett-burman试验结果分析Table 5 Plackett-burman test results analysis

通过表5可知，对极限断裂位移Y1的影响因素从大到小依次为X3、X1、X5、X2、X4，其中显著影响因素为X1、X3，由图7 帕累托图（Pareto Chart）的t-value检验不仅可得到各因素显著性排序，还可观察到因素正负效应，其中X1、X5、X2对目标值影响为负效应，X3、X4为正效应；对极限断裂载荷Y2的影响因素从大到小依次为X5、X3、X1、X4、X2，其中影响显著因素为X1、X3、X5，且X1、X4对目标值的影响为负效应，X3、X5对目标值的影响为正效应。

剔除对目标值影响不显著因素，各因素与极限断裂位移Y1和极限断裂载荷Y2的模型为：

3.2 最陡爬坡试验

根据Plackett-burman试验结果可知，X2、X4对目标值影响均不显著，因此X2值为1.44×109N·m-3、X4值为1.26×106Pa。选取因素X1、X3、X5初值分别为3.6×109N·m-3、2.1×106Pa、1.6，步长分别为-1.5×108N·m-3、1.8×105Pa、0.16 进行最陡爬坡试验，以求出更靠近标准值的参数组合，最陡爬坡试验方案设计及结果如表6所示。由前文可知，茎秆弯曲破坏时极限断裂位移和极限断裂载荷分别为3.559 mm、98.944 N，表6 仿真试验结果中位移及载荷与响应值的误差均先减后增。其中第4组试验结果与响应值误差最小，分别为0.67%、3.63%，因此最佳参数组合取值范围在试验4附近。继续取试验3、试验4、试验5 参数组合分别为低、中、高水平进行后续Box-Behnken数值模拟试验。

表6 最陡爬坡试验方案及结果Table 6 Program and results of the steepest climbing test

3.3 Box-Behnken试验

Plackett-burman试验筛选出具有显著效应的参数，最陡爬坡试验确定最佳范围，现以试验3、试验4、试验5参数组合分别为低（-1）、中（0）、高水平（1）进行Box-Behnken 试验，得出试验因素水平设计，如表7所示。试验选取各中心点对试验误差进行评估，试验方案及结果如表8所示。方差分析结果如表9所示。

表7 试验因素水平Table 7 Level of test factors

表8 Box-Behnken试验方案及结果Table 8 Box-behnken test protocol and results

表9 Box-Behnken试验方差分析Table 9 Box-behnken test analysis of variance

由表9可知，3个因素对极限断裂位移Y1和极限断裂载荷Y2的二次回归模型均显著且其失拟项不显著，表明能够准确预测两个目标值变化趋势。其中对Y1影响显著性由大到小为X3、X5、X1X5、X3、X1，其中X1不显著。

为得到更为精确的参数组合，需缩小3个因素取值范围。已知Y1和Y2目标值分别为3.559 mm、98.944 N，X2取值为1.44×109N·m-3，X4取值为1.26×106Pa，利用Design-expert 11.0 软件对试验因素进行优化[26-27]。得到最优参数组合X1、X3、X5分别为3.1277878×109N·m-3、2.657545×106Pa、2.1，此时Y1、Y2分别为3.575 mm、99.12 N，与目标值相对误差分别为0.45%、0.13%。

4 参数优化试验验证

4.1 试验条件

为进一步验证优化得到的大豆植株Bonding 模型黏结参数组合正确性，选取4LZ-1.6Z 型大豆收获机上常用纵轴流脱粒装置进行脱粒性能试验，搭建试验台架如图8a 所示，试验材料为提前在田间采收的南夏豆25 植株。每组试验设定喂入量为1.5 kg·s-1，总喂入时间2 s，在不同脱粒转速下试验，每组试验重复3 次取平均值。采用3D 建模软件Solidworks 对脱粒装置进行建模，如图8b 所示。试验选取脱出物夹带损失率为试验指标，脱出物收集盘位置如图9所示，具体计算方法如下：式中，Y3-夹带损失率（%）；m1-前收集盘中收集籽粒总质量（g）；m2-后收集盘中收集籽粒总质量（g）。

图8 脱粒台架试验Fig.8 Threshing bench test

图9 收集盘位置示意Fig.9 Collection disk position

4.2 试验结果

在滚筒转速分别为350、400、450 r·min-1时得到脱粒装置物理试验与仿真试验夹带损失率如图10 所示。可见，仿真试验夹带损失率与物理试验平均误差为7.07%，表明本文标定的Bonding 模型参数组合准确度高。

图10 仿真试验与台架试验得损失率对比Fig.10 Comparison of loss rate between simulation test and bench test

籽粒分布情况如图11所示。为更好统计大豆籽粒轴向分布规律，对前收集盘作轴向编号如图12所示，籽粒轴向分布数据如图13 所示。可见，大豆籽粒经脱粒装置脱粒分布大致相同，均呈递减趋势。

图11 籽粒分布情况Fig.11 Grain distribution

图12 前收集盘轴向编号Fig.12 Front collection disc axial number

图13 籽粒质量轴向分布Fig.13 Axial distribution of grain mass

大豆植株离散元模型在脱粒元件与凹板筛打击、搓擦作用下，跟随钉齿旋转同时进行轴向流动，当外加载荷超过黏结键强度时产生黏结键断裂，物料运动过程如图14 所示，表示物料在不同时刻状态，包括在凹板筛低处压缩及导流板处抛散。

图14 脱粒过程中大豆茎秆运动过程Fig.14 Movement of soybean stem during threshing

由图14 可见，物料随脱粒元件流动同时，易堆积在凹板筛右下部，此时物料受到打击、搓擦作用最强，最易破碎，造成物料破碎率高、脱粒功耗增大，导致含杂率增加、脱粒滚筒易堵塞。仿真情况与物理试验时情况一致，进一步验证本文标定的仿真参数准确性。

5 结 论

a.为测得所脱物料力学性能，利用质构仪与三点弯曲试验探头，通过弯曲试验得到大豆植株极限断裂位移和极限断裂载荷分别为3.559 mm、98.944 N。

b. 为对所建植株模型粘结参数进行标定，利用Placket-burman试验找出显著因子、最陡爬坡试验缩小参数区间、Box-Behnken试验确定最优参数组合。Placket-burman 试验结果表明，对大豆茎秆极限断裂位移有显著影响的因素为：单位面积法向刚度和临界法向应力；对极限断裂载荷有显著的影响因素为：单位面积法向刚度、临界法向应力和黏结半径比例。Box-Behnken 试验结果表明，单位面积法向刚度、临界法向应力和粘结半径比例对极限断裂位移影响显著，临界法向应力和粘结半径比例对极限断裂载荷影响显著。通过优化工具求解得到当单位面积法向刚度、切向刚度、临界法向应力、切向应力、粘结半径比例分别为3.1277878×109N·m-3、1.44×109N·m-3、2.657545×106Pa、1.26×106Pa、2.1时，极限断裂位移和极限断裂载荷分别为3.575 mm、99.12 N，与真实值相对误差分别为0.45%、0.13%。

c.为验证植株模型准确性，以相同条件进行大豆脱粒台架试验，对比仿真和物理试验损失率，两者平均误差为7.07%，证明所提出大豆植株Bonding模型参数可靠。