基于离散元法的茶园仿生铲减阻性能研究

姜嘉胤，董春旺，倪益华，徐家俊，李杨*，马蓉*

1. 浙江农林大学光机电工程学院，浙江 杭州 311300；2. 中国农业科学院茶叶研究所，浙江 杭州 310008；3. 浙江川崎茶业机械有限公司，浙江 杭州 311121

耕作具有清除杂草、提高土壤肥力等优点，但也是茶园管理中劳动力最为密集、机械化效率较低的环节之一。近年来，我国农业劳动力老龄化趋势明显，招工难、用工贵已成为制约浙江茶产业发展的瓶颈。面对茶产业生产“机器换人”发展的迫切需求[1]，急需针对浙江丘陵山区茶园耕作机具缺乏问题，研制自走式耕作机，提升茶园管理作业现代机械化水平。然而茶园多分布于丘陵地区，地势复杂且茶行间的间距较小，为防止茶树受到损伤，通常耕作设备以轻便型为主，因此动力也非常有限。此外，茶园土壤一般板结较为严重，耕作阻力较大，且杂草、茶树枝条较多，传统旋耕刀和深松趾易发生缠绕现象，极易对茶树的树根造成损伤，影响茶树生长，因此需要设计更加适合茶园的耕作机构来应对茶园复杂的耕作条件[2]。

仿生学的研究对解决工程实际中阻力大的问题具有重要作用[3]。文立阁[4]以蝼蛄挖掘爪趾为仿生原型，设计了一种仿生灭茬刀，经试验测得减阻率为 5.31%～10.02%。李默[5]仿螳螂前足结构，俞杰[6]仿家兔爪趾几何外形等设计了旋转耕作刀具，并通过试验验证了仿生刀具的减阻性能。李建桥等[7]模仿蜣螂体表触土部位结构特征设计了凸包型仿生犁，通过土槽试验测得其减阻率为6.6%～12.7%。Sun等[8]模仿熊爪设计了开沟机构，张思博[9]模仿蝼蛄爪趾形态设计了挖掘机斗齿，张智泓等[10]以砂鱼蜥头部为仿生原型，对深松趾尖进行优化，通过试验验证均能有效减小耕作阻力。此外，一些土壤洞穴动物的爪趾具有优异的挖掘能力，例如鼹鼠一个晚上可以挖掘出91 m长的洞穴[11]，鼢鼠平均挖掘距离为100 mm·min-1左右[12]等。王洪昌[13]基于鼢鼠爪趾内轮趾纵向曲线特征方程设计了一种苗间仿生除草趾，经土槽试验证明其在低速条件下除草的牵引阻力较小。汲文峰[14-15]和Tong等[16]基于鼹鼠前足爪趾轮趾曲线特征设计了仿生型旋耕-碎茬刀，试验表明在土壤破碎率和秸秆破碎率均满足我国国家标准要求的情况下，功耗小于传统刀片。上述研究表明，针对耕作机构的仿生设计，特别是仿土壤洞穴动物爪趾的机构设计，对部件减阻降耗具有重要借鉴意义。

因此，针对茶园土壤板结严重、行距小、耕作阻力大且易使耕作机构发生缠绕等问题，本研究对茶园土壤的参数进行测量与标定，利用鼹鼠爪趾优秀的挖掘能力，结合鼹鼠的挖掘动作，采用离散元的方法设计并优化出具有优秀减阻效果的耕作机构，旨在为解决茶园行间土壤耕作过程中阻力过大的问题提供方案和设计依据。

1 材料与方法

1.1 茶园土壤参数测算及标定

我国茶园管理方式大多为粗放式管理，对茶行间的土壤耕作较少，土壤板结较为严重，土壤的孔隙率较小。因此若使用普通的田间土壤参数进行仿真，可能导致误差较大，需对部分参数重新进行测定，对一些不易直接测得的参数使用响应曲面的方法进行标定。

1.1.1 茶园土壤物理参数测算

土壤是各种矿物质的几何体，在天然状态下，包含固体颗粒、水和空气，三者之间的相互作用及比例关系，反映出土壤的物理性质和物理状态，这些指标又都与土壤的力学性质相关[17]，因此测量茶园的土壤参数是研究茶园耕作刀具的基础。本研究检测土壤取自浙江农林大学东湖校区内的茶园，在茶园行间随机选取5个点，向下挖约5 cm除去表层土壤后使用环刀取样。根据土工试验方法标准[18]对土壤的含水率、密度等参数进行测量。

