何瑞银 段庆飞 陈信信 徐高明 丁启朔
(1.南京农业大学工学院,南京 210031;2.江苏省智能化农业装备实验室,南京 210031;3.江苏大学农业工程学院,镇江 212013)
秸秆还田作为作物秸秆最有效的利用方法之一[1],不仅可以避免秸秆焚烧造成的环境污染问题,而且还田后秸秆的腐解能够释放作物生长所必需的氮磷钾及其他微量元素,可以有效协调土壤的水肥气热关系,增强土壤的固碳能力。此外,秸秆还田对土壤结构、作物根系生长、化肥投入及作物产量等方面也有积极作用[2-4]。
目前机械化秸秆还田作业方式主要分为免耕覆盖还田、犁耕深翻还田、沟埋还田及旋耕还田等,其中免耕覆盖还田可蓄水保墒,改善土壤结构,但作业后秸秆覆盖在地表,腐解速率缓慢;犁耕深翻还田及沟埋还田能够将秸秆埋入较深土层,但会造成秸秆团聚,大幅降低秸秆的腐解速率,且犁耕及沟埋还田作业能耗大,作业后地表平整度较差;相比其他还田方式,旋耕还田作业后秸秆在土壤中的分布更加均匀,有效提升了秸秆的腐解速率,同时旋耕还田还具有能耗低、效率高等优点,故在我国长江中下游的稻麦轮作区,农户普遍采用旋耕的方式进行秸秆还田[5-7]。而经旋耕还田后秸秆在土壤中的空间分布质量会影响秸秆的腐解速率,秸秆在土壤中的空间分布质量越好则秸秆分布越均匀,其与土壤的接触面积越大,有利于土壤中微生物与秸秆的充分接触,进而大大加快了秸秆的腐解速率[5,8]。
近些年来,国内外学者针对还田作业后秸秆在土壤中的空间分布状况进行了相关研究[9-15]。目前主要集中研究了秸秆长度、作业方式及还田机具类型等对秸秆在土壤中空间分布状况的影响,但迄今为止,尚缺乏旋耕机具作业参数对秸秆空间分布状况影响的研究。离散元法适用于模拟不同颗粒在静态或动态条件下发生的接触变形,随着近些年来国内外学者的研究和试验发现,利用离散元仿真不仅可以模拟土壤颗粒间的相互作用,而且可以模拟农业机械的田间作业过程[16-21]。
因此,本研究利用离散元法对不同旋耕作业参数下秸秆的空间分布质量进行研究,进一步结合田间试验验证,通过评价分析不同作业参数下秸秆在土壤中的空间分布质量,以期为旋耕秸秆还田作业质量的快速预测研究提供支持,同时为旋耕机具的作业参数选择提供理论依据。
2020年11月在南京市六合区八百桥村试验田(118°55′E,32°25′N)进行田间旋耕试验,该地区为稻麦轮作区,土壤为壤质粘土,其中壤粒、砂粒、粘粒及有机物的质量分数分别为39.67%、21.20%、38.96%和3.02%。试验开始前,利用五点取样法对试验田块0~40 cm土层土壤的坚实度、含水率、孔隙度等进行测量,具体参数如表1所示。并利用50 cm×50 cm的钢制取样框对秸秆进行称量,测得田间秸秆量为749.62 g/m2。试验开始前对地表进行人工清茬。
表1 试验前0~40 cm土层土壤参数Tab.1 Soil parameters in 0~40 cm soil layer before test
1.2.1试验装备及材料
试验装备采用课题组自主研发的田间原位耕作试验台[22],其主要组成结构如图1所示,该试验台适用于牵引型和驱动型耕作部件的试验机理研究,其刀轴转速在150~400 r/min范围可调,前进速度在0.1~1.0 m/s范围可调。试验选用IT225C型旋耕刀,旋转半径为225 mm。试验材料为秋季水稻秸秆,为排除秸秆长度及秸秆质量等因素的影响,试验对象仅选择秸秆的茎秆部分,并将秸秆按30 mm长度切碎后进行喷漆处理。
图1 田间原位综合耕作试验台Fig.1 Field in situ integrated tillage test rig1.旋耕机具 2.计算机终端 3.伸缩立柱 4.地面导轨 5.牵引电机 6.控制箱 7.立柱升降电机 8.机具升降电机 9.悬架导轨
1.2.2试验方案
试验分别测试不同刀辊转速及前进速度对旋耕还田后秸秆空间分布质量的影响。试验开始前按照749.62 g/m2的秸秆量将处理好的秸秆均铺在大小为500 mm×2 400 mm的试验小区内。由田间原位综合耕作试验台带动旋耕刀组进行固定耕深为100 mm、幅宽为500 mm的旋耕作业。通过全面试验的方法,研究5种刀辊转速(240、260、280、300、320 r/min)及3种前进速度(0.25、0.50、0.75 m/s)旋耕作业后秸秆在土壤中的空间分布质量,每组试验重复3次。
1.2.3取样及测量方法
试验采用如图2所示特制取样框(500 mm×300 mm×50 mm)进行取样,依照水平划分的要求将样框内部均分为4部分区域。旋耕试验后,将取样框放置于作业后的试验小区内,将其压入土壤,将钢制隔板从侧面插入,之后将取样框周围土壤清除,筛选计算各区域内秸秆数量。分别取0~5 cm、5~10 cm、10~15 cm 3层土样。将取样框内秸秆-土壤筛分后计算各区域秸秆数量。
图2 特制取样框Fig.2 Special sampling frame
仿真旋耕装备选用江苏地区常规旋耕机,刀具类型为IT225C型弯刀。利用Pro/E 5.0对选用的旋耕装备的触土部件按照1∶1的比例进行三维建模,并以.igs格式将旋耕模型导入EDEM 2018仿真软件中,同时构建尺寸(长×宽×高)为5 000 mm×1 000 mm×400 mm的仿真虚拟土槽(图3)。通过更改运动参数来模拟5种刀辊转速(240、260、280、300、320 r/min)、3种前进速度(0.25、0.50、0.75 m/s)下的旋耕还田作业过程。
图3 旋耕秸秆还田离散元仿真模型Fig.3 Discrete element simulation model of rotary cultivation straw returning to field
2.2.1土壤-秸秆颗粒模型
本文选用球形颗粒模拟土壤。在离散元仿真中,颗粒模型尺寸越小,则仿真速度越慢,耗时越长,因此仿真中土壤颗粒尺寸一般远大于其实际尺寸,本文仿真土壤颗粒半径选为8 mm[23-25],同时根据表1测量的田间实际土壤类型与状态,用不同参数土壤模型来模拟实际耕层土壤(0~150 mm)与犁底层土壤(150~400 mm)。由于水稻秸秆具有易弯折性、中空性和含水率动态变化性,所以秸秆仿真一直是一项具有难度的研究课题。CHANDIO[26]在研究中发现,使用简化刚性秸秆模型进行仿真试验其结果与田间试验基本一致,故为了减少仿真工作时间,同时为了避免因秸秆弯曲或折断对田间试验结果造成的不利影响[27],本文采用秸秆刚性模型,用9个直径为6 mm、球心间距为3 mm组成的总长度为30 mm的线性模型来模拟实际水稻秸秆。
2.2.2土壤接触模型
选择合适的颗粒间接触模型是离散元仿真成功的首要条件,接触模型表示的是颗粒固体在准静止情况下的弹塑性分析结果。颗粒间所受力及力矩大小由接触模型的分析计算所决定。因此需要建立不同的接触模型来模拟不同的仿真对象[28]。本文田间试验土壤类型为壤质粘土,并且土壤间存在粘附力,故选用适用于粘性土壤的Hertz-Mindlin with Bonding模型,相比其他模型,该模型对颗粒间粘结作用及破碎程度的模拟较为可靠[29]。采用该模型仿真时,土壤颗粒间存在的抵抗法向及切向运动的作用力使得颗粒间产生粘结作用,而随着时间步长的增加,颗粒所受外界作用力逐渐增加,当外界作用力超过临界值后,粘结作用被破坏,此后土壤颗粒间不再受粘结作用的影响。
Hertz-Mindlin with Bonding模型的细观参数主要包括:颗粒粘结半径、粘结刚度和粘结临界应力。其中颗粒粘结半径可以反映湿颗粒含水率的大小,其在颗粒半径一定的情况下可以通过材料密度、含水率计算得到[30],本文耕层土壤与犁底层土壤颗粒粘结半径分别为9.5 mm和9.15 mm。而据已有研究可知,粘结刚度对颗粒的运动并不会产生明显影响,因此可取粘结刚度为固定值5×107N/m3[24]。临界应力决定了颗粒间的粘结强度,可以反映土壤的破碎程度以及耕作阻力的大小。本文引用丁启朔等[29]的方法,采用单轴土壤压缩试验获得田间不同土层土壤的最大应力级别,而后利用秸秆旋耕还田离散元仿真模型,进行不同临界应力下的旋耕作业仿真试验,将不同临界应力对应的耕作阻力与田间实测的耕作阻力进行误差分析对比,确定耕层与犁底层土壤颗粒模型的临界应力取值分别为3×105Pa与5×105Pa。
2.2.3颗粒模型参数
仿真所需要的模型参数主要分为材料参数与接触参数。其中材料参数包括土壤、秸秆、旋耕刀的泊松比、剪切模量和密度,而接触参数是指材料间的恢复系数、摩擦因数。材料参数可以通过查阅文献[21,31]以及实地测量获得,接触参数则主要通过引用文献[32-34]、实地测量及仿真标定获得,其中土壤与土壤间静、动摩擦因数通过仿真标定方法,利用离散元仿真进行休止角试验及贯入试验获得,土壤-旋耕刀静摩擦因数利用直剪试验测量获取[35];土壤-旋耕刀、秸秆-旋耕刀滚动摩擦因数通过斜板试验测量获取[36],其他接触参数则通过查阅相关文献获得。本文仿真试验所有参数如表2所示。
表2 仿真模型参数Tab.2 Simulation model parameters
仿真作业前0~15.0 s先在虚拟土槽深度150~400 mm 空间内生成497 359个犁底层土壤颗粒,犁底层土壤颗粒的土壤孔隙度约为40%,15.0~17.0 s在土槽深度0~150 mm区间内生成254 993个耕层土壤颗粒,耕层土壤颗粒的土壤孔隙度约为50%,17.0~23.1 s时间段内在耕层土壤表层生成均匀平铺的25 383个秸秆颗粒,秸秆颗粒生成区域大小为500 mm×2 400 mm,仿真秸秆量与田间实测秸秆量相同,控制在749.62 g/m2。后将三维旋耕装备模型导入,控制耕作深度为100 mm,耕后耕层深度约150 mm[11],23.1 s后改变旋耕装备刀辊转速及前进速度并设定仿真时间、步长、网格大小等值后开始旋耕还田仿真作业。仿真作业过程中,0~23.1 s为颗粒生成及稳定时间,23.1 s后为旋耕装备运动与颗粒沉降阶段。
仿真结束后,利用后处理模块中的Grid Bin Group设置相应计算区域测算耕后土壤中不同区域的秸秆数量,于土槽中心点设置整体计算区域大小为500 mm×2 400 mm,并以125 mm为单位长度沿旋耕装备幅宽方向进行如图4a所示的纵向划分,以300 mm为单位长度沿旋耕装备前进方向进行如图4b所示的横向划分,共计形成32个计算区域。
通过设置后处理模块Grid Bin Group中Z轴的坐标值,来对应耕后0~5 cm、5~10 cm、10~15 cm的垂直分层处理。并且通过设置显示秸秆颗粒数量来统计不同土层及区域内的秸秆数量。为便于区分,计算区域秸秆颜色为橙色,其他区域秸秆颜色为绿色,显示效果及各区域秸秆数量统计如图5所示。
图4 水平区域划分示意图Fig.4 Schematics of horizontal area division
将5种刀辊转速与3种前进速度下旋耕作业后各区域秸秆数量的仿真试验值与田间试验值进行测量统计,并以各区域的秸秆占比变异系数为评价指标评价秸秆在土壤中的垂直分布及水平分布质量,变异系数表示秸秆空间分布质量,数值越小,质量越优,反之则越差。
本文以5 cm为分割尺度对耕后0~15 cm土层进行垂直分层处理,结合0~5 cm、5~10 cm和10~15 cm土层的秸秆占比及其变异系数评价秸秆的垂直分布质量。
图5 秸秆显示效果及各层秸秆数量统计示意图Fig.5 Straw display effects and statistical schematics of number of straw in each layer
3.1.1刀辊转速对秸秆垂直分布的影响
通过测量统计垂直分层处理后的各层秸秆数量,可以确定各层秸秆占比。5种刀辊转速下旋耕还田作业后各层秸秆占比的仿真试验值与田间试验值如表3所示。由表3可知,随着转速的增加,仿真与试验值的变化趋势基本一致,0~5 cm土层的秸秆占比基本呈现递减的趋势,转速从240 r/min增加到320 r/min,0~5 cm土层的仿真试验秸秆占比减少7.3个百分点,田间试验秸秆占比减少4.3个百分点;5~10 cm土层的秸秆占比基本呈现递增的趋势,转速从240 r/min增加到320 r/min,仿真试验及田间试验秸秆占比分别增加6.1个百分点和6个百分点;而随着转速的增加,10~15 cm土层仿真试验秸秆占比增加1.2个百分点,田间试验秸秆占比减少1.7个百分点,秸秆占比变化量较小。总体趋势表明,转速的增加可以显著增加5~10 cm土层的秸秆数量,同时减少0~5 cm土层的秸秆数量。这可能是因为转速的增加会导致秸秆与土壤的切向加速度增加,使得秸秆沿垂直方向的位移增加,进而改变了各层秸秆的占比。并且由表3可知,秸秆主要集中分布在5~10 cm,这是因为正转旋耕工作过程中是将秸秆以挤压的形式埋入土壤的,这种形式会造成秸秆在垂直方向上的分布不均匀,大多数秸秆集中在5~10 cm的中层土壤中,这与陈青春等[11]得出的结论一致。
图6为5种刀辊转速作业后各层秸秆占比的变异系数,图中不同字母表示同一类型不同转速间的变异系数差异显著。由图6可知,随着刀辊转速的增加,旋耕作业后仿真与田间试验各层秸秆占比的变异系数基本呈现递增的趋势。其中仿真及田间试验在转速240 r/min时变异系数最小,分别为60.09%和80.65%,在转速320 r/min下变异系数最大,分别为74.11%和93.11%。这表明,无论是仿真还是田间试验,刀辊转速的增加总会使得秸秆在土壤中的垂直分布质量下降。其主要原因是由于增加转速会造成秸秆各层占比变化,导致各层秸秆占比相差较大,从而使得变异系数增大,秸秆在土壤中的垂直分布质量下降。5种刀辊转速下仿真值及试验值的差值分别为20.56%、24.12%、23.51%、23.46%和19%,平均值为22.13%。
表3 5种刀辊转速作业后各层秸秆占比Tab.3 Proportion of straw in each layer after operationat five kinds of knife roller speeds %
图6 5种刀辊转速作业后各层秸秆占比的变异系数Fig.6 Variation coefficient of proportion of straw in each layer after five kinds of knife roll speed operation
3.1.2前进速度对秸秆垂直分布的影响
3种前进速度(0.25、0.50、0.75 m/s)旋耕作业后仿真与田间试验的各层秸秆占比如表4所示。仿真与田间试验对比数据具有一定差别,但其规律趋势基本一致,即随着前进速度的不断增加,0~5 cm土层的秸秆占比基本呈现递增的趋势,前进速度从0.25 m/s增加到0.75 m/s,仿真与田间试验秸秆占比分别增加11.9个百分点和12.3个百分点;5~10 cm土层的秸秆占比呈现减少的趋势,随着前进速度的增加,仿真与田间试验秸秆占比分别减少4.2个百分点和1.3个百分点,秸秆占比变化不显著;10~15 cm土层的秸秆占比呈现递减趋势,随着前进速度的增加,仿真与田间试验秸秆占比分别减少7.7个百分点和11.0个百分点。这表明,前进速度的增加可以显著增加0~5 cm土层的秸秆占比,同时减少10~15 cm土层的秸秆占比,这是因为旋耕装备前进速度的增加会降低旋耕碎土性能,导致被挤压至下层土壤的秸秆量变少,更多的秸秆停留在上层土壤中。
表4 3种前进速度作业后各层秸秆占比Tab.4 Proportion of straw in each layer after operation at three forward speeds %
3种前进速度下旋耕还田后各层秸秆占比的变异系数如图7所示。由图7可知,3种前进速度旋耕作业后的仿真与田间试验各层秸秆占比的变异系数均呈先减小后增大的趋势,其中前进速度0.25 m/s条件下变异系数最大,仿真值与试验值分别为73.99%和88.89%,而前进速度0.50 m/s下变异系数最小,仿真值与试验值分别为61.00%和79.90%,前进速度增加至0.75 m/s后,仿真与试验各层秸秆占比的变异系数分别为66.74%和85.09%。这说明,随着前进速度的增加,秸秆在土壤中的垂直分布质量呈现先升高后下降的趋势,这一现象可能是因为前进速度的增加导致旋耕作业碎土效果降低,从而影响秸秆入土效果,使得各层秸秆比例发生了变化,从而影响了秸秆在垂直方向上的空间分布质量。3种前进速度下仿真值及试验值的差值分别为14.90%、18.90%和18.35%,平均值为17.38%。
图7 3种前进速度作业后各层秸秆占比的变异系数Fig.7 Variation coefficient of proportion of straw in each layer after operation at three forward speeds
将仿真值与试验值进行对比,得出不同转速下仿真与试验各层秸秆占比变异系数间的平均差值为22.13%,不同前进速度下仿真与试验各层秸秆占比变异系数间的平均差值为17.38%,这是由于EDEM离散元仿真的颗粒参数及作业参数值只能接近真实值,但无法与真实值保持一致,因此这些仿真参数会造成最终试验结果与田间实际试验结果有一定偏差,但产生的相对误差在可接受范围内,并且仿真与试验的数据变化趋势基本相符。
为了探究转速及前进速度对秸秆在土壤中水平分布质量的影响,本文将耕后区域按纵向与横向进行空间分割处理。
3.2.1刀辊转速对秸秆水平分布的影响
5种刀辊转速下各区域秸秆占比变异系数如图8所示。由图8可知,随着刀辊转速的变化,仿真与田间试验水平纵向划分及横向划分后的各区域秸秆占比变异系数没有明显变化规律,但刀辊转速对横向划分区域变异系数的影响大于其对纵向划分区域变异系数的影响,这是因为方会敏等[38]研究发现,随着转速的增加,秸秆沿机具前进方向的水平位移总是会大于同转速下的侧向位移。因此,刀辊转速的改变对秸秆在土壤中的水平分布质量并无规律性影响。5种刀辊转速下纵向划分后仿真值及试验值的差值分别为13.12%、13.91%、9.88%、10.38%和13.85%,平均值为12.23%;横向划分后仿真值及试验值的差值分别为11.47%、1.66%、15.71%、9.40%和15.86%,平均值为10.82%。
图8 5种刀辊转速作业后各区域秸秆占比变异系数Fig.8 Variation coefficient of proportion of straw in each area after five kinds of knife roll speed operation
图9 3种前进速度作业后各区域秸秆占比变异系数Fig.9 Variation coefficient of proportion of straw in each area after three forward speeds
3.2.2前进速度对秸秆水平分布的影响
图9为3种前进速度作业后各区域秸秆占比变异系数,由图9可知,随着前进速度的增加,纵向划分的各区域秸秆占比变异系数逐渐减小,仿真值与试验值分别由0.25 m/s时的17.26%和30.63%减小至0.75 m/s时的11.36%和20.12%,而与此同时横向划分的各区域秸秆占比变异系数整体增大,仿真值与试验值分别由0.25 m/s时的2.57%和11.37%增大至0.75 m/s时的3.79%和14.32%,但就整体而言,前进速度的增加可以适当优化秸秆在土壤中的水平分布质量。3种前进速度下纵向划分后仿真值及试验值的差值分别为13.37%、10.02%和8.76%,平均值为10.72%;横向划分后仿真值与试验值间的差值分别为11.37%、9.56%和14.32%,平均值为8.72%。
通过对仿真及田间试验数据分析发现,秸秆占比变异系数可以评价秸秆在土壤中的水平分布质量,并且不同转速下仿真与试验各区域秸秆占比变异系数间的差值平均为12.23%(纵向)和10.82%(横向);不同前进速度下仿真与试验各区域秸秆占比变异系数间的差值平均为10.72%和8.72%,产生的相对误差均在可接受范围内。
旋耕作业后秸秆的空间分布质量会对土壤的养分、结构以及作物的生长状况产生重要影响,而旋耕的作业方式、参数和机具类型以及秸秆的长度和类型等因素均会影响秸秆空间分布质量。目前现有的研究多集中在作业方式、机具类型、秸秆长度等方面,尚缺乏旋耕机具作业参数对秸秆空间分布质量影响的研究。本文利用离散元仿真结合田间试验的方法,研究评价了不同典型作业参数下秸秆的空间分布质量,为旋耕机械的作业参数优化和选择提供重要理论依据。
目前在利用离散元仿真模拟旋耕还田作业过程的研究中,大多采用刚性模型来模拟秸秆颗粒,虽然有部分学者研究了柔性秸秆颗粒的离散元仿真模型[37-39],但柔性秸秆颗粒模型建立较为困难且柔性秸秆弯曲断裂后其数量及分布状态都会发生相应变化,难以准确定位秸秆的空间位置,因此,本文选用刚性秸秆颗粒模型,并将仿真与田间试验秸秆长度定为30 mm,可以避免秸秆因田间旋耕作业导致的断裂弯曲现象,提高作业后秸秆位置的准确性,便于统计耕后秸秆数量。此外本文田间试验采用的取样方法为田间原位取样统计,旋耕试验后在田间现场统计秸秆在各区域分布数量,无需将样本带回实验室统计测量,避免了土壤样本因运输搬运过程中产生振动而引发的秸秆位置偏移。
国内外已有学者利用离散元仿真对还田作业后秸秆在土壤中的分布状况进行了相应研究,MARI等[40]利用离散元仿真建立了秸秆-土壤-圆盘犁模型,研究分析了圆盘犁作业后秸秆的位移变化。方会敏等[34]利用建立的秸秆-土壤-旋耕刀仿真模型,模拟研究了旋耕作业后秸秆的位移大小变化。但以上只是探究了耕后秸秆在水平方向和侧向的运动位移情况,没有从秸秆空间分布均匀性层面进行秸秆还田质量的评价。本文利用离散元法建立的旋耕秸秆还田仿真模型及区域划分方法,还可以模拟研究不同刀具类型、不同耕深等其他因素对旋耕还田秸秆空间分布的影响,能够为旋耕秸秆还田作业质量的快速预测研究提供支持。
(1)利用离散元法建立了秸秆旋耕还田仿真模型,并通过更改刀辊转速及前进速度参数值来模拟不同关键作业参数下的秸秆旋耕还田作业。利用田间原位综合耕作试验台进行了与仿真试验相对应的田间试验,将仿真及田间试验区域进行垂直分层划分及水平横纵向划分,计算各区域秸秆数量并将结果进行分析对比。
(2)在秸秆垂直分层处理中,仿真与田间试验结果表明刀辊转速的增加显著增加5~10 cm的秸秆占比,同时减少0~5 cm的秸秆占比,并且使得各层秸秆占比的变异系数总体呈现递增的趋势,总的来说,刀辊转速的增加降低了秸秆在土壤中的垂直分布质量;而前进速度的增加显著增加0~5 cm土层的秸秆占比,同时减少10~15 cm土层的秸秆占比,各层秸秆占比的变异系数呈先减小后增大的趋势,即前进速度的增加使得秸秆在土壤中的垂直分布质量呈现先提升后下降的趋势。在秸秆水平划分处理中,仿真与田间试验表明前进速度的增加可以有效减小纵向划分区域内的秸秆占比变异系数,同时前进速度的增加会小幅度增加横向划分区域内秸秆占比的变异系数,总体而言,前进速度的增加可以优化秸秆在土壤中的水平分布质量,而转速的增加对秸秆在土壤中的水平分布质量并没有规律性影响。
(3)通过将仿真值与田间试验值进行对比,发现垂直分布与水平分布的秸秆占比变异系数差值平均最大分别为22.13%和12.23%,仿真与试验误差在可接受范围内。研究发现,离散元仿真可以较好地模拟不同作业参数旋耕作业后秸秆的空间分布状态,可为旋耕秸秆还田作业质量的快速预测研究提供支持,也有助于为旋耕机械的作业参数选择提供理论依据。