湿粘水稻土深松过程离散元分析

丁启朔 任 骏 BELAL Eisa Adam,2 赵吉坤 葛双洋 李 杨

(1.南京农业大学工学院， 南京 210031； 2.扎林盖大学农业工程系， 扎林盖 06)

湿粘水稻土深松过程离散元分析

丁启朔1任 骏1BELAL Eisa Adam1,2赵吉坤1葛双洋1李 杨1

(1.南京农业大学工学院， 南京 210031； 2.扎林盖大学农业工程系， 扎林盖 06)

湿粘水稻土的深耕松对作业机具设计与作业参数的要求不同于旱作制，为探究水稻土条件下土壤深松扰动过程及其内在相互作用机理，本文结合田间实测土壤物理参数，借助离散元EDEM软件，建立适用于粘性水稻土的深松耕作离散元模型。利用粘性水稻土模型对水稻土的机械深松耕作过程进行离散元模拟仿真，并结合田间试验结果对机具耕作阻力、土壤宏观扰动进行对比验证，进一步从仿真的角度揭示深松土壤扰动的微观过程、土壤失效破碎的微观机理。结果表明，该模型下深松耕作阻力平均误差为6.63%，土壤扰动的起垄宽度平均误差为4.39%，起垄高度平均误差为19.22%；DEM仿真对土壤的微观扰动过程分析进一步论证了宏观试验测试结果及假说的正确，并且能够从土颗粒接触力学层面表达土壤扰动的边界生成、土体内部破碎等过程；结合DEM仿真角度提出的反映土体破碎程度指标——断裂系数，与传统试验的指标参数碎土系数对比两者误差为3.46%，该指标更有利于对土壤破碎过程微观机理的表达。

深松； 土壤扰动； 离散元； 水稻土

引言

稻麦轮作制是长江中下游地区的主导高产种植模式，然而该模式的机械化耕作问题日趋严重，因稻麦轮作周年干湿交替、稻季淹水以及长期浅旋代深耕翻造成土壤硬化、孔隙率低、耕层浅薄等问题日趋严重[1-5]，严重限制了稻麦作物的高产创建和品质提升，水稻土的合理机械化耕作已成为制约南方稻作制发展的瓶颈问题。为改善水稻土的结构性，近年来人们开始探索将北方已成熟应用的深耕深松技术应用于南方耕作的可能性[6-9]。不过相关研究表明现有的深松技术尚难于满足不同土壤条件下的深松要求[10]，因此对于特定土壤类型的深松性能差异性研究就显得尤为重要。黄玉祥等[11-12]在西北地区塿土条件进行了系统研究，分析了该土壤条件的深松扰动行为且提供了作业效果的方法；郑侃等[13]针对华北平原壤土区的特定土壤类型设计了一种新型的折线刃深松铲。鉴于水稻土土壤结构的复杂性和特殊性，适用于北方旱作区的深松机具及理论并不完全适合南方稻麦轮作区域[9]。因此须对南方水稻土条件下深松过程深入研究。

现有研究指出，土壤耕作(尤其是深耕松)的机理非常复杂，机械扰动过程受土壤空间变异、耕作机械动力学、土壤破碎等众多复杂因素影响，仅仅采用传统的试验方法难以深刻理解深松土壤的扰动行为[11]，因此需要进一步从内在机制上探讨与机械耕作密切相关的土壤物理响应及其表现[14]。但鉴于宏观测试方法的局限性，其定量化测试与描述并不容易实现，因此，近年来更强调借助离散元等数值模拟手段解析土壤的机械耕作过程[11,13,15-16]。实践表明，数值模拟方法的确为深刻阐释深耕松等机械耕作过程提供了有效的解决途径[10]，郑侃等[13]、李博等[17]、邓佳玉等[18]对不同工况下的深松铲耕作力学行为进行了研究，显示应用离散元法分析深松铲耕作过程的力学行为及效应的可行性和有效性，计算仿真已经成为揭示田间试验耕作力学过程的优势技术手段。

本文针对水稻土耕作力学行为的特殊性及其定量表达的困难性，以高粘结性为特征的水稻土离散元参数模型对深松土壤扰动过程进行机理分析，进一步结合田间试验验证，探讨湿粘性水稻土机械化深松的土壤失效机理及深松机具的力学行为与土壤结构效应，以期从微观扰动过程及破碎机理的仿真角度实现与田间宏观测试结果的配合，为南方耕作机具的研究提供借鉴及方法学的支撑。

1 田间试验

1.1 试验地点

试验地点位于南京市六合区八百桥试验田(118°59′E、31°98′N)，地势平坦，土壤类型为壤质粘土，该地区常年稻麦轮作。土壤pH值为7.6，土壤砂粒、壤粒、黏粒、有机物质量分数分别为24.06%、40.39%、35.55%和2.27%，液限为41.94%，塑限为25.18%。0～40 cm土层土壤容重、含水率、总孔隙度及硬度参数见表1。试验前进行人工地表清茬。

表1 0～40 cm土层基础物理参数Tab.1 Basic physical properties in 0～40 cm soil layer

1.2 试验方法

为保证田间试验的可靠性，本文在土壤耕作原位综合试验台[19](图1)上进行原位深松耕作试验。试验台内部台车上安装2个LKL-101/1T型拉压力传感器用于测量深松铲的水平耕作阻力。利用Excel 软件进行数据分析，获取耕作阻力波动曲线(图2a),对信号进行滤波处理后(图2b)截取波动稳定段并求其平均值，作为相应耕作深度下的水平耕作阻力。

图1 土壤耕作原位综合测试平台Fig.1 Schematic view of in-situ experimental platform1.牵引电动机 2.控制柜 3.链条 4.信号拖链 5.台车机构6.升降电动机 7.台车导轨 8.电源拖链 9.信号终端计算机10.移动导轨

