饲料油菜薹期收获茎秆破碎离散元仿真参数标定

廖宜涛 廖庆喜 周 宇 王在腾 蒋亚军 梁 方

(1.华中农业大学工学院， 武汉 430070； 2.农业农村部长江中下游农业装备重点实验室， 武汉 430070)

0 引言

油菜是我国区域分布最广、播种面积最大的油料作物，但受农村劳动力转移、生产过程机械化水平低、种植规模小、进口油料冲击等多重因素的影响，油菜纯油用经济效益不高[1]。油菜饲料化利用是油菜全价值链开发利用的新路径，主要是在冬春季节刈割抽薹或花期的油菜，全株用于鲜喂或青贮饲喂牲畜，其中薹期油菜鲜嫩多汁，特别适合鲜喂。在冬春季节，饲喂肉牛每日增加5 kg的新鲜油菜，日增重率可提高30%以上[2]；饲喂哺乳母羊，每日添加1 kg的新鲜油菜，能有效促进产后体重恢复，明显提高羊乳中的脂肪和总固含量[3]；以50%比例油菜，混合稻草、菌渣等为主原料发酵制成的全混合日粮饲喂湖羊，与经青贮玉米为主原料的全日混粮相比，湖羊生长速度和屠宰性状没有显著差异[4]。饲料油菜产量高、营养丰富、成本低，利用水稻、玉米、大豆等作物收获后的土地空闲期种植，可有效缓解牲畜在冬春季节的青饲料缺乏问题，具有良好的经济、生态和社会效益，极具发展潜力[5-6]。

饲料油菜生物量大，人工收获效率低、劳动强度大，机械化收获是促进饲料油菜快速发展的关键。前期研究表明，饲料油菜是含水率高、脆嫩易切断、易挤压破损的作物，使用现有饲草收获机械进行收获存在过度切碎成糊状、揉搓挤压营养流失、输送通道易堵塞等问题[7]。研制新型饲料油菜收获机械需要对其切碎、抛送、集料、卸料等多个部件进行结构参数与工作参数的优化设计。近年来，离散元法及其仿真软件EDEM在农业机械领域得到了广泛研究与应用[8-13]。应用离散元法研究饲料油菜茎秆与收获机械的作用机理，可为机具设计与优化提供理论依据，提高研发效率。

饲料油菜的茎秆脆嫩，表面蜡质层、角质层和表皮组织薄，茎秆内外未出现髓腔和木质化，可以用同一属性颗粒近似构建油菜茎秆物料模型。切碎后的油菜茎秆是散粒体颗粒，在EDEM软件中可以采用Hertz-Mindlin基本模型；可以通过添加粘结参数Hertz-Mindlin with bonding模型，模拟物料破坏过程[14-17]。离散元仿真模型参数包括材料泊松比、剪切模量和密度等本征参数，碰撞恢复系数、静摩擦因数和滚动摩擦因数等基本接触参数，以及模型的法向与切向接触刚度、法向与切向应力等粘结参数，参数的准确性直接影响仿真结果的可靠性[18-21]。

饲料油菜是一种新型青饲料作物，其物料特性与已有研究物料差异较大，相关研究较少，且缺乏类似物料参考。本文以薹期饲料油菜茎秆为研究对象，以EDEM仿真软件为平台，先采用Hertz-Mindlin基本模型，结合饲料油菜茎秆堆积角测试物理试验，通过因素显著性筛选和响应曲面试验，确定碰撞恢复系数、静摩擦因数和滚动摩擦因数等离散元仿真模型基本接触参数；再采用Hertz-Mindlin with bonding接触模型，结合饲料油菜茎秆弯曲破坏试验，通过响应曲面试验确定法向与切向接触刚度、法向与切向应力等粘结参数。以期为薹期饲料油菜机械化收获过程的离散元仿真分析提供基本参数，并为茎秆类作物切碎过程离散元仿真研究提供模型参数标定方法。

