联合收获油菜脱出物离散元仿真参数标定与试验

舒彩霞，杨 佳，万星宇，袁佳诚，廖宜涛，廖庆喜

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

0 引 言

油菜是中国主要油料作物，其生产种植面积占全国油菜总种植面积85%以上，全程全面高质高效推进油菜机械化发展是保障食用油安全的重要措施。油菜联合收获是油菜机械化收获的主要方式之一，可一次性完成切割、脱粒、分离、清选等工序，直接获得干净的油菜籽粒，生产效率高。在收获时，过早收获会产生脱粒不净、青籽多等问题；过晚收获容易造成落粒、损失率增加；因此，最适宜收获时期是田间油菜呈黄熟期后至完熟期之间，90%以上的油菜角果变成黄色和褐色，主分支向上收拢，果荚中黑籽粒在75%左右。

为提高油菜联合收获机清选装置的清选性能，国内外学者对清选装置进行了研究。油菜籽粒在进行清选作业时，脱出物中的籽粒、短茎秆及荚壳间的相互作用关系难以通过数学模型进行分析，采用离散元法（Discrete Element Method，DEM）和计算机流体力学（Computational Fluid Dynamics，CFD）对关键部件进行数值模拟可以有效模拟出物料运动轨迹，替代繁杂的台架试验，省时、省力，缩短研发周期，降低成本。在进行油菜脱出物清选仿真试验时，脱出物间精确的离散元接触参数和特征参数可以提高气固耦合清选仿真中颗粒间相互作用及运动规律的准确度，进而实现关键部件的参数优化。

为提高仿真模拟的准确性，使仿真结果更加接近实际作业过程，国内外学者基于堆积角试验针对不同物料的离散元接触参数进行了标定。Guo等通过物理试验与仿真试验相结合方法建立了香蕉茎秆力学模型并对其生物力学特性进行了研究。马彦华等采用6个接触参数测定试验，以此确定优化试验的参数区间，依次通过Plackett-Burman试验、最陡爬坡试验以及Box-Behnken 试验对苜蓿秸秆的离散元仿真参数进行标定。侯占峰等通过物理试验和仿真试验相结合的方法对冰草种子进行了参数标定。石林榕等通过堆积试验、查阅参考文献和计算了马铃薯基本物理参数和接触力学参数。综上所述，国内外学者为不同作物的数值模拟提供了可参考的离散元参数。由于油菜脱出物组分糅杂、成熟度不一致、物料间特性差异大，国内外学者围绕油菜脱出物的空气动力学特性、特征参数及组分配比进行了研究。陈立等对油菜脱出物中籽粒、茎秆、荚壳及轻杂的悬浮速度进行了测定。Zhan等对整株植物的抗裂荚性、抗主茎剪断性和共振频率等生物力学特性进行研究。雷小龙针对油菜籽粒的特征参数及接触参数进行了大量研究。袁佳诚等开展了脱出物质量比试验，发现随着切碎滚筒转速增加，脱出物中油菜茎秆平均长度逐渐降低，增加了脱出物中短茎秆的质量比，影响旋风分离清选装置的清选性能。相关研究可为仿真模拟的基础参数提供依据，但相关研究主要是对油菜单一组分的部分物理特性和油菜籽粒的部分特征参数及接触参数进行了研究，对于油菜脱出物（籽粒、短茎秆、荚壳）离散元接触模型参数鲜有报道。

本文采用物理试验与仿真试验相结合的方法进行相关参数标定，对油菜脱出物组分间的静摩擦系数，碰撞恢复系数进行标定，采用堆积试验对茎秆特征参数、茎秆间接触参数及茎秆-钢板接触参数进行标定，在优化标定后参数下，通过DEM-CFD气固耦合试验与台架试验验证了标定参数准确性，以期为油菜联合收获机械化仿真模拟提供离散元模型基础参数。

1 联合收获油菜脱出物物料特性测量与分析

1.1 试验材料