取15～30 g土壤样品放入铝盒中，使用电子秤进行称重，得到湿土质量m0，将铝盒与土壤放入干燥箱中，在 105℃下烘干，12 h后取出称重，得到干土质量md。

式中，θ0为土壤含水率；ρ0为土壤湿密度，g·cm-3；ρd为土壤干密度，g·cm-3；V为环刀内体积，cm3。

由公式（1）～（3）计算可得，土壤含水率θ0为 21.38%，ρ0土壤湿密度为 1.83 g·cm-3，干密度ρd为 1.45 g·cm-3。

式中，m1和m2分别为土壤颗粒和水分质量，g；ρ1和ρ2分别为土壤和水分的密度，g·cm-3；V1和V2分别为土壤颗粒和水分的体积，cm3；r1和r2分别为土壤颗粒半径与土壤颗粒黏结半径，cm。

对环刀法取得的土壤样品进行直剪试验，测得茶园土壤内摩擦角Ø为29.62°，黏聚力c为 13.84 kPa。通过文献资料可知[19]，土壤单位面积法向刚度σ、单位面积切向刚度τ和土壤泊松比μ可根据直剪试验得到的内摩擦角和黏聚力求得，计算公式如下：

式中，c为土壤黏聚力，kPa；Ø为土壤内摩擦角，（°）；K0为侧压系数；μ为土壤泊松比。

由试验测得的数据可知，板结土壤的黏聚力无法忽略，而在对非板结（松散）土壤进行相关研究时通常认为土壤黏聚力可忽略不计[20]，更大的黏聚力会使机构耕作时所受的阻力变大，导致小型机械耕作困难，因此需要对耕作机构进行仿生减阻的研究，使其更易打破土壤颗粒间的黏聚力。

1.1.2 茶园土壤参数标定

目前，土壤形状、密度、泊松比和剪切模量等本征参数可以通过仪器进行测量或从文献获取。土壤的接触参数很难进行实际测量，可借助离散元仿真优化标定来获得[22]。

土壤堆积角试验的前期处理阶段，需将土壤样品破碎并自然风干，风干后的土壤已从板结状态变为相对松散状态，土壤中的水分含量也大大降低，土壤颗粒之间的黏结键数量大大减少，使用 Hertz-Mindlin（noslip）模型，不考虑土壤颗粒之间的相互黏结情况，更适用于堆积角试验的标定。参数标定仿真中所使用的泊松比为0.336，剪切模量为1×106Pa[21]，土壤颗粒密度为 2 500 kg·m-3[21]。

通过查阅文献，得到土壤-土壤间的静摩擦系数、滚动摩擦系数和碰撞恢复系数分别在0.3～0.7、0.14～0.4 和 0.2～0.6 范围[22-29]。

采用二次回归正交旋转中心组合优化试验方法，试验因素水平如表1所示，通过EDEM 仿真得到各个水平下土壤堆积角的结果如表2所示，应用Design-Expert 8.0进行试验数据处理分析。

表1 仿真试验因素水平表Table 1 Simulation test factors

表2 土壤堆积试验仿真设计及结果Table 2 Design and results of soil accumulation simulation test

如表3所示，响应面回归模型中土壤堆积角模型P＜0.01，表明回归模型显著；失拟项P＞0.05，说明失拟不显著。决定系数R2为0.991，R2Adj为0.983，都接近于1，变异系数CV为2.64%，表明模型的回归方程可行度高。精确度为44.472，大于4，说明回归模型的精度较好。回归方程：

表3 回归模型的方差分析Table 3 Variance analysis of regression models for soil angle of repose

通过土壤堆积试验测定茶园土壤的实际堆积角，试验装置主要由漏斗和坡度仪组成，将茶园的土壤样品置于漏斗后使其自由落下，

在漏斗下形成小型土堆，取4个方向，每个方向之间间隔90°，使用坡度仪分别测量土堆的坡度，试验过程重复4次，取平均值为茶园土壤的堆积角，结果如表4所示。

表4 土壤堆积试验结果Table 4 Test for soil angle of repose

将试验所得的土壤堆积角代入回归方程，得到39组优化解，将各组优化解分别进行堆积试验仿真，得到其堆积角，取与实际堆积试验最接近的一组为 EDEM仿真所用的土壤参数。由表4可知，土壤堆积角与实际物理试验的误差为 0.044%，且仿真土壤堆的形状与真实土壤堆非常相似，如图1所示，最终标定得到的茶园土壤参数中，土壤-土壤恢复系数为0.28，土壤-土壤滚动摩擦系数为 0.37，土壤-土壤静摩擦系数为0.55。

