基于EDEM/FEM耦合算法平垄过程仿真

乔永昌，弋景刚，袁永伟

（河北农业大学 机电工程学院，河北 保定 071001）

根据农业部580个蔬菜重点县信息监测数据，2017年我国蔬菜累计种植面积1 313.25万a［2］，而蔬菜总需求量在2020年将达到5.9亿t［1］。温室起垄栽培工艺具有增加土壤透气性、提高水肥利用率、促进蔬菜根系生长、提高地表温度等优点［3-4］，所以温室多以起垄种植为主。目前温室整机机械主要以微耕机和大田旋耕机为主［5-7］，但是起垄种植模式使得地表高度差相对较大，旋耕机或微耕机在耕作时会因地表高低不平造成耕深一致性差、旋耕刀受力不均以及耕作后地表平整度差等问题，因此提高耕整土地耕深一致性、增加耕作后地表平整度以及降低旋耕刀波动受力是推广平垄作业的主要问题。常见的平垄方式主要是用推土板将垄推平，该种方式的缺点是消耗功率大，碎土质量差，平整效果差，且一次作业不能达到土地的待播状态。在平垄作业过程中，利用理论计算的方法很难准确分析出土壤对平垄部件的作用力，从而对平垄部件的强度计算及部件选型带来困难。离散元法是通过建立土壤颗粒模型，模拟仿真观察土壤颗粒运动情况，得出土壤颗粒间及土壤颗粒与工作部件之间的相互作用，从而对工作部件进行受力分析［8-9］。

本研究旨在基于EDEM/FEM耦合算法模拟平垄过程土壤对平垄部件产生的作用力及土壤粒子运动情况，对平垄部件工作参数选择并对关键部件进行强度校核，为平垄装置的设计提供可靠的依据。

1 整体结构设计

平垄作业是为了解决起垄种植土壤结构一次性旋耕作业不能达到旋耕作业耕深与作业后地表平整度要求的实际问题，平垄作业的主要形式与旋耕作业相同，都是靠刀具切削土壤实现松土的目的，但平垄作业与旋耕作业的区别在于平垄作业既要松土又要输土，具体来说就是要将一部分垄上土壤输送到垄沟处，实现耕作地面初步平整。平垄装置主要结构组成如图1所示。

图1 平垄部件结构示意图Fig.1 Diagram of ridge component structure

弯刀在刀轴上的布置形式不仅影响平垄作业质量、工作稳定性和功率消耗，而且也对刀轴的受力、机具所受阻力矩等影响较大，因此在同一作业区域内布置4把弯刀按照左弯刀、右弯刀周向间隔90°排列，同侧5组弯刀在平垄刀轴上按照4头螺旋线（螺旋线头数与同组内弯刀数量一致）布置，传动装置左右两侧弯刀对称式布置，该种布置方式较单头螺旋线方式工作时扭转力矩小，刀轴受力更为均匀。作业时传动装置左右两侧对称位置处左弯刀或右弯刀同时入土，减少刀轴所受轴向压力，保证整机左右两侧受力均匀一致。与弯刀交叉间隔式布置的带状螺旋叶片可将碎土输送至垄两侧的垄沟内，实现平垄作业。由图1可以知，在四头螺旋式弯刀与螺旋叶片组合布置的平垄装置结构中，带状螺旋叶片只输碎土而不切土，且带状螺旋叶片与平垄刀轴之间具有一定间隙，可有效防止平垄过程的壅土现象。

刀轴的左右两半轴上布置弯刀数Z1=Z/2=20把，其中Z为平垄部件弯刀总数，相邻两入土弯刀之间的转角为：

根据上述计算，将刀轴上的40把弯刀按照4头螺旋线、相邻入土转角18°进行布置，如图2所示。

图2 弯刀轴上布置平面展开图Fig.2 Planar expansion diagram of bending cutter arrangement on axis

2 离散元仿真模拟

应用离散元EDEM 2.7软件对弯刀与带状螺旋叶片组合布置的平垄部件进行耕作过程仿真模拟，既可以直观地观察平垄作业后土壤颗粒的运动状况和作业质量，又可以得到平垄过程工作部件的受力状况，为后期工作部件的受力分析提供可靠的依据。

2.1 离散元土壤接触模型与参数确定

