葡萄藤防寒土与清土部件相互作用的离散元仿真参数标定

马 帅，徐丽明，袁全春，牛 丛，曾 鉴，陈 晨，王烁烁，袁训腾

葡萄藤防寒土与清土部件相互作用的离散元仿真参数标定

马 帅，徐丽明※，袁全春，牛 丛，曾 鉴，陈 晨，王烁烁，袁训腾

（中国农业大学工学院，北京 100083）

为系统地研究中国北方地区沙壤土质地的葡萄藤防寒土及其与清土机清土部件常用材料（Q235钢、橡胶）相互作用的离散元仿真参数，以构建准确的土壤离散元仿真模型，该文选用整合延迟弹性模型（hysteretic spring contact model，HSCM）和线性粘附模型（liner cohesion model，LCM）作为土壤颗粒间的接触模型；基于土壤堆积试验，以土壤颗粒间恢复系数、静摩擦系数、滚动摩擦系数和土壤粘附能量密度为因素，以土壤堆积角为指标，利用EDEM进行通用旋转中心组合模拟试验，采用Design-Expert软件对试验数据进行回归分析，以实测的土壤堆积角作为优化目标值，获得土壤颗粒间的最佳接触参数组合；利用土壤屈服试验获得HSCM模型参数；基于斜面滑动法原理，利用倾斜板试验台测得土壤与Q235钢和橡胶之间的静摩擦系数，并以此为基础，采用土壤滑落试验，以滑动摩擦角为响应值，对土壤颗粒与Q235钢和橡胶之间的恢复系数和滚动摩擦系数进行寻优，得到最优解参数组合。为验证标定优化的离散元模型参数的准确性，采用刮土板土槽试验和仿真试验进行对比分析，获得刮土板在土槽试验和仿真试验中的水平前进阻力分别为228.36 N和213.79 N，两者之间的相对误差为6.38%，表明仿真模型中土壤的物理力学特性与实际土壤基本一致，验证了葡萄藤防寒土离散元仿真参数标定结果和研究方法准确可靠。研究结果可为基于离散元法研制适用于北方地区沙壤土质地的葡萄藤防寒土清土机提供理论基础和技术支撑。

土壤；离散元法；清土机；葡萄藤；防寒土；参数标定

0 引 言

中国葡萄种植区域广泛，其中北方葡萄产区由于冬季气候寒冷干燥，需在入冬前将修剪后的葡萄藤压倒在地面进行埋土防寒作业，等到第二年春天气候变暖，再把防寒土清除，葡萄藤上架[1-3]。其中，清土作业工作量大、劳动强度高、机械化程度低[4-5]，为促进葡萄产业的全程机械化发展，未来加大防寒土清土机的相关研制是必然趋势，而防寒土清土机作业时，清土部件与土壤直接接触，对于不同类型的土壤，其作业效果差异显著，故全面系统地研究葡萄藤防寒土离散元模型仿真参数，有助于为基于离散元法设计及优化防寒土清土机提供理论依据。

土壤离散元模型仿真参数主要包括材料本征参数、颗粒之间及颗粒与触土部件之间的接触参数、接触模型参数[6-8]。由于一般应用离散元法（discrete element method，DEM）建立的土壤颗粒模型尺寸和形状与实际土壤颗粒不同，因此，在运用离散元法对防寒土清土机与土壤相互作用仿真之前，需先根据实际土壤的质地类型，选择合适的土壤接触模型，并对其接触参数和接触模型参数进行标定[9-10]。

目前，国内外学者利用离散元法对土壤与农机部件的相互作用方面做了大量研究，对土壤等散粒体的物料特性参数标定做了大量工作。Ucgul等[11-12]通过结合Hertze-Mindlin及Hysteretic Spring接触模型，分析了土壤间有粘结力及无粘结力时的情况，同时考虑了土壤的塑性形变，为粘结力不大的土壤选择接触模型提供了参考；石林榕等[13]将线性内聚力模型添加到延迟弹性模型中，建立西北旱区农田土壤模型，并对不同含水率的土壤进行了离散元仿真参数的标定，该模型不仅可以体现土壤的塑性变形，而且能体现土壤间的粘结力特性；张锐等[14]应用默认的Hertz-Mindlin 接触模型，利用堆积角标定土壤颗粒在标准球和非标准球状态下相互作用的参数，建立了仿真模型，该模型只针对于沙壤土使用；王宪良等[15]基于代理模型的方法，利用The Edinburgh Elasto-Plastic Cohesion Model（ECM）弹塑性接触模型表示土壤颗粒接触模型，并对颗粒间的静摩擦系数和滚动摩擦系数进行了标定，该模型是考虑了土壤间弹塑性变形的一种非线性模型；李俊伟等[16]利用EDEM中的JKR模型对不同含水率黏重黑土颗粒之间的接触模型仿真参数及其与不同触土部件之间的接触参数进行了标定，该模型适应于含水率较高的土壤，一般水田中应用广泛；向伟等[17]基于土壤堆积试验对南方黏壤土的离散元仿真参数进行标定，构建了南方黏壤土精准的离散元仿真模型；丁启朔等[18]利用Hertz-Mindlin with Bonding模型建立了黏性水稻土的深松耕作离散元模型，该模型中颗粒间具有粘结键，适用于颗粒间粘结力特别大的土壤，尤其适用于混凝土和矿石等硬质颗粒。