南方水稻土的适宜深耕深松深度尚难明确，本文进行5个耕作深度(100、150、200、250、300 mm)的深松深耕试验，一方面可以完整反映耕深变化的耕作效应，另一方面则可观察水稻土条件深松铲不同部位的触土性能差异，耕作速度为0.1 m/s。一次深松作业完成之后，随机选取2 m耕作区域作为测量区域。沿着裂纹方向将测试区域内的不同尺度的土垡缓慢取出放于海绵垫上，自然风干后进行筛分处理，土垡尺度过大(大于64 mm)且难于筛分时采用手工测量。小尺度区间(32～64 mm、16～32 mm、8～16 mm、0～8 mm)的土壤破碎体使用筛分分级，分级尺度区间的上、下限算术平均值作为该尺度区间的平均土块径[20]，计算结果依次为48、24、12 mm，小于8 mm土壤粒径默认为4 mm。上述试验操作及测试方法见文献[9]，所有试验均重复3次。

图2 耕作阻力波动曲线Fig.2 Wave curves of tillage draught

2 EDEM仿真分析

2.1 深松铲建模

图3 试验用深松铲Fig.3 Subsoiler used in field test and its structural parameters

本试验选用凿型深松铲进行深松耕作，深松铲实物图如图3a所示，铲柄为圆弧形(仿制JB/T 9788—1999《深松铲和深松铲柄》所规定的圆弧形深松铲柄)，其结构简图如图3b所示。其中铲柄高度600 mm，铲柄宽度60 mm，铲尖长度160 mm，铲柄曲率半径300 mm，入土角为22°。为保证仿真分析的准确性，对试验用深松铲(图3)按照1∶1的比例，采用Pro/E建立深松铲的三维结构模型，并保存为.igs格式导入EDEM软件(图4)。

图4 EDEM仿真模型Fig.4 Simulation model of EDEM

2.2 土壤模型建立

本研究选用球形颗粒模拟田间土壤颗粒，离散元仿真中颗粒尺寸减小会导致仿真耗时呈几何级数增加，故仿真中的颗粒尺寸受计算时间和存储空间限制，总是比真实土壤尺寸大[15-16]，本文选定土壤颗粒半径为8 mm，处于可接受范围内。结合表1田间实测土壤状态，建立两种不同物理参数(耕层土壤、底层土壤)的土壤参数模型模拟田间水稻土分层结构，其中0～150 mm为耕层土壤，150～400 mm为底层土壤。

2.2.1 土壤接触模型

接触模型是离散元法的重要基础，其实质是准静态下颗粒固体的接触力学弹塑性分析结果[21]。接触模型的分析计算直接决定了颗粒受力和力矩的大小，对不同的仿真对象须建立不同的接触模型，用以提高仿真结果的准确性。试验地区为壤质粘土，颗粒表面粘附力较大(粘粒35.5%)，土颗粒接触模型选定为颗粒之间存在粘结作用的Hertz-Mindlin with Bonding模型，该模型可以很好地表达粘性土壤的破碎过程[15-16,22-23]。

该模型接触情况下，粘结作用既能传递作用力又能传递力矩，因此粘结可看作一系列具有恒定法向和切向刚度的弹性弹簧，这些弹簧均匀分布在两颗粒之间的球形区域内。粘结作用产生后，外界给予土颗粒的力Fn、Fτ和力矩Mn、Mτ随着时间步通过下面公式[24-25]从零开始增加。

δFn=-vnSnAδt

(1)

δFτ=-vτSτAδt

(2)

δMn=-ωnSnJδt

(3)

(4)

其中

式中A——接触区域面积Rb——粘结半径J——截面极惯性矩Sn、Sτ——粘结颗粒法向、切向刚度

vn、vτ——颗粒运动速度的法向、切向分量

ωn、ωτ——颗粒角速度的法向、切向分量

δt——时间步长

当外界作用力超过某个定义值时，粘结就被破坏。定义法向和切向临界应力为

(5)

当外界的作用力破坏了土壤颗粒之间的粘结后，土颗粒之间的相互作用将不再受到粘结作用的影响。

2.2.2 土壤模型参数

在离散元EDEM软件中，模型参数主要分为材料参数和接触参数，材料参数包括剪切模量、密度、泊松比，接触参数主要包括材料间的摩擦因数、恢复系数。材料参数和接触参数可以通过实地、实验测量及查阅已有本地区相同土壤类型的相关文献得到，其中接触参数通过查阅文献[13,16,26-27]可得：土壤-土壤及土壤-深松铲恢复系数均为0.6、土壤-土壤静摩擦因数0.5、土壤-深松铲静摩擦因数0.6、耕层土壤间动摩擦因数0.6、耕层土壤-深松铲动摩擦因数0.35、底层土壤间动摩擦因数0.24、底层土壤-深松铲动摩擦因数0.13。本文选用的针对粘性土的Hertz-Mindlin with Bonding土颗粒接触模型本身有5个微观参数，分别是：粘结法向刚度Sn、粘结切向刚度Sτ、粘结法向临界应力σmax、粘结切向临界应力τmax及颗粒粘结半径Rb。其中，前面4个参数反映颗粒之间的粘性大小，颗粒粘结半径则可以反映湿颗粒含水率的高低，这5个参数综合可以体现湿粘水稻土的高含水率、高粘性的特点。粘结半径Rb可通过公式计算出，在颗粒半径一定的情况下，湿颗粒的粘结半径可根据材料的密度、含水率计算得到[27]。其他微观参数由于不便通过试验测出，故本文选择虚拟仿真标定的方法得出。参照文献[28-29]的方法，为了减少标定参数的个数，可以取Sn=Sτ，σmax=τmax，已有的研究表明该模型下颗粒行为对粘结刚度参数不敏感，因而可取粘结刚度为5×107N/m3[30]。临界应力的取值关系到粘结强度，进而决定耕作过程中土壤的破碎程度及耕作阻力，因此本文在其他参数确定的情况下，以田间实测耕作阻力为标准校核临界应力的取值，建立适合本试验过程的深松铲-土壤离散元耕作模型。具体方法是在EDEM软件中建立如图4所示的模型，以颗粒接触模型中临界应力为变量(其他参数如前文所述)进行深松耕作的仿真试验，得出不同临界应力下土壤对深松铲的水平耕作阻力，然后参照田间试验相同工况下的水平耕作阻力进行误差分析，确定适宜的参数取值。