2021年5月于华中农业大学现代农业科技试验基地中采用自主研发的油菜联合收获机开展田间试验，试验材料为机直播华油杂62号油菜，平均种植密度为35～40 株/m，平均株高为1 655.44 mm，平均分支数为6。为准确反应茎秆的实际机械物理特性，采用五点取样法获得油菜茎秆，分别在1 m的测试区随机人工选取油菜植株，从离地350 mm开始，每隔100 mm截取油菜茎秆；收集脱出物（如图1所示）中荚壳及籽粒，将收集的材料分选后作为试验材料，按照其种类人工分选后并称量，其中籽粒、短茎秆和荚壳的质量比为2∶3∶5，含水率分别为30.41%～37.23%、40.92%～48.32%和66.67%～73.53%。

图1 油菜脱出物 Fig.1 The rape threshing mixture

1.2 茎秆特征参数测量与分析

茎秆材料特征参数包括密度、泊松比、弹性模量和剪切模量等。通过MNT951221型数显游标卡尺（测量误差为0.02 mm）对收集的油菜茎秆直径进行测量，测得油菜茎秆的平均直径为10.05 mm。根据文献[23]，通过体积计算和质量测定计算得到油菜茎秆密度为494 kg/m。

通过质构仪（TMS-Pro型，美国FTC公司）对随机选取的油菜茎秆进行径向单轴平板全压缩试验，试验采用直径75 mm 的平板压头，加载速度为50 mm/s，如图2a所示。茎秆出现破裂声时停止加载，采用定义法测量油菜茎秆泊松比；利用数显游标卡尺记录被挤压方向及与其垂直方向的直径变形量，使用质构仪记录载荷-位移数据（如图2b所示）。重复10次，通过式（1）～（3）计算得到油菜茎秆的弹性模量、泊松比和剪切模量平均值分别为4.44 MPa、0.45、1.53 MPa。

图2 油菜茎秆特征参数测量 Fig.2 Measurement of feature parameters of rape stalk

式中为物料的弹性模量，MPa；为最大压应力，Pa；为线应变；F为最大载荷，N；、为茎秆外径和内径，mm；为试样接触长度，75mm；Δ为茎秆直径变化量，mm；Δ为茎秆压缩前直径，mm。

式中为物料的泊松比；为被挤压方向应变；为与被挤压垂直方向应变；为油菜茎秆试验前被挤压方向初始直径，mm；Δ为油菜茎秆试验后被挤压方向直径变化量，mm；为油菜茎秆试验前与被挤压垂直方向的初始直径，mm；Δ为油菜茎秆试验后与被挤压垂直方向的直径变化量，mm。

式中为物料的剪切模量，MPa。

1.3 脱出物接触参数测量与分析

脱出物接触参数主要包括静摩擦系数、动摩擦系数及碰撞恢复系数等。油菜籽粒细小且脱出物生物学特性差异明显、形态各异。参数参考文献[14]，本文通过斜面法、自由跌落试验对脱出物组分间静摩擦系数和碰撞恢复系数进行测量，采用堆积角逼近的方式对茎秆与茎秆、钢板动摩擦系数进行标定，明确物料组分间的接触参数范围。

采用斜面法测量不同物料间的静摩擦系数。采用摩擦系数测定装置，待测板为300 mm×300 mm× 3 mm钢板，试验装置如图3所示。当物体静止于斜面上时：

图3 物料间静摩擦系数测量试验 Fig.3 Static friction coefficient test between materials

式中为物料的重力，N；为斜面与水平面的夹角，(°)；为待测物料与物料种群板的静摩擦系数。

通过不断增加斜面倾角，当物体刚好有滑动趋势时有：

此时斜面倾角正切值便是此时接触物体间的静摩擦系数。

在测定物料间静摩擦系数时，用镊子将油菜籽粒、荚壳、短茎秆等分别放至于物料种群板上，缓慢提高物料种群板，当物料出现向下运动时停止抬升物料种群板，记录此时的角度，通过式（5），即可得到物料间的静摩擦系数，每组试验重复20次。根据文献[14]，籽粒间静摩擦系数取0.5。在测量时，为了充分模拟物料与接触材料的接触特性，使试验数据有良好的适应性，物料每次放置位置尽量不重复，得到的油菜脱出物组分间的静摩擦系数如表1所示。