1 材料与方法

1.1 参数测定

1.1.1饲料油菜本征参数

试验选用华中农业大学现代农业试验基地种植的“华油杂62”优质双低油菜，该品种用于饲用，具有产量高、营养指数高的特点[22-23]。种植时采用2BFQ-6型油菜精量联合直播机播种，播种量5.25 kg/hm2，同步施复合肥675 kg/hm2；平均生长密度46.5株/m2。由于油菜茎秆底部粗壮、顶部细弱，测量统计其茎秆平均直径为15.18 mm；通过体积计算和质量测定计算薹期茎秆密度为1 049 kg/m3。

选取油菜茎秆植株中间段直径约15 mm的茎秆，制作长度为20 mm的圆柱颗粒，利用TMS-Pro型质构仪(FTC公司，美国)进行平板单轴压缩试验，设置压缩速度为5 mm/min，加载位移为4 mm，试验重复10次；通过单轴压缩试验前后高度和直径的变化，计算得到盛花期饲料油菜茎秆弹性模量、剪切模量和泊松比分别为16.61 MPa、5.89 MPa和0.41，变异系数分别为9.65%、9.21%和7.77%。

1.1.2茎秆堆积角

茎秆堆积角通过圆筒提升法试验[24]。获得根据饲料油菜鲜喂及青贮的切碎长度要求，将抽薹期饲料油菜的茎秆预制为长度30、40、50、60 mm共4种颗粒样品，按等质量比例混合均匀；试验时，将钢质圆筒(内径300 mm、高150 mm)固定在万能材料试验机上，使底面与试验台接触，向钢质圆筒内填充饲料油菜颗粒样品直至填满，将钢质圆筒以0.05 m/s的速度向上提升，使样品形成一个颗粒堆，如图1所示。试验重复5次，测量颗粒堆堆积角，统计平均值θ0为33.10°，变异系数为1.27%。

图1 饲料油菜茎秆颗粒堆积试验Fig.1 Angle of repose for fodder rape stalk particles

1.1.3茎秆破碎物理参数

为获得反映粘结模型的法向与切向刚度、法向与切向应力准确参考值，利用TMS-Pro质构仪(FTC公司，美国)以有支撑切割方式慢速对茎秆施加载荷，记录加载过程弯曲破坏最大载荷作为茎秆破坏力参考值。对薹期饲料油菜茎秆进行预制，获得长度60 mm、平均直径15 mm的饲料油菜茎秆进行试验，固定支撑两点间距为40 mm，以0.003 m/s速度进行加载，如图2所示。试验重复5次，统计最大破坏力平均值为49.94 N，变异系数为4.41%。

图2 饲料油菜茎秆破坏试验Fig.2 Bending broken test of fodder rape stalk

1.2 仿真模型

1.2.1茎秆颗粒堆积角仿真模型

在EDEM软件中，应用Hertz-Mindlin模型进行饲料油菜茎秆堆积仿真试验。该模型在计算饲料油菜茎秆颗粒与颗粒之间的运动与受力时，实质上是将饲料油菜茎秆颗粒之间的运动与受力过程进行分解，分别为饲料油菜茎秆颗粒间的法向运动、饲料油菜茎秆颗粒间的切向运动和饲料油菜茎秆颗粒间的滚动运动。饲料油菜颗粒近似圆柱体，为匹配实际圆筒提升试验中饲料油菜茎秆堆积状态，采用球颗粒组合的方式[25-26]，建立如图3所示不同长度饲料油菜茎秆模型，为简化模型，茎秆颗粒直径取15 mm。结合离散元仿真农业物料和藤茎类物料离散元仿真参数[27-28]，确定本研究中各仿真参数的变化范围如表1所示。

图3 饲料油菜茎秆颗粒几何模型Fig.3 Particle model of fodder rape stalk

表1 茎秆颗粒堆积角仿真模型参数Tab.1 Parameters of angle of repose simulation model for fodder rape stalk particles