图1 土壤堆积角试验与堆积角仿真对比Fig. 1 Comparison of soil accumulation angle test and accumulation angle simulation

1.2 基于鼹鼠爪趾的仿生耕作机构设计

茶园的耕作环境较为复杂，与传统的耕作环境有较大区别。茶行的间距较小，且行间的土壤中有较多茶树树枝堆积，树枝较硬，传统的旋耕刀和深松趾都会因树枝缠绕而无法进行正常耕作。鼹鼠的挖掘过程由工作行程和收臂行程两个部分组成（图2-A），基于四杆机构的茶园耕作机构可以很好地模仿鼹鼠的挖掘过程（图2-B），因其中具有收臂冲程，从而可以有效避免缠绕现象。耕作机构的装配图如图3所示。

图2 鼹鼠前足爪趾与耕作铲铲尖运动轨迹对比图Fig. 2 Comparison of the movement trajectory of the claw tip of mole forefoot and the tip of the tillage shovel

图3 基于四杆机构的茶园耕作机构Fig. 3 Tea garden tillage mechanism based on four-bar mechanism

1.2.1 三维点云提取和模型的建立

使用逆向工程技术对鼹鼠爪趾的几何特征进行提取，对鼹鼠爪趾进行处理后，使用3D扫描仪对其进行360°扫描，随后进行删除噪点、修复、平滑等操作得到鼹鼠爪趾的精确模型，如图4-A和图4-B所示。

图4 鼹鼠爪趾仿生特征获取Fig. 4 Mole's claw biomimetic feature acquisition

1.2.2 仿生曲线的提取

将重构后的模型导入 Solidworks2021进行特征选择，利用Matlab2019b进行坐标提取及特征拟合。鼹鼠前足五趾趾尖有较好的入土性能和滑切性能，其中第3趾的趾尖起着主要的入土作用，第2趾的趾尖起着主要的切土作用，其余爪趾的趾尖具有较好的滑切性能[12]。因此选择第3趾的特征进行拟合。

两个特征部位如图4所示，图4-C为趾尖特征曲线，图4-D为趾趾特征曲线，利用Matlab中的Curve fitting tool工具对其进行多项式拟合。

趾尖拟合方程：

趾趾拟合方程：

两条曲线特征的拟合曲线方程的R2均大于0.99，说明函数方程的拟合度很好。

1.2.3 仿生耕作趾的设计与建模

原型趾为浙江川崎茶业机械有限公司现用耕作趾（图5-A），耕作趾尺寸参数如图5-B、5-C所示，刃倾角α=30°，刃口长为16 mm，外圆弧半径R1=290 mm，内圆弧半径R2=310 mm，趾趾仿生特征用于替代圆弧R2。

图5 挖掘铲结构及参数示意图Fig. 5 Schematic diagram of the structure and parameters of the tillage shovel

结合所得到的拟合曲线，使用Solidworks进行曲面建模，将两个仿生特征分别与原型趾结合，并且除必要的仿生特征外，其余参数均与原型趾相同，并将两个仿生特征叠加设计，得到复合仿生趾。仿生趾结构示意图分别如图5-D—F所示。

1.3 离散元仿真

1.3.1 接触模型选择

采用Potyondy等[31]提出的Hertz-Mindlin with bonding模型。该模型用于模拟破碎、断裂等问题，采用小颗粒黏结成大块物料，外力作用下颗粒间黏结力会发生破坏，从而产生破碎及断裂效果。Bonding模型适用于模拟混凝土或岩石结构等的破碎效果，也可用于模拟土壤结构，通过分析断裂键的变化规律来研究农机具的碎土效果。锄齿与土壤之间的接触则采用默认的Hertz-Mindlin（no-slip）模型。

Hertz-Mindlin with bonding模型在EDEM中的应用参数，包括法向刚度Kn、剪切刚度Ks、单位面积法向刚度σ、单位面积切向刚度γ、黏结半径r2。各参数分别为Kn=1×108N·m-3[32]，Ks=5×107N·m-3[32]，由公式（4）计算得r2=0.005 9 m，由公式（5）、（6）计算得σ=24 340 Pa，γ=13 840 Pa。

1.3.2 挖掘机构与土壤相互作用离散元模型建立

为模拟鼹鼠挖掘的过程，耕作机构由四杆机构组成，使其运动过程与鼹鼠挖掘过程类似，而 EDEM2021中无法定义较为复杂的机构运动过程，采用Adams-EDEM耦合的方式，模型能够在 EDEM中完成四杆机构的运动过程，且在 Adams2020的后处理模块中能够直接导出扭矩、角速度、角加速度等，同时也能在 EDEM的后处理模块中得到土壤颗粒和耕作机构的运动和受力情况。

