麻山药种植区黏质土壤离散元参数标定与试验

刘月广，赵建国，齐浩凯，马志凯，李 浩，郝建军

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

麻山药是最早被认识和使用的中药材之一，有益脾、肺、肾的功效，主产地为河南、河北、山东、江苏、湖北等地［1］。平均直径为35 mm 左右，生长深度达0.8 ～1.2 m，平均含水量在80 % 左右，极限弯曲强度为4.7×105Pa，极限抗剪强度为1.4×105Pa［2］。麻山药物理特性脆，收获过程中极易折断，收获难度大。目前国内麻山药收获方式主要为机械开沟结合人工挖掘的半人工收获方式为主。该方式收获效率低、损伤率高，尤其在黏质土壤中，麻山药易折断现象更为严重，极大影响了麻山药收获后的品质及种植户的收益［3-4］。

麻山药采用深沟高畦种植，种植之前需深翻土壤（深度达1.5 m），破坏其犁底层和黏土层、改良其心土层和底土层性状［5］。麻山药与土壤固结部分主要为心土层和底土层，其生长环境土壤含水量高、黏性大，致使麻山药根-土结合机理比传统的农耕作物更为复杂，为麻山药机械化收获带来了困难［6］。因此，研究其生长过程中土壤状态和根土固结机理对根茎类作物的收获具有深远意义［7］。

为了研究不同的土壤特性，许多学者通过EDEM 软件中不同的接触模型来建立土壤模型［8-9］。杭程光等［10］采用Hertz-Mindlin（No slip），表征土壤的弹性变形以及流动性，分析了深松铲对土壤的扰动情况。Mustafa Ucgul 等［11］采用延迟弹性模型Hysteretic Spring Contact Model (HSCM) 模型，表征土壤的弹性和黏聚性，对所建立土壤模型进行直剪试验，得到不同垂直压力下土壤剪力线性增加的结果，证明该模型符合库仑定律，并且用该模型很好地模拟了土壤中的耕作力。文献［12-14］采用Hertz-Mindlin with JKR 用来表征土壤的吸附性和黏聚性，建立土壤模型，并验证了其准确性。石林榕等［15］采用延迟弹性模型Hysteretic Spring Contact Model（HSCM）和线性内聚力模型Liner Cohesion Model（LCM）表征土壤的弹性和塑性，针对西北旱区农田建立了土壤模型，并进行了鸭嘴插入土壤实验与仿真实验，通过试验与仿真的入土阻力对比，验证了其准确性。丁启朔等［16］用黏结模型（Hertz-Mindlin with Bonding）表征土壤的黏结破坏过程，建立了土壤模型，并通过对比仿真和田间试验的耕作阻力、土壤扰动、土壤破碎，验证了土壤模型的准确性。

参数标定对土壤模型的建立起到至关重要的作用［17-18］。许多学者采用土壤堆积角对离散元仿真过程中的颗粒参数和接触参数进行标定［13-14,19］。但以堆积角为响应值对参数进行标定，只能用于具有流动性的土壤模型。丁启朔［16］、祝英豪［20］等结合自己的研究，分别用耕作阻力和旋耕功耗对黏结模型（Hertz-Mindlin with Bonding）中的关键参数进行标定，得到了很好的效果。本文考虑到麻山药根土分离过程中是土和根的相对剪切运动，所以采用土壤剪切试验，以土壤极限剪力为响应值，对颗粒参数和接触模型参数进行标定。

1 参数标定过程

由文献［17］可知，颗粒材料的本征参数、颗粒之间的接触参数、接触模型的黏结参数对最终构建的土壤模型均有影响。通过筛选试验各项参数进行显著性分析，再对显著性参数进行优化可以快速地实现离散元参数标定［13,21］。本文参数标定过程如图1 所示。

图1 参数标定过程流程图Fig.1 Process of parameter calibration

1.1 土壤直切试验

由于心土层、底土层土壤一般位于地表300 mm以下，麻山药生长平均深度一般为600 ～1 000 mm。所以，将300 ～1 000 mm 按照100 mm 为单位均分7 层，每层土壤用3 个30 cm2×2 cm 的环刀取原状土，进行3 次直剪试验，3 次直剪试验的平均值即为该层土壤的剪切力。取土地点位于河北省高阳县麻山药种植区域。试验仪器为ZJ 型应变控制式直剪仪，压力级别在50 kPa 下以0.8 mm/min 的剪切速度在3 ～5 min 内剪碎。共进行24 组直剪试验得到8 组试验结果。直剪试验过程如图2 所示。

图2 土壤直剪试验过程Fig.2 Soil shear experiment process

