基于SPH算法的深松铲土壤切削过程仿真及试验研究

曹中华，崔晋波，湛小梅，李亚丽，王圆明，余雪源，刘汶树，宋树民

(重庆市农业科学院，重庆 401329)

0 引言

为避免土壤耕作过程中出现的地力下降问题，进行土壤深松是迫切的需要[1]。土壤深松技术对于耕层植被和土粒的结构基本不会有影响，有利于降低能耗及成本[2]。深松铲是实现深松技术的核心部件，深松机在田间作业的主要能耗来自于深松铲克服土壤阻力所做的功[3]。目前，深松作业面临的主要问题是作业时的切削阻力较大，为达到理想的耕作效果，研究深松铲耕作过程中的切削阻力问题刻不容缓。

本文以MASCHIO PINOCCHIO系列全方位深松机的核心部件深松铲为研究对象，该深松铲采用凿型铲加双翼铲结构，最大耕深分别达到450、550、650mm；深松铲通过铲尖入土深松、铲翼浅松侧翻土壤、破土铲对土壤进行破碎。由于工作过程中，耕深较深，土壤切削阻力较大，能耗较大，所以为了降低能耗，进行土壤切削阻力的研究必不可少。然而，目前主要是针对深松铲的结构参数进行研究，建立深松铲受力数学模型，对于研究深松铲土壤切削过程的研究还比较少，不能准确描述深松铲工作过程中的受力状态，更不能准确地建立深松铲土壤切削系统模型[4-8]。

传统方法是通过田间试验数据的理论修正来进行农业机具的优化研究，这种方法是直观、可靠的[9-10]。

然而，田间试验是一个复杂的系统工程，需要收到人、物、时间、天气等条件的制约，因此试验周期长，而且试验结果不能有效反应机具—土壤相互作用的过程，具有一定的局限[11-14]。20世纪末以来，研究者逐渐将数值模拟方法引入到农机具与土壤相互作用的研究中，可对农机具和土壤相互作用过程进行仿真，直观反映土壤的运动轨迹和刀具的受力大小[15-17]。因此，本文基于光滑粒子流体动力学(SPH)方法，对深松铲土壤切削过程进行仿真，分析深松铲在切削土壤过程中的运动特征和力学行为，并基于耕整地移动式田间动态试验台进行田间试验以验证仿真的准确性，揭示深松铲工作过程中切削阻力的变化规律，为未来进行深松铲的减阻耐磨特性研究和优化设计提供理论依据，也为其他耕作部件土壤切削过程的研究提供理论基础。深松铲主要有机械式和振动式两种类型，本文研究机械式深松铲。

1 深松铲切土动力学仿真流程

深松铲土壤切削动力学仿真流程如图1所示。

2 深松铲土壤切削过程仿真研究

本文以MASCHIO深松机的深松铲为研究对象，为了使模拟更符合实际情况，运用SPH法对深松铲土壤切削过程进行仿真分析。

2.1 SPH数值分析方法

SPH是一种基于插值理论的无网格的纯拉格朗日方法，实质是将连续体离散成相互作用的粒子，由各粒子分担连续体的质量、速度等物理量，求解时对各粒子分开求解进而获得整体信息。其计算不需要依赖网格，因此在解决大变形问题上具有独到优势。深松铲土壤切削过程是一个大变形过程，为了使仿真过程更接近实际情况，本文基于SPH法来仿真深松铲土壤切削过程。

图1 深松铲切土动力学仿真流程Fig.1 The flow chart of soil cutting process by subsoiler

2.2 建立仿真模型

2.2.1 土壤SPH模型

本文所选土壤为川渝地区较为普遍的紫色土壤，参照土工试验方法标准和相关研究对其进行参数测定[18-21]，测得土壤干密度为ρ=1 760kg/m3，内聚力c=5.2e+4Pa，内摩擦角φ=0.223rad相关参数如表1所示。

表1 土壤参数Table1 Parameters of soil

续表1

在LS-Prepost中建立SPH土壤模型，土壤大小为700mm×450mm×500mm，各方向粒子层数为70、45、50，构建完成后的土壤SPH模型如图2所示。

图2 土壤SPH模型Fig.2 The SPH model of soil

2.2.2 深松铲有限元模型

深松铲结构复杂，在ANSYS中建模困难，本文运用三维扫描仪逆向进行深松铲点云构建，然后在CATIA中进行三维实体建模。运用Hypermesh软件进行深松铲网格划分，网格单元尺寸为5mm，共452 759个单元，深松铲有限元模型如图3所示。深松铲的主要参数如表2所示。

图3 深松铲有限元模型Fig.3 Finite element model of subsoiler表2 深松铲参数Table 2 Parameters of subsoiler

