基于离散元法的油莎豆降阻挖掘装置设计与试验

何晓宁 张学军 赵 壮 尚书旗 王东伟 袁小伟

(1.新疆农业大学机电工程学院， 乌鲁木齐 830052； 2.青岛农业大学机电工程学院， 青岛 266109)

0 引言

油莎豆是原产于非洲及地中海沿岸国家的一种草本植物[1-2]，是集油、粮、牧、饲、绿化观赏于一体的新型经济作物，具有适应性广、含油量高、营养价值高等特性[3-4]。目前，我国大多数地区油莎豆收获以正旋旋耕和挖掘铲式作业方式为主，但是由于正旋挖掘方式存在埋果率高、壅土严重，挖掘铲式存在挖掘阻力大、土壤破碎率低等问题[5-6]，无法满足油莎豆规模化收获要求。国内外学者近年来对旋耕机理进行深入研究，ASAF等[7]研究了旋耕刀动力预测模型，推导了旋耕刀对土壤切削时的支反作用力和扭力方程，阐明了挖掘阻力产生的原理，ROSA等[8]利用有限元方法模拟出刀具耕作时的工作情况，在以不同耕速比为变量前提下，得出刀具不同结构参数对于性能指标的影响，LEE等[9]基于旱地水稻直播机的旋耕部件，研究旋耕刀在不同结构参数和工作参数下对旋耕作业扭矩、碎土率的影响。

文献[10-11]基于深松部件、旋耕部件作业之间的交互作用，设计了一种用于深耕的反旋深松联合作业耕整机，相比于普通深松旋耕联合整地机，牵引阻力降低了16.21%，方会敏等[12]建立了秸秆-土-旋耕刀离散元模型，进行EDEM仿真试验和田间试验对比，结果显示秸秆位移随着转速增加呈现增加的趋势，秸秆的水平运动位移总是大于同转速下的侧向位移，同时阐明了土壤的破碎机理，郭俊等[13]为研究正、反旋两种旋耕方式在不同转速下土壤和秸秆位移情况，运用示踪法原理，对铝块和秸秆进行标记，试验结果表明反转旋耕破碎率优于正旋。其他相关研究主要针对土壤破碎性能设计一种旋耕装置，未结合油莎豆种植模式和收获模式，适应性差，无法满足油莎豆收获基本要求。

本文基于离散元原理与旋耕刀-油莎豆团聚体力学模型，以反旋作业方式进行运动学以及旋耕刀结构参数分析，建立油莎豆团聚体离散元模型，通过EDEM仿真试验，分析反转旋耕刀结构参数对试验指标影响规律，结合试验结果，设计一种高效破土旋耕刀，以期解决油莎豆联合收获过程中挖掘阻力大、土壤破碎率低、埋果率高等问题。

1 油莎豆植株特性分析

由于油莎豆植株成簇生长，收获期时根系发达且含土量高[14]，如图1所示。本文针对油莎豆种植模式进行试验分析，发现河南民权地区油莎豆种植特点：行距15～20 cm，株距12～16 cm，生长深度为12～15 cm，为挖掘装置结构参数设计提供依据。

为保证后续仿真试验的精准性和科学性，结合河南民权地区油莎豆种植土壤特点，通过五点取样法进行土壤性质测定，得到其不同深度的土壤含水率、密度以及土壤坚实度，如表1所示。

表1 土壤物性参数Tab.1 Soil physical parameters

2 旋耕刀反旋工作原理

2.1 旋耕刀运动分析

当反转旋耕刀切削土壤时，旋耕刀在任何时刻的绝对速度都是由机器前进速度和刀轴转速矢量合成，运动轨迹模型如图2所示，设A为旋耕刀刀尖，其坐标为(X,Y)，其运动轨迹方程为

(1)

式中vm——前进速度，m/s

R——刀轴回转半径，mm

α——弯折角，(°)

ω——刀轴角速度，rad/s

t——刀轴运动时间，s

β——旋耕刀转动角，(°)

λ——旋耕速比

