土壤在耕作过程中土壤与机械间的相互作用是非常复杂的力学问题。土壤在耕作过程中由连续体向非连续离散体动态变化,耕作中实现了破碎、松散和运移的动力学过程演化。

早期学者采用三维非线性有限元方法研究土壤剪切,Chi等^[1]基于此方法分析窄齿耕作部件的土壤切削问题,通过有限元计算得到了耕作部件的耕作阻力;在此基础上,进一步分析了土壤及耕作部件表面的应力分布情况、土壤破坏区间及位移情况^[2]。然而由于耕作土壤是一种松散的、结构复杂的土壤,所以要准确研究土壤动态行为的影响因素和变化规律,一般使用离散单元法,从构成土壤的土粒细观角度进行模拟分析,许多研究者采用此法研究了土壤相关问题^{[3, 4, 5, 6, 7, 8, 9, 10]}。在研究农机具与土壤力学方面,Seidi等^[11]以深松铲为研究对象,对其在不同工作速度、不同耕作深度、不同入土角度下的工作过程进行模拟,得出铲尖在水平和垂直方向受到的阻力大小和变化趋势,与试验结果相一致,为减小阻力提供基础。为揭示波纹表面的降阻机制,刘国敏等^[12]利用离散元进行仿真,得出波纹表面与土壤间存在孔隙,从而使得接触面积减小、阻力降低。在此基础上,张锐等^[13]建立了颗粒非线性力学模型,并对波纹形复杂表面推土板前端土壤的动态扰动行为进行了离散元宏细观模拟。徐泳等^[14]提出了采用基于颗粒接触力学原理的离散元法对耕作土壤动力学问题进行数值模拟的新思路。

以上关于土壤耕作研究中没有考虑含秸秆土壤耕作非线性问题。笔者基于离散单元法,借助EDEM软件建立含秸秆土壤离散元模型,进行不同耕作速度下土壤耕作动态过程,分析秸秆、土壤颗粒和圆盘犁面非线性运动及动态演化情况。

土壤破碎是土壤在外力作用下从黏结的整体状态分离变成较小土块的过程,破碎过程是一个非常复杂的土块尺寸变化过程,与许多因素相关,主要包括:破碎过程的动力学参数,如土块的尺寸及破碎过程的运动状态;土颗粒的力学性质,如抗破碎能力(或者强度)、硬度、韧性、形状、尺寸、湿度、密度和均匀性等;外部条件,如土块破碎生成的破碎颗粒在破碎瞬间的相互作用、分布情况及操作条件等。

基于连续介质力学,土壤颗粒流变本构方程^[15]表示为:

式中:F为张量函数;D_ij为变形速度张量,,V_i和V_j分别为颗粒在i、j方向的速度矢量。式(1)进一步分解为:

基于以上公式,得到颗粒体系的本构方程为:

在单剪切土壤耕作切割作用条件下,无黏颗粒体系(c=0)的本构方程简化为:

在上述关系中,同时存在流速梯度的零次项、线性项和二次项。零次项是由于土壤颗粒之间静态支撑作用引起的,此时体系内部形成了相对稳定的力链网络;二次项是由于土壤颗粒之间的碰撞和扩散作用引起的;线性项是由于颗粒之间的相对滑动和挤压作用引起的,该项与普通流体的线性黏滞应力类似,但两者的机制有根本不同,颗粒体系的线性项是由于力链网络频繁形成和断裂引起的。

受黏塑性宾汉流体行为的启发,为复现颗粒在不同边界条件下的复杂流动,获得速度分布和颗粒浓度分布等重要参数,颗粒体系的应力分量写为:

影响土壤颗粒耕作过程中流动特性的主要因素包括颗粒体积分数、剪切速率、密度、颗粒大小、刚度、摩擦系数和恢复系数等。为更合理地描述土壤颗粒物质在耕作过程中由固态到流态的力学行为,可建立一个广义本构方程,即

考虑土壤颗粒在固态-流态转变过程中,应力张量与剪切速率的关系,上式可以进一步写成

建立土壤耕作试样,试样尺寸长度为200 mm,高度为100 mm,颗粒半径在0.8~1.3 mm范围内随机变化,按照孔隙比为0.1,在矩形区域内随机生成。力学参数见表 1所示。试样如图 1所示。

	图 1 土壤耕作试样 Fig. 1 Soil tillage sample
图选项