本文针对中国北方地区沙壤土质地的葡萄藤防寒土（以宁夏、甘肃地区沙壤土质地的葡萄园为主），选用整合延迟弹性模型（hysteretic spring contact model，HSCM）和线性粘附模型（liner cohesion model，LCM）作为土壤颗粒间的接触模型，采用仿真试验与物理试验相结合的方法，对葡萄藤防寒土离散元模型参数及其与清土机中清土部件常用材料（Q235钢、橡胶）的接触参数进行标定优化，并采用清土作业中常用的刮土板进行土槽试验与仿真试验的对比分析，对标定的优化仿真参数进行验证，以期获得较为准确的防寒土离散元模型参数。

1 土壤颗粒仿真模型参数标定

1.1 土壤颗粒间接触模型的选取

不同的接触模型使用不同的方程计算颗粒之间的接触力和阻尼力，故针对不同的土壤质地类型，须建立不同的接触模型。本文以中国北方地区（以宁夏、甘肃地区为主）沙壤土质地的葡萄园为主，入冬前，采用埋土机从行间取土，覆盖在葡萄藤上，形成一个横截面近似等腰梯形的土垄，由于经过整个冬季的沉积和雨雪，第二年春天清除时，防寒土具有一定的坚实度，平均含水率约为6%~8%，土壤颗粒之间既表现出明显的散粒体物料特性，又有一定的粘附力，同时当土壤颗粒受到外力压缩作用时会产生塑性形变。

在离散元软件EDEM内嵌的接触模型中，延迟弹性接触模型（hysteretic spring contact model，HSCM）可用于颗粒产生塑性形变的场合；线性粘附接触模型（liner cohesion model，LCM）可用于一般性粘结颗粒的快速计算，亦可用于具有一定含水率的颗[12]。为了真实反映沙壤土质地类型的葡萄藤防寒土宏观特性，本文选用整合延迟弹性模型和线性粘附模型作为土壤颗粒间的接触模型。由于假设摩擦力限制了HSCM方程中切向粒子的运动，因此将粘附力增加到法向接触力中[19]，粘附力F/a(N)的大小为

式中为颗粒间凝聚能力密度，J/m3；A为颗粒间接触面积，m2。A的计算公式如下：

式中r为颗粒间接触区域的半径，m，定义如下：

式中r为接触颗粒等效半径，m；E为颗粒模型等效杨氏模量，MPa；F为法向接触力，N。法向接触力F的计算公式如下：

式中1为加载刚度，N/m；2为卸载刚度，N/m；δ为颗粒间法向重叠量，m；0为颗粒间残余重叠量，m。此外，颗粒间法向阻尼力的计算公式为

接触颗粒间的总法向力为粘附力、法向接触力和阻尼力之和，即

式中F为总法向力，N；F为法向接触力，N；F为阻尼力，N；F为粘附力，N。

1.2 仿真模型参数标定

1.2.1 本征参数

土壤颗粒本征参数包括颗粒形状尺寸、密度、泊松比和剪切模量，这些参数是固定的，一般从测量结果和其他参考文献中获取。在离散元模拟中，当涉及大量土壤作业即生成的颗粒数量非常多时，建模时生成与实际土壤颗粒尺寸大小相同的颗粒模型是不切实际的，并且若采用任意形状颗粒进行建模，将会成倍增加仿真计算时间[20]，故为了建模的方便性和计算的快速性，本文采用离散元软件EDEM中默认的球形颗粒，颗粒基体半径设定为5 mm，球径倍率分布为0.95～1.05[13]；供试土壤取自宁夏贺兰山地区沙壤土质地的葡萄园内，其含水率在6%～8%之间，采用容积为2×10-4m3的环刀和量程为3 000 g、精度为0.01 g的电子天平测量供试土壤的密度，经10次重复测量，取其平均值1 389 kg/m3作为最终的仿真参数；泊松比和剪切模量依据文献[13,21]分别取0.35和2.89×1010Pa。

Q235钢与橡胶是防寒土清土机中常用的清土部件材料，根据文献[8]，Q235钢的密度、泊松比和剪切模量分别为7 850 kg/m3、0.28和8.2×1010Pa，橡胶的密度、泊松比和剪切模量分别为960 kg/m3、0.45和3 448 Pa。

1.2.2 接触参数标定

土壤颗粒间接触参数的标定基于堆积试验，标定之前，采用漏斗法测定供试土壤的堆积角，测试装置如图1a所示，主要由支架、漏斗（上开口直径为195 mm，下开口直径为80 mm，高度为100 mm）、土壤堆积钢板、钢尺、数字倾角仪（分辨率：0.05°，精度：±0.2°）等组成。土壤从漏斗中缓慢滑下，在钢板上堆积稳定后，采用钢尺和数字倾角仪配合测量此时的土壤堆积角，试验重复10次取平均值，最终测试结果为31.48°。

采用整合HSCM和LCM模型的土壤颗粒进行堆积仿真试验，如图1b所示，采用随机分布方式在漏斗上方开口位置生成土壤颗粒，待生成的所有颗粒从漏斗中落下并停止运动时，模拟结束，利用EDEM后处理中自带的Protractor配合Clipping工具分别从锥面的方向和方向测定土壤堆积角，取其平均值作为最终的仿真试验结果。

1.数字倾角仪2.钢板3.土壤4.钢尺5.漏斗6.支架

根据1.1节选取的土壤颗粒接触模型，在堆积角仿真试验中，以土壤颗粒间的恢复系数、静摩擦系数、滚动摩擦系数和LCM模型中的粘聚能量密度作为待标定参数[22]，以土壤堆积角1作为试验指标，进行4因素通用旋转中心组合模拟试验，根据试验设计原则，确定编码系数为2.0，结合仿真预试验结果和文献[23-24]，确定堆积角仿真试验因素的编码结果如表1所示，共实施31组试验（0水平试验重复7次），试验方案与结果如表2所示。