表1 接触物料间的静摩擦系数 Table 1 Static friction coefficient between contact materials

采用自由跌落试验测量不同油菜脱出物间碰撞恢复系数。将所测物料均简化为质点，物料从高度开始掉落，与地面接触板接触时（碰撞开始）质心速度为，受到固定接触面的碰撞冲量作用，物料变形增大，速度逐渐减小直至为0；此后弹性变形开始恢复，物料开始反向上升，物料离开接触面时的瞬间获得反向速度，并一直上升至高度处速度为0（碰撞结束）。碰撞过程如图4所示，由冲量定理可以求出碰撞恢复系数表达式，如式（6）所示。

图4 碰撞恢复系数测量原理图 Fig.4 Schematic diagram of collision recovery coefficient test

式中为物料质量，g；、为弹性变形前后的碰撞冲量，N·s；为碰撞恢复系数。

为清楚观察到物料的弹起高度，经过多次重复试验，随机选取茎秆制成20 mm的圆柱颗粒进行自由跌落试验，物料自地面接触板正上方200 mm处自由落体；分别以茎秆种群板、荚壳种群板以及钢板作为地面接触板进行测量。将网格纸粘于钢板上作为背景记录物料下降高度及弹起高度，采用美国VRI的Phantom系列高速摄影机正视放置于网格纸正前方，用镊子夹持物料自200 mm处释放，每组试验重复20次，将拍摄的照片导入电脑，通过式（6）计算出不同物料的碰撞恢复系数，如表2所示，根据文献[14]籽粒间碰撞恢复系数取0.6。

表2 接触物料间的碰撞恢复系数 Table 2 Collision recovery coefficient between contact materials

2 油菜联合收获茎秆颗粒离散元参数标定

由于油菜茎秆为各向异性材料，因此本文采用仿真逼近的方法对油菜茎秆接触参数进行标定。堆积角反应了颗粒物料流动及摩擦特性，与接触材料及自身物理特性相关。因此，堆积角物理试验常常被用来颗粒离散元参数标定。利用物理试验标定得到的参数范围开展油菜茎秆堆积试验，不断调整油菜茎秆间接触参数使仿真堆积角逼近实际堆积角，对应得到油菜茎秆间准确接触参数。

2.1 油菜茎秆实际堆积角测量

采用圆筒提升法获得油菜茎秆堆来测量堆积角。参考文献[32]，将收集的茎秆制成长度为20、30、40，50和60 mm茎秆颗料，经过多次重复试验，每种长度的茎秆各150 粒，共750 粒。试验所用的圆筒为内径200 mm、高度300 mm的钢制圆筒，将试样填满圆筒并放置于万能材料试验机平台，将圆筒以0.05 m/s匀速提升，试验重复3次。

正视拍照后将物理试验图片进行处理，应用Matlab中数字图像处理软件读取颗粒堆边缘图像，对图像进行灰度化、二值化，输出二值图，提取二值图边界轮廓，扫描轮廓图像上每一个像素，记录白点的坐标及个数，通过最小二乘法对所记录的白点进行线性拟合，在Matlab读取拟合直线的斜率，求得颗粒堆边缘轮廓的倾角。堆积角试验过程如图5所示，测得实际堆积角平均值为31.67°。

图5 堆积角左边缘拟合 Fig.5 Left edge fitting of actual repose angle

2.2 油菜茎秆模型建立及EDEM参数设置

在EDEM中采用Hertz-Mindlin无滑移模型进行油菜茎秆堆积试验（图6），为简化模型并提高仿真效率，采用圆形颗粒组合的方式形成与物理试验相同长度和直径的油菜茎秆，直径10 mm茎秆颗粒模型如图6a所示。通过多次预仿真试验并结合文献[14，34]及前文试验适当扩大试验参数取值范围，如表3所示。

