基于离散元法的铧式犁仿真优化分析及试验

2022-11-28刘明勇胡成龙谢柏林

浙江农业学报 2022年11期

刘明勇，胡成龙，*，谢柏林

(1.湖北工业大学，湖北武汉 430068； 2.湖北省农业机械工程研究设计院，湖北武汉 430068)

随着我国的快速发展，农业机械化水平不断提高，农业生产中的能耗问题也日益突出。我国现有耕地约1亿hm2[1]，在耕地耕作作业中，使用的农具种类并无太大的改变，铧式犁仍然被普遍的使用，但存在一些问题，为了保证耕作的效果，对犁体曲面进行再次设计和优化，减小犁体阻力。

离散元单元法(distinct element method，DEM)是首个多用途离散元素法建模软件，近年来被广泛地用来模拟土壤和农机具之间的相互作用，能够良好地反映出土壤颗粒的运动轨迹和颗粒之间的接触。Shmulevich[2]通过离散元软件，模拟不同曲面的叶片切割土壤，通过仿真与实验进行对比，通过离散元构建描述土壤-耕作相互作用的模型，其仿真结果与实验结果十分相近，表明DEM可以作为设计耕作形状过程中的预测模拟工具。

Godwin等[3]通过离散元软件建立了在不同工况下犁工作阻力的预测模型，探究犁在不同深度和速度下的阻力大小，并通过大量的实验验证其正确性，但对于犁体结构参数并没有进行过多的探究。Soni等[4]探究了犁面上不同形状的凸包对于土壤黏附和犁体阻力的影响，通过实验对比发现，凸包可以减少土壤黏附和阻力，但是其制造困难、加工昂贵，难以广泛应用。Bulgakov等[5]通过圆柱形的铧式犁研究了犁体阻力和工作速度的关系可通过以一个无量纲系数来估算，并进行实验验证其正确性。Ucgul等[6-8]使用EDEM中的滞回弹簧接触模型预测在凝聚力较低的土壤条件下的耕作工具阻力情况，其仿真结果和实验结果十分接近，并且还对两种模型下的土壤运动进行对比，得出此模型能更好地对土壤运动进行预测。此外，Ucgul等[9]通过仿真和实验进行对比研究铧式犁在不同工况下的土壤翻垡情况，得到了表层土壤翻垡深度与速度的关系，何永强等[10]使用离散元仅对犁体结构参数进行多因素分析，得出最佳参数。

上述的研究都证实离散元法能准确地预测土壤-机具的相互作用，并对铧式犁进行了理论和试验研究，对于铧式犁结构研究比较少，本文在此基础上对铧式犁的结构参数进行研究，讨论铧式犁的结构对于犁体-土壤相互作用和土壤扰动效果影响，并根据土壤实测参数建立了铧式犁与土壤的离散元力学模型，探究铧式犁在不同导曲线的情况下，对于犁体阻力的影响和扰动效果的情况，并对犁体结构进行参数分析，以期得到一个理想的铧式犁结构，为降低耕作阻力，提升效率提供理论依据。

1 犁体三维模型建立

1.1 犁体曲面成型方法

犁体作为铧式犁与土壤接触的主要工作部件，其不仅影响着土垡的翻转质量，而且对土壤耕作过程的动能消耗有很大影响。迄今为止，现阶段铧式犁曲面绘制方法主要分为水平元线法、倾斜直元线法、曲元线法等方法[7]。本文采用水平元线的方法绘制犁体曲面，以直线AB为元线，CD为导曲线，使直线AB沿着曲线CD运动改变其夹角θn并且始终与XOY平面平行，形成犁体曲面。其原理如图1所示。

图1 水平直元线法原理图Fig.1 Schematic diagram of horizontal straight line method

导曲线的形状由高度、开度、两端点切线夹角、铧刃起土角和直线段长度确定。高度取犁体曲面最大高度，当高度一定时，开度越小，则曲面越陡峭，碎土性能强，但阻力越大；开度值越大，则曲面平坦，起垡容易，阻力小，但碎土性能弱。为了探究导曲线对于铧式犁的工作阻力的影响，通过3种不同的曲线类型绘制导曲线，建立不同的犁体曲面，并保持其开度和高度不变，导曲线方程如表1所示，导曲线图形如图2所示。

表1 导曲线方程

