李霞+张思远+王维新+汤明军+郭金龙

摘要： 为降低耕作阻力，揭示深松铲与土壤之间的关系特性，根据深松铲切削土壤的工作特点，利用LS-DYNA模拟分析深松铲切削土壤的过程，获得切削土壤的应力变化规律。结果表明，当深松铲以0.49 m/s的初速度切削、深松深度为250 mm时，单个深松铲的最大切削阻力为2 230 N，土壤在1.75 s发生崩裂，达到深松效果。通过试验测试，单个深松铲受力约为2 332.5 N，与仿真误差为4.5%，验证了仿真的合理性，仿真数据结果可以为实际深松过程提供参考。

关键词： 土壤；切削；有限元法；深松铲；LS-DYNA模拟分析；阻力仿真

中图分类号： S222.19 文献标志码： A

文章编号：1002-1302（2017）22-0257-03

深松铲是深松的核心部件，通过深松可以改善土壤三相比，从而促进植物的生长[1]。然而，深松过程中有特别大的阻力。因此，须要探明在深松土壤切削过程中阻力的变化规律，达到减少阻力的目的。

在以往关于土壤切削的研究中，主要采用试验方法，但试验具有受外界环境因素影响大且经费多、周期长等缺点[2]。近50年来，随着数值分析模拟技术的发展，为研究土壤切削提供了一种新的方法。数学模拟切削土壤方法主要包括传统分析方法、有限单元法（finite element method，简称FEM）、离散单元法（distinct element method，简称DEM）[3]。近年来，多利用有限元法研究土壤切削问题，通过模拟仿真，可以方便清楚地观察到土壤和工作部件的情况，例如土壤的运移状况、破坏位置和工件部件各部分的受力状态。这些模拟仿真数据为耕作部件的设计优化提供了理论依据。

针对土壤非线性不均匀的特性，本研究采用动力分析方法探讨深松铲对土壤的切削过程，通过显式非线性动力分析软件ANSYS/LS_DYNA来进行深松铲切削土壤的数值模拟，建立深松铲深松工作过程的有限元仿真模型，对深松过程进行三维数值仿真研究。

1 有限元建模

1.1 土壤模型

ANSYS/LS-DYNA 12.1程序材料库中193号材料（MAT_DRUCKER_PRAGER）模型是一种弹塑性本构模型。由于这种材料本构模型能准确地体现土壤特性，简单实用而且能够得到准确的结果，所以被广泛地用于模拟岩土材料中，具体参数取值见表1[4]。Mohr-Coulomb屈服条件通常被使用在土壤切削仿真中，但其有不可忽略的缺点：它的屈服面在平面上表现为特殊的等边六边形（具有棱角奇异性的不等角），这大大增加了数值计算的困难。与Mohr-Coulomb屈服条件相比，Drucker-Prager屈服准则不仅计算简单，而且在选择恰当材料常数的前提下，就可以寻找到与之匹配的Coulomb模型。基于以上比较，本研究采用改进扩展后的Drucker-Prager准则[5]，即：

式中：C为偏应力第3应力不变量；I为单轴拉伸屈服应力与单轴压缩屈服应力比值；F1为凝聚力；为材料摩擦角；F2为Mises等效应力；F3为平均压应力；σ1、σ2、σ3为3个方向的主应力；t为应力。

选定土壤的实体模型为800×800×800的立方体，由于模型比较简单，所以直接在ANSYS中生成。土壤选用solid164单元，这种单元常用于三维显示结构实体，由8个节点构成而且支持所有许可的非线性特性，对土壤采用智能网格划分，选择土壤模型的12条边，设置网格数为20，单元形状采用六面体8个节点，划分为9 261个节点，8 000个单元[6]。

1.2 深松铲模型