对式(1)求导可得刀尖A(X，Y)处在x和y方向的分速度为

(2)

刀尖A处绝对速度为

(3)

联立式(2)、(3)可得

(4)

其中

λ直接影响旋耕刀的运动轨迹、切土节距以及机器功率损等工作性能。

切土节距计算式为

(5)

式中S——切土节距，mm

Z——刀轴上同一径向旋耕刀数量

n——旋耕刀转速，r/min

与正旋挖掘相比，机器以反旋的工作方式作业时，可以在低速下进行挖掘作业，保证机器在较低功耗下实现较高的土壤破碎率，提高机具的作业性能，同时土垡切削量小且稳定，作业后土壤颗粒层次分布，通透性好，利于后续耕整地作业[15]。

2.2 旋耕刀结构参数分析

在油莎豆反转挖掘收获过程中，旋耕刀正切刃和刀端弯曲部位首先开始切削土壤，然后正切刃由下至上切削土壤，正切面的主要功能为切土、碎土、抛土等，是反转旋耕刀重要组成结构，其结构形状和参数对于挖掘性能有着重要影响，其中影响正切面的结构参数主要有回转半径R、弯折角α以及单刀工作幅宽L[16]，如图3所示。本文结合油莎豆种植模式研究，将回转半径设为150 mm，通过改变弯折角α以及单刀工作幅宽L来改变正切面参数，进一步提高正切面的抛土性能，通过EDEM仿真探究其最佳结构参数组合。

3 离散元模型建立

3.1 土壤模型

结合前面对于土壤的特性分析，本文以参数平均值建立土壤颗粒床模型，在确保试验精度的前提下，将土壤颗粒进行适当放大，将土壤颗粒半径设为5 mm，同时为了确保土壤环境的真实性，采用颗粒组合的方式将土壤颗粒模型设置为以下4种模型[17-19]，即粒状颗粒、团粒状颗粒、块状颗粒、杆状颗粒，如图4所示。

3.2 油莎豆块茎建模

在实际油莎豆收获过程中，油莎豆块茎是挖掘过程中油莎豆团聚体主要组成部分，很多学者在模拟试验时，基本都忽略油莎豆块茎的影响[20]。为精准模拟在收获过程中油莎豆块茎与根系连接力破坏的影响规律，本文利用三维激光扫描技术确定油莎豆块茎的形状特征及三维尺寸，对油莎豆块茎的形状结构进行扫描、曲面重构，得出块茎的几何体模型，利用EDEM软件中的填充功能完成块茎离散元模型的建立，对油莎豆块茎三维尺寸进行试验分析，油莎豆块茎长为12 mm，宽为10 mm，高为8 mm，如图5所示。

3.3 油莎豆根系模型

为精准模拟油莎豆收获机在工作过程中旋耕刀与油莎豆团聚体的接触规律，建立油莎豆根系离散元模型。由于油莎豆根系具有柔软、细长以及方向不一致等特点[21]，本文对油莎豆根系建模进行简化处理，建立根系主要须茎，通过试验统计与分析，油莎豆根系长度范围为90～120 mm，半径为0.05～0.2 mm，结合对油莎豆根系长度的试验分析，通过颗粒组合的形式将油莎豆根系长度设为120 mm，如图6所示。

3.4 接触模型及参数设定

3.4.1接触模型

为精准模拟河南民权地区粘性土壤性质，土壤颗粒之间的接触模型选择Hertz-Mindlin with Bonding模型，如图7所示，该模型可以通过结合力将相邻2个土壤颗粒粘结在一起，结合力可以承受切向和法向位移[22-24]。在实际田间作业过程中，土壤与土壤之间存在粘附力，挖掘铲与土壤存在受力关系，该模型可以模拟土壤颗粒之间的粘结作用和土壤颗粒破碎的现象。

土壤颗粒间Hertz-Mindlin with Bonding接触模型计算式为

(6)

其中

(7)

式中Fn——法向接触力，N

Ft——切向接触力，N

Tn——法向力矩，N·m

Tt——切向力矩，N·m