注：a表示该项为物理试验测定，b表示该项为查阅文献获得，c表示该项为试验变量，其范围为该变量取值的上下限。下同。

在EDEM软件中，建立与圆筒一致的虚拟圆柱体作为颗粒工厂(内径300 mm、高150 mm)，颗粒工厂采用动态生成，生成速度为0.4 kg/s，生成时间为4 s，饲料油菜茎秆颗粒总量为1.6 kg。为了保证饲料油菜茎秆颗粒快速稳定，设定颗粒初始下落速度为1 m/s，颗粒生成完成且稳定后，以0.05 m/s的速度向上提升圆筒，仿真颗粒从圆筒底部流出，最终形成一个稳定的颗粒堆，模型如图4所示。

图4 茎秆颗粒堆积角仿真模型Fig.4 Simulation model of angle of repose for fodder rape stalk particles

1.2.2茎秆破碎仿真模型

饲料油菜在收获中涉及到切断过程，采用Hertz-Mindlin模型建立颗粒模型是一个整体，进行断裂过程数值计算无效，无法用于饲料油菜切碎过程模拟仿真。在Hertz-Mindlin模型基础上，加入了颗粒间粘结作用，即采用Hertz-Mindlin with bonding模型，形成具有一定机械强度的粘连颗粒体模型。以该颗粒体模拟生成秸秆，当颗粒间承受的应力达到设置的极限应力时，粘结就被破坏，破坏后的粒子间不再存在粘结力，可用于物料的破碎过程仿真研究。

建模过程如图5所示。先用小颗粒粘结法生成茎秆模型，替代球形颗粒组合法建立的茎秆模型，具体流程是在EDEM中建立直径2 mm的小颗粒模型和直径为15 mm、长度为60 mm、关于原点中心对称的圆柱体模型，使用小颗粒对圆柱体进行填充[29]；再通过三维软件Pro/E建立刀具模型，导入EDEM中；最后在EDEM中以刀具为中心建立两支撑板，支撑板之间间距为40 mm，设定刀具运动速度为0.003 m/s，方向垂直向下，建立茎秆弯曲破坏仿真模型，设置仿真计算固定时间步长为4×10-7s。仿真模型基本接触参数采用已确定的Hertz-Mindlin模型参数，粘结参数参考藤茎类物料仿真参数[25]，取值范围如表2所示。

图5 饲料油菜茎秆弯曲破坏仿真模型Fig.5 Bending broken models of fodder rape stalk

表2 茎秆弯曲破坏仿真模型粘结参数Tab.2 Bonding parameters of bending broken simulation model for fodder rape stalk particle

1.3 参数标定试验设计

1.3.1基本接触参数标定

(1)二水平因子试验设计：应用Design-Expert 10.0.4软件进行二水平因子试验设计，筛选显著影响堆积角的参数。根据表1给定的参数建立茎秆颗粒模型开展堆积角仿真试验，其中x1～x6分别选择上、下限两个数值作为高、低两个水平；设计16次试验。

(2)最陡爬坡试验：针对筛选出的主要参数，进行最陡爬坡试验，以确定最优值临近区域。其中非显著性基本接触参数取二水平因子试验中的中间水平，显著性参数按设定步长逐步增加，建立茎秆颗粒模型开展堆积角仿真分析，记录仿真试验得到的堆积角与物理试验相对误差变化，根据相对误差的变化趋势确定最优值临近区域。

(3)基本接触参数响应曲面试验：基于二水平因子试验和最陡爬坡试验结果，根据Box-Behnken Design(BBD)设计原理，取显著性参数最优值临近区域的高水平值、中心点和低水平值为响应曲面试验设计的高、中、低3个水平，开展堆积角仿真试验。仿真模型中的非显著性参数取值同最陡爬坡试验。

1.3.2粘结参数响应曲面试验