表3 油菜茎秆颗粒堆积角仿真模型参数 Table 3 Parameters of repose angle simulation model for rape stalk particles

在EDEM仿真试验中，建立与物理试验参数相同的钢制圆筒与底板模型，在圆筒上方建立Polygon平面作为颗粒工厂，经过多次重复试验，设定每种颗粒生成速率为75 粒/s，共生产时间为2 s，颗粒产生完全并静置稳定颗粒堆后采用与物理试验相同的速度进行提升；仿真总时间为4.5 s，时间步长为Rayleigh时间步长的20%，网格尺寸为颗粒半径的3倍，EDEM仿真试验结果如图6b所示，仿真试验图像处理如图6c所示。

图6 堆积角仿真试验 Fig.6 Simulation test of repose angle

2.3 仿真试验

为减少试验次数，应用Design Expert 10.0.7进行Plackett-Burman筛选试验，以茎秆堆积角为响应值，筛选出对堆积角影响显著的参数。将表3中的～进行最大、最小值分别编码水平+1、-1，编码结果如表4所示。分别测量左右两侧的堆积角数值，然后取平均值，Plackett-Burman试验结果如表5所示。

表4 Plackett-Burman试验方案 Table 4 Plackett-Burman test scheme

表5 Plackett-Burman试验结果 Table 5 Plackett-Burman test results

用Design Expert 10.0.7进行方差分析，得到各个参数的影响效果，如表6所示。根据贡献率进行显著性排序，在油菜茎秆堆积试验中，，对堆积角的影响较为显著，因此只考虑这3个因素进行最陡爬坡试验；其余显著性影响较小的取值分别为：取0.4，取1.5 MPa，取0.277，取0.790，取0.391。

表6 Plackett-Burman试验参数显著性分析 Table 6 Significance analysis of Plackett-Burman test parameters

基于Plackett-Burman堆积试验筛选出3个显著性影响参数（茎秆-茎秆静摩擦系数，茎秆-茎秆动摩擦系数，茎秆-钢板动摩擦系数），以仿真堆积角与实际堆积角的相对误差作为评价指标来确定试验参数的最佳范围，试验方案及结果如表7所示。

表7 最陡爬坡试验设计方案及结果 Table 7 The steepest climbing test design scheme and results

3个试验因素逐渐递增，堆积角相对误差呈先减小后增加的变化趋势，其中4号的堆积角相对误差最小，为1.89%，基于最陡爬坡试验确定以4号试验中的各个参数作为后期试验的中心点，3号、5号作为低水平与高水平，通过Box-Behnken响应面试验验证该模型的显著性关系，并通过优化得出最佳参数组合。由式（7）计算出仿真堆积角与实际堆积角的相对误差。

利用Design-Expert 10.0.7进行三因素三水平响应曲面试验设计，共进行17组仿真试验，试验设计方案与结果如表8所示。对Box-Behnken试验结果进行多元回归拟合，得到油菜茎秆堆积角与3个显著性参数二阶回归方程为

表8 Box-Behnken试验设计方案及结果 Table 8 Box-Behnken experiment design scheme and results

由该模型方差分析结果可知（表9所示），茎秆-茎秆动摩擦系数对堆积角影响极显著；茎秆-茎秆静摩擦系数，茎秆-钢板动摩擦系数对堆积角影响显著；、、对堆积角影响极显著；对堆积角影响显著；、对堆积角影响不显著；该模型的＜0.05表明该模型的各个参数与响应值之间的关系显著；失拟项=0.356 2＞0.05，表明该模型拟合良好；变异系数（CV）为1.62%，说明本试验具有较好的可靠性；决定系数=0.966 6，校正决定系数=0.923 8，表明拟合方程可靠度较高；精密度（Adeq Precisior）为15.726，表明该模型具有良好的精确度。

表9 Box-Behnken试验回归模型方差分析 Table 9 Variation analysis of Box-Behnken quadratic model