为了模拟种植草由于其根茎的生长,对周边土壤颗粒的凝聚作用,在此用较大颗粒,并施加竖向荷载5.0×10³ kN产生土壤颗粒间力链网,进而实现凝聚土体颗粒作用。推土板的倾斜角度为45°,推土板高度为20 mm,宽度为40 mm,入土深度为5 mm,给推土板施加水平向位移为每步2.5 mm。试验得出的推土板阻力与荷载步的关系曲线如图 2所示;种植草在边缘与中部的耕作过程见图 3所示。

	图 2 推土板前缘阻力随时间变化曲线 Fig. 2 Curve of compression force and step 散点为数值模拟结果;实线为拟合结果。 Scatter is the results of numerical simulation;the solid line is the fitting results.
图选项

	图 3 黏性土耕作模拟过程 Fig. 3 Cohesive soil tillage simulation process
图选项

通过图 2可以看到,推土板在没有遇到草颗粒时受阻基本均衡,推土板前缘翻耕到草颗粒固土范围内时,前缘阻力明显增大,峰值阻力是平均阻力的3.5倍。当翻耕通过草颗粒固土范围后,前缘阻力趋近于平均阻力。推土板阻力所表现出的波动现象主要是由于土体颗粒尺寸与固土的强度造成的推土板前缘阻力的差异性。草颗粒在中间与边缘处对于推土板的阻力基本一致。

综合图 2和图 3可以得出:在土壤耕作过程中,推土板前缘较后缘阻力大,前缘有明显的受阻区域,上层土壤颗粒在推土板的作用下发生明显的非线性运动过程,即稳固静态、准静-流动转变、流动变形的整个过程。前缘影响的区域土体颗粒有一定的松动现象,土壤经过推土板耕作以后呈现明显的松软状态,孔隙比有一定的增加。在此模拟过程中的草颗粒由于自身根茎的生长,至使周围土体颗粒的黏结力显著提高,推土板在遇到该部分时,前缘阻力明显增加,需要耗费较大的机械功;当耕作过后前缘阻力趋于前期平均阻力。由于土体细观颗粒尺寸不同以及力学参数的差异,细观颗粒单元的黏结强度的不断损伤,发生颗粒黏结失效进而滑移、翻滚和流动的物理过程。

按照试验平台与数值模拟平台10∶1等比尺方式,搭建秸秆覆盖情况下土壤耕作室内试验平台(图 4-a),研究3种尺寸秸秆在圆盘犁作用下运移及土壤耕作力学机制问题,在试验平台上测定了圆盘犁工作阻力及不同尺寸秸秆的位移。建立土壤耕作模型平台及秸秆模型(图 4-b和4-c),将3种尺寸分别为7、13和19 mm的秸秆分布在土壤颗粒表层。建立圆盘犁面几何模型,保持耕深为10 mm的前提下,并施加给该模型不同的耕作速度,对比分析结果参数。耕作参数如表 2所示。在圆盘犁面前缘选取以上3种尺寸的秸秆,进行秸秆覆盖情况下农田耕作过程非线性动态过程模拟。

	图 4 耕作平台及秸秆模型 Fig. 4 The tillage platform and straw model
图选项

综合图 5和图 6可以看出:圆盘犁在耕作过程中,壁面阻力逐步增大;牵引速度越大,阻力越大,试验结果与模拟结果基本吻合。土壤在圆盘犁切割翻耕作用下发生明显的位移,圆盘犁前缘土颗粒隆起现象明显;圆盘犁身后出现明显的浅沟;由于秸秆材料力学性能低于土壤颗粒,土壤颗粒带动秸秆发生更为明显的位移。从图 6可以看出:牵引速度越大,圆盘犁面前缘隆起土颗粒越多,在一定程度上前缘所受土壤的阻力越大,达到峰值阻力的时间越短。随着耕作的逐步完成,犁面前缘阻力逐步回落。

	图 5 耕作过程图 Fig. 5 Soil tillage simulation process
图选项

	图 6 圆盘犁面阻力变化曲线 Fig. 6 Soil tillage simulation process
图选项

由图 7可以看出:当犁面前缘接触到秸秆时,秸秆速度有明显的变化,7 mm秸秆相比其他2类秸秆速度变化较大,同时发生明显的相对位移。由于秸秆初始摆放位置不同,且自身质量较轻,因此不同尺寸秸秆发生了不同的分离现象:19与7 mm秸秆向两侧位移明显,而13 mm秸秆沿犁面前缘位移明显,并在^{2.5 s}时,13 mm秸秆又在犁面作用下重新发生二次接触和分离,此刻又发生相对明显的速度变化。

	图 7 圆盘犁前缘秸秆速度变化 Fig. 7 The curve of straw velocity