为了准确模拟出平垄部件的作业性能，需要建立适合的土壤接触模型并确定合理的仿真参数。接触模型的本质是固体颗粒在准静态条件下的关于接触力学的弹塑性分析结果，不同的接触模型对计算固体颗粒的力与力矩产生很大影响，因此，为了提高模拟结果的可靠性，针对不同作业对象，需要搭建正确的接触模型［10］。在离散元法颗粒运动过程中，颗粒i受到重力mig和接触力Fi两者的合力。其中，接触力Fi包括颗粒之间或平垄部件与颗粒之间的法向碰撞接触力Fn,ij、切向碰撞力Fτ,ij、法向阻尼、切向阻尼和法向结合力Fcon,ij［11］。

本实验组作业区在河北省保定市，该地区土壤属于轻壤土（颗粒表面的粘附力小，具有一定的压缩性），因此土壤间接触模型为Hysteretic Spring和Linear Cohesion，其中 Hysteretic Spring 接触模型能够很好地模拟土壤颗粒塑性变形，而Linear Cohesion接触模型可以很好地模拟法向粘聚力［12-13］。通过查阅相关文献［14-17］和实际测量得到主要仿真参数见表1。

表1 仿真参数表Table 1 Simulation parameter table

2.2 EDEM平垄仿真试验设计

平垄部件左右两侧对称式布置，其工作状况基本一致，为了减少仿真软件计算量与和缩短运行时间，现对平垄部件中一侧刀辊进行仿真模拟，根据上一节对垄结构尺寸的分析，利用EDEM软件分别建立垄层和垄底层的土壤模型，由于采用单刀辊仿真作业，故建立半条垄，垄层尺寸为（垄高×垄顶长×垄底长×垄长）150 mm×400 mm×500 mm×5 300 mm（土壤颗粒数为56 800），为方便观察土壤粒子的运动情况，将垄底层适当加宽，则垄底层的尺寸为（长×宽×高）5 300 mm×1 200 mm×60 mm（土壤颗粒数为84 500）。将Creo三维软件建立的平垄部件模型导入EDEM软件中，仿真土槽模型与平垄部件模型如图3所示。

图3 仿真试验结构图Fig.3 Structural diagram of simulation test

仿真试验以作业后耕作地表与非耕作地表高度差和输土距离为主要评价指标，以平垄作业的前进速度、平垄刀轴转速、螺旋叶片螺距为影响因素，设计1个3因素3水平的正交试验进行平垄作业仿真模拟。通过查阅农业机械设计手册及参照相关旋耕机作业参数，初步确定前进速度为0.42、0.69、0.97 m/s，平垄刀轴转速120、180、240 r/min，螺旋叶片螺距120、160、180 mm。试验过程每种因素组合重复3次，仿真时间均为16 s，时间步长1.4×10-5s。

2.3 平垄模拟结果分析

根据仿真试验设计过程，每种组合因素试验中每间隔0.5 m取土槽模型横断面进行测量；取耕作地表与非耕作地表高度差和螺旋叶片输土距离2个参数值，结果如表2所示。

表2 仿真结果Table 2 Simulation result

续表：

由表2中极差分析可知，3种因素对耕作地表与非耕作地表高度差影响程度由高到低为：螺旋叶片螺距、平垄刀轴转速、前进速度；3种因素对输土距离影响程度由高到低为：平垄刀轴转速、螺旋叶片螺距、前进速度。由2种评价指标中的均值，优选出平垄部件选用作业参数为：前进速度0.97 m/s，平垄刀轴转速180 r/min，螺旋叶片螺距160 mm，并以该参数组合在此进行试验，其仿真模拟过程如图4所示。

图4 平垄部件仿真过程图Fig.4 Simulation process diagram of ridge-flat parts

在图4中由左至右划分4个作业后区域，其中1号区域为未耕作区，由于平垄部件减速箱是中置结构，导致平垄作业时产生漏耕现象，因此该区域土壤颗粒仍保持原有状态；2号区域为平垄作业区域，该作业区域内对垄底层土壤颗粒扰动量较小，作业后垄层剩余土壤厚度很均匀；3号区域是由平垄刀具、螺旋叶片与挡土板共同作用将土壤颗粒输送而形成的，该区域内土壤颗粒堆积较多，但整体高度较1号区域内低且较为平坦；4号区域内垄层土壤颗粒较少，说明螺旋叶片输土过程对非作业区域内的影响较小，不会影响已经平整后地表的精度。对仿真模拟过的土壤模型进行实际测量，其结果为：耕作地表与非耕作地表高度差为16 mm，平均输土距离103 mm，满足平垄作业要求，因此将该组参数确定为平垄部件作业参数。

