新疆农田粉土离散元仿真参数标定

李清超，郑 炫，刘进宝，杨怀君，王子龙，张鲁云

（1.石河子大学机械电气工程学院，新疆石河子 832000；2.新疆农垦科学院机械装备研究所，新疆石河子 832000）

0 引言

【研究意义】新疆地区部分土壤属于粉土，土壤颗粒主要以粉粒为主，比重占到80%以上[1]，由于粒径较小，再加上新疆地区覆膜滴灌技术使用使得秋季农田土壤含有一定水分，具有一定黏性，不利于犁地作业[2～3]。新疆地区粉土土壤农田不适合免耕作业，耕翻成为基础性措施之一。目前对于新疆地区粉土土壤的犁地作业未能有准确描述，难以形成犁地翻垡作业的理论性指导，其主要的关键问题是缺少准确的参数模型进行仿真模拟，获取更优的土壤翻耕离散元模型，对机械设计人员有实际意义。【前人研究进展】离散元法被广泛应用在农业机械触土部件与土壤相互作用的研究上[4～8]，耕地模型及土壤与部件模型的准确性是保证结果有效的前提。建立一个系统、准确的新疆粉土土壤模型及粉土土壤与触土部件模型，是利用离散元分析方法对耕整地机械进行优化设计的基础。有研究通过物理试验与土壤仿真试验结合的方法，实现了对土壤离散元接触参数的标定[9～10]。Liang 等[11]以实测的土壤堆积角和土-膜的静摩擦系数为响应值，通过B-B试验，对土-土之间以及土-膜之间接触模型参数做了标定。戴飞等[12]用Hertz-Mindlin 无滑动模型，通过土壤堆积角、土壤与钢滑动摩擦角的试验方法，对全膜双茎沟覆膜土壤进行了标定。谢方平等[13]通过单轴密闭压缩试验和无侧限抗压强度试验，标定了基于Edinburgh Elasto-Plastic Adhesion（EEPA）模型的土壤离散元参数。【本研究切入点】新疆粉土土壤的待耕地，在常年滴水灌溉和耕翻作用下，其土壤耕性与其他地区不同土壤类型的耕地主要有阻力大、易与触土部件黏附等区别，仿真模型也会存在差异，难以通用。需要技术上重新标定，以校验用于耕地机械设计制造的依据。【拟解决的关键问题】以新疆地区粉土土壤的待耕地为研究对象，离散元仿真试验与物理试验相结合的方法，采用标定土壤与机具材料之间及土壤与土壤之间的离散元接触模型参数，为优化农业机械触土部件提供参考。

1 材料与方法

1.1 材料

1.1.1 土壤本构参数

土壤样品来源于2021年秋季作物收获后，采自新疆石河子市小麦屯村农场（44°21′31″N，85°55′59″E），前茬作物为棉花，土壤质地为粉土。采用5点取样法取耕翻作业0 ～30 cm深度的土壤，利用环刀（100 cm2）、YH-M50002电子天平（浙江东阳英衡智能设备厂，精度0.01 g）、DHG-9070A干燥机（上海一恒科技有限公司）、ZJ型应变控制式直剪仪（南京土壤仪器厂）、BT-2001 型激光粒度分布仪（湿法），对土壤的密度、含水率、泊松比、土壤粒径进行测量。密度、含水率、粒径的测量按照土壤试验标准进行；土壤泊松比的测量方法参照文献[14]。土壤本构参数：密度1 480 kg/m3，含水率14.5%，泊松比0.37，粒径分布＜0.001 mm、0.001 ～0.05mm和＞0.05mm分别占1.82%，85.55%和12.63%。计算泊松比。

表1 土壤本构参数Table 1 Soil constitutive parameters

1.1.2 土壤与65 Mn静摩擦系数

土壤与65 Mn之间的静摩擦系数是两者之间的特有属性，不会因形状改变而改变。土壤与65 Mn 之间的静摩擦系数由静摩擦试验台测得。将土壤用方刀制成1 cm×1 cm×1 cm的土块，放于65 Mn 斜面之上，旋转推杆来逐渐提高斜面角度，通过斜面仪来记录土块下滑瞬间的角度，由公式（2）计算得到静摩擦系数，重复10次，取平均值作为最终结果，得到土壤与65 Mn 之间的静摩擦系数为0.56。图1