图选项

由图 8可以看出:土壤颗粒在受到圆盘犁面耕作外力作用下,产生切向的位移滑动,当耕作力大于极限强度时,土壤颗粒产生分离断裂,表现为原来黏结在一起的土壤颗粒发生分离滑移,曲线达到极限峰值并出现下降段;随着耕作的不断推进,部分颗粒会受到二次切割分离,因此会出现二次峰值,可以看出二次峰值明显低于极限峰值强度。这是由于前期犁面切割作用,颗粒间的黏结力链发生损伤破坏,已经降低了该土壤颗粒相互黏结强度。不同的牵引速度下,颗粒间剪力变化基本一致,且数值模拟结果与理论结果吻合较好。土壤耕作过程中,圆盘犁面前缘较后缘触动影响区域不同,前缘有明显的触动影响区域,翻耕过后,由于土体颗粒间力链的断裂,呈现一定的松动现象和松软状态,孔隙比有一定的增加。耕作后的土壤表面呈现一定的起伏,主要是由于耕作过程改变了其原有土壤的表面形态,土颗粒间存在一定的摩擦效应,难以恢复起初形态。

	图 8 土壤颗粒剪力随时间变化 Fig. 8 The curve of soil particle shear
图选项

由图 9可以看出:不同长度的秸秆,在受到圆盘犁面作用下,发生不同的滚动旋转非线性运动;由于各自长度不同,旋转速率各有差异。研究发现,尺寸较大的19 mm秸秆发生旋转角速度最大,峰值角速度接近1100 (°)·s^-1;而质量较轻和尺寸较短的7和13 mm秸秆在与犁面接触后,翻滚的角速度波动明显,这是由于土层表面的平整性发生变化,以及和土壤颗粒的接触碰撞受影响程度较大;13 mm秸秆与犁面有一定的二次接触,这一点与图 7分析结果一致。

	图 9 秸秆角速度随时间变化 Fig. 9 The curve of straw angular velocity
图选项

从图 10和图 11综合来看:3种尺寸的秸秆在1.2 s之前,在土壤表层的X和Z向的位移分量随时间变化的趋势基本一致;然而从1.2 s开始,13 mm秸秆在2个方向的位移分量与7和19 mm秸秆的位移区别明显,这是由于耕作产生的土壤翻动以及土壤颗粒滚动与秸秆相互作用,致使不同尺寸的秸秆发生不同方向的位移。从数值模拟结果和试验结果对比分析发现,7 mm秸秆在受到翻耕土壤的影响下发生一定的埋置,不会发生较大的位移;19 mm的秸秆由于自身尺寸的影响,在犁面和翻耕土壤的作用下,产生部分埋置,扭转速度受到影响,位移不明显;而13 mm秸秆与犁面与翻耕土壤分离较为明显,该模拟结果与试验结果基本一致。

	图 10 秸秆沿X方向位移变化曲线 Fig. 10 The curve of straw X displacement
图选项

	图 11 秸秆沿Z方向位移变化曲线 Fig. 11 The curve of straw Z dispacement
图选项

搭建含秸秆土壤耕作试验平台,研究不同尺寸秸秆及耕作过程;基于离散单元法理论,建立了含秸秆土壤耕作模型,研究了不同耕作速度下土壤及秸秆非线性运动及力学特性规律,得出如下结论:

1) 耕作后的土壤表面呈现一定的起伏,主要是由于耕作过程改变了其原有土壤的表面形态。随着耕作的不断推进,部分颗粒会受到二次切割分离,因此会出现二次峰值,可以看出二次峰值明显低于峰值强度值,这是由于前期的损伤破坏,已经降低了该土壤颗粒黏结力。

2) 数值模拟分析了圆盘犁前缘受力变化情况,牵引速度越大,圆盘犁面前缘隆起土颗粒越多,在一定程度上前缘所受土壤的反作用力越大,达到峰值阻力的时间越短。随着耕作的逐步完成,前缘阻力降低明显。

3) 3种尺寸的秸秆在受到犁面作用下,7 mm秸秆和13 mm秸秆与犁面有二次接触,19 mm秸秆旋转速度较大;13 mm秸秆在1.2 s以后发生较大的位移,7 mm秸秆受翻耕土影响严重,产生较多埋置。