分层防寒土与接触式清土机具相互作用的离散元仿真参数标定

杨启志， 赫明胜， 施 雷， 朱梦岚， 李章彦， 何文兵

(江苏大学 农业工程学院， 江苏 镇江 212013)

中国西北地区地域辽阔，物种资源丰富，地貌类型多样，具有发展特色农业的潜力和优势[1].但是由于其特殊的地理环境，干旱少雨，人口稀少以及机械化程度较低等问题导致了农业发展停滞不前.例如宁夏地区的葡萄产业，由于冬季寒冷干燥，需要通过埋藤的方式保证葡萄藤的安全过冬，并且在春季到来之前清除其表面的覆土，其中春季的清土起藤作业工作量大，机械化程度低[2].因此，发展高效的专用清土起藤的机械装备将成为促进葡萄酒产业快速发展的关键因素.

现葡萄园内使用的清土机多为接触式清土，其清土部件势必会与土壤直接接触，土壤类型的不同将会直接导致清土质量的好坏，由于经过整个冬天的沉积以及春季土垄底部的返潮，导致了不同的土壤层有着不同的坚实度以及含水率[3]，这给清土机的设计及优化制造了难度，因此防寒土离散颗粒群与清土机械构件间的复杂动力学问题研究是关键.

采用EDEM仿真技术是目前土壤离散颗粒群与机械构件相互作用机理研究行之有效的方法，其核心是获取准确的土壤及与机具的接触参数，准确建立防寒土离散颗粒群与清土机械构件的EDEM模型.构建准确的土壤离散元仿真模型，主要包括土壤的本征参数、颗粒材料间以及与接触部件之间的接触参数[4-5].其中土壤本征参数可以通过一般的试验方法获得，而土壤颗粒间及与接触部件间的接触参数则很难通过试验获得.因此国内外很多学者对颗粒仿真参数标定方面做了大量的研究工作.文献[6]基于土壤的堆积角试验，通过EDEM进行通用旋转中心组合模拟试验，获取了土壤接触参数间的最佳接触参数；文献[7]采用Hertz-Mindlin with JKR接触模型，基于堆积角试验，确定了适用黏性土壤的离散元接触参数；文献[8]采用Hertz-Mindlin with JKR接触模型，基于堆积角试验对南方地区黏壤土仿真参数进行标定；文献[9]对西北旱区农田土壤进行离散元仿真参数标定；文献[10]通过Hertz-Mindlin (no slip)与Hertz-Mindlin with JKR接触模型相结合，计算了地面效应中泥沙颗粒的运动和分布；文献[11-12]通过结合Hertz-Mindlin(no slip)及Hysteretic spring接触模型，分别对有黏结力和无黏结力土壤的仿真参数进行标定，为黏性不大的土壤选择接触模型提供了参考依据.

综上所述，现在对于离散元土壤模型参数标定已经有了许多的研究，但是大部分都是基于堆积角试验，评价指标单一；且针对西北地区的葡萄藤防寒土研究也仅仅是局限于表面土壤，对于不同土壤层参数的研究更是没有出现，导致了EDEM仿真土垄模型建立困难，使得仿真的准确性不高.

因此文中根据西北地区土壤沙性特质，将土垄分为3层，测定不同层土壤颗粒的本征参数，选用Hertz-Mindlin (no slip)接触模型进行堆积角试验和直剪试验[13-14]，然后通过滑落试验和二因素通用旋转中心组合仿真试验获得不同层土壤颗粒与清土部件间的参数，在EDEM中建立3层土垄仿真模型.并通过实际田间清土试验与仿真试验进行对比分析，最终确定防寒土土壤的关键特性参数，为后续清土部件的设计和优化提供参考.

1 不同层防寒土本征参数的获取

土壤样品从宁夏回族自治区红寺堡葡萄酒种植基地获取，根据土垄厚度60 cm，将土垄分为3层，底层0～20 cm、中间层20～40 cm、上层40～60 cm.

1.1 土壤基本参数获取

土壤基本参数包括颗粒形状尺寸、密度以及含水率.粒径测量试验采用筛分法，确定土壤类型为沙壤土.土垄底层土壤含水率为11.65%，中间层为9.46%，上层为6.77%.通过环刀法测量不同土层的密度，土垄底层土壤平均密度为1.655 g·cm-3，中间层土壤平均密度为1.573 g·cm-3，上层土壤平均密度为1.425 g·cm-3.底层、中间层及上层土壤泊松比为0.31、0.33和0.34，剪切模量为2.97×1010、 2.93×1010和2.90×1010Pa.

1.2 土壤直剪试验

通过直剪法测量土壤内摩擦角和黏聚力[15-16]，主要仪器包括预压仪(图1a)和直剪仪(图1b).

图1 试验装置