ANSYS在构建复杂实体模型方面比较薄弱，所以本研究利用SolidWorks构建深松铲模型并通过接口导入到ANSYS中[7]。深松铲采用 65Mn材料，其密度为7.8×103 kg/m3，弹性模量为2.10万MPa，泊松比为0.3。由于深松铲是一个整体，所以采用计算精度较高的四面体单元solid45单元类型，采用智能自由网格划分方法对深松铲模型进行网格划分，设定尺寸级别为5级，为减小应力集中的影响，在深松铲铲刃部分进行网格细化，深松铲划分为12 724个节点， 19 046 个单元。在建立土壤-深松铲数值模拟模型时，对模型进行网格划分后如图1所示。

1.3 创建单元

为了给后面定义接触类型和施加边界条件作基准，要先创建单元（简称PART）。PARTS是1个单元集合体，它是由相同单元类型、材料号和实常数组成的。若要生成PARTS，必须要求选择所有实体。本案例生成的PARTS结果如图2所示，其中1代表深松铲模型，2代表土壤模型。网格划分数量相当是能够计算的必要条件，从图2可以看出，深松铲模型和土壤模型网格数量相当，这为后续计算提供了可能。

1.4 定义接触

ANSYS/LS_DYNA12.1接触类型多种多样，而且与接触类型对应的还有许多可选择接触的分析方式，所以在描述深松的过程中，要选择恰当的类型和参数。深松过程是利用深松铲来使坚硬的土壤破碎，在计算機模拟系统中，这种工作状况就是土壤材料单元的失效。在破碎土壤后，土壤材料的外部单元常常无法与深松铲接触，为了避免这一情况的发生，将深松铲与土壤之间的接触类型定义为面与面接触中的侵蚀接触（eroding surface-to-surface contact，简称ESTS）。深松铲与土壤之间静摩擦系数为0.60，动摩擦系数为0.32。

1.5 施加边界条件和载荷

由于实际土壤立方体材料尺寸足够大，所以为保证土壤结构在整个侵彻过程中不发生移动，定下表面为固定表面，施加全约束，土壤顶部为自由边界，其他面均设想为无限大区域，采用无反射边界约束[8]。为了节约计算时间，将土壤的左右侧面也固定起来，为土壤前后端面施加无反射边界条件。由于深松铲采用的是刚体材料，只须施加在其质心上。对深松施加Z方向的位移约束，限制深松铲的左右晃动，并约束深松铲3个方向的转动自由度。endprint

2 结果与分析

2.1 求解

上述操作完成后，進入求解分析阶段。首先设置求解时间，若设置值较大，则会延长计算时间，而且一般物理过程的计算终止时间多为毫秒级，因此，设置求解时间为0.4 s，设置输出步数为20步。设置完成后，选择输出文件类型，输出关键字文件递交给LS-DYNA求解，并通过后处理程序打开结果文件，进行分析。

2.2 结果分析

载荷大小影响土壤状态。随着外力的变化，土壤会依次经过线弹性—塑性—屈服—硬化这几个状态，在土壤硬化后继续加载，则土壤开始失效，直至最终破裂[9]。图3描述了深松深度为250 mm的土壤深松过程。在深松工作过程中，首先是深松铲铲尖接触土壤，由于受到深松铲铲尖的挤压和剪切作用，土壤发生变形和破坏。随着深松工作的继续，深松铲前进并且产生一定振动，此时深松铲铲柄也和土壤相互接触，导致土壤被剪切挤压的越发明显，受到的力也显著加大，被破坏的面积和进给方向的深度也增加，后面还未与深松铲接触的土壤由于受到前面被深松铲挤压破坏后土壤的挤压也逐渐开始变形，最终实现整个深松土壤的变形破坏，达到深松的目的[10]。观察1个周期深松铲对土壤的深松过程，整个过程结果输出20步。在1.25 s时，深松铲开始和土壤接触，并以049 m/s的初始速度开始深松土壤，由于受到深松铲的剪切力作用，土壤状态开始发生改变；在1.5 s时，因受深松铲的剪切力和挤压力作用，土壤在沿深松铲前刀面方向上发生变形失效，此时切削阻力最大，达到2 230 N。在1.75 s时，土壤发生崩裂。深松铲继续前进，后面的土壤也被深松。由此可见，通过仿真模拟可以较为直观地观察出深松铲深松工作过程中土壤的状态变化。

