基于EDEM 的粉垄和旋耕作业混肥效果研究

单发科， 康朔， 朱建锡， 王永维， 王俊*

（1.浙江大学生物系统工程与食品科学学院，杭州 310058； 2.浙江省农业机械研究院，浙江 金华 321051）

土壤施肥是提高作物产量的重要措施，化肥颗粒的分布范围和作物根系的吸收范围是否匹配及化肥分布的均匀性是影响化肥利用率的重要因素[1-2]。利用耕作机械将撒施在土壤表面的化肥颗粒与土壤混合是常用的施肥方法，不同耕作方式下的化肥分布存在差异，对作物生长产生不同影响[3-5]。石丽红等[6]研究发现，在水稻生产中，适当的肥料分布有利于提高水稻产量；张拥兵等[7]研究发现，不同施肥方式对葡萄果实品质影响较大，而施肥方式的不同也导致肥料分布的不同；雷金银等[8]研究发现，不同耕作方式改变了养分在土壤内的分布，提高了化肥利用率。化肥在土壤中的分布范围影响作物生长，而耕作方式则直接影响化肥颗粒在耕层内的分布[9-11]。因此，有必要研究不同耕作方式对化肥分布的影响。

离散元分析近年来被广泛应用于农业机械的研究中，目前针对施肥效果的仿真研究也多采用离散元法（discrete element method，DEM）。EDEM（expanded DEM）软件近年来被广泛应用于机械优化设计和作业过程仿真模拟等[12-13]研究，并取得一定进展。在EDEM 中，土壤由离散的颗粒组合而成，可以通过定义颗粒的大小、尺寸、位置组合以及结合方式等模拟不同的土壤特性[14-16]，DEM 是模拟耕作过程中机具与土壤相互作用的有效方式[17-18]。化肥颗粒为球状，在EDEM 中可以很方便地建立其仿真模型。

通过田间试验的探究不同耕作方式的混肥效果时步骤复杂且耗时耗力，为研究带来极大困难。因此，本研究通过分析化肥颗粒分布规律建立仿真预测模型并对旋耕、粉垄2 种耕作方式的适用性进行了验证。本研究旨在通过对仿真参数进行标定，基于EDEM 建立一种适用于不同耕作方式的离散元高精度混肥仿真模型。

1 材料与方法

1.1 刀具及整机结构

粉垄机具为浙江省农业机械研究院制造的五轴粉垄机，其主要工作部件为5 把螺旋粉垄刀具。粉垄刀具主体由单螺旋刀片和侧切削刃共同组成，刀轴直径（D）80 mm，螺旋刀片半径（a）80 mm，总高度（Hz）400 mm，螺旋角（α）28°，螺旋线圈数（n）为3。

旋耕机为河北耕耘农业机械制造有限公司生产的170 型旋耕机，用Handyscan 手持式3D 扫描仪获得170 型旋耕机的旋耕刀三维点云模型后进行逆向建模。

使用SolidWorks2018 建立三维模型，对实际机具进行1∶1 原尺寸还原，并将壳体、箱体、悬挂等非入土结构简化处理，保留耕作刀具。结合农艺要求，2种机具的耕作参数如表1所示。

表1 耕作参数Table 1 Tillage parameters

1.2 主要材料参数标定

土壤样本取自浙江省金华市婺城区徐村农业试验基地的试验田，取样方法参照《农业机械试验条件测定方法》[19]。使用与环刀同尺寸的光滑薄壁塑料管获取自然状态下不同土层深度（0—10、10—20、20—30、30—40 cm）的原样土壤样本。化肥样本取自实际耕作时使用的颗粒状复合肥（硫酸钾型复合化肥，火山能公司）。

土壤属性参数包括弹性模量、剪切模量、泊松比、含水率、密度等。使用TMS-PRO 质构仪（美国FTC 公司，检测精度>0.015%）于实验室标定土壤弹性模量，得到力-位移曲线。质构仪的压缩速度为15 mm·min-1，并设置在正压力到达500 N 时停止压缩。每层土壤深度下做4 组重复，共做16 组压缩试验，试验中一组土壤内部掺杂石块，所得曲线出现明显偏差，将其剔除后共获得15 条曲线。参照相关研究进展[13-15,23-25]选取土壤泊松比，通过计算获得土壤剪切模量。