参数名称单位取值深松铲前进速度m/s1深松铲旋转速度rad/s0

续表2

2.2.3 土壤—深松铲有限元模型

在建立土壤—开沟刀具模型时，为有效避免切削土壤过程中土壤大变形问题造成的网格畸变和负体积，本文采用SPH法建立土壤模型，在LS-PREPOST中直接生成SPH粒子，建立700mm×450mm×500mm的SPH土壤模型。图4为土壤—深松铲的有限元模型。

图4 土壤—深松铲有限元模型Fig.4 The finite element model of soil-subsoiler

2.3 深松铲切削土壤过程仿真与分析

运用LS-Prepost软件查看LS-DYNA971求解结束后生成的d3plot文件。本文通过LS-Prepost软件可以直观地看到深松铲对土壤进行深松的过程，同时可以看到深松铲的运动情况，以及土壤的运动、破碎和变形情况，深松铲和土壤的受力情况；而且，可以通过分析深松铲切削土壤过程中的应变和应力云图，了解土壤应力和应变的分布情况，可以清晰地看到应力集中的位置；通过提取切削力时程曲线分析深松铲的切削阻力。

2.3.1 切削过程仿真

从图4和图5可以看出深松铲土壤切削仿真过程的效果。图4为深松铲土壤切削有限元模型，也是初始状态，深松铲与土壤不接触，SPH粒子不产生任何变化。图5中，t=0.2s时刻，深松铲开始切削土壤，铲尖最先与土壤接触，土壤发生破坏，深松铲继续切削土壤运动，铲翼也对土壤进行切削；t=0.4s时刻，深松铲的铲尖和铲翼同时进行土壤切削；t=0.6s时刻，整个深松铲完全进入土壤，此刻也基本处于切削稳定状态，深松铲对土壤扰动范围也最大；t=0.8s时刻，深松铲土壤切削过程在设置时间内结束，深松铲铲尖切出土壤，深松铲通过对土壤的抬升、扭曲和翻转，土壤被切削部分的通透性加强，底层土壤的紧实度增加，该部分SPH粒子排列不均匀，与田间作业情况相吻合。

(a) t=0.2s (b) t=0.4s

(c) t=0.6s (d) t=0.8s图5 深松铲切削土壤过程Fig.5 The process of soil cutting by subsoiler

2.3.2 土壤的等效应力、应变情况分析

图6和图7反应了深松铲切削土壤过程土壤的等效应力和等效应变效果。

图6中，t=0.2s时刻深松铲铲尖进入土壤对土壤进行深松，此刻土壤受到最大等效应力为1.361MPa；t=0.4s时刻深松铲铲尖和铲翼同时对土壤进行深松和浅松，此刻最大等效应力为1.469MPa；t=0.6s时刻整个深松铲都进入土壤，对土壤进行深松、浅松和破土，此刻最大等效应力为3.184MPa；t=0.8s时刻深松铲铲尖切除土壤，切削过程结束，此刻土壤最大等效应力为3.025MPa。

通过图7得出：土壤所受的最大等效应变和最大等效应力分布基本一致，都分布在铲尖、铲翼和破土铲与土壤接触位置。深松铲土壤切削整个过程，SPH土壤粒子的扰动也在不断变化，刚开始接触时土壤粒子扰动较小。随着深松铲的铲尖、铲翼和破土铲逐渐进入土壤， SPH土壤粒子的扰动加大， 且内侧扰动明显大于外侧，说明深松铲内侧SPH土壤的应力、应变大于外侧。

(a) t=0.2s

(b) t=0.4s

(c) t=0.6s

(d) t=0.8s图6 土壤等效应力分布图Fig.6 Equivalent stress distribution map of soil

(a) t=0.2s

(b)t=0.4s

(c)t=0.6s

(d)t=0.8s图7 土壤等效应变分布图Fig.7 Equivalent strain distribution map of soil

2.3.3 深松铲切削阻力分析

深松铲土壤切削过程中深松铲受到土壤的反作用力即为切削阻力。深松铲的切削阻力可以通过在LS-Prepost 中载入rcforc文件进行查看。图8(a)为x、y、z3个方向的切削阻力；图8(b)为切削阻力合力Resultant Force的时程曲线。

(a) x、y、z方向切削阻力曲线

(b) 深松铲切削阻力合力曲线图8 深松铲切削阻力时程曲线Fig.8 Curve of tillage resistance of subsoiler

由图8(a)可以得出：深松铲切削阻力由x、y、z3个方向的阻力合成，值为负说明所受力为阻力。A曲线的走势说明了在深松铲运动方向受到最大的切削阻力；B曲线的走势说明了在深松铲切削运动过程中，在竖直方向也会受到较大的切削阻力；C曲线的走势说明在深松铲切削土壤运动过程中在两侧会受到土壤的挤压，产生较小的土壤挤压力，该方向的力对切削阻力的影响不大。

