基于离散元仿真参数的微型薯物料测定

余参参，段宏兵,，蔡兴奎，徐涛 ，姚飞虎， 陈志惠，严福勇

1.华中农业大学工学院，武汉 430070; 2.农业农村部马铃薯生物学与生物技术重点实验室，武汉 430070

我国马铃薯种植面积与产量均居世界首位，但单产与世界平均水平仍有差距。实践证明，采用脱毒微型马铃薯(简称“微型薯”)代替切块薯增产效果显著，单位面积产量一般可增加30%～50%，高的可增产3～4倍[1-2]。通过对脱毒种薯的扩繁而形成优质商品种薯，是降低种薯成本的有效手段。现阶段适用于脱毒种薯的机械化播种机具较少，因而造成脱毒种薯的种植过程机械化水平较低。同时由于脱毒种薯块茎的形状差异较大，农业生产中急需微型薯播种与分级技术与装备[3]。随着数字仿真技术的发展，离散元应用在作物与作物之间、作物与机械设备之间的仿真研究渐多。研究表明，机械设备的离散元仿真分析结果与种子物理特性存在较大联系[4-5]。目前部分学者针对种子的物理特征参数和部分接触力学参数进行了标定，结果均表明离散元仿真的平均误差值较小，基本与试验一致，可以缩短研究周期[6-7]。由于微型薯颗粒个体差异较大，品种、含水率等对物料参数的测定结果影响巨大，因此，本研究以3种不同品种的微型薯为研究对象，以离散元软件EDEM中微型薯机械设备研究所需的相关仿真参数为目标，对微型薯的物理力学特征参数和接触力学参数进行测定分析，旨在为微型薯相关机械设备在EDEM中的仿真参数的设置提供参考。

1 材料与方法

1.1 试验材料

由于不同品种的微型薯差异较大，因此，试验用薯选取种植相对广泛且具有一定代表性的“大西洋”“中薯5号”和“华薯1号”，3种种薯均取自华中农业大学农业农村部马铃薯生物学与生物技术重点实验室种植基地。采样方法为：选取种植基地某处大棚内长势较为稳定的区域，随机选取半径2.0 m范围采挖3个品种的种薯，直至各收集500个种薯为止，现场收集的种薯经标记后进行测量并记录数据，测量后置于冷库储存备用。测试结果显示，3种种薯含水率分别为71.75%～83.65%、75.45%～82.58%和74.49%～80.48%。

1.2 微型薯特征参数测定

1)形态参数与质量的测量。对采集的“大西洋”“中薯5号”“华薯1号”样品各500个进行三轴尺寸、曲率半径和质量测量。质量测定采用“苦竹”电子秤，精度0.01 g；采用“艾瑞泽”数显游标卡尺测量种薯三轴尺寸，可测量三轴极限尺寸，精度0.01 mm；曲率半径选用“天目”多功能数显半径测试仪，其可针对不同的圆弧形面和测量精度选用不同跨度的量爪，精确读出种薯的外圆弧半径，精度0.005 mm。三轴尺寸和曲率半径的测量均分别选取长、宽、高3个方向进行(图1)。

X为测量时长度方向，Y为宽度方向，Z为高度方向，测量三轴尺寸时分别选取每个方向上的2个顶点进行测量，顶点选取如1、2所示；测量曲率半径时，将量爪的中心分别与所测种薯3个方向的中心对齐后读取数据。X is the length direction, Y is the width direction, and Z is the height direction. When measuring triaxial dimensions, select two vertices in each direction for measurement. The vertices are selected as sign 1 and 2; when measuring the radius of curvature, measure data is recorded after the center of the paw is aligned with the center of the three directions.

微型薯的三轴尺寸分布统计见图2，曲率半径与质量概率分布统计情况见图3。图2显示，3个品种的微型薯尺寸多分布在宽度10～28 mm、厚度10～26 mm、长度13～40 mm，均占比90%以上，其中宽度主要集中在12～18 mm，占比60%以上，厚度集中在12～18 mm，占比70%以上，长度集中在16～25 mm，占比50%左右。