试验时，以0.8 mm·min-1的速度在不同垂直应力下进行剪切，记下量力环读数的峰值，每组试验重复3次，计算平均值.剪应力计算公式为

τ=kR,

(1)

式中：τ为剪切强度，kPa；k为量力环系数；R为量力环读数，mm.

根据试验结果拟合出每个土壤层试验的抗剪强度与垂直应力的一次方程，如图2所示.

图2 抗剪强度与垂直应力关系图

由图2可得，抗剪强度随着垂直应力的增加而增大，而在同一垂直应力下，土垄底层土壤的抗剪强度最大，上层土壤抗剪强度最小.抗剪强度与垂直应力回归方程符合摩尔-库伦理论[17-18]公式：

τmax=c+σtanφ，

(2)

式中：τmax为剪切应力最大值，kPa；φ为内摩擦角，(°)；c为黏聚力，kPa.

底层土壤的内摩擦角为33.052°，黏聚力为32.271 kPa；中间土壤的内摩擦角为30.233°，黏聚力为27.387 kPa；上层土壤的内摩擦角为28.443°，黏聚力为13.241 kPa.

1.3 土壤堆积角测定

堆积角也被称为休止角[19]，堆积角测量试验如图3所示.

图3 堆积角测量试验

土壤从漏斗中自由落下，在水平面上堆积形成夹角，使用数显倾角仪测量该倾角.每层土壤样本重复试验3次取平均值，最终获得底层土壤堆积角34.51°，中间层土壤堆积角32.28°，上层土壤堆积角30.12°.

2 土壤颗粒间接触参数与接触模型

目前对于离散元土壤模型参数标定已经有了许多研究，但大部分都是基于堆积角试验，将土壤视为松散颗粒，在自由落体下形成堆积角，评价指标单一.

针对清土过程中的土壤运动情况，文中将结合土壤堆积角试验与直剪试验来获取土壤参数.

2.1 参数标定试验的接触模型选取

西北地区酿酒葡萄种植产地普遍为沙壤土，为了真实的反映土壤特性，本研究选择Hertz-Mindlin (no slip)接触模型[20-23].

当两个球状颗粒发生弹性接触颗粒间的法向力Fn可由下式求得

(3)

式中：E*为等效弹性模量；R*为等效粒子半径；α为法向重叠量；r1、r2分别为颗粒球心位置矢量；R1、R2分别为颗粒半径；E1、E2和ν1、ν2分别为弹性模量和泊松比.

(4)

(5)

颗粒间切向力Ft可由下式计算：

(6)

式中：δ为切向重叠量；St为切向刚度；G*为等效剪切模量，G1、G2是两颗粒的剪切模量；R*为等效粒子半径.

(7)

2.2 颗粒间接触参数取值范围

根据EDEM软件中的颗粒材料数据库模块(GEMM)，获得所需的接触参数与JKR表面能的取值范围，如表1所示.

表1 仿真参数取值范围

2.3 仿真参数标定方法

根据颗粒间接触参数及取值范围，利用Design-expert软件设计Plackett-Burman试验，最陡爬坡试验法和Box-Behnken试验，以土壤堆积角、内摩擦角以及黏聚力为评价指标，筛选出在每个指标下影响显著的参数.

2.3.1试验设计

堆积角仿真试验如图4所示，建立与实际试验尺寸相同的漏斗，上口直径24 cm，下口直径5 cm.在漏斗上方生成土壤颗粒，在EDEM软件后处理中测量堆积角.在直剪仿真试验中，建立底面30 cm2、高2 cm的土样，对土壤颗粒进行加压，上部以0.8 mm·min-1的速度移动，获得最大横向剪切力F，结合公式获得土壤黏聚力及内摩擦角.

图4 仿真试验方法

2.4 结果与分析

2.4.1Plackett-Burman试验结果

试验结果如表2所示，利用Design-expert软件进行方差分析，结果如表3-5所示.由表可知以堆积角为指标的Plackett-Burman模型极显著，X1、X2和X3对堆积角、内摩擦角及黏聚力的影响均显著；以堆积角为评价指标，显著性大小排序为X2、X1、X3；以内摩擦角和黏聚力为评价指标，显著性大小排序为X2、X3、X1；其他因素均不显著.因此选择X1、X2、X3作为后续研究因素.

表2 Plackett-Burman试验设计与结果

表3 以堆积角为指标的参数显著性分析

表4 以内摩擦角为指标的参数显著性分析

表5 以黏聚力为指标的参数显著性分析

2.4.2最陡爬坡试验结果

基于Plackett-Burman试验，选择X1、X2和X3为自变量，每个因素设计5个水平，最陡爬坡试验结果见表6.由表可知参数值在水平4附近，因此选择水平4作为中心点，水平3和水平5分别为低、高水平.

表6 最陡爬坡试验设计与结果

2.4.3Box-Behnken试验结果及回归方程拟合