图2 导曲线图形Fig.2 Lead curve graph

1.2 犁体曲面绘制

通过上述方法绘制导曲线，其他犁体基本数都按照设计参数，如表2所示。确定各直元线面上元线角度[11-12]，如表3所示。用这3种曲线在solidworks中建立铧式犁模型如图3所示。

表2 犁体参数表

表3 直元线面元线角度

a，三维模型；b，正视图A 犁体幅宽；h，胫刃线高；H，顶边线最大高度；θ，翼边线夹角。a， Three-dimensional model； b， Front view of plough body A tillage width； h， Shin edge line height； H， Top edge line maximum height； θ， Wing edge line included angle.图3 犁体三维模型图Fig.3 3D model diagram of plough body

2 离散元仿真试验

2.1 土壤颗粒接触模型

土壤颗粒间的接触模型十分重要，不同的颗粒模型会导致不同的结果。由Walton等[13-16]提出的线性内凝聚滞回弹簧模型(HCSM)在滞回弹簧模型中添加了线性凝聚模型，Ucgul[17]首次使用HCSM模拟土壤-犁体相互作用，并发现此模型考虑了土壤的形变，相比于线性接触模型(HMCM)更好地反映了土壤与机具的运动接触。粒子之间的内聚力通过直接向正常接触力添加内聚力来定义。

图4 颗粒模型Fig.4 Particle model

在滞回弹簧模型中，法向接触力定义为：

(1)

式中：δ为土壤颗粒法向重叠量，mm；kn为法向刚度；cn为法向阻尼；kn和cn分别为：

(2)

(3)

式中：R=(1/R1+1/R2)-1为等效半径，mm；m=(1/m1+1/m2)-1为等效质量，kg；ν为土壤泊松比；E′为等效弹性模量；β=π/lne，e为土壤恢复系数。

将颗粒间凝聚力添加后，总法向力为：

FN=Fn+Fc。

(4)

式中：Fc=Ac·ξ为凝聚力，N；Ac为颗粒接触面积，mm2；ξ为凝聚力系数。

其土壤切向接触力表示为：

(5)

式中：kt为切向刚度；ct为切向阻尼；μn为土壤摩擦系数。

2.2 土壤参数确定

在EDEM软件中，土壤参数直接影响仿真结果的正确性。为了使仿真结果更加准确，以南方地区的黏性土壤为研究对象，对其进行含水率测定和剪切实验，使用SL-TSD型土壤测试仪测量土壤含水率并使用ZJ型无极调速应变控制式直剪仪进行剪切实验，其一次可以同时测量4组数据，其仪器如图5所示。

a，SL-TSD型土壤测试仪；b，ZJ型无极调速应变控制式直剪仪。a， SL-TSD soil tester； b， ZJ type stepless speed regulation strain control direct shear instrument.图5 实验测试仪器Fig.5 Experimental test instrument

对4组土壤样本加载50、100、200、400 kPa的垂直载荷，得出其抗剪强度和垂直载荷关系如图6所示，得到土壤的内摩擦角为15.4°，内凝聚力为13.06 kPa。其他仿真参数如表4所示[3-9，13-17]。

a，剪应力与剪切位移变化曲线；b，拟合曲线。a， Shear stress and shear displacement change curve； b， Fitting curve.图6 抗剪强度与垂直载荷关系图Fig.6 Diagram of shear strength and vertical load

表4 仿真参数

2.3 离散元模型建立

将solidworks建立的犁体仿真三维模型导入EDEM中，其土壤的各项参数照表4设置。模型中建立长方形土槽模型尺寸为1 800 mm×1 000 mm×300 mm；建立模型的时候，提前将已经生成并沉降后稳定的颗粒导入模型中，其耕深、耕宽如表4所示，土壤颗粒之间的接触模型采用线性内凝聚滞回弹簧模型(HCSM)，犁体的前进速度为0.8 m·s-1，在进行仿真对比时，仅改变犁体三维模型其他参数保持不变，仿真模型如图7所示。

图7 EDEM仿真模型Fig.7 EDEM simulation model