考虑到土壤模型的分层结构，分别以150 mm、300 mm为耕作深度校核耕层土壤和底层土壤的粘结强度。应先校核耕层土壤的粘结强度，得到确切的耕层土壤颗粒间临界应力取值后，再以此为基础校核底层土壤颗粒粘结的临界应力值。田间试验下测得150 mm、300 mm耕深下的耕作阻力分别为1.47 kN和5.21 kN。代表颗粒粘结强度的临界应力可以以宏观土壤剪切、压缩破坏强度为参考[27]，田间取样后做单轴压缩试验得耕层、底层土壤最大法向应力为105Pa级别，以此为基础设定一系列代表粘结强度的临界应力，进行仿真试验得到相应的耕作阻力并与田间试验真实阻力进行比较。仿真标定试验结果如表2所示，最终确定耕层、底层土壤颗粒粘结强度的临界应力为3×105Pa和5×105Pa，由于模型土壤颗粒与真实土壤颗粒粒径的差异，模型临界应力取值比真实土壤破坏强度稍大。最终仿真设计的所有参数如表3所示。

表2 土壤粘结临界应力校核参数Tab.2 Calibration process and effects of particle critical stiffness on tillage draught

表3 离散元模型基本参数Tab.3 Basic parameters of discrete element model

2.3 EDEM耕作模型

本文采用离散元EDEM软件进行整个深松耕作过程的仿真。为满足深松作业的要求，在模型中建立长×宽×高为1 000 mm×800 mm×400 mm的虚拟土槽[11]，深松铲模型导入后位于土体的一端，结合田间试验深松铲耕作速度为0.1 m/s，耕作深度可在模型建立时自主设定。结合田间实测土壤状态，土体模型分为上下两层。上层耕作层深度区间为0～150 mm，生成土颗粒按随机排列，土壤孔隙度约为50%；下层犁底层深度区间为150～400 mm，按简单立方晶格排列生成土颗粒，土壤孔隙度约为40%。仿真模型中共生成土壤颗粒79 824颗，仿真中涉及到所有参数如表3所示。对应表3在EDEM前处理器模块依次进行接触力学模型、土壤颗粒模型、几何模型和颗粒工厂等的设置。仿真试验为单因素试验，通过改变耕作深度来分析不同情况下对土壤扰动情况的影响。在EDEM求解器模块对仿真时步、仿真时间、数据保存间隔时间及网格大小等进行设置。仿真开始时生成土壤颗粒，待颗粒沉降稳定后深松铲开始运动，仿真时间总共为12 s，0～2 s为土壤颗粒生成及沉降稳定阶段，2～12 s为深松铲运动阶段。仿真结束后在EDEM后处理工具模块进行仿真结果数据的导出和分析。

3 结果与分析

3.1 耕作阻力与土壤扰动对比分析

图5a反映了不同耕作深度下田间实测水平耕作阻力与仿真水平耕作阻力的变化趋势。其中耕作深度为150 mm时误差最小，为0.68%，耕作深度200 mm时误差最大，为11.05%，5个耕作深度下耕作阻力的平均误差为6.63%。由图5a可知耕作阻力随着耕作深度的增加而增大，当耕作深度在200 mm以内时，耕作阻力平缓增加；当耕作深度大于200 mm时，耕作阻力急剧增加。

图5 不同耕作深度下实测与仿真结果对比Fig.5 Comparisons of simulated and measured results under different tillage depths

图5b反映了不同耕作深度下田间试验土壤宏观扰动状态与仿真下相对应值的对比。田间试验条件下起垄宽度随着耕作深度先增大后减小，在耕作深度200 mm时达到最大值，为422 mm，耕作深度100 mm时达到最小值，为264 mm；起垄高度随着耕作深度变化趋势相对平缓，总体也呈现出先增大后减小，在耕作深度200 mm时达到最大值(163 mm)，耕作深度300 mm时为最小值(103 mm)。仿真环境下起垄高度变化趋势与田间试验结果高度吻合，各耕作深度下起垄宽度平均误差为4.39%；起垄宽度变化总体趋势一致，但在耕作深度150 mm及300 mm时误差相对较大，分别为32.03%和39.81%，总体显示5个耕作深度的起垄高度平均误差为19.22%。

图6是仿真环境下不同耕作深度的深松铲垂直耕作阻力变化情况。其变化规律和水平耕作阻力变化一致，当耕作深度在200 mm以内时，增幅平缓；当耕作深度大于200 mm时，垂直耕作阻力急剧陡增。

图6 仿真环境下不同耕作深度的垂直耕作阻力 变化曲线Fig.6 Variation curve of vertical tillage draught in simulation under different tillage depths

本文以田间实测耕作阻力为校核标准，结合田间实测土壤物理参数建立水稻土离散元耕作模型。以深松阻力行为作校核显示仿真耕作阻力与田间实测阻力误差相对较小，在可接受范围内，表明DEM模型预测的匹配与适用性。但是鉴于田间试验条件地表虽经清茬处理仍然存在根茬和秸秆残留，而且虽然考虑到水稻土在深度方向的分层结构但并未描述其他方向的土壤力学空间变异，因此在土壤宏观扰动状态方面仿真与实测结果差异较大，尤其是起垄高度在某些点出现较大误差。总体上看，仿真结果与田间试验结果在耕作力学及土壤扰动两方面的效应基本一致，表明针对粘性水稻土的离散元耕作模型的土壤接触模型选定、仿真参数设定以及双层土壤模型建立的准确性。

田间试验和仿真结果都表明无论是耕作阻力(包括垂直耕作阻力)还是土壤宏观扰动，200 mm耕作深度是水稻土耕作深度的特征临界点，该临界深度对合理选择适宜的深松深度具有一定理论参考价值[9]，耕作的深度一定要因地制宜，尤其对于深松这种高能耗作业，相近工况下华北平原壤土区深松耕作水平阻力为1.80 kN[13]、东北平原黑土区土壤耕作阻力为3.86 kN[31]，然而在本文水稻土条件下，300 mm耕作深度时深松耕作牵引阻力已达到5.21 kN，这充分反映出湿粘水稻土深松作业与旱作土的差异性。因此，不同土壤条件的作业参数各不相同，且耕作深度等作业参数也需要进一步配合农艺要求。