图3显示，“大西洋”3个方向的曲率半径相差较大，符合程度较低，即3个方向的外圆弧半径差异较大，长度方向的尺寸大于宽度和高度，且高度方向和宽度方向尺寸较为接近，整体形状呈长椭圆形。“华薯1号”和“中薯5号”3个方向的曲率半径符合程度较高，即各个方向上的弯曲程度相近，三轴尺寸差别相对较小，相对来说外形偏圆，但“中薯5号”曲率半径较“华薯1号”符合程度更高，圆度更高，截面规则。“华薯1号”虽外形近圆，但其表面有小凹坑，截面形状有边角，近似椭圆形。不同品种的微型薯质量分布呈现出相似的变化规律，3个品种的种薯峰值均处于3～5 g，质量在10 g以下的种薯均占比90%以上；“中薯5号”质量分布范围最广，最小的不足1 g，最大的超过40 g。综上，同一品种的微型薯个体差异大，但基本符合正态分布；品种间形状差异较大，总体变化规律相似。

图2 三轴尺寸概率分布

A:“大西洋” “Daxiyang”; B:“中薯5号” “Zhongshu No.5”; C：“华薯1号” “Huashu No.1”; D：微型薯质量 Weight.

2)密度测量。由于微型薯表面较为光滑，体积较小，无疏松孔隙，液体难以浸润，可采用量筒法测量相关参数。随机选取3种微型薯各100个，选用量程1 500 mL的量筒和精度为0.01 g的电子秤，对3个品种的微型薯进行体积和质量的测量，试验重复5次，测量结果取平均值。 “大西洋”“中薯5号”“华薯1号”的密度分别为1 068、1 044、1 056 kg/m3。

3)剪切模量、泊松比的测定。本研究讨论在不超过屈服载荷下微型薯的弹性变形能力，采用万能试验机测量种薯破裂瞬间的相关数据，进而计算剪切模量和泊松比。试验仪器为瑞格尔生产制造的RGB-3005型万能试验机，可施加0～5 kN载荷。根据需要选用90 mm平板压头，通过该压缩仪器获得计算所需的横向变形量、纵向变形量、试验载荷力等数据参数，加载速度固定为10 mm/min。试验用种薯为测量物理力学特征参数时所用种薯，放置1个月后进行压缩试验测量(图4)。由于微型薯形状大小差异较大，考虑该差异可能对试验所测结果产生影响，试验时根据微型薯种植要求按质量小于3 g、3～5 g和6～10 g标准分组进行试验。

图4 压缩试验

剪切模量是剪切变形时剪切应力与剪切应变的比值[8]，其公式为：

(1)

式(1)中：E为弹性模量，MPa；μ为泊松比；G为剪切模量。由式(1)知，剪切模量可通过对弹性模量、泊松比的数值计算得到，泊松比、弹性模量参考文献[8]方法测定。

将3种微型薯按照质量差异分9组进行试验，每组测量50个，分别得9组剪切模量、泊松比和弹性模量。对各个品种的微型薯按质量分级段进行试验，结果发现弹性模量、泊松比与质量无关。如表1所示，仅列出了质量不一的“大西洋”的试验结果，试验结果统计仅取平均值。试验结果显示，同一品种间不同质量的微型薯，其弹性模量、泊松比及剪切模量无明显变化，但不同品种间差异较明显。

表1 微型薯的剪切模量和泊松比 Table 1 Shear modulus and Poisson’s ratio for potato minitubers

1.3 微型薯接触力学参数的测定

测定3种微型薯(“大西洋”“中薯5号”和“华薯1号”)分别与3种材料(钢板、ABS和各自品种微型薯薯群)的碰撞恢复系数、静摩擦系数和滚动摩擦系数。为消除微型薯的形状和质量影响，挑选大小相近、质量相当的种薯进行试验。当微型薯薯群作为试验材料时需要挑选大小相近，形状相同的若干较大微型薯粘结在一起。