式中，E为土壤弹性模量；G为土壤剪切模量；F为正压力，N；S为土样截面积，m3；dL为压缩距离，m；L为土样原长，m；v为土壤泊松比。

采用环刀法测定土壤容重，所用环刀为100 cm3定体积环刀（南京土壤仪器厂）。土壤含水率通过烘箱干燥法测定，将环刀取得的土样放入烘箱（GZX-9140ME 数显吹风干燥机，上海博迅实业有限公司）105 ℃烘干至恒重，土样做2 份平行测定。基于以下公式计算土壤容重和土壤含水率[21]。

式中，ρs为土壤容重，g·cm-3；m1为环刀的重量，g；m2为土壤和新鲜土壤的重量，g；V为环刀体积，cm3；Ws为样品的土壤含水量g∙kg-1；Ws为土壤含水率；Wf和Wd为土壤的鲜重和干重，g。

使用排油法测定化肥颗粒真实密度。供试材料的其余属性参数参考相关文献[22-24]。

不同材料间接触参数包括土-土间接触参数、土-肥间接触参数、土-钢间接触参数和肥-钢间接触参数；材料间接触参数包括碰撞恢复系数、滑动摩擦系数和滚动摩擦系数。采用下落法标定碰撞恢复系数，以土-肥间接触参数标定为例，将实地取得的原样土壤平铺在平板上，并进行压实以模拟真实土况；利用高速相机，记录从固定高度自由下落的化肥颗粒在与土壤平面发生碰撞后反弹的最高高度，并通过式（5）进行计算。

式中，ε为碰撞恢复系数；∆H为反弹高度，m；H为下落高度，mm。

每种材料间碰撞恢复系数做20 组重复，取平均值。采用自行搭建的摩擦系数测定实验平台（图1）测定摩擦系数，通过改变表面覆盖材料测定不同材料间的摩擦系数[25]。测定滑动摩擦系数时，将一种待测材料覆盖固定于测试平台之上，另一种待测材料利用环刀制成大小均匀的圆饼状置于其上，以保证材料间不发生滚动；测定滚动摩擦系数时，将一种待测材料覆盖固定于测试平台之上，另一种待测材料制成大小均匀的圆球状置于其上，以保证材料间不发生滑动。测定过程中记录所有材料质心的位置，并通过高速相机验证材料间的运动情况，通过数显式旋钮调速器调节转速，测量物体被抛离的距离和抛离高度，计算摩擦系数。每组实验通过多次重复试验以减小误差。

图1 摩擦系数测定平台及测定原理Fig. 1 Coefficient measuring platform and measuring principle

式中，μ为摩擦系数；r为旋转半径，m；h为坠落高度，m；l为抛离距离，m。

1.3 EDEM 仿真模型建立

使用离散元分析软件EDEM2018建立仿真模型，包括刀具模型和导入以及土槽模型。

依据刀具尺寸和耕深，设置土槽尺寸为400 mm×500 mm×1 200 mm。试验土槽每层土壤厚度为10 cm，每层生成15 万个土壤颗粒，共计60 万个[26]。结合试验田土壤特性选取合适的模型： 0—10 cm 土壤松散，呈粉壤土性质，选取Hertz-Mindlin 无滑动模型；10—40 cm 土层，土壤间存在明显黏结现象，呈黏性土性质，选取Hertz-Mindlin with Bonding 黏结模型对土壤颗粒进行约束，根据BPM 黏结模型，法向刚度和切向刚度由式（7）和式（8）计算[27-29]，黏结半径计算公式（9）可由土壤含水率的计算公式推导[18]。

式中，kn为法向刚度，N·m-1；ks为切向刚度，N·m-1；a 为常数，0—10 cm 土层取2/3，10—40 cm土层取1/2；ν1为颗粒1 的泊松比；ν2为颗粒2 的泊松比；E1为颗粒1 的弹性模量，Pa；E2为颗粒2 的弹性模量，Pa；r1为颗粒1 的半径，mm；r2为颗粒2的半径，mm；ω为土壤含水率；m1为土壤颗粒质量，kg；m2为水分质量，kg；ρ1为壤颗粒干密度，kg·m-3；ρ2为水密度，kg·m-3；R为土壤颗粒半径，m；RB为土壤颗粒黏结半径，m。

