基于EEPA接触模型的土壤耕作特性模拟及颗粒球型影响分析

万里鹏程 李永磊 苏 辰 赵 虎 董向前 宋建农 王继承

(中国农业大学 工学院，北京 100083)

随着农业现代化的发展，数字技术、仿真技术越来越多地应用于田间土壤环境的模拟中。离散元法(Discrete element method, DEM)是专门用来解决不连续介质问题的数值模拟方法。该方法适用于工作对象可转化为离散颗粒的研究[1-3]，能够模拟土壤多颗粒单元聚合体[4-6]及受力条件下表现出的动力学性质，广泛应用于土壤与机具相互作用的研究领域[7-12]。国内外已有研究对构建土壤离散元模型的颗粒球型与混配方式应用方法差异大，探究离散颗粒球型对土壤离散元模型仿真精度及计算效率的影响很有必要。

土壤耕作特性是土壤在耕作时反应出来的性质，包括耕作难易程度和耕后土壤表现出的状态[13-15]，离散元法模拟土壤耕作特性是目前广泛应用的研究方式。EEPA接触模型兼顾粘性、弹性与塑性，与实际作业区域的土壤加载力再卸载过程展现出宏观土壤力学特性一致，应用于粘性土壤模型研究，刘宏俊等[16]利用EEPA接触模型对稻麦周年地区黏重土壤进行相关参数标定；谢方平等[17]采用单轴密闭压缩试验和无侧限抗压强度试验标定了基于EEPA模型的土壤离散元参数。在模拟计算散体物料力学行为时，颗粒形状是重要的影响因素[18-21]，王宪良等[22]以粒径为10 mm的球形颗粒为基础，仿照土壤颗粒形状填充了核状、块体、条状、片状土壤颗粒单元，进行了0.95～1.05倍粒径随机分布的球形颗粒填充土槽，建立了土壤模型标定方法；Hang等[23]采用粒径12 mm的球形颗粒并填充了三球颗粒与四球颗粒为颗粒单元，分析了深松铲齿间距对土壤扰动的作用；Sun等[24]使用粒径8～12 mm随机分布单球颗粒填充土槽，分析了耕作深度大于40 cm的仿生深松铲的减阻、土壤扰动特性。国内外已有研究应用离散元法模拟土壤，选择不同颗粒球型建立虚拟模型[25-30]，却未阐明所选颗粒球型对土壤离散元模型仿真精度及计算效率的影响。

本研究拟以土槽实验室的土槽土壤为研究对象，采用静态堆积角试验、单轴密闭压缩试验测量土壤摩擦特性、粘性、弹性和塑性参数，建立挖掘犁铲-土壤离散元耦合仿真模型验证参数优化准确性，进一步分析比较颗粒球型、颗粒粒径配比对计算效率、仿真精度的影响，以期为土壤-机具互作离散元仿真的建模，EEPA接触模型的参数选择，颗粒球型的确定提供依据和参考。

1 EEPA接触模型

EEPA接触模型以7个模型参数表征法向粘结力加载、卸载、再加载、再卸载过程，其加载/卸载路径模型[30]见图1：最初加载时，颗粒间法向粘结力沿初始加载路径①，卸载后，切换到卸载/再加载路径②并在特定重叠处达到零，描述为塑性重叠区长度δp。若发生再加载，颗粒间法向粘结力加载方式将保持在路径②上，但若超过历史最大加载力，则切换到初始加载路径①。卸载超过塑料重叠区长度δp，将导致粘合力的变化，直到最大粘合力fmin达到f0。在这一点上进一步卸荷，将使法向重叠量和吸引力都减小，直至在重叠区长度δ=0处发生分离。如果再次发生加载，接触力将沿着一条与初始卸载/再加载路径平行的路径②，直至到达最初加载路径①重合点后，在初始加载路径①上进一步加载。