根据Hertz-Mindlin with bonding模型计算原理可知，颗粒间粘结键的断裂与法向接触刚度、切向接触刚度、临界法向应力、临界切向应力等4个粘结参数相关。为保证参数标定过程中参数范围的可靠性，避免参数取值超过范围造成不良影响，根据中心组合设计(Central composite design,CCD)原理，结合表2中参数上下限数值设计响应曲面试验；以标定的基本接触参数和设定的粘结参数建立模型，开展饲料油菜茎秆破碎仿真试验与分析，粘结参数编码如表3所示。

表3 粘结模型参数编码Tab.3 Coding of bonding parameters

2 结果与分析

2.1 基本接触参数标定

2.1.1堆积角影响参数筛选

对影响茎秆颗粒堆积角的基本接触参数：饲料油菜茎秆间碰撞恢复系数x1、饲料油菜茎秆间静摩擦因数x2、饲料油菜茎秆间滚动摩擦因数x3、饲料油菜茎秆-钢碰撞恢复系数x4、饲料油菜茎秆-钢静摩擦因数x5、饲料油菜茎秆-钢滚动摩擦因数x6作为试验因素，进行二水平部分因子试验设计，试验设计与仿真结果如表4所示，方差分析得到各参数的影响效果如表5所示。

表4 二水平因子试验设计与结果Tab.4 Regular two-level factorial design and results

表5 参数显著性分析Tab.5 Analysis of parameter significance

由表5可知，根据影响率进行显著性排序，在饲料油菜茎秆堆积试验中，x2、x3和x5对堆积角有显著性影响，其余参数影响极小。因此在后续最陡爬坡试验和基本接触参数响应曲面试验中只考虑这3个影响显著的参数，其余非显著性因素取中间水平，即x1取0.5、x4取0.5、x6取0.1。

2.1.2最陡爬坡试验

由于x2、x3和x5对堆积角的效应均是正值，即堆积角随参数增加而增加，故设计参数值递增，观察并记录各试验中堆积角和相对误差，试验设计与结果如表6所示。

在3个试验因素数值逐渐增大的过程中，堆积角随之增大，与实际物理试验的堆积角间的相对误差呈先减后增的趋势，3号试验中相对误差最小，故以3号试验中的各参数数值作为后期试验的中心点，2号试验和4号试验分别为低水平和高水平进行后续BBD响应面试验设计。

表6 最陡爬坡试验设计与结果Tab.6 Experimental results of the steepest ascent test

2.1.3响应曲面试验结果与分析

根据前面筛选试验得出的结果，基本接触参数取x2、x3和x5，利用Design-Expert 10.0.4进行三因素三水平响应曲面试验设计，中心水平设置5组重复，总共进行17组饲料油菜茎秆堆积角仿真试验。试验设计与仿真结果如表7所示。

表7 接触参数响应曲面试验设计与结果Tab.7 Response surface test design and results of basic contact parameters

应用Design-Expert 10.0.4软件对响应曲面试验结果进行拟合分析，不同数学模型拟合效果对比表明，使用二次全模型方程拟合效果较好，模型决定系数R2为0.974 7，模型显著(P=0.001 0)且失拟项不显著，方差分析结果如表8所示，方程为

(1)

由表8可知，二次全模型P小于0.01，且失拟项P大于0.05，证明模型可靠，其中x2、x3和x5以及x3和x5的平方项对堆积角的影响都显著，但各个因素之间的交互项的P值都比0.05大得多，为了保证响应曲面的生成和堆积角θ的准确性，保留所有因素建立二次多项式回归方程。

表8 堆积角响应曲面二次全模型方差分析Tab.8 Analysis of variance of angle of repose for response surface quadratic model

通过响应曲面分析所生成堆积角模型，在已知x2、x3和x5的情况下可以确定饲料油菜茎秆颗粒堆积角，但在已知堆积角的情况下，x2、x3和x5的具体数值是堆积角的等值面曲线，因此通过Design-Expert 10.0.4软件以实际试验堆积角与仿真试验堆积角之间相对误差δθ为响应值开展分析。结果表明二次全模型方程拟合较好，决定系数R2为0.913 7，方程为