根据回归模型方差分析可知，、这两个交互项对堆积角的影响不显著，交互项对堆积角的影响极显著。应用Design-Expert 10.0.7软件绘制三个显著性因素交互作用堆积角响应曲面，如图7所示。由图7a可知，当茎秆-茎秆间静摩擦系数一定时，随着茎秆-茎秆动摩擦系数的增加，堆积角变化趋势明显，对堆积角影响极显著，而茎秆-茎秆动摩擦系数一定时，随着茎秆-茎秆间静摩擦系数的增加，堆积角呈先增后减的趋势，对堆积角影响较小；由图7b可知，响应曲面为凸型曲面，茎秆-茎秆静摩擦系数和茎秆与钢板动摩擦系数对堆积角的影响均呈先增后减趋势，表明两因素对堆积角的影响基本相似；由图7c可知，茎秆与茎秆动摩擦系数对响应面趋势较陡，该因素比茎秆与钢板动摩擦系数影响显著。

图7 物料参数之间的交互作用响应曲面 Fig.7 Response surfaces of interaction between material parameters

利用Design-Expert 10.0.7软件的优化模块中以实际堆积角为目标值（31.67°）进行求解，对回归模型进行寻优，得到一组与物理试验测定结果较为相近的一组解为：茎秆-茎秆静摩擦系数为0.707，茎秆-茎秆动摩擦系数为0.015，茎秆-钢板动摩擦系数为0.012，其他非显著性参数取平均值。为验证离散元仿真的可靠性，以上述参数作为EDEM仿真参数，进行3组仿真试验，得到的休止角为30.90°、31.81°、31.80°，平均值为31.50°，与物理试验结果的相对误差为0.54%，验证了仿真试验的真实性与可靠性。

3 旋风分离清选仿真试验及验证

3.1 DEM-CFD标定参数验证

利用 Solidworks对旋风分离筒进行建模，在Workbench 17.0的mesh中采用四面体非结构化方法划分旋风分离器的网格，对进出口处网格加密处理，并进行边界层网格的划分，以保证其计算结果准确性，总体网格数约75万（如图8a）。对荚壳的模型进行三维扫描后在EDEM中进行实心颗粒填充；油菜籽粒简化为2 mm的球体；测量凹板筛分离出的油菜脱出物中短茎秆长度与直径，为简化计算量，结合文献[34-36]对茎秆颗粒进行简化，采用较为普遍应用的长度40 mm、直径3 mm的圆柱体短茎秆，如图8b所示。

图8 仿真网格及颗粒模型 Fig.8 Simulation grid and particle models

进料口为速度入口，吸杂口、出料口均为压力出口。EDEM和CFD的时间步长分别设置为5×10和1×10，气流介质为空气，其密度和黏性系数分别为1.225 kg/m和1.7894×10kg/(m·s)，重力加速度为9.81 m/s。物料颗粒力学参数及接触参数（非显著性参数取试验平均值）如表10和表11所示。旋风分离筒涉及复杂的密相气体-颗粒两相流动，在进行DEM-CFD耦合试验时，采用考虑颗粒对流体影响的Eulerian模型；颗粒相由EDEM 2018求解，为简化计算量、提高仿真效率，采用Hertz-Mindlin 无滑移接触模型，可以较好反应刚性物体间的力学行为；气相采用Fluent 17.0求解，Fluent将某一时间的流场计算至收敛后，通过曳力模型将流场信息转化为EDEM中作用在颗粒上的流体曳力，EDEM计算此时每个颗粒所受的外力（流体曳力、碰撞力等）从而更新颗粒此时的位置、速度、加速度等信息，最后这些信息以动量汇形式加载到CFD中，实现相间的耦合，由于清选筒内是三维强旋流流场，故选择RNG-模型，勾选适用于旋流的Swirl Dominated Flow选项。

表10 物料力学参数 Table 10 Mechanical parameters of materials

表11 物料接触参数 Table 11 Contact coefficient of materials