、 法向粘结力加载/卸载方向， 法向粘结力加载/卸载路径。 fn为法向粘结力，N；f0为恒定拖拉力， N；fmin为颗粒间最大粘结力，N；δ为重叠区长度；δp为颗粒塑性变形量；δmin为最大粘结力的颗粒重叠量；δmax为重叠区长度最大值；K1为原始加载刚度；K2为再加载刚度；Kadh为衰减刚度；x为黏附分支曲线幂指数；n为非线性曲线幂指数。 and are loading/unloading direction of normal bond force, is the loading/unloading path of normal adhesive force. fn is the normal bond force, N; f0 is constant drag force, N. fmin is the maximum bonding force between particles, N. δ is the length of overlap. δp is the plastic deformation of particles. δmin is the particle overlap of the maximum bonding force. δmax is the maximum length of overlap region. K1 is the original loading stiffness. K2 is reloading stiffness. Kadh is attenuation stiffness. x is the power exponent of adhesion branch curve. n is the power exponent of nonlinear curve.图1 EEPA-加载/卸载路径模型Fig.1 EEPA-Load/Unload path model

加载路径下的原始加载刚度为K1，表达式为：

(1)

式中：E*为杨氏模量；R*为颗粒间接触距离；n1为加载分支指数；塑性变形比λp由初始加载刚度、再加载刚度决定，表达式为：

(2)

式中：K2为再加载刚度。塑性重叠区变形量δp与λp、加载分支指数n1密切相关，表达式为：

(3)

式中：δ为重叠区长度。颗粒分离前抵抗拉脱的最大粘结力fmin的表达式为：

fmin=πΔγψr

(4)

式中：Δγ为粘附能；ψ为粘着常数；r为颗粒接触半径。

综上，颗粒间接触力主要由EEPA模型参数、颗粒间重叠区长度计算得出。耕作土壤的宏观力为定值，应用不同粒径的颗粒单元建立对应土壤模型时，模型颗粒间重叠区长度差异大，导致模型参数变化大，模型宏观表现差异显著；而等效粒径下的土壤颗粒单元将减弱颗粒间重叠区长度变化差异，该条件下的模型优化参数组合在表征土壤宏观力学性质具有通用性。

2 土壤力学特性参数测定

室内土槽试验能够有效模拟田间工况和较为准确便捷的获得试验结果，且具有较强的可重复性和可对比性，是研究土壤-机具-植物互相作用关系的重要方法。为便于获取土壤耕作特性参数，选择中国农业大学工学院土壤—机器—植物系统技术实验室土槽土壤为对象开展相关研究。

2.1 土壤物理参数

室内土槽土壤来源为北京地区棕壤土，质地比较黏重，为获得实际土壤的物理参数，取试验土样应用烘干法测得土壤实际含水率为16.09%，采用环刀法测得土壤实际密度为2 130 kg/m3，三轴剪切试验测得土壤剪切模量为0.96 MPa，取土壤泊松比为0.36。

2.2 土壤静态堆积角

静态堆积角试验[4]常用来表征土壤的流动性、粘性与摩擦特性。试验采用内径为200 mm，高度400 mm的有机玻璃材料圆筒，将盛满土壤的圆筒开口端向下置于底板上，沿竖直方向0.1 m/s提起圆筒得到土壤堆积角。

采用MATLAB 2017R软件对静态堆积角试验的原始图像进行二值化处理得到二值图像，借助Canny阶梯型边缘检测算子进行边缘检测，最小二乘法将边缘轮廓拟合成直线，添加水平线后通过霍夫变换函数提取二值图像中直线夹角(静态堆积角β)数值为33.52°。具体测量过程见图2。

图2 土壤静态堆积角(β)的测量过程Fig.2 Measurement process of soil static accumulation angle (β)

2.3 土壤单轴密闭压缩试验

