植烟沙壤土与触土部件相互作用的离散元仿真参数标定*

林恒矗，王琪，廖鹏，郑文鑫，何金成，张德晖

(1. 福建农林大学机电工程学院，福州市，350002； 2. 福建省农业信息感知技术重点实验室，福州市，350002； 3. 现代农业装备福建省高校工程研究中心，福州市，350002)

0 引言

近年来，基于离散元法的虚拟仿真技术及专业离散元软件EDEM在农业机械研究中得到广泛应用，在土壤工作部件、精量播种、收获机械等方面取得了一系列重要进展[1]。如王金武等[2]利用EDEM进行胡萝卜收获机松土铲的设计与优化，其设计的松土铲相较于凿式松土铲，前进阻力降低5.79%，胡萝卜拔取力降低了20.68%；鲍玉冬等[3]进行了机械采收蓝莓的数值模拟研究，通过能量变化值分析果实损伤程度，为国内机械采收蓝莓的推广提供了理论支持；宋少龙等[4]进行了分层施肥靴的参数优化，为精准对行分层施肥机的改造提供了参考；杜俊等[5]完成了对于旋埋刀辊功率分配特性的研究，得出不同工作条件下，旋埋刀辊中螺旋横刀功率占总功率比例在60%左右；赵淑红等[6]设计了玉米秸秆还田交互式分层深松铲并进行了离散元仿真，得出交互式分层深松铲破茬比率高，土壤扰动小，能耗低，满足深松作业的要求；石志鸣等[7]通过离散元仿真试验研究了翼型深松铲铲间距对土壤蓬松度、扰动系数的影响，仿真试验与田间试验平均误差为11.2%、13.5%。

本研究利用万能试验机进行土壤堆积角试验，而后设计显著性影响试验、最陡爬坡试验、响应曲面试验，依据试验设计利用离散元软件对土壤进行仿真堆积角试验，得出仿真堆积角数值，并以目标值进行求解，确定了土壤—土壤碰撞恢复系数、静摩擦系数及动摩擦系数等基本接触参数，以标定的参数进行仿真堆积角试验，对比实际堆积角试验值，验证标定结果可靠性，在此基础上，以标定出的数值构建植烟土壤土槽模型，并进行土壤—触土部件相互作用的离散元仿真试验，得出触土部件所受扭矩，后续进行土槽试验，验证土槽模型可靠性。

1 植烟土壤堆积角试验

离散元仿真模型参数主要包括材料本征参数和基本接触参数，本征参数包括物料密度、颗粒体几何尺寸、物料剪切模量和泊松比等，基本接触参数包括颗粒与颗粒、颗粒与材料之间的静摩擦因数、滚动摩擦因数、碰撞恢复系数等。大部分参数可通过试验测定得出，而部分参数(尤其是物理接触参数)很难通过试验获取，因而许多学者提出采用实测试验结合虚拟标定实现仿真参数的精准标定[8]。而堆积角是表征颗粒物料流动、摩擦等特性的宏观参数，且标定的参数是否具有可信度，可以用仿真堆积角与实际堆积角对照方法验证[1, 9]，故测量实际休止角值作为参数标定试验的参照值[8, 10]。

1.1 材料及试验设备

福建省烟草农业科学研究所烟草种植土壤(含水率29.9%，密度1 698 kg/m3)，岛津AG-Xplus微机控制电子万能试验机，碳钢圆筒(内径55 mm，高250 mm)，碳钢钢板(表面光滑平整)，相机。

1.2 土壤堆积角试验及结果

采用圆筒提升法获得土壤颗粒堆，测量实际堆积角[9]。试验前取钢板水平放置于万能试验机平台，将圆筒垂直放置在钢板上并倒入土壤样品320 g，将圆筒上端连接万能试验机，随后启动万能试验机，将圆筒以1 000 mm/min匀速提升，用相机拍下土壤堆积情况，试验重复10次。

将土壤堆积图像文件导入Matlab，对图像进行灰度化、二值化处理，利用Canny算子进行边缘检测并取反，得到边界图像如图1所示；将边界图像文件导入Origin，利用图像数字化工具得到边界图像的坐标值[11-13]，对单测坐标值进行线性拟合[14]，得到散点图与拟合直线以及线性拟合方程式，线性拟合度值R2=0.997 7，如图2所示。对斜率取反三角函数并转化为角度，即可得到堆积角数值；最终，得到堆积角平均值为34.43°。