表1 土壤堆积角仿真试验因素编码表

应用Design-expert 8.0软件对仿真试验结果进行显著性分析，结果如表3所示。由分析结果可知，在4个试验因素中，对土壤堆积角影响显著，对土壤堆积角影响不显著，和对土壤堆积角影响极显著；交互项因素中只有项对土壤堆积角影响极显著，其余因素的交互影响均不显著；二次方项因素中，2对土壤堆积角影响显著，2对土壤堆积角影响极显著，其余影响不显著。回归模型的值小于0.001，而失拟项的值大于0.05，说明回归模型极其显著且失拟不显著，回归有效。将不显著项删除后得到土壤堆积角1的回归方程如式（7）。

表3 回归模型显著性分析

注：**表示极显著（<0.01）；*表示显著（0.01<<0.05），下同。

Note: ** Means highly significant (<0.01); * means significant (0.01<<0.05), the same below.

1=-20.75+137.62+182.41+3.23×10-3-5.93×10-3

-129.172-189.672（7）

采用Design-expert 8.0软件Optimization-Numerical模块，以实测的土壤堆积角31.48°为目标，对堆积角回归模型进行寻优，得到多组优化解，经仿真验证试验，选取与物理试验得到的土壤堆积角最接近的一组，即土壤间恢复系数为0.51，静摩擦系数为0.65，滚动摩擦系数为0.06，粘聚能量密度为10 495 J/m3，此优化解下对应的堆积角仿真试验值为31.74°，与物理试验结果的相对误差为0.83%，无显著差异，说明所得的土壤颗粒间接触参数最优值准确可靠。从图2可以看出，采用优化参数进行仿真试验得到的土壤堆积锥形与实际物理试验得到的土壤锥形具有很高的相似性。

图2 土壤堆积角仿真与物理试验对比

1.2.3 接触模型参数测定

对于本文选用的HSCM和LCM模型，其中LCM模型的粘聚能量密度等于土壤内聚强度，已在1.2.2节中进行了标定；HSCM模型的屈服强度通过土壤屈服试验测得，测试装置如图3所示，主要由贯入圆轴（轴长150 mm、直径15 mm，底端带有厚度5 mm、直径20 mm的圆盘）、100 mm×100 mm×100 mm的空心铁盒、REGER万能试验机、底部支撑装置和顶部夹持装置等部分组成。颗粒间的阻尼系数和刚度系数采用EDEM中的默认值，分别为0.05和0.95[13,25-26]。

1.贯入圆轴 2.空心铁盒 3.底部支撑 4.土壤 5.顶部夹持装置 6.REGER万能试验机

试验前，根据空心铁盒容积大小，用电子天平测量1.39 kg供试土壤，装入空心铁盒并将顶部土壤压平，然后放在万能试验机上，贯入圆轴由顶部夹持装置固定，试验时，设置万能试验机向下运行速度为200 mm/min，在夹持装置的带动下，贯入圆轴匀速插入土壤中，在与万能试验机连接的计算机上保存得到贯入压力随插入位移的变化值，将压力值除以圆轴底端圆盘的水平横截面面积，即可得到贯入应力随插入位移的变化曲线，如图4所示。起始阶段，随着贯入圆轴插入深度的增加，贯入应力不断增大，当贯入圆轴的底端插入土壤中达到37 mm左右时，贯入应力开始下降，土壤发生失效，此时土壤受到的应力值为0.38 MPa，37～50 mm之间应力值趋于平缓，随着贯入圆轴继续下降，由于其底端不断接近铁盒底部，贯入圆轴本身开始受力变形，从而使应力值迅速上升。因此，HSCM模型中的屈服强度为0.38 MPa。

图4 贯入应力与位移之间的关系图

2 土壤与清土部件材料接触参数的标定

针对清土机清土部件的常用材料Q235钢和橡胶，需测量和标定土壤颗粒与其相互的接触参数。土壤与机械部件材料之间的接触模型为Hertz-Mindlin (no slip)，采用静摩擦试验测定土壤颗粒与2种材料之间的静摩擦系数，并基于静摩擦系数，采用土壤滑落试验标定土壤颗粒与2种材料之间的恢复系数和滚动摩擦系数。

2.1 土壤静摩擦试验

2.1.1 静摩擦试验原理

斜面滑动法是测量静摩擦系数的常用方法[16]，原理如图5所示，质量为的物体放置在倾角为的斜面上，其重力分解为2个力：沿斜面向下的分力×sin和垂直于斜面的分力×cos。当斜面倾角小于滑动临界角时，×sin小于物体与斜面间的静摩擦力，物体保持静止，随着斜面倾角的增加，×sin越来越大，当大于物体滑动临界角时，×sin＞，物体将开始沿着斜面下滑。基于图5中物体的受力平衡分析，通过公式（8）～（11），可计算得到两者之间的静摩擦系数。

式中为静摩擦系数。

注：为物体质量，kg；为斜面倾角，(°)；F为斜面对物体的支持力，N；为斜面对物体的摩擦力，N；为重力加速度，m·s2。

Note:is the mass of the object, kg;is the slope angle, (°);Fis the support force of slope on object, N;is the friction force of slope on object, N;is the acceleration of gravity, m·s2.

图5斜面滑动法测量原理

Fig.5 Measuring principle of slope sliding method