结合实际施肥，通过设置面颗粒工厂将化肥颗粒随机撒施到土壤表面。为提高仿真精度，同时控制计算时长，结合总计生成的颗粒数量及相关参考文献对颗粒半径进行放尺，放尺后土壤颗粒半径为4 mm，颗粒半径对仿真精度影响不大[30-31]。土、肥颗粒模型生成后，让其自然沉降3 s，生成黏结键[30]。

1.4 仿真模型准确性验证

通过耕作过程牵引力值验证仿真模型参数标定结果的准确性。设置仿真和田间试验时机具工作参数相同，粉垄耕作前进速度为0.3 m·s-1，刀具转速为300 r·min-1，刀具倾角0°；旋耕耕作前进速度0.2 m·s-1，刀具转速400 r·min-1。由EDEM软件输出刀具前进方向合力，即为仿真刀具水平方向牵引力。粉垄机械实际安装过程中，根据不同旋向分别按照120°和180°间隔安装刀具，因此仿真过程按5 把刀具排布角度进行，并将所得前进方向合力求和，得机具整体所受牵引力。

田间试验在浙江省金华市婺城区徐村农业试验基地进行，耕作过程中的牵引力信号通过板环式拉压力传感器（南京天光电气科技有限公司，量程2 t，电压输出范围0～20 mV）测得，通过数据采集系统（DH5902N 坚固型动态信号测试分析系统，江苏东华测试技术股份有限公司）采集电信号并按比例转化为待测力，所测得耕作力被用作验证仿真模型计算结果的准确性。试验所采用的板环式传感器可以测量沿传感器两销孔连线方向的力，当机具入土后保持匀速水平前进时，认为传感器测量处于平衡态，此时可根据机具与传感器形态位置计算牵引力。

1.5 混肥效果优化试验

结合仿真实验对2 种耕作方式下的混肥效果进行分析，以混合的均匀性为混肥效果评价指标，均匀性通过计算各土层化肥颗粒数变异系数（C）评价，C越小混肥更均匀。

式中，xi为第i层颗粒数，p为每层颗粒数平均值，n为层数。

结合实际农艺经验和耕作参数（表1），选取各因素取值范围。通过所建立的仿真模型进行试验，获取2 种耕作方式的分布规律预测模型，通过该预测模型计算获得最优工作参数。试验过程如图2 所示，在仿真模型中设置粉垄刀具运动为旋转下降和旋转前进2 部分，旋耕刀具运动为旋转前进，刀具工作参数根据正交试验表设置，因素的编码值和变异系数如表2 所示，其中，设置粉垄刀具倾斜角度向前进方向倾斜为负，反之则为正。

1.6 混肥效果田间验证

田间验证试验在浙江省金华市婺城区徐村农业试验基地进行。试验田长期进行农业活动，前茬作物为水稻，水稻收割完成后使用铧式犁对试验田进行翻耕，并施用有机肥对试验田理化性质进行调整改造，晾晒1 月后进行田间试验。试验用拖拉机型号为久保田M954KQ。由于实验所施用的化肥易溶于水，难以统计各层的数量分布，因此选用形状、大小、密度与化肥颗粒相似的聚丙烯颗粒以1∶1 的施用比例代替化肥。耕作前采用与撒肥相同的方法将聚丙烯颗粒均匀撒施到试验区域土壤表面，分别使用旋耕机和粉垄机对田地进行翻耕，将颗粒混入土壤。取样时，使用10 cm×10 cm×10 cm 合金钢方框以10 cm 为梯度在0—40 cm 土层的耕作扰动区域内取样，不同层的土壤分别装入不同的取样袋中。分别通过5 和2 mm 孔径筛子筛选，手动挑拣出聚丙烯颗粒，计数并统计颗粒分布情况。

2 结果与分析

2.1 参数标定结果

结合图3 的结果发现，本研究所取土壤样本在压缩过程中的力学特性表现类似于塑性材料，存在弹性阶段（OA）和屈服阶段，且进入屈服阶段以后，试样越压越扁，抗压能力继续增强，因此土壤弹性模量类比线弹性材料，通过A 处获得的轴向载荷计算得到。

图3 土壤力学分析Fig. 3 Soil mechanics analysis

最终得到的材料属性参数和接触参数如表3和表4 所示。可以看出，0—10、10—30 和30—40 cm 土层土壤物性存在明显差异，随着土层深度的增加，土壤被压实，土壤含水率也增加。最终确定黏结模型参数如表5 所示。化肥颗粒的半径通过实际取样，测量得平均值为1.94 mm，因此在EDEM 中设定模型半径为2 mm，根据实际施肥量计算得到试验区域面积内施用的颗粒数目为2 000个。