耕层土壤剪切模量小、颗粒粘结强度低，底层土壤剪切模量大、颗粒粘结强度高，耕作阻力变化趋势如图5a所示。土壤宏观扰动参数随着耕作深度先增大后减小，在200 mm耕作深度时宏观扰动参数达到最大，这显然是土壤自身属性和深松铲特点共同作用的结果。当耕作深度较浅(100 mm)时，深松铲的触土效应主要由铲尖和少部分曲柄段提供，此时耕作阻力小，对土壤的扰动也小；当耕作深度逐步增大到200 mm时，深松铲铲尖部分和曲柄弧段完全作用于耕层土壤和部分底层土壤，耕作阻力逐渐增大，土壤扰动量增大，曲柄弧形结构的抬土功能得到充分发挥，此时地表起垄参数达到最大值；当耕作深度继续加到300 mm时，深松铲作用于耕层土壤和底层土壤，由于底层土壤剪切模量大、粘结强度高，使得耕作阻力急剧增大，深松铲对土壤的扰动能力被大大削弱，地表起垄参数降低。由此可见，粘性水稻土的DEM模型仿真为土壤宏观扰动过程提供了机理性解释，DEM仿真也是实现不同耕作效应的耕作机具设计的重要方法学支撑。

3.2 深松土壤扰动微观机理

3.2.1 土壤扰动微观过程分析

在仿真环境下选取耕作方向上中间段厚度为100 mm的局部土体，以300 mm耕作深度为例分析深松铲对土壤的微观作用过程。如图7所示，分别为耕作方向及其侧向在不同时刻的截图，图8为土体内部粘结结构变化情况，图9为不同耕作深度下土体内部粘结断裂总数随时间的变化曲线。

图7 不同时刻深松铲对土壤扰动仿真结果Fig.7 Simulation result of soil disturbance of subsoiler at different time

图8 关键时刻土体内部粘结图示Fig.8 Diagrams of soil internal particle bond at important moment

图9 各耕作深度下粘结断裂总数随时间的变化曲线Fig.9 Variation curves of total broken bond under different tillage depths

在3.3 s时刻深松铲铲尖即将进入土体(图7)，此时在牵引力的作用下铲尖部分通过挤压土壤“凿”进土体，此时(图9)粘结断裂数逐步增加，土壤开始产生破碎；在5.4 s时刻深松铲曲柄刃段作用于土体，此时由于曲柄段的弧形结构使得作用区域内的土体被向上前方“抬”升(图7)，此时刻粘结图示(图8)深松铲作用下土体内部以抗剪强度为极限逐步形成深松扰动边界，扰动区域内土壤在深松铲作用下被抬升，达到松土效果；在6.6 s时刻深松铲直柄端作用于土体(图7)，此时随着深松铲的前进直柄端将前一阶段抬升起来的土体进行“剪”切，在深松铲两侧地表形成土垡及其裂纹，粘结图示(图8)表层下方土壤在直柄的剪切作用下进一步破碎，达到碎土效果；在5.4～6.6 s时段内(图9)粘结断裂总数剧增，这是由于深松铲曲柄段和直柄端对土壤的抬升、剪切作用土体内部发生着剧烈的剪切失效使土壤充分破碎；在7.5 s时刻深松铲完全离开土体(图7)，此时土体会形成深松后的沟槽，沟槽两侧土壤会有一部分回“落”，回填沟槽，土壤粘结断裂总数趋于一个稳定值(图9)，也说明深松铲对这一局部区域土体扰动完成。

传统田间试验条件下，传统的试验方法不便于深入理解认识深松铲-土壤的作用机理。已有的研究[10-11,32]表明凿型深松铲向前上方挤压土壤，随着深松铲的不断前进土壤产生剪切裂纹，破裂线从铲尖延伸到地表(图8)并发生破碎。本文从微观的DEM仿真角度进一步解释和论证了宏观试验测试结果及假说的正确。文献[11]在西北地区塿土土壤类型下对深松土壤扰动过程进行了试验与仿真研究，着重分析了深松铲对不同位置土壤的扰动机理，对不同层土壤的运动状态进行了对比分析。但是其并未结合深松铲自身结构特点分析机具-土壤作用过程。本文DEM仿真能够更深一步从微观层面充分结合深松铲的结构特点论证从铲尖到铲柄对土壤作用的各个细观过程，进一步明确了深松铲各部分结构的触土功能特点。同时结合仿真模型从土颗粒接触力学层面尝试表达了机具对土壤扰动的边界生成、土体内部破碎等过程，更加深入明确了深松铲对湿粘型水稻土的机械扰动过程，为优化深松铲结构设计及更好适应南方水稻土深耕松提供了理论参考。

图10 代表性土壤颗粒位置Fig.10 Location of representative soil particles

3.2.2 土壤扰动位移微观分析

图11 代表性颗粒不同方向位移随时间的变化曲线Fig.11 Changing curves of displacement of representative particles with time in different directions

为进一步定量分析深松对土壤的扰动过程，描述土壤在深松扰动过程的位置变化至关重要。参照文献[15]的方法，选取代表性位置土壤颗粒(即深松铲耕作路径上距地表0、100、200、300 mm处，如图10中橙色颗粒)分析其在深松铲扰动过程中X、Y、Z方向的位移变化趋势。由图11可知，在耕作方向(X)上0、300 mm处的表层和底层土颗粒运动位移较大，300 mm处颗粒较其他位置先被扰动；在侧向(Y)上0 mm处表层土颗粒侧向扰动位移明显大于其他位置，但总体位移变化不明显；在纵向(Z)上，开始扰动3.3～5.4 s阶段在深松铲的铲尖和曲柄刃的抬升作用下各位置处土颗粒都向上运动，5.4～7.5 s时段深松铲直柄端作用于这部分土壤，表层0、100 mm处土颗粒继续向上运动，而较底层200、300 mm处土颗粒则向下运动。