1)碰撞恢复系数测定。如果物体与固定平面碰撞,则碰撞前后固定平面的速度均为0，此时碰撞恢复系数可表示为：

Cr=v1/v0

(2)

式(2)中：v0为碰撞前微型薯的速率，m/s；v1为碰撞后微型薯的速率，m/s。碰撞试验装置如图5所示，以“华薯1号”与3种材料的碰撞为例。

试验全程采用高速摄影拍摄，频率为1 000 fps/s。微型薯从同一高度(400 mm)由静止开始做自由落体运动，与碰撞材料接触后碰撞产生回弹力使之向上运动，整个运动过程能量守恒，碰撞恢复系数可由公式(3)得：

(3)

式(3)中：g为重力加速度，9.8 m/s2；h2为碰撞后回弹高度，m；h1为碰撞前下落高度，m。当碰撞材料为微型薯薯群时，下落的种薯碰撞到微型薯薯群上任一种薯为有效试验，种薯下落到薯与薯之间的空隙时无效，不计入试验。每个品种的微型薯与3种碰撞材料均取20次有效试验，每次试验重复5次。试验结果显示：“大西洋”“中薯5号”“华薯1号”与钢板、ABS和各自品种微型薯薯群间的碰撞恢复系数分别为：0.577、0.482、0.449；0.564、0.480、0.460；0.562、0.459、0.440。

A:“华薯1号”与ABS碰撞 Collision between “Huashu No.1” and ABS; B:“华薯1号”与“华薯1号”碰撞 Collision between “Huashu No.1” and “Huashu No.1”; C:“华薯1号”与钢板碰撞 Collision between “Huashu No.1” and steel plate;1.碰撞材料 Collision material；2.刻度尺 Scale；3.碰撞的微型薯 Collision of minitubers; 4.试验台架 Experiment shelf.

2)静摩擦系数测定。微型薯置于倾斜材料上，重力分解为与斜面平行的力F1和与斜面垂直的力F2，最大静摩擦力f与F1大小相等，方向相反，即F1=f=mgsinα，此时斜面的倾角与最大静摩擦存在以下关系：

μ=f/F2=mgsinα/(mgcosα)=tanα

(4)

因此，可采用静摩擦测量装置进行测量。测量时将3种材料依次置于斜板上固定，微型薯(为消除滚动摩擦的影响，挑选各自品种代表且大小相近的4个种薯粘连在一起为1组)放置在固定好的材料上，放置位置相同，均为斜板300 mm处。缓慢调节斜板角度，直至微型薯和材料存在滑动趋势时停止角度调节，记录此时底板与材料之间高度h，由α=arctan(h/300)计算得到倾斜角度，即为静摩擦系数。每个品种的种薯选取10组进行试验，每组试验重复20次取平均值。结果显示：“大西洋”“中薯5号”“华薯1号”与钢板、ABS和各自品种微型薯薯群间的静摩擦系数分别为：0.469、0.493、0.434;0.386、0.390、0.356；0.426、0.411、0.38。

3)滚动摩擦系数测定。当点或线接触的两物体在法向和转矩作用下产生相对滚动时，其表面会产生阻碍相对运动趋势的形变，可用滚动摩擦力衡量[9]。滚动摩擦系数的测定可采用斜面滚动法。微型薯的滚动距离和滚动时间经高速摄像记录后，由式(5)计算相应的滚动摩擦系数。试验装置同静摩擦试验装置，材料固定在斜板上，微型薯放置在材料上同一位置由静止释放，以0 m/s开始在材料上运动。高速摄影记录滚动过程，利用帧数变化来计算滚动时间。设置相机帧数200 fps/s，从微型薯开始滚动的帧数算起，直至在斜面的滚动距离为200 mm时的帧数为止，帧数差即为滚动时间，滚动距离200 mm。微型薯滚动的过程中满足：

(5)

式(5)中：μ为滚动摩擦系数；θ为斜面夹角，(°)；S为滚动距离，m；g为重力加速度，9.8 m/s2；t为滚动时间，s。