土壤在实际作业过程中表现出松散性、可压实性，在受压条件下土壤颗粒将产生塑性变形并对施压部件产生反作用力(包含弹力)。单轴密闭压缩试验能够实现土壤颗粒群的单向压缩，测得土壤轴向压力与轴向应变的变化关系曲线，表征出土壤弹性与塑性。试验用圆筒内径为95 mm，注入土壤高度150 mm，使用美国FTC公司的TMS-Touch物性分析质构仪，设置圆盘压板加载速度为60 mm/min，向下压缩圆筒内土壤12 mm，测得加载压力与土壤的回弹高度。重复10次试验，计算得轴向压力均值为653.43 N(图3)。

图3 单轴密闭压缩试验结果Fig.3 Result of uniaxial closed compression test

2.4 土壤耕作阻力与耕作堆积角试验

土壤耕作阻力与耕作堆积角试验在中国农业大学土槽实验室进行。物理土槽(标记为TC0)中土槽土壤制备过程见图4：采用TCC电力四驱土槽试验台车的洒水功能对试验所需的土槽土壤进行均匀的洒水作业，重复3遍洒水作业，待土壤存水湿透，静置3天；检查土壤含水率，确保土壤含水湿透，保证土壤含水率在17%左右，并使用1GQN-125型旋耕机将土壤旋耕，将表面土壤全面松碎并打匀，保证旋耕深度稳定在20 cm；采用镇压辊保持额定前进速度对旋耕后的土壤进行压实，镇压3遍，使用坚实度仪测量0～200 mm平均土层坚实度为247 kPa，烘干法测得土壤的含水率为16.09%。

耕作阻力测定试验综合考虑了挖掘犁铲工作参数许用范围与土槽工作环境，设置挖掘铲耕作深度为150 mm，挖掘犁铲前进速度为0.2 m/s，单次试验前进距离为2 m。共进行6组试验，耕作阻力值取稳定区间数据点均值。耕作试验过程见图5。

耕作阻力表征土壤耕作特性的耕作难易程度，当挖掘铲前进速度为0.2 m/s，耕深150 mm时，耕作阻力数值为1275.05 N；耕作堆积角为动态堆积角，表征耕作中土壤状态，在上述工况下测得耕作堆积角数值为32.29°。

图4 土槽土壤制备过程Fig.4 Preparation process of soil

图5 耕作试验过程Fig.5 Tillage test process

3 土壤离散元模型构建

3.1 土壤离散元模型参数确定

根据项目团队前期研究结果并参考文献[16]-[24]，以粒径为13 mm的单球模型填充虚拟土槽(标记为TC1)，土槽土层外形尺寸为1 500 mm×1 000 mm×220 mm，颗粒数量约20万个。土壤本征参数：密度、剪切模量和泊松比。土壤颗粒接触参数：碰撞恢复系数е取值0.2～0.6，颗粒间静摩擦因数μs取值0.3～0.9，颗粒间滚动摩擦因数μr取值0.2～0.8。土壤颗粒EEPA接触模型参数：粘附能Δγ取值12～36 J/m2，塑性变形比λp取值0.2～0.6，黏结分支指数Χ取值2.5～5.5，切向刚度因子Ktm取值0.3～0.7，参考文献[17]确定恒定拖拉力f0为0，加载分支指数n1为1.5。以土壤静态堆积角、土壤轴向压力响应面分析确定土壤离散元模型仿真参数优化组合。

3.2 仿真参数显著性分析

以土壤静态堆积角和轴向压力为指标，开展Plackett-Burman试验，分析碰撞恢复系数е、静摩擦因数μs、滚动摩擦因数μr、粘附能Δγ、塑性变形比λp、 黏结分支指数Χ、切向刚度因子Ktm等7个主要仿真参数的显著性。Plackett-Burman试验结果见表1。

以α=0.05为显著性水平，静态堆积角方差分析见表2，土壤静态堆积角模拟显著性影响因素为塑性变形比λp、粘附能Δγ、滚动摩擦因数μr、切向刚度因子Ktm的P值为0.035 8，显著性较低；轴向压力方差分析见表3，轴向压力模拟显著性影响因素为碰撞恢复系数е、切向刚度因子Ktm和静摩擦因数μs。