图1 Matlab处理后的边界图像Fig. 1 Matlab’s processed boundary image

图2 线性拟合图像与方程Fig. 2 Linearly fits the image and equation

2 植烟土壤参数标定仿真试验

2.1 仿真模型构建

福建烟草种植土壤为沙壤土，土壤颗粒之间存在一定的黏附关系，故选用Hertz-Mindlin with JKR模型，在该模型中，法向弹性接触力的实现是基于Johnson-Kendall-Roberts理论[15-16]。

(1)

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

γ——相互作用参数的表面能，N/m2；

E*——等效弹性模量，Pa；

δτ——切向重叠量，m；

R*——等效接触半径，m。

(2)

(3)

式中：E——接触颗粒的弹性模量，Pa；

i、j——颗粒i、j；

v——接触颗粒的泊松比；

R——接触颗粒的接触半径，m。

切向重叠量δτ与法向重叠量δn的关系

(4)

Hertz-Mindlin with JKR切向弹力计算公式

(5)

(6)

式中：Fτ——切向弹力，N；

Sτ——切向刚度，N/m；

G*——当量剪切模量，Pa。

Hertz-Mindlin with JKR切向阻尼力计算公式

(7)

(8)

式中：Fτd——切向阻尼力，N；

m*——等效质量，kg；

β——常数；

为了简化运算的同时保证试验结果的准确性，土壤颗粒采用半径为1.0 mm的球形颗粒。结合颗粒本征参数的测定及查阅相关文献[6, 9, 17-22]等，确定了仿真模型参数范围，如表1所示。

表1 土壤颗粒堆积角仿真模型参数Tab. 1 Soil particle accumulation angle simulation model parameters

2.2 参数标定方法

在EDEM前处理界面中设置颗粒参数，建立颗粒工厂；建立与实际圆筒、钢板所对应的模型。颗粒生成的初始速度设置为Z轴-2 m/s；仿真试验的颗粒工厂生成总质量为0.32 kg的颗粒，颗粒生成速度为0.64 kg/s，生成时间为0.5 s，总计生成0.32 kg，颗粒稳定后圆筒以0.166 m/s的速度匀速向上运动，圆筒完全移开且土壤堆积稳定后，获取堆积角图像。

利用Design-Expert软件设计二水平部分因子试验，在EDEM软件中展开土壤堆积角仿真试验，筛选显著影响堆积角的参数；针对筛选出的参数，进行最陡爬坡试验，得到最优值所在区域；利用Design-Expert软件设计Box-Behnken Design试验，开展土壤堆积角仿真试验。

2.3 因素显著性筛选试验及结果分析

以土壤—土壤、土壤—碳钢之间的接触参数为试验因素，以堆积角仿真试验值为评价指标，开展二水平部分因子试验设计。试验结束后，获取仿真堆积角图像，将图像文件导入Matlab，对图像进行灰度化、二值化处理，利用Canny算子进行边缘检测，对检测后的边界进行开运算及取反处理，将边界图像文件导入Origin，利用图像数字化工具得到边界图像的坐标值，对单测坐标值进行线性拟合，得到边界散点图及其拟合直线，对拟合直线斜率取反三角函数并转化为角度，即可得到堆积角数值，如表2所示。

对试验结果进行显著性分析，如表3所示。由表3可知，根据试验结果中6个因素对于堆积角的影响程度进行排序，在土壤堆积角试验中X2、X3与X5对堆积角影响显著，其余因素影响较小。因此在后续试验中，非显著因素取中间水平，即X1取0.375、X4取0.45、X6取0.2。

表2 二水平因子试验设计与结果Tab. 2 Design and results of the two horizontal factor tests

表3 参数显著性分析Tab. 3 Parameter significance analysis

2.4 最陡爬坡试验及结果分析