对3个方向的位移分析表明，深松铲触土过程的土壤主要在耕作方向和纵向上运动，侧向扰动不明显。这与陈坤等[10]、 郭新荣[32]的研究结果一致，深松过程中土壤主要是向上前方运动。但是结合前文对扰动过程的定性和此处纵向位移定量分析表明，深松扰动过程中土壤有一个明显的抬升然后再回落的过程，黄玉祥等[11]的仿真研究也发现类似的现象。这说明深松过程土壤在纵向上先被深松铲抬起，随着深松铲的不断前进，土壤再向下回落，因此土壤不只单单的向上前方运动，而是耕作方向上向前运动而纵向上先向上抬升再向下回落的运动轨迹(图11c)。除此之外，粘性水稻土在整个过程中伴随剪切失效进而发生破碎，达到松碎土壤、调理耕层的效果。

3.2.3 土壤破碎程度对比分析

南方湿粘型水稻土的机械耕作特征显著不同于旱作制，古老东方文明有别于西方的一个典型的特征是稻作制的触土凸面设计强化了湿粘型水稻土的碎土效应[33]。然而鉴于深松铲并不具备犁的翻耕效应，因此强化深松铲在湿粘型水稻土条件的碎土效果就愈发成为现代保护性耕作理论的重要课题。耕后土壤破碎体的尺度分析是土壤宏观力学结构研究的基本内容[34]，传统的筛分法可以利用碎土系数从宏观上反映耕后土壤的破碎程度[35]，但这些都是基于耕后结果的后评价。DEM仿真的优势之一是有可能从触土过程中定量土壤的破碎程度。基于这一考虑，本文从仿真模拟手段在微观角度定义一个反映土壤破碎程度的指标：断裂系数I，其计算公式为

(6)

式中I——断裂系数Nb——扰动区域土壤颗粒粘结断裂数Ni——扰动区域土壤颗粒完整粘结数

扰动区域的边界划定为:以耕作方向上耕作稳定段中间的100 mm为区域长度，以耕作后起垄宽度作为区域宽度，耕作深度为区域深度。仿真结束后，提取每个耕作深度下扰动区域内的数据如表4所示。

表4 扰动区域粘结断裂数据Tab.4 Bond broken data in disturbed zone

图12是田间试验后测得不同耕作深度下土壤破碎体尺度分布情况，由图12可知本文水稻土深松耕作后，土壤破碎体粒径主要分布在100～400 mm区间(质量百分比大于50%)。耕作后土壤破碎体尺度分布能有效反映土壤的破碎程度，参照文献[35]的方法以小于200 mm的土块质量占总碎土质量的百分比定义为碎土系数。当耕作深度在200 mm以内时，粒径小于200 mm的土壤破碎体所占质量百分比处于一个相对较低的状态，分别为31.16%、28.62%、21.62%；当耕作深度大于200 mm时则迅速增加，耕作深度300 mm时达到最大值，为61.19%。图13为田间实测不同耕作深度下碎土系数和仿真情况下土壤颗粒粘结断裂系数的变化趋势，两者变化趋势基本一致，不同耕作深度下的两者平均误差为3.46%。仿真情况下断裂系数在耕作深度为200 mm以内时，基本保持在30%左右；当耕作深度大于200 mm时，土壤破碎率开始增大，耕作深度300 mm时达到最大值52.52%。

图12 田间实测土壤破碎体径级分布Fig.12 Relative size distribution of measured soil clods

图13 土壤破碎程度定量对比Fig.13 Comparison of measured and simulated soil crushing degrees

以上结果表明，田间试验下耕作深度较小时(200 mm以内)，深松耕作产生的土壤破碎体主要是粒径大于200 mm的大土垡结构(破碎土壤质量的60%以上)，土体整体破碎程度不高，此时深松扰动的效果是使被扰动土体破而不碎；仿真模拟环境下的结果显示，此时土体内部土壤颗粒粘结断裂系数相对较低，说明土颗粒之间的粘结断裂不充分、不完全，呈现出的结果就是土体破裂产生大土垡结构，土壤松动、破而不碎。随着耕作深度加大(超过200 mm)，土壤颗粒粘结断裂系数增大，此时对比田间试验土壤破碎体尺度分布结果也能印证这一点，耕作深度250、300 mm时碎土系数分别达到50.07%和61.19%，其中粒径小于100 mm小尺度碎块所占比例也分别达到23.02%、23.93%，整体破碎程度明显高于其他耕作深度情况。

以上对比分析说明，仿真环境下模型土壤颗粒粘结断裂系数指标也可以定量表达田间试验下土壤的破碎程度，其大小及变化趋势与传统试验下的碎土系数基本一致。这一指标也更有利于对深松土壤破碎过程微观机理的表达和理解。旱作制的深松铲优化设计强调触土曲面优化[36-37]、自激或受迫振

动[38-39]、仿生减阻[40-41]等措施，例如文献[13]针对华北平原壤土区不同土层特点设计新型折线破土刃深松铲，通过优化深松铲触土曲面强化破土碎土过程从而达到节能降阻的目的。其旨在通过强化触土过程达到减小耕作阻力的目的，并未过多关注耕后土壤破碎效果。然而南方水稻土的机械化深耕松更需要解决有效碎土及高效耕作层扰动的问题，本文针对水稻土耕作力学特性的离散元建模及仿真显示可以从深松铲的工作参数(耕作深度、耕作速度等)优化入手提高湿粘型土壤的深耕松效应。

4 结论

(1)基于田间实测的南方稻麦轮作区域水稻土耕层土壤、底层土壤的物理性质差异，选用带有粘结作用的离散元接触模型，结合田间试验及虚拟仿真标定的方法建立了离散元双层深松土壤模型。其中耕作阻力平均误差6.63%，土壤扰动垄宽平均误差4.39%、垄高平均误差19.22%。该模型的建立可为特定区域土壤耕作机具的研究提供理论参考模型。