2.1.2 土壤静摩擦试验

利用自制的倾斜试验台进行静摩擦试验，如图6所示，主要由底板、倾斜板、土壤挡板、铁丝和升降台等组成。其中4个土壤挡板安装在倾斜板上，形成一个无盖土槽容器，在土槽容器中装满供试土壤并将其顶部土壤压平，以消除由于表面起伏对被测清土部件材料滑动的影响，倾斜板通过铁丝连接在升降装置上，通过升降装置的上下运动带动倾斜板绕其末端转动轴旋转。

1.Q235钢板 2.土壤挡板 3.底板 4.倾斜板 5.数字倾角仪 6.土壤 7.铁丝 8.升降台9.橡胶片

由于土壤与不同材料之间的静摩擦系数是土壤与材料之间的特有属性，不会因为材料形状和大小的改变而改变，故试验时将预先制备好的被测清土部件材料样品（50 mm×50 mm×3 mm的Q235钢板、50 mm×50 mm× 3 mm的橡胶片）放在土壤表面顶端，将数字倾角仪校正后放置在倾斜板上，设定升降装置向上运行的速度为200 mm/min，在铁丝的牵引下倾斜板开始缓慢倾斜，当被测触土材料样品在土壤表面开始滑动的瞬间，关停升降装置，记下此时倾角仪的度数，每个试验重复10次取平均值，测量结果如表4所示，依据公式（11）计算静摩擦系数，即tan20.81°=0.38，tan25.44°=0.48，土壤与Q235钢和橡胶之间的静摩擦系数分别为0.38和0.48。

表4 Q235钢板和橡胶片在土壤表面滑动瞬间的倾角

2.2 土壤滑落试验

2.2.1 试验方法

针对沙壤土颗粒，土壤与清土部件材料（Q235钢和橡胶）之间的恢复系数和滚动摩擦系数不易直接测定，基于已测定的静摩擦系数，利用土壤滑落试验进行标定，为更精确地观察、控制试验状态并测定相应的试验结果，利用土壤颗粒在倾斜板上完全滑落时，对应滑动摩擦角的大小作为标定试验依据值。测定装置采用2.1.2节的自制倾斜试验台，当测量土壤颗粒在不同的触土材料上滑动时，只需在倾斜板上土槽容器底部更换放置对应的触土材料即可。经前期预试验可知，选用土壤质量的多少与其在清土部件材料上的滑动行为及完全滑落时对应滑动摩擦角大小无关，本文选择100 g供试土壤进行试验。试验时，将土壤堆放在触土材料表面的顶部，如图7a所示，然后设定升降装置向上运行的速度为200 mm/min，倾斜板缓慢倾斜，当土壤从材料表面的顶部完全滑落时，如图7b所示，关停升降装置，利用倾角仪测定此时倾斜板倾斜的角度即为滑动摩擦角的大小，每个试验重复10次取平均值。

图7 土壤在橡胶上的滑落试验

对应的土壤滑落仿真试验如图8所示，在EDEM软件中利用颗粒工厂生成同等质量的土壤颗粒并堆放在触土材料表面的顶部，如图8a所示，然后设定倾斜板绕其底部旋转轴开始转动，当土壤颗粒完全滑落时，如图8b所示，利用EDEM后处理中自带的Protractor工具测量倾斜板的倾角，此时倾斜板的倾角即为滑动摩擦角的仿真测定结果。

图8 土壤在橡胶上的滑落仿真试验

2.2.2 试验设计

以土壤颗粒与Q235钢板和橡胶之间的恢复系数和滚动摩擦系数为试验因素，以土壤在2种材料上完全滑落时对应的滑动摩擦角2和3为评价指标，分别进行二因素通用旋转中心组合仿真模拟试验。参考文献[24,27-28]，确定仿真试验因素水平的编码结果如表5所示，共实施13组响应面分析试验（零水平试验重复5次），试验方案与结果如表6所示。

表5 土壤滑落仿真试验因素编码表

表6 土壤滑落仿真试验方案与结果

2.2.3 试验结果分析与参数优化

应用Design-expert 8.0软件对仿真试验结果进行回归显著性分析，结果如表7所示。

由分析可知，在2个试验因素中，土壤与材料之间的滚动摩擦系数对土壤在钢板和橡胶上的滑动摩擦角的影响极显著，而土壤与材料之间的恢复系数对土壤在钢板上的滑动摩擦角影响不显著，对土壤在橡胶上的滑动摩擦角影响显著。土壤在钢板和橡胶上滑动摩擦角的回归模型的值均小于0.01，而失拟项的值均大于0.05，说明回归模型极其显著且失拟不显著，回归有效。将不显著项删除后得到土壤在钢板上的滑动摩擦角2、土壤在橡胶上的滑动摩擦角2的回归方程如式（12）～（13）。

表7 回归模型显著性分析

2=20.23+13.13（12）

3=5.08+11.37+119.99−214.402（13）

2.2.1节的试验测得土壤颗粒在钢板和橡胶上完全滑落时的滑动摩擦角分别为25.30°和28.10°。通过Design-expert 8.0软件Optimization- Numerical模块，分别以实测的2个滑动摩擦角为目标值进行寻优，得到多组优化解，经仿真验证试验，选取与实测试验得到的滑动摩擦角最接近的一组，即土壤与钢板间的恢复系数为0.60、滚动摩擦系数为0.37，结合2.1.2节测得的静摩擦系数0.38，此优化解下对应的土壤颗粒在钢板上完全滑落时的滑动摩擦角为25.40°，与实测结果的相对误差为0.40%；土壤与橡胶间的恢复系数为0.61、滚动摩擦系数为0.23，结合2.1.2节测得的静摩擦系数0.48，此优化解下对应的土壤在橡胶上完全滑落时的滑动摩擦角为28.80°，与实测结果的相对误差为2.49%；由此表明，标定优化后的土壤颗粒与钢板和橡胶之间的接触参数准确可靠。