根据最陡爬坡试验获得的3个水平，结果如表7所示，研究不同土壤-土壤恢复系数(A)、土壤-土壤静摩擦因数(B)和土壤-土壤动摩擦因数(C)对堆积角(Y1)、内摩擦角(Y2)及黏聚力(Y3)的综合评价指标(Y)[24-25]的影响.3个参数与综合评价指标Y值的二阶回归模型为

表7 Box-Behnken试验设计与结果

Y=-18.813 7+12.492 4A+27.076 3B+

73.008 3C-4.095 2AB+9.333 3AC+

5.714 3BC-9.255 6A2-
17.412 2B2-247.644 4C2.

(8)

方差分析如表8所示，由表可知，A、B、C、AB、A2、B2、C2对Y值有显著影响.在回归方程的一次项中，A、B、C都对Y值有非常显著的影响；在交互项中，只有AB为显著项，而AC、BC均为不显著项，表明土壤-土壤恢复系数与土壤-土壤静摩擦因数之间的交互影响非常明显；在二次项中，A2、B2和C2极显著，说明土壤-土壤恢复系数、土壤-土壤静摩擦因数和土壤-土壤动摩擦因数对Y值的影响是非线性的.

表8 回归模型方差分析

2.5 土壤参数优化与试验验证

根据上述试验测得的不同土壤层的堆积角、内摩擦角以及黏聚力，通过这些实际测量值的Y值作为目标值进行寻优处理，借助回归方程得到不同土壤层的最佳土壤参数，其中底层土壤颗粒、中间层土壤颗粒以及上层土壤颗粒的恢复系数分别为0.523、0.507和0.458，其静摩擦因数分别为0.676、0.518和0.661，其动摩擦因数分别为0.141、0.169和0.140.将所获得的最佳参数再次进行土壤的堆积角试验和直剪仿真试验，并与实际测量结果对比，如表9所示，相对误差最大在4%左右，表明结果真实可靠，可以接受.

表9 仿真试验与实际试验结果对比

3 土壤颗粒与清土机具接触参数

3.1 土壤静摩擦因数

测量静摩擦因数采用如图5所示的斜面滑动法.随着斜面倾角θ的增加，mgsinθ也随之增大，当mgsinθ大于f(最大静摩擦力)时，物体开始下滑.通过公式(9)-(12)，计算获得静摩擦系数.

图5 斜面滑动法原理

FN=mgcosθ，

(9)

f=mgsinθ，

(10)

f=μFN，

(11)

(12)

式中：μ为摩擦系数；m为物体质量，kg；θ为斜面倾角，(°)；FN为斜面对物体的支持力，N；f为最大静摩擦力，N；g为重力加速度，m·s-2.

将质量为60 g的 Q235钢板放置在装满土壤的容器表面的一侧，匀速提升摆放钢板的一侧，当钢板开始滑动的瞬间，停止上升，记录数显倾角仪的数值，进行5次试验取平均值，Q235在底层、中层及上层土壤表面滑动倾角为23.75°、21.31°和18.26°.根据公式(12)计算出静摩擦因数，底层土壤为0.44，中间层土壤为0.39，上层土壤为0.33.

3.2 土壤恢复系数和滚动摩擦因数

在实际试验中发现，防寒土土壤颗粒与刮板部件材料之间的弹性恢复系数和滚动摩擦因数很难直接测量，针对宁夏地区土壤颗粒的特殊性质，采用斜面滑动法对两个参数进行确定.通过二因素通用旋转中心组合仿真试验获得3层土壤模型通用回归方程，将防寒土土壤颗粒完全滑落时对应的Q235钢板倾角作为目标值进行寻优，最终获得恢复系数和滚动摩擦因数的值.

试验时，利用自制的摩擦角测量装置进行测量，将土壤样本堆放在钢板的一侧，如图6a所示，然后转动手柄缓慢匀速提升放置土壤的一侧，当全部土壤样本从钢板表面完全滑落时停止提升，如图6b所示.利用数显倾角仪测量此时钢板的倾斜角，重复5次试验取平均值，底层土壤的滑动摩擦角为25.36°，中间层土壤的滑动摩擦角为23.59°，上层土壤的滑动摩擦角为22.12°.

图6 滑落试验

以恢复系数D和滚动摩擦因数E为试验因素，以倾斜角Y4(滑动摩擦角)为评价指标，对防寒土滑落试验进行仿真，如图7所示.在EDEM后处理中测量土壤颗粒完全滑落时底板的倾角.

图7 滑落仿真试验

根据恢复系数和滚动摩擦因数的取值范围，设置13组试验，试验结果如表10所示，通过Design-expert软件进行显著性分析，结果如表11所示.由分析可知，恢复系数D对滑动摩擦角的影响显著，而滚动摩擦因数E极显著.土壤在Q235钢板上的滑动摩擦角Y4的回归方程为

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