直剪试验结果如表1 所示。由表1 可知，麻山药种植区域土壤心土层、底土层土壤质量和极限剪力并未有很大变化，这是因为在麻山药种植之前，对种植区域进行了深翻（深度达1.5 m），改良了心土层和底土层性状。由表1 可知，心土层、底土层土壤极限剪力均值为132.59 N。

表1 土壤直剪试验结果Table 1 Result of soil Diret Shear experiment

1.2 仿真参数标定方法

直剪仿真试验过程参照Freeman 单轴压缩试验［22］仿真标定方法。如图3 所示。在EDEM 中创建与实际土壤剪切盒尺寸相同的土壤剪切盒，先用直径为1 mm 的颗粒对其进行填充，再水平校准，之后用MCU 添加50 kPa 的压力，直至土壤颗粒稳定。对下剪切盒添加0.8 mm/min 的水平方向的剪切速度，直至上剪切盒水平方向受力不再增加为止。此时上剪切盒受到的水平方向的力即为土壤模型所承受的极限剪力。

图3 直剪仿真试验过程Fig.3 Simulation experiment process of diret Shear

1.2.1 Plackett-Burman 试验 本文采用Plackett-Burman 试验设计，以极限剪力为响应值，对仿真参数进行显著性分析。确定仿真参数11 个，分别用X1～X11表示。据文献［3, 10, 12, 13, 16, 19, 20,21］确定仿真参数的水平范围，如表2 所示。

表2 因素水平表Table 2 Table of factor level

1.2.2 最陡爬坡试验 通过Plackett-Burman 试验筛选出对极限剪力影响显著的因素，同时得到响应因子的显著性排序，在Plackett-Burman 试验结果的基础上，为使参数范围更小，针对显著的参数进行最陡爬坡试验。

1.2.3 Box-Behnken 试验 将显著性因素经过最陡爬坡试验优化后的参数范围作为Box-Behnken试验的高（+1）、低水平（-1），对该模型进行Box-Behnken 试验和响应面分析。同时为了试验的准确性，添加3 组中心水平试验，共进行15 组仿真试验。

2 参数标定结果与分析

2.1 Plackett-Burman 试验结果分析

Plackett-Burman 试验设计及仿真试验结果如表3 所示。对该试验结果进行分析，结果如表4 所示。由表4 可知，X2：土壤剪切模量(P＜0.01)；X9：黏结刚度(P＜0.01)，对极限剪力影响极显著；X5：土壤-土壤滚动摩擦系数（P＜0.05），对极限剪力影响显著，其显著性顺序为：X9＞X2＞X5。其余参数采用Plackett-Burman 试验中参数水平的中间水平。为了方便标记，后文中黏结刚度、土壤剪切模量、土-土滚动摩擦系数分别以A、B、C代指。

表3 Plackeet-Burman 试验设计及结果Table 3 Design and resule of Plackeet-Burman experiment

表4 Plackett-Burman 试验结果分析Table 4 Ansys of Plackett-Burman experiment results

2.2 最陡爬坡试验结果分析

最陡爬坡试验结果如表5 所示。由表中极限剪力结果可知，标定目标参数位于序号2、3 之间。所以，取2 为低水平、3 为高水平进行Box-Behnken 试验和参数优化。

表5 最陡爬坡试验设计及结果Table 5 Design and results of the steepest ascent experiment

2.3 Box-Behnken 试验结果分析

以黏结刚度(A)、土壤剪切模量（B）、土-土滚动摩擦系数(C)为响应因子进行Box-Behnken 试验。如表6 所示。

表6 Box-Behnken 试验设计及结果Table 6 Design and results of Box-Behnken experiment

该模型方差分析结果如表7 所示。结果表明：模型P＜0.000 1，表示该模型极显著；由各项P值可知，模型一次项A、B、C对极限剪力影响极显著（P＜0.01），交互作用项AB对极限剪力影响极显著（P＜0.01），二次项A2对极限剪力影响极显著（P＜0.01），二次项B2对极限剪力影响显著（P＜0.05），其余各项对极限剪力影响不显著。该模型二次多项式仿真模型拟合结果为：

表7 Box-Behnken 试验二次多项式模型方差分析Table 7 ANOVA of quadratic polynomial model of Box-Behnken experiment

模型决定系数R2=0.999 5，校正决定性系数adj-R2=0.998，非常接近于1，表示该回归模型拟合程度高。变异系数C.V.=0.5%，非常接近于0，表明该回归模型稳定性好，有良好的可信度和精度。

3 参数优化与土壤模型验证

在Box-Behnken 试验结果基础上，以实际测量极限剪力为目标值，对仿真参数进行优化求解，求得最优结果为：黏结刚度9.344×105N、剪切模量为5.976×106Pa、土-土滚动摩擦系数0.168。以最优参数（黏结刚度9.344×105N、土壤剪切模量为5.976×106Pa、土-土滚动摩擦系数0.168、其它不显著参数参考文献［3, 10, 12, 13, 16, 19, 20,21］）建立仿真直剪试验模型，求得极限剪力为131.16 N,与真实直剪试验所得结果基本一致，误差为1.07 %。