(2)对深松土壤的微观扰动过程进行了分析，进一步明确了深松铲与土壤相互作用机理。从仿真模型粘结断裂的角度表达了深松扰动边界的产生过程，土壤产生破碎的过程。对土壤运动位移的分析表明，土壤在耕作方向、垂直方向位移较大，侧向土壤扰动不明显。为南方因地制宜开展合理深松技术，例如新型深松铲的优化设计、合理的深松铲布局方式、最佳耕作深度的选择等提供了科学的决策依据。

(3)应用离散元EDEM软件，提出反映土壤破碎程度的指标“断裂系数”。本研究中当耕作深度在200 mm以内时，断裂系数稳定在30%左右；耕作深度大于200 mm时，断裂系数逐渐增大到52.52%。与传统田间试验下“碎土系数”指标无论是在变化趋势，还是相对误差都在可接受范围内。断裂系数可以作为衡量耕后土壤破碎程度的指标，断裂系数越大土壤破碎程度越高，反之亦然。

1 王恩姮，柴亚凡，陈祥伟. 大机械作业对黑土区耕地土壤结构性特征的影响[J]. 应用生态学报，2008，19(2): 351-356. WANG Enheng, CHAI Yafan, CHEN Xiangwei. Effects of heavy machinery operation the structural characters of cultivated soils in black soil region of Northeast China[J]. Chinese Journal of Applied Ecology, 2008, 19(2): 351-356. (in Chinese)

2 孙忠英，李宝筏. 农业机器行走装置对土壤压实作用的研究[J]. 农业机械学报，1998，29(3): 172-174. SUN Zhongying, LI Baofa. Research on the effect of soil compaction caused by the device of running agricultural machinery[J]. Transactions of the Chinese Society for Agricultural Machinery, 1998, 29(3): 172-174. (in Chinese)

3 雷友，曹国鑫，牛新胜，等. 土壤深耕对冬小麦根系在土壤剖面分布的影响[J]. 现代农业科技，2011(8): 272-273. LEI You, CAO Guoxin, NIU Xinsheng, et al. Influence of deep tillage on distribution of wheat root system in soil profile[J]. Modern Agricultural Sciences and Technology, 2011(8): 272-273. (in Chinese)

4 杨建昌，杜永，刘辉. 长江下游稻麦周年超高产栽培途径与技术[J]. 中国农业科学，2008，41(6): 1611-1621. YANG Jianchang, DU Yong, LIU Hui. Cultivation approaches and techniques for annual super-high-yielding of rice and wheat in the lower reaches of Yangtze River[J]. Scientia Agricultura Sinica, 2008, 41(6): 1611-1621. (in Chinese)

5 何进，李洪文，高焕文. 中国北方保护性耕作条件下深松效应与经济效益研究[J]. 农业工程学报，2006，22(10): 62-67. HE Jin, LI Hongwen, GAO Huanwen. Subsoiling effect and economic benefit under conservation tillage mode in Northern China[J]. Transactions of the CSAE, 2006, 22(10): 62-67. (in Chinese)

6 孙少华，张文斌. 机械化深松技术及苏州地区适应性研究[J]. 安徽农业科学，2012，40(20): 10701-10702. SUN Shaohua, ZHANG Wenbin. Research of subsoiling techniques and their adaptability in Suzhou[J]. Journal of Anhui Agricultural Sciences, 2012, 40(20): 10701-10702. (in Chinese)

7 张文斌，孙少华，李志刚. 苏州地区机械化深松技术试验研究与分析[J]. 农业装备技术，2012，38(5): 20-23.

8 张军，丁启朔，丁为民，等. 分层旋耕对水稻土耕层物理性状的影响[J/OL]. 南京农业大学学报，2015, 38(6): 1016-1022.http:∥nauxb.njau.edu.cn/oa/DArticle.aspx?type=view&id=201506022.DOI:10.7685/j.issn.1000-2030.2015.06.022. ZHANG Jun, DING Qishuo, DING Weimin, et al. Effect of stratified rotary tillage on paddy soil physical properties[J/OL]. Journal of Nanjing Agricultural University, 2015, 38(6): 1016-1022. (in Chinese)