微型薯的弹性模量较钢板、ABS小，微型薯运动时可忽略微观变形和弹性滞后，只有可能产生些许滑动影响试验结果。当倾斜角度较大时，滑动影响越大。为了降低这种影响，对3种微型薯进行预滚动，选取3种微型薯都能滚动的最小角度(16°)进行滚动试验；每个品种进行30次重复滚动试验取平均值。结果显示：“大西洋”“中薯5号”和“华薯1号”与钢板、ABS间的滚动摩擦系数分别为： 0.119、0.120；0.067、0.089；0.076、0.082。

1.4 仿真试验

1)仿真颗粒模型设置。根据各微型薯的物理力学特征参数建立微型薯颗粒模型，各品种按比例产生不同尺寸的种薯颗粒(图6)。

2)堆积角仿真参数。由于微型薯之间的滚动摩擦系数影响因素太多，难以直接测量，而堆积角与物料的摩擦特性有很大关系[10-11]，因此，目前通用的方法是采用离散元软件EDEM建立测量堆积角模型，并输入微型薯间的静摩擦系数、微型薯间的滚动摩擦系数以及微型薯和ABS间滚动摩擦系数进行仿真试验，改变EDEM软件中的各摩擦系数，得到不同的物料堆积模型，测得对应物料堆积模型的堆积角。选取当仿真结果的堆积角与实际试验中的堆积角近似相等时的滚动摩擦系数作为物料间的实际滚动摩擦系数。

由于微型薯间滚动摩擦系数较小，故前期仿真试验优先确定其范围为0.01～0.11，仿真时按照优先确定的范围将其等分设置为3个水平；微型薯间静摩擦系数、微型薯与ABS间滚动摩擦系数的设置前述试验测得的参数(表2)为基础值上下浮动各0.05。以堆积角θ为目标值，在Design Expert软件中分别对3种微型薯开展三因素三水平响应面设计，每组试验17组，3组试验共计51组。每组试验重复3次取平均值。其余各仿真参数设置如表2所示。

A:“大西洋” “Daxiyang”; B:“华薯1号” “Huashu No. 1”; C:“中薯5号” “Zhongshu No.5”.

表2 堆积角仿真参数设置 Table 2 Simulation parameter setting for stacking

3)堆积角试验。堆积角试验采用抽板法(图7)。仿真接触模型选择无滑移Hertz-Mindlin，材料为ABS，在软件Geometry页面中创建1个200 mm×150 mm×150 mm无盖长方体，一个虚拟的四边形盖板作为种薯下落到长方体内的基准平面。颗粒工厂设为动态，种薯的总质量3.0 kg，根据种薯的实际几何尺寸设置，种薯按照一定的数量比例以500个/s的速度填充长方体内部，填充时间为5 s。待种薯填满盒内后，创建1个四边形置于长方体底面充当承接掉落种薯的地板，设置长方体的一侧面以0.1 m/s的速度向上运动，使种薯能以较缓的速度从盒中自由散落在地板上，仿真时间设置为10 s，待种薯静止后测量堆积角度。真实试验与软件仿真模型相同，对3种微型薯进行5次抽板法的重复试验，得平均值θ，见式(6)～(8)。

A:软件仿真 Simulation; B:真实试验 Actual test.

(6)

(7)

(8)

2 结果与分析

2.1 响应面试验

响应面试验设计结果见表3(以“大西洋”为例)，对试验结果进行方差分析，发现自变量中“大西洋”间滚动摩擦系数(B1)极显著(P<0.000 1)，失拟项为0.038，该方程拟合较好，比较可靠，结果表明因素B1对堆积角的结果影响较大，“大西洋”间静摩擦系数(A1)、“大西洋”与ABS间滚动摩擦系数(C1)影响较小。为了得到合理的堆积角，仿真与试验结果均采用MATLAB对图片进行处理[12]，提取边界后进行直线拟合，获取堆积角的真实角度(图8)。

表3 “大西洋”微型薯仿真试验设计及结果 Table 3 Simulation test design and results of “Daxiyang”