3.3 土壤静态堆积角响应面分析

采用Box-Behnken试验设计方法，以塑性变形比λp、粘附能Δγ、滚动摩擦因数μr为试验因素开展圆筒提升土壤静态堆积角模拟试验。EEPA接触模型仿真参数中切向刚度因子Ktm表征模型弹性的参数，为简化响应面分析试验方案，确定切向刚度因子Ktm为区间中值0.5。碰撞恢复系数е、静摩擦因数μs、黏结分支指数Χ等参数分别取区间中值为0.4、0.6和4.0。圆筒提升模拟试验采用内径200 mm、高度400 mm圆筒模型，粒径13 mm单球颗粒自由填充形成土层。圆筒提升速度0.1 m/s，仿真时间设置为30 s。土壤静态堆积角模拟结果见表4。

表1 Plackett-Burman试验结果Table 1 Results Table of Plackett-Burman test

表2 静态堆积角试验方差分析Table 2 Variance analysis of static packing angle test

表3 轴向加载试验方差分析Table 3 Variance analysis of axial loading test

表4 土壤静态堆积角模拟结果Table 4 Simulation results of soil static accumulation angle

表5 静态堆积角响应面方差分析Table 5 Variance analysis of static packing angle response surface

3.4 土壤轴向压力响应面分析

采用Box-Behnken试验设计方法，以碰撞恢复系数е、切向刚度因子Ktm、静摩擦因数μs为试验因素开展土壤单轴密闭压缩仿真试验。EEPA接触模型仿真参数中塑性变形比λp、粘附能Δγ、滚动摩擦因数μr分别取优化值：0.31、26.30和0.52，黏结分支指数Χ等参数分别取值为4。

单轴密闭压缩模拟试验采用内径为95 mm、高度150 mm圆筒模型，粒径13 mm单球颗粒自由填充形成土层。下压盘实体下降速度为60 mm/min，设置仿真时间12 s。土壤轴向压缩仿真试验结果如表6所示。

轴向加载力响应面方差分析结果见表7：试验因素Ktm、μs达到极显著水平，因素е达到显著水平，其他因素均不显著。以实测轴向压力653.43 N为目标值，采用Design-Expert 12软件分析，优选显著影响弹塑性的因素参数组合为：碰撞恢复系数е为0.58、切向刚度因子Ktm为0.39、静摩擦因数μs为0.67。

表6 土壤轴向压缩仿真试验结果Table 6 Results of axial compression simulation test

表7 轴向压力响应面方差分析Table 7 Variance analysis of axial pressure response surface

3.5 土槽土壤-挖掘犁铲离散元仿真试验

采用优化后的仿真参数构建土槽土壤-挖掘犁铲离散元仿真模型，挖掘犁铲材料为65 Mn根据参考文献[4]确定本征参数：泊松比ν1为0.3，密度ρ1为7 865 kg/m3，剪切模量G为79 GPa；土壤与挖掘犁铲接触参数：碰撞恢复系数е1为0.5、静摩擦因数μs1为0.3、动摩擦因数μr1为0.1。设置挖掘犁铲作业速度0.2 m/s、挖掘深度150 mm，仿真时间5 s，采样间隔0.05 s。

挖掘犁铲模拟耕作阻力与实测耕作阻力对比见图6。在某一稳定作业区间内模拟耕作阻力均值为1 253.77 N，实测耕作阻力均值为1 275.05 N，误差约为1.70%，此离散元模型较好的模拟了挖掘犁铲的耕作阻力。

图6 实测耕作阻力-虚拟耕作阻力对比Fig.6 Comparison between measured and virtual tillage resistance

4 颗粒球型对仿真过程的影响分析

4.1 颗粒球型结构尺寸确定