9 丁启朔，葛双洋，任骏，等. 水稻土的深松阻力及土壤扰动效果研究[J/OL]. 农业机械学报，2017，48(1):47-56. http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?file_no=20170107&flag=1. DOI：10.6041/j.issn.1000-1298.2017.01.007. DING Qishuo, GE Shuangyang, REN Jun, et al. Characteristics of subsoiler traction and soil disturbance in paddy soil[J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2017,48(1):47-56. (in Chinese)

10 陈坤，胡晓丽，赵新子，等. 国内外深松铲研究现状与展望[J]. 农业与技术，2010，30(3): 30-34.

11 黄玉祥，杭程光，苑梦婵，等. 深松土壤扰动行为的离散元仿真与试验[J/OL]. 农业机械学报，2016,47(7):80-88. http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?file_no=20160712&flag=1. DOI:10.6041/j.issn.1000-1298.2016.07.012. HUANG Yuxiang，HANG Chengguang，YUAN Mengchan,et al. Discrete element simulation and experiment on disturbance behavior of subsoiling[J/OL].Transactions of the Chinese Society for Agricultural Machinery, 2016,47(7): 80-88. (in Chinese)

12 黄玉祥，杭程光，李伟，等. 深松作业效果试验及评价方法研究[J]. 西北农林科技大学学报:自然科学版，2015，43(11):228-234. HUANG Yuxiang，HANG Chengguang，LI Wei, et al.Subsoiling test and evaluation methodology of tillage quality[J].Journal of Northwest A&F University: Natural Science Edition,2015,43(11):228-234.(in Chinese)

13 郑侃，何进，李洪文，等. 基于离散元深松土壤模型的折线破土刃深松铲研究[J/OL]. 农业机械学报，2016，47(9):62-72. http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?file_no=20160910&flag=1. DOI:10.6041/j.issn.1000-1298.2016.09.010. ZHENG Kan，HE Jin，LI Hongwen，et al. Research on polyline soil breaking blade subsoiler based on subsoiling soil model using discrete element method [J/OL].Transactions of the Chinese Society for Agricultural Machinery,2016,47(9): 62-72.(in Chinese)

14 丁启朔，丁为民. 现代土壤机械耕作研究的综述[J]. 土壤通报， 2006，37(1): 149-153. DING Qishuo, DING Weimin. Modern soil tillage research: a review[J]. Chinese Journal of Soil Science, 2006, 37(1): 149-153. (in Chinese)

15 方会敏，姬长英，FARMAN Ali Chandio，等. 基于离散元法的旋耕过程土壤运动行为分析[J/OL]. 农业机械学报，2016，47(3): 22-28.http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?file_no=20160304&flag=1. DOI:10.6041/j.issn.1000-1298.2016.03.004. FANG Huimin，JI Changying，FARMAN Ali Chandio, et al. Analysis of soil dynamic behavior during rotary tillage based on distinct element method[J/OL].Transactions of the Chinese Society for Agricultural Machinery, 2016,47(3): 22-28.(in Chinese)

16 方会敏，姬长英，Ahmed Ali Tagar，等. 秸秆-土壤-旋耕刀系统中秸秆位移仿真分析[J/OL]. 农业机械学报，2016，47(1): 60-67. http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?file_no=20160109&flag=1. DOI:10.6041/j.issn.1000-1298.2016.01.009. FANG Huimin,JI Changying,AHMED Ali Tagar,et al. Simulation analysis of straw movement in straw-soil-rotary blade system [J/OL].Transactions of the Chinese Society for Agricultural Machinery,2016,47(1): 60-67.(in Chinese)

17 李博，刘凡一，陈军，等. 深松铲耕作阻力影响因素的离散元法仿真分析[J]. 农机化研究， 2015(2): 71-74. LI Bo, LIU Fanyi, CHEN Jun, et al. DEM Simulation on analyzing affecting factors of the soil resistance applied on the subsoiler [J].Journal of Agricultural Mechanization Research,2015(2): 71-74.(in Chinese)

18 邓佳玉，胡军，李庆达，等. 基于EDEM离散元法的深松铲仿真与试验研究[J]. 中国农机化学报，2016，37(4): 14-18. DENG Jiayu,HU Jun,LI Qingda,et al.Simulation and experimental study on the subsoiler based on EDEM discrete element method [J].Journal of Chinese Agricultural Mechanization, 2016, 37(4): 14-18.(in Chinese)

19 杨艳山，丁启朔，丁为民，等. 田间原位综合耕作试验台设计与应用[J/OL]. 农业机械学报，2016， 47(1): 68-74. http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20160110&journal_id=jcsam. DOI:10.6041/j.issn.1000-1298.2016.01.010. YANG Yanshan, DING Qishuo, DING Weimin, et al. Design and application of multi-purpose in-situ tillage tool testing platform[J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2016, 47(1): 68-74. (in Chinese)

20 WATTS C W, DEXTER A R. Traffic and seasonal influences on the energy required for cultivation and on the subsequent tilth[J]. Soil & Tillage Research,1994, 31(4): 303-322.

21 胡建平，周春健，侯冲，等. 磁吸板式排种器充种性能离散元仿真[J/OL]. 农业机械学报，2014， 45(2): 94-98.http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20140216&journal_id=jcsam. DOI：10.6041/j.issn.1000-1298.2014.02.016. HU Jianping, ZHOU Chunjian,HOU Chong, et al. Simulation analysis of seed-filling performance of magnetic plate seed-metering device by discrete element method[J/OL]. Transactions of the Chinese Society for Agricultural Machinery,2014, 45(2): 94-98.(in Chinese)

22 张绍军. 3种深松铲对深松作业效果影响研究[J]. 农业科技与装备，2012(6): 23-25. ZHANG Shaojun.Experimental analysis on the subsoiling effect of three subsoilers[J]. Agricultural Science & Technology and Equipment,2012(6): 23-25.(in Chinese)

24 EDEM. EDEM user guide[M]. Edinbuegh: EDEM-Solutions, 2010.

25 POTYONDY D O, CUNDALL P A. A bonded-particle model for rock[J]. International Journal of Rock Mechanics & Mining Sciences,2004, 41(8): 1329-1364.

26 AHMAD F. Performance evaluation of disc type furrow openers in no-till paddy soil environment[D]. Nanjing：Nanjing Agricultural University, 2015.

27 张永良. 基于离散元的逆转旋耕施肥播种机抛土性能仿真及试验研究[D].镇江：江苏大学, 2012.

28 ASAF Z, RUBINSTEIN D, SHMULEVICH I. Determination of discrete element model parameters required for soil tillage[J]. Soil & Tillage Research, 2007, 92(1-2): 227-242.

29 LINDE J V D. Discrete element modeling of a vibratory subsoiler[D]. Stellenbosch: Stellenbosch University, 2007.

30 MAK J, CHEN Y, SADEK M A. Determining parameters of a discrete element model for soil-tool interaction[J]. Soil & Tillage Research,2011, 118(5): 117-122.

31 齐关宇, 刘林, 赵艳忠, 等. 深松铲入土深度及铲形对耕作阻力影响[J]. 农机化研究，2015(11):161-165. QI Guanyu, LIU Lin, ZHAO Yanzhong, et al. Effects of subsoiler’s penetrating depth and spade shape on traction resistance[J]. Journal of Agricultural Mechanization Research, 2015(11):161-165. (in Chinese)

32 郭新荣. 土壤深松部件松土机理的分析[J]. 农机化研究，2005(3): 78-80. GUO Xinrong. Comparative study on mechanism of mutil-subsoiler deep loosening the soil[J]. Journal of Agricultural Mechanization Research,2005(3): 78-80.(in Chinese)

33 沈昌蒲，王秋华，温锦涛. 从犁的工作部件看中国和西欧设计思想的差异[J]. 农业机械学报， 2007，38(6): 178-183. SHEN Changpu,WANG Qiuhua,WEN Jintao. Difference of design idea of plow share between Chinese and West European[J]. Transactions of the Chinese Society for Agricultural Machinery,2007, 38(6): 178-183.(in Chinese)

34 丁启朔，丁为民，潘根兴，等. 土壤宏观力学结构与精准耕作[J]. 中国农业科学，2012，45(1): 26-33. DING Qishuo, DING Weimin, PAN Genxing, et al. Soil macro-structure and precision tillage[J]. China Agriculture Science, 2012, 45(1): 26-33. (in Chinese)

35 杨有刚，张宏，冯涛，等. 土壤浅深松联合松耕机设计和浅松土试验[J]. 机械工程学报， 2012，48(19): 163-168. YANG Yougang, ZHANG Hong, FENG Tao,et al. Topsoil and subsoil combined cultivator and top-soiling experiment[J]. Journal of Mechanical Engineering,2012,48(19): 163-168.(in Chinese)

36 郭志军，杜干，周志立，等. 土壤耕作部件宏观触土曲面减阻性能研究现状分析[J]. 农业机械学报，2011，42(6): 47-52. GUO Zhijun, DU Gan, ZHOU Zhili, et al. Actuality analysis of resistance-reducing properties on soil cultivating components with different macroscopic soil-engaging surfaces[J]. Transactions of the Chinese Society for Agricultural Machinery, 2011, 42(6): 47-52. (in Chinese)

37 张强，张璐，于海业，等. 复合形态深松铲耕作阻力有限元分析与试验[J/OL]. 农业机械学报， 2012，43(8): 61-65.http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20120812&journal_id=jcsam.DOI:10.6041/j.issn.1000-1298.2012.08.012. ZHANG Qiang,ZHANG Lu,YU Haiye,et al. Finite element analysis and experiment of soil resistance of multiplex-modality subsoiler[J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2012,43(8): 61-65.(in Chinese)

38 辛丽丽，李传琦，梁继辉，等. 考虑分段土壤作用力的振动减阻分析[J/OL]. 农业机械学报，2014, 45(2): 136-140. http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20140223&journal_id=jcsam. DOI:10.6041/j.issn.1000-1298.2014.02.023. XIN Lili, LI Chuanqi, LIANG Jihui, et al. Vibrating drag reduction considering acting force of piecewise soil[J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2014, 45(2): 136-140. (in Chinese)

39 蒋建东，高洁，赵颖娣，等. 土壤旋切振动减阻的有限元分析[J/OL]. 农业机械学报，2012，43(1):58-62.http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20120112&journal_id=jcsam. DOI:10.6041/j.issn.1000-1298.2012.01.012. JIANG Jiandong, GAO Jie, ZHAO Yingdi, et al. Finite element simulation and analysis on soil rotary tillage with external vibration excitation[J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2012, 43(1):58-62.(in Chinese)

40 张琰，黄河，任露泉. 挖掘机仿生斗齿土壤切削试验与减阻机理研究[J/OL]. 农业机械学报，2013, 44(1): 258-261.http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20130148&journal_id=jcsam. DOI:10.6041/j.issn.1000-1298.2013.01.048. ZHANG Yan,HUANG He,REN Luquan.Experiment and drag reduction mechanism of bionic excavator tooth during soil cutting[J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2013, 44(1): 258-261.(in Chinese)

41 张金波，佟金，马云海. 仿生减阻深松铲设计与试验[J/OL]. 农业机械学报，2014，45(4): 141-145.http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20140422&journal_id=jcsam. DOI:10.6041/j.issn.1000-1298.2014.04.022. ZHANG Jinbo, TONG Jin,MA Yunhai. Design and experiment of bionic anti-drag subsoiler [J/OL].Transactions of the Chinese Society for Agricultural Machinery, 2014, 45(4): 141-145.(in Chinese)

DEM Analysis of Subsoiling Process in Wet Clayey Paddy Soil

DING Qishuo1REN Jun1BELAL Eisa Adam1,2ZHAO Jikun1GE Shuangyang1LI Yang1

(1.CollegeofEngineering,NanjingAgriculturalUniversity,Nanjing210031,China2.DepartmentofAgriculturalEngineering,ZalingeiUniversity,Zalingei06,Sudan)

Tool design and working parameters for subsoiling in wet clayey paddy soil are different from those for dry land cropping system. Investigation on the subsoiling process and the governing mechanisms requires detailed description on the underlying influential factors. Field soil parameters were used as a reference to construct a DEM model for subsoiling which was suitable for wet clayey soil analysis and performed in a EDEM software. The proposed DEM model of clayey paddy soil was then implemented to simulate a subsoiling process, in which directly measured draft and macroscopic disturbance of field soil was compared. Microscopic process and mechanisms of soil fragmentation were assessed with simulated stages. Results showed that simulation error with the proposed DEM model was less than 6.63%. Mean error of the resulted micro-relief was 4.39%. Mean error of elevated ridge was 19.22%. The measured results from microscopic approach were re-evaluated and the correctness of related assumptions from microscopic approach was confirmed by DEM results, soil failure boundary evolution and soil failure with particle contacting models were depicted. A soil fragmentation index was proposed for describing soil fragmentation behavior from DEM perspectives. The error of fragmentation index with respect to measured soil fragmentation parameter was less than 3.46%, affirming that the proposed parameter was a suitable tool for microscopic description of soil failure imposed by subsoiling. This research provides a powerful technical basis for the interaction mechanism between soil and soil, and lays a theoretical foundation for the optimization design of the soil contact parts.

subsoiling; soil disturbance; DEM simulation; paddy soil

10.6041/j.issn.1000-1298.2017.03.005

2016-11-14

2016-12-16

国家重点研发计划项目(2016YFD0300900)、农业部现代农业装备重点实验室开放课题项目(201601002)和江苏省农机基金项目(201-051028)

丁启朔(1968—)，男，教授，主要从事土壤耕作力学研究，E-mail: qsding@njau.edu.cn

赵吉坤(1977—)，男，副教授，主要从事土壤数值模拟研究，E-mail: jikunzhao_2006@163.com

S222.12+9

A

1000-1298(2017)03-0038-11