为了更好地分析深度土壤的扰动状态，将土壤分为表层、中层和下层3层土壤，每层深度为100 mm，并将不同深度的土层颗粒设置为不同的颗粒，如图8所示。为了清晰地分析机具与土壤之间的相互运动，对模型沿着犁体运动方向和垂直运动方向进行剖视。图8-a为土壤分布纵向截面，分析犁体前进方向土壤的扰动行为；图8-b为土壤分布横向截面，分析犁体垂直于前进方向上的土壤扰动行为。

a，纵向分布；b，横向分布。a， Vertical distribution； b， Horizontal distribution.图8 犁体运动过程Fig.8 Plough body movement process

3 试验结果分析

3.1 扰动分析

3.1.1 土壤扰动分析

以抛物线曲面铧式犁为例，在仿真过程中选取几个时刻来观察耕作时土壤的扰动行为。为了清晰了解到土壤扰动随着犁体耕作过程的变化情况，从土壤的纵向剖面对不同层的土壤扰动进行分析，当犁体在不同位置的时候，土壤的扰动过程如图9所示，0.2 s时铧式犁铧的尖端进入了土壤，0.5 s时铧式犁绝大部分进入了土壤，并开始抛土，0.8 s铧式犁整个进入了土壤，进行持续性的抛土和挖土工作。

由图9可知，在整个工作过程中，表层土壤在整个工作过程中扰动范围最大，中层和下层次之，随着犁体不断地运动，其纵向和横向的扰动范围不断扩大。在0.2 s铧式犁的犁铧进入土壤时，下层土壤受到犁铧的挤压和剪切，沿着犁铧上升，并对表层和中层土壤造成扰动产生了轻微的位移，形成了轻隆起。0.5 s时随着铧式犁的不断切入，下层土壤逐渐被切开，随着犁面往上抬起，使表层和中层土随着下层土壤一同被抬起，使土壤开始堆积，在犁面挤压的作用下，土壤沿着犁面运动，横向土壤扰动范围加大，表层的土壤扰动范围增加最大，中层次之；0.8 s时铧式犁完全进入了土壤，由于犁面的挤压，越来越多的中层和下层土壤沿着犁面抬起堆积，表层土壤颗粒被挤压向侧边滚落，使中层和下层的土壤逐渐被翻垡到表层；1.2 s时铧式犁进行持续性的挖土和翻土过程，此时横向土壤扰动继续增大，表层土壤范围增加最大，此时中层土壤和深层土壤继续被抬起，中层和深层土壤被翻升暴露到表层，当土壤堆积高度超过了犁体高度时，少部分颗粒在犁体顶部被挤压掉落，大部分颗粒沿着犁面被抛翻向侧边运动。此时表层土壤一部分继续沿着犁体向前运动，另一部分被抛翻到侧边后，由于重力的作用和耕作过后的犁沟地形的影响回填形成了垄沟。

a，抛物线曲面纵向土壤扰动；b，椭圆曲面纵向土壤扰动；c，抛物线曲面横向土壤扰动；d，椭圆曲面横向土壤扰动；a中标注为高于表层的颗粒数；b中标注为抬升高度。a， Parabolic curved surface longitudinal soil disturbance； b， Elliptical curved surface longitudinal soil disturbance； c， Parabolic curved surface lateral soil disturbance； d， Elliptical curved surface lateral soil disturbance； The number of particles above the surface was marked in Fig. a； The height of uplift is marked in Fig.b.图9 不同类型犁体土壤扰动图Fig.9 Soil disturbance maps of different types of plough bodies

通过对比同一时刻不同类型犁体的土壤扰动状态和图中被抬升到表层土壤以上的颗粒数目可知，其土壤扰动过程大致相同。由图9对比可知，在0.5 s时由于犁体曲面不同，抛物线曲面抬升土壤颗粒到表层以上数目为19 697个，抬升高度为89 mm，而椭圆曲面抬升颗粒数为21 656个，抬升高度为97 mm，椭圆曲面抬升颗粒数明显多于抛物线曲面，通过对比横向扰动图即图9-c和d，椭圆曲面中颗粒堆积挤压的高度明显高于抛物线曲面中土壤堆积的高度；在0.8 s可明显看出，土壤沿着犁面运动抬升时土壤分布状态不同，相比较于抛物线曲面，椭圆曲面抬升土壤更加平直，抬升颗粒数目更多，抬升高度更高；1.2 s持续性耕作时，椭圆曲面从犁体上方掉落土壤明显多于抛物线曲面，说明椭圆曲面对于土壤具有更好的垂直方向上的抬升效果。