3 刮土板刮土验证试验

3.1 试验目的

标定的土壤离散元模型参数与土壤颗粒间的粘结力、内聚力、塑性形变等密切相关，选用不同的土壤模型参数建模，土壤所表现出来的宏观性质是不同的，其中表现较明显的是触土部件在土壤中运动时受到的阻力大小。采用刮土板进行刮土作业时，土壤在刮土板的作用下运动，不仅可以反映土壤颗粒模型间的相互作用，而且可以反映触土部件与土壤颗粒模型间的相互作用，故本文为验证标定优化后参数的可靠性和准确性，选用葡萄藤防寒土清除作业中常用的刮土板进行刮土试验，以作业时刮土板的水平前进阻力为响应值，将土槽试验中的实测值与EDEM中的仿真值进行对比分析，利用相对误差值判断所构建的葡萄藤防寒土离散元模型的有效性。

3.2 刮土板刮土土槽试验

3.2.1 试验条件准备

试验在中国农业大学工学院土槽中进行，为保证试验时土槽中土壤条件与实际葡萄园内的葡萄藤防寒土条件一致，试验前，通过人工制备土壤，包括旋耕、洒水、压实、翻运和起垄等处理，并根据春天清土时土垄的外形尺寸，制成了横截面近似于等腰梯形的2条土垄，每条长约22 m，土垄上底宽约25 cm，下底宽约70 cm，高约25 cm。用SC-900 型土壤紧实度仪测定土垄紧实度,从顶部到底部约在69～135 kPa之间，土垄底部紧实度大于顶部紧实度，用TDR 150型土壤含水率测定土垄的含水率约在6%～8%之间。

3.2.2 土槽试验

本文设计的刮土板机构主要由机架、刮土板、角度调节板和装在刮土板底部的橡胶等部分组成，通过三点悬挂系统挂接在土槽台车上，如图9a所示，刮土厚度通过升降系统调节，以土垄上表面为基准，设定刮土厚度为12 cm，试验时，通过操作控制台上的计算机，设定土槽台车的前进速度为1 km/h，数据保存频率为20 Hz，土槽台车带动刮土板前进并进行刮土作业，如图9b所示，将由六分力测力架测得的刮土板在作业时的水平前进阻力实时存储在计算机上。共进行2次平行刮土试验，并将2次平行试验测得的刮土板水平前进阻力的均值记为最终测试结果。

1.土槽 2.刮土板 3.角度调节板 4.橡胶 5.机架 6.土垄 7.六分力测力架 8.台车 9.控制柜 10.计算机

3.3 刮土板刮土仿真试验

运用EDEM软件建立土槽模型，设置其基本尺寸（长´宽´高）为1 800 mm´700 mm´270 mm，采用标定优化后的土壤离散元模型参数在土槽中生成总量为220 000个土壤颗粒模型，形成1 800 mm长的土垄模型，然后在土垄模型顶部设置一个向下运动的压板压实土壤，使建立的土垄离散元模型的紧实度和横截面尺寸与土槽试验的参数一致，然后将在SolidWorks中建立的刮土板三维模型按照1:1的比例导入EDEM中[29]，保证刮土板模型与土垄模型之间的相对位置与土槽试验时一致，试验装置材料为Q235钢和橡胶，其参数见1.2.1和2.2.3节，最终建立的EDEM仿真模型如图10a所示。

依据土槽试验条件，设定刮土板的刮土厚度为120 mm，水平前进速度为1 km/h，沿轴正方向。在EDEM求解器模块对仿真时间步长、仿真时间、数据保存间隔及网格大小进行设置[30]。为保证仿真的连续性，设置其仿真时间步长为1.40×10-6s（即Rayleigh 时间步长的40%），仿真时间为7 s，网格单元尺寸设置为颗粒平均半径的4倍，刮土板在EDEM中的仿真作业状态如图10b所示，待仿真作业完成后，在后处理中输出刮土板沿轴方向受到的阻力大小。

1.土垄模型 2.机架 3.刮土板 4.橡胶

3.4 试验结果与分析

土槽试验中，随着台车的前进，刮土板开始接触土壤，受到的水平阻力由0开始增大，以土垄起点为0点，当台车前进7.5 m左右时，刮土板与土壤稳定接触，阻力变化趋于平缓，如图11a所示，随着刮土板的前进，刮土板受到的阻力波动较大，分析原因可知，是由于人工起垄使土垄在不同位置的形状尺寸具有一定的误差以及土壤条件的差异造成。采用水平阻力的平均值作为最终的试验结果，通过计算其结果为228.36 N。

分析仿真试验可知图11b，刮土板作业时受到的水平前进阻力的整体变化趋势与土槽试验中大体一致，随着刮土板进入仿真计算区域与土壤颗粒模型开始接触，水平阻力不断增大，在1.5 s左右时，刮土板开始与土壤稳定接触，阻力变化趋于平缓，5.3 s左右以后，由于被刮土板壅积在最前方的土壤开始离开仿真计算区域，水平阻力值不断减小，水平阻力在1.5～5.3 s之间平缓变化，取其平缓阶段的阻力平均值作为最终的仿真试验结果，即为213.79 N。