表3 材料属性参数Table 3 Material property parameter

表4 材料接触参数Table 4 Material contact parameter

表5 黏结模型参数Table 5 Bonding model parameter

2.2 模型准确性分析

粉垄耕作验证结果如图4A 所示，将实际情况下机具受力信号分别放缩10%，作为仿真值的10%误差区间。粉垄耕作过程中，仿真所得刀具受力在入土初期波动较大，且呈现一定的周期性，波动范围为-100～300 N，总体呈上升趋势，这可能是由于机具入土初期，刀具入土深度较浅，与土壤相互作用部分较小，因此整体受力较小，刀具周期性与土壤的相互作用对受力的影响大于接触面积的影响。在入土中期，受力大幅提升，曲线波动不明显，这可能是由于随着耕作深度的增加，刀具与土壤接触面积对受力的影响大于周期性接触的影响。在入土后期，牵引力变化趋势和变化大小相似，但由于入土阶段总体受力较小，实际值的10%误差区间较小，仿真值仅有部分在误差区间内。在机具平稳前进阶段，仿真值基本全部位于10%的误差区间内，牵引力大小在1 100 N 波动，这是因为机具平稳前进时，相同深度土壤性质较为相似，刀具受土壤性质变化影响小，且仿真模型土壤性质参数与实际土壤性质参数的差异不大，因此预测值与真实值拟合程度较好。

图4 耕作牵引力曲线Fig. 4 Tillage traction force curve

旋耕机的牵引力验证结果如图4B所示，在同一时间长度内，仿真结果与实际信号同样呈现较大的波动，波动原因是工作过程中旋耕刀入土数量和入土深度呈周期性波动且旋耕机耕作过程中自身抖动较大。对信号数据进行平滑处理后，发现仿真结果和试验结果呈现相似的变化规律，波动周期大小相似但存在一定差异，差异存在原因可能是实际耕作时刀具转速难以精确控制。由平滑曲线发现，旋耕机实际受力略大于仿真受力，这可能是由于实际工作中旋耕机存在摩擦力产生的机械损耗[19]。

综合上述2 种模型的验证结果，认为该仿真模型具有一定的准确性，可以用于耕作过程的模拟分析。

2.3 混肥效果优化分析

模型方差结果如表6所示，响应曲面如图5所示。粉垄耕作模型显著性检验P=0.000 3，决定系数R2=0.91，失拟项P=0.092 1，回归模型极显著，失拟项不显著，拟合程度较高。由表6 可知，对均匀性的影响，A2、B2、C2极显著，显著性顺序为A2>C2>B2。试验因素与均匀性之间存在二次非线性关系，得到粉垄耕作（CF）的回归模型如下。

图5 标准差响应曲面Fig. 5 Response surface on standard deviation

表6 模型方差分析Table 6 Model analysis of variance

旋耕耕作模型显著性检验P=0.000 5，决定系数R2=0.93，失拟项P=0.063 9，回归模型显著，失拟项不显著，拟合程度较高。由表6 数据可知，对均匀性的影响，B、A2、B2极显著，AB显著，显著性顺序为A2>B2>B>AB。试验因素与均匀性之间存在二次非线性关系，得到旋耕耕作（CX）的回归模型如下。

结合所得回归模型，以获得最小变异系数为目标，优化方程求解得到最优混肥效果下的工作参数为，粉垄转速308 r·min-1，前进速度0.28 m·s-1，倾斜角度1.7°；旋耕转速420 r·min-1，前进速度0.25 m·s-1。根据该组参数，设计田间验证试验。

2.4 混肥效果田间验证结果分析

由EDEM 后处理模块分析化肥颗粒分布情况。分层建立土壤化肥颗粒取样的取样沟模型，统计每层统计模型中化肥颗粒的数量，最后按数量计算化肥颗粒在每层土层中的数量占比。仿真模型如图6所示，最终统计结果如表7所示。

图6 仿真结果Fig. 6 Simulation results

表7 化肥颗粒分布仿真和验证结果Tab. 7 Simulation and verification results of fertilizer particle distribution