由图8(b)可知：初始时刻，深松铲土壤未进行接触，切削阻力数值为0；紧接着，深松铲铲尖首先切入土壤，切削阻力增大，随着深松铲铲翼和破土铲逐渐切入，切削阻力逐渐增大；0.6s时刻，深松铲全部切入土壤，切削阻力达到最大值；0.8s时刻深松铲离开土壤，切削阻力逐渐降低。结合图6和图7土壤等效应力和等效应变变化情况，证明切削阻力和土壤等效应力和等效应变存在关联性。

0.5s时刻深松铲完全进入土壤，0.8s时刻深松铲开始切出土壤，提取图8中0.5～0.8s切削阻力数据，在Excel中进行计算，得出平均切削阻力为3.65kN。

3 田间试验验证

3.1 试验条件

试验在重庆市农业科学院农业机械研究所农机中试基地进行，以履带式拖拉机作为牵引动力。试验过程中，拖拉机前进速度为3.6km/h，耕深标定为450mm。深松铲切削阻力测试系统采用耕整地移动式田间动态试验台。耕作土壤为川渝地区较为普遍的紫色土壤，其物理特性参数如表3所示。

表3 试验区土壤物理特性参数Table 3 Soil physical parameters in the test area

该试验台由轻型履带式拖拉机、田间作业机具测试系统及相应测试传感器组成，主要适用于南方丘陵山区较平坦、开阔地块，可用于旋耕机、深松机等农机具在不同土壤工况(不同坚实度、不同含水率)下的田间作业性能试验，可直接测试深松铲工作过程中x、y、z3个方向的切削阻力。该系统由服务器端平台和远端平台两大部分组成。其中，远端平台能够实现传感器数据的实时采集与处理、试验参数配置、试验数据文件存储于上传及采集模块故障检测等多种功能。服务器端平台由计算机及相应测试应用软件组成，可以实现实时查看远端平台测试试验运行状态、试验结束后的数据文件远程存储等功能。

试验之前应进行传感器检查、参数设置、速度标定、耕深标定和模块检测等事项。田间试验过程如图9所示。

图9 田间试验过程Fig.9 Field test process

其他试验仪器包括TJSD-750坚实度仪、TZS-I水分仪、环刀。

该系统通过实时采集并处理农机具在田间作业时的工作性能参数和拖拉机常规性能参数，可用于研究分析农机具在不同作业工况下的作业性能情况及其与拖拉机能耗之间的关系，为农机具的研发、改进提供可靠的科研数据支撑。

3.2试验内容与结果分析

本文通过耕整地移动式田间动态试验台进行深松铲切削土壤田间试验，试验过程中深松铲耕作深度为450mm，土壤坚实度为1 021kPa，土壤含水率为21%，土壤容积密度为1 760kg/m3。同时，运用远端田间作业机具测试系统采集耕作过程中深松铲切削阻力数据，提取20组数据(见表4)，求取切削阻力平均值为3.542kN，与仿真结果相比误差为2.9%，验证了基于SPH法进行深松铲切削土壤过程的仿真是可行的。

表4 深松铲切削土壤试验前进速度-切削阻力表Table 4 Forward speed - cutting resistanceTable of soil cutting process by subsoiler

4 结论与展望

4.1 结论

1)基于SPH算法研究了深松铲土壤切削过程仿真，建立了深松铲-土壤有限元模型，直观展示了深松铲土壤切削过程，与实际情况基本吻合。

2)深松铲土壤切削仿真表明：不同时刻、不同方向等效应力变化规律不同。切削过程中，最大等效应力为3.184MPa，土壤所受的最大等效应变和最大等效应力分布基本一致，都分布在铲尖、铲翼和破土铲与土壤接触位置。

3)在耕深为450mm、深松铲前进速度为3.6km/h、土壤坚实度为1 021kPa、土壤含水率为21%情况下，进行深松铲土壤切削仿真，在深松铲完全进入土壤时求解土壤切削阻力平均值为3.65kN。

4)基于耕整地移动式田间动态试验台进行田间试验，测试深松铲切削阻力值。在耕作稳定后提取20组数据进行切削阻力平均值求解为3.542kN，与仿真结果对比误差为3.05%，验证了基于SPH法进行深松铲切削土壤过程的仿真是可行的。

4.2 展望

1)通过对深松铲的结构参数和运动参数的研究，结合数值模拟方法和田间试验，进行深松铲的优化设计，以达到进一步减阻降耗的目的。

2)基于流固耦合理论，采用数值模拟方法的耦合来进行深松铲切削土壤过程的仿真分析。