图11 刮土板水平前进阻力的变化曲线

由此可知，刮土板在土槽试验中受到的水平前进阻力与仿真试验结果的相对误差为6.38%，误差在可接受范围内，表明了仿真模型中土壤的物理力学特性与实际土壤一致，验证了葡萄藤防寒土离散元仿真参数标定结果和研究方法准确可靠。

4 结 论

针对中国北方地区沙壤土质地的葡萄藤防寒土，选用整合HSCM和LCM作为土壤颗粒间接触模型，采用实测试验与仿真试验相结合的方法，对土壤离散元模型参数及其与清土机清土部件常用材料（Q235钢、橡胶）的接触参数进行系统标定，确定了最优参数组合，并通过土槽试验和仿真试验的对比分析进行验证，得到如下主要结论：

1）在EDEM中，以土壤颗粒间恢复系数、静摩擦系数、滚动摩擦系数和土壤粘附能量密度为因素，以土壤堆积角为指标，设计了通用旋转中心组合模拟试验，对试验数据进行回归分析，表明只有土壤间静摩擦系数对土壤堆积角无影响，其他接触参数均有显著影响；以土壤堆积角实测试验结果作为优化目标值，利用响应面优化方法获得土壤间接触参数的最优组合：恢复系数0.51、静摩擦系数0.65、滚动摩擦系数0.06、粘聚能量密度10 495 J/m3；利用万能试验机，采用土壤屈服试验测得HSCM模型参数屈服强度为0.38 MPa.

2）基于斜面滑动法原理，利用自制的倾斜试验台直接测定土壤与Q235钢和橡胶之间的静摩擦系数分别为0.38和0.48，并以此为基础，采用土壤滑落试验，对土壤与Q235钢和橡胶间的恢复系数和滚动摩擦系数进行标定优化，得到土壤与Q235钢之间的恢复系数和滚动摩擦系数分别为0.60和0.37，土壤与橡胶之间的恢复系数和滚动摩擦系数分别为0.61和0.23.

3）为验证所标定优化的离散元模型参数的准确性，采用刮土板的土槽试验和仿真试验进行对比分析，获得刮土板在土槽试验和仿真试验中的水平前进阻力分别为228.36和213.79 N，两者之间的相对误差为6.38%，误差在可接受范围内，表明仿真土壤模型的物理力学特性与实际土壤基本一致，验证了葡萄藤防寒土离散元仿真参数标定结果和研究方法准确可靠。研究结果可为后期基于离散元法研制相关葡萄藤防寒土清土机提供理论基础和技术支撑。

[1]徐丽明，李超，王文斌，等. 葡萄生产机械化技术与装备[J]. 新疆农机化，2012(5)：24－26.

[2]袁全春，徐丽明，马帅，等. 我国北方葡萄冬季埋土和春季清土作业的机械化现状与思考[J]. 中外葡萄与葡萄酒，2017(6)：66－67.

[3]马帅，徐丽明，邢洁洁，等. 叶轮旋转式葡萄藤埋土单边清除机研制[J]. 农业工程学报，2018，34(23)：1－10.