3 试验验证

3.1 试验准备

为验证仿真结果的准确性，在新疆石河子市石河子大学教学试验场二连进行了试验。在试验田内选择1块地势比较平坦的未耕作土地，并在试验前按照振动深松机组试验方案的要求对试验区进行规划。确定长60 m（两侧各自留出 10 m 长的调整区，用于调整拖拉机的工作状态）、宽50 m的区域为振动深松机样机的试验场地，并在试验场地的两侧插上标杆，保证拖拉机在规划的区域内行驶。试验设备主要有振动深松机样机1台，雷沃M800-D型拖拉机1台，黑龙江省农业机械工程科学研究院研制的田间机械动力学参数遥测仪1套。

3.2 试验内容

将振动深松机样机与拖拉机3点悬挂连接，拖拉机的后置动力输出轴通过万向联轴器与样机的动力输入轴相连。拖拉机以慢二挡前进（秒表测得行驶40 m耗时82 s），当拖拉机进入测定区域时，记录遥测仪数据，并且重复试验3次。

3.3 试验结果

当拖拉机以慢二挡前进（秒表测得行驶40 m耗时82 s）时，重复3次试验，采集的数据如图4所示。

由图4可见，在深松作业过程中，3次试验的平均阻力分别为9.35、9.74、8.9 kN，则试验深松机阻力平均值为9.33 kN；在试验过程中，深松机共有4铲同时工作，且工况相同，所以单铲的力约为F=9.33/4=2.33 kN，试验力与仿真得到的力的误差为4.5%，因此利用仿真的方法来研究深松铲的土壤切削是合理可行的，同时也验证了该切削模型的正确性。

4 结论

通过建立深松铲-土壤的切削有限元模型，利用LS-DYNA进行数值分析，结果显示，当深松铲以0.56 m/s的初速度切削、深松深度250 mm时，单个深松铲最大切削阻力为2 230 N，土壤在1.75 s开始发生崩裂，达到深松效果，并通过试验验证了其准确性。仿真结果为优化深松铲结构、提高拖拉机功率的利用率提供了理论依据。

参考文献：

[1] 李洪文，陈君达，李问盈. 保护性耕作条件下深松技术研究[J]. 农业机械学报，2000，31（6）：42-45.

[2]王 敏，李萍萍，王纪章，等. 基于LS-DYNA茶园立式旋耕刀具土壤切削过程模拟[J]. 江苏农业科学，2015，43（9）：422-425.

[3]Wegner J L，Yao M M，Zhang X.Dynamic wave-soil-structure interaction analysis in the time domain[J]. Computers&Structures，2005，83（7）：2206-2214.

[4]姚仰平，张丙印，俊 高. 土的基本特性、本构关系及数值模拟研究综述[J]. 土木工程学报，2012，45（3）：127-150.

[5]苏翠侠，王燕群，蔡宗熙，等. 盾构刀盘掘进载荷的数值模拟[J]. 天津大学学报，2011，44（6）：522-528.

[6] 张 建，王 颖，王淑红. 基于LS-dyna的油菜靴鞋式开沟器工作阻力动态仿真[J]. 甘肃农业大学学报，2015，50（4）：156-159.

[7]蒋红旗，王繁生. 起重机吊臂结构有限元模态分析[J]. 农业机械学报，2006，37（3）：20-22.

[8]张冬梅，张 亚. 基于ANSY/SLS-DYNA头部形状对钻地斜侵彻土壤影响的数值分析[J]. 现代机械，2010（5）：21-22.

[9]刘 辉. 基于LS-DYNA的旱地旋耕刀工作机理研究[D]. 重庆：西南大学，2012.

[10] 辛丽丽，李传琦，梁继辉，等. 考虑分段土壤作用力的振动减阻分析[J]. 农业机械学报，2014，45（2）：136-140.endprint