基于SPH法的林地深松刀具仿真及试验

王恩民，李立君，汤刚车，欧阳益斌

(中南林业科技大学 机电工程学院，长沙 410004)

0 引言

林业机械的发展是机械化水平提高的重要保证，尤以深松机为重。目前，保护性耕作受到越来越多的重视[1-3]，土壤深松技术作为保护性耕作关键技术之一更是一种新型的土壤耕作方法。土壤深松技术是指在土壤不被翻转的情况下进行土壤疏松，对土壤中固、液、气3种成分比例进行调节，考虑到对林地作物情况影响，土壤耕深在200～300mm之间有利于根系生长，促进作物增产增收，对于保水保墒也发挥着重要作用[4]。

目前，对深松机研究多集中在旋耕刀或深松铲的研究，以及基于这两种刀具的受力分析及振动减阻原理研究，研究取得一系列成果[5-7]；但是，耕深不稳定、阻力大、功耗大、深松后土壤硬度及平整度达不到要求等一系列问题仍没有得到很好的解决。针对这一现状，设计了一种基于旋转切削原理的立式深松刀具，确定了刀具工作参数，并采用SPH法[8-11]建立刀具切削土壤模型，开展试验验证仿真可信度，为后续优化提供技术依据。

1 深松工作参数确定

1.1 深松刀结构设计

深松刀由法兰，加强肋、深松刀轴及深松刀片组成，结构如图1所示。

1.法兰 2.加强肋 3.深松刀轴 4.深松刀片图1 深松刀结构图Fig.1 Deep loosening tool structure

深松刀作业时，依靠刀轴带动刀片，刀片在刀轴上呈等螺距分布，工作头升降使刀片可以对不同深度的土壤进行深松，且土壤不被翻转，满足保护性耕作的要求。

切削过程中，深松刀切削阻力与切削角正切成正比[12]，但一定的切削角才能使土壤微粒向上运输，不至于堵塞造成阻力变大，同时保证切土效果。切削角取值范围为15°～25°。

1.2 深松刀运动分析

深松过程中，假设底盘以固定速度前进，深松刀轴以固定速度旋转，深松刀为匀速直线运动与匀速旋转运动的复合。其运动轨迹可用余摆线方程表示，即

式中R—深松刀旋转半径(m)；

ω—深松刀旋转速度(rad/s)；

v—深松刀前进速度(m/s)。

1.3 深松过程最低旋转速度确定

深松刀工作过程，为了使土壤微粒能随着深松刀旋转向上运输，深松刀旋转速度必须大于一个临界值。

1.3.1 土壤微粒速度分析

土壤微粒速度分析如图2所示。

(a) xoz平面内土壤微粒速度分析图

(b) xoy平面内土壤微粒速度分析图图2 土壤微粒速度分析图Fig.2 Analysis of soil particle velocity

设深松刀角速度为ω，前进速度为v0，对土壤微粒在刀边缘进行速度分布。

深松刀在旋转时，土壤微粒也在以ve的速度随刀旋转的同时，以vr的速度向上运动，va为vr与ve的合速度，与水平面夹角为β，α为刀片的倾角,则

(1)

水平方向上

vt=ve-vrx

(2)

vt=va·cosβ

(3)

ve=ω·R

(4)

vrx=vz·cotα

(5)

vz=va·sinβ

(6)

因为vt=ωa·R=va·cosβ，则

(7)

将式(3)～式(7)代入式(2)中得

(8)

式中ω—深松刀旋转速度；

α—深松刀片倾角；

β—土壤微粒绝对速度与水平面的夹角；

ωa—土壤微粒绝对角速度。

1.3.2 土壤微粒受力分析

土壤微粒受力分析如图3所示。

由式(3)得

dF1=N1·f1

(9)

dF2=N2·f2

(10)

式中dF1—深松刀片与土壤微粒间摩擦力，与vr方向相反；

dF2—沟壁与土壤微粒间摩擦力；

N1—刀片对土壤微粒支撑力；

N2—沟壁对土壤微粒支撑力；

f1—土壤与螺旋面摩擦因数；

f2—土壤内摩擦因数。

据图3可知：在N1的方向上有

N1=dF2·cos(90°-α-β)+dm·(g+az)cosα

(11)

式中dm—土壤微粒的质量；

az—土壤微粒的径向加速度。

在垂直FN1方向上有

(12)

(13)

式中FN1—dF1与N1的合力；

φ1—FN1与N1的夹角；

Fe—土壤微粒所受离心力；

N2—沟壁对土壤微粒的支持力。

(a) xoz平面内土壤微粒受力分析图

(b) xoy平面内土壤微粒受力分析图图3 土壤微粒受力分析图Fig.3 Analysis of soil particles by force

将式(10)、式(13)代入式(12)中得

代入式(8)得

令β=0得到深松刀需满足的最低角速度为

式中ω0—满足土粒升运最低角速度；

g—重力加速度；

R—深松刀片外径；

α—深松刀片倾角；

φ1—土壤质点与螺旋面摩擦角。

计算得：ω0=33.9rad/s。

1.4 深松前进速度确定

配套深松机前进速度如表1所示。

表1 前进速度表Table 1 Forward speedTable m/s

根据深松作业效率要求得深松机前进速度需大于0.5m/s。

2 刀具切削土壤过程模拟

2.1 SPH算法基本原理

刀具切削土壤属于高速切削，采用有限元网格法进行仿真易造成网格畸变导致结果不收敛。SPH算法是一种无网格算法，将土壤离散成一系列具有相互作用且带有质量和速度等物理量的点集合，通过粒子模拟流体运动规律[13]。其核心是插值理论，要求解各个质点在任意时刻能量和速度，需引入核函数；场函数经过核函数“光滑化”，再在整个求解域上积分，便得到了表征各个质点的动力学特征函数[14-16]，即