表11 回归模型显著性分析

Y4=-6.54+57.13D+115.57E-14.14DE-

49.99D2- 172.97E2.

(13)

以实测的滑动摩擦角为目标值进行寻优，得到不同土壤层的恢复系数为和滚动摩擦系数的值，底层、中层和上层土壤与Q235钢材之间的恢复系数0.622、0.588和0.549，动摩擦因数0.366、0.442和0.472.

4 试验验证

4.1 试验目的及方法

为了测试获得的防寒土仿真参数是否符合真实清土作业过程中土壤的特性，本试验结合田间试验与离散元仿真试验，验证仿真标定参数的准确性.在实际的刮板清土过程中，刮板对土垄两侧分别进行清土作业，在清土过程中土垄会在剪切、挤压等作用下坍塌.不同的土壤特性会呈现出不同的坍塌效果，因此以清土作业后剩余土垄的高度作为评价标准，将实测值与仿真获得值进行对比，利用两者的误差值判断标定的防寒土参数的正确性.

4.2 清土机田间试验

田间试验选在宁夏红寺堡葡萄酒种植地，为了保证在实际试验中土垄与仿真试验一致，在葡萄园区选取一条标准土垄，横截面为梯形，其上宽a约为70 cm，底宽b约为130 cm，高约60 cm，如图8a所示.刮土板总长为1 400 mm，高度600 mm，曲面半径400 mm，作业时侧边刮土厚度为300 mm，刮板的倾角45°.清土机的前进速度为1 m·s-1，刮土板对土垄两侧依次进行清土作业，作业过程如图8b所示.测量刮土板清土后的土垄高度，随机选5个测量点，取平均值作为最终的结果.

图8 田间试验

4.3 清土机仿真试验

在EDEM软件中建立土垄模型，其尺寸与实际土垄一致，土壤参数采用获得的最优参数.不同深度的土壤采用不同的颗粒形状，底层土壤采用块状三球体形状，中间层土壤采用双球体形状，顶层土壤采用默认的球形.生成的土垄如图9所示.然后将在Soildworks软件按照1 ∶1比例绘制刮土板模型导入到EDEM 软件中，其中刮土板的材料为Q235钢材，其密度、泊松比和剪切模量[26]分别为7 850 kg·m-3、0.28和8.2×1010Pa.

图9 仿真试验

设置刮板前进速度为1 m·s-1，入土深度、入土倾角保持与实际作业过程一致.在EDEM Simulator模块中设置时间步长为0.000 001 s，仿真总时间为10 s，其中0～7 s为土垄模型生成时间，7～10 s为刮板作业时间，数据保存间隔时间为0.05 s.随机选取稳定清土阶段的5个测量点，测量每个位置在刮土后的土垄高度，取平均值作为最终结果.

4.4 试验结果及分析

在田间试验中，清土作业后土垄高度的测量如图10所示.离散元仿真的土垄高度为52.89 mm，实际测量的土垄高度为55.82 mm，刮土板在田间试验中与仿真试验的相对误差为5.78%，在可接受的范围内，表明了通过试验与仿真结合测量的土壤力学特性参数及与清土部件的接触参数是准确的.

图10 实地测量

5 结 论

1) 依据接触参数取值范围，通过堆积角仿真试验以及直剪仿真试验，以土壤堆积角、内摩擦角以及黏聚力为评价指标，结合Plackett-Burman方法确定显著因素、最陡爬坡试验缩小范围、Box-Behnken响应面分析法确定回归方程，最终确定土壤-土壤恢复系数、土壤-土壤静摩擦因数和土壤-土壤动摩擦因数影响显著，其余参数均不显著.通过回归方程得到不同土壤层的最佳土壤参数，底层、中间层和上层的土壤颗粒间恢复系数为0.523、0.507和0.458，静摩擦因数为0.676、0.518和0.661，动摩擦因数为0.141、0.169和0.140.

2) 通过斜面滑动法确定3层土壤与Q235钢的静摩擦系数，其中底层土壤为0.44，中间层土壤为0.39，上层土壤为0.33，然后基于滑落试验和二因素通用旋转中心组合仿真试验获得不同土层的土壤颗粒与清土部件间的参数，其中底层土壤与Q235钢材之间的恢复系数为0.622、动摩擦因数为0.366，中间层土壤恢复系数为0.588，动摩擦因数为0.442，上层土壤恢复系数为0.549，动摩擦因数为0.472.

3) 为验证所标定优化的离散元土壤模型及与清土部件接触参数的准确性，通过田间试验与仿真试验进行对比分析，结果表明：仿真土壤模型与实际土壤误差较小，数值差异为5.78%，误差在可接受范围内，表明所标定的仿真参数准确可靠，为后期葡萄藤清土起藤机的设计和优化提供了参考.