3.1.2 土壤运动状态分析

为了探究在耕作范围内，犁体在不同位置时土壤的运动状态，选取0.2、0.8、1.2 s 3个时刻，对不同时刻土壤颗粒的速度大小和方向进行比较。

由图10可知，当犁铧进入土壤的时候(0.2 s)，下层土壤在犁铧的作用下沿着犁铧表面向上方运动，随着下层的土壤被抬升，中层和表层的土壤在下层的挤压下也向上抬升，此时深层土壤颗粒颜色最深其运动速度最大，中层和表层次之；当犁体完全进入土壤时(0.8 s)，由图10-a、b可知，表层土壤受到挤压，纵向扰动范围最大，具有速度的颗粒数目最多，下层的最少，根据颗粒的颜色可知，相同颜色的颗粒沿着犁体曲面曲线分布，下层和中层土壤主要沿着犁体曲面抬升逐渐运动到土壤表层，上层土壤在受到中层和下层的土壤挤压并在重力的作用下，在向前和向上抬升的同时，少部分被堆积超过了犁体高度，掉落到犁沟中，大部分沿着犁体运动方向向两边滑落；在犁体处于持续性的耕作状态时(1.2 s)，对比于0.8 s时犁体前端的土壤颗粒速度大小和方向变化不大，但由于耕作过后土层形成了犁沟，随着越来越多的下层和中层土壤被抬起暴露于表层，堆积高度超过了犁体高度，部分下层和中层土壤颗粒从犁体顶端被挤压掉落于犁沟中，另一部分沿着犁翼向侧边运动，抛翻在土壤表面。通过对比不同时刻的土壤分布和颗粒速度及其方向可知，土壤颗粒的速度与颗粒和犁体之间的距离成反比，与沿着犁面抬升高度成正比，并且相同速度的颗粒分布与犁体曲面曲线基本一致。

a，抛物线曲面土壤颗粒速度大小分布状态；b，抛物线曲面土壤颗粒速度方向分布状态；c，椭圆曲面土壤颗粒速度大小分布状态；d，椭圆曲面土壤颗粒速度方向分布状态。a， Parabolic surface soil particle velocity distribution status； b， Parabolic surface soil particle velocity direction distribution status； c， Elliptical surface soil particle velocity distribution status； d， Elliptical surface soil particle velocity direction distribution status.图10 不同类型犁体土壤运动状态Fig.10 Soil movement state of different types of plough bodies

通过对比不同犁体曲面的土壤运动状态，在犁体完全进入土壤(0.8 s)和犁体进行持续性耕作(1.2 s)时，与犁体顶端接触的土壤颗粒速度明显大于其他位置的土壤颗粒速度，方向与其他颗粒方向保持一致，均为垂直于犁面方向，由于抛物线顶端部分趋近于直线，抛物线曲面顶端速度较高的土壤颗粒数量明显多于椭圆曲面。

3.2 不同类型铧式犁阻力分析

为了分析不同的导曲线类型的铧式犁阻力情况，其所受阻力如图11所示。由图11可知，0 s犁体开始运动，犁铧开始切入土壤，阻力开始增加，随着犁体的不断切入，土壤逐渐堆积在犁体曲面上，此时犁体阻力出现快速增长，0.8 s当犁体完全进入土壤后，犁体进入持续性的耕作状态，土壤从犁铧切入，沿着犁胸向犁翼堆积挤压移动，最后从犁翼被抛出掉落，此时犁体阻力增加速度减慢，由于切入土壤是连续的，顶端土壤被挤压掉落是不规律的周期状态，因此在持续性的耕作状态时犁体所受阻力表现为稳定的波动状态。

由图11可知，在同一工况下，不同导曲线的铧式犁受到的犁体阻力由于犁体曲面的不同，阻力大小也不一样，但其总体趋势均表现为先快速增加，接着缓慢增长一段时间，最后稳定地在一定范围内波动。犁体从运动开始阻力逐渐增加，1 s时，开始达到相对稳定阶段。在相对稳定阶段，不同导曲线的犁体阻力顺序为：圆弧线>摆线>抛物线>椭圆线。抛物线阻力为1 664.415 N，以圆弧线为导曲线的平均阻力最大，为1 878.196 N，以椭圆线为导曲线的平均阻力最小，为1 591.551 N。仿真结果表明：当犁体完全进入土壤耕作，达到相对稳定阶段时，不同导曲线的铧式犁所受阻力不同，其圆弧段阻力最大，比寻常的抛物线段阻力增大了12.8%，摆线段也增加了4.4%；椭圆段阻力最小，比抛物线段阻力减小了8.5%。由此可知，不同导曲线其曲线的不同，导致铧式犁所受阻力不同。