图1 实际物理试验Fig.1 Actual physical test

式中：为静摩擦系数，为土块下滑瞬间角度。

1.1.3 土壤滚动距离

用65Mn 板搭建的斜坡试验台，标定土壤与65 Mn 之间的接触模型参数。土球颗粒由3D 打印的球形模具制得，半径为5 mm。在预试验的基础上，取斜面的高度为30 mm，斜面长度为100 mm，水平面长度为400 mm。将土球放于斜面上端中间位置处，使其沿斜面自然滚下，至土球停止运动。用直尺测量土球沿水平方向的滚动距离，重复10 次，平均值为227.29 mm，作为最终的目标值。

1.1.4 土壤堆积角

土壤堆积角试验装置，主要有万能试验机（AI-7000LAU30，高铁检测仪器有限公司）、65 Mn圆筒、65 Mn板组成。圆筒[15]内径为90 mm，高为270 mm，下方与65 Mn 板（500 mm × 500 mm）接触，内部装填土壤试样1 kg。采用万能试验机，以0.05 m/s 的速度向上垂直升起圆筒，土壤向下滑落在65 Mn板上，待土壤静置稳定后，用相机在垂直的方向上拍照，得到土样堆积角的正视图。后续使用matlab图像识别功能，提取边界，通过线性回归来获得土样堆积角的大小。重复10 次取平均值，得到土壤的堆积角为37.64°。图1

1.1.5 犁耕模型

选取颗粒形状为EDEM 软件自带的球状颗粒，半径为5 mm。采用EDEM 的Geometries 功能，生成2 500 mm×3 600 mm×400 mm 的土槽，颗粒总数为3617032个，将用solidworks绘制好的犁体曲面模型保存成x_t 格式并导入到EDEM中。前犁体和后犁体设置耕宽为450 mm，耕深为300 mm，水平前进速度分别设置为8、9 和10 km/h，总仿真时间为2.3 s，其中前0.3 s 生成土壤，后2s犁体模型完成在土壤中的运动。图2

图2 犁体仿真耕作过程Fig.2 Plow simulation tillage process.

1.2 方法

针对常用的耕翻机具材料65 Mn，采用静摩擦试验测定土壤与65 Mn 材料的静摩擦系数，运用斜坡物理试验与斜坡EDEM仿真试验相结合的方法来标定土壤与65 Mn钢之间的离散元接触模型参数。运用堆积角物理试验与堆积角EDEM仿真试验相结合的方法建立粉土土壤待耕地的离散元仿真模型。应用Design-Expert 软件，根据Central Composite 试验原理建立粉土土壤的斜坡滚动距离和堆积角与相关参数的回归模型，并找到最优参数组合。在最优参数组合下，比较堆积角的仿真值与实际值，以及耕翻仿真试验与耕翻田间试验的对照，验证模型的准确性。

1.2.1 接触模型选取

秋季待耕期间，由于膜下滴灌的影响，粉土土壤待耕地含有一定水分，颗粒间存在明显的黏附现象，且在耕翻时，土壤较容易黏附在触土部件上。Hertz-Mindlin with JKR 接触模型克服了经典的Hertz接触模型不涉及黏结力的缺点，通过表面能来表现颗粒间的粘结力，适用于有明显黏结和团聚的湿黏物料[16～19]。该模型的JKR 法向弹性接触力计算公式为：

式中FJKR：JKR法向弹性接触力，N；

E*：等效弹性模量（Pa）；

R*：等效接触半径（m）；

α：两颗粒切向重叠量（m）；

γ：表面能（J/m3）；

δ：两颗粒法向重叠量（m）。

其中，等效弹性模量和等效接触半径由式（5）和（6）定义：

式中E1、E2：两颗粒的弹性模量（Pa）；

V1、V2：两颗粒的泊松比；