为减少仿真所需时间又不影响仿真的精确性，将耕作机构进行一定程度的简化，取单侧的挖掘机构进行仿真且忽略轴承、螺栓等部件。将简化后的耕作机构材料参数设置为与Adams中steel材料参数相同。

采用长、宽、高分别为 0.9、0.5、0.4 m的模拟土槽，建立土壤颗粒模型时，为了减少仿真时间与运算量，设置颗粒半径为5 mm的球形颗粒，呈正态分布，标准差 0.05，使用EDEM2021中的Volume Packing功能生成土槽，该功能能在大批量生成土壤颗粒的同时进行预压，通过调节imposed solid fraction参数使其具有特定的孔隙率，从而大大缩短生成所需土槽的时间，imposed solid fraction设置为0.58，随机生成颗粒总数为19 708个。离散元仿真的具体参数如表5所示。

表5 离散元法仿真的基本参数Table 5 Basic parameters of discrete element method simulation

待土壤颗粒稳定后生成黏结键，将简化后的耕作机构分别导入到Adams和EDEM中进行耦合仿真求解。

2 结果与分析

2.1 仿生特征对扭矩的影响

为验证仿生趾的减阻效果，设置了多个不同的耕作条件，耕深为 10 cm和 8 cm，驱动臂转速为 180°·s-1、270°·s-1、360°·s-1，挖掘趾型为原型趾、趾尖仿生特征趾、趾趾仿生特征趾、复合仿生特征趾，进行全因子试验。仿真完成后从 Adams的结果文件中导出耕作过程中的扭矩数据，经过处理后得到的平均扭矩和最大扭矩分别如表6和表7所示，趾趾仿生特征在减小平均扭矩和最大扭矩方面都有明显的效果，平均扭矩和最大扭矩分别减小34.06%～39.29%和27.20%～37.62%；趾尖仿生特征可减小平均扭矩1.72%～5.04%，而在减小最大扭矩方面并没有较为明显效果；复合仿生特征可减小平均扭矩 36.61%～42.06%和减小最大扭矩20.67%～37.02%。平均扭矩减小比例结果如图6所示。

表6 平均扭矩Table 6 Average torque

表7 最大扭矩Table 7 Maximum torque

图6 平均扭矩减小比例统计图Fig. 6 Average torque reduction ratio statistics

2.2 仿生特征对水平阻力和垂直阻力的影响

过大的水平阻力和垂直阻力会使耕作机构产生振动，从而影响机构的耕作效果和稳定性，对于小型丘陵山地耕作机械的影响更为明显。由表8和表9可知，各仿生趾均能够起到减小耕作时的水平阻力和垂直阻力的效果，趾尖仿生特征能够减小0.49%～10.37%的最大水平阻力和3.46%～10.86%的最大垂直阻力；趾趾仿生特征能够减小 23.55%～38.89%的最大水平阻力和14.89%～36.63%的最大垂直阻力；复合仿生特征能够减小 32.11%～40.69%的最大水平阻力和18.59%～41.14%的最大垂直阻力。

表8 最大水平阻力Table 8 Maximum horizontal force

表9 最大垂直阻力Table 9 Maximum vertical force

使用第四强度理论对轴进行强度校核，公式如下：

式中，σr4为第四强度理论下轴所受到的最大应力，T为轴所受最大扭矩，M为轴在受到最大扭矩时刻所受的弯矩，Wz为轴的抗弯截面系数，L为力臂，Fx和Fy分别为轴所受的水平方向作用力和竖直方向作用力。

扭矩、水平力、竖直力均取最大值对轴进行强度校核，由表7～表9可知，最大扭矩为147.56 N·m，最大水平阻力为444.97 N，最大竖直阻力为 167.02 N耕作机构所用的轴直径为 30 mm，计算得σr4=48.47 Mpa，耕作轴材料为 45号钢，屈服极限为 355 Mpa，取安全系数 3，则[σ]=118.33 Mpa，得σr4=48.47 Mpa＜[σ]，完全符合要求。

2.3 仿生特征对能耗的影响

耕作机械产生能耗的主要部位是与土壤作用的刀具及传动机构，耕作所需能量减小能够提高耕作机械的持续工作能力。

式中，W为一次耕作过程所做的功，J；M为驱动臂的扭矩，N·m；ω为驱动臂的转速；η为功耗减少比例，%。