其中，W为核函数，使用辅助函数θ定义为

其中，v为空间维数；h为光滑长度。

2.2 刀具切削土壤模型建立

模型的前处理及数据的后处理均在Ls-Prepost中进行，仿真流程如图4所示。

土壤SPH模型材料采用MAT147，并采用了Drucke_Prager塑性模型，以Mohr_Cloulomb屈服准则进行修正[17]，考虑了粘塑性和切削应变失效等方面的因素，使土壤模型更加符合实际，使仿真结果更加可靠。土壤模型参数如下：

土壤密度/kg·m-3： 2.35

体积模量/Pa： 4.6×107

剪切模量/Pa：
1.86×107

土壤水分含量/%： 3.4

粘塑性参数：
1.1

土壤比重： 2.79

在Ls-Prepost前处理中，在土壤底部和两侧面添加固定约束，耕深为230mm，接触静摩擦因数为0.2，动摩擦因数为0.18，深松刀具与土壤接触方式为点面接触，最后通过定义曲线的方式来定义深松刀具初始速度。

图4 仿真流程图Fig.4 Simulation process

根据参数分析，切削角取值范围为15°～25°，刀具最低旋转速度为33.9rad/s，最低前进速度为0.5m/s，以牵引力为评价指标，因素水平表如表2所示。

表2 正交试验因素水平表Table 2 Simulation test factor level

2.3 仿真结果分析

前处理K文件在Ls-dyna中求解完毕后，深松过程牵引力随时间变化，如图5所示。求得其平均牵引力为1.092 3kN。

按表2进行9组仿真，试验安排及结果如表3所示。

由表3中极差分析可得：对牵引力影响显著性排序为切削角、旋转速度、前进速度，同时得出最优组合为A1B3C3，即切削角为15°、旋转速度为58rad/s、前进速度为0.7m/s时，深松牵引力最小。

图5 牵引力变化图Fig.5 Traction tendency chart表3 试验安排及结果Table 3 Experiment arrangement and result

试验序号A切削角/(°)B旋转速度/rad·s-1C前进速度/m·s-1牵引力/kN11111.092321221.176131330.855642121.295752231.164162311.095473111.356483231.264593321.1654K13.1243.7443.544K23.5553.6053.637K33.7863.1163.284k11.0411.2481.181k21.1851.2021.212k31.2621.0391.095R0.2210.2090.118最优组合A1B3C3

3 深松刀牵引力试验

3.1 试验目的及试验地概况

为验证仿真结果可信度，为后续优化提供依据，开展深松刀林地作业试验，验证仿真结果的正确性。

深松试验于2017年夏季在湖南省某林地进行，该地土壤常年未进行深松作业，含水量较低，地面硬度比较大。为了使试验场地土壤的含水量等情况基本一致，选择同一坡面及等高线平面进行试验。为了完成牵引力的测试，使用1台履带式除草机和深松机配合作业。

试验仪器设备为F0024艾德堡推拉力计，量程2 000kg，示值误差±0.5%；笔记本电脑；皮尺；土壤硬度计；秒表；标杆。

3.2 试验设计

3.2.1 性能指标

开展以刀具切削角、刀具旋转速度及刀具前进速度三水平三因素的正交试验，以深松过程牵引力大小作为评价指标，试验因素水平表如表2所示。

3.2.2 试验步骤

试验过程中，采用牵引法测试牵引力大小，在履带式除草机和深松机之间安装拉力传感器，试验时，牵引车辆从静止到匀速的过程尽可能缓慢，以减少冲击载荷影响；记录下在相同作业长度内深松机前进阻力及安装深松刀具后深松牵引力，二者相减即为深松刀具的牵引力。试验过程如图6所示。

(a) 试验前硬度测量 (b) 传感器安装调试图6 试验过程Fig.6 Experiment process

3.3 试验结果及分析

试验以牵引力为试验指标，传感器软件实时记录牵引力数值，如图7所示。

求得深松牵引力平均值为2.63kN，以及深松机前进阻力1.48kN，二者相减即为深松刀牵引力1.15kN。正交试验的安排及结果如表4所示。

3.4 数据分析及对比

1)得到与仿真结果相同的最优组合A1B3C3，对牵引力影响显著性排序为切削角、旋转速度及前进速度。

2)仿真误差按数值表示方法分为绝对误差和相对误差，即

绝对误差=仿真值-实测值

试验与仿真结果对比如表5所示。

图7 牵引力变化图Fig.7 Traction tendency chart表4 试验安排及结果Table 4 Experimental arrangement and results

试验序号A切削角/(°)B旋转速度/rad·s-1C前进速度/m·s-1牵引力/kN11111.1521221.2731330.8142121.3852231.3162311.0873111.4183231.3493321.29K13.233.943.64K23.773.923.94K34.043.183.46k11.081.311.21k21.261.311.31k31.351.061.15R0.270.250.16最优组合A1B3C3

表5 试验与仿真结果对比Table 5 Data analysis and contrastTable

续表5

4 结论

1)设计了一种立式深松刀具，对深松过程工作参数进行分析计算，得出了切削角范围为15°～25°，最低旋转速度为33.9rad/s，最低前进速度为0.5m/s。

2)分析切削土壤特点，基于SPH法对刀具切削土壤进行了仿真，采用3因素3水平正交试验表模拟仿真过程，得到切削角对阻力数值影响最大，且当切削角为15°、旋转速度为58rad/s、刀具前进速度为0.7m/s时，深松阻力最小。

3)开展了实地试验，采用牵引法收集牵引力数值，与仿真数据进行对比，误差在可接受的范围内，为下一步刀具优化提供技了术依据。