R1、R2：两颗粒的接触半径（m）。取表面能γ为0时，此时JKR模型与Hertz接触模型一致，即：

即使两颗粒没有直接接触，该模型也会在一定距离范围内，存在黏结力，颗粒间存在黏结力的最大间隙为：

式中，δC和αC分别为颗粒间存在粘结力的法向和切向最大间隙。

当颗粒间隙大于0，小于δC时，两者间黏结力达到最大值，计算公式为：

JKR接触模型中的表面能，能较接近地模拟土壤颗粒之间以及土壤与触土部件之间地黏结力，选用Hertz-Mindlin with JKR 接触模型来建立新疆粉土土壤类型待耕地的离散元仿真模型。

1.2.2 EDEM仿真试验

应用离散元软件EDEM 2018 版本建立仿真模型，所有仿真试验的土壤颗粒形状选为球状，半径设置为5 mm，生成方式为随机分布，三维结构模型的材料均为65 Mn。离散元仿真试验中的土壤和65 Mn 的本构参数，通过测量或查阅参考文献[20]获得。表2

表2 土壤与65Mn的本构参数Table 2 Constitutive parameters of soil and 65Mn

斜坡仿真试验模型与物理试验保持一致，将solidworks 绘制的斜坡模型保存成x_t 格式，导入到EDEM软件中，坐标原点位置设置在斜坡顶端，中间位置处，并在此位置设置颗粒工厂，颗粒生成方式为dynamic，颗粒数量为1，数据保存间隙为0.01 s，重力加速度为9.81 m/s2。颗粒球顺着斜坡滑下，到静止时，将水平坐标导出，再减去95.39 mm（斜坡水平长度），即为颗粒的水平滚动距离。表3

表3 土-65Mn仿真试验因素水平编码Table 3 Factor level coding table of soil-65Mn simulation test

参考文献[14～15，20～22]的接触参数选取范围，采用二次回归正交旋转中心组合试验的方法，确定土壤与65 Mn接触参数的范围。表4

表4 土-土仿真试验因素水平编码表因素水平Table 4 Factor level coding table of soil-soil simulation test

将solidworks 绘制的堆积角试验模型保存成x_t 格式，导入到EDEM 软件中，圆筒和65 Mn 板各尺寸与物理试验保持一致。在圆筒上方创建颗粒工厂，生成颗粒方式为dynamic，数量为1 kg，待落在圆筒中的颗粒稳定后，圆筒以0.05 m/s 的速度垂直向上运动，颗粒在重力作用下向下滑落，最终颗粒停止运动形成堆积角，测量。图3

图3 仿真试验Fig.3 simulation test.

参考文献[14～15，20～22]的接触参数选取范围，采用二次回归正交旋转中心组合试验的方法，确定土壤接触参数的范围。表5

表5 斜坡仿真试验设计及结果Table 5 Slope simulation test design and results

1.2.3 阻力对比

田间试验在新疆石河子市小麦屯村农场（44°21′31″N，85°55′59″E）进 行，土 壤 含 水 率为15.51%，土壤坚实度为2152.8 kPa。使用的仪器有：雷沃欧M904 拖拉机、翻转犁、遥测仪（Autobona，黑龙江省农业机械工程科学研究院制），直尺、卷尺、秒表等。量取地块中间100 m作为测试区间，在保证相应耕作速度的前提下，通过遥测仪测得耕作阻力，并对区间内的试验结果取平均值。

2 结果与分析

2.1 土壤与65Mn接触参数标定

研究表明，二次回归方程如下：

该回归模型P＜0.000 1，失拟项P=0.132 0 ＞0.05，决定系数R2=0.983 5，校正系数Radj2=0.968 1，均＞0.9，模型具有良好的准确度，该模型可较好地表达滚动距离与土-65 Mn 恢复系数、土-65 Mn 静摩擦系数、土-65 Mn 动摩擦系数、土-65 MnJKR表面能之间的关系。4个因素对滚动距离影响的主次顺序为土-65 Mn 动摩擦系数、土-65 MnJKR 表面能、土-65Mn 恢复系数、土-65 Mn 静摩擦系数。其中，C、D对滚动距离影响极显著，A、B对滚动距离影响不显著。一次交互作用均对滚动距离影响不显著。二次交互作用中，A2、B2、D2对滚动距离影响显著。C2对滚动距离影响极显著。选定静摩擦系数为0.56 的一组最优解。选择土球在65 Mn 板上滚动距离相近的优化解。表6，表7