图11 犁体阻力图Fig.11 Plough resistance diagram

3.3 耕作速度对犁体阻力的影响

为了提高工作效率，研究犁体的耕作速度对于犁体阻力的影响，分别计算犁体在0.8、1.2、1.6 m·s-1耕作速度情况下的犁体阻力并对比抛物线铧式犁的减阻率，结果如图12所示。由图12可看出，随着耕作速度的增加，犁体阻力也随之增加。通过选择合适的耕作速度，有利于降低犁体阻力提升效率。

图12 速度与阻力关系Fig.12 Velocity and drag relationship

3.4 安装角对犁体阻力的影响

为了研究犁体安装角对犁体阻力的影响，分别选取安装角为22.5°、25.0°、27.5°、30.0°、32.5°，其阻力如图13所示。由图13可以看出，从安装角22.5°到25.0°犁体阻力随着角度增大而减小，在25.0°达到最小值，从25.0°到32.5°，犁体阻力随着角度的增加而增加。当安装角较小时，犁铧部分较为平缓，此时犁胸抬升，土壤高度较高，做功大，故此时犁体阻力较大；在25.0°时，犁体阻力达到最小值，当安装角继续增大时，犁铧越陡峭，土壤受到的挤压增大，土壤被抬升和翻垡所做的功增大，因此犁体阻力较为明显地增大。从仿真结果可知，在安装角为25.0°时为最优值，犁体阻力最小。

图13 安装角与阻力关系Fig.13 Relationship between installation Angle and resistance

3.5 铧刃角对犁体阻力的影响

为了研究犁体铧刃角对犁体阻力的影响，分别选取铧刃角为40°、42°、44°、46°、48°，其阻力如图14所示。从图14可以看出，犁体所受阻力随着铧刃角的增加而增加，趋近于线性关系。这是因为铧刃角具有切土和横向推土的作用，铧刃角越大，其在前进方向上投影面积越大，与土壤的接触面积越大，土壤在横向推土中的位移越大，犁体曲面阻力也随之增加。由仿真结果可以看出，在一定的可选取范围内，铧刃角越小，犁体阻力越小。

图14 铧刃角与阻力关系Fig.14 Relationship between blade angle and resistance

综上，通过对犁体的前进速度、安装角、铧刃角进行仿真分析。前进速度和铧刃角与犁体阻力趋近于线性关系，随着前进速度和铧刃角的增加，犁体阻力也随之增加；但随着安装角的增大，犁体阻力却呈现出先减小后增加的趋势。故在选择犁体工况的时候，要选择较小的速度；在合适的铧刃角范围内，尽可能选择较小的铧刃角；安装角尽可能地选择在25°左右，降低犁体耕作阻力。

4 实验验证

4.1 试验方案

本次试验在湖北省农业机械工程设计研究院土槽实验台中进行，通过改变导曲线设计不同的犁体曲面可以减小犁体阻力，提升耕作效率，对仿真的结果进行实验验证，根据实验的要求，按7∶10的比例分别制作出抛物线曲面和椭圆曲面的犁体模型如图15所示。实验土槽尺寸为3 m×1 m×0.5 m，实验设备为交流伺服电机、驱动器、HP数显式拉推力计。

a，抛物线曲面模型；b，椭圆曲面模型。a， Parabolic surface model； b， Elliptical surface model.图15 犁体模型Fig.15 The plough body model

4.2 实验结果

通过电机和驱动器拖动行走架从而带动犁体运动，在土槽中分别以0.4、0.8、1.2 m·s-13种前进速度分别对两个犁体进行阻力测量。为了使仿真结果更为准确，先进行空载实验测出空载运行时行走架的阻力大小为132.3 N，在进行犁体阻力实验时减去空载时阻力大小即为犁体所受阻力大小。本实验测多组数据再取平均值得到结果如图16所示。