(2)

方差分析表明二次多项式回归模型P为0.000 4，且失拟项P大于0.05，模型可靠，其中除了x2和x5的交互项P远大于0.05之外，其余因素都具有较好的显著性。

通过Design-Expert 10.0.4软件约束求解工具对式(2)求误差最小极值点，得x2、x3和x5分别为0.43、0.05和0.29。将求解出的参数代入式(1)，可得出堆积角拟合值θ=34.28°，此时计算值与物理实现测定值θ0的相对误差为3.57%，表明建立的饲料油菜茎秆颗粒堆积角计算回归模型较好。

为验证求解的参数的可行性与准确性，利用计算所得的参数建立堆积角仿真模型开展仿真试验，5次重复试验测得饲料油菜茎秆平均堆积角为θs=32.35°，与θ0之间的相对误差δθ=2.27%，模拟结果与实际试验结果接近，表明确定的饲料油菜颗粒基本接触参数可行。

2.2 粘结参数标定

利用确定的基本接触模型参数值和CCD响应面试验设计确定的法向接触刚度x7、切向接触刚度x8、临界法向应力x9、临界切向应力x10等4个粘结参数值，开展饲料油菜茎秆弯曲破坏仿真试验，中心水平设置3组重复，总共进行27组仿真试验。试验设计与仿真结果如表9所示。

表9 粘结参数响应曲面试验设计与结果Tab.9 Response surface test design and results of bonding parameters

对响应曲面试验结果进行拟合分析，结果表明使用二次全模型方程拟合时，决定系数R2为0.995 4。方差分析(表10)可知，二次全模型P小于0.01，且失拟项P大于0.05；其中法向接触刚度x7、切向接触刚度x8及其二者平方项、交互项，临界法向应力x9和临界切向应力x10的平方项对弯曲破坏力有显著影响。

方程影响因子较多，在保证模型显著、失拟项不显著的情况下，剔除不显著项，对二阶回归模型进行优化调整，得到新的回归方程为

(3)

表10 弯曲破坏力响应曲面二次全模型方差分析Tab.10 Analysis of variance of bending broken force for response surface quadratic model

以破坏力与实测破坏力之间的误差Δf为响应值进行响应曲面分析，采用二次全模型方程进行拟合，方差分析结果如表11。

表11 弯曲破坏力误差响应面二次全模型方差分析Tab.11 Analysis of variance of error of bending broken for response surface quadratic model

各因素对破坏力相对误差Δf的影响显著性与对破坏力的影响显著性有较大差异，因素中临界法向应力x9、临界切向应力x10以及与这两项相关的交互项和平方项都不显著，因此临界法向应力x9和临界切向应力x10都取中心水平值，分别为40 MPa和5 MPa。

剔除不显著项，优化调整后建立破坏力误差Δf与各显著性因素之间的二次多项式回归方程为

(4)

对方程求最小值，得x7=4.60×109N/m，x8=3.55×108N/m。将参数代入式(4)中得Δf=1.61 N，相对误差为3.22%；将参数代入式(3)中得f=50.95 N，相对误差为1.98%；表明建立的饲料油菜茎秆颗粒破坏力计算回归方程可信。以计算所得的参数进行仿真试验，5次重复平均破坏力为50.64 N，与实际试验破坏力之间的相对误差为1.38%，模拟结果与实际试验结果无明显差异，标定参数有效。

3 茎秆破碎仿真试验验证

为了验证标定接触参数和粘结参数的可行性与准确性，确定建立的茎秆破碎离散元模型对于抽薹期饲料油菜茎秆具有普遍适用性，以抽薹期饲料油菜茎秆平均直径15 mm为中心，分别建立直径为13、14、16、17 mm的饲料油菜茎秆，模型参数按确定的参数(表12)进行设置，其他设置条件不变，建立离散元模型，与实际试验对比验证。