Ma Shuai, Xu Liming, Xing Jiejie, et al. Development of unilateral cleaning machine for grapevineburied by soil with rotary impeller[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2018, 34(23): 1－10. (in Chinese with English abstract)

[4]周伟彬. 葡萄起藤机的设计[D]. 银川：宁夏大学，2017.

Zhou Weibin. Design of Grape Vine Digging Machine[D]. Yinchuan: Ningxia University, 2017. (in Chinese with English abstract)

[5]郭家稳. 浅议哈密地区葡萄埋起藤全程机械化试验推广[J]. 农民致富之友，2017(8)：215－216.

[6]Janda A, Ooi J Y. DEM modeling of cone penetration and unconfined compression in cohesive solids[J]. Powder Technology, 2016, 293: 60－68.

[7]袁全春，徐丽明，邢洁洁，等. 机施有机肥散体颗粒离散元模型参数标定[J]. 农业工程学报，2018，34(18)：21－27.

Yuan Quanchun, Xu Liming, Xing Jiejie, et al. Parameter calibration of discrete element model of organic fertilizer particles for mechanical fertilization[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2018, 34(18): 21－27. (in Chinese with English abstract)

[8]刘文政，何进，李洪文，等. 基于离散元的微型马铃薯仿真参数标定[J]. 农业机械学报，2018，49(5)：125－135.

Liu Wenzheng, He Jin, Li Hongwen, et al. Calibration of simulation parameters for potato minituber based on EDEM[J]. Transactions of the Chinese Society for Agricultural Machinery, 2018, 49(5): 125－135. (in Chinese with English abstract)

[9]李博，陈军，黄玉祥. 机械与土壤相互作用的离散元仿真研究进展[J]. 农机化研究，2015(1)：217－222.

Li Bo, Chen Jun, Huang Yuxiang. Research progress of discrete element simulation of soil-machine Interaction[J]. Journal of Agricultural Mechanization Research, 2015(1): 217－222. (in Chinese with English abstract)

[10]Milkevycha V, Munkholmb L J, Chen Y, et al. Modelling approach for soil displacement in tillage using discrete element method[J].Soil and Tillage Research, 2018, 183: 60－71.

[11]Uchul M, Fielke J M, Saunders C. 3D DEM tillage simulation: Validation of a hysteretic spring (plastic) contact model for a sweep tool operating in a cohesionless soil[J]. Soil and Tillage Research, 2014, 144: 220－227.

[12]Ucgul M, Fielke J M, Saunders C. Three-dimensional discrete element modelling (DEM) of tillage: Accounting for soil cohesion and adhesion[J]. Biosystems Engineering, 2015, 129(1): 298－306.

[13]石林榕，赵武云，孙伟. 基于离散元的西北旱区农田土壤颗粒接触模型和参数标定[J]. 农业工程学报，2017，33(21)：181－187.

Shi Linrong, Zhao Wuyun, Sun Wei. Parameter calibration of soil particles contact model of farmland soil in northwest arid region based on discrete element method[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2017, 33(21): 181－187. (in Chinese with English abstract)

[14]张锐，韩佃雷，吉巧丽，等. 离散元模拟中沙土参数标定方法研究[J]. 农业机械学报，2017，48(3)：49－56.

Zhang Rui, Han Dianlei, Ji Qiaoli, et al. Calibration methods of sandy soil parameters in simulation of discrete element method[J]. Transactions of the Chinese Society for Agricultural Machinery, 2017, 48(3): 49－56. (in Chinese with English abstract)

[15]王宪良，胡红，王庆杰，等. 基于离散元的土壤模型参数标定方法[J]. 农业机械学报，2017，48(12)：78－85.

Wang Xianliang, Hu Hong, Wang Qingjie, et al. Calibration method of soil contact characteristic parameters based on DEM theory[J]. Transactions of the Chinese Society for Agricultural Machinery, 2017, 48(12): 78－85. (in Chinese with English abstract)

[16]李俊伟，佟金，胡斌，等. 不同含水率黏重黑土与触土部件互作的离散元仿真参数标定[J]. 农业工程学报，2019，35(6)：130－140.

Li Junwei, Tong Jin, Hu Bin, et al. Calibration of parameters of interaction between clayey black soil with different moisture content and soil-engaging component in northeast China[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2019, 35(6): 130－140. (in Chinese with English abstract)

[17]向伟，吴明亮，吕江南，等. 基于堆积试验的黏壤土仿真物理参数标定[J]. 农业工程学报，2019，35(12)：116－123.

Xiang Wei, Wu Mingliang, Lü Jiangnan, et al. Calibration of simulation physical parameters of clay loam based on soil accumulation test[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2019, 35(12): 116－123. (in Chinese with English abstract)

[18]丁启朔，任骏，Belal E，等. 湿粘水稻土深松过程离散元分析[J]. 农业机械学报，2017，48(3)：38－48.

Ding Qishuo, Ren Jun, Belal E, et al. DEM analysis of subsoiling process in wet clayey paddy soil[J]. Transactions of the Chinese Society for Agricultural Machinery, 2017, 48(3): 38－48. (in Chinese with English abstract)

[19]Zeng Z, Chen Y, Zhang X. Modelling the interaction of a deep tillage tool with heterogeneous soil[J]. Computers and Electronics in Agriculture, 2017, 143: 130－138.

[20]Shmulevich I. State of the art modeling of soil-tillage interaction using discrete element method[J]. Soil and Tillage Research, 2010, 111(1): 41－53.

[21]王金武，唐汉，王金峰，等. 悬挂式水田单侧修筑埂机数值模拟分析与性能优化[J]. 农业机械学报，2017，48(8)：72－80.

Wang Jinwu, Tang Han, Wang Jinfeng, et al. Numerical analysis and performance optimization experiment on hanging unilateral ridger for paddy field[J]. Transactions of the Chinese Society for Agricultural Machinery, 2017, 48(8): 72－80. (in Chinese with English abstract)

[22]郑侃. 深松旋耕作业次序可调式联合作业机研究[D]. 北京：中国农业大学，2018.

Zheng Kan. Study on Subsoiling and Rotary Tillage Operation Adjustable Combine Machine[D]. Beijing: China Agricultural University, 2018. (in Chinese with English abstract)

[23]张锐，韩佃雷，吉巧丽，等. 离散元模拟中沙土参数标定方法研究[J]. 农业机械学报，2017，48(3)：49－56.

Zhang Rui, Han Dianlei, Ji Qiaoli, et al. Calibration methods of sandy soil parameters in simulation of discrete element method[J]. Transactions of the Chinese society for Agricultural Machinery, 2017, 48(3): 49－56. (in Chinese with English abstract)

[24]戴飞，宋学锋，赵武云，等. 全膜双垄沟覆膜土壤离散元接触参数仿真标定[J]. 农业机械学报，2019，50(2)：49－56.

Dai Fei, Song Xuefeng, Zhao Wuyun, et al. Simulative calibration on contact parameters of discrete elements for covering soil on whole plastic film mulching on double ridges[J]. Transactions of the Chinese society for Agricultural Machinery, 2019, 50(2): 49－56. (in Chinese with English abstract)

[25]DEM Solution Ltd. EDEM 2.5 theory reference guide[R/OL]. 2014-12-05[2015-06-07].http://www.docin.com/p-980174717.html

[26]胡国明. 颗粒系统的离散元素法分析仿真[M]. 武汉：武汉理工大学出版社，2010.

[27]赵淑红，王加一，陈君执，等. 保护性耕作拟合曲线型深松铲设计与试验[J]. 农业机械学报，2018，49(2)：82－92.

Zhao Shuhong, Wang Jiayi, Chen Junzhi, et al. Design and experiment of fitting curve subsoiler of conservation tillage[J]. Transactions of the Chinese Society for Agricultural Machinery, 2018, 49(2): 82－92. (in Chinese with English abstract)

[28]孙伟，刘小龙，石林榕，等. 刮板升运带式膜上覆土装置覆土特性[J]. 机械工程学报，2016，52(7)：38－45.

Sun Wei, Liu Xiaolong, Shi Linrong, et al. Covering soil on plastic-film characteristics of scraper lifting belt mechanism[J]. Journal of Mechanical Engineering, 2016, 52(7): 38－45. (in Chinese with English abstract)

[29]方会敏，姬长英，Farman A C，等. 基于离散元法的旋耕过程土壤运动行为分析[J]. 农业机械学报，2016，47(3)：22－28.

Fang Huimin, Ji Changying, Farman A C, et al. Analysis of soil dynamic behavior during rotary tillage based on distinct element method[J]. Transactions of the Chinese Society for Agricultural Machinery, 2016, 47(3): 22－28. (in Chinese with English abstract)

[30]Di Renzo A, Di Maio F. An improved integral non-linear model for the contact of particles in distinct element simulations[J]. Chemical Engineering Science, 2005, 60(5): 1303－1312.

Calibration of discrete element simulation parameters of grapevine antifreezing soil and its interaction with soil-cleaning components

Ma Shuai, Xu Liming※, Yuan Quanchun, Niu Cong, Zeng Jian, Chen Chen, Wang Shuoshuo, Yuan Xunteng

(100083,)

Grapevine in seasonally frozen regions needs to be warm-insulated by soil in winter with the antifreezing soil removed in spring most mechanically by a soil removal machine. The purpose of this paper is to simulate the interaction between the insulating soil (with sandy loam texture) and the soil-cleaning materials (Q235 steel, rubber) commonly used in the soil removal machine, based on the discrete element method. The simulation model was constructed based on properties of the soil by integrating the hysteretic spring contact model (HSCM) and the linear cohesion model (LCM) as the contact model between soil particles. We took soil-soil restitution coefficient, soil-soil frictional coefficient, soil-soil rolling coefficient and soil cohesion energy density as the determinants and the soil accumulation angle as an evaluation index. The 4-factor universal rotation center combination simulation test, based on the EDEM, was used to regress the relationship between the determinants and the index using the Design-Expert software. The results showed that the soil-soil frictional coefficient was the only parameter that did not have significant effects on the soil accumulation angle. The best contact-parameter combination between soil particles was obtained by taking the physically measured soil accumulation angle as the optimization objective, which gave 0.51 for the soil-soil restitution coefficient, 0.65 for the soil-soil frictional coefficient, 0.06 for the soil-soil rolling frictional coefficient, and 10 495 J/m3for the soil cohesion energy density. The associated soil accumulation angle was 31.74o, with a relative error of 0.83% compared with the physically measured results. The universal testing machine for soil yield test was used to obtain the HSCM model parameters based on the change in penetration stress with the displacement, and the resultant soil yield strength was 0.38 MPa. The static frictional coefficient between soil and Q235 steel as well as the rubber measured by the inclination test bench was 0.38 and 0.48 respectively. These data and the EDEM were used to conduct the simulation test of the soil slip, with the restitution coefficient and the rolling frictional coefficient between soil and the materials taken as the determinants and the sliding frictional angle as evaluating index. Regressing the test date with the two-factor universal rotation center combination test showed that the rolling frictional coefficient between soil and the materials had a significant effect on the sliding frictional angle between soil and the steel plate and rubber. In contrast, the restitution coefficient between the soil and the materials impacted significantly on the sliding frictional angle between the soil and the rubber but not on the soil and the steel plate. The restitution coefficient and the rolling frictional coefficient between the soil and the Q235 steel as well as the rubber were optimized by using the measured sliding frictional angle as the optimization objective. The resultant restitution coefficient and the rolling frictional coefficient were 0.60 and 0.37 respectively for the soil and the Q235 steel, and 0.61 and 0.23 respectively for the soil and the rubber. Soil bin test and simulation test of the scraper were conducted to verify the calibrated discrete element model parameters. The horizontal forward resistance of the scraper in the soil bin test and simulation test was 228.36 and 213.79 N respectively, with a relative error of 6.38%. The results presented in this paper have important implications for using discrete element method to analyze the removal of grapevine-insulating soil.

soils; discrete element method; soil-cleaning machine; grapevine; antifreezing soil; parameters calibration

马 帅，徐丽明，袁全春，牛 丛，曾 鉴，陈 晨，王烁烁，袁训腾. 葡萄藤防寒土与清土部件相互作用的离散元仿真参数标定[J]. 农业工程学报，2020，36(1)：40－49.doi：10.11975/j.issn.1002-6819.2020.01.005 http://www.tcsae.org

Ma Shuai, Xu Liming, Yuan Quanchun, Niu Cong, Zeng Jian, Chen Chen, Wang Shuoshuo, Yuan Xunteng. Calibration of discrete element simulation parameters of grapevine antifreezing soil and its interaction with soil-cleaning components[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2020, 36(1): 40－49. (in Chinese with English abstract) doi：10.11975/j.issn.1002-6819.2020.01.005 http://www.tcsae.org

2019-09-04

2019-12-19

现代农业产业技术体系建设专项资金资助（CARS-29）

马 帅，博士生，主要从事生物生产自动化研究。Email：13637070719@163.com

徐丽明，教授，博士生导师，主要从事生物生产自动化技术和装备研究。Email：xlmoffice@126.com

10.11975/j.issn.1002-6819.2020.01.005

S152.9; S22

A

1002-6819(2020)-01-0040-10