3.1 田间试验

3.1.1 试验设备 为了测量土壤切削过程中心土层、底土层对切削部件的阻力大小，自行设计了1种切土铲（ZL 2019 2 1829103.3），如图4 所示。田间试验在河北省高阳县棒药种植区开展。试验设备主要包括东方红LX1804 型拖拉机，CYB-605ST型销轴测力传感器（单个销轴量程0 ～4.9 kN）、XSR40 系列彩色无纸记录仪。

图4 切土装置Fig.4 Soil cutting device

3.1.2 试验过程 先用链式开沟机开2 条长800 cm、相距约为450 mm 深1 m 的深沟，再通过人工修整以保证2 条沟的距离为（400±5 mm）；在“土墙”的一端挖出尺寸为长800 mm、宽400 mm、深1 m，用来放置切土铲。田间试验前期准备如图5 所示。

图5 田间试验前期准备Fig.5 Preliminary preparation for field trial

将切土装置通过销轴测力器与拖拉机后悬挂挂接，3 个销轴测力传感器的输出端分别与记录仪的3个通道相连接。将在拖拉机上挂接好的切土装置放置在“土墙”一端挖好的坑内，调整拖拉机后悬挂上拉杆，保证切土装置与地表垂直且切土装置底部切土铲上面距地表为650 mm。开启测力传感器，拖拉机以0.017 m/s 的速度使切土装置在土墙下方650 mm 位置匀速驶过，实验过程如图6 所示。

图6 田间试验过程Fig.6 Process of field trial

3.2 土壤切削仿真试验

在EDEM2018 中，创建与田间试验相同工况的仿真模型。为了节约仿真时间，参考相关学者的颗粒放尺方法［23-25］，在不与切土装置接触的土层对颗粒进行缩放，缩放之后颗粒直径为5 mm；切土铲接触的土层采用半径为1 mm 的球形颗粒，按照标定后最优参数堆积土槽模型。土槽堆积完成之后将在Creo5.0 中建立的切土装置三维模型按照1 ∶1 导入EDEM 中，调整切土铲位置，设置切土装置速度为0.017 m/s，设置仿真时间为240 s。仿真过程如图7 所示。

图7 仿真试验过程Fig.7 Process of the simulation experiment

3.3 结果与分析

田间试验阻力与仿真阻力如图8 所示。

图8 田间试验与仿真试验阻力对比Fig.8 Prison of resistance between field experiment and simulation experiment

从0 s 开始，田间试验阻力逐步上升，约75 s时趋于平稳；而仿真试验阻力快速上升，约30 s 时趋于平稳。这是由于挖掘放置切土装置的土坑过程中，对邻近该侧的“土墙”结构造成一定损伤，影响其土壤固结程度，进而影响其阻力大小。田间试验数据波动较大，而仿真试验数据曲线平缓、波动小。这是由于拖拉机在运行过程中速度会出现微小波动，且试验过程中拖拉机和切土装置存在振动。而仿真过程为理想状态，速度平稳且无振动存在。

通过比较田间试验与仿真阻力曲线可知，试验稳定阶段，田间试验阻力曲线以仿真阻力曲线为基准上下波动，与田间试验阻力曲线的3 次多项式拟合曲线趋势大体一致。通过计算之后可知，田间试验稳定阶段(75 ～240 s)和仿真试验稳定阶段(30 ～240 s)阻力均值分别为：8 413.70 N 和8 587.68 N。仿真试验稳定阶段与田间试验稳定阶段阻力相对误差为2.06 %。误差在可接受范围内。表明该仿真标定结果与实际土壤力学特性基本一致。

4 结论

（1）通过对山药生长区域土壤的直剪试验，得到其平均值为132.59 N。该结论可以为麻山药收获机械结构、材料、动力等部分研发工作提供数据支撑。

（2）针对影响土壤仿真模型极限剪力的11 组仿真参数，进行Plackeet-Burman 试验、Box-Behnken试验，并以直剪试验得到的平均极限剪力为目标值对其进行优化，优化后的结果为：黏结刚度9.344×105N、剪切模量为5.976×106Pa、土-土滚动摩擦系数0.168。最后通过田间试验与仿真试验阻力对比，验证了优化参数下土壤模型的准确性。该结论确定了麻山药种植田土壤离散元仿真参数，解决了麻山药种植田土壤离散元仿真中土壤基质间黏结特性等相关问题，可以在此基础上，结合麻山药特性建立麻山药根土复合体离散元仿真模型，从而促进麻山药收获离散元仿真的研究以及麻山药收获机械的研发。