表12 饲料油菜茎秆离散元仿真参数Tab.12 Parameters of fodder rape crop for discrete element simulation

试验效果如图6所示，抽薹期不同直径饲料油菜茎秆破碎状态无明显差异，仿真试验与物理试验破坏情况一致，表明建立的模型可行。虚拟仿真与物理试验测定的饲料油菜茎秆破坏力如表13所示。

图6 茎秆破坏仿真试验与实际试验对比Fig.6 Stalk broken simulation experiment and physical experiment

表13 破坏力结果对比Tab.13 Comparison results of broken force

不同直径饲料油菜茎秆破坏力，仿真试验和实际试验结果相对误差不大于4.21%，且破坏力与茎秆直径具有良好的线性关系，表明参数标定方法正确，建立的破碎模型准确。

4 结论

(1)试验测得薹期饲料油菜茎秆平均直径为15.18 mm，油菜茎秆密度为1 049 kg/m3，茎秆弹性模量、剪切模量和泊松比平均值分别为16.61 MPa、5.89 MPa和0.41；通过物理试验测得油菜茎秆颗粒的堆积角平均值为33.10°，茎秆弯曲最大破坏力平均值为49.94 N。

(2)饲料油菜茎秆间碰撞恢复系数、饲料油菜茎秆-钢碰撞恢复系数及饲料油菜茎秆-钢滚动摩擦因数对堆积角影响较小，其取值确定为0.5、0.5和0.1；饲料油菜茎秆间静摩擦因数、饲料油菜茎秆间滚动摩擦因数及其平方项、饲料油菜茎秆-钢静摩擦因数及其平方项对堆积角的影响均显著。以实际试验堆积角与仿真试验堆积角之间相对误差为响应值进行分析，建立相对误差的二次回归模型，以相对误差最小为目标，优化求解得出饲料油菜茎秆间静摩擦因数为0.43、饲料油菜茎秆间滚动摩擦因数为0.05、饲料油菜茎秆-钢静摩擦因数为0.29，利用基本接触参数建立模型，进行堆积角仿真试验，仿真结果与实测值相对误差为2.27%。

(3)通过颗粒替代建立饲料油菜茎秆破碎模型，根据CCD原理，建立了粘结参数与破坏力之间的二次回归模型，由方差分析可知，法向接触刚度、切向接触刚度及其二者平方项、交互项，临界法向应力和临界切向应力的平方项对破坏力有显著影响。以破坏力仿真值与实测值之间的误差为响应值进行分析，建立二次回归模型，优化求解得到法向接触刚度、切向接触刚度、临界法向应力和临界切向应力分别为4.60×109N/m、3.55×108N/m、40 MPa和5 MPa。利用粘结参数进行仿真试验，与实际破坏力之间的相对误差为1.38%。

(4)以标定的参数建立不同直径饲料油菜茎秆弯曲破坏仿真模型，仿真试验与实际试验破坏情况一致，两者相对误差不大于4.21%，且破坏力与茎秆直径具有良好的线性关系，说明标定的参数准确、可靠，建立的破碎模型正确、可行。

本研究以物理试验为参考，通过堆积角试验和对茎秆缓慢施加载荷的弯曲破坏试验，确定了饲料油菜薹期收获茎秆破碎离散元模型参数，为建立茎秆类作物切碎过程离散元仿真模型提供了一种方法。受油菜茎秆生物体物料的复杂性影响，茎秆弯曲破坏试验的取样仅以抽薹期茎秆中间段、直径约为15 mm的物料为代表，标定参数具有一定的适用范围。后续将针对不同生长期和不同部位茎秆进行系统性研究，并结合饲料油菜收获仿真分析，对模型参数进行修正，对饲料油菜收获切碎、输送关键装置进行优化改进。