由表3可知，X3、X5对堆积角的效应为正值，故设置仿真试验参数值递增；X2对堆积角的效应为负值，故设置仿真参数递减。再次进行堆积角仿真试验，记录仿真堆积角试验值、相对误差，试验参数及结果如表4所示。

表4 最陡爬坡试验设计与结果Tab. 4 Design and results of the steepest climb test

在2个试验因素逐渐增大、1个试验因素逐渐减小的过程中堆积角随之增大，与实际物理试验的堆积角的相对误差呈先减后增的趋势，仿真试验中4号试验相对误差最小，为3.96%，故以4号试验中的各参数数值作为后期试验的中心点，3号和5号试验分别作为低水平和高水平进行后续的BBD响应曲面试验设计。

2.5 响应曲面试验及结果分析

根据筛选的试验结果，取X2、X3与X5三个基本接触参数，利用Design-Expert软件进行三因素三水平的响应曲面试验设计，中心水平设置5组重复，共进行17组土壤堆积角仿真试验，试验设计及仿真结果如表5所示。

表5 接触参数响应曲面试验设计与结果Tab. 5 Contact parameters respond to surface test design and results

应用Design-Expert软件对响应曲面试验结果进行拟合分析，选取二次全模型，模型拟合度值R2=0.922，方差分析结果如表6所示，响应曲面分析如图3所示。由表6知，二次全模型回归模型极其显著(P=0.003 9<0.01)，且失拟项不显著(P=0.097 8>0.05)，证明模型可靠，其中X2、X3和X5的平方项对堆积角的影响都显著。堆积角仿真值θ与土壤—土壤静摩擦系数X2、土壤—土壤动摩擦系数X3、土壤—碳钢静摩擦系数X5的二次多项式关系为

θ=1.488-94.872X2+195.071X3+141.233X5+28.719X2X3+156.317X2X5+481.657X3X5+30.968X22-385.320X32-791.550X52

(9)

以实际试验堆积角数值为目标值，利用Design-Expert软件的目标点求取工具求得，X2、X3、X5的值分别为0.593、0.353、0.267。以此数值再次进行堆积角仿真试验，重复5次，得出堆积角平均值θ=34.403°，相对误差δθ=1.08%，堆积角仿真结果与试验结果接近，说明标定的参数准确可靠。

(a) X2、X3与堆积角曲面图

(b)X2、X5与堆积角曲面图

(c)X3、X5与堆积角曲面图图3 响应曲面分析图Fig. 3 Response surface analysis diagram

3 土壤—触土部件相互作用离散元仿真试验及结果分析

3.1 触土部件三维模型

除秆清根触土部件由耙齿与螺旋叶片两部分构成，螺旋叶片由一对满面式叶片与一对带式叶片组成，如图4所示。工作时，触土部件相对于前进方向逆时针旋转，耙齿将土壤击散，螺旋叶片将土壤推向两侧，从而减小作业后的土壤不平度。根据烟草根茬直径、深度及植烟土垄宽、垄高，设计耙齿梁长度为320 mm，带式叶片外径为200 mm，内径为125 mm，满面式叶片半径为120 mm，耙齿长度为70 mm。

(a) 触土部件正视图 (b) 触土部件斜视图图4 触土部件三维模型Fig. 4 3D model of soil contact parts1.带式叶片 2.满面式叶片 3.耙齿 4.耙齿梁固定座 5.触土部件轴 6.耙齿梁

3.2 仿真模型

为了更直观地了解触土部件对土壤的作业效果，并得到触土部件所受扭矩，进行土壤—触土部件相互作用的离散元仿真。由Solidworks构建触土部件三维模型，导入EDEM软件中；在EDEM软件中建立3 000 mm×1 000 mm×700 mm大小的土槽模型，设置两个颗粒工厂，均生成半径为5 mm土壤颗粒；颗粒工厂1以质量生成颗粒，共计600 kg，生成速度为300 kg/s；颗粒工厂2以质量生成颗粒，共计146 kg，生成速度为146 kg/s，作为土垄；仿真时间步长设置为Rayleigh时间步长的10%[22]，触土部件行进速度为0.3 m/s，转速为6.5 rad/s。

3.3 仿真结果与分析

仿真试验前后作业效果对比如图5所示。作业后，中间的土垄被挖掘并被螺旋叶片推向两侧。