旋风分离筒内的气流分布规律及物料间的接触参数决定了筒内物料的运动状态，调节旋风分离筒的透明度，完成可视化的清选分离仿真过程显示。由图9可以看出物料在筒内的运动状态，物料从切向进入旋风分离筒内并绕着筒内壁向下运动，茎秆和荚壳在运动过程中被吸杂口处的负压吸走，籽粒沿着内壁不断向下运动从出料口处流出或者部分积聚在出料口的底部，总体来看数值模拟可以较为准确的反应籽粒及杂余的运动轨迹。在EDEM后处理模块Grid Bin Group统计出料口和吸杂口处的籽粒、荚壳、短茎杆的质量，通过式（9）计算得到旋风分离清选装置油菜籽粒清洁率为94.42%，损失率为3.96%。

图9 筒内物料运动状态及运动轨迹 Fig.9 Movement state and motion trajectory of the material in the cylinder

式中Y为油菜籽粒清洁率，%；为出料口油菜籽粒质量，g；为出料口油菜脱出物总质量，g；Y为损失率，%；为吸杂口油菜籽粒质量，g。

3.2 台架验证试验

本试验依托华中农业大学工学院自主研发的收获关键部件试验台开展台架验证试验，试验台主要包括输送带、割台、纵轴流脱粒分离装置和旋风分离清选装置等组成，如图10所示。收获关键部件试验台相关参数实时可调并且相关参数可以在线监测。

图10 台架试验 Fig.10 Bench test

试验结束后，收集出料口及吸杂口处物料进行人工分选并称量，通过式（9）计算油菜籽粒清洁率与损失率。入料口速度为6 m/s，吸杂口速度为13 m/s时，旋风分离清选装置油菜籽粒清洁率为91.84%，损失率为4.28%。仿真试验与台架试验的油菜籽粒清洁率与损失率相对误差分别为2.81%和7.48%。

4 结 论

1）研究得出了油菜联合收获脱出物的相关基础参数。开展了油菜脱出物特征参数、接触参数测量，测得油菜茎秆的泊松比、剪切模量分别为0.4、1.5 MPa；钢板与籽粒、茎秆、荚壳的静摩擦系数与碰撞恢复系数范围分别为0.214～0.390、0.711～0.869、0.647～0.781；0.543～0.667、0.212～0.570、0.141～0.283。茎秆与茎秆、荚壳的静摩擦系数与碰撞恢复系数范围分别为0.650～0.754、0.637～0.751；0.187～0.367、0.158～0.265。茎秆与籽粒的静摩擦系数与碰撞恢复系数分别为0.292～0.516；0.515～0.632；籽粒与荚壳的静摩擦系数与碰撞恢复系数范围分别为0.635～0.907；0.566～0.612，荚壳与荚壳的静摩擦系数与碰撞恢复系数范围分别为0.625～0.724；0.141～0.265。

2）开展了油联合收获油菜茎秆离散元参数标定，通过Plackett-Burman试验，筛选出对堆积角影响显著的参数为：茎秆-茎秆静摩擦系数、茎秆-茎秆动摩擦系数、茎秆-钢板动摩擦系数，通过Design-Expert 10.0.7软件优化模块得到显著性参数最佳组合参数为0.707、0.015、0.012，在最优组合参数下仿真堆积角平均值为31.50°，与实际堆积角的相对误差为0.54%。

3）结合标定的油菜脱出物不同组分间接触参数，开展了基于DEM-CFD旋风分离清选过程数值分析和台架验证试验。结果表明：旋风分离清选仿真试验油菜籽粒清洁率为94.42%，损失率为3.96%；仿真试验与台架试验籽粒清洁率与损失率相对误差分别为2.81%和7.48%。表明标定的油菜脱出物接触参数可为油菜联合收获机具的离散元仿真分析提供基础参数。

联合收获油菜茎秆取材和参数标定还需进一步扩大样本，后续将持续开展不同地区、不同部位油菜茎秆力学特性的具体性研究；荚壳易碎、呈长扁内凹壁结构，常规力学测试方法难以测出其特征参数和接触参数；后续将结合高速摄影技术和能量守恒定律对荚壳及脱出物不同组分间接触参数进一步标定。