在前述研究基础上，以13 mm为等效粒径确定双球、三球、四球颗粒结构尺寸，简图见图7。双球模型颗粒半径为4.5 mm，(中心点坐标长度单位为mm，下同)中心点坐标分别为(2，0，0)、(-2，0，0)；三球模型颗粒半径为4 mm，中心点坐标分别为(0，0，2.5)、(-2.165，0，-1.25)、(2.165，0，-1.25)；四球模型的颗粒半径为3.5 mm，中心点坐标分别为(2.121，0，2.121)、(-2.121，0，2.121)、(2.121，0，-2.121)、(-2.121，0，-2.121)。

4.2 不同球型颗粒填充条件下虚拟土槽仿真试验

4.2.1虚拟土槽建立与仿真平台配置

此外，为了使新的生产系统能够灵活、可靠地工作，并且能够顺利运转物流过程，Leadec既使用现有的技术，同时也研发了自己的数字化解决方案，从而优化自身以及客户的内部流程。例如，Leadec的专业人员正在为生产、物流和原材料的控制开发个性化的软件解决方案。因为无论是在产品中还是在工厂基础设施中的传感器和执行器相互联系起来，都可以大大改变维护保养的措施。通过个性化的软件解决方案，Leadec可以更好地预见、规划这些改善措施，最终为客户创造更好的经济效益。

采用单球、单球变粒径、双球、三球、四球、多球混合等填充方式建立6种虚拟土槽，标记为TCi，各虚拟土槽的填充方式表示为：i=1，单球；i=2，单球变粒径；i=3，双球；i=4，三球；i=5，四球；i=6，多球混合。单球变粒径填充时以粒径13 mm为基准，尺寸缩放比为0.5∶1.0∶1.5，填充数量比为3∶4∶3；多球混合填充时单球、双球、三球、四球四种颗粒球型的尺寸缩放比和填充数量比为1∶1∶1∶1。

仿真平台硬件配置：处理器为Intel(R) Xeon(R) Gold 6226R CPU @ 2.90 GHz(双CPU，64核)，显卡(GPU)为GeForce RTX 3080。在EDEM2020软件中设置固定时间步1.2×10-5s、仿真时间5 s、写出时间间隔0.01 s、Grid Size(网格大小)为3Rmin。为有效评估颗粒球型对计算效率的影响，仿真过程不使用动态计算域。

图7 颗粒球型结构尺寸简图Fig.7 Size diagram of granular spherical structure

4.2.2评价指标

为合理评估离散元仿真模型仿真精度与计算效率，采用耕作堆积角相对误差φi、耕作阻力相对误差σi评价模型仿真精度；采用计算时间增长率τi评价模型仿真计算效率。

耕作堆积角相对误差φi(i=1，2，…，6)由式(5)计算：

(5)

式中：αi为虚拟土槽TCi中耕作堆积角的模拟值，(°)；α0为土槽土壤耕作堆积角实测值，(°)。耕作阻力相对误差σi由式(6)计算：

(6)

式中：Fi为虚拟土槽TCi中挖掘犁铲耕作阻力模拟值，N；F0为挖掘犁铲耕作阻力实测值。计算时间增长率τi由式(7)计算：

(7)

式中：Ti为虚拟土槽TCi仿真计算时间，s；T1为单球填充虚拟土槽TC1仿真计算时间，s。

4.3 试验结果分析

4.3.1模型仿真精度

土壤耕作阻力与耕作堆积角试验结果表明，物理土槽(TC0)土壤的耕作堆积角为32.29°、耕作阻力1 275.05 N，各虚拟土槽类型TCi中耕作堆积角模拟值αi测定见图8。

图8 不同土槽类型虚拟耕作堆积角模拟值的测定Fig.8 Determination of simulation values of virtual tillage accumulation angle for different soil groove types