(a) 试验前

(b) 试验后图5 仿真试验Fig. 5 Simulated test

仿真试验中触土部件所受到的扭矩值在14.38～63.82 N·m间震荡，扭矩最大值为63.82 N·m，如图6所示。

图6 仿真试验扭矩值Fig. 6 Torque value of simulation test

4 土槽试验

4.1 试验条件

试验地点：福建农林大学现代农业装备福建省高校工程研究中心土槽试验室(土槽长×宽：7.0 m×2.4 m)，试验设备：土槽车(1.8 m×1.4 m，电机3 kW×2)，变频器(JF1800G；MITSUBISHI，A700)，扭矩传感器(CYT-302，精度为0.5%)，采集软件为M400数据采集管理软件。

将触土部件安装在土槽车，以变频器调节土槽车行进速度与触土部件转速，设定行进速度为0.3 m/s，转速为6.5 rad/s，传动方式采用齿数比为25∶60的同步带轮传动，传动效率取0.97。

4.2 试验方法及指标

调整土槽土壤状态使其接近于植烟沙壤土状态，使用变频器调整刀具转速及土槽车行进速度，使其与仿真试验一致，而后进行刀具试验。作业完成后，获取扭矩传感器数据与仿真试验扭矩数据进行对比分析。

由式(10)～式(12)计算触土部件所受扭矩。

(10)

(11)

故

(12)

式中：Tr——传感器所受扭矩，N·m；

T——触土部件所受扭矩，N·m；

FP——皮带受力，N；

r——小带轮半径，m；

R——大带轮半径，m；

η——传动效率。

4.3 结果分析

土槽试验前后作业效果对比如图7所示。作业后，中间的土垄被挖掘并被螺旋叶片推向两侧。对比图5与图7可知，仿真试验与土槽试验的作业效果一致。

土槽试验扭矩值与仿真试验扭矩值对比，如图12所示。土槽试验扭矩值图像与仿真试验扭矩值图像存在差异，这是由于土槽试验土壤并非均匀分布且人工起垄在不同位置的形状尺寸具有一定的差异性。土槽试验扭矩值在33.40～56.49 N·m间震荡，扭矩最大值为56.49 N·m，与仿真试验值(63.82 N·m)误差为11.48%，两者扭矩曲线趋势大致相同。对比仿真试验中15个扭矩波峰平均值(45.08 N·m)与土槽试验中15个扭矩波峰平均值(49.34 N·m)，误差为9.45%，误差在可接受范围内，表明仿真模型中土壤与实际土壤的物理力学特性基本一致，土壤模型可靠。

(a) 试验前

(b) 试验后图7 土槽试验Fig. 7 Soil tank test

图8 土槽试验与仿真试验扭矩值Fig. 8 Torque value of soil tank test and simulation test

5 结论

1) 通过显著性仿真试验、最陡爬坡仿真试验、响应曲面仿真试验，得出土壤—土壤碰撞恢复系数、土壤—碳钢碰撞恢复系数、土壤—碳钢动摩擦系数对堆积角影响较小，其取值分别为0.375、0.45、0.2；响应曲面仿真试验优化求解得到：土壤—土壤静摩擦系数、土壤—土壤动摩擦系数、土壤—碳钢静摩擦系数分别为0.593、0.353、0.267；利用标定的参数建立模型，进行5次堆积角仿真试验，均值θ=34.403°，相对误差δθ=1.08%，仿真结果与试验结果接近，说明标定的参数准确可靠。

2) 触土部件离散元仿真试验与土槽试验结果表明，触土部件所受扭矩仿真值与实测值曲线趋势大致相同，仿真试验中触土部件所受到的扭矩最大值(63.82 N·m)与实际试验值(56.49 N·m)的误差为11.48%。仿真试验中15个扭矩波峰平均值(45.08 N·m)与土槽试验中15个扭矩波峰平均值(49.34 N·m)的误差为9.45%，说明土壤模型准确可靠。

本研究以物理试验为参考，通过堆积角试验和土槽试验，确定了准确的南方植烟沙壤土离散元模型，为该类型土壤与触土部件相互作用的动力学研究提供了技术支持。