vn——法向速度，m/s

vt——切向速度，m/s

kn——法向刚度，N/m

kt——切向刚度，N/m

ωn——法向角速度，rad/s

ωt——切向角速度，rad/s

A——单位接触面积，mm2

J——惯性矩，mm4

δt——时间步长，s

Ra——颗粒粘结半径，mm

土壤颗粒间法向和切向应力达到一定极值时粘结键将会断裂，设断裂时法向应力和切向应力分别为σmax、τmax，当粘结键断裂时，有

(8)

3.4.2离散元参数设定

为了确保仿真试验的准确性，本文通过参数试验以及查阅文献[25]，确定仿真试验所涉及的接触参数和基本物理参数，如表2所示。

3.5 虚拟仿真试验过程

为模拟油莎豆挖掘装置实际作业状态，根据河南民权地区油莎豆种植模式，运用EDEM建立长1 200 mm、高400 mm虚拟土槽，其中油莎豆株距120 mm，行距150 mm，如图8所示。

在虚拟仿真试验过程中，将旋耕刀前进速度设置为0.5 m/s，刀轴转速30 r/min，旋耕刀弯折角110°，旋耕刀工作幅宽30 mm；为保证仿真过程精确性和连续性，设置仿真步长为3.14×10-6[26-28]，网格单元尺寸9 mm，为最小颗粒半径的3倍，颗粒粘结半径为5.5 mm，颗粒生成时间0.4 s，仿真时长2 s，通过仿真试验可得到旋耕刀对土壤的破碎情况，从EDEM后处理工具中导出旋耕刀所受阻力以及土壤破碎率的实时数据。

表2 离散元参数Tab.2 Discrete element parameter values

4 仿真试验

4.1 挖掘过程动态分析

工作阻力是衡量联合收获机综合性能的主要技术参数，在旋耕刀工作过程中，通过EDEM后处理分析不同时刻下挖掘铲与油莎豆团聚体破坏动态过程，如图9、10所示。

油莎豆反转旋耕作业过程中(图9、10)，旋耕刀正切刃和刀端弯曲部位首先开始切削底层土壤，土壤垂直于旋耕刀尖上向上方及两侧运动，深层土壤被打破，随着旋耕刀不断对中部土层进行剪切、撕裂、翻转、抛掷等作用，土壤失效区域范围变大，呈松散状态，油莎豆块茎和油莎豆根系受到扰动，如图10所示，土壤呈红色状态，逐渐沿旋耕刀尖切线方向运动，在旋耕刀往复作用下油莎豆团聚体在被破坏的同时被抛起，完成反转挖掘过程。

4.2 单因素试验设计

为探究影响挖掘作业质量的主要因素取值范围，以弯折角α以及单刀工作幅宽L为试验因素进行单因素试验，以工作阻力和刀轴扭矩为试验指标，由上述旋耕刀结构参数分析可知，油莎豆反转作业质量与弯折角α以及单刀工作幅宽L因素密切相关，根据河南油莎豆种植农艺要求，选取弯折角为90°～170°，工作幅宽为25～65 mm研究其对工作阻力和刀轴扭矩的影响规律，为后续多因素试验提供理论基础。

4.2.1单因素试验

根据EDEM仿真数据，导出工作阻力和刀轴扭矩在不同仿真时间段的变化规律，如图11、12所示。在旋耕刀工作过程中，当旋耕刀未接触土壤时工作阻力和刀轴扭矩为0，随着刀轴的转动，旋耕刀开始切削土壤，所受工作阻力和刀轴扭矩开始逐渐增大，随着切削运动的进行，旋耕刀切削土壤的体积不断增加，当旋耕刀开始接触油莎豆团聚体，对油莎豆团聚体进行破坏，旋耕刀所受阻力瞬间增大，在旋耕刀与油莎豆团聚体破坏过程中耕深达到最大值时，所受阻力达到最大值；当旋耕刀越过最大耕深时，切土体积逐渐减少，所受阻力相应减少；由于旋耕刀在切削土壤时每转过一定的角度都有相同的4把旋耕刀进入土壤，因此刀轴在工作时受到的总阻力呈现周期性变化。