从表7 可以看出，在预测准确性方面，旋耕后的化肥分布预测在0—10 和20—30 cm 土层较精准，预测误差分别为4.47%和4.18%，而10—20 cm土层的预测误差为8.65%。旋耕耕作方式下的化肥分布总体呈由上而下逐渐减少趋势，20—30 cm土层化肥颗粒数量占比较少。对其结果进行分析，原因可能是仿真过程中旋耕耕作不能到达20 cm 以下土层深度，但存在少量化肥颗粒沿耕作产生的空隙滑落至耕作区底部；实际取样过程中，耕作机具耕深不稳定，时常大于设定耕深，因此实际取样的颗粒数占比比仿真更高，从而导致颗粒数量波动较大。粉垄耕作过程的预测结果在10—20和20—30 cm土层较为精准，分别为2.32%和0.59%，但0—10 和30—40 cm 土层误差较大，分别为8.03%和6.31%。化肥颗粒分布在数量上大体呈现由上到下逐层减少的规律。对出现的较大误差进行分析，0—10 cm 土层范围内实际测得数量小于仿真数量，原因可能是该土层内化肥颗粒数量较多，田间取样时，少量颗粒滑落至10—20 cm 土层，而10—20 cm 土层实际取样所得颗粒数大于仿真颗粒数也印证了这一点。30—40 cm土层误差较大原因可能为，一方面，深层土层内，化肥颗粒数量较少，取样时数量很少的几粒化肥颗粒数量波动就会导致较大误差的产生；另一方面，同样存在耕深不稳定现象导致实际取样中30—40 cm土层化肥颗粒数量占比更高。

综合表7 数据，混肥效果仿真预测值和测量值最大误差为8.65%，最小误差为0.59%，且仿真与实际测得的化肥颗粒分布呈现相似的规律，理论分析与实际结果非常接近，说明该预测模型具有可靠性。

3 讨论

相较于粉垄耕作方式，旋耕在耕层内能够近似实现化肥的均匀分布，适用于根系浅、根数量多的须根系作物如水稻、小麦、葱等。粉垄耕作的影响区域较大，在每个钻削部件的工作直径范围内都产生影响，其土壤运动范围在深度上较大，因此便于实现地表撒肥深施。粉垄耕作下耕层内的化肥数量近似倒三角分布，相较于其他耕作方式，能够实现不同土壤深度下分布不同化肥颗粒量，更加贴合部分根系较深作物的生长过程，如茄子、青花菜、萝卜等。因此，结合不同作物的生长习性和种植农艺，利用预测模型对特定耕作条件下的化肥分布进行模拟预测，根据预测结果选取最优的耕作方式，可以提高农业施肥效率，改善施肥效果。

目前已有研究利用离散元方法对耕作机具进行仿真优化，陈新予等[32]应用EDEM 软件对旋耕机抛土效果进行了验证；鲍洋清[33]研制了新型仿生深松铲，基于EDEM 软件探究了深松铲的减阻作用和土壤颗粒的运动规律；高亮[34]建立了耕刀切割土壤的离散元模型，并对耕刀形状进行优化，进一步探究了耕刀形状与功耗间的关系。现有研究多从减小耕作阻力和降低能耗角度出发，模拟耕作机具与土壤间相互作用，优化工作部件设计参数，从而实现对机具结构的改进，鲜有从化肥混合效果角度出发，对机具结构进行优化的研究。

本文使用离散元法模拟了粉垄机具和旋耕机具对地表撒施化肥颗粒的混合过程和耕作过程，通过优化计算获得最优混肥效果下的工作参数。对比田间试验结果和仿真试验结果，发现每层化肥颗粒占比的仿真预测值和测量值最大误差为8.65%，最小误差为0.59%，证明了预测模型的准确性。研究结果表明，不同耕作机具作用下的化肥颗粒空间分布大不相同。

综上所述，一方面，本研究结果在旋耕机具设计时可以以化肥颗粒在耕层内混合更均匀为目标对刀具结构等进行优化，在粉垄机具设计时可以以自下而上，由少到多、等差分布为目标进行优化，或结合作物根系分布范围，以化肥分布与根系分布相匹配为目标，对刀具长度、螺旋叶片宽度等进行优化；另一方面，在机具应用时，可以通过重新标定、设置不同的仿真参数，实现不同土壤条件下的工况模拟，通过虚拟试验方法，以达到所需的化肥混合效果为目的，对机具工作时的工作参数进行优化，提高作业质量，从改善化肥颗粒空间分布角度出发，对旋耕、粉垄等耕作机械设计和应用提出新的要求。