3 耦合仿真分析

3.1 Hypermesh前处理

将导入到EDEM软件中的平垄部件中带状螺旋叶片模型（.stp）导入到HyperMesh软件中，首先以拓扑状态并进行渲染的方式来检查模型表面是否存在缺失或重叠，根据检查结果对几何模型进行清理修复［18-19］，最终成为闭合的几何体；然后进行网格划分，网格的类型和大小都会对仿真分析结果产生很大影响，在该软件中以六面体网格的收敛性最好且计算精度高［20］，故对螺旋叶片划分大小为5 mm的六面体网格，并对划分完成的网格进行质量检查，保证划分网格时的网格变形控制在不影响仿真分析质量的范围内。最终对带状螺旋叶片模型划分的网格模型得到104 121个节点和50 302个单元，网格划分模型如图5所示。

图5 网格划分模型Fig.5 Mesh generation model

3.2 Workbench静力学分析

静力学分析之前，需要将前期仿真分析得到的数据模型导入到ANSYS Workbench中，在EDEM软件中将仿真得到的土壤对螺旋叶片产生的应力结果以.axdt格式导出，然后再导入至ANSYS Workbench中，在HyperMesh中将划分好网格模型的螺旋叶片以.inp文件格式导出后也再导入至ANSYS Workbench中，至此完成软件间的数据交换。

工作过程的载荷主要由弯刀切削土壤时土壤对刀具的切削阻力和带状螺旋叶片排土过程受到的土壤阻力两部分组成，本研究仅针对带状螺旋叶片进行分析。在静力学分析时，将土壤对带状螺旋叶片部件产生的阻力设置成压力，且该压力大小数值与EDEM仿真过程导出的应力数据一致，将固定约束施加在带状螺旋叶片与刀轴连接位置处。最后，在ANSYS Workbench后处理模块中进行求解分析，得到带状螺旋叶片的总变形云图、等效应力云图和等效应变云图如图6、图7和图8所示。

图6 总变形云图Fig.6 Total Deformation Nephogram

图7 等效应变云图Fig.7 Equivalent strain nephogram

图8 等效应力云图Fig.8 Equivalent stress nephogram

由带状螺旋叶片的静力学分析结果可以得出：

（1）平垄作业过程中，带状螺旋叶片最大总体变形产生在其螺旋外边缘处，且靠近平垄部件中间位置的带状螺旋叶片外边缘处变形量最大。带状螺旋叶片输土过程越靠近其外侧输土量越大，且带状螺旋叶片外边缘处距离其固定端最远，所以造成该处变形量最大，最大变形为1.46 mm，在带状螺旋叶片使用许用范围内。

（2）带状螺旋叶片最大应变和最大应力均产生在其与刀轴连接位置处，且在两者固定连接位置处变形与应力最严重，最大应变为6.16×10-4mm，最大应力为122.22 MPa，等效应力、应变值均在带状螺旋叶片许用值范围内，即带状螺旋叶片的刚度和强度满足平垄作业的使用要求。

4 结论

（1）本文以EDEM、Hypermesh以及ANSYS Workbench三软件联合仿真的方法对平垄部件进行平垄过程仿真模拟以及关键部件的后处理分析，根据河北保定地区土壤类型选取土壤颗粒接触类型并建立了离散元仿真结构模型。

（2）建立离散元仿真正交试验，根据仿真模拟作业后耕作地表与非耕作地表高度差和输土距离2个评价指标确定平垄部件最优工作参数为：前进速度0.97 m/s，刀轴转速180 r/min，带状螺旋叶片螺距160 mm，并在该参数下再次进行平垄模拟，其作业结果为耕作地表与非耕作地表高度差为16 mm，平均输土距离103 mm，满足平垄作业要求。

（3）根据有限元后处理分析得出，带状螺旋叶片最大变形为1.46 mm，最大应变为6.16×10-4mm，最大应力为122.22 MPa，其等效应力、应变值均在带状螺旋叶片许用值范围内，表明平垄部件设计合理。