A:堆积角 Stacking angle; B:提取边界 Extraction boundary; C:直线拟合 Straight line fitting.

使用软件的优化功能分析确定符合的参数，以自变量A1、B1、C1的变化范围为基础，将实际堆积角作为响应值(目标)对回归模型求解(式6)。式(6)中，A、B、C分别为微型薯间静摩擦系数、微型薯间滚动摩擦系数、微型薯与ABS间滚动摩擦系数；1、2、3分别为“大西洋”“中薯5号”和“华薯1号”，如A1为“大西洋”间静摩擦系数。确定 “大西洋”“中薯5号”和“华薯1号”与各自品种间的滚动摩擦系数分别为0.031、0.020和0.024。

2.2 试验结果验证

在几何体模块新建一个仿真模型(图9A)，创建一个500 mm×500 mm的四边形底面置于圆筒底部，一个无底无盖圆筒200 mm(半径)×350 mm(高)，将前述试验测得的物理力学特征参数、接触力学参数和仿真得到的微型薯间的滚动摩擦系数输入EDEM软件参数设置中，其余相关设置参见表2。设置微型薯颗粒的总质量为2.0 kg，仿真时间10 s。待种薯填满圆筒且稳定后，设置圆筒以0.1 m/s的速度向上做线性运动，直至仿真结束，测得“大西洋”“中薯5号”“华薯1号”的仿真堆积角分别为23.377°、19.824°和20.913°。

实际试验(图9B)与仿真试验方法一致，进行5次重复试验后得“大西洋”“中薯5号”“华薯1号”的堆积角分别为23.778°、19.423°和21.314°，与仿真的堆积角试验结果相比，相对误差分别为1.69%、2.06%、1.88%，表明采用此方法可得到较为相似的结果，与实际情况相符，即微型薯离散元颗粒参数可为EDEM模拟仿真微型薯颗粒与机械设备优化提供参考。

A:仿真试验 Simulation test; B 真实试验 Actual test.

3 讨 论

本研究测量了湖北地区常见3种微型薯(“大西洋”“中薯5号”“华薯1号”)的物理力学特征参数，包括其基本物理特征参数(三轴尺寸、质量、曲率半径)、密度、泊松比和剪切模量。测试结果表明，密度、剪切模量和泊松比与微型薯的基本物理特征参数相关性较小，但基本物理特征参数差异较大，品种差异较大。在设计、仿真计算时可以根据用户需求选择所需参数。3种微型薯和各自品种微型薯薯群之间的接触力学参数较钢板、ABS均低，3种微型薯和ABS的静、动摩擦系数均高于钢板，但碰撞恢复系数则呈相反趋势，ABS均低于钢板。“大西洋”的接触力学参数均高于“中薯5号”和“华薯1号”，3种材料中“大西洋”与钢板的接触力学参数最大，与ABS和微型薯薯群之间的接触力学参数依次降低，表明其弹性恢复性能最好，摩擦力较大，相对不容易损伤，但流动性较差，机械化播种性能差。“中薯5号”和“华薯1号”的接触力学参数差异较小，“中薯5号”略大于“华薯1号”。不同品种对接触力学参数影响大，在设计、仿真时要充分考虑品种差异，不能忽略。采用响应曲面法设计试验，以静摩擦系数、滚动摩擦系数为试验变量，堆积角为试验指标，根据试验因素与变量的二阶回归方程进行方差分析，结果表明，各品种微型薯之间的滚动摩擦系数差异极显著。利用Design Expert软件中的优化模块将试验测得的实际堆积角作为响应值对得到的回归模型求解，确定了“大西洋”“中薯5号”“华薯1号”与各自品种间的滚动摩擦系数分别为0.031、0.020和0.024。通过圆筒提升试验对仿真得到的滚动摩擦系数进行验证，结果显示仿真结果与试验结果的差异均小于3%，说明本研究所测定的微型薯颗粒实际参数可作为EDEM软件中仿真参数使用，仿真结果误差较小。