各虚拟土槽条件下耕作堆积角模拟值与实测值的相对误差值由大到小依次为：TC1>TC3>TC4>TC2>TC6>TC5(表8)。虚拟土槽TC2、TC6、TC5耕作堆积角模拟值相对误差小于10%，能够较好模拟土槽土壤的耕作堆积角。单球、双球、三球颗粒模型结构相对简单，颗粒间互锁作用较弱，虚拟土槽TC1、TC3及TC4的耕作堆积角模拟值相对误差分别为20.32%、16.44%、13.04%，耕作堆积角仿真精度较低。

表8 耕作堆积角模拟值与实测值误差分析Table 8 Error analysis between simulated and measured values of tillage accumulation angle

各虚拟土槽条件下耕作阻力模拟值与实测值的相对误差值由大到小依次为：TC3>TC4>TC2>TC5>TC6>TC1(表9)，各虚拟土槽挖掘犁铲耕作阻力模拟值的相对误差值均小于8%，能较好模拟土槽土壤的耕作阻力。

表9 耕作阻力模拟值与实测值误差分析Table 9 Error analysis between simulated and measured values of tillage resistance

4.3.2模型仿真计算效率

采用计算时间增长率τi评价模型仿真计算效率，不同虚拟土槽类型(TCi)仿真计算效率分析见表10。可知：同等尺寸虚拟土槽颗粒填充数量由大到小依次为：TC3>TC4>TC6>TC1>TC5>TC2。随着颗粒表面凸起数量增加，接触频次增多，计算时间增长率由大到小依次为：TC5>TC6>TC4>TC2>TC3>TC1。双球填充方式仿真时间增加率为140.01%、颗粒填充数量约为31万个，四球填充方式仿真时间增长率为328.68%、颗粒填充数量增

表10 不同虚拟土槽类型仿真计算效率分析Table 10 Simulation calculation efficiency analysis of different virtual soil bin types

长率为4.23%。仿真效率受土壤颗粒填充数量、接触关系及颗粒球型等因素影响，大计算规模时颗粒球型对仿真效率影响更显著。

综合考虑仿真精度及计算效率，土壤耕作阻力模拟宜选用单球颗粒填充方式、耕作堆积角模拟宜选用四球或多球混合填充方式建立虚拟土槽。

5 结 论

本研究以室内土槽土壤为研究对象，采用EEPA接触模型模拟土壤耕作特性，建立了6种虚拟土槽，研究颗粒球型对土壤离散元模型仿真精度及计算效率的影响，得到结论如下：

1)通过试验建立起耕层土壤离散元模型，确定并优化了基于EEPA接触模型的单球颗粒土槽土壤离散元模型仿真参数：塑性变形比λp为0.31、粘附能Δγ为26.30 J/m2、滚动摩擦因数μr为0.52、碰撞恢复系数е为0.58、切向刚度因子Ktm为0.39、静摩擦因数μs为0.67、黏结分支指数Χ为4。

2)依次进行了6种虚拟土槽土壤-挖掘犁铲仿真试验，耕作堆积角模拟值与实测值的相对误差值由大到小为：TC1>TC3>TC4>TC2>TC6>TC5，其中TC2、TC6、TC5耕作堆积角模拟值相对误差小于10%，适用于土壤耕作堆积角等耕作土壤表现状态的研究；耕作阻力模拟值与实测值的相对误差值由大到小为：TC3>TC4>TC2>TC5>TC6>TC1，6种虚拟土槽的耕作阻力模拟值与实测值相对误差均小于8%，均可用于土壤耕作阻力的研究。

3)填充同等尺寸虚拟土槽所需颗粒填充数量由大到小依次为：TC3>TC4>TC6>TC1>TC5>TC2；土槽土壤-挖掘犁铲仿真计算时间增长率由大到小依次为：TC5>TC6>TC4>TC2>TC3>TC1。颗粒填充方式显著影响仿真效率，其中单球土槽仿真计算时间短，具有较优的计算效率。