耕作机构进行一次耕作动作时扭矩所做的功W可由公式（14）求得，功耗减少比例η由公式（15）可知与驱动臂转速无关，对扭矩-时间曲线求积分，积分结果如表10所示，代入公式求得功耗减少比例，如图7所示。各个耕作条件下，趾尖仿生特征趾的功耗减少比例为1.58%～4.84%，趾趾仿生特征趾的功耗减少比例为29.02%～34.73%，复合仿生特征趾的功耗减少比例为30.84%～38.15%，功耗减少比例与平均扭矩减小比例的结果基本一致。从功耗减少比例与平均扭矩减小比例来看，趾尖仿生特征对此机构的减阻效果并不显著。

表10 扭矩-时间乘积结果Table 10 Torque time product result

图7 功耗减少比例统计图Fig. 7 Statistical chart of the power reduction ratio

2.4 仿生特征对铲面接触力的影响

以转速180°·s-1、耕深8 cm为例，在EDEM2021的后处理界面中依据不同速度对颗粒进行着色处理，得到耕作机构在4个典型时刻的土壤速度云图，如图8所示。从0.65～0.95 s，耕作机构对土壤造成了严重干扰，耕作土壤和未耕作土壤之间产生了强大的内部剪切力，使土层与耕作趾的内侧趾面发生挤压，从而产生摩擦阻力。土壤颗粒与趾面的接触力情况如图9所示，趾面与土壤颗粒的接触力都分布于内侧趾面，因此改变内侧趾面的曲线能够达到减阻效果。将趾面所受法向接触力和切向接触力导出后绘制成图，如图10和图11所示，趾趾仿生特征能够有效减小土壤对趾面的法向与切向接触力，趾尖仿生特征在减小法向与切向接触力方面并不明显且不稳定，因此趾趾仿生特征在减阻方面有明显的效果，趾尖仿生特征在减阻方面的效果较小。

图8 耕作过程的土壤速度云图Fig. 8 Cloud diagram of soil particle velocity during the tillage process

图9 原型铲铲面与土壤颗粒的接触力图Fig. 9 Contact force diagram of prototype shovel blade surface and soil particles

图10 耕作过程中铲面所受法向接触力图Fig. 10 Normal contact force on the shovel surface during tillage

图11 耕作过程中铲面所受切向接触力图Fig. 11 Tangential contact force on the shovel surface during tillage

3 结论

本研究通过试验测得了茶园土壤的实际物理参数，利用 EDEM软件模拟土壤堆积角试验，以茶园土壤间的恢复系数、静摩擦系数、滚动摩擦系数为试验因素，进行二次回归正交旋转中心组合试验，得出恢复系数、静摩擦系数、动摩擦系数对土壤堆积角的影响都非常显著，利用响应曲面优化方法得到土壤间接触参数的最优解为恢复系数 0.28、静摩擦系数0.58、滚动摩擦系数0.37。

根据鼹鼠爪趾优秀的破土效果，并结合鼹鼠的挖掘过程设计基于四杆机构的耕作机构，取其趾尖和趾趾仿生特征，在原型耕作趾的基础上设计出趾尖仿生特征趾、趾趾仿生特征趾、复合仿生特征趾。使用Adams-EDEM耦合仿真，模拟耕作机构单次挖掘过程，仿真结果表明，各个仿生特征趾都具有减阻效果，其中趾尖仿生特征的平均扭矩减小比例为1.72%～5.04%，功耗减少比例为1.58%～4.84%，能够减小 0.49%～10.37%的最大水平阻力和3.46%～10.86%的最大垂直阻力，但在减小最大扭矩方面无明显效果；趾趾仿生特征的平均扭矩减小比例为34.06%～39.29%，功耗减少比例为29.02%～34.73%，同时能够减小最大扭矩比例为27.20%～37.62%，能够减小23.55%～38.89%的最大水平阻力和 14.89%～36.63%的最大垂直阻力。从仿真结果来看，趾趾仿生特征可以有效减小耕作过程中趾面所受的法向接触力和切向接触力，从而达到减阻的效果，趾尖仿生特征在这一方面的效果较小且不稳定，因此减阻效果较小。复合仿生特征的平均扭矩和功耗减小比例分别为36.61%～42.06%、30.84%～38.15%，能够减小 32.11%～40.69%的最大水平阻力和18.59%～41.14%的最大垂直阻力，整体减阻效果约为两个仿生特征之和，同时具有与趾趾仿生特征相似的减小最大扭矩的能力，趾尖仿生特征和趾趾仿生特征具有一定的可叠加性，两个特征结合后设计得到的复合仿生趾具有更好的减阻效果。