4.2.2弯折角α对试验指标的影响

固定旋耕刀前进速度0.5 m/s，刀轴转速30 r/min，旋耕刀工作幅宽30 mm，弯折角α分别取90°、110°、130°、150°和170°，进行5组仿真试验，以工作阻力和刀轴扭矩作为试验指标，观察旋耕刀所受三向阻力最大值的变化规律，如图13、14所示。在弯折角90°～150°范围内，弯折角的变化对于工作阻力影响较为显著，旋耕刀所受三向阻力逐渐增加，在150°～170°之间，弯折角的变化对于工作阻力影响较低，旋耕刀所受三向阻力变化较为平稳；当弯折角在90°～170°范围内，对于x方向刀轴扭矩先增加至51 N，然后逐渐降低趋于平稳，当弯折角在90°～150°范围内，对于y、z方向刀轴扭矩变化影响显著，在150°～170°区间范围内其刀轴扭矩变化趋于平稳。

4.2.3工作幅宽L对工作阻力的影响

固定旋耕刀前进速度0.5 m/s，刀轴转速30 r/min，弯折角90°，旋耕刀工作幅宽L分别取25、35、45、55、65 mm，进行5组仿真试验，旋耕刀所受三向阻力最大值随转速的变化情况如图15、16所示。当工作幅宽在25～55 mm区间范围内，x、z方向工作阻力逐渐增加，工作幅宽的变化对于x方向工作阻力影响较为显著，在55～65 mm之间，x方向工作阻力变化较为平稳；工作幅宽在25～55 mm区间范围内，对于y方向工作阻力变化影响显著，在55～65 mm范围其工作阻力有略微下降趋势，但阻力波动范围较小；当工作幅宽在25～55 mm区间范围内，对于x方向刀轴扭矩变化影响显著，但工作幅宽在55～65 mm内其刀轴扭矩略微下降，当工作幅宽在25～55 mm区间范围内，y、z方向刀轴扭矩逐渐增加，在55～65 mm之间，刀轴扭矩变化较为平稳。

根据以上分析，为得到最佳旋耕刀结构参数，选取弯折角110°～150°，工作幅宽35～55 mm为试验因素取值范围。

4.3 多因素仿真试验

在前期油莎豆团聚体离散元模型建立及旋耕刀动力学分析的基础上，为得到理想的旋耕刀结构参数组合，进行EDEM正交旋转虚拟仿真试验研究，结合前期理论分析和单因素试验结果，设定试验因素编码如表3所示。

表3 试验因素编码Tab.3 Test factors coding

以弯折角α、工作幅宽L为试验因素，采用二因素五水平二次正交旋转组合设计进行旋耕刀虚拟仿真试验，根据所得数据进行结果分析，并对影响指标的主要因素进行显著性分析，试验结果如表4所示，表中A、B为编码值。

表4 仿真试验结果Tab.4 Simulation test results

对表4数据进行二次多元回归拟合，得到工作阻力和刀轴扭矩对编码自变量的二次多元回归方程

y1=703.22+103.48A+33.72B+5.42AB+
222.50A2+47.92B2

(9)

y2=31.06+2.09A+3.19B+0.15AB-
1.37A2+1.96B2

(10)

同时经过Design-Expert软件处理后得到关于工作阻力和刀轴扭矩的方差分析结果如表5、6所示。

表5 工作阻力方差分析Tab.5 Working resistance analysis of variance

表6 刀轴扭矩方差分析Tab.6 Analysis of variance of cutter shaft torque

由表5可知，模拟的一次项弯折角A对工作阻力的影响显著，二次项A2对工作阻力影响极其显著。从单因素水平分析，各项因素对工作阻力的影响由大到小为弯折角A、工作幅宽B。响应面分析如图17所示。

由表6可知，模拟的一次项工作幅宽B对刀轴扭矩的影响极其显著，一次项弯折角A和二次项B2对刀轴扭矩影响显著。从单因素水平分析，各项因素对工作阻力的影响由大到小为工作幅宽B、弯折角A，响应面分析如图18所示。