表6 斜坡仿真试验方差Table 6 Variance analysis of slope simulation test

表7 土-65Mn接触模型参数Table 7 Parameters of soil-65mn contact model

2.2 土壤与土壤接触参数标定

研究表明，建立堆积角Y2与土-土恢复系数X1、土-土静摩擦系数X2、土-土动摩擦系数X3、土-土JKR表面能X4之间的二阶回归模型。表8

表8 土-土参数标定试验设计与结果Table 8 Design and results of soil-soil parameter calibration test

模型P值小于0.000 1，失拟项P值大于0.05，决定系数R2均大于0.9，具有良好的准确度，此模型可较好地表达堆积角与土-土恢复系数、土-土静摩擦系数、土-土动摩擦系数、土-土JKR 表面能之间的关系。其中，有4个回归项影响极显著（P＜0.01），分别为X4、X1X4、X3X4、X22，2个回归项对模型影响显著（P＜0.05），分别为X2、X42。表9

表9 土-土参数标定试验方差Table 9 Analysis of variance of soil-soil parameter calibration test

取堆积角37.64°为目标，选择与目标值相近的优化解。选择软件给出的多组最优解中的第一个组解作为土-土接触模型的参数，最终确定土-土接触模型参数。表10

表10 土-土接触模型参数Table 10 Parameters of soil-soil contact model

2.3 物理试验验证

研究表明，土壤参数标定后的滚动距离、堆积角的相对误差分别为6.05%、1.28%。标定后堆积角的形状与实际形状更为接近。表11，图4

图4 土壤堆积角仿真试验与物理试验对比Fig.4 Comparison between soil accumulation angle simulation test and physical test.

表11 试验结果对比Table 11 comparison of test results

2.4 犁体阻力试验验证

研究表明，仿真的滚动距离比实际滚动距离大了6.05%，堆积角小了1.28%，在耕速10 km/h情况下，犁体耕作阻力小了6.49%；土壤含水率在田间试验时，与取土阶段相比发生了轻微的变化，也是造成误差的原因之一。所建立的新疆粉土土壤待耕地离散元仿真模型，在误差范围内是可以替代土壤的实际参数的，适合秋季进行耕翻的干旱地区。表12，图5

表12 阻力对比Table 12 Resistance comparison

图5 田间试验Fig.5 field experiment

3 讨论

研究标定的是新疆地区易于粘附的粉土土壤离散元模型，虽然与东北地区黑土的离散元模型选用的接触模型是一致的[20]，但具体接触参数标定结果不同，新疆农田粉土的JKR 表面能要明显小于东北黑土；北方地区的葡萄藤防寒土所选用的基础模型为整合延迟弹性模型和线性粘附模型，与研究的接触模型有一定区别[21]；与新疆棉田耕后土壤离散元模型相比，研究所选的接触模型更适合有一定粘附特性的土壤[22]。研究与黑土区玉米秸秆-土壤混料离散元模型相比，只考虑土-土之间以及土-65Mn 之间的接触参数，并未考虑残茬的相关作用[15]。

4 结论

4.1 利用二次回归拟合方法获得土-65 Mn之间最优接触参数组合：恢复系数0.51、静摩擦系数0.56、动摩擦系数0.08、JKR表面能4.12。

4.2 土-土之间最优接触参数组合为：恢复系数0.57、静摩擦系数0.65、动摩擦系数0.23、JKR表面能4.49。

4.3 对比仿真试验和实际试验的滚动距离被堆积角，误差分别为：6.05%、1.25%。

4.4 犁体8、9、10 km/h 作业速度下，仿真试验和实际试验的阻力相对误差为：5.69%、5.95%、2.74%，误差在可接受范围内，验证了标定的新疆粉土土壤离散元接触参数的准确性。