4.4 试验优化

为得到试验因素最佳旋耕刀结构参数组合，对试验进行优化设计。建立参数化数学模型，结合试验因素的边界条件，以降低工作阻力和刀轴扭矩为目标，对旋耕刀工作阻力和刀轴扭矩进行回归分析，建立非线性规划参数模型为

(11)

利用Design-Expert对数学模型进行分析求解，弯折角130°、工作幅宽45 mm为最优参数组合，此时工作阻力为703.22 N，刀轴扭矩32.06 N·m，将最优参数进行EDEM虚拟仿真试验验证，其工作阻力为928.63 N，刀轴扭矩34.56 N·m，与优化结果基本一致。在相同参数设置下与目前正旋旋耕挖掘装置进行试验对比，结果表明，正旋挖掘装置工作阻力为1 015.52 N，刀轴扭矩为36.72 N·m，工作阻力降低了8.6%，刀轴扭矩降低了5.9%，验证了反旋式旋耕刀片具有一定的降阻作用。

5 田间试验

5.1 田间试验条件与方案

旋耕刀挖掘阻力和刀轴扭矩并不是衡量作业效果主要标准，由于油莎豆目前无相关收获作业标准，本文参考NY/T 502—2016中的收获标准以及GB/T 5668—2017中的土壤破碎率计算方法，将埋果率和土壤破碎率作为田间试验的指标，目的是检验优化设计的旋耕刀作业性能，结合理论分析与仿真模拟结果进行田间验证，试验地点选择河南省民权县油莎豆种植基地，该地区油莎豆种植模式为平作，行距为140 mm，株距为120 mm，油莎豆生长深度160 mm左右，如图19所示。

试验划为5个区域，每个试验区域长50 m，将新型旋耕刀和IT245旋耕刀分别安装在油莎豆联合收获机上，进行油莎豆收获性能对比分析。机器作业速度设置为0.5 m/s，将埋果率和土壤破碎率作为试验指标，进行试验分析，如图20所示。

5.2 田间试验结果与分析

试验结果如表7所示。

表7 田间性能试验对比Tab.7 Field performance test comparison %

由表7可知，新型旋耕刀作业时埋果率最高值为1.23%，最低值为0.85%，平均值为1.07%，土壤破碎率最高值为94.61%，最低值为92.67%，平均值为93.48%，标准旋耕刀埋果率最高为2.85%，最低为1.75%，平均值为2.27%，土壤破碎率最高为92.87%，最低为91.54%，平均值为92.18%，与标准旋耕刀相比，新型旋耕刀油莎豆块茎埋果率减低了1.20个百分点，土壤破碎率提高了1.3个百分点，参考NY/T 502—2016的规定，进一步提高了油莎豆的收获效率。

6 结论

(1)通过对旋耕刀进行运动学分析和结构分析，结合虚拟仿真试验，得到旋耕刀结构参数弯折角取值范围为110°～150°，工作幅宽取值范围为35～55 mm。

(2)建立了旋耕刀-土壤-油莎豆团聚体离散元模型，进行了油莎豆虚拟仿真挖掘试验，从微观角度分析了土壤破碎的动态过程和油莎豆团聚体的速度关系，根据Design-Expert试验分析，确定了旋耕刀的最佳结构参数：弯折角130°，工作幅宽45 mm，在相同参数设置下与正旋旋耕方式进行对比试验，工作阻力降低了8.6%，刀轴扭矩降低了5.9%，证明反旋挖掘具有降阻作用。

(3)为验证反旋刀片设计的合理性，以埋果率和土壤破碎率进行油莎豆收获性能试验，同时与标准旋耕刀进行试验对比，结果表明：反旋作业方式埋果率1.07%，土壤破碎率93.48%，与标准旋耕刀相比，新型旋耕刀油莎豆块茎埋果率减低了1.2个百分点，土壤破碎率提高了1.3个百分点。